一种修正的GFRP本构模型及其双层板结构抗冲击特性预测

张 萌，杨 扬，张岳青，柳 柳，王 慧

(1.长安大学 理学院，西安 710064；2.西北工业大学 航空学院，西安 710072；3.中国船舶重工集团公司第七〇五研究所，西安 710075)

玻璃纤维增强聚合物基材料(Glass Fiber Reinforced Polymer，GFRP)是一种典型的纤维增强型复合材料，具有轻质高强度、耐腐蚀、可设计性好等诸多优点，在航空航天飞行器减重设计、船舶外壳制造等领域具有重要应用。在实际使用中，GFRP会不可避免地遭受外来物的撞击，因此对材料的动态响应计算模型及抗冲击特性进行深入研究具有重要意义[1-5]。

目前，对于各向同性材料冲击动力学响应计算方法的研究已较为成熟，简单来说就是分别利用状态方程和本构关系计算材料的球应力和偏应力。但是GFRP是典型的各向异性复合材料，其在冲击载荷作用下的力学响应十分复杂，状态方程和本构关系在计算过程中高度耦合，通常情况下球应力在各个方向上产生的应力不同，可能会带来形状改变；同时，偏应力也会对体应变产生一定的影响[6]。目前已有的针对GFRP动态材料模型的研究大致可以分为两类：其一是宏观模型，这类模型将纤维-基体看作统一的整体，忽略二者之间的相互作用，赋予统一的材料参数，计算效率较高，但是各个应力分别计算，忽略了各向异性材料中固有的耦合效应(如Autodyn商用软件中提供的GFRP材料模型)；其二是微观模型，这类模型将纤维和基体分别建模，更注重对材料损伤机理的细节剖析，此类模型考虑材料的各向异性，但通常情况下模型计算量极大，效率较低。本文主要研究GFRP在冲击载荷下的宏观力学响应，因此重点对宏观动态材料模型进行改进。

宏观各向异性材料模型的核心问题是如何对材料内部复杂的力学响应进行解耦。Anderson等[7]提出了一种适用于求解各向异性材料中耦合响应问题的计算方法，继而Hayhurst等[8-10]运用此方法提出了一种适用于各向异性材料冲击仿真的材料模型框架，并在Kevlar-Epoxy、CFRP等经典复合材料的高速冲击数值仿真中予以了应用，取得了较好的仿真结果。本文也将基于上述模型框架，对Autodyn软件中现有的GFRP材料模型中的状态方程和强化模型进行修正，继而推导并获取相关材料参数，最终实现对GFRP-Al双层板结构抗冲击特性的仿真预测，为工程设计领域提供一定的理论和技术支撑。

1 修正的GFRP动态本构模型

对于GFRP复合材料来说，准确描述其在冲击载荷作用下的力学响应特性，有两点至关重要：① 固体材料在高速或超高速冲击载荷作用下的瞬态变形极大，局部作用区域内将呈现出明显的类流体行为，因此状态方程的准确描述与计算尤为重要。② GFRP是典型的塑性基复合材料，材料在屈服与后屈服阶段力学性能的准确刻画，对其抗冲击特性的预测也是至关重要的。Autodyn软件中现有的GFRP材料模型(记为GFRP-old)采用的是Puff状态方程和Von Mises屈服准则，但是首先Puff状态方程表征的是体应变与球应力的单一关系，忽略了各向异性材料中偏应变对球应力的影响；其次，Von Mises屈服准则多用于各向同性材料，它表征的实际上是单位体积形状改变的弹性位能，并没有完全体现出各向异性材料中各个方向力学响应的差异性。本文将基于文献[7]中的解耦方法，重点对GFRP-old模型中的状态方程和屈服模型进行修正。

1.1 状态方程的修正

由于正应变为体应变ε0和偏应变εd之和，此时各向异性材料基本应力-应变关系可记为：

(1)

由式(1)得出σii(i=1,2,3)的表达式，同时球应力P为平均正应力(σ11+σ22+σ33)/3的相反数，由此可推导出修正的状态方程为：

(2)

