杨海波，刘 枫，李凡珠，卢咏来，张立群

（北京化工大学 有机无机复合材料国家重点实验室，北京 100029）

采用有限元方法研究橡胶材料/制品的力学响应行为已成为一种方便有效的方法[1-5]。有限元分析中，精确的计算结果与可靠的材料模型和准确的材料参数密切相关，即本构方程的选取非常重要。在有限元分析中[6-10]较少提及本构方程选取的过程、依据以及准确性。因此系统地研究和总结不同应力-应变试验数据及不同本构方程对橡胶材料/制品有限元分析结果的影响很有必要。

本工作以炭黑填充天然橡胶（NR）制作的圆柱形橡胶试样的压缩变形过程为分析对象，分别采用不同的本构方程对三类应力-应变试验曲线进行拟合，并用有限元方法预测其压缩变形。通过对压缩变形试验数据的对比，判断不同应力-应变曲线及本构方程的计算精度。

1 实验

1.1 试样制备

胶料配方如下：NR 100，炭黑N234 40，氧化锌 3，硬脂酸 1，硫黄 1，促进剂CBS 1.5。使用密炼机将各组分按一定加料顺序在NR中混合和分散均匀，然后将混炼胶停放一定时间，使用圆盘式硫化仪测定硫化曲线，得到t90，再经平板硫化机硫化和气动裁压机冲压等步骤分别得到单轴拉伸（UT）试样、平面拉伸（PT）试样、等双轴拉伸（ET）试样[11]以及圆柱形单轴压缩试样。其中UT试样尺寸为50 mm×5.6 mm×2 mm，PT试样尺寸为60 mm×10 mm×2 mm，圆柱形单轴压缩试样[12]直径为17.8 mm，高为25 mm，如图1所示。

图1 各类试样示意

1.2 本构（应力-应变）曲线测试

使用电子万能试验机对橡胶试样做多次拉伸回复试验，以消除橡胶材料的应力软化效应（Mullins效应），得到稳定的力-位移数据。再将力-位移数据转为应力-应变数据。在进行有限元分析时，选取以下3类应力-应变试验曲线来拟合橡胶材料的超弹性本构模型：（1）UT应力-应变曲线；（2）单轴压缩（UC）应力-应变曲线；（3）UT，PT，ET组合的应力-应变曲线。对于PT试验，由于试样垂直于拉伸方向上的横截面积比单轴拉伸大得多，导致试样与夹具之间易产生滑脱，尤其在大变形加载时。为解决这一问题，采用平面拉伸试验结合光学测试技术，在试样表面涂上银粉，即可轻松地追踪试样的真实变形情况，从而修正由于试样的滑脱而导致力学数据不准确的问题。该方法类似于激光测距仪的功能。ET试验数据测定使用易瑞博科技（北京）有限公司的精密测试设备[11]。

1.3 试验验证数据测试

在有限元分析中，计算的是圆柱形橡胶试样单轴下压4.45 mm过程中对应的反作用力。为验证应力-应变曲线以及本构方程的计算精度，使用拉力机得到圆柱形橡胶试样的力-位移曲线。为保证试验数据的可靠性，单轴压缩试验进行4次，测试结果曲线高度重合。

2 有限元模型

2.1 几何模型

以炭黑填充NR材料制作的圆柱形单轴压缩试样为建模对象，采用轴对称模型建模，如图2所示，其中虚线为对称轴。

图2 圆柱形试样轴对称几何模型

2.2 材料模型

在材料模型的选取上，对3类应力-应变数据采用不同超弹性本构方程拟合，超弹性本构模型包括：一阶和二阶多项式模型（P_N1，P_N2）、一阶至六阶减缩多项式模型（RP_N1，RP_N2，RP_N3，RP_N4，RP_N5，RP_N6）、一阶至六阶Ogden模型（O_N1，O_N2，O_N3，O_N4，O_N5，O_N6）以及Arruda Boyce（AB），Marlow（Mar），Van Der Waals（VdW）模型共6类17种[7，13]，共有51个组合计算结果。

2.3 边界条件

在边界条件的施加上，有限元模型中载荷的施加与圆柱形橡胶试样的压缩变形试验完全相同，试样的上下表面建立两个刚性面来模拟试验设备的上下两个夹持面，载荷施加在两个刚性面上。刚性面与试样间建立接触。其中下端刚性面固定，上端刚性面受到Y轴负方向的4.45 mm的位移加载。整个模型共包含1 150个（CAX8RH）单元、1 224个节点。

3 结果与讨论

3.1 基于UT应力-应变数据的超弹性本构方程

