随着纤维增强树脂基复合材料(Fiber reinforced plastic, FRP)结构设计水平和制造能力的提升，FRP在结构上的应用由少及多、由次要变主要、由非承力到主承力、由静力覆盖到考虑冲击损伤及疲劳问题。FRP已成为“先进结构”不可或缺的部分。

科学基本原理大多源于试验、观察于试验、论证于试验，但试验费时、费力、费财。数值仿真尽管不能完全代替试验，但可为FRP结构性能评估和优化设计提供依据，从而大幅减少试验次数，提高试验的经济和时间效率。

诸多学者在FRP数值仿真领域已经获得积极的科研结果，Rose^[1]较为全面地综述了FRP结构渐进损伤数值分析的若干问题及其研究进展。目前仿真领域的研究热点将趋向集中于两个方面，一是研究复杂载荷工况下的数值模拟。例如文献[2~4]数值分析了复材的冲击响应历程及其冲击后损伤。文献[5, 6]中仿真了复材的冲击后压缩响应。二是着重研究机理更为复杂(相比层内响应)的层间分层扩展响应。例如文献[7, 8]中提出一种新的层间本构，用于模拟混合模式下的分层扩展。Harper^[9]研究粘聚区长度对分层模拟效果的影响。Turon^[10]数值分析了疲劳载荷下复材分层损伤的起扩展情况。

本文在多年研究工作的基础上，提出一种适用于模拟FRP综合力学性能的数值仿真架构，并总结其关键要素，最后给出多种载荷工况下的数值分析算例。

1 材料损伤的定义

FRP力学性能仿真本质上与载荷型式无关，其核心在于构造一个数值仿真架构来定量描述FRP的力学行为响应。力学行为响应是能量传递的过程，材料在吸收外载输入的能量时，其内部便伴随着损伤的起始与演化。FRP力学性能分析亦即等价于确定材料内部损伤的形成机理及其累积信息。

物理机理上，损伤的定义形式十分多样。宏观上有刚度、强度、寿命，细观上有局部硬度、透波性，微观上有裂纹密度等。损伤状态下各种性能参数与完好状态下相应参数的比值，皆可定义为损伤D

$ D = 1 - \frac{{X\left( t \right)}}{{{X_0}}} $

(1)

式中:X(t)为t时刻的材料性能参数; X₀为对应的初始性能参数。通常地，我们多用材料的刚度性能作为X来定义损伤。

目前的数值仿真中，该性能参数通常选择为刚度(损伤前)或者能量释放率(损伤后)。如可以将损伤定义为

$ D = 1 - \frac{{{G^{\rm{d}}}}}{{{G_{\rm{c}}}}} $

(2)

式中: G^d为损伤状态时的能量释放率; G_c为临界能量释放率。

2 FRP仿真的物理架构

FRP仿真的物理架构流程图如图 1所示，其涉及5个物理场，包含4类核心模块。

2.1 物理场

物理场分为两类，一类是载荷场(代表外部信息)，另一类是材料场(代表内部信息)。

载荷场是多种载荷的线性叠加，其作用于FRP结构后，材料内部产生相应顺序的损伤累积。损伤不可逆原则意味着先期载荷相当于给后期载荷引入了不同程度的初始损伤。本征场主要包含FRP初始几何参数，它代表结构材料最基本的分布信息。损伤场与属性场都是本征场的另一类表征，属性场基于力学性能在本征场的基础上进一步丰富材料信息，损伤场基于损伤定义在本征场的基础上建立了无量纲信息场。应力应变场是整个循环的运行依据，仿真架构的核心模块集中于此。

通过力学分析，基于载荷场和力学属性场，可导出应变场。进而基于材料本构，推得应力场。材料完好本构导出的应力场配合损伤起始准则判断损伤萌生，材料损伤本构导出的应力场配合损伤扩展准则判断损伤质变。若未满足起始准则，损伤状态为0；若首次满足起始准则，损伤状态由式(1)确定；若已有损伤(包含初始损伤)且未满足损伤扩展准则，损伤状态由式(2)确定；若满足损伤扩展准则，则损伤质变为完全失效，即损伤状态为1。损伤状态的确定更新了损伤场，若损伤场出现贯穿性的失效材料分布，则整个架构停止分析；反之，则进行材料信息的更新与传递。对于完全失效的材料，将其等效为空洞，并由本征场进行拓扑信息的更新；对于部分失效的材料，由属性场结合损伤定义，对力学参数进行更新，更新完成后直接执行当前载荷步下的新一轮内循环分析。若所有材料均未损伤萌生，则由载荷场增加载荷增量步，进行新一轮的外循环分析。

2.2 核心模块

