蘆冠達(dá)，黃爭(zhēng)鳴

同濟(jì)大學(xué) 航空與力學(xué)學(xué)院,上海 200092

復(fù)合材料因高比強(qiáng)度、高比模量、可設(shè)計(jì)性好等優(yōu)點(diǎn)，在航空航天領(lǐng)域得到廣泛應(yīng)用，其極限承載能力是工程設(shè)計(jì)人員最為關(guān)注的一項(xiàng)性能指標(biāo)。Huang的橋聯(lián)模型[1-2]，僅需組份材料的性能參數(shù)，即可完成對(duì)復(fù)合材料彈塑性應(yīng)力場(chǎng)的計(jì)算，而且所有公式皆為解析式，十分方便于工程應(yīng)用。橋聯(lián)模型誕生不久，受邀參加了第一屆世界范圍破壞評(píng)估(WWFE-I)[3]，預(yù)報(bào)精度雖超過(guò)了其他參評(píng)的細(xì)觀(guān)力學(xué)理論，但總體評(píng)分上尚不如一些宏觀(guān)理論[4-5]。

近年，橋聯(lián)模型在復(fù)合材料破壞分析和強(qiáng)度預(yù)報(bào)方面取得重要進(jìn)展，導(dǎo)出了添加纖維后基體中的各應(yīng)力集中系數(shù)[6]，修正了層合板最終破壞判據(jù)[7]，完善了纖維和基體界面開(kāi)裂后對(duì)材料承載能力的影響[8]，添加了層合板非致命破壞后的部分剛度衰減策略[9]。這些工作使橋聯(lián)理論預(yù)報(bào)WEEI-I的精度顯著提升，整體精度比賽會(huì)第一名Zinoviev理論的精度更高[10]，更重要的是，組分材料的輸入?yún)?shù)與其原始性能完全相同，無(wú)需由任何復(fù)合材料的實(shí)驗(yàn)數(shù)據(jù)來(lái)反演確定輸入?yún)?shù)。除了組分材料的原始性能，僅需提供單向復(fù)合材料的橫向拉伸強(qiáng)度，以確定界面粘接強(qiáng)弱，這給實(shí)際應(yīng)用帶來(lái)極大方便。

有限元(FE)數(shù)值模擬已成為工程結(jié)構(gòu)尤其是復(fù)雜結(jié)構(gòu)分析不可或缺的工具，在復(fù)合材料領(lǐng)域也如此。但目前的商用有限元軟件，主要只提供基于宏觀(guān)破壞判據(jù)的復(fù)合材料強(qiáng)度計(jì)算模塊，需要使用者提供眾多復(fù)合材料測(cè)試數(shù)據(jù)，這給用戶(hù)帶來(lái)很大不便。歐洲航天局在1992年開(kāi)始開(kāi)發(fā)ESAComp軟件，使用經(jīng)典層合板理論計(jì)算各層應(yīng)力值，并把初始層破壞(First Ply Failure)作為最終破壞判據(jù)[11]。2002年，F(xiàn)riebel等開(kāi)發(fā)了Digimat材料用戶(hù)子程序[12]，主要基于Mori-Tanaka細(xì)觀(guān)力學(xué)方法，但對(duì)復(fù)合材料的破壞分析，依然需要提供復(fù)合材料的實(shí)驗(yàn)數(shù)據(jù)。2015年，韓國(guó)研究人員Lee等運(yùn)用宏觀(guān)理論P(yáng)uck失效準(zhǔn)則與損傷有限單元法，通過(guò)子程序開(kāi)發(fā)，成功將該理論添加入ABAQUS的計(jì)算中，得到相關(guān)層合板算例的強(qiáng)度譜，并與實(shí)驗(yàn)結(jié)果作了對(duì)比[13]。這些工作豐富了復(fù)合材料結(jié)構(gòu)分析功能，但如何減少對(duì)復(fù)合材料實(shí)驗(yàn)數(shù)據(jù)的依賴(lài)尚有不足。

