熱塑性預浸料賦形褶皺缺陷的全局靈敏度分析

張佳晨,郭淵,呂柄熠,校金友*,文立華*,侯曉

（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 賦形褶皺缺陷的全局靈敏度指標計算方法

熱塑性預浸料賦形褶皺缺陷受溫度和壓強參數影響較大，但由于影響機制復雜，且工藝參數間存在耦合作用，現有方法無法量化工藝參數對褶皺的影響程度和工藝參數之間的耦合程度?；赟obol指標，本文發展一種能夠量化賦形工藝參數對褶皺缺陷影響程度的全局靈敏度分析方法，通過比較輸入參數的改變對模型輸出參數方差的貢獻大小評估輸入參數的影響程度。

本文有兩個輸入參數，分別為賦形溫度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)預浸料單穹頂賦形褶皺缺陷的影響程度，驗證了該方法的有效性。