0 引言

C/SiC复合材料具有耐高温、高比强、高比模、抗氧化、抗烧蚀，对裂纹不敏感，不发生灾难性损毁等特点^[1]。

国内外对纤维增韧陶瓷基复合材料的力学行为的实验研究颇为重视，然而力学行为模拟计算上大多集中于树脂基复合材料。从T. K.Jacobsen等^[2]得到的C/SiC机织复合材料单轴拉伸应力-应变曲线中可以看出：由于C/SiC复合材料加载过程中裂纹、空洞等缺陷的逐渐扩展，其应力-应变曲线呈现伪塑性，即模量逐渐下降的非线性关系。针对树脂基复合材料，文献[3-5]基于数值方法研究了三维编织复合材料的渐进损伤，并进行了强度预测，然而，陶瓷基复合材料基体模量大、强度低，造成基体率先破坏，且基体呈现脆性，这些与树脂基复合材料呈现不同的特点。陶永强等^[6]将2D机织结构简化为0°/90°正交铺层结构，对复合材料在单轴拉伸载荷作用下的渐进损伤进行了分析。然而这一宏观平均方法并不能体现材料细观上缺陷萌生到逐渐扩展的过程，也无法体现不同材料不同的损伤特性。此外，许多模型都未能体现陶瓷基复合材料因制备工艺在纤维束内外基体中形成的缺陷，或是仅仅采用Weibull分布在外部基体中预制一些较小缺陷，这些都与实际结构有很大出入，影响了最后结果的分析。

本文基于考虑纤维束相互挤压的卵圆形纤维束截面单胞模型，引入周期性位移边界条件，采用细观多尺度有限元方法，建立了二维C/SiC平纹机织复合材料的渐进损伤多尺度模型。该模型在两个尺度上对其单胞的力学性能进行考虑：在第一尺度上建立了含预制缺陷的由SiC基体包裹的T300单丝碳纤维各向异性模型，以模拟纤维束的制备缺陷和性能退化；在第二尺度上依据第一尺度均匀化方法得到的力学性能，建立平纹机织C/SiC复合材料单胞模型，结合周期性边界条件，针对不同部分的破坏特点，分别使用Hashin准则和莫尔强度理论，通过ANSYS有限元软件对平纹机织复合材料的渐进损伤过程进行数值模拟。该模型基于纤维束和纯基体相应损伤类型所造成的材料性能退化，模拟了单调拉伸状态下损伤产生、发展的过程，并最终得到了单轴拉伸损伤状况下材料参数的变化规律。

1 能量法原理

由均匀化方法^[7]，复合材料的力学性质可由代表性体积单元体现。在弹性范围内，代表性体积单元的宏观力学行为可以由应力张量σ和应变张量ε体现，两者之间由等效刚度矩阵C^H相关联：

$ \mathit{\boldsymbol{\bar \sigma }} = {\mathit{\boldsymbol{C}}^\mathit{\boldsymbol{H}}}\mathit{\boldsymbol{\bar \varepsilon }} $

(1)

式中，$ \mathit{\boldsymbol{\bar \sigma = }}\frac{1}{V}\int\limits_\mathit{\Omega } {\mathit{\boldsymbol{\sigma }}{\rm{d}}\mathit{\Omega }}, \mathit{\boldsymbol{\bar \varepsilon = }}\frac{1}{V}\int\limits_\mathit{\Omega } {\mathit{\boldsymbol{\varepsilon }}{\rm{d}}\mathit{\Omega }} $，V为代表性体积单元的体积。Ω为代表性体积单元所占用的空间。

考虑三维各向异性材料的情形，式(1) 可表示为：

$ \left[{\begin{array}{*{20}{c}} {{{\bar \sigma }_{11}}}\\ {{{\bar \sigma }_{22}}}\\ {{{\bar \sigma }_{33}}}\\ {{{\bar \sigma }_{12}}}\\ {{{\bar \sigma }_{23}}}\\ {{{\bar \sigma }_{31}}} \end{array}} \right] = \left[{\begin{array}{*{20}{c}} {C_{1111}^H}&{C_{1122}^H}&{C_{1133}^H}&0&0&0\\ {C_{1122}^H}&{C_{2222}^H}&{C_{2233}^H}&0&0&0\\ {C_{1133}^H}&{C_{2233}^H}&{C_{3333}^H}&0&0&0\\ 0&0&0&{C_{1212}^H}&0&0\\ 0&0&0&0&{C_{2323}^H}&0\\ 0&0&0&0&0&{C_{3131}^H} \end{array}} \right]\left[{\begin{array}{*{20}{c}} {{{\bar \varepsilon }_{11}}}\\ {{{\bar \varepsilon }_{22}}}\\ {{{\bar \varepsilon }_{33}}}\\ {{{\bar \varepsilon }_{12}}}\\ {{{\bar \varepsilon }_{23}}}\\ {{{\bar \varepsilon }_{31}}} \end{array}} \right] $

