雙材料V型切口應(yīng)力強度因子的加料有限元分析*

楊軍輝，韓珺禮，雷勇軍，蒙上陽

(1．國防科技大學 航天科學與工程學院， 湖南 長沙　410073； 2．北京特種機電技術(shù)研究所， 北京　100012)

?

雙材料V型切口應(yīng)力強度因子的加料有限元分析*

楊軍輝1,2，韓珺禮2，雷勇軍1，蒙上陽2

(1．國防科技大學 航天科學與工程學院， 湖南 長沙410073； 2．北京特種機電技術(shù)研究所， 北京100012)

摘要：應(yīng)用Williams本征函數(shù)展開和線性變換求解V型切口端部漸進位移場。將該位移場加入常規(guī)等參單元位移模式中，構(gòu)造雙材料V型切口加料單元和過渡單元的位移模式，推導加料有限元方程。建立帶V型缺口雙材料三點彎曲梁試件和直角界面端平面問題的加料有限元模型，求解有限元方程可直接得到應(yīng)力強度因子。計算結(jié)果與用其他方法得到的數(shù)據(jù)吻合，驗證了方法的正確性，可用于雙材料V型切口結(jié)構(gòu)斷裂特性計算分析。

關(guān)鍵詞：雙材料V型切口；漸進位移場；加料單元；過渡單元；應(yīng)力強度因子

隨著結(jié)合材料工程應(yīng)用范圍的不斷擴大，因結(jié)合材料界面缺陷導致的結(jié)構(gòu)失效問題日益引起人們的廣泛關(guān)注。界面端部在結(jié)合材料中不可避免，而界面端部由于材料幾何結(jié)構(gòu)或工藝制作上的原因出現(xiàn)的V型切口，則是一種常見的界面缺陷表現(xiàn)形式。

Pageau等[1]應(yīng)用Williams本征函數(shù)方法研究了結(jié)合材料V型切口和界面折點奇異性平面問題，給出了特征值為復數(shù)時位移和應(yīng)力場角函數(shù)求解方法。Tan等[2]發(fā)展了一種計算應(yīng)力強度因子的一維線性有限元方法。Gu等[3]給出了一種求解雙材料V型切口奇異性指數(shù)的有限元方法，該方法通過設(shè)定切口端點漸進位移場或應(yīng)力函數(shù)的分離變量解，應(yīng)用變分原理將求解的問題轉(zhuǎn)化為一個特征值問題，求解特征值問題可直接得到奇異性指數(shù)。隨后，Pageau等[4]用類似的方法，通過給定極坐標下奇異位移場的位移模式，利用虛功原理將求解奇異性指數(shù)和角函數(shù)轉(zhuǎn)化為一個特征值問題，求解特征方程可直接得到奇異性指數(shù)和位移角函數(shù)值，并將該方法推廣到三維雙材料V型切口問題[5]，Zhang等[6]用同樣的方法提出了一種適用于冪硬化材料界面切口端點奇異場的非線性有限元方法。在上述研究基礎(chǔ)上，平學成等[7-8]提出了基于位移的V型切口非協(xié)調(diào)元法，程長征等[9]用邊界元方法對雙材料V型切口應(yīng)力奇異性進行了分析。王海濤[10]提出了一種分析不同材料V型切口應(yīng)力奇異性的一維雜交有限元方法，通過扇形區(qū)域在角度方向上離散得到特征矩陣方程。

上述方法多以一維有限元特征矩陣計算特征值，應(yīng)力強度因子的獲取以應(yīng)力雜交元和邊界元方法為主，在計算網(wǎng)格前處理和擴展性方面，不便于推廣到主流有限元程序中。為此，楊軍輝等提出了一種基于加料有限元方法的雙材料V型切口應(yīng)力奇異性分析方法。目前，加料有限元方法在單相介質(zhì)斷裂[11]和界面裂紋問題中已得到應(yīng)用[12]，但在雙材料V型切口問題中還未見相關(guān)文獻報道。

1雙材料V型切口漸進位移場

