非等溫樹脂傳遞模塑成型仿真建模及應用

趙 亮, 高勝暉*, 段躍新

（1.中國航發(fā)北京航空材料研究院，北京 100095；2.北京航空航天大學，北京 100191）

樹脂基復合材料有著良好的力學性能、質量輕、耐腐蝕等特點，在航空、交通、電子電力等領域應用廣泛，其中一種高效的成型技術是樹脂傳遞模塑成型（ resin transfer molding， RTM） 。 采用RTM 工藝模擬仿真，可以有效的預測樹脂填充時間，優(yōu)化注膠口位置和溢料口位置，避免干斑缺陷，從而改善傳統(tǒng)試錯的工藝驗證方式，有效縮短工藝設計時間，節(jié)約成本。在實際生產(chǎn)過程中，非等溫RTM 過程是更為接近真實生產(chǎn)過程的物理描述。相比等溫RTM 的仿真建模，非等溫RTM 需要同時考慮溫度和流場的耦合計算，對算法的可靠性和收斂性均提出了較高的要求。國外商業(yè)化軟件如法國的PAM-RTM、比利時的RTM-Worx 等均實現(xiàn)了非等溫RTM 工藝的仿真求解，在國內企業(yè)得到廣泛應用。長期以來國內針對RTM 仿真求解器的研發(fā)主要以等溫求解為主，針對非等溫過程的熱流耦合求解仍然處于基礎研究階段，目前國內尚缺少成熟的針對工程應用的非等溫RTM 工藝仿真模型，在實際生產(chǎn)過程中長期依賴國外商業(yè)化模型。隨著我國對基礎工業(yè)仿真重視程度的進一步提升，工藝制造仿真作為基礎工業(yè)仿真的組成部分，也受到廣泛的關注和重視。

RTM 工藝的充模過程是一種典型的非定常流動現(xiàn)象，對于該過程的模擬仿真國內外學者分別采用有限元、有限元控制體積法、有限差分法等進行了研究。有限元方法可以用來模擬樹脂充模的工藝過程,但傳統(tǒng)有限元方法固有的特點使得求解過程較繁瑣，內存消耗較大，計算速度較慢，且在追蹤流動前鋒上存在一定不足[1]。達西定律適用于多孔介質中低速穩(wěn)態(tài)流動的描述，當采用達西定律描述樹脂在纖維中的流動時，需要假定在短時間內樹脂的流動是一種穩(wěn)態(tài)流動，然后可采用有限元方法求解，最終得到樹脂流動的速度和壓力分布[2-4]。

控制體/有限元法綜合了有限體積方法（FVM）和有限元方法（FEM）的優(yōu)點，逐漸成為主流的RTM 仿真工藝數(shù)值算法。Bruschke 等[5]采用CV/FEM 方法模擬了樹脂在滲透率各異的纖維鋪層中的等溫充模過程，經(jīng)過實驗驗證，在充填時間和流動前鋒方面數(shù)值結果與實驗結果吻合較好。Bruschke 等較早考慮了RTM 工藝中的非等溫的現(xiàn)象，采用CV/FEM 模擬了熱傳導和固化反應等對工藝過程的影響。區(qū)別于傳統(tǒng)的控制體/有限元方法，Trochu 等[6]和趙亮等[7]采用了類CV/FEM方法，將網(wǎng)格單元看作一個控制單元，將壓力場結果存儲到單元的重心位置，通過有限元方法得到了充模過程中的壓力值，通過比較，與實驗結果吻合較好。北京航空材料研究院2000 年與北京航空航天大學梁志勇等[8]采用CV/FEM 方法合作開發(fā)了具備自主知識產(chǎn)權的二維RTM 工藝仿真軟件BHRTM，該軟件可在二維空間迅速簡便實現(xiàn)RTM 工藝幾何創(chuàng)建和網(wǎng)格的自動劃分，在后處理結果展示中可對樹脂流動和壓力場演化過程進行偽三維顯示。尹明仁等[9]基于CV/FEM 方法開發(fā)了模擬軟件平臺BHRTM-2，該平臺可實現(xiàn)充模過程的可視化顯示，對工藝設計中可能出現(xiàn)的干斑等問題起到了較好的指導作用。施飛[10]基于CV/FEM 方法實現(xiàn)了在非結構四面體網(wǎng)格中RTM 工藝的數(shù)值模擬，包括等溫和非等溫工藝過程，模擬結果與實驗結果吻合較好。

