热固性树脂基复合材料固化变形数值模拟

Numerical Simulation of Curing Deformation in Thermosetting Resin Matrix Composites

摘要

了解复合材料的固化变形机理并有效预测和控制其变形是确保复合材料结构产品成功应用的关键。为准确预测热固性树脂基复合材料的固化变形，准确表征复合材料固化过程中的非机械应变是关键前提。文章设计一种用于直接测量复合材料的非机械应变，并结合黏弹性本构模型，对复合材料层合板及C型构件的固化过程进行数值模拟的实验装置。结果显示，复合材料固化过程中残余应力和固化度的模拟结果与文献中的实验数据吻合性较好，C型复合材料构件固化变形的模拟结果与实验值的误差约为0.08%，验证了该方法的准确性与可靠性。

Abstract

Understanding the curing deformation mechanism of composite materials and effectively predicting and controlling their deformation is crucial for the successful application of composite material structures. To accurately predict the curing deformation of thermosetting resin-based composites, one of the key prerequisites is the precise characterization of the non-mechanical strain during the curing process. The paper presents an experimental setup designed to directly measure the non-mechanical strain of composite materials. By combining a viscoelastic constitutive model, numerical simulations of the curing process for composite laminates and C-shaped components were conducted. The results indicate that the simulated residual stresses and curing degrees during the curing process are in good agreement with experimental data from the literature. The simulation error for the curing deformation of the C-shaped composite component is approximately 0.08%, thus verifying the accuracy and reliability of the proposed method.

Graphical abstract

关键词

热固性树脂基复合材料 / 固化变形 / 非机械应变 / 黏弹性本构模型 / 数值模拟

Key words

Thermosetting resin matrix composites / Curing deformation / Non-mechanical strain / Viscoelastic constitutive model / Numerical simulation

热固性树脂基复合材料因其比强度高、比模量高、耐腐蚀、抗疲劳等显著优点，在航空航天、汽车、电子等高性能应用领域得到广泛应用^[1]。复合材料的固化过程通常伴随着一系列复杂的物理与化学反应。由于材料的各向异性、热胀冷缩效应、树脂基体的化学收缩反应及温度梯度等因素^[2]，固化后内部存在的残余应力会导致固化后变形^[3]。因此，建立复合材料固化模型对精确预测构件性能具有重要意义^[4-5]。利用有限元法预测复合材料固化成型并补偿工装型面，可减少实验迭代，降低固化变形影响^[6]。

早期，研究人员提出用于预测复合材料在固化成型过程中的变形和残余应力的多种模型。WHITE等^[7]利用动态热机械分析(DMA)实验研究3051-6树脂的黏弹特性，并总结3051-6树脂的权重因子与松弛时间之间的关系。ZOCHER等^[8]提出一种黏弹性本构模型，利用广义Maxwell模型来表征固化过程中材料的松弛模量。SVANBERG等^[9]提出Path-dependent本构模型，该模型在简化传统黏弹性本构模型的基础上，假设高弹态阶段的材料转换因子为0，而在玻璃态阶段，转换因子则趋向于无穷大。近年来，国内关于黏弹性本构模型的研究逐渐增多。闵荣等^[10]基于广义Maxwell模型，构建一种能够预测固化变形和残余应力的三维数值模型。丁安心^[11]则对黏弹性、CHILE和Path-dependent三种固化力学本构模型进行比较。结果表明，在预测固化回弹角和残余应力时，黏弹性模型与Path-dependent模型在数值上的差异可以忽略不计。DAI等^[12]研究材料力学性能随温度或固化程度变化的本构模型，从简单的一维线弹性模型扩展到复杂的三维非线性模型。

在确定复合材料固化模型之后，准确表征材料的非机械应变是确保固化变形预测精确性的关键前提之一。对于复合材料热膨胀系数的研究，NELSON等^[13]考虑不同方向的热膨胀系数，准确预测L型和U型构件的固化变形。赵明等^[14]研究复合材料的宏观收缩率，发现固化过程中最大收缩率出现在厚度方向。王乾等^[15]提出，复合材料单层板在纤维方向的热膨胀系数明显小于纤维法向方向的热膨胀系数。在固化完成后的降温阶段，沿厚度方向的收缩量远大于纤维方向，这种差异会导致复合材料V型结构发生角度回弹。

针对复合材料的预测，以往文献大多从理论公式推导非机械应变，即热应变和化学应变的解析式，再结合黏弹性本构模型，这大大增加了复合材料预测的难度。为了简化复合材料固化变形过程的预测，本文设计一种测量复合材料厚度方向变形的装置，用于记录固化反应全过程中的厚度变化。将厚度变化转化为非机械应变，并进一步构建非机械应变与固化时间的关系，取代传统非机械应变解析式的推导，从而显著简化了复合材料固化变形的预测难度，并结合黏弹性本构模型，建立有限元模型，对C型试件的固化变形进行预测，并通过与实际固化变形结果的对比分析，验证了基于黏弹性本构模型的简化预测方法的准确性。

