雙材料界面裂紋的加料有限元方法*

楊軍輝，雷勇軍，蒙上陽(yáng)

(1.國(guó)防科技大學(xué)航天科學(xué)與工程學(xué)院，湖南長(zhǎng)沙410073;2.中國(guó)人民解放軍63961部隊(duì)，北京100012)

由多相物質(zhì)組成的新材料在工程上的應(yīng)用日益增多，不同材料的宏觀結(jié)構(gòu)之間存在著界面，因此，工程結(jié)構(gòu)中的界面是廣泛存在的:如固體火箭發(fā)動(dòng)機(jī)殼體/絕熱層、絕熱層/襯層(包覆層)、襯層/推進(jìn)劑藥柱之間的界面，復(fù)合材料、多相材料中的異質(zhì)界面等。界面是各種組合結(jié)構(gòu)最薄弱的地方，往往容易產(chǎn)生裂紋，裂紋在一定條件下的不穩(wěn)定擴(kuò)展將導(dǎo)致組合結(jié)構(gòu)的解體，因此，界面問(wèn)題變得越來(lái)越重要［1］，而如何精確求解應(yīng)力強(qiáng)度因子等斷裂參量是界面斷裂問(wèn)題的重要內(nèi)容。計(jì)算界面裂紋應(yīng)力強(qiáng)度因子常用的方法主要有權(quán)函數(shù)法［2］、邊界元方法［3］、奇異有限元方法［4］和擴(kuò)展有限元法［5］等數(shù)值方法。加料有限元方法最早由Benzly［6］提出用于解決單相材料的斷裂力學(xué)問(wèn)題。1985年Chen［7］首先采用加料有限元法解決雙材料界面裂紋問(wèn)題，通過(guò)在12節(jié)點(diǎn)四邊形單元引入加料項(xiàng)構(gòu)造了加料界面裂紋元，但加料單元與常規(guī)單元之間沒有采用過(guò)渡單元。此后，Bayram［8］將加料單元與罰函數(shù)結(jié)合解決基于接觸模型的界面裂紋問(wèn)題，段靜波［9］等將加料單元應(yīng)用于粘彈性材料，研究了單相線粘彈性材料中的斷裂問(wèn)題。Chen等的方法較好地求解了雙材料界面裂紋的應(yīng)力強(qiáng)度因子，但工程應(yīng)用不便，且不能直接應(yīng)用于斜界面裂紋。

1 界面裂紋加料單元位移模式

1.1 加料裂尖單元

加料單元包括加料裂尖單元以及為了保證加料單元和常規(guī)單元位移協(xié)調(diào)一致而引入的過(guò)渡單元。圖1所示雙線彈性材料界面裂紋ν2分別表示材料1和材料2的彈性模量和泊松比)，在裂尖直角坐標(biāo)系x'o'y'中，裂紋尖端位移場(chǎng)表達(dá)式如下［7］:

圖1 界面裂紋坐標(biāo)系Fig.1 Interfacial crack coordinate

其中，ux'，uy'為x'和y'方向的位移，r，θ為裂尖局部極坐標(biāo)下的徑向和周向坐標(biāo)，下標(biāo)j表示兩種不同的材料(j=1，2)，μ為剪切模量，K1，K2為界面裂紋張開型和滑開型模態(tài)的應(yīng)力強(qiáng)度因子(Stress Intensity Factor，SIF)，ζ為結(jié)合材料的振蕩指數(shù)，定義如下:

Dundurs參數(shù)為:

其中，κ為Kolosov常數(shù)，對(duì)于平面應(yīng)力κ=(3-ν)/(1+ν)，對(duì)于平面應(yīng)變?chǔ)?3-4ν，式(1)中其他參數(shù)表達(dá)式如下:

將裂尖漸進(jìn)位移場(chǎng)加入到相應(yīng)的常規(guī)單元位移場(chǎng)中，需要進(jìn)行坐標(biāo)變換，把裂尖直角坐標(biāo)系x'o'y'中的位移ux'，uy'變換為整體坐標(biāo)系xoy中的位移ux，uy，變換關(guān)系式如下:

