试验所用的离子型中间膜为SG(SentryGlas^®，简称SG)，试件切割于同一卷SG材料，名义厚度为1.52 mm。考虑到试件尺寸对动态拉伸试验结果的可信度有较大影响，试验前结合数值模拟和预试验，对GB/T528 − 2009^[13]中的1型~4型四种试件尺寸以及Zhang等^[8]采用的试件尺寸(下文称为ZXH型)进行考察。根据试件在中高应变率下单轴拉伸试验中的断裂位置和工程应变率恒定性，确定1 s⁻¹、10 s⁻¹和100 s⁻¹下采用GB/T528 − 2009^[13]中的2型试件，300 s⁻¹、500 s⁻¹和800 s⁻¹下采用ZXH型试件。试件尺寸如图1所示。

图  1  试件尺寸　/mm

Figure 1.  Dimensions of specimens

试验加载装置为INSTRON VHS 高速液压伺服试验机，最大拉伸速率为20 m/s。试验过程中采用PHOTRON FASTCAM SA-X 高速摄影仪记录试件标距段的变形，试验装置如图2所示。试件标记如图3所示。

图  2  试验设备

Figure 2.  Test machine

图  3  试件的标记

Figure 3.  Marks on specimens

试验前，将试件放入温度箱中加热至预设温度并持温3 h以上。为记录试件的实际试验温度，试验开始时，从温度箱中取出两个相同的试件，一个用于加载，另一个用于测定加载完毕时试件的实际温度。

采用位移模式加载，加载速率根据预设应变率和试件标距段长度(见图1)确定。试验加载原理如图4所示，试件通过夹具与INSTRON拉伸试验机相连，下端固定，上端由试验机的驱动器施加恒定的拉伸速率V，从而给试件提供一个恒定的预设应变率。

图  4  加载原理

Figure 4.  Loading principle

试验在不同温度和应变率下进行，共42组工况：7种预设温度−40 ℃、−20 ℃、0 ℃、20 ℃、40 ℃、50 ℃、60 ℃，以及6种预设应变率1 s⁻¹、10 s⁻¹、100 s⁻¹、300 s⁻¹、500 s⁻¹、800 s⁻¹。每个工况进行三次重复试验。

图5所示为高速摄像仪记录的20 ℃和10 s⁻¹下试件加载至破化的全过程：加载初期，试件就发生明显的轴向拉伸变形和横向收缩，但直至被拉断试件也未出现明显的颈缩现象。

为确保试验量测结果能表征标矩段的材料力学行为，试件的断裂位置应位于标距段中。试验结果显示：在−40 ℃、−20 ℃、0 ℃和20 ℃下，试件的断裂位置基本都在标距段以内。而在40 ℃、50 ℃和60 ℃下，部分试件的断裂位置在标距段以外，但仍处在试件的直线段内并靠近标矩段，可以认为应力状态与标距段内基本一致。图6给出了20 ℃和40 ℃下试件的断裂位置。

图  5  试件典型加载至破坏过程 (20 ℃-10 s⁻¹-2)

Figure 5.  Typical deformation process of specimen 20 ℃-10 s⁻¹-2

图  6  试件断裂位置

Figure 6.  Fracture locations in specimens

试验过程中材料表现出明显的大变形特性，因此，采用真实应力和真实应变描述离子型中间膜材的力学性能，如式(1)和式(2)所示：

式中：${\sigma _{{\rm{t}}}}$和 ${\varepsilon _{{\rm{t}}}}$分别为真实应力和真实应变；${\sigma _{{\rm{e}}}}$为工程应力，即荷载和试件标矩段的初始横截面面积之比； ${\varepsilon _{{\rm{e}}}}$为工程应变，即试件标距段的变形与初始长度之比。

为考虑试验时间内的试件温度变化，试验结束时记录了试件的试验温度。表1给出了预设温度值与试验温度的对比。

图7为预设温度为20 ℃、预设应变率为1 s⁻¹下三个试件的真实应变率和工程应变率时程曲线，分别由真实应变时程和工程应变时程求导得：

图  7  试验应变率的确定

Figure 7.  Determination of test strain rate

如图7(a)所示，随着轴向变形的增大，试件的横截面面积减小，进而导致真实应变率明显减小。恒速率加载模式不能实现真实应变率的恒定，只能保持工程应变率的相对平稳，如图7(b)所示。因材料的拉伸力学性能主要用于描述大变形阶段的贡献，取工程应变率时程曲线中峰值点和破坏点之间的平均值来表征工况，称为试验应变率。表2给出了各工况的试验应变率。

