热塑性预浸料赋形褶皱缺陷的全局灵敏度分析

张佳晨,郭渊,吕柄熠,校金友*,文立华*,侯晓

（1.西北工业大学 航天学院，西安 710072；2.航天海鹰(镇江)特种材料有限公司，镇江 212001；3.中国航天科技集团有限公司，北京 100048）

热塑性复合材料具有成型快、断裂韧性高、抗冲击性能好和可回收等优点，已成为航空航天和国防尖端装备的理想材料[1-2]。赋形是热塑性复合材料结构成型的第一个环节，通过设置一定的温度和压力，使预浸料贴合模具，形成复杂型面[2-4]。然而，在复杂型面赋形过程中预浸料纱线存在弯曲、滑动等行为，导致面内剪切和面外弯曲等多种变形相互耦合，进而形成宏观褶皱缺陷[5-6]，显著影响结构的成型质量和力学性能。赋形温度、压力等多工艺参数的相互耦合极大程度上影响了褶皱等缺陷的有效调控。因此，如何明确工艺参数对赋形褶皱缺陷的影响程度、实现工艺参数的精准控制成为制约热塑性复合材料高质量成型的关键[7-9]。

在赋形过程中，温度、压力等工艺参数对赋形质量的影响规律复杂，工艺参数不合适会导致预浸料产生褶皱缺陷，因此国内外学者针对热塑性预浸料成型进行了实验和数值研究。Liang等[10]在不同温度下对热塑性预浸料进行悬臂梁弯曲实验，发现温度变化会改变预浸料的弯曲刚度，进而影响结构褶皱的尺寸和大小。在此基础上，Boisse等[11-12]分析了弯曲刚度对纺织复合材料起皱的影响，结果表明弯曲刚度的增加会导致褶皱尺寸的增加。Wang等[13]和吕柄熠等[14]通过研究碳纤维增强复合材料赋形过程发现赋形温度强烈影响材料起皱行为，当温度低于树脂的熔融温度时，预浸料刚度较大，成型过程中会发生面内剪切，从而产生褶皱，当温度较高时褶皱不明显。Khan等[15]和孔令国等[6]研究了压边力对玻璃纤维平纹机织物成型的影响，结果表明纤维剪切角随压边力增加而增大。但是，目前针对工艺参数对赋形褶皱缺陷的影响主要是规律性研究，定量表征多种工艺参数对赋形褶皱影响程度的工作较少。

为了定量表征工艺参数对热塑性预浸料赋形褶皱缺陷的影响程度，首先需要针对褶皱缺陷的严重程度建立评价指标[16-17]，针对该问题研究人员已经发展了多种缺陷表征方法。Lee等[18]提出了平均半径法，利用边缘到圆心的平均距离及其标准差评估半球形样品上的褶皱程度，然而这种方法只适用于表征半球形型面的褶皱缺陷。随着研究的深入，发展出能够表征复杂型面上赋形褶皱缺陷的距离法，通过计算成型结构表面到理想表面的距离识别褶皱[19]。Viisainen等[20]利用三维数字图像相关方法，通过褶皱面积比和最大褶皱幅度表征褶皱行为。Shen等[21]使用图像分析方法对褶皱进行三维重建，通过计算褶皱形状与半球形和方盒形模具形状的偏差表征褶皱严重程度。距离法为褶皱幅度提供了简单的评价方法，但由于忽略了表面曲率，无法评估局部褶皱的严重程度。Dörr[22]提出了曲率法，将三维模型的曲面三角化，通过计算曲面的表面平均曲率衡量褶皱程度。然而，现有的缺陷表征方法多关注褶皱的单一特点，难以准确表征各种典型褶皱缺陷的形状和大小[16-17]。