其中，(xo'，yo')為裂尖在整體坐標(biāo)系中的坐標(biāo)。對(duì)于4節(jié)點(diǎn)四邊形等參元，在常規(guī)位移模式中加入界面裂紋裂尖漸進(jìn)位移項(xiàng)后，可得:

其中，ui(i=1，2)分別表示總體坐標(biāo)系下加料界面裂紋單元x，y方向的位移，ξ，η分別為單元局部坐標(biāo)系的坐標(biāo)軸，fij(r，θ)為裂尖角函數(shù)，根據(jù)界面裂紋的裂尖漸進(jìn)位移場(chǎng)表達(dá)式，有:

將4節(jié)點(diǎn)四邊形等參元的節(jié)點(diǎn)局部坐標(biāo)(ξi，ηi)代入式(6)中，可得到廣義坐標(biāo)αij的表達(dá)式，再將得到的廣義坐標(biāo)回代到式(6)中，得到加料界面裂尖單元的位移模式如下:

其中，Nm(ξ，η)為常規(guī)等參元形函數(shù)，mk為單元節(jié)點(diǎn)數(shù)為第m個(gè)節(jié)點(diǎn)處的節(jié)點(diǎn)位移值θ)為角函數(shù)fij(r，θ)在第m個(gè)節(jié)點(diǎn)處的函數(shù)值，Kj為附加節(jié)點(diǎn)自由度，即界面裂紋應(yīng)力強(qiáng)度因子。

1.2 過(guò)渡單元

在加料裂尖單元和常規(guī)單元之間采用過(guò)渡單元，以消除兩種單元因位移模式引起的位移不協(xié)調(diào)問(wèn)題，從而保證有限元解得收斂，提高計(jì)算精度。在加料裂尖單元位移模式的基礎(chǔ)上，引入一個(gè)調(diào)整函數(shù)Z(ξ，η)來(lái)構(gòu)造過(guò)渡單元內(nèi)位移模式，即:

其中，調(diào)整函數(shù)Z(ξ，η)須滿足:Z(ξ，η)在過(guò)渡單元與加料裂尖單元交界邊上為1，在過(guò)渡單元與常規(guī)單元交界邊上為0，實(shí)現(xiàn)從加料裂尖單元到常規(guī)單元的協(xié)調(diào)過(guò)渡。采用相對(duì)簡(jiǎn)單的線性調(diào)整函數(shù):

在有限元分析過(guò)程中，需要根據(jù)過(guò)渡單元與加料裂尖單元的連接方式以及過(guò)渡單元的局部坐標(biāo)系，選擇適當(dāng)?shù)谋磉_(dá)式。例如，當(dāng)過(guò)渡單元的邊ξ=1與加料裂尖單元連接時(shí)，滿足條件的調(diào)整函數(shù)Z(ξ，η)表達(dá)式為Z(ξ，η)=0.5(1+ξ)。

2 有限元方程

為了便于推導(dǎo)有限元方程，將加料單元、過(guò)渡單元位移模式統(tǒng)一寫成如下向量形式:

其中:N是常規(guī)單元形函數(shù)，是節(jié)點(diǎn)位移列向量;N k加料單元或過(guò)渡單元的附加形函數(shù)附加自由度向量。N k的具體形式如下:

對(duì)于加料單元Z(ξ，η)≡1，過(guò)渡單元按式(9)取值，將位移向量式(11)代入位移應(yīng)變關(guān)系中，得到:

其中，B表示單元常規(guī)應(yīng)變矩陣，B k則是由于加料界面裂紋單元位移模式中引入裂紋尖端位移項(xiàng)而產(chǎn)生的附加項(xiàng)，稱之為加料裂紋元附加應(yīng)變矩陣。將式(13)代入單元應(yīng)力應(yīng)變關(guān)系式式中，應(yīng)用總體勢(shì)能泛函為:

其中，b為體力，P為面力，Ω表示求解域，Γ為面力積分邊界，D為材料矩陣，取值依據(jù)界面兩側(cè)材料具體參數(shù)確定，可得到:

