基于橋聯模型預測層合板分層萌生載荷

顧嘉杰, 黃爭鳴

（同濟大學 航空航天與力學學院，上海 200092）

復合材料因比強度和比剛度高、可設計性好等特點，在航空結構中廣泛使用，用量可占飛機總重的50%以上[1]。而層合板的分層破壞是復合材料中最常見的破壞形式之一[2-4]，其原因是層合板受載時，各單層板的泊松比和拉剪耦合系數不匹配，自由邊會產生很高的層間應力，導致層合板提前破壞[5-6]，因此，對層合板的分層研究具有重要意義。

預測層合板的分層萌生的方法目前主要有兩種：強度方法[7-12]和能量方法[13-17]。強度方法是計算出層間應力，與層間強度進行比較，該方法需引入特征長度來解決層間應力的奇異性問題[18-19]，但對于特征長度還沒有統一確定方法。Kim等[7]假設特征長度為一個單層板厚度，層間強度為單層板的橫向強度，這種確定方法也被其他學者使用[8]。Sun等[9]則選擇了6個特征長度來預測分層萌生，確定特征長度為兩個單層板厚度最合適。Brewer等[10]認為特征長度應通過實驗測得，與單層板厚度無關。Lorriot[11-12]提出了一個確定特征長度和層間強度的方法，通過一種層合板的分層實驗數據得到的臨界長度和層間強度，再應用到其他角度層合板的分層預測中。能量方法是根據層間能量釋放率來判斷分層是否發生，需要臨界能量釋放率和有效裂紋長度[13]。Martin等[14]將強度和能量兩種方法結合起來，認為只有兩個條件同時滿足分層才會發生。此外內聚力單元（cohesive element）也用來模擬層合板分層的萌生和擴展，需要層間的剛度、強度和臨界能量釋放率等參數。Turon[15]等對層間剛度提出了一個顯示公式，認為當層間上某一點不再承受應力時，則分層發生。Uguen[16]等分析了這三種參數對分層萌生的影響，預測了層合板[±θn]s的分層載荷，最大誤差27.3%。Mohammadi等[17]用有厚度的界面單元來模擬分層萌生，認為當界面上某一點的應力達到層間強度時，則分層出現。

層合板本質上由基體樹脂將單層板粘在一起，其分層萌生一定由基體樹脂控制，基體樹脂性能對層合板的分層有很大影響[20]。本工作對層合板的分層萌生提出一個基于細觀力學橋聯模型的分析方法。由組分材料預測出單層板的性能，然后通過有限元計算層合板的應力場，用橋聯模型計算纖維和基體樹脂的應力，并通過Mohr判據來判斷基體樹脂的破壞，最后對T800/914層合板[±θn]s進行預測。

1 有限元模型

1.1 材料參數

Lorriot[11]對T800/914層合板[±θn]s進行了分層萌生的實驗測定，文中只給出了單層板的彈性常數，但根據文獻[21]由橋聯模型預測的IM7/914-C單層板的性能與T800/914單層板相匹配，見表1，因此用IM7/914-C中組分材料的性能來代替T800/914，見表2。

表1 單層板的力學性能Table1 Mechanical properties of UD composites

表2 組分材料IM7和914-C的力學性能[21]Table2 Original mechanical properties on IM7 fiber and 914-C matrix[21]

通過表2中組分材料預測單層板的性能時，纖維體積含量Vf調整為0.55來使兩種單層板的性能一致。從分層實驗[22]中看出，分層萌生之前層合板的應力應變曲線保持線性，因此本工作也沒考慮組分材料的非線性。

1.2 層合板離散

層合板[±θ]s的尺寸如下：長度L = 40 mm、寬度b = 15 mm、總厚度h = 4t，單層板厚度t = 0.125 mm。如圖1所示，坐標X、Y、Z分別對應層合板的長度、寬度和厚度。一個均勻的拉伸應力σ作用在表面X = L。由于對稱性，實際建模只取層合板的上半部分，邊界條件為：

式中，u，v，w為對應著X、Y、Z方向的位移。

圖1 層合板[±θ]sFig.1 Geometry of a [±θ]s laminate

以ABAQUS軟件為平臺，對每個單層賦予預測的單層板性能和相應的角度，用8節點實體單元（C3D8R）劃分網格，在X-Y平面，寬度方向劃分136個單元，長度方向107個單元，并且越靠近自由邊網格越密，如圖2所示。每個單元在局部坐標下可看成是單層板的一個特征體元（RVE）。

