三维正交机织复合材料在重复低速冲击下的有限元仿真

Finite Element Simulation of 3D Orthogonal Woven Composites under Repeated Low-Velocity Impact

摘要: 为了建立碳纤维三维正交机织复合材料(3DOWCs)在重复低速冲击下的有限元仿真，本工作将仿真与实验进行比较，以验证仿真的有效性。通过在自制的3D机织机上制造了3DOWCs板件，并对其进行50 J重复低速冲击试验，获得断裂形貌与力学曲线。使用多步骤方法模拟重复冲击的试验，冲击过程中引入重启动技术与人工阻尼。建立了多尺度的有限元模型，并将增强丝束材料属性与损伤模式集成至ABAQUS/Explicit用户定义材料子程序(VUMAT)中，3DOWCs损伤模型包含基于位移的失效模型与三维Hashin失效准则，以模拟3DOWCs在重复低速冲击载荷的损伤演化。使用对称边界条件，建立了1/4有限元模型，以减少计算成本。与实验结果相比，首次冲击的力学曲线、三次冲击的峰值力、损伤分布的连续性验证了该模型的有效性。

Abstract: In order to establish the finite element simulation of carbon fiber 3D orthogonal woven composites (3DOWCs) under repeated low-velocity impact, this work compares the simulation with the exper-iment to verify the effectiveness of the simulation. 3DOWCs plates were manufactured on a home-made 3D woven machine and subjected to 50 J repeated low-velocity impact tests, and the fracture morphology and mechanical curves were obtained. The multi-step method is used to simu-late the repeated impact test, and the restart technology and artificial damping are introduced in the impact process. A multi-scale finite element model is established, and the material properties and damage modes of the reinforced tow are integrated into the ABAQUS/Explicit user defined material subroutine (VUMAT). The 3DOWCs damage model includes a displacement-based failure model and a 3D Hashin failure criterion to simulate the damage evolution of 3DOWCs under the re-peated low-velocity impact load. Using symmetric boundary conditions, a 1/4 finite element model is established to reduce the calculation cost. Compared with the experimental results, the validity of the model is verified by the mechanical curve of the first impact, the peak force of the third impact and the continuity of the damage distribution.

1. 引言

三维正交机织复合材料(3DOWCs)与传统二维(2D)层合复合材料相比，由于其Z向的捆绑纱存在，抑制了分层的发展与裂纹的传播，已广泛运用于民用运动防护装备、军事领域承重部件。出于设计产品与优化的目的，3DOWCs其结构的动态响应与损伤发展应得到充分研究 [1] [2] [3] 。因此，3DOWCs在低速冲击载荷下的力学行为研究变得有意义。近些年，研究者对3DOWCs在低速冲击载荷下开展了研究 [4] [5] 。但在生活、制造与日常维护中，结构部件往往承受的冲击是重复或多次的，而不是仅受一次单独的冲击 [6] 。复合材料受到低能级重复冲击会因损伤累计导致分层、树脂脱粘和纱线断裂等损伤，而受到低能级单一冲击可能不会有损伤 [7] 。因此，研究复合材料结构在重复冲击下的损伤发展具有重要意义。然而目前复合材料冲击响应与损伤累计研究方法通过实验进行 [8] [9] [10] ，但为了更好理解复合材料在重复冲击下的损伤累计与发展，还需数值模拟的帮助 [11] 。目前研究者对与重复冲击建模仿真主要通过LS Dyna 3D或Abaqus/explicit软件进行 [12] [13] [14] 。Tian和Zhou [12] 建立了一种新的三维渐进退化有限元模型，并将其编写为LSDyna用户子程序供GLARE冲击模拟调用，以研究GLARE在单次与多次冲击下的损伤机理。在最近Zhou等人的工作中 [13] ，将三维Hashin失效准则、连续损伤模型和基于等效位移的损伤演化模型进行编码载入用户定义的子程序(VUMAT)，基于CZM双线性牵引分离法的界面单元被用于模拟复合材料层合板的分层。在Reazasefat等人的工作中 [14] ，Puck模型在VUMAT中得以编码，采用多步骤形式模拟重复冲击。然而，据作者了解，现在对于重复冲击的仿真还仅限于层合板，这项工作将对3DOWCs开发一种重复低速冲击仿真的方法，使用多尺度计算方法并在ABAQUS/Explicit的用户定义材料子程序(VUMAT)进行实现，仿真结果将与实验原始数据进行比对。

