李 晨，黃 甲，陳 程，高麗敏

（1.中國(guó)商飛北京民用飛機(jī)技術(shù)研究中心民用飛機(jī)結(jié)構(gòu)與復(fù)合材料北京市重點(diǎn)實(shí)驗(yàn)室，北京 102211；2.上海理工大學(xué)機(jī)械工程學(xué)院，上海 200093）

先進(jìn)復(fù)合材料由于具備比模量和比強(qiáng)度高、耐腐蝕性能和抗疲勞性能好、設(shè)計(jì)和制造一體化成型等諸多優(yōu)點(diǎn)，而被廣泛應(yīng)用于航空、航天、交通運(yùn)輸、能源和建筑等行業(yè)[1]。當(dāng)前，先進(jìn)復(fù)合材料已由最初作為產(chǎn)品次承力部件和內(nèi)飾材料逐步發(fā)展成為主承力部件的材料，但復(fù)合材料較高的使用成本在一定層面上限制了其在工業(yè)領(lǐng)域的大規(guī)模應(yīng)用。復(fù)合材料制品的成本主要來自原材料和制造工藝兩部分，其中制造工藝方面的成本占到70%左右[2]。

作為一種典型的先進(jìn)復(fù)合材料低成本成型工藝，復(fù)合材料液體成型技術(shù)（Liquid Composites Molding，LCM）因具有生產(chǎn)效率高、制件表面光潔度好、尺寸精度高、閉模成型環(huán)境污染小的諸多優(yōu)點(diǎn)得到越來越多的關(guān)注。VARTM 工藝作為液體成型工藝中較常見的方法，近年來已被廣泛應(yīng)用于工業(yè)生產(chǎn)制造。例如，俄羅斯MC21 支線客機(jī)和龐巴迪C 系列飛機(jī)都將液體成型技術(shù)用于飛機(jī)機(jī)翼制造[3]。圖1[4]為一個(gè)典型的VARTM工藝。

傳統(tǒng)的對(duì)VARTM 工藝進(jìn)行優(yōu)化的方法往往需要依靠經(jīng)驗(yàn)并通過大量試驗(yàn)對(duì)樹脂流動(dòng)情況進(jìn)行預(yù)測(cè)。這種方法需要耗費(fèi)大量人員精力、資源和材料，拉長(zhǎng)研制周期。而采用基于流體力學(xué)方程的有限元數(shù)值模擬方法對(duì)模具內(nèi)樹脂流動(dòng)過程進(jìn)行模擬可以有效降低VARTM 工藝設(shè)計(jì)成本，有利于優(yōu)化工藝參數(shù)，從而提高產(chǎn)品開發(fā)效率。國(guó)外有很多學(xué)者對(duì)此展開研究，Drapier[5]和Blais 等[6]利用Darcy 定理耦合Stokes 方程的方法研究了樹脂在預(yù)成型體內(nèi)外的單相不飽和流動(dòng)情況，但是在數(shù)值收斂性上仍然存在很多問題。Gantois 等[7–8]采用邊界元的方法在宏觀部件尺度上模擬樹脂流動(dòng)進(jìn)行了比較，雖然計(jì)算速度較快但是只對(duì)邊界預(yù)測(cè)精度較低。近年來，國(guó)內(nèi)學(xué)者也針對(duì)RTM 工藝仿真開展了廣泛研究[9–11]，楊波[12]、Li 等[13]采用有限元方法對(duì)樹脂流動(dòng)和缺陷在微觀和細(xì)觀層面進(jìn)行了預(yù)測(cè)，但沒有探討宏觀尺度情況；戴福洪等[14]采用控制體積和等效滲透系數(shù)方法模擬邊緣效應(yīng)，以避免干斑工藝缺陷的產(chǎn)生，然而控制體積方法計(jì)算精度較差，影響了預(yù)報(bào)結(jié)果。目前，國(guó)內(nèi)采用有限元方法對(duì)樹脂單相（不含空氣相）流動(dòng)的VARTM 工藝樹脂浸潤(rùn)過程進(jìn)行的研究相對(duì)較少，在宏觀尺度上模擬的主要控制方程還是Darcy 定理，對(duì)樹脂流動(dòng)前沿輪廓一般采用水平集或控制體積等傳統(tǒng)方法。本文采用理查德方程作為流動(dòng)控制方程，應(yīng)用一種輪廓函數(shù)描述樹脂瞬時(shí)流動(dòng)前沿，基于COMSOL Multiphysics 多物理場(chǎng)仿真軟件的偏微分方程（Partial Differential Equation，PDE）模塊，通過有限元方法實(shí)現(xiàn)了對(duì)樹脂充模過程中樹脂流動(dòng)前沿和流場(chǎng)壓強(qiáng)的預(yù)報(bào)，并將該方法得到的數(shù)值仿真結(jié)果與所設(shè)計(jì)的試驗(yàn)進(jìn)行對(duì)比，證明了該方法的正確性。