式中：P为球应力，P(ε0)为普通的状态方程(此处为了与GFRP-old模型进行对比，选为如下式(3)所示的Puff状态方程)、Cij为刚度系数。由式(2)可以清晰地看出，偏应变εd对球应力的计算结果亦具有重要影响，可以看作是对传统状态方程的修正。

(3)

式中：A1、A2、A3和T1、T2为体积常数、Γ0为Gruneisen系数、H为膨胀系数、e为物质的比内能、es为升华能。

1.2 强化模型的修正

GFRP是典型的塑性基复合材料，随着冲击变形过程的发展，材料达到屈服应力后会进入后屈服响应阶段。首先，确定屈服准则：Chen等[11]提出了一种普适的二次屈服函数，该屈服函数充分考虑了各向应力的差异性，并且证明了Tsai-Hill、Von Mises等常用屈服准则均为其特例，因此，本文也选用此屈服函数。其次，后屈服响应的描述：材料进入屈服状态后，由于内部纤维和基体的相互作用，会经历一个后屈服的强化阶段，之后才会达到最终失效，已有的研究中已证明后屈服阶段对于力学响应的准确刻画至关重要。在文献[11]中二次屈服函数的基础上，可通过材料的等效应力和等效塑性应变曲线表征材料的后屈服阶段[12]：

(4)

至此，GFRP材料修正的动态本构模型已经建立完成，将其记为GFRP-new模型，接下来将利用典型材料性能试验获取上述模型中的材料参数。

1.3 模型参数

文献[13]针对一种GFRP试样完成了材料性能测试实验和两种厚度(1.95 mm和3.9 mm)的GFRP薄板冲击实验。利用文献[13]基本材料性能参数推导本文修正模型中需要的模型参数[14]，如表1～3所示。其中，本文中约定11方向为面外方向，22和33方向为面内方向。刚度矩阵通过柔度矩阵求逆计算得到；为了与GFRP-old模型进行对比，Puff状态方程参数直接选用Autodyn中现有数据；屈服模型参数可参照文献[11-12]及式(4)计算得出；失效模型参数直接由文献[13]得出。

表1 GFRP刚度系数及状态方程参数

表2 GFRP的屈服模型参数

表3 GFRP的失效模型参数

2 修正模型合理性验证

为了验证本文修正模型的正确性与合理性，利用Autodyn软件建立球形弹丸-薄板冲击二维轴对称计算模型，并与已有结果进行对比分析。

在下面的所有验证算例中，弹丸选用Autodyn材料库中的Steel1006钢，直径5 mm；GFRP分别选用Autodyn材料库中的GLASS-EPXY材料(记为GFRP-old模型)和前文修正的动态本构模型及参数(记为GFRP-new模型)进行计算。计算时选用SPH求解器模拟冲击中的动态大变形现象，粒子间距设置为0.1 mm，采用变光滑长度计算。

2.1 冲击实验结果对比分析

参照文献[13]中的两种实验工况，弹丸及薄板的计算参数如表4所示。图1(a)～(b)分别给出两种工况下的计算结果与实验结果的对比图，从图中可以看出：① 相比于Autodyn中自带的GFRP-old模型，本文修正模型的计算结果与实验结果更为接近，具体来说：1.95 mm工况时均方误差由4 255.2减小为1 506.7，降低了64.5%；3.9 mm工况时均方误差由6 296.8减小为2 384.5，降低了62.1%，验证了修正模型的有效性。② 单独来看GFRP-new模型，发现随着弹丸冲击速度的增大，计算结果与实验值的误差逐渐减小。这是由于当冲击速度逐渐增大时，弹丸撞击点局部区域变形愈发剧烈，材料内部的类流体行为将更为显著，这就使得状态方程的准确计算变得尤为重要。本文修正模型正是在常规状态方程中引入了偏应变的修正项，因此在高速冲击下的计算精度得到了提升，这也进一步验证了修正模型的合理性。

(a)1.95 mm

表4 两种实验工况

2.2 碎片云形貌对比分析