隨著Huang的橋連理論不斷發(fā)展與完善，如何將該理論與CAE(Computer Added Engineering)軟件結(jié)合分析工程結(jié)構(gòu)的破壞與強(qiáng)度問(wèn)題，成為該理論發(fā)展的一大方向。2007年，Huang首次將橋聯(lián)模型結(jié)合商業(yè)有限元軟件ABAQUS，開(kāi)發(fā)出用戶(hù)子程序UGENS[7,14]。本文在該工作基礎(chǔ)上，將UGENS程序改編，添加了上述橋聯(lián)模型的最新進(jìn)展，并應(yīng)用于WWEF-I考題分析和復(fù)雜層合板結(jié)構(gòu)T型接頭的有限元模擬與強(qiáng)度預(yù)報(bào)，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好。

1 理論概述

橋聯(lián)模型UGENS程序基礎(chǔ)詳見(jiàn)文獻(xiàn)[2]，為節(jié)省篇幅，本文只簡(jiǎn)要列出各有關(guān)公式。

1.1 均值內(nèi)應(yīng)力計(jì)算

(1)

(2)

(3)

(4)

(5)

(6)

上述各式中Aij為橋聯(lián)矩陣A中的元素，A為[2]

(7)

其中：

(8)

應(yīng)用于層合板分析，有[2]：

(9)

(10)

(11)

(12)

(13)

進(jìn)一步將式(12)代入式(1)和式(2)中，可求出纖維和基體的內(nèi)應(yīng)力增量。

1.2 真實(shí)內(nèi)應(yīng)力

(14)

(15a)

(15b)

(16)

(17)

(18)

式中：

(19)

纖維中的點(diǎn)應(yīng)力均勻，其真實(shí)應(yīng)力與均值應(yīng)力相同。

1.3 破壞判據(jù)

當(dāng)前步纖維和基體中的真實(shí)應(yīng)力以及單層板分擔(dān)的外力按如下更新：

(20)

(21)

(σ)l=(σ)l-1+dσl

(22)

基于當(dāng)前真實(shí)應(yīng)力，若纖維或基體已經(jīng)破壞，則認(rèn)為復(fù)合材料達(dá)到了破壞。程序中，纖維采用修正的最大應(yīng)力破壞判據(jù)[15]：

(23)

(24)

一旦纖維破壞即式(23)成立，認(rèn)為層合板達(dá)到了致命破壞，對(duì)應(yīng)的外載定義為層合板的極限強(qiáng)度，程序終止計(jì)算?；w的破壞式(24)成立只是一種非致命破壞，需折減單層板的剛度。本文采用的剛度衰減格式為

(25)

存在纖維始終不破壞的情況，持續(xù)加載及每一層的非致命破壞，必然會(huì)導(dǎo)致層合板變形過(guò)大。為此，在程序中還施加了應(yīng)變控制條件：當(dāng)層合板最大應(yīng)變超過(guò)極限應(yīng)變時(shí)，也視為發(fā)生了致命破壞，即需滿(mǎn)足

(26)

式中：h為粘結(jié)層厚度。

聚合物復(fù)合材料最大應(yīng)變一般不超過(guò)15%[10]，一般層板結(jié)構(gòu)破壞時(shí)的極限應(yīng)變不超過(guò)10%。以下算例中均取10%作為極限應(yīng)變條件。

1.4 用戶(hù)子程序

上述理論發(fā)展均添加到UGENS子程序，與2007版本相比[14]，只需額外多輸入一個(gè)單向復(fù)合材料的橫向拉伸強(qiáng)度，以便確定纖維和基體界面何時(shí)開(kāi)裂，若缺省，認(rèn)為界面始終理想?；w的各應(yīng)力集中系數(shù)，根據(jù)纖維和基體的原始性能參數(shù)計(jì)算得到。

2 單層板與層合板強(qiáng)度預(yù)報(bào)——WWFE分析

2.1 算例介紹

