基于自适应方法的橡胶刚度求解研究

廖小平 李达 韩雁青

【摘 要】针对工程分析中橡胶大变形无法计算的问题进行研究，采用自适应方法求解橡胶刚度曲线。建立与试验原型机一致的有限元模型，包含橡胶、橡胶底座及试验机冲头的系统。当橡胶发生大变形无法计算时，采用自适应方法继续求解。该方法通过提取变形后的网格并进行重划分，再进行结果映射，最后经过多次结果映射并重启动分析，得到完整的橡胶刚度曲线。应用实例表明：自适应求解方法能够解决橡胶大变形无法计算的问题，为橡胶参数反求提供条件，为橡胶新结构的设计阶段提供验证方法，缩短开发橡胶新结构的试验周期。

【关键词】橡胶；大变形；自适应方法；刚度曲线

【中图分类号】U463.33 【文献标识码】A 【文章编号】1674-0688（2016）08-0067-04

0 前言

在汽车工业中，橡胶材料的元件主要集中在发动机、变速器、传动轴与车架或车身的连接处，以及悬架系统中。橡胶的主要作用是衰减发动机、传动系及路面传递到车身的振动，可以在很大程度上改善车辆的NVH性能。

橡胶材料的变形是一个非常复杂的过程，其伴随着大位移、大应变。近年来，橡胶本构模型一直是国内外学者研究的热门课题。首先在橡胶有限元分析中，橡胶本构模型类型的选取，对分析精度尤为重要。刘萌[1]等人采用Marc非线性有限元分析程序，用工程实用的测试方法，对橡胶材料单轴拉伸力学行为进行测试，拟合处理得到Mooney Rivlin模型材料常數，并与实际测试结果进行对比。结果表明，有限元模型计算所得力与位移关系曲线与实测曲线吻合较好，相对误差小于6%。殷闻[2]等人则对比分析了2种橡胶本构模型Mooney Rivlin和Yeoh的参数，并进行了数学计算。他们采用非线性有限元软件对三维哑铃试样进行不同载荷水平下的单轴拉伸试验仿真，得出Mooney Rivlin模型能更好地模拟橡胶的中小变形行为的结论。

对于复杂结构的橡胶，用常规的分析方法无法得到完整的橡胶刚度曲线，以致橡胶制品的可靠性及疲劳性能一直停留在试验验证的阶段。如果采用自适应方法重划网格达到完整的橡胶刚度曲线，那么有限元技术的应用将缩短橡胶的设计周期，降低开发成本。因此，当已知橡胶衬套的设计刚度曲线时，根据工程经验由具体的结构形式确定橡胶本构模型参数，通过自适应方法对原始设计模型进行橡胶有限元仿真，得到完整的橡胶刚度曲线后，再对比与设计刚度曲线的差距，以此确定橡胶是否满足刚度要求，最后阶段才通过试验验证，减少试验周期。

1 橡胶有限元网格自适应方法

对于橡胶大变形结构件，为提高分析精度得到完整的橡胶刚度曲线，可以采用有限元自适应分析方法。自适应分析有3种方法：第一种是ALE自适应网格；第二种是自适应网格重划；第三种是网格间的求解变换[2-3]。基于3种自适应分析方法的各自优缺点，虽然都可以解决分析过程中网格畸变问题，但只有网格间的求解变换这种方法最适用于橡胶刚度的求解。它是用一个新的网格替代因变形过大而严重扭曲的原有网格，把原来的分析结果自动映射到新网格上，然后继续分析[3-4]。

自适应分析流程如图1所示：在多次网格重划中，直接选取最大位移为初始分析步的最大压缩量，但由于橡胶压缩到一定程度时，网格会严重扭曲而使分析失败；这时再根据最后迭代的分析步中计算出初始分析步的压缩量，然后再重新计算；接着，通过有限元软件把计算得到的结果文件转化为有限元模型，并生成几何进行重划，再把上一步计算得到的结果通过网格间的变换映射到新的网格中，同时要启动重启动分析步。重复上述过程，直至完成所需的压缩量为止。最后，把设计刚度曲线与仿真分析刚度曲线作对比，判定新的设计结构是否满足设计要求。两次重划分刚度曲线的连接点如果出现刚度点偏移的现象，必须在刚度偏移点的上一步骤重新选取更小的压缩量来加载，再重新分析。

2 橡胶本构模型