2. 试样制备与试验过程

2.1. 三维正交机织复合材料的制备

2.1.1. 三维正交纤维预制体织造

三维正交织物结构具有相互垂直排列的经纱和纬纱，捆绑纱与经纱同向，垂直贯穿织物厚度方向，跨越一纬后反向垂直穿过织物厚度，具有相反运动规律的两组捆绑纱以1:1上下交织，实现经纬纱的紧密交织。三维正交机织预制件由自研多剑杆织机平台织造，其经密(包括经纱方向的捆绑纱)为3根/cm，纬密为3根/cm，经纱3层，且纬纱4层。碳纤维型号为12 K东丽T700SC，并且经纱每束12 K，而纬纱每束24 K。

2.1.2. 复合材料制备

本工作通过真空树脂转移法(VARTM)将树脂导入织物中进行固化成型，环氧树脂其型号为法国AXSON公司的EPOLAM 2040，与之搭配使用的固化剂型号为EPOLAM 2042。在室温20℃，树脂与固化剂质量配比为100:32，混合后用搅拌棒充分搅拌5分钟以得到均质的混合物，将其置于真空状态20分钟，以剔除树脂中的气泡，再经过负压(0.9 bar)使环氧树脂缓慢渗透预制体，待完全渗透后，在室温环境中静置24 h达到初步固化后，将其置于75℃恒温箱中16小时，完成复合材料固化。通过金刚石线切割机，将复合材料切割成尺寸为150 mm × 100 mm的试样，长边沿纬纱方向，短边沿经纱方向。将试样切割成尺寸为20 mm × 20 mm，放置坩埚中，通过马弗炉加热到650˚，待坩埚冷却测定纤维体积含量，每种试样重复三次，得到参数如表1所示。

Table 1. 3DOWCs sample parameters

表1. 3DOWCs试样参数

2.2. 重复低速冲击试验

重复低速冲击试验在低速冲击落锤试验装置上进行。每次冲击按ASTM-D7136测试协议下进行，所使用的夹具为测试协议所规定的，试验夹具上有四个夹子，夹子末端有橡胶阻尼器吸收振动，以避免冲击过程中夹具产生振动所带来的影响，测试夹具中心有长宽为100 mm × 75 mm的矩形空槽。落锤试验机上配备光电传感器与电磁开关，用于避免二次冲击对实验造成的影响。冲头采用直径为16 mm的半球形冲头，冲击能量选取50 J进行试验。落锤试验机设置光电传感器和电磁开关用于防二次冲击，冲击头安装有压电传感器用于反馈冲击瞬间的冲击力，并且数据采集系统采样频率为20 kHz，采用高速摄影仪记录冲头位移校准传感器输出冲击力的大小。

3. 数值模型

3.1. 有限元建模

本节建立多尺度有限元模型，模拟3DOWCs在重复低速冲击中的力学响应。如图1所示，为了减小计算成本，建立了1/4纱线尺度的复合材料，材料截面使用对称边界条件进行约束，如图1(a)所示。建立了三个刚性冲击器按照多步骤方式模拟试验过程中的三次冲击。冲击过程引入重启动技术，将前一次的冲击结果作为下一次冲击前的初始状态，此外，为了减小复合材料在冲击后的残余振动加入人工阻尼。冲击能量由预定义初速度与质量进行分配。3DOWCs材料包含树脂基体与三种增强丝束(经纱、纬纱和捆绑纱)。为了降低几何建模复杂性，假设经纱与纬纱平直且垂直，而捆绑纱简化为“方波”形状。按照3DOWCs截面中纱线尺寸，绘制三种纱线截面，通过Solidworks软件拉伸成实体纱线，再导入Abaqus软件中通过布尔求差得到所需树脂，最后将实体纱线与树脂进行布尔求和成完整复合材料板。经过网格划分后，一共394,800个C3D8R单元，其中树脂260,980个。三种增强丝束纱线被认为横向各向同性材料，其纤维取向如图1(b)所示。

