基于擴展二相雜交應力有限元法的兩相材料界面裂紋斷裂力學分析

江海濤, 李歡, 楊文彩

(云南農業大學機電工程學院, 昆明650201)

隨著科學技術的發展,層狀復合材料或結構如復合材料層合板,因其比強度及比模量高、抗疲勞斷裂性能好等良好的力學性能在航空和電子等一些重要工業領域得到廣泛的應用。實際應用過程中,由于復合材料界面兩側各相材料性能各異,以及制造過程中存在的缺陷,在外載荷的作用下,結構斷裂產生的裂紋通常發生在界面處,界面裂紋的擴展也是降低復合材料性能的主要因素之一。例如,基巖和混凝土壩之間在外力的作用下,界面處產生大量初始微裂紋,進而擴展相交形成宏觀裂紋,最后擴展貫通導致結構失效。因此,對于界面裂紋的分析研究對優化復合材料,提高其應用性能及可靠性是十分重要的。

研究裂紋擴展的方法分為理論、實驗以及數值模擬,對于實際環境中材料的應用,為更準確地觀察測量裂紋的產生和裂紋尖端附近應力場以及應變場,在邊界條件和幾何條件下數值計算方法更加有效,具有更高的適應性和精確性。有限元法是數值計算中最可靠、應用最普遍的方法。在有限元方法中,對于裂紋的研究依賴于計算網格的劃分,在裂紋擴展過程中,計算網格要與裂紋的幾何形狀一致,并隨裂紋的每一步擴展重新劃分網格。Williams[1]研究了各向同性雙材料界面裂紋,結果表明界面裂紋尖端附近位移場和應力場具有振蕩奇異性。由于其振蕩行為,對各向同性雙材料界面裂紋進行建模采用針對均勻材料的傳統求解程序已經不充分,即使在裂紋尖端附近使用非常細的網格來捕捉界面附近的高梯度和奇異性的應力應變場,單元公式中沒有包含高階奇異項,精度沒有太大提高,導致收斂性極差。

為解決傳統有限元方法中的問題,Gao等[2]采用擴展有限元法(extended finite element method, XFEM)研究了各向異性斷裂韌性材料中的裂紋擴展行為。該方法在研究復合材料層合板時,成功地用于預測復合材料中典型的復雜裂紋場景。例如,分層模擬[3-6],微尺度模擬[7],多尺度模擬[8]和基體裂紋/分層相互作用[9-10]。在數值分析中裂紋面不必再與單元界面重合,但由于裂紋尖端附近節點和自由度的增多,計算成本提高,而且基于位移的擴展有限元法其單元模式還是低階位移單元,收斂性還是不好。1960年初期,對于固體結構的有限元分析都是以單元節點位移為未知數,根據勢能原理協調元或余能原理平衡元構建剛度矩陣。Pian[11]根據余能原理假定單元內部平衡的應力和相鄰單元間協調的邊界條件創建的一個單元定義為雜交元。雜交元分為雜交應力元和雜交應變元,雜交應力元根據位移和應力變分原理進行推導,最后求解是以節點位移為未知數的有限元。Voronoi單元有限元法(voronoi cell finite element method, VCFEM)是Ghosh等[12]基于Pian[11]提出的假設應力雜交元法建立的單元格式。Guo等[13]提出了一種新的考慮界面裂紋和基體裂紋的Voronoi單元,在文獻[14]的基礎上新加入了裂紋裂尖附近奇異應力場,準確地刻畫了裂尖周圍的應力奇異,并計算出準確的裂尖應力強度因子,而且裂紋的不連續性在單元中能被準確描述,并構建反映裂紋從界面裂紋向基體擴展和裂紋的跨單元的擴展的網格重劃分算法,模擬了含大量微結構和裂紋的顆粒增強復合材料的界面裂紋和基體裂紋的擴展貫穿過程[15]。徐越等[16]以擴展有限元法的仿真結果為訓練數據,葉片裂紋位置和裂紋長度為輸入參數,建立了裂紋尖端應力強度因子的多層反向傳播人工神經網絡(back propagation artificial neural network, BPANN)。