下图2给出了采用两种GFRP材料模型计算的薄板背部碎片云形状，弹丸速度1 000 m/s，方形薄板厚度为5 mm，边长依然为100 mm，整个模型粒子总数为51 916个。可以看出：① 采用两种材料模型计算所得的薄板背部碎片云扩展区域尺寸基本相同，分布直径D约为30 mm，这说明两种材料模型在描述GFRP抗冲击响应的本质特征方面基本一致；② 从碎片云形貌来看，利用本文提出的改进的材料模型计算的碎片云细节更加丰富，这主要是由于球应力在碎片云扩展和成形的计算中具有重要作用，而本文通过对状态方程进行修正，获得了更为合理的球应力计算结果。

(a)GFRP-old模型

2.3 塑性波传播对比分析

图3给出了采用两种材料模型计算所得的对应时刻的塑性波传播过程，其中深黑色表示材料处于弹性状态，浅灰色表示材料处于塑性状态。GFRP-old模型中采用了von Mises屈服准则，在弹丸侵彻过程中，内部材料状态呈现了如图3(a)中箭头所示的繁杂的分布。本文选用了更为普适的二次屈服函数，计算结果显示材料内部的塑性波沿弹靶接触点逐渐向薄板四周均匀传播，计算结果更为可信。

(a)GFRP-old模型

3 GFRP-Al双层防护结构抗冲击特性预测

3.1 不同弹丸冲击下破坏模式分析

上节通过对冲击实验结果、碎片云形貌和塑性波传播三个方面的对比，说明了改进的本构模型的合理性和有效性。本节将利用改进的动态本构模型对GFRP薄板的抗冲击特性进行仿真预测。利用Autodyn软件建立GFRP-Al双层板冲击计算模型：弹丸分别选用Steel1006钢球和Al球，直径5 mm，冲击速度为1 000 m/s；前板选用GFRP板，后板选用Al板，二者厚度均为5 mm，边长均为100 mm，板间距50 mm。计算时仍选用SPH求解器，粒子间距设置为0.1 mm，整个模型粒子总数为111 916个。

图4(a)、(c)依次给出了GFRP-Al双层结构在钢球冲击作用下3个时刻的仿真计算结果，可以看出：① 当弹丸穿过第一层GFRP薄板后，会嵌于碎片云中，并随碎片云一起冲击后板。② 铝板此时抵挡了所有GFRP碎片，并发生了局部的塑性变形；而钢弹丸此时即将撞击铝板，剩余速度约为599 m/s。③ 钢弹丸二次撞击铝板，并造成了铝板穿孔，飞行速度进一步降低为约243 m/s。

图4(d)、(f)对应给出了GFRP-Al双层结构在铝球冲击作用下3个典型时刻的仿真计算结果。可以看出：① 对比图4(a)和(d)可知，相同冲击速度下，铝球撞击GFRP薄板后形成的碎片云扩散区域更小，分布直径约为21 mm，但是两种弹丸的穿孔直径基本相当，也就是说二者的碎片含量基本相当，因此铝球弹丸撞击后形成的碎片云团密度更大。② 在铝球即将接触后板时，如图4(g)、(h)所示，后板整体变形较钢球撞击情形更大，同时会形成局部凸起。这是因为与钢弹丸相比，一方面铝球质量较小，撞击前板后剩余能量较小，穿甲性能较弱，因此导致的整体变形更大；另一方面，铝球撞击面自身变形较大，凹陷较为明显，对碎片颗粒产生一定的聚集效应，从而对后板产生较强的局部作用形成凸起。③ 当铝弹丸二次撞击后板时，会与碎片颗粒一起对后板局部产生集中作用，在后板中心形成塞块沿冲击轴线飞出；而钢弹丸对后板的二次撞击，由于其剩余速度较铝弹丸更大，冲击能量更大，更易形成绝热剪切，产生两个碎块向两侧斜向方向飞出，如图4(i)、(j)。

(a)0.04 ms

3.2 不同冲击速度下破坏模式分析

