基于随机分布代表性体积的颗粒增强复合材料拉伸断裂破坏模拟

2021-11-16丁桥华宁志华

材料科学与工程学报 2021年5期

丁桥华，宁志华

(暨南大学力学与建筑工程学院，重大灾害与控制教育部重点实验室，广东广州 510632)

1 引言

颗粒复合材料具有各组分的优点，同时通过各组分相的合理配比，又能削弱各组分的缺点，具有传统材料难以比拟的可设计性，但其复杂的细观结构和界面相会使得影响材料力学响应的因素变得更加复杂。颗粒增强复合材料的主要失效形式有颗粒-基体界面脱黏、基体断裂及颗粒断裂等。为明确其细观组分与力学性能之间的关系，需建立合理的细观力学模型来探讨其有效性能和损伤破坏机理。

在理论分析方面，ERDOGAN等[1]给出了基体裂纹位于颗粒外部、内部及贯穿时应力强度因子的精确解。九十年代初，LI等[2-3]获得了裂纹靠近直至穿透颗粒时裂尖能量释放率的精确解。魏金币[4]引入了颗粒复合材料中带有弱化表面相，建立了由基体、完整颗粒以及带有弱化表面相颗粒的理论分析模型。ROBERTA SBURLATI等[5]重点研究了由梯度界面区包围的固体球形颗粒增强的复合材料的有效体积和剪切模量。在试验研究方面，NANDY等[6]研究了氧化铝颗粒对碳化硅基体的断裂韧性的影响，并重点分析了氧化铝颗粒按照梯度形式分布对裂纹路径的影响。RANADE等[7]应用扫描电子显微镜观察并研究了不同的界面增强相对复合材料断裂模式以及断裂韧性的影响。LEGGOE等[8]采用高速光学装置研究了颗粒增强金属基复合材料中裂纹动态扩展问题。TIPPUR等[9-10]也采用光学测试法探讨动载荷作用下裂纹和刚性夹杂间的相互作用，着重讨论了颗粒的尺寸以及界面粘接强度对复合材料断裂行为的影响。郭晓磊[11]通过纳米压痕和高空间分辨率研究了体积分数为15% SiC颗粒增强的2009 Al基复合材料在不同粒径下的界面的硬度演变。

基于有限元分析的数值模拟已成为计算细观力学中的主要研究手段。WILLIAMS[12]等通过有限元方法研究了颗粒对基体内部裂纹扩展的影响，考虑了单颗粒和多颗粒规则分布两种情形。AYYAR和CHAWLA[13]对含初始裂纹的SiC-Al复合材料进行裂纹扩展分析，采用真实的颗粒分布模型。杨慧[14]基于Mori-Tanaka方法，提出一个可研究规则分布颗粒增强复合材料有效弹性性能的力学模型。翁琳[15]和焦文文[16]建立了颗粒随机分布的三维模型对颗粒增强复合材料的宏观力学行为进行研究。沈日麟[17]基于内聚力单元，建立区分颗粒复合材料各组分并能够显示模拟材料损伤演化的有限元模型。江上洋[18]基于单颗粒模型研究了含初始裂纹的颗粒增强复合材料断裂问题，分析了基体与颗粒模量比值对裂纹扩展路径和能量释放率的影响。

在上述文献的研究中，所采用的细观模型大多为单胞模型或考虑颗粒规则分布，个别研究者考虑了颗粒的随机分布，但未对其随机性进行评估或验证。为对颗粒增强复合材料的有效性能及其细观力学行为进行更为有效的分析，本研究基于计算细观力学，建立颗粒随机分布的二维代表性体积单元(Representative volume element, RVE)，并运用不同的统计分析方法对随机分布进行评估，以选取最具随机性和各向同性的细观结构。在此基础上，运用扩展有限元方法对细观尺度下陶瓷复合材料的拉伸断裂破坏进行模拟，探讨颗粒及界面相参数对裂纹扩展和材料断裂韧性的影响。

2 RVE单元的选取

在进行材料设计及构件性能优化前，通常需要对颗粒增强复合材料的有效性能进行评估。RVE和周期性单元是有效性能评估中常用的两种模型，周期性单元适用于颗粒周期性分布的情形。由于实际应用的复合材料，其颗粒分布一般不具备周期性，尤其当涉及材料塑性屈服或者损伤演化时，更不适宜采用周期性。实际的颗粒增强复合材料，其增强相通常是非均匀分布的。因此，采用颗粒随机分布的RVE才能对复合材料力学行为进行高保真度的评估。

