未硫化橡胶黏弹塑性本构模型及有限元实现

图1 未硫化橡胶试件

Fig.1 Specimen of uncured rubber

1.2 橡胶实验及结果分析

实验均是在Instron E3000试验机上开展,包括单调拉伸和循环拉伸加卸载两种变形模式,通过单调拉伸实验,可以观察到未硫化橡胶的弹性水平和屈服特性,而循环拉伸加卸载实验则可综合反映出未硫化橡胶的超弹性、黏性、塑性等力学行为特征.作者对轮胎生产中压延、成型、定型等工艺过程进行了初步仿真(此部分工作还未整理发表),结果表明上述工艺过程中未硫化橡胶应变率最大不超过 0.1 s^-1,故本文在每种变形模式下均开展4种不同应变率( $\overset{\cdot}{ε}$ =0.001,0.003,0.03,0.1 s^-1)下的拉伸力学实验,以观察其应变率相关性.每次实验开展3次重复性实验以保证实验结果的可靠性.

单调拉伸工况下未硫化橡胶的拉伸实验结果如图2(a)所示.图中:ε_eng为工程应变;σ_eng为工程应力.从图中可以看出,不同应变率下未硫化橡胶具有相似的非线性力学特征,在拉伸初始时具有很小的线性段,随后发生软化行为,最终应力趋于平缓.胶料拉伸应力水平表现出明显的应变率相关性,随着应变率增加,初始刚度显著增加,应力水平提高.相比于硫化橡胶,未硫化橡胶在100%的应变范围内未出现反“S”现象,可能由其内部未形成交联网络所致.

图2

图2 未硫化橡胶不同应变率下拉伸实验结果

Fig.2 Mechanical behaviors of uncured rubber at different strain rates

循环拉伸加卸载实验中,未硫化橡胶同样表现出明显的应变率相关性,且在加载段的应变率相关性略大于卸载段的应变率相关性,如图2(b)所示.此外,还可以看出,未硫化橡胶存在较为明显的Mullins软化现象,以及显著的迟滞效应和残余变形.Mullins软化效应反映了试件在首次拉伸时的损伤,为了更清晰地观察未硫化橡胶的软化现象,图3给出了其在不同应变率下随拉伸循环次数(C)增加,拉伸上升段在前20%应变范围内的切变模量(G).可见,随着循环次数增加,胶料切变模量明显降低,且在较高应变率时降低速度更快.迟滞效应反映了循环加卸载实验中的能量耗散,未硫化橡胶迟滞效应同样受到应变率的影响,随着应变率增加,迟滞效应有所增强.由于未硫化橡胶黏塑性较强,每次卸载后均存在一定的塑性残余变形,且残余变形量随着循环次数的增加而累计增加,随着应变率的增加有所降低.

图3

图3 不同应变率下4次上升段切变模量

Fig.3 Shear modulus in four loading steps at different strain rates

2 未硫化橡胶黏弹塑性本构模型

基于对未硫化橡胶实验结果的分析,提出了TN模型对未硫化橡胶的力学行为进行表征,模型的黏塑性响应基于对变形梯度的乘法分解,具体推导过程介绍如下.

2.1 变形梯度分解

对于不可压缩材料,变形梯度F通过乘法分解为形变分量(Distortional Part)与体积分量(Dilatational Part)^[17]:

$\begin{matrix} F = F_{d i s} F_{d i l} \end{matrix}$

(1)

形变分量与体积分量分别具有以下形式^[18]:

$\begin{matrix} F_{d i s} = J^{- \frac{1}{3}} F, F_{d i l} = J^{\frac{1}{3}} I \end{matrix}$

(2)

式中: $J$ 为变形梯度的行列式; $I$ 为单位张量;det(F_dis)=1.

此时左Cauchy-Green变形张量可表示为

$\begin{matrix} \begin{matrix} B = F F^{T} = J^{2 / 3} B^{*} \end{matrix} \end{matrix}$

(3)

$\begin{matrix} B^{*} = F_{d i s} F_{d i s}^{T} = J^{- 2 / 3} B \end{matrix}$

(4)

式中:B^*为B的形变分量.

对于近似不可压缩橡胶类材料黏塑性变形,形变分量也可分解为弹性变形分量和黏塑性分量^[19]:

$\begin{matrix} F_{d i s} = F_{e} F_{v} \end{matrix}$

(5)

式中:F_e为变形梯度的弹性分量;F_v为变形梯度的黏塑性分量.

此时可获得速度梯度L_dis:

$L_{d i s} = {\overset{\cdot}{F}}_{d i s} F_{d i s}^{- 1} = {\overset{\cdot}{F}}_{e} F_{e}^{- 1} + {\overset{\cdot}{F}}_{e} {\overset{\cdot}{F}}_{v} F_{v}^{- 1} F_{e}^{- 1} = L_{e} + {\tilde{L}}_{v}$

(6)

$\begin{matrix} {\tilde{L}}_{v} = {\overset{\cdot}{F}}_{e} {\overset{\cdot}{F}}_{v} F_{v}^{- 1} F_{e}^{- 1} = {\tilde{D}}_{v} + {\tilde{W}}_{v} \end{matrix}$

(7)

式中:L_e为速度梯度的弹性分量; ${\tilde{L}}_{v}$ 为速度梯度的黏塑性分量; ${\tilde{D}}_{v}$ 为黏塑性伸长率分量; ${\tilde{W}}_{v}$ 为旋转率分量.为了使变形梯度中的刚体旋转得到确定,可令旋转率分量 ${\tilde{W}}_{v}$ ≡0,该约束不会改变本构模型的表征结果^[20],此时 ${\tilde{L}}_{v}$ 可写为

