考虑缠结效应的超弹性本构模型1)

肖 锐 向玉海 钟旦明 曲绍兴

(浙江大学工程力学系,杭州 310027)

引言

建立超弹性模型对于橡胶和其他软材料的工程应用具有重要的意义[1-2].当前关于软材料超弹性理论模型主要包括唯象模型和分子统计理论模型[3-4].唯象模型的建立主要是基于连续介质力学以及软材料的实验所观察到的宏观力学行为.常见的唯象模型包括Neo-Hooken 模型、Mooney-Rivlin 模型、Ogden模型、Gent 模型、Yeoh 模型等[3].其中Neo-Hookean模型、Gent 模型和Yeoh 模型都是柯西−格林应变张量第一不变量I1的函数,Mooney-Rivlin 模型同时依赖于柯西−格林应变张量的第一和第二不变量(I1,I2),而Ogden 模型是3 个主方向拉伸比的函数.这些唯象模型在描述橡胶材料的超弹性变形行为取得了一定的成功[5-7],例如,张希润等[5]采用压入法研究软材料的力学性能,结果显示基于Mooney-Rivlin 模型的有限元结果能与压入实验结果吻合.Sasso 等[6]发现Ogden 模型和二阶的Mooney-Rivlin 模型可以描述橡胶在单轴、等双轴和鼓泡实验中的力学响应.

针对唯象超弹性模型的不足之处,部分学者[8-9]对上述模型作了扩展,例如,李雪冰等[8]在Yeoh 模型的自由能基础上增加了依赖于应变张量第二不变量的部分,克服了经典Yeoh 模型在预测等双轴曲线时“偏软”的特性.近年来,国内外相关学者[10-13]也提出了多种新型的超弹性唯象模型,例如,Lopez-Pamies[10]提出了基于应变第一不变量的指数型自由能形式,Mangan 等[11]提出了同时依赖于柯西−格林应变第一和第二不变量的新型超弹性自由能形式,该模型能准确表征大变形条件下简单剪切变形切应力和切应变的线性关系,以及单轴变形下的应变强化行为.

另一类被广泛用来描述软材料超弹性行为的模型是基于分子统计理论模型,其主要建模思路是将统计理论应用于分子链网络中,进而得出变形自由能函数.常见的分子统计理论模型包括Treloar 高斯链模型[14],Wang-Guth 三链模型[15]、Arruda-Boyce 八链模型[16]以及全链网络模型[17].Treloar[14]提出的高斯链分布模型,与唯象的Neo-Hookean 模型,在数学形式上是完全等价的,该模型只能描述轻度交联橡胶材料的应力响应,而对于交联度高的橡胶网络,这类模型仅能描述小应变范围内的应力响应[18].Boyce等[19]发现唯象的Gent 模型和Arruda-Boyce 模型在多种加载条件下的应力−应变关系是接近的,这两个模型均可较为准确地模拟不同交联度橡胶材料的单轴加载条件下的应力响应.Wu 和van der Giessen[17]建立了全链模型,全链模型需要对各个方向进行积分,计算起来不够方便.上述全链模型是基于仿射假设,而Diani 和Tallec[20]通过能量最小化原则建立了非仿射的全链超弹性模型.Hossain 等[21]针对八链模型和全链模型的不足,详细评估了几类扩展的八链和全链模型在模拟复杂加载情况下的表现.

针对Treloar[22]关于硫化橡胶的单轴、纯剪切和等双轴加载的经典实验数据,Steinmann 等[23]详细评估了14 种唯象模型和分子统计理论模型的预测能力,发现这些广泛应用的模型并不能很好地预测硫化橡胶在多种变形模式下的力学响应.针对经典超弹性模型的不足之处,部分学者[24-28]提出了引入缠结约束效应,其主要的建模思路是将缠结约束效应用Doi和Edwards[29]所提出的管模型(tube model)来近似表征.相关模拟结果显示,引入缠结约束效应的超弹性模型能更为准确地预测橡胶的力学行为,但上述模型也有相应的不足之处,例如,Miche 等[25]所建立的考虑缠结效应的超弹性模型是基于非仿射假设,需开展全方向积分计算,而Davodson 和Goubourne[27]以及Xiang 等[28]所建立的缠结约束模型的自由能是主方向拉伸比的函数,这导致这些模型有限元植入方面不够方便.

基于微观变形机制,本文建立了只依赖于应变不变量的超弹性模型,该模型同时包括熵弹性和缠结约束效应.并通过与文献中的实验数据进行对比验证了本文所建立模型的预测能力.最后,本文比较了多个缠结约束模型在不同加载模式下自由能的变化.

1 本构理论

1.1 单链变形自由能

针对单根分子链的变形行为,这里采用Miehe等[25]提出的微观力学模型来进行描述.该模型的自由能包括两部分,第一部分是源于分子链自由变形所引起的熵弹性自由能,第二部分是由于周边分子链缠结约束引起的自由能