(2)

和此细观结构相关的应变能可表示为：

$ \begin{array}{l} E = \int\limits_\mathit{\Omega } {\frac{1}{2}({\sigma _{11}}{\varepsilon _{11}} + {\sigma _{22}}{\varepsilon _{22}} + {\sigma _{33}}{\varepsilon _{33}} + {\sigma _{12}}{\varepsilon _{12}} + {\sigma _{23}}{\varepsilon _{23}} + {\sigma _{31}}{\varepsilon _{31}})} {\rm{d}}\mathit{\Omega = }\\ \frac{1}{2}({{\bar \sigma }_{11}}{{\bar \varepsilon }_{11}} + {{\bar \sigma }_{22}}{{\bar \varepsilon }_{22}} + {{\bar \sigma }_{33}}{{\bar \varepsilon }_{33}} + {{\bar \sigma }_{12}}{{\bar \varepsilon }_{12}} + {{\bar \sigma }_{23}}{{\bar \varepsilon }_{23}} + {{\bar \sigma }_{31}}{{\bar \varepsilon }_{31}}) \end{array} $

(3)

在如表 1所示特殊边界条件下，联立(2)、(3) 两式可以帮助推断^[8]代表性体积单元等效刚度矩阵C^H中的系数，并由此可导出材料单胞整体本构方程。同理可得其他矩阵系数。

2 纤维束模型的建立 2.1 几何模型的建立

根据图 1(a)所示^[1]的纤维束内部的电镜扫描结果，碳纤维增韧SiC复合材料(CFCC-SiC) 的内部孔隙源自于预制体结构。总体而言，CFCC-SiC内存在两个尺度的孔隙，纤维束间孔隙和纤维束内孔隙(同一纤维束内单丝纤维之间的孔隙)。据此，在纤维束尺度上，可以按照图 1(b)所示选取代表性体积元，并建立了如图 2(a)所示的纤维束模型，以此进行对纤维束初始材料常数的预测。因此，可以在第一步模型的基础上，于SiC材料中预制一如图 2(b)所示的初始缺陷。由此，可以模拟含初始缺陷的纤维束，在单调拉伸过程中，纤维束中的SiC材料逐渐损伤，以至最后失效的过程。

如图 2(b)所示，中间为半径R的模拟单丝纤维的圆柱体，外层模拟在纤维束制备过程中分布于纤维束内部的长宽高分别为a、a、H的SiC材料，并在中间预制一直径为d的缺陷，以模拟前文所述纤维束内部孔隙, 得到的数值如表 2所示。

2.2 有限元模型的建立

依据推算得到表 2的几何参数和表 3的材料常数，并通过特殊边界条件下的加载，依照前述能量法原理，可确定如表 4所示的整体纤维束单胞的材料常数。

由于T300单丝碳纤维在轴线方向上具有远大于SiC的拉伸强度，而有较小的刚度，因此在渐进损伤过程中，SiC材料会首先发生破坏。且SiC为高温陶瓷脆性材料，在此采用最大拉应力理论判定其失效。得到的z向(和单丝纤维轴线方向相同) 应力-应变结果如图 3所示。

从应力-应变曲线中可以看出，与树脂基复合材料破坏主要由C纤维体现不同，单向陶瓷基复合材料纵向破坏呈现两级破坏特点，当应变加载到0.1%左右时，SiC基体开始破坏，应力略微下降，单胞部分损伤，随后由于C纤维的承载应力继续上升。当应变加载到1.6%左右时，C纤维失去承载能力，平均应力开始下降。由此可获得纤维束的起始破坏强度为218 MPa，最终轴向拉伸强度为2 376 MPa。另外施加周期性边界条件^[9]可获得，纤维束的横向拉伸强度以及剪切强度，在表 5中列出。由表 3可以看出C纤维和SiC的压缩强度均较高，因此，本文计算中不考虑纤维束的压缩破坏。

