1 概述

本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过

1）   基础理论

2）   商软操作

3）   自编程序

三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。

有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。

一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元求解器，通过自研求解器和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。

==几何非线性的T.L.和U.L.描述方法==

物理空间的很多量和方程都是相对某个参照系或者时刻点的，几何非线性的方程也不例外。K.J.Bathe教授1979年首次提出了应用于工程问题的几何非线性理论，他将几何非线性理论公式分为完全拉格朗格式(T.L. Total Lagrangrian格式)和更新拉格朗日格式(U.L. Updated Lagrangian格式)两种公式描述方式，物理方程中的所有的物理量统一在一种描述方式下表示，最终建立几何非线性的求解方程。此后，这两种描述方式称为商用结构有限元最广泛采用的几何非线性描述方法方法，此文将简单介绍一下物理量、网格的Lagrangina和Euler描述，虚功原理的T.L.和U.L.描述方法，然后通过和自主结构求解器iSolver的比对，验证Abaqus的几何非线性的描述方法。

1.1 初始构型和当前构型

我们先从构型的含义说起，几何非线性理论由三维连续体的虚功方程，提出了适用大位移、大转动几何非线性的解法。若把物体视为由无数质点所构成，并把这些质点称为材料点，一个物体从初始状态由于受到外部载荷运动，运动到另一个状态的过程中，每一个时刻，物体中所有材料点的位置的集合，定义为该物体的一个构型（configuration），且材料点不会凭空消失。

物体的材料点在0时刻占有的区域称为初始构型C0（initial configuration），材料点的位置在在当前t时刻占有的区域称为当前构型Ct（current configuration）。譬如下面的Abaqus和iSolver计算的几何大变形下的悬臂梁，初始时刻，悬臂梁在一个水平平面上，由未变形的体单元组成初始构型。而下压的任意时刻t，悬臂梁将弯曲，由变形的体单元组成当前构型。

1.2 物理量的Lagrangian和Euler描述方式

为描述各个物理量，可以为初始构型和当前构型选取不同的参考坐标系。当然，一般情况下，无论是初始构型还是当前构型，都以空间中固定不变的全局直角坐标系O来参考，此时，C0下任意一个材料点a在该全局坐标系O下记为X，因为材料点不会凭空消失，在t时刻C1构型下，必然会移动到某一位置A，该点在同一个全局坐标系下的坐标为x。

在数学上，X和x就相当于两种自变量，对于参考坐标系O中任意的变量y，就有两种表示方法。

1.2.1 Lagrangian描述

第一种：表示为X的函数：y=f(X)。此时我们称为Lagrangian（Lagrangian英文后面统一简称为L）坐标描述方法。典型物理量：Green应变。

• 注：Green应变也称为Green-Lagrange 应变，猜测可能就是因为是Lagrange描述的原因。

1.2.2 Eulerian描述

第二种：表示为x的函数：y=f(x)。此时我们称为Eulerian（Eulerian英文名后面统一简称为E）坐标描述方法。典型物理量：真实应变。

1.2.3 两种描述分析

这两者表述的是同一个模型，举个简单的例子，譬如材料点随时间t的运动为x=2*X*t。那么在坐标系O下，在t时刻的速度的Lagrangian和Eulerian描述方法分别为：

v(X)=2*X和v(x)=2*x/2t=x/t

这两种描述方法都是相对同一坐标系O的物理量的表示方法，那么得到的数值结果完全一致，只不过写成函数的时候自变量不一致一样。

1.3 网格的Lagrangian和Eulerian描述

在讲T.L.和U.L.之前，我们先解释一下L.代表的Lagrangian网格的含义。前面所述物理量的坐标描述也有L和E之分，网格也有L和E之分。这两者是不同的概念。在L和E网格定义中，表示的是两种网格随时间变化的位置定位方式，而物理量的坐标描述仅仅是变量的表示方式，这两者没有直接关系。为了说明L和E网格的区别，举一个简单的例子，譬如二维平面的运动如下：x=X+tY, y=Y，那么一个长方形运动后的变形为一个平行四边形，具体的材料点a(在全局坐标系O中坐标为X，Y)变形后为A（在O中坐标为x，y）。

网格主要由节点组成，而节点在全局空间O中的坐标决定了网格的形状。节点的坐标有两种取法：

1.3.1 网格的Lagrangian描述

无论在哪个时刻，取同一个材料点变形后的xi，yi坐标，此时如果一开始网格点在a点，那么变形后必然在A点，而且，如果用L描述，那么可以得到，Xi，Yi=常数。那么得到的网格将类似下方，随构型变形而自动更新变化。这就是Lagrangian网格。

1.3.2 网格的Eulerian描述

无论在哪个时刻，取同一个材料点变形前的Xi，Yi坐标，此时如果一开始网格点在a点，那么变形后还是在a点。那么得到的网格将类似下方，固定在空间不变。这就是Eulerian网格。

结构有限元中一般都是Lagrangian网格，这也是为何T.L.和U.L.方法称为Lagrangian方法的原因。

1.4 虚功原理的T.L.和U.L.描述方式

有限元中，T.L.和U.L.用于非线性问题中的虚功原理的描述方式，是整个有限元求解方程的描述方式，而前面物理量的描述方式是虚功原理中各个量的描述方式。