Figure 1. (a) Simulation model of repeated low-velocity impact; (b) Schematic diagram of 3DOWCs simulation model

图1. (a) 重复低速冲击仿真模型；(b) 3DOWCs仿真模型示意图

3.2. 材料模型

3.2.1. 树脂材料属性

树脂基体是一种弹塑性材料，其受力过程可分为弹性阶段与塑性阶段。树脂基体弹性阶段的材料属性采用生产厂家提供的参数如表2所示，而树脂塑性阶段理论模型采用J-2各向同性塑性硬化理论，通过对固化后的树脂基体进行准静态压缩试验获得名义应力应变曲线。通过下式转换成真实应力应变，截选塑性阶段的参数输入Abaqus软件中，以代表树脂基体的塑性行为，如表3所示。

${\sigma }_T={\sigma }_N\left(1+{\epsilon }_N\right)$ (1)

${\epsilon }_T=\ln \left(1+{\epsilon }_N\right)$ (2)

其中 ${\sigma }_T$ ， ${\epsilon }_T$ 分别表示真实应力与真实应变，而 ${\sigma }_N$ ， ${\epsilon }_N$ 分别表示名义应力与名义应变。

Table 2. Mechanical parameters of resin matrix

表2. 树脂基体力学参数

Table 3. Plastic parameters of resin matrix

表3. 树脂基体塑性参数

对于环氧树脂具有应变率相关性，可使用Abaqus软件中Cowper-Symonds准则进行如下表示：

${\stackrel{\dot{}}{\epsilon }}^{pl}=d{\left(R-1\right)}^n$ (3)

$R=\frac{\sigma }{{\sigma }_0}$ (4)

其中 ${\stackrel{\dot{}}{\epsilon }}^{pl}$ 为塑性应变率，R是屈服应力比， $\sigma$ 与 ${\sigma }_0$ 分别是动态及静态加载下的屈服应力；n与d为材料常熟，可通过在不同应变率压缩试验下，拟合应力应变曲线获得，本文中环氧树脂基体的材料常数n与d分别为2.9与4880。

3.2.2. 纱线材料属性

增强丝束主要由两部分组成，一部分是纤维，另一部分是树脂，这两部分组成使增强丝束成为横向各向同性的材料。其纵向力学性能主要来源碳纤维束，而横向力学性能主要是环氧树脂基体为主导，其本构模型如下式表示：

$\sigma =\left[\begin{array}{cc}{C}^0& \\ & G^0\end{array}\right]\cdot \epsilon$ (5)

$C^0=\left[\begin{array}{ccc}{E}_{11}\left(1-v_{23}{v}_{32}\right)\Delta & E_{11}\left(v_{21}+v_{31}{v}_{23}\right)\Delta & E_{11}\left(v_{31}+v_{21}{v}_{32}\right)\Delta \\ & E_{22}\left(1-v_{13}{v}_{31}\right)\Delta & E_{22}\left(v_{32}+v_{12}{v}_{31}\right)\Delta \\ & & E_{33}\left(1-v_{12}{v}_{21}\right)\Delta \end{array}\right]$ (6)

$G^0=\left[\begin{array}{ccc}{G}_{11}\Delta & & \\ & G_{22}\Delta & \\ & & G_{33}\Delta \end{array}\right]$ (7)

$\Delta ={\left(1-v_{12}{v}_{21}-v_{23}{v}_{32}-v_{13}{v}_{31}-2v_{21}{v}_{32}{v}_{13}\right)}^{-1}$ (8)

式中 $\sigma$ 与 $\epsilon$ 为应力与应变，而 $C^0$ 与 $G^0$ 为代表性单元(RUC)的弹性矩阵与剪切刚度矩阵。

