超弹性本构模型对橡胶隔振器静态特性预测影响的研究

超弹性本构模型对橡胶隔振器静态特性预测影响的研究
王文涛;上官文斌;段小成
【摘要】The stress-strain data in uniaxial, planar and biaxial tension tests of standard rubber specimen are measured, based on which the parameters of 5 common hyperelastic constitutive models are identified with least square principle, and the fitting accuracies of different models are compared. Then taking the rubber isolator of a real vehicle as an example, its force-displacement curves under three typical engineering strains with different constitutive models are calculated and compared with measured data. The results show that in dealing with smaller strain values, Mooney-Rivlin model has the highest accuracy and computational stability among all the models, while in dealing with large strain values, Van der Waals model and 0gden3 model have better accuracies than other models.%通过对某橡胶材料标准试件的单轴拉伸、平面拉伸和等双轴拉伸的应力-应变数据的测试,并利用最小二乘原理识别了5种常见橡胶超弹性本构模型的参数,对比了不同模型间的拟合精度.以某款车型的橡胶隔振器为实例,计算不同本构模型在3种典型工程应变下的力-位移曲线,并与实测数据进行对比.结果表
明,MooneyRivlin模型在处理较小应变数据时具有较好的计算精度和计算稳定性,Van der Waals模型和3阶Ogden模型在处理较大应变数据时具有较好的计算精度.
【期刊名称】《汽车工程》
【年(卷),期】2012(034)006
【总页数】8页(P544-550,539)
【关键词】橡胶;超弹性本构模型;力-位移特性;计算精度与稳定性;模型参数识别【作者】王文涛;上官文斌;段小成
【作者单位】华南理工大学机械与汽车工程学院,广州510640;广州城市职业学院机电工程系,广州510405;华南理工大学机械与汽车工程学院,广州510640;宁波拓普集团股份有限公司,宁波315800;华南理工大学机械与汽车工程学院,广州510640;宁波拓普集团股份有限公司,宁波315800
【正文语种】中文
前言
橡胶材料是由长链、大分子和网状交连结构所构成的超弹性材料，具有良好的减振和隔振作用，因此广泛应用于汽车隔振器、轮胎和密封件等［1］。

橡胶材料具有很强的非线性特征，准确描述其力-位移的非线性性质并建立相应的数学模型有较大的难度，因此橡胶材料超弹性本构模型的恰当选取及其参数的识别引起了学术界和工业界的广泛关注。

目前，可用于橡胶隔振器力-位移分析的本构模型较多，常见的有Mooney-Rivlin(MR)模型、Van der Waals(VDW)模型等［2］。

在不同的激振幅值、激振频率和初始条件下，这些本构模型对橡胶件(如汽车隔振器等)的力-位移计算值与实测值的误差较大。

因此，研究橡胶材料本构模型在常见承载工况下的计算精度及其模型参数的识别方法对汽车隔振器的应力-应变分析具有重要的工程意义。

国内外对橡胶本构模型的研究主要集中于模型的推导或在某一应变载荷下的计算精度［3－4］。

文献［5］和文献［6］中研究了MR模型处理橡胶试件大变形的精
度问题，发现其具有较高的精度。

文献［7］中在轮胎滚动阻力研究中得出Yeoh
模型适应范围较广的结论。

文献［8］中基于单轴试验材料参数对比研究了Arruda-Boyce(A-B)模型、VDW模型和Yeoh模型的相对计算精度。

文献［9］
中研究了多项式模型及1阶、2阶和3阶Ogden模型的稳定计算区间问题。

上述研究尚未对不同本构模型在不同材料应变组合下的拟合精度做出探讨;在不同加载
条件下，不同橡胶本构模型在汽车橡胶隔振器静态预测计算中的计算精度研究及其材料应变大小组合选择的研究尚不多见。

本文中选取汽车橡胶隔振器静态预测计算广泛使用的MR模型、Yeoh模型、3阶Ogden(Ogden3)模型、A-B模型和VDW模型5种橡胶本构模型进行静态预测计算精度的研究。

为保证汽车橡胶隔振器静态预测计算精度，获取较为精确的橡胶材料本构模型参数是关键。

因此，进行了橡胶材料的单轴拉伸、平面拉伸和等双轴拉伸的工程应力-工程应变测试工作，利用最小二乘法识别不同本构模型的模型参数。

选取某车型橡胶隔振器进行力-位移静态特性和材料应变组合选择原则的研究。