其中，U為總體位移列陣，K為應(yīng)力強(qiáng)度因子列陣，其他符號(hào)表達(dá)式如下:

其中，ns表示加料單元數(shù)量，nt表示過(guò)渡單元數(shù)量，no表示常規(guī)單元數(shù)量，上標(biāo)e表示單元，根據(jù)最小勢(shì)能原理，式(15)分別對(duì)U，K變分，得到有限元方程如下:

3 加料單元在界面裂紋中的應(yīng)用

為便于對(duì)比分析，本節(jié)不涉及物理單位。

3.1 中心界面裂紋

對(duì)含中心界面裂紋的彈性體承受均布拉力的問(wèn)題，按平面應(yīng)變進(jìn)行了計(jì)算。其中，y軸單向受拉如圖2所示，x，y軸雙向受拉如圖3所示。裂紋長(zhǎng)度a=2，有限元建模時(shí)，取平板邊長(zhǎng)b=10a，近似滿足無(wú)限大條件。根據(jù)對(duì)稱性，取1/2模型建模，裂尖局部網(wǎng)格適當(dāng)加密(網(wǎng)格尺度l/a=1/40)，共劃分為800個(gè)4節(jié)點(diǎn)四邊形單元，850個(gè)節(jié)點(diǎn)。為了考察加料單元和過(guò)渡單元的配置區(qū)域?qū)τ?jì)算精度的影響，采用兩種配置方案:第1種方案是在上半平面(y〉0，對(duì)應(yīng)材料1)配置2個(gè)加料單元，4個(gè)過(guò)渡單元(如圖4所示);第2種方案是在上半平面配置6個(gè)加料單元，4個(gè)過(guò)渡單元(如圖5所示)，兩種方案下半平面(對(duì)應(yīng)材料2)均采用與上半平面同樣的配置。計(jì)算單元?jiǎng)偠染仃嚂r(shí)，加料單元和過(guò)渡單元用8×8高斯積分，常規(guī)單元采用2×2高斯積分。

圖2 單向受拉工況Fig.2 Uniaxial tensile load case

圖3 雙向受拉工況Fig.3 Biaxial tensile load case

圖4 加料單元和過(guò)渡單元配置方案1Fig.4 Scheme one of enriched element and transition element

雙向受拉工況材料參數(shù)和載荷取值見表2，應(yīng)力強(qiáng)度因子計(jì)算結(jié)果以及和文獻(xiàn)［7］對(duì)比結(jié)果見表3。

圖5 加料單元和過(guò)渡單元配置方案2Fig.5 Scheme two of enriched element and transition element

表2 材料參數(shù)及載荷工況Tab.2 Material parameters and load cases

表1 單向受拉應(yīng)力強(qiáng)度因子Tab.1 Stress intensity factor values of uniaxial tensile load

表3 雙向受拉應(yīng)力強(qiáng)度因子Tab.3 Stress intensity factor values of biaxial tensile load

文獻(xiàn)［7］加料單元和常規(guī)單元均為12節(jié)點(diǎn)四邊形單元，但加料單元和常規(guī)單元之間沒有采用過(guò)渡單元。由對(duì)比結(jié)果可見，在各種工況下，K1，K2計(jì)算值與文獻(xiàn)［7］相差均不超過(guò)5%，K0的相對(duì)誤差不超過(guò)1%。與K1相比，K2誤差相對(duì)較大，經(jīng)分析，文獻(xiàn)［7］K2計(jì)算值與準(zhǔn)確結(jié)果相比誤差較大，而加料有限元法K2計(jì)算值與準(zhǔn)確結(jié)果更接近。K2準(zhǔn)確值可參考文獻(xiàn)［10］，該文獻(xiàn)應(yīng)用17節(jié)點(diǎn)應(yīng)力雜交元，保留了應(yīng)力場(chǎng)漸進(jìn)解的前16項(xiàng)，其計(jì)算結(jié)果是相當(dāng)準(zhǔn)確的，與加料有限元法對(duì)比結(jié)果見表4。因此，盡管有限元模型沒有采用高階單元，但由于應(yīng)用了加料單元和過(guò)渡單元，并在裂紋局部網(wǎng)格適當(dāng)加密，仍達(dá)到了令人滿意的計(jì)算精度。以上兩個(gè)算例表明，有限元模型及計(jì)算方法對(duì)界面裂紋問(wèn)題是有效的，能方便地給出應(yīng)力強(qiáng)度因子的正確結(jié)果。

