預成形體滲透率預測及剪切變形的影響

楊 波，金天國，畢鳳陽，2，李建廣

（1哈爾濱工業(yè)大學 機電工程學院，哈爾濱150001；2黑龍江工程學院 機電工程學院，哈爾濱150050）

復合材料液體模塑成型工藝（Liquid Composite Molding，LCM）如今正廣泛應用于制造復雜結構的復合材料構件，然而以試驗和經(jīng)驗為主要手段的模具設計、工藝參數(shù)的確定方法成本高、效率低，很大程度上推升了該工藝制件的價格，限制了LCM工藝的應用推廣。因此采用數(shù)值模擬代替實驗來預測模具填充過程中可能產(chǎn)生的缺陷，進而優(yōu)化模具設計和工藝參數(shù)的方法越來越受到重視。

樹脂在預成型體中的流動可以看做是不可壓縮黏性流體在多孔介質(zhì)中的滲流問題，通常采用達西定律描述，聯(lián)立求解達西定律與連續(xù)性方程即可獲得流場參數(shù)，進而預測工藝過程中可能產(chǎn)生的缺陷。在此過程中，滲透率是進行填充仿真過程中最重要的輸入?yún)?shù)，對仿真結果精度有極大的影響。

滲透率的研究是為樹脂充模仿真服務的，然而僅僅研究單胞在經(jīng)緯紗線正交狀態(tài)下的滲透率是不足的，因為織物鋪覆在模具上時會在約束下產(chǎn)生變形，這樣滲透率會發(fā)生較大變化?？椢镤伕策^程中可能發(fā)生的變形有剪切、矯直、起皺、拉伸和滑移等，其中剪切是織物變形的主要形式。通常忽略其他變形模式，認為織物鋪覆在細觀上的體現(xiàn)就是單胞剪切變形［1，2］。以某雷達天線屏蔽器的LCM工藝模具填充仿真為例，其仿真過程如圖1所示。首先，建立模具型腔的幾何模型，并進行離散；然后利用鋪覆算法計算出織物在模具內(nèi)受力變形導致的局部剪切角，根據(jù)滲透率與剪切角的關系獲得預成形體的滲透率分布；最后根據(jù)滲透率分布求解達西定律獲得填充過程參數(shù)。可以看出，研究預成形體在剪切變形后的滲透率變化有重要意義。

圖1 填充仿真過程Fig.1 Illustration of filling simulation for LCM

目前主要通過實驗測定滲透率，然而實驗成本高、時間長、對環(huán)境敏感，尤其是剪切滲透率很難準確測定。因此，通過數(shù)值方法模擬流體流過預成形體單胞，獲得流場參數(shù)后回帶入達西定律預測復雜預成型體滲透率的方法正逐漸顯出其優(yōu)勢。

Verleye等［3，4］針對無皺褶織物和平紋織物建立了細觀幾何模型，采用有限差分法求解Stokes方程預測單胞滲透率，并且研究了剪切的影響，最后通過實驗進行了驗證。Naoki Takano等［5］采用漸進均勻理論研究了織物細觀和微觀滲透率，采用有限元法求解流動方程，并對單胞剪切滲透率進行了深入探討。Chen等［6］采用FLOTRAN CFD軟件仿真了平紋織物單胞中的細觀流動，預測了細觀滲透率，并通過大量仿真研究了幾何參數(shù)對滲透率的影響。國內(nèi)研究者目前大都采用實驗法測定滲透率［7－10］，對數(shù)值模擬方法預測滲透率的研究很少。其中，戴福洪等［11］采用均勻化法預測了單胞滲透率，并通過與文獻中的實驗結果比較，驗證了其方法的合理性。倪愛清等［12］采用均勻化理論分析了流體在多孔介質(zhì)中的流動問題，推導出廣義達西定律，分別對單向纖維織物和二維平面織物的滲透率進行了驗證，考察了單胞微觀結構對滲透率的影響。吳炎等［13］采用有限元法針對雙向縫合氈單胞進行了流動分析，預測了滲透率，并通過一組實驗進行了驗證。