1.1雙材料V型切口特征值

雙材料V型切口端點O(如圖1所示)附近的應(yīng)力和位移場可用Williams本征函數(shù)展開表示，在極坐標系中，Williams給出的應(yīng)力函數(shù)為：

Φ(r,θ)=rλ+1f(θ)

(1)

圖1　雙材料V型切口Fig.1　Bi-material V-notch

Φ(r,θ)滿足雙調(diào)和方程▽4Φ(r,θ)=0，由此可得

f(θ)=VRVck

(2)

其中:VR=[sin(λ+1)θcos(λ+1)θsin(λ-1)θcos(λ-1)θ]，Vck=Q[Ck,1Ck,2Ck,3Ck,4]T，Ck,1、Ck,2、Ck,3、Ck,4為常系數(shù)項，由邊界條件決定，λ為特征值，下標k代表第k種材料，Q為標準化系數(shù)。相應(yīng)的應(yīng)力和位移場可表示為：

(3)

將位移表示成矩陣向量的形式。

(4)

其中：μk為剪切模量；平面應(yīng)變狀態(tài)下κk=νk，平面應(yīng)力狀態(tài)下κk=νk/(1+νk)，νk為泊松比。在極坐標系下，應(yīng)力自由邊界條件和界面連續(xù)條件為：

(5)

將式(3)應(yīng)力、位移表達式代入邊界條件(5)中，得到線性齊次特征方程。

MVc=0

(6)

其中：M為各向同性彈性材料界面端V型切口特征矩陣，它與特征值λ、剪切模量μ、切口角θ1、θ2有關(guān)；Vc=[Vc1Vc2]T為解向量。式(6)有非零解的充要條件是特征矩陣M的行列式等于零，即

detM=0

(7)

式(7)即為各向同性彈性材料界面端V型切口問題的特征方程，求解該方程可得到特征值λ。圖2給出了平面應(yīng)變狀態(tài)下，E1=1.0，ν1=0.3，E2=10.0，ν2=0.3，θ1從90°變化到180°(θ1=-θ2)特征值變化規(guī)律。從計算結(jié)果可見，依據(jù)界面端切口材料參數(shù)、切口張開角的不同，在0

圖2　特征值隨切口角度的變化規(guī)律Fig.2　Variation of eigenvalue with V-notch surface angle

1.2特征值為實數(shù)時位移場角函數(shù)

當有兩個實特征根時，漸進應(yīng)力、位移場可表示為：

(8)

當特征值λ1，λ2為實數(shù)時，其對應(yīng)的應(yīng)力強度因子K1，K2是解耦的，根據(jù)通常的V型切口應(yīng)力強度因子定義，可得：

(9)

將應(yīng)力表達式(3)、解向量Vc代入式(9)中，可求得標準化系數(shù)Q與應(yīng)力強度因子Kk之間的關(guān)系。

(10)

式(10)表示解向量標準化系數(shù)和應(yīng)力強度因子之間的關(guān)系，Q1，Q2為：

(11)

將式(10)代入式(4)中，與式(8)比較后可得到極坐標下位移場角函數(shù)表達式。

(12)

圖3　雙材料V型切口位移場角函數(shù)(實特征值)Fig.3　Displacement field angular functions of bi-material V-notch (real roots)

圖3為E1=1.0，ν1=0.3，E2=10.0，ν2=0.3，θ1=-θ2=140°時平面應(yīng)變狀態(tài)下位移場角函數(shù)曲線，此時特征值λ1=0.578 4，λ2=0.771 1。退化為均質(zhì)材料V型切口時(E1=1.0，ν1=0.2，E2=E1，ν2=ν1)位移場角函數(shù)曲線如圖4所示，計算結(jié)果和解析解[13]一致，表明本文方法得到的位移角函數(shù)是正確的。

圖4　均質(zhì)材料V型切口位移場角函數(shù)Fig.4　Displacement field angular functions of homogeneous material V-notch

1.3特征值為復數(shù)時的位移場角函數(shù)