各物理场之间存在相互交错的个体关系及内外两个循环的整体关系。众多关系都会最终汇集到若干核心模块上：材料点本构、损伤起始准则、损伤演化及其扩展准则、初始材料属性场的建立(尺度效应与就位效应)。

2.2.1 材料点本构

材料区分层内本构及层间本构，层内(层间)本构主要通过应力应变(应力位移)关系进行描述。试验测得的本构曲线往往包括3个阶段：线性上升段、非线性上升段及非线性下降段。仿真模型常用的本构曲线为双线性本构，即包含线性上升段及线性下降段。本文认为这两种本构都不是真实的材料点本构，综合材料物理机理及定量描述的可行性，材料本构应该包含非线性上升段及线性下降段^[11]，如图 2所示。

材料点本构中的非线性与微损伤无关，它主要取决于材料本身的剪切非线性程度。因为材料点的尺度基本与微损伤相当，其本构并不包含微损伤非线性累积引起的非线性上升段，而是直接在损伤后进入下降段。至于下降段是线性还是非线性，目前尚不易通过试验验证，其中线性下降假设使用频次最多，模型参数的确定也更为便捷。

2.2.2 损伤起始准则

损伤起始准则分为层内和层间两部分。层内损伤起始准则又细分为单向单模准则(最大应力/应变准则)、多向单模准则(Tsai-Hill^[12]及Tsai-Wu^[13]准则)、多向多模准则(Puck^[14]及LaRc^[15]准则)。层间损伤准则也可细分为层间单向准则(最大层间正/剪应力准则)、层间多向准则(Brewer-Lagace^[16]准则及Yeh^{[4, 17]}准则)、层内层间耦合准则。

文献[18]对众多准则进行了测评，最终优选Puck准则(式3)及LaRC准则为层内损伤起始准则、优选改进的Yeh准则^[19](式4)为层间损伤起始准则。式(3, 4)参数解释具体参见文献[18, 19]。

