基于FEA的復合材料結(jié)構(gòu)極限承載失效預測

章 凌，趙優(yōu)存，李 祎，楊 帆，崔 浩

(1.北京宇航系統(tǒng)工程研究所，北京 100076；2.西北工業(yè)大學民航學院，太倉 215513)

0 引言

連續(xù)纖維增強樹脂基復合材料是以樹脂聚合物為基體、連續(xù)纖維等為增強材料，通過復合工藝制備而成的高性能材料，具有高比強度、高比剛度和耐腐蝕等諸多優(yōu)點[1]，廣泛應用于C929寬體客機、長江系列渦扇發(fā)動機以及新一代運載火箭等國產(chǎn)重大裝備[2]。

對于纖維增強復合材料的損傷分析，由于試驗成本較高，業(yè)內(nèi)廣泛使用有限元對試驗進行模擬仿真，并且通過漸進損傷模型預測復合材料結(jié)構(gòu)的失效過程[3]?，F(xiàn)有有限元軟件的內(nèi)置損傷材料模型多數(shù)基于二維層合板理論，有限元分析 (Finite Element Analysis，F(xiàn)EA)軟件中雖然內(nèi)嵌三維正交各向異性模型，但仍然缺乏三維失效判據(jù)[4]，這導致當使用內(nèi)嵌Hashin、LaRC03等失效判據(jù)進行失效判斷或通過用戶自定義子程序引入其他失效判據(jù)時，常需要針對層合板每一層分別建立單層模型[5-6]。這種較為精細的模型計算效率低，計算時間長，為此，Chou等[7]提出一種關(guān)于復合材料層合板的三維等效彈性常數(shù)計算方法，通過各單層板的屬性可計算出層合板的三維等效彈性常數(shù)，從而顯著提高復合材料層合板，特別是厚板問題的計算效率，但Chou等僅是通過該方法對一簡單的層合板問題進行研究，得出各層對應的應力、應變場，仍未解決層合板的失效問題。在此方法基礎上，Bogetti等[8]針對厚板的失效問題，在多尺度計算方法中加入損傷準則，并基于FORTRAN語言編寫有限元子程序求解厚板問題；Staniszewski等[9]在Bogetti的工作基礎上增加了材料的非線性，完善了材料的響應。上述兩種針對厚板的多尺度損傷計算方法，需要在計算過程中不斷獲取單元的信息，包括材料以及鋪層信息，這種讀取會造成計算效率的降低，賈利勇等[4]基于上述研究成果針對厚板問題提出了一種新的非線性多尺度分析方法，該方法通過剛度等效和應力-應變分解建立了一種多尺度分析模型，實現(xiàn)了復合材料厚板結(jié)構(gòu)在子層壓板水平的計算以及在鋪層水平的失效判斷，并通過編寫子程序?qū)秃喜牧鲜讓邮нM行判斷。這種方法僅能判斷復合材料首層的失效，無法解決復合材料層合板的漸進失效問題，同時，該方法未考慮層合板處于小應變階段時各鋪層未進入損傷狀態(tài)，而是在層合板每個應變階段均對各鋪層進行損傷判斷，這會導致計算效率降低。

本文采用并行多尺度的方法，基于經(jīng)典層合板理論根據(jù)各單層板的鋪層信息以及材料剛度計算層合板的等效剛度。在層合板小應變階段不考慮各鋪層損傷，超出一定應變閾值后，引入多種失效判據(jù)對各單層板進行漸進損傷判斷，將仿真結(jié)果與FEA軟件內(nèi)嵌判據(jù)結(jié)果以及試驗結(jié)果進行對比，以驗證該方法的效率與準確性。

1 數(shù)值模型

1.1 并行多尺度計算流程

