清华大学计算固体力学第四次课件_Lagrangian网格.ppt

1、非线性有限元第4章 Lagrangian网格,计算固体力学,第4章 Lagrangian网格,引言 UL控制方程，弱形式 UL有限元离散 编制程序 旋转公式 TL格式，弱形式，有限元半离散化,1 引言,在Lagrangian网格中，节点和单元随着材料移动，边界和接触面与单元边缘保持一致，处理较为简单。积分点也随着材料移动，本构方程总是在相同材料点处赋值，这对于历史相关材料是有利的。基于这些原因，在固体力学中广泛地应用Lagrangian网格。 UL格式，Eulerian（空间）坐标和Cauchy应力； TL格式，名义应力，PK2应力，Green应变张量。,2 UL控制方程 弱形式,考虑一个物体

2、，占有域，边界为。 连续体力学行为的控制方程是： 1 质量（或物质）守恒，标量方程； 2 线动量和角动量守恒，张量方程，包含n个偏微分方程（n维数）； 3 能量守恒，通常称作热力学第一定律，标量方程； 4 本构方程，应力应变或应变率的关系，对称张量； 5 应变位移方程。,2 UL控制方程 弱形式,边界条件：在二维问题中，面力或速度的每个分量都必须预先指定在整个边界上；但是，面力和速度的同一个分量不能指定在边界上同一点处。其分量可以指定在不同于总体坐标系的局部坐标系上。 速度边界条件等价于位移边界条件：如果给定了位移，可以通过时间微分得到速度；给定了速度，可以通过时间积分得到位移。,2 UL控制

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

4、程，力边界条件和力连续性条件。微分方程一般称为强形式 。,2 UL控制方程 弱形式,3 UL有限元离散,有限元近似,在有限元方法中，运动,近似地表示为,当前构形中的节点坐标,小写的下标表示分量，如三维 大写的下标表示节点值,默认对重复的指标求和；在小写指标的情况下，对空间的维数进行求和，而在大写指标的情况下，对节点的编号进行求和。 在求和中的节点数目取决于所考虑的域：当考虑整个域时，对整个域中的所有节点求和；当考虑一个单元时，对这个单元的所有节点求和。,3 UL有限元离散,有限元近似,当一个节点具有初始位置,有,节点J总是对应于相同的材料点XJ ，在L网格中，节点总是和材料点保持一致,定义节点

5、位移：,位移场：,取位移的材料时间导数得到速度：,速度是位移的材料时间导数，即当材料坐标固定，对时间求偏导数。由于形状函数不随时间改变，因此速度是由相同形状函数给出的。节点位移上面的点表示普通导数，因为它仅是时间的函数。,3 UL有限元离散,有限元近似,加速度是速度的材料时间导数,速度梯度为,变形率给出为,变分函数或变量不是时间的函数，因此将变分函数近似为,虚拟节点速度,3 UL有限元离散,有限元近似,作为构造离散有限元方程的第一步，将变分函数代入虚功率原理中，得到,利用除,以外的节点上虚节点速度的任意性，则动量方程的弱形式为,在任何指定速度的地方，虚速度必须为零。,离散运动（动量）方程为,3

6、 UL有限元离散,离散运动（动量）方程,内部节点力,外部节点力,惯性节点力,内部节点力代表着物体的应力。这些表达式既可以应用于整体网格，也可以应用于任意单元或单元集。这些表达式包含形状函数对应于空间坐标的导数和在当前构形上的积分。在非线性有限元方法中，对于更新的Lagrangian网格，这是一个关键的方程；它也应用于Eulerian和ALE网格。,离散方程,3 UL有限元离散,是关于节点速度的,个常微分方程系统。,半离散运动（动量）方程,是不受约束的节点速度分量的数目,称作自由度的数目。,为了完成这个方程系统，要附加上单元积分点处的本构方程 和以节点速度形式表示的变形率。,积分点与材料点是一致

7、的。,离散方程,在网格中nQ个积分点表示为,为应力张量的独立分量数目，,在二维平面应力问题中，由于应力张量对称， ;,在三维问题中， 。,3 UL有限元离散,通过任何常微分方程的积分方法，如RungeKutta法或中心差分法，可以对这个常微分方程系统进行时间积分；见第6章。,离散方程,式中,这是一个标准的初值问题，包括含有速度,和应力,的一阶常微分方程。消去变形率，所有未知量的个数就变为,半离散运动（动量）方程为关于时间的常微分方程,3 UL有限元离散,对于平衡问题，加速度为零，控制方程成为,这是真正意义的离散平衡方程。如果本构方程是率无关的，那么离散平衡方程是关于应力和节点位移的非线性代数方