$\small{ \left\{\begin{array}{*{20}{l}} \frac{1}{{{\varepsilon _{1T}}}}\left( {{\varepsilon _1} + \frac{{{\nu _{f12}}}}{{{E_{f1}}}}{m_{\sigma f}}{\sigma _2}} \right) = 1&{\rm{FF\;tensile}}\\ \frac{1}{{{\varepsilon _{1C}}}}\left| {\left( {{\varepsilon _1} + \frac{{{\nu _{f12}}}}{{{E_{f1}}}}{m_{\sigma f}}{\sigma _2}} \right)} \right| = 1 - {\left( {10{\gamma _{21}}} \right)^2}&{\rm{FF\;compression(kinking)}}\\ \sqrt {{{\left( {\frac{{{\tau _{21}}}}{{{S_{21}}}}} \right)}^2} + {{\left( {1 - p_{ \bot //}^{\left( + \right)}\frac{{{Y_T}}}{{{S_{21}}}}} \right)}^2}{{\left( {\frac{{{\sigma _2}}}{{{Y_T}}}} \right)}^2}} + p_{ \bot //}^{\left( + \right)}\frac{{{\sigma _2}}}{{{S_{21}}}} = 1 - \left| {\frac{{{\sigma _1}}}{{{\sigma _{1D}}}}} \right|&{\rm{IFF}}\;{\sigma _2} \ge 0\\ \frac{1}{{{S_{21}}}}\left( {{{\sqrt {\tau _{21}^2 + \left( {p_{ \bot //}^{\left( - \right)}{\sigma _2}} \right)^2} }} + p_{ \bot //}^{\left( - \right)}{\sigma _2}} \right) = 1 - \left| {\frac{{{\sigma _1}}}{{{\sigma _{1D}}}}} \right|&{\rm{IFF}}\;{\sigma _2} < 0\;{\rm{and}}\;0 \le \left| {\frac{{{\sigma _2}}}{{{\tau _{21}}}}} \right| \le \frac{{R_{ \bot \bot }^A}}{\left|\tau _{21c}\right|}\\ \left[ {{{\left( {\frac{{{\tau _{21}}}}{{2\left( 1 + p_{ \bot \bot }^{\left( - \right)}\right){S_{21}} }}} \right)}^2} + {\rm{ }}{{\left( {\frac{{{\sigma _2}}}{{{Y_C}}}} \right)}^2}} \right]\frac{{{Y_C}}}{{ -( {\sigma _2})}} = 1-\left| {\frac{{{\sigma _1}}}{{{\sigma _{1D}}}}} \right|&{\rm{IFF}}\;{\sigma _2} < 0\;\;{\rm{and}}\;\;0 \le |\frac{{{\tau _{21}}}}{{{\sigma _2}}} |\le \frac{{\left| {{\tau _{21c}}} \right|}}{{R_{ \bot \bot }^A}} \end{array} \right.} $

(3)

$ \left\{ \begin{array}{*{20}{l}} {\left( {\frac{{{\hat t_n}}}{{{R_n}}}} \right)^2} + {\left( {\frac{{{\hat t_s}}}{{{R_s}}}} \right)^2} + {\left( {\frac{{{\hat t_t}}}{{{R_t}}}} \right)^2} = 1&{\hat t_n} \ge 0\\ {\left( {\frac{{{\hat t_s}}}{{{R_s} - p_{ns}^C{\hat t_n}}}} \right)^2} + {\left( {\frac{{{\hat t_t}}}{{{R_t} - p_{nt}^C{\hat t_n}}}} \right)^2}&{\hat t_n} < 0 \end{array} \right. $

(4)

Puck准则与改进的Yeh准则均考虑了层内(间)各向应力的耦合作用，更贴近复合材料力学各向异性的行为机理。Puck准则针对材料组分的不同，区分纤维损伤(Fiber fracture, FF)与纤维间损伤(Inter fiber fracture, IFF)。使得材料损伤起始点得以具象化，损伤准则中的应力不再是复材层板的宏观应力，而是具象到纤维或者纤维间的当地应力。损伤准则不再是宏观应力与宏观强度之间的不等式关系，转而基于材料损伤起始点的当地应力与当地强度建立临界包线，更符合层内材料真实的失效机理。改进的Yeh准则区分了层间正(负)法向应力对层间损伤起始造成的不同影响。法向正应力促进损伤起始，法向负应力不仅无促进作用，反而起到抑制损伤起始的作用。故而改进的Yeh准则并没有剔除法向负应力，而是将其作为损伤阻抗引入到表达式之中。此外，由于层内层间的耦合作用，可通过数值节点间的传力进行表征，理论上行之有效的层内层间耦合应力并未引入至层内(层间)损伤起始准则。

2.2.3 损伤演化及其扩展准则

损伤演化很大程度上取决于材料点本构的假设以及损伤定义的方式，如图 3所示。整个损伤演化历程中有两个临界点，第一个临界点对应于损伤起始，材料满足损伤起始准则后，需基于式(1)量化起始损伤。第二个临界点对应于损伤扩展，在两个临界点之间，基于式(2)的损伤定义以及材料点损伤本构的假设，确定损伤演化曲线。

基于能量法的损伤演化依靠临界能量释放率判定损伤扩展，混合模式下的临界能量释放率一般由纯Ⅰ型、Ⅱ型临界能量释放率间接推出。常用的方法按照相应的待定参数可分为三参数模型(如B-K^[20]准则、Exp-Hackle^[21]准则)及四参数模型(Power-Law^[22]准则、Bilinear^[23]准则)。

文献[19]基于试验数据对上述准则进行了综合测评，并优选B-K准则(式(5))作为损伤演化终止准则

$ {G_{{\rm{TC}}}} = G{_{{\rm{IC}}}} + \left( {{G_{{\rm{II}}C}} - {G_{{\rm{IC}}}}} \right){\left( {\frac{{{G_{{\rm{II}}}}}}{{{G_{\rm{T}}}}}} \right)^m} $

(5)

式中:G_ⅠC与G_ⅡC分别为Ⅰ、Ⅱ型临界能量释放率；G_ⅡC/G_T为混合模式比。

B-K准则为半经验公式，虽然为三参数模型，但预测准确度接近四参数模型。且其表达式简单，参数便于确定，目前分层分析多采用B-K准则。

2.2.4 尺度效应与就位效应

材料初始属性场是当前研究领域较为忽视的问题，大部分研究直接将试验获得的层压板宏观性能作为材料初始属性，构建了一个单值化的均质场。该方法忽略了初始缺陷造成的性能差异，这种属性场内不同材料点之间的差异性与材料点尺度相关(尺度效应)，也与材料点位置相关(就位效应)。

若极度放大材料点的尺度，将整个层压板看作一个材料点，若干层压板的宏观性能可构成随机场Ξ(Λ, Ω)，Λ和Ω分别代表随机变量Ξ的均值与方差。若缩小材料点尺度，将其视作层压板中的组分，若干组分的细观性能也构成随机场ξ(λ, ω)，λ和ω分别代表随机变量ξ的均值与方差。尺度效应即在材料点尺度改变前后，定量描述随机场(Ξ与ξ)特征参数之间的演化关系。针对材料刚度，层压板弹性性能K是组分弹性性能k的均质化参量，因此有

$ {\lambda _k} = {\varLambda _K} $

(6)

但随着尺度的降低，分散性应呈增大趋势

$ {\omega _k} > {\mathit{\Omega} _K} $

(7)

针对材料强度，层压板强度性能S的试验值是组分强度性能s极小值参量^[24]，进而有

$ {\lambda _s} = \frac{{{\varLambda _S}}}{{\left( {1 + {a_{\min }}{w_s}} \right)}} $

(8)

$ {{\omega _s} > {\mathit{\Omega} _S}} $

(9)

式中，a_min为极小值分布系数。

就位效应最早由Parvizi^[25]于1978年提出，我们认为应将这个问题分为两个部分：一是相同子层位于不同位置，子层的表观强度不一致，这与子层间的约束作用相关；二是相同夹杂位于不同位置，层压板宏观力学性能不一致，这与初始缺陷的分布相关。初始属性的确定需综合考虑这两方面要素，在材料点位置改变前后，定量描述随机场(Ξ与ξ)特征参数之间的转化关系。基于反演理论^[26]，给出了Ξ(Λ, Ω)与ξ(λ, ω)之间的关系表达式

$ \left\{ \begin{array}{l} \sigma \left( {\varXi ,\mathit{\xi} } \right) = k\frac{{\rho \left( {\varXi ,\mathit{\xi} } \right)\varTheta \left( {\varXi ,\mathit{\xi} } \right)}}{{\eta \left( {\varXi ,\mathit{\xi} } \right)}}\\ \sigma \left( \mathit{\xi} \right) = \int {\sigma \left( {\varXi ,\mathit{\xi} } \right)} {\rm{d}}\varXi \end{array} \right. $

(10)

式中:ρ(Ξ, ξ)代表层压板宏观性能试验信息; Θ(Ξ, ξ)代表层压板与其内部材料点之间的物理信息; k及η(Ξ, ξ)分别为正则化常数及均质化概率分布函数。

3 FRP仿真架构的数值实现 3.1 FRP静力加载工况的力学性能预测 3.1.1 静力加载仿真模型

在静力加载工况下，图 1的物理架构可以具体化到图 4的数值仿真模型^[19]。

3.1.2 算例对比

对T700碳纤维/环氧树脂开孔试件进行缺口强度仿真分析^[19]：试件长220 mm，宽36 mm，共16铺层，名义厚度3 mm。其中，试件开孔直径分为6, 9, 12 mm三组。各组的数值模拟结果与试验结果对比如图 5~7所示，具体的仿真数据见表 1。

经对比可以发现，3类不同孔径的开孔件，其数值仿真结果均与试验数据分散带相贴合。表明静力加载工况下，仿真模型搭建合理，证明仿真架构贴近实际的物理损伤机理。

3.2 FRP冲击后压缩工况的力学性能预测 3.2.1 冲击后压缩仿真模型

在冲击及冲击后压缩的载荷工况下，图 1的物理架构又可以具化到图 8的数值仿真模型，它相当于内嵌了部分图 4所示的数值仿真模型。具体的仿真数据见表 2。

3.2.2 算例对比

对T300碳纤维/环氧树脂试件进行冲击后压缩分析^[11]。试件长150mm，宽100mm，共30铺层，名义厚度5mm。冲头采用直径16mm半球形钢冲头。冲击损伤及冲击后压缩响应的预测结果如图 9, 10所示。

经对比发现：冲击损伤面积及长短轴大小与试验结果接近，满足精度要求。冲击后的压缩响应曲线也较为合理地预测了试验结果，压缩极值误差满足精度要求。表明冲击后压缩工况下，仿真模型搭建正确，再次证明仿真架构贴近实际损伤机理。

3.3 FRP疲劳加载工况的力学性能预测 3.3.1 疲劳寿命仿真模型

在疲劳加载工况下，图 1的仿真架构可以具化为图 11所示的数值仿真模型^[24]。

3.3.2 算例对比

对玻璃纤维/环氧树脂复合材料疲劳寿命进行预测^[24]。试件长130mm，宽20mm，共8个铺层，名义厚度3mm。疲劳预测结果如图 12所示，具体的疲劳寿命数据见表 3。

经对比可以发现，FRP层压板的疲劳寿命预测值与试验结果接近，满足精度要求。这表明疲劳载荷工况下，仿真模型搭建正确，仿真架构贴近实际的疲劳损伤机理。

4 结论

本文提出一种适用于模拟FRP综合力学性能的仿真架构，对其包含的核心模块进行了讨论，给出了最佳的模块搭建方案。根据3种不同工况下的算例结果，证明该仿真架构合理且有效，并得到以下结论：

(1) 对于材料点本构，损伤前最好为非线性本构(考虑剪切非线性)，损伤后为线性软化本构。

(2) 层内损伤起始准则优选Puck准则或LaRC系列准则，层间损伤起始准则优选改进Yeh准则。

(3) 损伤扩展准则优选B-K准则。

(4) 网格尺度与单元本构及材料属性相互耦合，如何准确找到仿真精度与计算效率之间的平衡，还需要进一步研究。