1 數(shù)值模擬方法

對(duì)樹脂流動(dòng)進(jìn)行數(shù)值模擬的過程可以歸屬為計(jì)算流體力學(xué)（Computational Fluid Dynamics，CFD）的范疇，一個(gè)完整的計(jì)算流體力學(xué)過程如圖2所示。

控制方程是整個(gè)計(jì)算過程的核心，是對(duì)物理過程的公式化表述；初始值和邊界條件規(guī)定了物理過程的起始情況和邊際范圍；對(duì)計(jì)算域的網(wǎng)格劃分、方程和初始條件邊界條件的離散是有限元的基礎(chǔ)，通過循環(huán)求解以達(dá)到數(shù)值結(jié)果滿足收斂條件。

1.1 控制方程

為滿足公式的使用條件，需要進(jìn)行如下假設(shè)：

（1）纖維預(yù)成型體為連續(xù)的多孔介質(zhì)結(jié)構(gòu)，并且在樹脂浸潤(rùn)過程中不會(huì)發(fā)生形變。

（2）樹脂在纖維預(yù)成型體中的流動(dòng)為多孔介質(zhì)滲流并滿足達(dá)西流動(dòng)特性。

（3）樹脂的流動(dòng)為恒溫流動(dòng)，充模過程中樹脂黏度保持恒定，且忽略樹脂的化學(xué)反應(yīng)。

1.1.1 滲流方程Darcy定理

Darcy 定理：

其中，稱=–為 達(dá)·西ΔP速度（m/s），表示平均體速度；μ是樹脂的黏度（Pa·s）；P為壓強(qiáng)（Pa）；K代表多孔介質(zhì)纖維增強(qiáng)材料的滲透率張量（m2）。K為二階張量，對(duì)于笛卡爾坐標(biāo)系下的二維多孔介質(zhì)材料的滲透率張量可以表示為：

三維滲透率張量可表示為：

1.1.2 質(zhì)量守恒定律

圖2 計(jì)算流體力學(xué)過程Fig.2 Process of computational fluid dynamics

質(zhì)量守恒定律或稱連續(xù)性方程，是任何流體流動(dòng)過程中需要遵守的普遍性方程。其中，ρ表示流體密度；v是速度向量；ε表示多孔介質(zhì)的預(yù)成型體的孔隙率。

當(dāng)式（2）中流體流動(dòng)達(dá)到穩(wěn)態(tài)平衡時(shí)，流體不可被壓縮，此時(shí)為常數(shù)，式（2）可化簡(jiǎn)為：

1.1.3 理查德方程（Richard方程）

將經(jīng)典Darcy 方程（式（1））代入式（2），同時(shí)引入無量綱的飽和度變量S∈[0,1]，可得到式（4）即為簡(jiǎn)化的Richard 方程。

1.2 輪廓方程

理查德方程建立了瞬時(shí)條件下的飽和度S和壓強(qiáng)P之間的關(guān)系。對(duì)樹脂流動(dòng)前沿的輪廓描述，本文采用Kojima 等[15]提出的輪廓方程，如式（5）所示，它通過建立S（P）方程確立了變量壓強(qiáng)P和飽和度S的直接關(guān)系。這里α定義為數(shù)值形狀因子（1/Pa），不同α選取值對(duì)飽和度S和壓強(qiáng)P關(guān)系的影響如圖3所示。

