1 引言

复合材料具有比强度高、比刚度大、耐腐蚀性强和抗冲击性良好等优点，广泛应用于在航空航天、汽车、船舶等高科技制造领域^[1]。其中，三维四向编织复合材料，由于具有多向纤维束构成的空间网状结构，消除了传统“层”的概念，面外性能大大提升^[2]，具有承受复杂载荷能力的优势，因此广泛应用于承受复杂热/机械载荷的航空航天部件。

三维四向编织复合材料由纤维束和基体组成，空间结构十分复杂。微观尺度上，一束纤维束包含上千根纤维丝以及与之相结合的基体。在细观编织尺度上，一般通过代表性体积单元（RVE）来预测其宏观等效力学性能。根据周期性结构分析和实验观测，可以建立三维编织复合材料的细观尺度RVE模型，又称单胞^[3~6]。目前，国内外已经对多尺度分析方法已经展开了一系列的研究，Hassani等^[7]总结介绍了双尺度渐进均匀化分析方法，Lions等^[8]在均匀化理论的基础上提出了小参数渐进展开双尺度分析方法，吴世平等^[9]利用小参数摄动展开技术引入了细观和宏观两种尺度坐标描述细观单胞的非均匀化构造和宏观均匀化结构。郑晓霞等^[10]详细地论述了多尺度分析方法在纤维增强复合材料弹塑地等力学性能中的研究进展。

现阶段关于三维四向编织复合材料刚度与强度的计算模型，大多是假定材料没有内部缺陷或并未考虑孔隙随机分布等因素对力学性能的影响。然而，不论是基于树脂转移塑模（RTM）工艺的树脂基复合材料，还是基于化学气相沉积（CVI）工艺的陶瓷基复合材料，在制造过程中都不可避免的会产生孔隙缺陷^[11]。树脂在复合材料中渗透不够充分时，会分别在纤维束和基体中产生孔隙缺陷，纤维束中的孔隙缺陷称为干斑，基体中的孔隙缺陷称为孔穴^[12]，两种不同的孔隙缺陷都会对复合材料的力学性能产生不同程度的影响。任明法等^[13]基于两尺度代表体元对含孔隙复合材料单层板的弹性常数进行了预测，但这仅限于对单层板的研究，并且其所建立的含孔隙RVE模型为均匀化后的模型，并未对纤维束和基体进行区分。同时，由于孔隙随机分布于均匀化的复合材料内部，而不是仅限于基体中，其假设与实际相悖。Shen等^[14]和石多奇等^[15]分别通过在基体中添加孔穴的RVE模型预测了二维编织和三维编织陶瓷基复合材料的宏观等效力学性能，但均未考虑干斑对纤维束性能的影响以及由纤维束尺度的不均匀性所产生的效应。

基于此，本文采用双尺度分析方法研究了干斑和孔穴对三维四向编织复合材料宏观等效力学性能的影响。在纤维束尺度上，采用通用单胞法来预测纤维束等效力学性能参数；在编织结构细观尺度上，利用RVE模型和细观力学有限元法预测得到其宏观等效弹性常数。将计算结果与实验结果进行对比，验证了双尺度分析方法的正确性。分别在纤维束和基体中引入干斑和孔穴，讨论了两种孔隙缺陷对三维四向编织复合材料力学性能的影响规律。

2 双尺度分析方法

尽管三维四向编织复合材料的空间结构十分复杂，但是由于其具有周期性排布特点，所以在细观编织尺度上可以将其分成许多的RVE模型^[6]。本文采用由底至上的双尺度分析方法：在微观尺度上由基体和纤维丝属性推导纤维束材料性能；然后在细观尺度上基于由基体和纤维束构成的单胞模型得到宏观材料性能，三维四向编织复合材料的双尺度分析模型如图 1所示。

2.1 纤维束通用单胞模型

三维四向编织复合材料的纤维束与单向纤维增强复合材料的性质基本相同。通常情况下，一束理想的纤维束是由500~2000根纤维丝和与之相关联的基体组成^[16]，纤维束的微观结构如图 2所示。

本文采用通用单胞法（GMC方法）获取纤维束等效弹性常数。GMC方法是一种基于均匀化理论，用于计算周期性复合材料代表性体积单元力学性能的高效算法，可以直接获得单胞的等效弹性参数。GMC方法将复合材料纤维束中的单胞划分为多个方形的子胞，本文采用的二维通用单胞模型如图 3所示。

通用单胞法基本过程主要由以下步骤组成^[17]：确定代表性体积单元，根据图 3所示的二维胞元结构确定胞元宏观应力（或应变）与微观应力（或应变）关系如下