當特征方程(7)有一對共軛復根時，應(yīng)力強度因子K1，K2之間是耦合在一起的，此時界面端V型切口應(yīng)力強度因子定義為：

(13)

應(yīng)力和位移表示[1]為：

(14a)

(14b)

其中：β+iζ=λ，i為虛數(shù)單位；φ，α為標準化參數(shù)；角函數(shù)hijk(θ)，gik(θ)由特征值、特征方程基礎(chǔ)解向量和應(yīng)力、位移表達式得到。圖5給出了E1=1.0，ν1=0.3，E2=10.0，ν2=0.3，θ1=-θ2=30°，平面應(yīng)變狀態(tài)下特征值λ=0.548 5±0.072 4i時位移場角函數(shù)曲線。當θ1=-θ2=180°時，雙材料V型切口退化為界面裂紋，文獻[14]給出了直角坐標系下界面裂紋尖端位移場角函數(shù)解析解，通過坐標變換可得到極坐標系下的解析解。由圖6可見本文解與文獻[14]解析解一致，表明本文的結(jié)果是正確的。

圖5　雙材料V型切口位移場角函數(shù)(復特征值)Fig.5　Displacement field angular functions of bi-material V-notch (complex roots)

圖6　雙材料界面裂紋尖端位移場角函數(shù)Fig.6　Displacement field angular functions of bi-material interfacial crack

對于雙材料V型界面切口，當特征值為復數(shù)時，與雙材料界面裂紋相似，應(yīng)力和位移場角函數(shù)中仍包含riζ項，因此，其理論解同樣會出現(xiàn)應(yīng)力振蕩與位移相互嵌入的不合理現(xiàn)象。作為V型界面切口的一個特例，Shih等[15]和Rice[16]對界面裂紋中這一現(xiàn)象進行研究，許金泉[17]對其做了比較全面的總結(jié)分析，此處不再贅述。

2加料有限元方程

2.1加料切口單元位移函數(shù)

以位移表達式(8)為例，對其進行坐標變換，得到直角坐標系o′x′y′下的位移場，與位移式(14a)除表達形式略有不同外，變換過程完全相同。

(15)

其中，fij(r,θ)(i,j=1,2)表示直角坐標系下的切口尖端角函數(shù)，形式為：

(16)

對8節(jié)點四邊形等參元，在常規(guī)單元位移模式中加入式(15)后，可得加料切口單元位移模式。

(17)

其中，ui(i=1,2)分別表示總體坐標系下加料切口單元x，y方向的位移，αij為廣義坐標，ξ，η分別為單元局部坐標系的坐標軸，fij(r,θ)如式(16)所示。將8節(jié)點四邊形等參元的節(jié)點局部坐標(ξi,ηi)代入式(17)中，可得到廣義坐標αij的表達式，再將得到的廣義坐標回代到式(17)中，得到以節(jié)點位移和應(yīng)力強度因子表示的加料切口單元位移函數(shù)。

(18)

2.2過渡單元位移函數(shù)

在加料切口單元位移模式的基礎(chǔ)上，引入一個調(diào)整函數(shù)Z(ξ,η)來構(gòu)造過渡單元，以消除加料單元和常規(guī)單元之間因位移模式引起的位移不協(xié)調(diào)問題，從而保證有限元解的收斂，提高計算精度。

(19)

式(19)中調(diào)整函數(shù)Z(ξ,η)須滿足：Z(ξ,η)在過渡單元與加料切口單元交界邊上為1，在過渡單元與常規(guī)單元交界邊上為0，實現(xiàn)從加料切口單元到常規(guī)單元的協(xié)調(diào)過渡。采用相對簡單的線性調(diào)整函數(shù)。

(20)

在有限元分析過程中，需要根據(jù)過渡單元與加料切口單元的連接方式以及過渡單元的局部坐標系，選擇適當?shù)谋磉_式。

2.3有限元方程

為了便于推導有限元方程，將加料單元、過渡單元位移模式統(tǒng)一寫為：

(21)

(22)