基于并行多尺度方法，本文編寫的子程序其計算過程如圖 1所示。與通常理解的介于宏觀-細觀-微觀體系的多尺度方法不同，本文所指多尺度參考賈利勇等[4]、Bogetti等[8]以及Stanis-zewski等[9]關(guān)于厚板問題所提出的多尺度理解，即對于層合板從“層合板整體”到“其中每一層”兩個尺度進行分析?；谠摲椒ǖ膶雍习逵邢拊Ｐ?，一個單元中包含N個鋪層，該單元即為一個層合板單元。計算起始時輸入各單層板剛度，從而計算得出層合板的等效剛度。層合板在載荷較低時處于線彈性階段，此時無需對各單層板進行損傷判斷；隨著加載載荷的增大，當層合板的應變大于預設安全閾值時，開始對每一層單層板進行損傷判斷。對損傷后的單層板的剛度進行折減，并反饋計算得出損傷后層合板的等效剛度，最后更新應力矩陣，進行下一個的加載步。

1.2 層合板彈性響應

在并行多尺度方法中，對于層合板有限元模型，一個單元中包含N個鋪層，該單元即為一個層合板單元，通過經(jīng)典層合板理論層合板的等效剛度可由下式計算得出

(1)

(2)

子程序中，通過輸入層合板的等效剛度可從FEA有限元軟件中反饋計算得到層合板的應變以及應力。當層合板應變值大于預設應變閾值時，進入單層板的損傷判斷，基于正交各向異性復合材料的應力應變關(guān)系[10]，相對于材料主軸坐標系，單層板面內(nèi)應力應變關(guān)系為

(3)

其中

(4)

式中，σ11，σ22，τ12是單層板的面內(nèi)應力分量；ε11，ε22，γ12是上述應力對應應變。E1與E2分別是平行于纖維方向和垂直于纖維方向的楊氏彈性模量；G12是面內(nèi)剪切模量；ν12與ν21為面內(nèi)泊松比。

單層板層間應力應變關(guān)系為

(5)

Q33=E3,Q44=G23,Q55=G13

(6)

式中：σ33，τ23，τ13是單層板的層間應力分量；ε33，γ23，γ13是上述應力對應應變。E3是垂直于單層板平面方向的楊氏彈性模量；G23與G13是層間剪切模量。

1.3 漸進損傷模型

本文針對單層板纖維、基體以及層間3個維度的損傷判斷，分別采用基于應力描述的最大應力準則、二維Puck失效準則[11]和Ye分層失效準則[12]來判斷各單層板的初始失效。相較于以往學者僅采用一種失效準則對材料的損傷進行分析，此處所界定的3種準則同時對單層板進行3個維度的失效判斷仍然合理，原因如下：

以下是對材料各方向損傷起始的定義，有效應力被用來判定材料的損傷起始，當纖維出現(xiàn)拉伸破壞(σ11≥0)時

(7)

當纖維出現(xiàn)拉伸破壞(σ11<0)時

(8)

相較于二維Puck失效準則，三維Puck失效準則具有更好的預測結(jié)果，但后者的理論需要進行額外的算法來搜索潛在的斷裂平面，這一步需要消耗大量的計算時間，故此處選擇Puck理論的平面應力格式[13]對基體損傷進行簡單討論，二維Puck失效準則認為基體失效有3種模式，即基體拉伸失效(模式A)、基體壓縮失效(模式B)以及基體由于更大的壓縮應力而失效(模式C)。

當基體出現(xiàn)拉伸破壞(σ22≥0)時，基體破壞為模式A

(9)

(10)

(11)

當層間出現(xiàn)破壞時

(12)

(13)

(14)

(15)

在損傷演化方面，采用雙線性漸進損傷模式對材料剛度進行折減[15-16]，定義損傷變量為

(16)

式中，δ0為破壞起始點的等效位移，有效位移δ的值為有效應變εeff與單元特征長度LC的乘積，δC為完全破壞時的等效位移，其值由混合模態(tài)斷裂韌性GT與完全斷裂時破壞的應力計算得出。

2 驗證試驗