表4 K2對(duì)比結(jié)果Tab.4 Compare results of K2

3.2 斜界面裂紋

受均勻拉伸的具有斜單邊界面裂紋的矩形板如圖6所示，b=2，a/b=0.5，裂紋傾斜角為θ。材料1和材料2的泊松比固定為0.3，則結(jié)合材料參數(shù)只與彈性模量E1，E2有關(guān)，將E2固定為1，E1取不同值，可得到不同的結(jié)合材料參數(shù)，按平面應(yīng)變計(jì)算斜界面裂紋應(yīng)力強(qiáng)度因子。

當(dāng)E1=E2時(shí)，退化為單相材料矩形板單邊斜裂紋問(wèn)題，單向拉伸即可產(chǎn)生Ⅰ，Ⅱ復(fù)合型應(yīng)力強(qiáng)度因子，當(dāng)θ=45°，a/b=0.5時(shí)，計(jì)算結(jié)果為K1=1.208，K2=-0.572，與文獻(xiàn)［11］結(jié)果一致。

當(dāng)E1取值從1變化到1000時(shí)，振蕩指數(shù)變化區(qū)間為(0，0.0934)，分別計(jì)算了a/b=0.5，裂紋傾角分別為45°，67.5°以及裂紋傾角為45°，a/b分別等于0.3，0.5時(shí)K0和模態(tài)比(K2/K1，Modal Ratio，MR)的值。

為了便于比較分析，將K0變形為K0=F0σ為橫坐標(biāo)，F(xiàn)0為縱坐標(biāo)做出變化曲線。圖7是θ=45°，a/b=0.5時(shí)的有限元模型，采用4節(jié)點(diǎn)四邊形單元，單元數(shù)量為1250個(gè)，節(jié)點(diǎn)數(shù)量為1227個(gè)，裂紋傾角為67.5°以及裂紋長(zhǎng)度為0.3b時(shí)單元數(shù)量分別為1255個(gè)和1327個(gè)。根據(jù)計(jì)算結(jié)果繪制的曲線如圖8～11所示。

圖6 雙材料矩形板單邊斜界面裂紋Fig.6 Bimaterial rectangular plate unilateral oblique interface crack

圖7 斜界面裂紋加料有限元模型Fig.7 Unilateral oblique interface crack EFE model

由圖8～11可見，在各種不同裂紋長(zhǎng)度和傾角下，隨著結(jié)合材料振蕩指數(shù)的增大，復(fù)合應(yīng)力強(qiáng)度因子的模逐漸減小，而模態(tài)比的絕對(duì)值則逐漸增大，說(shuō)明界面兩側(cè)材料的差異性增強(qiáng)了剪切效應(yīng)，并且裂紋傾角越大，模態(tài)比曲線變化越明顯，表明界面導(dǎo)致的KⅡ變化越大，這與馬開平等［12］用半權(quán)函數(shù)法得到的結(jié)果是一致的。由圖10、圖11可見，界面裂紋傾角相同的情況下，裂紋長(zhǎng)度越小，復(fù)合應(yīng)力強(qiáng)度因子越小;并且與裂紋傾角相比，由裂紋長(zhǎng)度變化而引起的剪切模態(tài)KⅡ變化較小。

圖8 不同裂紋傾角下SIF計(jì)算結(jié)果Fig.8 SIF for different crack angles

圖9 不同裂紋傾角下MR計(jì)算結(jié)果Fig.9 MR for different crack angles

圖10 不同裂紋長(zhǎng)度下SIF計(jì)算結(jié)果Fig.10 SIF for different crack lengths

圖11 不同裂紋長(zhǎng)度下MR計(jì)算結(jié)果Fig.11 MR for different crack lengths

4 結(jié)論