現提出一種擴展二相雜交應力有限元法,通過構建包含邊界裂紋、水平中心裂紋以及傾斜裂紋的兩相材料的擴展二相雜交應力有限元單元模型對界面裂紋進行斷裂力學分析。通過該方法得到的數值模擬結果和傳統有限元分析的數值結果進行比較,驗證該方法的有效性和準確性。

1 界面裂紋尖端附近應力場

如圖1所示,兩種材料的楊氏模量、泊松比和剪切模量分別為Ei、υi、ui(i=1,2),界面兩側材料中各存在一點Q1、Q2。在以界面裂紋尖端為坐標原點的局部坐標系(X1O1Y1)和(X2O2Y2)下極坐標分別為(L1,θ1)和(L2,θ2)。上層材料極角θ1在[0,π]范圍內沿X1軸正方向逆時針測量,極角θ2在[-π,0]范圍內沿X2軸正方向順時針測量。下層材料中極角θ1在[-π,0]范圍內沿X1軸正方向順時針測量,極角θ2在[0,π]范圍內沿X2軸正方向逆時針測量。包含水平中心界面裂紋的兩相材料的界面裂紋近端應力σ場可表示為

圖1 裂紋尖端附近應力場Fig.1 Stress field near the crack tip

(1)

Liε=cos(εlnL)+isin(εlnL)表現出應力場的振蕩奇異性;ε為振蕩參數或雙材料系數,其表達式為

(2)

式(2)中:βd為連接界面裂紋兩側彈性材料參數的第二個Dundur參數,第一個Dundur參數為αd。兩個Dundur參數的計算公式為

(3)

式(3)中:μ1、μ2為彈性常數;k1、k2為Kolosov常數,計算式為

(4)

從式(1)、式(2)中可知,第二個Dundur參數對界面裂紋尖端應力場分布的影響較大,似乎不受第一個Dundur參數的影響。如圖2所示為不同泊松

圖2 不同泊松比下雙材料Dundur參數βd 隨彈性模量比變化曲線Fig.2 Dundur parameterβdchanges curve with elastic modulus ratio under different Poisson’s ratios

比下第二個Dundur參數隨楊氏模量比的變化規律。由圖2可知,當兩種材料的彈性常數相同時,第二個Dundur參數βd=0。另外,從圖2可以看出,βd可能取正值,也可能取負值,這取決于(Ei,υi)的匹配。當兩種材料的泊松比取值相同時,βd取正數表示上層比下層軟,取負數表示下層比上層軟。

(5)

(6)

2 描述包含界面裂紋的兩相材料的雜交應力有限元公式

以圖3所示包含界面中心裂紋單元Ωe為例,其由上層材料Ω1和下層材料Ω2組成。即Ωe=Ω1∪Ω2單元外邊界由四種類型:給定位移邊界?Ωr、單元

圖3 包含中心界面裂紋的兩相材料單元模型特征邊界Fig.3 Characteristic boundary of two-phase material element model with central interfacial crack

通過最小余能原理,在滿足平衡方程和邊界條件情況下,余能泛函表達式為

(7)

(8)

利用Lagrange乘子將以上幾種約束條件引入最小余能泛函,得到修正最小余能泛函,表達式為

(9)

式(9)中:下標1指上層第一相材料,下標2指下層第二相材料,兩相材料界面連接處變量下標為12;u為位移;n為邊界外法線;σ1、σ2分別為兩相材料內部平衡應力場。雜交應力單元模型中應力場的構建選取對X-THSFEM的收斂性和效率有重要影響。X-THSFEM中的應力場由三個分量組成,應力函數導出的前兩部分合并起來表示為Pprβpr,界面裂紋尖端周圍的奇異項表示為Picrβicr,單元內各相材料應力表達式為