$\left\{ \begin{array}{l} \overline \sigma = \frac{1}{{mn}}\sum\limits_{\beta {\rm{ = }}1}^{{N_\beta }} {\sum\limits_{\alpha = 1}^{{N_\alpha }} {{m_\alpha }{n_\beta }{{\overline \sigma }^{\left( {\alpha \beta } \right)}}} } \\ \overline \varepsilon = \frac{1}{{mn}}\sum\limits_{\beta {\rm{ = }}1}^{{N_\beta }} {\sum\limits_{\alpha = 1}^{{N_\alpha }} {{m_\alpha }{n_\beta }{{\overline \varepsilon }^{\left( {\alpha \beta } \right)}}} } \end{array} \right. $

(1)

式中$\bar \sigma $和$\bar \varepsilon $分别代表胞元的宏观应力和宏观应变；$m={\sum\limits_{\alpha = 1}^{{N_\alpha }} {{m_\alpha }} }$；$n=\sum\limits_{\beta {\rm{ = }}1}^{{N_\beta }} {{n_\beta }} $；${\bar \sigma ^{\left( {\alpha \beta } \right)}}$和${\bar \varepsilon ^{\left( {\alpha \beta } \right)}}$分别是子胞(αβ)的平均应力和平均应变。

利用单胞界面之间以及RVE中子胞界面之间的平均位移和平均应力的连续性条件建立方程。最早的GMC方法采用子胞平均应变作为未知量，而后Pindera和Bednarcyk^[18]提出了利用子胞平均应力作为未知量的高效GMC算法，大大减少了计算过程中所需方程的数量，本文即采用此种方法。子胞中的应变可以由宏观应变以及子胞内塑性应变确定。

${\bar \varepsilon ^{\left( {\alpha \beta } \right)}} = {A^{\left( {\alpha \beta } \right)}}\bar \varepsilon + {D^{\left( {\alpha \beta } \right)}}\varepsilon _s^I$

(2)

式中A^(αβ)和D^(αβ)分别为相应的集中系数矩阵；$\varepsilon _s^I$为所有子胞塑性应变组成的列向量。

复合材料宏观本构方程可以表示为

$\bar \sigma = {B^{\rm{*}}}\left( {\bar \varepsilon - {{\bar \varepsilon }^I}} \right)$

(3)

式中ε^I表示胞元塑性应变组成的列向量，B^*代表复合材料的等效刚度矩阵，表达式为

$ {B^{\rm{*}}} = \frac{1}{{mn}}\sum\limits_{\beta = 1}^{{N_\beta }} {\sum\limits_{\alpha = 1}^{{N_\alpha }} {{m_\alpha }{C^{\left( {\alpha \beta } \right)}}{A^{\left( {\alpha \beta } \right)}}} } $

(4)

本文根据上述方法编写通用单胞模型计算程序，进行后续计算分析。

2.2 编织结构RVE模型

选取两种具有不同编织角的RVE模型如表 1所示^[19]。通过分析编织机床的周期性运动，记录携纱器的运动轨迹，以表 1中第一种编织结构（No.1）为例，分析预成型编织件的空间拓扑结构^[20]，得到纱线在单胞中的空间位置如图 4（a）所示，图 4（b），（c），（d）分别展示了纤维束、基体以及单胞整体的几何结构模型。

选取模型材料为三维四向碳/环氧树脂编织复合材料，各组分弹性性能参数如表 2所示^[19]。（E₁₁E₂₂分别代表面外弹性模量与面内弹性模量，G₁₂G₁₃分别代表面外剪切模量与面内剪切模量，μ₁₂表示泊松比）。采用Hypermesh软件对上述单胞模型进行网格划分如图 5所示，其中的分网过程需要保证单胞相对两个面的节点具有一致性，这是为了确保在后续计算过程中周期性边界条件的顺利施加，从而有效预测单胞模型的力学性能。

2.3 预测结果

依据上述所建立的微观-细观双尺度分析方法，对三维四向编织复合材料进行宏观弹性常数预测，采用2.1中所述方法，建立纤维束的通用单胞模型，应用MATLAB自编程序，计算纤维束的力学性能参数如下：E₁₁=208.150GPa，E₂₂=E₃₃=25.966GPa，μ₁₂=μ₁₃=0.273，μ₂₃=0.378，G₁₂=G₁₃=12.942GPa，G₂₃=9.289GPa。

细观尺度上，将纤维束的力学性能参数作为均质化参数，应用到三维四向编织复合材料RVE有限元模型中。利用有限元软件ABAQUS/Standard对单胞模型施加周期性边界条件^{[21, 22]}，对其进行x，y，z方向纯拉伸和xy，yz，zx方向纯剪切的应力分析，计算过程中需要保证每种情况的应变均为单胞尺寸大小的1%。图 6为单胞模型沿z方向单轴拉伸的应力分布云图，图 7为沿yz方向纯剪切的应力分布云图。从两幅图中可以看出，载荷主要由纤维束承担，基体只承担了较少部分的应力，同时，在纤维束与基体的交界面处出现了局部应力集中现象。从图 7中可以看出，在剪切过程中，尽管单胞的横截面发生了翘曲，但是整个单胞模型的相对表面却始终保持着同样的变形，证明周期性边界条件的成功施加，达到了“周期性变形”的目的。

整体小变形情况下，应变-应力关系可表示为${\bar \zeta _{ij}} = {S_{ijkl}}{\bar \sigma _{kl}}$。其中，${\bar \zeta _{ij}}$和${\bar \sigma _{kl}}$分别是RVE单元的平均应变和平均应力，${S_{ijkl}}$是不同情况下对应的柔度矩阵。根据柔度矩阵，得到三维四向编织复合材料的力学性能参数。

依据本文提出的双尺度分析方法，三维四向碳/环氧树脂编织复合材料的弹性常数的预测结果及误差对比如表 3所示，可以发现预测值与实验值^[19]基本相同，其误差来源在于本文所建立的RVE模型并未考虑纤维束的扭曲变形等因素，但总体上数值精度拟合较好，这说明本文的双尺度分析方法可以较好地预测三维四向编织复合材料的宏观弹性常数。

3 孔隙缺陷对力学性能的影响规律 3.1 孔隙缺陷随机分布

目前工程中通常采用无损检测技术^[12]可以有效测得纤维束中孔隙率P_fv与基体中孔隙率P_mv的大小。复合材料成型工艺过程中，孔隙缺陷的大小和位置随机排布，导致在仿真过程中难以量化表征孔隙的分布形式，传统的多孔结构会使得网格划分过程非常复杂，网格划分质量较差。本文则是采用气孔单元法^[23]来模拟纤维束和基体中的孔隙在单胞模型中的存在形式。具体操作方法即在单胞模型中，通过Monte-Carlo仿真技术，按照无损检测实测得到的孔隙率，在纤维束及基体中随机投放气孔单元。其中，气孔单元视为无刚度，但为了保证计算的收敛性，不能使得每个气孔单元的刚度均为0，而是选取一个很小的正值（本文取10^-5Pa）。

本节为了消除气孔单元位置随机性所造成的影响，以表 1中第一种编织结构（No.1）为例进行了40组气孔单元的随机投放，按孔隙率P_mv=P_fv=4%，在纤维束和基体中分别投放气孔单元如图 8所示。

运用上文中所建立的双尺度分析模型，分别计算出每次投入气孔单元时三维四向碳/环氧树脂编织复合材料的弹性常数。其中，E₁₁的概率分布如图 9所示，可以看出，E₁₁的分布全部集中在13.6~14.2GPa的范围内。这说明只要孔隙率保持不变，气孔位置的随机分布对复合材料沿编织方向的宏观弹性常数影响不大。

3.2 孔隙缺陷百分比

对于复合材料而言，孔隙缺陷的体积分数对其力学性能有着重要的影响，通过本文中所建立的双尺度分析模型，对表 1中第一种编织结构（No.1）的三维四向编织复合材料进行分析，纤维束中不同孔隙率所对应的纤维束等效力学参数计算结果如表 4所示。从表中可以看出，随着P_fv的不断增加，所有的弹性常数均呈现下降的趋势；其中，除沿编织方向弹性模量E₁₁外，其余参数均随着P_fv的增加而显著减小，结果表明纤维束中孔隙的存在对于纤维束力学性能的影响很大。

图 10展示了表 1中第一种编织结构（No.1）的三维四向编织复合材料在不同纤维束孔隙率以及不同基体孔隙率下的力学性能计算结果。从图 10（a），（b）中可以看出，单胞纵向方向的弹性模量大于横向方向的弹性模量，说明三维四向编织复合材料纵向拉伸性能强于横向拉伸性能，当P_mv增加时，纵向弹性模量E₁₁与横向弹性模量E₂₂，E₃₃均减小。图 10（c），（d）中，随着P_mv的增加，横向泊松比μ₂₃逐渐减小且趋势明显，纵向两个方向的泊松比μ₁₂和μ₁₃呈现缓慢上升趋势，表明孔隙缺陷的存在，促进了三维四向编织复合材料在纵向载荷下的横向变形。图 10（e）中，随着P_mv的增加，G₁₂，G₁₃呈现与μ₁₂，μ₁₃相反的趋势。

当P_fv逐渐增大时，除了横向泊松比μ₂₃外，其它宏观力学性能参数均随着P_fv的增加而减小。对比各力学性能参数随着P_fv和P_mv的变化趋势，可以发现，P_fv对三维四向编织复合材料力学性能参数的影响作用要大于P_mv。对于弹性模量和剪切模量而言，P_fv的影响效果则更加明显。

4 结论

通过本文研究，得到如下结论：

（1）当孔隙率一定时（P_mv=P_fv=4%），沿编织方向弹性模量仅从13.6GPa变化到14.2GPa，说明孔隙位置的随机分布对沿编织方向弹性模量的影响几乎可以忽略。

（2）孔隙率对三维四向编织复合材料力学性能参数有一定的影响。无论是随着纤维束中孔隙率的增加还是基体中孔隙率的增加，除了横向泊松比外，其余的力学性能参数均呈现下降的趋势。

（3）纤维束中孔隙率对于三维四向编织复合材料力学性能的影响作用要大于基体中孔隙率对其的影响作用。