分
析了橡胶隔振器在3种加载条件下实测与计算的力-位移数据，结果发现:橡胶材料本构模型的拟合精度与其计算精度的变化趋势基本一致;橡胶隔振器在力-位移计算中，为提高其计算精度，橡胶材料应变组合最大工程应变水平应≥计算试件最大工程应变水平;在所选工程应变范围内，MR模型的计算稳定性最好，计算值与实测
值误差最小，VDW模型与Ogden3模型在较大工程应变范围内也具有较好的计
算精度。

1 橡胶超弹性本构模型
就橡胶材料力学特性而言，其应力-应变关系较金属材料具有更强的非线性。

为描
述橡胶材料的应力-应变非线性特性，一般假设其外部载荷所做的功全部存储于弹
性体内，通常将反映其变形梯度与应变势能函数关系的模型称为超弹性本构模型［10］。

通过对应变能密度函数的应变不变量求导，可得橡胶材料的工程应力与
工程应变之间的本构关系［11］。

目前橡胶材料本构模型主要有基于分子统计理论和基于唯象理论两大类［12］。

上述两种理论的本构模型种类较多，在不同应变范围内的计算精度差别较大。

汽车橡胶隔振器受载工况复杂，故仅讨论在工程计算中应用较多、且在各应变区间计算精度不同的5种本构模型:基于唯象理论的MR模型、Yeoh模型与Ogden3模型和基于热力学或分子统计理论结构的A-B模型与VDW模型。

应变势能函数U 的一般形式［13］为
式中:I1、I2、J分别为1阶、2阶、3阶应变不变量，它们为3个主拉伸比的函数;C1、C2、…、Cm为m个超弹性材料的剪切特性常数;d1、d2、…、dn为n个超弹性材料的压缩特性常数。

I1、I2、J与超弹性材料的3个主拉伸比λ1、λ2、λ3的关系为
1.1 MR模型
MR模型的应变势能函数［14］为
式中:W为单位体积的应变势能;C10、C01和D1均为依赖于温度的模型参数，当D1=0时，该材料为完全不可压缩;J el为弹性体积应变比。

1.2 Yeoh模型
Yeoh模型的应变势能函数［15］为
式中:Ci0为依赖于温度的模型参数，i=1，2，3。

1.3 Ogden3模型
Ogden模型的应变势能函数［16］为
式中μi、αi均为描述温度函数的模型参数。

1.4 A-B模型
A-B模型又称八链模型［17］，其应变势能函数为
式中:C1=，C5=，这些模型参数通过热力学统计方法得到，具有明确的物理含义;λm为锁死伸长率;D=2/K，K为体积模量。

1.5 VDW模型
VDW模型的应变势能函数［6］为
式中为初始剪切模量;a为全局相互作用参数;β为不变量混合参数。

2 橡胶材料试验和模型参数识别方法
2.1 橡胶材料试验
若要精确的识别橡胶材料本构模型的参数，通常须进行单轴拉伸、平面拉伸和等双轴拉伸材料试验，平面拉伸对应纯剪切变形、等双轴拉伸对应单轴压缩［18］。

图1为3种材料应力-应变测试装置和加载方向。

根据汽车橡胶减振器的实际变形和选材情况，选取某种橡胶材料作为测试对象，测量标准试件在0.5、1和2的3种典型工程应变下的材料数据。

标准试件工程应变的计算式为
式中:εC为工程应变;l0为试件初始长度;l为试件承载后长度;Δl为试件伸长量。

因此0.5、1和2的工程应变分别相当于标准试件伸长量为原有试件长度的0.5倍、1倍和2倍。

图1 材料实验装置和加载方向
试片在初始几个加载循环时，其应力-应变曲线与稳定时区别较大，随着加载次数
的增加，同一应变下的应力最大值减小，这种应力减少现象称为马林斯效应(Mullin effect)，如图2所示。

为保证计算精度，工程计算中一般选用3次加载循环以后加载段的材料数据进行模型参数拟合［19］。

图2 实测橡胶材料工程应力-工程应变关系
2.2 本构模型参数识别方法
假设橡胶为不可压缩材料，对其势能函数在单轴拉伸、等双轴拉伸和平面拉伸3
种状态下的拉伸率求导，可求得对应拉伸状态理论工程应力［20］:
式中:λU、λB和λP分别为单轴、等双轴和平面拉伸的拉伸率;σthi为理论工程应力。

式(9)中不同拉伸状态下的应变不变量I1和I2分别为
利用最小二乘法拟合实测材料的应力-应变数据，进而求取唯一的本构模型参数。

本构模型的拟合误差为
式中:σtesti为实验测试应力;E为实测应力与计算工程应力的相对误差平方和。