(10)

(11)

單元各邊界位移由相應邊界段節點位移插值得到,表達式為

(12)

(13)

(14)

式(14)中:

即β=H-1Gq,

(15)

其中單元剛度矩陣表達式為

Ke=GTH-1G

(16)

3 應力函數的構造

X-THSFEM單元模型中應力場函數的構造應考慮遠場應力和界面形狀對應力場的影響以及裂紋尖端附近應力奇異性。圖1中X-THSFEM單元模型中應力函數由3個部分組成,表達式為

Φ=Φpoly+Φrec+Φicr

(17)

式(17)中:Φpoly為Airy高階多項應力函數刻畫遠場應力;Φrec為互作用應力函數描述界面的形狀;Φicr為構造奇異應力場函數,描述界面裂紋尖端附近的應力集中。

3.1 多項式Airy應力函數的構造

為刻畫遠場應力,利用Airy應力函數構造高階多項應力函數。Airy高階多項應力函數在局部坐標中的表達式為

(18)

式(18)中:(xc,yc)為單元中心坐標,利用縮放系數Li=max|(xi,yi)-(xc,yc)|∕10,將全局坐標系(x,y)轉變為局部坐標系(ξ,η),其中:ξ=(xi-xc)∕Li,η=(yi-yc)∕Li。由于局部坐標系與全局坐標系互相平行,故可直接得到全局坐標系下的應力分量為

(19)

3.2 互作用應力函數的構造

f(x,y)=1所在界面呈直線形狀,當(x,y)→∞時,1/f(x,y)→0可以作為映射函數來構造基于界面形狀的互作用應力函數,直界面中由點到線的距離公式為

(20)

基于上述映射函數,可以構造局部坐標系(ξ,η)中基于界面形狀的互作用應力函數。各階段的互作用應力函數表示為

(21)

兩相材料第二部分應力分量為

(22)

3.3 界面裂紋尖端附近奇異性應力場函數的構造

以包含水平中心界面裂紋的兩相材料的X-

THSFEM模型單元為例,界面裂紋長度為2a,不同各向同性材料分別位于界面上下兩側。

如圖4所示,兩種材料的楊氏模量、泊松比和剪切模量分別為Ei、υi、ui(i=1,2),界面兩側材料中各存在一點,在以界面裂紋尖端為坐標原點的局部坐標系(X1O1Y1)和(X2O2Y2)下,L1和L2分別為點Q1和Q2與對應裂紋尖端的徑向距離。

將式(5)、式(6)代入式(1)中可得到點Q1和Q1在界面裂紋尖端1和2附近應力分量表達式為

圖4 包含水平中心界面裂紋的兩相材料坐標系Fig.4 Coordinate system of two-phase material with horizontal central interfacial crack

(23)

式(23)中:f1、f2、g1、g2、h1和h2為雙材料系數,插值函數和坐標L、θ已知,復合應力強度因子的實部K1和虛部K2相當于第2節中提到的應力系數βicr。同時考慮兩個界面裂紋尖端的影響,計算X-THSFEM中任意積分點的裂紋尖端奇異應力。由坐標變換矩陣T可以得到全局奇異應力分量為

(24)

全局應力場應充分反映遠離裂紋尖端和界面的應力、界面形狀對應力場的影響以及界面裂紋尖端附近的應力奇異性。因此,兩相材料中的全場應力分量可以由上述3個部分疊加而成,表達式為

(25)

4 最小二乘法確定尖端應力強度因子

應力強度因子可以通過應力和位移等參數進行求解,在數值方法求解中,可以通過路徑獨立積分比如相互獨立積分和J積分。將采用基于點的應力的最小二乘法計算應力強度因子。

如圖5所示,在裂紋長度為2a的兩相材料模型單元中,在裂紋尖端附近選取大量點形成半徑為r的環形曲面,將點的坐標及全局應力代入式(23)中得到線性方程