本工作使用数值计算方法，通过构造四边形代表性单元(RUC)几何模型以代表增强丝束，如图2所示。纤维丝束周围由树脂浸润包裹，碳纤维丝束与环氧树脂基体材料属性由供应商提供。通过Abaqus软件对RUC单元模型设置边界条件并在显式分析中进行载荷求解，计算获得浸润纱线的弹性参数，而强度属性通过Chamis细观力学公式进行计算，碳纤维浸润纱线材料参数如表4所示。

Figure 2. Schematic diagram of RUC model

图2. RUC模型示意图

Table 4. Reinforced tow material parameters

表4. 增强丝束材料参数

3.3. 损伤起始演化

3.3.1. 树脂损伤准则

采用Abaqus软件显式求解器中内嵌的延展破坏准则与剪切破坏准则对树脂损伤演化进行模拟。此准则的原理是通过判断等效塑性应变积累是否达到阈值，开始生成延性或剪切损伤。之后，树脂基体刚度退化依照基于能量耗散的损伤演化法则，假设损伤变量与塑性位移成正比。依照Hooputra提出的简化函数解析式 [15] ，对于延展破坏准则，可用应力三轴度η与等效塑性应变率 ${\stackrel{\dot{}}{\bar{\epsilon }}}_D^{pl}$ 的函数表示等效塑性应变 ${\bar{\epsilon }}_D^{pl}$ ：

${\bar{\epsilon }}_D^{pl}\left(\eta ,{\stackrel{\dot{}}{\bar{\epsilon }}}_D^{pl}\right)=\frac{\left({\epsilon }_D^{+}\sinh \left[k_0\left({\eta }^{-}-\eta \right)\right]+{\epsilon }_D^{-}\sinh \left[k_0\left(\eta -{\eta }^{+}\right)\right]\right)}{\sinh \left[k_0\left({\eta }^{-}-{\eta }^{+}\right)\right]}$ (9)

式中 ${\epsilon }_D^{+}$ 与 ${\epsilon }_D^{-}$ 是材料发生延展损伤时，由轴向拉伸与压缩引起的等效塑性应变； $k_0$ 为材料参数； ${\eta }^{-}=-1/3$ 与 ${\eta }^{+}=1/3$ 是轴向拉伸与压缩的应力三轴度；此破坏准则，基于塑性应变的状态量 ${\omega }_D$ 作为材料延展损伤启始的标准，当 ${\omega }_D$ 为1时，延展破坏开始发生：

${\omega }_D\text{=}{\displaystyle \int \frac{\text{d}{\bar{\epsilon }}^{pl}}{{\bar{\epsilon }}_D^{pl}\left(\eta ,{\stackrel{\dot{}}{\bar{\epsilon }}}^{pl}\right)}}=1$ (10)

而关于剪切破坏准则， ${\bar{\epsilon }}_S^{pl}$ 为材料损伤启始时的等效剪切应变函数，其变量 ${\tau }_S$ 与 ${\stackrel{\dot{}}{\bar{\epsilon }}}^{pl}$ 分别为剪切应力与应变率：

${\bar{\epsilon }}_S^{pl}\left({\tau }_S,{\stackrel{\dot{}}{\bar{\epsilon }}}^{pl}\right)=\frac{{\epsilon }_S^{+}\sinh \left[f_S\left({\tau }_S-{\tau }_S^{-}\right)\right]+{\epsilon }_S^{-}\sinh \left[f_S\left({\tau }_S^{+}-{\tau }_S\right)\right]}{\sinh \left[f_S\left({\tau }_S^{+}-{\tau }_S^{-}\right)\right]}$ (11)

${\tau }_S\text{=}\frac{\left(1-k_S\eta \right){\sigma }_{eq}}{{\tau }_{\max }}$ (12)

式中 $k_S$ 与 $f_S$ 均为材料参数， ${\tau }_S^{+}$ 与 ${\tau }_S^{-}$ 是材料发生剪切损伤时，由轴向拉伸与压缩引起的等效塑性应变，此破坏准则，基于塑性应变的状态量 ${\omega }_S$ 作为材料剪切破坏启始的标准，当 ${\omega }_S$ 为1时，剪切破坏开始发生：

