各向同性超弹性本构模型数值计算及验证

周建雄， 魏志刚，毛 欢， 韦雅宁， 刘迎松

(安徽工业大学 机械工程学院，马鞍山 243032)

1 引 言

橡胶材料由于其独特的物理和化学性质，作为减震和密封等部件广泛应用于工业工程中，建立真实可靠的本构关系来描述橡胶材料的力学特性变得尤为重要[1]。部分模型已经写入有限元材料库中用于仿真设计，如Arruda-Boyce[2]、Mooney-Rivlin[3]、Yeoh[4]和Ogden[5]等模型在ABAQUS中只需给出材料参数就可用于仿真模拟，而本构模型的选取也会影响仿真结果的准确性。近些年涌现较多拟合性和可靠性较好的本构模型，尚未开发相应模型本构子程序用于有限元做仿真模拟以更好地满足工程应用和科学研究[6]。

橡胶材料的弹性本构方程可表示为应变不变量的函数，如Mooney-Rivlin模型和Yeoh模型，也可表示为主伸长率 的函数形式，如Ogden模型和Valanis-Landel[7]模型。在ABAQUS中，不变量表示的超弹性模型，可用UHYPER编写材料子程序，Franus等[8]对MCIZ模型编写UHYPER子程序用于有限元实现对微可压缩超弹性薄壳的膨胀分析，Jemiolo等[9]编写Mooney-Rivlin和Blatz-Ko模型的UHYPER子程序，研究材料模型泊松比的取值对壳体结构的稳定性影响。由伸长率表示的超弹性本构模型需要用更具一般性的UMAT子程序编写。Asim等[10]开发了用于表征弹性材料线粘弹性和非线性超弹性的UMAT子程序，Suchocki等[11]开发了幂指数模型不可压缩超弹性UMAT子程序。

本文开发了几个新提出的橡胶弹模型的UHYPER和UMAT子程序。基于开发的UHYPER子程序进行了多孔橡胶板单轴拉伸有限元仿真，验证了子程序的正确性以及本构模型对复杂应变场预测的准确性；基于开发的UMAT子程序进行单轴、双轴和剪切变形实验模拟，验证了子程序的可靠性。

2 本构关系的数值计算

2.1 相关物理量的计算方法

在UMAT中编写主伸长率形式的弹性模型，需要定义应力张量和弹性张量，下面简要介绍这些量的得出方法。

2.1.1 变形张量

在连续介质力学变形理论中[12]，对于任何detF≠0的张量F，可极分解为F=RU=VR。式中R为正交张量，代表纯转动；右和左伸长张量U和V为对称正定张量，代表纯粹的变形。在有限应变状况下，常用右和左Cauchy-Green变形张量C和B来描述物体的变形。

C=FTF=(URT)(RU)=U2

B=FFT=(VR)(RTV)=V2

(1)

体积变化率J=V/V0=|F|，通过求解变形梯度F的行列式获得。

(2)

在数值计算中，所需特征值和特征向量可用迭代雅可比方法显式计算。Kopp[13]利用Jacobi算法、QL算法和Cuppen算法完成了三阶矩阵对角化的数值计算，并提供了C语言和Fortran语言开源代码，其中Jacobi算法最为精确，适合求解本文所需的特征值和特征向量，可将该代码写成子程序编入UMAT中。

2.1.2 应力张量和弹性张量

Connolly等[14]给出了主伸长率表示的超弹性本构模型数值实现方法，并基于Fortran77编写UMAT子程序用于ABAQUS，但其存在一定的局限性，其编写的UMAT子程序中未加入体积应变能部分。本文给出了应力张量和弹性张量详细推导过程，并将体积应变能函数代入公式推导，基于Fortran77编写UMAT子程序用于有限元模拟。

(1) 应力张量

材料的变形应变能可分解为等容变形和体积变形两部分，即

(3)

对应变能函数求导可得第二类Piola-Kirchhoff应力

(4)

式中 利用链式求导得

(a=1,2,3) (5)

橡胶材料变形一般为大变形，有限元仿真中使用Cauchy应力σ，其可由第二类Piola-Kirchhoff应力S获得

(6)

第一类Piola-Kirchhoff主应力Pa(名义应力)在工程中普遍应用，Cauchy应力与其变换关系为

(7)

可根据本构模型的函数形式选择应力变换，本文用于编写UMAT的Wei模型是用名义应力P表示的柯西弹性模型，则需将P应力变换为σ应力。

(2) 弹性张量

由于超弹性具有率无关性，物质弹性张量可以通过第二类Piola-Kirchhoff应力和变形张量等价定义为

Q=2∂S/∂C

