碳纤维增强树脂（CFRP）基复合材料由于高比强度、高比刚度、抗疲劳、耐腐蚀和设计性强等优点，被广泛应用于航空航天、汽车工业、能源等领域，在工业生产中通常采用固化工艺进行成型^{[ [1]  HUI X Y, XU Y J, ZHANG W H. An integrated modeling of the curing process and transverse tensile damage of unidirectional CFRP composites[J]. Composite Structures, 2021, 263: 113681.
1 ]}。CFRP复合材料固化成型过程受到热、化学、力等多场耦合作用，同时受到树脂基体的固化收缩、材料与模具之间热膨胀系数不匹配等多种因素的影响，其内部会产生非均匀的固化残余应力场，影响成型后复合材料的承载性能^{[ [2]  ZHAO L G, WARRIOR N A, LONG A C. A micromechanical study of residual stress and its effect on transverse failure in polymer-matrix composites[J]. International Journal of Solids and Structures, 2006, 43(18–19): 5449-5467.
[3]  ANTUNES M B, ALMEIDA J H S JR, AMICO S C. Curing and seawater aging effects on mechanical and physical properties of glass/epoxy filament wound cylinders[J]. Composites Communications, 2020, 22: 100517.
2-3 ]}。因此，在CFRP复合材料的力学性能评估过程中，需要合理考虑固化残余应力的影响。

目前，围绕复合材料固化成型分析^{[ [4]  DING A X, LI S X, SUN J X, et al. A thermo-viscoelastic model of process-induced residual stresses in composite structures with considering thermal dependence[J]. Composite Structures, 2016, 136: 34-43.
[5]  LI D N, LI X D, DAI J F, et al. A comparison of curing process-induced residual stresses and cure shrinkage in micro-scale composite structures with different constitutive laws[J]. Applied Composite Materials, 2018, 25(1): 67-84.
[6]  孙勇毅, 许英杰, 唐闻远, 等. 共固化成型复合材料加筋壁板的固化变形仿真技术研究[J]. 航空制造技术, 2022, 65(4): 107-114, 120.SUN Yongyi, XU Yingjie, TANG Wenyuan, et al. Simulation of curing deformation for co-cured stiffened composite panel[J]. Aeronautical Manufacturing Technology, 2022, 65(4): 107-114, 120.
4-6 ]}和力学行为分析^{[ [7]  XU Y N, GAO Y K, WU C, et al. On design of carbon fiber reinforced plastic (CFRP) laminated structure with different failure criteria[J]. International Journal of Mechanical Sciences, 2021, 196: 106251.
[8]  GU F W, YUAN X Y, ZHU X L, et al. Numerical study of composite laminates subjected to low-velocity impact using a localized damage algorithm of Puck’s 3D IFF criterion[J]. Engineering Fracture Mechanics, 2020, 228: 106901.
[9]  CAO D F, DUAN Q F, HU H X, et al. Computational investigation of both intra-laminar matrix cracking and inter-laminar delamination of curved composite components with cohesive elements[J]. Composite Structures, 2018, 192: 300-309.
7-9 ]}方面的研究取得了一定进展，但是二者大多处于割裂状态，CFRP复合材料的力学行为分析中鲜有系统性考虑固化残余应力影响的研究。如孙亮亮等^{[ [10]  孙亮亮, 丁安心, 祖磊, 等. 三维粘弹性复合材料层合板固化残余应力研究[J]. 应用力学学报, 2016, 33(5): 779-785, 933.SUN Liangliang, DING Anxin, ZU Lei, et al. Research on process-induced residual stress for composite laminates based on three-dimensional viscoelasticity[J]. Chinese Journal of Applied Mechanics, 2016, 33(5): 779-785, 933.
10 ]}采用三维粘弹性模型对固化残余应力进行了计算，预测了复合材料层合板可能出现分层的位置。Wen等^{[ [11]  WEN Z M, GONG X L. Influence of residual stresses on the damage of composite laminates under tensile loading[J]. Applied Mechanics and Materials, 2015, 784: 361-368.
11 ]}利用增量打孔法对不同固化成型工艺下[0₂/θ₂]_s复合材料的残余应力分布进行了测量，并从理论上预测了承载条件下的复合材料损伤形貌。Bondarchuk等^{[ [12]  BONDARCHUK D A, FEDULOV B N, FEDORENKO A N. The effect of residual stress induced by manufacturing on strength on free edge of carbon–epoxy composite with [00/900]n layup[J]. Procedia Structural Integrity, 2019, 18: 353-367.
12 ]}通过数值模拟方法对残余应力进行了预测，并比较了固化残余应力对复合材料单轴加载下损伤分布的影响。Cinar等^{[ [13]  CINAR K, GUVEN I, ERSOY N. Effect of residual stress on the bending response of L-shaped composite laminates[J]. Composite Structures, 2020, 246: 112425.
13 ]}建立了固化过程和受载损伤行为分析的耦合模型，总结了残余应力对层间分层和初始破坏载荷的影响规律。然而，上述工作仍有待进一步完善。首先固化成型分析与力学性能分析的割裂，难以实现固化成型到力学性能分析的全过程预测；其次固化残余应力的精确计算是研究其影响的重要基础，而对材料性能时变演化规律的忽略难以保证固化残余应力场的计算精度；最后固化残余应力对损伤失效形式的影响规律尚需深化。

