橡胶材料的基本实验及本构关系模型

第3章：橡胶材料的基础实验及本构模型
作为一种具有良好弹性性能的工程材料，硫化橡胶早在19世纪就被广泛应用于密封、承载、减振降噪等工业领域。

而橡胶轨道减振器的使用则是最近20年来的事情，然而，不同于金属材料仅需要几个参数描述其材料特性，橡胶的行为复杂，材料本构关系是非线性的。

它的力学行为对温度，环境，应变历史，加载的速率都非常敏感，这样使得描述橡胶的行为变得更为复杂。

而橡胶的制造工艺和成分也对橡胶力学性能有显著的影响。

简单依赖单向拉伸性能实验并不能完全描述材料包括压缩及剪切在内的所有力学行为，这也意味着对橡胶轨道减振器进行有限元分析和结构模拟，必须对橡胶材料进行包括拉伸、压缩，剪切及体积实验等在内的全部基础实验。

3.1 橡胶基础实验简介
描述橡胶材料的基础实验有8种（如图3-1）：单轴拉伸和压缩实验，双轴拉伸和压缩实验，平面拉伸和压缩（纯剪）实验以及测定体积变化的实验（拉或压）。

在长期的研究和实验，发现从单轴拉伸，双轴拉伸，平面拉伸及体积压缩实验中能够获得足够精确的实验数据。

因此，目前国际上定义橡胶材料力学行为的实验为：单向拉伸、双向拉伸、平面剪切及体积压缩。

图3-1 橡胶材料的8种基础实验
对有限元分析所用的实验数据，一个重要的要求是，实验时实验试样应能达到“纯”的应变状态，这样得到的应力应变曲线是我们期望的能代表橡胶的行为特性的状态。

有限元程序通常需要输入的应力应变实验数据范围应大于要分析结构的预期的最大应力应变范围。

通常，理想状态应该是测得在几种准静态荷载模式下的应力应变曲线，这样可以选择出最合适的材料的本构模型以及反映这种模型的参数。

图3-2是本课题研究工作中所用到的一组橡胶材料数据，该实验在美国AXEL实验室完成，材料是公司生产轨道减振器产品所用配方。

图3-2 橡胶基础实验数据
3.2 橡胶材料的基础实验
3.2.1单轴拉伸实验
单轴拉伸实验是最常用到的一种实验，有很多种橡胶拉伸的实验标准。

但是为有限元分析的实验要求比标准的实验方法还要高些，最为明显的是实验要达到一个纯的拉伸状态，也就是实验应该尽量减小对试样侧面的约束。

可用有限元技术分析确定出满足设计要求的合适的试样的长和宽度的比例。

图3-3 单轴拉伸实验
有限元计算的结果显示长度方向的尺寸至少是厚度和宽度的10倍。

并且试样也没有必要制作成哑铃形状，简单的正方形长条试样就可获得很好的实验效果，实验时，为了测得纯的单向拉伸的应力应变曲线，测量的位置应远离夹持部分，使用不接触式视频或激光应变传感器能够获得精确的满足要求的实验数据。

3.2.2平面拉伸实验
图3-4 纯剪实验
如图3-4所示的纯剪实验与拉伸实验类似，实际是一个比较宽的拉伸实验，
有限元分析试件的几何特征表明，实验精度对试件长和宽的比例非常敏感，当试件的宽度大于其长度的10倍时，可以获得精确的实验数据。

由于橡胶材料几乎不可压缩，因此实验时，在与拉伸方向成45度的地方出现了纯剪状态。

试件宽度要远大于拉伸方向的长度，试件的厚度方向约束尽量减小，让试件在厚度方向可以缩小。

有限元分析结果表明：在远离夹具的地方，特别是试件中心(如a区)，形成了“纯”的剪切状态，因此，测量应变的部位应该远离夹具并且尽量避免在试样边缘，应变的测量一定的使用非接触的传感器，视频或激光应变传感器都是很好的测试设备。

图3-5：剪切实验的应变分布
3.2.3等轴拉伸实验
不同于金属材料，橡胶的拉伸和压缩的应力应变曲线差别很大，仅有拉伸实验是不足以很好地确定橡胶本构模型。

因此，有限元分析的实验数据必须包括能够反映压缩性能的实验数据，然而，在单轴压缩实验（图3-7）中“纯”的单向压缩状态很难实现，一方面是因为在压缩时，试件和压缩装置之间的摩擦对应力状态的影响非常大。