$\begin{matrix} {\tilde{L}}_{v} = {\tilde{D}}_{v} = {\overset{\cdot}{F}}_{e} {\overset{\cdot}{F}}_{v} F_{v}^{- 1} F_{e}^{- 1} \end{matrix}$

(8)

2.2 TN模型

本文提出的TN模型包含3个并联网络,如图4所示.其中一个网络(网络A)为超弹性网络,用于表征未硫化橡胶的超弹性应力响应;另外两个网络为黏塑性网络(网络B、C),用于表征其黏塑性响应.不同于通常所采用的广义Maxwell模型,黏塑性网络B、C采用了基于高分子链微观运动学机制的Bergström-Boyce流动模型,能更好地表征未硫化橡胶的非线性黏塑性流动特征.同时,模型引入了基于分子链运动演化机制的损伤模型来表征未硫化橡胶的Mullins损伤软化效应.

图4

图4 TN模型框架

Fig.4 Frame of TN model

模型Cauchy应力σ_tol为3个网络Cauchy应力的和:

$\begin{matrix} σ_{t o l} = σ_{A} + σ_{B} + σ_{C} \end{matrix}$

(9)

2.2.1 TN模型应力响应

橡胶类材料的本构行为常通过应变能密度函数进行表征^[1,21],目前针对橡胶类高分子材料已经发展出多种不同形式的应变能密度函数,包括唯象应变能密度函数和基于分子网络的应变能密度函数^[22⇓-24].其中,Arruda等^[25]基于非高斯链分子网络所提出的八链模型参数简单,表征效果优异,应用较为广泛.考虑到未硫化橡胶具有相似的分子网络构型,在此采用八链模型应变能密度函数来表征未硫化橡胶的超弹应力响应,应变能密度表达式为

$\begin{matrix} W = n k_{B} T (\sqrt{N} λ_{c h a} β + N l n \frac{β}{s i n h β}) \end{matrix}$

(10)

$\begin{matrix} β = L^{- 1} (\frac{λ_{c h a}}{\sqrt{N}}) \end{matrix}$

(11)

式中: W为应变能密度;L为朗之万函数,L(x)=coth x-1/x;n为分子链密度(单位体积分子链的数量);N为分子链联接点(包括分子链的缠结点,分子链与填充颗粒的物理吸附点等)之间的分子链段数目; 为八链网络中的链伸长; $k_B$ 为玻尔兹曼常数;T为参考温度.

对于超弹性力学响应,Cauchy应力可用应变能密度函数^[18]表示为

$\begin{matrix} σ = \frac{2}{J} F \frac{\partial W}{\partial C} F^{T} \end{matrix}$

(12)

式中:C为右Cauchy-Green变形张量.将应变能密度函数表达式(10)和(11)代入上式,并将C以不变量的形式写入,即可获得八链模型的Cauchy应力表达式^[25]:

$\begin{matrix} σ = \frac{n k_{B} T}{3 J} \frac{\sqrt{N}}{λ_{c h a}} L^{- 1} (\frac{λ_{c h a}}{\sqrt{N}}) d e v (B^{*}) + K (J - 1) I \end{matrix}$

(13)

式中: $d e v (B^{*}) = B^{*} - \frac{1}{3} t r (B^{*}) I$ ;K为体积模量,用于表征材料的可压缩性.

对于本文中三网络框架,3个网络的Cauchy应力可通过八链模型分别表示为

$\begin{matrix} σ_{A} = \frac{n_{A} k_{B} T}{3 J} \frac{\sqrt{N_{A}}}{λ_{c h a}^{A}} L^{- 1} (\frac{λ_{c h a}^{A}}{\sqrt{N_{A}}}) d e v (B^{*}) + K (J - 1) I \end{matrix}$

(14)

$\begin{matrix} σ_{B} = \frac{n_{B} k_{B} T}{3 J} \frac{\sqrt{N_{B}}}{λ_{c h a}^{B}} L^{- 1} (\frac{λ_{c h a}^{B}}{\sqrt{N_{B}}}) d e v (B_{e B}^{*}) + K (J - 1) I \end{matrix}$

(15)

$\begin{matrix} σ_{C} = \frac{n_{C} k_{B} T}{3 J} \frac{\sqrt{N_{C}}}{λ_{c h a}^{C}} L^{- 1} (\frac{λ_{c h a}^{C}}{\sqrt{N_{C}}}) d e v (B_{e C}^{*}) + K (J - 1) I \end{matrix}$

(16)

超弹性网络(网络A)的Cauchy应力由整体变形梯度决定;黏塑性网络(网络B、C)的变形梯度可进行乘法分解为弹性分量和黏塑性分量:F=F_eBF_vB,F=F_eCF_vC,其Cauchy应力仅由弹性分量提供,即

$\begin{matrix} \begin{array}{l} d e v (B_{e B}^{*}) = B_{e B}^{*} - \frac{1}{3} t r (B_{e B}^{*}) I \\ B_{e B}^{*} = J^{- 2 / 3} B_{e B} = J^{- 2 / 3} F_{e B} F_{e B}^{T} \end{array}\} \end{matrix}$

(17)

$\begin{matrix} \begin{array}{l} d e v (B_{e C}^{*}) = B_{e C}^{*} - \frac{1}{3} t r (B_{e C}^{*}) I \\ B_{e C}^{*} = J^{- 2 / 3} B_{e C} = J^{- 2 / 3} F_{e C} F_{e C}^{T} \end{array}\} \end{matrix}$

(18)