3.1.1 方程拟合

不同本构方程对UT应力-应变试验数据的拟合结果见图3。由图3可见：P_N2，Mar以及二阶以上的O_Ni方程对于UT试验数据的拟合精度较高；P_N1，VdW，AB，O_N1以及三阶以下的RP_Ni方程的拟合数据小应变下拟合值低于试验值，中等应变下拟合值高于试验值，而大应变下拟合值远低于试验值；三阶及以上的RP_Ni方程拟合效果较好，但拟合曲线不平滑，存在波动，稳定性不佳。

图3 不同本构方程对UT应力-应变试验数据的拟合结果

3.1.2 有限元分析结果

UT试验数据拟合的超弹性本构方程对应的圆柱形试样的变形计算结果见图4。由图4可见：在使用UT应力-应变试验数据拟合本构方程时，拟合精度极高的O_Ni方程（二阶及以上）和P_N2方程的预测结果与试验数据差别较大，尤其在大位移变形情况下；而拟合精度极高的Mar方程在整个变形范围内的计算值与试验值接近；P_N1，VdW和O_N1方程的预测结果类似，在较小位移下（3 mm以内）的计算值与试验值接近，在较大位移下的计算值高于试验值；RP_N1和AB方程的预测结果类似，在小位移下的计算值低于试验值，而当位移增大时（5～7 mm），计算值与试验值接近；RP_N2方程在小位移下（3 mm以内）的预测值略低于试验值，随着加载位移的增大，其计算值与试验值吻合；三阶及以上的RP_Ni方程的整体预测结果与试验数据接近，但随着方程阶数的提高，材料模型的稳定性下降，计算出现收敛问题。

图4 UT试验数据拟合的超弹性本构方程对应的圆柱形试样的变形计算结果

3.2 基于UC应力-应变数据的超弹性本构方程

3.2.1 方程拟合

不同本构方程对UC应力-应变试验数据的拟合结果见图5。由图5可见：P_N1，P_N2，VdW，Mar和二阶及以上的O_Ni方程对UC试验数据的拟合精度较高；AB，O_N1和RP_N1方程在较小应变下（压缩应变小于0.3）拟合值低于试验值，较大应变下拟合值高于试验值；二阶及以上的RP_Ni方程拟合效果较好，但拟合曲线同样出现不平滑的现象，拟合曲线存在波动，稳定性不佳。

图5 不同本构方程对UC应力-应变试验数据的拟合结果

3.2.2 有限元分析结果

UC试验数据拟合的超弹性本构方程对应的圆柱形试样的变形计算结果见图6。由图6可见：使用UC应力-应变数据拟合本构方程时，拟合精度较高的P_N1，P_N2，VdW，Mar和二阶及以上的O_Ni方程的预测结果与试验数据接近，尤其是Mar方程；AB，O_N1和RP_N1方程在加载位移范围内计算值均低于试验值；二阶及以上的RP_Ni方程预测精度较高，尤其是RP_N5方程；RP_N2，RP_N4和RP_N6方程由于材料模型的不稳定性导致计算出现收敛问题。由此可见，相对于UT试验数据，采用UC试验数据拟合的超弹性本构方程预测圆柱形橡胶试样的单轴压缩变形行为计算精度大幅度提高。

图6 UC试验数据拟合的超弹性本构方程对应的圆柱形试样的变形计算结果

3.3 基于3种应力-应变数据的超弹性本构方程

3.3.1 方程拟合

不同本构方程对UT，PT，ET组合的应力-应变试验数据的拟合结果见图7—9。由图7可知：P_N1和AB方程与RP_N1方程的拟合结果类似；P_N2方程对ET试验数据的拟合效果较佳，但却不能很好地描述UT和PT数据；VdW方程在应变小于0.7时的整体拟合精度较高；Mar方程在应变小于0.4时可以很好地描述3类应力-应变曲线，当应变大于0.4时，ET的拟合数据开始高于相应的试验数据，但对于UT和PT数据，拟合数据与试验数据完全重合。由图8可知，RP_N1和RP_N2方程的整体拟合精度不理想，三阶以上的RP_Ni方程在较大应变时出现不稳定问题，RP_N3方程的拟合精度相对较高，值得一提的是不同阶数的RP_Ni方程在较小应变下（应变小于0.2）的拟合精度均比较理想。由图9可见，三阶及以下的O_Ni方程对UT，PT及ET组合的应力-应变试验数据的拟合精度较差，三阶以上的O_Ni方程可以很好地描述UT和ET试验数据，但当应变大于0.6时，PT的拟合数据与UT拟合数据重合，总体而言，O_N5和O_N6方程的拟合精度较高。