1 测量复合材料非机械应变

1.1 测量装置与方法

图1为试验装置。本装置是基于热压罐成型技术^[16]设计的。真空泵管连接外部的真空泵，可以为复合材料的固化成型提供稳定的压力；风扇可以促进内部空气流动，从而使装置内的温度分布更加均匀；保温层可以隔绝外部环境，从而实现保温效果；千分表显示复合材料的形变量；联轴器连接千分表与石墨棒；石墨棒能够将复合材料的形变量传递至千分表；多块石墨板可以起到支撑作用，选用石墨材料是由于其具有耐高温、低热膨胀系数和高强度等特性；加热片提供热源；密封袋提供密闭环境；热电偶连接外部温控表，可以调节和控制装置内部温度。

在试验过程中，首先通过温控仪设定复合材料的固化温度曲线，同时使用真空泵抽取密封袋内的空气，确保真空状态。固化过程中，风扇需持续转动，以保证装置内温度分布均匀。复合材料固化过程中的形变量通过石墨棒传递至千分表，记录千分表的读数变化，以获取相关数据。

1.2 非机械应变测量结果与分析

本文测试复合材料为碳纤维/树脂(AS4/3501-6)。将整个固化过程中多个阶段近似看成多个线性关系，将非机械应变与固化时间的关系定义为复合材料的非机械应变系数，则有：

∆{\epsilon }_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}j}={\alpha }_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}}∆t(1)

式(1)中：∆{\epsilon }_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}j}为第j个样本的非机械应变，%；{\alpha }_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}}为非机械应变系数，\mathrm{m}\mathrm{i}{\mathrm{n}}^{-\mathrm{1}}；∆t为固化时间，min。

图2为固化成型过程中的非机械应变及温度变化。升温阶段(0~50 min)：温度上升，树脂仍未固化(固化度低)，此时非机械应变以热膨胀为主，几乎没有固化收缩。固化阶段(50~150 min)：50~100 min时为保温阶段，此时固化度增加缓慢；100~150 min时，固化反应加速体积收缩明显，导致非机械应变迅速下降。恒温固化完成阶段(150~250 min)：150~200 min时，固化度接近1，非机械应变下降缓慢；200~250 min时，固化度为1，且处于保温阶段，{\alpha }_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}}趋于稳定。冷却阶段(250~300 min)：温度下降，固化完成后的材料产生的热收缩，非机械应变下降。

表1为{\alpha }_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}}与固化温度的关系。

2 理论模型

2.1 热-化学模型

固化度场和温度场作为热固性树脂基复合材料固化过程中的基本物理量，二者之间具有强耦合关系^[17]。为研究复合材料的固化行为，依据傅里叶定律建立数学模型:

\rho C_{\mathrm{p}}\frac{\partial T}{\partial t}=k_{xx}\frac{{\partial }^{\mathrm{2}}T}{\partial x^{\mathrm{2}}}+k_{yy}\frac{{\partial }^{\mathrm{2}}T}{\partial y^{\mathrm{2}}}+k_{zz}\frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}+\dot{Q}(2)

式(2)中：ρ为复合材料的密度，kg/m³；C_{\mathrm{p}}为复合材料比热容，J/(kg·K)；T为热力学温度，K；t为时间，min；k_xx 、k_yy 、k_zz 分别为复合材料在3个不同方向上的热传导系数，W/(m\bulletK)；\dot{Q}为复合材料内部产生的热量，J。

\dot{Q}由树脂固化放热产生，其表达式为：

\dot{Q}={\rho }_{\mathrm{r}}(\mathrm{1}-V_{\mathrm{f}})H_{\mathrm{r}}\frac{\mathrm{d}\alpha }{\mathrm{d}T}(3)

式(3)中：{\rho }_{\mathrm{r}}为树脂的密度，kg/m³；V_{\mathrm{f}}为纤维的体积分数，%；H_{\mathrm{r}}为单位质量树脂完全固化时释放的总热量，J/kg；\frac{\mathrm{d}\alpha }{\mathrm{d}T}为固化反应速率。

\frac{\mathrm{d}\alpha }{\mathrm{d}T}的表达式为：

\frac{\mathrm{d}\alpha }{\mathrm{d}T}=f(\alpha ,T)(4)

式(4)中：α为固化度；T为热力学温度，K。

f(\alpha ,T)的具体形式由固化动力学的性质决定。一般来说，固化动力学模型主要分为机理模型和唯象模型两类^[18]，其中唯象模型被广泛采用。3501-6树脂的固化动力学方程为：

