有限矩形板基于權(quán)函數(shù)法的應(yīng)力強(qiáng)度因子計算

陳 昊，董登科

(中國航空工業(yè)集團(tuán)有限公司 中國飛機(jī)強(qiáng)度研究所，陜西 西安 710065)

1 引 言

飛機(jī)在使用過程中，其材料和結(jié)構(gòu)會因為各種原因而產(chǎn)生裂紋缺陷。裂紋在載荷和環(huán)境條件的共同作用下不斷擴(kuò)展，最終會導(dǎo)致結(jié)構(gòu)斷裂，引發(fā)嚴(yán)重的后果。因此，對含裂紋材料和結(jié)構(gòu)進(jìn)行研究對于航空工業(yè)具有重要的意義。

基于斷裂力學(xué)的基本理論對線彈性裂紋尖端場的分析表明，所有影響裂紋尖端力學(xué)狀態(tài)的變量，僅通過應(yīng)力強(qiáng)度因子K來影響裂紋體的行為。一旦K確定，裂紋尖端區(qū)內(nèi)的應(yīng)力應(yīng)變和位移分量也隨之確定。鑒于K在裂紋問題分析中的核心地位，各種裂紋體在任意載荷下K的求解成為斷裂力學(xué)中非常重要的研究內(nèi)容。

許多學(xué)者對應(yīng)力強(qiáng)度因子的求解方法進(jìn)行了深入的研究。Wu等建立了基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法[1]和Glinka-Shen基于兩種(或三種)參考載荷和一個幾何條件的通用權(quán)函數(shù)法[4]，并結(jié)合有限元軟件ANSYS計算應(yīng)力強(qiáng)度因子的權(quán)函數(shù)。

2 權(quán)函數(shù)法

權(quán)函數(shù)法是一種通用性很強(qiáng)、效率很高的應(yīng)力強(qiáng)度因子求解方法。因為權(quán)函數(shù)把載荷條件和幾何條件分開考慮，只包含裂紋體幾何特征，不受載荷影響。因此，一旦確定了權(quán)函數(shù)m(a,x)，就可以計算不同裂紋尺寸在不同載荷下的應(yīng)力強(qiáng)度因子，并且一次計算就能夠得到一種載荷情況下一條完整的K(a)曲線，而有限元等數(shù)值方法一次計算只能得到一個特定裂紋長度a下的K值。

含裂紋結(jié)構(gòu)的應(yīng)力強(qiáng)度因子K的計算是對權(quán)函數(shù)m(a,x)和應(yīng)力分布σ(x)的乘積進(jìn)行積分。此處的應(yīng)力分布是未含裂紋結(jié)構(gòu)在假想裂紋處的應(yīng)力分布。

2.1 裂紋的權(quán)函數(shù)法

對于二維裂紋結(jié)構(gòu)的應(yīng)力分布形式，應(yīng)力強(qiáng)度因子計算公式為[3]：

(1)

式中，α和ξ分別是無量綱的裂紋長度和坐標(biāo)(a和x分別是實際裂紋長度和坐標(biāo)，W是裂紋體的特征尺寸)，σ(ξ)是不考慮裂紋的存在、假想裂紋位置的應(yīng)力。

斷裂力學(xué)中，權(quán)函數(shù)求解法最早由Bueckner提出，Rice對其做了進(jìn)一步的研究。高精度權(quán)函數(shù)的推導(dǎo)與驗證是采用權(quán)函數(shù)法處理裂紋問題的關(guān)鍵。本文二維裂紋的權(quán)函數(shù)推導(dǎo)方法主要為：基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法和基于多個參考載荷情況的權(quán)函數(shù)法。

2.1.1 基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法

Wu等建立了基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法，對裂紋在一種參考載荷情況下的裂紋面位移分別提出了規(guī)范化的級數(shù)展開表達(dá)式[3]。利用裂紋尖端場的特性、自洽條件、裂紋嘴位移和合理的幾何條件及其各種組合，推導(dǎo)位移級數(shù)的相關(guān)函數(shù)，進(jìn)而由裂紋面位移對裂紋長度求偏導(dǎo)數(shù)，確定裂紋權(quán)函數(shù)m(α,ξ/α)中的系數(shù)βi(α)。

對于中心裂紋，該解析權(quán)函數(shù)的通式為：

(2)

(3)

對于邊緣裂紋，該解析權(quán)函數(shù)的通式為：

(4)

(5)