2.1 RVE的生成

为分析数值产生的微结构，必须使用足够大的RVE尺寸来代表整体材料。RVE的特征尺寸可用变量δ表示，该变量与RVE的边长及颗粒的半径有关，可表示为：

δ=LRVE/R

(1)

式中：LRVE为RVE的边长，R为颗粒的半径。PATHAN等[19]认为，对于所有的体积分数和弹性性质来说，δ≥24已足够大，能提供稳定的力学性能并且分散性低于5%。因此本研究选择δ=40作为RVE的特征尺寸。

运用Python编写随机方法程序来生成具有周期性几何的2D RVE 的模型，如图1所示，RVE的边长L为20 mm，颗粒的半径R为0.5 mm，界面厚度为0.05 mm，即为颗粒半径的0.1倍，颗粒的面积分数为0.50。

图1 生成恒定颗粒直径的2D RVE(Vf =0.5)

对于颗粒随机分布的情形，需要对其进行统计分析，以便选取更优的分布模型作为RVE。

2.2 RVE的评估

采用颗粒随机分布的RVE对复合材料的有效性能及破坏进行数值模拟，需要对RVE进行统计评估，以选取最具随机性和各向同性的细观结构[20]。如果只对复合材料进行线弹性分析，并不一定需要对细观结构进行选择，但若涉及非线性分析或内聚力界面，则非常有必要选择适当的RVE。另一方面，如果要对材料进行有限变形、界面损伤或粘弹性等方面的有限元数值模拟，由于耗时太长，因此有必要通过统计分析，选择适当的代表性细观结构，以便既能减少数值分析的细观结构数量，又能满足材料力学响应评估的精确性[21]。

2.2.1径向分布函数径向分布函数g(r)是描述颗粒密度随到参考粒子的距离变化的函数。径向分布函数适用于评估颗粒直径为常数、颗粒体积分数介于10%～80%的RVE单元颗粒不均匀程度。g(r)的定义如下[22]：

×θ(r+Δr-rij)

(2)

其中，r为到参考粒子的距离，rij是颗粒i和j之间的距离，N是粒子数，N(N-1)/2是加权项。V是区域的体积(2D情形下为面积)，Vr=π(2r+Δr)Δr是内半径为r、外半径为r+Δr的圆环面积。θ(x)函数的表达式如下：

(3)

当颗粒i和j位于r和r+Δr之间时，θ(rij-r)θ(r+Δr-rij)等于1，否则等于0。当r充分大时，g(r)随着N的增大趋近于1。

图2显示了不同颗粒体积分数的径向分布函数。从图可见，当r较小时(r=1)，随着体积分数的增加，g(r)的峰值明显增大；随着r增加，不同体积分数的g(r)的值都趋近于1；体积分数越大，g(r)越快收敛于1，且曲线更平滑，这是由于随着体积分数的增大，基体富集区域越小，颗粒的连续性越好。由此可见，g(r)可用于评估颗粒分布的非均匀程度。

图2 不同体积分数的随机分布2D RVE及其径向分布函数

在相同体积分数下，g(r)曲线的波动能反映颗粒分布的不均匀程度。给定50%的颗粒体积分数，生成不同随机分布的RVE模型，绘制出g(r)曲线见图3。从图可见，不同模型的g(r)曲线其波动性有差异，Model 3的g(r)曲线更加平滑和低波动频率，表明Model 3的颗粒随机分布的连续性和均匀性更好。

图3 不同模型的径向分布函数(Vf=0.5)

2.2.2平均近邻距离变异系数为获得更多的颗粒随机分布统计评估数据，采用MATLAB计算颗粒的Voronoi镶嵌。同样，给定体积分数为50%，生成不同的颗粒随机分布RVE，分别计算出其Voronoi镶嵌如图4所示。图中所有坐标轴均进行了标准化，其刻度范围为(0, 1)。从图可见，相同体积分数下，不同RVE的Voronoi镶嵌并没有显著差别。