为此，本文搭建了关联固化过程与承载力学行为的一体化分析计算框架，并以L型试件受四点弯准静态加载为例进行分析和验证，实现了固化过程和力学性能分析的全过程准确预测。首先综合考虑了固化过程热–化–力多场耦合作用及纤维和树脂的化学收缩、热膨胀、时变演化性能，建立了考虑材料时变特性的固化多场耦合模型，准确表征固化残余应力的产生和演化过程。随后，将固化残余应力场作为材料力学行为计算的预定义场，分析比较了固化残余应力对材料损伤行为的影响。

CFRP复合材料固化过程的温度场计算取决于两个方面：（1）包括热压罐温度、预浸料和模具传热在内的外部热量；（2）预浸料固化反应释放的内部热量。因此，温度场计算需要综合考虑传热和固化反应耦合过程，可以采用如下所示的耦合固化放热项的热传导方程^{[ [14]  ZHU Q, GEUBELLE P H, LI M, et al. Dimensional accuracy of thermoset composites: Simulation of process-induced residual stresses[J]. Journal of Composite Materials, 2001, 35(24): 2171-2205.
14 ]}。

式中，ρ_c和C_c分别为CFRP复合材料的密度和比热容；T为温度；$\dot{q}$为树脂基体交联反应放出的热量。$\dot{q}$和C_c可表示为

式中，ρ_r和H_r分别为树脂密度和树脂单位质量固化反应释放的总热量；dα/dt为树脂固化速率；C_r为树脂的比热容；α为固化度；t为时间；V_f为纤维体积分数；C_f为纤维的比热容；K_xx、K_yy、K_zz分别为CFRP复合材料沿x、y、z方向的导热系数。其中，K_xx通过下式计算。

对于横向各向同性材料的单向层，导热系数K_yy、K_zz在垂直于纤维方向的两个方向上相等，所以得到^{[ [15]  LI X Y, WANG J H, LI S X, et al. Cure-induced temperature gradient in laminated composite plate: Numerical simulation and experimental measurement[J]. Composite Structures, 2020, 253: 112822.
15 ]}

式中，$B=2（K_r/K_{tf}-1）$；K_lf和K_tf分别为纤维纵向和横向的导热系数；K_r为树脂基体的导热系数。本文所用材料为AS4/8552复合材料，其中AS4纤维和8552环氧树脂的性能参数如表1所示^{[ [16]  WANG Q, LI T, YANG X F, et al. Multiscale numerical and experimental investigation into the evolution of process-induced residual strain/stress in 3D woven composite[J]. Composites Part A: Applied Science and Manufacturing, 2020, 135: 105913.
16 ]}。

