纤维增强复合材料具有比强度高、比刚度大、可设计强以及良好的抗疲劳性能等优点, 在航空航天等众多领域得到广泛应用。然而复合材料复杂的损伤失效过程以及破坏模式给力学性能仿真分析以及结构应用带来巨大的难度^[1-2]。

近年来基于代表性体积单元(representative volume element, RVE)的微观力学损伤模型在预测复合材料损伤过程和强度性能中得到广泛应用^[3-5]。微观力学模型中对纤维分布形貌以及组分材料本构关系的准确描述是预测单向复合材料力学性能的关键。González等^[6]对比了6组纤维随机分布模型仿真得到的复合材料横向压缩应力-应变曲线, 结果表明纤维分布形式对损伤后的应力-应变曲线变化有一定影响。Totry等^[7-8]进一步研究表明不同随机分布形式在预测复合材料C/epoxy和C/peek横向压缩强度的偏差分别为3%和7%。上述研究表明纤维分布形式会对单向复合材料力学性能的预测精度产生影响, 需要在数值分析中加以考虑。

组分材料性能参数是影响单向复合材料力学性能的一个重要因素。一般来说组分材料中纤维与基体的性能可以通过强度试验获得, 然而界面强度很难通过传统试验方法获得, 并且对于同一材料, 界面脱黏起始点判定、试验方法以及试样几何差异等因素的不同, 同样会导致测出的数据分散性很大^[9]。当前数值分析中界面强度的选取主要有3种方案：①采用Mohr-Coulomb(M-C)弹塑性准则对基体力学行为进行描述时, 选取M-C准则中基体的黏聚力作为界面强度^{[6, 10-11]}；②采用基体强度作为界面强度^[12]；③通过分析界面强度对仿真结果的影响规律, 反演出与试验结果较为吻合的界面参数值^{[5, 13]}。然而3种方案中方案①和方案②只对于特定材料性能的仿真具有很好的精度，方案③需要进行复合材料力学性能试验, 花费较高且时间历程过长。

基于上述原因, 本文通过C++编写纤维随机分布算法, 使用python在Abaqus中参数化生成包含纤维、基体和界面的复合材料代表性体积单元模型, 采用扩展的线性Drucker-Prager(D-P)弹塑性准则^[14]和双线性内聚力模型分别描述基体和界面的本构关系, 最终建立了复合材料二维损伤微观力学模型。在此基础上研究分析纤维分布形式与界面强度对单向复合材料横向压缩强度的影响。

1 微观力学数值分析模型 1.1 微观结构中纤维分布的生成

在单向复合材料固化过程中, 纤维单丝浸润在流动性树脂中, 铺设出的单向层纤维呈随机分布。为了研究纤维随机分布对单向层力学性能的影响, 采用随机扩张算法^[15]生成纤维随机分布模型, 生成步骤如下所示:

1) 输入纤维体积含量和纤维半径, 计算出目标区域尺寸；

2) 利用随机函数生成第一根纤维的圆心坐标(x₁, y₁), 设为当前纤维；

3) 生成第二根纤维(x₂, y₂)。此点位置由(x₂, y₂)与当前纤维(x₁, y₁)之间的距离l和夹角θ决定, l和夹角θ为随机数, 满足l_min+2r < l < l_max+2r, 0≤θ≤2π。其中r为纤维半径, l_min为纤维之间的最短距离, 由经验公式或微观分布结构中获取。l_max是控制纤维体积含量的主要参数；

4) 重复步骤3), 直到在当前纤维周围不能再生成新的纤维。判断新生成纤维与已生成纤维之间距离是否小于l_min, 若小于则删去此纤维尝试重新生成。若尝试次数达到规定, 则进行步骤5)；

5) 选取第二根纤维为当前纤维, 重复步骤3)和步骤4)；

6) 对于每一根纤维都重复步骤3)~步骤5), 直至纤维体积含量满足或目标区域已填满, 结束。

本文纤维随机分布算法采用C++编程, 再通过python语言在Abaqus中参数化生成包含纤维、基体和界面的复合材料RVE模型。

1.2 FEM有限元模型

为了准确预测复合材料的宏观力学性能, RVE模型必须包含足够多的纤维。文献[6]指出RVE模型包含至少30根纤维时足以表征材料的性能。本文采用上述算法生成5组不同的纤维随机分布模型, 模型包含36根纤维, 纤维半径为5 μm, 体积含量为60%, RVE模型的尺寸为68.7 μm×68.7 μm。其中一组纤维随机分布RVE模型如图 1所示。