全局灵敏度分析方法能够量化单变量或多变量相互作用对输出参数的贡献程度，稳定性较高，已经广泛应用于热塑性复合材料结构力学性能的定量分析。Pitchai等[23]研究了单向复合材料纤维和树脂的杨氏模量和泊松比变化对复合材料整体等效模量和泊松比的灵敏度，揭示了组分材料性质的不确定性对整体结构有效性质的影响。Balokas等[24]对碳纤维增强复合材料的等效性质进行了灵敏度分析，结果表明复合材料体积分数对结构等效刚度影响最大，基体抗拉强度对结构等效强度影响最大。Sharma等[25]研究了材料属性、纤维取向和热膨胀系数等属性对变刚度复合材料层合板屈曲温度的影响进行了灵敏度分析，结果表明横向热膨胀系数对层合板屈曲温度的贡献最大。但是，目前的灵敏度分析方法多应用于复合材料力学性能的研究，热塑性复合材料赋形工艺参数对褶皱缺陷的影响机制复杂，因此相关研究较少。

为了实现热塑性复合材料低缺陷成型，本文基于非正交本构模型建立热塑性机织物预浸料赋形仿真模型，探究赋形工艺参数对褶皱缺陷的影响规律。针对褶皱缺陷表征问题，将织物到模具表面的距离和织物表面的曲率相结合，提出一种更有效的褶皱缺陷量化表征方法。引入基于方差的灵敏度指标，研究赋形工艺参数和褶皱缺陷程度的定量关系。通过量化平纹碳纤维/聚碳酸酯(CF/PC)机织物预浸料赋形工艺参数对褶皱缺陷的影响程度，证明了灵敏度分析方法的有效性。

1 赋形褶皱缺陷量化表征及灵敏度分析

本文基于非正交本构模型进行热塑性预浸料赋形仿真模拟，获得复合材料赋形褶皱缺陷，针对赋形褶皱缺陷提出一种更合理的量化表征方法。引入灵敏度指标表征赋形工艺参数对褶皱缺陷的影响程度，提出赋形褶皱缺陷的灵敏度分析方法。

1.1 热塑性预浸料宽温域赋形有限元仿真方法

本文研究对象为平纹CF/PC预浸料，具体参数如表1所示。机织复合材料在赋形过程中的主要变形模式是经纬纱线之间的角度变化，纱线的再定向导致机织复合材料赋形过程中呈现显著的非线性和非正交各向异性材料特性。对于机织复合材料仿真有两种主流有限元方法，分别是离散法和连续法。离散法将复合材料看作一系列单胞，对单胞进行分析，这种方法不适用于整个赋形过程的分析。连续法是基于连续介质力学理论建立的方法，根据不同的研究对象可以选择次弹性、超弹性和非正交等本构模型进行赋形仿真分析。为了准确表征复合材料宽温域赋形过程中的非线性大变形行为，本文在团队前期建立的基于连续介质力学的非正交本构模型的基础上，在200～250℃温度范围内讨论了温度和压强工艺参数对褶皱缺陷的灵敏度。

表1 平纹碳纤维/聚碳酸酯(CF/PC)机织物预浸料参数[14]Table 1 Plain weave carbon fiber/polycarbonate (CF/PC)woven fabric prepreg parameters[14]

非正交本构模型的坐标系定义如图1所示，其中e1-e2是局部正交坐标系。为了追踪纱线在赋形过程中的方向变化，引入形f1-f2局部非正交坐标系和m1-m2正交坐标系。坐标轴f1和f2分别表示预浸料的经纬纱线方向，用于研究纱线拉伸行为；m1是变形过程中e1和e2的平分线；m2是m1沿顺时针方向的恒定正交坐标轴，用于研究纱线旋转行为。

图1 非正交本构模型坐标系定义Fig.1 Definition of coordinate systems of non-orthogonal models

非正交本构模型将在非正交坐标系中获得纱线的应力分量通过坐标变换转化到全局X-Y坐标下，推导出机织复合材料全局坐标系下的本构关系[14]，计算公式如下：

其中：σf1和σf2为纱线拉伸引起的应力分量；σm1和σm2为纱线旋转引起的应力分量；σXX、σXY和σYY为全局坐标系下的应力张量；α为f1和e1之间的夹角；β为f1和f2之间的夹角。

与正交本构模型相比，非正交本构模型可以追踪机织复合材料在赋形过程中纱线的走向和纤维剪切角的变化，能够更精确地描述机织物复合材料的赋形变形。