對于加料單元Z(ξ,η)≡1，過渡單元按式(17)取值，將位移向量式(21)代入位移應(yīng)變關(guān)系ε=Lu中，得到：

(23)

其中：B表示單元常規(guī)應(yīng)變矩陣；Bk則是由于加料界面切口單元位移模式中引入切口尖端位移項而產(chǎn)生的附加項，稱之為附加應(yīng)變矩陣。將式(23)代入單元應(yīng)力應(yīng)變關(guān)系式σ=Dε中，應(yīng)用總體勢能泛函，得到：

(24)

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

(25)

其中：U為總體位移列陣，K為應(yīng)力強度因子列陣。其他符號表達式為：

(26)

其中：下標ns表示加料單元數(shù)量，nt表示過渡單元數(shù)量，no表示常規(guī)單元數(shù)量，上標e表示單元。根據(jù)最小勢能原理，式(25)分別對U，K變分，得到有限元方程為：

(27)

3算例分析

3.1算例1

具有V型切口的三點彎曲梁試件，由兩種材料構(gòu)成，如圖7所示。試件厚度B=1 mm，切口深度d=5 mm,寬度w=10 mm，長度H=4w, 切口角θ1=-θ2=135°，集中力P=1 N。構(gòu)成試件的兩種材料具有相同的泊松比，ν1=ν2=0.3，彈性模量E2=1 MPa,兩種材料彈性模量比E1/E2分別為1，3，5，7，10，按平面應(yīng)力問題計算不同彈性模量比下的應(yīng)力強度因子。

圖7　算例1模型和有限元網(wǎng)格示意圖Fig.7　Schematic of geometry and mesh for example 1

有限元模型共劃分為856個8節(jié)點四邊形單元，2713個節(jié)點，切口端點加料單元尺度為5.4×10-2w，在切口尖端設(shè)置4個加料單元、8個過渡單元。有限元網(wǎng)格劃分如圖7所示，切口尖端加料單元和過渡單元配置如圖8所示。計算單元剛度矩陣時，加料單元和過渡單元采用10×10高斯積分，常規(guī)單元采用3×3高斯積分。

圖8　算例1網(wǎng)格劃分和加料單元配置方案Fig.8　Meshes and enriched element scheme of example 1

特征值計算結(jié)果見表1，由表1可見各工況下兩個特征值均為實數(shù)，隨著E1/E2比值的增大，λ1逐漸增大而λ2逐漸減小。由表2應(yīng)力強度因子計算結(jié)果可知，隨著界面兩側(cè)材料差異性的增強，兩種模態(tài)的應(yīng)力強度因子的絕對值是增大的。

表1　算例1特征值計算結(jié)果

表2算例1應(yīng)力強度因子計算結(jié)果

Tab.2Stress intensity factors of example 1

E1/E2F1F2本文解文獻[9]解本文解文獻[9]解12.1402.10100032.44382.3935-0.56581-0.647052.83912.7857-0.92665-1.008973.24113.1793-1.2425-1.2938103.87663.7993-1.7044-1.6886

文獻[9]采用邊界元方法對相同的問題進行了計算，得到了應(yīng)力強度因子數(shù)值解，并標準化處理為Fi=KiBwλi/(6P)。為了便于對比，本文應(yīng)力強度因子采用同樣的處理方法，各種工況下應(yīng)力強度因子計算結(jié)果列于表2?？梢姼鞴r下本文結(jié)果與邊界元方法計算結(jié)果吻合，表明本文方法的正確性與有效性。

3.2算例2

由兩種彈性材料構(gòu)成的直角界面端平板，平板高為2H，寬W，受單向均勻分布拉伸載荷σ0作用，如圖9所示。為便于和相關(guān)文獻結(jié)果對比，結(jié)構(gòu)、載荷和材料屬性與文獻[18]一致，取H=16 mm，W=8 mm,選用鉛、軋制鋅、碳鋼、軋制純銅、鋁、灰口鑄鐵和白口鑄鐵等材料形成5種材料組合工況，材料屬性見表3。