目前的研究主要有以下局限：（1）計算大都將紗線視為不可滲透的，忽略了紗線內(nèi)部的樹脂流動，這樣的預測結果不能體現(xiàn)預成形體單胞的雙尺度特性；（2）受幾何建模或計算方法的限制，預成形體單胞的幾何模型被大大簡化：紗線截面簡化為矩形，紗線的卷曲路徑也采用分段直線來代替。這樣的簡化損害了滲透率預測的精度；（3）研究大都針對正交單胞，對剪切變形單胞的滲透率求解方法研究較少，滲透率預測結果對實際應用的幫助有限。本工作考慮紗線的內(nèi)部流動建立樹脂流動模型，建立逼近真實的預成形體正交單胞幾何模型；基于Chorin投影法和Adams－Bashforth差分格式構造求解樹脂流動控制方程的高分辨率TVD數(shù)值格式，預報預成形體滲透率；在此基礎上，針對剪切單胞，采用貼體坐標法實現(xiàn)其滲透率的預測，研究單胞剪切變形對其滲透率的影響。

1 流動數(shù)學建模

根據(jù)Darcy定律，預成形體滲透率可以用下式表示：

其中，K表示多孔介質(zhì)的滲透率張量；u＝（u，v，w）表示樹脂流動速度，上劃線表示體積平均；μ表示樹脂的黏度；p表示壓力。流場速度和壓力分布為未知量，需要通過數(shù)值求解預成形體中流場控制方程獲得。

宏觀預成形體由紗線編織成纖維布，再通過鋪覆過程獲得，結構復雜度非常高，不可能一次性仿真構件預成形體內(nèi)的流動，但是根據(jù)其周期性可知，預成形體是由大量重復排列的編織結構成，稱之為細觀單胞，以下簡稱單胞，圖2所示為平紋織物單胞?？紤]紗線的可滲透性，單胞內(nèi)部的樹脂流動由兩個尺度的流動組成，一是紗線間區(qū)域的樹脂流動，稱之為細觀流動；二是樹脂在紗線內(nèi)部纖維單絲間的流動，稱之為微觀流動。

圖2 平紋織物單胞內(nèi)的樹脂流動Fig.2 The resin flow of plain fabric unit－cell

將樹脂看做不可壓縮黏性牛頓流體，填充過程視為層流、恒溫的流動過程，這樣，紗線間的樹脂細觀流動可以用不可壓縮Navier－Stokes方程表示：

其中，f表示液體的單位質(zhì)量力，ρ表示流體密度，Δ表示Laplacian算子。

由于樹脂流動速度很慢，忽略非線性對流項，另外，樹脂的密度對滲透率的預測沒有影響，因此樹脂在紗線間的細觀流動可以用下式表達：

對于紗線內(nèi)部纖維間的流動，可以將紗線看做多孔介質(zhì)，采用Brinkman方程描述流動：

其中，Kyarn表示紗線內(nèi)部的微觀滲透率張量，用來表達紗線內(nèi)部纖維間孔隙的滲透性能，受到纖維分布狀態(tài)的影響，Kyarn通常通過經(jīng)驗公式（如Gebart模型［14］）獲得。對比式（3）和式（4）可以發(fā)現(xiàn)，Brinkman方程僅細觀流動控制方程多一個余項μu／Kyarn，因此，可以將單胞內(nèi)的整個流動全部采用Brinkman方程來描述，紗線外部的細觀流動可看做是K－1yarn＝0的Brinkman流動。

聯(lián)立求解式（4）和連續(xù)性方程，即可獲得單胞內(nèi)流場的壓力分布和速度分布，回帶入式（1），即可求得織物單胞的滲透率值K。

2 正交單胞滲透率預測

2.1 幾何建模及邊界條件