8、程组。对于率相关材料，为了获得非线性代数方程组，任何率形式都必须在时间上离散。,对于线性问题，控制方程也可以写成 KUF 矩阵位移法的刚度方程形式。,离散方程,3 UL有限元离散,单元坐标,通常建立有限元是采用以母单元坐标形式表示形状函数，简称为单元坐标。单元坐标的例子有三角形坐标和等参坐标。 下面说明以单元坐标形式表示形状函数的用法。证明在L网格中，单元坐标可以考虑为材料坐标的另一种形式。这样，在L网格中将形状函数表示为单元坐标的形式，在本质上等价于把它们表示为材料坐标的形式。,一个单元的三个域： 母单元域 当前单元域,3. 初始（参考）单元域,3 UL有限元离散,单元坐标,相关的映射,通过

9、映射的合成描述每一单元的运动,3 UL有限元离散,单元坐标,运动近似给出为,形状函数仅是母单元坐标的函数；运动的时间相关性完全反映在节点坐标上。上式代表了在单元的母域和当前构形之间的一个时间相关映射。,在t0时写出这个映射，得到,在一个L单元中，材料坐标和单元坐标之间的映射是时间不变的。如果这个映射是一对一的，则在L网格中可以将单元坐标看作是材料坐标的代用品，因为在一个单元中的每一材料点具有唯一的单元坐标编号。为了在0中在单元坐标和材料坐标之间建立唯一的对应关系，单元数目必须成为编号的一部分。如果单元坐标不能代替材料坐标，则网格不是L格式(见第7章)。事实上，应用初始坐标X作为材料坐标主要源于

10、解析过程；在有限元方法中，应用单元坐标作为材料编号是更自然的。,3 UL有限元离散,单元坐标,单元坐标是时间不变的，可以将位移、速度和加速度表示为形状函数的形式：,4 编制程序,在有限元方程的程序编制中，通常采用两种方法： 1将指标表示直接处理为矩阵方程。 2. 使用Voigt标记，将矩形应力和应变矩阵转换为列矩阵。 每种方法都有其优点。在框4.3中总结了两种形式的离散方程。,从指标表示到矩阵形式的转换是比较任意的，并取决于个人的偏爱。在大多数情况下，将单指标的变量解释为列矩阵；当解释为行矩阵时，其过程就会有所不同。,在 Voigt 标记中，将应力和变形率表示为列向量的形式。,4 编制程序,4

11、 编制程序,数值积分,节点力、质量矩阵和其它单元矩阵的积分不是由解析计算的，而是应用数值解答，称为数值积分。最广泛应用的是Gauss积分,式中nQ个积分点的权重wQ和坐标值Q有表可查；见附录3。,指定方程在母单元域上进行积分，其积分区间为1，1。,一个二维单元的Gauss积分为,在非线性分析中，采用积分点数的规则一般基于在线性分析中的相同规则；对于一个规则的单元，积分点数目的选择是能恰好积分内部节点力。一个单元的规则形式，是指仅通过母单元的拉伸而不是剪切能得到的形式；例如，二维等参单元的一个矩形。 对于一个4节点四边形单元，如何选择内部节点力的积分点数目？ 由于速度是双线性的，单元中的变形率和

12、B矩阵是线性的。如果应力是与变形率线性相关，那么它将在单元内线性变化。内部节点力的被积函数是近似为二次的，因为它是B矩阵和应力的乘积。在Gauss积分中，对于一个二次函数的精确求解在每一方向上需要两个积分点，所以对于线性材料，需要22个点的积分得到内部节点力的精确解。 对于线性本构方程的积分，几乎得到内部节点力精确解的积分公式，称为完全积分。,4 编制程序,完全积分,4 编制程序,局部减缩积分,对于完全不可压缩或接近不可压缩的材料，运动必须是等体积的,0,J = 1,式中K是体积模量， 是剪切模量。在任意的等体积运动中，单元的整个体积将保持常数，即在整个单元中的运动必须是等体积的， 否则，当K