圖9　算例2模型和有限元網(wǎng)格示意圖Fig.9　Schematic of geometry and mesh for example 2

根據(jù)載荷和幾何對稱性，取1/2建立有限元模型，共劃分為720個8節(jié)點四邊形單元，2287節(jié)點，界面端部設(shè)置2個加料單元，6個過渡單元,單元尺度為1.5×10-2w,界面端部加料單元和過渡單元配置如圖9所示。計算單元剛度矩陣時，加料單元和過渡單元采用10×10高斯積分，常規(guī)單元采用3×3高斯積分。

表3　材料參數(shù)

文獻[18]采用邊界元方法，通過應(yīng)力外插求取應(yīng)力強度因子，并對應(yīng)力強度因子做無量綱處理F=K/(σ0w1-λ)。本文應(yīng)力強度因子計算結(jié)果以及和文獻[15]對比情況見表4，由表4可知，應(yīng)力強度因子隨著特征值增大而增大，兩種方法得到的結(jié)果變化趨勢是一致的，并且數(shù)值吻合良好，差別很小，表明本文的計算方法是正確有效的。

表4　算例2應(yīng)力強度因子計算結(jié)果

4結(jié)論

應(yīng)用Willams本征函數(shù)展開和線性變換相結(jié)合的方法得到雙材料V型切口漸進位移場。與復勢函數(shù)法等方法相比，該方法推導過程及最終表達形式采用矩陣向量的形式，相對而言較為簡潔直觀，適于數(shù)值計算。

將雙材料V型切口漸進位移場加入到常規(guī)等參元位移模式中構(gòu)造了加料單元和過渡單元的位移表達式，推導了加料有限元方程，求解有限元方程得到應(yīng)力強度因子。通過帶V型缺口的雙材料三點彎曲梁試件和直角界面端平板受拉兩個算例，表明本文方法的正確性。該方法不僅能得到雙材料V型切口端點附近的應(yīng)力場和位移場角函數(shù)值，并且能通過求解加料有限元方程直接得到廣義應(yīng)力強度因子，避免了外插法需要對應(yīng)力場二次處理才能得到應(yīng)力強度因子的不便以及由此帶來的精度損失，是對雙材料V型切口進行應(yīng)力奇異性分析的一種有效方法。

參考文獻(References)

[1]Pageau S S, Gadi K S,Biggers S B Jr,et al. Standardized complex and logarithmic eigensolutions forn-material wedges and junctions[J]. International Journal of Fracture, 1996, 77(1):51-76.

[2]Tan M A, Meguid S A.Analysis of bimaterial wedges using a new singular finite element[J]. International Journal of Fracture, 1997, 88(4):373-391.

[3]Gu L, Belytschko T.A numerical study of stress singularities in a two-material wedge [J]. International Journal of Solids and Structures, 1994, 31(6):865-889.

[4]Pageau S S, Joseph P F,Biggers S B Jr.Finite element analysis of anisotropic materials with singular inplane stress fields [J].International Journal of Solids and Structures, 1995, 32(5):571-591.

[5]Pageau S S, Biggers S B Jr.A finite element approach to three dimensional singular stress states in anisotropic multi-material wedges and junctions [J].International Journal of Solids and Structures, 1996, 33(1):33-47.

[6]Zhang N S, Joseph P F.A nonlinear finite element eigen analysis of singular plane stress fields in bimaterial wedges including complex eigenvalues[J]. International Journal of Fracture, 1998, 90(3):175-207.

[7]平學成, 謝基龍, 陳夢成, 等.各向異性兩相材料尖劈奇性場的非協(xié)調(diào)元分析[J]. 力學學報, 2005, 37(1): 24-31.

PING Xuecheng, XIE Jilong, CHEN Mengcheng, et al.Anon-confirming finite element analysis of singular fields in prismatic anisotropic bimaterial wedges[J]. Acta Mechanica Sinica, 2005, 37(1): 24-31.(in Chinese)