預成形體單胞幾何模型通常采用通用CAD（Computer Aided Design，計算機輔助設計）軟件（如Pro／E、CATIA等）建立，受其限制，很難建立逼真的單胞幾何模型，因此需要將紗線截面近似為多邊形，紗線中心線路經(jīng)簡化為分段直線，這樣的簡化模型與單胞的真實結構差距很大，雖然方便了建模和計算，但是損害了滲透率預報精度。本工作利用ACIS幾何建模內(nèi)核，采用曲線曲面建模方法建立該逼真的織物正交單胞幾何模型如圖3所示。

圖3 平紋織物正交單胞幾何模型Fig.3 Geometry model of plain fabric orthogonal unit－cell

預測單胞細觀X向（Y向類似）滲透率時的求解區(qū)域及邊界條件設置如圖4所示，在入口Γ3和出口Γ6施加恒壓邊界條件分別為P0和0；由于預成形體是由多層織物堆疊而成，而且單胞具有周期性，在其余邊界施加周期性邊界條件，即Γ1和Γ4，Γ2和Γ5分別成對施加周期邊界條件。

圖4 求解區(qū)域及邊界條件設置Fig.4 Calculation region and boundary conditions

基于有限差分法求解單胞中樹脂流動控制方程。對流動區(qū)域幾何模型采用矩形網(wǎng)格離散，為了避免計算產(chǎn)生沒有物理意義的鋸齒形流動速度或壓力分布，本工作采用交錯網(wǎng)格布置方式。以二維網(wǎng)格布置為例，交錯網(wǎng)格布置方式如圖5所示。壓力在由實心點表示的網(wǎng)格點上計算，水平和垂直流動速度分別在由灰色實心點和空心點表示的網(wǎng)格點上計算。這樣壓力和速度就在不同的網(wǎng)格點上計算，可以消除產(chǎn)生鋸齒形流速和壓力分布的可能性。三維交錯網(wǎng)格與二維相似，不再贅述。

圖5 交錯網(wǎng)格布置Fig.5 The staggered grid

2.2 細觀流動數(shù)值求解

由于本工作的流動計算考慮了紗線的可滲透性，而且紗線的滲透率較低，紗線內(nèi)部的流動阻力會在紗線和孔隙界面處造成明顯的速度不連續(xù)現(xiàn)象，這樣在采用普通差分格式計算流動時，可能會呈現(xiàn)出震蕩現(xiàn)象，影響計算的穩(wěn)定性，這種震蕩完全是因為數(shù)值的原因。為了避免上述現(xiàn)象的發(fā)生，應該選擇一種滿足總變差減?。═VD）條件的高分辨率時間離散格式，本工作采用滿足該條件的Adams－Bashforth離散格式來進行時間離散。

Adams－Bashforth格式用于對流方程時是一個一步，三層和時間前差的差分格式，它的截斷誤差為O（Δt2，Δx2），它可以被看成是對二階時間導數(shù)有限差分近似。Adams－Bashforth的基本形式為：

其中，Γ表示所計算的物理量，n表示時間步，Δt表示時間步長，i，j，k表示節(jié)點序號。

采用Chorin投影法解耦單胞內(nèi)流動的動量守恒方程（4）和連續(xù)性方程，它是一種預估校正法。第一步利用動量守恒方程求解獲得預測速度場，這個預測速度場不滿足連續(xù)性方程。第二步通過連續(xù)性方程推導出用于求解壓力場的泊松類方程，通過求解該方程獲得下一時間步的壓力分布。最后利用第二步求得的壓力場將預測速度場映射到無散度的速度場中，該速度場即為下一時間步的速度場分布。將Brinkman方程的時間項利用上述Adams－Bashforth差分格式離散，并應用到該映射法中，則整個計算過程可以表示為：

（1）已知第n時間步的速度分布un，求解下式獲得預測速度場u＊：

（2）聯(lián)立：