13、是一个非常大的数时(一个接近于不可压缩材料)，任何非零体积应变将吸收几乎全部的能量。,内部节点力的完全积分可能引起单元的自锁，即出现很小的位移而不收敛或收敛得非常慢。考虑一种线性材料，如果将线性弹性应变能分解为静水和偏量部分，可以写为,4 编制程序,局部减缩积分,为了克服这个困难，最容易的方法是使用局部减缩积分。在局部减缩积分中，压力为不完全积分，而应力矩阵的其余部分为完全积分。为此，将应力张量分解为偏斜部分和静水部分,体积自锁源于单元没有能力准确地表示一个等体积运动。为了消除自锁，必须设计应变场，这样在假设的应变场中整个单元的膨胀为零：为了避免自锁，对于任意保持单元体积的速度场，整个单元的应

14、变场必须是等体积的。对于等体积运动的四边形单元，为了克服沙漏模式，在整个单元中膨胀必须为零。,4 编制程序,局部减缩积分,在局部减缩积分中，压力为不完全积分，以此确定与外力平衡的应力场。如杂交单元,4 编制程序,局部减缩积分,注意到膨胀部分和偏斜部分是彼此正交的，因此内部虚功率为,将变形率表示为形状函数的形式，膨胀和偏斜部分的被积函数分别为,局部减缩积分包含在偏斜功率上的完全积分和在膨胀功率上的减缩积分。对于一个4节点四边形单元的局部减缩积分为,关于内力的局部减缩积分表达式为,4 编制程序,局部减缩积分,只有四边形和六面体单元才能采用减缩积分；而所有楔形、四面体和三角形实体单元采用完全积分。,

15、关于减缩积分的单个积分点是母单元的质心。偏斜部分是通过完全积分方式积分的。这种方法类似于应用于不可压缩材料线性分析的算法。关于其它单元的局部减缩积分算法，可以采取对于线性有限元类似的修正方法建立起来。,单元力和矩阵转换,例题：三角形3节点单元,4 编制程序,例题,四边形单元和其它二维等参单元,建立二维等参单元的变形梯度、变形率、节点力和质量矩阵的表达式。对于4节点四边形单元给出详细的表达式。以矩阵形式给出内部节点力的表达式。,是母单元的节点坐标,形状函数和节点变量,节点速度,节点坐标,变形率和内部节点力,4 编制程序,对于形状函数公式，映射到空间坐标不是可逆的。不能直接将,表达成 x 和 y

16、的形式，所以形状函数的导数是通过隐式微分计算出来,当前构形相对于单元坐标的Jacobian矩阵给出为,对于4节点四边形单元，上式成为,4 编制程序,Jacobian矩阵,考虑平衡方程,定义残差,通过Newton迭代求解和线性化,内力Jacobian切线刚度,外力Jacobian载荷刚度,整体Jacobian为二者之差,ABAQUS的UMAT就是求内力的Jacobian矩阵切线刚度,Jacobian矩阵是时间的函数，,的逆给出为,4 编制程序,4节点四边形单元采用单元坐标表示的形状函数的梯度为,采用空间坐标表示的形状函数的梯度为,出现在分母中,4 编制程序,速度梯度给出为,内部节点力,上式适用于

17、任意的二维等参单元。因为,是关于单元坐标的有理函数，所以不能对上式进行解析积分，通常使用数值积分。 对于4节点四边形单元，22 Gauss积分为完全积分。但是在平面应变问题中，对于不可压缩或几乎不可压缩的材料，完全积分会出现单元自锁，因此必须应用局部减缩积分。 对于4节点四边形单元，沿着每条边的位移是线性的。所以，其外部节点力与3节点三角形单元外部节点力是一致的。,出现在分母中，被积函数,4 编制程序,一致质量矩阵,根据Lobatto积分，使积分点与节点重合，或者将单元的整个质量平均分配到4个节点上可以得到一个集中的对角质量矩阵，分配到4个节点上给出,集中质量矩阵,例题,主从连接线,4 编制程

