橡胶材料作为一种工程减振材料，广泛应用于汽车减振系统(如汽车悬置、衬套等)。相比于传统的燃油汽车，电动汽车具有高频激励的特点，因此，橡胶隔振器的高频动态特性逐渐成为人们研究的重点。目前，对于橡胶隔振器，其静态性能计算方法已较为成熟[1-3]，而对橡胶隔振器动态特性计算方法的研究则相对较少。
橡胶隔振器的动态特性与预载、激振振幅和激振频率等相关[4-5]。为了描述橡胶的动态特性，前人提出了很多不同的模型，如Maxwell模型、Kelvin-Voigt模型、BERG模型、DZIERZEK模型、修正的边界曲面模型[6-9]等，或者以它们的不同组合而得到的模型来表征橡胶隔振器的动态特性[6-9]，但是少有研究提出可行的办法将材料的力学行为特性模型转化为有限元软件中的材料本构模型。另外，新发展的分数导数模型能以较少的参数数量有效反映橡胶衬套动态特性的频率影响[10-15]，吴杰、上官文斌[16]将粘弹性分数导数模型用于橡胶隔振器动态特性的建模及应用上，取得了一定的成果。但分数导数模型当前时刻计算结果依赖于前面所有结果，不利于模型的实时仿真。因此，有必要研究橡胶的材料特性，求出其本构参数，并将其转换为有限元软件中的材料参数。
本文将基于橡胶材料实验，揭示橡胶材料的本构关系。采用由Mooney-Rivlin和多个Maxwell模型叠加组成的非线性粘弹性模型，使用多种实验方法对模型的本构参数进行拟合，并将其转换为有限元软件中的材料参数，使用跨点动刚度有限元模型，对某款橡胶悬置的动态特性进行有限元分析。此外，文中还建立了橡胶隔振器等效力学模型，分析了原点动刚度和跨点动刚度的区别；同时，文中搭建了橡胶隔振器有限元模型，分别分析其跨点动刚度和原点动刚度，并与实验结果进行对比。最后，文中还简单对比了三种拟合粘弹性参数方法的优缺点。
1 橡胶隔振器静动态特性
1.1 原点动特性与跨点动特性
对于橡胶隔振器，常使用动刚度和滞后角来描述其动特性；根据测试时力传感器布置位置不同，将其测得的动特性分为原点动特性和跨点动特性。若力传感器布置在输入端，即力取值点与位移输入点为同一点，则计算所得结果为原点动特性，通常所说的动刚度和滞后角为原点动刚度和原点滞后角；若力传感器安装在输出端，则此时测得的动刚度和滞后角为跨点动刚度与跨点滞后角。原点动特性与跨点动特性是从不同的角度描述橡胶隔振器的动态特性，它们之间的关系取决于橡胶隔振器的惯性质量和试验激振频率，即橡胶隔振器所产生的惯性力。
橡胶隔振器的等效力学模型如图1所示，当试验激振频率较低时，由于橡胶隔振器的质量较小，此时惯性力较小，因此，此时原点动特性与跨点动特性的区别可忽略不计。但是，在高频振动激励下，惯性力的影响不能忽略不计，此时原点动特性与跨点动特性差别较大，不能再将二者等同。
图1 橡胶隔振器等效力学模型
Fig.1 Equivalent mechanical model of rubber isolator
1.1.1 动特性测试中惯性力影响分析
在原点动特性测试中，其力传感器布置在输入端，力传感器随激励端的运动而运动，同时输入端连接件、测试件工装等也随激励端的运动而运动，因此，输入力中既包含激励端的作动力F，又包含传感器、连接件、测试件工装等附加质量产生的附加惯性力Fm。因附加惯性力Fm并不引起弹性隔振元件的弹性变形，所以必须在输入力Fi中扣除附加惯性力Fm才能获得作动力F。
图2是原点动特性试验中输入力Fi、作动力F和附加惯性力Fm的矢量关系图，其中X=Asin(ωt)为位移输入。其中，作动力幅值
可表示为：
式中：α为输入力与输入位移间的夹角；θ为作动力与输入位移的夹角，即原点滞后角。
图2 输入力、作动力与附加惯性力的矢量关系图
Fig.2 Vector diagram of input force, motive force and additional inertia force
附加惯性力Fm可由附加惯性质量和振动加速度求得。附加惯性质量包括传感器、连接件、测试件工装等的质量。由于原点动特性可通过作动力计算得到，因此，需要根据式(1)对输入力进行修正，扣除附加惯性力Fm成分，得到作动力，但由于附加惯性力与附加惯性质量和测试频率有关，且对于不同的测试件，其工装不同，即附加惯性质量也会不同，此时，附加惯性力的扣除将非常复杂繁琐。因此，对于橡胶隔振器的高频动特性试验，建议使用跨点动刚度测试法，将力传感器布置在输出端，可消除惯性力对测试结果的影响。
1.2 橡胶隔振器有限元模型
如图3所示的橡胶悬置有限元模型，其橡胶外管刚性连接至结点1，约束此结点作为边界条件；橡胶内管关联至中间结点2，以该结点作为加载点；在加载点施加沿z方向的简谐位移激励，若采集加载点(结点2)处反力来计算动刚度，则可得到原点动刚度；若取橡胶外管连接点(结点1)处反力来计算动刚度，则可得到跨点动刚度。图4为动刚度模型等效力学示意图。
图3 橡胶悬置有限元模型
Fig.3 Finite element model of rubber mount
图4 橡胶悬置动刚度模型等效示意图
Fig.4 Equivalent schematic diagram of dynamic stiffness model of rubber mount
对于简单简谐动态实验，施加一小振幅简谐位移激励，则加载位移、加速度可表示如下：
式中：A为激励振幅；ω为激励频率；对图4中的等效质量m进行受力分析有：
式中：
均为复刚度；m为橡胶悬置等效质量；
分别为弹簧1、弹簧2对应的形变量。
在低频阶段，随频率ω逐渐增大，惯性力
快速增大；但在低频阶段，随频率ω逐渐增大，橡胶悬置的刚度
变化不大；因此，由式(5)、式(6)可知，随频率ω逐渐增大，
减小，
增大，即如图5所示，原点动刚度随频率(0 Hz~200 Hz)的增大而减小，相反，跨点动刚度随频率(0 Hz~200 Hz)的增大而增大。
图5 原点动刚度与跨点动刚度仿真结果对比
Fig.5 Comparison of analysis results between origin and cross-point dynamic stiffnesses
1.3 橡胶隔振器静动态实验结果分析
对图3所示的橡胶悬置进行静动态特性测试，其静动态实验在MTS 831液压伺服振动试验台上进行，实验采用固定衬套外圈、激励内圈的方式。图6为±8000 N静态加载-卸载过程得到的力-位移曲线，图7为不同振幅激励下，橡胶悬置动刚度随频率的变化关系。由图6可知，静态加载时，橡胶材料的应力和应变不再是线性对应关系；由图7可知，当激励振幅较小时，橡胶悬置的动刚度几乎不随振幅的变化而变化；振幅一定时，橡胶悬置的动刚度随频率的变化而变化。
图6 ±8000 N静态加载-卸载过程力-位移曲线
Fig.6 Force-displacement curve during ±8000 N static loading and unloading
图7 不同激励振幅下动刚度随频率的变化关系
Fig.7 Dynamic stiffness under different excitation amplitudes
综上所述，在激励振幅较小时，激励频率是影响橡胶悬置动态特性的决定性因素，即橡胶悬置的动态特性表现出频率相关性；而激励振幅对橡胶悬置动态特性的影响可以忽略。同时，橡胶的静态加载具有强的非线性。
2 橡胶材料的叠加非线性本构模型及本构参数识别
为了更准确地描述橡胶悬置的静动态特性，根据橡胶悬置在试验中表现出的非线性和频率相关性，本文采用多个Maxwell模型来描述其频率相关性，即粘弹性；每个Maxwell模型由弹性单元(
)与粘性单元(tr_i)串联组成。由于橡胶表现出了很强的非线性，因此，Mooney-Rivlin模型被用于模拟橡胶的弹性变形部分，其中C10、C01均为材料模型参数；Mooney-Rivlin模型和多个Maxwell模型叠加组成了如图8所示的非线性粘弹性本构模型。
图8 MooneyRivlin-Maxwell叠加非线性粘弹性模型
Fig.8 MooneyRivlin-Maxwell nonlinear viscoelastic model
2.1 Mooney-Rivlin单元
对于Mooney-Rivlin模型，其应变势能可表示为：
式中：D1为材料模型参数；Jel为弹性体积比，假设橡胶为体积不可压缩材料；
分别为第1阶、第2阶应变不变量，其表达式为：
式中，偏伸长率
为总体积比，λi为主伸长率。初始剪切模量μ0、体积模量K0由式(10)给出：
材料参数C10、C01可由单轴拉伸、等双轴拉伸和平面拉伸试验拟合得到，具体试验拟合方法见文献[17-18]；本文拟合得到的Mooney-Rivlin本构参数如表1所示。
表1 拟合得到的Mooney-Rivlin本构参数
Table 1 Mooney-Rivlin constitutive parameters obtained by fitting
Mooney-Rivlin本构参数 C10/MPa C01/MPa拟合数值 0.36163 0.099851
2.2 Maxwell粘弹性单元
在t时刻，粘弹性单元的总应力为：
式中：γ (s)为剪切应变；GR(t)为剪切模量；τve(t)为粘弹性应力，归一化剪切松弛函数gR(t)可表示为：
式中：
为需要进行识别的本构参数，分别表示各个Maxwell粘弹性单元对应的松弛时间和弹性模量比例系数；GR(0)为瞬时剪切模量；
为各个Maxwell单元对应的弹性模量。
2.2.1 在时域范围进行粘弹性参数识别：通过橡胶松弛实验识别
橡胶试片松弛试验的试验工装和试验装置如图9所示。其中，橡胶试片长a =20 mm、宽b =20 mm、厚度δ=5 mm，试验装置为英斯特朗材料试验机，试验时，对试样施加γ= 0.2的剪切变形，保持形变γ不变，记录下工装处反力F(t)随时间的变化关系，实验结果如图10所示。
图9 松弛/蠕变试验装置及工装
Fig.9 Relaxation/creep experimental device and test tooling
图10 松弛试验结果
Fig.10 Relaxation test results
在t时刻，粘弹性单元的总应力τve(t)为：
式中：γ(s)为松弛应变；GR(t)为松弛剪切模量。
其中，归一化剪切松弛函数gR(t)可表示为：
式中，GR(0)为瞬时剪切模量，其中：
式中：F(t)为松弛实验测得的力随时间的变化关系；2A为橡胶试片总剪切面积。则归一化剪切松弛函数gR(t)可表示为：
因此，可得到归一化剪切松弛函数gR(t)随时间的变化关系如图11所示；另外，归一化剪切松弛函数gR(t)可表示为：
对式(19)进行拟合可得到粘弹性本构参数(参数1)，拟合得到的本构参数具体数值如表2所示。
图11 归一化剪切松弛函数随时间变化的关系
Fig.11 Relation of normalized shear relaxation function with time
表2 松弛实验拟合得到的粘弹性本构参数(参数1)
Table 2 Viscoelastic constitutive parameters obtained by relaxation experiment fitting (parameter 1)
本构参数ve 1/MPa gg2/MPa ve g tr_1/s tr_2/str_3/s 3/MPa ve数值 0.0472 0.0379 0.0675 430.4 994276294
tr_i、
分别为各个Maxwell粘弹性单元对应的松弛时间和弹性模量比例系数，则各个Maxwell单元对应的弹性模量
可表示为：
2.2.2 在时域范围进行粘弹性参数识别：通过橡胶蠕变实验识别
橡胶试片蠕变试验的试验工装和试验装置如图9所示，试验方法为对试样施加F=65 N的剪切力，保持剪切力不变，测试样的剪切形变γ(t)随时间t的变化关系，实验结果如图12所示。
图12 蠕变试验结果
Fig.12 Creep test results
其中，在t时刻，剪切柔度Js(t)为：
式中：τ0为蠕变实验所加的剪切应力常量；γ(t)为剪切形变随时间的关系。则归一化剪切柔度js(t)可表示为：
因此，可得到归一化剪切柔度js(t)随时间变化的关系，如图13所示。
图13 归一化剪切柔度随时间变化的关系
Fig.13 Relation of normalized shear compliance with time
通过卷积积分，将蠕变数据转换为松弛数据：
根据式(19)、式(22)、式(23)，可拟合得到粘弹性本构参数(参数2)，具体数值如表3所示。
表3 蠕变实验拟合得到的Maxwell单元本构参数(参数2)
Table 3 Constitutive parameters of Maxwell elements obtained by creep experiment fitting (parameter 2)
本构参数 gve 1/MPa ge 2v/MPa gve 3/MPa 1/s t_ 2/s r t_3/s r t_r数值 0.03850 0.03054 0.07258 421.8 10 42678 177
2.2.3 在频域范围进行粘弹性参数识别：通过橡胶试片动刚度和滞后角识别
在进行简单剪切动态实验时，对试件施加一小振幅简谐位移激励：
式中：x为激振振幅；ω为激励频率，则剪切应力为：
式中：Gs(ω)、Gl(ω)分别为储存剪切模量和损失剪切模量。通过动态剪切试验(试验装置和试验工装如图14所示)可得到不同频率下的动刚度Kd和滞后角φ，则存储刚度Ks和损失刚度Kl为：
存储力(Fs(ω))和损失力(Fl(ω))可由式(27)计算得到：
式中，x为激振振幅。则存储剪切应力τs(ω)和损失剪切应力τl(ω)为：
式中，4A为橡胶试片总剪切面积，剪切应变可由式(29)计算得到：
式中，δ为橡胶试片厚度，存储剪切模量Gs(ω)和损失剪切模量Gl(ω)为：
图14 动态实验装置及工装
Fig.14 Dynamic experimental device and tooling
可以使用参数ωRe(g*)和ωIm(g*)来表示材料的粘弹性，其中，参数ωRe(g*)和ωIm(g*)分别为归一化剪切松弛函数在频域的实部和虚部，可分别表示为：
式中，G∞为准静态剪切模量，可以直接用粘弹性参数ωRe(g*)和ωIm(g*)来表示材料的粘弹性(参数3)，参数ωRe(g*)和ωIm(g*)随频率的变化关系如图15和图16所示。
图15 ωRe(g*)随频率的变化关系(参数3)
Fig.15 Relation between ωRe(g*) and frequency(parameter 3)
图16 ωIm(g*)随频率的变化关系(参数3)
Fig.16 Relation between ωIm(g*) and frequency(parameter 3)
2.2.4 时域、频域粘弹性参数相互转换关系
频域的剪切存储模量(Gs(ω))和剪切损失模量(Gl(ω))可以写成如下形式：
式中：GR(0)为瞬时模量；Gs(ω)为存储模量；Gl(ω)为损失模量；ω为角频率；tr_i、
分别为各个Maxwell粘弹性单元对应的松弛时间和弹性模量比例系数。由式(34)~式(37)可实现频域和时域参数的相互转换。
3 计算结果与分析
3.1 计算结果
分别使用由松弛实验(参数1)、蠕变实验(参数 2)和简谐动态实验(参数3)获得的粘弹性参数，使用跨点动刚度有限元模型，对某款橡胶悬置进行动刚度分析，具体的橡胶悬置动刚度CAE分析结果与实验结果如图17所示。
图17 不同粘弹性参数进行CAE分析结果与测试结果对比
Fig.17 Comparison of CAE analysis results and test results with different viscoelastic parameters
由图17可知，使用由松弛实验(参数1)、蠕变实验(参数2)和简谐动态实验(参数3)获得的粘弹性参数对橡胶悬置进行有限元动刚度仿真得到的结果与试验值均较吻合，经计算可知相对误差均小于10%，满足工程应用要求。
综上所述，采用文中提出的由Mooney-Rivlin模型和多个Maxwell模型叠加组成的非线性粘弹性模型可以较好地表征橡胶材料的高频动态特性，同时验证了文中三种拟合Maxwell粘弹性参数方法的有效性。
3.2 三种粘弹性参数识别方法比较
通过蠕变实验拟合粘弹性参数时，需将蠕变数据转换为松弛数据(式(23))，二者的转化精度取决于橡胶材料中Boltzmann叠加原理的有效性；因此，对于粘弹性参数拟合方法，松弛实验拟合的粘弹性参数较蠕变实验的更精确。使用简谐动态实验进行粘弹性参数拟合时，由于需要在高频条件下测试，因此，对实验工装和橡胶试验试片的要求较高，其要求在参数拟合的频率范围内不能发生共振，否则会引入结构形状因素，导致识别结果不准确；对于频变特性比较明显的粘弹性材料，使用频域(简谐动态实验)方法进行粘弹性参数拟合较松弛实验和蠕变实验的更精确，因为，对于频变特性比较明显的橡胶，需要使用不同的prony级数组来描述不同频率段的粘弹性特性，即时域(松弛、蠕变)识别方法需要分段进行，因此，该情况下使用简谐动态实验进行粘弹性参数拟合更方便精确。
4 结论
(1) 文中首先搭建了橡胶隔振器的等效力学模型和有限元模型，分析了原点动刚度和跨点动刚度的区别，分析结果表明，对于橡胶隔振器，分析高频动态特性时，应使用跨点动特性模型。对于橡胶隔振器，由于其质量较小，所以当试验激振频率较低时，其产生的惯性力较小，此时原点动特性与跨点动特性的区别可忽略不计。但是，在分析其高频动态特性时，由于振动加速度较大，其惯性力的影响不能忽略，此时原点动特性与跨点动特性存在着较大的差别，不能再将二者等同；该结论为实际工程应用提供了参考。
(2) 测橡胶隔振器的高频动态特性时，为了消除附加惯性力对测试结果的影响，应使用跨点动态特性测试法。在原点动特性测试中，其力传感器布置在输入端；测试时，力传感器随作动头的运动而运动，同时输入端连接件、测试件工装等也随作动头的运动而运动，此时，输入力中既包含作动头的作动力，又包含传感器、连接件、测试件工装等附加质量产生的附加惯性力；因此，计算原点动特性时，需要对输入力进行修正，扣除附加惯性力成分，得到作动力，但由于附加惯性力与附加质量和测试频率有关，且对于不同的测试件，其工装不同，即附加惯性质量也会不同，此时附加惯性力的计算将非常复杂繁琐。因此对于橡胶隔振器的高频动特性试验，建议使用跨点动刚度测试法，将力传感器布置在输出端，消除附加惯性力对测试结果的影响。
(3) 文中采用一种由Mooney-Rivlin模型和多个Maxwell模型叠加组成的非线性粘弹性本构模型来描述橡胶的动态特性，分析结果验证了该模型的有效性。同时，文中给出了在时域范围(松弛实验、蠕变实验)和频域范围(简谐动态实验)拟合Maxwell粘弹性参数的方法，并分别使用由这三种方法拟合得到的粘弹性参数对某款橡胶悬置进行动刚度有限元分析，验证了三种粘弹性参数拟合方法的有效性。这些参数拟合方法为实际工程应用提供了参考，在实际工程中，可根据实际实验条件，选择一种参数拟合方法进行粘弹性参数识别。
(4) 文中简单对比了三种拟合粘弹性参数方法的优缺点，对于粘弹性的频变特性比较明显的材料，使用简谐动态实验进行粘弹性参数拟合较松弛实验和蠕变实验的更精确；由于橡胶材料中Boltzmann叠加原理的有效性的限制，松弛实验拟合的粘弹性较蠕变实验的更精确。使用简谐动态实验进行粘弹性参数拟合时，由于需要在高频条件下测试，因此，对实验工装和橡胶试验试片的要求较高，其要求在参数拟合的频率范围内不能发生共振，否则会引入结构形状因素，导致识别结果不准确。
(5) 文中设计的本构模型参数拟合实验均是基于简单的橡胶试片，即以简单的橡胶试片为材料级别研究对象，设计时域(松弛、蠕变)和频域(动刚度)三种参数识别实验，分别拟合得到橡胶材料的本构模型参数，并将拟合得到的本构模型参数用于具体橡胶隔振件的动刚度仿真中；仿真结果验证了以橡胶试片为材料级别研究对象的三种拟合粘弹性参数方法的有效性；即从橡胶试片拟合得到的本构模型参数可用于复杂零件级别的性能预测。因此，在工程实践中，可用试验成本较低、更容易获取的橡胶试片来进行材料级别试验，得到橡胶材料的本构模型参数，然后将其用于复杂零件的性能预测。
[1]何小静, 上官文斌.橡胶隔振器静态力-位移关系计算方法的研究[J].振动与冲击, 2012, 31(11): 91-97.He Xiaojing, Shangguan Wenbin.Calculating methods for force versus displacement relation of a rubber isolator[J].Journal of Vibration and Shock, 2012, 31(11): 91-97.(in Chinese)
[2]Kaya N, Erkek M Y, Caner Güven.Hyperelastic modelling and shape optimisation of vehicle rubber bushings [J].International Journal of Vehicle Design, 2016, 71(1/2/3/4): 212-225.
[3]魏志刚, 陈效华, 吴沈荣, 等.汽车控制臂衬套在整车工况下的有限元仿真[J].工程力学, 2016, 33(5): 234-240.Wei Zhigang, Chen Xiaohua, Wu Shenrong, et al.Finite element simulation of rubber bushing of automobile control arm in under full vehicle loading condition [J].Engineering Mechanics, 2016, 33(5): 234-240.(in Chinese)
[4]Shoyama T, Fujimoto K.Direct measurement of high-frequency viscoelastic properties of pre-deformed rubber [J].Polymer Testing, 2018, 67: 399-408.
[5]Lee Y H, Kim J S, Kim K J, et al.Prediction of dynamic stiffness on rubber components considering preloads [J].Materialwissenschaft und Werkstofftechnik, 2013, 44(5):372-379.
[6]Park S W.Analytical modeling of viscoelastic dampers for structural and vibration control [J].International Journal of Solids and Structures, 2001, 38(44/45): 8065-8092.
[7]Lee H S, Shin J K, Msolli S, et al.Prediction of the dynamic equivalent stiffness for a rubber bushing using the finite element method and empirical modeling [J].International Journal of Mechanics and Materials in Design, 2017, 15(1): 77-91.
[8]Khajehsaeid H, Baghani M, Naghdabadi R.Finite strain numerical analysis of elastomeric bushings under multi-axial loadings: A compressible visco-hyperelastic approach [J].International Journal of Mechanics & Materials in Design, 2013, 9(4): 385-399.
[9]李占龙, 孙大刚, 燕碧娟, 等.考虑形状参数的分数黏弹性振子频响特性[J].工程力学, 2017, 34(2): 232-240.Li Zhanlong, Sun Dagang, Yan Bijuan, et al.Frequency response of fractional viscoelastic oscillator considering geometric factor [J].Engineering Mechanics, 2017,34(2): 232-240.(in Chinese)
[10]Rickard Ö, Wentzel H, Kari L, et al.Constitutive modelling of the amplitude and frequency dependency of filled elastomers utilizing a modified boundary surface model [J].International Journal of Solids and Structures,2014, 51(19/20): 3431-3438.
[11]Wollscheid D, Lion A.The benefit of fractional derivatives in modelling the dynamics of filler-reinforced rubber under large strains: a comparison with the Maxwell- element approach [J].Computational Mechanics, 2014, 53(5): 1015-1031.
[12]Wollscheid D, Lion A.Predeformation- and frequencydependent material behaviour of filler-reinforced rubber:Experiments, constitutive modelling and parameter identification [J].International Journal of Solids and Structures, 2013, 50(9): 1217-1225.
[13]Pritz T.Five-parameter fractional derivative model for polymeric damping materials [J].Journal of Sound and Vibration, 2003, 265(5): 935-952.
[14]Shi H, Wu P.A nonlinear rubber spring model containing fractional derivatives for use in railroad vehicle dynamic analysis [J].Proceedings of the Institution of Mechanical Engineers, Part F: Journal of Rail and Rapid Transit,2016, 230(7): 1745-1759.
[15]Wang B, Kari L.A nonlinear constitutive model by spring,fractional derivative and modified bounding surface model to represent the amplitude, frequency and the magnetic dependency for Magneto-sensitive rubber [J].Journal of Sound and Vibration, 2019, 438: 344-352.
[16]吴杰, 上官文斌.采用粘弹性分数导数模型的橡胶隔振器动态特性的建模及应用[J].工程力学, 2008, 25(1):162-166.Wu Jie, Shangguan Wenbin.Modeling and applications of dynamic characteristics for rubber isolators using viscoelastic fractional derivative model [J].Engineering Mechanics, 2008, 25(1): 162-166.(in Chinese)
[17]胡小玲, 刘秀, 李明, 等.炭黑填充橡胶超弹性本构模的选取策略[J].工程力学, 2014, 31(5): 34-42.Hu Xiaoling, Liu Xiu, Li Ming, et al.Selectionstrategies of hyperelastic constitutive models for carbon black filled rubber [J].Engineering Mechanics, 2014, 31(5): 34-42.(in Chinese)
[18]Devendra M O, Abhijit L.CAE simulation approach to predict behavior of hyper-elastic (rubber) material [J].SAE Technical Paper, 2016, 9(1): 336-345.