與連續(xù)性方程

獲得泊松類方程：

求解式（9）獲得n＋1時間步的壓力分布pn＋1。

（3）將第一步求得的u＊和第二步求得的pn＋1帶入式（7），獲得n＋1時間步的速度分布。

上述公式中Δt表示時間步長，對于顯式計算的時間推進法，時間步長有如下穩(wěn)定性條件：

其中，Δd表示最小空間網(wǎng)格步長（包含三個方向）。

對于壓力項和擴散項的離散，均采用具有二階精度的中心差分格式：

在以上算法中，每一步迭代均需要且僅需要使用前兩個時間步的流場信息，因此啟動該迭代需要提供起始兩個時間步的流場信息，這兩個時間步的流場信息可以通過很多方法獲得，本工作前兩步采用隱式算法獲得。上述顯式差分計算算法流程如圖6所示。

圖6 算法流程Fig.6 The algorithm flowchart

2.3 細觀流動求解結果分析

為了驗證以上算法的正確性，本節(jié)預測Syncoglass R420織物（表1）的滲透率，其正交單胞幾何建模結果如圖7所示。

表1 Syncoglass R420織物幾何參數(shù)Table1 Geometry parameters of Syncoglass R420fabric

圖7 Syncoglass R420織物單胞模型Fig.7 Meso－geometry model of Syncoglass R420fabric

根據(jù)經(jīng)驗公式求得Syncoglass R420織物紗線軸向和徑向滲透率分別為kalong＝2.71×10－12m2，ktrans＝5.88×10－13m2。采用本文所述算法計算該織物單胞內(nèi)的流場參數(shù)，邊界條件為入口壓力P0＝1Pa，流體黏度μ＝1Pa·s。X方向流動仿真結果的壓力分布和速度分布分別如圖8，9所示。

圖8 X向流動壓力分布Fig.8 Pressure distribution of X－direction flow

圖9 X向流動速度分布Fig.9 Velocity distribution of X－direction flow

將計算所得壓力和速度分布回帶入式（1）獲得單胞滲透率，本工作算法預報的滲透率為Kxx＝1.657×10－10m2，Kyy＝1.873×10－10m2。戴福洪等［11］針對同一織物，采用均勻化方法，不考慮紗線的可滲透性而預報的結果為Kxx＝1.29×10－10m2，Kyy＝1.228×10－10m2。而實驗測定的該織物滲透率為Kxx＝1.43×10－10m2，Kyy＝1.79×10－10m2。比較這幾個結果可以看出，在忽略紗線可滲透性情況下預測的滲透率值小于實驗值，而本工作的滲透率預測值略大于實驗值?？紤]到實驗過程中模具對預成形體的擠壓作用，本工作的預測結果較為準確，證明了該算法的可行性。

3 剪切變形對滲透率的影響

3.1 剪切單胞滲透率預測

預成型體在發(fā)生剪切變形時，不可避免地引起厚度以及紗線內(nèi)部纖維排布和狀態(tài)的變化，但是由于這些變形量很小，而且相對于剪切角度，這些因素對單胞內(nèi)流動的影響非常小，因此本工作忽略這些因素，僅考慮剪切角對細觀滲透率的影響。單胞剪切鎖合角度（最大剪切角度）一般為15～35°，很少超過40°，與織物拓撲和材料相關。對正交單胞模型進行簡單處理即可獲得剪切單胞模型，剪切角為θ的平紋織物單胞如圖10所示。

圖10 剪切角為θ的平紋織物單胞模型Fig.10 Plain fabric unit－cell with shear angleθ

對剪切單胞樹脂流動區(qū)域做貼體坐標系如圖11所示（僅剪切平面），將有剪切的物理域流動求解問題轉換到無剪切的計算域進行求解。（x，y）和（ε，η）分別為物理域和計算域的自變量。

對于三維剪切單胞的貼體坐標的物理域和計算域自變量有如下關系：

