基于多线性本构与损伤耦合的叠层陶瓷基复合材料数值预测方法.pdf

航空科学技术Aeronautical Science&TechnologyJun.25 2023 Vol.34 No.06 54-65基于多线性本构与损伤耦合的叠层陶瓷基复合材料数值预测方法刘斌1，曹立阳2，王波1，杨腾飞1，刘永胜2，司源31.西北工业大学，陕西 西安 7100722.超高温结构复合材料重点实验室，陕西 西安 7100723.航空工业济南特种结构研究所 高性能电磁窗航空科技重点实验室，山东 济南 250023摘 要：本文提出一种用于预测复杂应力状态下叠层C/SiC损伤破坏的有限元渐进损伤方法（FE-PDM），并提出与试验载荷位移曲线、损伤等进行对比、迭代并最终确定模型参数的反演方法。FE-PDM方法包括基于应变控制的3D失效准则，正交各向异性多线性本构关系，基于刚度矩阵的损伤因子耦合方法，以及预测分层的内聚力方法。通过2D叠层C/SiC的面内拉伸、面内剪切、三点弯曲试验分别验证了FE-PDM方法，并分别采用扫描电子显微镜及X射线技术对试样断口形貌和内部损伤进行了分析。结果表明，通过较少的控制参数，FE-PDM方法可精确地预测2D叠层C/SiC各阶段的应力应变曲线拐点、损伤起始与演化过程及载荷位移响应曲线等，并与试验结果吻合良好。关键词：陶瓷基复合材料；FE-PDM；本构关系；损伤演化中图分类号：V214.8 文献标识码：A DOI：10.19452/j.issn1007-5453.2023.06.008 陶瓷基复合材料（CMCs）具有强度高、耐高温、密度低、硬度大等优点，广泛用于航空航天领域，如高超声速飞行器、航天飞机、航空发动机等1-2。二维（2D）织物叠层及其演化型是CMCs常见的细观结构类型3。但CMCs在微/细观尺度上存在的孔洞缺陷4、复杂的裂纹扩展规律、各向异性以及伪塑性行为给其数值模拟带来了诸多困难，特别是如何快速准确地预测复杂受力状态下陶瓷基复合材料的损伤情况。现阶段，陶瓷基复合材料微、细观结构的多尺度建模技术与计算方法得到了快速的发展。C.Chateau 等5建立了SiC/SiC复合材料微观模型，主要研究了化学气相渗透（CVI）工艺下材料的孔隙率对材料弹性性能的影响。L.Borkowski等6提出了一种微观力学与热弹性渐进损伤相耦合的多尺度方法，该方法可以精确预测热-力耦合条件下平纹编织C/SiC复合材料的弹性模量和损伤行为，并且在模型中引入了制造过程中由热失配引起的初始损伤。Borkowski认为传统的分析方法不能体现CMCs力学行为的多尺度现象，但是也没有考虑到该多尺度方法用于分析CMCs宏观力学问题需要付出的计算代价与时间成本，并且引入的初始损伤的位置效应还需要在研究中进一步讨论。Dong Hongnian等7指出目前很少有学者对编织CMCs的疲劳寿命进行多尺度分析预测，研究多集中在预测CMCs的模量、应力/应变响应方面，因此提出了对CMCs疲劳寿命的多尺度预测方法，但是其微观力学模型是理论假设推导而非有限元模型，其适用性还需进一步验证。Liu Bin等8从CMCs（如SiCf/SiC）的微观尺度出发，建立了带周期性边界条件的代表性体积单元模型（3D-RVE），系统地分析了SiCf/SiC复合材料的刚度和热残余应力（TRS），并在制造过程中加入高温条件，对环向、径向和轴向的纤维、界面和基体的热残余应力进行了详细的计算和分析。收稿日期：2023-03-01；退修日期：2023-04-15；录用日期：2023-05-16基金项目：国家自然科学基金（51902256）；航空科学基金（2020Z057053002）；工业装备结构分析国家重点实验室开放基金（GZ21115）引用格式：Liu Bin，Cao Liyang，Wang Bo，et al.Numerical prediction method of laminated ceramic matrix composite based on the multilinear constitutive and the coupling of damageJ.Aeronautical Science&Technology，2023，34（06）：54-65.刘斌，曹立阳，王波，等.基于多线性本构与损伤耦合的叠层陶瓷基复合材料数值预测方法 J.航空科学技术，2023，34（06）：54-65.刘斌 等：基于多线性本构与损伤耦合的叠层陶瓷基复合材料数值预测方法随着声发射技术（AET）、数字图像技术（DIC）及电子计算机断层扫描（CT）技术的快速发展，CMCs的力学试验研究得到了重大突破。微、纳米CT无损检测技术可以真实反映材料的初始缺陷，并将其引入多尺度模型中，从而建立含缺陷的多尺度模型。L.Glbart等9将数字图像技术引入编织SiC/SiC复合材料的细观尺度研究中，提出了一种新的计算方法用于预测裂纹张开位移、裂纹密度、轴向应力及应变。Liu Yu等10采用CT原位测试技术对平纹SiC/SiC复合材料进行了拉伸损伤演化研究，得到了材料试验过程中真实的微观结构和损伤演化图像，并使用深度学习方法进行损伤识别。Chen Yanfei等11基于X射线微断层扫描数字图像（IB-FEM）技术，捕获了CVI工艺下C/SiC复合材料的细观结构图像。2D IB-FEM结果表明，完整模型与含缺陷模型之间的应力应变曲线偏差随着应变的增大而增大，3D IB-FEM计算得到的损伤演化、等效应力/应变曲线均与原位X射线试验结果吻合。试验发现，损伤一开始从缺陷处缓慢扩展，然后到达某一应力水平后迅速扩展，这说明材料内部的原始缺陷是损伤演化的主要因素和起始位置。因此，从CMCs宏观本构关系的角度出发，把握其本构关系、损伤模式、刚度退化是本文所采取的仿真策略。目前对CMCs宏观本构关系的研究主要集中在单应力状态上。Jiao Guiqiong等12针对2D C/SiC复合材料进行面内剪切及偏轴拉伸试验研究，并结合材料宏观应力应变本构关系，提出了面内经、纬向损伤系数及剪切损伤系数对材料刚度折减的影响。有限元方法适用于复杂应力状态的预测，模型中需考虑不同的材料失效准则，如Puck13准则、Hashin14准则，利用刚度退化的方法进行应力等效折减。这些准则从20世纪发展至今已经十分成熟，而且预测精度和计算效率都很高，但大都是针对连续纤维增强树脂基复合材料以及单向带铺层结构提出和验证的，涉及编织CMCs结构相对较少。Zhang Yi等15重点研究了CVI工艺下 z-pin加强2D C/SiC复合材料的剪切力学行为，提出采用Chang-Lessard准则对2D C/SiC复合材料进行失效判断，为了拟合拉伸和剪切应力/应变的损伤系数，采用具有双线性本构关系的内聚力单元模型对z-pin剪切面进行描述。Gao Xiguang16应用渐进式失效分析（PFA）与DIC技术对带圆孔C/SiC板进行拉伸试验，研究其应变分布与损伤演化，提出了连续折减策略（NCRCDS）的非线性本构关系，应力应变曲线通过试验进行了指数型拟合，拟合效果良好，但是拟合公式中需要确定的参数过多，如面内拉伸需要确定 5 个参数，面内剪切需要确定三个参数。CMCs的应力应变数据通常具有较大的分散性，为了反映真实的力学行为，本文利用CMCs简化的非线性本构关系，同时考虑到正应力与剪应力的耦合情况，提出了基于应变控制失效的6种损伤准则。基于这6种失效模式和CMCs的多线性本构关系，考虑到不同应力之间的耦合效应，采用刚度退化的方法，利用Abaqus/Explicit中的材料用户子程序VUMAT进行编写，将其应用于复杂应力状态下CMCs的损伤预测。1 材料2D CMCs的结构由多层织物堆叠而成，堆叠方式可分为同相态（IPM）、异相态（OPM）和混合随机态（RPM）。本文提出的FE-PDM方法包含CMCs宏观本构关系及其渐进损伤演化，适用于这三种状态的 CMCs 结构。2D CMCs的层间不是理想的平面，而是凹凸不平的曲面，同时纤维束间存在细观孔洞以及纤维之间存在微观孔隙，如图1所示（以C/SiC为例）。在单层织物中，定义经纱方向为1方向，纬纱方向为2方向，垂直于单层织物的面外方向为3方向。在本文的第3节中，建立了CMCs的宏观本构关系，并基于这三个方向建立了CMCs失效模式和损伤演化行为。本文算例采用的是C/SiC复合材料，该材料由T300碳纤维正交编织成碳布预制体，采用化学渗透工艺沉积热解碳界面层（PyC），基于化学气相渗透（CVI）技术沉积SiC基体形成。图12D叠层陶瓷基复合材料（以C/SiC为例）:细观结构、微观结构、缺陷及方向定义 Fig.12D ceramic matrix composites(taking C/SiC as an example)explanations:meso-structure,micro-structure,defects and direction definition 55航空科学技术Jun.25 2023 Vol.34 No.062 FE-PDM方法2.1 损伤模式和破坏准则3D Hashin准则适用于复合材料层压板损伤17-18。Y.J.Lee和Huang19提出，刚度系数在局部点退化后，应力可能会产生较大波动，无法连续增大，从而在下一增量步时导致计算失败，但应变比应力更加连续、平滑，用其作为失效判据更加合适。本文开发的 FE-PDM 方法用来精确模拟CMCs宏观结构受复杂应力下的损伤破坏，如图2所示。本研究采用基于应变控制的三维失效准则和三维Hashin准则的修正形式，失效模式包括经向拉伸和压缩、纬向拉伸和压缩以及面外拉伸和压缩共6种。这6种损伤模式可以表示为式（1）式（6）的基于应变控制的失效准则，用于描述层内单元的起始破坏。模式1：经向拉伸 110f1=(11XT)2+(12S12)2+(13S13)2（1）模式2：经向压缩 110f1=(11XC)2+(12S12)2+(13S13)2（2）模式3：纬向拉伸 220f2=(22YT)2+(12S12)2+(23S23)2（3）模式4：纬向压缩 220f2=(22YC)2+(12S12)2+(23S23)2（4）模式5：面外拉伸 330f3=(33ZT)2+(13S13)2+(23S23)2（5）模式6：面外压缩 330f3=33ZC（6）式 中，11,22,33,12,23,13表 示 应 变 的 6 个 分 量，XT,XC,YT,YC,ZT,ZC,S12,S23,S13表示损伤起始时的应变分量，f1，f2，f3表示CMCs材料是否进入损伤状态，如果f1，f2，f3等于1，则将相应的应变分量记为损伤起始应变，FE-PDM方法将降低CMCs的原位刚度。2.2 层内损伤演化通过降低M.Kachanov等20提出的刚度系数来考虑材料损伤。如式（7）所示，采用A.Matzenmiller21和I.Lapczyk等22提出的模型来计算正交各向异性材料刚度矩阵系数的退化。有效应力 与标称应力之间的关系为23-25=M=Cd（7）式中，M、Cd和分别表示损伤因子、损伤刚度矩阵和应变矩阵。=112233233112T（8）=112233233112T（9）M=diag 11-d1 11-d211-d3 11-d4 11-d5 11-d6（10）=S=Sd=(MS)（11）式中，S、Sd分别为原始弹性矩阵和损伤弹性矩阵。Sd=MS（12）S=1E1-21E2-31E312E11E2-32E313E1-23E2-1E21G231G311G12（13）Cd可以表示为Cd=Sd-1=MS-1（14）bzyyzyxzxxzxyy?int?f?+/?y?z?1,2?fczyxint23ftOtcmcft+?fc?fc?int?3intftfcO+int23ftOfc?1,2tcmcft O图22D叠层CMCs层内损伤演化策略Fig.2Damage evolution strategy of intralaminar 2D CMCs56刘斌 等：基于多线性本构与损伤耦合的叠层陶瓷基复合材料数值预测方法C d=C11 C12 C13C21 C22 C23C31 C32 C33 C44 C55 C66（15）=1-(1-d1)(1-d2)v12v21-(1-d1)(1-d3)v13v31-(1-d2)(1-d3)v23v32-2(1-d1)(1-d2)(1-d13)v12v23v31其中，Cij=Cji。C11=(1-d1)E1-(1-d1)(1-d2)(1-d3)E12332（16）C22=(1-d2)E2-(1-d1)(1-d2)(1-d3)E21331（17）C33=(1-d3)E3-(1-d1)(1-d2)(1-d3)E31221（18）C23=(1-d2)(1-d3)E323+(1-d1)(1-d2)(1-d3)E31321（19）C13=(1-d1)(1-d3)E313+(1-d1)(1-d2)(1-d3)E31223（20）C12=(1-d1)(1-d2)E212+(1-d1)(1-d2)(1-d3)E21332（21）C44=(1-d4)G23（22）C55=(1-d5)G31（23）C66=(1-d6)G12（24）从式（7）、式（15）式（21）可知，1、2、3除与自身应变相关外，还与另两个方向的应变相关，将其定义为耦合关系。如1与C11、C12和C13相关，分别受d1、d2和d3影响。采用B.G.Falzon等26提出的方法计算损伤系数如下di()=fifi-oi(1-oii)（25）（1）d1和d2的计算方法图3（a）给出了CMCs经向和纬向上拉伸与压缩的本构关系。对于拉伸过程，根据CMCs正交各向异性的本构关系，应力应变曲线共有4个阶段：第一个拐点表示产生了隧道裂纹(tc,int)，第二个拐点表示基体开裂(mc,2)，第三个拐点表示纤维断裂(ft,3)。其中，在第三阶段，基体裂纹不断增大，直至达到饱和状态，基体裂纹发生偏转，出现纤维基体界面脱黏。在第四阶段，当应力达到ft后存在纤维拉出现象。对于压缩过程，当应力和应变达到阈值(fc,fc)时，CMCs发生脆性断裂。对于拉伸和压缩试验，通过利用式（25），d1、d2可分别表示为式（26）、式（27）。1，220时d1 or d2=0 int1-int-int(sa-mc)(-int)(2-int)mc int 2 2-2-2(ft-sa)(-2)(3-2)sa-int2-int(sa-mc)2mc 2 3ftft-3(1-3)+2-23-2(ft-sa)3sa-int2-int(sa-mc)2mc 3220时（26）d1，d2=0 fc1 fc（27）图3（b）显示了当材料同时承受拉伸和切应力时，裂纹的萌生、扩展和损伤演化规律。在相对较低的应力水平下，垂直于拉伸应力方向的纤维束中出现隧道裂纹；随着应力的不断提高，基体裂纹沿着拉伸应力方向在纤维束中相继出现。在基体裂纹发生偏转和纤维基体界面脱黏后，纤维束中的纤维最终承受拉伸应力和切应力直至断裂。但对于1方向或2方向上的压缩与剪切组合情况，不会出现隧道裂纹和基体开裂的现象，最终损伤模式为压缩剪切断裂，如图3（c）所示。根据损伤模式的6种状态，CMCs的失效标准可设计如下：当CMCs进入失效状态时，三个法向应力和三个切应力将通过刚度退化的方式减小，这将在3.2节中进行详细阐述。（2）d3的计算方法图4（a）给出了CMCs面外方向上的拉压本构关系。对于CMCs的拉伸过程，其应力应变曲线为双线性阶段：第一个拐点(in,int)表示隧道裂纹、基体裂纹和分层的产生。通过式（25），d3可表示为式（28）、式（29），分别应用于拉伸和压缩过程。此外，压缩过程中系数d3的计算方法还适用于低能量冲击以及三点弯曲等情况。对于带剪切的面外拉压力，图4（b）、图4（c）基本可以说明隧道裂纹和分层的破坏。面外压缩试验没有用示意图的形式表示，而是以三点弯曲为示例在第4节中进行试验与仿真模拟验证。330时d3=0 fcftft-int()1-int int （28）330时d3=fc1-fc int（29）（3）d4、d5、d6的计算方法图5给出了CMCs面内、面外剪切的本构关系，应力应变曲线为双线性阶段：第一个拐点表示CMCs的屈服点57航空科学技术Jun.25 2023 Vol.34 No.06(y,int)，第二个拐点表示最终剪切破坏点(b,f)。通过利用式（25），d4、d5和d6表示为式（30）0 or 0d4，d5，d6=0 int1-int-int(b-y)(-int)(f-int)y int(30)2.3 分层内聚力模型基于Dugdale-Barenblatt模型的内聚力模型（CZM），假设在裂纹尖端周围存在一个小范围的塑性区。CZM单元用于表示塑性区，不需要预制初始裂纹。图6显示了弹塑性裂纹、CZM单元几何结构和应力状态27-28。如图6（a）所示，CZM单元承受法向拉力和滑动剪切力。CZM单元有以下三个方向（n、s、t）的应变n=nTo,s=sTo,t=tTo（31）m=2s+21+n2=2shear+n2（32）O/ybintf图5面外剪切本构方程示意图及其损伤模式Fig.5The constitutive relation of out-of-plane shearOaCMCs?int23ft?fc?fc?mc?fc?1,2fc图3经纬向拉伸与压缩本构方程以及损伤模式1和3Fig.3The constitutive relation of warp or weft tension andcompression and damage mode 1&3 aCMCs?Ointftfcint3fc图4面外拉压本构方程示意图及其损伤模式Fig.4 The constitutive relation of out-of-plane tension andcompression 58刘斌 等：基于多线性本构与损伤耦合的叠层陶瓷基复合材料数值预测方法式中，m为等效位移，To表示CZM单元厚度。CZM单元的本构关系如式（33）所示t=tntstt=KnnKnsKntKnsKssKstKntKstKtt nst=K（33）CZM单元损伤过程包括损伤起始和扩展两个阶段，式（34）式（38）描述了分离过程，D是退化参数。tn=(1-D)tn tn 0 tn tn 0 0shear n 0(fS)2+(ft)2 n 0（42）3 FE-PDM方法验证为了验证基于应变控制的渐进损伤方法，对常温下三个不同的载荷情况进行了试验和模拟仿真，研究了2D C/SiC复合材料的面内拉伸、面内剪切和三点弯曲问题。3.1 面内拉伸与剪切表1和表2为面内拉伸和面内剪切两组试验各自的材料属性。在表1中，通过试验得到的应力应变曲线与模拟仿真进行对比，反演了弹性常数和损伤控制参数，如图7所示。表1为2D C/SiC复合材料的弹性模量、泊松比、应力、应变等参数。在图7（a）中，可以发现，通过FE-PDM图6内聚力单元与双线性内聚力法则Fig.6Cohesive force element and bilinear cohesion rule59航空科学技术Jun.25 2023 Vol.34 No.06方法来模拟2D C/SiC复合材料拉伸试验，试验与仿真的应力应变曲线在线性阶段吻合良好，但在非线性阶段具有较大的分散性。此外，通过数值模拟可以清楚地了解损伤系数d1的变化趋势，d1在不同阶段呈非线性增加。表2包括2D C/SiC复合材料的切变模量、泊松比、应力和应变等参数。如图7（b）所示，试验结果一致性良好，呈平滑的非线性趋势。通过FE-PDM方法，将面内剪应力应变曲线设置为双线性曲线，其中有两个拐点需要控制，并与试验进行比对。结果表明，仿真和试验的一致性很好，这表明双线性本构关系可以很简单和准确地代表试验的非线性关系，同时也不需要太多的控制点和参数。此外，模拟还给出了损伤系数d6的变化趋势，在第二阶段呈非线性增加。3.2 三点弯曲试验为了验证组合应力，采用三点弯曲试验进行综合应力示例。三点弯曲试验具有拉伸和压缩的法向应力，以及平面外的切应力。因此，该示例能够验证本文FE-PDM方法的适用性。试验装置与试样如图8所示。样品的尺寸为40mm5.2mm3.8mm（长度、宽度和高度），支架跨度为30mm。2D C/SiC复合材料的弹性常数见表332，2D C/SiC的损伤控制参数见表4。由图9可知，试样上半部分承受压缩应力，下半部分承受拉伸应力。其中，最大的拉伸应力位于试样底部中心线上。针对拉伸试验，拉伸应力和损伤系数d1与应变的关系如图10（a）所示，整体呈非线性变化；对于压缩试验，压缩应力和损伤系数d1与应变的关系如图10（b）所示，应力应变曲线呈线性趋势，损伤系数d1=0表明材料直到最终断裂才进入破坏状态。图11（a）显示了面外的切应力分布，其中试件的左半部分与右半部分受力对称。从试件的z方向来表1 2D C/SiC复合材料（拉伸）弹性常数及损伤控制参数Table 1 Elastic constants(tension)and damage-controllingparameters of 2D C/SiC 参数数值参数数值E11/GPa70mc/MPa160E22/GPa70ft/MPa250120.02int/2000130.022/2500230.023/6500tc/MPa140ft/7000表2 2D C/SiC复合材料（剪切）弹性常数及损伤控制参数Table 2 Elastic constants(shear)and damage-controllingparameters of 2D C/SiC 参数数值G12/GPa60y/MPa80b/MPa190int/1300f/133003002502001501005020406080100120140020004000600080000200040006000800010000120001400000.20.40.60.81.0?MPaMPa00.20.40.60.81.01234d1FEPDM12345d6FEPDMab图7试验与数值模拟结果对比Fig.7Comparison between experimental and simulated result表3 2D C/SiC复合材料的弹性常数32Table 3 Elastic constants of 2D C/SiC32参数数值E11/GPa33E22/GPa33E33/GPa33120.02130.02230.02G12/GPa35G13/GPa35G23/GPa35图8试验装置及试样Fig.8Experimental setup and specimen60刘斌 等：基于多线性本构与损伤耦合的叠层陶瓷基复合材料数值预测方法看，应力集中逐渐远离中心线，即裂纹不位于中心平面上，这与试验结果一致，如图12（b）所示。图11（b）显示了面外压缩应力分布，可以看出，试样压头下方的区域承受了与损伤模式六相对应的法向压力。然而，面外剪切和法向压力都并未达到破坏标准，因此，应力与应变呈线性变化，且损伤系数为零。如图12所示，通过X射线技术对三点弯曲试样中心截面进行无损CT扫描，精度为2m。CMCs的坐标如图12（a）中所示，其中轴1表示经纱方向，轴2表示纬纱方向，轴3是面外的法线方向。裂纹主要出现在垂直平面13和水平平面12。如图12（b）所示，垂直面13的CT图像显示了裂纹从试样底部向顶部的扩展过程，水平面12的CT图像表明试样下表面受拉，裂纹出现在试样的下表面。这表明试样的破坏模式表现为包含拉伸、压缩、剪切模式以及这些模式的组合，与方程（1）和图3（b）中所示的FE-PDM模型的失效假设一致。基于上述分析，通过FE-PDM方法计算得到了三点弯曲试验的损伤演化过程如图14所示；同时将仿真模拟得到的载荷位移曲线与试验曲线进行对比（如图13（b）所示），虽然试验得到的载荷位移曲线的斜率和极限荷载具有一定的分散性，但模拟的载荷位移曲线位于多条试图11面外剪切/压缩应力云图Fig.11Out-of-plane shear/compression stress nephogram图9沿经向的应力分布Fig.9Stress distribution along warp direction表4 2D C/SiC复合材料的损伤控制参数Table 4 Damage-controlling parameters of 2D C/SiC参数数值参数数值tc/MPa50fc/MPa-300mc/MPa65fc/-9000ft/MPa170y/MPa180int/1500b/MPa2002/2000int/50003/7000f/9000ft/80001801601401201008060402000.050.100.150.200.250.300.350501001502002503000.050.100.150.200.250.300.3500.51.0MPaMPad1FEPDMd1FEPDM?ab0.500.5图10拉伸/压缩应力和损伤系数与应变Fig.10Tensile/compressive stress and damage coefficient vs strain61航空科学技术Jun.25 2023 Vol.34 No.06验曲线的中间位置，并且与试验曲线趋势一致，因此试验与仿真结果吻合良好。如图14所示，三点弯曲试验主要破坏模式表现为包含拉伸、压缩、剪切模式以及这些模式的组合。FE-PDM 方法模拟了 4 个不同的破坏阶段，SDV28（损伤状态变量）表示试样处于哪一破坏阶段：无损伤（SDV28=0）、隧道裂纹（SDV28=1）、基体开裂和界面脱黏（SDV28=2）以及纤维断裂和拔出（SDV28=34）。第二阶段中（SDV28=1），隧道裂纹发生在试样中心截面的底部，相对而言，裂纹向两侧的扩展速度比向顶部的扩展速度更快。然后，损伤转变为基体开裂、裂纹饱和、裂纹偏转和界面脱黏。到达第三阶段，裂纹不断向两侧和顶部扩展，拉伸应力与应变保持线性关系，直至纤维最终断裂。仿真模拟的断裂位置位于CMCs试样的中心平面，与图14所示位置一致。因此，本文提出的FE-PDM方法能够适应载荷位移的非线性变化，并通过应力与应变的相关性进行控制。4 结论针对陶瓷基复合材料（CMCs）正交各向异性的本构关图13SEM断口图像与载荷位移曲线图Fig.13SEM fracture image and load displacement curve图12三点弯曲试样的显微CT图像Fig.12Micro-CT image for the middle part of the three-pointbending specimen 图14FE-PDM计算的损伤起始与演化过程Fig.14Damage initiation and evolution computed by FE-PDM62刘斌 等：基于多线性本构与损伤耦合的叠层陶瓷基复合材料数值预测方法系，本文提出了一种基于应变控制的有限元渐进损伤计算方法（FE-PDM），并通过面内拉伸、面内剪切和三点弯曲试验进行了验证，结论如下：（1）FE-PDM 方法包含 6 种失效模式：经向拉伸和压缩、纬向拉伸和压缩以及面外拉伸和压缩。对于这些失效模式，本文提出了基于应变控制的适用于各正交方向损伤因子相互耦合的刚度折减策略，即通过 6 个损伤系数（d1d6）来表征应力的退化，并在刚度矩阵中将各向正应力通过非对角线元素进行相互关联。（2）基于失效准则，根据CMCs的多线性本构关系，利用刚度退化的方法处理进入失效模式的单元。损伤系数方程通过应变进行控制，并考虑了不同方向上损伤系数之间的耦合关系。此外，还建立了CMCs的层间内聚力模型，用于模拟CMCs的分层损伤。（3）提出了基于试验结果与模拟结果进行对比迭代的参数反演方法。通过FE-PDM方法数值模拟得到的应力应变曲线、载荷位移曲线、损伤起始与演化过程与试验结果吻合较好，最终确定了模型的控制参数，并验证了FE-PDM方法的有效性和准确性。（4）FE-PDM方法可在商用有限元软件Abaqus及其二次开发vumat用户子程序下实现，具有控制参数少、计算速度快等优势。参考文献1 Aveston J,Cooper G A,Kelly A.The properties of fibre compositesM.Guildford,England:IPC Science and Technology Press,1971.2 Naslain R R.The design of the fibre-matrix interfacial zone in ceramic matrix compositesJ.Composites Part A:Applied Science and Manufacturing,1998,29(9-10):1145-1155.3 Cao Liyang,Wang Jing,Liu Yongsheng,et al.Effect of heat transfer channels on thermal conductivity of silicon carbide composites reinforced with pitch-based carbon fibersJ.Journal of the European Ceramic Society,2022,42(2):420-431.4 郭广达，成来飞，叶昉.航空发动机热结构部件的RMI工艺研究进展J.航空科学技术，2022，33(8):1-8.Guo Guangda,Cheng Laifei,Ye Fang.Research progress on RMI technology for thermal structural components of aircraft enginesJ.Aeronautical Science&Technology,2022,33(8):1-8.(in Chinese)5 Chateau C,Glbart L,Bornert M,et al.Micromechanical modeling of the elastic behavior of unidirectional CVI SiC/SiC compositesJ.International Journal of Solids and Structures,2015,58:322-334.6 Borkowski L,Chattopadhyay A.Multiscale model of woven ceramic matrix composites considering manufacturing induced damageJ.Composite Structures,2015,126:62-71.7 Dong Hongnian,Gao Xiguang,Zhang Sheng,et al.Multi-scale modeling and experimental study of fatigue of plain-woven SiC/SiC compositesJ.Aerospace Science and Technology,2021,114:106725.8Li Fei,Liu Bin,Zhong XiaoPing,et al.Computational prediction of SiCf/SiC stiffness and thermal residual stress by 3D micro-scale FEA methodsC/Proceedings of the 5th China Aeronautical Science and Technology Conference,2022:177-183.9 Bernachy-Barbe F,Glbart L,Bornert M,et al.Characterization of SiC/SiC composites damage mechanisms using digital image correlation at the tow scaleJ.Composites:Part A,2015,68:101-109.10Zhang Daxu,Liu Yu,Liu Hailong,et al.Characterisation of damage evolution in plain weave SiC/SiC composites using in situ X-ray micro-computed tomographyJ.Composite Structures,2021,275(8):114447.11Shiganga A,Song Weili,Chen Yanfei,et al.Stress field and damage evolution in C/SiC woven composites:Image-based finite element analysis and in situ X-ray computed tomography testsJ.Journal of the European Ceramic Society,2021,41:2323-2334.12Li Jun,Jiao Guiqiong,Wang Bo,et al.Damage characteristics and constitutive modeling of the 2D C/SiC composite:Part I-Experiment and analysisJ.Chinese Journal of Aeronautics,2014,27(6):1586-1597.13Puck A,Schurmann H.Failure analysis of FRP laminates by means of physically based phenomenological modelsJ.Composites Science and Technology,1998,62(12-13):1633-1662.14Hashin Z.Failure criteria for unidirectional fiber compositesJ.Jourmal of Applied Mechanics,1980,47:329-334.15Zhang Yi,Zhang Litong,He Jiangyi,et al.Modelling shear behaviors of 2D C/SiC z-pinned joint prepared by chemical 63航空科学技术Jun.25 2023 Vol.34 No.06vapor infiltrationJ.Ceramics International,2018,44(6):6433-6442.16Gao Xiguang,Yu Guoqiang,Xue Jiangang,et al.Failure analysis of C/SiC composites plate with a hole by the PFA and DIC methodJ.Ceramics International,2017,43(6):5255-5266.17Hashin Z,Rotem A.A fatigue criterion for fiber-reinforced materialsJ.Journal of Composite Materials,1973,7(4):448-464.18Huang C H,Lee Y J.Experiments and simulation of the static contact crush of