采用Kuhn 和Gr¨un[30]提出的Langevin 统计模型来表示熵弹性所引起的自由能,其具体形式如下

其中,N是Kuhn 链段数量,λ 是单链的伸长比,£−1是Langevin 逆函数,其中Langevin 函数的定义为£(x)=cothx−1/x,Ψf0是常数.

采用Doi 和Edwards[29]提出的管模型(图1)来表征周边分子链缠结约束导致的变形自由能

其中,α 是与管的截面积形状相关的系数,b是Kuhn长度,d0是管的初始直径,υ 是管的截面积收缩比,Ψc0是常数.

图1 高分子链变形的管模型Fig.1 Tube model for polymer chains

1.2 微观−宏观关系

在单链变形自由能的基础上建立宏观变形自由能,需首先确定宏观和微观变形之间的映射关系,这里采用微球模型和仿射假设来获得相应的关系.如图2 所示,认为初始未变形的分子链均匀分布在一个微球表面上,施加变形后,微球在主应变方向上变形为一椭球,其变形梯度张量可表示为

其中λ1,λ2,λ3为3 个主方向上的拉伸比.变形后,分子链的长度变为

其中φ 是在初始构型中分子链在XOY平面的投影与X轴的夹角,而θ 是分子链与Z轴的夹角(见图2).

图2 基于仿射假设的微球模型Fig.2 Microsphere model based on affine assumption

这里假设微球表面所有分子链长度的均方值为宏观变形条件下分子链等效长度

采用同样方式,将变形后微球微单元的面积与初始面积比的均方值假定为宏观变形管面积收缩比[26],其中微球上无限小单元变形后面积dS与初始未变形面积dS0的比ξ 可表示为[31]

相应的面积比均方值为

1.3 宏观变形本构关系

结合单链自由能函数和微观−宏观映射关系,采用式(6)中的代替式(2)中的伸长比λ,采用式(8)的代替式(3)中的微管截面积收缩比υ,并假定单位体积内的分子链数目为n.则宏观的Helmholtz 自由能可表达为

其中Gc=nkT是分子链特征刚度,而Ge=αkTN(b/d0)2/3 是由于缠结约束所引起的剪切模量,Ψ0为常数,其作用是保证在初始未变形条件下自由能为0.这里在缠结约束的自由能部分添加了一个系数3,其主要原因是使得在小变形纯剪切条件下,缠结约束效应的剪切模量为Ge,而当N≫1 时,在小变形纯剪切变形时,熵弹性部分的剪切模量近似为Gc.

在主方向上,柯西应力可表示为

其中p是静水压部分,以保证体积不可压.将式(9)代入式(10)可得

在单轴和双轴变形过程中,可假定变形是加载在e1和e2方向上,在e3方向上应力为0,可依据这一条件将静水压p消除,进而得到在加载方向上的应力

虽然Arruda 和Boyce[16]在推导八链模型时采用了不同的方法,但最终得到的自由能形式与本文的熵弹性部分是完全一致的,因此本文的模型可认为是在八链模型的基础上进一步引入了分子链的缠结约束效应.当分子链的Kuhn 链段长度很大时,八链模型会退化为Neo-Hookean 模型.因此,这里也给出以Neo-Hookean 模型作为熵弹性部分,在单轴和双轴加载条件下的应力表达式

相应的名义应力可表示为

2 结果

采用上述建立的本构理论来模拟文献中关于橡胶材料在单轴和双轴变形条件下的应力响应,并与相应的Arruda-Boyce 八链模型和Neo-Hookean 模型进行对比分析.首先模拟Treloar[22]开展的关于硫化橡胶和Zhao[32]开展的关于TPE 橡胶的实验,上述工作对同一橡胶均同时开展了单轴拉伸、纯剪切和等双轴拉伸实验.其中纯剪切实验的实验设置与单轴拉伸实验是一致的,其主要区别在于所采用的试样的宽度远大于长度(一般来说,宽度/长度>10),在该条件下,可认为宽度方向的拉伸比为1,长度方向拉伸比为λ1,厚度方向拉伸比为1/λ1.值得注意的是,纯剪切变形和简单剪切变形是不一样的,具体可参加Thiel 等[33]的相关工作.

在单轴拉伸、纯剪切和等双轴变形情况下,3 个主方向的拉伸比可依次表示为

这里采用直接对实验数据进行拟合来获得模型参数,具体是获得以下方程的最小值