图8为20 ℃、不同试验应变率下的工程应力-应变曲线和真实应力-应变曲线。可以看出：加载初期，工程应力-应变曲线和真实应力-应变曲线类似，均呈现基本线性关系；当曲线接近屈服点时，两者逐渐表现出明显的差别：工程应力随着工程应变的增加先减小后保持一定的应力水平直至断裂，而真实应力-应变曲线则相继进入强化阶段直至试件断裂，部分曲线还在强化前表现出软化现象。

图  8  20 ℃、不同试验应变率下应力-应变曲线(虚线为工程应力-应变曲线，实线为真实应力-真实应变曲线)

Figure 8.  Strain-stress curves for different test strain rates at 20 ℃ (dashed lines represent engineering stress-strain curves, solid lines represent true stress-strain curves)

图9(a)~图9(d)给出了−33 ℃、20 ℃、47 ℃和54 ℃、所有试验应变率下的真实应力-应变曲线。所有曲线均经历了初始线性、屈服和强化阶段。对于应变率为1 s⁻¹~100 s⁻¹、试验温度为−20 ℃~35 ℃的部分工况，曲线还会在屈服后进入应变软化阶段。

各温度下，随着应变率提高，相同应变水平的应力都有提高，且该效应在1 s⁻¹~100 s⁻¹内更明显。此外，温度低于20 ℃时，应变率对线性阶段斜率的影响较小。而温度高于20 ℃时，随应变率的提高，线性段的斜率显著增大。

图  9  不同试验温度和试验应变率下真实应力-应变曲线

Figure 9.  True strain-stress curves at different test temperatures and test strain rates

为量化分析温度和应变率对离子型中间膜力学性能的影响，根据真实应力-应变曲线确定了所有工况的关键力学参数，包括初始弹性模量、名义屈服强度和名义屈服应变，以及抗拉强度和相对应的极限应变。其中，初始弹性模量取应变范围为0%~1%的割线模量，名义屈服强度取残余应变为0%~2%时的应力值，屈服应变为相应的应变值。需要指明的是，在讨论关键力学参数的应变率效应时，为描述一致，采用预设应变率，但该值与实际的应变率有一定差别。

图10给出了初始弹性模量随试验温度和预设应变率的变化曲线。结果表明：离子型中间膜的初始弹性模量随温度的升高而减小，且减小幅度与温度区间密切相关：当温度低于20 ℃时，初始弹性模量模量随温度的变化不明显，基本在1000 MPa左右；当温度高于20 ℃时， 初始弹性模量随温度的升高而显著减小，40 ℃时，初始弹性模量减小至500 MPa左右，仅为20 ℃时的一半。

图  10  初始弹性模量随试验温度和预设应变率变化的对比

Figure 10.  Comparison of initial modulus with test temperature and preset strain rate

此外，应变率对初始弹性模量的影响也与温度区间相关：当温度低于20 ℃时，应变率对初始弹性模量的影响较小且没有明显规律，而当温度高于20 ℃时，初始弹性模量随应变率的提高而略有增大。

图11给出了屈服强度随试验温度和预设应变率的变化，结果表明：无论在何种应变率下，离子型中间膜的屈服强度均随温度的升高而减小；而随应变率的提高，屈服强度表现出增大的趋势。

图  11  屈服强度随试验温度和预设应变率变化的对比

Figure 11.  Comparison of yield stress with test temperature and preset strain rate

图12给出了屈服应变随试验温度和预设应变率的变化。当温度低于20 ℃时，屈服应变随温度的升高而变化较小；当温度高于20 ℃时，屈服应变随温度的升高而增大。预设应变率100 s⁻¹下，当温度由47 ℃升高至54 ℃时，屈服应变反而减小。这与试验曲线特征和屈服应变的确定方法有关，如图9所示，47 ℃的曲线有明显的线弹性段，而54 ℃的曲线更趋于超弹性特性。

图12同时表明应变率对屈服应变的影响规律并不明显。

图  12  屈服应变随试验温度和预设应变率变化的对比

Figure 12.  Comparison of yield strain with test temperature and preset strain rate

图13给出了抗拉强度随试验温度和预设应变率的变化，结果表明：抗拉强度随温度升高而显著减小。但应变率的影响并不显著，这可能是不同工况下极限应变的差异导致的，如图14所示，极限应变随应变率的提高而减小。而在相同的应变水平下，应力均随应变率的提高而增大，如图9所示。