有限元模型采用平面应变单元CPE3划分基体和纤维网格, 以及二维单元COH2D4划分界面网格, 每个模型平均包含单元数约19 000。为了克服非线性数值求解中的收敛性问题, 基于Abaqus/explicit进行RVE模型的准静态分析, 设定模型质量放大系数为1 000, 以保证准静态求解的同时减小计算时间。

RVE模型施加周期性边界条件, 在相对平行的2条边以及角点上满足

(1)

式中，u表示边界上节点的位移向量, 下标L, R, T, B代表 4条边, LB, LT, RB, RT代表 4个角点, 如图 1所示。在Abaqus中通过“equation constraint”来实现公式(1)。同时分别约束角点RB和LT的水平自由度和竖直自由度, 固定角点LB, 以限制RVE模型的刚体运动。角点RB处施加水平位移载荷, 如图 1所示。

1.3 基体与界面本构关系

树脂基体在横向压缩载荷作用下表现出明显的弹塑性力学行为, 并且其屈服行为对静水压力敏感。结合扩展的线性D-P屈服准则和塑性准则模拟基体材料在静水压力下的塑性以及损伤退化力学行为, D-P屈服准则控制方程为

(2)

式中, p为静水压力, q为Mises等效应力, r为偏应力第三不变量, β和d分别是D-P准则的内摩擦角和黏聚力。k是三轴拉伸强度与三轴压缩强度之比, 引入了基体材料在拉伸和压缩情况下不同的屈服行为。k为0.78, β为19.11°, d为83.1 MPa。

M-C屈服准则同样对基体的弹塑性行为进行了描述, 非相关联流动法则中M-C准则参数可以通过与D-P准则的相互转换关系(3)式获得

(3)

式中，ϕ, c分别为M-C准则的内摩擦角和黏聚力。

采用塑性准则定义材料初始损伤的等效塑性应变ε^pl, 它表示为应力三轴度η的函数。η刻画了不同的三维应力状态, 取0, 1/3, -1/3分别表示纯剪切、单轴拉/压等应力状态, 对应的ε^pl分别为0.5, 0.05和0.5^[10]。初始损伤发生后, 基体表现出脆性断裂, 取断裂能G_m=0.5 kJ/m²。

纤维/基体的界面采用经典的双线性内聚力单元表征其损伤起始与退化的力学行为

(4)

(5)

(4) 式为初始损伤控制方程, (5)式为损伤扩展控制方程。T和G分别表示为界面应力与能量释放率, 下标n, t, s分别表示正向和2个切向方向。选取内聚力单元的初始刚度K=10⁸ GPa^[5]。界面正向强度T_n⁰=120 MPa, 与基体压缩强度相同, 切向强度T_s⁰=180 MPa, 正向断裂能G_n^c=3 kJ/m², 切向断裂能G_s^c=6 kJ/m²。

2 横向压缩性能计算及结果分析 2.1 RVE模型验证

为了验证所建立的RVE有限元模型在预测复合材料横向压缩强度的准确性, 选取试验中复合材料E-glass/MY750作为研究对象^[16]。文献[5]对该材料进行了有限元数值仿真, 与本文仿真结果进行对比验证。纤维E-glass弹性模量E_f=74.0 GPa, 泊松比v_f=0.2, 在横向压缩载荷作用下认为不发生破坏^{[5, 8]}。基体弹性模量E_m=3.35 GPa, 泊松比v_m=0.35, 拉伸强度σ_mt和压缩强度σ_mc分别为80 MPa和120 MPa。

图 2所示为RVE有限元模型预测的横向压缩应力-应变曲线与文献[5]仿真结果的对比。可以看出, RVE模型预测的复合材料横向压缩强度为138.8 MPa, 与文献[16]中试验值145 MPa的误差在4.3%以内, 比文献[5]仿真结果更接近试验值, 说明本文所建立的RVE模型可以很好地预测单向复合材料横向压缩强度。

本文RVE模型与文献[5]预测的失效应变都略高于试验的失效应变1.2%, 这主要是由于数值仿真假设基体的硬化模量为零, 而实际上在横向压缩载荷作用下, 单向复合材料中基体会产生塑性硬化, 延性减弱, 进而使得失效应变下降。

2.2 纤维分布对复合材料横向压缩力学性能的影响

为了研究纤维随机分布对单向复合材料横向压缩力学性能的影响, 建立5组不同纤维分布的RVE模型进行分析, 计算获得在横向压缩载荷作用下复合材料的应力-应变曲线如图 3所示。

