王 波 李 昂 楊振宇

（1 中國航空工業(yè)集團有限公司防務工程部，北京 100022）

（2 北京航空航天大學固體力學所，北京 100083）

文 摘 針對平面編織復合材料預制體在曲面鋪覆過程中的局部變形，建立了三維有限元模型，利用商業(yè)有限元軟件Abaqus 模擬了預制體在鋪覆成型過程中的纖維束變形規(guī)律。研究了0°和45°碳纖維織物在球面鋪覆過程的變形過程和局部剪切變形。結果表明，纖維束之間的滑動和纖維束起皺是該預制體在球面鋪覆過程中的典型變形模式。在0°織物的球形鋪覆變形中，0°和90°纖維束的剪切變形角最小，45°方向纖維束的剪切變形角最大；在45°織物的鋪覆變形中，0°和90°纖維束的剪切變形角最大，45°方向纖維束的剪切變形角最小。

0 引言

編織復合材料在航空航天、船舶和汽車等領域具有非常廣泛的應用。隨著結構的復雜程度不斷提高，編織復合材料常常需要制備具有復雜曲率的三維結構，比如球面、錐面等[1]。平面編織復合材料在纖維鋪覆過程中，需要根據(jù)結構的形狀對預制體施加預先的變形，因此鋪覆過程中纖維的位置、纖維取向和局部纖維體積分數(shù)會重新分布。復合材料內部的纖維分布和取向決定了各類載荷作用下的復合材料內部的傳力路徑，紗線束之間的相互作用形式也決定了復合材料的損傷和破壞機制[2]。因此，獲取預制體鋪覆成型后纖維束的變形規(guī)律，對于認識最終復合材料結構的整體力學性能至關重要。

關于纖維預制體的在鋪覆過程中的變形研究，受到國內外學者的廣泛關注。C.MACK 等[3]提出了漁網(wǎng)模型對織物的鋪覆過程進行了預測，建立了可以適用于大部分曲面的通用微分方程，并且利用半球鋪覆實驗驗證了其可行性。R.E.ROBERTSON 等[4]在漁網(wǎng)模型基礎上提出了更接近實際的纖維假設，可以更加精確的預測連續(xù)纖維織物的無褶皺鋪覆。此后，研究者對漁網(wǎng)模型的映射算法、假設和初始條件等進行了完善[5]，并將漁網(wǎng)法應用到了更加復雜的曲面上，如雙曲率曲面[6-7]。為了準確的仿真織物鋪覆過程中的變形過程，有限元方法開始在編織復合材料鋪覆的仿真預測中被廣泛的開發(fā)和采用[8]。M.A.KHAN等[9]建立了一種以3D 桁架單元代替單向纖維、以3D 膜單元代替樹脂的預浸料模型。M.DUHOVIC等[10]通過模擬纖維的實際制造過程開始著重于模擬纖維中復雜的接觸關系，并且以3D 梁單元來模擬單向纖維。彭雄奇等[11]基于各向異性超彈性本構模型，研究了碳纖維編織復合材料在沖壓成型過程中的纖維重新排布和重新取向，并通過實驗對數(shù)值模擬進行了驗證。

本文通過建立平面編織復合材料球形鋪覆的有限元模型，模擬了碳纖維平面編織預制體在成型過程中纖維束的變形規(guī)律，詳細分析了各個部位的纖維位置和取向的變化。本文的研究為評估預制體質量以及后續(xù)預測結構的整體力學行為提供了有意義的參考。

1 有限元模型

1.1 平面編織復合材料預制體單胞結構

圖1為單胞模型的幾何尺寸示意圖，其中紗線束的橫截面幾何尺寸a，b 是紗線的橫截面長、短軸的長度。單胞模型中的單胞經(jīng)紗/緯紗方向長度Wx與單層平紋織物厚度方向上的長度Wy決定了單元胞體的大小。單元胞體模型中幾何參數(shù)之間的關系可由如下的公式計算而得。

圖1 平面編織預制體的單元胞體模型Fig.1 Unit cell model for the plane woven composites preforms

紗線束曲率半徑大小R 為

紗線束之間厚度方向的間隙大小f 為

根據(jù)纖維間的集合關系，紗線束直線段的長度大小L 為

紗線束橫截面的面積大小S 為

紗線束中直段的傾角φ 的大小為

由于經(jīng)紗和緯紗方向上，單胞模型各包含兩條紗線束長度的纖維，因此單胞模型中單方向的完整紗線束長度L0為

纖維體積分數(shù)V 為

