CN106202683A - 一种考虑相关效应的橡胶材料粘超弹本构模型的建模方法 - Google Patents

Abstract

本发明公开了一种考虑相关效应的橡胶材料粘超弹本构模型的建模方法，方法步骤如下：通过准静态单轴拉伸实验，获得橡胶材料的应力‑应变，通过SHTB装置进行冲击单轴拉伸实验，获得高应变率下橡胶材料的应力‑应变；根据准静态单轴拉伸实验和冲击单轴拉伸实验获得的橡胶材料的应力‑应变，获得其对应的响应规律和变形特性；再构建橡胶材料粘超弹本构模型。本发明针对橡胶材料拉伸力学性能的研究，建立橡胶材料粘超弹本构模型，为橡胶材料拉伸力学性能研究提供参考，解决了不同应变率下，橡胶材料力学表征问题。

Description

一种考虑相关效应的橡胶材料粘超弹本构模型的建模方法

技术领域

本发明属于橡胶材料拉伸力学性能研究的技术领域，具体涉及一种考虑相关效应的橡胶材料粘超弹本构模型的建模方法。

背景技术

橡胶材料的本构研究应包括超弹性和粘弹性。目前，关于橡胶拉伸力学性能的研究往往局限于超弹性，有限的关于橡胶率相关本构模型的研究是将超弹性本构模型和粘弹性本构模型进行组合，如H.Pouriayevali等将本构方程表示为与坐标系无关的形式，提出了描述橡胶类材料粘超弹性的本构模型，然而由于建模过程中引入了由应变不变量组成的多项式，参数较多且物理含义不明确，可能会出现多解状况。

发明内容

本发明的目的在于提供一种考虑相关效应的橡胶材料粘超弹本构模型的建模方法，针对橡胶材料拉伸力学性能的研究，建立橡胶材料粘超弹本构模型，为橡胶材料拉伸力学性能研究提供参考，解决了不同应变率下，橡胶材料力学表征问题。

实现本发明目的的技术解决方案为：

一种考虑相关效应的橡胶材料粘超弹本构模型的建模方法，方法步骤如下：

步骤一，对不同应变率下橡胶材料进行单轴拉伸实验，获得橡胶材料的响应规律和变形特性：

步骤1-1，通过准静态单轴拉伸实验，获得橡胶材料的应力-应变；

步骤1-2、通过SHTB装置进行冲击单轴拉伸实验，获得高应变率下橡胶材料的应力-应变；

步骤1-3、根据准静态单轴拉伸实验和冲击单轴拉伸实验获得的橡胶材料的应力-应变，获得其对应的响应规律和变形特性；

步骤二、构建橡胶材料粘超弹本构模型：

步骤2-1、构建Exp-ln超弹性本构模型；

步骤2-2、构建基于非线性粘弹性方法的粘超弹本构模型。

步骤2-1、构建Exp-ln超弹性本构模型，具体方式如下：

在准静态加载条件下应变率橡胶表现为超弹性，设单轴加载方向的伸长比为λ，橡胶为不可压缩材料，则三向主伸长比为：λ₁＝λ，λ₂＝λ₃＝λ^-1/2，其中λ₁为X轴方向主伸长比，λ₂为Y轴方向主伸长比，λ₃为Z轴方向主伸长比；因此，变形梯度张量F及左Cauchy应变张量B为

$F=\left|\begin{array}{ccc}\mathrm{\λ}& 0& 0\\ 0& {\mathrm{\λ}}^{-1/2}& 0\\ 0& 0& {\mathrm{\λ}}^{-1/2}\end{array}\right|,B=\left|\begin{array}{ccc}{\mathrm{\λ}}^2& 0& 0\\ 0& {\mathrm{\λ}}^{-1}& 0\\ 0& 0& {\mathrm{\λ}}^{-1}\end{array}\right|---\left(1\right)$

超弹性材料的应力应变关系用应变能函数W来度量，根据能量守恒可得，各向同性不可压缩橡胶材料的本构方程为

${\mathrm{\σ}}_e=-P_{e}I+2(\frac{\∂W_e}{\∂I_1^B}+I_1^B\frac{\∂W_e}{\∂I_2^B})B-2\frac{\∂W_e}{\∂I_2^B}B.B---\left(2\right)$

式中，σ_e为Cauchy应力，P_e为静水压力，左Cauchy应变张量B的第一不变量左Cauchy应变张量B的第二不变量I为单位矩阵；

Exp-ln应变能函数W_e为

$W_e=A\[\frac{1}{a}\exp \left(a\right(I_1^B-3\left)\right)+b(I_1^B-2)\×(1-ln(I_1^B-2\left)\right)-\frac{1}{a}-b\]---\left(3\right)$

式中，超弹性参数A＝G/2，G为橡胶材料剪切模量，超弹性参数a与链节伸长极限有关，超弹性参数b与中等应变量值有关；

联立式(2)和式(3)，则加载方向超弹性应力及垂直于加载方向的横向超弹性应力为