由应力-应变曲线可以看出, 5组曲线在弹性段几乎重合, 随着初始损伤的发生, 曲线开始分离, 且随着应变的增大差别越来越明显, 最终导致压缩强度出现差异。这主要是因为各组RVE模型中纤维分布位置、尤其是纤维间最近距离分布不同, 使得初始损伤产生区域以及最终破坏路径出现不同, 从而引起复合材料应力-应变曲线以及横向压缩强度的差别。

表 1列出了5组纤维分布形式预测的复合材料横向压缩强度与试验值的误差。可以看出5组纤维分布形式的预测结果与试验值较为吻合, 最大误差为6.2%, 说明本文RVE模型尺寸以及界面强度的取值可以很好地预测单向复合材料的宏观强度性能。但不同纤维分布之间预测的横向压缩强度存在偏差, 模型RVE-2与RVE-4之间的偏差达到8.5%。由此可见为了更准确地预测单向复合材料力学性能以及损伤破坏形貌, 需要建立符合真实微观纤维分布特征的RVE模型。

2.3 界面强度对复合材料横向压缩力学性能的影响

为了对比现有研究中界面强度的3种选取方案对复合材料横向压缩力学性能的影响, 选取界面强度T_n⁰分别为:

a) M-C准则的黏聚力

根据关系(3)式可以获得为49 MPa;

b) 基体压缩强度

基体压缩强度为120 MPa;

c) 根据试验数据反演界面强度

通过选取不同界面强度, 仿真获得相应的复合材料横向压缩强度。比较发现当界面强度为140 MPa时, 复合材料横向压缩强度与试验值最为吻合, 因此选择140 MPa作为反演法得到的界面强度。

3种界面强度选取方案预测的复合材料横向压缩强度以及与试验值的误差如表 2所示。

可以看出, 界面强度为M-C准则黏聚力时对复合材料横向压缩强度的估值明显偏低。随着界面强度提高, 复合材料横向压缩强度升高, 并且逐渐逼近试验值。这说明此时在横向压缩载荷作用下, 界面损伤对复合材料横向压缩强度预测精度产生了影响。

方案c中界面强度140 MPa时预测的横向压缩强度显然与试验值的误差最小。此时复合材料压缩强度基本不再随界面强度增大而变化, 这说明实际上基体塑性损伤, 而不是界面损伤在复合材料横向压缩破坏过程中起主导作用, 这一结论也与文献[6]得到的结果相符。

图 4进一步给出了界面强度为140 MPa时, 复合材料在横向压缩载荷作用下的损伤断裂形貌。可以看出在横向压缩载荷作用下, 复合材料最终形成与加载垂直方向成48°的剪切断裂带。这与M-C弹塑性理论计算得到的基体断裂角50.8°相吻合^[5], 说明基体剪切塑性断裂是引起复合材料横向压缩破坏的主要因素, 而界面损伤只是在基体出现损伤的局部区域内由于应力集中引起的。

由上述分析可知, 复合材料在实际压缩破坏过程中, 界面损伤不是导致材料破坏主要因素, 此时继续增大界面强度仿真得到的复合材料横向压缩强度基本不再升高。因此本文采用一种忽略界面单元的有限元模型, 用此简化的无界面单元模型仿真得到的复合材料横向压缩应力-应变曲线以及损伤破坏形貌如图 5所示。可以看出, 无界面单元模型预测的复合材料横向压缩强度为142.9 MPa, 与试验值最为吻合, 误差仅为1.4%。同时在横向压缩载荷作用下复合材料最终断裂形貌与界面强度为140 MPa时的情况相近。这说明了通过准确描述基体的塑性损伤退化行为, 无界面单元可以更好地预测单向复合材料的横向压缩力学性能, 避免了数值仿真中界面参数选取的不确定性问题。

3 结论

1) 本文建立的纤维随机分布代表性体积单元微观力学模型可以准确预测单向复合材料的力学性能与损伤形貌。

2) 纤维随机分布会引起复合材料损伤初始位置与断裂路径改变, 是引起复合材料横向压缩强度差异性的一个因素。

3) 基体塑性损伤而不是界面损伤在复合材料横向压缩破坏过程中起主导作用。通过准确描述基体的塑性损伤退化行为, 无界面单元可以很好地预测单向复合材料的横向压缩力学性能, 避免了数值仿真中界面参数选取的不确定性问题。