式中:F_eB和F_eC分别为B、C网络的弹性变形梯度;F_vB和F_vC为B、C网络的黏塑性变形梯度.

$λ_{c h a}^{A}$ 、 $λ_{c h a}^{B}$ 、 $λ_{c h a}^{C}$ 为各自的弹性伸长量:

$\begin{matrix} λ_{c h a}^{A} = \sqrt{t r (B^{*}) / 3} \end{matrix}$

(19)

$λ_{c h a}^{B} = \sqrt{t r (B_{e B}^{*}) / 3}$

(20)

$\begin{matrix} λ_{c h a}^{C} = \sqrt{t r (B_{e C}^{*}) / 3} \end{matrix}$

(21)

2.2.2 黏塑性流动模型

要获得黏塑性网络B、C的弹性变形梯度,首先需要求解其黏塑性变形梯度,故需要合适的模型来对黏塑性流动进行表征.由1.2节实验结果可知,未硫化橡胶力学行为较为复杂,通常的Maxwell模型不足以较好地表征其非线性黏塑性,故本文采用了基于高分子链微观动力学机制的Bergström-Boyce流动模型.该流动模型认为橡胶类聚合物材料的黏塑性变形是分子链的布朗运动(Brownian Motion)和“爬行”运动(Reptation Motion)联合作用的结果^[18],基于此假设,模型具有以下的表达形式:

$\begin{matrix} {\overset{\cdot}{γ}}_{v} = (λ_{c h a} {- 1 + ξ)}^{η} {(\frac{‖τ‖}{\hat{τ}})}^{m} \end{matrix}$

(22)

式中: ${\overset{\cdot}{γ}}_{v}$ 为有效流动速率; $\hat{τ}$ 为剪切流动阻力;m为剪切流动指数;η为应变指数;ξ为应变修正因子;τ为施加的Kirchhoff应力张量;‖τ‖为驱动黏塑性流动的Kirchhoff应力张量的Frobenius范数, $‖ τ ‖ = ‖ d e v (σ) ‖ = \sqrt{t r (σ σ^{T})}$ 代数式 $(λ_{c h a} {- 1 + ξ)}^{η}$ 是由“爬行”运动驱动,用于表征与伸长量相关的有效黏度; ${(‖ τ ‖ / \hat{τ})}^{m}$ 则是用于表征历史相关的非线性响应,由布朗运动和“爬行”运动联合驱动.

黏塑性伸长率由 ${\overset{\cdot}{γ}}_{v}$ 和流动方向N_v共同决定:

$D_{v} = {\overset{\cdot}{γ}}_{v} N_{v}$

(23)

$\begin{matrix} N_{v} = \frac{τ}{‖τ‖} = \frac{d e v (σ)}{‖ d e v (σ) ‖} \end{matrix}$

(24)

将式(22)和(23)代入黏塑性变形速度梯度表达式(8),即可获得黏塑性变形梯度的时间导数:

$\begin{matrix} {\overset{\cdot}{F}}_{v} = {\overset{\cdot}{γ}}_{v} (F_{v})^{- 1} \frac{d e v (σ)}{‖ d e v (σ) ‖} F_{e} F_{v} \end{matrix}$

(25)

网络B、C中黏塑性变形梯度的时间导数分别为

$\begin{matrix} {\overset{\cdot}{F}}_{v B} = {\overset{\cdot}{γ}}_{v B} (F_{v B})^{- 1} \frac{d e v (σ)}{‖ d e v (σ) ‖} F_{e B} F_{v B} \end{matrix}$

(26)

$\begin{matrix} {\overset{\cdot}{F}}_{v C} = {\overset{\cdot}{γ}}_{v C} (F_{v C})^{- 1} \frac{d e v (σ)}{‖ d e v (σ) ‖} F_{e C} F_{v C} \end{matrix}$

(27)

对式(26)和(27)进行时间积分即可获得B、C网络的黏塑性变形梯度.

2.2.3 Mullins损伤效应

此TN模型从橡胶高分子链的分子演化理论^[26-27]出发来考虑未硫化橡胶的Mullins软化效应.对于颗粒填充聚合物材料,填充颗粒和高分子链之间的联接包括强化学联接(共价键)及弱物理联接^[28].在拉伸过程中,物理联接(填充颗粒与高分子链之间),以及部分化学联接(高分子链之间或填充颗粒与高分子链之间)达到极限伸长而发生断裂,引起软化现象.对于未硫化炭黑填充橡胶,物理联接的强度远远小于化学共价键联接,故可认为拉伸过程中发生断裂的化学共价键数量相对于物理联接点很少,可以忽略不计^[27].

TN模型中,A网络代表了由物理联接构成的炭黑颗粒填充高分子链网络,用于表征未硫化橡胶的超弹性响应,而B、C黏塑性网络代表了由化学共价键联接构成的自由分子链网络,用于表征未硫化橡胶的黏塑性响应.Mullins损伤主要来自于A网络中填充颗粒与高分子链之间物理联接的破坏,该破坏一方面减少了填充颗粒与高分子链之间的联接点数目,引起分子网络中N的增加,另一方面导致n减少^[26].物理联接的破坏程度与λ_cha有关,故n的减少由λ_cha决定,其变化情况可用A网络中切变模量G_A=n_Ak_BT的整体变化来表示^[27]:

$\begin{matrix} G_{A} = n_{A} k_{B} T = \frac{1 + C_{4} (λ_{c h a}^{A} - 1)}{C_{5} + C_{6} (λ_{c h a}^{A} - 1)} \end{matrix}$

(28)

N的变化同样与λ_cha有关:

