豐艷霞，徐 旺，楊震霆，周震寰

（1. 大連理工大學(xué) 工程力學(xué)系，遼寧 大連 116024；2. 武漢科技大學(xué) 冶金工業(yè)過程系統(tǒng)科學(xué)湖北省重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430065）

0 引言

裂紋是工程結(jié)構(gòu)中的主要缺陷形式之一，其失穩(wěn)擴(kuò)展將導(dǎo)致安全事故[1].在工程實(shí)際中，多裂紋情況更為常見，應(yīng)力強(qiáng)度因子是定量評估裂紋是否穩(wěn)定的重要依據(jù)，因此多裂紋應(yīng)力強(qiáng)度因子（SIFs）的計(jì)算是斷裂力學(xué)中的一個(gè)重要問題.

國內(nèi)外學(xué)者已經(jīng)在多裂紋問題的應(yīng)力強(qiáng)度因子求解方面開展了大量研究.奇異積分方程法和復(fù)變函數(shù)法是常用的兩類解析方法.基于奇異積分方程法，CHEN Y Z等[2]求解了有限尺寸板中多裂紋的SIFs.MONFARED M M等[3]分析了非均勻正交各向異性平面中多裂紋問題，并計(jì)算了對應(yīng)的SIFs. HAGHIRI H等[4]求解了功能梯度正交各向異性半平面邊界處多平行或垂直裂紋的SIFs.ZHAO J F等[5]計(jì)算了無限大板中多孔邊裂紋的SIFs.解析方法雖然可以獲得解析表達(dá)式，但是還受到無限域問題、簡單邊界條件等限制.為克服上述問題，大量數(shù)值方法被應(yīng)用于多裂紋問題的求解，如有限元法（FEM）[6]、邊界元法（BEM）[7]、無網(wǎng)格法[8].FEM雖然更簡單，但是針對斷裂問題的計(jì)算精度不高.BEM和無網(wǎng)格方法更精確，但是依賴于基本解，無網(wǎng)格方法的計(jì)算效率尚待解決.

基于解析辛方法[9-11]和FEM，提出一種適用于計(jì)算平面多裂紋問題的辛離散有限元方法（FEDSM）.對含裂紋結(jié)構(gòu)進(jìn)行網(wǎng)格劃分，并將其劃分為裂紋尖端的奇應(yīng)力區(qū)域和剩余的無奇異性區(qū)域.在近場內(nèi)建立哈密頓系統(tǒng)，利用獲得的辛本征函數(shù)進(jìn)行節(jié)點(diǎn)未知量變換，得到FEDSM的計(jì)算列式，獲得多裂紋對應(yīng)的SIFs及其裂尖附近奇異物理場的解析表達(dá)式.

1 近場內(nèi)哈密頓體系與解析解

1.1 問題的提出

整體含裂紋結(jié)構(gòu)的區(qū)域劃分見圖1（a）.為了方便表示，近場選取半徑為RN的圓形.在近場內(nèi)，選取極坐標(biāo)(r，θ)，r軸沿徑向方向，原點(diǎn)位于裂紋尖端，對應(yīng)的笛卡爾直角坐標(biāo)為(x，y)，ГFN為兩類區(qū)域的分界線，見圖1（b）.

圖1 多裂紋結(jié)構(gòu)示意Fig.1 representation of multi-crack structure

1.2 近場內(nèi)哈密頓控制方程

為在近場內(nèi)建立哈密頓體系，定義廣義坐標(biāo)變量ξ，并定義變量對其導(dǎo)數(shù)為=?f/?ξ（f為任意變量）.對于平面應(yīng)力問題，哈密頓體系下的基本未知量[9]可以表示為

式中，q和p分別為原變量及其對偶變量向量，Pa·m；u和ν分別為徑向和環(huán)向位移，m；sr和srθ為廣義應(yīng)力，N，其中sr=rσr和srθ=rτrθ；σr和τrθ分別為徑向正應(yīng)力和剪應(yīng)力，N；r為極坐標(biāo)半徑，m；θ為極坐標(biāo)旋轉(zhuǎn)角度，°.忽略裂紋面載荷的影響，極坐標(biāo)下平面應(yīng)力問題的哈密頓控制方程可以表示為

式中，H為哈密頓算子矩陣；Ψ為哈密頓體系下的基本未知量.

式中，E為彈性模量，Pa；υ為待定系數(shù).

對應(yīng)的裂紋表面邊界條件為

1.3 近場內(nèi)解析解

在哈密頓體系下，哈密頓控制方程式（2）可以歸結(jié)為哈密頓算子矩陣H的本征問題，即