3 C/SiC细观单胞有限元模型 3.1 单胞几何模型

图 4为C/SiC复合材料^[1]纤维束之间电镜扫描图，明显可见位于板材中间部分由于制作工艺而产生的孔隙，而正是这部分孔隙更容易成为材料破坏的起始位置，有必要对这部分单胞进行建模分析。由此，可按照图 5(a)选取代表性体积单元，并建立类似如图 5(b)所示的单胞模型，其横截面尺寸见图 6。

取复合材料中典型机织结构为一个单胞，如图 6所示为机织结构剖面图，假设机织结构中，纤维束截面形状为如图经纱所示的压扁形状，由两段相交圆弧组成，圆弧半径为R₁。纤维束形状由两段同心圆弧段组成，两段圆弧半径分别为R₁, R₁+h/2。模型中相关参数可由电镜照片中经纱纬纱图形测量得到。

3.2 单胞有限元模型

图 7(b)为二维C/SiC复合材料单胞有限元模型，浅灰色部分为基体，其余部分为四根交叠的纤维束。其中的纤维束单元为在第一步计算中得到了其初始模量及渐进损伤演化规律的平均化模型，且由于纤维束的分布特点，四根交叠纤维束由八个有不同方向材料参数的部分组成，同时单元坐标系的方向随纤维束的轴线旋转，如图 8所示。

为了顺利施加周期性边界条件^[9]，在有限元网格划分时，必须保证主从面上相应网格节点一一对应，采用网格扫掠的方法来实现，模型主要采用20节点六面体solid 187单元离散，纤维束与纤维束分界面及纤维束与基体分界面均采用共节点方式连接，以保证界面处变形的连续性。

3.3 材料性能退化方案及分析流程

针对C/SiC二维机织复合材料的特点建立了相应的损伤判据。对于纤维束单元，可建立^[4]基于损伤状态变量的受损材料的刚度矩阵，如图 7所示，其中，D_i代表不同损伤模式的损伤度。其定义为：

$ {D_i} = 1-\frac{{{E_i}}}{{{E_{io}}}}(i = 1, 2, \cdots, 6) $

(4)

式中，E_io代表i方向初始模量，E_i代表i方向实际模量。i=1, 2, 3分别代表纤维束轴向(L向)、横向(T向) 和法向(Z向)；i=4, 5, 6为TZ向、LZ向和LT向。达到损伤准则后，通过设定不同损伤模式的损伤度值，更新该损伤积分点的刚度矩阵，随着载荷增加，在以后增量步模拟中可以达到不同模式的损伤准则，从而继续更新刚度矩阵，实现渐进损伤模拟。对于C/SiC平纹机织复合材料纤维束，由于陶瓷和碳纤维的压缩强度均较大，且纵向弹性模量和横向弹性模量相差较大，因此仅考虑纵向拉剪和横向拉剪两种破坏形式。在此针对纤维束采用Hashin准则^{[4, 12]}，三维Hashin失效准则具体形式如下。

纤维束轴向拉剪损伤(M1, M2)：

$ {\left( {\frac{{{\sigma _1}}}{{{X_T}}}} \right)^2} + {\left( {\frac{{{\sigma _{12}}}}{{{S_{12}}}}} \right)^2} + {\left( {\frac{{{\sigma _{13}}}}{{{S_{13}}}}} \right)^2} \ge 1 $

(5)

纤维束横向拉剪损伤(M3)：

$ {\left( {\frac{{{\sigma _2} + {\sigma _3}}}{{{Y_T}}}} \right)^2} + {\left( {\frac{1}{{S_{23}^2}}} \right)^2}(\sigma _{23}^2-{\sigma _{22}}{\sigma _{33}}) + {\left( {\frac{{{\sigma _{12}}}}{{{S_{12}}}}} \right)^2} + {\left( {\frac{{{\sigma _{13}}}}{{{S_{13}}}}} \right)^2} \ge 1 $

(6)

损伤准则中X_T为纤维束轴向拉伸强度，Y_T为纤维束横向拉伸强度，S₁₂、S₂₃、S₁₂为剪切强度。

当纤维束达到损伤准则后，对各个方向设定不同损伤度，与到达模拟渐进损伤的方法不同，针对基体材料为SiC陶瓷，拉伸与压缩强度相差较大，因此采用莫尔强度准则判定其失效。

