0 引言
橡胶应用广泛，在汽车、航空、船舶及机械行业中常作为密封、冲击吸能和抗震材料，具有重要的价值[1-2]。橡胶在周期运转下的弹性变形可以恢复，黏性变形不能恢复，会被吸收进而转变成热量，成为引起橡胶温度升高的滞后热源[3-4]。滞后生热引起的温度升高会影响橡胶的力学性能，严重时会造成橡胶损坏，因此研究橡胶材料的滞后生热具有重要的实际意义。
文献[5-7]采用解耦的方法将滚动轮胎的热力耦合分析分解为变形分析、损耗计算和热传导分析，为橡胶温度场的有限元分析建立了基本框架。初红艳等[8]在分析橡胶材料生热机理的基础上，对挤压接触对滚的一对钢-橡胶辊在转动过程中的生热进行研究，分析了随时间变化的胶料温度升高的规律。吴福麟等[9]根据解耦的分析思想进行了滚动轮胎稳态温度场的有限元分析，用实测的轮胎内部点的稳态温度场对分析结果进行了验证。但在实际生产中，并非所有的橡胶结构变形产生的应力应变关系均为传统三角函数关系。
本文以一对挤压对滚的钢-橡胶辊为研究对象，采用解耦方法进行瞬态力学场分析，再通过损耗计算得出节点生热率，最后将滞后生热作为内热源进行热分析，得到其稳态温度场分布。实验验证了仿真分析的有效性。
1 橡胶材料生热机理
橡胶辊的生热来源主要有两个：一是在转动过程中，胶料内部的交变变形产生的热量，二是运动过程中的摩擦生热。滞后生热是由胶料的交变变形引起的。橡胶材料具有黏弹性，在循环挤压旋转过程中，应力与应变的不同步导致机械功损失转化为热能，即滞后生热[10]。生热率是橡胶材料在单位体积、单位时间内产生的滞后生热。研究橡胶辊的生热特性，计算出橡胶辊生热率，在温度场分析过程中可以用有限元法将其生热率离散到各个节点上，从而计算得到整个橡胶辊的温度场。
橡胶辊在旋转过程中的变形引起了应力和应变，但橡胶的黏弹特性使应力与应变不同步，呈现相位滞后，形成损耗角。KABE等[11]用三角函数表述应力σ(t)与应变ε(t)的变化规律(用于研究橡胶材料滞后生热)：
(1)
式中，t为时间，s；ω为角速度，rad/s；ε0为应变振幅；σ0为应力振幅；δ为滞后损失相位差。
应力与应变关系如图1所示，滞后损失回路即滞后圈如图2所示。
图1 应力应变随时间变化关系
Fig.1 Stress and strain relationship over time
图2 滞后圈
Fig.2 Graph of hysteresis loss
材料出现滞后现象时，每一个周期均会有热能损失——力学损耗。对橡胶辊内任意点而言，在一个转动周期内，橡胶材料损耗的能量即所产生的热量为 [12]
(2)
式中，σ、ε分别为瞬时应力和应变。
由式(2)可知，要计算生热量Q，需知道应力和应变随时间的变化关系。把式(1)中的应力展开成σ=σ0sinωtcosδ+σ0cosωtsinδ，可见应力由两部分组成：与应变同相位的σ0cosδ和与应变相位相差90°的σ0sinδ。于是应力和应变的关系可用一个与应变同相位的量E′和一个与应变相差90°的量E″表示[13]：
σ=ε0(E′sinωt+E″cosωt)
(3)
E′=σ0cosδ/ε0 E″=σ0sinδ/ε0
式中，E′为储存弹性模量；E″为损失弹性模量。
把式(1)、式(3)代入式(2)，得到每个旋转周期橡胶辊料所产生的热量
(4)
其中，tanδ为损耗因子，是用来表征橡胶材料能量损失的一个重要参数。则橡胶辊的生热率为
(5)
式(5)是用橡胶辊应力应变三角函数关系算得的生热率公式，本文将针对一对挤压对滚的钢-橡胶辊运用解耦的方法，先进行力学场仿真，找到橡胶辊的应力应变的变化规律，并进行损耗计算，将计算出的生热率对橡胶辊进行稳态热分析，再进行试验验证。
2 挤压对滚的钢-橡胶辊的力学场仿真
运用仿真软件ANSYS对图3中的一对挤压对滚的钢-橡胶辊进行瞬态挤压对滚的力学场仿真分析，得出模型的变形情况及应力应变大小。根据实际情况将模型简化，取两辊的一段进行建模，建立两辊初始时刻相切的模型。模型上部为钢辊，直径为65 mm，长度为30 mm；下部为橡胶辊，中心轴为钢轴，外面包裹着厚度为10 mm的橡胶层。为缩短ANSYS仿真计算时间、提高计算的准确性，橡胶辊只建立橡胶层部分的模型，内部由一个有质量的节点代替钢轴，并在此节点上加载钢轴的质量和转动惯量。橡胶辊模型为内径45 mm、外径65 mm、长度30 mm的空心圆柱体。
图3 挤压对滚的钢-橡胶结构
Fig.3 Squeezing and rolling steel-rubber structure
为简便计算，钢辊建模选择SOLID281壳体单元，弹性模量为206 GPa，泊松比为0.3，密度为7 800 kg/m3。橡胶辊建模选择SOLID186实体单元，密度为1 300 kg/m3，橡胶邵氏硬度为50。力学场分析关注橡胶的变形情况及应力应变大小，故采用能体现大变形的超弹性M-R本构关系来代替黏弹性本构关系，ANSYS中输入的待定力学性能常数C10=0.436，C01=0.03。钢辊、橡胶辊的网格划方法分别为自由网格划分和扫掠网格划分。对两辊进行接触对的设置，并给钢辊施加0.5 MPa的下压力和21 rad/s的角速度，钢辊向下挤压时间为1 s，然后带动橡胶辊旋转1 s进行瞬态分析。
图4 胶辊应力云图
Fig.4 Stress nephogram of rubber roller
图4为邵氏硬度为50的橡胶辊旋转1 s之后的应力云图，可看出，由于钢辊与橡胶辊挤压旋转，橡胶辊在接触区域产生形变，在沿轴挤压力的作用下向外膨胀，胶辊轴向两端的应力最大；沿周向方向，越靠近挤压区域的应力越大，挤压区域应力最大，应变规律与应力一致。以胶辊中心截面径向方向为研究对象，沿径向由内向外依次提取5个节点的最大应力应变，如图5所示，应力应变由内向外先增大后减小，中间位置应力应变最大。
图5 沿径向由内向外节点应力应变曲线
Fig.5 Stress and strain of nodes along the radial direction from inside to out
3 损耗与生热率计算
损耗计算中，基于超弹性本构关系的变形分析结果，用黏弹性性质来计算损耗生热，采用节点等效法计算得到节点生热率。
提取胶辊上任一节点在旋转过程的应力应变，如图6所示，由于力学场仿真中橡胶材料采用超弹模型进行变形分析，故其应力应变同步变化；随着胶辊进出挤压区，应力先逐渐增大再逐渐减小至趋于零。实验模型中，橡胶应力应变的变化不符合正弦曲线，而是符合高斯曲线，即在挤压区中心为应力应变峰值，MATLAB拟合出的二次高斯曲线如图7所示。
图6 任一节点在旋转过程中应力应变图
Fig.6 Stress-strain graph of arbitrary point in rotating process
图7 应力拟合图
Fig.7 Fitting figure of stress
力学场仿真中，橡胶采用的是超弹性模型，故其应力应变同步变化，不存在黏弹滞后，但实际实验中的橡胶具有黏弹特性，故在损耗计算中，根据德国GAOBA EPLECOR 500N DMA动态热机械分析仪实验所得的丁腈橡胶(邵氏硬度为50)黏弹性参数，计算出滞后时间，并代入到应力应变关系式，求解出橡胶材料的生热率。计算过程如下：
根据动态热机械分析仪实验，在室温为20 ℃的条件下，硬度为50的胶辊所用橡胶材料的损耗因子tanδ=0.126，故损耗角δ=0.125 34 rad。
根据式(1)，可计算出滞后时间
t=δ/ω=0.125 34/21=0.005 969 s
(6)
运动周期
T=2π/ω=2π/21=0.299 2 s
(7)
将滞后时间加到应力应变随时间的变化数据中，使应变滞后于应力5.969 ms，可得出应力应变关系，如图8所示。
图8 应力应变滞后图
Fig.8 Lag figure of stress-strain
图9 滞后圈
Fig.9 The graph of hysteresis loss
根据滞后的应力应变数据可以得出该节点的滞后圈，如图9所示。节点的生热量即为滞后圈的面积，采用MATLAB自适应Simpson法进行数值积分求解，可以计算出节点生热率。以橡胶辊某节点为研究对象，提取这些节点的应力应变随时间变化的数据，并代入滞后时间，得出每个节点的滞后圈，按照上述损耗计算方法，可计算出这5个节点的高斯生热率。根据这5个节点在旋转过程中的最大应力应变，结合式(5)进行计算，可得转速为21 rad/s，sinδ=0.125 012条件下节点的三角函数生热率与高斯函数生热率，如表1所示。
表1 节点生热率
Tab.1 Heat generation rate of nodes
节点序号最大应力(kPa)最大应变三角函数生热率J/(m3·s)高斯函数生热率J/(m3·s)1170.3610.095 7390.021 4090.085 4132209.6100.120 0670.033 0350.185 4553288.7460.166 9240.063 2670.403 5684283.7570.162 4910.060 5230.375 6415159.5730.092 5580.019 3870.098 530
三角函数生热率只需将节点最大应力应变与固定常数相乘即可得出，而高斯生热率需要提取旋转过程中的所有应力应变，加入滞后时间再进行积分运算才能得出结果。由于橡胶辊与钢辊对滚为周期性旋转运动，故在理想状态下，不同节点在相同位置产生的应力应变应相等。橡胶辊共有4 444个任意分布的节点，若将所有节点一一进行数据提取和积分计算，则计算量太大且大部分是重复数据。对比表1中各节点的三角函数生热率和高斯生热率可看出，2种方法计算出来的生热率具有相同的变化规律，所以所有节点的高斯生热率不必一一进行拟合积分，只需将这5个节点的高斯生热率与三角函数生热率的数值进行对比拟合，得出三角函数生热率的修正公式，即提取出橡胶辊上全部节点的最大应力应变，计算出所有节点的三角函数生热率，再利用修正公式即可求出所有节点适于本实验模型的节点高斯生热率。这种计算方法既大大减小了计算量，又可保证数值的准确性，避免重复计算产生不必要的误差。运用MATLAB对这5个点的两种生热率进行拟合，得出的修正公式为
q=x2+6.532x-0.042 68
(8)
其中，x为三角函数生热率；q为对应的高斯函数生热率。5个节点修正之后的生热率以及修正误差如表2所示，修正误差均没有超过15%，故修正公式可用。将所有节点在转动过程中的最大应力及应变代入到式(5)，再利用式(8)进行修正，得出所有节点的高斯函数生热率，4 444个点的数据在此不一一列出。
表2 不同方法节点生热率
Tab.2 Heat generation rate of nodes by different methods
节点序号三角函数生热率J/(m3·s)高斯函数生热率J/(m3·s)修正生热率J/(m3·s)绝对修正误差(%)10.021 4090.085 4130.097 57914.24420.033 0350.185 4550.175 2215.51830.063 2670.403 5680.378 4596.22240.060 5230.375 6410.359 8614.20150.019 3870.098 5300.084 66914.068
4 橡胶辊温度场仿真
4.1 稳态温度场分析
由于本文主要研究橡胶的黏弹性损耗产生的滞后生热，故在仿真模型中只对橡胶辊进行温度场分析,其中，橡胶材料的参数如表3所示。
表3 材料参数
Tab.3 Material parameters
弹性模量(GPa)0.008密度(kg·m3)1 300热导率(W/(m·K))0.25质量热容(J/(kg·K))1 700热膨胀系数(K-1)65×10-5泊松比0.47
ANSYS稳态温度场仿真分析中，橡胶辊模型的单元为三维稳态热分析单元SOLID90，该单元等效于力学场分析中的橡胶结构单元SOLID186，则可划分出与力学场分析中的橡胶辊相同的网格，且可保证热分析中橡胶辊上所有节点的编号与力学场分析中的节点编号对应一致。热边界条件设置为：橡胶外表面与空气的传热系数设为280 W/(m2·K)，内表面与钢质辊芯的传热系数设为2 000 W/(m2·K)；环境温度设置为20 ℃。内热源为计算所得的节点生热率，利用ANSYS APDL BF生热率载荷命令将上一节中求出的所有节点的高斯生热率施加到对应的节点上，得到橡胶辊模型的节点生热率，如图10所示。热边界条件设置全部完成后，对橡胶辊进行稳态分析设置并求解。
图10 橡胶辊生热率图
Fig.10 Heat generation rate figure of rubber roller
图11 橡胶辊温度云图
Fig.11 Temperaturenephogram of rubber roller
4.2 温度场的仿真结果分析
图11为橡胶辊在达到稳定状态的温度场云图，橡胶辊整体温度最高的地方为胶辊轴向的两端、径向的中间位置，最高温度为39.882 6 ℃；沿橡胶辊中心截面处径向方向，内部温度最低，为26.003 1 ℃，温度由内向外先增大后减小；沿橡胶辊外表面轴向方向，两端温度最高，为36.798 2 ℃，中间温度最低，为35.256 1 ℃。沿径向由内向外依次提取橡胶辊中心截面的5个节点，其温度由内向外先增大后减小，最高温度出现在橡胶层中间位置，如图12所示，与对应节点的应力应变(图5)规律一致，说明橡胶材料滞后生热与应力应变的大小成正相关，且生热率越高，产生温度越高。
图12 径向温度数据
Fig.12 Temperature data along radial direction
将用三角函数计算出的全部节点的生热率按照4.1节的稳态温度场分析方法和条件进行仿真，得出的温度规律与图12所示的规律基本一致，但温度不同，最高温度仅24.530 5 ℃，最低温度为21.007 5 ℃。
5 挤压对滚的钢-橡胶辊的温升试验
5.1 试验设备及材料
主要试验设备为：钢-橡胶辊挤压对滚试验平台、红外温度枪以及LX-A邵氏硬度计。试验所用橡胶辊为钢质辊芯外包覆一层橡胶，橡胶层厚度为10 mm、硬度为50，橡胶层材料为加入炭黑补强剂的硫化丁腈橡胶，胶辊和钢辊直径均为65 mm。
5.2 试验方案
将试验平台的转速设定为200 r/min，两挤压辊间的压力为0.5 MPa。对橡胶辊进行试验，试验开始前，记录初始环境温度。开始试验后，每隔5 min进行一次温度测量：使用红外温度枪测量两辊中部辊间挤压区域的温度(图3)和挤压区两端温度，当挤压区域温度上升至趋于稳定时，终止试验。
5.3 试验数据分析
试验过程中，橡胶辊与钢辊挤压对滚，橡胶材料由于周期性运转而产生滞后损失，导致内部热量积累，引起温度升高。每5 min记录一个辊间温度，随着橡胶辊与钢辊的循环挤压运转，两辊挤压区两端的温度最高，中部的温度最低，且中部辊间温度先骤然升高，再持续缓慢升高直至最大值，中心辊间温度数据如表4所示。
表4 中心辊间温度数据
Tab.4 Temperature data of roller center
时间(min)051015202530温度(℃)2024.927.329.931.833.334.2时间(min)35404550556065温度(℃)35.536.837.538.239.540.441.1
5.4 仿真与试验结果的对比
温升试验中，硬度为50的橡胶辊中心辊间最高温度为41.1 ℃，升高了21.1 ℃；温度仿真结果中，相同初始温度下，橡胶辊外表面中间位置最高温度为35.256 1 ℃，升高了15.256 1 ℃，可见仿真的结果与实际试验结果相仿，验证了以高斯函数生热率作为热源进行生热分析的正确性。仿真中，橡胶的温度完全由橡胶滞后生热产生。试验中，钢辊与橡胶辊有相对滚动，故存在摩擦生热，但滚动状态下摩擦产生的生热量较少，且摩擦生热主要导致橡胶辊表面的温升，所以试验与仿真的温差部分来源于钢辊与橡胶辊对滚产生的摩擦生热。与橡胶滞后生热相比，摩擦生热影响较小。温升试验与温度场仿真中，中心辊间温度都没有两端温度高，这是由于胶辊两端受钢辊挤压变形量大、产生的应力应变大。
6 结论
(1)两辊旋转挤压过程中，橡胶层的黏弹特性引起形变滞后，损失的机械能转变成热能，导致旋转过程中的钢辊和橡胶辊的温度升高；胶辊的应力应变随时间的变化关系并非传统的三角函数关系，而是符合高斯函数的变化关系。
(2)由仿真分析可知，橡胶辊中心截面处的温度在径向由内向外先增大后减小，中间位置的值最大；在轴向，由于两端变形大，应力应变较大，所以温度最高，试验测得的数据也具有相同的规律。
(3)由仿真和试验分析可知，胶辊与钢辊挤压对滚产生的热量大部分源于橡胶材料滞后生热，摩擦生热的影响较小。
(4)采用高斯函数的仿真和试验分析得出了相同的温度分布规律，数值上也非常接近；采用传统三角函数进行仿真得出的温度数据与实际情况相差太大，证明了采用高斯函数进行节点生热率的计算是符合该橡胶结构的，也证明了仿真模型和解耦方法的有效性。
[1] 王宝珍,胡时胜,周相荣.不同温度下橡胶的动态力学性能及本构模型研究[J].实验力学,2007,22(1)：1-6.
WANG Baozhen,HU Shisheng,ZHOU Xiangrong.Research of Dynamic Mechanical Behavior and Constitutive Model of Rubber under Different Temperatures[J].Journal of Experimental Mechanics,2007,22(1)：1-6.
[2] SARKAWI S S,DIERKES W K,NOORDERMEER J W M.The Influence of Non-rubber Constituents on Performance of Silica Reinforced Natural Rubber Compounds[J].European Polymer Journal,2013,49 (10)：3199-3209.
[3] 何燕,马连湘,黄素逸.轮胎橡胶材料生热率的测定及分析[J].橡胶工业,2004,51 (1)：53-55.
HE Yan,MA Lianxiang,HUANG Suyi.Measurement and Analysis of Heating Rate of Tire Rubber[J].Rubber Industry,2004,51 (1)：53-55.
[4] 何燕.轮胎非稳态温度场的研究[D].武汉:华中科技大学,2003.
HE Yan．Research on Unsteady Temperature Field of Tire[D].Wuhan：Huazhong University of Science and Technology,2003.
[5] SEGALMAN D J.Modeling Tire Energy Dissipation for Power Loss Calculations[J].SAE Technical Paper,1981：810162.
[6] BROWNE A L,ARAMBAGES A.Modeling the Thermal State of Tires for Powers Loss Calculations [J].SAE Technical Paper,1981：810163.
[7] EBBOTT T G,HOHMAN R L,JEUSETTE J P,et al.Tire Temperature and Rolling Resistance Prediction with Finite Element Analysis[J].Tire Science and Technology,1999,27(1)：2-21.
[8] 初红艳,许康健,黄为,等.橡胶材料生热及其对力学性能的影响[J].北京工业大学学报,2017(11):32-37.
CHU Hongyan,XU Kangjian,HUANG Wei,et al.Heat-generating of Rubber and Its Influence on Mechanical Properties[J].Journal of Beijing University of Technology,2017(11):32-37.
[9] 吴福麟，李子然，夏源明.不同载荷和初始气压下滚动轮胎稳态温度场的测试与有限元分析[J].工程力学,2008,25(1)：54-60.
WU Fulin，LI Ziran，XIA Yuanming.Experiment and Finite Element Analysis of Steady-state Temperature Distribution in Rolling Tire under Different Loads and Initial Inflation Pressures[J].Engineering Mechanics,2008,25(1)：54-60.
[10] 曾玉铧.橡胶压缩屈挠温升与生热和导热关系的研究[D].广州:华南理工大学,2015.
ZENG Yuhua.Relationship between Temperature Rise with Heat Generation and Thermal Conduction of the Rubber in the Compression Flexometer Test[D].Guangzhou：South China University of Technology,2015.
[11] KABE K,SHIMOMURA I.Application of Tyre Stress Analysis and Thermal Analysis to Tyre Durability[J].International Polymer Science and Technology,1984,11(8)：32-39.
[12] 张强,蒋豹,刘昱良.定子橡胶衬套生热影响因素分析[J].化工机械,2017,44 (3)：314-318.
ZHANG Qiang,JIANG Bao,LIU Yuliang.Analysis of Factors Affecting Heat Generation of Stator Rubber Lining[J].Chemical Engineering &Machinery,2017,44 (3)：314-318.
[13] 宋喜岗,崔淑华,郑雪莲.汽车轮胎稳态温升机理的研究[J].机电产品开发与创新,2011,24 (3)：44-46.
SONG Xigang,CUI Shuhua,ZHENG Xuelian.The Research on Vehicle Tires Thermal Mechanism[J].Development &Innovation of Machinery &Electrical Products,2011,24 (3)：44-46.