${\mathrm{\σ}}_{11}^e=-P_e+2{\mathrm{A\λ}}^2\[\exp \left(a\right({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\left)\right)-bln({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\]---\left(4\right)$

${\mathrm{\σ}}_{22}^e={\mathrm{\σ}}_{33}^e=-P_e+2{\mathrm{A\λ}}^{-1}\[\exp \left(a\right({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\left)\right)-bln({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\]---\left(5\right)$

联立式(4)和式(5)，结合获得加载方向的超弹性本构模型为

${\mathrm{\σ}}_{11}^e=2A({\mathrm{\λ}}^2-{\mathrm{\λ}}^{-1})\[\exp \left(a\right({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\left)\right)-bln({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\]---\left(6\right)$

式中，λ＝1+ε₁₁，其中P₁₁为工程应力，ε₁₁为工程应变，由式(6)结合橡胶材料准静态单轴拉伸实验数据，利用最小二乘法拟合确定两个超弹性参数a、b的值。

步骤2-2中，构建基于非线性粘弹性方法的粘超弹本构模型，具体方法如下：

设粘超弹本构模型σ(ε,t)为

$\mathrm{\σ}(\mathrm{\ϵ},t)={\mathrm{\σ}}_0\left(\mathrm{\ϵ}\right)g\left(t\right)={\mathrm{\σ}}_e+{\mathrm{\σ}}_v={\mathrm{\σ}}_0\left(\mathrm{\ϵ}\right)\left(g_{\mathrm{\∞}}\right(\mathrm{\ϵ})+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{g}_i(\mathrm{\ϵ}\left)\exp \right(-\frac{t}{{\mathrm{\θ}}_i}\left)\right)---\left(7\right)$

式中，σ₀(ε)为橡胶材料的瞬态应力响应，σ_e＝σ₀(ε)g_∞(ε)，g_∞(ε)、g_i(ε)均为无量纲参数，表征超弹性和粘弹性的权重，θ_i为松弛时间，N为松弛时间的组数。σ_v为粘弹本构模型，g(t)为时间相关系数；

结合Leaderman卷积积分和式(7)，橡胶材料非线性粘超弹性本构σ为

$\mathrm{\σ}={\∫}_{-\mathrm{\∞}}^t\[\frac{\∂}{\∂\mathrm{\τ}}\left({\mathrm{\σ}}_0\right(\mathrm{\ϵ}\left)\right(g_{\mathrm{\∞}}\left(\mathrm{\ϵ}\right)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{g}_i\left(\mathrm{\ϵ}\right)\exp \left(-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_i}\right)\left)\right)\]d\mathrm{\τ}---\left(8\right)$

其中τ为特定时间；

将式(8)分解为超弹性应力部分和粘弹性应力部分，三维全量格式粘超弹本构模型为

$\mathrm{\σ}={\mathrm{\σ}}_e+{\mathrm{\σ}}_v={\mathrm{\σ}}_0\left(\mathrm{\ϵ}\right)\·g_{\mathrm{\∞}}\left(\mathrm{\ϵ}\right)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\∫}_{-\mathrm{\∞}}^t\frac{\∂}{\∂\mathrm{\τ}}\left({\mathrm{\σ}}_0\right(\mathrm{\ϵ})\·g_i(\mathrm{\ϵ}\left)\right)\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_i})d\mathrm{\τ}---\left(9\right)$

式(9)中 为瞬态应力的粘弹性部分；Y.Anani认为瞬态粘弹性应力和超弹性应力具有相同形式的应变能函数，因此，对于各向同性不可压缩的橡胶材料表示为

${\mathrm{\σ}}_0^v\left(\mathrm{\ϵ}\right)=-P_{v}I+2(\frac{\∂W_v}{\∂I_1^B}+I_1^B\frac{\∂W_v}{\∂I_2^B})B-2\frac{\∂W_v}{\∂I_2^B}B.B---\left(10\right)$

式中，A'、a'、b'均为粘弹性参数，P_v为静水压力；W_v为粘弹性应变能；

假设加载前材料的应力状态不影响加载后的应力状态，即时间积分下限为零，为减少模型参数数量，取N＝1时的松弛时间θ_i，联立式(9)和式(10)得

${\mathrm{\σ}}_v=-P_{v}I+{\∫}_0^t\frac{\∂}{\∂\mathrm{\τ}}\[2(\frac{\∂W_v}{\∂I_1^B}+I_1^B\frac{\∂W_v}{\∂I_2^B})B-2\frac{\∂W_v}{\∂I_2^B}B.B\]\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_1})d\mathrm{\τ}---\left(11\right)$

由式(11)，获得加载方向粘弹性应力及垂直于加载方向的横向粘弹性应力为

${\mathrm{\σ}}_{11}^v=-P_v+{\∫}_0^t\frac{\∂}{\∂t}\[2A^{\′}{\mathrm{\λ}}^2\left(\exp \right(a^{\′}\left(I_1^B-3\right))-b^{\′}ln(I_1^B-2\left)\right)\]\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_1})d\mathrm{\τ}---\left(12\right)$

${\mathrm{\σ}}_{22}^v={\mathrm{\σ}}_{33}^v=-P_v+{\∫}_0^t\frac{\∂}{\∂t}\[2A^{\′}{\mathrm{\λ}}^{-1}\left(\exp \right(a^{\′}\left(I_1^B-3\right))-b^{\′}ln(I_1^B-2\left)\right)\]\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_1})d\mathrm{\τ}---\left(13\right)$

联立式(12)和式(13)，结合得到在加载方向的非线性粘弹性本构关系为

${\mathrm{\σ}}_{11}^v={\∫}_0^t\frac{\∂}{\∂\mathrm{\τ}}\[2A^{\′}({\mathrm{\λ}}^2-{\mathrm{\λ}}^{-1})\left(\exp \right(a^{\′}\left(I_1^B-3\right))-b^{\′}ln(I_1^B-2\left)\right)\]\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_1})d\mathrm{\τ}---\left(14\right)$

在定应变率单轴拉伸作用下，应变率加载时间将上述关系代入式(14)并结合式(6)得橡胶材料加载方向的粘超弹本构模型σ₁₁为

