模拟线弹性材料能有多难？

线弹性模型是结构力学分析中最基础的材料模型。虽然听上去微不足道，但模型中却包含不少难以一眼看出的重要细节。在本篇文章中，我们将深入讨论线弹性材料模型的相关理论和应用，并且大致介绍其各向同性和各向异性、材料数据的容许值、不可压缩性，以及与几何非线性之间的相互作用。

各向同性线弹性

在绝大多数涉及线弹性材料的仿真中，我们都要模拟完全不具有方向敏感性的各向同性材料。描述这种材料只需要两个独立材料参数。有很多方法来选择合适的参数，不过其中某些参数则更为常用。

杨氏模量、剪切模量和泊松比

杨氏模量、剪切模量和泊松比是材料数据表中最常见的参数。它们不是独立参数，因为剪切模量 $G$ 可以用杨氏模量 $E$ 和泊松比 $\nu$ 推算出来，计算公式为

$G=\frac{E}{2(1+\nu)}$

杨氏模量可以通过单轴拉伸试验直接测量，而剪切模量可以通过纯扭试验测量。

在单轴试验中，泊松比用于确定材料的横向收缩（或拉升）程度。其容许范围为 -1 < $\nu$ < 0.5，当其为正值时，表示材料受拉时会在厚度方向上发生收缩。有一些材料具有负泊松比，它们被称为 Auxetics。葡萄酒瓶中软木塞的泊松比接近于零，所以无论是拔出或推入，它的直径都不受影响。

对于大部分金属和合金而言， $\nu\approx$ 1/3，所以剪切模量大约等于杨氏模量的 40%。

在给定 $\nu$ 的可能取值后，剪切模量和杨氏模量的比值范围为

$\frac{1}{3}<\frac{G}{E}<\infty$

当 $\nu$ 接近 0.5 时，材料变得不可压缩，此类材料会给分析带来特定问题，我们将在下文说明。

体积模量

体积模量 $K$ 反映一定均匀压力下的体积变化。它可以用杨氏模量 $E$ 和泊松比 $\nu$ 推导出来，表示为：

$K=\frac{E}{3(1-2\nu)}$

当 $\nu$ = 1/3 时，体积模量与杨氏模量的值相等，但是对于不可压缩材料（ $\nu$ $\rightarrow$ 0.5）， $K$ 趋近无穷。

体积模量往往与剪切模量被一起指定。从某种意义上讲，这两个参数是最独立的参数选择。体积变化仅仅取决于体积模量，扭曲则完全受剪切模量决定。

拉梅常数

拉梅常数 $\mu$ 和 $\lambda$ 常见于论述弹性的数学论文。我们可以用拉梅常数方便地表述应力张量 $\sigma$ 和应变张量 $\varepsilon$ 间的全三维本构关系： $\sigma=2\mu\varepsilon+\lambda trace(\varepsilon)I$ 常数 $\mu$ 即为剪切模量， $\lambda$ 则可以写为

$\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)}$

点击此链接，可查看各类弹性参数之间的完整转换表。

线弹性材料的不可压缩性

橡胶一类的材料几乎不可压缩。从数学角度来说，完全不可压缩材料与可压缩材料具有本质的区别。因为不会发生体积变化，因此无法确定其平均应力。平均应力（压力）p 是体积变化 $\Delta$ $V$ 的函数，所以其状态方程 $p=f（\Delta V）$ 不成立，而必须用一个约束说明进行替代 $\Delta V=0$ 不可压缩性还有另一个角度需要注意， $（1-2\nu）$ 项会出现在本构方程的分母中，即如果 $\nu$ = 0.5 ，则方程将除以零。那么我们是否可以设 $\nu$ = 0.499 ，从而对不可压缩材料实现近似模拟吗？

想法可行，但在这种情况中，基于标准位移的有限元公式可能得出不理想的结果，这是由“锁定”现象引起的。造成的后果包括：

模型过于生硬。
应力呈棋盘式分布。
方程的病态导致求解器发生错误或发出警告。

补救方法是使用“混合方程”，将压力作为额外自由度引入。在 COMSOL Multiphysics 中，勾选材料模型设置窗口中的“几乎不可压缩材料”复选框，即可启用混合方程。

当泊松比约大于或等于 0.45 时，体积模量比剪切模量大超过一个数量级，因此使用混合公式是明智的做法。其效果示例如下图所示。

一个简单的平面应变模型中的应力分布，v= 0.499。上方的图表示基于标准位移的方程，下方的图表示混合方程。