其中，z和ζ分別為單胞厚度方向在物理域和計算域的自變量；θ為剪切角。

可以看出，求解區(qū)域從物理域到計算域的轉換可以大大簡化數(shù)值計算的難度，但是樹脂流動控制方程式（4）僅適用于物理域，需要將其向計算域轉化，可以說求解區(qū)域的簡化是以控制方程的復雜化為代價的。下面給出轉換到計算域中控制方程的具體形式。

圖11 貼體坐標系轉換示意圖（a）物理域；（b）計算域Fig.11 Coordinate transformation of body－fitted coordinates（a）physical domain；（b）computational domain

對于Brinkman方程式（4），在計算域上的表達式為：其中，S＝K－1yarn。連續(xù)性方程式（8）在計算域上的具體形式為：

即：

由于單胞剪切造成的紗線微觀滲透率主方向的與坐標系方向產(chǎn)生夾角，在計算之前需要進行滲透率方向變換。

設滲透率主方向與坐標系方向夾角角度為θ，為了進行數(shù)值計算，需將該處紗線滲透率Kyarn轉換到坐標系方向K′上，如圖12所示。根據(jù)張量理論，有如下轉換關系：

式中，kalong和ktrans分別為紗線的微觀軸向和徑向滲透率；且有k21＝k12。

圖12 滲透率方向轉換Fig.12 The transform of permeability direction

將式（15）和（16）作為控制方程，利用上一節(jié)所述顯式計算方法進行剪切單胞下的流動參數(shù)的求解，即可獲得流場參數(shù)，進而預測滲透率。

3.2 滲透率與剪切角的關系

單胞剪切角為20°時，以Y向為壓力差方向的流場仿真結果如圖13，14，15所示，分別表示壓力、Y向速度和X向速度分布。

圖13 剪切單胞壓力分布Fig.13 Pressure distribution of sheared unit－cell

圖14 剪切單胞Y向流動速度分布Fig.14 Y－direction velocity distribution of sheared unit－cell

圖15 剪切單胞X向流動速度分布Fig.15 X－direction velocity distribution of sheared unit－cell

從上述數(shù)值計算結果可以看出，單胞發(fā)生剪切變形后，Y方向的壓力差會產(chǎn)生較大的X 方向的樹脂流動，因此坐標軸X，Y向不再是單胞滲透率的主方向。根據(jù)仿真結果可以分別計算Kxx，Kyy，Kxy，進而通過張量的標準型變換可獲得剪切單胞主滲透率值K11和K22。剪切角每隔5°進行流動仿真，獲得剪切單胞主滲透率比值與剪切角的關系如圖16所示。通過與文獻［3］，［8］比較發(fā)現(xiàn)，該關系曲線與實驗和理論計算所得主滲透率值－剪切角關系相吻合，證明了以上剪切單胞滲透率預測算法的正確性。

圖16 單胞主滲透率比值與剪切角的關系Fig.16 Relationship between the ratio of principal permeability and shear angle

4 結論

（1）考慮紗線的可滲透性，利用Brinkman方程建立了表達織物單胞內(nèi)部樹脂流動的控制方程，達到了將紗線內(nèi)部和外部流動統(tǒng)一表達的目的。

（2）建立逼近真實的正交單胞幾何模型，基于Adams－Bashforth差分格式和Chorin投影法構造了求解樹脂流動控制方程的高分辨率TVD數(shù)值格式，利用達西定律預測了單胞的滲透率，算例表明該算法預測值與實驗值更接近，驗證了算法的準確性。

（3）采用貼體坐標法完成了剪切變形后流動控制方程從物理域到計算域的轉換，進而實現(xiàn)了剪切單胞滲透率的數(shù)值預測，考察了剪切單胞主滲透率比與剪切角的關系，比較表明與文獻中的實驗和理論計算結果相吻合，證明了該剪切單胞滲透率預測算法的合理性。