式中，F(xiàn)i(α)為裂紋面位移展開級數(shù)的系數(shù)。

2.1.2 基于多個參考載荷情況的權(quán)函數(shù)法

考慮到早期基于裂紋張開位移的P-A方法在求解精度和適用范圍方面的缺陷以及在權(quán)函數(shù)的獲取過程中所需數(shù)值積分的麻煩，一些研究者改用其他方法，即利用多個參考載荷情況的已知K解來推導(dǎo)權(quán)函數(shù)。

Glinka和Shen確立了一種適用于所有裂紋幾何的通用權(quán)函數(shù)法[5]，該權(quán)函數(shù)的確定是基于兩個或三個參考載荷情況及相關(guān)幾何條件，其通用權(quán)函數(shù)基本公式為：

(6)

式中，系數(shù)M1、M2、M3需要基于兩種(或三種)載荷情況下的已知K解(附加幾何條件)，通過求解一組聯(lián)立方程確定。

上述兩種權(quán)函數(shù)法都需要一個或多個參考載荷情況下的K值(及裂紋嘴位移)[7]。獲得這些已知條件的途徑，一是從文獻(xiàn)中查找理論精確解或其他方法得到的高精度數(shù)值解，二是采用經(jīng)過反復(fù)驗證的數(shù)值方法(有限元、邊界元等)計算確定。本文采用有限元計算參考載荷下的應(yīng)力強(qiáng)度因子。

2.2 格林函數(shù)

由參考載荷情況應(yīng)力強(qiáng)度因子導(dǎo)出權(quán)函數(shù)之后，結(jié)合給定應(yīng)力分布下的應(yīng)力強(qiáng)度因子，并與其他方法得出的應(yīng)力強(qiáng)度因子進(jìn)行對比，是目前常用的評估應(yīng)力強(qiáng)度因子算法準(zhǔn)確性的一種方式。但是，這種方法實際上并不能很準(zhǔn)確地反映權(quán)函數(shù)的精度，因為由權(quán)函數(shù)法計算應(yīng)力強(qiáng)度因子，在導(dǎo)出權(quán)函數(shù)后，需要對權(quán)函數(shù)m(a，x)和無裂紋體假想裂紋處的應(yīng)力σ(x)的乘積沿整個裂紋面做積分，而積分具有平均效應(yīng)，所以由積分得到的應(yīng)力強(qiáng)度因子不能真實反映權(quán)函數(shù)這一參數(shù)本身的精確性。當(dāng)應(yīng)力分布沿裂紋面劇烈變化，尤其當(dāng)應(yīng)力方向改變時，應(yīng)力強(qiáng)度因子的誤差甚至可能會遠(yuǎn)大于權(quán)函數(shù)本身的誤差。因此，評價權(quán)函數(shù)算法精度的更好的方法是針對不同裂紋長度，沿著裂紋面逐點對權(quán)函數(shù)進(jìn)行比較，即比較格林函數(shù)[8]。這種比較方式不會受權(quán)函數(shù)與應(yīng)力乘積的積分平均效應(yīng)的影響。

格林函數(shù)G(a，x)又稱為影響函數(shù)[9]，表示裂紋面在任意位置受到一對單位集中力P作用時的無量綱應(yīng)力強(qiáng)度因子。格林函數(shù)G(a，x)與裂紋面受一對單位集中力P作用下的應(yīng)力強(qiáng)度因子K沿著集中力作用的位置逐點對應(yīng)，其與權(quán)函數(shù)的關(guān)系如下：

(7)

3 有限矩形板裂紋應(yīng)力強(qiáng)度因子計算

對于平板結(jié)構(gòu)，工程上常見的裂紋形式為中心裂紋和邊緣裂紋。很多學(xué)者提出了各種求解權(quán)函數(shù)的方法，本文分別采用基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法和基于多個參考載荷條件的Glinka-Shen通用權(quán)函數(shù)法計算有限矩形板中心裂紋和邊緣裂紋的應(yīng)力強(qiáng)度因子。

3.1 基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法

3.1.1 有限矩形板中心裂紋問題

考慮有限矩形板中心裂紋問題，板寬度為2W、長度為2H，裂紋長度為2a(有限矩形板中心裂紋如圖1所示)。選取板半寬W作為該裂紋幾何特征尺寸，無量綱裂紋長度α和坐標(biāo)ξ分別為：α=a/W,ξ=x/W。

圖1 有限矩形板中心裂紋

