帶裂紋功能梯度材料薄板SIFs分析的廣義參數(shù)Williams單元

徐 華， 楊 濤， 韓林君， 楊綠峰,3

(1. 廣西大學(xué) 土木建筑工程學(xué)院/工程防災(zāi)與結(jié)構(gòu)安全教育部重點(diǎn)實(shí)驗(yàn)室/廣西防災(zāi)減災(zāi)與工程安全重點(diǎn)實(shí)驗(yàn)室， 廣西 南寧 530004；2. 中國(guó)電建集團(tuán)華中電力設(shè)計(jì)研究院有限公司， 河南 鄭州 450007； 3. 廣西壯族自治區(qū) 住房和城鄉(xiāng)建設(shè)廳， 廣西 南寧 530028)

功能梯度材料[1](functionally graded material， FGM)是一種材料屬性呈連續(xù)變化的新型功能材料，因其具有緩解熱應(yīng)力、釋放殘余應(yīng)力及提高黏結(jié)的能力，被廣泛應(yīng)用于航空航天、能源、建筑與生物醫(yī)學(xué)等領(lǐng)域.受生產(chǎn)過(guò)程及工作環(huán)境等因素影響，F(xiàn)GM與均勻材料一樣不可避免地產(chǎn)生裂紋，且其材料屬性的非均勻性導(dǎo)致FGM結(jié)構(gòu)中裂紋的斷裂行為更為復(fù)雜，通常呈現(xiàn)為混合形式.在外荷載作用下，這些裂紋易進(jìn)一步擴(kuò)展，嚴(yán)重時(shí)會(huì)導(dǎo)致FGM結(jié)構(gòu)的整體破壞.因此，為提高帶裂紋FGM結(jié)構(gòu)的安全性，需深入研究裂尖斷裂力學(xué)參數(shù)的變化規(guī)律，而應(yīng)力強(qiáng)度因子(stress intensity factors，SIFs)作為結(jié)構(gòu)斷裂分析的重要參數(shù)之一，表征了裂尖應(yīng)力-應(yīng)變場(chǎng)的強(qiáng)弱程度，對(duì)深入研究帶裂紋FGM結(jié)構(gòu)的破壞機(jī)理及提高結(jié)構(gòu)的安全性具有重要意義.

近年來(lái)，眾多學(xué)者圍繞帶裂紋FGM結(jié)構(gòu)中裂尖SIFs問(wèn)題進(jìn)行了大量研究.Guo等[2-3]建立了分段指數(shù)模型(PE模型)，通過(guò)求解奇異積分方程組，研究了功能梯度分層結(jié)構(gòu)中穿越界面任意方向的裂紋，發(fā)現(xiàn)其夾角對(duì)混合型裂紋尖端SIFs的影響較大.Pan等[4]在Guo等[2-3]建立的PE模型基礎(chǔ)上，同樣采用求解奇異積分方程組研究了含共線裂紋的任意熱機(jī)械屬性功能梯度條I型熱應(yīng)力，分析了非均勻常數(shù)和裂紋幾何參數(shù)對(duì)裂尖熱應(yīng)力強(qiáng)度因子(TSIF)的影響.Wang等[5]基于開(kāi)爾文基本解建立了FGM結(jié)構(gòu)斷裂分析的邊界-域積分方程，采用徑向積分邊界元法分析含圓盤(pán)狀裂紋的FGM三維體中裂尖SIFs，研究表明材料梯度平行于裂紋面方向變化對(duì)I型和II型SIFs均有顯著影響.Zhang等[6]應(yīng)用徑向積分邊界元法對(duì)三維連續(xù)非均勻各向同性線彈性FGM結(jié)構(gòu)進(jìn)行了斷裂分析，研究了材料梯度對(duì)裂紋張開(kāi)位移和SIFs的影響.前述文獻(xiàn)采用的奇異積分方程組和徑向積分邊界元法均是半解析半數(shù)值方法，需經(jīng)大量繁瑣的理論推導(dǎo)過(guò)程，難以適用于非規(guī)則邊界和復(fù)雜荷載條件的結(jié)構(gòu).隨著計(jì)算機(jī)運(yùn)算能力的提升，數(shù)值方法得到了迅速發(fā)展，Abdollahifar等[7]首次將無(wú)網(wǎng)格局部彼得洛夫伽遼金法用于分析帶邊界裂紋功能梯度條的I型裂紋裂尖SIFs；隨后，又采用此方法對(duì)帶裂紋FGM板的I型和II型裂紋SIFs、最大能量釋放率、裂紋擴(kuò)展角和非均勻載荷作用下功能梯度板的動(dòng)態(tài)SIFs進(jìn)行了研究.結(jié)果表明，此方法具有較高的精度，相比傳統(tǒng)伽遼金法可節(jié)省大量計(jì)算時(shí)間[8-9]，但計(jì)算效率仍不高.滕子浩等[10]對(duì)含有間斷的非均勻材料的斷裂問(wèn)題，將虛節(jié)點(diǎn)多邊形單元的形函數(shù)引入到擴(kuò)展有限元(extended finite element method，XFEM)中，提出了一種基于四叉樹(shù)結(jié)構(gòu)的動(dòng)態(tài)網(wǎng)格細(xì)化方法，并對(duì)非均質(zhì)材料中的裂紋擴(kuò)展問(wèn)題進(jìn)行研究，結(jié)果表明該方法相比傳統(tǒng)XFEM具有更好的精度、收斂性以及計(jì)算效率，但需對(duì)裂紋附近的網(wǎng)格進(jìn)行細(xì)化，增加了工作量.由此可見(jiàn)，F(xiàn)GM結(jié)構(gòu)的斷裂行為備受關(guān)注，而現(xiàn)行的分析方法還存在諸多問(wèn)題，半解析法應(yīng)用范圍有限，數(shù)值法雖有廣泛的適用性和較好的精度及計(jì)算效率，但尚有改進(jìn)空間.