式（2）中树脂固化速率dα/dt由以下固化动力学模型（式（6）和（7））^{[ [15]  LI X Y, WANG J H, LI S X, et al. Cure-induced temperature gradient in laminated composite plate: Numerical simulation and experimental measurement[J]. Composite Structures, 2020, 253: 112822.
15 ， [17]  ERSOY N, GARSTKA T, POTTER K, et al. Modelling of the spring-in phenomenon in curved parts made of a thermosetting composite[J]. Composites Part A: Applied Science and Manufacturing, 2010, 41(3): 410-418.
[18]  ERSOY N, GARSTKA T, POTTER K, et al. Development of the properties of a carbon fibre reinforced thermosetting composite through cure[J]. Composites Part A: Applied Science and Manufacturing, 2010, 41(3): 401-409.
17-18 ]}和材料类型确定。其中，m和n为反应级数，通过DSC试验结果拟合得到；α_C为与温度相关的临界固化度。

式中，K为固化反应速率常数；R、T和A分别为普适气体常数、绝对温度和频率因子；ΔE为反应活化能；α_C=α_C0+α_CTT，C_α、α_C0和α_CT均为相关拟合系数。与式（6）相关的固化动力学参数如表2所示^{[ [16]  WANG Q, LI T, YANG X F, et al. Multiscale numerical and experimental investigation into the evolution of process-induced residual strain/stress in 3D woven composite[J]. Composites Part A: Applied Science and Manufacturing, 2020, 135: 105913.
16 ]}。

在弹性范围内，CFRP复合材料的固化残余应力σ通过胡克定律求解。其中，弹性应变ε^e、热膨胀应变ε^th及化学收缩应变ε^sh共同组成复合材料固化成型过程中的总应变ε。

式中，C为刚度矩阵。

考虑到树脂基体的时变特性，引入CHILE 模型描述固化过程中树脂基体模量E_r随温度和固化度的演化过程，如式（10）所示^{[ [19]  JOHNSTON A, VAZIRI R, POURSARTIP A. A plane strain model for process-induced deformation of laminated composite structures[J]. Journal of Composite Materials, 2001, 35(16): 1435-1469.
19 ]}。通过式（11）~（12）计算相应温度下的泊松比υ_r和剪切模量G_r。

式中，$E_r^0$和$E_r^{\infty }$分别为树脂初始和固化完成后的弹性模量；T₁和T₂分别为树脂模量开始和停止变化时的温度；T^*=T_g–T，T_g=268+220α为树脂基体的玻璃化转变温度。

考虑树脂基体和纤维的热膨胀与化学收缩效应，引入Bogetti等^{[ [20]  BOGETTI T A, GILLESPIE J W. Process-induced stress and deformation in thick-section thermoset composite laminates[J]. Journal of Composite Materials, 1992, 26(5): 626-660.
20 ]}提出的计算方法，采用下式分别计算复合材料的化学收缩应变和有效热膨胀系数。

化学收缩应变：

热膨胀系数：

式中，1和2分别为纵向和横向方向；β₂和β₃均为垂直于纤维方向的热膨胀系数；β_1f和β_2f分别为纤维的纵向、横向热膨胀系数；β_r为树脂的热膨胀系数。

利用Abaqus用户子程序实现复合材料固化成型热–化–力多场耦合过程的计算，建立如图1所示的固化成型计算流程。首先通过热–化耦合分析模拟复合材料的传热过程，计算得到单元每一节点的温度和固化度；其次，将计算得到的各节点温度值和固化度值作为已知条件，通过热–力耦合分析计算得到固化残余应力。

将热–化耦合模型、热–力耦合模型和准静态加载模型进行集成，研究固化残余应力对准静态受载的CFRP复合材料损伤行为的影响。L型结构是飞机翼梁、翼肋等构件的常用结构形式，因此如何准确表征L型复合材料结构的准静态响应和损伤行为对结构设计具有重要意义。建立考虑固化残余应力的L型试件四点弯加载模型，如图2所示。参照标准ASTM D6415，试件采用[0]₁₆单向铺层，长度为95 mm，宽度为25 mm，上下加载滚轮的直径均为25 mm，滚轮跨距分别为50 mm和75 mm。