有限元建模均采用薄殼單元。只考慮雙軸向正應(yīng)力加載時(shí)，在殼體結(jié)構(gòu)左下兩邊施以對(duì)稱(chēng)邊界條件，如圖1所示；存在剪應(yīng)力的算例，則用左下節(jié)點(diǎn)鉸支的邊界條件，如圖2所示。以此達(dá)到整個(gè)殼體受均勻載荷的目的。

表1 4種復(fù)合材料的組份材料性能參數(shù)[17]Table 1 Constituent properties of four composites[17]

注：F和M分別表示纖維與基體。

表2 4種復(fù)合材料中基體材料的塑性參數(shù)[17]Table 2 Plastic parameters of matrices in four composites[17]

表3 4種基體材料的應(yīng)力集中系數(shù)及界面開(kāi)裂參數(shù)Table 3 Stress concentration factors and cracking parameters of matrices in four composites

圖1 雙軸加載邊界條件[12]Fig.1 Boundary condition for biaxial loading[12]

圖2 受剪應(yīng)力邊界條件[12]Fig.2 Boundary condition for shear stress[12]

2.2 計(jì)算結(jié)果

根據(jù)表1和表2的輸入數(shù)據(jù)并應(yīng)用表3中的參數(shù)，對(duì)WWFE的14個(gè)考題進(jìn)行了分析計(jì)算。為省篇幅，這里只給出最終結(jié)果的得分情況，如表4所示。評(píng)估組織者共設(shè)置了5大類(lèi)125個(gè)評(píng)分點(diǎn)。當(dāng)理論預(yù)報(bào)值與實(shí)驗(yàn)值的偏差在±10%之內(nèi)時(shí)得分為A，預(yù)報(bào)值與實(shí)驗(yàn)值偏差在±10%～±50%之間得分為B，超出±50%時(shí)得分為C[5]。根據(jù)評(píng)分標(biāo)準(zhǔn)，本文有限元求解A+B之和的得分比例為86.4%，遠(yuǎn)高于WWEF評(píng)估中得分最高的Zinoviev理論76.8%的A+B得分比例。證實(shí)了二維橋聯(lián)理論UGENS子程序預(yù)報(bào)復(fù)合材料強(qiáng)度特性的有效性。

表4 WWFE評(píng)估匯總表Table 4 Summary of results for WWFE

注：RC1～RC5依次表示單層板破壞強(qiáng)度分析、層合板初始破壞分析、層合板最終破壞分析、層合板應(yīng)力應(yīng)變曲線(xiàn)分析、定性特征分析。

3 T型接頭受彎失效分析

3.1 材料參數(shù)

T型接頭是復(fù)合材料結(jié)構(gòu)連接的一種常見(jiàn)形式，在航空領(lǐng)域應(yīng)用廣泛，也受到相關(guān)學(xué)者的廣泛研究。為進(jìn)一步驗(yàn)證二維橋聯(lián)理論UGENS子程序的功效，這里對(duì)文獻(xiàn)[18]報(bào)告的T型接頭在彎曲載荷下的極限破壞進(jìn)行了有限元模擬，并與實(shí)驗(yàn)數(shù)據(jù)對(duì)比。

文獻(xiàn)[18]給出了單層板的力學(xué)性能參數(shù)和層合板幾何數(shù)據(jù)，但未能提供組分材料的性能參數(shù)。本文分析中，纖維材料參數(shù)取自文獻(xiàn)[19]，基體材料參數(shù)取自文獻(xiàn)[20]，不計(jì)基體塑性變形，詳見(jiàn)表5。

根據(jù)文獻(xiàn)[18]，T型接頭由不同結(jié)構(gòu)層合板共固化而成，本文對(duì)兩個(gè)層合板的結(jié)合處以及曲率劇變(即90°彎角)處導(dǎo)致的富樹(shù)脂區(qū)均添加了一層ABAQUS自帶的Cohesive單元。

幾何上，Cohesive單元類(lèi)似無(wú)厚度的平面單元，該單元上可承受3種作用力，即正應(yīng)力tn、剪應(yīng)力ts和tt,分別定義[21]為

(27)

表5 T700/QY9611組分材料的力學(xué)性能[19-20]Table 5 Mechanical properties of T700 fiber and QY9611 matrix[19-20]

