纤维增强复材C 形圆台壳件固化变形的预测方法

陶裕梅,郑子君,邵家儒

（重庆理工大学 机械工程学院,重庆 400054）

复合材料因其较高的比强度，在航空、车辆、船舶领域的应用越来越广泛。热压罐工艺是纤维增强复合材料的一种常见成形工艺，温度、压力场均匀可控，可同时处理多个模具，具有成形质量好、效率高的优势。

纤维增强复合材料的固化是非常复杂的物理化学过程，由于材料化学收缩、热变形、模具坯料相互作用等因素，脱模后工件必定会产生和模具形状的偏差，即固化变形。预估固化变形的趋势和大小，对优化工艺设计、提高制造质量有着重要的意义。在实验研究方面，肖光明等在大量C 形和L 形工件的固化实验的基础上，建立了变形量的响应面[1]。王雪明等对大量曲面构件曲率与成型缺陷的统计分析发现，分层等拐角缺陷与曲率半径成负相关[2]。王乾等通过实验研究了纤维含量、富树脂厚度对V 型零件成形回弹的影响，并建立了回弹预测有限元模型[3]。李桂洋等通过实验验证了经过优化成形温度和压力等工艺参数可提高制件内部质量和承载能力[4]。胡海晓通过预埋光纤布拉格光栅（fiber Bragg grating，FBG）传感器，选取层板、L 型和C 型试样，追踪固化过程中的固化度和模量变化[5]。数值模拟方面，邹尧等考虑了化学收缩和热应变对长桁类复合材料零件固化变形的影响，并提出型面补偿方法，将成形误差降低到4%以内[6]。朱海华等采用热结构耦合分析的方式模拟了固化过程，研究热导率各向异性对温度分布的影响[7]。董丰路等考虑了热载荷下复合材料杆件的固化变形，提出杆件截面设计应具有对称性以避免过度弯曲[8]。Bellini 等采用线弹性模型计算了不同铺层顺序下，热变形导致的层合板翘曲[9]。Johnston 等指出固化过程中树脂力学性质与固化度和温度有关，提出了CHILE 本构模型[10]。梁群等采用CHILE 模型计算了大尺寸壳体的变形梁，并估算得出残余应力可达兆帕级[11]。Ersoy 等根据固化变形时间相关性不明显的现象，提出了属性跃变的path-dependent 模型[12]，Ding 等在此基础上进行了改进[13]。黄尚洪等结合了热化学模型和黏弹性本构关系对固化过程进行模拟，并与光纤光栅传感器的检测结果进行了对比验证[14]。有学者还对比了树脂固化模拟中，线弹性、黏弹性、CHILE 模型和path-dependent 模型等本构模型的适用性，表明黏弹性模型的计算结果最好，但相应的参数测定和建模难度增加[13,15-16]。Jiao 等对近年来固化变形机理、模拟方案和控制方法上的认识进行了总结[17]。

研究表明，导致固化变形的应力大部分是在玻璃化后降温阶段产生的热应力[18]。在变形预测的理论研究方面，常采用热弹性力学模型来推导典型形状的热致固化变形量理论解。由于层合平板和C 形圆壳的几何简单，工程中常见的L 形和U 形等工件可以视为两者的组合。Kim 等[19]和Abouhamzeh等[20]分别推导了层合板在固化后的翘曲变形。Radford 等提出了C 形圆壳的热回弹角度预测公式[21-22]。Radford 等得到的公式基于弹性本构且只考虑残余正应力，Wisnom 等[23]和Ersoy 等引入了剪切应力和材料相变的影响[12]。Ding 采用广义平面应变模型，进一步考虑了面外应力的影响[24]。上述理论研究文献考虑了常见定截面L 形、U 形和C 形零件，预测效果良好，但实际生产中截面渐变的工件更常见，以小型固定翼飞机为例，机翼前缘圆弧、机头罩、机身后部覆盖件等均可以看作是截面渐变的C 形工件。这类工件的固化变形主要借助有限元软件进行数值分析，并往往需要开发自定义材料库，实现较为繁琐。

C 形圆台壳件可以认为是C 形圆壳件的推广，其沿母线方向上，截面半径线性改变，常可作为宽度渐变结构的覆盖件，使用范围比后者更广。在设计初期若能对此类工件的固化变形量有合理估计，对通用航空领域提高制造精度、缩短工艺开发周期有着重要的意义。本工作研究纤维增强复合材料C 形圆台壳件的固化变形规律，以得到能快速预测此类零件变形量的方法。

1 C 形圆台壳固化变形量理论分析

1.1 理论公式推导