(8)

物质弹性张量分解为等容部分和体积变化部分,Q=Qiso+Qvol。

等容弹性张量的谱形式为[15]

Na⊗Nb⊗Nb⊗Na)

(9)

式(9)第二项，当λa和λb数值相近时，分子分母趋于零，由洛必达法则得

(10)

(11)

体积变化部分为

(12)

上述本构模型中，材料是假定为不可压缩的，体积应变能为零。但完全不可压缩会在有限元仿真变形中导致数值奇异，因此有限元仿真实际上认为材料是近似不可压缩的。Abaqus中采用的体积应变能表达式为

(13)

式中Da取值决定了材料的可压缩性，Da取零则认为材料是完全不可压缩的。初始体积模量为k0=1/D1，体积模量过大会造成模型体积自锁，根据有限元模拟结果对Da取值进行不断调整，Da≈10-3时,有限元结果与模型理论值一致，且求解较为迅速。

在ABAQUS中，空间弹性张量用Cauchy应力的Jaumann率来定义[16]，表示为QJ。

QJ=J-1FFQFTFT=

J-1(Qiso+Qvol)+[(σ⊙I)+(I⊙σ)]=

[(na⊗nb)⊗(na⊗nb+nb⊗na)]+

[(na⊗nb)⊗(na⊗nb+nb⊗na)]+

[(σ⊙I)+(I⊙σ)]

(14)

最终的表达式中要求的量很多，且都与第二类Piola-Kirchhoff应力S有关，在应力张量部分有提到，Cauchy应力σ、第二类Piola-Kirchhoff应力S和第一类Piola-Kirchhoff应力P(名义应力)存在变换关系，则不管本构模型由何种应力描述，都可变换为第二Piola-Kirchhoff应力S进行推导。本文基于Fortran77完成对Wei模型和Ogden模型UMAT子程序开发，并将其用于有限元模拟来验证子程序的正确性。

2.2 使用模型

2.2.1 伸长率形式模型

(1) Wei模型[17]

在不考虑材料不可压缩特性时，Wei模型名义主应力为

(15)

式中i,j,k=1,2,3,且取值互不相等，a,c,d,e和m为材料常数，拟合橡胶三类简单试验数据[18]得,a=0.001，c=2.04，d=0.42，e=0.17,m=0.85。

对于不可压缩材料，该模型的柯西主应力偏量表达式为

(16)

将柯西主应力转变为编程所需名义应力偏量

(17)

名义主应力偏量关于修正伸长率的偏导为

(18)

(2) Ogden模型

(19)

式中μi和αi为材料参数，本文选择拟合性较好的三阶Ogden模型，拟合橡胶三类简单试验数据得[19],μ1=0.006，α1=4.17，μ2=0.709，α2=1.108，μ3=-0.006，α3=-2.125。

名义主应力偏量可由等容应变能对修正伸长率求偏导得

(20)

对于所需的名义主应力偏量及其关于修正伸长率偏导可参考式(17,18)求解。

2.2.2 不变量形式模型

使用UHYPER开发子程序遵循固定格式，需定义应变能函数，并计算应变能关于不变量和体积变化率的一阶、二阶和三阶偏导数，将其存储到相应数组中，常用变量有

U(1)=W

(21)

选取三个近来提出的超弹性模型编写UHYPER子程序来探究其预测复杂应变场的能力。模型如下。

(1) Exp-Ln[20]模型

(22)

式中A，a和b为材料常数，拟合未填充多孔硅橡胶板拉伸试验数据[21]得A=0.166，a=0.202，b=0.386。

(2) Mansouri[22]模型

W=A1{exp[m1(I1-3)]-1}+

B1{exp[n1(I1-3)]-1}

(23)

式中A1,B1,m1和n1为材料常数，拟合未填充多孔硅橡胶板拉伸试验数据得A1=1.54，B1=5.57，m1=0.089，n1=0.004。

(3) New neo -Hookean[23]模型

(24)

3 数值模拟

3.1 非均匀拉伸变形有限元模拟

Meunier对未填充多孔硅橡胶板进行拉伸试验，描述橡胶材料在复杂应变状态下的力学性质，且拉伸试验数据可用于评估超弹性本构模型预测复杂应变场的准确性。本文使用UHYPER子程序编写上述本构模型的子程序进行多孔橡胶板的拉伸仿真，评估这些模型预测结果的精度。

图1为采用的多孔橡胶板的形状尺寸，孔中心C3和C5之间沿圆心线做了切割。仿真中将底端固定，限制横纵方向上的自由度，在上顶面沿竖直方向施加57.3 mm的位移载荷，限制水平自由度。