试件的侧面在压缩时不可能是自由的膨胀。

即使摩擦系数很小也会引起侧面出现显著的剪切应变，而且最大的剪应变往往可能超过其最大
的压缩应变。

另一方面则是由于无法确定摩擦系数，实验数据无法修正，因此利用单向压缩实验无法得到合适的实验数据。

对于不可压缩材料，把静水压力（拉力）叠加在另外一种单轴压缩应力状态上，得到一种新的应力状态，即等轴拉伸应力状态，但是材料的应变状态不会发生变化，如下图：
图3-6：应力状态的叠加
所以，等轴拉伸实验和单轴压缩实验是等价的。

他们之间的应力应变关系为：
1)1/(1)1(23
−+=+=b c b c c εεεσσ （3-1） 其中b ε为等轴拉伸状态工程应变，b σ为等轴拉伸状态工程应力,c ε为压缩实验工种应变，c σ为压缩状态工程应力。

图3-7： 压缩实验 图3-8：等轴拉伸 因此，对于不可压缩或者几乎不可压缩橡胶等轴拉伸实验等价于单轴压缩实验。

尽管这种实验要比单轴压缩实验复杂得多，但是这种实验可以得到很“纯”的应力状态，实验的结果可以让我们得到精确的橡胶实验数据。

双轴应变状态可以通过对一个橡胶圆盘进行径向拉伸得到，有限元分析结果表明：在圆盘中心区域形成了均匀的双向应变云图（如图3-8），因此利用橡胶圆盘试样获得双向拉伸应变是可行的。

实际可行的装置如图3-9，应变测量同样使用的是非接触的传感器。

图3-9： 双轴实验装置
3.2.4体积压缩实验
橡胶材料几乎不可压缩，若计算中考虑橡胶材料的可压缩性时，需要进行体积压缩实验。

实验时把一个圆柱形试件放置在固定的刚性容器中进行压缩（如图3-10），从实验中可直接获得工程应力应变曲线，其初始斜率即为体积模量，一般情况下，这个值比剪切模量高2～3个数量级。

本课题所使用的橡胶材料的体积压缩实验数据见图3-10所示，其初始体积模量为2.61GPa。

图3-10 体积压缩实验及实验曲线
3.3 橡胶材料的本构模型及系数定义
常用的对橡胶力学性能的描述方法主要分为两类：一类是认为橡胶为连续介质的现象学描述；另一类是基于热力学统计的方法。

现象学的描述方法假设在未变形状态下橡胶为各向同性材料，即长分子链方向在橡胶中是随机分布的。

这种各向同性的假设是用单位体积（弹性）应变能密度来描述橡胶特性的。

而基于统计热力学方法的理论则认为：观察到橡胶中的弹性恢复力主要来自橡胶中熵的减少，熵的减少是由于橡胶的伸长使得橡胶结构由高度的无序变得有序。

由对橡胶中分子链的长度，方向以及结构的统计得到橡胶的本构关系。

经过长期的理论研究，到目前为止，较成熟的本构模型，有基于连续介质力学理论的多项式形式模型和Ogden 形式模型以及后来基于热力学统计理论的模型。

按此标准划分的两类模型如下表所示：
z 基于热力学统计的本构模型
1. Arruda-Boyce 模型
2. Van der Waals 模型
z 现象学模型：（Phenomenological models）
1. N 次多项式（Polynomial）模型
2. Ogden 模型 材料系数 2 4 2N 2N
表3-1：有限元分析中的两类橡胶本构模型
3.3.1多项式形式
对于各向同性材料，应变能加法分解成应变偏量能和体积应变能两部分，形式如下：
)1()3,3(21−+−−=J g I I f U
（3-2） 令∑=−=N i i i
J D g 12)1(1，并且进行泰勒展开，可以得到下式： ∑∑=+=−+−−=
N j i N i i i j i ij J D I I C U 11221)1(1)3()3( （3-3）
这种形式为多项式表示的应变能，参数N 为我们选择的多项式阶数。

D i 的
值决定材料是否可压：如所有D i 都为0，说明材料是完全不可压的。

对于多项式
模型，无论N 值如何，初始的剪切膜量0µ，初始的体积膜量0k ，都决定于多项式一阶（N=1）系数：
10011002 ),(2D k C C =+=µ （3-4）
对于完全多项式，如果1=N 则只有线性部分的应变能量保留下来，即Mooney-Rivlin 形式：
21201110)1(1)3()3(−+−+−=J D I C I C U （3-5）
这种形式是最简单的超弹模型，未知精确参数的材料通常使用这种形式。

这种函数不能够表示应力应变曲线的大应变部分的“陡升(upturn)” 行为，但是在小应变和中等应变时是可以较好的模拟材料特性的。

图3-11是ABAQUS 用Mooney-Rivlin
本构模型拟合的单轴拉伸实验数据的结果。

图3-11：M-R 模型拟合 图3-12：Yeoh 模型拟合
多项式模型的特殊形式可以由设定参数为0来得到。

如果所有0=ij C （0≠j ），则得到减缩的多项式模型：
i N i N i i i
i J D I C U 21010)1(1)3(−+−=∑∑== （3-6）
当3N =，则减缩多项式为Yeoh 形式[5]，它是减缩多项式的特殊形式：
i i i i i
i J D I C U 2313010)1(1)3(−+−=∑∑== （3-7）
此时，如果)(1O C 10=，第二个系数将为负，并且比第一个系数小1－2个数量级；例如， )01.0()1.0(20O O C −−～的值为，30C 的值将再小1－2个数量级，
并且为正，)40.1()20.1(30−+−+E O E O C ～的值就为，这种量级关系将产生典型的S 形的橡胶的应力－应变曲线。

在小变形的时候，10C 代表初始剪切膜量；在中等变形情况时，由于第二个系数20C 为负值，曲线软化；在大变形情况时，由于第三个系数30C 为正，曲线又变硬。

因此Yeoh 形式具有良好的拟合橡胶大变形的能力。

利用给出的这组实验数据，材料的参数可由最小二乘法决定，这样可以使得误差最小。

对于n 组应力应变的实验数据，相对误差E 取最小值。

()211∑=−=n i test i th i T E
（3-8）
其中，test i T 是实验数据中的应力值，th i T 是按照本构关系这种伸长率对应的应力表达式。

多项式形式的势能关于ij C 是线性的；因此使用线性的最小二乘法（对于Ogden，Arruda-Boyce 以及Van der Waals 形式的势能关于ij C 是非线性的，需要用到非线性最小二乘法）。

对于完全的多项式模型我们可以把th
i T 写成如下的形式：
,...1 ,)(1n k X C T N j i k ij ij th i ==∑=+λ （3-9）
其中)(k ij X λ的函数表达与应力状态有关（单向拉伸，双轴，纯剪），前面讨论过，1=N 是一阶多项式（即Mooney-Rivlin 形式）。

当橡胶材料视为完全不可压时， N=3为缩减三阶多项式（即Yeoh 形式）为了使相对误差最小，令：
0=∂∂ij
C E （3-10） 然后得到下面的一组方程，共有)(3N N 2
1M +=个： ....1 )()()
()(1112N m l T X C T X X n k test k
k lm ij n k N j i test k k lm k ij =+=∑∑∑===+λλλ （3-11） 解这个线性方程组可以得到确定的系数ij C 。

本构方程
弹性系数：C ij Mooney-Rivlin C 01= 4.826E-02
C 10= 0.579
Yeoh C 10= 0.696
C 20=-9.147E-02
C 30= 1.724E-02
表3-2：拟合的材料参数
3.3.2 Ogden 形式
Ogden 应变能以三个主伸长率1λ，2λ，3λ为变量。

在有限元软件中应变能的形式如下：
i N i i N i i i J D U i i i 2113212
)1(1)3(2∑∑==−+−++=αααλλλαµ （3-12） 其中132131=→=−λλλλλi i J 。

Ogden 应变能函数的第一部分只与1I 和2I 有关。

Ogden 应变能函数不能写成仅由1I 和2I 表示的形式。

如果2a 21N 21−===,,α，得到Mooney-Rivlin 模型，在Ogden 模型中，0µ由全部系数决定：
∑==N i i 1
0µµ （3-13）
对于Ogden，Arruda-Boyce 和 Van der Waals，其形式势能关于他们的系数是非线性的；因此需要一个非线性的最小二乘法。

我们使用Marquard-Levenberg 运算法则[13]。

m a i i ....1,=，是超弹性材料的系数，其中m 是本构模型中待定系数的个数。

例如，在Ogden 模型中N m 2=（在Arruda-Boyce 模型中2m =，在Van der Waals 模型中4m =）。

系数在如下的反复迭代中得到：
[])()(111)()()()1(r k r jk m j n k ij
r jk r ik r i r i E P P P a a −==+∑∑+−=γδ （3-14）
其中r 是反复迭代次数，n 是实验数据的组数。

k E 是相对误差：
test
k
th
k test k k T T T E −= （3-15） ik P 是相对误差关于系数i a 的导数。

i
th
k test
k i k ik a T T a E P ∂∂−=∂∂=1 （3-16） 用线性最小二乘法计算出)0(i µ，初始化系数)0(i α后，用上面的推导方法，并用到下面的公式：
)(2
11−−−=∂∂i i c i
i th k T ααλλαµ （3-17） λλλαµ
λλαµαααααln )(2)(21111−−−−−−−−=∂∂i i i i c i i c i
i i th k c T （3-18）
⎪⎩
⎪
⎨⎧−−−=planar if ,1biaxia if ,2uniaxial if ,21c （3-19）
即可求出Ogden 模型的系数，本课题的橡胶实验数据用Ogden 形式N=3，N=4进行了拟合，其拟合的橡胶参数见表3-3所示。

本构方程 系数mu 系数alpha Ogden,N=3
mu 1=1.257 mu 2=4.081E-04 mu 3=5.611E-05
alpha 1= 1.0985 alpha 2= 11.983 alpha 3=-10.055 Ogden,N=4
mu 1= 47.646 mu 2= 4.883e-03 mu 3=-82.364 mu 4= 36.078
alpha 1=-1.550 alpha 2= 9.306 alpha 3=-1.825 Alpha 4=-2.082
Ogden 形式拟合的参数
图3-13、3-14是ABAQUS 用3次Ogden 本构模型、4次Ogden 本构模型拟合实验数据的结果。

从拟合结果可知，实验曲线与预测曲线基础吻合，但对平面曲线拟合出现了一定的误差，因此，如果用这组数据去分析受拉、压工况下的橡胶产品是合适的，但用此参数对轨道减振器产品进行受力分析会产生较大误差，应选择更好的模型。

图1-13： Ogden N=3 图1-14： Ogden N=4
3.3.3 Arruda-Boyce 形式
Arruda-Boyce 形式[6]应变能定义如下：
⎟⎟⎠
⎞⎜⎜⎝⎛−−+−=∑
=−J J D I C U i i
i
i m
i
ln 211)3(25
1
1
2
2λ
µ （3-10）
其中：673750
519C 705019C 105011C 201C 21C 54321=====
,,,, Arruda-Boyce 模型也叫做eight-chain 模型。

51...C C 的值由热力学统计方法得到，是具有物理意义的。

系数µ代表初始的剪切膜量，即0µ。

系数m λ表示锁死应变（locking stretch），位置大约在应力应变曲线最陡的地方。

初始的体积膜量为D 2K 0
=。

因为只有两个参数，这样哪怕只有很少的材
料行为已知，材料的本构关系也可以得到。

Kaliske 和Rothert(1997)成功地证明这种体积应变能表达式应用于大部分工程弹性材料都足够精确。

图1-15是ABAQUS 用Arruda-Boyce 模型拟和的一组实验数据。

Arruda-Boyce 模型的材料参数只有两个，改变这两个参数的值只能改变应力应变的比例而无法改变曲线的形状，由于试验数据和Arruda-Boyce 模型预测的曲线形状不同，因此使用Arruda-Boyce 本构模型就不能很好地模拟这组材料的特性。

图1-5 Arruda-Boyce 形式本构模型
3.3.4 V an der Waals 形式
Van der Waals 模型定义的应变能为：
[]⎟⎟⎠⎞⎜⎜⎝⎛−−+⎪⎭
⎪⎬⎫⎪⎩⎪⎨⎧⎟⎟⎠⎞⎜⎜⎝⎛−−+−−−=J J D I U m ln 21123~32)1ln()3(223
2αηηλµ （3-11） 其中：21)1(~I I I ββ+−=，参数β是把1I 和2I 混合成I ~
用到的线性的参数，
3
3
~2
−−=
m I λη。

可以看到这种应变能有4个独立的系数。

由图1-6可以看到，Van der Waals 模型比Arruda-Boyce 模型有更好地拟
合实验数据的能力。

它不但可以改变应力应变的比例，而且可以改变拟合曲线的形状。

而且对压剪应变状态具有良好的拟合预测能力，因此本课题所使用的本构模型为Van der Waals 形式，其拟合的材料参数见下表3-4所示。

mu Lambda_m A Beta 1.651
3.743
1.103
7.59e-02
Van der Waals 的拟合参数
图1-6 Van der Waals形式本构模型。