按照图3（a）所示的固化工艺曲线对AS4/8552复合材料进行固化成型^{[ [21]  CHEN J L, WANG J H, LI X Y, et al. Monitoring of temperature and cure-induced strain gradient in laminated composite plate with FBG sensors[J]. Composite Structures, 2020, 242: 112168.
21 ]}，树脂在固化过程中经历了复杂的化学和物理变化，性能具有显著的时变特性。第1阶段树脂基体受热从固态向黏流态转变，此阶段应力松弛时间短，无法积累应力；第2阶段以凝胶点（固化度α=0.3）为起始，树脂基体交联反应产生的热膨胀和固化收缩效应对残余应力积累产生一定贡献，但该阶段树脂基体仍处于黏流态，与纤维之间力的传递微弱，所以只积累了小部分应力；第3阶段为玻璃化转变后的保温阶段，此阶段树脂基体固化度提高较小，但固化收缩和热膨胀效应引起非机械应变的变化，因此在此阶段积累了一定的残余应力；第4阶段为降温阶段，树脂基体此时呈玻璃态，固化残余应力在降温阶段积累显著（图3（b））。

为了减少有限元网格数量，提高计算效率，根据L型试件的对称性，建立如图4所示的1/4分析模型。热–化耦合模型中，给定297 K的初始温度，并将试件的温度边界条件写入子程序，计算固化过程的温度场及固化度场。在热–力耦合分析中建立了两个分析步模拟固化过程，首先，对试件外部施加0.7 MPa的均匀压强并给定对称约束条件；然后撤去所有约束和外部载荷，对试件一侧自由度进行固定，模拟试件完成固化后的脱模过程。热–化耦合模型和热–力耦合模型中定义的单元类型分别为DC3D8R和C3D8R，均采用Abaqus/Standard进行分析计算。

通过Abaqus/Explicit分析计算L型试件四点弯准静态加载过程，在相邻层间插入厚度为0.001 mm的Cohesive单元层，如图5所示，用以表征铺层层间损伤。有限元模型中上滚轮完全固定，下滚轮以10 mm/s的速度上移。

3D Hashin损伤准则被广泛应用于复合材料层内损伤失效预测，包括纤维拉伸失效、纤维压缩失效、基体拉伸失效和基体压缩失效4种失效模式，具体形式如式（17）~（20）所示^{[ [22]  CHANG F K, CHANG K Y. A progressive damage model for laminated composites containing stress concentrations[J]. Journal of Composite Materials, 1987, 21(9): 834-855.
[23]  MOHAMMED M A, TARFAOUI M. A progressive damage modelling of glass/epoxy cylindrical structure subjected to low-velocity impact[J]. Engineering Failure Analysis, 2022, 134: 106036.
22-23 ]}。本文使用Fortran语言在Vumat子程序中实现对材料损伤起始的判断。

纤维拉伸失效（σ₁≥0）：

纤维压缩失效（σ₁<0）：

基体拉伸失效（σ₂+σ₃≥0）：

基体压缩失效（σ₂+σ₃<0）：

式中，X_T、X_C、Y_T、Y_C分别为层合板纵向拉伸、纵向压缩、横向拉伸和横向压缩强度；S_ij（i，j=1，2，3）为层合板的剪切强度；σ_i（i=1，2，3）和σ_ij（i，j=1，2，3）分别为有效法向应力和剪切应力。参考文献[ [13]  CINAR K, GUVEN I, ERSOY N. Effect of residual stress on the bending response of L-shaped composite laminates[J]. Composite Structures, 2020, 246: 112425.
13 ]中的AS4/8552复合材料相关力学参数，具体如表3所示^{[ [13]  CINAR K, GUVEN I, ERSOY N. Effect of residual stress on the bending response of L-shaped composite laminates[J]. Composite Structures, 2020, 246: 112425.
13 ]}。