式中，μ和Ψ為辛特征對.

式（2）的本征解可以分為兩類[9].

① 零本征值本征解（ 0μ= ）為

② 非零本征值本征解（μ=n/2，n= 1，2，3，…）為

式中，上標(biāo)s和a分別為對稱和反對稱本征解；A1～A4為待定系數(shù).α11=α12=α13=α14=-α21=-α23=1，α12= -α24=(-3+υ-μ-υμ)/(3-υ-μ-υμ)，α32=α34=Eμ(3-μ)/(3-υ-μ-υμ)，α42=-α44=Eμ(1-μ)/，α31=α33=-α41=-α43=Eμ(1+υ).

由式（6）～式（8）可以表示式（2）的通解，其中位移函數(shù)可以表示為

式中，uc為位移函數(shù)向量，m；q1為徑向位移，m，q1=u；q2為環(huán)向位移，m，q2=v；ri為有限元節(jié)點(diǎn)坐標(biāo)半徑，m；θi為有限元節(jié)點(diǎn)坐標(biāo)旋轉(zhuǎn)角度，°；i為第i個(gè)節(jié)點(diǎn)，i= 1 -ZN（ZN為近場內(nèi)節(jié)點(diǎn)數(shù)）；j為第j階辛本征解，j= 1 -ZM（ZM為選取的本征解項(xiàng)數(shù)）；c為待定系數(shù)向量；Φ為變換矩陣.

2 辛離散有限元法

2.1 多裂紋結(jié)構(gòu)的有限元列式

對于整體結(jié)構(gòu)，其有限元列式可以表示為

式中，u為節(jié)點(diǎn)位移向量，m；f為外力向量，N；K為有限元?jiǎng)偠染仃?；下?biāo)F和N分別代表遠(yuǎn)場和近場區(qū)域.

為了分析多裂紋問題，將整個(gè)裂紋體劃分為幾個(gè)子域，每個(gè)子域僅包含一個(gè)裂紋，見圖2.

圖2 子域劃分示意Fig.2 signal of sub-domains divisions

圖2 中Mi為第i個(gè)子域，Гi為相鄰域分界線，i{1，2，∈ …，N}.各子域內(nèi)的有限元列式如下.

第1個(gè)子域?yàn)?/p>

第i個(gè)子域?yàn)?/p>

第N個(gè)子域?yàn)?/p>

由式（11）～式（13）可知，整體結(jié)構(gòu)的有限元列式可以表示為

式中，

2.2 多裂紋結(jié)構(gòu)FEDSM列式及應(yīng)力強(qiáng)度因子

通過求解辛離散有限元列式（14），可獲得問題在近場內(nèi)的解析解為

在近場內(nèi)，廣義應(yīng)力的奇異項(xiàng)為

根據(jù)斷裂力學(xué)定義，I型和II型SIFs為

式中，KI為型應(yīng)力強(qiáng)度因子；KII為型應(yīng)力強(qiáng)度因子.

3 數(shù)值算例

數(shù)值算例將對工程中常見的I型和I、II混合型平面斷裂問題展開研究，驗(yàn)證所提方法的正確性和精確性.

3.1 算例1

為驗(yàn)證所提出多裂紋問題FEDSM方法的準(zhǔn)確性，選取含有一對變裂紋長度的彈性圓盤，見圖3，圓盤半徑為R，裂紋長度為a，在該圓盤外邊界作用有沿徑向向外的均勻應(yīng)力σ.在計(jì)算中，選取ANSYS中的PLANE183單元對圓盤建模，近場半徑RN= 0.1R，辛本征解取20 項(xiàng).由于結(jié)構(gòu)具有對稱性，因此只需考慮其中一條裂紋的SIFs即可.圖4為SIFs隨裂紋長度的變化，F(xiàn)EDSM的計(jì)算結(jié)果與文獻(xiàn)[12]的結(jié)果吻合良好，直接證明了所提方法是準(zhǔn)確的.此外，為了說明近場尺寸的影響，取a/R為 0.5，圖5為不同RN/R對應(yīng)的SIFs隨本征解展開項(xiàng)數(shù)的變化，當(dāng)RN/R≥10%時(shí)，近場尺寸對計(jì)算結(jié)果已無影響.

圖3 含邊裂紋圓盤Fig.3 edge-cracked circular disk

圖4 SIFs隨裂紋長度的變化Fig.4 variations of SIFs with crack length

圖5 近場尺寸對應(yīng)力強(qiáng)度因子的影響 Fig.5 effect of the Near field size on the stress strength factor