橡胶分子链由许多链节组成，其间多通过链节节点处化学交联而形成交联网络结构。链节一端节点到另一端节点的距离向量称为末端距向量。从分子或原子运动原理出发，采用统计法，通过对长链分子弹性性质的研究，可确定橡胶的宏观本构关系。由于组成原子的微布朗运动，橡胶长链分子可能有许多不同的构象。当没有外力作用，分子链的卷曲构象熵通常趋于最大值。当有外力作用时，分子链的构象改变，构象熵也发生变化。若分子链由3个长为t的链节组成，链节末端距向量为γ0，如果γ0=γ0≤nt，则可采用Gauss统计理论，建立材料本构模型。橡胶材料为超弹性材料，反映其应力、应变关系的模型称为本构模型。描述橡胶的本构模型种类相当多，均基于连续介质力学理论。19世纪以来，橡胶本构关系的理论模型很多且基本趋于成熟，主要分为3类：第一种是基于分子统计学理论的本构模型，如Arruda-Boyce模型、Kuhn-Grun模型；第二种是以应变不变量表示的应变能密度函数，其中比较有代表性的有Mooney-Rivlin，Yeoh；第三种是以主伸长率表示的应变能函数。不同的本构模型适用于不同的橡胶材料，所反映的精度也有一些差异[5-8]，因此选择合理的橡胶材料的本构模型尤为重要。

2.1 Mooney-Rivlin本构模型

对于各向同性材料，假如I1、I2和I3分别为右Cauchy Green变形张量C的第一、第二、第三基本不变量，对于初始无应力构形的超弹性材料，应变能函数W可表示如下：

W=W（I1，I2，I3）（1）

I=trC=C∶I=C（2）

I2=■[（trC）2-tr（C）2]=■（■-CijCji） i，j=1，2，3（3）

I3=detC（4）

C=F TF（5）

C=α？字÷？坠？锥（6）

J=deft（7）

式中，I为二阶张量不变量，F为变形梯度，х和Х分别为同一点在变形前后的坐标，J为变形后与变形前的体积比。

Tschoegl N W[11]认为，含高阶项的Mooney-Rivlin模型能更好地适应填充与非材料填充，Mooney-Rivlin[8]通过物质相变理论和大量试验，探讨了不可压缩各向同性超弹性材料有限变形弹性理论，假设单位体积的储能函数是右Cauchy Green变形张量的第一和第二基本不变量函数，建立的橡胶材料应变能函数如下：

W=C10（I1-3）+C01（I2-3）+■（U-1）2

式中，C10、C01、D1为待确定的描述橡胶力学特性的参数、U是应变能密度。

由于橡胶材料的几何和材料的双重非线性、超弹性及体积的不可压缩性，使得对任何橡胶制品进行有限元分析时都会遇到同样的难题：即大变形时计算的收敛性及接触问题的收敛性。这是目前在橡胶制品的有限元分析中存在的一个共同的难题。如果按常规的方法，即使画再小的网格，由于网格畸变，也不能完成所需的压缩量。橡胶本构模型的材料参数对有限元的计算结果影响很大，即使对于同一个本构模型，选择不同的材料参数也会得到不同的计算结果。参考文献[6]的实验结论：两参数Mooney-Rivlin模型能更好地模拟橡胶材料中等应变范围的应变能。Mooney-Rivlin模型能很好地描述橡胶变形在150%内的特性，研究文献广泛采用，因为它是两参数本构模型，系数拟合方面更方便。

3 应用实例

在某车的后悬缓冲块设计初期，已知后悬缓冲块的设计刚度曲线时，根据工程经验由具体的结构形式确定橡胶材料本构模型Mooney-Rivlin的2个参数。但如果用试验方法验证橡胶刚度，周期较长。通过对原始设计模型进行有限元分析，对比仿真刚度曲线与设计刚度曲线的差距，以此确定橡胶是否满足刚度要求，减少开发周期。

3.1 橡胶刚度试验原型机

橡胶刚度试验原型机主要由3个部分组成：橡胶底座、橡胶本体、试验机冲头（如图2所示）。

3.2 橡胶刚度试验物理模型的建立

建立橡胶材料的本构模型时，必须同时考虑橡胶材料的非线性和几何非线性。在ABAQUS有限元分析软件中，根据图3中的模型结构，假设橡胶底座及试验冲头不发生变形，建立缓冲块的三维有限元模型，网格大小为2～3 mm，橡胶底座及试验冲头建二维刚体单元，网格尺寸为5 mm（如图3所示）。检查橡胶和橡胶上下端面的固定和加载装置的干涉情况，以便计算结果更好地收敛。