$\begin{matrix} \sqrt{N_{A}} = \sqrt{\frac{1 + C_{7} (λ_{c h a}^{A} - 1)}{C_{8} + C_{9} (λ_{c h a}^{A} - 1)}} \end{matrix}$

(29)

式中:C₄~C₉均为材料常数.

将式(28)和(29)代入Cauchy应力表达式(14),可得包含应力软化效应的炭黑填充分子链网络的Cauchy应力表达式:

$\begin{array}{l} σ_{A} = \frac{1}{3 J} \frac{1 + C_{4} (λ_{c h a}^{A} - 1)}{C_{5} + C_{6} (λ_{c h a}^{A} - 1)} \times \\ \sqrt{\frac{1 + C_{7} (λ_{c h a}^{A} - 1)}{C_{8} + C_{9} (λ_{c h a}^{A} - 1)}} L^{- 1} (\frac{λ_{c h a}^{A}}{\sqrt{N_{A}}}) \times \\ d e v (B^{*}) + K (J - 1) I \end{array}$

(30)

化学共价键联接构成的自由链网络(网络B、C)对Mullins效应的影响主要来自于分子链的运动演化回复,并表现为胶料模量的增大.在此采用微观泡沫模型(MFM)^[29]的形式来计及分子链运动的影响,并将原模型中的缩减密度改写为应力回复因子R(取值应大于1,本文取为1.2)以体现分子链的运动回复对应力的增强效果,从而获得了切变模量与体积模量的表达式:

$\begin{matrix} G (G_{0}, R) = G_{0} {(\frac{R + α}{1 + α})}^{h} \end{matrix}$

(31)

$\begin{matrix} K (K_{0}, R) = K_{0} {(\frac{R + α}{1 + α})}^{h} \end{matrix}$

(32)

式中: $G_{0} 、 K_{0} 为网络初始切变模量和体积模量, G_{0} = n k_{B} T; α 和 h$ 为模量缩放因子.

将式(31)和(32)代入Cauchy应力表达式(15)和(16),可获得包含分子运动演化效应的B、C网络Cauchy应力表达式:

$\begin{matrix} σ_{B} = \frac{G_{0 B}}{3 J} {(\frac{R + α}{1 + α})}^{h} \frac{\sqrt{N_{B}}}{λ_{c h a}^{B}} L^{- 1} (\frac{λ_{c h a}^{B}}{\sqrt{N_{B}}}) \times d e v (B_{e B}^{*}) + K_{0 B} {(\frac{R + α}{1 + α})}^{h} (J - 1) I \end{matrix}$

(33)

$\begin{matrix} σ_{C} = \frac{G_{0 C}}{3 J} {(\frac{R + α}{1 + α})}^{h} \frac{\sqrt{N_{C}}}{λ_{c h a}^{C}} L^{- 1} (\frac{λ_{c h a}^{C}}{\sqrt{N_{C}}}) \times d e v (B_{e C}^{*}) + K_{0 C} {(\frac{R + α}{1 + α})}^{h} (J - 1) I \end{matrix}$

(34)

式中:G_0B、G_0C和K_0B、K_0C分别为B、C网络的初始切变模量和体积模量.

2.3 本构模型表征效果

曲线拟合即实验曲线与表征曲线之间误差最小化的过程,本文采用最小二乘法进行误差计算:

$\begin{matrix} m i n \overset{P}{\sum_{i = 1}} f (σ_{i}^{m o d} - σ_{i}^{e x p}) \end{matrix}$

(35)

式中:为本构模型的预测应力; 为实验应力值;P为实验数据点的数量.拟合时材料参数的初值是根据参数的实际物理意义以及相关文献报道综合确定.对于G、K及,根据材料实际力学性质,取值均需大于0.本文中给定G和的初值在 0~10 MPa 范围内, 考虑到未硫化橡胶是近似不可压缩的,K取值可相对较大,给定其初值为150 MPa.对于N、m及α、h等唯象参数,为防止拟合过程发散,借鉴了文献[30-31]中的初值选择范围;模量因子 $C_4$ ~ $C_9$ 的参数初值均取为1.表1中给出了拟合中各参数的大致取值区间.

表1 三网络模型拟合参数

Tab.1 Ultimate optimization parameters of NT model

参数名称	参数符号	参数值	参数取值区间
切变模量	G_0B, G_0C /MPa	1.562 2, 0.403 87	>0
链段密度	N_B, N_C	9.610 1	(1, 10)
体积模量	K/MPa	180.656	>0
流动阻力	${\hat{τ}}_{B}$ , ${\hat{τ}}_{C}$ /MPa	1.311 8, 5.810 4	>0
流动指数	m_B, m_C	2.340 1, 5.610 7	(1, 20)
B、C网络模量缩放因子	α, h	1.254 2, 2.681 4	(1, 10)
A网络模量因子	C₄, C₅, C₆	0.367 4, 17.185 9, 2.292 1	—
A网络模量因子	C₇, C₈, C₉	1.048 1, 0.354 1,-0.052 4	—

新窗口打开| 下载CSV

模型参数拟合过程采用梯度下降法,拟合结果如图5所示.图中:ε_true为真应变;σ_true为真应力.可以看出,不同应变率下的拟合曲线与实验曲线吻合良好, 未硫化橡胶的复杂力学行为特征, 包括迟滞效应、残余变形及软化现象均得到了较好的表征.4种不同应变率下的表征结果在同一套拟合参数下完成,拟合参数如表1所示.

图5

图5 TN模型表征结果

Fig.5 Fitting result using TN model

3 TN模型的有限元实现及验证

获得的本构模型若要应用于工程实践,则需要将其进行数值实现.本文依托于有限元软件Abaqus,完成了TN模型用户材料子程序(UMAT)的开发,并进行了仿真计算以验证模型的有效性和数值稳定性.

