EI、Scopus 收录
中文核心期刊

显式模拟类橡胶材料应力软化引起的不可恢复变形及其各向异性特征

王晓明, 吴荣兴, 蒋义, 肖衡

王晓明, 吴荣兴, 蒋义, 肖衡. 显式模拟类橡胶材料应力软化引起的不可恢复变形及其各向异性特征. 力学学报, 2021, 53(7): 1999-2009. DOI: 10.6052/0459-1879-21-060
引用本文: 王晓明, 吴荣兴, 蒋义, 肖衡. 显式模拟类橡胶材料应力软化引起的不可恢复变形及其各向异性特征. 力学学报, 2021, 53(7): 1999-2009. DOI: 10.6052/0459-1879-21-060
Wang Xiaoming, Wu Rongxing, Jiang Yi, Xiao Heng. Explicitly modeling permanent set and anisotropy property induced by stress softening for rubber-like materials. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(7): 1999-2009. DOI: 10.6052/0459-1879-21-060
Citation: Wang Xiaoming, Wu Rongxing, Jiang Yi, Xiao Heng. Explicitly modeling permanent set and anisotropy property induced by stress softening for rubber-like materials. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(7): 1999-2009. DOI: 10.6052/0459-1879-21-060
王晓明, 吴荣兴, 蒋义, 肖衡. 显式模拟类橡胶材料应力软化引起的不可恢复变形及其各向异性特征. 力学学报, 2021, 53(7): 1999-2009. CSTR: 32045.14.0459-1879-21-060
引用本文: 王晓明, 吴荣兴, 蒋义, 肖衡. 显式模拟类橡胶材料应力软化引起的不可恢复变形及其各向异性特征. 力学学报, 2021, 53(7): 1999-2009. CSTR: 32045.14.0459-1879-21-060
Wang Xiaoming, Wu Rongxing, Jiang Yi, Xiao Heng. Explicitly modeling permanent set and anisotropy property induced by stress softening for rubber-like materials. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(7): 1999-2009. CSTR: 32045.14.0459-1879-21-060
Citation: Wang Xiaoming, Wu Rongxing, Jiang Yi, Xiao Heng. Explicitly modeling permanent set and anisotropy property induced by stress softening for rubber-like materials. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(7): 1999-2009. CSTR: 32045.14.0459-1879-21-060

显式模拟类橡胶材料应力软化引起的不可恢复变形及其各向异性特征

基金项目: 2020访问工程师项目校企合作资助项目(FG2020038)
详细信息
    作者简介:

    王晓明, 讲师, 智能材料本构建模. E-mail: wangxiaoming.g@163.com

  • 中图分类号: O343, O345

