减振橡胶疲劳黏滞生热的仿真分析-源文件与子程序详解
文章来源：仿真侠
00
—
论文链接
论文链接：
减振橡胶疲劳黏滞生热的仿真分析[J]
http://dx.chinadoi.cn/10.13465/j.cnki.jvs.2021.12.026
http://dx.chinadoi.cn/10.13465/j.cnki.jvs.2021.12.026
减振橡胶疲劳黏滞生热的仿真分析[D]
http://dx.chinadoi.cn/10.27426/d.cnki.gxtdu.2020.000389
http://dx.chinadoi.cn/10.27426/d.cnki.gxtdu.2020.000389
01
—
研究背景
02
—
03
—
炭黑填充橡胶的超弹性本构模型
04
—
减振橡胶疲劳黏滞生热的试验与理论分析
4.4 黏滞耗能率与频率幅值温度关系
05
—
仿真分析
5.1 有限元分析流程
5.2 子程序设计
5.3圆柱试样黏滞生热仿真分析
5.4沙漏试样黏滞生热仿真分析
06
—
6.2基于自热温升的橡胶疲劳寿命仿真流程
07
—
总结
采用单轴拉伸(ST)、平面拉伸(PT)和等双轴拉伸(ET)三种橡胶测试试验，拟合相关试验参数，得到了橡胶材料的Ogden本构模型及相关参数。
采用一种修正的Kraus模型定量描述了橡胶材料动态损耗模量随温度、载荷频率和应变幅值的变化规律。得到了生热率与温度、载荷频率和应变幅值的函数关系式。
利用依黏弹性理论得出的黏滞生热率与温度、载荷频率和应变幅值的函数关系式，编制了相应的计算程序。建立了减振橡胶疲劳黏滞生热的有限元分析方法。
通过将经典疲劳模型中用作疲劳寿命预测指标的最大主应变替换为稳态温升，在幂律模型的基础上开发了一种方法来快速评估橡胶结构的疲劳寿命。
08
—
源文件与操作步骤（沙漏试样为例）
8.1分析流程
仿真分析主要包括三个环节：变形分析、热源计算与热分析。（1）在变形分析环节，对材料和减振元件施加设定的载荷历史，采用超弹性本构描述橡胶材料的力学行为，求解每个加载时刻有限元模型中各积分点的应变状态；（2）在热源计算环节，对应每一加载时刻，将变形分析中对应的载荷频率、应变状态（动态应变幅值）以及热分析中得到的温度作为输入变量，通过自编的Fortran语言子程序，计算得到各积分点的黏滞生热率；（3）依已知的材料参数和问题的热边界条件进行Abaqus热分析，得出温度分布后再将温度场数据返回到自编子程序，对黏滞生热强度和温度场进行迭代计算，从而得出橡胶材料和减振元件各位置的温升历程。
8.2建模
abaqus/cae操作
8.3赋材料属性
钢：
CAE操作：
inp文件：
*Material, name=ste*Conductivity43.,*Density 7.83e-09,*Elastic210000., 0.3*Expansion 8e-06,*Specific Heat 4.4e+08,
橡胶：
CAE操作：
inp文件：
*Material, name=rub*Conductivity 0.2,*Density 1.2e-09,*Expansion 8e-05,**为hetval子程序定义生热率做准备*Heat Generation*Hyperelastic, n=3, ogden 1.9042, 1.0625, -1.924e-10, -17.7, 0.0003185, 12.3795, 1e-05, 1e-05 1e-05,*Specific Heat 1.7e+09,**以下两个参数均为usdfld子程序所需，定义预定义场与三个过程变量。*User Defined Field*Depvar 3,
8.4网格划分
CAE操作：
8.5定义分析步
静力分析：
CAE操作：
瞬态热力耦合分析：
CAE操作：
8.6定义相互作用
CAE操作：
8.7定义载荷边界
试验加载方式为应力比R=0.1频率5Hz最大载荷400N正弦循环加载，即4N-400N-4N循环。故第一个静力分析步加4N，通过usfld子程序提取此时的应变场；第二个分析步加400N，同样通过usfld子程序提取此时的应变场，（应变场-4N➖应变场-400N）➗2 即为应变幅值，传入hetval子程序定义生热率场在第三个瞬态热力耦合分析步生效。
CAE操作：
inp文件：
*Step, name=Step-1, nlgeom=YES*Static0.1, 1., 1e-05, 1.*BoundarySet-3, 1, 1Set-3, 2, 2Set-3, 6, 6*CloadSet-2, 2, 40.*End Step*Step, name=Step-1_2, nlgeom=YES*Static0.1, 1., 1e-05, 1.*CloadSet-2, 2, 400.*End Step*Step, name=Step-2, nlgeom=YES*Coupled Temperature-displacement, creep=none, deltmx=10.1., 4000., 0.04, 4000.*SfilmSurf-1, F, 20., 0.025*End Step
最后需要在inp文件model level添加*Initial Conditions关键字，初始化场变量，（usdfld子程序所需，/CAE中使用edit key words操作。）
*Initial Conditions, type=FIELD, VARIABLE=1*step
8.8提交计算
批处理提交
call abaqus job=*** user=***.f
09
—
子程序
C =====================================================C This ABAQUS user subroutine was produced byC Xia JiangC (E-mail: jiangxia127@hotmail.com)CC Xiangtan UniversityC =====================================================CCCSUBROUTINE HETVAL(CMNAME,TEMP,TIME,DTIME,STATEV,FLUX,1 PREDEF,DPRED)CINCLUDE 'ABA_PARAM.INC'CC Set frequencyCPARAMETER(F=10.,PI=3.1415926)CCHARACTER*80 CMNAMECREAL*8 LM,ratio_LMCDIMENSION TEMP(2),STATEV(3),PREDEF(*),TIME(2),FLUX(2),1 DPRED(*)CC value of max strain,pass in from usdfld.CEMAX = STATEV(1)CC value of current temperature,pass in from usdfld.CT = STATEV(2)CC value of min value of strain,pass in from usdfld.CEMIN = STATEV(3)CC increment of strain,pass in from usdfld.CDE = abs(EMAX-EMIN)/2CC step timeCTi = time(1)CC value of Loss ModulusCLM = exp(-0.014*(T-23.))*(1.967+0.013*F)*(DE/0.01)**0.434/(1.+(DE/0.01)**0.868)C LM = exp(-0.014*(T-23.))*2*0.9206*(DE/1.24)**0.1845/(1+(DE/1.24)**(2*0.7845))CC Modification based on dynamic property softening effectCc ratio_LM = 0.1619*exp(-Ti/163.7594)+0.7997CC Heat flux calculationCC FLUX(1) = F*PI*DE**2.*LM*ratio_LMFLUX(1) = F*PI*DE**2.*LMCC Write Heat flux,Time,Temprature,StrainMAX,StrainMIN,DStrain,Loss Modulus to msg fileCc write(7,*)'--------------------------------------------------------------------'c write(7,*)'Time ', time(2),' Temprature',Tc write(7,*)'StrainMAX',EMAX,'StrainMIN ',STATEV(3)c write(7,*)'HeatFlux ',FLUX(1),'DStrain ',DEc write(7,*)'LossModulus',LMc write(7,*)'--------------------------------------------------------------------'CCRETURNENDCCCSUBROUTINE usdfld(field,statev,pnewdt,direct,t,celent,time,dtime,1 cmname,orname,nfield,nstatv,noel,npt,layer,kspt,kstep,kinc,2 ndi,nshr,coord,jmac,jmtyp,matlayo,laccflg)Cinclude 'aba_param.inc'Ccharacter*80 cmname,ornameCcharacter*3 flgray(15)Cdimension field(nfield),statev(nstatv),direct(3,3),t(3,3),time(2),1 coord(*),jmac(*),jmtyp(*)Cdimension array(15),jarray(15)CC value of current temperature:Ccall getvrm('TEMP',array,jarray,flgray,jrcd,1 jmac, jmtyp, matlayo, laccflg)CTEM = ARRAY(1)CC Store temperature as a solution dependent state variableCSTATEV(2) = TEMCC Value of current strain :Ccall getvrm('NEP',array,jarray,flgray,jrcd,1 jmac, jmtyp, matlayo, laccflg)CDSTRAIN = abs(ARRAY(1))CC Maximum value of strain up to this point in time:Ccall getvrm('SDV',array,jarray,flgray,jrcd,1 jmac, jmtyp, matlayo, laccflg)CEPSMAX = ARRAY(1)CC Store maximum value of strain as a solution dependent state variableCSTATEV(1) = MAX(DSTRAIN,EPSMAX)CC MINrage value of strainCif(time(2).eq.1) thenCC Store MINrage value of strain as a solution dependent state variableCSTATEV(3) = STATEV(1)CC Write StrainMIN to msg fileCC write(7,*)'MINrage value of strain',STATEV(3)Cend ifCRETURNEND