3.2 算例2

圖6 為含有一對邊裂紋的矩形板，矩形板尺寸為2h×2b，裂紋長度分別為a1和a2，上下兩端受到均勻拉伸載荷作用σ.有限元單元類型、近場半徑和本征解項(xiàng)數(shù)與算例1相同.算例中的SIFs數(shù)值均為無量綱數(shù)值，其中無量綱特征參數(shù) 0Kbσ= .

圖6 含雙裂紋矩形板Fig.6 rectangular plate with double edge-cracked

圖7 為左端裂紋SIF隨裂紋長度a1的變化.由圖7可知，本計(jì)算結(jié)果與文獻(xiàn)[13]結(jié)果一致，再一次證明所提出的FEDSM適用于多裂紋問題的SIFs求解.此外，圖7表明非對稱邊裂紋的SIFs隨裂紋長度的增加而逐漸增加.該現(xiàn)象表明整體結(jié)構(gòu)和裂紋幾何參數(shù)是影響多裂紋問題的重要影響因素.

圖7 非對稱邊裂紋SIFs隨長度a2的變化Fig.7 variations of SIFs for non-symmetrical edge cracks with the crack length a2

3.3 算例3

算例3為較復(fù)雜的I、II混合型平面斷裂問題.考慮含有1對共線斜裂紋的矩形板，其幾何參數(shù)見圖8，長寬比h/b為2，裂紋長度均為a.有限元單元類型、近場半徑、本征解項(xiàng)數(shù)和無量綱特征參數(shù)與算例2相同.

圖8 含雙斜裂紋的矩形板Fig.8 rectangular plate with two inclined cracks

表1和表2分別列出了不同裂紋角度和長度對應(yīng)的I型和II型SIFs.為驗(yàn)證本計(jì)算結(jié)果的正確性，表格中同時(shí)給出了ANSYS的計(jì)算結(jié)果，選用單元類型與FEDSM相同，并且裂紋尖端單元采用1/4節(jié)點(diǎn)單元.從表1和表2中的數(shù)據(jù)可以看出，F(xiàn)EDSM和ANSYS獲得的I型和II型SIFs數(shù)據(jù)結(jié)果一致.I型和II型SIFs都隨裂紋長度的增加而增加；隨著裂紋角β的增加，I型SIFs減小，而II型SIFs繼續(xù)增大.該現(xiàn)象與算例2類似，表明裂紋幾何參數(shù)是影響其SIFs數(shù)值的重要因素.

表1 不同裂紋角度和長度雙斜裂紋的無量綱I型SIFs（h/b = 2）Tab.1 normalized mode I SIFs of two inclined cracks for various crack angles and lengths (h/b = 2)

表2 不同裂紋角度和長度雙斜裂紋的無量綱II型SIFs（h/b = 2）Tab.2 normalized mode II SIFs of two inclined cracks for various crack angles and lengths (h/b = 2)

除SIFs外，裂紋尖端附近的奇異應(yīng)力場分布也是結(jié)構(gòu)斷裂分析的關(guān)鍵.選取在裂紋尖端周圍的長和寬為l1=l2= 0.1b的正方形路徑，見圖9（a），對應(yīng)的奇異應(yīng)力變化規(guī)律見圖9（b）～圖9（d），所獲應(yīng)力分量與FEM結(jié)果非常吻合.

圖9 沿正方形路徑l1 = l2 =0.1b的無量綱應(yīng)力分量Fig.9 dimensionless stress components along the square path with l1=l2=0.1b

4 結(jié)論

提出一種適用于多裂紋結(jié)構(gòu)應(yīng)力強(qiáng)度因子計(jì)算的辛離散有限元方法.該方法通過將整體結(jié)構(gòu)劃分為多個(gè)子域?qū)崿F(xiàn)多裂紋問題簡化，并在各個(gè)子域內(nèi)將解析辛方法與傳統(tǒng)有限元方法相結(jié)合，直接獲得應(yīng)力強(qiáng)度因子和該區(qū)域內(nèi)物理場的解析表達(dá)式.與其他數(shù)值方法相比，辛離散有限元方法具有3個(gè)優(yōu)勢：① 通過在近場內(nèi)進(jìn)行未知變量變換，將未知量的數(shù)量大幅降低，直接減少斷裂分析的計(jì)算時(shí)間和計(jì)算機(jī)的存儲需求.② 該方法無需引入新的特殊單元，并且無需后處理程序即可求解SIFs.③ 該方法可以直接獲得近場中的位移和應(yīng)力的解析解.