圖3 Cohesive單元的線(xiàn)彈性-線(xiàn)軟化本構(gòu)模型Fig.3 Linearly elastic-softening constitutive model for a cohesive element

在圖3所示的線(xiàn)彈性-線(xiàn)軟化本構(gòu)模型中，當(dāng)δ=δ0時(shí)材料屈服，當(dāng)δ=δmax時(shí)材料開(kāi)裂。假設(shè)粘結(jié)層厚度為h,對(duì)應(yīng)的3個(gè)應(yīng)變?yōu)?/p>

(28)

當(dāng)δ<δ0時(shí)，材料為線(xiàn)彈性，此時(shí)有

(29)

當(dāng)δ0≤δ<δmax時(shí)，材料處于軟化區(qū)，此時(shí)有

(30)

式中：D為損傷因子，其大小由圖3確定。

由于粘結(jié)層破壞往往并非單一開(kāi)裂破壞模式所致，應(yīng)考慮混合模式的開(kāi)裂準(zhǔn)則。本文采用ABAQUS提供的二次應(yīng)力判據(jù)和B-K準(zhǔn)則，二次應(yīng)力判據(jù)將正應(yīng)力引起的Ⅰ型斷裂與剪應(yīng)力引起的Ⅱ型斷裂結(jié)合起來(lái)，簡(jiǎn)化為一個(gè)應(yīng)力-位移曲線(xiàn)，其中損傷起始的臨界點(diǎn)由式(31)確定，當(dāng)滿(mǎn)足式(31)時(shí)裂紋開(kāi)始產(chǎn)生[21]。

(31)

鑒于該T型板結(jié)構(gòu)為整體共固化而成，本文將粘結(jié)層取為與純基體的力學(xué)性能參數(shù)相同更為科學(xué)，但文獻(xiàn)[18]采用了另一組粘結(jié)層參數(shù)，一并列于表6中。分別應(yīng)用兩組粘結(jié)層參數(shù)進(jìn)行了有限元模擬。

圖4 COHESIVE單元B-K準(zhǔn)則圖解[21]Fig.4 Schematic of B-K criterion used for COHESIVE elements[21]

表6 粘結(jié)層的力學(xué)性能[18]Table 6 Mechanical properties of cohesive layer[18]

3.2 結(jié)構(gòu)參數(shù)

T型構(gòu)件由兩個(gè)L型腹板、底板、中心富樹(shù)脂區(qū)3部分構(gòu)成，各部尺寸詳見(jiàn)圖5。底板3的厚度為4 mm，含32層鋪層，L型腹板1和2的厚度都為1.5 mm，由13層單層板制成。各部分鋪層角詳見(jiàn)表7。

圖5 T型構(gòu)件結(jié)構(gòu)圖[18]Fig.5 Structural diagram for T-joint[18]

表7 T型構(gòu)件鋪層參數(shù)[18]Table 7 Lamination parameters of T-joint[18]

3.3 有限元建模

結(jié)構(gòu)示意如圖6，在離底板下底面80 mm處施加500 N的水平力，在底部?jī)蓚?cè)距中心90 mm處設(shè)置固支邊界。

層板用S4R殼單元，中央三角富樹(shù)脂區(qū)用C3D8R實(shí)體單元。在曲率巨變或應(yīng)力集中區(qū)域進(jìn)行網(wǎng)格加密。表7中的各子結(jié)構(gòu)彼此接觸面上插入COH3D8 Cohesive單元。

圖6 T型構(gòu)件有限元模型Fig.6 FE modelling for T-joint

3.4 結(jié)果與討論

圖7為T(mén)型結(jié)構(gòu)受載343 N時(shí)的位移云圖。達(dá)到182 N時(shí)，右側(cè)L型層合板在與富樹(shù)脂層連接的部分單元發(fā)生基體的壓縮破壞，程序?qū)Πl(fā)生非致命破壞的單元進(jìn)行剛度折減；當(dāng)荷載為343 N時(shí)，由于局部應(yīng)變過(guò)大，左側(cè)L型腹板與底板連接區(qū)域發(fā)生致命破壞，層合板單元(圖8)的應(yīng)變超過(guò)臨界值，此時(shí)計(jì)算終止。