本課題組研究建立的廣義參數(shù)Williams單元[11](簡(jiǎn)記為W單元)可直接求解裂尖SIFs，無(wú)需網(wǎng)格加密和后處理，精度很高，且已開(kāi)發(fā)了相應(yīng)的快速解法程序塊[12].目前，W單元的研究成果主要集中于均勻材料的斷裂分析，鑒于當(dāng)前FGM結(jié)構(gòu)的斷裂行為研究尚需深入，本文建立了帶裂紋FGM薄板I-II復(fù)合型裂紋尖端SIFs分析的W單元新格式，并分析了彈性模量E的變化形式、裂紋長(zhǎng)度及裂紋傾角對(duì)裂尖SIFs的影響，為帶裂紋的FGM薄板裂尖SIFs求解提供了新思路.

1 FGM薄板的廣義參數(shù)有限元模型

含中心斜裂紋均勻拉伸FGM薄板如圖1所示，板寬為W，高為H，厚為t，裂紋長(zhǎng)度為2a，傾斜角為γ，薄板上下邊界承受均布拉應(yīng)力σ0.以裂紋中心作為整體坐標(biāo)原點(diǎn)o，水平方向?yàn)閤軸，逆時(shí)針旋轉(zhuǎn)90°為y軸，建立直角坐標(biāo)系xoy.以任意裂尖為坐標(biāo)原點(diǎn)oi，沿裂紋擴(kuò)展方向?yàn)閤i軸，逆時(shí)針旋轉(zhuǎn)90°為yi軸，建立裂尖局部坐標(biāo)系xioiyi(圖中i取1或2).文獻(xiàn)[2]對(duì)含裂紋FGM薄板研究表明:泊松比μ對(duì)裂尖SIFs的影響較小，可忽略，因此本文研究取μ為常數(shù)，僅考慮彈性模量E隨坐標(biāo)軸分布，并假設(shè)其為整體坐標(biāo)的函數(shù).

因裂尖應(yīng)力場(chǎng)具有奇異性，將整個(gè)計(jì)算模型劃分為裂尖奇異區(qū)域和外圍常規(guī)區(qū)域.FGM矩形薄板常規(guī)區(qū)網(wǎng)格可利用ANSYS有限元軟件中四邊形8結(jié)點(diǎn)等參單元自動(dòng)離散，如圖2所示，因此類(lèi)單元中不同位置E值不同，稱(chēng)之為非均勻單元.裂尖奇異區(qū)網(wǎng)格離散如圖3所示，其離散過(guò)程如下:選取以裂尖oi為中心點(diǎn)，邊長(zhǎng)為l的正方形微小區(qū)域?yàn)槠娈悈^(qū)，將oi與各邊中點(diǎn)和角點(diǎn)連接形成8個(gè)相同三角形條元，利用一組環(huán)繞裂尖oi且相互平行的折線將每個(gè)條元?jiǎng)澐譃閚個(gè)梯形微單元和1個(gè)裂尖三角形微單元(紅色區(qū)域)，任意相鄰折線Γn，Γn-1到裂尖oi的距離之比為常數(shù)α(α<1)，α稱(chēng)為奇異區(qū)徑向離散比例因子，n稱(chēng)為奇異區(qū)徑向離散單元層數(shù).選取條元oiL1L2作為典型單元對(duì)其組成進(jìn)行詳細(xì)說(shuō)明:最外層梯形微單元L1L2L3L4的邊L1L2為常規(guī)區(qū)與奇異區(qū)的公共邊界，其上3個(gè)黑色實(shí)心圓點(diǎn)表示常規(guī)單元結(jié)點(diǎn)，另外的5個(gè)空心圓點(diǎn)位于奇異區(qū)內(nèi)，為虛結(jié)點(diǎn)，故稱(chēng)微單元L1L2L3L4為過(guò)渡單元；最內(nèi)層三角形微單元因面積過(guò)小，可忽略其剛度貢獻(xiàn)，因此，可用條元L1L2L5L6的剛度代替oiL1L2整個(gè)條元的剛度，因條元內(nèi)各微單元的位移場(chǎng)受Williams級(jí)數(shù)控制，故稱(chēng)之為W單元.