本文中測試層合板材料為T800級碳纖維增強復合材料，其材料參數(shù)如表 1與表 2所示。

表1 T800級碳纖維增強復合材料模量與泊松比

表2 T800級碳纖維增強復合材料強度值

如圖 2所示，開孔拉伸試驗試樣長度L=250 mm，寬度W=36 mm，開孔直徑D=6 mm，層合板鋪層方式為[0°]8、[±45°]4S，每一層單層板厚度均為0.15 mm。

3 算例分析

3.1 單個單元仿真

為驗證在FEA軟件中損傷演化的準確性，采用單個單元仿真方法對子程序進行測試。如圖 4(a)建立1 mm×1 mm單個單元加載模型，在一個端面的4個節(jié)點處進行相同大小位移加載，另一端面的4個節(jié)點處進行固支。此處對以纖維拉伸方向損傷演化為例，圖 4(b)顯示了單元雙線性損傷演化的過程，紅色線為單元應力，藍色線為單元損傷變量，單元應力不超過2 067 MPa 時單元為線彈性、損傷為0，單元應力達到2 067 MPa時開始發(fā)生損傷，此后單元剛度降低、損傷積累，單元損傷達到1.0時完全失效。單元應力峰值為2 067 MPa，與纖維方向拉伸強度值2 071 MPa近似相等。

3.2 開孔板損傷分析

圖 5與圖 6分別是[0°]8、[±45°]4S鋪層層合板拉伸試驗的并行多尺度計算方法仿真結(jié)果以及試驗數(shù)字圖像相關(guān)技術(shù)(Digital Image Correlation，DIC)結(jié)果?？梢钥闯觯簩τ陂_孔層合板拉伸試驗，其損傷主要發(fā)生在孔邊，[0°]8鋪層層合板纖維從孔邊起裂并沿著孔板的寬度方向擴展，基體從孔邊起裂沿著孔板長度方向擴展，這符合試驗完全破壞前一幀的DIC結(jié)果，與試驗結(jié)果相近；[±45°]4S鋪層層合板的損傷主要以基體拉伸破壞為主，基體在孔邊起裂并沿著45°方向擴展，最后形成X形裂紋，仿真結(jié)果同樣符合試驗完全破壞前一幀的DIC結(jié)果。

圖 7、圖 8分別對應[0°]8、[±45°]4S鋪層層合板拉伸的載荷-位移曲線，比較了在相同網(wǎng)格尺寸下的模擬結(jié)果與試驗結(jié)果以及使用FEA軟件中內(nèi)嵌的Hashin損傷準則預測結(jié)果。

表 3中為兩種不同鋪層開孔板拉伸失效載荷的對比，結(jié)果表明，在同樣網(wǎng)格尺寸的情況下使用本文計算方法得出的失效載荷與試驗結(jié)果相近，使用本文計算方法能夠較為準確地預測開孔板拉伸試驗的結(jié)果，誤差能夠保持在5%以下，但在本文計算方法中未考慮剪切非線性的影響，故從圖 8曲線中可以看出，針對[±45°]4S鋪層層合板拉伸試驗的預測，本文方法不能很好地預測損傷前的非線性趨勢。

表3 開孔板拉伸失效載荷誤差對比

3.3 夾筋板損傷分析

圖 9是橫跨兩筋條四點彎試驗的并行多尺度計算方法仿真結(jié)果以及試樣破壞結(jié)果。從仿真結(jié)果可以看出，在筋條與蒙皮的連接處首先出現(xiàn)初始裂紋，并且該裂紋沿著筋條連接處擴展至橫縱筋條交叉點處，與試驗結(jié)果相近。

圖10為橫跨兩筋條四點彎試驗的載荷-位移曲線，比較了在不同網(wǎng)格尺寸下的模擬結(jié)果與試驗結(jié)果以及使用FEA軟件中內(nèi)嵌的Hashin損傷準則預測結(jié)果。