[8]平學成, 陳夢成,謝基龍.復合材料尖劈和接頭端部奇性場的反平面問題研究[J]. 固體力學學報, 2005, 26(2):193-198.

PING Xuecheng, CHEN Mengcheng,XIE Jilong. Singular stress states in tips of anisotropic multi-material wedges and junction subjected to antiplane shear[J]. Acta Mechanica Solid Sinica, 2005, 26(2):193-198.(in Chinese)

[9]Cheng C Z, Niu Z R, Recho N.Analysis of the stress singularity for a bi-material V-notch by the boundary element method[J]. Applied Mathematical Modelling, 2013, 37(22): 9398-9408.

[10]王海濤.分析不同材料界面應(yīng)力奇異性的一維雜交有限元方法[J]. 工程力學, 2009, 26(2):21-26.

WANG Haitao.A one-dimensional hybrid finite element method for the analysis of stress singularities at bimaterial interface [J]. Engineering Mechanics,2009, 26(2):21-26.(in Chinese)

[11]段靜波,雷勇軍.線粘彈性材料中三維裂紋問題的加料有限元法[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)

[12]Ayhan A O, Kaya A C, Nied H F.Analysis of three-dimensional interface cracks using enriched finite elements[J]. International Journal of Fracture, 2006, 142(3):255-276.

[13]Carpenter W C.Mode I and Mode II stress intensities for plates with cracks of finite opening[J]. International Journal of Fracture, 1984, 26(3): 201-214.

[14]Chen E P.Finite element analysis of a bimaterial interface crack[J]. Theoretical and Applied Fracture Mechanics, 1985, 3(3): 257-262.

[15]Shih C F,Asaro R J.Elastic-plastic analysis of cracks in biomaterial interface[J].Journal of Applied Mechanics, 1988, 55(2): 299-316.

[16]Rice J R.Elastic fracture mechanics concepts for interfacial cracks[J]. Journal of Applied Mechanics,1988, 55(1): 98-103.

[17]許金泉.界面力學[M].北京:科學出版社,2006.

XU Jinquan.The mechanics of interface[M]. Beijing:Science Press, 2006. (in Chinese)

[18]唐亮,許金泉.直角結(jié)合異材界面端應(yīng)力強度系數(shù)的經(jīng)驗公式[J]. 力學季刊, 2005, 26(1): 96-101.

TANG Liang, XU Jinquan.An empirical formula for stress intensity coefficient of orthogonal bonded materials near interface end[J]. Chinese Quarterly of Mechanics,2005, 26(1): 96-101.(in Chinese)

Enriched finite element analysis of stress intensity factors of bi-material V-notch

YANGJunhui1,2,HANJunli2,LEIYongjun1,MENGShangyang2

(1.College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China；2.Beijing Institute of Special Electromechanical Technology, Beijing 100012, China)

Abstract：The V-notch asymptotic displacement field was derived through an approach based on the Williams’ series expansion and linear algebraic transforms. By incorporating the displacement expressions to the common isoparametric elements, the enriched and transition element displacement model were obtained, and then the enriched finite element equation was derived consequently. The enriched finite element model for a V-notched bi-material three-point bending beam and an orthogonal bonded materials interface end plane problem were constructed. The stress intensity factors can be solved directly from the finite element equation. Comparisons between the results and the published data computed with other algorithm indicate that the present method is correct and can be used to analyze the fracture property of the V-notched bi-material structure.

Key words：bi-material V-notch; asymptotic displacement field; enriched element; transition element; stress intensity factor

中圖分類號：O346.1

文獻標志碼：A

文章編號：1001-2486(2016)01-156-07

作者簡介：楊軍輝(1979—)，男，河北深澤人，博士研究生，E-mail：yangjun_hui@126.com;雷勇軍(通信作者)，男，教授，博士，博士生導師，E-mail：leiyj108@nudt.edu.cn

基金項目：國家自然科學基金資助項目(11272348)

*收稿日期：2015-02-02

doi:10.11887/j.cn.201601025

http://journal.nudt.edu.cn