3 橡胶本构模型参数识别与拟合精度
3.1 应力-应变组合
工程计算中一般利用橡胶材料的单轴、等双轴和平面拉伸实验数据识别其模型参数［4，21］。

对3种应变水平0.5、1和2的材料数据进行组合，如表1所示。

汽
车橡胶隔振器在实际应用中单轴压缩变形和纯剪切变形(相当于等双轴拉伸变形)相对较小，因此，在材料应变组合中以单轴拉伸工程应变最大，平面拉伸工程应变次之，等双轴拉伸工程应变最小。

表1 材料应变组合方式组合编号单轴拉伸最大工程应变等双轴拉伸最大工程应变
平面拉伸最大工程应变1# 0.5 0.5 0.5 2# 1 0.5 1 3# 2 1 2
表1中1#、2#和3#应变组合分别对应橡胶材料承载较小、中等和较大的应变状况。

3.2 模型参数及其拟合误差
表2～表6为利用实测工程应力-工程应变数据识别的模型参数和不同模型的拟合误差。

表2 Yeoh模型参数及拟合误差拟合参数1#材料应变组合拟合误差2#材料应变组合拟合误差3#材料应变组合拟合误差C10 0.404 0.351 0.422 C20 3.15%0.053 5.96%0.017 11.73%0.082 C30 0.037 0.014 0.002
表3 Ogden3模型参数及拟合误差拟合参数1#材料应变组合拟合误差2#材料应变组合拟合误差3#材料应变组合拟合误差μ1 2.861 0.698 3.817 μ2 0.001 0.005 2.976 μ3 2.072 0.000 0.000 α1 9.98%0.265 10.68%0.680 9.82%3.188 α2 12.878 7.262 3.706 α3 11.979 0.635 7.619
表4 A-B模型参数及拟合误差拟合参数1#材料应变组合拟合误差2#材料应变组合拟合误差3#材料应变组合拟合误差μ 0.7 0.6 0.786 λm 0.8 μ0 11.78%0.723 14.24%0.647 16.06%1099.32 2111.23 1211.98
表5 VDW模型参数及拟合误差常数 1#材料应变组合拟合误差2#材料应变组合拟合误差3#材料应变组合拟合误差μ 0.877 0.776 0.914 λm 3.741 5.015 3.842 a 9.19%0.941 10.15%0.665 10.37%0.196 D 0.951 β 0.190 0.133 0 0 0
表6 MR模型不同应变组合下拟合参数及拟合误差拟合参数1#材料应变组合拟合误差2#材料应变组合拟合误差3#材料应变组合拟合误差C10 0.352 0.290 0.278 C01 0.027 2.56%0.060 3.39%0.040 9.15%
综合分析各表可知，就拟合误差而言，最大工程应变为0.5的1#材料应变组合中MR和Yeoh模型的拟合精度较高，其拟合误差约为3%;其余3种本构模型的拟合误差均在10%左右。

最大工程应变为1的2#材料应变组合中，MR模型的拟合误差最小，约为
3.39%;Yeoh模型的拟合误差次之，约为6%;Ogden3模型和VDW模型的拟合误差约为10%。

最大工程应变为2的3#材料应变组合中，MR模型的拟合误差最小，为
9.15%;Ogden3模型和VDW模型的拟合误差均在10%左右，在3种应变组合中A-B模型的拟合误差都最大，约在12% ～16%的范围。

4 橡胶隔振器力-位移关系的分析
4.1 试验设备及加载条件
试验设备为MTS 831测试系统，可测量橡胶件的静态与动态的力-位移曲线，还可测量橡胶件的滞后角与蠕变等性能参数，如图3所示。

图3 MTS831测试系统
橡胶试件静态加载条件如表7所示。

试验加载形式为轴向位移，其最大工程应变0.38、1和1.56分别对应于车辆在较为平顺、较为颠簸和极为颠簸的路况行驶。

加载位置为其端面，为消除马林斯效应的影响，3次加载循环后通过力、位移传感器记录其端面响应试验数据。

表7 橡胶试件静态测试条件加载位移/mm最大工程应变测试数量预载次数加载速率/(mm/min)加载方式0～12 0.38 8 3 10轴向0～30 1 8 3 10轴向0～50 1.56 8 3 10轴向
4.2 橡胶隔振器力-位移计算和有限元分析
被测橡胶隔振器结构和尺寸如图4所示。

该试件两端的圆形金属块硫化在橡胶体内，可进行拉伸和轻微压缩测试。

其两端面均可作为加载端。

首先使用CATIA建立3D数模。

两端金属块视为刚度无穷大，利用hypermesh 软件绘制出如图5所示的有限元网格模型，共包含12 386个六面体单元。

