0 引言 连续纤维增韧陶瓷基复合材料能从根本上克服传统结构陶瓷的脆性缺陷。C/SiC复合材料以其比强度、比模量高，耐高温等优良的性能广泛应用于航空发动机热端部件及高超声速飞行器热防护系统等方面，是重要的国防装备材料[1-2]。 CFCC-SiC本征结构异常复杂，制备工艺特点和基体与纤维线胀系数失配产生的热残余应力造成数量众多、尺度不一的细观裂纹和孔隙，这些原始缺陷在较低应力时就导致基体开裂、纤维脱粘等现象，细观损伤的不断积累导致材料失效。通过建立单胞物理模型将复杂的细观结构解耦，能有效分析复合材料的损伤、断裂、强度问题及各相间的应力传递机制、破坏机制。 目前，复合材料热残余应力的研究主要集中在粗纤维、金属基和树脂基复合材料，对C/SiC复合材料研究报道甚少，缺乏系统深入的研究。刘徵等[3]利用直径30 μm的微小X射线束衍射法研究了SiC/Al模型复合材料Al基体残余应力分布情况；谢薇等[4]研究了处理工艺对C/Mg复合材料热残余应力的影响；Ward等[5]用拉曼光谱法测量了两种SiC纤维的轴向热残余应力。C/SiC复合材料中T300直径微小，SiC基体耐腐蚀及脆硬性等因素对实验研究带来很大困难，同时难以测得具体应力分布且实验结果具有离散性，数值模拟方法为C/SiC复合材料热残余应力研究提供了有力的工具。 本文通过CVI工艺[6-7]制备单向C/SiC复合材料，采用细观力学分析方法对纤维、基体中热残余应力分布规律进行预测，使用ABAQUS6.11有限元软件建立“纤维-界面-基体”单胞物理分析模型，对制备过程进行仿真，研究界面厚度、纤维体积分数、制备温度等因素对复合材料横截面内热残余应力分布的影响规律。研究所揭示规律有助于指导材料的设计、制备及微纳力学性能试验，为探明CFCC加工的材料去除机理奠定理论基础。 1 单向C/SiC复合材料 1.1 材料制备 CVI具体制备过程如下：根据体积分数进行质量换算称重所需T300碳纤维并编织单向纤维预制体，采用带气孔的专用石墨夹具固定预制体 (图 1)。以丙烯为反应气体经数小时在纤维表面沉积PyC涂层，再以过量H2为载气通过鼓泡方式将CH3SiCl3引入反应室 (摩尔比10:1)，1 000℃条件下经30 h CVI致密化后获得成型材料。 1.2 材料表征 采用美国Quanta-FEG-250扫描电子显微镜和日本VHX-1000光学显微镜观察材料显微形貌，成型材料密度2.1 g/cm3、孔隙率10%，体积分数40%、界面层厚度100~200 nm。图 2为材料样品及显微形貌，纤维、界面、基体、孔隙等结构单元使得复合材料呈尺度跨度的多层次性、细观结构的非连续性、空间分布的非均匀性特征，如原始孔隙的存在及基体分布不均匀，这是导致材料热残余应力复杂分布的重要原因。 2 细观力学建模及分析 考虑三相同心圆模型如图 3所示。 极坐标下由弹性力学理论并考虑对称性有： 平衡方程 $\frac{\text{d}{{\sigma }_{r}}}{\text{d}r}+\frac{{{\sigma }_{r}}-{{\sigma }_{\theta }}}{r}=0$ (1) 几何方程 ${{\varepsilon }_{r}}=\frac{\text{d}u}{\text{d}r},{{\varepsilon }_{\theta }}=\frac{u}{r}\text{ }$ (2) 物理方程 $\left\{ \begin{align} & {{\varepsilon }_{r}}=\frac{1}{E}\text{ }({{\sigma }_{r}}-\nu {{\sigma }_{\theta }}) \\ & {{\varepsilon }_{\theta }}=\frac{1}{E}\text{ }({{\sigma }_{\theta }}-\nu {{\sigma }_{r}}) \\ \end{align} \right.$ (3) 式中,ν、E为材料泊松比、弹性模量，联立各式并考虑对称性得通解表达式： $\left\{ \begin{align} & {{\sigma }_{ri}}=\frac{{{A}_{i}}}{{{r}^{2}}}+2{{C}_{i}} \\ & {{\sigma }_{\theta i}}=-\frac{{{A}_{i}}}{{{r}^{2}}}+2{{C}_{i}} \\ \end{align} \right.$ (4) ${{u}_{i}}=\frac{1}{{{E}_{i}}}[-(1+{{\nu }_{i}})\text{ }\frac{{{A}_{i}}}{r}+2(1-{{\nu }_{i}})\text{ }{{C}_{i}}r]$ (5) 式中，下标i=1、2、3与基体、界面、纤维对应，Ai、Ci作为待定系数由如下边界条件确定： $\begin{matrix} {{A}_{3}}=0 & \left( r\to 0 \right) \\ \end{matrix}$ (6) $\begin{matrix} {{\sigma }_{r1}}=0 & (r={{r}_{1}}) \\ \end{matrix}$ (7) $\begin{matrix} {{\sigma }_{r1}}={{\sigma }_{r2}} & ~{{u}_{1}}+{{\alpha }_{1}}\Delta Tr={{u}_{2}}+{{\alpha }_{2}}\Delta Tr & (r={{r}_{2}}) \\ \end{matrix}$ (8) $\begin{matrix} {{\sigma }_{r2}}={{\sigma }_{r3}} & {{u}_{2}}+{{\alpha }_{2}}\Delta Tr={{u}_{3}}+{{\alpha }_{3}}\Delta Tr & (r={{r}_{3}}) \\ \end{matrix}$ (9) 联立各式得： $\left[ {\matrix{ {} \hfill & {1/{r_1}^2} \hfill & 0 \hfill & 2 \hfill & 0 \hfill & 0 \hfill \cr {} \hfill & {1/{r_2}^2} \hfill & { - 1/{r_2}^2} \hfill & 2 \hfill & { - 2} \hfill & 0 \hfill \cr {} \hfill & 0 \hfill & {1/{r_3}^2} \hfill & 0 \hfill & 2 \hfill & { - 2} \hfill \cr {} \hfill & {{K_{41}}} \hfill & {{K_{42}}} \hfill & {{K_{43}}} \hfill & {{K_{44}}} \hfill & 0 \hfill \cr {} \hfill & 0 \hfill & {{K_{52}}} \hfill & 0 \hfill & {{K_{54}}} \hfill & {{K_{55}}} \hfill \cr } } \right]\left[ \matrix{ {A_1} \hfill \cr {A_2} \hfill \cr {C_1} \hfill \cr {C_2} \hfill \cr {C_3} \hfill \cr} \right] = \left[ \matrix{ 0 \hfill \cr 0 \hfill \cr 0 \hfill \cr ({\alpha _2} - {\alpha _1})\Delta T{r_2} \hfill \cr ({\alpha _3} - {\alpha _2})\Delta T{r_3} \hfill \cr} \right]$ (10) 式中,ΔT为制备温差，α为线胀系数，K为与ν、E有关系数。修改边界条件，在两相同心圆模型中i=1、2代表基体、纤维，则有： $\left\{ \begin{align} & {{\sigma }_{r2}}={{\sigma }_{\theta 2}}=\frac{\left( 1-m \right)({{\alpha }_{2}}-{{\alpha }_{1}})\Delta T{{E}_{1}}{{E}_{2}}}{{{E}_{2}}(1+m+m{{\nu }_{1}}-{{\nu }_{1}})-{{E}_{1}}(1-m-{{\nu }_{2}}-m{{\nu }_{2}})} \\ & {{\sigma }_{r1}}=\frac{({{\alpha }_{2}}-{{\alpha }_{1}})\Delta T{{E}_{1}}{{E}_{2}}}{{{E}_{2}}(1+m+m{{\nu }_{1}}-{{\nu }_{1}})-{{E}_{1}}(1-m-{{\nu }_{2}}-m{{\nu }_{2}})}\cdot \left( 1-\frac{r_{1}^{2}}{{{r}^{2}}} \right) \\ & {{\sigma }_{\theta 1}}=\frac{({{\alpha }_{2}}-{{\alpha }_{1}})\Delta T{{E}_{1}}{{E}_{2}}}{{{E}_{2}}(1+m+m{{\nu }_{1}}-{{\nu }_{1}})-{{E}_{1}}(1-m-{{\nu }_{2}}-m{{\nu }_{2}})}\cdot \left( 1+\frac{r_{1}^{2}}{{{r}^{2}}} \right) \\ \end{align} \right.$ (11) 该式表明纤维径向和周向残余应力等值同号且为定常应力场，而基体内径向和周向残余应力随与界面距离增加而快速减小。平面应变问题只需将式中ν、E、α替换为${\nu \over {1 - \nu }}、{E \over {1 - {\nu ^2}}}$、(1+ν)α即可，式中$m = {{r_1^2} \over {{r^2}}}$。 3 有限元模型建立 3.1 界面内聚力模型数值模拟方法 内聚力模型 (Cohesive Zone Model) 较传统建模方法能更好的反映界面热残余应力传递规律。如图 4所示，本构模型采用Mi等[8]提出的双线性内聚力模型，材料的强度、裂纹的形成和扩展都由本构关系决定。δinit、δfail分别为材料初始损伤及完全失效等效位移，本构关系见式 (12)。 $\sigma = \left[ \matrix{ {\sigma _n} \hfill \cr {\sigma _s} \hfill \cr {\sigma _t} \hfill \cr} \right] = \left[ \matrix{ \left( {1 - D} \right){\rm{ }}{K_{nn}} & 0 & 0 \hfill \cr 0 & \left( {1 - D} \right){\rm{ }}{K_{ss}} & 0 \hfill \cr 0 & 0 & \left( {1 - D} \right){\rm{ }}{K_{tt}} \hfill \cr} \right]\left[ \matrix{ {\varepsilon _n} \hfill \cr {\varepsilon _s} \hfill \cr {\varepsilon _t} \hfill \cr} \right]$ (12) 选取基于能量、线性软化、独立模式的损伤演化规律。课题组已有实验[9]和数值模拟结果对比修正得纯Ⅰ型、纯Ⅱ型及纯Ⅲ型破坏的最大名义应力分别为10、30、30 MPa，断裂能取6 N/m。 3.2 数值模拟建模及实验验证 图 5为C/SiC复合材料抛光横截面显微形貌图。依此建立图 6所示二维六方单胞物理模型。模型包括T300、PyC界面、SiC基体，T300直径7 μm，基体尺寸由体积分数算出。根据对称性取1/4模型，界面附近纤维和基体网格加密以提高计算精度，远离界面区域为不均匀稀疏网格以降低计算成本。纤维和基体单元类型为四节点平面应变单元 (CPE4R)，界面为内聚力单元 (COH2D4)，通过“tie”约束保证接触的位移、应力协调。 数值模拟制备过程中界面粘结完好、没有材料缺陷，且最高温度时无热应力，不考虑纤维与纤维之间应力场的相互影响，对分析模型左边及底边施加对称约束边界条件。定义各相材料性能见表 1。 鉴于与刘徵等[3]研究对象结构相似且残余应力产生原因相同，将此建模方法用于SiC/Al复合材料热残余应力数值模拟，图 7所示基体径向应力场理论计算、数值模拟、试验结果吻合良好,表明此数值建模方法合理可行。 4 数值模拟结果与讨论 4.1 复合材料热残余应力数值模拟研究方法 研究主要考察C/SiC复合材料纤维和基体中径向、周向热残余应力的分布规律。结合CVI工艺及工程应用实际情况，采用控制单因素变量的方法研究界面厚度、纤维体积分数、制备温度对热残余应力分布的影响规律，方案见表 2。规定模型中纤维水平径向为0°方向，逆时针角度为正。 4.2 复合材料横截面内热残余应力分布规律 纤维、基体中残余应力分布规律与公式 (11) 基本一致，但数值模拟能获得应力场具体分布规律。 图 8为径向热残余应力分布图，图 9为基体残余应力云图 (H=200 nm,Vf=40%,T=1 000℃)。 数值模拟结果表明，纤维内近似定常力场，但表面应力有微小差别，基体中σr变化梯度与角度有关，且σθ出现应力集中，最大拉应力位置及数值见图 9。 基体中径向压应力对纤维形成箍筋作用，有利于延迟界面横向脱粘；过大周向拉应力使制备工艺形成的原始孔隙成为裂纹源，未达到基体断裂强度时就出现开裂。同时，细观结构中脆性陶瓷基体的应力集中现象将一直保持在复合材料整个承载过程中直至达到强度极限[12]，这是导致复合材料横向强度低于基体强度的重要原因。 图 10为纤维表面径向、周向热残余应力分布 (H=200 nm,Vf=40%,T=1 000℃)。可知，纤维表面σr、σθ是周期为π/3的压应力函数，纤维表面最大σr位于θ=30°+n·60°方向，最小σr位于θ=n·60°方向；而σθ的分布规律与σr相反，最大σθ位于θ=n·60°方向，最小σθ位于θ=30°+n·60°方向。 制备降温过程中多基体区收缩程度大于少基体区，对周围区域形成牵拉作用,导致基体收缩方向由少基体区指向多基体区，如图 11箭头所示。这是导致图 8~图 10中纤维、基体热应力场复杂分布规律的主要原因。 4.3 界面层厚度对热残余应力分布规律影响 图 12为界面厚度对纤维轴心热残余应力值影响趋势 (Vf=40%,T=1 000℃)。 纤维中σr、σθ等值同号，界面越厚压应力值越大。由图 13及图 9可知随界面厚度增加基体中σr分布状态变化不大，靠近界面基体σr压应力显著增大导致基体中应力梯度增大；基体中σθ为拉应力，应力集中区域相同，随H增加应力集中现象更严重，最大值由148.54 MPa升至530.01 MPa，这是因为较厚的界面层不易变形，热应力的释放更加困难。因此制备过程中应控制沉积时间使界面层厚度适当。 4.4 纤维体积分数对热残余应力分布规律影响 图 14为体积分数对纤维轴心热残余应力值影响趋势 (H=100 nm，T=1 000℃)。可知，Vf与残余压应力值成负相关。 由图 15及图 13(a)(b) 可知，Vf对基体σr影响较小，但对σθ分布规律影响明显。基体中拉应力梯度增大，应力集中现象严重，最大σθ由111.49 MPa增加至315.05 MPa，且位置转移导致基体承载危险点变化。Vf=30%时少基体区没有应力集中，多基体区及与界面粘结处已出现应力集中，Vf＞40%后，应力集中现象更加严重，使役过程中应力集中区成为细观结构危险区。故复合材料Vf不宜过高，但过低又难以实现增韧和补强作用，工程实际常用的体积分数为35%~45%。 4.5 制备温度对热残余应力分布规律影响 图 16为制备温度对纤维轴心热残余应力值影响趋势 (H=100 nm，Vf=40%)。可知，纤维σr、σθ与制备温度成线性关系，与公式 (7) 理论计算一致。由图 17及图 11(a)(b) 可知，不同制备温度条件下基体中σr、σθ分布规律基本相同，随制备温度的升高各相残余应力都有增大，但σθ变化幅度较大。 5 结论 (1) 细观力学模型计算表明纤维内径向和环向残余应力等值同号且为定常应力场，基体内径向和环向残余应力随与界面距离增加而快速减小。 (2) 成型材料具有多层次、非连续、非均匀的结构特征，数值模拟六方模型中纤维表面应力是周期为π/3的周期性函数，多基体区界面附近和少基体区会出现应力集中。 (3) 界面厚度、纤维体积分数、制备温度对复合材料热残余应力分布规律有影响。实际制备工艺中界面厚度100 nm左右为宜，体积分数在35%~45%为宜，制备温度1 000℃为宜。 (4) 细观力学模型分析、数值模拟分析与试验结果吻合较好。基于数值模拟的方法研究复合材料热残余应力是可行的，通过对制备过程的仿真能够弥补力学分析和实验手段的不足，有效揭示复合材料热残余应力的分布规律。