${\omega }_S={\displaystyle \int \frac{\text{d}{\bar{\epsilon }}^{pl}}{{\bar{\epsilon }}_S^{pl}\left(\eta ,{\stackrel{\dot{}}{\bar{\epsilon }}}^{pl}\right)}}=1$ (13)

3.3.2. 增强丝束的损伤准则

增强丝束损伤起始与演化使用三维Hashin准则 [16] [17] ，Hashin准则中包含着纤维方向与树脂方向的拉伸损伤与压缩损伤。通过VUMAT用户子程序建立增强丝束的损伤演化模型，并且当增强丝束的计算积分点满足损伤判据时会进行损伤模式输出，其损伤起始判据列于表5 (下标数字1表示增强丝束的轴向，数字2、3表示增强丝束的横向)。 $F_f^t\left(\sigma \right)$ ， $F_f^c\left(\sigma \right)$ ， $F_m^t\left(\sigma \right)$ 和 $F_m^c\left(\sigma \right)$ 为四种损伤模式失效指标，其数值分别以状态变量SDV18至SDV21进行输出。

Table 5. Damage initiation criterion of reinforcing tow

表5. 增强丝束损伤起始判据

当损伤起始判据成立，增强丝束出现失效情况，其材料刚度将会降低。本文通过引入损伤变量( $d_f$ 与 $d_m$ )来表示纤维与树脂基体的损伤程度，并通过函数 $\alpha =\alpha \left(F\right)$ 表示损伤演化程度：