從圖3可知，隨α升高數(shù)值模擬的收斂性逐漸升高，綜合考慮數(shù)值模擬的收斂性、飽和度S和壓強(qiáng)P計(jì)算的準(zhǔn)確性，本文選取α=0.001。將式（5）運(yùn)用鏈?zhǔn)椒▌t進(jìn)行變換可得式（6）。這樣式（4）就轉(zhuǎn)化為關(guān)于壓強(qiáng)P與時(shí)間t的方程。

1.3 網(wǎng)格劃分

圖3 考慮變化的飽和度S和壓強(qiáng)P關(guān)系Fig.3 Relation of S and P considering variation of α

合理的網(wǎng)格劃分是后續(xù)采用有限元方法高效計(jì)算的基礎(chǔ)，在COMSOL Multiphysic 軟件中，PDE 方程采用Galerkin 方法[16]進(jìn)行離散。在采用CFD 進(jìn)行分析和邊界層網(wǎng)格劃分過程中，將控制方程離散到基函數(shù)上時(shí)，三角形單元和四邊形單元是最常用的兩種線性網(wǎng)格單元類型，如圖4所示。

為驗(yàn)證本文方法的正確性，建立了一個(gè)幾何結(jié)構(gòu)模型（圖5），在矩形區(qū)域內(nèi)預(yù)設(shè)了兩個(gè)小矩形區(qū)域阻礙樹脂流動(dòng)，在臨近下端的位置預(yù)設(shè)阻礙區(qū)域與整個(gè)矩形邊緣區(qū)域相聯(lián)通。采用三角形單元對(duì)幾何模型進(jìn)行有限元?jiǎng)澐趾蟮慕Y(jié)果如圖5（b）所示，共包含5760 個(gè)域單元和270 個(gè)邊界單元。

圖5 幾何結(jié)構(gòu)模型Fig.5 Geometrical model

1.4 邊界條件和初始值

如圖5所示，除了兩個(gè)位于幾何模型中部的注入口，其他邊界條件均設(shè)為不可滲透。

模擬開始前，模具內(nèi)壓強(qiáng)與出口壓強(qiáng)都為0；樹脂開始注入時(shí)，入口壓強(qiáng)(Pin)和出口壓強(qiáng)(Pout)的初始值分別為1×105Pa 和0。

2 試驗(yàn)設(shè)計(jì)和參數(shù)設(shè)定

試驗(yàn)設(shè)備由樹脂灌注系統(tǒng)和位于樹脂灌注系統(tǒng)上部的檢測(cè)系統(tǒng)兩部分構(gòu)成。樹脂灌注系統(tǒng)由密封鋪放在柔性真空袋內(nèi)的單層編織纖維布、導(dǎo)管和真空泵組成；檢測(cè)系統(tǒng)由CCD（Camera Coupled Device）相機(jī)構(gòu)成，負(fù)責(zé)記錄試驗(yàn)中不同時(shí)刻樹脂流動(dòng)前沿的位置，試驗(yàn)裝置如圖6所示。

試驗(yàn)和數(shù)值模擬選用的材料參數(shù)值如表1所示，纖維預(yù)成型體的滲透率沿水平方向（Kh）和豎直方向（Kv）各異。

3 結(jié)果對(duì)比驗(yàn)證

圖6 試驗(yàn)設(shè)計(jì)Fig.6 Experiment design

表1 材料參數(shù)Table 1 Material parameters

在不同時(shí)刻，樹脂流動(dòng)前沿位置和壓力場(chǎng)的試驗(yàn)結(jié)果和數(shù)值模擬結(jié)果如圖7所示，其中圖7（a）、（e）、（i）、（m）為不同時(shí)刻試驗(yàn)所得樹脂流動(dòng)前沿的記錄；圖7（c）、（g）、（k）、（o）為文獻(xiàn)[8]中采用邊界元法獲得的壓強(qiáng)場(chǎng)結(jié)果；圖7（b）、（f）、（j）、（n）和圖7（d）、（h）、（l）、（p）分別為采用本文方法所獲得的樹脂流動(dòng)前沿和壓強(qiáng)場(chǎng)的結(jié)果。