圖7 T型接頭破壞時(shí)的位移云圖Fig.7 Displacement contour of T-joint under failure load

模擬的T型接頭在加載點(diǎn)處的荷載-位移曲線(xiàn)如圖9所示，極限荷載為343 N，而原文獻(xiàn)中的實(shí)驗(yàn)結(jié)果為288.5 N，兩者相差54.5 N，預(yù)報(bào)誤差為18.9%，低于原模擬結(jié)果19.8%的誤差率。由于極限破壞由臨界應(yīng)變條件控制，進(jìn)一步降低臨界應(yīng)變至5%，預(yù)報(bào)的極限載荷為307.2 N，與實(shí)驗(yàn)對(duì)比誤差更小。因?yàn)楸纠欣w維和基體皆為脆性材料，取5%的臨界應(yīng)變似更合理。這說(shuō)明，將橋聯(lián)理論通過(guò)UGENS用于復(fù)雜復(fù)合材料結(jié)構(gòu)的破壞與強(qiáng)度模擬十分有效，剛度與實(shí)驗(yàn)結(jié)果相比略小，相對(duì)誤差為8.21%。原因可能是取自其他文獻(xiàn)的組分材料彈性性能并非與實(shí)際的完全一致。不出所料，若將粘結(jié)層的強(qiáng)度參數(shù)取自文獻(xiàn)[18]，預(yù)報(bào)的復(fù)合材料強(qiáng)度會(huì)明顯降低，見(jiàn)圖9，這也表明，在接觸面引入Cohesive單元是恰當(dāng)?shù)?。若各結(jié)構(gòu)通過(guò)粘接而成，粘接面接觸不良以及非一致性，將導(dǎo)致粘接強(qiáng)度明顯降低，對(duì)復(fù)合材料結(jié)構(gòu)的承載能力也會(huì)產(chǎn)生顯著影響。

圖8 T型接頭破壞位置及形式Fig.8 Failure modes and locations of T-joint

圖9 載荷-位移曲線(xiàn)對(duì)比Fig.9 Comparison of load-displacement curves

3.5 T-型接頭設(shè)計(jì)

3.5.1 不同組分材料的影響

由于發(fā)生致命破壞的是L型腹板，只需對(duì)該構(gòu)件進(jìn)一步考慮。首先，比較3種不同組分材料體系,其力學(xué)性能參數(shù)如表8所示。

將L型腹板用上述3種材料替換，底板性能依然由表5提供，由此得到的T-型接頭載荷-位移曲線(xiàn)如圖10所示，不同材料體系(用于L腹板)的極限破壞模式如表9所示。

表8 3種復(fù)合材料的組分材料性能參數(shù)Table 8 Constituent properties of three composites

圖10 4種組分材料載荷-位移曲線(xiàn)對(duì)比Fig.10 Comparison of load-sdisplacement curves of T-joints made of 4 constituent systems

表9 4種組分材料體系的極限破壞模式Table 9 Ultimate failure modes of four constituent systems

圖10表明，兩種碳纖維材料體系預(yù)報(bào)的結(jié)構(gòu)剛度相同，但強(qiáng)度明顯有異，這是由于AS4纖維的拉壓強(qiáng)度遠(yuǎn)低于T700纖維的水平。而兩種玻璃纖維/環(huán)氧樹(shù)脂材料預(yù)報(bào)的模量與強(qiáng)度均相近，源于兩者采用了相同性能的纖維。綜合而言，碳纖維樹(shù)脂基體組合比玻璃纖維/樹(shù)脂基體組合有更高的模量與強(qiáng)度，而原結(jié)構(gòu)所用的T700/QY9611材料腹板相比其他3種材料有更強(qiáng)的抗彎曲載荷能力。

3.5.2 不同鋪層角