有限差分法對控制方程直接差分化，程序實現(xiàn)相對較為簡單，最有代表性的方法是貼體坐標法[11-12]，該方法首先需要對網(wǎng)格進行雅克比空間變換，得到規(guī)則形狀的計算域,隨后采用有限差分法求解。Coulter 等[13]根據(jù)達西定律,采用有限差分法模擬了樹脂在充模過程的非等溫流動,考慮了充模過程中樹脂因熱傳遞、固化反應等引起的溫度變化最終對充模過程的影響。Friedrichs 等[14]和施飛等[15]采用有限差分方法中的貼體坐標法模擬了樹脂流過形接頭的情形,考慮了不同部件間網(wǎng)格并不一定完全一致的情形，最終得到了不同時刻的樹脂流動前沿曲線以及最終的壓力分布情況。李海晨等[16]采用有限差分算法模擬了樹脂傳遞模塑成型工藝的等溫過程，得到了充模過程中前鋒形狀的變化和不同時刻下壓力場的分布情況，仿真結果與其他程序仿真結果吻合較好。

控制體有限元方法相比于傳統(tǒng)的有限元方法和有限體積方法具有編程簡單的優(yōu)勢，同時相對于有限差分和邊界元等方法具有更好的守恒特性，適用于RTM 工藝中樹脂流動過程的仿真建模。從仿真經(jīng)驗來講，在相同數(shù)值精度下，六面體網(wǎng)格相比于四面體網(wǎng)格計算結果要更精確。然而，在RTM工藝仿真中，將控制體有限元方法應用于六面體單元的仿真模型并不多見，本工作將基于六面體網(wǎng)格單元構建RTM 工藝的數(shù)值模型。

1 模型

1.1 物理模型

在真實的 RTM 工藝中，樹脂在流動過程中與纖維之間不斷發(fā)生熱交換，使得樹脂本身的溫度、黏度、固化度不斷變化。固化反應的效應又反過來影響樹脂的溫度和黏度，各因素之間相互影響，最終對樹脂的流動模式產(chǎn)生改變，建立非等溫模型是目前數(shù)值刻畫這些因素的主要方法。

（1） Darcy 方程

在樹脂傳遞模塑成型工藝中，由于樹脂充填速度一般很小，流動雷諾數(shù)較低，樹脂在多孔介質中的流動服從Darcy 定律，當假定樹脂黏度在短時間內不發(fā)生變化時，樹脂流動速度與壓力梯度成正比，Darcy 方程為：

式中：V是 Darcy 速率；[K] 是纖維增強體滲透率的二階張量；μ是樹脂黏度； ?P是壓力梯度。

在非等溫 RTM 模型中，樹脂的流動與等溫情況有所差異，卻依然遵循 Darcy 定律，但樹脂黏度μ不再是一個常數(shù)，而是隨著溫度和固化度變化不斷變化的量?？坍嫎渲髯冃阅艿哪Ｐ陀泻芏啵缫浑A等溫模型、凝膠模型、雙阿累尼烏斯黏度模型等。

（2）連續(xù)性方程

連續(xù)性方程是流動的守恒方程，在流體微小體積單元內，質量隨時間的變化率應等于在這段時間內流入該微小體積的凈質量

式中：ρ是樹脂密度；t是時間；?是纖維的孔隙率。由于樹脂流動速度較低，黏度較大，在短時間內可假定樹脂的密度不發(fā)生變化，即 ?ρ/?t =0。連續(xù)性方程可簡化為：

數(shù)值求解中，速度和壓力的邊界條件為：（1）在壁面處，速度為無穿透邊界條件，壓力取零法向梯度邊界條件；（2） 注膠口處，可根據(jù)輸入的工藝條件選擇恒壓注射和恒流注射，若采用恒壓注射，p = p0，p0為注射壓力；如果采用恒流注射，V = V0，V0是注射速度。

（3） 能量守恒方程

在 RTM 數(shù)值仿真的非等溫模型中，溫度隨時間和空間不斷變化，變化的溫度影響了樹脂的黏度和固化度，而固化反應的發(fā)生又反過來影響了樹脂的溫度和黏度，黏度的變化通過影響流動速度也必然引起溫度和固化度的變化，這三個量的綜合作用下最終影響了樹脂的流動模式。在此過程中，影響因素眾多，需建立相應的模型和方程來定量刻畫諸因素的影響，首先考慮關于溫度的能量方程。

