专业论文｜一种基于数值模拟的复合材料损伤参数标定方法

曹鹏  张伟  原中晋  雷在林

达索析统（上海）信息技术有限公司北京分公司

纤维增强复合材料性能仿真中，Hashin渐进损伤准则被广泛应用于含损伤结构强度的预测。对于不同材质和结构的复合材料，Hashin损伤参数是不同的。准确获得这些参数，进而得到有意义的计算结果，对于前期仿真十分重要。对于多层复合材料结构，直接测试获得仿真需要的参数可能并不太容易。本文尝试一种方法，利用有限元软件Abaqus和参数优化工具Isight，通过结合测试数据的拟合，确定Hashin损伤参数。

Hashin损伤模型预测弹性-脆性材料的各向异性损伤。它主要应用于纤维增强复合材料，并考虑四种不同的破坏模式：（1）纤维张力破坏；（2）纤维压缩破坏；（3）拉伸和剪切下基体开裂；（4）压缩和剪切下基体压溃。

有限元软件Abaqus包含Hashin损伤模型，属于材料本构模型的一部分。其定义包括损伤起始（Damage Initiation）和损伤扩展（Damage Evolution）两部分。损伤起始是指在应力张量满足特定准则时，开始发生损伤，此时结构的承载能力会下降。损伤扩展描述在损伤发生后，结构刚度的下降规律。Abaqus的Hashin损伤起始准则如下所列。定义Hashin损伤起始准则需要输入(1)-(4)式中的六个参数：纵向拉伸强度（X^T）、纵向压缩强度（X^C）、横向拉伸强度（Y^T）、横向压缩强度（Y^C）、纵向剪切强度（S^L）、横向剪切强度（S^T），分别对应四种不同失效模式。^[1]

Abaqus的Hashin损伤扩展模型需要指定四种失效模式的断裂能（也称为能量释放率）：纵向（纤维）拉伸断裂能（G_ft）、纵向（纤维）压缩断裂能（G_fc）、横向（基体）拉伸断裂能（G_mt）、横向（基体）压缩断裂能（G_mc）。Abaqus仿真模型中材料属性定义在单个铺层上，因此Hashin损伤参数是针对单个铺层的，但由于多层复合材料中各层铺设角不同，中间铺层的性能会受周边铺层的影响，即原位效应。因此输入的单层参数也需要包含原位效应的影响。

原位效应（In-situ effect）最早发现于Parvizi的试验中^[2]，其含义是当复合材料是交叉铺层时（例如[0°，90°，0°]），中间铺层的横向拉伸和纵向剪切强度会相对于单独铺层增大，这是由于中间层周围存在取向不同的铺层，这会对中间层的基体拉伸破坏有一定抑制作用^[3]，如图1所示。由于原位效应的存在，在仿真模型中，基体拉伸相关的Hashin损伤参数不能直接输入单层试样测试的结果，而需要包含原位效应的修正，才能得到更合理的结果。在Hashin损伤模型的10个参数中，受原位效应影响的参数主要为Y^T、S^L和G_mt。

通过直接测试获得包含原位效应修正的参数是比较困难的，这是由于如果测试单层试样，则一定无法包含相邻层的影响，而如果测试多层试样，得到的是整体性能而不是单层性能，也无法直接用于仿真。因此本文尝试一种结合仿真和测试的拟合方法，通过Abaqus + Isight联合仿真，结合测试数据标定仿真模型中的单层Hashin损伤参数，并包含原位效应的影响。

图1  原位效应导致中间层横向属性发生变化

本节主要描述包含原位效应的Hashin损伤参数的确定方法，Hashin损伤准则的10个参数的标定方法如表1所列。这里重点描述需要考虑原位效应的参数。

表1  Hashin损伤参数的标定方法

首先标定包含原位效应的横向（基体）拉伸强度（Y^T）。图2为[0°，90°，0°]多层试样的拉伸测试数据^[4]，其加载方向为0°以确保中间层所受载荷为横向拉伸。可以看到当拉伸变形量达到一定特定值时，材料中出现损伤导致刚度下降，拐点ε₀为损伤发生时的应变量。在有限元软件Abaqus中建立和测试一致的多层试样模型，并复现此加载状态。有限元模型可以使用一个单元的简单模型。模型和测试试样具有相同的尺寸和铺层顺序。此时在材料属性中可以不包含Hashin损伤模型。边界条件采用和测试相同的设置，施加损伤发生时的应变ε₀（根据单元尺寸和应变换算成位移边界条件后施加，如图3所示）。在求解模型后，我们可以获得90°铺层中的Y方向的应力分量S22，此应力即损伤发生时中间层的横向应力，可以认为此应力分量就是中间层的Y^T。由于试样是多层试样，该应力强度已经包含了原位效应的影响。

图2  多层试样的测试曲线

图3  仿真模型获得Y^T