图3 给出了本文所建立的模型和八链模型对Treloar[22]所开展的关于硫化橡胶的实验的模拟结果,相应的模型参数和模拟误差见表1.由图3(a)可见,考虑了缠结效应的超弹性模型能同时准确的模拟单轴、纯剪切和等双轴拉伸的实验结果,与之相比,八链模型虽然能较好地预测硫化橡胶在单轴拉伸的应力响应,但对等双轴拉伸的结果的预测明显偏低(图3(b)).此外,考虑缠结效应的模型也能更好地体现单轴和纯剪切实验中应力的区别.图4(a)和图4(b)给出了关于TPE 橡胶的实验和相应的模拟结果,其主要结论与硫化橡胶是一致的.表1 的实验结果和模拟结果的误差值也进一步证实本文所建立的本构模型的预测能力明显好于八链模型.

针对Treloar 的实验数据,Treloar[34]和Steinmann等[23]系统地评估了经典唯象模型(Neo-Hookean 模型,Ogden 模型,Mooney-Rivlin 模型,Gent 模型,Yeoh模型等)和基于分子统计模型(三链模型、四链模型、八链模型、仿射模型等)的预测能力,结果显示这些模型均不能同时较好地预测单轴、纯剪切和等双轴的实验结果.为了与本文所建立的缠结模型进行对比,进一步采用以下两个三参数模型来模拟TPE 橡胶的实验结果.第一个模型是Yeoh 模型,其自由能形式为

其中c1,c2以及c31是Yeoh 模型的3 个参数.第二个模型是广义Rivlin 模型,其自由能形式如下

图3 关于硫化橡胶的实验和模拟结果对比Fig.3 Comparison between the experimentally measured and simulated results of vulcanized rubber

其中,c1,c2以及c32是广义Rivlin 模型的3 个参数.当c2=0,该模型退化为Mooney-Rivlin 模型.

表1 基于分子统计理论模型参数和模拟误差Table 1 Parameters and simulation error for models based on statistical mechanics

表2 给出了Yeoh 模型和广义Rivlin 模型的参数以及模拟误差,其误差均明显大于本文所建立的考虑缠结效应的超弹性模型.由图4(c)和图4(d)可见,Yeoh 模型能较好地预测单轴实验数据,但与八链模型类似,不能较好地模拟等双轴加载条件下应力响应.在所有模型中,广义Rivlin 模型的模拟效果最差,在单轴变形的大应变区间和双轴变形的小应变区间与实验结果都有较大偏差.由于Mooney-Rivlin 模型是该模型的一个特例,可预见Mooney-Rivlin 模型也不能较好的预测相关实验结果.

图4 关于TPE 橡胶的实验和模拟结果对比Fig.4 Comparison between the experimentally measured and simulated results of TPE

表2 唯象模型参数和模拟误差Table 2 Parameters and simulation error for phenomenological models

本文所建立的模型被进一步用来模拟Jones 和Treloar[35]所开展的关于另一类硫化橡胶(记为硫化橡胶B) 的实验结果.其实验设计如下:橡胶试样先在e2方向上被拉伸至一固定伸长比,然后在e1方向上开展拉伸实验,在e3方向上试样不受力.图5 给出了相应的实验结果,其中,当X轴的伸长比的数值大于时,X轴对应的是λ1,Y轴对应σ1−σ2,而当X轴的数值小于时,X轴对应λ3,而Y轴对应−σ2.在拟合过程中,我们发现相对误差值随着拉伸极限N的增大而降低,因此这里选择在Neo-Hookean 的模型的基础上考虑缠结效应,即采用式(13) 来拟合实验结果,相关的模型参数见表1.如图5(a) 所示,考虑了缠结效应的本构模型能准确地预测实验结果,而不考虑缠结效应的Neo-Hookean模型只能在一定的变形范围内与实验结果相吻合,而总体的模拟效果较差.在本文中,缠结约束部分是基于柯西−格林应变张量的第二应变量I2,在体积不可压条件下,I2可表示为当λ3值较小时,I2随着λ3的减小而迅速增大,进而导致一个较大的应力响应,这与实验结果是一致的.Neo-Hookean 模型仅依赖于柯西−格林应变张量第一不变量,因此不能够预测当λ3值较小时,随着λ3减小应力响应的迅速增大.

图5 关于硫化橡胶B 的实验和模拟结果对比Fig.5 Comparison between the experimentally measured and simulated results of vulcanized rubber

3 讨论

大量的实验数据显示,当用八链模型来模拟单轴变形实验时,八链模型能很好地预测实验结果,但本文和其他相关论文[8,23]的研究结果都显示当需要同时考虑橡胶的多种加载模式时,八链模型不能较好地描述实验结果,而考虑缠结效应的超弹性模型能克服八链模型的不足之处.

考虑缠结约束效应的超弹性模型的应力响应同时包含熵弹性部分和缠结约束部分,在本文中,前者依赖于柯西−格林应变张量第一不变量,而后者是柯西−格林应变张量第二不变量的函数,因此,这两部分对不同加载情况下的响应是不同的.