1.2 基于距离和曲率结合的量化表征方法

热塑性复合材料在赋形过程中纤维发生面内剪切和弯曲，诱发褶皱缺陷[14,26]，常见的局部典型褶皱缺陷形貌如图2所示。其中，图2(a)和图2(c)包含较大的平台型褶皱，赋形温度较低时在压边区域产生，图中1位置和2位置显示两种褶皱不同的弯曲方向；当赋形温度升高时，预浸料刚度降低，压边区域的褶皱尺寸逐渐变小，形成如图2(d)所示的波纹型褶皱；在赋形过程中，穹顶部分会产生如图2(b)所示的褶皱，尺寸小但表面变化剧烈。通常采用距离法和曲率法表征褶皱缺陷程度。距离法通过计算织物节点到模具表面的距离表征褶皱缺陷严重程度，计算方法如图3所示。假设模具表面点P'是距离节点P最近的点，P'点不一定与模具网格上的节点重合。针对预浸料每一个P点首先寻找模具表面与节点P距离最近的单元ABCD，通过式(10)定义单元ABCD的法向量nP，式中AB和AC均为向量；再通过式(11)获得节点P到模具表面点P'的距离dP，式中PP'和AP均为向量。距离法能够快速识别一般褶皱，但对于如图2(b)所示的小距离大曲率类型的褶皱，使用距离法表征褶皱程度很小，然而此类型褶皱曲率很大，无法准确识别。

图2 热塑性复合材料局部典型褶皱缺陷形貌Fig.2 Shape of local typical wrinkle defects of thermoplastic composites

图3 织物到模具的距离和织物表面曲率计算Fig.3 Calculation of fabric-to-mold distance and fabric surface curvature

为了通过表面曲率评价褶皱程度，Dörr[22]提出了曲率法。曲率法通过织物表面的曲率大小表征褶皱缺陷程度，计算方法如图3所示。第一步，如果共用节点P的单元中存在四边形单元，将其拆分为三角形单元，以便计算单元法向量。第二步，通过式(12)计算节点P处的节点法向量，先计算共用节点P的每个三角形单元计算法向量，再将共用节点P的N个三角形单元的法向量加权平均得到节点P处的法向量nP。第三步，将nP定义的平面作为节点P处的切平面，在节点P的切平面上建立局部坐标系{e1,e2,e3}。根据式(13)和(14)，在局部坐标系中采用二次曲面多项式对共享节点P的单元节点进行曲面拟合，式中a、b、c、d和e均为多项式系数，之后根据拟合曲面得到节点P处的法向量(式(15))。第四步，根据nn

Pew建立新的局部坐标系{e1,e2,e3}进行曲面拟合，经过多次拟合修正二次曲面。第五步，通过式(16)～(19)计算得到节点P的曲面曲率κP，式中KP和HP分别为节点P的高斯曲率和平均曲率，κ1和κ2分别为节点P的第一和第二主曲率。曲率法能够有效捕捉织物表面的弯曲变化，但在识别局部褶皱形状时无法识别褶皱的方向。例如，使用曲率法表征图2(a)和图2(c)中的缺陷时，1位置和2位置的褶皱曲率值相同，但是两种褶皱的走势并不相同，属于不同的褶皱类型。

为了更准确地表征赋形褶皱，本文提出一种新的褶皱缺陷表征方法。该方法综合考虑褶皱形貌的高度和曲率典型特征，结合距离法和曲率法的优势，实现任意褶皱形貌的综合表征。首先，通过式(20)和(21)对节点P处的距离和曲率进行归一化。由于不同位置的网格密度具有一定差距，若直接计算各节点的距离和曲率的平均值，在网格较密的区域距离和曲率对平均值的贡献更大，为了消除网格密度的影响，对每个节点的数据进行单元面积加权。通过式(22)得到预浸料褶皱距离大小，通过式(23)获得预浸料褶皱曲率大小。如果直接通过dP和κP进行衡量，只能说明预浸料自身的褶皱缺陷程度大小，没有建立不同预浸料褶皱缺陷之间的关系，无法比较不同工艺参数之间的预浸料褶皱缺陷程度。由于式(20)和(21)使用最大值对距离和曲率进行了归一化，导致无法通过对比不同预浸料的dP和κP衡量其褶皱大小。因此，以不同工艺参数下褶皱的距离和曲率最大值为依据，通过式(24)进行褶皱距离和曲率的全局归一化处理，通过kP定量表征预浸料褶皱缺陷的严重程度。