图1 多孔橡胶板几何尺寸

有限元仿真得出的变形图轮廓(红色轮廓)与Meunier试验变形图(灰色图片)重合比较如图2所示。

图2 试验与仿真变形对比

在单轴拉伸变形状态下，伸长率为1.7时，图2中红色轮廓线与灰色图片重合度非常好，表明三个模型仿真结果和试验结果十分契合，所编子程序能够按照本构方程模拟计算，且所建橡胶本构模型可以准确预测橡胶在复杂应变场下的力学性质。由于此拉伸试验伸长率仅为1.7，是小变形试验，很难分辨出三个模型预测复杂应变场的优劣，可进一步增加伸长率，图3为伸长率为2.5时各模型预测多孔橡胶板的应力分布及变形仿真。

图3 应力分布及变形轮廓对比

可以看出，Exp-Ln模型和New neo-Hookean模型预测多孔橡胶板在伸长率2.5时的应力分布及大小较为相近，且变形轮廓重合度较好，相比之下，Mansouri模型预测的应力要小近30%，且孔洞附近的变形较前两个模型要稍大些。

3.2 单轴、等双轴和纯剪变形有限元模拟

在ABAQUS中使用由Wei模型及Ogden模型编写的UMAT子程序对三类简单拉伸实验做有限元模拟，获得柯西应力σ关于名义应变ε的仿真结果，由P=σ/(1+ε)将其转换为名义应力关于伸长率的关系，并将仿真结果与本构模型求解所得主变形方向下名义应力解析解进行对比。

3.2.1 Wei模型

(25)

图4 单轴拉伸状态下理论值与应用UMAT子程序获得有限元仿真结果对比

对比上述单轴、等双轴和纯剪变形有限元仿真结果和理论结果，直观上发现两者的曲线较为吻合，进一步绘制三类简单变形下仿真结果与理论结果间的相对误差，如图7所示。

图5 等双轴拉伸状态下理论值与应用UMAT子程序获得有限元仿真结果对比

图6 剪切变形状态下理论值与应用UMAT子程序获得的有限元仿真结果对比

图7 有限元模拟与Wei本构方程理论结果相对误差

由图7可知，对Wei模型所编UMAT用于橡胶单轴与双轴变形模拟的相对误差随伸长率逐渐增大，但最大误差不超过3%，而UMAT对剪切变形计算的相对误差随伸长率逐渐减小，最大相对误差出现在初始变形阶段，数值接近7%。表明所编UMAT与本构方程解析解存在偏差，该偏差可通过在ABAQUS中减小 取值来减小相对误差。

3.2.2 Ogden模型

与Wei模型的处理方法相同，结合式(20,25)得Ogden模型名义主应力理论值

(26)

图8 单轴拉伸状态下理论值与应用UMAT子程序获得的有限元仿真结果对比

图9 等双轴拉伸状态下理论值与应用UMAT子程序获得的有限元仿真结果对比

图10 剪切变形状态下理论值与应用UMAT子程序获得的有限元仿真结果对比

图11 有限元模拟与Ogden本构方程理论结果相对误差

对Ogden模型所编UMAT用于橡胶三类简单变形模拟的相对误差随伸长率逐渐增大，其相对误差较小，最大误差不超过2.5%。结合上文由Wei模型编写UMAT对三类简单变形实验模拟结果可知，所编子程序能够按照本构模型给出的力学表达式正确模拟计算，其相对误差来源于体积应变能函数Da的取值无法趋近于0，在模拟计算中无法避免。

4 结 论

针对各向同性超弹性本构模型数值计算，本文主要推导了主伸长率表示的超弹性本构模型数值实现方法，并基于Fortran77编写UMAT用于有限元模拟，所得结果如下。

(1) 结合Exp-Ln模型、Mansouri模型和New neo-Hookean模型开发UHYPER子程序用于多孔橡胶板拉伸仿真，模拟结果与实验数据十分契合，其中Exp-Ln模型和New neo-Hookean模型能更准确地预测较高伸长率下的复杂应变场。

(2) 由Wei模型和Ogden模型开发UMAT子程序用于单轴、双轴和剪切变形有限元仿真，仿真结果与本构方程解析解十分接近，表明所编子程序能够按照本构模型给出的力学表达式正确模拟计算，其中的偏差来源于ABAQUS中为避免出现体积自锁，体积应变能函数中Di的取值不能完全趋于0。

本文所提方法可用于主伸长率表示的超弹性本构模型数值计算，并编写UMAT子程序，可供后续用户开发使用。改进之处在于可编写相关求导子程序来替代用户求解主应力对伸长率的偏导。