\frac{\mathrm{d}\alpha }{\mathrm{d}T}=\left\{\begin{array}{c}(K_{\mathrm{1}}+K_{\mathrm{2}}\alpha )(\mathrm{1}-\alpha )(\mathrm{0.47}-\alpha ),\alpha \le \mathrm{0.3}\\ K_{\mathrm{3}}(\mathrm{1}-\alpha ),\alpha \mathrm{0.3}\end{array}\right.(5)

K_i=A_i\mathrm{e}\mathrm{x}\mathrm{p}(\frac{-∆E_i}{RT}),(i=\mathrm{1,2},\mathrm{3})(6)

式(5)~式(6)中：K_i为自催化模型的反应速率常数；A_i为自催化模型的频率因子，\mathrm{m}\mathrm{i}{\mathrm{n}}^{-\mathrm{1}}；∆E_i为自催化模型的活化能，J/mol；R为理想气体常数，J/(mol\bulletK)。

表2为AS4/3501-6预浸料固化动力学参数^[19]。

2.2 复合材料细观力学

在预测复合材料固化变形时，通常需要采用细观力学方法。目前，针对复合材料的工程弹性常数，已有多种解析解，包括Reuse和Voigt模型^[20]、Hashin-Shtrikman界限、周期分布微观结构模型(PMM)以及自洽场模型(SCFM)等。Reuss和Voigt模型是较为简单的模型。Voigt模型假设复合材料的应变在各向同性条件下相同，给出了弹性常数的上界；而Reuss模型假设各相应力相等，提供了弹性常数的下界。这两个模型属于极限理论，适用于均匀分布的简单复合材料^[21]。Hashin-Shtrikman界限模型只能提供一个范围，不能精确预测^[20]。PMM模型的计算公式较为复杂，且计算资源消耗较大^[22-23]。相比之下，SCFM在处理周期性微观结构时具有较高的计算效率，并能有效处理多尺度问题，因此在复合材料固化变形预测中得到广泛应用。

根据SCFM模型^[24]，可得出：

\left\{\begin{array}{l}{E}_{\mathrm{1}}=E_{\mathrm{1}\mathrm{f}}+E_{\mathrm{r}}\left(\mathrm{1}-V_{\mathrm{f}}\right)+\left[\frac{\mathrm{4}\left(V_{\mathrm{r}}-v_{\mathrm{12}\mathrm{f}}^{\mathrm{2}}\right)k_{\mathrm{f}}{k}_{\mathrm{r}}{G}_{\mathrm{r}}\left(\mathrm{1}-V_{\mathrm{f}}\right)V_{\mathrm{f}}}{\left(k_{\mathrm{f}}+G_{\mathrm{f}}\right)k_{\mathrm{r}}+\left(k_{\mathrm{f}}-k_{\mathrm{r}}\right)G_{\mathrm{r}}{V}_{\mathrm{f}}}\right]\\ E_{\mathrm{2}}=E_{\mathrm{3}}=\frac{\mathrm{1}}{{\left(\mathrm{4}{k}_{\mathrm{t}}\right)}^{-\mathrm{1}}+{\left(\mathrm{4}{G}_{\mathrm{23}}\right)}^{-\mathrm{1}}+{\left(v_{\mathrm{12}}^{\mathrm{2}}/E_{\mathrm{11}}\right)}^{-\mathrm{1}}}\\ G_{\mathrm{12}}=G_{\mathrm{13}}=G_{\mathrm{r}}\left[\frac{G_{\mathrm{12}\mathrm{f}}+G_{\mathrm{r}}+\left(G_{\mathrm{12}\mathrm{f}}-G_{\mathrm{r}}\right)V_{\mathrm{f}}}{G_{\mathrm{12}\mathrm{f}}+G_{\mathrm{r}}-\left(G_{\mathrm{12}\mathrm{f}}-G_{\mathrm{r}}\right)V_{\mathrm{f}}}\right]\\ G_{\mathrm{23}}=\frac{G_{\mathrm{r}}[k_{\mathrm{r}}(G_{\mathrm{r}}+G_{\mathrm{23}\mathrm{f}})+\mathrm{2}{G}_{\mathrm{23}\mathrm{f}}{G}_{\mathrm{r}}+k_{\mathrm{r}}(G_{\mathrm{23}\mathrm{f}}-G_{\mathrm{r}})V_{\mathrm{f}}]}{k_{\mathrm{r}}(G_{\mathrm{r}}+G_{\mathrm{23}\mathrm{f}})+\mathrm{2}{G}_{\mathrm{23}\mathrm{f}}{G}_{\mathrm{r}}-(k_{\mathrm{r}}+\mathrm{2}{G}_{\mathrm{r}})(G_{\mathrm{23}\mathrm{f}}-G_{\mathrm{r}})V_{\mathrm{f}}}\\ v_{\mathrm{12}}=v_{\mathrm{13}}=\\ v_{\mathrm{12}}{V}_{\mathrm{f}}+v_{\mathrm{r}}(\mathrm{1}-V_{\mathrm{f}})+\left[\frac{(v_r-V_{\mathrm{12}f})(k_r-k_f)G_r(\mathrm{1}-V_f)V_f}{(k_f+G_r)k_r+(k_f-k_r)G_r{V}_f}\right]\\ v_{\mathrm{23}}=\frac{\mathrm{2}{E}_{\mathrm{1}}{k}_{\mathrm{t}}-E_{\mathrm{1}}{E}_{\mathrm{2}}-\mathrm{4}{v}_{\mathrm{12}\mathrm{f}}^{\mathrm{2}}{k}_{\mathrm{t}}{E}_{\mathrm{2}}}{\mathrm{2}{E}_{\mathrm{1}}{k}_{\mathrm{t}}}\\ k_{\mathrm{f}}=\frac{E_{\mathrm{1}\mathrm{f}}}{\mathrm{2}(\mathrm{1}-v_{\mathrm{12}\mathrm{f}}-\mathrm{2}{v}_{\mathrm{12}\mathrm{f}}^{\mathrm{2}})}\\ k_{\mathrm{r}}=\frac{E_r}{\mathrm{2}(\mathrm{1}-v_{\mathrm{r}}-\mathrm{2}{\mathrm{v}}_{\mathrm{r}}^{\mathrm{2}})}\\ k_{\mathrm{t}}=\frac{(k_f+G_r)k_r+(k_f-k_r)G_r{V}_f}{(k_f+G_r)-(k_f-k_r)V_f}\end{array}\right.(7)

式(7)中：E、G分别为弹性模量和剪切模量，MPa；v为泊松比；G_{\mathrm{r}}为瞬时剪切模量，MPa；k_{\mathrm{f}}为平行于纤维方向体积模量，MPa；k_{\mathrm{r}}为树脂体积模量，MPa；k_{\mathrm{t}}为中间变量，MPa；r、f代表树脂和碳纤维；1、2、3代表方向，1是平行于纤维方向。

表3为AS4碳纤维性能。

复合材料为正交各向异性，其本构方程为：

\left(\begin{array}{c}\genfrac{}{}{0pt}{}{{\mathit{\sigma }}_{\mathrm{1}}}{{\mathit{\sigma }}_{\mathrm{2}}}\\ {\mathit{\sigma }}_{\mathrm{3}}\\ {\mathit{\sigma }}_{\mathrm{4}}\\ {\mathit{\sigma }}_{\mathrm{5}}\\ {\mathit{\sigma }}_{\mathrm{6}}\end{array}\right)=\left(\begin{array}{cccccc}{\mathit{C}}_{\mathrm{11}}& {\mathit{C}}_{\mathrm{12}}& {\mathit{C}}_{\mathrm{13}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\mathit{C}}_{\mathrm{21}}& {\mathit{C}}_{\mathrm{22}}& {\mathit{C}}_{\mathrm{23}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\mathit{C}}_{\mathrm{31}}& {\mathit{C}}_{\mathrm{32}}& {\mathit{C}}_{\mathrm{33}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{C}}_{\mathrm{44}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{C}}_{\mathrm{55}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{C}}_{\mathrm{66}}\end{array}\right)\left(\begin{array}{c}\genfrac{}{}{0pt}{}{{\mathit{\epsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{1}}^{}}{{\mathit{\epsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{2}}^{}}\\ {\mathit{\epsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{3}}^{}\\ {\mathit{\epsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{4}}^{}\\ {\mathit{\epsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{5}}^{}\\ {\mathit{\epsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}\mathrm{6}}^{}\end{array}\right)(8)

其中：

\left\{\begin{array}{l}{\mathit{C}}_{\mathrm{11}}=\frac{\mathrm{1}-v_{\mathrm{23}}{v}_{\mathrm{32}}}{\nabla }{E}_{\mathrm{1}}\\ {\mathit{C}}_{\mathrm{12}}=\frac{v_{\mathrm{21}}+v_{\mathrm{31}}{v}_{\mathrm{23}}}{\nabla }{E}_{\mathrm{1}}\\ {\mathit{C}}_{\mathrm{13}}=\frac{v_{\mathrm{31}}+v_{\mathrm{21}}{v}_{\mathrm{23}}}{\nabla }{E}_{\mathrm{1}}\\ {\mathit{C}}_{\mathrm{22}}=\frac{\mathrm{1}-v_{\mathrm{13}}{v}_{\mathrm{31}}}{\nabla }{E}_{\mathrm{2}}\\ {\mathit{C}}_{\mathrm{23}}=\frac{v_{\mathrm{32}}+v_{\mathrm{12}}{v}_{\mathrm{31}}}{\nabla }{E}_{\mathrm{2}}\\ {\mathit{C}}_{\mathrm{33}}=\frac{v_{\mathrm{32}}+v_{\mathrm{12}}{v}_{\mathrm{31}}}{\nabla }{E}_{\mathrm{2}}\\ {\mathit{C}}_{\mathrm{21}}={\mathit{C}}_{\mathrm{12}}{\mathit{C}}_{\mathrm{31}}={\mathit{C}}_{\mathrm{13}}{\mathit{C}}_{\mathrm{32}}={\mathit{C}}_{\mathrm{23}}\\ {\mathit{C}}_{\mathrm{44}}=G_{\mathrm{23}}{\mathit{C}}_{\mathrm{55}}=G_{\mathrm{13}}{\mathit{C}}_{\mathrm{66}}=G_{\mathrm{12}}\\ \nabla =\mathrm{1}-v_{\mathrm{12}}{v}_{\mathrm{21}}-v_{\mathrm{23}}{v}_{\mathrm{32}}-v_{\mathrm{31}}{v}_{\mathrm{13}}-\mathrm{2}{v}_{\mathrm{21}}{v}_{\mathrm{32}}{v}_{\mathrm{13}}\end{array}\right.(9)

式(8)~式(9)中：{\mathit{\sigma }}_i为复合材料的应力矩阵，MPa；{\mathit{C}}_{ij}为复合材料的模量矩阵，MPa；{\mathit{\epsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}j}^{}为复合材料的有效应变矩阵，其表达式为^[25]：

{\mathit{\epsilon }}_{\mathrm{e}\mathrm{f}\mathrm{f}j}^{}={\mathit{\epsilon }}_j-{\mathit{\epsilon }}_{\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{t}j}^{}(10)

式(10)中：{\mathit{\epsilon }}_j表示总应变张量。

2.3 黏弹性本构模型

线性黏弹性本构方程用于预测固化复合材料结构中残余应力松弛，其表达式为：

{\sigma }_i\left(t\right)={\int }_{-\mathrm{\infty }}^t{C}_{ij}(\alpha ,T,t-\tau )\frac{\partial {\epsilon }_{\mathrm{e}\mathrm{f}\mathrm{f}j}^{}}{\partial \tau }\mathrm{d}\tau(11)

式(11)中：{\sigma }_i为应力，MPa；\alpha为固化度；t和\tau分别为当前时间和过去时间，min。

当固化度\alpha达到参考固化度{\alpha }_{\mathrm{r}}时，考虑到高分子材料热变流特性，假设t\le0时没有发生应变，式(11)可改写为：

{\sigma }_i={\int }_{\mathrm{0}}^t{C}_{ij}({\alpha }_{\mathrm{r}},T_{\mathrm{r}},{\xi }^t-\xi \text{'})\frac{\partial {\epsilon }_{\mathrm{e}\mathrm{f}\mathrm{f}j}^{}}{\partial \tau }\mathrm{d}\tau(12)

式(12)中：{\alpha }_{\mathrm{r}}为参考固化度，AS4/3501-6的{\alpha }_{\mathrm{r}}为0.98^[26]；T_{\mathrm{r}}为参考温度，一般取值为298 K；{\xi }_t为从初始时间到时间t的累积缩减时间；{\xi }^{\text{'}}为从初始时间到时间τ的累积缩减时间，min，二者与{\alpha }_{\mathrm{T}}的关系为：

{\xi }_t={\int }_{\mathrm{0}}^t\frac{\mathrm{d}s}{{\alpha }_{\mathrm{T}}({\alpha }_{\mathrm{r}},T)}(13)

{\xi }^{\mathrm{\text{'}}}={\int }_{\mathrm{0}}^{\tau }\frac{\mathrm{d}s}{{\alpha }_{\mathrm{T}}({\alpha }_{\mathrm{r}},T)}(14)

式(13)~式(14)中：{\alpha }_{\mathrm{T}}为转换因子。

AS4/3501-6的{\alpha }_{\mathrm{T}}可表示为^[27]：

{\alpha }_{\mathrm{T}}={\mathrm{10}}^{[-a_{\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}(\frac{\mathrm{1}}{\alpha -\mathrm{1}})-a_{\mathrm{2}}](T-T_{\mathrm{r}})}(15)

式(15)中：a_{\mathrm{1}}、a_{\mathrm{2}}为常数，对于复合材料AS4/3501-6来说，a_{\mathrm{1}}=1.400 0，a_{\mathrm{2}}=0.071 2。

在本构模型中采用广义Maxwell模型。图3为广义麦克斯韦模型。在该模型中，通常以模量作为衡量麦克斯韦弹簧性能的指标，以松弛时间来评估麦克斯韦黏壶的性能^[28]。

则应力松弛模量可表示为：

C_{ij}\left(\xi \right)=C_{ij}^{\mathrm{\infty }}+(C_{ij}^{\mathrm{0}}-C_{ij}^{\mathrm{\infty }}){\sum }_m^n{W}_m\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\xi }{{\tau }_m}\right)

=C_{ij}^{\mathrm{\infty }}+\overline{C_{ij}}{\sum }_m^n{W}_{m}exp(-\frac{\xi }{{\tau }_m})(16)

式(16)中：W_m为Maxwell模型第m支的权重系数；{\tau }_m为Maxwell模型第m支的松弛时间，min；\overline{C_{ij}}=C_{ij}^{\mathrm{0}}-C_{ij}^{\mathrm{\infty }}，C_{ij}^{\mathrm{\infty }}为平衡模量，C_{ij}^{\mathrm{0}}为初始模量，MPa，与温度和固化度无关。

在CHILE模型中，BOGETTI等^[29]认为复合材料模量与α成线性关系。

C_{ij}=\left\{\begin{array}{c}{C}_{ij}^{\mathrm{\infty }}(\alpha {\alpha }_{\mathrm{c}\mathrm{1}})\\ (\mathrm{1}-{\alpha }_{\mathrm{m}\mathrm{o}\mathrm{d}})C_{ij}^{\mathrm{\infty }}+{\alpha }_{\mathrm{m}\mathrm{o}\mathrm{d}}{C}_{ij}^{\mathrm{0}}({\alpha }_{\mathrm{c}\mathrm{1}}\alpha {\alpha }_{\mathrm{c}\mathrm{2}})\\ C_{ij}^{\mathrm{0}}(\alpha {\alpha }_{\mathrm{c}\mathrm{1}})\end{array}\right.(17)

{\alpha }_{\mathrm{m}\mathrm{o}\mathrm{d}}=\frac{\alpha -{\alpha }_{\mathrm{c}\mathrm{1}}}{{\alpha }_{\mathrm{c}\mathrm{2}}-{\alpha }_{\mathrm{c}\mathrm{1}}}(18)

式(17)~式(18)中：{\alpha }_{\mathrm{c}\mathrm{1}}、{\mathrm{\alpha }}_{\mathrm{c}\mathrm{2}}分别为玻璃态点、凝胶点时固化度；C_{ij}^{\mathrm{0}}为3.2 GPa，C_{ij}^{\mathrm{\infty }}一般取值为C_{ij}^{\mathrm{0}}/\mathrm{1}\mathrm{000}。

在CHILE模型中，树脂在固化过程中的瞬时剪切模量为^[30]：

G_{\mathrm{r}}=\frac{C_{ij}}{\mathrm{2}(\mathrm{1}+v_{\mathrm{r}})}(19)

v_{\mathrm{r}}=\frac{C_{ij}^{\mathrm{0}}-C_{ij}(\mathrm{1}-\mathrm{2}{\mathrm{v}}_{\mathrm{r}}^{\mathrm{0}})}{\mathrm{2}(\mathrm{1}+C_{ij}^{\mathrm{0}})}(20)

式(20)中：v_{\mathrm{r}}^{\mathrm{0}}为玻璃态时的泊松比，一般取0.35。

3501-6树脂的松弛时间的计算公式^[27]为：

\left\{\begin{array}{l}{\tau }_m(\alpha )={\mathrm{10}}^{\{\mathrm{l}\mathrm{g}[{\tau }_m({\alpha }_{\mathrm{r}})]+[f(\alpha )-(\alpha -{\alpha }_{\mathrm{r}})\mathrm{l}\mathrm{g}(\frac{{\mathrm{10}}^{\mathrm{9.9}}}{{\tau }_{m(a_{\mathrm{r}})}})]\}}\\ f(\alpha )=\mathrm{9.134}\mathrm{7}{\alpha }^{\mathrm{2}}+\mathrm{0.608}\mathrm{9}\alpha -\mathrm{9.369}\mathrm{4}\end{array}\right.(21)

式(21)中：{\tau }_m({\alpha }_{\mathrm{r}})为松弛时间，min。

表4为3501-6树脂的松弛时间和权重系数^[27]。

将式(12)转换成应力增量方程。

∆{\sigma }_i^{t+∆t}={\sigma }_i^{t+∆t}-{\sigma }_i^t(22)

假设转换因子在很小的时间内呈线性变化，则：

\frac{\partial {\epsilon }_{\mathrm{e}\mathrm{f}\mathrm{f}j}^{}(\tau )}{\partial \tau }\approx \frac{\partial {\epsilon }_{\mathrm{e}\mathrm{f}\mathrm{f}j}^{}(\tau )}{∆t}(23)

则应力增量可表示为：

∆{\sigma }_i^{t+∆t}=C_{ij}^{\mathrm{*}}∆{\epsilon }_{\mathrm{e}\mathrm{f}\mathrm{f}j}^{t+∆t}(\tau )+∆{\sigma }_i^{\mathrm{R}}(24)

式(22)~式(24)中：∆t为极短的时间增量；{\sigma }_i^{t+∆t}为t+\mathrm{\Delta }t时刻的应力增量，MPa；{\sigma }_i^t为t时刻的应力增量，MPa；C_{ij}^{\mathrm{*}}和∆{\sigma }_i^{\mathrm{R}}为自定义符号^[10]，MPa。

其中，C_{ij}^{\mathrm{*}}、∆{\sigma }_i^{\mathrm{R}}分别可表示为：

C_{ij}^{\mathrm{*}}=C_{ij}^{\mathrm{\infty }}+\overline{C_{ij}}{\sum }_m^n{W}_m\frac{{\alpha }_{\mathrm{T}}{\tau }_m}{∆t}\left[\mathrm{1}\right.-exp(-\frac{∆{\xi }^{t+∆t}}{{\tau }_m})](25)

∆{\sigma }_i^{\mathrm{R}}={\sum }_m^n{S}_{im}^t[exp(-\frac{∆{\xi }_t}{{\tau }_m})-\mathrm{1}](26)

式(25)~式(26)中：∆{\xi }^{t+∆t}为t+\mathrm{\Delta }t时刻的缩减时间，min；n为maxwell模型总支数；S_{im}^t为历史状态变量，MPa。

式(26)中的S_{im}^t可表示为：

S_{im}^t=\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{∆{\xi }_t}{{\tau }_m})S_i^{t-∆t}+\overline{C_{ij}}{W}_m\frac{{\alpha }_{\mathrm{T}}{\tau }_m∆{\epsilon }_{\mathrm{e}\mathrm{f}\mathrm{f}j}^t}{∆t}[\mathrm{1}-\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{∆{\xi }^t}{{\tau }_m})](27)

3 模型验证

3.1 固化变形模拟过程

由于复合材料在固化过程中会经历热、化学和力学等多物理场的耦合变化，导致变形、应力及残余应力的产生，有限元软件的标准材料库通常无法直接处理这种复杂的多物理场耦合过程。通过二次开发可以自定义材料属性及其耦合关系，从而更精确地模拟复合材料的固化行为。

本文在模拟复合材料固化变形时，使用USDFLD、FILM、DISP、HETVAL、UMAT和UEXPAN等子程序。USDFLD、FILM、DISP和HETVAL子程序主要用于热传导-固化分析，而UMAT和UEXPAN则用于应力-位移分析。具体来说，UMAT定义了应力及材料的雅可比矩阵，其中应变为机械应变；而UEXPAN则用于定义复合材料的非机械应变。

3.2 模型的验证

为了准确预测复合材料的固化变形，本文建立非机械应变与固化时间之间的关系，并结合黏弹性本构模型，对材料的固化成型过程进行数值模拟分析。为验证该方法的准确性，采用WHITE等^[7]的研究结果作为对比。

图4为WHITE的有限元模型。其中a=4b，采用[0,90]_S铺层。材料AS4/3501-6的固化动力学参数及力学性能参数如表4所示。

图5为AS4/3501-6固化度与残余应力计算结果比较。从图5可以看出，本文提出的方法所得的固化度-固化时间曲线和残余应力-固化时间曲线与文献中的计算结果高度吻合，进一步验证了该方法的准确性。

3.3 预测C型构件

为了进一步验证所提方法的可靠性，本文模拟了复合材料C型构件与模具相互作用下的成型过程。图6为构件和模具的有限元模型。在该模拟中，构件的内表面与模具外表面接触，并在接触面上定义了相互作用属性，模具底部施加温度载荷，构件外表面施加0.1 MPa压力载荷^[31]。

图7为固化结束后C型构件位移云图。从图7可以看出，C型构件在固化完成脱模后发生显著变形，最大变形出现在构件上表面最右端。具体来说，最大相对位移为0.735 1 mm，而y轴方向上的最大位移为0.731 5 mm，这归因于固化结束后脱模时，构件内的残余应力释放导致形变。

表5为实验与数值模拟变形参数对比。从表5可以看出，两者具有良好的一致性，夹角误差约为0.05°，构件上下端点间的距离偏差约为0.08 mm，对应的相对误差约为0.08%。这些结果充分验证了模型的准确性和可靠性。

4 结论

为简化复合材料固化变形的预测过程，本文设计一种装置，用于测量复合材料的非机械应变，并建立了黏弹性本构模型。在此基础上，对复合材料层合板和C型构件的固化过程进行数值模拟，得出以下结论：通过模拟复合材料层合板的成型过程，并与WHITE等^[7]实验结果进行对比，发现模拟得到的固化度和残余应力曲线与实验数据高度一致，验证了该方法的准确性和可靠性。为进一步验证该方法的合理性，本文模拟了复合材料C型构件的固化过程，比较实验中构件夹角和上下端点间距离的变化。结果显示，预测值与实际值的误差仅为0.08%，充分证明该方法在复合材料固化变形预测中的有效性。

参考文献

[1] 李军, 刘燕峰, 倪洪江, 等. 航空发动机用树脂基复合材料应用进展与发展趋势[J]. 材料工程, 2022, 50(6): 49-60.

[2] 元振毅, 魏方见, 孔令飞, 等. 模具-复合材料相互作用的复合材料固化变形仿真方法研究[J]. 复合材料科学与工程, 2023(10): 10-16.

[3] 杨智勇, 左小彪, 范春浩, 等. 典型复合材料结构固化变形仿真预测与控制验证[J]. 宇航材料工艺, 2024, 54(2): 54-61.

[4] PATE S B, KORYCKI R. Numerical optimization of stress concentration in composite structures for different material arrangement[J]. Materials, 2021, 14(11): 2957.

[5] CHEN S Y, JING X S, XIE X, et al. Numerical simulations and optimization of foam filled free-curved surface composite part[J]. Materials Research Express, 2022, 9(5): 055602.

[6] 黄岗领, 盛毅, 陈志霞, 等. 复合材料构件固化变形预测及其工装型面补偿设计[J]. 纤维复合材料, 2023, 40(3): 87-93.

[7] WHITE S R, KIM Y. Process-induced residual stress analysis of AS4/3501-6 composite material[J]. Mechanics of Composite Materials and Structures, 1998, 5(2):153-186.

[8] ZOCHER M A, GROVES S E, ALLEN D H. A three-dimensional finite element formulation for thermoviscoelastic orthotropic media[J]. International Journal for Numerical Methods in Engineering, 1997, 40(12): 2267-2288.

[9] SVANBERG J M, HOLMBERG J A. Prediction of shape distortions. Part Ⅱ. Experimental validation and analysis of boundary conditions[J]. Composites Part A: Applied Science and Manufacturing, 2004, 35(6): 723-734.

[10] 闵荣, 元振毅, 王永军, 等. 基于黏弹性本构模型的热固性树脂基复合材料固化变形数值仿真模型[J]. 复合材料学报, 2017, 34(10): 2254-2262.

[11] 丁安心. 热固性树脂基复合材料固化变形数值模拟和理论研究[D]. 武汉: 武汉理工大学, 2016.

[12] DAI J F, XI S B, LI D N. Numerical analysis of curing residual stress and deformation in thermosetting composite laminates with comparison between different constitutive models[J]. Materials, 2019, 12(4): 572.

[13] NELSON R H, CAIRNS D S. Prediction of dimensional changes in composite laminates during cure[J]. International ASMPE Symposium and Exhibition, 1989, 34(2): 2397-2410.

[14] 赵明, 王怡敏, 刘伟, 等. 基于宏观线收缩率的复合材料固化变形仿真预测[J]. 工业技术创新, 2022, 9(1): 30-35.

[15] 王乾, 关志东, 王仁宇, 等. 结构参数对复合材料V型构件固化变形影响试验及解析分析[J]. 复合材料学报, 2017, 34(12): 2722-2733.

[16] DONG C C, ZHOU J X, JI X Y, et al. Study of the curing process of carbon fiber reinforced resin matrix composites in autoclave processing[J]. Procedia Manufacturing, 2019, 37: 450-458.

[17] 龙习坤, 李树健, 陈蓉, 等. 树脂基复合材料制件微波固化数值模拟[J]. 复合材料学报, 2023, 40(4): 2415-2426.

[18] 陈忠丽. 纤维复合材料结构固化变形过程的有限元模拟及在线监测[D]. 济南: 山东大学, 2020.

[19] BOGETTI T A, GILLESPIE J. Two-dimensional cure simulation of thick thermosetting composites[J]. Journal of Composite Materials, 1991, 25(3): 239-273.

[20] NEMAT-NASSER S, HORI M. Micromechanics: Overall properties of heterogenous materials[M]. Amsterdam: Elsevier, 1998.

[21] 李录贤, 王铁军. 复合材料等效弹性参数预估方法间的关系[C]//中国科协第五届青年学术年会——陕西卫星会议论文集. 上海: 同济大学, 2005.

[22] LUCIANO R, BARBERO E J. Formulas for the stiffness of composites with periodic microstructure[J]. International Journal of Solids and Structures, 1994, 31(21): 2933-2944.

[23] LUCIANO R, BARBERO E J. Analytical expressions for the relaxation moduli of linear viscoelastic composites with periodic microstructure [J]. Journal of Applied of Applied Mechanics, 1995, 62(3): 786-93.

[24] 张艺澄. 变厚度碳纤维复合材料构件固化变形预测与控制[D]. 南京: 南京航空航天大学, 2022.

[25] 姚文标. 二维机织复合材料结构件固化变形研究[D]. 天津: 天津工业大学, 2023.

[26] KIM Y K, WHITE S R. Stress relaxation behavior of 3501-6 epoxy resin during cure[J]. Polymer Engineering Science, 1996, 36(23): 2852-2862.

[27] LOOS A C, SPRINGER G S. Curing of epoxy matrix composites[J]. Journal of Composite Materials, 1983, 17(2): 135-169.

[28] 王晓霞. 热固性树脂基复合材料的固化变形数值模拟[D]. 济南: 山东大学, 2012.

[29] 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.

[30] 30] O'BRIEN D J, MATHER P T, WHITE S R. Viscoelastic properties of an epoxy resin during cure[J]. Journal of Composite Materials, 2001, 35(10): 883-904.

[31] 包正弢. 热固性树脂基复合材料构件的固化变形仿真研究[D]. 广汉: 中国民用航空飞行学院, 2019.

基金资助

国家自然科学基金(52371097)