下面考慮L型腹板不同鋪層角對(duì)結(jié)構(gòu)承載能力的影響。假定層數(shù)、厚度與材料性能不變，設(shè)計(jì)了3種不同鋪層角，如表10所示。

依然取材料T700/QY9611，結(jié)構(gòu)致命破壞模式如表11所示,預(yù)報(bào)的不同腹板鋪層T型結(jié)構(gòu)載荷-位移曲線(xiàn)如圖11所示。

從圖11中可以看出，4種不同鋪層的T型接頭在彎曲載荷下的力學(xué)響應(yīng)差異很大，說(shuō)明鋪層方案對(duì)復(fù)合材料結(jié)構(gòu)優(yōu)化影響顯著。按剛度大小排序依次是：鋪層3>鋪層1>原鋪層>鋪層2；按極限承載能力排序依次是：鋪層3>鋪層1>原鋪層>鋪層2。由此可見(jiàn)，鋪層3的抗彎能力較原鋪層的抗彎載荷與剛度都有顯著提高，為4種鋪層設(shè)計(jì)中的最優(yōu)方案。

表10 T型接頭中3種L型腹板鋪層參數(shù)Table 10 Parameters of three lay-ups of L-web in T-joint

表11 4種鋪層結(jié)構(gòu)極限破壞模式Table 11 Ultimate failure modes of four lay-ups

圖11 4種鋪層的載荷-位移曲線(xiàn)對(duì)比Fig.11 Comparison of load-displacement curves with 4 kinds of lay-ups

4 結(jié) 論

本文將最新發(fā)展的橋聯(lián)理論與有限元軟件ABAQUS結(jié)合，通過(guò)編寫(xiě)UGENS殼單元子程序，實(shí)現(xiàn)了對(duì)復(fù)合材料層合板結(jié)構(gòu)在復(fù)雜載荷下的漸進(jìn)破壞模擬和強(qiáng)度預(yù)報(bào)。通過(guò)對(duì)WWFE中14道考題的模擬預(yù)報(bào)，其結(jié)果相較實(shí)驗(yàn)的精確度提高到86.4%，超過(guò)該競(jìng)賽評(píng)分最高的Zinoviev理論(76.8%)。在此基礎(chǔ)上，對(duì)復(fù)合材料復(fù)雜結(jié)構(gòu)——T型接頭進(jìn)行了有限元極限分析，結(jié)果與有關(guān)文獻(xiàn)中的測(cè)試報(bào)告吻合良好(誤差8.21%)，并進(jìn)一步對(duì)T型接頭中腹板的組分材料與鋪層方案進(jìn)行了優(yōu)化設(shè)計(jì)。應(yīng)用本文介紹的程序和方法，可以對(duì)更復(fù)雜復(fù)合材料或復(fù)合材料與其他材料組合結(jié)構(gòu)進(jìn)行破壞分析、強(qiáng)度預(yù)報(bào)和設(shè)計(jì)選型。

z=x2+ix3=reiφ,N2=aFe-iψ+akeiψ

N3=Faeiψ+e-iψak

H=a(cosψ+2λsinψ)(0.5-F)

C=(k-1)(cosψ-2λsinψ)a2exp[2λ(ψ-π)]

D=(1-k)a3exp[2λ(ψ-π)]

式中：a為纖維的半徑，可取單位值；開(kāi)裂角ψ由下述方程求解：

R(eiφ)=(eiφ-eiψ)0.5+iλ(eiφ-e-iψ)0.5-iλe-iφ

H0=(1-k)(1+4λ2)sin2ψ

J1=kG0-1-2(1-k)ξexp(2λψ)cosψ

J2=2(1-k)ξexp(2λψ)sinψ

J3=2(1-k)ξ·

exp(2λψ)[J1cosψ-J2sinψ]/J2

需要指出的是，ξ=1對(duì)應(yīng)奇異開(kāi)裂，此時(shí)開(kāi)裂角ψ無(wú)解。然而，總可以微調(diào)纖維/基體性能(其測(cè)試存在離散性)使得ξ≠1。

基體的臨界von Mises等效應(yīng)力計(jì)算公式為