计算固体力学网格

计算固体力学网格第一页，共四十四页，2022年，8月28日第4章Lagrangian网格

引言UL控制方程，弱形式UL有限元离散编制程序旋转公式TL格式，弱形式，有限元半离散化第二页，共四十四页，2022年，8月28日1引言

在Lagrangian网格中，节点和单元随着材料移动，边界和接触面与单元边缘保持一致，处理较为简单。积分点也随着材料移动，本构方程总是在相同材料点处赋值，这对于历史相关材料是有利的。基于这些原因，在固体力学中广泛地应用Lagrangian网格。UL格式，Eulerian（空间）坐标和Cauchy应力；TL格式，名义应力，PK2应力，Green应变张量。第三页，共四十四页，2022年，8月28日2UL控制方程

弱形式考虑一个物体，占有域Ω，边界为Г。连续体力学行为的控制方程是：1质量（或物质）守恒，标量方程；2线动量和角动量守恒，张量方程，包含n个偏微分方程（n－维数）；3能量守恒，通常称作热力学第一定律，标量方程；4本构方程，应力－应变或应变率的关系，对称张量；5应变－位移方程。第四页，共四十四页，2022年，8月28日2UL控制方程

弱形式第五页，共四十四页，2022年，8月28日边界条件：在二维问题中，面力或速度的每个分量都必须预先指定在整个边界上；但是，面力和速度的同一个分量不能指定在边界上同一点处。其分量可以指定在不同于总体坐标系的局部坐标系上。速度边界条件等价于位移边界条件：如果给定了位移，可以通过时间微分得到速度；给定了速度，可以通过时间积分得到位移。2UL控制方程

弱形式初始条件：可以是速度和应力，或者是位移和速度。第一组初始条件更适合于大多数工程问题，因为确定一个物体的初始位移通常是很困难的。初始应力通常为已知的残余应力，有时候可以测量或者通过平衡解答估算。例如，当一个钢件经过铸锭成型后确定其位移几乎是不可能的。对于在工程部件中的残余应力场，经常能够给出较准确的估计。类似地，在埋置管道中，靠近管道周围的土壤或岩石的初始位移的概念是毫无意义的，而初始应力场可以通过平衡分析估计出来。因此，以应力形式的初始条件更加实用。第六页，共四十四页，2022年，8月28日虚功率原理是动量方程，面力边界条件和内部力连续性条件的弱形式。微分方程的积分形式一般称为弱形式。强形式或广义动量平衡，包括动量方程，力边界条件和力连续性条件。微分方程一般称为强形式。2UL控制方程

弱形式第七页，共四十四页，2022年，8月28日3UL有限元离散有限元近似

在有限元方法中，运动近似地表示为当前构形中的节点坐标

小写的下标表示分量，如三维大写的下标表示节点值

默认对重复的指标求和；在小写指标的情况下，对空间的维数进行求和，而在大写指标的情况下，对节点的编号进行求和。在求和中的节点数目取决于所考虑的域：当考虑整个域时，对整个域中的所有节点求和；当考虑一个单元时，对这个单元的所有节点求和。第八页，共四十四页，2022年，8月28日3UL有限元离散有限元近似

当一个节点具有初始位置有

节点J总是对应于相同的材料点XJ，在L网格中，节点总是和材料点保持一致定义节点位移：

位移场：

取位移的材料时间导数得到速度：

速度是位移的材料时间导数，即当材料坐标固定，对时间求偏导数。由于形状函数不随时间改变，因此速度是由相同形状函数给出的。节点位移上面的点表示普通导数，因为它仅是时间的函数。第九页，共四十四页，2022年，8月28日3UL有限元离散有限元近似

加速度是速度的材料时间导数

速度梯度为

变形率给出为

变分函数或变量不是时间的函数，因此将变分函数近似为

虚拟节点速度

第十页，共四十四页，2022年，8月28日3UL有限元离散有限元近似

作为构造离散有限元方程的第一步，将变分函数代入虚功率原理中，得到

利用除以外的节点上虚节点速度的任意性，则动量方程的弱形式为

在任何指定速度的地方，虚速度必须为零。

离散运动（动量）方程为第十一页，共四十四页，2022年，8月28日3UL有限元离散离散运动（动量）方程内部节点力

外部节点力

惯性节点力

