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

楊軍輝，雷勇軍，蒙上陽

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

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

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

1.1 加料裂尖單元

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

圖1 界面裂紋坐標系Fig.1 Interfacial crack coordinate

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

Dundurs參數為:

其中，κ為Kolosov常數，對于平面應力κ=(3-ν)/(1+ν)，對于平面應變κ=3-4ν，式(1)中其他參數表達式如下:

將裂尖漸進位移場加入到相應的常規單元位移場中，需要進行坐標變換，把裂尖直角坐標系x'o'y'中的位移ux'，uy'變換為整體坐標系xoy中的位移ux，uy，變換關系式如下:

其中，(xo'，yo')為裂尖在整體坐標系中的坐標。對于4節點四邊形等參元，在常規位移模式中加入界面裂紋裂尖漸進位移項后，可得:

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

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

其中，Nm(ξ，η)為常規等參元形函數，mk為單元節點數為第m個節點處的節點位移值θ)為角函數fij(r，θ)在第m個節點處的函數值，Kj為附加節點自由度，即界面裂紋應力強度因子。

1.2 過渡單元

在加料裂尖單元和常規單元之間采用過渡單元，以消除兩種單元因位移模式引起的位移不協調問題，從而保證有限元解得收斂，提高計算精度。在加料裂尖單元位移模式的基礎上，引入一個調整函數Z(ξ，η)來構造過渡單元內位移模式，即:

其中，調整函數Z(ξ，η)須滿足:Z(ξ，η)在過渡單元與加料裂尖單元交界邊上為1，在過渡單元與常規單元交界邊上為0，實現從加料裂尖單元到常規單元的協調過渡。采用相對簡單的線性調整函數:

在有限元分析過程中，需要根據過渡單元與加料裂尖單元的連接方式以及過渡單元的局部坐標系，選擇適當的表達式。例如，當過渡單元的邊ξ=1與加料裂尖單元連接時，滿足條件的調整函數Z(ξ，η)表達式為Z(ξ，η)=0.5(1+ξ)。

2 有限元方程

為了便于推導有限元方程，將加料單元、過渡單元位移模式統一寫成如下向量形式:

其中:N是常規單元形函數，是節點位移列向量;N k加料單元或過渡單元的附加形函數附加自由度向量。N k的具體形式如下:

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

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

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

其中，U為總體位移列陣，K為應力強度因子列陣，其他符號表達式如下:

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

3 加料單元在界面裂紋中的應用

為便于對比分析，本節不涉及物理單位。

3.1 中心界面裂紋

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

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

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

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

雙向受拉工況材料參數和載荷取值見表2，應力強度因子計算結果以及和文獻［7］對比結果見表3。

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

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

表1 單向受拉應力強度因子Tab.1 Stress intensity factor values of uniaxial tensile load

表3 雙向受拉應力強度因子Tab.3 Stress intensity factor values of biaxial tensile load

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

表4 K2對比結果Tab.4 Compare results of K2

3.2 斜界面裂紋

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

當E1=E2時，退化為單相材料矩形板單邊斜裂紋問題，單向拉伸即可產生Ⅰ，Ⅱ復合型應力強度因子，當θ=45°，a/b=0.5時，計算結果為K1=1.208，K2=-0.572，與文獻［11］結果一致。

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

為了便于比較分析，將K0變形為K0=F0σ為橫坐標，F0為縱坐標做出變化曲線。圖7是θ=45°，a/b=0.5時的有限元模型，采用4節點四邊形單元，單元數量為1250個，節點數量為1227個，裂紋傾角為67.5°以及裂紋長度為0.3b時單元數量分別為1255個和1327個。根據計算結果繪制的曲線如圖8～11所示。

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

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

由圖8～11可見，在各種不同裂紋長度和傾角下，隨著結合材料振蕩指數的增大，復合應力強度因子的模逐漸減小，而模態比的絕對值則逐漸增大，說明界面兩側材料的差異性增強了剪切效應，并且裂紋傾角越大，模態比曲線變化越明顯，表明界面導致的KⅡ變化越大，這與馬開平等［12］用半權函數法得到的結果是一致的。由圖10、圖11可見，界面裂紋傾角相同的情況下，裂紋長度越小，復合應力強度因子越小;并且與裂紋傾角相比，由裂紋長度變化而引起的剪切模態KⅡ變化較小。

圖8 不同裂紋傾角下SIF計算結果Fig.8 SIF for different crack angles

圖9 不同裂紋傾角下MR計算結果Fig.9 MR for different crack angles

圖10 不同裂紋長度下SIF計算結果Fig.10 SIF for different crack lengths

圖11 不同裂紋長度下MR計算結果Fig.11 MR for different crack lengths

4 結論

1)通過在常規單元位移模式中引入界面裂紋漸進位移場，得到了界面裂紋加料單元和過渡單元的位移模式，推導了界面裂紋加料有限元方程，對一般界面裂紋平面問題進行了計算，結果表明應力強度因子數值解的精度滿足工程要求。

2)應用加料單元和過渡單元后，即使采用低階線性單元，網格劃分規模相對較小的情況下，也能得到較為精確的應力強度因子數值解。一般而言，應力強度因子主要模態的精度要好于次要模態的精度，加料單元和過渡單元的配置方式對應力強度因子計算精度有一定的影響，但對于主要模態和復合應力強度因子精度的影響不大。

3)加料有限元方法可直接從有限元方程求解得到應力強度因子，不再需要通過單元節點位移或單元應力插值的方式獲取，是處理界面斷裂問題的一種快速、有效的計算方法。

References)

［1］許金泉.界面力學［M］.北京:科學出版社，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］段靜波，雷勇軍.線粘彈性材料中三維裂紋問題的加料有限元法［J］.國防科技大學學報，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］中國航空研究院.應力強度因子手冊［M］.北京:科學出版社，1981.China Aeronautical Establishment.Stress intensity factor handbook［M］.Beijing:Science Press，1981.(in Chinese)

［12］馬開平，柳春圖.雙材料界面裂紋平面問題的半權函數法［J］.應用數學和力學，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)