$ {\sigma _{{\rm{TM}}}} = {\sigma _1}- \frac{{[{\sigma _{\rm{t}}}]}}{{[{\sigma _{\rm{e}}}]}}{\sigma _3} \le [{\sigma _{\rm{t}}}] $

(7)

式中，[σ_t]为拉伸强度，[σ_c]为压缩强度。基于Hashin准则和莫尔准则的两种材料的刚度退化方案如表 6所示。其中，M1代表达到纤维束内部基体破坏强度的拉剪耦合破坏，M2代表达到纤维束整体破坏强度的拉剪耦合破坏。

通过ANSYS用户二次开发程序APDL，在材料本构关系中嵌入损伤准则和材料性能退化方案。通过不同损伤准则对组分材料单元进行判断，由此在损伤过程中退化相应的材料性能；利用Newton-Raphson方法计算各增量步单胞中单元应力应变，根据损伤准则更新单元本构矩阵，非线性求解达到收敛，再继续增加载荷步，如此循环完成分析。

$ \begin{array}{l} \left[{\begin{array}{*{20}{c}} {{\sigma _{11}}}\\ {{\sigma _{22}}}\\ {{\sigma _{33}}}\\ {{\tau _{23}}}\\ {{\sigma _{31}}}\\ {{\sigma _{12}}} \end{array}} \right] = \left\{ {\begin{array}{*{20}{c}} {(1 -{D_1}){C_{11}}}&{}&{}\\ {\left( {1 -\frac{{{D_1} + {D_2}}}{2}} \right){C_{12}}}&{(1 -{D_2}){C_{22}}}&{}\\ {\left( {1 - \frac{{{D_1} + {D_3}}}{2}} \right){C_{13}}}&{\left( {1 - \frac{{{D_2} + {D_3}}}{2}} \right){C_{23}}}&{(1 - {D_3}){C_{33}}}\\ {}&{}&{}\\ {}&{}&{}\\ {}&{}&{} \end{array}} \right.\\ \left. {\begin{array}{*{20}{c}} {}&{}&{}\\ {}&{}&{}\\ {}&{}&{}\\ {(1 - {D_4}){C_{44}}}&{}&{}\\ {}&{(1 - {D_5}){C_{55}}}&{}\\ {}&{}&{(1 - {D_6}){C_{66}}} \end{array}} \right\}\left\{ {\begin{array}{*{20}{c}} {{\varepsilon _{11}}}\\ {{\varepsilon _{22}}}\\ {{\varepsilon _{33}}}\\ {{\gamma _{23}}}\\ {{\gamma _{31}}}\\ {{\gamma _{12}}} \end{array}} \right\}\\ \end{array} $

(8)

4 数值分析与讨论 4.1 单胞几何尺寸的确定

由于单胞选取的是孔隙率较大的中间部分，且其孔洞由两部分组成，分别分布在纤维束内部的SiC和纤维束外部的基体中。因此，在此依据文献[13]给出的电镜扫描图对纤维束尺寸进行测量如图 9所示。其中，文献[13]中给出基体厚度为0.01 mm，其他尺寸依此比例确定，如表 8所示。其中g可由几何关系求得，满足关系：

$ {R_1} = \frac{{16{L^2} + 16Lg + 4{g^2}-{h^2}}}{{8h}} $

(9)

据此，可反解求得g。

同时，查得C/SiC平纹机织复合材料内部密度^[14]ρ_B为1.56 g/cm³，依据(10) 式及C纤维和SiC的密度ρ_f和ρ_m，以及假设C纤维的体积分数V_m和复合材料外部致密部分一致的前提下，可推知此时的总体孔隙率V_p：

$ {\rho _B} = {\rho _{\rm{f}}} \times {V_{\rm{f}}} + {\rho _{\rm{m}}} \times {V_{\rm{m}}} + {\rho _{\rm{p}}} \times {V_{\rm{p}}} $

(10)

且有：

$ {V_{\rm{m}}} = {V_{{\rm{m, i}}}} + {V_{{\rm{m, o}}}}, {V_{\rm{p}}} = {V_{{\rm{p, i}}}} + {V_{{\rm{p, o}}}} $

(11)

式中，V_{m, i}和V_{m, o}分别为纤维束内部和外部SiC的体积分数；V_{p, i}和V_{p, o}分别为纤维束内部和外部孔隙率。