最后，将网格模型导入Abaqus有限元分析软件，约束其任一端面的6个方向自由度，
在另一端轴向加载，测得圆弧凹点的拉伸工程应力和工程应变值。

图4 橡胶隔振器结构
图5 橡胶隔振器网格模型
首先将橡胶元件离散化;随之选择位移插值函数;然后进行单元分析，单元应变、单元应力关于位移参数的函数分别为
式中:f为单元体内任意一点的位移列阵;N为形函数矩阵;δe为单元节点列阵;ε为单元应变;B为几何矩阵;σ为单元应力;D为弹性矩阵;S为单元应力矩阵。

根据虚功原理求取单元体节点的力与位移的关系［14］式为
式中:F e为节点载荷;K e为刚度矩阵。

由于橡胶元件有很强的非线性特性，故不能简单使用式(16)求取其单元的应力-位移的关系。

通常使用式(3)～式(6)所示的超弹性体本构模型来求取它们之间的应力-应变关系。

进而使用式(17)和式(18)求橡胶元件加载载荷。

Abaqus软件中具有较多的非线性超弹性体本构模型，如MR模型和Ogden模型等，试验中橡胶隔振器的应变、位移采用该软件进行。

对橡胶隔振器加载3种工程应变0.38、1和1.56时，通过 Abaqus软件计算3种应变所对应的加载位移如表7所示。

4.3 本构模型计算精度分析
试验共测取8个试件的力-位移数据，间隔为0.1mm。

由于在较小位移间隔内橡胶隔振器的力-位移可视为线性关系，故可在整个加载位移范围内选取100个点进行线性插值，求取各点力-位移平均值，然后计算各本构模型的相对误差，即
式中:E Mean为平均相对误差;N为插值点数，N=100;为第i个插值点处的实测拉伸力为第i个插值点处的计算拉伸力。

4.3.1 不同材料应变组合对橡胶元件静态力-位移的计算精度影响
图6 各本构模型计算值与实测值的力-位移曲线
选取非线性较强的0～50mm位移加载下的实验数据进行分析。

图6为各本构模
型在0～50mm位移加载时的计算值与实测值的力-位移关系曲线。

由图6可知，在3种不同材料应变组合方式中，5种橡胶本构模型的力-位移曲线均不能重叠，
说明各模型力-位移计算值均受材料应变组合影响。

图6中3#应变组合的曲线均与实测曲线较为接近;1#应变组合的曲线与实测曲线偏离最远;2#应变组合的曲线居中。

说明在橡胶隔振器静态预测计算时材料应变应与测试对象的最大工程应变基本接近。

图6中，VDW模型、Yeoh模型和Ogden3模型的力-位移曲线较为分散，而
MR模型在3种材料应变组合中计算偏差均不很大，说明MR模型具有较为稳定
的计算能力，VDW模型、Yeoh模型、Ogden3模型的计算精度受材料应变组合
的影响较大。

因此，在下面各模型的计算精度对比研究中，选择材料应变组合最大工程应变与试件最大工程应变较为接近的组合进行橡胶材料的静态力-位移预测计算。

4.3.2 本构模型在某橡胶元件的计算精度分析
图7为3种加载位移下各模型的力-位移关系曲线。

由图7(a)可见，在0～12mm 的加载范围内(最大工程应变为0.38)，MR模型的力-位移曲线与实测曲线最为接近，其余模型均与实测曲线偏离较大，Yeoh模型偏离最远。

由图7(b)可见，在0～30mm的加载范围内(最大工程应变为1)，MR模型、VDW模型和Ogden3模型的力-位移曲线与实测曲线最为接近，Yeoh模型与A-
B模型的力-位移曲线与实测曲线偏离较远。

图7 3种加载位移下各模型的力-位移曲线
由图7(c)可见，在0～50mm的加载范围内(最大工程应变为1.56)，各模型的力-
位移曲线与实测曲线接近趋势与图7(b)中基本一致;其中MR模型、VDW模型和
Ogden3模型的计算曲线与实测曲线更接近。

另外，使用3种材料应变组合计算5种模型在0～0.38、0～1和0～1.56工程应变范围内橡胶隔振器力-位移的计算误差分别为:MR模型的相对计算误差为6.31%;在静态计算时MR模型所选材料应变组合的工程应变应大于或等于试件承载最大工程应变;其余4种模型在橡胶隔振器力-位移计算中，在0～0.38工程应变范围内均具有10%以上的相对误差;在0～1、0～1.56工程应变范围内，VDW模型和Ogden3模型的相对误差较小。