1)通過(guò)在常規(guī)單元位移模式中引入界面裂紋漸進(jìn)位移場(chǎng)，得到了界面裂紋加料單元和過(guò)渡單元的位移模式，推導(dǎo)了界面裂紋加料有限元方程，對(duì)一般界面裂紋平面問(wèn)題進(jìn)行了計(jì)算，結(jié)果表明應(yīng)力強(qiáng)度因子數(shù)值解的精度滿足工程要求。

2)應(yīng)用加料單元和過(guò)渡單元后，即使采用低階線性單元，網(wǎng)格劃分規(guī)模相對(duì)較小的情況下，也能得到較為精確的應(yīng)力強(qiáng)度因子數(shù)值解。一般而言，應(yīng)力強(qiáng)度因子主要模態(tài)的精度要好于次要模態(tài)的精度，加料單元和過(guò)渡單元的配置方式對(duì)應(yīng)力強(qiáng)度因子計(jì)算精度有一定的影響，但對(duì)于主要模態(tài)和復(fù)合應(yīng)力強(qiáng)度因子精度的影響不大。

3)加料有限元方法可直接從有限元方程求解得到應(yīng)力強(qiáng)度因子，不再需要通過(guò)單元節(jié)點(diǎn)位移或單元應(yīng)力插值的方式獲取，是處理界面斷裂問(wèn)題的一種快速、有效的計(jì)算方法。

References)

［1］許金泉.界面力學(xué)［M］.北京:科學(xué)出版社，2006.XU Jinquan.The mechanics of interface［M］.Beijing:Science Press，2006.(in Chinese)

［2］Gao H J.Weight function method for interface crack in anisotropic bimaterials［J］.International Journal of Fracture，1992，56(2):139-158.

［3］He W J，Bolander J E J，Lin D S，et al.A boundary element for crack analysis at a bimaterial interface［J］.Engineering Fracture Mechanics，1994，49(3):405-410.

［4］Wu Y L.Variable power singular interface elements for a crack normal to the bimaterial interface［J］.Engineer Ingfracture Mechanics，1994，48(6):763-772.

［5］Sukumar N，Huang Z Y，Prévost J H，et al.Partition of unity enrichment for bimaterial interface cracks［J］.International Journal for Numerical Methods in Engineering，2004，59(8):1075-1102.

［6］Benzley S E.Representation of singularities with isoparametric finite elements［J］.International Journal for Numerical Methods in Engineering，1974，8(3):537-545.

［7］Chen E P.Finite element analysis of a bimaterial interface crack［J］.Theoretical and Applied Fracture Mechanics，1985，3(3):257-262.

［8］Bayram Y B，Nied H F.Enriched finite element-penalty function method for modeling interface cracks with contact［J］.Engineering Fracture Mechanics，2000，65(1):541-557.

［9］段靜波，雷勇軍.線粘彈性材料中三維裂紋問(wèn)題的加料有限元法［J］.國(guó)防科技大學(xué)學(xué)報(bào)，2012，34(3):6-11.DUAN Jingbo，LEI Yongjun.The enriched finite element method for 3-D fracture problems in viscoelastic materials［J］.Journal of National University of Defense Technology，2012，34(3):6-11.(in Chinese)

［10］Lin K Y，Mar J W.Finite element analysis of stress intensity factors for cracks at a bi-material interface［J］.International Journal of Fracture，1976，12(4):521-531.

［11］中國(guó)航空研究院.應(yīng)力強(qiáng)度因子手冊(cè)［M］.北京:科學(xué)出版社，1981.China Aeronautical Establishment.Stress intensity factor handbook［M］.Beijing:Science Press，1981.(in Chinese)

［12］馬開平，柳春圖.雙材料界面裂紋平面問(wèn)題的半權(quán)函數(shù)法［J］.應(yīng)用數(shù)學(xué)和力學(xué)，2004，25(11):1135-1142.MA Kaiping，LIU Chuntu.Semi-weight function method on computation of stress intensity factors in dissimilar materials［J］.Applied Mathematics and Mechanics，2004，25(11):1135-1142.(in Chinese)