通过牵引分离定律的Cohesive zone模型表征层间损伤行为，开始受载时界面表现为线弹性，当载荷增加到临界强度后，界面发生损伤，承载能力持续减弱，直至完全消失，此时界面发生脱粘失效，材料出现层间分层损伤。内聚力单元参数如表4所示^{[ [13]  CINAR K, GUVEN I, ERSOY N. Effect of residual stress on the bending response of L-shaped composite laminates[J]. Composite Structures, 2020, 246: 112425.
13 ]}，其中σ_13C和σ_33C分别表示剪切载荷下的名义应力和法向载荷下的名义应力，G_IC和G_IIC表示应变能释放率。内聚力单元的本构模型可参考文献[ [9]  CAO D F, DUAN Q F, HU H X, et al. Computational investigation of both intra-laminar matrix cracking and inter-laminar delamination of curved composite components with cohesive elements[J]. Composite Structures, 2018, 192: 300-309.
9 ， [24]  CAO D F, HU H X, DUAN Q F, et al. Experimental and three-dimensional numerical investigation of matrix cracking and delamination interaction with edge effect of curved composite laminates[J]. Composite Structures, 2019, 225: 111154.
24 ]中的描述。

复合材料固化度和固化速率随温度的变化情况如图6所示。第1次保温之前（0~50 min），材料的固化速率较低，固化度在第1次保温阶段缓慢上升至0.1（50~110 min），当再次进入加热阶段（110~140 min），固化速率显著升高并达到最大值，固化度也显著提高。在第2次保温阶段（140~260 min）和最后的冷却阶段（260~340 min），固化速率逐渐降低，材料完成固化。

L型试件固化成型脱模前后（应力释放前后）固化残余应力的分布情况如图7所示，其中，Mises、S₁₁、S₃₃和S₁₃分别为米塞斯应力、切向应力、径向应力和剪应力；1方向表示沿着L型试件R角和平板长边的方向； 2方向表示沿着宽度方向； 3方向表示沿着厚度方向。脱模后，试件内部分残余应力得到释放，应力值数量级减小，且应力释放前后，残余应力均主要分布在试件的弯曲区域，因此本文将L型试件弯曲区域的应力分布作为研究重点。

应力释放前（脱模前），L型试件弯曲部分的上下表面分别受到切向和径向的拉压应力。撤去模具约束和所受罐压外载，部分应力得到释放，切向压应力和切向拉应力（S₁₁）分别分布在弯曲部分上表面和厚度方向中性轴部位，最大值分别为12.86 MPa和27.85 MPa；试件弯曲部分与平板部分交界处受到剪应力作用，在应力释放前后该区域的内侧始终受到负向剪应力，外侧始终受到一定的正向剪应力，应力释放后（脱模后），内外侧剪应力（S₁₃）值均降低，且应力值均在1 MPa左右。

CFRP复合材料L型试件弯曲部分失效形式包括层间分层和横向基体开裂。造成层间分层的主要原因是层合板厚度方向的径向应力过大，而较大的弯曲应力则会造成横向基体开裂，因此须对CFRP复合材料层合板中的径向应力，以及由切向应力和剪应力所组成的弯曲应力进行着重考虑和分析。

L型试件损伤失效起始前弯曲部位的应力分布如图8所示。在四点弯载荷作用下，试件上下表面分别受到切向压应力和切向拉应力，中间层受径向拉应力作用，而剪应力主要集中在试件的弯曲部分与平板部分的交界处，靠近下表面部分受到反向剪应力作用。由图8可知，剪应力值（S₁₃）远小于径向应力值（S₁₁）。因此，根据应力云图结果，并结合8552环氧树脂抗拉强度的定量分析，可以评估试件的初始破坏模式。

通过选取试件弯曲区域上下表面单元积分点并进行计算后发现，引入固化残余应力使得该区域上表面的平均和最大切向压缩应力值有所下降，下表面的平均和最大切向拉应力值也降低。产生该现象的原因是，在切向拉应力预应力场的作用下，四点弯加载中的部分载荷需要承担中间层的预拉应力，因此在相同的加载条件下，固化残余应力场会分担一部分四点弯作用力，使试件上下表面的应力值相对降低。值得注意的是，图8中的最大、最小切向应力值均分布在L型试件的直板部分。如图7所示，直板部分的残余应力场分布与受四点弯载荷试件的应力场分布一致，均为上表面受径向压应力，下表面受径向拉应力，因此固化残余应力场的作用会增大试件平板部分的最大和最小径向应力值。由于弯曲部分径向残余应力值较小，拉应力和压应力均在10 MPa以内，因此固化残余应力对径向应力的分布影响较小，但是会增大受载下的径向应力值。就剪应力而言，在弯曲部分与直板连接部分的交界处，中性轴上半部分受反向剪应力，下半部分受正向剪应力，相比于未考虑固化效应的四点弯结果，考虑固化效应情况下的剪应力值更大。