内部节点力代表着物体的应力。这些表达式既可以应用于整体网格，也可以应用于任意单元或单元集。这些表达式包含形状函数对应于空间坐标的导数和在当前构形上的积分。在非线性有限元方法中，对于更新的Lagrangian网格，这是一个关键的方程；它也应用于Eulerian和ALE网格。离散方程

第十二页，共四十四页，2022年，8月28日3UL有限元离散是关于节点速度的个常微分方程系统。

半离散运动（动量）方程是不受约束的节点速度分量的数目,

称作自由度的数目。

为了完成这个方程系统，要附加上单元积分点处的本构方程和以节点速度形式表示的变形率。

积分点与材料点是一致的。离散方程

在网格中nQ个积分点表示为为应力张量的独立分量数目，在二维平面应力问题中，由于应力张量对称，;在三维问题中，

。第十三页，共四十四页，2022年，8月28日3UL有限元离散通过任何常微分方程的积分方法，如Runge－Kutta法或中心差分法，可以对这个常微分方程系统进行时间积分；见第6章。离散方程

式中这是一个标准的初值问题，包括含有速度和应力的一阶常微分方程。消去变形率，所有未知量的个数就变为半离散运动（动量）方程为关于时间的常微分方程第十四页，共四十四页，2022年，8月28日3UL有限元离散对于平衡问题，加速度为零，控制方程成为

这是真正意义的离散平衡方程。如果本构方程是率无关的，那么离散平衡方程是关于应力和节点位移的非线性代数方程组。对于率相关材料，为了获得非线性代数方程组，任何率形式都必须在时间上离散。对于线性问题，控制方程也可以写成KU＝F矩阵位移法的刚度方程形式。离散方程

第十五页，共四十四页，2022年，8月28日3UL有限元离散单元坐标

通常建立有限元是采用以母单元坐标形式表示形状函数，简称为单元坐标。单元坐标的例子有三角形坐标和等参坐标。下面说明以单元坐标形式表示形状函数的用法。证明在L网格中，单元坐标可以考虑为材料坐标的另一种形式。这样，在L网格中将形状函数表示为单元坐标的形式，在本质上等价于把它们表示为材料坐标的形式。一个单元的三个域：母单元域□当前单元域3.初始（参考）单元域第十六页，共四十四页，2022年，8月28日3UL有限元离散单元坐标

相关的映射

通过映射的合成描述每一单元的运动第十七页，共四十四页，2022年，8月28日3UL有限元离散单元坐标

运动近似给出为

形状函数仅是母单元坐标的函数；运动的时间相关性完全反映在节点坐标上。上式代表了在单元的母域和当前构形之间的一个时间相关映射。在t＝0时写出这个映射，得到

在一个L单元中，材料坐标和单元坐标之间的映射是时间不变的。如果这个映射是一对一的，则在L网格中可以将单元坐标看作是材料坐标的代用品，因为在一个单元中的每一材料点具有唯一的单元坐标编号。为了在Ω0中在单元坐标和材料坐标之间建立唯一的对应关系，单元数目必须成为编号的一部分。如果单元坐标不能代替材料坐标，则网格不是L格式(见第7章)。事实上，应用初始坐标X作为材料坐标主要源于解析过程；在有限元方法中，应用单元坐标作为材料编号是更自然的。第十八页，共四十四页，2022年，8月28日3UL有限元离散单元坐标

单元坐标是时间不变的，可以将位移、速度和加速度表示为形状函数的形式：第十九页，共四十四页，2022年，8月28日4编制程序在有限元方程的程序编制中，通常采用两种方法：

1．将指标表示直接处理为矩阵方程。2.使用Voigt标记，将矩形应力和应变矩阵转换为列矩阵。每种方法都有其优点。在框4.3中总结了两种形式的离散方程。

从指标表示到矩阵形式的转换是比较任意的，并取决于个人的偏爱。在大多数情况下，将单指标的变量解释为列矩阵；当解释为行矩阵时，其过程就会有所不同。

在Voigt标记中，将应力和变形率表示为列向量的形式。

第二十页，共四十四页，2022年，8月28日4编制程序第二十一页，共四十四页，2022年，8月28日4编制程序数值积分

节点力、质量矩阵和其它单元矩阵的积分不是由解析计算的，而是应用数值解答，称为数值积分。最广泛应用的是Gauss积分

式中nQ个积分点的权重wQ和坐标值ξQ有表可查；见附录3。指定方程在母单元域上进行积分，其积分区间为[－1，1]。