$\begin{array}{c}{\mathrm{\σ}}_{11}={\mathrm{\σ}}_{11}^e+{\mathrm{\σ}}_{11}^v=2A({\mathrm{\λ}}^2-{\mathrm{\λ}}^{-1}\[\exp (a\left({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\right))-b\ln ({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\]+\\ 2A^{\′}{\mathrm{\θ}}_1(1-\exp (-\frac{\mathrm{\λ}-1}{\stackrel{\·}{\mathrm{\λ}}{\mathrm{\θ}}_1}\left)\right)(2\mathrm{\λ}\stackrel{\·}{\mathrm{\λ}}+\frac{\stackrel{\·}{\mathrm{\λ}}}{{\mathrm{\λ}}^2})\·\[\exp \left(a^{\′}\right({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\left)\right)-b^{\′}\ln ({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\end{array}---\left(15\right)$

通过式(15)结合动态单轴拉伸实验数据，在求得两个超弹性参数a、b的值的基础上，利用最小二乘法拟合确定三个粘弹性参数A'、a'、b'和松弛时间θ₁的值。

本发明与现有技术相比，其显著优点在于：相对传统的橡胶类材料粘超弹性的本构模型，本发明将Exp-ln超弹性本构模型与基于广义粘弹性方法的粘弹性本构模型进行组合，建立了一种预测应变范围大、参数少且物理意义明确的橡胶材料粘超弹本构模型。

附图说明

图1为本发明橡胶材料准静态单轴拉伸实验装置示意图。

图2为本发明分离式霍普金森拉伸装置(SHTB)示意图。

图3为本发明一维弹性波传播示意图。

图4为本发明橡胶试件8示意图。

图5为本发明带螺纹孔的铝合金入射杆示意图。

图6为本发明胶粘夹具示意图。

图7为本发明试件胶粘示意图。

图8为本发明橡胶试件SHTB实验连接示意图。

图9为本发明准静态下拉伸应力-应变实验值。

图10为本发明350s^-1下入射、反射、透射信号曲线图。

图11为本发明550s^-1下入射、反射、透射信号曲线图。

图12为本发明900s^-1下入射、反射、透射信号曲线图。

图13为本发明1200s^-1下入射、反射、透射信号曲线图。

图14为本发明1600s^-1下入射、反射、透射信号曲线图。

图15为本发明不同应变率下拉伸应力-应变实验值示意图。

图16为本发明橡胶材料霍普金森拉杆实验数值模型示意图。

图17为本发明准静态下应力-应变曲线对比图。

图18为本发明不同应变率下应力-应变曲线对比图。

图19为本发明考虑相关效应的橡胶材料粘超弹本构模型的建模方法的流程图。

具体实施方式

下面结合附图对本发明作进一步详细描述。

本发明的原理如下：

建立橡胶材料粘超弹本构模型，并利用自主搭建的准静态单轴拉伸实验装置和分离式霍普金森拉伸装置(SHTB)对其进行实验验证，为橡胶材料拉伸力学性能研究提供参考。首先，为了使本构模型具有预测应变范围大、参数少且物理意义明确的特点，将超弹性本构模型与粘弹性本构模型进行组合，建立了橡胶材料粘超弹本构模型；其次，常见Instron实验机在测试小尺寸试件时存在的夹持困难、量程过大等问题，采用自主搭建了一种准静态单轴拉伸实验装置，该实验装置结构简单，能够精确测量小尺寸试件在较小拉力作用下的应力应变响应；再次，针对SHTB装置中试件连接方式存在的问题，设计一种用于SHTB装置的新型试件连接方式，结合橡胶材料为板材这一实际情况，避免试件因连接受到破坏，增强了连接强度，提高了实验效率和实验精度；最后，为验证用户子程序和橡胶材料粘超弹本构模型的有效性，根据单轴拉伸实验获得的应力-应变曲线拟合了粘超弹本构模型材料参数，推导了本构模型的三维增量格式，根据ABAQUS软件子程序的相关约定，编写了本构模型用户子程序VUMAT，验证了用户子程序和橡胶材料粘超弹本构模型的有效性。

以H.Pouriayevali等关于橡胶本构模型的研究思路为例，结合图1至图19，本发明所述的一种考虑相关效应的橡胶材料粘超弹本构模型的建模方法，方法步骤如下：

步骤一，对不同应变率下橡胶材料进行单轴拉伸实验

针对橡胶材料试件，开展不同应变率下单轴拉伸实验，获得应力-应变响应规律和变形特性。

步骤1-1，通过准静态单轴拉伸实验，获得橡胶材料的应力-应变。

搭建橡胶材料准静态单轴拉伸实验装置，包括静态应变仪1、拉力传感器2、数显螺旋千分尺3、直线导轨4、支架5、试件固定角铝7。旋动数显螺旋千分尺3即可拉动直线导轨4和支架5，进而驱动拉力传感器2对橡胶试件8进行单轴加载；整个装置中，支架5起到支撑数显螺旋千分尺3、连接直线导轨4的滑块和拉力传感器2、支撑试件固定角铝7的作用；橡胶试件8两端通过502胶与拉力传感器2端部和试件固定角铝7粘结。拉力传感器2用来测量试件作用力F并通过静态应变仪1的0通道读取，试件8的相对位移ΔL通过数显螺旋千分尺3读取。力F-位移ΔL与试件工程应力P₁₁-工程应变ε₁₁的转换关系式如下：

$P_{11}=\frac{F}{A_0},{\mathrm{\ϵ}}_{11}=\frac{\mathrm{\Δ}L}{L}---\left(16\right)$

式中，A₀、L为试件的初始实验段受力面积和初始实验段长度，其中P₁₁为工程应力。

步骤1-2、通过SHTB装置进行冲击单轴拉伸实验，获得高应变率下橡胶材料的应力-应变。

杆系直径为12.7mm，入射杆9、透射杆10和吸收杆11长分别为600mm、1200mm和1200mm，管状子弹内径、外径和长度分别为30mm、40mm和300mm，橡胶和帘线/橡胶复合材料均为低阻抗材料，为保证阻抗匹配，入射杆9和透射杆10均采用铝合金材料。(SHTB装置包括气压控制系统12、数据采集系统13、高速摄影系统14、入射杆9、透射杆10、吸收杆11、管状子弹15。)

在铝合金入射杆9和透射杆10两端攻螺纹孔，加工用于胶粘的夹具，其端部与入射杆9和透射杆10螺纹孔配合，与试件8胶粘处为矩形槽。为提高实验效率，夹具加工十套，使用AB胶同时进行胶粘；为保证试件8与拉杆轴线的同轴度，在AB胶固结时，使用带槽支架支承夹具。将橡胶试件8固结后胶粘夹具与入射杆9和透射杆10装配。

入射杆9和透射杆10分别贴有两组电阻应变片，管状子弹15以一定速度撞击法兰盘时，在入射杆9端产生拉伸入射波形向杆的另一端(图中向右)传播，入射杆9上应变片G₁记录下入射应变波形，这一波形传入试件8后在试件8中经试件前后界面多次透射、反射，在试件中基本达到应力平衡。试件被看成为没有厚度的试件界面，应变波形在这个界面上，一部分被反射回入射杆9向左传播，再次由应变片G₁记录；其余部分则透射到透射杆10内继续向右传播，透射杆10上应变片G₂记录下透射应变波形。在实验中实测的是由数据采集系统13获得的电压-时间曲线，通过电压与应变的对应关系，即可将电压波形转变应变波形。

根据一维弹性应力波的传播理论，利用在入射杆9和透射杆10上记录的信号ε_I、ε_R、ε_T，可以根据式(17)-式(20)计算出试件两个端面的力F₁、F₂和位移u₁、u₂。

F₁＝E_oA_o(ε_I+ε_R) (17)

F₂＝E_oA_oε_T (18)

u₁＝C_o∫(ε_I-ε_R)dt (19)

u₂＝C_o∫ε_Idt (20)

式中，ε_I是入射应变；ε_R是反射应变；ε_T是透射应变；E_o是杆的弹性模量；A_o是杆的横截面积；C_o是杆内弹性波的速度；t是时间。

在试件8内部满足应力平衡的条件下，试件8的应力σ_s、应变ε_s、应变率可用式(21)-式(23)表示：

${\mathrm{\σ}}_s=\frac{F_1+F_2}{2A_s}---\left(21\right)$

${\mathrm{\ϵ}}_s=\frac{u_1-u_2}{2L_s}---\left(22\right)$

${\stackrel{\·}{\mathrm{\ϵ}}}_s=\frac{{\mathrm{d\ϵ}}_s}{dt}---\left(23\right)$

式中，A_s、L_s分别为试件的横截面积；t为时间。

在动态应力平衡条件下，根据以下式(24)-式(26)可以得到材料在不同应变率下的工程应力σ_s(t)-工程应变ε_s(t)关系：

${\mathrm{\σ}}_s\left(t\right)=\frac{E_o{A}_o}{2A_s}{\mathrm{\ϵ}}_T\left(t\right)---\left(24\right)$

${\mathrm{\ϵ}}_s\left(t\right)=-\frac{2C_o}{L_s}{\∫}_0^t{\mathrm{\ϵ}}_R\left(t\right)dt---\left(25\right)$

${\stackrel{\·}{\mathrm{\ϵ}}}_s\left(t\right)=-\frac{2C_o}{L_s}{\mathrm{\ϵ}}_R\left(t\right)---\left(26\right)$

步骤1-3、根据准静态单轴拉伸实验和冲击单轴拉伸实验获得的橡胶材料的应力-应变，获得其对应的响应规律和变形特性。

通过静态应变仪1的0通道读取试件作用力F，通过数显螺旋千分尺3读取试件的相对位移ΔL，根据式(16)可以获得准静态下的工程应力P₁₁-工程应变ε₁₁实验值。

调节管状子弹15的撞击速度，得到橡胶试件8在350s^-1、550s^-1、900s^-1、1200s^-1、1600s^-1应变率下的冲击拉伸响应(考虑橡胶材料力学性能的离散性，每组应变率完成4组平行实验)。结合图10-图14，入射杆9、透射杆10上应变片G₁、G₂通过数据采集系统获得入射、反射、透射信号，实测的各应变率下入射、反射、透射信号。

由实测电压信号波形获得橡胶材料冲击拉伸规律的方法如下：将反射波和透射波电压信号去除零漂、截取有效段；根据电压与应变的对应关系，将有效电压信号转变为入射应变ε_I、反射应变ε_R和透射应变ε_T；由透射应变ε_T利用式(24)可以得到试件8的应力时间曲线，由反射波ε_R利用式(26)可以得到试件8的应变率时间曲线；利用式(25)即反射波ε_R对时间积分可以得到试件8的应变时间曲线，最后利用时间对应关系就可以得到试件的工程应力-工程应变关系。

利用上述方法，获得橡胶试件在350s^-1、550s^-1、900s^-1、1200s^-1、1600s^-1应变率下的工程应力P₁₁-工程应变ε₁₁实验值。结合图9和图15，得出橡胶试件的应力随应变的增大而增大；随应变率的增大，相同应变下，橡胶试件的应力增大。橡胶材料的力学响应表现出明显的率相关性。

步骤二、构建橡胶材料粘超弹本构模型

将Exp-ln超弹性本构模型与基于广义粘弹性方法的粘弹性本构模型进行组合，建立了一种预测应变范围大、参数少且物理意义明确的橡胶材料粘超弹本构模型。

步骤2-1、构建Exp-ln超弹性本构模型

准静态加载条件下橡胶表现为超弹性，设单轴加载方向的伸长比为λ，橡胶为不可压缩材料，则三向主伸长比为：λ₁＝λ，λ₂＝λ₃＝λ^-1/2，其中λ₁为X轴方向主伸长比，λ₂为Y轴方向主伸长比，λ₃为Z轴方向主伸长比；因此，变形梯度张量F及左Cauchy应变张量B为

$F=\left|\begin{array}{ccc}\mathrm{\λ}& 0& 0\\ 0& {\mathrm{\λ}}^{-1/2}& 0\\ 0& 0& {\mathrm{\λ}}^{-1/2}\end{array}\right|,B=\left|\begin{array}{ccc}{\mathrm{\λ}}^2& 0& 0\\ 0& {\mathrm{\λ}}^{-1}& 0\\ 0& 0& {\mathrm{\λ}}^{-1}\end{array}\right|---\left(1\right)$

超弹性材料的应力应变关系用应变能函数W来度量，根据能量守恒可得，各向同性不可压缩橡胶材料的本构方程为

${\mathrm{\σ}}_e=-P_{e}I+2(\frac{\∂W_e}{\∂I_1^B}+I_1^B\frac{\∂W_e}{\∂I_2^B})B-2\frac{\∂W_e}{\∂I_2^B}B.B---\left(2\right)$

式中，σ_e为Cauchy应力，P_e为静水压力，左Cauchy应变张量B的第一不变量左Cauchy应变张量B的第二不变量I为单位矩阵。

采用参数具有实际物理意义且适用应变范围较广的Exp-ln应变能函数，所述Exp-ln应变能函数W_e为

$W_e=A\[\frac{1}{a}\exp \left(a\right(I_1^B-3\left)\right)+b(I_1^B-2)\×(1-ln(I_1^B-2\left)\right)-\frac{1}{a}-b\]---\left(3\right)$

式中，超弹性参数A＝G/2，G为橡胶材料剪切模量，超弹性参数a与链节伸长极限有关，超弹性参数b与中等应变量值有关。A、a、b分别决定了橡胶材料小、中、大应变阶段的力学响应，因此，该应变能函数的应变预测范围较广且需拟合的参数较少。

联立式(2)和式(3)，则加载方向超弹性应力及垂直于加载方向的横向超弹性应力为

${\mathrm{\σ}}_{11}^e=-P_e+2{\mathrm{A\λ}}^2\[\exp \left(a\right({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\left)\right)-bln({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\]---\left(4\right)$

${\mathrm{\σ}}_{22}^e={\mathrm{\σ}}_{33}^e=-P_e+2{\mathrm{A\λ}}^{-1}\[\exp \left(a\right({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\left)\right)-bln({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\]---\left(5\right)$

联立式(4)和式(5)，结合可得加载方向的超弹性本构模型为

${\mathrm{\σ}}_{11}^e=2A({\mathrm{\λ}}^2-{\mathrm{\λ}}^{-1})\[\exp \left(a\right({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\left)\right)-bln({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\]---\left(6\right)$

式中，λ＝1+ε₁₁，其中P₁₁为工程应力，ε₁₁为工程应变，由式(6)结合橡胶材料准静态单轴拉伸实验数据，利用最小二乘法拟合确定两个超弹性参数a、b的值。

步骤2-2、构建基于非线性粘弹性方法的粘超弹本构模型

结合C.Gamonpilas给出的非线性粘弹性本构模型的描述方法，将应变ε和时间t对应力的影响效应进行解耦，避免了引入多项式，粘超弹本构模型σ(ε,t)为

$\mathrm{\σ}(\mathrm{\ϵ},t)={\mathrm{\σ}}_0\left(\mathrm{\ϵ}\right)g\left(t\right)={\mathrm{\σ}}_e+{\mathrm{\σ}}_v={\mathrm{\σ}}_0\left(\mathrm{\ϵ}\right)\left(g_{\mathrm{\∞}}\right(\mathrm{\ϵ})+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{g}_i(\mathrm{\ϵ}\left)\exp \right(-\frac{t}{{\mathrm{\θ}}_i}\left)\right)---\left(7\right)$

式中，σ₀(ε)为橡胶材料的瞬态应力响应，σ_e＝σ₀(ε)g_∞(ε)，g_∞(ε)、g_i(ε)均为无量纲参数，表征超弹性和粘弹性的权重，θ_i为松弛时间，N为松弛时间的组数。σ_v为粘弹本构模型，g(t)为时间相关参数；

积分形式的本构方程能直接反映粘弹性材料的记忆特性，结合Leaderman卷积积分和式(7)，橡胶材料非线性粘超弹性本构σ为

$\mathrm{\σ}={\∫}_{-\mathrm{\∞}}^t\[\frac{\∂}{\∂\mathrm{\τ}}\left({\mathrm{\σ}}_0\right(\mathrm{\ϵ}\left)\right(g_{\mathrm{\∞}}\left(\mathrm{\ϵ}\right)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{g}_i\left(\mathrm{\ϵ}\right)\exp \left(-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_i}\right)\left)\right)\]d\mathrm{\τ}---\left(8\right)$

其中τ为特定时间

将式(8)分解为超弹性应力部分和粘弹性应力部分，三维全量格式粘超弹本构模型为

$\mathrm{\σ}={\mathrm{\σ}}_e+{\mathrm{\σ}}_v={\mathrm{\σ}}_0\left(\mathrm{\ϵ}\right)\·g_{\mathrm{\∞}}\left(\mathrm{\ϵ}\right)+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\∫}_{-\mathrm{\∞}}^t\frac{\∂}{\∂\mathrm{\τ}}\left({\mathrm{\σ}}_0\right(\mathrm{\ϵ})\·g_i(\mathrm{\ϵ}\left)\right)\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_i})d\mathrm{\τ}---\left(9\right)$

式(9)中 为瞬态应力的粘弹性部分。Y.Anani认为瞬态粘弹性应力和超弹性应力具有相同形式的应变能函数，因此，对于各向同性不可压缩的橡胶材料可表示为

${\mathrm{\σ}}_0^v\left(\mathrm{\ϵ}\right)=-P_{v}I+2(\frac{\∂W_v}{\∂I_1^B}+I_1^B\frac{\∂W_v}{\∂I_2^B})B-2\frac{\∂W_v}{\∂I_2^B}B.B---\left(10\right)$

式中，A'、a'、b'均为粘弹性参数，P_v为静水压力。W_v为粘弹性应变能；

假设加载前材料的应力状态不影响加载后的应力状态，即时间积分下限为零，为减少模型参数数量，取N＝1时的松弛时间θ_i，联立式(9)和式(10)可得

${\mathrm{\σ}}_v=-P_{v}I+{\∫}_0^t\frac{\∂}{\∂\mathrm{\τ}}\[2(\frac{\∂W_v}{\∂I_1^B}+I_1^B\frac{\∂W_v}{\∂I_2^B})B-2\frac{\∂W_v}{\∂I_2^B}B.B\]\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_1})d\mathrm{\τ}---\left(11\right)$

由式(11)可得加载方向粘弹性应力及垂直于加载方向的横向粘弹性应力为

${\mathrm{\σ}}_{11}^v=-P_v+{\∫}_0^t\frac{\∂}{\∂t}\[2A^{\′}{\mathrm{\λ}}^2\left(\exp \right(a^{\′}\left(I_1^B-3\right))-b^{\′}ln(I_1^B-2\left)\right)\]\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_1})d\mathrm{\τ}---\left(12\right)$

${\mathrm{\σ}}_{22}^v={\mathrm{\σ}}_{33}^v=-P_v+{\∫}_0^t\frac{\∂}{\∂t}\[2A^{\′}{\mathrm{\λ}}^{-1}\left(\exp \right(a^{\′}\left(I_1^B-3\right))-b^{\′}ln(I_1^B-2\left)\right)\]\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_1})d\mathrm{\τ}---\left(13\right)$

联立式(12)和式(13)，结合可得在加载方向的非线性粘弹性本构关系为

${\mathrm{\σ}}_{11}^v={\∫}_0^t\frac{\∂}{\∂\mathrm{\τ}}\[2A^{\′}({\mathrm{\λ}}^2-{\mathrm{\λ}}^{-1})\left(\exp \right(a^{\′}\left(I_1^B-3\right))-b^{\′}ln(I_1^B-2\left)\right)\]\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_1})d\mathrm{\τ}---\left(14\right)$

在定应变率单轴拉伸作用下，应变率加载时间将上述关系代入式(14)并结合式(6)得橡胶材料加载方向的粘超弹本构模型为

$\begin{array}{c}{\mathrm{\σ}}_{11}={\mathrm{\σ}}_{11}^e+{\mathrm{\σ}}_{11}^v=2A({\mathrm{\λ}}^2-{\mathrm{\λ}}^{-1}\[\exp (a\left({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\right))-b\ln ({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\]+\\ 2A^{\′}{\mathrm{\θ}}_1(1-\exp (-\frac{\mathrm{\λ}-1}{\stackrel{\·}{\mathrm{\λ}}{\mathrm{\θ}}_1}\left)\right)(2\mathrm{\λ}\stackrel{\·}{\mathrm{\λ}}+\frac{\stackrel{\·}{\mathrm{\λ}}}{{\mathrm{\λ}}^2})\·\[\exp \left(a^{\′}\right({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-3\left)\right)-b^{\′}\ln ({\mathrm{\λ}}^2+2{\mathrm{\λ}}^{-1}-2)\end{array}---\left(15\right)$

通过式(15)结合动态单轴拉伸实验数据，在求得两个超弹性参数a、b的值的基础上，利用最小二乘法拟合确定三个粘弹性参数A'、a'、b'和松弛时间θ₁的值，如表1。由于采用了参数少且物理意义明确的Exp-ln超弹性本构模型和建立在应变、时间效应解耦基础上的非线性粘弹性本构模型，提出的粘超弹本构模型具有预测应变范围大、参数少且物理意义明确的特点。表1为本发明粘超弹本构模型参数值。

表1

步骤三、对粘超弹本构模型进行实验验证。

步骤3-1、对粘超弹本构模型进行参数拟合

将橡胶材料单轴拉伸实验数据代入粘超弹本构模型中，进行最小二乘法拟合，即可获得本构模型的6个待定参数，参数拟合步骤为：

(1)根据准静态加载条件下橡胶材料单轴拉伸实验数据，结合式(6)并采用最小二乘法确定两个超弹性参数a、b的值。

(2)将a、b代入式(15)中，选取应变率为900s^-1的曲线进行非线性拟合，确定三个粘弹性参数A'、a'、b'和松弛时间θ₁的值。至此橡胶材料粘超弹本构模型中所有参数全部求得。

步骤3-2、粘超弹本构模型三维增量格式与VUMAT开发

为进行橡胶材料粘超弹本构模型VUMAT子程序开发，对式(9)三维全量格式本构进行有限元离散，建立三维应力增量与应变增量之间的关系。

针对橡胶粘超弹本构(9)中的超弹本构部分，采用Kirchhoff应力张量S与Green应变张量E获得超弹本构模型σ_e的应力更新。设变量G＝2E+I，E、G用分量形式表示为(i,j＝1,2,3)，e_i为正交曲线坐标的基矢，ξ为应变。则三个应变不变量为即有

$\frac{\∂I_1^B}{\∂G}=I,\frac{\∂I_2^B}{\∂G}=2G,\frac{\∂I_3^B}{\∂G}=3G^2---\left(27\right)$

$\left\{\begin{array}{c}\frac{{\∂}^2{I}_1^B}{\∂G^2}=0e_i\⊗e_j\⊗e_k\⊗e_l\\ \frac{{\∂}^2{I}_2^B}{\∂G^2}=2{\mathrm{\δ}}_{ki}{\mathrm{\δ}}_{lj}{e}_i\⊗e_j\⊗e_k\⊗e_l\\ \frac{{\∂}^2{I}_3^B}{\∂G^2}=3({\mathrm{\δ}}_{ki}{g}_{ij}+g_{jk}{\mathrm{\δ}}_{lj})e_i\⊗e_j\⊗e_k\⊗e_l\end{array}\right.,k,l=1,2,3---\left(28\right)$

设由应变能函数定义S对E偏导数

又

$\frac{\∂}{\∂G}\left(W_{ei}\frac{\∂I_i}{\∂G}\right)=\frac{\∂I_i}{\∂G}\⊗\frac{\∂W_{ei}}{\∂G}+\frac{\∂}{\∂G}\left(\frac{\∂I_i}{\∂G}\right)W_{ei}=\frac{\∂I_i}{\∂G}\⊗\left(W_{eij}\frac{\∂I_j}{\∂G}\right)+\frac{{\∂}^2{I}_i}{\∂G^2}{W}_e---\left(29\right)$

由式(27)得

$\left\{\begin{array}{c}\frac{\∂I_1}{\∂G}\⊗\left(W_{e1j}\frac{\∂I_j}{\∂G}\right)=I\⊗(W_{e11}I+2W_{e12}G+3W_{e13}{G}^2)\\ \frac{\∂I_2}{\∂G}\⊗\left(W_{e2j}\frac{\∂I_j}{\∂G}\right)=2G\⊗(W_{e21}I+2W_{e22}G+3W_{e23}{G}^2)\\ \frac{\∂I_3}{\∂G}\⊗\left(W_{e3j}\frac{\∂I_j}{\∂G}\right)=3G^2\⊗(W_{e31}I+2W_{e32}G+3W_{e33}{G}^2)\end{array}\right.---\left(30\right)$

令D_ijkle_i为切线本构张量，采用Kirchhoff应力张量S与Green应变张量E，超弹本构模型的三维增量格式为

dS＝D_ijkldE (31)

结合式(28)和式(29)，即可给出

D_ijkl＝4[W_e11δ_jiδ_lk+2W_e12δ_jig_kl+3W_e13δ_jig_kmg_ml+2W_e21g_ijδ_lk+4W_e22g_ijg_kl+6W_e23g_ijg_kmg_ml+3W_e31g_img_mjδ_lk+6W_e32g_img_mjδ_lk+6W_e32g_img_mjg_kl+9W_e33g_img_mjg_kng_nl+2W_e2δ_kiδ_lj+3W_e3(δ_kig_lj+g_ikδ_li)] (32)

将式(2)应变能函数W_e代入式(32)求得切线本构张量D_ijkle_i，通过编制VUMAT子程序获得Kirchhoff应力张量S和Green应变张量E的更新，对于不可压缩橡胶有σ_e＝FSF^T、E＝1/2(F^TF-ξ)，即可获得超弹本构模型σ_e的应力更新。

针对橡胶粘超弹本构中的粘弹本构部分，由于粘弹本构模型σ_v中粘性积分项与变形历史相关，采用递推格式进行求解获得粘弹本构模型σ_v的三维增量格式。取N＝1时的松弛时间，设时刻t的粘弹性应力为

${\mathrm{\σ}}_v^t=-P_{v}I+{\∫}_0^t\frac{\∂}{\∂\mathrm{\τ}}\left({\mathrm{\σ}}_0^v\right(\mathrm{\ξ}\left)\right)\exp (-\frac{t-\mathrm{\τ}}{{\mathrm{\θ}}_i})d\mathrm{\τ}---\left(33\right)$

则时刻t+Δt的粘弹性应力为

${\mathrm{\σ}}_v^{t+\mathrm{\Δ}t}=-P_{v}I+{\∫}_0^{t+\mathrm{\Δ}t}\frac{\∂}{\∂\mathrm{\τ}}\left({\mathrm{\σ}}_0^v\right(\mathrm{\ξ}\left)\right)\exp (-\frac{t+\mathrm{\Δ}t-\mathrm{\τ}}{{\mathrm{\θ}}_i})d\mathrm{\τ}---\left(34\right)$

对于式(34)采用分步积分，第一部分积分起止时间为(0,t)，第二部分积分起止时间为(t,Δt)，当Δt很小时，近似认为t和Δt时刻的应变率相等，式(34)可化简为

${\mathrm{\σ}}_v^{t+\mathrm{\Δ}t}={\mathrm{\θ}}_i\frac{{\mathrm{d\σ}}_0^v\left(\mathrm{\ξ}\right)}{d\mathrm{\τ}}(1-\exp (-\frac{\mathrm{\Δ}t}{{\mathrm{\θ}}_i}\left)\right)+{\mathrm{\σ}}_v^t\exp (-\frac{\mathrm{\Δ}t}{{\mathrm{\θ}}_i})---\left(35\right)$

式(35)即为实现粘弹本构模型σ_v应力更新的三维增量格式，其中项采用和超弹本构部分相同的更新方式。

至此，将超弹本构模型、粘弹本构模型的应力更新进行叠加，即获得了进行子程序开发所需的本构方程三维增量格式，子程序开发的任务就是在材料积分点上完成上式的应力计算，并对主程序进行应力更新。针对式(31)、式(35)和VUMAT所提供的数据接口，制定粘超弹本构模型子程序开发实施步骤为：

(1)从ABAQUS主程序中读入增量步开始时刻材料属性：a、b、A'、a'、b'、θ₁；变形梯度F^t、F^t+Δt；时间增量Δt等；

(2)从用户自定义状态变量载入Green应变张量E，根据E＝1/2(F^TF-ξ)求解Green应变增量dE；

(3)由式(31)计算Kirchhoff应力增量dS，根据σ_e＝FSF^T求解获得超弹本构模型σ_e的应力更新；

(4)由式(35)计算获得粘弹本构模型σ_v应力更新；

(5)将超弹本构模型σ_e的应力更新和粘弹本构模型σ_v应力更新进行叠加，并返回增量步结束时刻的应力；

(6)更新状态变量。

根据以上步骤，用FORTRAN语言编写子程序，从而将橡胶粘超弹本构模型嵌入到ABAQUS中。

步骤3-3，对粘超弹本构模型的有效性进行实验验证

为验证本构模型的有效性，将开发的橡胶粘超弹本构模型的子程序VUMAT用于模拟霍普金森拉杆实验。结合图16，橡胶试件8采用C3D8R单元进行网格划分，共计2148个单元，为了减小数值模型规模、简化计算，并没有建立实际的霍普金森拉伸装置(如法兰盘，吸收杆等)，而是利用了数值模拟的便捷性，直接将管状子弹15与入射杆9粘接在一起，采用显式动态方法进行计算，冲击加载以对管状子弹15施加瞬间速度实现，通过改变子弹的撞击速度来实现不同应变率的加载。

将霍普金森拉杆实验数值模型中显式动态算法修改为隐式静态算法，在管状子弹15上施加位移边界条件，即可模拟橡胶试件8的准静态单轴拉伸实验。

在霍普金森拉杆实验数值模型的入射杆9和透射杆10上按应变片粘贴位置选取单元，获得入射波、反射波和透射波的应变波形，根据不同应变率下应变波形的数值解，结合式(24)-式(26)，即可获得不同应变率下橡胶材料试件8拉伸应力-应变数值解。结合图17-图18，在工程应变ε₁₁<0.02范围内，工程应力数值解与实验值符合较好；提出的粘超弹本构模型能够预测橡胶材料的单轴加载方向力学响应，验证了粘超弹本构模型的有效性。

Claims (3)

1.一种考虑相关效应的橡胶材料粘超弹本构模型的建模方法，其特征在于，方法步骤如下：

步骤一，对不同应变率下橡胶材料进行单轴拉伸实验，获得橡胶材料的响应规律和变形特性：

步骤1-1，通过准静态单轴拉伸实验，获得橡胶材料的应力-应变；

步骤1-2、通过SHTB装置进行冲击单轴拉伸实验，获得高应变率下橡胶材料的应力-应变；

步骤1-3、根据准静态单轴拉伸实验和冲击单轴拉伸实验获得的橡胶材料的应力-应变，获得其对应的响应规律和变形特性；

步骤二、构建橡胶材料粘超弹本构模型：

步骤2-1、构建Exp-ln超弹性本构模型；

步骤2-2、构建基于非线性粘弹性方法的粘超弹本构模型。

2.根据权利要求1所述的考虑相关效应的橡胶材料粘超弹本构模型的建模方法，其特征在于，步骤2-1、构建Exp-ln超弹性本构模型，具体方式如下：

在准静态加载条件下应变率橡胶表现为超弹性，设单轴加载方向的伸长比为λ，橡胶为不可压缩材料，则三向主伸长比为：λ₁＝λ，λ₂＝λ₃＝λ^-1/2，其中λ₁为X轴方向主伸长比，λ₂为Y轴方向主伸长比，λ₃为Z轴方向主伸长比；因此，变形梯度张量F及左Cauchy应变张量B为

$F=\left|\begin{array}{ccc}\mathrm{\&lambda;}& 0& 0\\ 0& {\mathrm{\&lambda;}}^{-1/2}& 0\\ 0& 0& {\mathrm{\&lambda;}}^{-1/2}\end{array}\right|,B=\left|\begin{array}{ccc}{\mathrm{\&lambda;}}^2& 0& 0\\ 0& {\mathrm{\&lambda;}}^{-1}& 0\\ 0& 0& {\mathrm{\&lambda;}}^{-1}\end{array}\right|---\left(1\right)$

超弹性材料的应力应变关系用应变能函数W来度量，根据能量守恒可得，各向同性不可压缩橡胶材料的本构方程为

${\mathrm{\&sigma;}}_e=-P_{e}I+2(\frac{\&part;W_e}{\&part;I_1^B}+I_1^B\frac{\&part;W_e}{\&part;I_2^B})B-2\frac{\&part;W_e}{\&part;I_2^B}B.B---\left(2\right)$

式中，σ_e为Cauchy应力，P_e为静水压力，左Cauchy应变张量B的第一不变量左Cauchy应变张量B的第二不变量I为单位矩阵；

Exp-ln应变能函数W_e为

$W_e=A\&lsqb;\frac{1}{a}\exp \left(a\right(I_1^B-3\left)\right)+b(I_1^B-2)\&times;(1-ln(I_1^B-2\left)\right)-\frac{1}{a}-b\&rsqb;---\left(3\right)$

式中，超弹性参数A＝G/2，G为橡胶材料剪切模量，超弹性参数a与链节伸长极限有关，超弹性参数b与中等应变量值有关；

联立式(2)和式(3)，则加载方向超弹性应力及垂直于加载方向的横向超弹性应力为

${\mathrm{\&sigma;}}_{11}^e=-P_e+2{\mathrm{A\&lambda;}}^2\&lsqb;\exp \left(a\right({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-3\left)\right)-bln\left({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-2\right)\&rsqb;---\left(4\right)$

${\mathrm{\&sigma;}}_{22}^e={\mathrm{\&sigma;}}_{33}^e=-P_e+2{\mathrm{A\&lambda;}}^{-1}\&lsqb;\exp \left(a\right({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-3\left)\right)-bln({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-2)\&rsqb;---\left(5\right)$

联立式(4)和式(5)，结合获得加载方向的超弹性本构模型为

${\mathrm{\&sigma;}}_{11}^e=2A({\mathrm{\&lambda;}}^2-{\mathrm{\&lambda;}}^{-1})\&lsqb;\exp \left(a\right({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-3\left)\right)-bln({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-2)\&rsqb;---\left(6\right)$

式中，λ＝1+ε₁₁，其中P₁₁为工程应力，ε₁₁为工程应变，由式(6)结合橡胶材料准静态单轴拉伸实验数据，利用最小二乘法拟合确定两个超弹性参数a、b的值。

3.根据权利要求1所述的考虑相关效应的橡胶材料粘超弹本构模型的建模方法，其特征在于，步骤2-2中，构建基于非线性粘弹性方法的粘超弹本构模型，具体方法如下：

设粘超弹本构模型σ(ε,t)为

$\mathrm{\&sigma;}(\mathrm{\&epsiv;},t)={\mathrm{\&sigma;}}_0\left(\mathrm{\&epsiv;}\right)g\left(t\right)={\mathrm{\&sigma;}}_e+{\mathrm{\&sigma;}}_v={\mathrm{\&sigma;}}_0\left(\mathrm{\&epsiv;}\right)\left(g_{\mathrm{\&infin;}}\right(\mathrm{\&epsiv;})+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_i(\mathrm{\&epsiv;}\left)\exp \right(-\frac{t}{{\mathrm{\&theta;}}_i}\left)\right)---\left(7\right)$

式中，σ₀(ε)为橡胶材料的瞬态应力响应，σ_e＝σ₀(ε)g_∞(ε)，g_∞(ε)、g_i(ε)均为无量纲参数，表征超弹性和粘弹性的权重，θ_i为松弛时间，N为松弛时间的组数。σ_v为粘弹本构模型，g(t)为时间相关系数；

结合Leaderman卷积积分和式(7)，橡胶材料非线性粘超弹性本构σ为

$\mathrm{\&sigma;}={\&Integral;}_{-\mathrm{\&infin;}}^t\&lsqb;\frac{\&part;}{\&part;\mathrm{\&tau;}}\left({\mathrm{\&sigma;}}_0\left(\mathrm{\&epsiv;}\right)\left(g_{\mathrm{\&infin;}}\left(\mathrm{\&epsiv;}\right)+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_i\left(\mathrm{\&epsiv;}\right)\exp \left(-\frac{t-\mathrm{\&tau;}}{{\mathrm{\&theta;}}_i}\right)\right)\right)\&rsqb;d\mathrm{\&tau;}---\left(8\right)$

其中τ为特定时间；

将式(8)分解为超弹性应力部分和粘弹性应力部分，三维全量格式粘超弹本构模型为

$\mathrm{\&sigma;}={\mathrm{\&sigma;}}_e+{\mathrm{\&sigma;}}_v={\mathrm{\&sigma;}}_0\left(\mathrm{\&epsiv;}\right)\&CenterDot;g_{\mathrm{\&infin;}}\left(\mathrm{\&epsiv;}\right)+\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{\&Integral;}_{-\mathrm{\&infin;}}^t\frac{\&part;}{\&part;\mathrm{\&tau;}}\left({\mathrm{\&sigma;}}_0\left(\mathrm{\&epsiv;}\right).g_i\left(\mathrm{\&epsiv;}\right)\right)\exp \left(-\frac{t-\mathrm{\&tau;}}{{\mathrm{\&theta;}}_i}\right)d\mathrm{\&tau;}---\left(9\right)$

式(9)中 为瞬态应力的粘弹性部分；Y.Anani认为瞬态粘弹性应力和超弹性应力具有相同形式的应变能函数，因此，对于各向同性不可压缩的橡胶材料表示为

${\mathrm{\&sigma;}}_0^v\left(\mathrm{\&epsiv;}\right)=-P_{v}I+2(\frac{\&part;W_v}{\&part;I_1^B}+I_1^B\frac{\&part;W_v}{\&part;I_2^B})B-2\frac{\&part;W_v}{\&part;I_2^B}B.B---\left(10\right)$

式中，A'、a'、b'均为粘弹性参数，P_v为静水压力；W_v为粘弹性应变能；

假设加载前材料的应力状态不影响加载后的应力状态，即时间积分下限为零，为减少模型参数数量，取N＝1时的松弛时间θ_i，联立式(9)和式(10)得

${\mathrm{\&sigma;}}_v=-P_{v}I+{\&Integral;}_0^t\frac{\&part;}{\&part;\mathrm{\&tau;}}\&lsqb;2(\frac{\&part;W_v}{\&part;I_1^B}+I_1^B\frac{\&part;W_v}{\&part;I_2^B})B-2\frac{\&part;W_v}{\&part;I_2^B}B.B\&rsqb;\exp (-\frac{t-\mathrm{\&tau;}}{{\mathrm{\&theta;}}_1})d\mathrm{\&tau;}---\left(11\right)$

由式(11)，获得加载方向粘弹性应力及垂直于加载方向的横向粘弹性应力为

${\mathrm{\&sigma;}}_{11}^v=-P_v+{\&Integral;}_0^t\frac{\&part;}{\&part;t}\&lsqb;2A^{\&prime;}{\mathrm{\&lambda;}}^2\left(\exp \right(a^{\&prime;}\left(I_1^B-3\right))-b^{\&prime;}ln(I_1^B-2\left)\right)\&rsqb;\exp (-\frac{t-\mathrm{\&tau;}}{{\mathrm{\&theta;}}_1})d\mathrm{\&tau;}---\left(12\right)$

${\mathrm{\&sigma;}}_{22}^v={\mathrm{\&sigma;}}_{33}^v=-P_v+{\&Integral;}_0^t\frac{\&part;}{\&part;t}\&lsqb;2A^{\&prime;}{\mathrm{\&lambda;}}^{-1}\left(\exp \right(a^{\&prime;}\left(I_1^B-3\right))-b^{\&prime;}ln(I_1^B-2\left)\right)\&rsqb;\exp (-\frac{t-\mathrm{\&tau;}}{{\mathrm{\&theta;}}_1})d\mathrm{\&tau;}---\left(13\right)$

联立式(12)和式(13)，结合得到在加载方向的非线性粘弹性本构关系为

${\mathrm{\&sigma;}}_{11}^v={\&Integral;}_0^t\frac{\&part;}{\&part;\mathrm{\&tau;}}\&lsqb;2A^{\&prime;}({\mathrm{\&lambda;}}^2-{\mathrm{\&lambda;}}^{-1})\left(\exp \right(a^{\&prime;}\left(I_1^B-3\right))-b^{\&prime;}ln(I_1^B-2\left)\right)\&rsqb;\exp (-\frac{t-\mathrm{\&tau;}}{{\mathrm{\&theta;}}_1})d\mathrm{\&tau;}---\left(14\right)$

在定应变率单轴拉伸作用下，应变率加载时间将上述关系代入式(14)并结合式(6)得橡胶材料加载方向的粘超弹本构模型σ₁₁为

$\begin{array}{c}{\mathrm{\&sigma;}}_{11}={\mathrm{\&sigma;}}_{11}^e+{\mathrm{\&sigma;}}_{11}^v=2A({\mathrm{\&lambda;}}^2-{\mathrm{\&lambda;}}^{-1}\&lsqb;\exp \left(a\left({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-3\right)\right)-b\ln \left({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-2\right)\&rsqb;+\\ 2A^{\&prime;}{\mathrm{\&theta;}}_1(1-\exp (-\frac{\mathrm{\&lambda;}-1}{\stackrel{\&CenterDot;}{\mathrm{\&lambda;}}{\mathrm{\&theta;}}_1}\left)\right)\left(2\mathrm{\&lambda;}\stackrel{\&CenterDot;}{\mathrm{\&lambda;}}+\frac{\stackrel{\&CenterDot;}{\mathrm{\&lambda;}}}{{\mathrm{\&lambda;}}^2}\right).\&lsqb;\exp \left(a^{\&prime;}\right({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-3\left)\right)-b^{\&prime;}ln({\mathrm{\&lambda;}}^2+2{\mathrm{\&lambda;}}^{-1}-2)\end{array}---\left(15\right)$

通过式(15)结合动态单轴拉伸实验数据，在求得两个超弹性参数a、b的值的基础上，利用最小二乘法拟合确定三个粘弹性参数A'、a'、b'和松弛时间θ₁的值。