L型试件损伤形貌如图9所示，其中图9（a）为不考虑固化残余应力的仿真结果，该情况下复合材料的损伤形式为层间分层，分层位置位于第9层和第10层之间；图9（b）为考虑固化残余应力的仿真结果，损伤形式为横向基体开裂和层间分层，基体开裂出现在材料的第8层和第10层，且第10层的Cohesive单元在沿厚度方向的外侧部分和直板部分均达到极限承载应力值，即其QUADSCRT值（二次牵引破坏准则系数）等于1，Cohesive单元失效被删除，试件在这两处均出现了局部分层现象。

如前所述，径向应力是影响复合材料层合板分层的重要因素，过高的弯曲应力则会引起层合板内基体开裂。考虑固化残余应力会在L型试件内部引入初始的径向拉应力，因此在相同加载条件下，试件更易出现层间分层，在引入固化残余应力后，不仅在原有的第9层和第10层出现了层间分层，在第10层和第11层之间也出现了局部分层的现象。同时，弯曲部分切向应力值减小和直板部分切向应力值增大的共同作用，使得L型试件平板部分第8层和第10层的前端出现层内横向基体开裂，并逐步扩展。

将仿真结果与试验结果进行对比，需要说明的是，本文根据文献[ [13]  CINAR K, GUVEN I, ERSOY N. Effect of residual stress on the bending response of L-shaped composite laminates[J]. Composite Structures, 2020, 246: 112425.
13 ]中的固化工艺和力学性能测试条件进行仿真建模，因此可将该试验结果（图10）作为本文计算方法的试验验证。由图10可知，损伤出现在第8层和第10层，仿真结果（图9）与试验结果基本吻合。由此可见，引入固化残余应力将直接影响复合材料受载的损伤形式和损伤位置。图11对比了四点弯试验仿真和试验（文献[ [13]  CINAR K, GUVEN I, ERSOY N. Effect of residual stress on the bending response of L-shaped composite laminates[J]. Composite Structures, 2020, 246: 112425.
13 ]中的3次试验）的载荷–位移曲线，可以看出，考虑固化残余应力的有限元模型在较小的位移下就能引起材料的损伤，失效载荷略有降低。

（1）基于固化多场耦合模型，分析了固化过程对碳纤维增强树脂（CFRP）基复合材料后续加载力学行为的影响，固化过程分析中考虑了各组分材料的时变特性，并根据固化过程中树脂和纤维热膨胀和化学收缩效应，预测了固化CFRP复合材料的非均匀固化残余应力场。应力释放后，残余应力值数量级减小，切向压拉应力分别分布在弯曲部分上表面和厚度方向中性轴部位，其最大值分别为12.86 MPa和27.85 MPa；弯曲与平板部分交界处受到剪应力作用，在应力释放前后，该区域的内侧和外侧始终分别受到负向和正向剪应力，但应力释放后，应力值降低为1 MPa左右。

（2）建立耦合固化过程和力学加载过程模型，仿真结果表明，考虑固化残余应力模型计算得到的结果为第8层和第10层内的横向基体开裂，与不考虑固化残余应力条件下的结果（出现在第9层和第10层间的层间分层）相异，但与试验结果一致。因此，引入固化残余应力可以获得更高的计算精度，且残余应力会对复合材料损伤位置和损伤模式产生显著影响。此外，引入残余应力降低了最大加载位移和失效载荷。仿真与试验结果的一致性进一步验证了本文所提计算方法的有效性和将固化过程引入加载行为计算的必要性。