在模具內，當假定纖維溫度與樹脂的溫度相同時，即Tf= Tr= T，該假設忽略了樹脂與纖維之間的熱傳遞，關于樹脂的能量方程和纖維的能量方程可不必分開求解，而是共同遵循平衡模型下的能量方程：

式中：等式左邊第一項是樹脂和纖維增強體溫度隨時間的變化率，偏導數(shù)前的系數(shù)是樹脂和纖維的綜合熱效應；第二項是樹脂隨充模過程引起的熱對流效應；等式右邊第一項是熱擴散項，表征由不同方向上溫度的不均勻性引起的熱傳導，其中krf代表xyz三個方向上的熱傳導系數(shù)；右邊最后一項是固化熱源項，表征因固化反應放熱帶來的熱量的產(chǎn)生。

（4） 固化動力學方程

在 RTM 成型工藝中，當達到一定溫度時，樹脂分子間便會發(fā)生聚合，即樹脂的固化反應。固化反應使得樹脂流動性變差，黏度增加，但固化反應中釋放的熱量又促使樹脂黏度變小。樹脂固化反應的控制方程為：

式中：α為樹脂固化度，取值范圍為0～1；f（α，Tr）為樹脂的固化反應速率，由固化反應的模型刻畫，在數(shù)值求解時，采用了與能量方程中固化熱源項類似的處理方式，對該項進行了線化處理。

（5） 流變模型

在RTM 成型工藝中，由于溫度在空間和時間尺度上不斷變化，樹脂的黏度也在不斷發(fā)生變化，黏度的變化最終對填充結果有著較大的改變。樹脂黏度模型主要包括理論模型、經(jīng)驗模型和半經(jīng)驗模型，在數(shù)值求解中引入黏度模型可用來定量刻畫黏度隨溫度時間等的變化。黏度模型也常常表述為流變模型，其中經(jīng)驗/半經(jīng)驗模型因其經(jīng)濟實用性在 RTM 工藝具有更廣泛的應用。經(jīng)驗模型包括工程黏度模型、雙阿累尼烏斯黏度模型、（williamslandel-ferry， WLF）方程等[17-18]。其中雙阿累尼烏斯模型方程是適用性較廣的經(jīng)驗流變模型，該模型是在研究環(huán)氧樹脂體系黏度變化規(guī)律的基礎上提出的半經(jīng)驗公式，需假設樹脂體系的固化反應為一級反應或某些總體為非一級動力學反應的樹脂體系；工程黏度模型是另一種適用較廣的經(jīng)驗方程，其基本原理是樹脂體系物理和化學反應相互作用最終導致了樹脂固化過程中黏度的變化；WLF 方程適用于研究熱固性樹脂體系的黏度變化。

1.2 計算模型

控制體/有限元方法并不限定網(wǎng)格類型，二維中的多邊形，三維中的多面體也均適用。如圖1 所示，以規(guī)則的四邊形單元為例，控制體是由包含該節(jié)點的四邊形單元各邊的中點與四邊形重心的連線圍成，每個控制體歸屬中心節(jié)點所有，數(shù)值計算的結果也均存儲于節(jié)點處。如圖1 所示，Ni,j是陰影構成的控制的中心節(jié)點，Q1～Q8 是構成該控制體的小面。在RTM 工藝中，樹脂的流動遵循物質守恒定律和達西定律，根據(jù)這兩個方程可求得關于壓力的拉普拉斯方程。CV/FEM 方法將壓力定義在每個節(jié)點上，在每個小控制體上對壓力方程進行離散可得到一些列的代數(shù)方程，求解后便可得到計算域內壓力的分布，再結合達西定律，就可求得樹脂的流動情況即樹脂速度場的分布情況。

圖1 CV/FEM 單元示意圖Fig. 1 Diagram of the CV/FEM cell

可采用控制體填充系數(shù)f來描述樹脂流動的前鋒。充填系數(shù)f代表了樹脂填充控制體的情況，當f=1 時，樹脂已完全充滿控制體；當f<1 時，樹脂未充滿控制體；當f=0 時，樹脂尚未流動到該控制體，屬于未充填區(qū)域。在流動的前鋒，充填系數(shù)f介于0 和1 之間，統(tǒng)計所有滿足0

在算法設計中，為保證數(shù)值解在時間和空間方向互不影響，在能量方程求解中采用了半離散方法；為減小因分離式求解多個方程帶來的數(shù)值誤差，采用了一種三收斂標準[10]。下面對半離散方法和三收斂標準做簡要概述。