2 FGM薄板裂尖奇異區(qū)位移場(chǎng)與SIFs

裂紋尖端奇異區(qū)位移場(chǎng)采用Williams級(jí)數(shù)表示，并截取前m+1項(xiàng)，即

(1)

(2)

3 建立FGM薄板整體控制方程

3.1 FGM薄板常規(guī)區(qū)單元控制方程

在FGM薄板中，8結(jié)點(diǎn)等參單元各結(jié)點(diǎn)處的彈性模量值不相等，根據(jù)有限單元法可將單元位移函數(shù)表示成結(jié)點(diǎn)位移函數(shù)，將FGM薄板中常規(guī)區(qū)非均勻單元內(nèi)部任意位置處的彈性模量表示為結(jié)點(diǎn)彈性模量的函數(shù)，即

(3)

根據(jù)有限單元法，建立FGM薄板中常規(guī)區(qū)非均勻單元控制方程:

(4)

(5)

常規(guī)區(qū)中部分結(jié)點(diǎn)位于常規(guī)區(qū)與奇異區(qū)公共邊界，即公共結(jié)點(diǎn)；其余結(jié)點(diǎn)位于外圍常規(guī)區(qū)，則FGM薄板常規(guī)區(qū)單元整體控制方程可分塊表示:

(6)

3.2 FGM薄板裂尖奇異區(qū)W單元控制方程

以圖3中的條元L1L2L5L6為例，說(shuō)明i裂尖1個(gè)W單元控制方程的建立過(guò)程.

將四邊形8結(jié)點(diǎn)等參單元位移向量表示為

(7)

令W單元中所有子單元的8個(gè)結(jié)點(diǎn)按相同順序排列，則相鄰兩層梯形微單元內(nèi)層與外層相同編號(hào)結(jié)點(diǎn)到裂尖的距離之比為α，由各層子單元的相似性可知其極坐標(biāo)關(guān)系為.

(8)

根據(jù)廣義參數(shù)有限元法，結(jié)點(diǎn)的局部位移場(chǎng)服從整體位移場(chǎng)，由式(1)和式(8)可知W單元中第k層四邊形8結(jié)點(diǎn)等參單元結(jié)點(diǎn)位移向量為

(9)

φi=[ai,0,bi,0,ai,1,bi,1,…,ai,m,bi,m]T.

(10)

將式(10)帶入式(9)可得

(11)

將式(11)代入非均勻單元控制方程式(4)可得

(12)

(13)

Ki,sφi=fi,s.

(14)

因過(guò)渡單元部分結(jié)點(diǎn)位于公共邊界，有結(jié)點(diǎn)號(hào)，而內(nèi)部5個(gè)結(jié)點(diǎn)位于奇異區(qū)，為虛結(jié)點(diǎn)，其位移場(chǎng)需進(jìn)行廣義變換.根據(jù)結(jié)點(diǎn)所處區(qū)域不同，將過(guò)渡單元控制方程進(jìn)行分塊:

(15)

(16)

3.3 整體控制方程

同一裂尖，周?chē)蠾單元均可按照式(16)形成控制方程，則i裂尖8個(gè)W單元控制方程可集成為

(17)

對(duì)于含有λ個(gè)裂尖的FGM薄板，根據(jù)常規(guī)區(qū)控制方程式(6)與單個(gè)奇異區(qū)所有W單元控制方程式(17)，將所有局部坐標(biāo)下的裂尖奇異區(qū)控制方程轉(zhuǎn)換為整體坐標(biāo)后與常規(guī)區(qū)控制方程進(jìn)行集成，可得模型整體控制方程:

(18)

4 算例分析

4.1 中心斜裂紋