為了更直觀準(zhǔn)確地量化對(duì)比試驗(yàn)和數(shù)值模擬得到的樹脂流動(dòng)前沿結(jié)果，本文定義浸潤(rùn)比（λ），λ=S樹脂/S纖維，其中S樹脂代表對(duì)應(yīng)時(shí)刻下平面內(nèi)樹脂流過區(qū)域面積；S纖維代表整個(gè)纖維鋪放區(qū)域的面積。不同時(shí)刻模擬值和試驗(yàn)值獲得的浸潤(rùn)比（λ）如圖8所示。試驗(yàn)和數(shù)值模擬的差異可以理解為，數(shù)值模擬中只考慮樹脂的單相流體流動(dòng)，不包含空氣氣體相對(duì)流動(dòng)的阻礙。而在實(shí)際試驗(yàn)中樹脂灌注初始階段，由于真空袋密封區(qū)域內(nèi)還有部分空氣未排出，且殘余空氣沒有均勻分布在纖維預(yù)成型體內(nèi)，所以樹脂流動(dòng)初始階段試驗(yàn)獲得的浸潤(rùn)比λ試驗(yàn)<λ模擬；而伴隨樹脂灌注繼續(xù)進(jìn)行，空氣逐漸被排出密封區(qū)域，樹脂流動(dòng)加快，則λ試驗(yàn)>λ模擬。

不同時(shí)刻數(shù)值模擬與試驗(yàn)的λ絕對(duì)值差，即|λ試驗(yàn)–λ模擬|小于10%，具有良好的準(zhǔn)確性。將不同時(shí)刻壓強(qiáng)場(chǎng)結(jié)果的極值與文獻(xiàn)中結(jié)果進(jìn)行了對(duì)比，亦得到近似結(jié)果。本文方法在與文獻(xiàn)[7]中所列相似的硬件平臺(tái)上運(yùn)行相同案例和相近計(jì)算規(guī)模時(shí)，計(jì)算時(shí)間僅為57s，遠(yuǎn)低于該文獻(xiàn)中采用邊界元法計(jì)算的200s，這證明了本方法在計(jì)算效率上的優(yōu)越性。因此，可以得出本文所用數(shù)值模擬方法對(duì)壓強(qiáng)和數(shù)值流動(dòng)前沿的預(yù)報(bào)都較為準(zhǔn)確，證明了該方法的正確性。

4 結(jié)論

本文基于理查德方程和一種輪廓函數(shù)，利用COMSOL Multiphysics 多物理場(chǎng)仿真軟件，實(shí)現(xiàn)了對(duì)樹脂浸潤(rùn)纖維預(yù)成型體過程的數(shù)值模擬預(yù)測(cè)。預(yù)測(cè)結(jié)果經(jīng)與試驗(yàn)結(jié)果相比較，證明該方法具有準(zhǔn)確的預(yù)報(bào)結(jié)果。

本文所實(shí)現(xiàn)的方法，目前僅用于二維模型的流動(dòng)仿真，由于實(shí)際生產(chǎn)制造中許多復(fù)材部件具有一定的厚度尺寸，因此在后續(xù)研究中可以探索將該方法擴(kuò)展到三維尺度模型上。其次，本文中所用的纖維增強(qiáng)材料的滲透率數(shù)值為材料本身標(biāo)定，在后續(xù)研究中可以探究基于本文建立的宏觀尺度方法，將基于數(shù)值模擬對(duì)細(xì)觀結(jié)構(gòu)進(jìn)行分析獲得的滲透率數(shù)值迭代入本文方法中，進(jìn)行多尺度滲流研究。同時(shí)，溫度變化對(duì)流體流動(dòng)特性的影響諸如樹脂黏度變化等，也可以納入后續(xù)數(shù)值模擬的研究范圍。

圖7 驗(yàn)證結(jié)果對(duì)比Fig.7 Comparison of results

圖8 不同時(shí)刻浸潤(rùn)比模擬值和試驗(yàn)值Fig.8 Infusion ratio λ of results from simulation and experiments at variant time