選取裂紋面均布應(yīng)力作為推導(dǎo)中心裂紋權(quán)函數(shù)的參考載荷。對于裂紋面受均布應(yīng)力的中心裂紋，α=0時，應(yīng)力強(qiáng)度因子和裂紋嘴位移都存在理論精確解；α=1時，應(yīng)力強(qiáng)度因子存在理論極限值。

α=0時：

fr(0)=1Vr(0)=2.0

α=1時：

(8)

對于有限矩形板中心裂紋，本文利用有限元軟件ANSYS19.0建立模型，模擬不同裂紋長度的中心裂紋，并進(jìn)行應(yīng)力分析。模擬裂紋長度α包括0.01、0.05、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9共11個模型。

為了簡化計算，提高運(yùn)行效率，本文建立1/4模型，并施加對稱邊界條件。裂紋尖端處采用1/4節(jié)點奇異單元，網(wǎng)格類型設(shè)置單元PLANE183，裂紋尖端網(wǎng)格如圖2所示。材料屬性：楊氏模量E=7.17×105MPa，泊松比v=0.3。

圖2 裂紋尖端網(wǎng)格

利用ANSYS進(jìn)行應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur計算，結(jié)果如表1所示。

表1 有限元計算結(jié)果

根據(jù)有限元計算得到的Kr和Ur值，采用式(9)對其進(jìn)行正則化處理，得到無量綱應(yīng)力強(qiáng)度因子fr(α)和無量綱裂紋嘴張開位移Vr(α)，然后再通過多項式擬合，得到fr(α)和Vr(α)的高精度擬合多項式。

Vr(α)=ur(α,ξ=0)E′/(σ0α)

(9)

(10)

式中的多項式系數(shù)λi和γn如表2所示。

表2 均布應(yīng)力下fr(α)和Vr(α)的多項式系數(shù)

把fr(α)和Vr(α)的高精度擬合多項式代入式(11)，可以方便地求得?(α)的積分式。將其代入式(3)，就能確定裂紋面位移函數(shù)中的Fj(α)，j=1～3。

(11)

再將Fj(α)代入式(2)，便能確定βi(α)的值。對βi(α)的計算結(jié)果進(jìn)行多項式擬合，得到有限矩形板中心裂紋βi(α)的多項式擬合表達(dá)式：

β1=2

β2=(-26.32α7+90.14α6-126.2α5+91.83α4-38.16α3+9.297α2-0.5675α-0.0001263)/(1-α)＾1.5)

β3=(17.49α7-90.07α6+164.4α5-144.4α4+65.71α3-14.97α2+1.378α-0.0002538)/(1-α)＾1.5)

β4=(19.82α7-33.51α6+1.365α5+29.83α4-23.12α3+6.86α2-0.7436α+0.000396)/(1-α)＾1.5)

(12)

把所求的βi(α)代入下式，就能確定有限矩形板中心裂紋的解析權(quán)函數(shù)。

3.1.2 有限矩形板邊緣裂紋問題

考慮有限矩形板邊緣裂紋問題，板寬度為W、長度為H，裂紋長度為a(有限矩形板邊緣裂紋如圖3所示)。選取板寬W作為該裂紋幾何特征尺寸，無量綱裂紋長度α和坐標(biāo)ξ分別為：α=a/W，ξ=x/W。

圖3 有限矩形板邊緣裂紋

選取裂紋面均布應(yīng)力作為推導(dǎo)邊緣裂紋權(quán)函數(shù)的參考載荷。對于裂紋面受均布應(yīng)力的邊緣裂紋，α=0時，應(yīng)力強(qiáng)度因子和裂紋嘴位移都存在理論精確解；α=1時，應(yīng)力強(qiáng)度因子存在理論極限值。

α=0時：

fr(0)=1.1215

Vr(0)=2.9086

α=1時：

(14)

對于有限矩形板邊緣裂紋，本文利用有限元軟件ANSYS19.0建立模型，模擬不同裂紋長度的邊緣裂紋，并進(jìn)行應(yīng)力分析。模擬裂紋長度α包括0.01、0.05、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9共11個模型。

為了簡化計算，提高運(yùn)行效率，本文建立1/2模型，并施加對稱邊界條件。裂紋尖端處采用1/4節(jié)點奇異單元，網(wǎng)格類型設(shè)置單元PLANE183。材料屬性：楊氏模量E=7.17×105MPa，泊松比v=0.3。

利用ANSYS進(jìn)行應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur計算，結(jié)果如表3所示。