取一含中心斜裂紋的FGM矩形薄板，如圖1所示.板寬W=0.4 m，高H=0.6 m，厚度t取單位厚度，裂紋長(zhǎng)度2a=0.2 m，傾斜角為γ，薄板上下邊界承受均布拉應(yīng)力σ0，以裂紋中心點(diǎn)為坐標(biāo)原點(diǎn)o，水平方向?yàn)閤軸，逆時(shí)針旋轉(zhuǎn)90°為y軸，建立整體坐標(biāo)系xoy，泊松比μ取0.3.分別取彈性模量E呈指數(shù)型和線性型分布進(jìn)行分析:

(19)

式中:β1和β2為材料非均勻參數(shù)；E0表示坐標(biāo)原點(diǎn)處彈性模量.另外，板左、右邊界處的彈性模量分別用EL，ER表示，上、下邊界處彈性模量分別用ET，ED表示.

將FGM板進(jìn)行有限元網(wǎng)格離散，奇異區(qū)尺寸取0.02 m，每個(gè)奇異區(qū)離散為8個(gè)扇形條元，外圍常規(guī)區(qū)采用大型通用有限元分析軟件ANSYS離 散為684個(gè)四邊形8結(jié)點(diǎn)等參單元，整個(gè)計(jì)算模型共有2 172個(gè)結(jié)點(diǎn).計(jì)算結(jié)果如下:

1) 彈性模量E沿x軸呈指數(shù)型分布:取板的左邊界彈性模量EL=380 GPa，右邊界彈性模量ER=116 GPa，式(19)中β2=0，經(jīng)計(jì)算可得材料非均勻參數(shù)β1=-0.029，選取不同的裂紋傾角γ，研究其與ki,g的變化規(guī)律，此模型與文獻(xiàn)[13]一致，并比較其結(jié)果，如圖4所示.

文獻(xiàn)[13]通過(guò)大型通用有限元軟件ANSYS中的1/4奇異單元對(duì)模型進(jìn)行求解，由圖4可知:本文解與文獻(xiàn)[13]解吻合很好，證明了本文方法用于求解FGM裂紋問(wèn)題的正確性.當(dāng)彈性模量E沿x軸呈指數(shù)型分布時(shí)，左右裂尖的ki,Ⅰ及ki,Ⅱ變化趨勢(shì)一致；隨著裂紋傾角γ的增大，左右裂尖ki,Ⅰ均逐漸減小，且彈性模量大的一側(cè)ki,Ⅰ越大，ki,Ⅱ則均呈現(xiàn)先增大后減小的趨勢(shì)，左右裂尖值較接近.

2) 彈性模量E呈線性型分布:取板的左邊界彈性模量EL=380 GPa，整體坐標(biāo)原點(diǎn)處的彈性模量E0=E(0,0)=210 GPa，研究彈性模量E沿x軸呈線性型分布時(shí)對(duì)ki,g的影響規(guī)律，計(jì)算結(jié)果如圖5所示.

對(duì)比圖5與圖4的計(jì)算結(jié)果發(fā)現(xiàn):二者只在值的大小上有細(xì)微差別，變化趨勢(shì)一致.由此可知:彈性模量E沿x軸的任意單調(diào)變化形式均不改變ki,g的變化趨勢(shì)，且當(dāng)彈性模量E沿x軸呈指數(shù)型或線性型分布時(shí)，僅對(duì)ki,I有較為顯著的影響，且彈性模量較大一側(cè)ki,I較大，彈性模量較小一側(cè)ki,I較小.

3) 彈性模量E沿x，y軸同時(shí)呈指數(shù)型分布:取E(-W/2,-H/2)=380 GPa，E(W/2,H/2)=116 GPa，E0=210 GPa，經(jīng)計(jì)算可得材料非均勻參數(shù)β1=β2=-0.029.研究彈性模量E滿足式(19)，沿x，y軸同時(shí)呈指數(shù)型分布時(shí)，ki,g的變化規(guī)律，計(jì)算結(jié)果如圖6所示.

將圖4與圖6進(jìn)行比較可知:二者ki,Ⅰ的值相同，而ki,Ⅱ的值出現(xiàn)了較大差異，說(shuō)明沿y軸方向分布的彈性模量E改變時(shí)，主要對(duì)兩裂尖的ki,Ⅱ值產(chǎn)生影響.進(jìn)一步分析可知:當(dāng)彈性模量E分布形式呈單調(diào)變化且梯度與荷載方向平行或垂直時(shí)，分別使中心斜裂紋兩裂尖的ki,Ⅰ或ki,Ⅱ值產(chǎn)生差異.

4.2 邊界裂紋