同时依据电镜扫描图获得的复合材料单胞尺寸建立单胞模型，利用CATIA分析工具获取此时外部孔隙率V_{p, o}，基体和纤维束的体积分数V_{m, o}和V_f。由式(11)，可以最终推知，纤维束内部SiC和孔隙的体积分数V_{m, i}和V_{p, i}，如表 7所示。复合材料单胞几何参数如表 8所示。

依照前述能量法原理，计算得到C/SiC平纹机织复合材料的等效弹性常数以及文献[1]中查得的材料弹性常数见表 9。对比表 9中的数据可以看到，计算的面外剪切弹性模量结果略小于实验值。这是由于预制件在CVI过程中，在外部则由于沉积比较充分，填充比较完整。而计算模型只考虑复合材料中的这部分残留孔隙，并未考虑外部完整填充型单胞，导致了计算结果稍小于实验值。

4.2 应力-应变曲线

图 10分别给出了C/SiC机织复合材料试件在损伤准则下的应力-应变预测曲线和试验曲线。从图 10可以看出，二维机织复合材料试件试验曲线和模拟曲线在初始阶段均呈现线性关系。模拟曲线逐步偏离线性是由于局部损伤的发生对宏观应力-应变曲线产生影响，在损伤累积到一定程度以后，曲线开始平缓下降，损伤的逐渐累积将导致材料的最终破坏。

二维机织C/SiC复合材料的力学行为主要由纤维束决定。根据最大拉应力准则判定SiC材料损伤，并引入了纤维束的渐进损伤规律。由于纤维束轴向拉伸强度较大，数值模拟中纤维束损伤的单元不多。在ε=0.06%之前，单胞未发生损伤，应力-应变曲线呈现线性；在ε=0.06%之后，基体单元开始发生损伤，随着损伤单元的增多，应力-应变曲线开始呈现非线性；在ε=0.3%时，继续加载，纤维束各种损伤相对稳定扩展，材料应力-应变曲线逐渐下降，并保持一定的承载能力，材料具有一定的延展性。如图 11(a)和11(b)分别为C/SiC平纹机织复合材料纱线和基体的损伤扩展图。由ε=0.06%的纱线损伤图可以看出由于纤维束横向强度较低，因此和纵向纤维交会部分会先产生损伤。这与实验中观察到的结果一致。由图 11中ε=0.6%和ε=1.2%可以看出, 沿加载方向纤维束受损严重，且中部薄弱处最先破坏。由基体损伤场扩展图可知，基体损伤由中部缺陷开始，并逐渐扩展至整个基体，且在应变较小时损伤较大。

4.3 渐进损伤过程中的材料参数变化

C/SiC单胞单轴拉伸受损后，其本构方程可写为：

$ \mathit{\boldsymbol{\bar \sigma }} = \mathit{\boldsymbol{C}}_D^H\mathit{\boldsymbol{\bar \varepsilon }} $

(12)

式中，σ, ε为单胞整体的平均应力和平均应变，C_D^H为随应力大小增加而变化的损伤本构矩阵；当处于卸载过程，材料未达继续损伤条件时，此本构矩阵保持不变。在某一载荷水平下，依第1节中方法施加边界条件后，可以得到发生损伤后单胞的平均应力矩阵σ和平均应变矩阵ε。由此可计算得到此载荷水平下损伤本构矩阵C_D^H中的系数，并由此推算出单轴作用下发生损伤后的材料参数变化如表 10所示，据此可以获得材料的损伤本构矩阵。

5 结论

针对平纹机织C/SiC复合材料，利用电镜扫描图，分别建立了内部和外部缺陷模型。同时采取"两步走"的策略，将前一个尺度的计算结果作为后一个尺度性能计算的基础，引入破坏准则，验证了本文所提模型以及计算方法的有效性。

(1) 初始阶段C/SiC二维机织复合材料单胞的材料常数与文献中多个实验给出的数据区间一致，证明了所提出的多尺度复合材料单胞的合理性。

(2) 应力-应变曲线和实验相差较小以及损伤场的变化趋势与电镜观察一致，证明了基于莫尔强度理论和Hashin准则的分别针对于基体和纤维束的破坏准则的合理性。

(3) 可以进一步研究含界面多层陶瓷基复合材料板件中，因制备工艺和技术造成的孔隙大小及其材料组分对复合材料等效性能的影响，进而展开材料性能的优化设计工作。