图  13  抗拉强度随试验温度和预设应变率变化的对比

Figure 13.  Comparison of failure stress with test temperature and preset strain rate

图  14  极限应变随试验温度和预设应变率变化的对比

Figure 14.  Comparison of failure strain with test temperature and preset strain rate

G’SELL模型是适用于固相聚合物的经验型动态本构关系^[14]。该模型基于对固相聚合物试验数据的分析，将其力学响应归结为粘弹性和塑性共同作用的结果，且这两个因素是相互独立的。此外，G’SELL模型假设材料的粘弹性只受应变影响，而材料的塑性发展是应变和应变率的共同作用的结果。G’SELL本构关系表达为：

式中：比例因子为$K$；粘性参数为$\varpi $；应变强化系数为$h$；应变率系数为$m$；温度系数为$a$，共5个模型参数。基于本文试验数据，标定模型参数列于表3。

为验证表3所给本构模型参数的有效性，图 15 分别对比了本文及已有文献[3, 8]中的试验应力-应变曲线和按表3所给参数预测的本构关系曲线，可以发现吻合良好。

图  15  模型预测结果与试验结果的比较

Figure 15.  Comparison between predicted results and test results

考虑到通用软件的材料模型库没有包括G’SELL模型，为将第3节给出的本构关系应用于数值分析，本文基于LS-DYNA平台开发了材料本构关系的用户子程序。

本文采用基于弹性张量的应力补偿法对主程序的每一个增量步的应力进行更新^[15]，算法原理如图16。首先根据Mises准则判断当前增量步材料是否屈服，若屈服则实际应变增量为${\varDelta }{{\varepsilon }_{t + \varDelta t}} = {\varDelta \varepsilon }_{t + \varDelta t}^{\rm{e}} + {\varDelta \varepsilon }_{t + \varDelta t}^{\rm{p}}$。先假设应变增量都是弹性应变即${\varDelta }{{\varepsilon }_{t + \varDelta t}} = {\varDelta \varepsilon }_{t + \varDelta t}^{\rm{e}}$，由广义胡克定律计算试探应力${\sigma }_{t + \varDelta t}^{{\rm{tr}}} = {\sigma }_t^{{\rm{tr}}} + {{D}^{\rm{e}}}:{\varDelta \varepsilon }_{t + \varDelta t}^{\rm{e}}$，其中${{D}^{\rm{e}}}$是弹性矩阵。然后计算塑性应变增量，并对试探应力进行负补偿，即从试探应力减去塑性应变增量对应的应力增量部分，${\sigma }_{t + \varDelta t}^{{\rm{tr}}} = {\sigma }_t^{{\rm{tr}}} - {{D}^{\rm{e}}}:{\varDelta \varepsilon }_{t + \varDelta t}^{\rm{p}}$。

图  16  应力补偿更新算法

Figure 16.  Stress compensation method

基于应力补偿更新算法，应力更新步骤为：1)计算试探应力；2)根据Mises准则判断是否屈服；3)若屈服，则计算塑性应变增量，然后在试探应力中减去塑性应变增量对应的应力增量部分，对试探应力进行负补偿。具体步骤阐述如下：根据广义胡克定律，各向同性材料的应力-应变关系表示为：

式中，$\lambda $和$\mu $为Lame常数。

试探应力表示为：

根据Mises屈服准则，屈服方程为：

式中：${{{S}}_{{{ij}}}}$为试探应力的偏应力；${\sigma _{\rm{y}}} = f\left( {\varepsilon ,\dot \varepsilon } \right)$为动态屈服应力。若未屈服，以线弹性方法更新应力。若屈服，先求解等效塑性应变增量$\Delta \varepsilon _{}^{\rm{p}}$：

式中：${S^{{\rm{tr}}}}$为试探应力的等效偏应力；$H \!=\! \partial {\sigma _{\rm{y}}}/\partial \Delta \varepsilon _{}^{\rm{p}}$为增量步开始时的塑性硬化率。根据${J_2}$流动法则，塑性应变增量的张量与偏应力张量${{{S}}_{{{ij}}}}$成正比例，得塑性应变增量张量：

由于塑性变形部分体积不变，即 ${\rm{tr}}\left( {{{\varDelta \varepsilon }}_{}^{{{\rm{p}}}}} \right) = 0$，所以应力更新为：

即：