其中：κP,max和κP,min分别表示预浸料节点曲率的最大值和最小值；dP,max和dP,min分别表示预浸料节点距离的最大值和最小值；κmax和dmax分别表示所有工艺参数下预浸料节点曲率和距离的最大值；N为节点的数量；AP为节点P的面积，定义为共用节点P的单元面积之和。

分别用距离法、曲率法和综合表征方法量化图2中的4种局部典型褶皱，结果如图4所示。使用距离法可以表示局部褶皱的走势，但是其数值无法表征小距离大曲率的褶皱。例如图4(b) 和图4(d)中虚线框处的褶皱，虽然距离数值很小，但是此处褶皱曲率很大，此时数值无法表征褶皱的严重程度。使用曲率法可以很好地捕捉织物表面的细小褶皱，但是对于平台形状的褶皱，曲率法计算的数值无法表征平台部分。例如图4(a)和图4(c)中褶皱的平台部分，使用曲率法计算的数值均为0，但实际上存在褶皱缺陷。此外，对于图4(c)虚线框部分，曲率法只能捕捉到织物表面的弯曲，但是无法识别弯曲的方向。综合表征方法将两种方法相结合，在曲率较小的区域其数值可以表征褶皱的走势，在曲率较大区域可以通过数值变化表征织物弯曲程度，能够准确表征4种局部典型缺陷。

图4 使用距离法、曲率法和综合表征方法表征典型褶皱缺陷Fig.4 Characterization of typical wrinkle defects using distance,curvature and combined methods,respectively

1.3 赋形褶皱缺陷的全局灵敏度指标计算方法

热塑性预浸料赋形褶皱缺陷受温度和压强参数影响较大，但由于影响机制复杂，且工艺参数间存在耦合作用，现有方法无法量化工艺参数对褶皱的影响程度和工艺参数之间的耦合程度。基于Sobol指标，本文发展一种能够量化赋形工艺参数对褶皱缺陷影响程度的全局灵敏度分析方法，通过比较输入参数的改变对模型输出参数方差的贡献大小评估输入参数的影响程度。

本文有两个输入参数，分别为赋形温度T和上压边压强p，输出参数为式(24)中的因子kP。ST是赋形温度T的总Sobol指数，表示T变化对kP的影响程度。

式中：S1表示T单独变化对kP的影响程度，称为一阶Sobol指数；S12为T和p之间的相互作用对kP的影响程度，称为二阶Sobol指数；用Sp表示赋形压强p的总Sobol指数，表示p变化对kP的影响程度。

式中，S2表示p单独变化对kP的影响程度。各阶Sobol指数均介于0和1之间，值越大，表示输入参数对输出参数的影响越大。

一阶Sobol指数S1由式(27)计算得到。式中，V(kP)为kP的总方差；E(kP|X1=T)为kP给定T的条件期望，通过给定T后计算kP的均值得到；V[E(kP|X1=T)]为条件期望的方差；两个方差的比值表征了kP的方差中能被T预测和解释的百分比，即T对kP的影响程度。一阶Sobol指数S2和二阶Sobol指数S12分别由式(28)和式(29)计算得到，计算方法与S1类似。

式中：E(kP|X2=p)为kP给定p的条件期望；E(kP|X1=T,X2=p)为kP给定T和p的条件期望；V[E(kP|X2=p)]和V[E(kP|X1=T,X2=p)]分别为两个条件期望的方差。

2 结果与讨论

本节首先将仿真结果与已有的实验结果对比，验证本文所采用的赋形有限元仿真方法的正确性。然后，在不同赋形温度和压强下进行预浸料成型模拟仿真，分析赋形工艺参数对机织物预浸料赋形褶皱缺陷的影响程度，验证定量分析方法的有效性。