图1 为纤维增强复合材料的典型热压罐固化工艺曲线。由图1 可见，在固化度X低于凝胶点固化度Xc之前，复合材料的树脂基体处于黏流态，不能形成明显的残余应力；在凝胶点和玻璃化点之间，材料处于橡胶态，弹性模量较小，但能够因化学收缩和热膨胀在内部产生一定的应力；随着固化进行，材料弹性模量逐渐增大，当玻璃温度Tg高于工艺温度时，材料完全玻璃化，模量达到最大，在此阶段构件内部能够产生明显的应力。因此导致固化变形的应力大部分是在玻璃化后降温阶段产生的热应力[18]。

图1 固化温度工艺曲线以及固化参数随时间的变化Fig.1 Curing temperature process curves and changes of curing parameters over time

基于上述论断，对于C 形圆壳（图2（a）），Radford等提出了经典的热回弹角度预测公式[21]：

图2 两种C 形构件模型示意图（a）圆壳模型；（b）圆台壳模型Fig.2 Schematic diagrams of two C-shaped component models（a）round shell model；（b）round pedestal shell model

式中：ατ、αn分别是C 形圆壳模型环向和法向的热膨胀系数；θ为圆心角（见图2（a））。

文献[22]认为，固化反应导致的化学收缩也服从相同的机制，于是式（1）可以改写为：

式中：βτ、βn分别是C 形圆壳模型环向和法向的化学收缩应变。

由于式(1) 和式(2) 高度的简洁性，至今仍然是估算C 形圆壳角度变形最重要的“拇指法则”。

本工作推导式（1）、（2）的更一般情形，即半径随母线长度变化时的情形。此时模型可称为C 形圆台壳，如图2(b) 所示。建立正交坐标系，x、y、s轴分别沿着母线、厚度和水平截面内弧长方向。参考文献[20，22]，对C 形圆台壳件的固化变形做如下假设：（1）变形主要源于降温阶段的热变形，此阶段工件变形是线弹性的；（2）壳体是沿各方向均匀铺层的复合材料，宏观上是横观各向同性（transversely isotropic）的；（3）顶底半径差和厚度相对于圆台半高处半径而言是小量，变形量相对于原始尺寸是小量；（4）工件变形前后均为C 形圆台壳，且各材料点在xy截面内的应变状态相同；（5）忽略剪切变形的影响。

显然，此类工件的几何形状可以用5 个参数描述，即半高处半径γ、半顶角ϕ、母线长度h、厚度b以及圆心角θ。根据假设4，描述固化变形只需给出这5 个参数的相对变化，这里采用虚功原理推导。

由假设4，工件厚向和母线方向的应变分别为：

又在xy平面上归一化坐标为(x,y)的环向纤维(其中x,y∈(−0.5,0.5))，原始长度为：

利用全微分公式，可得变形后的环向应变为：

本问题中总应变由热变形以及应力造成，在忽略泊松比的影响时有：

式中：Eτ、En分别是横观各向同性材料中横向和法向的杨氏模量。由式(6) 可以解得正应力分量 σx、σy、σx。另一方面，虚应变分量δεx、δεy、δεs可通过在式（3）、（5）中用虚增量δb、δh、δr、δθ、δϕ代替参数增量∆b、∆h、∆r、∆θ、∆ϕ得到：

由假设5 不考虑切应变分量。于是内力虚功为：

由式(5)可见式(8)被积函数的分母中出现了变量x,y，不利于积分。根据假设3，hsinϕ ≪r,b≪r，在计算环向应力 σs时，可以将式(5)的分母部分替换为其一阶近似：

此时式(8)的被积函数已是x,y的多项式，积分得到简化。可以算得：

其中：

根据虚功原理，内力虚功δWin等于外力虚功δWex。本问题中只考虑热应变导致的变形，因此外力虚功为零。对比式(10)中各虚增量的系数，可以得到式(11)中的5 个式子均为零。求解该代数方程组即得：

可以看到式(12)的结果并不包含材料的弹性参数。特别指出：虽然在式(6)中假设了泊松比为零，但这仅仅是本工作为了式(11)的简洁性而进行的简化。如果在式(6)中考虑泊松比，式(11)将非常复杂，但是最终得到的结果仍然是式(12)，即固化变形量也不包含弹性参数的影响。

当ϕ=0时，工件退化为了C 形圆壳件，此时的式(12)指出半顶角不再改变，母线、半径和厚度的相对变化量分别为对应方向的热应变值，而圆弧角度回弹量即为式(1)。因此C 形圆壳件的固化预测式可以视作本工作公式的特殊情形。当需要进一步考虑化学收缩应变的影响时，也可模仿式(2)修改预测公式：