σ=A+K

(26)

式(26)中:σ包含式(10)、式(11)兩個部分;A包含與每個相點坐標有關的插值函數;界面裂紋應力強度因子K可以通過式(26)兩邊乘以A-1得到。

圖5 用于計算應力強度因子的域Fig.5 Fields used to calculate SIF

5 數值算例

通過X-THSFEM單元模型數值模擬結果與商用計算軟件ABAQUS的有限元結果或已有文獻結果進行比較加以驗證X-THSFEM的有效性和準確性。

算例1以包含水平邊界裂紋的雙材料矩形板為例。

如圖6所示,包含水平邊界裂紋的雙材料矩形板長3L=15 mm,寬度L=5 mm,邊界裂紋長度a=1 mm,上下層材料的楊氏模量和泊松比分別為E1=4.8×104MPa,E2=480 MPa,υ1=υ2=0.3,板材上下邊界施加1 MPa的拉伸載荷。相同幾何形態下,通過對X-THSFEM單元模型的分析和ABAQUS軟件所得到結果進行比較,如圖7所示為應力分量σx、σy、σxy云圖。

為比較包含邊界裂紋的兩相材料在E1/E2不同的情況下,無量綱應力強度因子的取值變化情況,如表1所示。從表1可知,X-THSFEM得到的裂紋尖端附近的無量綱SIF將與Gu等[17]給出的邊界元方法計算的參考結果和Jiang等[18]基于廣義有限差分法計算的參考結果非常接近。進一步說明本文所提方法的準確性。

算例2以包含水平中心界面裂紋的雙材料矩形板為例。

如圖8所示,板長2L=10 mm,寬度L=5 mm,中心界面裂紋長度a=2 mm,上下層材料的楊氏模量和泊松比分別為E1=4.8×104MPa,E2=480 MPa,υ1=υ2=0.3板材上下邊界施加1 MPa的拉伸載荷P。在相同幾何形態下,通過如圖9所示X-THSFEM模型與ABAQUS軟件所得到的水平應力分量σx云圖和垂直應力分量σy云圖,可以觀察到,X-THSFEM模型得到的這些應力云圖與ABAQUS軟件得到的應力云圖中應力分布趨于一致。

當E1/E2=100時,為保證應力強度因子的收斂性,如表3所示通過對應力系數的個數(多項式應力函數的項數M和互作用應力函數的項數N)不同取值,應力強度因子K1的變化。

圖8 包含水平中心界面裂紋的雙材料矩形板Fig.8 A double-material rectangular plate with a horizontal central interface crack

圖9 包含水平中心裂紋的雙材料矩形板X-THSFEM 模型和ABAQUS模型應力云圖Fig.9 X-THSFEM model and ABAQUS model of bi-material rectangular plate with horizontal central crack stress nephogram

表2 包含邊界裂紋的兩相材料無量綱應力強度因子Table 2 Dimensionless stress intensity factors for two-phase materials containing boundary cracks

如表3所示,當多項式應力函數的項數M=117時,互作用應力函數的項數N取7～25,應力強度因子K1保持不變,為保證應力強度因子的收斂性,最終可確定M=117,N=18。

當裂紋面與水平方向呈一定傾斜角度時,如圖10所示,其中L=5 mm,裂紋面傾斜角度為θ。上層材料的楊氏模量和泊松比分別為E1=4.8×104MPa,υ1=0.3;下層材料楊氏模量和泊松比分別為:E2=480 MPa,υ2=0.3,板材上下邊界施加1 MPa的拉伸載荷。

為證明本文方法的準確性,如圖11所示,當裂紋傾斜角度θ分別為15、30、45、60時,通過對X-THSFEM模型與ABAQUS軟件得到的水平應力分量包含不同傾斜角度界面裂紋的兩相材料在E1/E2不同的情況下,無量綱應力強度因子的取值變化情況,如表4所示,X-THSFEM得到的裂紋尖端附近的無量綱SIF與ABAQUS軟件得到的無量綱SIF趨于一致,進一步說明本文所提方法的準確性。

