复合材料层压板低速冲击刚度退化方案仿真研究
伊鹏跃，于哲峰，汪海
（上海交通大学 航空航天学院，上海 200240 ）
摘要：针对复合材料层压板低速冲击仿真，根据应力更新、材料弹性参数折减和基于应变渐进失效的不同刚度退化思路，改进传统损伤刚度折减方法，通过ABAQUS分别编写了三种刚度退化方案的VUMAT子程序，引入粘滞规律克服刚度退化的数值计算收敛困难，结合实验进行有限元仿真，研究比较了不同刚度退化方案下冲击响应的异同，结果表明：改进的三种刚度退化方案都可较准确地描述低速冲击下复合材料失效过程；应力更新方案，思路简单清晰，但失效过程应力变化剧烈，增量步数多，计算效率低；弹性参数折减方案中，根据失效模式调整折减系数，结合粘滞规律，响应平稳；前两种方案对冲击损伤形式只能定性，无法定量表征；而渐进失效方案引入合理的损伤变量，不但冲击响应连续而且能较好地表征材料损伤形式与程度。
关键词：复合材料 低速冲击 刚度退化 VUMAT 有限元
Stiffness degradation simulation methods for composite laminate
subjected to low-velocity impact
YI Pengyue, YU Zhefeng, WANG Hai
(School of Aeronautics and Astronautics, Shanghai Jiao Tong University, Shanghai 200240,China) Abstract: To study simulation of composite laminate under low-velocity impact, stiffness reduction strategies based on stress updating , degradation of elastic constant and progressive damage have been modified and implemented in the user material subroutine (VUAMT) in ABAQUS. The difficulty of convergence in numerical calculation is overcome by viscous regularization. Differences between the three methods have been studied according to the finite element analysis and experimental data. Conclusions drawn from the comparison are that, the analytical results are ideal compared to test results. The three stiffness degradation process can describe the process of material damage exactly. Solution based on stress updating is simple and clear，but the number of increment and iteration is large due to stress varying greatly during material damage and computer calculation is inefficient. Elastic constant degradation is able to adjust the degraded coefficient based on damage mode. The impact response combining viscous regularization is ideal. Material damage mode is predicted but damage degree is unknown in this two ways. Progressive damage degradation is based on material progressive damage model and the response is continuous. Material damage mode and its degree can be well simulated. Key words: Composite; low-velocity impact; Stiffness degradation; VUMAT; Finite element analysis.
引言
复合材料由于具有高的比刚度和比强度、疲劳性能好等优点，在航空航天领域和其他现代工业中得到越来越广泛地应用。复合材料结构在制造和使用过程中，常常会受到外物低速冲击载荷的作用，在其内部产生损伤，如基体开裂、分层、纤维断裂等，严重降低了结构的强度和使用寿命，对结构安全构成潜在威胁。因此复合材料层合板低速冲击损伤性能受到广泛关注。
由于复合材料损伤失效机理的复杂性，目前暂无精确解析解，国内外许多学者通常借助有限元方法基于不同的失效准则及刚度退化方案对复合材料的损伤失效进
行研究。J.P.Hou[1]
等通过对失效后相应的应力分量置零来实现刚度的等效折减；滕锦[2]
等也利用类似的应力更新方案来研究了
z-pin增韧复合材料的冲击损伤；C.F.Li[3]
等在上述思路的基础上提出了对失效后材料弹性参数进行折减的方案并成功运用到
有限元的分析预测中；贾建东[4]
等基于Hashin失效准则对失效后材料弹性参数乘以相应的折减系数来实现材料性能的衰减并成功预测了复合材料层合板冲击后剩余
压缩强度；张彦[5]
等推导了基于Hashin失效准则，建立了基于应变损伤并引入损伤变量的复合材料层合板逐渐累积损伤模型，并成功预测了层合板低速冲击损伤；Yuequan
Wang[6]
等根据Hashin失效准则及损伤演化规律，构建损伤变量实现对损伤材料刚度阵的折减，并成功预测了带孔层合板的拉伸与
压缩极限强度；王仁鹏[7]
等基于Hashin失效准则通过引入损伤变量分别对材料的弹性参数进行折减，并成功分析了层合板在静压痕作用下的损伤阻抗。
B.G.Faozon[8]
等发展了复合材料损伤后的非线性弹塑性本构关系，并考虑了材料失效后裂纹闭合对材料刚度阵的复原作用；
L.M.Daniel[9]
等提出了基于应变率的复合材料失效准则，给出了基于应变率的材料参数的变化规律。虽然关于复合材料损伤演化的研究在不断进行中，然而总结目前有限元分析中广泛应用的复合材料层合板失效后刚度退化方案，基本思路一般有以下三类：
(1)通过将失效后相应应力分量降为零来更新应力从而实现刚度的等效退化。(2)材料失效后对相应弹性参数进行一定程度的折减来实现刚度的退化。(3)利用基于应变的失效判据并引入损伤变量将应变与材料刚度折减程度进行关联来实现材料的损伤演化。本文基于以上三种基本思路分别改进并提出了三种相应的刚度折减方案，通过ABAQUS的VUMAT子程序二次开发接口，实现了复合材料层合板三种刚度折减方案的材料本构模型，将三种具体方案分别应用到复合材料低速冲击有限元分析中，来研究不同方案下层合板低速冲击响应和损伤情况，并参照文献实验数据对各结果进行分析比较，总结不同折减方案对仿真结果的影响。
1.3D Hashin失效准则
由于Hashin[10]
失效准则已被成功地发展并较好地应用于复合材料层合板的损伤破坏分析中，所以本文在判断复合材料失效时选用三维Hashin失效准则，表达式如下：纤维断裂（?11?0）： 222e???11???????f?X???12???13? (1)
?t??S12??S13? 纤维屈曲（?11?0）： 2e????11?fX? (2)?c? 基体拉伸失效（?22??33?0）： 2e???22??33??1?Y????2m?2???23??22?33??t??S23??22(3)
??12???13?
?S????12??S13?基体压缩失效（?22??33?0）：
2e1??Yc?????222?33?m?Y????1???22??33?????c???2S12????2S12?(4)
1???2???2
?2??12
S2?2322?33????S???13?2312??S13?拉伸分层失效（?33?0）： ??222e?33??Z?????13??????13?d??(5)t??S13??S13??0 剪切分层失效（?33）：
22??????(6)ed??13???13?SS?13??13? 式(1)-(6)中：1为纤维纵向方向；2为横向方向；3为厚度方向，?11、?22、?33、?12、?13、
?23分别为1、2、3方向正应力和1-2、1-3、在子程序中实施时，为每一种失效模式定义一个状态变量，一旦发生该种失效，则相应的状态变量置为1。每个增量步开始时，首先判断该单元各状态变量是否为1，若为1，则对相应的应力分量置零，然后利用主2-3方向剪应力。Xt、Xc分别为纤维纵向拉伸和压缩强度；Yt、Yc分别为横向拉伸和压缩强度；Zt、Zc分别为厚度方向拉伸和压缩强度 ，S12、S13、S23分别为1-2,1-3,2-3方向剪切强度。当ei ?1，i=f，m，d时，材料发生失效。
2.刚度退化方案
2.1应力更新刚度退化方案
当复合材料层合板中单元应力水平满足基于应力描述的失效准则后，材料将发生破坏，其承载能力与破坏前相比会发生显著的变化。复合材料层合板失效前其性能一般表现为线弹性，材料破坏后，单元的相应应力可视为零，即认为破坏的单元对结构的承载不再做贡献。根据此原则，通过相应的应力分量置零实现对破坏后的单元进行等效刚度退化，具体应力更新方案（下称刚度退化方案1）如表1所示。
表1 应力更新方案
Table.1 Stress updating strategy
失效形式 应力更新
纤维断裂 ?11?0?22?0?33?0?12?0?13?0?23?0 纤维屈曲 ?11?0?22?0?33?0?12?0?13?0?23?0
基体拉伸 ?22?0?12?0 基体压缩 ?22?0
拉伸分层 ?33?0?13?0?23?0 剪切分层
?33?0?13?0?23?0
根据上述思路，当单元发生失效后相应方向应力置为零，等效于该单元相应的刚度阵分量为零，如此单元在该方向不再对冲击头产生反力作用，在有限元中单元易过度变形扭曲，接触反力-时间历程曲线有大振幅振动，难以观察其变化趋势，且接触力幅值整体偏低。为了保证有限元计算，须对其进行改进，程序中通过应力更新对失效单元保留一定的残余应力来防止单元刚度的过度软化，如此取得了较理想的仿真结果，证明了其可行性。
程序所传入的应变增量，与刚度阵求积得到应力增量，通过与上一增量步该单元应力相加实现应力更新，将该步应力值代入失效准则，再判断是否有新的失效模式发生，若有，则对相应的状态变量置1同时将该模式对应的应力分量置0；对于已经产生某种失效的单元，在增量步开始时已对该失效模式对应的应力置零，所以该步其应力更新后保留的残余应力大小即为该增量步中的应力增量。至此该单元的此增量步结束，然后计算下一个单元的此增量步，直至该增量步下的所有单元的应力更新完毕后，开始进入下一增量步，如此循环。
2.2参数折减刚度退化方案 在计算过程中，材料积分点应力满足失效准则后，则对材料的弹性参数乘以某一折减系数来实现失效后材料的刚度降，该种退化方案认为虽然单元已经发生损伤破坏，但只要该层仍旧埋在层合板中，它就继续对层合板的刚度有贡献。对于拉伸破坏与压缩破坏，材料刚度的退化量是不同的，材料发生压缩破坏时，将保留更多的刚度。由于参数退化方法选择是否得当对求解结果有很大影响，基于传统思路，结合有限元分析比较，改进了不同失效模式所对应折减的弹性参数及其折减系数，提出如下具体折减方案（下称刚度退化方案2），如表2所示。
表2 弹性参数折减方案
Table.2 Elastic constant degradation 失效形式 弹性参数折减
纤维断裂
Qd?0.01Q?Q?E11,G12,G13,?12,?13? 纤维屈曲
Qd?0.03Q?Q?E11,G12,G13,?12,?13? 基体拉伸
Qd?0.1Q?Q?E22,G12,G23,?12,?23? 基体压缩
Qd?0.15Q?Q?E22,G12,G23,?12,?23? 拉伸分层
Qd?0.05Q?Q?E33,G13,G23,?13,?23? 剪切分层
Qd?0.07Q?Q?E33,G13,G23,?13,?23? 表2中，E11、E22、E33、G12、G13、G23、v12、v13、v23、分别为1、2、3方向的弹性模量和1-2、
1-3、2-3方向的剪切模量与泊松比。
VUMAT子程序中同样为各种失效模式定义对应的状态变量，失效发生时置
1，然后根据表2为相应的折减系数赋值。增量步开始时，判断单元各失效模式对应的状态变量值，根据弹性参数折减系数还原上一增量步时该单元的最终材料参数值，然后施加应变增量，更新应力，将该步应力值带入失效准则，判断是否有新的失效发生，若有，则对相应的状态变量置1，并为相应的折减系数赋值，利用新的材料参数重新计算该增量步下的应力值，至此该单元的此增量步结束。
同一单元发生两种或两种以上的失效模式时，若需要对同一参数同时进行折减，由于对应各种失效模式的折减系数各不相同，基于材料损伤的不可逆性，此时取折减系数的最小值。
2.3渐进损伤刚度退化方案
当复合材料失效准则满足后，其应力和刚度并不是急剧下降，而是一个逐渐退化的过程，当结构内部发生损伤后，损伤区域的应力分布变化很剧烈，不适合再使用上述基于应力描述的失效准则。而应变在复合材料结构发生损伤前后变化比较平缓，因此更适合被用作复合材料结构中损伤演化的判据。
使用应力、应变关系式(7)，可以将基于应力描述的失效表达式转变为基于应变描述的表达式。
XEf,cf,cf,tf,cc?11?11,Xc?E11?11,Yt?E22?22,Yc?E22?22?Z?Ef,tf,cff?t33?33,Zc?E33?33,S12?G12?12,S13?G13?13?Sf? (7)
23?G23?23,?11?E11?11,?22?E22?22,?33?E33?33??12?G12?12,?13?G13?13,?23?G23?23??根据公式(1)-(7)推导得基于应变描述的3D
Hashin失效准则：
纤维断裂（?11?0）： e????22211???12???12?f??f,t???f??(8)??f?11??12???12? 纤维屈曲（? 11?0）： e??211?f????f,c?(9)11? 基体拉伸失效（?22??33?0）：
e????22E2233?33???23?m??f,t?f,t???f??22E22?22???23?E22(10)
22E33?22?33????12???13?f?
?G?2???f?2323???12??13?基体压缩失效（?22??33?0）：
?eE?f,c221?m??22?2?f,??E22?22?E33?33????2G12?f12?Ec22?22???22?E22?22?E33?33?23?2Gf??????E22?22E33?3312?12???f??(11)
23?G?f?2323?22??12???f?????
13?
12???f?13?拉伸分层失效（?33?0）：
222e??33?????f,t???13???23?d??f???f?(12)33???13???23? ?0 ）：
剪切分层失效（?3322e???13???23?d???f???f?(13) 13???23?式中：e、ee
fm、d分别为纤维失效因子、基体失效因子、分层失效因子；?11、?22、?33、
?12、?13、?23分别为1、2、3方向正应变与1-2、1-3、2-3方向工程剪应变（VUMAT中
剪应变为张量剪应变），?f,tf,c11、?11分别为纤维
纵向拉伸与压缩强度对应的极限应变，?f,t22、?f,c22分别为横向拉伸与压缩强度对应的极限应变，?f,tf,c33、?33分别为沿厚度方向拉伸与压
缩强度对应的极限应变，?ff、?f12、?1323分别为1-2、1-3、2-3方向剪切强度对应的极限剪
切应变。 基于上述转换，现引入与失效因子相关的损伤变量df、dm、dd来表征纤维、基体和分层损伤的程度，它们的定义域为[0,1],当
di=0时，表示未有损伤产生，当di=1时，
表示材料完全破坏失效，当0 材料内部有损伤存在，该部分材料的刚度需 要进行折减退化。 为了表征由应变描述的Hashin准则中失效因子与损伤变量的关系，本文引入了由 Linde[11]等提出的复合材料的损伤演化模型, 通过下面的关系式将失效因子与损伤变量关联起来。 ?d??E11??11f,i?2?ff?1?Lc/G? f??f?1?e?/ff，ff?ef(14) ?f,i??E22?22?(15) ??f,i2c??E33??33??fd?1?L/Gd??(16) dd?1?e?/fd，f?edm?1?e????f2m?1?Lc/Gm??/fm，fm?em对于刚度退化材料，其本构模型可能会 由于其刚度瞬间的大幅度折减易引起数值计算收敛困难，在显式分析步中往往表现为dd式中：Lc 为单元的特征长度，G f、Gm、Gd为材料沿纤维纵向、横向及厚度方向的断裂韧性，i的值根据材料应力状态受拉或受压分别赋值为t或c。 该种损伤演化方法基于断裂韧性并引入了单元特征长度，如此能够尽量减小应变软化的本构模型中网格密度对结果精确度的影响。 在损伤变量的基础上，需要对损伤材料的刚度阵进行折减，由于其刚度退化与应变相关，大变形时为防止材料过度软化，分析比较有限元结果，对传统方法改进，通过引入与纤维、基体与分层对应的折减系数d1、d2、d3，为材料保留一定的刚度防止其过度软化，具体方案（下称刚度退化方案3）如下，各分量表达式如式(17)-(25)，其余分量为0。 Cd11??1?d1df??1?d1df?C11 (17) Cd22??1?d2dm??1?d2dm?C22(18)? Cd33?1?d3dd??1?d3dd?C33 (19) Cd12??1?d1df??1?d2dm?C12 (20) Cd13??1?d1df??1?d3dd?C13 (21) Cd23??1?d2dm??1?d3dd?C23 (22) Cd44??1?d1df??1?d2dm?C44 (23) Cd55??1?d1df??1?d2dm?C55 (24) Cd66??1?d1df??1?d3dd?C66 (25) CddddCdd21?C12C31?C1323?C23式中： d (26) 1、d2、d3定义域为[0,1],可通过合理赋值防止材料刚度的过度软化，此处取d1=d2=d3=0.8。 增量步开始时，首先更新应变值，根据应变描述的失效准则，判断是否有失效发生，若有，计算损伤变量值与折减后的刚度阵，最后通过折减后的刚度阵与应变值求积来更新该增量步的应力，至此该单元的此增量步结束。当没有失效发生时，损伤变量df、dm、dd一直为0。由于材料损伤的不可逆性与不可愈合性，损伤变量值是历史相关的，同一单元在每一增量步下都计算一次损伤变量值，损伤变量的值取至此增量步为止的所有值中的最大值。 3.粘滞规律 单元的扭曲而导致分析进程的终止。为克服 此问题，本文引入由Duvaut等[12] 提出的粘滞规律技术，对于方案2中参数的折减系数和方案3中的损伤变量不直接使用，而是用与其对应的粘滞变量替代并应用到刚度退化的过程中，粘滞变量的定义如下： d?1v???d?dv?(27) 其中：d?v为粘滞变量，dv为dv对时间的导数，d为折减系数或损伤变量，?为粘滞系数。粘滞系数的选择取决于对结果精度的要求，当粘滞系数选取较大时，其会明显地延缓材料的刚度折减进程，参照时间增量步长，在收敛的前提下，应保证该系数尽量小，为方便各方案间比较分析统一取?=2.5X 10-4。 在数值计算中，t0??t时刻的粘滞变量表示为： d?t?v?t0??t?????td?t0??t?????tdv?t0?(28) 4.仿真模型 本文参照文献[13] 中的试验数据，材料参数如表3所示： 表3 材料参数 Table.3 Material property 弹性参数 强度参数 E11 127GPa Xt 1400 MPa E22=E33 10 GPa Xc 930 MPa G12=G13 5.4 GPa Yt 47 MPa G23 3.05 GPa Yc 130 MPa V12=V13 0.34 GPa Zt 62.3 MPa V23 0.36 Zc 130MPa S12=S13 53 MPa 密度 1540kg/m3 S23 89 MPa 复合材料层合板直径为80mm，单层厚度0.18mm，铺层数为10，共有[0]10、[0/90/0/90/0]s和[+45/-45/+45/0/90]s三种铺层顺序不同的试验件（下称试件1、2、3），周边固支，冲头直径为16mm,质量为1.205kg，冲击速度为3.13m/s，冲击能量为5.91J。 在ABAQUS中建立有限元模型，对于分层损伤的模拟，本文没有使用Cohesive单
低速冲击下复合材料刚度退化方案仿真研究12-15