表3 有限元計算結(jié)果

根據(jù)有限元計算得到的Kr和Ur值，根據(jù)式(15)對其進(jìn)行正則化處理，得到無量綱應(yīng)力強(qiáng)度因子fr(α)和無量綱裂紋嘴張開位移Vr(α)，然后再通過多項式擬合，得到fr(α)和Vr(α)的高精度擬合多項式。

(15)

式中的多項式系數(shù)λi和γn如表4所示。

表4 均布應(yīng)力下fr(α)和Vr(α)的多項式系數(shù)

把fr(α)和Vr(α)的高精度擬合多項式代入公式，可以方便地求得?(α)的積分式。將其代入式(4)，就能確定裂紋面位移函數(shù)中的Fj(α),j=1～3。再將Fj(α)代入式(5)，便能確定βi(α)的值。對βi(α)的計算結(jié)果進(jìn)行多項式擬合，得到有限矩形板邊緣裂紋βi(α)的多項式擬合表達(dá)式：

β1=2
β2=-68.01α7+208.4α6-274.4α5+205α4-105.1α3+32.42α2+1.568α+0.845
β3=246α7-668.3α6+697α5-316.4α4+44.01α3+23.32α2-9.025α+1.573
β4=-364.4α7+1065α6-1230α5+675.1α4-166.7α3+3.357α2+5.667α-0.786
β5=133α7-392.7α6+459.1α5-256.2α4+64.59α3-2.57α2-1.353α+0.05477

(16)

把所求的βi(α)代入下式，就能確定有限矩形板邊緣裂紋的解析權(quán)函數(shù)。

3.2 基于多個參考載荷條件的Glinka-Shen通用權(quán)函數(shù)法

本文求解權(quán)函數(shù)的另一種方法是基于多個參考載荷條件的Glinka-Shen通用權(quán)函數(shù)法。

3.2.1 有限矩形板中心裂紋問題

對于有限矩形板中心裂紋，選取均布應(yīng)力、線性應(yīng)力、二次應(yīng)力作為推導(dǎo)中心裂紋權(quán)函數(shù)的參考載荷。

均布應(yīng)力：σ(x)=1

線性應(yīng)力：σ(x)=(1-x/a)

二次應(yīng)力：σ(x)=(1-x/a)2

(18)

本文利用有限元軟件ANSYS19.0建立矩形板有限元模型，模擬不同應(yīng)力條件下不同裂紋長度的中心裂紋，并進(jìn)行應(yīng)力分析。模擬裂紋長度α包括0.01、0.05、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9共11個模型。

為了簡化計算，提高運(yùn)行效率，本文建立1/4模型，并施加對稱邊界條件。在裂紋尖端處采用1/4節(jié)點奇異單元，其余網(wǎng)格類型設(shè)置為單元PLANE183。材料屬性：楊氏模量E=7.17×105MPa，泊松比v=0.3。

利用ANSYS進(jìn)行應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur計算，結(jié)果如表5所示。

表5 應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur的有限元計算結(jié)果

將通過ANSYS求解得到的不同應(yīng)力分布下的應(yīng)力強(qiáng)度因子作為參考應(yīng)力強(qiáng)度因子和對應(yīng)的應(yīng)力分布函數(shù)，代入式(19)，并聯(lián)立方程組，可求得不同裂紋長度對應(yīng)的M1、M2、M3的值，如表6所示。

(19)

表6 不同裂紋長度對應(yīng)的M1、M2、M3的數(shù)值

續(xù)表6

根據(jù)計算得到的M1、M2、M3的值，對其進(jìn)行正則化處理，然后再通過多項式擬合，得到M1、M2、M3的高精度擬合多項式：

M1=0.7752α8+0.1064α7-0.7361α6-1.112α5-0.6818α4+1.542α3+0.06072α2+0.6471α-0.2981
M2=-4.281α8+1.314α7+8.549α6-1.274α5-4.031α4-0.2923α3+2.151α2-2.121α+1.52
M3=4.24α8-1.959α7-10.1α6+4.009α5+7.003α4-2.499α3-2.462α2+2.46α-0.4551

(20)

把所求得的M1、M2、M3的擬合多項式代入下式，就能確定有限矩形板中心裂紋的解析權(quán)函數(shù)，可求解對應(yīng)載荷作用下的應(yīng)力強(qiáng)度因子。

(21)

3.2.2 有限矩形板邊緣裂紋問題