$\alpha \left(F\right)=\{\begin{array}{c}{d}_{\max },F\ge 1\\ 0,F<1\end{array}$ (14)

对于纤维损伤变量 $d_f$ ，可分为拉伸损伤 $d_f^t$ 主导与压缩损伤 $d_f^c$ 主导，表示如下：

$d_f=\{\begin{array}{c}{d}_f^t=\alpha \left(F_f^t\right),F_f^t>F_f^c\\ d_f^c=\alpha \left(F_f^c\right),F_f^t\le F_f^c\end{array}$ (15)

而对于树脂基体损伤变量 $d_m$ ，可分为拉伸损伤 $d_m^t$ 主导与压缩损伤 $d_m^c$ 主导，表示如下：

$d_m=\{\begin{array}{c}{d}_m^t=\alpha \left(F_m^t\right),F_m^t>F_m^c\\ d_m^c=\alpha \left(F_m^c\right),F_m^t\le F_m^c\end{array}$ (16)

基于Murakami-Ohno损伤理论 [16] [18] ，构造出对于增强丝束发生损伤后的单元等效应力 $\tilde{\sigma }$ ：

$\tilde{\sigma }=\left[\begin{array}{cc}C\left(D\right)& \\ & G\left(D\right)\end{array}\right]\cdot \epsilon$ (17)

将损伤变量 $d_f$ 与 $d_m$ 引入至原有的增强丝束刚度矩阵中，得出损伤后具有刚度退化矩阵，其各分量如下表示：

$C_{11}=\left(1-d_f\right)C_{11}^0$ (18)

$C_{22}=\left(1-d_f\right)\left(1-d_m\right)C_{22}^0$ (19)

$C_{12}=\left(1-d_f\right)\left(1-d_m\right)C_{12}^0$ (20)

$C_{13}=\left(1-d_f\right)\left(1-d_m\right)C_{13}^0$ (21)

$C_{23}=\left(1-d_f\right)\left(1-d_m\right)C_{23}^0$ (22)

$G_{12}=\left(1-d_f\right)\left(1-s_{mt}{d}_m^t\right)\left(1-s_{mc}{d}_m^c\right)G_{12}^0$ (23)

$G_{23}=\left(1-d_f\right)\left(1-s_{mt}{d}_m^t\right)\left(1-s_{mc}{d}_m^c\right)G_{23}^0$ (24)

$G_{13}=\left(1-d_f\right)\left(1-s_{mt}{d}_m^t\right)\left(1-s_{mc}{d}_m^c\right)G_{13}^0$ (25)

为了控制树脂基体因拉伸损伤与压缩损伤而导致增强丝束剪切刚度的退化，引入了剪切变量因子 $s_{mt}$ , $s_{mc}$ ，本文中设定 $s_{mt}$ 与 $s_{mc}$ 分别为0.9与0.5。当增强丝束出现损伤时，刚度退化模型的计算导致材料刚度降低，并且单元应力同步更新，随后当前增量步节点计算完成后会返回损伤变量与应力至主程序进行迭代运算。

4. 仿真结果有效性分析

4.1. 仿真与实验数据对比

如图3(a)所示为实验与仿真首次冲击力–时间曲线比较，从图中可以看出曲线趋势总体相似。冲击力首先随时间逐渐上升，到达峰值后开始震荡一段时间，接着逐渐下降至零。图3(b)为三次冲击峰值力实验与仿真对比，可以发现仿真峰值力与实验峰值力数值相近，并且首次冲击力均小于后两次冲击峰值力。但其中首次冲击仿真冲击接触时间与实验有一定误差和仿真峰值力相较于实验峰值力略高。这可归结于仿真模型的简化，同时在建模过程中还忽略了试样制备时的系统误差，比如固化时会存在些许气泡等。

Figure 3. (a) Comparison of first impact value and experimental force-time; (b) Comparison between numerical value and experimental peak force

图3. (a) 首次冲击数值与实验力–时间比较；(b) 数值与实验峰值力比较

4.2. 仿真与实验损伤形貌对比

如图4所示为实验形貌与仿真树脂损伤变量图冲击正面的对比，可以看出首次冲击后，冲击正面存在树脂破损区域，而仿真中恰好也存在局部树脂损伤区域。随着冲击次数的增加，3DOWCs遭受第3次冲击后，其冲击正面树脂损伤区域扩大并且损伤加重了，观察仿真中的树脂损伤区域也明显的扩大了，如图4(b)所示。这表明对于树脂损伤的预测，仿真是可行的。

Figure 4. (a) Comparison of the first impact frontal experiment and simulation topography; (b) Comparison of the third impact frontal test and simulation topography

图4. (a) 首次冲击正面实验与仿真形貌比较；(b) 第3次冲击正面实验与仿真形貌比较比较

图5(a)为实验形貌与仿真树脂损伤变量图冲击背面的对比，实验背面形貌发现存在沿纬纱方向的树脂裂纹，而对应仿真冲击背面观察到沿纬纱方向的树脂损伤区域是明显大于沿经纱方向。而图5(b)为第3次冲击后实验背面形貌与增强丝束中的树脂拉伸损伤变量图，当损伤达到破裂能时，网格会被删除。观察第3次冲击背面形貌有捆绑纱的拔出与沿纬纱的断裂纹，而在仿真中也存在捆版纱的损伤与沿纬纱方向的纱线翘起，此外，这也再次证明该仿真模型可以预测3DOWCs的损伤。

Figure 5. (a) Comparison of the first impact back experiment and simulation topography; (b) Comparison of yarn damage in the third impact

图5. (a) 首次冲击背面实验与仿真形貌比较；(b) 第3次冲击纱线损伤比较

5. 结论

本工作制作了碳纤维三维正交机织复合材料板件，并进行了重复低速冲击实验。在ABAQUS/Explicit使用用户材料子程序VUMAT对3DOWCs重复低速冲击进行仿真模拟。本文旨在为需通过有限元分析去理解3DOWCs在重复低速冲击的动态力学响应与损伤机理提供仿真的思路与方法。通过建立多尺度模型，使用基于位移的失效准则模拟树脂的损伤演化，三维Hashin准则模拟增强丝束的损伤演化。仿真结果与实验数据和断裂的形貌进行比较，证明了该仿真方法用于建立3DOWCs重复低速冲击仿真是有效的。

致谢

感谢浙江理工大学为本文的工作提供了有效帮助，同时还感谢编辑老师与评审老师为本文的工作提出宝贵批评与建议。