图4 不同随机分布的颗粒中心Voronoi镶嵌(Vf=0.5)

为进一步对颗粒随机分布进行量化描述，采用 YANG等[23]提出的平均近邻距离变异系数(Coefficient of variation of mean near-neighbor distance, COVd)来进行统计评估。COVd在表征颗粒分布的团聚程度及非均匀性方面具有较强的有效性和较高的灵敏度，其数值与颗粒的体积分数、颗粒的尺寸和形态基本无关。

COVd的计算表达式如下[23]：

COVd=σd/d

(4)

式中：d是所有颗粒平均近邻距离的均值，σd是平均近邻距离的标准差。

COVd值的范围是0～1，一般而言，其数值越高，颗粒分布越聚集，颗粒的分布就越不均匀；当COVd值越接近零，颗粒的分布越均匀[24]。

基于图4中三种RVE模型的Voronoi镶嵌，计算出其相应的COVd值，分别为0.5429，0.6000和0.4423。可见在图4的所有模型中，Model 3的COVd最小，表明该模型的颗粒团聚程度最小。结合三种模型的g(r)曲线(见图3)，可认为Model 3的颗粒随机分布均匀性最好，因此选取Model 3作为RVE。

2.2.3有效弹性常数分析生成颗粒随机分布的RVE后，可通过分析材料的有效弹性常数验证所选取的RVE是否具有代表性。为便于对比，选取文献[25]的单向纤维增强玻璃/环氧复合材料，其纤维和基体的弹性模量及泊松比分别为：E1f=72 GPa，E1m=3.35 GPa，v12f=0.21和v12m=0.35。采用图4中的Model 3(Vf=0.5，不含界面层)作为RVE，其尺寸为20 mm×20 mm，纤维半径R=0.5 mm。

假设复合材料沿垂直纤维方向受拉，并处于线弹性平面应变状态，纤维和基体理想粘结。对于单向长纤维增强复合材料，在平面应变状态下，可采用颗粒随机分布的2D模型作为其代表性体积单元。施加周期性边界条件，采用三节点线性平面应变单元(CPE3)。对于2D模型，由于增强纤维是连续的，在线弹性情况下，体积平均模量可用体积平均应力和相应的体积平均应变来表示，因此有效弹性模量和泊松比可分别用下式计算：

(5)

为预测弹性模量E22和泊松比ν23，可在RVE的上侧施加0.5%的应变，其正应力σ22分布如图5(a)所示；为预测E33和ν32，则在RVE的右侧施加0.5%的应变，其正应力σ33分布如图5(b)所示。从图可见，两种加载方式下的应力分布明显不同，位于纤维与基体界面的应力集中分布与加载方向一致，高值正应力出现在纤维聚集的局部区域。

图5 2D RVE沿加载方向的正应力分布云图

根据式(5)可计算出复合材料的有效弹性常数，如表1所示。从表中数据可见，预测的弹性常数与LI等[25]的结果较为吻合。需要说明的是，预测的弹性常数与HALPIN和KARDOS[26]采用Halpin-Tsai公式的预测结果相比有一定的偏差，这是由于HALPIN 和KARDOS[26]中所采用的模型颗粒是结构化排列，而非随机分布。PATHAN等[27]也曾报道过颗粒随机分布与结构化排列模型预测的有效弹性模量存在类似的偏差。

表1 随机分布2D RVE的有效弹性常数预测

此外，单向纤维增强复合材料通常表现出横观各向同性，即其弹性常数满足以下关系式：

(6)

由表1中数据易知，所预测的有效弹性常数满足式(6)，表明纤维在RVE两个横向方向上的分布均较为理想，所选取的RVE是有效的。

3 拉伸断裂破坏分析

颗粒增强复合材料有多种断裂机制，包括基体开裂，界面剥离和颗粒破坏等。微观结构和组分性质通过激活不同的断裂机制来确定整体的断裂行为，因此建立材料的抗断裂能力与细观属性之间的联系非常重要。另一方面，经典有限元方法在模拟裂纹扩展时，往往首先计算裂纹尖端的应力强度因子和应变能密度因子等，然后与实验获得的数值进行比较，进而判断裂纹是否扩展。但如果材料系统比较复杂，其断裂韧性受到多种因素影响而变得离散，通过实验来获取相关参数就较为困难。