增强型层合纤维增强复合材料渐进损伤模型：基于3D Hashin 失效准则的多层级分析与验证
增强型层合纤维增强复合材料渐进损伤模型：基于 3D Hashin 失效准则的多层级分析与验证
Source: Zhang Y.; Van Paepegem W.; De Corte W. An Enhanced Progressive Damage Model for Laminated Fiber-Reinforced Composites Using the 3D Hashin Failure Criterion: A Multi-Level An alysis and Validation. Materials 2024, 17, 5176. DOI: 10.3390/ma17215176（MDPI 开放获取，CC BY 4.0；PMCID: PMC11547127；原文 35 页）
Collected: 2026-10-18
Published: 2024-10-24（Received: 2024-09-24；Revised: 2024-10-17；Accepted: 2024-10-19）
Note： 本文件为原文全文中文译本，专有名词、人名与公式符号保留英文。全部插图（图 1～图 24，其中图 24 跨页拆为 a/b 两张）与 Table 7 的嵌图均提取自原 PDF，存于 attachs/hashin-figNN.png、attachs/hashin-table7-pN.png。公式由原文 MathML 还原为 LaTeX。文末「译者注：Hashin 与损伤演化的实现要点」一节为译者归纳，非原文内容。
作者：Yichen Zhang（根特大学 结构工程与建筑材料系）、Wim Van Paepegem（根特大学 材料、纺织与化学工程系）、Wouter De Corte（根特大学，通讯作者）
摘要：本文在 ABAQUS/Explicit™ 2021 框架内通过 VUMAT 子程序实现了一个基于 3D Hashin 失效准则的渐进损伤模型（progressive damage model, PDM），借助对损伤变量的精确计算以及对非单调加载条件的兼容，增强了对应力（应变）复合材料弹性段与软化段力学性能及损伤演化的刻画能力。在随后的多层级验证中发现：该模型能够准确模拟单元层面的主要失效模式，并削弱单元尺寸的影响，保证非单调加载下可靠的力学行为表达；在层合板层面，模型还能准确预测开孔单层板与层合板的弹性行为和损伤演化，并给出最终失效时的裂纹带。本文同时强调了正确选择特征长度（characteristic length）的重要性，以及如何通过选取合适的取值来削弱网格尺寸的影响。与 ABAQUS 内置 2D 模型相比，该 3D VUMAT 子程序表现出更高的精度与有效性，证明了其在刻画纤维增强聚合物（FRP）材料力学行为与损伤机制方面的价值。所增强的 3D PDM 能够准确描述复合材料在单调或非单调加载、简单或复杂应力状态下的软化过程，有效降低网格依赖性，并合理捕捉失效裂纹带，适合日后用于复合材料模型在复杂三维应力状态下的模拟与数值问题求解。
关键词：FRP；渐进损伤模型（PDM）；多层级分析；非单调加载；特征长度；数值模拟
1. 引言
在材料工程领域，纤维增强聚合物（FRP）作为现代结构应用的轻质材料备受瞩目，其在比强度、耐腐蚀性与设计灵活性方面优势突出 [1]。这类复合材料通常由聚合物基体与玻璃、碳或芳纶纤维增强体构成 [2]，已经深刻改变了从航空航天、汽车、医疗到建筑与运动器材等众多行业 [3, 4, 5, 6, 7]。FRP 的独特性能使构件既能轻量化又具备极高强度，因而在高性能应用中不可或缺 [8]。
模拟 FRP 在各种应力状态下的行为，对提升其性能并保证实际应用的可靠性至关重要 [9, 10]。在预测单向纤维增强塑料（UD-FRP）层合板损伤起始的众多理论中，Hashin [11]、Puck [12] 与 LaRC05 [13] 失效准则因精度高而获得业界广泛认可。其中 Hashin 失效准则是最成熟的、用于理解与预测单向复合材料失效行为的框架之一，其普适性与有效性已在工业应用中得到大量验证。该准则为评估复合材料纤维与基体组分的不同失效模式提供了稳健的方法，从而支持对多种载荷条件下力学失效的综合分析。与此同时，渐进损伤模型（PDM）[14–16] 在 FRP 复合构件的失效分析中扮演关键角色：通过把材料刚度退化与纤维、基体及界面的微观光损伤机制关联起来，这类模型能有效预测复杂载荷条件下损伤的起始及其扩展 [14, 17–24]。
在实际应用层面，基于 Hashin 失效准则的 2D PDM 已被纳入 Abaqus/Explicit™ 框架。此外，许多研究者发展的 PDM 并不考虑厚度方向的应力-应变分量。Maimí 等 [25] 提出了基于 LaRC04 失效准则的连续损伤模型，用于预测复合材料层合板在平面应力状态下的力学行为。Fakoor 等 [26] 对其建立的二维损伤模型开展了与材料软化相关刚度退化方法的敏感性分析。Ridha 等 [27] 与 Su 等 [28] 用各自发展的 2D PDM 预测了含缺口层合板的强度与失效模式。然而，尽管 ABAQUS 材料库自带的标准 2D Hashin 模型（本文下称 2D Hashin 模型）应用广泛，其在冲击、压入（indentation）、钻孔等三维场景中的使用却相当有限。
为解决这一问题，若干三维渐进损伤模型相继被提出。Donadon 等 [20] 与 Li 等 [29] 用这类模型预测复合层合板受低速冲击损伤的力学响应，但其分析未包含厚度方向的基体失效。Huang 等 [19] 与 Xin 等 [21] 利用发展的三维损伤模型预测 FRP 复合材料在冲击载荷下的力学响应。此外，Divse 等 [16] 关注包含钻孔损伤的层合板力学行为建模。Maimí 等 [30] 提出的模型虽考虑了厚度方向的刚度退化，却忽略了剪切模量 G23 的退化。最后，Lee 等 [31] 采用三维模型模拟了开孔层合板的拉伸试验。
综上，使用三维模型刻画 FRP 复合材料的力学行为总体上能获得相对准确的结果。但这些研究在计算等效损伤初始应变时，普遍直接用极限强度除以弹性模量 [14, 16, 25, 29, 32–34]；在复杂应力条件下这种做法并不准确。而且，目前尚缺同时考虑面外刚度退化与非单调加载、又提供完整多层级验证的三维渐进损伤模型文献。
针对上述不足，本文在 ABAQUS VUMAT 框架内实现了一个融合 3D Hashin 失效准则的渐进损伤模型。该模型显式考虑厚度方向的刚度退化以及非单调加载带来的影响；此外，通过引入一种确定等效初始损伤应变与应力的新方法，显著提高了有限元模拟的预测精度。随后，模型经历了从单元级到层合板级的多尺度充分验证。结果表明，该模型能够准确预测复合材料在包括非单调条件在内的多样复杂应力状态下的力学行为。
2. 基于 3D Hashin 准则的连续损伤框架
本节详细阐述所发展的基于 3D Hashin 失效准则的 PDM 的理论基础与实现流程。该模型的一项显著创新在于采用了确定等效损伤初始应变 与等效损伤初始应力 的新方法（见 2.3 节）。该方法能够在复杂应力状态下更准确地计算等效损伤起始参数，从而更精确地刻画 Hashin 损伤的起始与演化。此外，该 PDM 引入了厚度方向的刚度退化以及非单调加载的影响，并利用单元特征长度降低模拟结果对网格尺寸的依赖。本节将展开模型的理论基础以及子程序中关键参数的确定方法，最后给出流程图，说明在 ABAQUS/Explicit™ 中经由 VUMAT 子程序实现的材料点损伤演化全过程。
2.1. 刚度列式
存在损伤时，单向层合板（UD）的刚度矩阵会发生相应变化。本研究将 Kachanov [22] 与 Matzenmiller 等 [23] 针对平面应力条件提出的损伤刚度矩阵推广到三维应力状态 [15]。损伤张量 、有效应力（effective stress） 与名义应力 之间的关系由式 (1) 给出，其中 是由损伤参数构成的对角矩阵，见式 (2)：
式中 、 、 （取值 0～1）分别为纤维、基体与剪切模式的全局损伤变量，可按式 (3) 表达：
式中 、 、 、 （取值 0～1）分别为纤维拉伸、纤维压缩、基体拉伸与基体压缩的局部损伤变量，由损伤演化方法确定；当某变量达到 1 时，表示该方向的组分完全失效。这四个变量的计算方法见 2.2 节。另外，为控制基体拉、压失效导致的剪切刚度退化，引入系数 与 ，取值分别为 0.9 与 0.5 [31, 32, 35, 36]。
随后，损伤后层合板的柔度矩阵由式 (4) 得到，式 (5) 即材料损伤后的柔度矩阵 ：
进一步，对含损伤的柔度矩阵 求逆即得损伤刚度矩阵 （式 (6)）。求得的 用于本研究所开发的 VUMAT，并按式 (7) 更新应力：
其中
2.2. 损伤起始准则
单向复合材料的失效准则通常用以失效模式相关的应力多项式函数表示。为描述 UD 单层板在三维应力条件下的损伤起始状态，本文纳入了 3D Hashin 准则 [11]。该渐进损伤起始准则考虑四种损伤模式：纤维拉伸（fiber tension, FT）、纤维压缩（fiber compression, FC）、基体拉伸（matrix tension, MT）与基体压缩（matrix compression, MC）。通过编制用户材料子程序（VUMAT），这四个三维失效准则在 ABAQUS 中得以实现。
三维 Hashin 损伤起始准则 [11] 的数学表达为式 (11)–(14)：
• （1）纤维拉伸模式（ ）：
• （2）纤维压缩模式（ ）：
• （3）基体拉伸模式（ ）：
• （4）基体压缩模式（ ）：
、 、 、 为各模式的损伤起始参数； 为相应方向的应力分量； 、 为 UD 单层板的纵向拉伸、压缩强度， 、 为横向拉伸、压缩强度； 、 、 为三个方向的剪切强度。Hashin 准则中涉及的所有失效模式极限强度均取正值。每个损伤起始参数取值 0～1：任一参数达到 1，即表示相应模式的损伤开始演化。
2.3. 损伤演化
确定 UD-FRP 单层板损伤后的力学行为至关重要。损伤发生前材料呈线弹性，因此令全局损伤变量（ 、 、 ）为 0，即可按式 (8) 的本构关系确定损伤前刚度矩阵。损伤发生后继续加载会使层合结构刚度下降，式 (8) 的刚度矩阵随之更新，此时全局损伤变量从 0 开始增大。这些全局损伤变量由四个局部损伤变量 （i = ft, fc, mt, mc）按式 (3) 计算得到，各局部变量取值从 0 增至 1。局部损伤变量 用式 (15) 计算。
需要说明的是，全局与局部损伤变量都以 0～1 标度：0 表示无损伤，1 表示完全损伤。
式 (15) 中， 为某一损伤模式的等效应变，其计算详见式 (18)–(21)。此外，借助等效应变 与特征长度 ，可得到单元的等效位移 与等效应力 （见表 1）。
表 1 各损伤模式的等效位移与等效应力
式 (15) 中， 与 分别对应各损伤模式损伤起始时刻与损伤最终时刻的等效应变；上标表示四种损伤模式 FT、FC、MT、MC。准确确定 与 对正确描述刚度退化阶段的力学行为至关重要。需要澄清的是：当损伤起始变量 时，该损伤模式下材料点损伤起始时刻的等效应变即记为 （如基体拉伸模式记为 ）。
在迄今查阅的全部文献 [14, 16, 25, 29, 32–34] 中， 都按式 (16) 确定；当法向或剪切分量之一不显著时，由此得到的损伤初始应变是准确的。但当法向与剪切分量都不可忽略时，算得的等效初始应变 就不再准确。这种偏差出现在单元处于复杂应力状态时——此时不应只用式 (16) 而忽略剪切分量来计算 。
例如，对含 ±45° 铺层的层合板施加拉伸时，面内剪应力 对材料损伤起始参数 影响显著。若用式 (16) 计算初始损伤应变（ ），就完全忽略了剪切分量对 的贡献，显然不合理。
因此，本文采用新方法确定单元在各损伤模式下的 ：不再用固定的材料极限强度与弹性模量预先计算等效初始损伤应变，而是严格遵循 Hashin 失效准则，在损伤起始变量 首次等于 1 的时刻确定该损伤模式的等效应变，并以该值作为单元的初始等效应变。该方法把损伤起始阶段与刚度退化阶段衔接起来，保证每个单元都能依据其真实应力状态获得初始等效应变 。
类似地，等效最终应变 也不再按式 (17) 计算。而是当损伤起始变量 首次等于 1 时，先用表 1 中的公式计算等效初始损伤应力 ；再根据单元所处的损伤模式，从式 (18)–(21) 中选取相应公式计算该应力状态下的等效最终损伤应变 。式 (18)–(21) 中包含各失效模式的断裂能 与特征长度 。
本文采用的局部刚度退化遵循线性软化模型（见图 1）。该模型能有效捕捉损伤演化过程中刚度的渐进折减，从而在多种载荷条件下对材料行为给出更准确、更真实的模拟。
图 1 本文 PDM 理想线性软化模型的应力-应变响应。
• （1）纤维拉伸模式（ ）：
• （2）纤维压缩模式（ ）：
• （3）基体拉伸模式（ ）：
• （4）基体压缩模式（ ）：
式中 为 Macaulay 括号算子，定义为 。
2.4. 特征长度
应变软化会因应变局部化而产生显著的网格依赖性，网格加密时耗散能随之减小 [16]。这意味着断裂过程中耗散的能量与单元实际尺寸直接相关 [37]。为此引入特征长度 来控制损伤过程，以准确描述应变软化。要在损伤过程中尽量削弱网格依赖性，必须基于单元特征长度 对其断裂能进行归一化。
在 ABAQUS 中，实体单元默认的 为单元体积的立方根 [35]。本研究实现的软化模型中，单元厚度方向的 取实际厚度，面内方向取面积的平方根。这一配置旨在尽量减小计算量的同时，尽可能消除较薄六面体单元在面内方向特征长度上的误差。
此外，由于单元在局部化起始时刻的弹性能必须大于其断裂能，损伤起始应变必须小于最终损伤应变 [25]，以保证单元在到达局部化与进一步软化之前具有足够的弹性能驱动该过程：
因此， 的最大取值为
该处理提升了计算模型的精度与可靠性，保证模拟结果的实用与准确。特征长度对软化段力学行为的影响将在 4.2 节专门分析。
2.5. 3D Hashin 渐进损伤模型的实现
依据本节前述 3D Hashin 损伤准则并基于断裂能，该渐进损伤模型借助 VUMAT 用户子程序在 Abaqus/Explicit 框架内得以实现。子程序的运行逻辑简述如下：
首先，子程序从有限元模型中读取材料属性以及上一增量步各材料点（积分点）的「旧」解相关状态变量（SDV）；随后按给定的位移增量更新单元应变；接着检查各材料点的局部损伤变量 ：若 （无退化），则用弹性刚度矩阵 计算应力；否则基于损伤刚度矩阵 （式 (8)）更新应力。
接下来，子程序检查上一增量步的损伤起始变量 是否等于 1。
若 ，说明材料刚度退化已经发生。子程序随后判断当前等效应变 是否超过历史最大等效应变 ：若 ，则计算新的局部损伤变量 ，并把 更新为 。检查 这一步对处理非单调加载下等效应变 回落的情形至关重要。
若 ，说明该材料点尚未起始损伤：子程序按 3D Hashin 理论（式 (11)–(14)）计算各失效模式的损伤起始参数 （i = ft, fc, mt, mc）。当 首次等于 1 时：赋 ；按式 (18)–(21) 与表 1 计算等效应变 、等效应力 、等效初始应变 与等效最终损伤应变 ；并更新 。
最后，基于新算得的局部损伤变量 ，更新材料点的损伤刚度矩阵 与应力状态。整个流程逻辑见图 2。
图 2 基于 3D Hashin 准则的 PDM 实现流程图。
通过该实现，本文 VUMAT 不仅能在材料点任意应力状态下准确刻画损伤起始与退化过程，还考虑了单元厚度方向的刚度退化，并可适用于非单调加载条件。下一节将对该 VUMAT 子程序的多个方面开展全面验证。
2.6. 层间分层建模
本研究在 ABAQUS 框架内采用内聚区模型（cohesive zone model, CZM）结合通用接触算法模拟 FRP 层合板的层间分层。模型利用基于表面的内聚接触（surface-based cohesive interaction），通过双线性牵引-分离（traction-separation）关系描述界面粘接接触应力行为 [38, 39]。层间损伤起始采用二次名义应力准则（式 (24)）预测，式中 、 、 为法向与两个切向的接触应力， 、 、 为相应模式的临界接触应力。分层扩展则选用基于混合模式能量释放率的 B-K 准则 [40] 描述界面分离：
3. 数值模拟验证
本节目的是验证所建立的基于 3D Hashin 失效准则的 PDM 的准确性与普适性。验证思路由简单（单元级）到复杂（层合板级）递进。需要说明的是，ABAQUS 材料库中的 2D Hashin 准则只能用于壳单元（S8R）或连续壳单元（SC8R）。因此在单元级，通过将本文 3D VUMAT 的结果与 ABAQUS 内置 2D 损伤模型对比，证明其精度与优越性；在层合板级，进一步将模拟结果与文献 [41] 的试验数据对比，印证子程序的准确性。本研究所有实体单元均采用带沙漏控制的八节点线性减缩积分六面体单元（C3D8R）；2D 模型采用同样带减缩积分与沙漏控制的八节点四边形面内通用连续壳单元（SC8R）。
3.1. 单元级
本部分检验所开发的 3D VUMAT 子程序在单元级能否按输入设定给出预期的力学响应。研究对象为边长 1 mm 的 3D 实体单元（C3D8R），材料为 E-glass/Polyester（玻纤/聚酯），其材料参数为假定值（详见表 2）。作为对照，建立尺寸与材料参数相同的 2D 单元（SC8R）。
表 2 E-glass/Polyester 材料参数（假定值）
3.1.1. 单单元验证
纤维拉、压失效与基体拉、压失效构成 UD 层合板的四种基本失效模式。因此有必要针对每种失效模式分别开展单元验证，各工况的边界条件见图 3。除非另有说明，下文中 1 方向为纤维方向，2 方向为面内基体方向，3 方向为面外厚度方向。由于施加于单元的拉、压力要么平行要么垂直于纤维方向，单元在垂直于加载方向的横向保持平衡，故不约束垂直于加载方向的两个自由度——这一边界设置旨在避免加载中因泊松比效应约束变形，使模型保持单轴应力状态。鉴于 VUMAT 采用 3D Hashin 损伤准则，另对 3 方向的基体失效以及面外面内剪切失效（2–3、1–3 平面）进行了补充考核。各工况结果在 4.1.1 节详细讨论。
图 3 单单元试验的边界条件。
3.1.2. 非单调拉-压加载验证
数值模拟中，由于材料与接触的非线性，结构内单元的应变状态往往是非单调的。因此有必要考核单元在非单调加载下的刚度退化行为。共进行两类非单调加载模拟：一沿纤维方向（1 方向），一沿基体方向（2 方向）。加载边界条件与 3.1.1 节相同（图 3a–d）。加载序列见表 3：第 1、3 步的设置使单元发生损伤起始但不至最终失效；最后两步（第 5、6 步）使单元发生极限拉伸失效。结果在 4.1.3 节展示讨论。
表 3 非单调加载序列（给定位移）
3.1.3. 网格敏感性试验
复合材料单元的网格尺寸对裂纹带的生成与扩展 [24, 42] 以及单元自身的力学行为影响显著。为评估损伤模型能否降低单元尺寸对模拟结果的影响，本研究采用与 Falzon 和 Apruzzese [43] 及 Reinoso 等 [44] 相同的验证方法：构造边长 1 mm 的立方体模型，分别包含 1、27、125、343 个三维实体单元（C3D8R）与二维连续壳单元（SC8R）。为简化，取横向拉伸为基本加载条件：沿 2 方向对所有模型施加 0.007 mm 拉伸位移，边界条件与材料参数保持不变。相关结果在 4.1.3 节详细讨论。
3.2. 特征长度评估
合理选取单元特征长度 是获得准确模拟结果的关键。本节通过对单位立方体模型进行不同网格划分，研究 取值对模拟结果的影响，以期确定 3D 复合材料单元的最优特征长度设置。不同网格划分方式见图 4：Case 1、2、3 与 Case 4、5、6 保证加载方向所有内部单元的几何长度一致，前者含单层单元、后者含两层单元；Case 7 亦为两层单元，但其纵向与横向几何长度不同。模拟采用两种 设置方式：（a）使用 ABAQUS 实体单元的默认 ；（b）使用按加载方向各单元实际几何尺寸确定的 。随后比较两种方式的结果。同样为简化，只关注基体方向（2 方向）拉伸加载，边界条件与材料参数与前述一致。本节结果将在 4.2 节讨论。
图 4 不同网格划分方式的模型。
3.3. 层合板级
在单元级验证 VUMAT 子程序之后，本节进一步考察该 PDM 能否准确刻画层合板尺度的力学响应。采用文献 [41] 的开孔板拉伸试验结果对比验证模拟结果。试验用 IM7/8552 碳纤维复合材料的材料参数取自 [44]，见表 4。
表 4 IM7/8552 材料参数
3.3.1. 单层开孔板拉伸试验
为在层合板尺度上逐步验证 3D VUMAT 子程序的准确性，从文献 [41] 试验的多铺层层合板中选取 0°、90°、45° 三个单层开孔板。该初步检验旨在确认子程序能否预测 FRP 在单层板级别的弹性模量、极限强度与失效模式。各开孔单层板几何尺寸为 （长 × 宽 × 厚），中央圆孔直径 2 mm，与试验多层试件一致。尺寸见图 5。沿 1 方向施加拉伸位移，不约束单层板 2 方向的自由度。
图 5 开孔层合板有限元模型的几何与边界条件。
3.3.2. 开孔层合板 [90/0/±45]3S 拉伸试验
在单元级与单层板级均验证 3D 模型精度之后，进一步考察其在层合板内的表现。本节按 [41] 的材料参数与几何尺寸建立开孔板模型并进行拉伸模拟。单层厚度 0.125 mm，铺层顺序 [90/0/±45]3S。为对比，另建五个采用 ABAQUS 材料库 2D Hashin 准则（SC8R 单元）的模型，与试验及 3D 数值结果比较。各试件几何布局见图 5，尺寸见表 5。IM7/8552 的层间性能参数 [44] 见表 6。
表 5 开孔层合板试件几何尺寸
表 6 IM7/8552 层间性能参数 [44]
4. 结果与讨论
4.1. 单元验证
4.1.1. 单单元验证
单单元验证结果见表 7：四种失效模式均可被准确刻画，输入参数与数值模拟结果完全一致。此外可以观察到，除剪切损伤外，损伤变量的演化平滑，且与应力-应变关系直接对应。值得注意的是，剪切损伤模式下最终失效时剪切刚度（ 、 、 ）并不降至零，单元仍保留 10% 的残余剪切刚度——这是因为本子程序中剪切损伤系数 取 0.9（见 2.1 节）。
表 7 单单元验证的应力-应变关系与损伤变量演化（原表跨页，嵌图整体提取为下面两张图；每行左列为应力-应变关系，右列为损伤变量演化；行序自上而下：纤维拉伸损伤（1 方向）、纤维压缩损伤（1 方向）／基体拉伸损伤（2、3 方向）、基体压缩损伤（2、3 方向）、剪切损伤（1–2、1–3、2–3 平面））
由此可确认：单单元在简单应力条件下，该 PDM 能够准确反映其力学行为与极限失效。
另外值得一提：用 ABAQUS 内置 2D Hashin 材料模型也做了单单元测试。在平面应力条件（1–2 平面）下它能给出正确结果，但一旦存在面外应力分量（S33、S23、S13）就不适用了。
4.1.2. 非单调拉-压加载验证
图 6 给出采用增强 3D VUMAT 的实体单元沿纤维方向（1 方向）经历六阶段单轴非单调加载的力-位移曲线。第一阶段为单轴拉伸段（Stage 0→1），拉伸位移使刚度在本阶段结束前开始退化， 减小：至点 1， 从初始值 33,457 MPa 降到约 29,150 MPa。继续加载，应力状态由拉（Stage 0→2）转为压（Stage 2→3）。进入压缩段后 保持在 29,150 MPa 左右，说明第一阶段拉伸引起的刚度退化在压缩中被延续计入。随压缩继续， 降至约 9300 MPa（Stage 3→4）。当单元再次回到拉伸（Stage 4→5）时， 仍与 Stage 3→4 相同，约 9300 MPa——说明压缩造成的刚度退化同样被计入其后的拉伸加载。
图 6 3D PDM 在单轴纵向（1 方向）拉-压非单调载荷下的载荷-位移曲线。
而 ABAQUS 内置的 2D Hashin 损伤模型给不出同样的结果。图 7 对比了相同非单调加载下 2D 与 3D PDM 的应力-应变曲线：两条曲线明显不重合。2D 模型中，1 方向初始拉伸造成的刚度退化没有在其后的压缩段得到体现：初始拉伸结束时模量 约 29,100 MPa，但其后的压缩段模量仍保持为初始值 33,457 MPa；同样，由压转拉时，压缩造成的刚度退化也没有在拉伸段中被计入。图中虚线为单元在相应单轴单调载荷下直到失效的应力-应变曲线。
图 7 单轴纵向（1 方向）拉-压非单调载荷下 3D PDM 与 ABAQUS 内置 2D PDM 的应力-应变关系对比。
同样地，对 3D VUMAT 在横向（2 方向）非单调加载下的力学行为进行了刻画，力-位移曲线见图 8。第一拉伸段结束（Stage 0→2）时，拉伸导致的 退化使其从初始值 9883 MPa 降至约 7480 MPa；随后的压缩段（Stage 2→3）在发生压缩损伤前 保持约 7480 MPa；从压缩段（Stage 3→4）转入第二拉伸段（Stage 4→5）时， 在两段中保持一致，约 2880 MPa。作为对比，相同载荷下 3D 与 2D 模型的应力-应变曲线见图 9：与 1 方向的结果类似，2D 模型未能同时计入拉伸与压缩加载造成的刚度退化。
图 8 3D PDM 在单轴横向（2 方向）拉-压非单调载荷下的载荷-位移曲线。
图 9 单轴横向（2 方向）拉-压非单调载荷下 3D PDM 与 ABAQUS 内置 2D PDM 的应力-应变关系对比。
结论：本文发展的 3D PDM 通过在同一单元内同时考虑同方向的拉伸与压缩损伤，能有效刻画材料在非单调加载下的损伤行为，从而更准确地表达非单调载荷下材料的力学行为。其意义在于避免在非单调应力条件下高估材料的刚度与强度，得到更合理准确的结果 [25, 44–46]。反之，用 ABAQUS 材料库 2D PDM 建立的模型在非单调加载条件下倾向于高估材料的刚度与强度。
4.1.3. 网格敏感性试验
网格敏感性分析结果见图 10～图 12，加载方向为横观方向（2 方向）。由图 10 可见，对 3D 模型，不同网格密度下的力-位移曲线保持一致；而且尽管未像 [44] 那样预先折减中部区域单元的强度，裂纹带仍出现在垂直于加载方向的中部区域（图 11a，SDV3 为基体拉伸损伤变量）。这表明模型能把失效限制在一条单元宽度的裂纹带内。相反，2D 模型达不到预期效果：损伤段刚度退化速率随网格密度增大而减小（图 10），且失效单元并不集中在同一层的裂纹带内（图 11b）。
图 10 单轴横向（2 方向）载荷下 3D 与 2D 单元的载荷-位移曲线对比。
图 11 面内基体方向（2 方向）单轴拉伸下（a）3D 与（b）2D 模型的极限失效形态。
图 12 2 方向单轴拉伸下 7 个工况的载荷-位移曲线：（a）真实特征长度 ；（b）默认特征长度 。
译注：原文图 12 的 (a)/(b) 标题如此书写（ 为 default、 为 true），与 4.2 节的记号 （默认）、 （真实）一致，疑为原文笔误处按原文保留。
上述结果证实：3D 模型的力学行为几乎不受网格尺寸影响，并能在模型中部形成合理的裂纹带。这体现了所提 PDM 在保证失效机制预测准确可靠方面的优越性。
4.2. 特征长度评估
ABAQUS 中实体单元的默认特征长度（ ）按单元体积的立方根计算；真实特征长度记为 ，即单元在相应方向（本试验为 2 方向）的实际几何长度。单元在相应方向的真实特征长度与默认特征长度之比记为 ，见表 8。
表 8 各工况的 、 与
图 12a 给出七个模型采用默认特征长度 的拉伸力-位移曲线：随 增大，模型的失效位移也增大。
这一现象源于单元断裂能按式 (25) 计算：
当所选特征长度（ ）大于加载方向的特征长度（ ）时，单元获得更高的断裂能，失效位移更大；反之若 小于 ，单元失效位移更小。只有当赋予的特征长度 等于 时，单元才呈现预期的断裂能与失效位移。单元在给定加载方向的真实特征长度 由其在该方向的几何长度定义。
为验证这一结论，七个模型均把各单元的特征长度赋为 。结果见图 12b：赋予正确的 后，所有工况都得到准确的载荷-位移曲线。此外，Case 7 进一步说明：即使不同层单元在加载方向的几何长度不同，只要每个单元被赋予与其几何尺寸匹配的特征长度，模型的损伤演化过程仍能被准确表达。
4.3. 层合板级验证
4.3.1. 开孔单层板拉伸试验
本节分别讨论 0°、45°、90° 三种开孔单层板在拉伸载荷下的应力分布与损伤演化。
• 0° 开孔单层板拉伸试验
单层板的初始损伤由基体拉伸失效（SDV3）引起，失效部位位于孔边、沿板轴对称分布（见图 13）；此阶段未观察到纤维拉伸损伤。这与无孔 0° 单层板在纵向拉伸下由纤维拉伸损伤（SDV1）主导的情形形成对比。差异的原因在于：开孔单层板孔边因应力集中处于复杂应力状态。由图 13 可观察纵向正应力（S11）、横向正应力（S22）与面内剪应力（S12）的分布：它们都在孔边达到最大值，分别为 889.6 MPa、30.8 MPa 与 69.7 MPa。但按 2.2 节的损伤起始准则，这些应力水平尚不足以起始纤维损伤，却足以造成基体拉伸失效。虽然横向正应力最大值只有 30.8 MPa，但该处的高剪应力（69.7 MPa）与横向正应力（S22）、面内剪应力（S12）共同作用触发了基体拉伸损伤。这进一步佐证了 2.5 节（译注：原文如此，实为 2.3 节）所提增强型等效损伤初始应变计算方法的合理性。
图 13 0° 单层板损伤演化过程中的拉伸基体/纤维损伤变量及相应应力分量分布。
随拉伸载荷增大，基体拉伸损伤持续演化，同处也开始出现纤维拉伸损伤，只是纤维拉伸损伤的扩展滞后于基体拉伸损伤。此时板内面内应力分布与损伤起始时相似，但应力水平更高。
如图 13 所示，开孔板的最终失效由纤维拉伸与基体拉伸两种失效模式共同造成，裂纹平行于纤维方向即加载方向，这与复合材料裂纹沿纤维方向扩展的一般认识一致。因此结论是：开孔 0° 单层板受拉失效由基体与纤维失效共同控制。
另外，0° 单层板受拉的应力-应变曲线见图 16：计算得拉伸方向弹性模量 168,116.21 MPa，略低于按经典层合板理论（CLT）对无孔 0° 单层板算得的 171,420 MPa。这证实了 0° 单层板弹性段的准确性。
• 45° 开孔单层板拉伸试验
对 45° 开孔单层板，加载中最显著的应力分量是横向正应力（S22）与面内剪应力（S12），初始损伤由基体拉伸损伤（SDV3）控制。损伤起始时这两个应力分量的分布见图 14：S22 与 S12 的最大值分别约为 58.9 MPa 与 50.3 MPa。虽然任一分量单独都未达到损伤起始临界值，但依据 2.2 节，这些应力的联合效应可诱发基体拉伸失效（式 (13)）。因此基体拉伸损伤（SDV3）出现在 S22 尚低于临界值（62.3 MPa）的位置。
图 14 45° 单层板损伤演化过程中的拉伸基体损伤变量及横向正应力、面内剪应力分布。
此外，图 14 也给出了损伤演化过程中 SDV3 与 S22、S12 的分布：可见裂纹沿纤维方向扩展，应力分布与损伤起始时刻基本相似。
最后，图 14 给出了 45° 开孔单层板失效时的裂纹扩展方向：与 0° 单层板一样平行于纤维方向。还可观察到单层板两部分在 2 方向发生相对位移——这是因为加载过程中 2 方向的自由度未被约束。
另外还验证了 45° 单层板弹性段的刚度：绘制加载过程的应力-应变曲线（图 16），得 45° 单层板加载方向弹性模量约 13,328.5 MPa，略低于按 CLT 对无孔单层板算得的 13,433.6 MPa。
• 90° 开孔单层板拉伸试验
与 0°、45° 铺层相比，90° 单层板受拉时的应力分布相对简单（图 15）：整个加载过程中横向正应力（S22）显著高于其他分量。损伤起始由基体拉伸损伤（SDV3）控制。横向应力分布见图 15，最大值 62.07 MPa，位于孔上下边缘。
图 15 90° 单层板损伤演化过程中的拉伸基体损伤变量及横向正应力分布。
损伤发生后，横向正应力峰值由初始最大值（62.07 MPa）因刚度折减降至 56.35 MPa；期间裂纹继续沿纤维方向扩展。
完全断裂时，裂纹仍沿纤维方向（垂直于加载方向）扩展，说明 90° 开孔单层板的拉伸失效由基体拉伸损伤主导。
类似地，受拉应力-应变曲线（图 16）表明 90° 单层板沿板纵向的弹性模量为 9011.24 MPa，略低于按 CLT 对无孔条件算得的 9080 MPa。
图 16 0°、90°、45° 开孔单层板受拉的载荷-位移曲线。
在 0°、45°、90° 三种单向开孔单层板的拉伸模拟中可观察到：0° 案例的损伤由基体拉伸（SDV3）与纤维拉伸（SDV1）共同控制，而另外两个案例的失效始终由基体拉伸失效（SDV3）主导；且基体损伤诱发的裂纹总是沿纤维方向扩展。这些观察与 Prabhakar 和 Waas [47]、Khechai 等 [48]、Zhang 等 [34] 的试验结果，以及 Prabhakar 和 Waas [47]、Forghani 等 [49]、Zhang 等 [34]、Divse 等 [16] 的有限元研究结果一致。
因此可以得出结论：该 VUMAT 子程序能准确刻画 0°、45°、90° 开孔单层板受拉时的损伤演化与失效模式。
4.3.2. 开孔层合板 [90/0/±45]3S 拉伸试验的数值模拟结果与试验结果对比
为验证 3D PDM 在层合板级的准确性，本节先用 3D VUMAT 深入分析讨论各铺层在弹性段与软化段的应力分布状态与最终失效模式；随后将 3D 层合板模型预测的整体极限强度与失效模式和 Camanho 等 [41] 的试验结果比较；同时给出采用 ABAQUS 材料库 2D Hashin 准则获得的开孔层合板极限强度，以突出本文 3D VUMAT 的优越性。
试件为五个不同宽度（W）与孔径（D）的开孔层合板，铺层 [90/0/±45]3S，宽径比 W/D 均为 6。几何尺寸见表 5，几何构型见图 5。
首先对孔径 2 mm 的 3D 模型做网格敏感性分析：随网格尺寸按 1.5、1.0、0.75、0.5 mm 递减，极限强度 依次升高，分别为 466.6、520.2、529.0、537.0 MPa。为保证模拟精度，后续模拟网格尺寸均取 0.5 mm。
此外，为提高计算效率启用了质量缩放（mass scaling）：缩放因子取 100、500、1000 时结果非常接近，故为保证计算效率取 1000。
• 开孔层合板各铺层的损伤预测
先讨论不同单层的应力与损伤状态。为简化，本部分选取孔径 6 mm 的模型（中等尺寸）作详细分析。虽然层合板共 20 层，但相同纤维方向的铺层观察到相似的损伤模式，因此分别在 0°、±45°、90° 铺层中各选四层展开讨论。另外，为清晰展示各铺层的损伤过程，选取四个关键时刻观察各铺层的应力分布与损伤状态：0.5 倍极限载荷、0.8 倍极限载荷、极限载荷以及失效后。
• （a）0° 铺层
整个拉伸过程中，纵向正应力（S11）、横向正应力（S22）与面内剪应力（S12）都显著，损伤主要由纤维拉伸（SDV1）与基体拉伸（SDV3）控制。
损伤起始时（0.5 倍极限载荷），由于应力集中，S11、S12、S22 的最大值都出现在孔边（图 17），分别为 1665 MPa、86.26 MPa 与 28.1 MPa。虽然这些值单独都未达到损伤的极限强度，但它们的联合效应在该区域引发了纤维拉伸与基体拉伸损伤。
图 17 [90/0/±45]3S 层合板 0° 铺层在不同载荷水平下的拉伸基体/纤维损伤变量及相应应力分量分布。
随加载进行，至 0.8 倍极限载荷，纤维与基体拉伸损伤都继续演化，且基体拉伸损伤（SDV3）的演化快于纤维拉伸损伤（SDV1），局部区域达到 1，即这些部位发生基体拉伸失效（图 17）。这与 4.3.1 节 0° 单层板观察到的现象一致；此阶段应力分布变化不大，最大应力仍出现在孔边。
达到极限载荷时，损伤与应力分布仍与前一阶段相似。此时圆孔边缘已发生显著刚度退化，孔上下边缘 S11 达到 2240 MPa，明显高于其他铺层沿拉伸方向的应力水平。这说明层合板拉伸试验中大部分载荷由 0° 铺层承担。
失效后，0° 铺层的基体拉伸与纤维拉伸损伤分布见图 21：最终断裂发生在板中部并穿过圆孔，裂纹方向垂直于加载方向；最终断裂处基体与纤维拉伸损伤都接近 1，说明最终失效由两种损伤模式共同造成。
本节的观察与 4.3.1 节类似：无论单独承载还是处于层合板的不同铺层，0° 铺层的最终失效都由纤维拉伸（SDV1）与基体拉伸（SDV3）损伤共同造成。整个过程中 S11、S22、S12 都显著，且孔边附近的纵向正应力在失效前接近临界值。
• （b）±45° 铺层
±45° 铺层的应力状态与损伤模式非常相似，故一并讨论。拉伸过程中 ±45° 单层板的损伤主要由基体拉伸损伤（SDV3）控制，应力状态则由面内剪应力（S12）与横向正应力（S22）主导。
以 45° 铺层为例：损伤起始时（0.5 倍极限载荷），S22 与 S12 在圆孔上下边缘最显著（图 18），分别为 44.93 MPa 与 76.6 MPa。由于横向应力 S22 与剪应力 S12 的联合作用，该区域发生基体拉伸损伤（图 18）。类似地，−45° 铺层在此时刻也表现出相近的损伤与应力分布（图 19）。
图 18 [90/0/±45]3S 层合板 45° 铺层在不同载荷水平下的拉伸基体/纤维损伤变量及相应应力分量分布。
图 19 [90/0/±45]3S 层合板 −45° 铺层在不同载荷水平下的拉伸基体/纤维损伤变量及相应应力分量分布。
随加载继续，损伤区域扩大、损伤程度加深（图 18、图 19）。此外可观察到 S12 呈对称分布，而 S22 呈反对称分布。
极限失效时，±45° 铺层圆孔附近局部区域因应力集中发生基体拉伸失效。但由于层合板的极限承载能力主要由 0° 铺层提供，±45° 铺层的基体拉伸损伤对层合板整体刚度影响有限，不会触发层合板的最终失效。
±45° 铺层的最终失效形态见图 21：断裂区仍位于板中部、穿过孔区并垂直于加载方向；断裂带两侧的三角区域也发生基体拉伸失效。可观察到两侧三角区的裂纹不连续、呈间歇分布。该现象可用文献 [50] 提到的基体裂纹附近应力分布模式解释：在某一临界长度（ ）内，越靠近裂纹的区域应力水平越低，因此新裂纹只能出现在临界长度边界或之外。
• （c）90° 铺层
90° 铺层的损伤过程主要由基体拉伸损伤驱动，最显著的应力分量为 S22（横向正应力）与 S12（面内剪应力）。
0.5 倍极限载荷时（图 20），横向正应力 S22 的最大值出现在紧邻孔边的区域（60.8 MPa），该处发生基体拉伸损伤；最大剪应力沿板中心线反对称分布，最大值 48.5 MPa。这说明基体拉伸损伤由这些应力共同作用导致，但以 S22 分量为主。
图 20 [90/0/±45]3S 层合板 90° 铺层在不同载荷水平下的拉伸基体/纤维损伤变量及相应应力分量分布。
随损伤演化，0.8 倍极限载荷时的应力分布与损伤起始时（0.5 倍极限载荷）相似。此阶段孔周上下区域出现数条由基体拉伸损伤（SDV3）造成的小横向裂纹，并沿板宽方向（平行于纤维方向）横向扩展。
当层合板达到极限强度时，孔上下出现等腰三角形损伤区，裂纹沿三角形边缘间歇分布，与 ±45° 铺层讨论过的损伤形态类似。此外可观察到，未发生刚度退化区域的应力水平显著高于已损伤区域。
90° 铺层的极限失效同样由基体拉伸损伤（SDV3）造成。与前面讨论的 0°、±45° 铺层类似，极限失效时的最终断裂区垂直于加载方向（图 21）。
图 21 [90/0/±45]3S 层合板各铺层的极限失效：（a）0° 铺层极限基体拉伸损伤（SDV3）；（b）45° 铺层极限基体拉伸损伤（SDV3）；（c）−45° 铺层极限基体拉伸损伤（SDV3）；（d）90° 铺层极限基体拉伸损伤（SDV3）；（e）0° 铺层极限纤维拉伸损伤（SDV1）。
小结：最终失效前，[90/0/±45]3S 层合板的大多数铺层（90° 铺层除外）裂纹不会大面积扩展，而是集中在孔周。最终失效时，四种铺层方向都产生穿过孔区并垂直于加载方向的横向断裂带。此外，除 0° 铺层外，其余铺层沿断裂带两侧出现由基体拉伸损伤造成的三角形损伤区。除 0° 铺层的极限失效由基体拉伸与纤维拉伸损伤共同造成外，其余铺层的极限损伤都主要由基体拉伸损伤（SDV3）驱动。拉伸载荷主要由 0° 铺层承担，层合板的最终失效由该层断裂引起。
• 层合板整体极限载荷与失效模式
层合板的纵向应力水平 （MPa）定义为 ，其中 P（N）为拉伸载荷，W（mm）为宽度，T（mm）为层合板总厚度。按此定义可生成不同孔径试件受拉的应力-应变曲线，见图 22。
图 22 数值模拟预测的不同开孔尺寸 [90/0/±45]3S 层合板的应力-应变曲线。
基于图 22 的计算发现：孔径 6 mm 试件的数值模拟弹性模量为 62,774.5 MPa，略低于按经典层合板理论（CLT）对无孔构件算得的 64,598.8 MPa，误差 2.8%，属合理范围；其余四个试件的弹性模量亦然。
五个试件受拉失效应力（试验与数值）列于表 9。3D 模型相对试验的最大误差为 9.75%，出现在孔径 4 mm 试件；最小误差 −3.37%，出现在孔径 2 mm 试件。2D 模型的最大误差为 10.89%（D = 8 mm），最小误差 5.88%（D = 6 mm）。图 23 给出 2D、3D 数值与试验结果中极限应力随开孔直径（D）的关系。尽管五个试件的极限强度预测存在偏差，但 3D 与 2D 模型得到的强度随层合板尺寸变化的趋势与试验结果完全一致：试件尺寸增大，极限强度单调下降。数值模拟（2D、3D）与试验都表明，极限强度差异最大的是孔径 2 mm 与 4 mm 的层合板；且极限强度的变化量随构件尺寸增大而减小，说明尺寸效应随尺寸增大而减弱；当 D = 8 mm 时，层合板基本不受尺寸效应影响。
表 9 五个试件的拉伸失效应力（试验与数值对比）
图 23 试验与数值结果的尺寸效应对比。
值得注意的是：2D 模型预测极限强度的平均绝对误差为 9.28%，高于 3D 模型的 6.00%；且 2D 模型的极限强度预测一致高于试验值。原因有二：（a）2D 模型不考虑面外应力分量，使软化段出现得偏「晚」；（b）2D 模型无法同时计入同方向因非单调加载引起的拉、压损伤，导致模型极限承载力偏高（对比讨论见 4.1.2 节）。
因此，与 ABAQUS 材料库的 2D Hashin 损伤模型相比，本文 3D PDM 通过考虑厚度方向应力分量以及非单调加载引起的刚度退化，能更准确地预测开孔板受拉的极限强度与尺寸效应。
此外，图 24 给出开孔层合板在试验（图 24a）与数值模拟（图 24b）中的极限失效形态。对比可见：两者在最终失效时都出现穿过中央圆孔的横向裂纹带。说明本文 3D PDM 能合理预测开孔层合板的极限失效模式。
图 24 [90/0/±45]3S 层合板的极限失效：（a）试验结果 [41]；（b）3D 数值结果。（原图 24 跨页，此处拆为 a/b 两张）
基于数值模拟结果与试验发现的对比讨论可以得出结论：3D PDM 能准确预测含缺口层合板受拉时的力学性能与损伤机制，并反映尺寸效应。
5. 结论
本研究通过 VUMAT 子程序在 ABAQUS 中集成了基于 3D Hashin 失效准则的渐进损伤模型（PDM），并开展了从简单单单元级到复杂层合板级的全面验证。主要结论如下。
子程序引入了一种获取等效损伤初始应变与应力的新方法，能更准确地刻画复合材料在复杂应力状态下的力学性能与损伤演化。
在单单元尺度，模型精确表达了单轴应力状态诱发的四种失效模式：纤维拉伸损伤、纤维压缩损伤、基体拉伸损伤与基体压缩损伤。此外，在单轴非单调加载下，该 PDM 有效关联了同方向的拉伸与压缩损伤，结果合理。在多单元尺度，模型消除了单元尺寸对模拟结果的影响，准确反映了弹性段与损伤演化段的力学行为。
特征长度对准确表达复合材料力学行为影响显著。对三维实体单元，特征长度应取单元在损伤模式对应方向上的几何长度。选对特征长度可消除网格尺寸对模拟结果的影响。
模型在层合板尺度同样达到预期效果：准确预测了 0°、±45°、90° 单层板的力学行为。最终失效时所有铺层方向的裂纹均沿纤维方向；90° 与 ±45° 铺层的最终失效主要由基体拉伸失效造成，而 0° 铺层由基体与纤维拉伸失效共同造成。
增强型 3D PDM 能准确预测不同尺寸、[90/0/±45]3S 铺层开孔层合板受拉时的力学行为与最终失效模式。3D 模型模拟极限强度的平均误差（6.00%）小于采用 ABAQUS 内置 2D Hashin 损伤准则的 2D 模型（9.28%）。该 3D 模型还能准确反映开孔层合板的尺寸效应：随层合板尺寸增大，尺寸效应逐渐减弱。整体失效模式的观察显示：失效前裂纹扩展有限，除 90° 铺层外，远离孔区的铺层无可见裂纹；加载过程中 0° 铺层主要承担拉伸载荷，其失效最终导致层合板失效；最终失效时所有铺层都出现垂直于加载方向的横向断裂带，断裂带两侧形成近似三角形的损伤区，与试验观察一致。
与 ABAQUS 材料库内置 2D 模型相比，本文提出的 3D PDM 能更准确、更真实地表达 FRP 材料在空间应力状态下的力学行为与损伤机制。但受篇幅限制，本研究未把 3D PDM 应用于空间载荷条件下的低速冲击等场景。后续工作将把 3D PDM 应用于更广泛的空间力学场景，如冲击损伤与空间构型复合材料模拟，以进一步验证和完善该 PDM 的精度与优越性，期待其在今后的复合材料有限元模拟中发挥重要作用。
致谢与其他声明（节译）
• 致谢： 作者感谢 VUB（布鲁塞尔自由大学）P. De Baets 与 W. Devillé 先生在试验方面的支持与协助。
• 作者贡献： 概念化，Yichen Zhang 与 Wouter De Corte；方法学、软件、验证、数据分析与论文初稿，Yichen Zhang；指导与论文修改，Wim Van Paepegem 与 Wouter De Corte。所有作者已阅读并同意稿件出版版本。
• 数据可用性： 支撑本研究所报结果的数据可向通讯作者索取。
• 利益冲突： 作者声明无利益冲突。
• 基金： 本研究未获得外部资助。Yichen Zhang 感谢中国留学基金委员会（CSC）的博士生奖学金资助（编号 202108440181）。
译者注：Hashin 准则与损伤演化的实现要点（译者归纳，非原文内容）
以下为译者面向「在 Abaqus VUMAT 里实现 3D Hashin 渐进损伤」这一主题提炼的要点清单，便于日后与 wiki/constitutive/ 下 MML、UMMDp 等塑性本构文章互参：
四种模式一把抓： 3D Hashin 把损伤拆成 FT / FC / MT / MC 四个局部变量 （式 (11)–(14) 判起始，式 (15) 管演化），再按式 (3) 组装成 、 、 三个全局变量进入刚度矩阵（式 (2)、(8)、(9)）。与 2D Hashin 的本质差异：MT/MC 用 的正负判方向，天然覆盖厚度方向分量；MC 判据里 一项不可漏。
等效初始应变不要预计算（本文核心创新）： 通行做法是 （式 (16)），复杂应力状态下会漏掉剪切分量的贡献（如 ±45° 层受拉时 主导）。正确做法：在 首次等于 1 的那一刻，现场记录 与 （式 (18)–(21)），把损伤起始与软化段无缝衔接，使每个单元的初始等效应变随真实应力状态而定。
损伤演化用断裂能标定： 线性软化下 是等效应变的线性函数（式 (15)）；最终等效应变由 反推（用 而非 ，见式 (18)–(21)），这是把网格依赖性从解中挤出去的关键——软化段曲线下面积恒等于 。
特征长度必须自己给对： ABAQUS 实体单元默认 ，与本构假设不匹配。应按加载（损伤）方向的几何长度赋值；本文取厚度方向=实际厚度、面内=√面积。取值上限受式 (23) 约束 （保证起始应变 < 最终应变，否则单元软化前弹性能不足，会出现 snap）。4.2 节的七个 Case 证明： 给错，失效位移随 系统性漂移；给对则各网格密度曲线重合。
非单调加载靠历史最大等效应变： 每个材料点存 （SDV），只有 时才更新 并刷新历史最大值；卸载/反向加载时 冻结。这保证同方向的拉损伤与压损伤互不覆盖（拉→压→拉全过程刚度退化单调累积）。ABAQUS 内置 2D Hashin 做不到这一点，会系统性高估非单调载荷下的刚度与强度（4.1.2 节的量化对比：同一单元反向加载后 2D 模型模量回弹到初始值，3D VUMAT 保持退化值）。
剪切残余刚度由 、 控制： 取 0.9 / 0.5，意味着基体拉/压失效不会把剪切刚度打成 0（4.1.1 节：失效后仍保留 10% 残余剪切刚度），对数值稳定性很重要。
SDV 布局约定： 文中 SDV1=纤维拉伸损伤变量、SDV3=基体拉伸损伤变量（表 7 的四行损伤演化曲线即按此观察），实现时可据此规划 STATEV/SDV 排布（与 MML/UMMDp 的 SDV 索引表同属一类工程问题）。
单元与分层： 实体用 C3D8R（减缩积分+沙漏控制），2D 对照只能用 SC8R；层间分层不在 VUMAT 内实现，而是用 surface-based cohesive（双线性 traction-separation）+ 二次名义应力起始准则（式 (24)）+ B-K 混合模量扩展准则，配合通用接触。
数值实践参数： 开孔层合板网格收敛到 0.5 mm（1.5→0.5 mm 极限强度从 466.6 升到 537.0 MPa）；质量缩放因子取 1000 已与 100/500 无差别；W/D=6 的五个开孔试件 3D 模型平均绝对误差 6.00%（2D 为 9.28%，且系统性偏高——因为它漏了面外分量与非单调退化）。
适用边界： 模型未涉及低速冲击等强动态场景（作者列为后续工作）；文中 E-glass/Polyester 的单元级材料参数为假定值，IM7/8552 参数取自文献 [44]，复现时注意区分。
参考文献
（共 50 条，保留英文原文以便检索；与原文编号一致）
Elchalakani M., Higgoda T.M., Kimiaei M., Yang B. Pultruded FRP Tubular Section Members, Connections, and Frames: A Review on Experimental Studies. Structures. 2023;56:104986. doi: 10.1016/j.istruc.2023.104986. Al-Furjan M.S.H., Shan L., Shen X., Zarei M.S., Hajmohammad M.H., Kolahchi R. A Review on Fabrication Techniques and Tensile Properties of Glass, Carbon, and Kevlar Fiber Reinforced Rolymer Composites. J. Mater. Res. Technol. 2022;19:2930–2959. doi: 10.1016/j.jmrt.2022.06.008. Jarrett W., Jeffs S.P., Korkees F., Rawson M. The Opportunities and Challenges of Hybrid Composite Driveshafts and Their Couplings in the Aerospace Industry: A Review. Compos. Struct. 2023;320:117203. doi: 10.1016/j.compstruct.2023.117203. Monaldo E., Nerilli F., Vairo G. Basalt-Based Fiber-Reinforced Materials and Structural Applications in Civil Engineering. Compos. Struct. 2019;214:246–263. doi: 10.1016/j.compstruct.2019.02.002. Velea M.N., Wennhage P., Zenkert D. Multi-Objective Optimisation of Vehicle Bodies Made of FRP Sandwich Structures. Compos. Struct. 2014;111:75–84. doi: 10.1016/j.compstruct.2013.12.030. Sreejith M., Rajeev R.S. Fiber Reinforced Composites. Elsevier; Amsterdam, The Netherlands: 2021. Fiber Reinforced Composites for Aerospace and Sports Applications; pp. 821–859. Topol H., Al-Chlaihawi M.J., Demirkoparan H., Merodio J. Bifurcation of Fiber-Reinforced Cylindrical Membranes under Extension, Inflation, and Swelling. J. Appl. Comput. Mech. 2023;9:113–128. doi: 10.22055/jacm.2022.40949.3677. Karim M.A., Abdullah M.Z., Deifalla A.F., Azab M., Waqar A. An Assessment of the Processing Parameters and Application of Fibre-Reinforced Polymers (FRPs) in the Petroleum and Natural Gas Industries: A Review. Results Eng. 2023;18:101091. doi: 10.1016/j.rineng.2023.101091. Oya N., Hamada H. Mechanical Properties and Failure Mechanisms of Carbon Fibre Reinforced Thermoplastic Laminates. Compos. Part A Appl. Sci. Manuf. 1997;28:823–832. doi: 10.1016/S1359-835X(97)00035-3. Devarajan B., Lakshminarasimhan R., Murugan A., Rangappa S.M., Siengchin S., Marinkovic D. Recent Developments in Natural Fiber Hybrid Composites for Ballistic Applications: A Comprehensive Review of Mechanisms and Failure Criteria. Facta Univ. Ser. Mech. Eng. 2024:343–383. doi: 10.22190/FUME240216037D. Hashin Z. Failure Criteria for Unidirectional Fiber Composites. J. Appl. Mech. 1980;47:329–334. doi: 10.1115/1.3153664. Puck A., Schürmann H. Failure Analysis of FRP Laminates by Means of Physically Based Phenomenological Models. Compos. Sci. Technol. 2002;62:1633–1662. doi: 10.1016/S0266-3538(01)00208-1. Pinho S., Darvizeh R., Robinson P., Schuecker C., Camanho P. Material and Structural Response of Polymer-Matrix Fibre-Reinforced Composites. J. Compos. Mater. 2012;46:2313–2341. doi: 10.1177/0021998312454478. Rahul K., Yudhanto A., Ravindranath P.K., Jack D.A. Progressive Damage Modeling in Open Hole Composite Laminates with Ultrasound-Informed Drilling-Induced Delamination. Compos. Part A Appl. Sci. Manuf. 2024;184:108262. doi: 10.1016/j.compositesa.2024.108262. Divse V., Marla D., Joshi S.S. Progressive Damage Analysis in an Open Hole Compression of FRP Laminates Including Fiber Kinking and Pre-Existing Damage. Compos. Part A Appl. Sci. Manuf. 2023;169:107523. doi: 10.1016/j.compositesa.2023.107523. Divse V., Marla D., Joshi S.S. 3D Progressive Damage Modeling of Fiber Reinforced Plastics Laminates Including Drilling-Induced Damage. Compos. Part A Appl. Sci. Manuf. 2022;163:107230. doi: 10.1016/j.compositesa.2022.107230. Trombini M., Enea M., Arruda M.R.T., Pagani A., Petrolo M., Carrera E. 1D Higher-Order Theories for Quasi-Static Progressive Failure Analysis of Composites Based on a Full 3D Hashin Orthotropic Damage Model. Compos. Part B Eng. 2024;270:111120. doi: 10.1016/j.composites b.2023.111120. Pham D.C., Cui X., Ren X., Lua J. A Discrete Crack Informed 3D Continuum Damage Model and Its Application for Delamination Migration in Composite Laminates. Compos. Part B Eng. 2019;165:554–562. doi: 10.1016/j.composites b.2019.02.045. Huang L., Tao Y., Sun J., Zhang D., Zhao J. Assessment of Numerical Modeling Approaches for Thin Composite Laminates under Low-Velocity Impact. Thin-Walled Struct. 2023;191:111053. doi: 10.1016/j.tws.2023.111053. Donadon M.V., Iannucci L., Falzon B.G., Hodgkinson J.M., De Almeida S.F.M. A Progressive Failure Model for Composite Laminates Subjected to Low Velocity Impact Damage. Comput. Struct. 2008;86:1232–1252. doi: 10.1016/j.compstruc.2007.11.004. Xin S.H., Wen H.M. A Progressive Damage Model for Fiber Reinforced Plastic Composites Subjected to Impact Loading. Int. J. Impact Eng. 2015;75:40–52. doi: 10.1016/j.ijimpeng.2014.07.014. Kachanov M. On the Concept of Damage in Creep and in the Brittle-Elastic Range. Int. J. Damage Mech. 1994;3:329–337. doi: 10.1177/105678959400300402. Matzenmiller A., Lubliner J., Taylor R.L. A Constitutive Model for Anisotropic Damage in Fiber-Composites. Mech. Mater. 1995;20:125–152. doi: 10.1016/0167-6636(94)00053-0. Lapczyk I., Hurtado J.A. Progressive Damage Modeling in Fiber-Reinforced Materials. Compos. Part A Appl. Sci. Manuf. 2007;38:2333–2341. doi: 10.1016/j.compositesa.2007.01.017. Maimí P., Camanho P.P., Mayugo J.A., Dávila C.G. A Continuum Damage Model for Composite Laminates: Part II—Computational Implementation and Validation. Mech. Mater. 2007;39:909–919. doi: 10.1016/j.mechmat.2007.03.006. Fakoor M., Mohammad Navid Ghoreishi S. Experimental and Numerical Investigation of Progressive Damage in Composite Laminates Based on Continuum Damage Mechanics. Polym. Test. 2018;70:533–543. doi: 10.1016/j.polymertesting.2018.08.013. Ridha M., Wang C.H., Chen B.Y., Tay T.E. Modelling Complex Progressive Failure in Notched Composite Laminates with Varying Sizes and Stacking Sequences. Compos. Part A Appl. Sci. Manuf. 2014;58:16–23. doi: 10.1016/j.compositesa.2013.11.012. Su Z.C., Tay T.E., Ridha M., Chen B.Y. Progressive Damage Modeling of Open-Hole Composite Laminates under Compression. Compos. Struct. 2015;122:507–517. doi: 10.1016/j.compstruct.2014.12.022. Li X., Ma D., Liu H., Tan W., Gong X., Zhang C., Li Y. Assessment of Failure Criteria and Damage Evolution Methods for Composite Laminates under Low-Velocity Impact. Compos. Struct. 2019;207:727–739. doi: 10.1016/j.compstruct.2018.09.093. Maimi P., Mayugo J.A., Camanho P.P. A Three-Dimensional Damage Model for Transversely Isotropic Composite Laminates. J. Compos. Mater. 2008;42:2717–2745. doi: 10.1177/0021998308094965. Lee C.-S., Kim J.-H., Kim S., Ryu D.-M., Lee J.-M. Initial and Progressive Failure An alyses for Composite Laminates Using Puck Failure Criterion and Damage-Coupled Finite Element Method. Compos. Struct. 2015;121:406–419. doi: 10.1016/j.compstruct.2014.11.011. Liao B.B., Liu P.F. Finite Element An alysis of Dynamic Progressive Failure of Plastic Composite Laminates under Low Velocity Impact. Compos. Struct. 2017;159:567–578. doi: 10.1016/j.compstruct.2016.09.099. Liu Y., Hou Y., Sapanathan T., Meng L., Xu Y. Multiscale Modeling of the Mechanical Behavior of 3D Braided CFRP Composites under Uniaxial Tension. Compos. Struct. 2023;306:116601. doi: 10.1016/j.compstruct.2022.116601. Zhang D., Zhou J., Wang J., Zhang W., Guan Z. A Comparative Study on Failure Mechanisms of Open-Hole and Filled-Hole Composite Laminates: Experiment and Numerical Simulation. Thin-Walled Struct. 2024;198:111730. doi: 10.1016/j.tws.2024.111730. Ait Mohammed M., Tarfaoui M. A Progressive Damage Modelling of Glass/Epoxy Cylindrical Structure Subjected to Low-Velocity Impact. Eng. Fail. A nal. 2022;134:106036. doi: 10.1016/j.engfailana l.2022.106036. Ni K., Chen Q., Wen J., Cai Y., Zhu Z., Li X. Low-Velocity Impact and Post-Impact Compression Properties of Carbon/Glass Hybrid Yacht Composite Materials. Ocean Eng. 2024;292:116448. doi: 10.1016/j.oceaneng.2023.116448. Baant P., Oh B.H. Crack Band Theory for Fracture of Concrete. Matériaux Constr. 1983;16:155–177. doi: 10.1007/BF02486267. De Maio U., Greco F., Lonetti P., Pranno A. A Combined ALE-Cohesive Fracture Approach for the Arbitrary Crack Growth A nalysis. Eng. Fract. Mech. 2024;301:109996. doi: 10.1016/j.engfracmech.2024.109996. Kuhn C., Müller R. A Continuum Phase Field Model for Fracture. Eng. Fract. Mech. 2010;77:3625–3634. doi: 10.1016/j.engfracmech.2010.08.009. Benzeggagh M.L., Kenane M. Measurement of Mixed-Mode Delamination Fracture Toughness of Unidirectional Glass/Epoxy Composites with Mixed-Mode Bending Apparatus. Compos. Sci. Technol. 1996;56:439–449. doi: 10.1016/0266-3538(96)00005-X. Camanho P.P., Maimí P., Dávila C.G. Prediction of Size Effects in Notched Laminates Using Continuum Damage Mechanics. Compos. Sci. Technol. 2007;67:2715–2727. doi: 10.1016/j.compscitech.2007.02.005. Maimí P., Camanho P.P., Mayugo J.A., Dávila C.G. A Continuum Damage Model for Composite Laminates: Part I—Constitutive Model. Mech. Mater. 2007;39:897–908. doi: 10.1016/j.mechmat.2007.03.005. Falzon B.G., Apruzzese P. Numerical An alysis of Intralaminar Failure Mechanisms in Composite Structures. Part II: Applications. Compos. Struct. 2011;93:1047–1053. doi: 10.1016/j.compstruct.2010.06.022. Reinoso J., Catalanotti G., Blázquez A., Areias P., Camanho P.P., París F. A Consistent Anisotropic Damage Model for Laminated Fiber-Reinforced Composites Using the 3D-Version of the Puck Failure Criterion. Int. J. Solids Struct. 2017;126–127:37–53. doi: 10.1016/j.ijsolstr.2017.07.023. Völkerink O., Petersen E., Koord J., Hühne C. A Pragmatic Approach for a 3D Material Model Considering Elasto-Plastic Behaviour, Damage Initiation by Puck or Cuntze and Progressive Failure of Fibre-Reinforced Plastics. Comput. Struct. 2020;236:106280. doi: 10.1016/j.compstruc.2020.106280. Divse V., Marla D., Joshi S.S. Finite Element An alysis of Tensile Notched Strength of Composite Laminates. Compos. Struct. 2021;255:112880. doi: 10.1016/j.compstruct.2020.112880. Prabhakar P., Waas A.M. A Novel Continuum-Decohesive Finite Element for Modeling in-Plane Fracture in Fiber Reinforced Composites. Compos. Sci. Technol. 2013;83:1–10. doi: 10.1016/j.compscitech.2013.03.022. Khechai A., Tati A., Guerira B., Guettala A., Mohite P.M. Strength Degradation and Stress An alysis of Composite Plates with Circular, Square and Rectangular Notches Using Digital Image Correlation. Compos. Struct. 2018;185:699–715. doi: 10.1016/j.compstruct.2017.11.060. Forghani A., Poursartip A., Vaziri R. An Orthotropic Non-Local Approach to Modeling Intra-Laminar Damage Progression in Laminated Composites. Int. J. Solids Struct. 2019;180–181:160–175. doi: 10.1016/j.ijsolstr.2019.07.015. Jalalvand M., Czél G., Wisnom M.R. Damage An alysis of Pseudo-Ductile Thin-Ply UD Hybrid Composites—A New An alytical Method. Compos. Part A Appl. Sci. Manuf. 2015;69:83–93. doi: 10.1016/j.compositesa.2014.11.006.