为了进一步探究弹丸以不同冲击速度撞击作用下CFRP-Al双层结构的损伤破坏模式，建立同3.1节相同的仿真模型，其中弹丸选用Steel1006钢球，冲击速度分别为200 m/s、500 m/s和1 000 m/s，其余参数保持不变。计算结果如图5所示，可以看出：① 随着冲击速度的升高，GFRP薄板均发生了穿孔，其破坏形貌并无明显差异，但是速度较高时，冲击能量较大，靶板内部破坏更为严重，在薄板的撞击侧会出现碎片颗粒飞溅。② 随着冲击速度的升高，弹丸穿透首层GFRP薄板后的剩余速度逐渐升高，后板Al板的破坏模式出现差异。当v=200 m/s时，钢弹丸穿透GFRP薄板后的剩余速度为126 m/s，在与Al板和碎片云团相互挤压作用后，速度进一步降为42 m/s，无法穿透第二层Al板；当冲击速度增大为v=500 m/s时，弹丸剩余速度较高，会造成Al板出现局部塞块；而当冲击速度进一步升高为v=1 000 m/s时，后板的破坏形貌再次发生变化，出现如图所示向两侧飞行的剥落碎块(三维模型中应是一个向外侧飞行的圆环)。③ 进一步观察弹丸与Al板撞击的局部图，可以看出当冲击速度较低时(图5(d)、(e))，Al板整体变形较大，此时在弹丸与Al板变形形成的凹坑内聚集了大量的碎片云颗粒，这些颗粒与弹丸的速度方向一致，它们共同作用于第二层Al板，造成Al板局部出现冲塞破坏；而当冲击速度增大为1 000 m/s时(图5(f))，弹丸穿透GFRP薄板后的剩余速度较高，达到了707 m/s，剩余动能较大，此时弹丸会迅速排开碎片云颗粒直接撞击第二层Al板；同时由于初始撞击速度较大，弹丸与GFRP薄板作用后，其形状由球形近似变为椭球形，撞击后板时接触面更大。这种双重的影响会使得Al板撞击轴线附近易形成绝热剪切，产生如图所示向两侧飞行的剥落碎块(三维模型中为一个剥落的圆环)。

(a)v=200 m/s

3.3 抗冲击特性规律分析

为了全面评估GFRP-Al双层结构的抗冲击特性，分别改变弹丸冲击速度、弹丸直径和靶板厚度，对弹丸依次穿透双层靶板结构后的剩余速度进行规律性分析。该算例中，弹丸均选用Steel1006钢球，第一层靶板为GFRP薄板，第二层靶板为Al薄板。具体设置如表5，计算结果分别如图6～8所示。

图6 弹丸剩余速度随冲击速度变化图

表5 计算参数

从计算结果可以看出：① 弹丸穿透第一层GFRP板的临界穿透速度约为100 m/s，继续穿透第二层Al板的临界穿透速度约为600 m/s，也就是说当弹丸初始冲击速度小于600 m/s时，它们将无法穿透第二层靶板，而是出现回弹，在图中将其剩余速度记为0；② 弹丸穿透第一层GFRP板的临界穿透直径约为1 mm，继续穿透第二层Al板的临界穿透直径约为4 mm；③ 当GFRP板厚度不超过10 mm时，1 000 m/s的钢弹丸可以穿透所有双层板结构。④ 无论第一层还是第二层靶板，一旦弹丸穿透，弹丸剩余速度随冲击速度基本呈现线性变化趋势，随弹丸直径和GFRP板厚度呈现非线性变化趋势。

图7 弹丸剩余速度随弹丸直径变化图

图8 弹丸剩余速度随GFRP薄板厚度变化图

4 结 论

本文针对GFRP各向异性材料，采用经典的Anderson解耦方法对Autodyn中传统的GFRP材料模型中的状态方程和强化模型进行了修正。与已有实验结果相比，修正模型的计算均方误差降低了60%、碎片云细节更为丰富、塑性波传播计算结果更为合理。进一步，将改进模型应用于GFRP-Al双层板结构的破坏机理和抗冲击特性分析，结果显示后板的破坏模式与弹丸撞击的初始能量及穿透首层GFRP薄板后形成的碎片云形貌密切相关，当初始动能较小时，后板撞击点附近形成塞块沿轴线飞出；而当初始动能较大时，在后板撞击点附近易发生绝热剪切，出现向两侧飞行的剥落碎块(三维模型中为一个剥落的圆环)。