［1］LOUIS M，HUBER U.Investigation of shearing effects on the permeability of woven fabrics and implementation into LCM simulation［J］.Composites Science and Technology，2003，63：2081－2088.

［2］HAMILA N，BOISSE P.Simulations of textile composite reinforcement draping using a new semi－discrete three node finite element［J］.Composites：Part B，2008，39：999－1010.

［3］VERLEYE B，CROCE R，GRIEBEL M.Permeability of textile reinforcements：simulation，influence of shear and validation［J］.Composites Science and Technology，2008，68：2804－2810.

［4］VERLEYE B，LOMOV S V，LONG A.Permeability prediction for the meso－macro coupling in the simulation of the impregnation stage of resin transfer moulding［J］.Composites：Part A，2010，41：29－35.

［5］TAKANO N，ZAKO M，OKAZAKI T，et al.Microstructurebased evaluation of the influence of woven architecture on permeability by asymptotic homogenization theory［J］.Composites Science and Technology，2002，62：1347－1356.

［6］CHEN Z R，YE L，LU M.Permeability predictions for woven fabric preforms［J］.Journal of Composite Materials，2010，44（13）：1569－1586.

［7］李海晨，張明福，王彪.RTM工藝增強纖維滲透率測量方法研究［J］.航空材料學報，2001，21（1）：51－54.LI H C，ZHANG M F，WANG B.Method on measuring fibre permeabilities in resin transfer molding［J］.Journal of Aeronautical Materials，2001，21（1）：51－54.

［8］田正剛，祝穎丹，張垣，等.剪切效應對纖維增強材料滲透率的影響［J］.武漢理工大學學報，2005，27（2）：4－12.TIAN Z G，ZHU Y D，ZHANG Y，et al.Influence of shearing effects on the permeability of fiber reinforcement［J］.Journal of Wuhan University of Technology，2005，27（2）：4－12.

［9］祝君軍，段躍新，陳吉平，等.碳纖維經(jīng)編織物定型參數(shù)及滲透特性［J］.復合材料學報，2012，29（3）：42－48.ZHU J J，DUAN Y X，CHEN J P，et al.Packifier parameters and permeability characteristics of non－crimp stitched carbon fabric［J］.Acta Materiae Compositae Sinica，2012，29（3）：42－48.

［10］羅云烽，段躍新，王偉，等.VIP工藝中L型構件滲透規(guī)律的實驗研究［J］.材料工程，2009，（增刊2）：15－19.LUO Y F，DUAN Y X，WANG W，et al.Experimental study on permeability of L－type preform for vacuum infusion process［J］.Journal of Materials Engineering，2009，（Suppl 2）：15－19.

［11］戴福洪，張博明，杜善義.用均勻化方法預報平紋織物的滲透率［J］.復合材料學報，2009，26（2）：90－93.DAI F H，ZHANG B M，DU S Y.Permeability prediction of fabric preform using homogenization method［J］.Acta Materiae Compositae Sinica，2009，26（2）：90－93.

［12］倪愛清，王繼輝，朱以文.復合材料液體模塑成型工藝中預成型體滲透率張量的數(shù)值預測［J］.復合材料學報，2007，24（6）：50－56.NI A Q，WANG J H，ZHU Y W.Numerical prediction of saturated permeability tensor of a woven fabric for use in the fluid simulation of liquid composite molding［J］.Acta Materiae Compositae Sinica，2007，24（6）：50－56.

［13］吳炎，晏石林，譚華.具有空隙尺寸雙尺度特性的纖維織物滲透率的預測［J］.固體力學學報，2008，29（S）：80－84.WU Y，YAN S L，TAN H.Predication of permeability in a glass－fiber mat with dual－scale pore－size［J］.Chinese Journal of Solid Mechanics，2008，29（S）：80－84.

［14］GEBART B A.Permeability of unidirectional reinforcements for RTM［J］.Journal of Composite Materials，1992，26（8），1100－1133.