图7 多项式和AB，Mar，VdW模型本构方程对UT，PT及ET组合的应力-应变试验数据的拟合结果

图8 减缩多项式模型本构方程对UT，PT及ET组合的应力-应变试验数据的拟合结果

图9 Ogden模型本构方程对UT，PT及ET组合的应力-应变试验数据的拟合结果

3.3.2 有限元分析结果

采用UT，PT及ET组合的应力-应变数据拟合本构方程进行有限元分析时，整体计算结果的精度进一步提高。UC试验数据拟合的超弹性本构方程对应的圆柱形试样的变形计算结果见图10。图10（a）表明，AB，P_N1和P_N2方程的计算结果低于试验数据，VdW方程除了在小位移（小于1 mm）之外计算结果与试验数据相吻合，Mar方程在加载位移范围内的计算精度均较高。由图10（b）可知，RP_N1和RP_N2方程的计算结果低于试验数据，三阶及以上的RP_Ni方程计算精度较高，但仍存在小位移（小于1 mm）下计算精度较低的问题。由图10（c）可知，O_N1，O_N2和O_N3方程的计算结果低于试验数据，而O_N5和O_N6方程的计算精度高于较低阶数的O_Ni方程，且其计算结果与VdW方程类似，即小位移（小于1 mm）下的计算精度较低。

图10 UT，PT，EC组合的试验数据拟合的超弹性本构方程对应的圆柱形试样的变形计算结果

3.4 不同试验数据及超弹性本构方程的有限元分析结果对比

通过上述3类试验数据拟合橡胶材料的超弹性本构方程计算圆柱形橡胶试样的压缩变形结果可知，仅提供UT应力-应变数据拟合本构方程时，不能根据拟合精度判断各个方程的优劣，如拟合精度极高的O_Ni方程（二阶及以上）和P_N2方程的计算结果与试验数据差别极大。同样，高拟合精度的Mar方程在整个变形范围内的计算值与试验值接近，而拟合精度较差的RP_N2方程对应的计算结果却与试验数据接近。通过UC应力-应变数据拟合本构方程计算圆柱形橡胶试样的压缩变形行为时，其整体计算精度远高于通过UT数据拟合本构方程，而且拟合精度较高的本构方程对应的计算结果的精度亦较高。

相对于仅采用UT试验数据或UC试验数据，同时使用UT，PT以及ET试验数据拟合本构方程预测圆柱形橡胶试样的单轴压缩变形的整体计算精度进一步大大提高，如Mar，VdW以及高阶的O_Ni和RP_Ni方程的计算结果与试验数据吻合度极高。这主要是因为UT对应纯拉伸状态，PT对应纯剪切状态，ET对应纯压缩状态。本工作中，虽然圆柱形橡胶试样处于单轴压缩的加载工况，但其并非处在纯压缩状态，如图11所示。由图11可见，除了纵向压缩变形外[图11（b）]，该试样亦存在横向拉伸[图11（a）]和周向拉伸变形[图11（c）]，而图11（d）表明由于橡胶试样上下表面与上下加载面之间存在约束，区域A处的剪切变形较为明显。

图11 圆柱形橡胶试样的名义应变各分量

4 结论

对于应力-应变试验数据，仅采用单轴拉伸试验数据拟合本构方程时，不能单纯依赖拟合精度判断本构方程的优劣；而采用单轴压缩试验数据拟合本构方程时，其整体计算精度大大提升；同时采用单轴拉伸、平面拉伸和等双轴拉伸试验数据拟合本构方程时，计算精度进一步提升，且可依据拟合精度判断相应超弹性本构方程的优劣。

若仅提供一种变形状态下的应力-应变数据拟合本构方程，不能只根据拟合精度判断各方程的优劣，如采用单轴压缩试验数据拟合本构方程预测橡胶试样单轴压缩力学行为的整体表现均较好；若提供单轴拉伸、平面拉伸和等双轴拉伸试验数据拟合本构方程时，可以选用拟合精度高的O_Ni，P_Ni，VdW，Mar等方程作为橡胶材料的超弹性本构方程；而在仅有单轴拉伸数据的情况下，选用RP_Ni，AB，Mar方程较好，不要选用拟合精度高的O_Ni和P_N2方程。应变加载条件不同，同种本构方程在不同应变区间内的计算精度亦有所差别。

要准确预测处于复杂变形条件下的橡胶材料及制品的变形行为，需提供相应橡胶材料单轴拉伸、平面拉伸和等双轴拉伸试验数据，并选用拟合精度高的本构方程作为橡胶材料模型。