模型有限元实现过程中采用的应力更新算法为二阶精度的显式龙格库塔算法(Explicit Midpoint Rule)^[32],该算法流程清晰,灵活性强,能够较好地平衡计算精度、数值稳定性以及易数值植入性.采用t代表前一个增量步,t+1代表后一个增量步,则该应力更新算法步骤如下.

(1) t时刻从Abaqus中读取各变形梯度张量:F_t、 $F_{t}^{v B}$ 、 $F_{t}^{v C}$ 、F_t₊₁.

(2) 由应力计算公式(14)~(16)和黏塑性流动法则式(22)、(23),计算t+1/2时刻黏塑性网络B、C的黏塑性变形梯度 $F_{t + 1 / 2}^{v B}$ 、 $F_{t + 1 / 2}^{v C}$ .

(3) 对t时刻和t+1时刻总变形梯度平均获得t+1/2时刻总变形梯度 $F_{t + 1 / 2}$ .

(4) 根据步骤(2)、(3)的计算结果,计算t+1/2时刻B、C网络的弹性变形梯度 $F_{t + 1 / 2}^{e B}$ 、 $F_{t + 1 / 2}^{e C}$ .

(5) 根据应力公式(14)~(16)和式(22)、(23),以及F_t₊₁、 $F_{t + 1 / 2}^{e B}$ 、 $F_{t + 1 / 2}^{e C}$ ,计算黏塑性网络B、C在 t+1 时刻的黏塑性变形梯度 $F_{t + 1}^{v B}$ 、 $F_{t + 1}^{v C}$ .

(6) 检验B、C网络的黏塑性变形梯度是否满足误差准则,若满足则进入第(7)步,否则更新时间步长后返回第(2)步.

(7) 计算t+1时刻B、C网络的弹性变形梯度并根据式(14)~(16)计算t+1时刻Cauchy应力.

(8) 储存自定义状态变量并计算一致切线刚度矩阵(Jacobian矩阵).

获得TN模型用户材料子程序后,为验证其有效性,首先采用此子程序对一个单元在不同工况下(拉伸、压缩、剪切)的力学行为进行了数值计算,各工况下单元的Mises应力(S)分布情况如图6所示.图7给出了在不同应变率下拉伸工况的仿真结果与测试结果的对比情况.可以看出,仿真结果能准确再现不同应变率下的拉伸应力应变曲线.

图6

图6 一个单元的仿真结果

Fig.6 Simulation results of a single element

图7

图7 拉伸仿真结果与单调拉伸实验结果对比

Fig.7 Comparison of simulation results with experiments under uniaxial tensile conditions

为了进一步验证数值仿真精度,开展了未硫化橡胶在 $\overset{\cdot}{ε}$ =0.03 s^-1下的多段拉伸松弛加卸载实验,即首先将试件拉伸至40%应变,松弛30 s后将试件拉伸至60%应变,再松弛30 s后拉伸至80%应变,再逐步将应变卸载至60%、40%.为了对比仿真结果与实验结果,按照标准试件尺寸建立了试件的有限元模型,如图8(a)所示.试件包括 2 079 个节点,1 216 个单元,单元类型为C3D8H.仿真过程中固定一端的加持段,对另一端施加位移边界条件.最终获得的仿真结果与测试结果对比情况如图8(b)所示,可以看出,仿真结果与测试结果吻合非常好.这表明TN模型能较好表征未硫化橡胶材料拉伸过程中的各种非线性黏弹塑性变形特征,其残余变形、迟滞效应以及松弛段内的应力恢复情况都得到了较好的体现.

图8

图8 多段松弛加卸载实验与仿真结果对比

Fig.8 Comparison of multistep tensile relaxation test and simulation result

另一方面,未硫化橡胶多用于成型工艺,此过程中胶料需经历较为复杂的变形,故本构模型在复杂变形工况下的计算稳定性尤为重要.为了验证TN模型的数值稳定性,对轮胎胎面胶压入模具的成型过程进行了数值仿真,如图9所示.此模型实际上是对轮胎胎面胶纵向花纹沟成型工艺的简化,如图9(a)所示,模型轴对称分布,且仅包含一条花纹沟.采用轴对称有限元模型对此过程进行仿真,如图9(b)所示,将胎面胶分为5层以便于观察胶料流动,总厚度为20 mm,宽度为80 mm,采用的单元类型为CAX4H,节点数为 2 721 个,单元数为 2 594 个.模具采用离散刚体,深度为15 mm,中间的拱形突起高8.75 mm,宽7.5 mm.对胎面胶上层节点施加径向向下的位移边界条件,速度为0.5 mm/s,橡胶与模具之间的接触摩擦因数设为0.5.胎面胶压入模具后的各层胶料分布如图9(c)所示,由于未硫化橡胶具有近似不可压缩性,胶料充满模具后由两侧溢出.图9(d)和9(e)给出了此时刻胎面胶的应力和应变分布情况,可以看出,中间的拱形突起处胶料流动较为剧烈,存在一定的应力应变集中,靠近模具边缘处的应力应变分布高于底部平坦部位.胶料流动情况以及应力应变分布情况与文献[6,8]中实验及仿真结果基本一致.综上可见,TN模型在复杂变形工况下仍能给出合理的计算结果,表明TN模型具有良好的数值稳定性.

图9

图9 胎面胶压入模具过程仿真

Fig.9 Simplified molding process of a tire tread

4 结论