式中，η 為纖維斂集率。根據(jù)以上參數(shù)之間的聯(lián)系[12]，即可在有限元軟件中構建出平面編織復合材料的單元胞體模型。

1.2 預制體成型的有限元模型

平紋織物預制體的球形鋪覆實驗如圖2所示，采用一個球形壓頭對纖維布作為一個整體進行沖壓。

圖2 鋪覆成型的有限元模型Fig.2 Finite element model for the draping

為了保證模型網(wǎng)格的計算精度，紗線束均采用結構化網(wǎng)格，并選擇Abaqus 中的線性減縮積分單元C3D8R。由于經(jīng)紗和緯紗的邊界處由于截面變化較大，需要對截面進行網(wǎng)格細分來控制邊界處網(wǎng)格的質量。同時，為了減小求解器的計算量，我們將裝配網(wǎng)格模型分別導出為經(jīng)紗緯紗的孤立網(wǎng)格，這樣在不進行質量縮放的情況下大大提高模型計算效率。

由于纖維預制體在受到?jīng)_壓作用而變形的過程中涉及到大量的接觸作用，本文中均采用罰函數(shù)的方法來處理接觸問題。懲罰接觸力可以表示為

式中，k 是懲罰剛度，d0是未解決過程的初始補償距離，dcur是當前補償距離。如果dcur減小至d0，則將d0重置為dcur。標準的庫侖摩擦模型認為，當?shù)刃Σ翍Ζ觘q小于與接觸壓力p 成比例的臨界應力τcrit=μp 時不產生滑移。

因此，可以根據(jù)臨界剪應力提出滑移條件如下

式中，τmax為用戶自定義的最大剪應力值。在等效應力大于等于臨界應力時，會發(fā)生滑移。如果此時摩擦是各向同性的，則滑移的方向和摩擦應力方向一致，則有

根據(jù)庫侖模型，摩擦系數(shù)被定義為一個與等效滑移速率和接觸壓力有關的公式

因此，本文中摩擦系數(shù)只與滑移速率和接觸壓力有關。根據(jù)已有實驗，本文模擬中采用的纖維束間的摩擦系數(shù)為0.22[13]，摩擦性質的方向性為各向同性。

2 結果與討論

2.1 織物的鋪覆變形

由圖3可見，隨著半球模具的沖壓作用，0°纖維預制體不斷向半球表面鋪覆。同時由于壓邊圈的作用，大部分纖維呈拉伸狀態(tài)。對于纖維取向為0°和90°的纖維來說，其在軸向的彈性模量遠大于其橫向上的彈性模量，因此纖維預制體的0°和90°方向上的纖維拉伸變形較為明顯，邊界中點附近產生了很大的向內縮進，邊角處幾乎無縮進。45°方向上相鄰經(jīng)緯紗之間的角度變化較為明顯，可以看出這部分纖維主要靠互相之間的角度變化來適應整體結構的大變形過程。圖4為45°纖維預制體在半球表面鋪覆過程的變形情況，邊角附近產生了很大的向內縮進，而邊界中點縮進并不明顯。0°和90°方向上相鄰經(jīng)/緯紗之間的角度變化較為顯著，表明存在較大的局部剪切變形。

圖3 0°織物的變形過程Fig.3 Deformation process of the 0° fabric

圖4 45°織物的變形過程Fig.4 Deformation process of the 45° fabric

2.2 纖維束的局部變形

圖5給出了0°平紋編織碳纖維預制體在半球曲面上的纖維局部剪切變形角分布情況，其中圖5（a）橫坐標為測量點處的半徑Ri與半球半徑R 之比。由圖5（a）可以看出，纖維預制體在0°和90°兩個方向上的剪切變形角十分接近于0，變化趨勢趨于平緩。很顯然這是由于0°織物的纖維取向為0°和90°，碳纖維在這個方向的鋪覆成型過程中主要承受軸向拉伸作用，變形主要為軸向變形。同時碳纖維的軸向拉伸模量較大，因此相鄰紗線間并沒有產生明顯的角度變化；相反，在±45°上的纖維隨著數(shù)據(jù)點的位置Ri/R 的增大，剪切變形角不斷增大，且變化趨勢越來越快，最終在臨近半球赤道處達到了剪切變形角的最大值40.66°，并產生了局部起皺現(xiàn)象。這主要是由于在±45°方向，纖維軸向和橫向均受力，而碳纖維的橫向彈性模量遠低于其拉伸模量，因此在橫向產生的變形遠大于軸向，導致紗線間產生的角度變化較大，剪切變形明顯。同時相鄰紗線間的縫隙也相應增大，局部產生明顯的間隙缺陷。在半球面上，緯度60°處，纖維預制體的纖維取向的變化較小，剪切變形不明顯，如圖5（b）所示。而在緯度為30°處，纖維預制體的纖維取向的變化較大，基本規(guī)律為0°、90°、180°經(jīng)線處剪切變形角趨近于0；在45°和135°經(jīng)線處剪切變形角分別達到最大值27.6°和27.8°，同時0°～90°區(qū)域和90°～180°經(jīng)線區(qū)域的變化規(guī)律基本呈對稱趨勢。