在仅涉及位移自由度的解中，其应力分布图在左端（即存在约束的位置）呈现出扭曲的状态。使用混合公式后，这些扭曲几乎全部消失。

正交各向异性和各向异性

在一般情况下，线弹性材料的材料属性都具有方向敏感性。其中最普遍的特性是各向异性，这意味着全部六个应力分量都取决于各自不同的应变分量。完整表示这些分量需要 21 个材料参数，很明显，获取全部数据是一项艰巨的任务。如果将应力 $\sigma$ 和应变 $\varepsilon$ 当作矢量，那么二者的关系就可以用“6×6”本构对称矩阵 $D$ ，通过如下方程式进行表示 $\sigma=D\varepsilon$ 幸运的是，各向异性材料通常会表现一定的对称性。在正交各向异性材料中，有三个正交方向上的剪切作用和轴向动作实现解耦。也就是说，当材料沿着其中一个主方向拉伸时，它只会在两个正交方向上收缩，而不会受剪切力的作用。完整描述正交各向异性材料需要九个独立材料参数。

当以柔度形式记录时，正交各向异性材料的本构关系会更加清晰明了，其中 $\varepsilon=C\sigma$ ：

由于柔度矩阵必然是对称的，因此使用的十二个常数可通过符合下方形式的三个对称关系减少为九个

$\frac{\nu_{YX}}{E_{Y}}=\frac{\nu_{XY}}{E_{X}}$

请注意 $\nu yx\ne\nu xy$ ，因此在处理正交各向异性的数据时，一定要确保使用预期泊松比值。并不是所有来源中的数值都是相同的。

各向异性和正交各向异性常见于非均质材料。其材料属性通常不是由测量得到的，人们更倾向于通过从微观到宏观尺度的均质化过程计算这些属性。关于这种均质化作用在完全不同的研究背景中的讨论，请访问这篇文章进行查看。

对于非各向同性材料来说，使用描述各向同性材料的类似材料参数进行计算，在取值上可能会受到一些限制。我们虽难以立即发现这些限制，但有两个注意事项会对我们找到限制有所帮助：

本构矩阵 $D$ 必须是正定的。

a. 对于各向异性材料，唯一的选择是检查是否所有特征值都为正。

b. 对于正交各向异性材料，适用条件为：全部六个弹性模量皆为正，且

$\nu_{XY}\nu_{YX}+\nu_{YZ}\nu_{ZY}+\nu_{ZX}\nu_{XZ}+\nu_{YX}\nu_{ZY}\nu_{XZ}<1$

2. 如果材料的压缩率低，则必须使用混合方程。

a. 我们可以估计等效体积模量和剪切模量的值。

b. 在不确定的情况下，为了避免可能出现的误差，最好在混合方程中引入额外自由度。

几何非线性

在解决几何非线性问题时，“线弹性”的含义实际上是一个常规问题。这里的问题是，我们有多种应力和应变的数学表示方式。如希望了解有关应力和应变的不同测量方式，请阅读之前的一篇文章。

因为在 COMSOL Multiphysics 中，主应力和应变量分别为第二 Piola-Kirchhoff 应力和 Green-Lagrange 应变，因此线弹性自然地被解释为两个量之间的线性关系。人们有时将这种材料称作 St. Venant 材料。

人们通常凭直觉认为“线弹性”指的是简单拉伸试验中力与位移的线性关系。事实并非如此，因为应力和应变都取决于变形。为了理解这一点，我们来看一看横截面为正方形的条块。

条块的初始长度为 $L_{0}$ ，初始横截面积 $A_{0}=a_{0}^{2}$ ，其中 $a_{0}$ 表示横截面的初始边长。假设条块的延伸距离为 $\Delta$ ，故现在的长度 $L=L_{0}+\Delta=L_{0}(1+\xi)$ 在这里， $1+\xi$ 表示轴向拉伸， $\xi$ 可解释为工程应变。横截面的新边长为 $a=a_{0}+d=a_{0}(1+\eta)$ ，其中 $\eta$ 表示横向工程应变。力可表示为轴向 Cauchy 应力 $\sigma_{x}$ ，与当前横截面积的乘积：
$F=\sigma_{x}A=\sigma_{x}A_{0}(1+\eta)^{2}$
为了使用线弹性关系，Cauchy 应力 $\sigma$ 必须表示为第二 Piola-Kirchoff 应力 $S$ 。变换规则为 $\sigma=J^{-1}FSF^{T}$ 其中 $F$ 表示变形梯度张量，体积比例定义为 $J=det（F）$ 。不考虑细节时，单轴的情况可表示为