一个二维单元的Gauss积分为第二十二页，共四十四页，2022年，8月28日在非线性分析中，采用积分点数的规则一般基于在线性分析中的相同规则；对于一个规则的单元，积分点数目的选择是能恰好积分内部节点力。一个单元的规则形式，是指仅通过母单元的拉伸而不是剪切能得到的形式；例如，二维等参单元的一个矩形。对于一个4节点四边形单元，如何选择内部节点力的积分点数目？由于速度是双线性的，单元中的变形率和B矩阵是线性的。如果应力是与变形率线性相关，那么它将在单元内线性变化。内部节点力的被积函数是近似为二次的，因为它是B矩阵和应力的乘积。在Gauss积分中，对于一个二次函数的精确求解在每一方向上需要两个积分点，所以对于线性材料，需要2×2个点的积分得到内部节点力的精确解。对于线性本构方程的积分，几乎得到内部节点力精确解的积分公式，称为完全积分。4编制程序完全积分

第二十三页，共四十四页，2022年，8月28日4编制程序局部减缩积分

对于完全不可压缩或接近不可压缩的材料，运动必须是等体积的

0J=1式中K是体积模量，

是剪切模量。在任意的等体积运动中，单元的整个体积将保持常数，即在整个单元中的运动必须是等体积的，否则，当K是一个非常大的数时(一个接近于不可压缩材料)，任何非零体积应变将吸收几乎全部的能量。内部节点力的完全积分可能引起单元的自锁，即出现很小的位移而不收敛或收敛得非常慢。考虑一种线性材料，如果将线性弹性应变能分解为静水和偏量部分，可以写为第二十四页，共四十四页，2022年，8月28日4编制程序局部减缩积分

为了克服这个困难，最容易的方法是使用局部减缩积分。在局部减缩积分中，压力为不完全积分，而应力矩阵的其余部分为完全积分。为此，将应力张量分解为偏斜部分和静水部分

体积自锁源于单元没有能力准确地表示一个等体积运动。为了消除自锁，必须设计应变场，这样在假设的应变场中整个单元的膨胀为零：为了避免自锁，对于任意保持单元体积的速度场，整个单元的应变场必须是等体积的。对于等体积运动的四边形单元，为了克服沙漏模式，在整个单元中膨胀必须为零。

第二十五页，共四十四页，2022年，8月28日4编制程序局部减缩积分

在局部减缩积分中，压力为不完全积分，以此确定与外力平衡的应力场。如杂交单元第二十六页，共四十四页，2022年，8月28日4编制程序局部减缩积分

注意到膨胀部分和偏斜部分是彼此正交的，因此内部虚功率为将变形率表示为形状函数的形式，膨胀和偏斜部分的被积函数分别为局部减缩积分包含在偏斜功率上的完全积分和在膨胀功率上的减缩积分。对于一个4节点四边形单元的局部减缩积分为关于内力的局部减缩积分表达式为第二十七页，共四十四页，2022年，8月28日4编制程序局部减缩积分

只有四边形和六面体单元才能采用减缩积分；而所有楔形、四面体和三角形实体单元采用完全积分。

关于减缩积分的单个积分点是母单元的质心。偏斜部分是通过完全积分方式积分的。这种方法类似于应用于不可压缩材料线性分析的算法。关于其它单元的局部减缩积分算法，可以采取对于线性有限元类似的修正方法建立起来。单元力和矩阵转换

例题：三角形3节点单元

第二十八页，共四十四页，2022年，8月28日4编制程序例题四边形单元和其它二维等参单元

建立二维等参单元的变形梯度、变形率、节点力和质量矩阵的表达式。对于4节点四边形单元给出详细的表达式。以矩阵形式给出内部节点力的表达式。

是母单元的节点坐标

形状函数和节点变量

节点速度

节点坐标

第二十九页，共四十四页，2022年，8月28日变形率和内部节点力

4编制程序对于形状函数公式，映射到空间坐标不是可逆的。不能直接将表达成x和y的形式，所以形状函数的导数是通过隐式微分计算出来

当前构形相对于单元坐标的Jacobian矩阵给出为对于4节点四边形单元，上式成为第三十页，共四十四页，2022年，8月28日4编制程序Jacobian矩阵考虑平衡方程

定义残差

通过Newton迭代求解和线性化