取一含邊界裂紋的FGM矩形薄板，如圖7所示.板寬W=1 m，高H=2 m，厚度t為單位厚度，裂紋長(zhǎng)度為a，傾斜角度γ=0°，薄板上下邊界承受均布拉應(yīng)力σ0，以裂紋與左邊界的交點(diǎn)為坐標(biāo)原點(diǎn)o，水平方向?yàn)閤軸，逆時(shí)針旋轉(zhuǎn)90°為y軸，建立整體坐標(biāo)系xoy.泊松比μ=0.23，分別取E為指數(shù)型和線性型進(jìn)行分析，其表達(dá)式為

(20)

將FGM板進(jìn)行有限元網(wǎng)格離散，奇異區(qū)尺寸為0.025 m，奇異區(qū)離散為8個(gè)扇形條元，外圍常規(guī)區(qū)采用大型通用有限元軟件ANSYS離散為817個(gè)四邊形8結(jié)點(diǎn)等參單元，整個(gè)計(jì)算模型共有2 588個(gè)結(jié)點(diǎn).計(jì)算結(jié)果如圖7所示.

1)令E1=ER=dE(0, 0)，分別取d=0.1，0.2，0.5，2.0，5.0，10.0，通過(guò)設(shè)置不同的裂紋長(zhǎng)度a，研究不同裂紋長(zhǎng)度下彈性模量梯度對(duì)裂尖SIFs的影響，此模型與文獻(xiàn)[14-15]一致，并將計(jì)算結(jié)果分別與文獻(xiàn)[14-15]解進(jìn)行比較，如圖8及圖9所示.

文獻(xiàn)[14-15]分別采用奇異積分方程組及1/4奇異單元求解含邊界裂紋FGM薄板裂尖SIFs，兩種方法的計(jì)算精度均相對(duì)較高.由圖8、圖9可知:本文解與文獻(xiàn)[14-15]解吻合得很好，證明了本文方法用于求解邊界裂紋FGM薄板結(jié)果的正確性.當(dāng)彈性模量E沿x軸呈指數(shù)型或線性型分布時(shí)，kI的變化趨勢(shì)均一致；當(dāng)裂紋長(zhǎng)度a為0.4，0.6 m時(shí)，kI隨著E1/E0的增大而逐漸減??；而當(dāng)裂紋長(zhǎng)度a=0.2 m時(shí)，kI隨E1/E0的增大先增大后減小，此時(shí)是由于裂紋尺寸過(guò)小，裂尖過(guò)于靠近邊界，且E1/E0<0.5時(shí)梯度過(guò)大共同作用造成此現(xiàn)象.

2) 令E1=E(W,H/2)=dE(0, 0)，分別取d=0.1，0.2，0.5，2.0，5.0，10.0，研究當(dāng)彈性模量E沿x，y軸同時(shí)呈指數(shù)型分布時(shí)，由于彈性模量梯度與裂紋面不平行，導(dǎo)致kI，kII同時(shí)存在，本文計(jì)算結(jié)果如圖10所示.當(dāng)彈性模量E沿x，y軸同時(shí)呈指數(shù)型分布時(shí)，kI的變化趨勢(shì)與彈性模量E沿x軸呈指數(shù)型或線性型分布時(shí)一致，kI隨著彈性模量比E1/E0的增大而逐漸減小，但當(dāng)裂紋長(zhǎng)度a=0.2 m時(shí)，kI在E1/E0接近于0時(shí)，會(huì)隨著E1/E0的增大而增大；而kII則隨著彈性模量比E1/E0的增大而逐漸增大.

5 結(jié) 論

1) 結(jié)合材料屬性變化的W單元能有效地解決FGM薄板平面斷裂問(wèn)題，以裂尖SIFs作為基本未知量在剛度方程中直接求解，得到高精度結(jié)果.

2) FGM結(jié)構(gòu)中裂尖位置處彈性模量越大的一側(cè)裂尖SIFs越大.分布形式呈單調(diào)變化的彈性模量E，當(dāng)其梯度與荷載方向平行或垂直時(shí)，將使中心斜裂紋兩個(gè)裂尖的I型或II型SIFs值產(chǎn)生差異.

3) 彈性模量E沿水平方向呈線性型分布、指數(shù)型分布或水平及豎直雙向呈指數(shù)分布時(shí)，SIFs隨裂紋傾角γ的變化趨勢(shì)相同，即在含裂紋的FGM薄板斷裂分析中，當(dāng)彈性模量E分布形式呈單調(diào)變化時(shí)，裂尖SIFs隨裂紋傾角γ的變化規(guī)律不受彈性模量E分布形式的影響.