非线性问题在有限元中都是以增量法来求解，增量法的基本思想就是将施加外载荷的过程分割为若干个时间增量步，假定这个时间增量步是固定的，dt,2dt,3dt,….，每一步只施加一个比较小的载荷增量，总载荷引起的位移就是所有增量步位移增量的累加。

具体的理论和Abaqus实现过程可参考我们以前系列文章4：非线性问题的求解。

由虚功原理在(n+1)dt时刻应变能为：

上面每个量都是都是参考(n+1)dt时刻的坐标的，此时将相对参考时刻的坐标写在左下方，将要求的时刻点写在左上方，那么就是：

1.4.1 T.L.描述

注意，有些读者认为T.L.只能用于小应变，但这种说法个人认为并不正确，这里的推导没有涉及是否是大应变还是小应变，同时，后面文章我们将证明U.L.和T.L描述是可以等价转换的，既然U.L.可用于大应变，那么T.L.也适用于大应变、大转动、大位移情况。

1.4.2 U.L.描述

以当前时刻ndt构型为参考构型，用Eulerian坐标x描述物理量，那么就是U.L.方法（Updated Lagrangian方法）

• 上式的转换数学上并不严格，还差一个应变对时间的导数，严格推导和系列文章22：几何非线性的刚度矩阵求解类似，需要将上述表达式先做线性化。

1.5 Abaqus的内部描述方式验证方法

在Abaqus所有的帮助手册中，找不到任何T.L.还是U.L.的信息，由上节所述，我们只能反过来，通过查看Abaqus的变量的度量方式，来猜测Abaqus用了T.L.还是U.L.方式，Abaqus文档中多处提到了各种变量的度量，但并不一定真实，而且前后明显矛盾，最好的方式还是直接查Abaqus程序的结果。

注意，此处不一定能采用后处理模块中显示的应变量。为了和试验对应，Abaqus几何非线性后处理显示是自然应变和Cauchy应力，它将程序中的相关度量在显式时做了转换，与程序中的真正用于计算的度量并不一致。详细说明可以查看本系列文章27： Abaqus内部计算和显示的应变。

我们通过UMAT等子程序来查看Abaqus子程序接口中的应变应力度量，同时和iSolver采用同样度量的结果比对，大体猜测Abaqus几何非线性的描述方式。

1.5.1 Abaqus中T.L.的描述

Abaqus中S4R5、STRI3、STR65、S4RS、S8R5、S3RS、B33，B23等单元都采用T.L.的描述方式。采用一个简单的算例，证明如下：

1.5.1.1 算例介绍

参数如下：

尺寸：5X1，厚度0.1。

材料：Young’s Modulus 1e8， Poisson Ratio 0.3。

左侧两个节点固支。

右侧两个节点每个加集中力1e5，x方向。

划分为一个单元。

1.5.1.2 Abaqus应变结果

Abaqus中采用S4R5单元，几何非线性单元NLGeom打开，得到一个迭代步后结果：

1.5.1.3 iSolver应变结果

iSolver中也采用S4R5的度量方式，同时采用T.L.方式，得到的应变度量如下，可发现，和Abaqus完全一致。

1.5.2 Abaqus中U.L.的描述

Abaqus中C3D8/C3D8R、S4/S4R、S3/S3R、B31，B21等非线性单元都采用T.L.的描述方式。采用一个上述同样的算例，证明如下：

1.5.2.1 Abaqus应变结果

几何非线性开关NLGeom=On，同时单元类型改为S4R。得到应变：

1.5.2.2 iSolver应变结果

iSolver中也采用S4R的度量方式，同时采用U.L.方式，得到的应变度量如下，可发现，和Abaqus完全一致。

1.6 T.L.和U.L.描述方式的本质

既然Bathe教授在开辟几何非线性时就创立了T.L.和U.L.两种虚功原理的描述方式，而且有限元商软也是按这两种方式实现的，那Abaqus中为何从头到尾没提T.L.还是U.L.的描述方式呢？

我们猜测这是Abaqus有意为之。这两种虚功原理的描述本来就是描述的同一个物理对象，所以理论上你采取任何一种都可以得到最终的正确解，但无论哪种描述都需要它用到物理量的描述方式结合起来满足一定的要求，譬如功的共轭、能量守恒、动量守恒等条件。同时，因为有限元中物理量的描述方式可以相互转换的，为了有限元的结果更能和实际试验的度量对应，Abaqus也会在T.L.中采用Euler描述的物理量的度量方式或者U.L.中采用Lagrangian描述的物理量的度量方式，譬如它本身提供的UMAT是真实应变和真实应力，但你完全可以不用管UMAT输入的DSTRAN等应变值，而用它提供的另一个参数变形梯度DFGRD1按Lagrangian描述方式实现Green应变（当然，这样得到的Jacobian矩阵DDSDDE和应力STRESS都没问题，但你改不了Abaqus内部利用STRESS来计算初始应力刚度和本身的材料刚度矩阵的算法，所以还需要Abaqus后面提供更多的接口才行），或者在UEL采用Green应变和2nd PK力来实现刚度阵AMATRX和内力RHS，也就是说已经突破了传统的T.L.和U.L.的界限，所以干脆Abaqus就不再提T.L.还是U.L.方式了。

==总结==

本文共介绍了三种近似相关的描述方式，Abaqus和iSolver的对应实现方式分别如下：

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删