(1) 对未硫化橡胶开展了不同应变率下( $\overset{\cdot}{ε}$ =0.001, 0.003, 0.03, 0.1 s^-1)的单调拉伸和循环拉伸加卸载实验.结果表明,未硫化橡胶的拉伸力学行为具有明显的应变率相关性,拉伸应力水平随着应变率增加而显著提高.循环加卸载工况下,胶料表现出明显的迟滞效应、残余变形及Mullins损伤效应,且随着应变率增加,迟滞效应增强,残余变形减小,反映出了较为复杂的非线性黏塑性特征.

(2) 根据实验结果,提出一个用于表征未硫化橡胶黏弹塑性力学行为的TN本构模型,模型由一个超弹性网络和两个非线性黏塑性网络构成.3个网络的应力响应均基于超弹八链模型,同时在黏塑性网络中采用了Bergström-Boyce流动法则描述胶料的非线性黏塑性流动.此外,模型从橡胶高分子链的分子演化理论出发考虑了未硫化橡胶的 Mullins 损伤软化效应.采用TN模型对未硫化橡胶循环拉伸加卸载实验曲线进行表征,表征结果与实验结果吻合良好,胶料的各种拉伸力学特征均得到了较好的体现.

(3) 依托有限元仿真软件Abaqus,完成了TN模型用户材料子程序(UMAT)的开发,并对试件单调拉伸和多段松弛拉伸加卸载实验进行了仿真,仿真结果与实验结果吻合较好,验证了此TN模型的有效性.最后,采用TN模型对轮胎胎面胶压入模具的过程进行数值计算,进一步验证了模型在复杂变形工况下良好的数值稳定性.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

MARK

J E

, ERMAN

, ROLAND

C M

, et al.

The science and technology of rubber[M]. 4th ed. Waltham: Academic Press, 2013: 337-340.

[2]

董义军.

农业子午线轮胎成型胎坯质量缺陷原因分析与解决措施

[J]. 轮胎工业, 2021, 41(5): 327-330.

DONG

Yijun

Cause analysis and solutions of green tire quality defects of agricultural radial tire

[J]. Tire Industry, 2021, 41(5): 327-330.

DOI:10.1016/j.trpro.2019.07.020 URL [本文引用: 1]

[3]

ANDRZEJ

, MICHAL

, SEWERYN

, et al.

Characteristics and investigation of selected manufacturing defects of passenger car tires

[J]. Transportation Research Procedia, 2019, 40: 119-126.

[4]

KALISKE

, ZOPF

, BRÜGGEMANN

Experimental characterization and constitutive modeling of the mechanical properties of uncured rubber

[J]. Rubber Chemistry and Technology, 2010, 83(1): 1-15.

DOI:10.5254/1.3548264 URL [本文引用: 4]

The properties of uncured rubber are characterized by a viscoelastic material formulation in order to develop a finite element method (FEM) to predict the deformation behavior of this material during processing. As material formulation, a viscoelastic material model is considered, which consists of a nonlinear Hooke spring connected in parallel to a finite number of Maxwell elements. To identify the material parameters, various tests have to be conducted. The chosen test procedure and data analysis strategy are presented. An evolutionary optimization procedure is used to fit the material parameters to the measurements by minimizing the mean square error of the approximation. At the end, the suitability of the chosen material model and the identified material parameters is shown. Finally, the result of a molding test is compared to the corresponding FEM simulation.

[5]

FENG

X J

, WEI

Y T

, LI

Z C

, et al.

Research on nonlinear viscoelastic constitutive model for uncured rubber

[J]. Engineering Mechanics, 2016, 33: 212-219.

[本文引用: 3]

[6]

DAL

, ZOPF

, KALISKE

Micro-sphere based viscoplastic constitutive model for uncured green rubber

[J]. International Journal of Solids and Structures, 2018, 132-133: 201-217.

[本文引用: 4]

[7]

ZOPF

, KALISKE

Numerical characterisation of uncured elastomers by a neural network based approach

[J]. Computers & Structures, 2017, 182: 504-525.

DOI:10.1016/j.compstruc.2016.12.012 URL [本文引用: 2]

[8]

ZOPF

, HOQUE

S E

, KALISKE

Comparison of approaches to model viscoelasticity based on fractional time derivatives

[J]. Computational Materials Science, 2015, 98: 287-296.

DOI:10.1016/j.commatsci.2014.11.012 URL [本文引用: 4]

[9]

CHEN

, ZHOU

, et al.

Radial tire construction design method based on finite element simulation

[J]. Journal of Donghua University, 2017, 34(03): 150-157.

[10]

王国林, 周伟, 周海超, 等.

子午线轮胎二次法成型过程仿真

[J]. 机械设计与制造, 2017(7): 179-182.

WANG

Guolin

, ZHOU

Wei

, ZHOU

Haichao

, et al.

Simulation of the radial tire second-stage building process

[J]. Machinery Design & Manufacture, 2017(7): 179-182.

DOI:10.2346/tire.18.460405 URL [本文引用: 2]

[11]

ZHOU

, WANG

Wide-base tire-building process and design optimization using finite element analysis

[J]. Tire Science and Technology, 2018, 46(4): 242-259.

The wide-base tire is a relatively new design that originated to replace dual tires because of its potential for improved performance. However, during the construction process, the wide-base tire is more likely to experience tread deformation and uneven stress distribution. The goal of this study is to incorporate numeric techniques for the construction and design optimization of a wide-base, heavy vehicle, pneumatic tire. First, four conditions of the tire (385/55R22.5)–building process, including gluing of components on the main drum, gluing of components on the auxiliary drum, green tire, and finalizing the capsule vulcanizing machine, were simulated using finite element analysis. Second, to solve the difference in the tire's (435/50R19.5) material distribution between the real manufactured structure and the theoretical structure, the curved surface drum-building method and the parameters of the curved surface drum were determined by tire construction simulation. In this article, we present the method for collecting tire material, the measurement process, the analysis method, some general results, and statistics on the wide-base tire. Finally, validation of results of the simulation and measurement are given.