内力Jacobian－切线刚度外力Jacobian－载荷刚度整体Jacobian为二者之差ABAQUS的UMAT就是求内力的Jacobian矩阵－切线刚度第三十一页，共四十四页，2022年，8月28日Jacobian矩阵是时间的函数，的逆给出为

4编制程序4节点四边形单元采用单元坐标表示的形状函数的梯度为

采用空间坐标表示的形状函数的梯度为出现在分母中第三十二页，共四十四页，2022年，8月28日4编制程序速度梯度给出为

内部节点力

上式适用于任意的二维等参单元。因为是关于单元坐标的有理函数，所以不能对上式进行解析积分，通常使用数值积分。对于4节点四边形单元，2×2Gauss积分为完全积分。但是在平面应变问题中，对于不可压缩或几乎不可压缩的材料，完全积分会出现单元自锁，因此必须应用局部减缩积分。对于4节点四边形单元，沿着每条边的位移是线性的。所以，其外部节点力与3节点三角形单元外部节点力是一致的。出现在分母中，被积函数第三十三页，共四十四页，2022年，8月28日4编制程序一致质量矩阵

根据Lobatto积分，使积分点与节点重合，或者将单元的整个质量平均分配到4个节点上可以得到一个集中的对角质量矩阵，分配到4个节点上给出集中质量矩阵

第三十四页，共四十四页，2022年，8月28日例题主－从连接线

4编制程序连接线经常用于连接采用不同单元尺度的网格部分，因为它们比使用三角形或四面体单元连接不同尺度的单元更方便。通过约束从属节点的运动，使其与连接主控节点的附近边界的场一致，保证跨过连接线的运动的连续性。根据转换规则，导出节点力和质量矩阵。第三十五页，共四十四页，2022年，8月28日例题主－从连接线

4编制程序通过运动约束给出从属节点的速度，使沿连接线两侧的速度保持协调，即C0。这个约束可以表达为节点速度的一个线性关系，可以写成为连接后模型的节点力为

矩阵A是从线性约束得到的，“^”表示在两侧连接到一起之前的分离模型的速度

主控节点力是分离模型的主控节点力和转换的从属节点力的和。这些公式对于外部和内部节点力都适用

一致质量矩阵为

第三十六页，共四十四页，2022年，8月28日4编制程序4节点四边形单元，沿任何边界上的速度是线性的。从属节点3和5与主控节点1和2重合，而从属节点4与节点1的距离为节点速度可以写成为节点力则给出为主控节点1的力是节点力两个分量的转换是一致的；这种转换对于内部和外部节点力都适用。如果两条边只是在法向相连接，则需要在节点上建立局部坐标系。通过公式，联系节点力的法向分量，而切向分量保持各自独立。例题主－从连接线

第三十七页，共四十四页，2022年，8月28日5旋转公式在结构单元中，如杆、梁和壳，处理固定坐标系是很棘手的。例如，考虑一个旋转杆件。最初，唯一的非零应力是σx，而σy等于零。当杆件旋转时，以应力张量整体分量的形式表示单轴应力的状态就麻烦了。克服这种困难的一般途径是在杆中嵌入一个随杆旋转的坐标系－旋转坐标系考虑一个坐标系始终连接着节点1和2，在单轴应力状态下

杆的变形率类似地应用分量描述

第三十八页，共四十四页，2022年，8月28日在旋转格式中考虑的关键问题是材料旋转的定义，有两种方法：1在每一个积分点嵌入一个坐标系，在某种意义上它随着材料旋转。对于任意大的应变和旋转都有效。在大多数情况下，应用极分解原理定义旋转。然而，当材料的指定方向有一个较大的刚度并需要准确地表现出来时，由极分解提供的旋转在Cartesian坐标系下未必能提供最好的旋转；见第9章。5旋转公式对于一些单元，如杆或常应变三角形单元，整个单元的刚体转动是相同的，那么在单元中嵌入一个单一坐标系就足够了。对于高阶单元，如果应变比较小，不需要坐标系准确地随材料旋转。例如定义旋转坐标系，使它与单元的一边重合。如果与嵌入坐标系相关的旋转为θ阶，则在应变中的误差具有θ2阶。这样，只要θ2与应变相比很小，单个嵌入坐标系是合适的。称作小应变、大转动问题。2在一个单元中嵌入一个坐标系并随着单元旋转。第三十九页，