在纤维增强复合材料固化过程的升温阶段，树脂处于黏流态或橡胶态，模量较小，不易形成较大的残余应力，对固化变形的影响很小；而在保温—降温阶段，化学应变和热应变均是收缩的，且法向由于缺乏纤维的支持，变形量大于横向，即整体上有βτ<βn<0,ατ∆T<αn∆T<0。代入式(13)可知，固化后C 形圆台壳工件的母线缩短、厚度变小、半径变小、半顶角减小、而圆心角增大。

1.2 模拟正交实验验证

由式(12)、(13) 可知，C 形圆台壳的固化变形量与多个几何参数有关。为验证本理论计算的合理性，采用有限元方法对C 形圆台壳件的热致固化变形进行正交模拟实验（图3）。注意到模型尺寸的等比例放缩不会影响结果，可固定壳厚为1 mm，从而只需考虑4 个变量：圆台半高处半径r、半顶角ϕ、母线长度h以及圆心角θ。采用四变量三水平正交实验设计，参数取值如表1。复合材料类型为玻璃纤维增强树脂，采用的正交织物材料可视为横观各向同性的，其等效材料参数如表2[25]。根据此类复合材料的典型工艺曲线（图1），取整个固化降温阶段的总温差为107 ℃。

表1 正交实验变量及对应水平选择Table 1 Orthogonal test variables and corresponding level selections

表2 玻璃纤维增强树脂基复合材料（Hexcel7781/LY5052/ HY5052,Vf=0.49）的等效力学性能Table 2 Equivalent mechanical properties of glass fiber reinforced resin-based composite material（Hexcel7781/LY5052/HY5052,Vf=0.49）

图3 C 形圆台壳的1/2 有限元模型Fig.3 Symmetrical finite element model of C-shaped circular pedestal shell

根据对称性，建立圆台壳的1/2 对称模型，并采用六面体实体单元，结构化网格剖分。在对称面上额外施加简支撑。旋转单元坐标系使之与图2(b)中的xys坐标系重合。对结构施加均匀温度载荷，进行考虑几何非线性的弹性计算。

计算完成后取6 个关键点（图3），折算有限元模拟的工件几何参数改变量，计算结果与预测公式(12)对比，结果如图4 所示。可以看出，在不同几何参数组合和相当大的参数范围内，理论解和数值解基本吻合，由图4 还可以看到，有限元计算的变形量总是略大于本工作计算结果，这是因为本工作对工件变形后仍为圆台面的假定增加了结构的刚性。由于半顶角的增量 ∆ϕ的绝对值较小，且由多个几何参数间接计算得到，与有限元模拟的差距最大。从正交实验还可以看出，当圆心角和半径较大（因素C、B 取水平3）时，对变形模式的假设与真实情况偏差越大，误差更加明显。

图4 C 形圆台壳件几何参数变量的公式预测结果与有限元对比Fig.4 Comparison of formula prediction results of geometric parameter variables of C-shaped circular truncated shell and finite element

总体而言，本工作提出的理论解与有限元模拟结果吻合良好，表明了其在预测C 形圆台壳件热致固化变形时的正确性和适用性。

2 基于path-dependent 模型的固化模拟方案

前述讨论考虑了热压罐固化工艺中已固化的复合材料由于降温阶段的温差导致变形（图1）。这是因为铺层较厚时，玻璃化前的材料模量较低，其化学收缩和热应变不会形成较大的应力；然而文献表明，当铺层较薄时，工件在玻璃化前受到模具的约束较强，也能够形成明显的应力，从而影响固化变形[12]。若要进一步考虑玻璃化前橡胶态阶段的化学-热变形的影响，可以基于path-dependent 本构模型进行有限元模拟[13-15]。如图5 所示，该模型常忽略材料在黏流态的弹性模量，而认为在橡胶态和玻璃态时的切线弹性模量分别为常数；当前温度T低于玻璃化温度Tg时，材料立即玻璃化，切线模量发生突变，但应力仍连续变化，该本构关系的增量表达式为：

图5 path-dependent 本构模型示意图Fig.5 Schematic diagram of path-dependent constitutive model

式中：Dr,Dg分别为橡胶态和玻璃态下复合材料的四阶弹性张量。

通常基于有限元软件二次开发自定义材料的方式来实现上述本构模型[12-13]，增大了建模分析的难度。这里通过将固化度场与力学场解耦，使用带初应力的热弹性分析来等效实现path-dependent 模型。