18、序,连接线经常用于连接采用不同单元尺度的网格部分，因为它们比使用三角形或四面体单元连接不同尺度的单元更方便。通过约束从属节点的运动，使其与连接主控节点的附近边界的场一致，保证跨过连接线的运动的连续性。 根据转换规则，导出节点力和质量矩阵。,例题,主从连接线,4 编制程序,通过运动约束给出从属节点的速度，使沿连接线两侧的速度保持协调，即C0。这个约束可以表达为节点速度的一个线性关系，可以写成为,连接后模型的节点力为,矩阵A是从线性约束得到的，“”表示在两侧连接到一起之前的分离模型的速度,主控节点力是分离模型的主控节点力和转换的从属节点力的和。这些公式对于外部和内部节点力都适用,一致质量矩阵为,4

19、 编制程序,4节点四边形单元，沿任何边界上的速度是线性的。从属节点3和5与 主控节点1和2重合，而从属节点4与节点1的距离为,节点速度可以写成为,节点力则给出为,主控节点1的力是,节点力两个分量的转换是一致的；这种转换对于内部和外部节点力都适用。如果两条边只是在法向相连接，则需要在节点上建立局部坐标系。通过公式，联系节点力的法向分量，而切向分量保持各自独立。,例题,主从连接线,5 旋转公式,在结构单元中，如杆、梁和壳，处理固定坐标系是很棘手的。例如，考虑一个旋转杆件。最初，唯一的非零应力是x，而y等于零。当杆件旋转时，以应力张量整体分量的形式表示单轴应力的状态就麻烦了。克服这种困难的一般途径是

20、在杆中嵌入一个随杆旋转的坐标系旋转坐标系,考虑一个坐标系,始终连接着节点1和2，在单轴应力状态下,杆的变形率类似地应用分量,描述,在旋转格式中考虑的关键问题是材料旋转的定义，有两种方法： 1 在每一个积分点嵌入一个坐标系，在某种意义上它随着材料旋转。对于任意大的应变和旋转都有效。在大多数情况下，应用极分解原理定义旋转。然而，当材料的指定方向有一个较大的刚度并需要准确地表现出来时，由极分解提供的旋转在Cartesian坐标系下未必能提供最好的旋转；见第9章。,5 旋转公式,对于一些单元，如杆或常应变三角形单元，整个单元的刚体转动是相同的，那么在单元中嵌入一个单一坐标系就足够了。对于高阶单元，如果

21、应变比较小，不需要坐标系准确地随材料旋转。例如定义旋转坐标系，使它与单元的一边重合。如果与嵌入坐标系相关的旋转为阶，则在应变中的误差具有2阶。这样，只要2与应变相比很小，单个嵌入坐标系是合适的。称作小应变、大转动问题。,2 在一个单元中嵌入一个坐标系并随着单元旋转。,5 旋转公式,在旋转坐标系下，矢量v的分量与整体分量之间的关系为,对于速度场，有限元近似的旋转分量可以写成为,速度梯度张量的旋转分量给出为,旋转变形率张量给出为,旋转公式仅用来计算内部节点力。外部节点力和质量矩阵常常在整体坐标系下计算。运动的半离散方程以整体分量的形式进行处理。,5 旋转公式,以旋转分量的形式建立,的表达式。内部节

22、点力的标准表达式为,根据链规则,从Cauchy应力变成为旋转应力的转换,应用R的正交性,表达式可以转化为矩阵形式,5 旋转公式,旋转Cauchy应力率是客观的(框架不变性)，因此，以旋转Cauchy应力率和旋转变形率之间的关系直接表示本构方程,对于次弹性材料,这里的弹性响应矩阵也表达为旋转分量的形式。上式中一个 令人感兴趣的特征就是对于各向异性材料，不需要改变,矩阵而能够反映转动。,另一方面，如果C矩阵的分量是在固定坐标系中表达的，对于各向异性材料，当材料转动时C矩阵随之变化。,5 旋转公式,应用旋转方法建立3节点三角形单元的速度应变和内部节点力的表达式。 单元的初始和当前构形如图。在初始时旋转坐标系与整体坐标系之间有 一个角度0；通常选择0为零，但是对于各向异性材料，将初始,轴定位在各向异性的一个方向可能是理想的；例如，在复合材料中，将,定位在纤维方向可能是有用的。旋转坐标系的当前角度是。,例题,三角形单元,5 旋转公式,例题,三角形单元,运动可以表示为三角形面积坐标的形式,在单元中的位移和速度场则