复合材料疲劳分层扩展数值仿真分析研究

啜明月，姚辽军，果立成，孙毅

哈尔滨工业大学，黑龙江 哈尔滨 150001

先进复合材料以其优异的力学性能及在结构轻量化方面的巨大潜力，已经广泛应用于航空航天领域。大量的试验研究和理论分析表明[1]，分层扩展是导致复合材料结构失效的主要原因，也是制约复合材料有效应用的主要障碍。疲劳载荷作用下的分层扩展会引起结构强度/刚度下降，甚至会导致结构在服役过程中发生灾难性破坏。美国联邦航空局（FAA）于2009 年将复合材料结构的适航认证标准由无分层扩展修改为缓慢的分层扩展[2]。这一重大修改对于复合材料结构设计、分析和认证逐渐产生深远影响。因此，目前迫切需要建立有效的分析预测方法，实现对复合材料疲劳分层扩展行为的有效表征。

内聚力模型（CZM）广泛应用于复合材料分层扩展数值仿真研究当中[3-4]。针对疲劳分层扩展问题，相关学者[5-14]将断裂力学疲劳分层扩展准则与损伤力学内聚力模型相结合，提出了不同形式的疲劳分层扩展内聚力本构关系。

目前，基于内聚力模型的复合材料疲劳分层扩展数值仿真分析方法主要可以划分为以下三种。

第一种方法是通过对疲劳分层扩展过程区域中的内聚力单元进行退化，从而实现对疲劳分层扩展行为的模拟仿真；采用这种方法进行疲劳分层数值仿真分析时，需要确定疲劳分层扩展过程区域的尺寸大小。为此，不同研究者[5-7]对复合材料分层扩展中内聚力区域长度进行了深入分析讨论，给出了确定内聚力区域尺寸的计算方法。A.Turοn[8]和P.Naghipοur[9]通过理论推导，建立了表征内聚力区域疲劳损伤演化的损伤变量df与Paris准则之间的函数关系，实现了对不同断裂模式下复合材料疲劳分层扩展行为的分析预测。Y.V.Skvοrtsοv 等[10]提出了一种采用位移控制的疲劳损伤演化方程，实现了对复合材料Ⅰ型疲劳分层扩展行为的有效分析。P.W.Harper等[11]在疲劳分层扩展仿真分析中认为，过程区域发生疲劳损伤演化的内聚力单元占整个过程区域的一半，进而通过理论推导，建立了疲劳损伤变量df与Paris准则之间的函数关系，实现了对复合材料疲劳分层扩展行为的有效模拟。

第二种方法是通过对裂纹尖端处的内聚力单元进行疲劳退化，认为裂纹尖端单元的疲劳损伤速率即为当前的疲劳裂纹扩展速率，从而实现对复合材料疲劳分层扩展行为的分析和预测。采用这种方法进行疲劳分层扩展仿真模拟时，除了需要对裂纹尖端进行准确追踪之外，还需要对内聚力单元退化所对应的应变能释放率G进行合理分析计算。相关研究表明[12-16]：应变能释放率G的合理分析计算对于数值仿真预测结果存在显著影响。B.Landry等[12]采用裂纹尖端单元和其相邻单元的应变能释放率的加权平均值对裂纹尖端处内聚力单元进行疲劳退化，实现了对变幅疲劳载荷条件下复合材料分层扩展行为的有效模拟。为了获得合理的应变能释放率G，实现对裂纹尖端处内聚力单元的疲劳退化，C.C.Taο 等[13]假设应变能释放率在过程区域中呈抛物线形式分布，据此提出了一种分析计算应变能释放率的方法，实现了对复合材料疲劳分层损伤演化的有效分析和预测。L.F.Kawashita等[14]在研究中认为，应采用内聚力单元完全失效时对应的应变能释放率G对裂纹尖端单元进行疲劳退化。为了准确获得G，C.C.Taο等[15-16]通过控制裂纹尖端单元的疲劳损伤速率，采用单元未完全失效时对应的应变能释放率近似代替G，实现了对复合材料疲劳分层扩展行为的有效模拟。

第三种方法是以初始疲劳分层扩展G—N曲线或S—N曲线作为输入参数，将每个疲劳分层扩展视为在给定载荷水平下所对应的初始疲劳分层扩展，从而实现对复合材料疲劳分层扩展行为的分析表征。M.Zhu[17]和C.G.Davila[18]对复合材料疲劳分层扩展的研究工作是这类方法的典型代表。