$\sigma_{x}=\frac{F_{xX}}{F_{yY}F_{zZ}}S_{X}=\frac{(1+\xi)}{(1+\eta)^{2}}S_{X}$

由于对受到单轴拉伸作用的 St. Venant 材料来说，其轴向应力与轴向应变的关系为 $S_{x}=E\epsilon_{x}$ ，由此得出

$F=S_{x}A_{0}(1+\xi)=EA_{0}(1+\xi)\varepsilon_{ X}$

已知 Green-Lagrange 应变张量在轴向上的项定义为

$\varepsilon_{X}=\frac{\partial_{u}}{\partial X}+\frac{1}{2}(\frac{\partial_{u}}{\partial X})^{2}=\xi+\frac{1}{2}\xi^{2}$

所以力与位移的关系为

$F=EA_{0}(1+\xi)(\xi+\frac{1}{2}\xi^{2})=EA_{0}(\xi+\frac{3}{2}\xi^{2}+\frac{1}{2}\xi^{3})$

具有几何非线性的线弹性材料实际上意味着力与工程应变（或者力与位移之间，因为 $\Delta=L_{0}\xi$ ）之间存在立方关系，如下图所示。

如图所示，受压侧的材料刚度接近于零： $\xi=\sqrt{1/3}-1\approx-0.42$ 。在实际操作中，这意味着在这一应变水平下，仿真就会失败。有人认为，现实中不存在大应变下呈线性的材料，所以在实践中不会出现这样的问题。然而线弹性材料通常在远远超出应力合理范围的情况下使用，例如：

通常，在引入更复杂的材料模型之前，您可能想快速检查一遍“数量级”。
模型中存在奇异点，并且致使某一点上产生了极高的应变。

点击此处，阅读更多关于奇异点的知识。

接触问题中的研究总是围绕几何非线性问题。

高压缩应变通常在分析过程中出现于某一时刻的局部接触区域。

对于以上所有情况来说，如果压缩应力过大，求解器也许会无法求解。如果您怀疑我们的示例属于这种情况，绘制最小主应变是一个很好的检测方法。如果它小于 -0.3 左右，我们就预测到会发生类似故障。由 Green-Lagrange 应变得到的临界值结果为 -1/3，如果这是一个问题，您应该考虑更换一个合适的超弹性材料模型。

压缩或许不是唯一的问题。在上述分析里，泊松比没有出现在方程中。那么横截面的情况如何呢？

根据单轴情况中的定义，横向应变与轴向应变的关系为 $\varepsilon_{Y}=-\nu\varepsilon_{X}$ 当这些应变是 Green-Lagrange 应变时，这便是一个非线性关系，可表述为

$\frac{\partial\upsilon}{\partial Y}+\frac{1}{2}(\frac{\partial\upsilon}{\partial Y})^{2}=-\nu(\frac{\partial\upsilon}{\partial X}+\frac{1}{2}(\frac{\partial u}{\partial X})^{2})$

因此横截面的变化具有很强的非线性。求解这个二次方程可得出如下的工程应变之间的关系

$\eta=\sqrt{1-\nu(\xi^{2}+2\xi)}-1$

结果如下图所示。

St. Venant 材料受到单轴拉伸作用时，其横向位移随轴向位移而变化。图中显示了五个不同的泊松比的值。

如您所见，当泊松比的值更高时，横截面在大拉伸应力的作用下的塌陷更为迅速。

如果选择另一种应力和应变数学表述方式，例如 Cauchy 应力与对数或“真实的”应变成正比，那么我们将得到完全不同的响应。当力-位移响应取决于泊松比的值时，这种材料的刚度反而会随着拉伸而下降。虽然在两个不同的仿真平台中，利用大应变弹性计算出的结果存在巨大差异，但两种材料仍旧皆可称为“线弹性”材料。

关于线弹性材料的结束语

本文阐释了线弹性材料在使用上的一些限制，重点强调了不可压缩性，以及线弹性与大应变的结合可能产生的使用误区。

如果您有兴趣了解更多涉及结构力学问题的材料建模内容，请阅读下列文章：

经授权转载自cn.comsol.com/blogs，原作者 Henrik Sönnerlind 。

来源：知乎 www.zhihu.com

作者：COMSOL 中国

【知乎日报】千万用户的选择，做朋友圈里的新鲜事分享大牛。
点击下载