EXPLICITLY MODELING PERMANENT SET AND ANISOTROPY PROPERTY INDUCED BY STRESS SOFTENING FOR RUBBER-LIKE MATERIALS

  • 摘要: 类橡胶材料在经过初次加载后会产生应力软化现象, 也就是Mullins效应. 实验证明应力软化现象会导致材料产生不可恢复变形, 同时引入各向异性特征. 本文基于对数应变构造一个多轴可压缩应变能函数, 先引入耗散来表征应力软化现象, 再引入依赖耗散大小的不可恢复变形量以及各向异性特征量, 使得新模型既可以表征Mullins效应, 又能模拟应力软化作用下产生的不可恢复变形和各向异性特征. 本文在各向同性形函数的基础上, 通过球坐标系的思想, 进一步发展并提出了一个任意方向适用的各向异性形函数. 新模型在材料尚未发生软化(耗散为0)的情况下, 表现出各向同性; 一旦发生应力软化(耗散大于0), 则变为各向异性. 随着加载−卸载循环的累积, 耗散逐渐变大, 不可恢复变形也随之变大直到达到一个稳定的值, 各向异性特性也逐渐变得明显. 新方法得到的结果可以精确匹配经典的实验数据, 并预测不同方向的应力软化现象以及由此产生的不可恢复变形和各向异性特征.
    Abstract: Stress softening, known as the Mullins effect, is observed in rubber-like materials after initial loading cycles. Experimental observations have shown that the Mullins effect leads to a permanent set and induces anisotropic properties. A multi-axial, compressible strain energy function based on the Hencky strain is proposed to account for the stress softening and simulate the permanent set and anisotropic properties by introducing two variables, separately characterizing the permanent set and anisotropic properties, which are dependent on dissipation. A comprehensive, explicit shape function is proposed by using the spherical coordinate system to make the model effective in arbitrary direction. The new model exhibits isotropic properties when it not yet induced stress softening, and becomes anisotropy when the Mullins effect occurs. The residual strain increases and reach a fixed value as the loading-unloading cycles proceeds, and the anisotropy becomes more obvious. The simulation results are shown good accord with the classical experimental data and successfully forecast the permanent set and anisotropic properties induced by stress softening.
  • 经过初始加载−卸载后的类橡胶材料再一次加载后, 初始卸载应力以下部分会将发生显著的应力软化现象, 该现象也被称为Mullins效应. 随着加载−卸载循环的进行, 应力软化现象逐渐变得显著, 并趋于一个稳定的状态, 而应力−应变曲线也逐渐变化并发展到一个稳定的状态[1-3]. 天然橡胶和填充橡胶都有可能发生Mullins效应, 而填充物的种类和比重则会影响应力软化的程度.

    过去几十年来, 为了从微观尺度解释Mullins效应, 研究人员们进行了大量的理论研究. Bueche[4-5]认为当材料发生变形后, 填充颗粒之间附着的部分橡胶链发生了断裂和松弛, 使得填充颗粒间的距离发生仿射变形, 从宏观上表现为应力软化现象; Lee和Williams[6]通过结构重组以及化学键永久性断裂来说明Mullins效应的产生原因; Mullins[3]列举了两种情况. 对于天然橡胶, 橡胶长链分子之间在剪切过程中的解链和定向导致了应力软化. 对于填充橡胶, 应力软化的发生则是因为填充颗粒之间以及填充颗粒和橡胶分子之间的相互作用被断开而导致. 填充颗粒的存在使得应力软化效应被进一步增强; 以上观点基本都认为应力软化是因为橡胶分子和填充颗粒之间的相互作用而产生, 而Harwood等[7]提出了不同的见解, 他们认为橡胶微结构中的网络节点位移导致了分子网状结构的重组, 从而产生了应力软化效应. 除此之外, 还有其他的微结构理论[8-10]对Mullins效应的产生做了理论解释.

    在前面微结构理论的基础上, 科研工作者们先构造微观尺度下的本构关系, 再通过一系列平均方法, 建立宏观尺度下的本构模型. Blanchard和Parkinson[11]提出了一个半经验公式去描述填充硫化橡胶初始加载后的“应力−伸长比”关系, 其中引入了二个微结构相关的参数, 分别表示网状结构数量和分子链受限延伸率. 通过分析应力软化的特征, 他们给出了网状结构数量变化的公式, 并最终得到匹配实验数据的“应力−应变”结果; Govindjee和Simo[8]在Bueche[4-5]理论研究的基础上进一步发展, 通过对统计代表样本体积(SRSV)进行宏观平均, 最终得到一个三维的本构模型; Marckmann等[9]发展了Arruda和Boyce[12]的八链理论, 将该模型中的两个参数进行改进, 使他们依赖于分子链伸长比的极限值, 得到了能模拟应力软化效应的新模型.

    另一部分学者不考虑橡胶分子的内部微观结构, 从宏观现象出发构建唯象的本构模型. Mullins和Tobin[13]将橡胶内部区域分为两大类, 分别是“硬区”和“软区”, 材料的变形主要是在“软区”产生, 而“硬区”对变形几乎没有贡献. 加载过程中“硬区”部分逐渐转为“软区”, 卸载后重新加载, 当加载应力达到初始卸载应力以后, “转化”过程才会继续. 此外, 为了模拟填充橡胶的应力−应变关系, Mullins和Tobin[14]又引入了应变放大因子; 文献[15-17]在前面研究的基础上进一步优化和完善; Simo[18]另辟蹊径, 在基于变形梯度的弹性势中加入了一个表征损伤的量, 用来解释类橡胶材料应力软化过程中出现的分子链破坏, 微结构损伤. 在此基础上, 根据引入量演化方程的不同, 提出了各种不同的结果[19-22]; Ogden和Roxburgh[23]在伪弹性模型的基础上只引入了一个变量去表示理想化的应力软化效应.

    Mullins效应通常伴随着不可恢复变形的产生[3](随着时间的流逝会慢慢变小, 在退火后静置足够长时间基本消失), 其大小和橡胶中的填充物含量百分比以及加载应力大小均有关系[24-25]. 通过大量实验证明, 应力软化会使得材料性质变成各向异性, 先前伸长过的方向上应力软化会更加明显[1-3,26]. 各向异性特征也会直接导致材料发生不可恢复变形. 为了模拟Mullins效应产生的不可恢复变形和各向异性特征, Dorfmann和Ogden[27]在伪弹性理论[23]的基础上额外引入了一个表征残余应变的量, 可以模拟单轴拉伸−压缩循环加载下产生的应力软化和不可恢复变形. 他们的模型可以证明材料循环加载下发生了各向异性, 但是无法通过模型进行量化模拟; Göktepe和Miehe[10]将分子内部网络结构分解成CC (crosslink-to-crosslink)结构和PP (particle-to-particle)结构. Mullins效应引起的结构损伤和各向异性主要发生在PP结构中, 且会自然而然地使得模型产生不可恢复变形; Diani等[26]根据Marckmann等[9]的方法, 定义了不同方向上的耗散, 然后将其引入到Pawelski[28]的本构法则中, 构建了新的本构模型. 他们的方法可以预测应力软化、残余变形以及各向异性特征; 谈炳东等[29-30]根据纤维增强复合材料连续介质力学理论, 提出了各向异性超弹性本构模型, 并能很好地预测0° ~ 45°各向异性力学性能.

    关于Mullins效应产生不可恢复变形和各向异性特征的最新研究成果主要是在经典的模型上进行改进[31-34]. 以上提到的模型尽管可以模拟Mullins效应引起的不可恢复变形和各向异性特征, 但是都存在着不足. 微结构模型的参数往往需要迭代求解, 计算量大; 唯象模型则需要引入没有物理意义的隐式参数. 近几年来, Xiao等[35-36]和Wang等[37]用显式的方法提出了类橡胶材料的多轴弹性势, 可以精确匹配4个基准实验, 同时预测非等双轴拉伸的变形趋势. 王晓明等[38]在此基础上引入了耗散量, 并通过分析应力−应变关系, 给出耗散表达, 最终使得模型可以模拟理想的Mullins效应滞回圈. 通过进一步改造形函数的表达, 该方法还能模拟类橡胶材料过载破坏的情况[39]. 在弹性势中加入温度的项, 可以模拟形状记忆聚合物的热力学行为[40].

    本文以前面的研究为基础, 在形函数中新引入表征不可恢复变形大小和各向异性特征的两个变量, 以上两个量均依赖于依赖耗散大小. 其中, 不可恢复变形量随着耗散的增大而增大, 直到趋于稳定; 各向异性特征量在应力软化之前不起任何作用(材料各向同性), 一旦在某一方向发生Mullins效应, 其值将发生变化, 在任意其他方向加载−卸载后的应力−应变关系将和初始加载方向上的应力−应变关系发生明显的偏离(引入各向异性), 在达到相同的变形情况下, 应力值小于初始方向的应力值. 当加载方向和初始加载方向呈90°时, 这样的差异最为明显.

    类橡胶材料本构模型的建立, 通常需要提出合适的弹性势. 本章通过以下4个步骤构造弹性势: 首先, 通过对数应变张量和耗散标量构造弹性势; 其次, 通过对数应变偏量构造3个不变量, 并将它们引入弹性势, 使其多轴可压缩且能匹配多个变形模式; 再次, 将3种变形模式下考虑耗散的单轴形函数积分得到各自的单轴弹性势; 最后, 通过Hermite插值方法, 利用3个不变量和单轴弹性势, 构造考虑耗散的多轴可压缩弹性势.

    理想的超弹性模型没有能量的耗散, 而Mullins效应的存在表明材料在加载−卸载过程中出现了能量的损耗. 为了使模型具备模拟应力软化的功能, 需要引入一个标量$\kappa $来表示耗散. 此外, 相比于其他模式的应变张量, 基于对数应变${\boldsymbol{h}}$的弹性势$W$推导得到应力−应变关系的过程更加简单和直接[41]

    $${\boldsymbol{\tau }} = J{\boldsymbol{\sigma }} = \frac{{\partial W}}{{\partial {\boldsymbol{h}}}}$$ (1)

    式中, $W, J, {\boldsymbol{\tau }}, {\boldsymbol{\sigma }}$分别表示弹性势、体积比、基尔霍夫应力以及柯西应力.

    对数应变与变形梯度的关系为

    $${\boldsymbol{h}} = \frac{1}{2}\ln ({\boldsymbol{F}} \cdot {{\boldsymbol{F}}^{\rm{T}}}) = \frac{1}{2}\sum\limits_{i = 1}^3 {(\ln \lambda _i^2)} {{\boldsymbol{n}}_i} \otimes {{\boldsymbol{n}}_i}$$ (2)

    其中${\boldsymbol{B}} = {\boldsymbol{F}} \cdot {{\boldsymbol{F}}^{\rm{T}}}$, 也称为左柯西−格林张量. 方程(2)中, ${\lambda _i},{{\boldsymbol{n}}_i}$, i = 1,2,3, 是材料伸长比(特征值)和特征向量. 方程(1)中的弹性势$W$是对数应变张量和耗散标量的函数

    $$ W=W({\boldsymbol{h}}, \kappa )$$ (3)

    基于对数应变张量的3个基本不变量表达为

    $${i_s} = \operatorname{tr} {{\boldsymbol{h}}^s} = \sum\limits_{i = {{1}}}^{{3}} {{{(\ln {\lambda _i})}^s}} $$ (4)

    其中s = 1, 2, 3. 对数应变偏量的表达为

    $${\boldsymbol{\tilde h}} = {\boldsymbol{h}} - \frac{1}{3}(\operatorname{tr} {\boldsymbol{h}}){\boldsymbol{I}}$$ (5)

    根据方程(5), 当trh = 0时, i1 = 0, 表示变形前后体积变化忽略不计(材料不可压缩), 使得${\boldsymbol{h}} = {{\tilde {\boldsymbol{h}}}}$. 本文考虑可压缩情况, 所以对数应变偏量和对数应变可以不相等.

    类似于方程(4)的表达, 基于对数应变偏量的3个基本不变量为

    $${j_s} = \operatorname{tr} {{\tilde{\boldsymbol{h}}}^s}$$ (6)

    其中s = 1, 2, 3. 新的3个基本不变量和初始的3个基本不变量之间存在着以下关系

    $$\left. {\begin{array}{*{20}{c}} {{j_1} = 0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \\ {{j_2} = {i_2} - i_1^2{\kern 1pt} {\rm{/3}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \\ {{j_3} = {i_3} - {i_1}{i_2} + {\rm{2}}i_1^3{\rm{/9}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } \end{array}} \right\}$$ (7)

    从以上3个基本不变量出发构造3个新的不变量. 第1个不变量称为压缩不变量, 用${\gamma _{\rm{1}}}$表示, 能控制材料是否可压缩的条件. 当材料不可压缩, 其值为0, 当材料可压缩, 其值不为0; 第2个称为桥联不变量, 用${\gamma _{\rm{2}}}$表示, 能实现模型的多轴扩展. 在单轴情况下会自动退化到相应的对数应变; 第3个不变量称为匹配不变量, 用${\gamma _{\rm{3}}}$表示, 能将3种不同变形模式(单轴拉伸和压缩、等双轴拉伸和压缩、以及平面应变)整合到一个统一的模型当中. 通过分析我们可以得到满足条件的3个不变量表达为

    $$\left. {\begin{array}{*{20}{l}} {{\gamma _{\rm{1}}} = {i_1} = \ln J } \\ { {\gamma _2} = \sqrt {{\rm{2}}{j_2}{\rm{/3}}} } \\ {{\gamma _3} = \sqrt 6 {j_3}{\rm{\Biggr/}}\sqrt {j_2^3} } \end{array}} \right\}$$ (8)

    理想情况下加载−卸载循环作用下的应力−应变曲线都是重合的, 引入了耗散也就意味着每一次循环都会发生能量损耗, 应力−应变曲线随着循环的进行而逐渐变化, 因此其形函数的给定必须考虑耗散.

    单轴拉伸压缩变形模式下, 假设受力方向应力为${\tau _{\rm{u}}}$, 对数应变为$h$, 耗散用$\kappa $表示, 则应力−应变形函数可以表示为

    $$ {\tau }_{\rm{u}}={f}_{\rm{u}}(h, \kappa )$$ (9)

    对方程(9)积分, 就可以得到单轴拉伸情况下的弹性势

    $$ {w}_{\rm{u}}(h, \kappa )={\displaystyle \underset{0}{\overset{h}{\int }}{\tau }_{\rm{u}} {\rm{d}}h=}{\displaystyle \underset{0}{\overset{h}{\int }}{f}_{\rm{u}}(h, \kappa ){\rm{d}}h}$$ (10)

    等双轴拉伸的应力−应变关系可以通过单轴实验的关系推导出来(参考文献[35]中的Theorem A), 再通过积分可以得到其弹性势.

    平面应变是对一个长条形橡胶块(z方向比xy方向尺寸大很多), 在z方向进行固定, x方向和y方向, 一端拉伸, 另一端自由. 假设加载方向应力为${\tau _{\rm{P}}}$, 对数应变为$h$, 固定方向应力为${\tau _{{\rm{pf}}}}$, 我们给定其两个方向的应力−应变关系为

    $$ \left. {\begin{array}{*{20}{l}} {{\tau _{\rm{p}}} = g(h,\kappa )}\\ {{\tau _{\rm{pf}}} = {g_{\rm{f}}}(h,\kappa )} \end{array}} \right\}$$ (11)

    在该变形模式下, 只有在载荷方向做功, 因此我们得到应变能弹性势为

    $$ {w}_{\rm{p}}(h, \kappa )={\displaystyle \underset{0}{\overset{h}{\int }}{\tau }_{\rm{p}}{\rm{d}}h}={\displaystyle \underset{0}{\overset{h}{\int }}g(h, \kappa ){\rm{d}}h}$$ (12)

    方程(3)给定的弹性势依赖于对数应变张量和耗散, 其中前者可以通过1.2节中给定的3个不变量来替换, 即

    $$ W=W({\gamma }_{1},{\gamma }_{2},{\gamma }_{3}, \kappa )$$ (13)

    将方程(13)代入到方程(1), 得到

    $${\boldsymbol{\tau }} = \frac{{\partial W}}{{\partial {\gamma _1}}}{\boldsymbol{I}} + \frac{2}{3}\frac{{\partial W}}{{\partial {\gamma _2}}}\gamma _2^{ - 1}{\tilde{\boldsymbol{h}}} + \frac{{\partial W}}{{\partial {\gamma _3}}}{\mathop {\boldsymbol{h}} \limits^ \smallsmile}$$ (14)

    其中

    $${\mathop {\boldsymbol{h}} \limits^ \smallsmile} = 4\gamma _2^{ - 3}{{\tilde{\boldsymbol{h}}}^2} - 2\gamma _2^{ - 2}{\gamma _3}{\tilde{\boldsymbol{h}}} - 2\gamma _2^{ - 1}{\boldsymbol{I}}$$ (15)

    式中${\tilde{\boldsymbol{h}}}$是对数应变偏量, ${\boldsymbol{I}}$是单位二阶张量. 张量${\boldsymbol{I}},{\tilde{\boldsymbol{h}}},{\mathop {\boldsymbol{h}} \limits^ \smallsmile}$之间相互正交

    $$\left. {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\boldsymbol{I}}:{\tilde{\boldsymbol{h}}} = {\bf{0}}} \\ {{\boldsymbol{I}}:{\mathop {\boldsymbol{h}} \limits^ \smallsmile} = {\bf{0}}} \end{array}} \\ {{\tilde{\boldsymbol{h}}}:{\mathop {\boldsymbol{h}} \limits^ \smallsmile} = {\bf{0}}} \end{array}} \right\}$$ (16)

    可以作为对数应变张量的3个基.

    接下来利用Hermite插值方法, 我们可以得到多轴可压缩弹性势为

    $$W = \frac{{1 - 2\upsilon }}{3}{w_{\rm{u}}}\left( {\frac{{{\gamma _1}}}{{1 - 2\upsilon }},\kappa } \right) + \frac{{1 + \upsilon }}{6}\left[{{{Z}}^ + }{(1 + {\gamma _3})^2} + {{{Z}}^{\rm{ - }}}{(1 - {\gamma _3})^2}\right]$$ (17)

    其中Z +Z 表达为

    $$\left. {\begin{array}{*{20}{c}} {{{{Z}}^ + } = (2 - {\gamma _3}){w_{\rm{u}}}\left( {\dfrac{{1.5{\gamma _2}}}{{1 + \upsilon }},\kappa } \right) + ({\gamma _3} - 1){Y_ + }}\\ {{{{Z}}^{\rm{ - }}} = (2 + {\gamma _3}){w_{\rm{u}}}\left( { - \dfrac{{1.5{\gamma _2}}}{{1 + \upsilon }},\kappa } \right) + ({\gamma _3} + 1){Y_ - }} \end{array}} \right\}$$ (18)

    ${Y_ + },{Y_ - }$分别等于${\left. {\dfrac{{\partial W}}{{\partial {\gamma _3}}}} \right|_{{\gamma _{\rm{3}}}{\rm{ = 1}}}}$${\left. {\dfrac{{\partial W}}{{\partial {\gamma _3}}}} \right|_{{\gamma _{\rm{3}}}{\rm{ = - 1}}}}$, 形式为

    $$ \left. {\begin{array}{*{20}{l}} {\begin{array}{*{20}{l}} {{Y_ + } = \dfrac{5}{2}{w_{\rm{u}}}\left( {\dfrac{{1.5{\gamma _2}}}{{1 + \upsilon }},\kappa } \right) - \dfrac{1}{2}{w_{\rm{u}}}\left( { - \dfrac{{1.5{\gamma _2}}}{{1 + \upsilon }},\kappa } \right) - }\\ {\;\;\;\;2{w_{\rm{p}}}\left( {\dfrac{{3\sqrt 3 {\gamma _2}}}{{4(1 + \upsilon )}},\kappa } \right) - \dfrac{{\sqrt 3 {\gamma _2}}}{{4(1 + \upsilon )}}\left[ {g\left( {\dfrac{{3\sqrt 3 {\gamma _2}}}{{4(1 + \upsilon )}},\kappa } \right) - } \right.}\\ {\left. {\;\;\;\;2{g_{\rm{f}}}\left( {\dfrac{{3\sqrt 3 {\gamma _2}}}{{4(1 + \upsilon )}},\kappa } \right)} \right]} \end{array}}\\ {\begin{array}{*{20}{l}} {{Y_{\rm{ - }}} = \dfrac{{\rm{1}}}{2}{w_{\rm{u}}}\left( {\dfrac{{1.5{\gamma _2}}}{{1 + \upsilon }},\kappa } \right) - \dfrac{{\rm{5}}}{2}{w_{\rm{u}}}\left( { - \dfrac{{1.5{\gamma _2}}}{{1 + \upsilon }},\kappa } \right) + }\\ {\;\;\;\;2{w_{\rm{p}}}\left( {\dfrac{{3\sqrt 3 {\gamma _2}}}{{4(1 + \upsilon )}},\kappa } \right) - \dfrac{{\sqrt 3 {\gamma _2}}}{{4(1 + \upsilon )}}\left[ {g\left( {\dfrac{{3\sqrt 3 {\gamma _2}}}{{4(1 + \upsilon )}},\kappa } \right) - } \right.}\\ {\left. {\;\;\;\;2{g_{\rm{f}}}\left( {\dfrac{{3\sqrt 3 {\gamma _2}}}{{4(1 + \upsilon )}},\kappa } \right)} \right]} \end{array}} \end{array}} \right\}$$ (19)

    式中$\upsilon $表示泊松比, 等于侧向应变${h_{\rm{l}}}$比上受力方向应变$h$的相反数

    $$\upsilon {\rm{ = - }}\frac{{{h_{\rm{l}}}}}{h}$$ (20)

    当材料变形前后体积不变, ${h_l}{\rm{ = - 0}}{\rm{.5}}\;h$, 从式(20)可以得到泊松比$\upsilon $=0.5.

    从式(17)出发, 可以得到单轴实验下的应力应变关系为方程(9), 平面应变下的应力−应变关系为方程(11), 而等双轴下的应力−应变关系[38]

    $${\tau _{\rm{e}}} = - \frac{1}{{2\upsilon }}f\left( { - \frac{1}{\upsilon }h,\kappa } \right)$$ (21)

    本节首先通过应变−应力图中面积所代表的能量关系推导出耗散$\kappa $的显式表达; 其次, 随着耗散的增加, 分析不可恢复变形大小的演变规律以及各向异性引入对应力−应变关系的影响, 分别通过${h_P}(\kappa )$$ \varphi (\kappa )$来表示不可恢复变形大小以及各向异性对应力−应变关系的影响程度; 最后, 利用耗散$\kappa $, 不可恢复变形大小${h_P}(\kappa )$, 以及各向异性影响项$ \varphi (\kappa )$, 结合球坐标, 构造新的形函数表达形式.

    图1所示, 橡胶材料在初始加载下, 应力−应变沿着(a)前进, 在达到P1点后开始卸载, 卸载路径沿着(b)运行, 最后产生不可恢复变形hP1, 完成一个循环, 产生的能量耗散可以用图1中(1)的面积来表示; 再次加载后, 初始卸载应力点以下的部分发生明显的应力软化, 而之后这样的软化就不是很明显[1-3]. 因此, 第2次加载路径应该是从(b)变成(c), 到达P2点后再次卸载, 卸载路径沿着(d)走(此时软化现象更加严重, 卸载曲线和初始加载曲线之间的偏离更加严重), 最后产生不可恢复变形hP2, 产生的能量耗散可以用图中(2)的面积来表示; 第3次循环的路径可以通过类似的规律推导.

    图  1  Mullins效应产生不可恢复变形示意图
    Figure  1.  Schematic of Mullins effects with permanent set

    通过大量应力软化的实验数据分析, 循环加载−卸载作用下能量耗散的变化规律有以下3点: (1)随着循环加载−卸载的进行, 能量耗散逐渐累积; (2)从单个加载−卸载循环来看, 能量耗散最小的情况是加载和卸载路径重合(不产生应力软化), $\kappa {\rm{ = 0}}$; (3)从单个加载−卸载循环来看, 能量损耗最大的情况是卸载曲线极度弯曲(应力软化最大化), $\kappa = {\kappa _{\rm{m}}}$. 其中

    $${\kappa _{\rm{m}}} = \int\limits_0^{{h_{\rm{m}}}} {\bar f} (h){\rm{d}}h$$ (22)

    式中, $\bar f(h)$表示单个循环中加载阶段的应力−应变曲线, ${h_{\rm{m}}}$表示卸载点应变.

    滞回圈面积大小无法超过${\kappa _{\rm{m}}}$, 因此

    $${\rm{0}} \leqslant \kappa {\rm{/}}{\kappa _{\rm{m}}} \leqslant {\rm{1}}$$ (23)

    耗散随着卸载应力${\tau _{\rm{m}}}$的变化规律可以通过双曲正切函数得到

    $$\kappa = \frac{{{\kappa _{\rm{m}}}}}{2}\Big\{\tanh [{\alpha _1}({\tau _{\rm{m}}} - {\tau _{\rm{r}}})]{\rm{ + 1}}\Big\}$$ (24)

    式中, ${\tau _{\rm{r}}}$${\alpha _1}$是可调参数. 当${\tau _{\rm{m}}} = {\tau _{\rm{r}}}$时, $\kappa = \dfrac{{{\kappa _m}}}{2}$; 参数$\alpha $表示$\kappa $随着${\tau _{\rm{m}}}$变化的速度.

    耗散的引入会产生不可恢复变形以及各向异性特征. 本节讨论不可恢复变形以及各向异性特征随着耗散累积的变化规律, 分别通过${h_P}(\kappa )$$ \phi (\kappa )$来表征不可恢复变形大小和各向异性对应力−应变关系的影响程度.

    在保持最大应变的情况下, 不可恢复变形的变化规律如下[27]: (1)不可恢复变形量随着循环的进行会逐渐变大; (2)第一个循环前后产生的变形最大, 后面依次减少; (3)在循环次数足够多的情况下, 不可恢复变形不再变化, 达到一个稳定值. 假设有n次循环, hPi表示每一次产生的不可恢复变形, i = 1, 2, ···, n. 则有

    $${h_{{{Pn}}}} > {h_{{{P}}(n - 1)}} > \cdots > {h_{{{P2}}}} > {h_{{{P1}}}}\qquad\qquad\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$$ (25)
    $${h_{{{P1}}}} > ({h_{{{P2}}}} - {h_{{{P1}}}}) > ({h_{{{P3}}}} - {h_{{{P2}}}}) > \cdots > ({h_{{{Pn}}}} - {h_{{{P(n - 1)}}}})$$ (26)
    $$\mathop {\lim }\limits_{{{n}} \to \infty } {h_{{{Pn}}}} = {h_{{{P(n - 1)}}}} = {h_{{P}}}\qquad\qquad\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\qquad$$ (27)

    耗散的引入导致了不可恢复变形的产生. 因此, 假设

    $${h_{{P}}} = {h_{{P}}}(\kappa ) = {{{p}}_1}\tan({{{p}}_2}\kappa )$$ (28)

    ${h_{{P}}}$随着耗散的增大而增大, 直到趋于一个稳定的数值. ${p_1}$${p_2}$是可调参数.

    图2所示, 假设六面体橡胶块初始情况下为各向同性. 从任意方向加载, 其性质都是一样的. 一旦在某一方向(比如x轴方向)加载然后卸载, Mullins效应的产生将导致材料引入了各向异性特征, 此时, 再沿着初始方向(x轴方向)加载−卸载的应力−应变关系将和沿着其他方向(比如y轴或者z轴)的发生明显的差异.

    图  2  六面体橡胶块三维受力示意图
    Figure  2.  Schematic of hexahedral rubber block with three-dimensional force

    为了证明并量化Mullins效应引入的各向异性性质. Diani等[26]做了相关实验. 首先, 将A和B两块相同的橡胶六面体同时在x轴方向进行同样的拉伸加载, 然后卸载. 此时两块试样均发生了Mullins效应且都产生了相同的不可恢复变形. 接下来, 对A进行x轴方向的拉伸加载和卸载, 对B进行y轴方向(垂直于初始方向)的拉伸加载和卸载. 将两块试样的第2个循环的应力−应变曲线(Mullins效应已经饱和)进行比较, 结果如图3所示, 其中, 应变类型采用对数应变. 从图3可以得到二个结论: (1) A和B两块试样的应力−应变曲线不重合且有明显的差异, 说明Mullins效应确实引入了各向异性特征; (2) A的曲线明显高于B曲线(同样应变情况下, A的应力更大), 说明B的应力软化更加严重. 假设在x轴上的应力−应变形函数为$f(h,\kappa )$, 那么与其垂直方向加载−卸载的形函数为$\tilde f(h,\kappa ) = \varphi (\kappa )f(h,\kappa )$, 其中$\varphi (\kappa )$表征了各向异性程度, 其满足以下条件

    $$\left. {\begin{array}{*{20}{c}} {\mathop {{\rm{lim}}}\limits_{\kappa \to 0} \varphi (\kappa ){\rm{ = 1}}} \\ {\mathop {{\rm{lim}}}\limits_{\kappa \to \infty } \varphi (\kappa ){\rm{ = 0}}} \end{array}} \right\}$$ (29)

    式(29)表示如果没有发生耗散($ \kappa =0$), 那么垂直方向的应力−应变关系和初始方向一致, 材料没有发生各向异性特征; 如果耗散很大, 那么材料发生极度的各向异性, 垂直方向几乎没有抵抗变形的能力. 基于以上条件, 可以假设$ \varphi (\kappa )$的形式为

    $$\varphi (\kappa ){\rm{ = - }}\frac{1}{2}{\rm{tanh[}}{\alpha _2}(\kappa - {\kappa _{\rm{r}}}){\rm{] + }}\frac{1}{2}$$ (30)

    式中, ${\alpha _2}$${\kappa _{\rm{r}}}$是可调参数. 在初始方向和垂直方向中间的形函数, 其各向异性特征项的值理论上介于0 ~ 1之间.

    图  3  Diani等[26]证明引入各向异性的实验数据
    Figure  3.  An experimental evidence of induced anisotropy by Diani et al.[26]

    首先给出各向同性的形函数形式; 然后利用球坐标转换和2.2.2节中给定的各向异性项$ \varphi (\kappa )$进行扩展, 得到能模拟各向异性特征的形函数; 最后, 在单轴情况下, 利用有理插值的方法给定带有不可恢复变形${h_P}(\kappa )$的形函数表达. 将结果代入到前面的各向异性形函数, 从而构造统一的形函数形式.

    形函数分为加载曲线和卸载曲线, 如图1所示. 初始加载曲线沿着路径(a)前进, P1点卸载以后发生应力软化(耗散的产生). 重新加载直到上一次卸载应力之前的部分将沿着上一次卸载曲线路径(b)行走. 一旦应力超过P1点应力, 应力软化降低, 应力−应变曲线沿着(c)前进. 再次在P2点卸载后重新加载, 将沿着从(d)到(e)的方向前进; 卸载曲线的变化和耗散紧密相关, 耗散越大, 软化越严重, 和同一滞回圈内的加载曲线偏离就越大. 满足以上条件的形函数构造形式如下

    $$\begin{split} \widehat{f}(h,\kappa )=&\frac{1}{4}(1+\chi )(1+\pi )\overline{f}(h)+\\ &[(1+\chi )(1-{\pi})+2(1-\chi )]\tilde{f}(h,\kappa ) \end{split}$$ (31)

    其中

    $$\left. {\begin{array}{*{20}{l}} {\tilde f(h,0) = \bar f(h)} \\ {\pi = \tanh [{\alpha _3}(\tau - {\tau _{\rm{m}}})]} \\ {\chi = \dot h/\left| {\dot h} \right|} \end{array}} \right\}$$ (32)

    式中${\tau _{\rm{m}}}$表示上一次卸载应力, ${\alpha _3}$是可调参数, $\dot h$表示应变率. 参数$\chi $$\pi$具有如下性质

    $$\chi = \left\{ {\begin{array}{*{20}{c}} {{\rm{1}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{loading}}} \\ {{\rm{ - 1}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{unloading}}} \end{array}} \right.$$ (33)
    $$\pi = \left\{ {\begin{array}{*{20}{c}} {{\rm{1}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \tau > {\tau _{\rm{m}}}} \\ {{\rm{ - 1}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \tau < {\tau _{\rm{m}}}} \end{array}} \right.$$ (34)

    下面分析加载曲线和卸载曲线. 初次加载时, 将$\kappa = 0$, $\chi = 1$, $\pi = 1$代入到方程(31), 得到加载曲线方程为

    $$\hat f(h,\kappa ) = \bar f(h)$$ (35)

    有了加载历史以后, $\chi = 1$不变, $\kappa {\rm{ > }}0$, 加载曲线方程变为

    $$\hat f(h,\kappa ) = \left\{ {\begin{array}{*{20}{c}} {\tilde f(h,\kappa ),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \tau < {\tau _{\rm{m}}}} \\ {\bar f(h),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \tau > {\tau _{\rm{m}}}} \end{array}} \right.$$ (36)

    然后分析卸载曲线, 此时$\chi {\rm{ = - }}1$, $\kappa {\rm{ > }}0$, 得到曲线方程为

    $$\hat f(h,\kappa ) = \tilde f(h,\kappa ){\kern 1pt} $$ (37)

    方程(31)得到的各向同性形函数还需要进一步改进, 使其能够产生各向异性性质. 改进形式为

    $$f(h,\kappa ) = \hat f(h,\kappa ){\kern 1pt} ({\cos ^{\rm{2}}}\phi {\rm{ + }}{\kern 1pt} \varphi (\kappa ){\kern 1pt} {\rm{si}}{{\rm{n}}^{\rm{2}}}\phi {\rm{si}}{{\rm{n}}^{\rm{2}}}\theta {\rm{ + }}{\kern 1pt} \varphi '(\kappa ){\rm{si}}{{\rm{n}}^{\rm{2}}}\phi {\rm{co}}{{\rm{s}}^{\rm{2}}}\theta )$$ (38)

    图4所示, r表示受力方向, $\phi $表示受力方向在x轴和y轴组成的平面投影和x轴的夹角, $\theta $表示受力方向和z轴的夹角. 假设初始受力方向为x轴, $\varphi (\kappa )$表示与其垂直的y轴方向各向异性程度, 而$\varphi '(\kappa )$表示与其垂直的z轴方向各向异性程度, 他们均可用方程(30)的形式给出.

    图  4  任意方向受力的球坐标示意图
    Figure  4.  Schematic of loading in arbitrary direction by Spherical coordinate system

    在初始x轴方向加载−卸载的基础上, 如果继续在同样的方向进行加载, 那么r方向和x轴重合, 此时$\phi = {{\rm{0}}^ \circ }$, $\theta {\rm{ = 9}}{{\rm{0}}^ \circ }$, 得到

    $$f(h,\kappa ) = \hat f(h,\kappa )$$ (39)

    如果再次加载方向为y轴, 则$\phi {\rm{ = 9}}{{\rm{0}}^ \circ }$, $\theta {\rm{ = 9}}{{\rm{0}}^ \circ }$, 得到

    $$f(h,\kappa ) = \hat f(h,\kappa ){\kern 1pt} \varphi (\kappa )$$ (40)

    如果再次加载方向为z轴, 则$\phi {\rm{ = 9}}{{\rm{0}}^ \circ }$, $\theta = {{\rm{0}}^ \circ }$, 得到

    $$f(h,\kappa ) = \hat f(h,\kappa )\varphi '(\kappa )$$ (41)

    本节需要给出方程(31)中函数$\bar f(\kappa )$$\tilde f(h{\rm{,}}\kappa )$具体形式, 然后代入到方程(38), 得到统一的形函数表达. $\bar f(\kappa )$$\tilde f(h{\rm{,}}\kappa )$表示的是单轴情况下基准实验形函数, 可以通过有理插值的方法给出. 对于单轴拉伸−压缩变形模式下的形函数给定为

    $$\bar f(h) = {\bar f_{\rm{u}}}(h) = {E_0}h\left[ {\dfrac{{{\alpha _0}}}{{\left( {1 - \dfrac{h}{{{h_{{\rm{10}}}}}}} \right)\left( {1 + \dfrac{h}{{{h_{{\rm{20}}}}}}} \right)}} + 1 - {\alpha _0}} \right]$$ (42)
    $$\begin{split} \tilde f(h,\kappa ) =& {\tilde f_{\rm{u}}}(h,\kappa ) = E(\kappa )(h - {h_{\rm{P}}}(\kappa ))\\ &\left\{ {\dfrac{{\alpha (\kappa )}}{{\left[ {1 - \dfrac{{h - {h_{\rm{P}}}(\kappa )}}{{{h_{\rm{1}}}(\kappa )}}} \right]\left[ {1 + \dfrac{{h - {h_{\rm{P}}}(\kappa )}}{{{h_{\rm{2}}}(\kappa )}}} \right]}} + 1 - \alpha (\kappa )} \right\} \end{split}$$ (43)

    关于以上两个函数具备的性质在前面的研究工作[38]中已经阐明. 其中方程(43)中的$E(\kappa )$, ${h_1}(\kappa )$, ${h_2}(\kappa )$$\alpha (\kappa )$是依赖耗散的函数, 当耗散$\kappa {\rm{ = 0}}$, 他们的初始值即为方程(42)中的${E_0}$, ${h_{10}}$, ${h_{20}}$${\alpha _0}$. 和前面研究不一样的地方在于本文的方程中加入了依赖耗散的不可恢复变形${h_{{P}}}(\kappa )$; 等双轴实验的形函数结合方程(21)和(38)得到平面应变情况下的形函数形式, 其加载方向为

    $$\left. \begin{array}{l} f(h) = \bar g(h) = \dfrac{2}{3}{E_0}h\left( {\dfrac{{{\alpha _{{\rm{q0}}}}}}{{ {1 - \dfrac{{{h^2}}}{{h_{{\rm{q0}}}^2}}} }} + 1 - {\alpha _{{\rm{q0}}}}} \right)\\ \tilde f(h,\kappa ) = \tilde g(h,\kappa ) = \dfrac{2}{3}E(\kappa )[h - {h_P}(\kappa )]\cdot \\\qquad \left\{ {\dfrac{{{\alpha _{\rm{q}}}(\kappa )}}{{ {1 - \dfrac{{{{[h - {h_P}(\kappa )]}^2}}}{{h_{\rm{q}}^2(\kappa )}}} }} + 1 - {\alpha _{\rm{q}}}(\kappa )} \right\} \end{array} \right\}$$ (44)

    固定方向为

    $$\left. {\begin{array}{*{20}{c}} {\bar f(h) = {{\bar g}_{\rm{f}}}(h) = \dfrac{{\rm{1}}}{3}{E_0}h\left( {\dfrac{{{\alpha _{{\rm{qf0}}}}}}{{ {1 - \dfrac{{{h^2}}}{{h_{{\rm{q}}0}^2}}} }} + 1 - {\alpha _{{\rm{qf0}}}}} \right)}\\ \tilde f(h,\kappa ) = {{\tilde g}_{\rm{f}}}(h,\kappa ) = \dfrac{{\rm{1}}}{3}E(\kappa )[h - {h_P}(\kappa )]\cdot\\ \qquad \left\{ {\dfrac{{{\alpha _{{\rm{qf}}}}(\kappa )}}{{{1 - \dfrac{{{{[h - {h_P}(\kappa )]}^2}}}{{h_{\rm{q}}^2(\kappa )}}} }} + 1 - {\alpha _{{\rm{qf}}}}(\kappa )} \right\} \end{array}} \right\}$$ (45)

    式中${h_{\rm{q}}}(\kappa )$, ${\alpha _{\rm{q}}}(\kappa )$${\alpha _{\rm{qf}}}(\kappa )$是依赖耗散$\kappa $的函数, 他们的初始值分别为${h_{{\rm{q0}}}}$, ${\alpha _{{\rm{q0}}}}$${\alpha _{{\rm{qf0}}}}$.

    第2节构造的形函数代入到第1节的本构方程中, 能够自动得到对应变形模式的应力−应变关系. 3个基准实验和Mullins效应产生的理想滞回圈数据都能精确匹配[38]. 本文重新优化了形函数, 加入了考虑不可恢复变形和各向异性的项, 本章给出模型结果和经典实验数据的对比.

    单轴情况下对数应变$h = \ln \lambda $, 代入到方程(42)和(43), 再结合方程(31)和(38), 即可得到单轴拉伸−压缩情况下的应力−应变形函数.

    Diani等[26]做了单轴拉伸加载−卸载循环实验. 分别保持伸长比$\lambda $为1.5, 2, 2.5和3不变, 进行10次加载−卸载, 产生了相应的4个稳定的不可恢复变形.

    为了得到模型结果, 给出方程(43)中的$E(\kappa )$ (MPa), ${h_1}(\kappa )$, ${h_2}(\kappa )$$\alpha (\kappa )$的具体表达为

    $$\begin{split} E(\kappa ) = {\rm{9}}\exp ({\rm{ - 0}}{\rm{.5}}\kappa ) + 1 \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\\[-10pt] \end{split}$$ (46)
    $$\begin{split} {h_1}(\kappa ) =& \{{\rm{exp}}( - 8\kappa ) + 0.25 + 0.45\{\tanh [0.5(\kappa - 1.6)]{\rm{ + }}1\}\}\cdot {\kern 1pt} \\ &[ - 0.06\tanh (10\kappa - 35) + 0.9] \\[-10pt]\end{split} \;\;\;\;\;$$ (47)
    $$\begin{split} {h_2}(\kappa ){\rm{ = 10}}\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\\[-10pt] \end{split}$$ (48)
    $$\alpha (\kappa ) = {\rm{0}}{\rm{.1 + 0}}{\rm{.06}}\exp ( - 10\kappa ) + 0.13\{{\rm{tanh}}[1.8(\kappa - 2.6)] + 1\}$$ (49)

    通过方程(46) ~ 式(49)得到, 当$\kappa {\rm{ = 0}}$时初始值(方程(42)中的参数)分别为${E_0}{\rm{ = 10\;MPa}}$, ${h_{10}}{\rm{ = 1}}{\rm{.4}}$, ${h_{20}}{\rm{ = }} $$ 10$${\alpha _0}{\rm{ = 0}}{\rm{.16}}$.

    方程(28)中, 取${p_1} = 0.25$, ${p_2} = 0.27$. 方程(24)中控制耗散变化规律的两个参数${\tau _{\rm{r}}}$${\alpha _1}$分别给定为10.15 MPa和0.056. 实验中的4次卸载应力分别为4.49 MPa, 7.68 MPa, 11.23 MPa和15.98 MPa. 再根据方程(42)确定${\kappa _{\rm{m}}}$的大小分别为1.01, 2.51, 4.68和7.05, 由此可以确定耗散功大小. 将耗散代入到方程(46) ~ 式(49)确定每一次循环参数值, 方程(32)中的${\alpha _3} = 100$. 最后根据方程(31)计算每一次循环的加载曲线和卸载曲线. 与实验结果的对比如图5所示.

    图  5  模型结果和实验结果[26]对比图. 横坐标表示对数应变, 纵坐标表示基尔霍夫应力
    Figure  5.  Comparison between model result and the experimental data[26]. “x” axis represents the Hencky strain, and “y” axis represents the Kirchhoff stress

    实线表示模型模拟的结果, 其中黑线表示第1次加载路径, 红色表示第1次卸载然后再次加载的路径, 蓝色表示第2次卸载然后再次加载的路径, 洋红色表示第3次卸载然后再次加载的路径, 褐色表示第4次卸载然后再次加载的路径, 紫色表示第5次卸载然后再次加载的路径. 空心圆点表示实验数据. 从图中可以看出模型可以精确模拟和匹配前四次加载−卸载的实验数据, 而对于第5次卸载(假设卸载应力为21.08 MPa, 塑性功为9.17)的曲线, 能够进行合理的预测.

    为了证明Mullins效应会引入各向异性特征, Diani等[26]做了相关实验, 结果如图3所示. 假设方程(38)中的$\varphi \left( \kappa \right)$$\varphi '\left( \kappa \right)$的形式一致, 都通过方程(30)来确定, 其中参数${\alpha _2} = 0.{\rm{52}}$, ${\kappa _{\rm{r}}}{\rm{ = 3}}{\rm{.11}}$. 初次加载到$\lambda {\rm{ = 2}}$, 对数应变为$\ln \lambda {\rm{ = 0}}{\rm{.6931}}$, 通过前面的方法可得耗散为2.51. 代入到方程(30)得到$\varphi \left( \kappa \right) = 0.{\rm{6}}5$, $\varphi '\left( \kappa \right) = 0.{\rm{6}}5$. 假设初始加载方向为x轴, 后面再次加载分别在x轴和y轴, 通过方程(38)可以得到模型结果, 和实验结果的对比如图6所示.

    图  6  Diani等[26]的实验数据和模型结果对比
    Figure  6.  Comparison between experimental data of Diani et al.[26] and model result

    黑色圆点和空心圆点分别表示x方向和y方向的实验数据, 黑色实线和红色实线分别表示x方向和y方向的模型结果. 从图6可以看出, 模型结果和实验结果可以精确匹配.

    除了匹配实验结果, 模型还能预测其他方向受力的应力−应变滞回圈关系. 在$\varphi \left( \kappa \right)$$\varphi '\left( \kappa \right)$形式一致的前提下, 方程(38)转化为

    $$f(h,\kappa ) = \hat f(h,\kappa ){\kern 1pt} [{\cos ^{\rm{2}}}\phi {\rm{ + }}{\kern 1pt} \varphi (\kappa ){\kern 1pt} {\rm{si}}{{\rm{n}}^{\rm{2}}}\phi ]$$ (51)

    $ \phi ={0}^{\circ }, {30}^{\circ }, {60}^{\circ }$${90}^{\circ } $的情况下, 结果对比如图7所示. 其中$ \phi ={0}^{\circ }$${90}^{\circ } $分别表示x方向和y方向, 图6结果证明这两个方向受力得到的模型结果和实验数据可以精确匹配. 而$ \phi ={30}^{\circ }$${60}^{\circ } $方向上的应力−应变关系可以通过模型进行合理的预测.

    图  7  不同方向上模型结果对比
    Figure  7.  Comparison of model results in different directions

    本文在前面研究[38]的基础上, 改进了形函数的构造形式, 其中加入了依赖耗散的不可恢复变形项和控制各向异性的项. 通过图5 ~ 图7的结果, 模型结果不仅可以精确匹配实验数据, 同时可以对结果做合理的预测. 从而证明新的模型可以模拟Mullins效应产生的不可恢复变形以及由此引入各向异性特征.

    本文创新点在于以下3点: (1)改进了加载−卸载形函数的统一模式(方程(31)), 引入了新的量${\text{π}} $, 实现了加载曲线从应力显著软化($\tau < {\tau _{\rm{m}}}$)到不显著软化($\tau > {\tau _{\rm{m}}}$)的平滑过度; (2)形函数中引入了不可恢复变形部分${h_{\rm{P}}}$, 并给出其具体形式, 且符合2.2节给出的3个条件, 最后能精确匹配和预测实验结果(图5); (3)通过球坐标构建了任意方向受力的形函数改进形式(方程(38)), 并在其中加入了控制各向异性的项$\varphi \left( \kappa \right)$$\varphi '\left( \kappa \right)$及其具体表达(方程(30)). 引入的各向异性满足2.2.2节中提出的条件且能精确匹配和预测实验结果(图6).

  • 图  1   Mullins效应产生不可恢复变形示意图

    Figure  1.   Schematic of Mullins effects with permanent set

    图  2   六面体橡胶块三维受力示意图

    Figure  2.   Schematic of hexahedral rubber block with three-dimensional force

    图  3   Diani等[26]证明引入各向异性的实验数据

    Figure  3.   An experimental evidence of induced anisotropy by Diani et al.[26]

    图  4   任意方向受力的球坐标示意图

    Figure  4.   Schematic of loading in arbitrary direction by Spherical coordinate system

    图  5   模型结果和实验结果[26]对比图. 横坐标表示对数应变, 纵坐标表示基尔霍夫应力

    Figure  5.   Comparison between model result and the experimental data[26]. “x” axis represents the Hencky strain, and “y” axis represents the Kirchhoff stress

    图  6   Diani等[26]的实验数据和模型结果对比

    Figure  6.   Comparison between experimental data of Diani et al.[26] and model result

    图  7   不同方向上模型结果对比

    Figure  7.   Comparison of model results in different directions

  • [1]

    Mullins L. Softening of rubber by deformation. Rubber Chemistry and Technology, 1969, 42(1): 339-362 doi: 10.5254/1.3539210

    [2]

    Mullins L. Effect of stretching on the properties of rubber. Rubber Chemistry and Technology, 1948, 21(2): 281-300 doi: 10.5254/1.3546914

    [3]

    Mullins L. Permanent set in vulcanized rubber. Rubber Chemistry and Technology, 1949, 22(4): 1036-1044 doi: 10.5254/1.3543010

    [4]

    Bueche F. Molecular basis of the Mullins effect. Journal of Appli-ed Polymer Science, 1960, 4: 107-114 doi: 10.1002/app.1960.070041017

    [5]

    Bueche F. Mullins effect and rubber-filler interaction. Journal of Applied Polymer Science, 1961, 5: 271-281 doi: 10.1002/app.1961.070051504

    [6]

    Lee M, William M. Application of a new structural theory to pol-ymers. A. Uniaxial stress in crosslinked rubbers. Journal of the Mechanics and Physics of Solids, 1985, 45: 1805-1834

    [7]

    Harwood JAC, Mullins L, Payne AR. Stress softening in natural rubber vulcanizates. Part II. Stress softening effects in pure gum and filler load rubbers. Journal of Applied Polymer Science, 1965, 9: 3011-3021.

    [8]

    Govindjee S, Simo JC. A mirco-mechanically based continuum damage model for carbon black-filled rubbers incorporating Mullin-s’ effect. Journal of the Mechanics and Physics of Solids, 1991, 39(1): 87-112 doi: 10.1016/0022-5096(91)90032-J

    [9]

    Marckmann G, Verron E, Gornet L, et al. A theory of network alteration for the Mullins effect. Journal of the Mechanics and Physics of Solids, 2002, 50: 2011-2028 doi: 10.1016/S0022-5096(01)00136-3

    [10]

    Göktepe S, Miehee C. A micro-macro approach to rubber-like mat-erials. Part III: The micro-sphere model of anisotropic Mullins-typedamage. Journal of the Mechanics and Physics of Solids, 2005, 53: 2259-2283 doi: 10.1016/j.jmps.2005.04.010

    [11]

    Blanchard AF, Parkinson D. Breakage of carbon-rubber networks by applied stress. Rubber Chemistry and Technology, 1952, 44(4): 799-812

    [12]

    Arrude EM, Boyce MC. A three-dimensional constitutive model forthe large stretch behavior of rubber elastic materials. Journal of the Mechanics and Physics of Solids, 1993, 41(2): 389-412 doi: 10.1016/0022-5096(93)90013-6

    [13]

    Mullins L, Tobin, NR. Theoretical model for the elastic behavior of filler reinforced vulcanized rubbers. Rubber Chemistry and Technology, 1957, 30(2): 555-571 doi: 10.5254/1.3542705

    [14]

    Mullins L, Tobin NR. Stress softening rubber vulcanizates Part I. Use of a strain amplification factor to describe the elastic behavior of filler-reinforced vulcanized rubber. Journal of Applied polymer Science, 1965, 9: 2993-3009 doi: 10.1002/app.1965.070090906

    [15]

    Klüppel M, Schramm J. An advanced micro-mechanical model of hyper-elasticity and stress softening of reinforced rubbers // Dorf-mann, Eds. Constitutive Models for Rubber. 1999, 211-218

    [16]

    Govindjee S. An evaluation of strain amplification concepts via Monte Carlo simulations of an ideal composite. Rubber Chemistry and Technology, 1997, 70: 25-37 doi: 10.5254/1.3538416

    [17]

    Bergström JS, Boyce MC. Mechanical behavior of particle filled elastomers. Rubber Chemistry and Technology, 1999, 72: 633-656 doi: 10.5254/1.3538823

    [18]

    Simo JC. On a fully three-dimensional finite-strain viscoelastic da-mage model: Formulation and computational aspects. Computer Methods in Applied Mechanics & Engineering, 1987, 60: 153-173

    [19]

    Miehe C. Discontinuous and continuous damage evolution in Ogde-n-type large-strain elastic materials. European Journal of Mechanics A/Solids, 1995, 14: 697-720

    [20]

    Miehe C, Keck J. Superimposed finite elastic-viscoelastic-plastoelas-tic stress response with damage in filled rubbery polymers. Experi-ments, modeling and algorithmic implementation. Journal of Mech-anics and Physicals of Solids, 2000, 48(2): 323-365 doi: 10.1016/S0022-5096(99)00017-4

    [21]

    Kaliske M, Nasdala L, Rothert H. On damage modeling for elasticand viscoelastic materials at large strain. Computer & Structures, 2001, 79: 2133-2141

    [22]

    Lin RC, Schomburg U. A finite elastic-viscoelastic-elastoplastic material low with damage: theoretical and numerical aspects. Com-puter Methods in Applied Mechanics & Engineering, 2003, 192: 1591-1627

    [23]

    Ogden RW, Roxburgh DG. A pseudo-elastic model for the Mullins effect in filled rubber. Proceedings of the Royal Society, 1999, 455: 2861-2878 doi: 10.1098/rspa.1999.0431

    [24]

    Li on, A. A constitutive model for carbon black filled rubber: Exp-erimental investigations and mathematical representation. Continuum Mechanics and Thermodynamics, 1996, 8(3): 153-169 doi: 10.1007/BF01181853

    [25]

    Li on, A. On the large deformation behaviour of reinforced rubber at different temperatures. Journal of the Mechanics and Physics of Solids, 1997, 45(11-12): 1805-1834 doi: 10.1016/S0022-5096(97)00028-8

    [26]

    Diani J, Brieu M, Vacherand JM. A damage directional constitutive model for Mullins effect with permanent set and induced anisotropy. European Journal of Mechanics A/Solids, 2006, 25: 483-496 doi: 10.1016/j.euromechsol.2005.09.011

    [27]

    Dorfmann A, Ogden RW. A constitutive model for the Mullins eff-ect with permanent set in particle-reinforced rubber. International Journal of Solids and Structures, 2004, 41: 1855-1878 doi: 10.1016/j.ijsolstr.2003.11.014

    [28]

    Pawelski H. Softening behavior of elastomeric media after loading in changing directions//Besdo D, Schuster RH, Ilheman J, Eds. Proceedings of the Second European Conference of the Constitutive Models for Rubber. Balkema, 1999, 27-34.

    [29] 谈炳东, 许进升, 孙朝翔, 等. 短纤维增强三元乙丙橡胶横观各向同性黏—超弹性本构模型. 力学学报, 2017, 49(3): 677-684 (Tan Bingdong, Xu Jinsheng, Sun Chaoxiang, et al. A transversely isotropic visco-hyperelastic constitutive model for short fiber reinforced EPDM. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(3): 677-684 (in Chinese) doi: 10.6052/0459-1879-16-380
    [30] 谈炳东, 许进升, 贾云飞, 等. 短纤维增强 EPDM 包覆薄膜超弹性本构模型. 力学学报, 2017, 49(2): 317-323 (Tan Bingdong, Xu Jin-sheng, Jia Yunfei, et al. Hyperelastic constitutive model for short fiber reinforced EPDM inhibitor film. Chinese Journal of Theoreti-cal and Applied Mechanics, 2017, 49(2): 317-323 (in Chinese) doi: 10.6052/0459-1879-16-324
    [31]

    Mitsuteru A, Yoshiyuki K, Yoshimi S, et al. Constitutive modeling for texture reinforced rubber by using an anisotropic viscohypere-lastic model. Doboku Gakkai Ronbunshuu A, 2010, 66(2): 194-205 doi: 10.2208/jsceja.66.194

    [32]

    Rebouah M, Chagnon G. Permanent set and stress softening consti-tutive equation applied to rubber like materials and soft tissues. Acta Mechanica, 2013, 225(6): 1685-1698

    [33]

    Nasdala, L, Kempe A, Rolfes, R. An elastic molecular model for rubber inelasticity. Computational Materials Science, 2015, 106: 83-99 doi: 10.1016/j.commatsci.2015.04.036

    [34]

    Külcü ID, Tanrverdi HB. A constitutive model for hysteresis: the continuum damage approach for filled rubber-like materials. Archive of Applied Mechanics, 2020, 90: 1771-1781 doi: 10.1007/s00419-020-01695-2

    [35]

    Xiao H. An explicit, direct approach to obtaining multi-axial elasti-cpotentials which exactly match data of four benchmark tests for incompressible rubberlike materials-Part 1: incompressible deformat-ions. Acta Mechanica, 2012, 223(9): 2309-2063

    [36]

    Xiao H, Ding XF, Cao J, et al. New multi-axial constitutive mod-dels for large elastic deformation behaviors of soft solids up to breaking. International Journal of Solids and Structures, 2017, 109: 123-130 doi: 10.1016/j.ijsolstr.2017.01.013

    [37]

    Wang XM, Li H, Yin ZN, et al. Multiaxial strain energy functions of rubberlike materials: An explicit approach based on polynomial interpolation. Rubber Chemistry and Technology, 2014, 87(1): 168-183 doi: 10.5254/rct.13.86960

    [38] 王晓明, 吴荣兴, 肖衡. 显式模拟类橡胶材料Mullins效应滞回圈. 力学学报, 2019, 51(2): 484-493 (Wang Xiaoming, Wu Rongxing, Xiao Heng. Explicit modeling the hysteresis loops of the mullins effect for rubber-like materials. Chinses Journal of Theoretical and Applied Mechanics, 2019, 51(2): 484-493 (in Chinese) doi: 10.6052/0459-1879-18-334
    [39] 王晓明, 郑东, 吴荣兴等. 类橡胶材料变形直到破坏的显式本构模型. 力学季刊, 2019, 40(2): 252-264 (Wang Xiaoming, Zheng Dong, Wu Rongxing, et al. Explicit constitutive model of rubber-like materials up to failure. Chinese Quarterly of Mechanics, 2019, 40(2): 252-264 (in Chinese)
    [40] 王晓明, 肖衡. 显式方法精确模拟形状记忆聚合物热力学行为. 固体力学学报, 2020, 41(4): 366-378 (Wang Xiaoming, Xiao Heng. Accurately modeling thermomechancial behavior of shape memory polymer with explicit Method. Chinese Journal of Solid Mechani-cs, 2020, 41(4): 366-378 (in Chinese)
    [41]

    Xiao H, Chen LS. Hencky’s logarithmic strain and dual stress-stra-in and strain-stress relations isotropic finite hyperelasticity. Internat-ional Journal of Solids and Structures, 2003, 40(6): 1455-146 doi: 10.1016/S0020-7683(02)00653-4

  • 期刊类型引用(7)

    1. 郑刚,刘新福,刘广胜,魏韦,邓泽鲲,刘春花. 井口简易防喷器密封胶芯形变及其密封特性分析. 机床与液压. 2025(01): 140-145 . 百度学术
    2. 宋义虎,李承宇. 橡胶复合材料的Mullins效应. 工程塑料应用. 2024(06): 70-81 . 百度学术
    3. 王晓明,田兴兴,张振,肖衡. 显式方法模拟类橡胶材料率相关Mullins效应. 力学学报. 2024(12): 3553-3563 . 本站查看
    4. 龚臣成 ,陈艳 ,戴兰宏 . 聚脲弹性体力学性能与本构关系研究进展. 力学学报. 2023(01): 1-23 . 本站查看
    5. 郭鑫,陈素文. 考虑Mullins效应的硅酮胶本构模型. 力学学报. 2023(06): 1308-1318 . 本站查看
    6. 王晓明,肖衡. 基于有理插值方法模拟SMAs循环加载下的变形行为. 应用数学和力学. 2023(06): 694-707 . 百度学术
    7. 王晓明,王译仪,肖衡,金天付. 显式方法精确模拟形状记忆合金循环荷载下的变形行为. 力学季刊. 2022(03): 614-628 . 百度学术

    其他类型引用(0)

图(7)
计量
  • 文章访问数:  903
  • HTML全文浏览量:  396
  • PDF下载量:  65
  • 被引次数: 7
出版历程
  • 收稿日期:  2021-02-02
  • 录用日期:  2021-06-07
  • 网络出版日期:  2021-06-07
  • 刊出日期:  2021-07-17

目录

/

返回文章
返回