当复合材料铺层较少时，可以忽略温度梯度，假设整个工件内的温度分布均匀。这在传热学上称为集总热容法，其适用条件为毕渥数Bi满足：

式中：hf是对流换热系数；λ是固体的导热系数；δ0是固体中心到表面的特征长度，这里取复合材料厚度的一半。例如按热压罐内对流换热系数h=85.2 W/(m2•K)[26]，玻纤复合材料的平均导热系数0.6 W/(m•K)估算[27]，则当铺层厚度低于1.4 mm时，采用上述假设是合理的。

由于固化度的变化率只是当前固化度X和温度T的函数[28]：

可知工件各处固化的程度相同。根据DiBenedetto 公式[27]，玻璃化温度只与固化度有关：

式中：Tg0,Tg∞分别是尚未固化和完全固化后的玻璃化温度，与 λ一样是材料常数。由式(17) 易知，工件各处将同时玻璃化，因此有限元模拟时可以根据复合材料所处的状态分步处理，流程如下：

首先，在给定工艺温度曲线T(t)后，利用式(16)、(17)求解固化度和玻璃化温度曲线，找到凝胶点(XC=X0)和玻璃化点(T=Tg)，如图1 所示；

其次，将工件一面固定在模具上，施加从凝胶点到玻璃化点之间由于升温和固化导致的本征应变，并在自由面施加热压罐内压力，使用橡胶态的弹性参数计算变形和应力（表1 和表2）；

再次，以初始应力和初始缺陷的形式继承橡胶态阶段的应力和变形，进一步施加从玻璃化点到工艺结束由于降温和固化导致的本征应变，改用玻璃态的弹性参数计算变形和应力；

最后，移除模具的约束，使用玻璃态的弹性参数计算回弹量和残余应力。

上述模拟方案避免了引入自定义材料二次开发，更加容易实施，并且当忽略固化度和温度在工件内的空间梯度时，与常用的path-dependent 本构模型[12-13,15]是完全等效的。

为验证本工作提出的简化模拟方案的可行性以及计算效率上的优势，建立文献[12-13]中C 形模型的固化变形分析，并与实验测量的结果进行对比（图6）。结果表明，本工作提出的简化方案的回弹角与文献和实验测量值吻合良好。其中厚度为3 mm 的C 形圆壳的计算时间与文献[13]对比如表3 所示。由于本工作的简化实现方案避免了使用非线性自定义材料，在网格数同为2400，时间步数同为218 步的情况下，相比原始path-dependent方案减少时间代价约80%。

图6 C 型构件回弹预测结果与实验、文献对比（a）（90）n（n=4、8、12、16）；（b）（90/0）ns（n=1、2、3、4）Fig.6 Comparison of spring-in prediction results of C-shaped model with experiments and literatures（a）（90）n（n=4,8,12,16）;（b）（90/0）ns（n=1,2,3,4）

表3 3 mm 厚C 形模型的不同本构模型运行时间对比Table 3 Comparison of time costs of different constitutive laws for C-section with 3 mm thickness

3 工程案例与实验

某小型固定翼飞机的机头罩安装后存在明显的缝隙，为验证本工作提出的公式和有限元模拟方案的可行性，对机头罩进行成型预测及装配约束的变形分析。图7 为机头罩模型，模型长775 mm，截面近似为C 形圆槽，截面尺寸沿长度方向逐渐变化，前端槽宽900 mm，后端宽1125 mm，可以近似为C 形圆台壳模型。零件采用玻璃纤维环氧织物复合材料LY5052/7781，性能参数如表2，单层正交织物预浸料厚0.25 mm，铺层方案为[45◦/0◦/−45◦/90/40◦]，其中0°纤维方向为垂直于机头罩横轴线方向。成形时热压罐内压力为0.7 MPa，温度工艺曲线如图1 所示。在常温下脱模完毕后，将零件对称轴紧靠在模具上，此时零件侧沿与模具间有明显的间隙。取侧沿上的三点（如图8 所示），采用游标卡尺测量了水平位移量。再分别通过热弹性有限元模型、本工作简化的path-dependent 模型以及本工作理论预测公式对机头罩的固化变形量进行估算，并与实际测量结果进行对比。

图7 某小型固定翼飞机机头罩几何模型Fig.7 Geometric model of nose cover of a small fixed-wing aircraft

图8 机头罩1/2 模型固化变形云图及实验测量点Fig.8 Curing deformation contour and experimental measurement points of symmetrical model of nose cover