圖2 層合板X-Y平面的網格劃分Fig.2 Meshing of X-Y plane of a laminated plate

2 內應力計算

2.1 坐標轉換

求特征體元中纖維和基體的應力之前，有限元計算的各單元應力{σ}G= {σxx，σyy，σzz，σyz，σxz，σxy}T需通過式（2）轉化成局部坐標下的應力{σ} = {σ11，σ22，σ33，σ23，σ13，σ12}T，如圖3所示，x、y、z對應整體坐標，x1、x2、x3對應局部坐標。

式中，[Tij]c是坐標變換矩陣，可見參考文獻[23]。

圖3 層合板局部坐標和整體坐標示意圖Fig.3 Global coordinates and local coordinates

2.2 橋聯模型

橋聯模型[24]基本思路是：復合材料受外載作用且沒有破壞時，其組分材料（纖維和基體樹脂）內應力之間必存在一個系數矩陣，使兩者之間互為函數關系。設{σm}為樹脂應力，{σf}為纖維應力，則有：

且任何復合材料細觀力學理論模型都必須滿足基本方程：

結合式（3）和式（4）確定纖維和基體樹脂的應力：

三維橋聯矩陣[A]可用以下形式表達：

其中[I]為6階單位矩陣；{σij}為單元的應力向量；α和β為橋聯參數。文獻[24]通過對18種單向板的剛度預測，發現橋聯參數α和β均取0.3時與測量值綜合吻合最好，因此，將橋聯參數均取0.3。Em為基體樹脂的彈性模量；Gm為基體樹脂的剪切模量；νm是基體樹脂的泊松比； Ef,11為纖維軸向的彈性模量； Ef,22為纖維橫向的彈性模量， νf,12是纖維面內的泊松比；Gf,12為纖維的面內剪切模量。

2.3 分層判據

如前所述，分層萌生源自基體樹脂破壞，本研究采用Mohr判據[25]來檢測基體樹脂的破壞情況。Mohr認為在各種應力狀態下材料發生破壞時，由主應力σ1和σ3所作的所有應力圓有一條公共的包絡線，如果材料某一應力狀態下的最大應力圓與該包絡線內切，材料就破壞。本研究用基體樹脂的單軸拉伸實驗和純剪切實驗破壞時的應力圓的二次曲線公切線，來近似該包絡線，如圖4所示，σu,t是樹脂的拉伸強度，σu,s是樹脂的剪切強度。

基體的應力狀態{σ} = {σ11，σ22，σ33，σ23，σ13，σ12}T用{，，，，，}T= δ{σ11，σ22，σ33，σ23，σ13，σ12}T替代，假設二次曲線和應力圓的方程分別為：

圖4 由純剪切和單軸拉伸實驗近似的Mohr包絡線Fig.4 Mohr envelope approximated by pure shear and uniaxial tensile tests

當材料破壞時，二次曲線和應力圓相切且只有一個切點，其充分必要條件是：

當從式（14)中解出0 ＜ δ ≤ 1，意味著材料發生破壞，否則材料將不會破壞。若基體樹脂產生了拉伸破壞（σn＞ 0），認為層合板萌生了分層。

3 結果與討論

3.1 有限元模擬結果

層合板[±θn]s受軸向拉應力時，（+θ/-θ）層間應力σzz和σyz非常小。若層合板發生分層， 必然由層間應力σxz引起[11]。通過在各個單層板厚度方向上劃分不同數目的單元來研究自由邊層間應力分布。

以層合板[±10]s為例，圖5為單層板厚度方向上劃分不同單元數（m = 2、4、6、8、10）計算出的層間應力σxz在層間（+10/-10）上沿寬度方向的分布，橫坐標是自由邊沿內寬度方向的距離，縱坐標是層間應力大小，m代表單層板厚度方向的單元數，Lorriot[11]計算的層間應力也列于圖5中。從圖5中可知，層間應力σxz只存在于自由邊很小區域中，其寬度是一個層合板厚度h。當單層板厚度方向上劃分不同數目單元時，在自由邊y = 0處有限元計算的層間應力σxz不同，即自由邊緣處的層間應力σxz不收斂，這種情況早已發現[18-19]。這是由單層板性能的不連續造成的，因為應力奇異性只出現在層間（+θ/-θ）上，在性能連續的中面則沒有出現[19]。

