学力学出身，我们也未必能全面认识非线性力学

weixin

导言：搞力学软件开发这么些年，个人经验感觉即使是学力学出身的人也未必能全面地了解非线性力学。比如：

  · 看到过有人写论文把全量Lagrangian（以下简写为TL）形式的公式强行变换为更新Lagrangian（以下简写为UL）公式。

  · 若干年前我开发的开源软件FrontISTR公开后，有人通过该项目开发的负责人找到我，指责其基于UL的弹塑性解析部分的程序是错误的。因为在他的测试计算中，即使仍处在弹性状态的加载后完全卸载，变形量不能回零。

上面两例涉及的人员无疑是专业人士，但是第一例中所做工作有点无事生非的意思，第二例则反映了对UL法性质的不理解。本文试图对对非线性力学做一个简单概括。由于不想把问题摊得太开，本文将仅限于弹性材料。另外，线弹性太简单，本文就不赘述了。

01、什么是弹性

变形体在经历加载变形，完全卸载后，完全恢复到未变形前的原始形状。本文将具有这一变形特征的变形称为弹性变形。如果你有一定的数学功底，要定义出一套满足这一条件的数学公式想必不难。但是，现实往往充满着种种变数......

02、非线性弹性力学模型

（1）Cauchy弹性

Cauchy弹性本构表述Cauchy应力σ 为左右伸长张量U、V 或其他不受刚体旋转影响的变形度量，如C=FTF 的函数

参考Cayley-Hamilton理论，当σ 和C 是各向同性张量时，可以写为下述代数形式

很明显，采用该模型时，当变形为零时，应力也为零。这一弹性的基本特征可以保证。但是这一公式不能保证封闭变形路径下的能量守恒。Cauchy弹性模型好在简单，坏在粗暴。个人认为，这是一个应该湮灭于历史尘埃的的模型。也就不再多说了。当然更主要的是，下面的模型堪称完美。

（2）Green弹性，超弹性 (Hyperelasticity)

如果弹性变形能函数是一个状态函数，或者说是可以写成恰当微分或全微分形式的函数，那么它必将是能量守恒的。而且，可以保证在零变形下的状态不变。我们可以利用这一性质来构造弹性模型。

记变形体的弹性势能为某变形度量A 的函数为Ψ(A)。那么与其对偶的应力就可以写成：

那么，我们应该选择哪个变形度量？很多现代的连续体力学教材都用不小的篇幅来论述有限变形度量方法。其实如果我们退出力学领域，从数学的角度远眺一下这个问题，会觉得它是很简单的。对于表述尺寸和角度变化这样的变形问题，微分几何学告诉我们可以使用度规张量。记未变形状态下的黎曼度规为G，变形后的黎曼度规为g，那么

就可以完美地表现该变形过程的特征。而这个定义就是连续体力学教材中出现的Green-Lagrangian应变。

从上式出发，我们可以得到E 的共轭应力。关于共轭的概念，在我的这篇文章中有一些简单的介绍，可以参考。从共轭的角度出发，你可以自由地调换上述能量泛函中的应力，应变定义。在实际应用中，Cauchy-Green变形张量和第2PK应力似乎最为常见。

超弹性模型是一个完美地表述了弹性变形特征的模型。由于使用了状态函数，其最终变形状态与变形路径无关。所以说，如果使用两个力学计算软件来计算同一超弹性变形问题，其计算结果必须是基本相同了。如果不同，那么，我们可以合理地怀疑其中至少一个软件实装是错误的。

由于超弹性材料本构要求计算与变形前度规的比较，这就是所谓全量Lagrangian算法的由来。这种算法是以变形前构形为基准来进行计算的。如果你非要用更新Lagrangian法来计算超弹性（正如在导言中举得例子那样），你需要把上述本构关系通过一系列坐标变换转换到当前构形。要干这样的活儿，为了推公式多死几个脑细胞到没什么，怕的是计算机也觉得心塞：凭什么增加我的计算量，要加工资云云。

到这里，弹性计算问题似乎是解决了。但是，你要知道，这个世界往往是不完美的！

（3）亚弹性 (hypoelasticity) 模型

还存在一种微分或率形式的弹性本构dσ=f(dε)，称为亚弹性模型。导入这类模型的原因是这个世界上除了有与变形路径无关而变形的弹性材料，还有与变形履历相关的材料。塑性变形完全取决于变形体当前的变形状态。比如，如果你使用associated flow rule，那么屈服面的法线方向就是当前的变形方向。这时，采用上面的率形式的本构方程就变得是理所当然的了。对于弹塑性材料需要做弹塑性分解dε=dεe+dεp，然后把弹塑性本构写成dσ=f (dε-dεp) 的形式。在这种情况下，如果变形还处在弹性变形状态，我们使用的就是亚弹性模型。

与上一节相同，我们仍然从几何学的角度来选择应力，应变的定义。首先考虑当前构形下黎曼度规的变化，那么用黎曼度规g 沿位移u 方向的李微分就非常合适。度规的李微分是一种比较特殊的李微分，把它展开可以得到([3]p.14)

这个就是Cauchy小应变张量ε 。由于与其共轭的应力张量是Cauchy应力张量σ ，因此我们可以亚弹性公式写为

（注：虽然数学上对二阶对称应力张量的李微分的定义是严格的，但是经过连续体力学研究者多年的努力，在连续体力学界存在多种客观性应力率。这些客观性应力率实际上只有一个是符合严格的李微分的定义。其他的都是丐版（[4], Chapter 1, Box 6.）