2.1 典型构件赋形仿真与验证

单穹顶结构赋形实验装置如图5(a)所示。模型由冲头、上压边、模具(下压边)和预浸料组成，预浸料尺寸为250 mm×250 mm，冲头直径120 mm。实验中通过环境箱控制赋形温度，通过紧固螺栓控制赋形压强。将预浸料和模具温度升高至目标温度后，控制冲头下移，到达指定位置后赋形完成[14]。根据实验建立单穹顶预浸料赋形的三维有限元模型，几何模型如图5(b)所示。ABAQUS®模型中预浸料一共有2 244个S4 R线性四边形单元和288个S3 R线性三角形单元，迭代收敛性良好，冲头、上压边和模具使用R3 D4和R3 D3刚性单元。为了避免预浸料尺寸影响，通过压边力和上压边与预浸料的接触面积计算得到上压边压强p，通过改变预浸料的纤维剪切角变化规律控制赋形温度T。将赋形温度T和上压边压强p作为模型输入，预浸料成型后，将预浸料的穹顶和压边部分的褶皱形状作为模型输出。

图5 单层平纹CF/PC预浸料的单穹顶结构赋形实验及模型示意图Fig.5 Schematic diagram of the preforming model of the single dome structure of CF/PC monolayer prepreg

为了评价有限元模型的准确性，在温度210℃、压强1.0 kPa下进行预浸料单穹顶实验仿真，得到的褶皱与实验结果对比，如图6(b)所示。在圆顶区域设置7个等距采样点如图6(a)所示。对各位置的剪切角进行采样，将模型与实验结果对比发现纤维剪切角变化模拟准确，如图6(c)所示。对比证明所建有限元模型与实验吻合良好，可以计算出预浸料的褶皱缺陷。

图6 210℃下CF/PC预浸料赋形褶皱缺陷仿真与实验对比Fig.6 Comparison of wrinkle defect in simulation and experimental results of CF/PC prepreg at 210℃

2.2 赋形褶皱缺陷的全局灵敏度分析

热塑性织物预浸料赋形过程中的褶皱是压缩和剪切变形相互作用的结果。在复杂型面成型过程中，预浸料纱线方向产生压应力，织物通过面内剪切变形进行补偿，当织物经纬向剪切角达到临界值时会发生瞬时屈曲。此时如果织物继续变形，局部剪切力和面内压应力都会增加，需要通过面外变形进行补偿，导致织物产生褶皱缺陷。选取9组不同的温度和压强作为赋形工艺参数，对CF/PC预浸料进行单穹顶仿真，赋形褶皱形貌如图7所示。从赋形褶皱形状可以看出，当赋形温度相同时，随着上压边压强的增加，褶皱逐渐减少；当上压边压强相同时，赋形温度越高，褶皱缺陷越少。

图7 不同赋形温度T和压强p下的褶皱缺陷形貌Fig.7 Shape of wrinkle defects at different forming temperatures T and pressures p

在热塑性复合材料赋形过程中，赋形温度和压强相互耦合，影响预浸料的褶皱形貌。当赋形温度较低时，树脂黏度大，预浸料面内剪切刚度大，在赋形过程中不易产生面内剪切，此时预浸料变形行为主要有压缩行为主导。在预浸料上表面施加压强能够增大预浸料层内摩擦和预浸料与模具之间的摩擦，在赋形过程中能够产生张力补偿压应力。随着压边力逐渐增大，对预浸料施加的张力也增大，面内压应力降低，面外变形产生的褶皱减少。随着温度升高，树脂黏度降低，引起剪切变形阻力降低，预浸料在赋形过程中更容易发生面内剪切变形，对压应力的补偿更充分，此时较小的压边力就能有效抑制褶皱缺陷。