表3 應力系數不同的情況下應力強度因子取值Table 3 Values of SIF under different stress coefficients

圖10 包含傾斜界面裂紋的雙材料矩形板Fig.10 Bimaterial rectangular plate with inclined interface crack

σx云圖和垂直應力分量σy云圖進行比較,從圖11中可以觀察到,X-THSFEM模型得到的應力云圖與ABAQUS軟件得到的應力云圖中應力分布趨于一致。

算例3包含多個界面裂紋的層合板。

如圖12所示,對于由多種各向同性材料組成,包含不同類型界面裂紋的層合板,通過構建包含多個單元的擴展二相雜交應力有限元模型對界面裂紋進行斷裂力學分析,其中每個二相雜交應力單元包含邊界裂紋、中心界面裂紋以及同時包含中心界面裂紋和邊界裂紋。構建的擴展二相雜交應力有限元模型包含9個單元。

如圖13所示,每層材料彈性模量及泊松比:E1=E4=730 GPa,E2=E3=73 GPa,E5=E6=7.3 GPa,υ1=υ2=υ3=υ4=υ5=υ6=0.2。第1號單元兩側同時包含邊界裂紋,第2、3、5、7、8號單元只一側包含邊界裂紋,第4、9號單元包含水平中心裂紋,第6號單元包含水平中心裂紋以及兩側同時包含邊界裂紋。通過對包含多個單元的擴展二相雜交應力有限元模型進行分析得到應力分量σx、σy、σxy云圖,如圖14所示。為進一步說明X-THSFEM對于分析多種材料組成,包含不同類型界面裂紋的層合板具有很大優勢,如表5所示,研究裂紋尖端SIFs隨彈性模量比值的變化趨勢,選取兩組上下兩層材料彈性模量不同比值的有限元模型,當隨著單元內上下層材料彈性模量之比的增大,K1相應增大,K2的絕對值也相應增大。

表4 包含不同傾斜角度界面裂紋的兩相材料 無量綱應力強度因子Table 4 Dimensionless stress intensity factors for two-phase materials with interfacial cracks at different inclined angles

圖11 包含傾斜界面裂紋的雙材料矩形板X-THSFEM模型和ABAQUS模型的應力分量云圖Fig.11 Stress component nebulae of double-material rectangular plates with inclined interface crack of X-THSFEM model and ABAQUS model

圖12 包含不同種界面裂紋的多相材料矩形板Fig.12 Rectangular plate of polyphase material containing different kinds of interfacial cracks

圖13 擴展二相雜交應力有限元模型Fig.13 Extended two-phase hybrid stress finite element model

表5 彈性模量不同比值時界面裂紋應力強度因子Table 5 Stress intensity factors of interfacial crack with different ratio of elastic modulus

6 結論

提出了一種擴展二相雜交應力有限元法,通過構建包含邊界裂紋、水平中心裂紋以及傾斜裂紋的兩相材料的擴展二相雜交應力有限元單元模型對界面裂紋進行斷裂力學分析。應用Lagrange乘子法將各個界面裂紋面面力為零邊界條件和未斷裂界面的面力連續條件引入修正余能泛函中,在每一相材料域內假設高階多項式應力函數,充分考慮界面形狀的影響,構建互作用應力函數,引入雙材料界面裂紋裂尖附近奇異應力場函數。基于Fortran程序求解后得到單元應力場,采用最小二乘法計算了界面裂紋的應力強度因子。通過該方法得到的數值模擬結果和商業有限元軟件ABAQUS模擬結果進行比較,結果趨于一致,進而對由不同種各向同性材料組成,包含多種界面裂紋的層合板進行分析,說明該方法的有效性和準確性。為后期對界面裂紋的擴展研究提供理論和方案基礎。