如以较小应变材料应变组合计算较大试件变形的力-位移，Ogden3模型和A-B模型将产生20%左右的相对误差。

5 结论
(1)不同橡胶材料应变组合对橡胶本构模型的计算精度影响较大。

在5种模型中，MR模型比其他模型具有更为稳定的计算能力，在汽车橡胶隔振器的力-位移计算中推荐优先使用。

(2)不同位移加载条件，各本构模型的力-位移计算精度差异较大。

在0～0.38工程应变范围内，MR模型的计算精度最高，其相对计算误差约为6.5%。

在0～1和0～1.56工程应变范围Yeoh模型、VDW模型的计算相对误差较小;在0～1工程应变范围内，Ogden3模型的力-位移计算精度最高;Ogden3模型、A-B模型用较小应变材料组合计算试件较大工程应变时将产生很大的相对计算误差。

(3)各橡胶本构模型的静态预测计算精度与实测计算精度排序趋势基本一致。

参考文献
［1］ Shangguan Wen-Bin．Engine Mounts and Power Train Mounting Systems:a Review［J］．International Journal of Vehicle Design，2009，49(4):237 －258．
［2］ Wadham-Gagnon M．Hyperelastic Modelling of Rubber Behaviour in Finite Element Software［D］:McGill University，2006:117 －118．
［3］吴杰，上官文斌，潘孝勇．采用超弹性-粘弹性-弹塑性本构模型的橡胶隔振器动态特性计算方法［J］．机械工程学报，2010，46(14):109－114．
［4］ Charlton D J，Yang J．A Review of Methods to Characterize Rubber Elastic Behaviour for Use in Finite Element Analysis［J］．Rubber Chemistry and Technology，1994，67:481 －503．
［5］ Mooney M A．A Theory of Large Elastic Deformation［J］．Journal of Applied Physics，1940(11):582 －592．
［6］ Arruda E M，Boyce M C．A Three-Dimensional Constitutive Model for the Large Stretch Behavior of Rubber Elastic Materials［J］．Journal of the Mechanics and Physics of Solids，1993，41(2):389－412．
［7］ Ghosh P，Saha A，Mukhopadhyay R．Prediction of Tyre Rolling Resistance Using FEA［J］．Constitutive Models for Rubber，2003(3):141－145．
［8］ Seibert D J，Schoche N．Direct Comparison of Some Recent Rubber Elasticity Models［J］．Rubber Chem Technol，2000，73．
［9］ Li-Rong W，Zhen-Hua L．Modeling Method of Constitutive Law of Rubber Hyperelasticity Based on Finite Element Simulations［M］．Vol 76 Akron，OH，ETATS-UNIS，American Chemical Society，2003．
［10］徐明．橡胶类超弹性材料非线性本构关系研究和有限元分析［D］．北京:北京航空航天大学，2002:78－79．
［11］郭仲衡．非线性弹性理论［M］．北京:科学出版社，1980:68－69．［12］ Guo Z，Sluys L J．Constitutive Modelling of Hyperelastic Rubberlike Materials［J］．HERON，2008，53:109 －132．
［13］ Treloar L R G．The Physics of Rubber Elasticity［M］．Oxford U-
niversity Press，1975:150 －153．
［14］ Manual．Inc A T．Theory Manual［M］．Hibbit，Karlsson ＆Sorensen(HKS)2009:326－327．
［15］ Yeoh O H．Some Forms of Strain Energy Function for Rubber ［J］．Rubber Chemistry and Technology，1993，66:754 －771．
［16］ Ogden R W．Large Deformation Isotropic Elasticity-on the Correlation of Theory and Experiment for Incompressible Rubberlike Solids ［J］．Proceedings of the Royal Society of London Series A，Mathematical and Physical Sciences，1972，326(1567):565 －584．
［17］ Ling Y．Constitutive Models，Load-deflection Relations，and Instabilities of Elastomeric Materials and Structures［D］:State University of New York at Binghamto，1992:132 －134．
［18］杨晓翔．非线性橡胶材料的有限单元法［M］．北京:石油化工出版社，1999:121－123．
［19］ Dolwichai P，Limtragool J．Hyperelastic Material Models for Finite Element Analysis with Commercial Rubber［C］．Technology and Innovation for Sustainable Development Conference，2006:770－774．［20］ Shangguan W-B，et al．Finite Element Analysis of Static Elastic Characteristics of the Rubber Isolators in Automotive Dynamic Systems ［J］．SAE，2003，112:185 － 193．
［21］ Gent A．A New Constitutive Relation for Rubber［J］．Rubber Chemistry and Technology Journal，1996，69:59 －61．。