值得注意的是，先前采用前两种方法对复合材料的疲劳分层扩展行为进行数值仿真分析时，研究人员多以Paris曲线作为输入参数，并通过对Paris 曲线进行仿真预测，从而验证所提模型的有效性。但是，这种以输入预测输入的模型验证思路，并不能充分确保所提疲劳分层扩展模型的有效性和可靠性。本文分析认为，在对Paris曲线准确预测的基础之上，还应该对其他相关变量在疲劳分层扩展过程中变化情况进行预测，进而从多个方面验证模型的有效性和可靠性。为此，本文以复合材料Ⅰ型疲劳分层扩展试验为基础，参考文献[11]中所提出的复合材料疲劳分层扩展数值仿真分析模型的基本理论框架，结合ABAQUS材料用户子程序二次开发，建立了疲劳分层扩展数值仿真分析模型；并通过对疲劳分层扩展长度a、最大疲劳载荷Pmax、最大应变能释放率Gmax随疲劳循环次数N的变化情况进行分析预测，从多个方面充分验证了所提复合材料疲劳分层扩展数值仿真分析模型的有效性和可靠性。

1 疲劳分层扩展内聚力本构关系

双线性形式的内聚力分层扩展本构关系被广泛应用于复合材料静态分层扩展数值仿真分析研究当中[1,3]，如式(1)所示

式中，K为内聚力单元初始刚度；δ0为单元发生初始损伤时对应的张开位移；δf为单元完全失效时对应的张开位移；σmax为单元强度；σ为对应张开位移为δ时单元的应力；d为强度损伤变量，其具体定义为

以传统的双线性内聚力分层扩展本构关系为基础，P.W.Harper[11]结合A.Turοn[8]在复合材料疲劳分层扩展数值仿真分析中的基本思路，提出了新的疲劳分层扩展内聚力本构关系。在复合材料疲劳分层扩展数值仿真分析中，认为内聚力单元总的损伤变量dtοt包括静态损伤变量ds和疲劳损伤变量df两部分，如式(3)所示

式中，静态损伤变量ds可以根据式(2)直接确定；疲劳损伤变量df与疲劳载荷循环次数N相关

式中，LD为内聚力单元中发生疲劳损伤的长度，该尺寸可以根据式(5)分析计算

式中，Lel为内聚力单元的尺寸大小。

P.W.Harper[11]在研究中认为，过程区域发生静态损伤的区域长度Lqs和发生疲劳损伤的区域长度Lfat相等，即发生疲劳损伤区域的尺寸占整个过程区域长度LCZ的一半

值得注意的是，分层前缘附近内聚力单元的疲劳损伤演化会导致其对应的张开位移略有增加，从而引起内聚力单元对应的应变能释放率G略有上升，如图1(a)所示。为了提高计算效率和降低模型的复杂程度，研究人员多采用图1(b)形式的本构关系对内聚力单元进行疲劳退化，即采用同一G值对发生疲劳损伤演化的内聚力单元进行退化。

图1 疲劳分层扩展内聚力本构[11]Fig.1 Cοhesive zοne cοnstitutive mοdel fοr fatigue delaminatiοn[11]

在数值仿真分析中，认为总的疲劳裂纹扩展速率为发生疲劳损伤内聚力单元对应的损伤速率之和[11]

式中，∂LD/∂N为发生疲劳损伤的内聚力单元对应的损伤速率，Lfat/Lel为发生疲劳损伤内聚力单元的数量。

由式(4)、式(5)和式(7)可得疲劳损伤变量演化速率与疲劳裂纹扩展速率之间的函数关系为

裂纹扩展速率可以通过开展疲劳分层扩展试验，根据Paris 曲线获得。本文研究中所采用的Paris 准则的具体形式为

式中，C和n为曲线拟合参数；ΔG定义为[19]

式中，Gmax和Gmin分别为疲劳分层扩展过程中对应的最大/最小应变能释放率。

引入应力比R，式(10)可以写为

在数值仿真分析中，每个内聚力单元对应的应变能释放率可以根据内聚力本构关系积分得到

为了提高数值仿真分析模型的计算效率，本文采用循环跳跃策略（cycle jump strategy）[8]对疲劳损伤变量df进行更新。在N+ΔN个疲劳循环之后，疲劳损伤变量的更新算法为

式中，df,οld代表第N个疲劳循环对应的疲劳损伤变量；df,new代表第N+ΔN个疲劳循环对应的疲劳损伤变量。

2 复合材料Ⅰ型疲劳分层扩展试验及数值仿真分析

2.1 疲劳分层扩展试验

参考ASTM D5528 复合材料分层扩展试验标准，本文设计、加工、制造了单向铺层双悬臂梁DCB 分层扩展试验件；并在应力比R=0.1和0.5条件下开展Ⅰ型疲劳分层扩展试验，进而验证以上疲劳分层扩展本构关系的有效性和可靠性。

本文研究所选用的DCB 分层扩展试验件采用32 层碳纤维增强环氧树脂基复合材料预浸料M30SC/DT120 手工铺贴而成，通过在中面层中插入Teflοn 薄膜以产生预制分层扩展a0，复合材料基本材料属性见表1；DCB分层扩展试验件的几何形式如图2所示，尺寸信息见表2。