看起来一切顺利？但是微分形式的弹性本构关系没有任何机制保证弹性变形下的能量守恒和变形的复原特性。要证明这一点很简单：找到一个反例即可。比如下面的文献[5]~[7]。

结合上面的论述，回到引言中提到的故事2：指责FrontISTR的弹塑性解析不能再现弹性变形的特征是正确的。但是，指责软件的实装有误就没道理了，因为这本来就是亚弹性模型特征之一。

与超弹性材料本构要求全量Lagrangian算法相对应，微分形式的本构用更新Lagrangian算法来处理最为自然。这就是导入更新Lagrangian算法的基本原因。

微分形式的材料本构也给数值计算提出难题：积分运算怎么进行？我们知道，单靠增量的累加的计算精度是可疑的，需要优秀的算法去弥补。这方面的问题说起来话太长，可以参见文献[8]~[10]。

也有研究者，如[11]，关心微分形式弹性本构的不完备问题，试图找到一个可积的本构关系。不过这种打补丁式的研究方向缺乏一种简洁优雅的气质，非理论研究的正道。更重要的是，下面的方法要清晰明了得多。

03、乘式弹塑性分解

把变形梯度F 按弹塑性成分做如下分解F=F e+F p。它对应与下面的物理图景。

变形梯度的弹塑性分解

采用这样的方法原则上可以对弹塑性变形分别处理，得到在弹性变形范围内能量守恒的计算结果。遗憾的是，乘式分解并不适应于数值计算。从这个起点出发还需要大量的公式变换，并不轻松。

4. 在开源软件FrontISTR中的实装

FrontISTR具有处理各种材料本构关系的功能。首先，材料定义有一个NLGEOM开关，如果关闭，所有的计算都采用小变形计算。如果打开，那么

  · 弹性材料，粘弹性材料 ——》采用Cauchy-Green变形张量和第2PK应力描述变形，采用TL计算公式

  · 弹塑性材料，粘塑性材料 ——》采用个Green应变增量和Cauchy应力描述变形，采用UL计算公式

（本文上文并未论及粘性，粘弹性可以理解为如果时间足够长，卸载后的变形体可以恢复原状的变形体）

用户在使用过程中并不需要知道其中的差异，原则上也不允许对算法进行修改。但是这个留有一个后门（在用户手册上没有记载），允许力学高手自行调整算法。比如在输入文件中写下

！HYPERELASTICITY， CAUCHY

这里的CAUCHY关键字将告诉计算程序用UL法亚弹性增量计算。而输入

！PLASTICITY，KIRCHHOFF

中的关键字KIRCHHOFF将告诉计算程序用TL法全量的Cauchy-Green变形张量和第2PK应力来进行计算。

这原本是我用来调试程序用的功能。虽然考虑到避免滥用而隐藏了起来，但是并非完全无用。比如说你的材料是弹塑性体，但是在你的问题中大部分变形都处于弹性变形状态，塑性变形领域很小（比如火车车轮，只有它与轨道接触的一小块区域才有大变形）。此时，采用TL法就不失为一个明智的选择。

04、一点感言

敲了这么多文字难免带一些发散的思考。感觉现在的亚弹性，弹塑性计算仍然问题多多。需要某个数学高手来一个总括。另外，放弃旧有思路，从李微分，李代数的方向下手未必不是一个可行的选择。本文也在尝试从度量空间的角度切入非线性连续体力学这个主题，衷心希望能给读者一些有益的提示。

参考文献：
[1] Y.B.Hu, R.W.Odgen: Nonlinear Elasticity, Theory and Applications, Cambridge University Press, 2001
[2] Stuart S. Antman: Nonlinear Problems of Elasticity, Springer, 2005
[3] Nick Sheridan :　Hamilton's Ricci Flow
[4] J. E. Marsden and T. J. R. Hughes. Mathematical Foundations of Elasticity. Dover Publications, Inc., New York, 1994
[5] Koji, M., Bathe, K. J., "Studies of finite element procedures-stress solution of a closed elastic strain path with stretching and shearing using the updated Lagrangian Jaumann formulation", Computers and Structures, Vol. 26, 1987, pp. 175–179.
[6] Truesdell, C., Noll, W., “The non-linear field theories of mechanics”, third edition. Springer: Berlin, 2003
[7] Lin, R. C., Schomburg, U., Kletschkowski, T., "Analytical stress solutions of a closed deformation path with stretching and shearing using the hypoelastic formulations", European Journal of Mechanics–A/Solids, Vol. 22, 2003, pp. 443–461.
[8] Jamshid Ghaboussi、 David A. Pecknold、 Xiping Steven Wu : Nonlinear Computational Solid Mechanics, CRC, 2017
[9] Allen C etal : A Comparison of Lagrangian-Eulerian approaches for tracking the kinematics of high deformation solid motion, SAND2009-5154
[10] Miklos Zoller: Error Analysis on Numerical Integration Algorithms in a Hypoelasticity Framework
[11]  Xiao, H., Bruhns, O. T., Meyers, A., “The integrability criterion in finite elastoplasticity and its constitutive implications”, Acta Mechanica, Vol. 188, 2007, pp. 227-244.

来源：仿真秀App微信公众号（ID：fangzhenxiu2018），作者：hillyuan 力学博士  仿真软件开发者。