（1） 半離散方法

能量方程的求解采用了半離散的控制體/有限元方法，半離散方法將空間偏導數(shù)項的離散和時間方向的推進完全分開，可使空間的誤差和精度等與時間推進的穩(wěn)定性、收斂的加速性等無關。經(jīng)過空間離散，能量方程就轉化為了如下形式：

式中：P(T)為空間離散后關于不同控制體上溫度T的線性組合，方程由原來的偏微分方程轉化成了形式上的常微分方程，在時間方向推進即可求得方程的解。本工作在能量方程的時間方向采用了四階Runge-Kutta（龍格-庫塔）方法，該方法可以構造高精度來求解初值問題，至今仍為實際應用的重要方法。四階 Runge-Kutta 方法的具體計算格式如下：

（2） 三收斂標準

在與非等溫相關的方程如能量方程、固化動力學方程的求解中，采用了分離式求解算法。雖然方程的求解得到了較大的簡化，但也容易因為更新的滯后造成誤差的累計從而引起數(shù)值解的不準確。基于以上考慮，本工作采用了一種三收斂的標準，在算法中設計了一個內循環(huán)，具體表述為：在一個時間步內，壓力和速度計算完成后，就進入到非等溫場的計算內循環(huán)中，只有當溫度場、黏度場、固化度場均循環(huán)收斂時，才會跳出該循環(huán)進入下一個時間步。

1.3 驗證模型

仿真模型采用C++語言進行了求解器的開發(fā)，基于QT 及python 語言進行了前后處理器的開發(fā)，仿真模型命名為BIAM Composites-RTM（簡稱BIAM-RTM）。仿真模型可模擬各種復雜的構件和多種工藝工況，最終得到壓力場、速度場、溫度場等仿真結果。驗證模型選取了平板實驗件和工程大型零部件機匣，現(xiàn)對兩個驗證模型如下說明。

北京航空航天大學段躍新團隊對復合材料平板實驗件進行了工藝實驗。平板實驗件幾何模型如圖2（a）所示，其尺寸為32 cm×32 cm，樹脂黏度為0.062 Pa?s，玻璃布滲透率為9.7092×10?10m2，纖維含量為29%，模腔厚度5 mm，底邊線恒壓注射，注射壓力為0.02 MPa，實驗充填時間為114 s。

圖2 驗證模型的幾何 （a）平板實驗件；（b）機匣件Fig. 2 Geometry of verification models （a）flat test piece；（b）receiver piece

隨后，選取航空發(fā)動機零件中典型機匣件進行了應用驗證。機匣復合材料實驗件幾何模型如下圖2（b）所示，其半徑尺寸為0.42 m，高0.25 m，樹脂黏度為0.19 Pa?s，預制件滲透率為徑向4.94×10?13m2，纖維體積含量為52%，底邊線恒壓注射，注射壓力為1.2 MPa，實驗充填時間為7200 s。

2 結果與討論

2.1 平板實驗件應用驗證

應用BIAM-RTM 對平板件進行了數(shù)值模擬，工藝參數(shù)與實驗條件保持一致。BIAM-RTM 仿真結果如圖3 所示。圖中分別給出了四個不同時刻（5.76、19.19、53.73、105.56 s）的填充度f和對應的壓力分布情況。底邊線同時進膠，前鋒點位于同一高度上，隨著不斷注膠，前鋒點向上平直推進，在同一高度上的壓力值也基本平直。由于平板件實驗中無法得到某點處的速度和壓力值隨時間變化曲線，為了驗證仿真結果的可靠性，本工作采用PAMRTM 對平板件進行了相同工況的數(shù)值模擬。

圖3 平板件填充度f 和壓力P 隨時間的變化 （a），（e）t=5.76 s；（b），（f）t=19.19 s；（c），（g）t=53.73 s；（d），（h） t=105.56 sFig. 3 Variations in fraction and pressure of flat piece over time （a），（e）t=5.76 s；（b），（f）t=19.19 s；（c），（g）t=53.73 s；（d），（h） t=105.56 s

PAM-RTM 在RTM 工藝模擬方面非常強大，具有較高可靠性而受到業(yè)內的廣泛認可和好評。PAMRTM 充填時間為115.586 s，與實驗結果114 s 接近。表1 給出了兩款軟件與實驗充填時間的對比誤差。