同时，下一增量步的等效塑性应变和等效塑性应变率均可得出。

当LS-DYNA软件材料库的材料模型不能满足用户需求时，可采用用户自定义子程序自定义材料的本构，并在K文件中使用关键字*MAT_USER_DEFINED_MATERIAL_ MODELS进行调用。LS-DYNA主程序在$t + \varDelta t$时，传递给用户子程序t时刻的应变张量${{{\sigma }}_{{t}}}$，应变增量${{\varDelta }}{{{\varepsilon }}_{{{t + \varDelta t}}}}$，以及状态变量，要求用户子程序给出$t + \varDelta t$时刻的应力增量${{\sigma }}_{{{t + \varDelta t}}}^{}$，即用户子程序需要对本构方程进行求解，完成应力张量的更新。用户子程序的计算流程如图17所示。

根据上述应力更新算法，按照LS-DYNA用户子程序接口规范，基于FORTRAN语言编写离子型中间膜本构关系的用户子程序。

图  17  用户子程序流程图

Figure 17.  Flowchart of user subroutine

首先，考量用户子程序的可行性，采用对单个单元在单轴加载下的动态响应数值模拟来进行验证。这是因为材料行为应当与结构的形态无关，且单个单元模型几何简单，加载明确，计算积分点少，节约时间且便于检查子程序合理性。

建立SOLID164单元的立方体单元如图18，边长为1000 mm。固定单元下表面的1、2、3和4四个节点的${\textit{z}}$方向自由度、1节点$x$方向自由度和2节点$y$方向自由度。固定上表面5节点$x$方向自由度，7节点$y$方向自由度。在上表面的5、6、7和8四个节点施加${\textit{z}}$方向速度荷载V=0.1 m/s、1 m/s和300 m/s。材料参数取值同表3。

图  18  单元模拟的边界和荷载条件　/mm

Figure 18.  Boundary and loading condition of element simulation

图19为20 ℃和0.1 s⁻¹、1 s⁻¹以及300 s⁻¹下单元数值模拟结果和本文试验结果的对比，由图可见两个结果吻合非常好，进而验证了子程序开发的可行性和正确性。

图  19  单元数值试验结果和本文试验结果的对比

Figure 19.  Comparison between testing results and element simulation results

其次，考量该用户子程序在构件模拟应用中的可行性和准确性。对本文进行的单调拉伸试验，调用用户子程序进行数值模拟。 试件的左端固支约束，右端为一沿试件轴向施加的速度荷载3 m/s，即标距段的预设应变率为300 s⁻¹，荷载施加时间为0.006 s。采用实体单元SOLID164建立有限元模型，网格划分时，沿宽度方向将整个试件进行32等分，沿长度方向的划分为：两端夹持段20等分、中间标距段200等分、圆弧段25等分、厚度方向3等分。建立的有限元模型及网格划分如图20。

图  20  有限元模型及网格划分

Figure 20.  Finite element model

单调拉伸数值模拟与试验力-位移曲线对比见图21，可见对于荷载的预测较准确。但数值模拟对于荷载峰值点时间的预测晚于试验值，这可能是数值和试验中加载位置的差异导致的，数值模拟中约束和加载均位于端部，试件的夹持段也参与了变形，而试验中仅夹具中间的试件参与变形。

图  21  用户子程序数值模拟与本文试验测得的力位移曲线的对比

Figure 21.  Comparison of force-displacement curve of the tensile tests and simulation

本文对离子型中间膜材料进行了−40 ℃ ~60 ℃和1 s⁻¹~800 s⁻¹下的动态拉伸试验，分析了温度和应变率的影响，基于G’SELL模型标定了模型参数，开发了用户子程序，主要结论如下：

(1)离子型中间膜表现出显著的温度效应。随温度的升高，材料的初始弹性模量和强度均明显减小。当温度从20 ℃升高至40 ℃，初始弹性模量大致减小为原来的一半。

(2)离子型中间膜表现出一定的应变率效应，在各试验温度下，相同应变水平处的应力均随应变率增加而增大，当材料进入强化阶段后，应变率效应更显著。因此，模拟夹层玻璃破碎后力学行为时，有必要考虑应变率效应。

(3)基于G’SELL模型，标定了模型参数，并通过对比模型预测曲线与已有研究试验曲线进行了验证。

(4)基于应力补偿算法开发了用户子程序，将本文所提出本构模型引入LS-DYNA软件，通过数值结果与试验结果的对比，证明用户子程序的有效性。