目前成型工艺参数对褶皱形貌的影响多为规律性探究，难以准确分析各工艺参数对赋形褶皱缺陷的影响程度，需要对其进行深入表征和量化分析。因此，本文针对赋形工艺参数对褶皱程度的影响进行灵敏度分析。PC的熔融温度为220～230℃，温度过低时预浸料无法成型，温度过高时经纬纱线容易引起间隙，因此设置温度区间为200～250℃，设置压强区间为0.2～2.0 kPa；采用Sobol方法对T和p进行随机均匀采样，获得不同的工艺参数组合。

预浸料的赋形温度与剪切角变化一一对应，因此可采用剪切角变化规律表征模型温度输入。进行褶皱缺陷的灵敏度分析需要获得采样点处预浸料的褶皱，涉及多个温度，因此对210℃、220℃、230℃、240℃和250℃温度下CF/PC预浸料偏轴拉伸实验获得的纤维剪切角变化曲线进行拟合，得到拟合曲面[27]，如图8所示。

图8 不同温度下CF/PC预浸料的纤维剪切角变化拟合曲面Fig.8 Fitting surface for fiber shear angle variation of CF/PC prepreg at different temperature

将赋形温度和压强作为输入变量，将成型后CF/PC预浸料的褶皱程度kP作为输出变量，计算得到不同工艺参数下kP的变化如图9所示。通过灵敏度分析得到的Sobol指数如图10所示。其中First-order对应变量的一阶Sobol指数，表征单个变量对输出响应方差的独立贡献；Total对应变量的总Sobol指数，表征单个变量对输出响应方差的独立贡献及与其他变量间的共同贡献之和。结果表明：在CF/PC复合材料成型过程中，赋形温度对褶皱缺陷的一阶Sobol指数为0.50，总Sobol指数为0.74；赋形压强对褶皱缺陷距离的一阶Sobol指数为0.20，总Sobol指数为0.45。在热塑性复合材料赋形过程中，赋形温度和压强对赋形褶皱缺陷都有较大影响，其中温度的贡献更大。此外，两个参数的总Sobol指数与一阶Sobol指数相差较大，说明在复合材料赋形过程中存在较强的双参数耦合作用。

图9 不同赋形工艺参数下CF/PC预浸料的褶皱缺陷程度Fig.9 Degree of wrinkle defects in CF/PC prepreg under different forming process parameters

图10 CF/PC预浸料赋形工艺参数对褶皱缺陷的影响程度Fig.10 Influence of forming process parameters of CF/PC prepreg on wrinkle defects

为了验证提出的灵敏度分析方法的有效性，探究赋形温度和压强变化时褶皱程度kP的变化幅度。以220℃、1.0 MPa为基准，赋形温度和压强分别变化±5%和±10%，预浸料褶皱程度kP的变化幅度如图11所示。可以看出，在4种情况下，赋形温度引起kP的变化幅度均大于赋形压强，与分析结果吻合。

图11 褶皱缺陷程度随赋形温度和压强的变化幅度Fig.11 Variation amplitude of the degree of wrinkle defects with the forming temperature and pressure

3 结 论

本文基于非正交本构模型构建热塑性预浸料赋形仿真模型，获得复合材料赋形褶皱缺陷。针对热塑性预浸料赋形褶皱缺陷，将织物到模具表面的距离和织物表面的曲率相结合，提出了一种定量评估褶皱缺陷严重程度的量化表征方法。基于灵敏度指标建立了赋形工艺参数与褶皱缺陷之间的对应关系，发展了一种热塑性预浸料赋形褶皱缺陷全局灵敏度分析方法，探究了工艺参数对赋形褶皱缺陷的影响程度。得到的结论如下：

(1)提出的褶皱缺陷综合量化表征方法能够准确表征预浸料赋形褶皱缺陷。使用不同表征方法对宏观典型褶皱缺陷进行评价，结果表明综合表征方法对各种典型缺陷评价具有更好的效果；

(2)提出的赋形褶皱缺陷全局灵敏度分析方法能够定量分析不同赋形工艺参数对褶皱缺陷的影响程度及不同工艺参数之间的耦合作用，为进一步优化赋形工艺提供理论依据。通过探究赋形温度和压强对平纹碳纤维/聚碳酸酯(CF/PC)预浸料单穹顶赋形褶皱缺陷的影响程度，验证了该方法的有效性。