同樣，圖6給出了45°平紋編織碳纖維預制體在半球曲面上的局部剪切變形角分布情況。由圖6（a）可以看出，碳纖維預制體在45°和-45°兩個方向上的剪切變形角十分接近于0，變化并不顯著。相反，在0°和90°上纖維的剪切變形角隨著半徑比Ri/R 不斷增大，最終在臨近半球赤道處達到了剪切變形角的最大值41.29°。圖6（b）可見，在緯度60°處，纖維預制體的纖維取向的變化較小，剪切變形不明顯；而在緯度30°處，纖維預制體的纖維取向的變化較大。表現(xiàn)為45°、135°處纖維的剪切變形角趨近于0°，而在0°、180°和90°經(jīng)線處剪切變形角分別達到最大值23.5°、22.4°和23.5°，同時0°～90°區(qū)域和90°～180°經(jīng)線區(qū)域的變化規(guī)律呈對稱趨勢。

圖5 0°織物半球曲面上不同方向的纖維局部剪切變形角分布Fig.5 Local shear deformation angle distribution of fibers in different directions on the hemispherical surface of 0° fabric

圖6 45°織物半球曲面上不同方向的纖維局部剪切變形角分布Fig.6 Local shear deformation angle distribution of fibers in different directions on a hemispherical surface of 45° fabric

通過數(shù)值模擬的結果與已有實驗之間的對比（圖7），可以發(fā)現(xiàn)本文的有限元模擬很好的預測了球形曲面鋪覆中平面編織纖維網(wǎng)絡的變形。

圖7 數(shù)值模擬結果與實驗[14]的對比Fig.7 Comparison between the numerical simulation and experiments

從成型過程中纖維整體的變形規(guī)律來看，預制體中紗線都會受到的拉伸、壓縮和剪切，其中拉伸導致了邊界的縮進；壓縮（主要由壓邊圈限制作用）導致了部分纖維的屈曲；剪切導致了纖維層間的大角度變化。這三者的共同作用之外還會產生層間滑動和起皺等缺陷，如圖8所示。

圖8 纖維束之間的滑動和彎曲變形Fig.8 Sliding between fiber bundles and bending deformation

起皺是由于纖維預制體在鋪覆過程中需要由平面轉變?yōu)槿S曲面，因此在平面向曲面過渡的區(qū)域內預制體局部的幾何形狀會發(fā)生突變，單元會產生較大的變形，即紗線間產生更大的接觸壓力。此時在過渡區(qū)域會產生局部的較大的壓縮變形，纖維會發(fā)生類似梁的“翹曲”，使得半球赤道處的過渡區(qū)域的經(jīng)紗和緯紗發(fā)生起皺現(xiàn)象，與BOISSE 等[15]提出的由于紗線壓縮變形導致的起皺現(xiàn)象吻合。

3 結論

本文采用有限元法分別對0°和45°碳纖維平面編織布在球形壓頭沖壓成型的過程進行了仿真，獲得了與實驗一致的變形模式，詳細分析了纖維束在成型過程中的剪切變形與分布情況。這些數(shù)據(jù)對不同取向編織物鋪覆成型的結構設計和編織布的剪裁具有重要的意義，對于該碳纖維布的成型性能的準確預測，也具有很強的實際應用價值。通過對比分析，得到了以下結論:

（1）滑動和起皺是編織復合材料預制體在球面鋪覆過程中的典型變形模式，纖維的局部剪切變形角可以用作碳纖維編織物鋪覆成型性能的一種評判標準；

（2）鋪覆過程中纖維的局部變形與織物中纖維的取向角密切相關，通過改變織物的鋪層角度，可以調控成型后預制體在不同方向上纖維的變形程度；

（3）本文的有限元方法很好的預測了球形鋪覆的成型過程，并且對平紋編織復合材料預制體的變形做出合理的預測，該方法還可望推廣應用于其他復雜曲率的不可展曲面鋪覆結構的預測。