圖5 厚度上劃分不同數目單元時T800/914層合板[±10]s在界面（+10/-10）z = t上層間應力 σxz沿y方向的應力分布Fig.5 Interlaminar stress σxz at interface (+10/-10) z = t in T800/914 [±10]s laminate versus y for different numbers of elements along thickness of each layer

圖6 為單層板厚度方向上劃分不同單元數（m = 2、4、6、8、10）計算出的層間應力 σxz在自由邊y = 0沿厚度方向的分布，并與文獻[12]進行了對比，橫坐標是對稱面沿厚度方向的距離，縱坐標是層間應力大小，在層合板上表面z = 0.5h和對稱面z = 0處，層間應力σxz等于0，而在層間（+10/-10）z= 0.25h處，層間應力σxz達到最大值，與實驗觀察到的分層位置一致[12]，同樣能觀察到應力奇異性。

圖6 厚度上劃分不同數目單元時層合板[±10]s在自由邊y = 0處層間應力σxz沿z方向的應力分布Fig.6 Interlaminar stress σxz at y = 0 in T800/914 [±10]s laminate versus z for different number of elements through thickness of each layer

對于應力奇異性，可通過沿單層板厚度方向對層間應力σxz平均來消除。假定單層板沿厚度方向上劃分m個單元，則在單層板中將x、y坐標相同而z（厚度）坐標不同的m個單元的應力進行算術平均。圖7為沿單層板厚度劃分不同單元數（m =2、4、6、8、10）時，厚度方向平均后的層間應力σxz沿寬度方向的應力分布。從圖7可知，當單層板厚度方向上單元數m大于2時，平均后的層間應力σxz在自由邊y = 0處是一個常數，與m無關。因此將單層板厚度方向平均后的應力作為橋聯模型的輸入參數，計算該點處基體樹脂和纖維的應力。

圖7 平均后層合板T800/914 [±10]s的層間應力σxz沿y方向的應力分布Fig.7 Average interlaminar stress σxz in T800/914 [±10]s laminate versus y

3.2 分層萌生預測

利用有限元方法、橋聯理論和基體樹脂的破壞判據，只需提供基體樹脂的原始強度即可預測層合板分層萌生載荷，T800/914層合板[±10n]s, [±20n]s和 [±30n]s（n = 1～5）的分層萌生載荷的實驗值[11]和預測值見表3。

從表3可知，除了層合板[±304]s的誤差是21.5%，其他層合板的預測值與實驗值的誤差都在20%之內，平均誤差是10.3%，與文獻[14]用內聚力單元預測出的分層載荷的最大誤差27.3%相比要低。預測值表明：相同角度的單層板的層數n越多，分層萌生的載荷就越小；層數n相同時，角度θ越大，分層萌生的載荷越小，這與實驗結果一致。除了層合板[±105]s外，其他層合板的預測值都偏大，這可能是沒有考慮纖維和基體樹脂的界面脫粘造成的。雖然分層由基體樹脂的破壞所控制，但纖維和基體樹脂之間界面脫粘往往更早出現[26]，一旦界面發生脫粘，就會導致基體樹脂提前破壞。

結果顯示，橋聯理論預測的單層板性能與文獻[11]一致，結合有限元計算出層合板的應力場也與文獻[11]吻合。在此基礎上，橋聯模型計算出的基體樹脂應力場在自由邊附近很復雜，以層合板[±10]s為例，如圖8所示，在局部坐標下，除基體樹脂的面外正應力很小外，其他應力都存在，用Mohr判據來判斷基體樹脂的破壞能夠綜合所有應力的影響，而且只需基體樹脂的原始強度就能預測出層合板[±θn]s分層萌生載荷。

表3 T800/914層合板[±θn]s的預測值與實驗值[11]對比Table3 Predicted and measured delamination initiation loads[11] for T800/914 angle-plied laminates

圖8 橋聯模型得到層合板[±10]s的基體應力沿寬度方向的應力分布Fig.8 Stresses of matrix determined by bridging model in T800/914 [±10]s versus y

4 結論

（1）有限元計算層合板[±θn]s的應力時，自由邊層間應力存在奇異性，可通過沿單層板厚度方向平均予以解決。

（2）研究層合板[±θn]s分層萌生問題時，通過基體樹脂的單軸拉伸強度和純剪切強度近似的破壞包絡線，可足夠精確地判斷基體樹脂在復雜應力狀態下的破壞情況。

（3）相比于依賴大量分層實驗數據的宏觀方法，橋聯模型僅根據組分材料性能，就能預測出T800/914層合板[±θn]s的分層萌生載荷。