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

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

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

0 引言

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

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

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

1 數(shù)值模型

1.1 并行多尺度計(jì)算流程

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

1.2 層合板彈性響應(yīng)

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

(1)

(2)

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

(3)

其中

(4)

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

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

(5)

Q33=E3,Q44=G23,Q55=G13

(6)

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

1.3 漸進(jìn)損傷模型

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

以下是對(duì)材料各方向損傷起始的定義，有效應(yīng)力被用來(lái)判定材料的損傷起始，當(dāng)纖維出現(xiàn)拉伸破壞(σ11≥0)時(shí)

(7)

當(dāng)纖維出現(xiàn)拉伸破壞(σ11<0)時(shí)

(8)

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

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

(9)

(10)

(11)

當(dāng)層間出現(xiàn)破壞時(shí)

(12)

(13)

(14)

(15)

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

(16)

式中，δ0為破壞起始點(diǎn)的等效位移，有效位移δ的值為有效應(yīng)變?chǔ)舉ff與單元特征長(zhǎng)度LC的乘積，δC為完全破壞時(shí)的等效位移，其值由混合模態(tài)斷裂韌性GT與完全斷裂時(shí)破壞的應(yīng)力計(jì)算得出。

2 驗(yàn)證試驗(yàn)

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

表1 T800級(jí)碳纖維增強(qiáng)復(fù)合材料模量與泊松比

表2 T800級(jí)碳纖維增強(qiáng)復(fù)合材料強(qiáng)度值

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

3 算例分析

3.1 單個(gè)單元仿真

為驗(yàn)證在FEA軟件中損傷演化的準(zhǔn)確性，采用單個(gè)單元仿真方法對(duì)子程序進(jìn)行測(cè)試。如圖 4(a)建立1 mm×1 mm單個(gè)單元加載模型，在一個(gè)端面的4個(gè)節(jié)點(diǎn)處進(jìn)行相同大小位移加載，另一端面的4個(gè)節(jié)點(diǎn)處進(jìn)行固支。此處對(duì)以纖維拉伸方向損傷演化為例，圖 4(b)顯示了單元雙線性損傷演化的過(guò)程，紅色線為單元應(yīng)力，藍(lán)色線為單元損傷變量，單元應(yīng)力不超過(guò)2 067 MPa 時(shí)單元為線彈性、損傷為0，單元應(yīng)力達(dá)到2 067 MPa時(shí)開(kāi)始發(fā)生損傷，此后單元?jiǎng)偠冉档汀p傷積累，單元損傷達(dá)到1.0時(shí)完全失效。單元應(yīng)力峰值為2 067 MPa，與纖維方向拉伸強(qiáng)度值2 071 MPa近似相等。

3.2 開(kāi)孔板損傷分析

圖 5與圖 6分別是[0°]8、[±45°]4S鋪層層合板拉伸試驗(yàn)的并行多尺度計(jì)算方法仿真結(jié)果以及試驗(yàn)數(shù)字圖像相關(guān)技術(shù)(Digital Image Correlation，DIC)結(jié)果。可以看出：對(duì)于開(kāi)孔層合板拉伸試驗(yàn)，其損傷主要發(fā)生在孔邊，[0°]8鋪層層合板纖維從孔邊起裂并沿著孔板的寬度方向擴(kuò)展，基體從孔邊起裂沿著孔板長(zhǎng)度方向擴(kuò)展，這符合試驗(yàn)完全破壞前一幀的DIC結(jié)果，與試驗(yàn)結(jié)果相近；[±45°]4S鋪層層合板的損傷主要以基體拉伸破壞為主，基體在孔邊起裂并沿著45°方向擴(kuò)展，最后形成X形裂紋，仿真結(jié)果同樣符合試驗(yàn)完全破壞前一幀的DIC結(jié)果。

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

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

表3 開(kāi)孔板拉伸失效載荷誤差對(duì)比

3.3 夾筋板損傷分析

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

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

圖1 并行多尺度方法計(jì)算流程圖

圖2 開(kāi)孔板拉伸幾何形狀及試驗(yàn)條件

(a)試驗(yàn)條件

(a)纖維方向拉伸加載

(a)纖維損傷示意圖

(a)基體損傷示意圖

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

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

(a)層間破壞

圖10 橫跨兩筋條四點(diǎn)彎試驗(yàn)載荷-位移曲線

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

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

表4 橫跨兩筋條四點(diǎn)彎試驗(yàn)失效載荷誤差對(duì)比

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

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

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

表5 計(jì)算效率對(duì)比

4 結(jié)束語(yǔ)

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

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