圖1 并行多尺度方法計算流程圖

圖2 開孔板拉伸幾何形狀及試驗條件

(a)試驗條件

(a)纖維方向拉伸加載

(a)纖維損傷示意圖

(a)基體損傷示意圖

圖7 [0°]8鋪層層合板拉伸載荷-位移曲線

圖8 [±45°]4S鋪層層合板拉伸載荷-位移曲線

(a)層間破壞

圖10 橫跨兩筋條四點彎試驗載荷-位移曲線

圖11 不同網(wǎng)格密度下四點彎有限元模型的載荷-位移曲線

表 4為三橫跨兩筋條四點彎試驗的失效載荷的對比，結(jié)果表明，在同樣網(wǎng)格尺寸的情況下使用本文計算方法得出的失效載荷與試驗結(jié)果相近，誤差不超過14%，相較于開孔板試驗的預測誤差，預測四點彎試驗產(chǎn)生較大誤差的原因是在四點彎試驗過程中層合板發(fā)生層間破壞，層合板產(chǎn)生了較大的翹曲，但由于子程序中無法獲取變形曲率值，故利用該并行多尺度方法進行仿真模擬時，難以獲取層合板的彎矩與曲率，導致預測結(jié)果比試驗結(jié)果偏大。

表4 橫跨兩筋條四點彎試驗失效載荷誤差對比

3.4 網(wǎng)格相關(guān)性驗證及計算效率對比

為檢驗本文所提出方法的網(wǎng)格相關(guān)性，在不同的網(wǎng)格密度下對所建立的有限元模型進行計算，所有有限元模型均通過網(wǎng)格相關(guān)性檢驗。圖 11為四點彎有限元模型在不同網(wǎng)格密度下的載荷-位移曲線，其中ele為單元網(wǎng)格尺寸，單位為mm，該結(jié)果證明了本文所提出的方法不會因網(wǎng)格尺寸的變化而造成結(jié)果的較大變化。

本文對計算效率進行橫向?qū)τ诓⒂枰栽u估，使用相同的計算內(nèi)核數(shù)，對相同網(wǎng)格尺寸情況下并行多尺度方法耗時與內(nèi)嵌Hashin損傷判定方法耗時進行對比，計算時間從每個算例的.msg文件中獲取，評估結(jié)果如表 5所示，可得知，相較于使用FEA軟件內(nèi)嵌Hashin損傷判定方法，采用并行多尺度計算方法能夠減少計算總耗時并且提高每一步計算的速度，能夠縮短43%～85%的總計算時長，并且將每步計算效率提高1倍。

表5 計算效率對比

4 結(jié)束語

本文基于經(jīng)典層合板理論通過建立層合板與層合板中每一層單層板的剛度矩陣，利用元軟件分別對層合板、單層板兩個尺度進行損傷判斷。采用基于應力描述的最大應力準則、二維Puck失效準則以及Ye分層失效準則對每一層進行損傷初始判斷。采用雙線性漸進損傷模式對建立了復合材料層合板的損傷演化，并通過折減單層板的剛度反饋給層合板實現(xiàn)層合板的等效剛度折減，完成多尺度的損傷判斷。

本文所使用并行多尺度計算方法通過T800級碳纖維增強復合材料開孔板拉伸試驗以及橫跨兩筋條四點彎試驗的算例驗證，該方法進行了網(wǎng)格相關(guān)性檢查，將仿真結(jié)果與試驗結(jié)果進行對比。結(jié)果表明：通過該方法建立的模型能夠較好地預測T800級碳纖維增強復合材料開孔板拉伸試驗以及橫跨兩筋條四點彎試驗的失效載荷，預測值與試驗值平均誤差低于10%；與使用內(nèi)嵌Hashin損傷判定方法相比，該方法能夠縮短43%～85%的總計算時長，并且將每步計算效率提高1倍；本文所使用并行多尺度方法受網(wǎng)格尺寸影響較小。