表1 平板件BIAM-RTM 與PAM-RTM 充填時間與實驗對比誤差Table 1 Relative error of filling time compared to the experimental result of BIAM-RTM and PAM-RTM on flat piece

取平板件中心點的結果與PAM-RTM 仿真結果進行比對，速度及壓力的對比曲線如圖4 所示，兩款軟件仿真結果較為接近，趨勢吻合較好。取峰值點進行定量比對，BIAM-RTM 的速度在峰值為0.00282 m/s，PAM-RTM 為0.00273 m/s，誤差為（ 0.00282-0.00273） /0.00273=3.30%。 BIAM-RTM的壓力在峰值為0.01090 MPa ，PAM-RTM 為0.01122 MPa，誤差為（0.01122-0.0109）/0.01122=2.85%。

圖4 平板實驗件中心點處兩款軟件速度和壓力隨時間變化曲線 （a）速度；（b）壓力Fig. 4 Curves of velocity and pressure versus time between two softwares for the center point of flat test piece （a）velocity；（b）pressure

由以上分析可得出結論，BIAM-RTM 在平板件的模擬對比中，充填時間與實驗結果和PAMRTM 仿真結果基本相當；中心點處的速度和壓力與PAM-RTM 結果吻合較好，速度峰值點和壓力峰值點與PAM-RTM 結果較為接近。

2.2 機匣件應用驗證

采用BIAM-RTM 對機匣件進行了模擬，模擬結果如圖5 所示，圖中分別給出了四個不同時刻（58.64、1257.38、2255.89、5846.90 s）下填充系數(shù)f和壓力的分布情況。底邊線同時進膠，在相同高度下填充系數(shù)和壓力值近乎相等。BIAM-RTM 模擬充填時間為7076 s，與實驗結果7200 s 較為接近。同樣采用PAM-RTM 在同樣工藝條件下對機匣件進行了數(shù)值仿真，PAM-RTM 得到的充填時間為7206 s。表2 給出了兩款軟件模擬充填時間與實驗值的對比誤差，誤差在2%范圍內。

表2 機匣件BIAM-RTM 與PAM-RTM 充填時間與實驗對比誤差Table 2 Comparison error of filling time between BIAMRTM and PAM-RTM

圖5 機匣件填充度f 和壓力P 隨時間的變化 （a）填充度；（b）壓力；（1）t=58.64 s；（2）t=1257.38 s；（3）t=2255.89 s；（4）t=5846.90 sFig. 5 Variations in fraction and pressure of receiver piece over time （a）f；（b）p；（1） t =58.64 s；（2） t =1257.38 s；（3） t =2255.89 s；（4）t =5846.90 s

在同一高度上的控制體具有相同的模擬結果，取機匣軸向中間位置任意點結果進行比對，速度及壓力的對比曲線如圖6 所示。觀察圖6，從整體來講，BIAM-RTM 與PAM-RTM 模擬結果吻合較好，壓力和速度在上升段的斜率不同，但在工藝末期，速度和壓力值基本吻合。取峰值點進行比對，BIAM-RTM 的壓力在峰值為533230 Pa ，PAMRTM 軟件為503812 Pa，誤差為（533230?503812）/503812=5.84%； BIAM-RTM 的速度在峰值為0.04094 mm/s，PAM-RTM 為0.0393 mm/s，誤差為（0.04094-0.0393）/0.0393=4.17%。

圖6 機匣件中間點兩款軟件速度和壓力隨時間變化曲線 （a）速度；（b）壓力Fig. 6 Curves of velocity and pressure versus time between two softwares for the midpoint of the receiver piece （a）velocity；（b）pressure

由以上分析可得出結論，針對機匣件的模擬結果，自研軟件BIAM-RTM 與對標軟件PAM-RTM在充填時間、速度及壓力的仿真結果吻合較好，相對誤差在6% 范圍以內。

3 結 論

（1）BIAM-RTM 分別在復合材料平板實驗件和大型機匣結構件上進行了模擬仿真和應用驗證。在充填時間、充填速度、壓力分布等關鍵工藝指標的仿真結果上與國外商業(yè)模型精度一致，同時仿真結果與實驗結果吻合較好，表明了仿真模型與算法的合理性。

（2）大型機匣類環(huán)形結構件的仿真結果與PAM-RTM 在充填時間、速度及壓力等方面誤差在6%以內，這表明BIAM-RTM 基本具備了針對實際復雜結構零件開展工程仿真的能力。