[12]

KIM

, BERGER

, KALISKE

Strain rate-dependent behavior of uncured rubber: Experimental investigation and constitutive modeling

[DB/OL]. (2022-01-04)[2022-04-05]. https//doi.org/10.5254/rct.21.78981.

[本文引用: 3]

[13]

, WANG

, LI

, et al.

Experimental investigation and constitutive modeling of uncured carbon black filled rubber at different strain rates

[J]. Polymer Testing, 2019, 75: 117-126.

DOI:10.1016/j.polymertesting.2019.02.005 URL [本文引用: 3]

[14]

WANG

, LI

, et al.

Effect of the temperature and strain rate on the tension response of uncured rubber: Experiments and modeling

[J]. Mechanics of Materials, 2020, 148: 103480.

DOI:10.1016/j.mechmat.2020.103480 URL [本文引用: 2]

[15]

GUO

, WANG

High-rate tensile behavior of silicone rubber at various temperatures

[J]. Rubber Chemistry and Technology, 2020, 93(1): 183-194.

DOI:10.5254/rct.19.81562 URL [本文引用: 1]

The effects of strain rate and temperature on the tensile behavior of silicone rubber were investigated. The quasi-static uniaxial tensile experiments were conducted using an electromechanical testing system, and the high-rate uniaxial tensile tests were performed employing a modified split Hopkinson tension bar technique for low-strength and low-impedance materials. The tensile responses were obtained at strain rates of 0.001–1400 s−1 and temperatures ranging from −50 to 50 °C. The experiments reveal that the tensile stress–strain behavior of silicone rubber is nonlinear and highly dependent on strain rate and temperature. The values of stiffness and nominal stress at a given elongation increase with increased strain rate and decrease with increasing temperature. It is appropriate to postulate that the tensile response at high strain rates arises from the combination of hyperelasticity and viscoelasticity. According to the incompressibility assumption, a phenomenologically inspired visco-hyperelastic model was proposed to describe the constitutive behavior of silicone rubber over wide ranges of strain rates and temperatures.

[16]

魏明杰, 王银龙, 刘敏, 等.

温度影响下炭黑增强橡胶复合材料的循环拉伸力学行为

[J]. 实验力学, 2020, 35(6): 1030-1040.

WEI

Mingjie

, WANG

Yinlong

, LIU

Min

, et al.

Temperature effect on the mechanical behavior of carbon black reinforced rubber composites under cyclic tensile loadings

[J]. Journal of Experimental Mechanics, 2020, 35(6): 1030-1040.

DOI:10.1016/j.ijplas.2017.10.011 URL [本文引用: 1]

[17]

GUO

, ZAÏRI

, GUO

A thermo-viscoelastic-damage constitutive model for cyclically loaded rubbers. Part I: Model formulation and numerical examples

[J]. International Journal of Plasticity, 2018, 101: 106-124.

[18]

BERGSTRÖM

J S

Mechanics of solid polymers: Theory and computational modeling[M]. Norwich: William Andrew Publishing, 2015: 141-150.

[本文引用: 3]

[19]

JARRAH

H R

, ZOLFAGHARIAN

, HEDAYATI

, et al.

Nonlinear finite element modelling of thermo-visco-plastic styrene and polyurethane shape memory polymer foams

[J]. Actuators, 2021, 10(3): 46.

DOI:10.3390/act10030046 URL [本文引用: 1]

This paper presents nonlinear finite element (FE) models to predict time- and temperature-dependent responses of shape memory polymer (SMP) foams in the large deformation regime. For the first time, an A SMP foam constitutive model is implemented in the ABAQUS FE package with the aid of a VUMAT subroutine to predict thermo-visco-plastic behaviors. A phenomenological constitutive model is reformulated adopting a multiplicative decomposition of the deformation gradient into thermal and mechanical parts considering visco-plastic SMP matrix and glass microsphere inclusions. The stress split scheme is considered by a Maxwell element in parallel with a hyper-elastic rubbery spring. The Eyring dashpot is used for modelling the isotropic resistance to the local molecular rearrangement such as chain rotation. A viscous flow rule is adopted to prescribe shear viscosity and stress. An evolution rule is also considered for the athermal shear strengths to simulate macroscopic post-yield strain-softening behavior. In order to validate the accuracy of the model as well as the solution procedure, the numerical results are compared to experimental responses of Styrene and Polyurethane SMP foams at different temperatures and under different strain rates. The results show that the introduced FE modelling procedure is capable of capturing the major phenomena observed in experiments such as elastic and elastic-plastic behaviors, softening plateau regime, and densification.

[20]

BOYCE

M C

, WEBER

G G

, PARKS

D M

On the kinematics of finite strain plasticity

[J]. Journal of the Mechanics and Physics of Solids, 1989, 37(5): 647-665.

DOI:10.1016/0022-5096(89)90033-1 URL [本文引用: 1]

[21]

CARROLL

M M

Molecular chain networks and strain energy functions in rubber elasticity

[J]. Philosophical Transactions of the Royal Society A, 2019, 377(2144): 20180067.

DOI:10.1002/msd2.v1.1 URL [本文引用: 1]

[22]

MELLY

S K

, LIU

, et al.

A review on material models for isotropic hyperelasticity

[J]. International Journal of Mechanical System Dynamics, 2021, 1(1): 71-88.