對于有限矩形板邊緣裂紋，選取均布應(yīng)力、線性應(yīng)力、二次應(yīng)力作為推導(dǎo)邊緣裂紋權(quán)函數(shù)的參考載荷。本文利用有限元軟件ANSYS19.0建立有限矩形板的有限元模型，模擬不同應(yīng)力條件下不同裂紋長度的邊緣裂紋，并進(jìn)行應(yīng)力分析。模擬的裂紋長度α包括0.01、0.05、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9共11個模型。

為了簡化計算，提高運(yùn)行效率，本文建立1/2模型，并施加對稱邊界條件。在裂紋尖端處采用1/2節(jié)點奇異單元，其余網(wǎng)格類型設(shè)置為單元PLANE183。有限元模型的材料屬性：楊氏模量E=7.17×105MPa，泊松比v=0.3。

利用ANSYS進(jìn)行應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur計算，結(jié)果如表7所示。

表7 應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur的有限元計算結(jié)果

續(xù)表7

將通過ANSYS求解得到的不同應(yīng)力分布下的應(yīng)力強(qiáng)度因子作為參考應(yīng)力強(qiáng)度因子和對應(yīng)的應(yīng)力分布函數(shù)，代入式(19)，并聯(lián)立方程組，可求得不同裂紋長度對應(yīng)的M1、M2、M3的值，如表8所示。

表8 不同裂紋長度對應(yīng)的M1、M2、M3的數(shù)值

根據(jù)計算得到的M1、M2、M3的值，對其進(jìn)行正則化處理，然后再通過多項式擬合，得到M1、M2、M3的高精度擬合多項式。

M1=-0.1514α8+0.005979α7+1.767α6+1.647α5-1.11α4-0.04461α3+1.221α2+2.335α+1.619

M2=0.868α8+0.05433α7-7.228α6-5.851α5+4.987α4+0.8334α3-3.057α2-4.74α-2.548

M3=2.115α8+4.88α7+1.841α6-2.226α5+1.432α4+7.806α3+7.94α2+10.53α+6.98

(22)

把所求得的M1、M2、M3的擬合多項式代入公式，就能確定有限矩形板邊緣裂紋的解析權(quán)函數(shù)，可求解對應(yīng)載荷作用下的應(yīng)力強(qiáng)度因子。

4 權(quán)函數(shù)精度評價

格林函數(shù)G(a,ξ/α)又稱為影響函數(shù)，表示裂紋面在任意位置受到一對集中力P作用時的無量綱應(yīng)力強(qiáng)度因子。

4.1 中心裂紋的格林函數(shù)

在有限矩形板中心裂紋的裂紋面任意位置處施加一對集中力P(單位厚度)，計算格林函數(shù)并繪制函數(shù)的圖像，如圖4所示。

(23)

圖4 中心裂紋的格林函數(shù)

4.2 邊緣裂紋的格林函數(shù)

在有限矩形板邊緣裂紋的裂紋面任意位置處施加一對集中力P(單位厚度)，計算格林函數(shù)并繪制函數(shù)的圖像，如圖5所示。

(24)

圖5 邊緣裂紋的格林函數(shù)

5 結(jié) 論

本文通過不同的權(quán)函數(shù)法計算有限矩形板中心裂紋和邊緣裂紋的權(quán)函數(shù)，主要內(nèi)容如下：

(1)基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法，是通過參考載荷下的裂紋面位移Ur(a,x)對裂紋長度a求偏導(dǎo)數(shù)推導(dǎo)確定。

(2)使用ANSYS計算參考應(yīng)力強(qiáng)度因子，確定權(quán)函數(shù)的關(guān)鍵一步是求得參考載荷情況下的高精度裂紋張開位移Ur(a,x)，然后得到權(quán)函數(shù)的解析表達(dá)式m(α,ξ)。

(3)基于Glinka-Shen通用權(quán)函數(shù)法采用三個參考應(yīng)力計算相應(yīng)的參考應(yīng)力強(qiáng)度因子，計算權(quán)函數(shù)中的待定參數(shù)M1、M2、M3，最終確定應(yīng)力強(qiáng)度因子的權(quán)函數(shù)。

(4)采用格林函數(shù)驗證了所求權(quán)函數(shù)的正確性，為飛機(jī)結(jié)構(gòu)和材料的應(yīng)力強(qiáng)度因子確定、疲勞裂紋擴(kuò)展預(yù)測提供了一種簡便、高效的方法。