本模型包含节点21 379个，单元18 476个，其中包括2 222个刚体单元、16 254个六面体单元。刚体单元的单元类为R3D3/R3D4，六面体采用杂交单元C3D8H（如图3所示）。根据实际情况约束底座全部自由度及冲头除压缩位置自由度，并在沖头控制点加载，同时主要部位的橡胶自接触模拟，与其他部位用绑定接触模拟。

确定Mooney-Rivlin本构模型参数。首先是对已有车型的多个同类缓冲块进行分析，对比仿真刚度与测设刚度的匹配情况，应用参数反求，通过设定目标值，反求出最佳材料参数，然后根据经验确定Mooney-Rivlin本构模型的2个参数。

3.3 自适应分析过程与重启动过程

初始阶段设定一个最大的橡胶压缩量A，提交abaqus计算分析，当网格发生畸变时，无法收敛而终止分析，这时通过计算最后迭代步与压缩总量的乘积来确定初始分析步的初始压缩量a1，再重新提交abaqus计算分析。查看初始分析步橡胶刚度曲线是否正常，变形网格如图4所示。建立好以上模型后，必须要设置重启动，为第二次分析做准备。

二次分析：提取初始分析最后分析步的变形后的结果网格（如图4所示），生成几何模型后重划，再重新划质量好的三维网格，并将模型的不同区域赋予与初始模型相同的材料及单元属性，也要设置重启动参数，为第三次重划网格分析做准备。然后设定一个剩余压缩量（A-a1），提交abaqus计算分析。当网格发生畸变时，无法收敛而终止分析，这时通过计算最后迭代步与压缩总量的乘积来确定二次分析步的初始压缩量a2，同时要让step1分析结果映射到二次分析step2中，使用关键字*MAP SOLUTION，再重新提交abaqus重启动分析。查看初始分析步橡胶刚度曲线和变形图是否正常（如图5所示），如果前后的刚度点基本重合，则说明二次分析完成，如刚度点出现大的偏移（如图6所示），则说明初始分析时单元出现了过大的网格畸变，需要施加一个更小的压缩量，接着重新进行计算，这是自适应方法的关键技术之一。如果橡胶变形过大，则需要3次甚至多次网格重划来完成整个分析，以得到一个较真实的计算结果（如图7、图8所示）。

3.4 结果讨论

本例通过3次网格自适应方法，即3次网格间的变化和重启动求解得到一条完整的橡胶刚度曲线。分析刚度曲线是由初始分析步刚度曲线、二次分析刚度曲线、三次分析步刚度曲线及四次分析步刚度曲线，总共4个分析步首位相连得到完整的刚度曲线（如图9所示）。由最终的分析刚度曲线与设计刚度曲线的对比结果看，本例的缓冲块基本满足设计刚度要求，可以采用此结构。

4 结论

（1）针对结构复杂的橡胶，在压缩到一定程度后无法计算的问题，采用自适应方法重划分网格继续计算，分析得到完整的橡胶刚度曲线，为橡胶新结构的设计阶段提供验证方法。

（2）本文简化橡胶台架实验模型，能有效地提高计算效率。

（3）该方法具有较强的工程实用性，可以进一步应用到橡胶的参数反求中，以缩短产品开发周期。

参 考 文 献

[1]刘萌，王青春，王国权.橡胶Mooney—Rivlin模型中材料常数的确定[J].橡胶工业，2011，58（4）：241-245.

[2]殷闻，勒晓雄，仝光.两种常用橡胶本构模型的有限元分析及其仿真[J].上海机电学院学报，2010（4）：215-218.

[3]赵腾伦，姚新军.abaqus在机械工程中的应用[M].北京：中国水利水电出版社，2007.

[4]刘展，祖景平，钱英莉.Abaqus基础教程与实例详解[M].北京：中国水利水电出版社，2007.

[5]朱艳峰，刘锋，黄小清，等.橡胶材料的本构模型[J].橡胶工业，2005，53（2）：119-125.

[6]李晓芳，杨晓翔.橡胶材料的超弹性本构模型[J].弹性体，200，15（1）：50-58.

[7]危银涛，杨挺.橡胶类材料大变形本构关系及其有限元方法[J].汽车技术，2011（11）：27-30.

[8]纪奎江，邓本诚，李俊山.丁腈橡胶加工和应用[M].北京：石油化学工业出版社，1978：19-21.

[9]Mooney MJ.A theory of large elastic deformation

[J].Journal of Applied Physics，1940，11（6）：582-592.

[10]Tschoel N W.Constitutive equations for elastom-ers[J].Journal of Polymer Science，1971（12）：959-970.

[责任编辑：陈泽琦]