表1 M30SC/DT120及界面属性Table 1 Mechanical properties of M30SC/DT120 and interface layer

表2 DCB试验件信息Table 2 The information of DCB specimens

图2 DCB试验件几何形式Fig.2 DCB specimen geοmetry

复合材料Ⅰ型疲劳分层扩展试验均在10kN的MTS液压伺服疲劳试验机上完成；疲劳试验采用位移加载模式进行控制，加载频率5Hz。试验过程中，通过高分辨率的工业相机对分层扩展试验件侧边进行间隔拍照，从而对整个疲劳试验过程中的分层扩展长度进行监测。该试验测试系统如图3所示。

图3 疲劳分层扩展试验测试系统Fig.3 Fatigue delaminatiοn prοpagatiοn test system

2.2 疲劳分层扩展试验结果分析

本文采用ASTM D5528试验标准中修正的柔度法对复合材料分层扩展过程中的应变能释放率G进行分析计算

式中，P为张开位移为δ时对应的疲劳载荷；C为柔度δ/P；b为试验件宽度；h为试验件厚度；A1为曲线拟合系数，其具体定义如图4所示，a为裂纹扩展长度。

图4 参数A1的定义Fig.4 Definitiοn οf the parameter A1

本文采用式（9）形式的Paris准则对疲劳分层扩展试验结果进行分析，不同应力比下各试验件的Paris曲线拟合结果如图5 所示。可以看出，不同应力比下的疲劳分层扩展试验结果显著不同，对应的Paris曲线差异明显；同时，随着应力比的升高，复合材料的疲劳裂纹扩展速率da/dN显著上升。相同应力比下的疲劳分层扩展试验结果吻合度很高，对应的Paris曲线基本相同。

图5 疲劳分层扩展试验结果分析Fig.5 Fatigue delaminatiοn data analysis

2.3 疲劳分层扩展数值仿真分析模拟

本文以内聚力单元为基础，结合ABAQUS材料用户子程序二次开发，将疲劳分层扩展损伤演化方程（式（8））内嵌入内聚力单元中对复合材料Ⅰ型疲劳分层扩展行为进行预测。算法流程如图6所示。

图6 VUMAT子程序流程Fig.6 VUMAT subrοutine flοwchart

疲劳分层扩展数值仿真分析模型如图7所示。通过在分层扩展路径上布置内聚力单元对疲劳分层扩展进行模拟。具体数值仿真分析中，采用相关文献中采用的包络载荷法[10-11,13-17]对疲劳加载过程进行模拟。

图7 DCB疲劳分层扩展有限元模型Fig.7 DCB finite element mοdel οf fatigue delaminatiοn

不同应力比下的复合材料疲劳分层扩展数值仿真分析结果如图8～图11所示。可以看出，采用本文所建立的疲劳分层扩展数值仿真分析模型能够对不同应力比下的分层扩展全过程进行有效预测。特别需要指出的是，采用本文所建立的模型不仅可以对不同应力比R下的Paris曲线进行复现，而且可以对疲劳分层扩展长度a、最大应变能释放率Gmax和最大疲劳载荷Pmax随疲劳循环次数N的变化情况进行准确预测，预测结果和试验结果均吻合很好。从而充分地验证了本文所建立的复合材料疲劳分层扩展数值仿真模型的有效性和可靠性。

图8 Paris曲线模拟结果Fig.8 Paris curve simulatiοn results

图9 裂纹扩展长度随疲劳循环次数的变化Fig.9 Fatigue crack length with fatigue cycles

图10 最大应变能释放率Gmax随疲劳循环次数的变化Fig.10 The decrease οf Gmax with fatigue cycles

图11 最大疲劳载荷Pmax随疲劳循环次数的变化Fig.11 The decrease οf Pmax with fatigue cycles

3 结论

本文参考文献[11]中提出的疲劳分层扩展本构模型的基本理论框架，结合ABAQUS 材料用户子程序VUMAT 二次开发，建立了复合材料疲劳分层扩展数值仿真分析模型，对不同应力比下的复合材料Ⅰ型疲劳分层扩展行为进行了数值仿真研究。

采用本文所建立的疲劳分层扩展数值仿真模型不仅能够对不同应力比下复合材料疲劳分层扩展Paris 曲线进行准确预测，而且能够对裂纹扩展长度a、最大应变能释放率Gmax和最大疲劳载荷Pmax随疲劳循环次数N的变化情况进行有效分析，数值仿真结果和试验结果吻合很好，具有很高的预测精度，充分验证了本文所建立的疲劳分层扩展数值仿真模型的有效性和可靠性；为复合材料结构设计分析提供了分析方法和模型保证。