由于结构的对称性，有限元分析时只取右半部分。根据固化温度工艺曲线，求解玻璃化转变温度Tg以及固化度X随时间的变化如图1 所示。取XC=0.5[28]，可算得橡胶阶段和玻璃阶段的固化温差分别为∆Tr=−10oC,∆Tg=−97oC。基于热弹性模型的有限元模拟采用玻璃态弹性参数，将橡胶态和玻璃态下的温度应变和化学收缩之和作为本征应变施加；基于path-dependent 本构模型的有限元模拟方案则如前所述。为了移除刚体位移，回弹阶段在对称轴的两端施加简支撑。

采用理论解公式进行预测时，由于机头罩模型的截面并非标准的C 形圆弧，将模型用两个顶角相同、并排相切的圆台壳来拟合，拟合参数如表4。再采用式（12）或式（13）分别计算两段圆壳的变形量，拼接后得到整体的变形。最后用模具形状叠加上变形量得到零件的最终形状。

表4 机头罩简化几何模型参数Table 4 Simplified geometric model parameters of nose cover

4 结果分析与讨论

在各预测变形模型中提取侧沿的最终形状，并将侧沿在xy平面上的坐标与实验结果进行对比，结果如图9 所示。由图9 可以看到，各种模型预测的机头罩变形形式相近，都表明轮廓将内移，变形量沿x正方向逐渐减小的特点，理论解模型、热弹性模型和path-dependent 模型预测的机头罩半跨长平均缩小量分别是8.1 mm、7.6 mm、6.1 mm，均与实测值7.7 mm 基本接近。总体而言，基于pathdependent 模型的有限元解考虑了相变导致的力学性质改变，实现复杂度相对较高，但只在测量点1 明显优于其他解。基于热弹性模型的有限元解未考虑橡胶态与玻璃态的力学性质差异，高估了固化时的残余应力，变形量更大。理论解式(13)与热弹性模型基于相同的假设，因此两者的结果非常相近，在模型前端的差异主要是因为工件并非标准的圆台壳。而理论解式(12)没有考虑橡胶态阶段化学收缩对固化变形的影响，这在铺层较少的情况下将明显低估变形量[12]。

图9 不同变形轮廓与实验测量结果的对比Fig.9 Comparison of different deformation profiles with experimental measurement results

壳罩类零件起到保护内部控制系统和电器装置，维持飞机整体的空气动力学外形的作用，其外观形状非常重要。为进一步研究固化变形对壳罩装配后形状的影响，模拟装配过程，提取用pathdependent 模型算出的已有固化变形的零件，将其对称轴与侧沿固定在原设计位置处，此时的零件形状与原设计曲线对比如图10 所示。

由图10 可以看出，在此种约束下机头罩模型的构型与设计曲面有一定的偏差。图中标记处为铆钉布置位置，可见AB 段在设计曲面下方，而BC 段在设计曲面上方。实际装配完成后的照片如图11 所示。考虑到与机头罩铆接的发动机罩是较矮的穹顶形，固化变形较小，其轮廓可以作为原始设计的参照。可以看到机头罩的AB 段确实略低于发动机罩，而BC 段明显更高，两者间有明显的缝隙。这从侧面验证了本工作有限元模拟方案的正确性，也说明复合材料的固化变形对装配后的飞机外观有着明显的影响。

图10 机头罩1/2 模型装配对称轴与侧沿后的变形云图Fig.10 Deformation contour of half of nose cover model after assembling symmetry axis and side edges

图11 某小型固定翼飞机机头罩安装后照片Fig.11 Photo of a small fixed-wing aircraft after the installation of the hood

5 结论

（1）采用虚功原理推导出C 形圆台壳热致固化变形的理论预测公式。公式表明：固化后C 形圆台壳工件的母线变短、厚度变小、半径变小、半顶角减小、圆心角增大。

（2）为验证提出公式在各几何模型参数组合下的适用性，选择4 个关键变量进行三水平的正交实验，并建立相应的有限元模型进行对比，两种方法的计算结果一致性较好，表明公式可以用于初步估计复合材料的固化变形量。

（3）对于铺层厚度较小的复合材料，提出基于path-dependent 模型的简化模拟方案，相比原始path-dependent 模型实现更简单，计算时间可减少约80%。以近似C 形圆台壳的某飞机机头罩为例，分别采用本工作公式、热弹性有限元模型、简化的path-dependent 模型进行固化变形模拟，均与实测变形量基本吻合。

（4）采用有限元方法模拟了零件固化变形对装配后最终形状的影响，模拟结果能够解释实际制造中发现的装配缝隙，这可能对飞机的气动特性造成影响，因此分析复合材料固化变形是非常必要的。