[23]

XIANG

, ZHONG

, RUDYKH

, et al.

A review of physically based and thermodynamically based constitutive models for soft materials

[J]. Journal of Applied Mechanics, 2020, 87(11): 110801.

DOI:10.1115/1.4047776 URL [本文引用: 1]

In this paper, we review constitutive models for soft materials. We specifically focus on physically based models accounting for hyperelasticity, visco-hyperelasticity, and damage phenomena. For completeness, we include the thermodynamically based viscohyperelastic and damage models as well as the so-called mixed models. The models are put in the frame of statistical mechanics and thermodynamics. Based on the available experimental data, we provide a quantitative comparison of the hyperelastic models. This information can be used as guidance in the selection of suitable constitutive models. Next, we consider visco-hyperelasticity in the frame of the thermodynamic theory and molecular chain dynamics. We provide a concise summary of the viscohyperelastic models including specific strain energy density function, the evolution laws of internal variables, and applicable conditions. Finally, we review the models accounting for damage phenomenon in soft materials. Various proposed damage criteria are summarized and discussed in connection with the physical interpretations that can be drawn from physically based damage models. The discussed mechanisms include the breakage of polymer chains, debonding between polymer chains and fillers, disentanglement, and so on.

[24]

GILLES

, ERWAN

Comparison of hyperelastic models for rubber-like materials

[J]. Rubber Chemistry and Technology, 2006, 79(5): 835-858.

DOI:10.5254/1.3547969 URL [本文引用: 1]

The present paper proposes a thorough comparison of twenty hyperelastic models for rubber-like materials. The ability of these models to reproduce different types of loading conditions is analyzed thanks to two classical sets of experimental data. Both material parameters and the stretch range of validity of each model are determined by an efficient fitting procedure. Then, a ranking of these twenty models is established, highlighting new efficient constitutive equations that could advantageously replace well-known models, which are widely used by engineers for finite element simulation of rubber parts.

[25]

ARRUDA

E M

, BOYCE

M C

A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials

[J]. Journal of the Mechanics and Physics of Solids, 1993, 41(2): 389-412.

DOI:10.1016/0022-5096(93)90013-6 URL [本文引用: 2]

[26]

MARCKMANN

, VERRON

, GORENT

, et al.

A theory of network alteration for the Mullins effect

[J]. Journal of the Mechanics and Physics of Solids, 2002, 50(9): 2011-2028.

DOI:10.1016/S0022-5096(01)00136-3 URL [本文引用: 2]

[27]

KHAJEHSAEID

Development of a network alteration theory for the Mullins-softening of filled elastomers based on the morphology of filler-chain interactions

[J]. International Journal of Solids and Structures, 2016, 80: 158-167.

DOI:10.1016/j.ijsolstr.2015.10.032 URL [本文引用: 3]

[28]

VYAZOVKIN

, SBIRRAZZUOLI

Isoconversional kinetic analysis of thermally stimulated processes in polymers

[J]. Macromolecular Rapid Communications, 2010, 27(18): 1515-1532.

DOI:10.1002/(ISSN)1521-3927 URL [本文引用: 1]

[29]

ROBERTS

A P

, GARBOCZI

E J

Elastic properties of model random three-dimensional open-cell solids

[J]. Journal of the Mechanics and Physics of Solids, 2002, 50(1): 33-55.

DOI:10.1016/S0022-5096(01)00056-4 URL [本文引用: 1]

[30]

BERGSTRÖM

J S

A library of advanced user materials: Version 5.0.0

[EB/OL]. (2018-04-17) [2022-04-05]. http//PolyUMod.com/.

URL [本文引用: 1]

[31]

DREHER

M L

, NAGARAJA

, BERGSTRÖM

J S

, et al.

Development of a flow evolution network model for the stress-strain behavior of poly(L-lactide)

[J]. Journal of Biomechanical Engineering, 2017, 139(9): 091002.

DOI:10.1115/1.4037071 URL [本文引用: 1]

Computational modeling is critical to medical device development and has grown in its utility for predicting device performance. Additionally, there is an increasing trend to use absorbable polymers for the manufacturing of medical devices. However, computational modeling of absorbable devices is hampered by a lack of appropriate constitutive models that capture their viscoelasticity and postyield behavior. The objective of this study was to develop a constitutive model that incorporated viscoplasticity for a common medical absorbable polymer. Microtensile bars of poly(L-lactide) (PLLA) were studied experimentally to evaluate their monotonic, cyclic, unloading, and relaxation behavior as well as rate dependencies under physiological conditions. The data were then fit to a viscoplastic flow evolution network (FEN) constitutive model. PLLA exhibited rate-dependent stress–strain behavior with significant postyield softening and stress relaxation. The FEN model was able to capture these relevant mechanical behaviors well with high accuracy. In addition, the suitability of the FEN model for predicting the stress–strain behavior of PLLA medical devices was investigated using finite element (FE) simulations of nonstandard geometries. The nonstandard geometries chosen were representative of generic PLLA cardiovascular stent subunits. These finite element simulations demonstrated that modeling PLLA using the FEN constitutive relationship accurately reproduced the specimen’s force–displacement curve, and therefore, is a suitable relationship to use when simulating stress distribution in PLLA medical devices. This study demonstrates the utility of an advanced constitutive model that incorporates viscoplasticity for simulating PLLA mechanical behavior.

[32]

TOMAS

, CISILINO

A P

, FRONTINI

P M

Accurate, efficient and robust explicit and implicit integration schemes for the Arruda-Boyce viscoplastic model

[J]. Mecnica Computacional, 2008, 14: 1003-1042.