这里给出在不同加载条件下两部分的归一化应力,即相应的应力除以剪切模量.如图6 所示,在小应变条件下,来自熵弹性和缠结效应的归一化应力是相同的.当变形逐步增大时,两者之间逐渐偏离.在单轴拉伸和纯剪切变形情况下,缠结效应的归一化应力小于相应的熵弹性归一化应力,且会逐步达到平台并发生软化行为,这与实验结果相一致[27].与之不同,在等双轴拉伸情况下,缠结效应的归一化应力不出现软化行为,且在一定的应变范围,其数值大于熵弹性部分,直至接近材料的拉伸极限时,熵弹性归一化应力急剧增加,其数值将高于缠结约束引起的归一化应力.

经典八链模型能较好地模拟橡胶单轴加载的应力响应,但在模拟橡胶等双轴拉伸时,会低估实验所测得的应力响应,而图6 的结果显示缠结效应在单轴加载条件下对应力的贡献较小,而在等双轴加载条件下,缠结效应具有显著效应,这解释了考虑缠结约束效应的超弹性模型能同时预测单轴、纯剪切和等双轴实验结果的原因.

图6 熵弹性和缠结约束效应归一化应力的对比Fig.6 Comparison between the normalized stress from entropy and entanglement effects

当前,国内外相关学者[24-28]也建立了相应的包含缠结约束效应的超弹性模型,其主要建模思路均是基于管模型,本文将这些模型进行对比分析.Davidson 和Goulbourne[27]采用非仿射的管模型给出了如下的缠结约束自由能

Xiang 等[28]基于三链模型,提出了如下的自由能方程来表示缠结约束效应

为了与本文的模型进行对比,在原模型的自由能基础上增加了系数1/2,以保证上述模型在小应变条件下纯剪切的剪切模量均为Ge.本文相应的缠结约束自由能可表达为

图7 对比了上述模型在单轴、纯剪切和等双轴拉伸情况下的归一化自由能,即相应的自由能除以缠结剪切模量.在单轴拉伸情况下,本文所采用的缠结约束自由能最小,而基于非仿射约束的自由能最大;在纯剪切情况下,3 个自由能的大小比较接近,且基于非仿射假设的自由能和基于三链模型的缠结约束自由能是完全一致的;与单轴拉伸情况不同,在等双轴拉伸情况下,基于非仿射假设的缠结约束自由能最小,而相应的基于三链模型的缠结约束自由能最大.

图7 3 种缠结约束自由能的对比Fig.7 Comparison between three different forms of free energy density of entanglement effects

图7 3 种缠结约束自由能的对比(续)Fig.7 Comparison between three different forms of free energy density of entanglement effects(continued)

在熵弹性模型的基础上,引入上述3 种缠结约束自由能,均能较好地模拟文献中的实验数据[27,28].因此,基于当前文献中的实验数据不足以全面评估哪类缠结约束模型能更好地预测橡胶材料在不同加载情况下的变形行为.亟需开展针对单一橡胶的全面的实验表征,其实验表征可同时考虑单轴拉伸、压缩、平面应变压缩、纯剪切、非等双轴以及等双轴拉伸变形等.这样的实验数据将有益于评估各类模型的适用性.

4 结论

本文建立了包含缠结约束效应的超弹性模型,其主要的建模思路是在微观尺度上,对于单根分子链采用了Langevin 统计模型来表示熵弹性所引起的自由能,采用了管模型来表征缠结约束引起的自由能.采用仿射假设,将各个方向分子链伸长比的均方根作为宏观变形分子链的伸长比,把微球表面微单元收缩比的均方根作为宏观变形管横截面积的收缩比.所建立的模型含有3 个参数:分子链特征刚度、分子链拉伸极限和缠结约束剪切模量.当不考虑缠结效应时,该模型退化为Arruda-Boyce 八链模型.

与经典熵弹性模型相比,本文所建立的模型能更好地表征文献中橡胶在单轴、纯剪切和等双轴多个加载模式下的实验数据,也能准确地预测不同预拉伸比双轴拉伸条件下橡胶的力学响应.上述结果充分揭示了缠结约束效应对橡胶等软材料宏观力学行为的重要影响,在理论建模中需加以考虑.

本文所建立的宏观本构关系仅依赖于柯西−格林应变张量的不变量,因此相比于文献中考虑缠结约束的超弹性模型在有限元植入上更为方便,例如,常见的商业软件ABAQUS 的UHYPER 子程序可通过直接输入基于应变不变量的自由能进行有限元分析,而无需编写相应的基于刚度矩阵的用户自定义子程序,这为软材料的超弹性行为的分析提供了方便[36].目前我们正在开展橡胶材料在多种变形条件下的变形行为和具有复杂几何尺寸试样的力学表征,以更为全面地评估各类超弹性本构模型.