断裂能描述材料出现损伤后，刚度下降的规律，可以认为是应力应变曲线中刚度下降段围成的面积，如图4。Hashin损伤演化准则中的横向（基体）拉伸断裂能（G_mt）也存在原位效应的影响。这里仍然使用图2的试验数据进行G_mt的标定，但此时无法通过有限元计算直接获得G_mt值，因此需要结合Isight的Data Matching功能，对仿真和测试的加载曲线进行拟合确定G_mt。

现有测试曲线是E/E₀-ε曲线，E/E₀是相对于未损伤模型的归一化刚度，ε为应变，而仿真获得的是力-位移曲线（F-δ曲线），因此需要将测试数据进行转换，获得每个测点的力和位移，以便和仿真结果进行对比。标定仍然使用图3所示的模型进行分析，但是与之前模型有区别，需要做如下设置。

首先需要在材料属性中添加Hashin损伤模型，将Hashin损伤起始准则中的横向拉伸强度设置为上一节确定的Y^T，其他方向的参数设置一个大值（例如1e15），这样确保损伤仅以横向（基体）拉伸的模式发生，如图5。

图4  断裂能是应力应变曲线刚度下降段围成的面积

图5  设置Hashin损伤起始准则

在Hashin损伤演化准则中，将除G_mt之外的其他值都设置成大值，保证损伤仅以横向拉伸模式演化（图6），而G_mt作为变量后续在Isight模型中自动调整。为了方便测试和仿真结果的对比，在Abaqus中定义与测试数据应变量对应的载荷步历程，并在这些时刻输出支反力和位移，从而获得仿真的F-δ曲线。

利用Isight集成仿真模型并进行参数寻优，设计变量为G_mt值，可以利用测试数据的F-δ曲线按图4中围成的面积来估算的G_mt范围，作为DOE和优化的上下限。利用Isight的Data Matching组件进行曲线拟合（图7），最终得到最优的G_mt，此时测试和仿真的F-δ曲线基本重叠。（图8）

图6  设置Hashin损伤演化准则

图7  Isight Data Matching用于曲线拟合

图8  初始和最优G_mt下仿真和测试曲线对比

S^L描述平面内纵向剪切（剪应力方向为横向）的损伤起始强度，受原位效应影响。此时可以利用具有40-50°铺层试样的测试数据进行标定。图9是具有[0/±40₄/0_1/2]_s铺层试样的测试数据，其中40°铺层在拉伸状态下其应力状态近似于纯剪切，这样可以用标定Y^T类似的方法标定S^L。仍然建立一个单元的模型，具有和测试试样相同尺寸和铺层，在材料属性中不设置Hashin损伤定义。施加变形量为损伤起始时的应变ε₀。计算完成后提取40°铺层的剪应力S12作为S^L，如图10。

图9  [0/±40₄/0_1/2]_s铺层试样测试数据

图10  有限元模型和结果提取

Hashin损伤模型中的X^T，X^C和Y^C不包含原位效应，因此可以利用单层试样测试数据直接获得其强度值。G_ft, G_fc和G_mc也不包含明显原位效应，可以首先获得单层试样测试的加载曲线，一般如图4所示形式，G值可以通过计算曲线围成的面积获得。S^T可以取S^T=Y^C/2。

前面描述了利用测试数据和Abaqus + Isight仿真进行Hashin损伤参数拟合的方法。将此方法应用于Hexcel 8552 AS4 Uni Directional Tape材料^[5]，获得Hashin损伤参数。在复合材料加筋板模型上进行了包含目视勉强可见冲击损伤（BVID）损伤的壁板屈曲分析，并和测试数据进行了对比，结果具有一致性^[6]。

本文尝试一种结合测试数据和仿真进行复合材料Hashin损伤参数的拟合方法。这种方法需要一定量的复合材料单层和多层试样拉伸、压缩测试数据作为支撑，其优势在于获得的Hashin损伤参数能够考虑原位效应的影响，而这种影响通常不容易从测试直接获得。因此，方法得到的Hashin损伤参数更加合理。同时，仿真流程本身并不复杂，比较容易实现对大量测试数据的批处理。当形成了仿真参数数据积累之后，能够提升前期设计中，仿真预测复合材料损伤以及包含损伤结构强度的可信度。

[1] 达索系统SIMULIA Abaqus用户手册

[2] A. Parvizi, K. Garrett, J. Bailey, Constrained cracking in glass fibre-reinforced epoxy cross-ply laminates, Journal of Material Science 13 (1978) 195–201.

[3] Pinho, Silvestre Taveira. Modelling failure of laminated composites using physically-based failure models. Diss. Imperial College London (University of London) (2005).

[4] Barbero, E. J., F. A. Cosso, R. Roman, and T. L. Weadon. "Determination of material parameters for Abaqus progressive damage analysis of E-glass epoxy laminates." Composites Part B: Engineering 46 (2013): 211-220.

[5] Hexcel 8552 AS4 Unidirectional Prepreg at 190 gsm &35% RC Qualification Material Property Data Report, National Institute for Aviation Research

[6] 曹鹏，张伟等，含BVID复合材料的数值模拟和优化方法，SAMPE会议论文，2020

-END-