陳 昊,董登科
(中國(guó)航空工業(yè)集團(tuán)有限公司 中國(guó)飛機(jī)強(qiáng)度研究所,陜西 西安 710065)
飛機(jī)在使用過(guò)程中,其材料和結(jié)構(gòu)會(huì)因?yàn)楦鞣N原因而產(chǎn)生裂紋缺陷。裂紋在載荷和環(huán)境條件的共同作用下不斷擴(kuò)展,最終會(huì)導(dǎo)致結(jié)構(gòu)斷裂,引發(fā)嚴(yán)重的后果。因此,對(duì)含裂紋材料和結(jié)構(gòu)進(jìn)行研究對(duì)于航空工業(yè)具有重要的意義。
基于斷裂力學(xué)的基本理論對(duì)線彈性裂紋尖端場(chǎng)的分析表明,所有影響裂紋尖端力學(xué)狀態(tài)的變量,僅通過(guò)應(yīng)力強(qiáng)度因子K來(lái)影響裂紋體的行為。一旦K確定,裂紋尖端區(qū)內(nèi)的應(yīng)力應(yīng)變和位移分量也隨之確定。鑒于K在裂紋問(wèn)題分析中的核心地位,各種裂紋體在任意載荷下K的求解成為斷裂力學(xué)中非常重要的研究?jī)?nèi)容。
許多學(xué)者對(duì)應(yīng)力強(qiáng)度因子的求解方法進(jìn)行了深入的研究。Wu等建立了基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法[1]和Glinka-Shen基于兩種(或三種)參考載荷和一個(gè)幾何條件的通用權(quán)函數(shù)法[4],并結(jié)合有限元軟件ANSYS計(jì)算應(yīng)力強(qiáng)度因子的權(quán)函數(shù)。
權(quán)函數(shù)法是一種通用性很強(qiáng)、效率很高的應(yīng)力強(qiáng)度因子求解方法。因?yàn)闄?quán)函數(shù)把載荷條件和幾何條件分開(kāi)考慮,只包含裂紋體幾何特征,不受載荷影響。因此,一旦確定了權(quán)函數(shù)m(a,x),就可以計(jì)算不同裂紋尺寸在不同載荷下的應(yīng)力強(qiáng)度因子,并且一次計(jì)算就能夠得到一種載荷情況下一條完整的K(a)曲線,而有限元等數(shù)值方法一次計(jì)算只能得到一個(gè)特定裂紋長(zhǎng)度a下的K值。
含裂紋結(jié)構(gòu)的應(yīng)力強(qiáng)度因子K的計(jì)算是對(duì)權(quán)函數(shù)m(a,x)和應(yīng)力分布σ(x)的乘積進(jìn)行積分。此處的應(yīng)力分布是未含裂紋結(jié)構(gòu)在假想裂紋處的應(yīng)力分布。
對(duì)于二維裂紋結(jié)構(gòu)的應(yīng)力分布形式,應(yīng)力強(qiáng)度因子計(jì)算公式為[3]:
(1)
式中,α和ξ分別是無(wú)量綱的裂紋長(zhǎng)度和坐標(biāo)(a和x分別是實(shí)際裂紋長(zhǎng)度和坐標(biāo),W是裂紋體的特征尺寸),σ(ξ)是不考慮裂紋的存在、假想裂紋位置的應(yīng)力。
斷裂力學(xué)中,權(quán)函數(shù)求解法最早由Bueckner提出,Rice對(duì)其做了進(jìn)一步的研究。高精度權(quán)函數(shù)的推導(dǎo)與驗(yàn)證是采用權(quán)函數(shù)法處理裂紋問(wèn)題的關(guān)鍵。本文二維裂紋的權(quán)函數(shù)推導(dǎo)方法主要為:基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法和基于多個(gè)參考載荷情況的權(quán)函數(shù)法。
2.1.1 基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法
Wu等建立了基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法,對(duì)裂紋在一種參考載荷情況下的裂紋面位移分別提出了規(guī)范化的級(jí)數(shù)展開(kāi)表達(dá)式[3]。利用裂紋尖端場(chǎng)的特性、自洽條件、裂紋嘴位移和合理的幾何條件及其各種組合,推導(dǎo)位移級(jí)數(shù)的相關(guān)函數(shù),進(jìn)而由裂紋面位移對(duì)裂紋長(zhǎng)度求偏導(dǎo)數(shù),確定裂紋權(quán)函數(shù)m(α,ξ/α)中的系數(shù)βi(α)。
對(duì)于中心裂紋,該解析權(quán)函數(shù)的通式為:
(2)
(3)
對(duì)于邊緣裂紋,該解析權(quán)函數(shù)的通式為:
(4)
(5)
式中,F(xiàn)i(α)為裂紋面位移展開(kāi)級(jí)數(shù)的系數(shù)。
2.1.2 基于多個(gè)參考載荷情況的權(quán)函數(shù)法
考慮到早期基于裂紋張開(kāi)位移的P-A方法在求解精度和適用范圍方面的缺陷以及在權(quán)函數(shù)的獲取過(guò)程中所需數(shù)值積分的麻煩,一些研究者改用其他方法,即利用多個(gè)參考載荷情況的已知K解來(lái)推導(dǎo)權(quán)函數(shù)。
Glinka和Shen確立了一種適用于所有裂紋幾何的通用權(quán)函數(shù)法[5],該權(quán)函數(shù)的確定是基于兩個(gè)或三個(gè)參考載荷情況及相關(guān)幾何條件,其通用權(quán)函數(shù)基本公式為:
(6)
式中,系數(shù)M1、M2、M3需要基于兩種(或三種)載荷情況下的已知K解(附加幾何條件),通過(guò)求解一組聯(lián)立方程確定。
上述兩種權(quán)函數(shù)法都需要一個(gè)或多個(gè)參考載荷情況下的K值(及裂紋嘴位移)[7]。獲得這些已知條件的途徑,一是從文獻(xiàn)中查找理論精確解或其他方法得到的高精度數(shù)值解,二是采用經(jīng)過(guò)反復(fù)驗(yàn)證的數(shù)值方法(有限元、邊界元等)計(jì)算確定。本文采用有限元計(jì)算參考載荷下的應(yīng)力強(qiáng)度因子。
由參考載荷情況應(yīng)力強(qiáng)度因子導(dǎo)出權(quán)函數(shù)之后,結(jié)合給定應(yīng)力分布下的應(yīng)力強(qiáng)度因子,并與其他方法得出的應(yīng)力強(qiáng)度因子進(jìn)行對(duì)比,是目前常用的評(píng)估應(yīng)力強(qiáng)度因子算法準(zhǔn)確性的一種方式。但是,這種方法實(shí)際上并不能很準(zhǔn)確地反映權(quán)函數(shù)的精度,因?yàn)橛蓹?quán)函數(shù)法計(jì)算應(yīng)力強(qiáng)度因子,在導(dǎo)出權(quán)函數(shù)后,需要對(duì)權(quán)函數(shù)m(a,x)和無(wú)裂紋體假想裂紋處的應(yīng)力σ(x)的乘積沿整個(gè)裂紋面做積分,而積分具有平均效應(yīng),所以由積分得到的應(yīng)力強(qiáng)度因子不能真實(shí)反映權(quán)函數(shù)這一參數(shù)本身的精確性。當(dāng)應(yīng)力分布沿裂紋面劇烈變化,尤其當(dāng)應(yīng)力方向改變時(shí),應(yīng)力強(qiáng)度因子的誤差甚至可能會(huì)遠(yuǎn)大于權(quán)函數(shù)本身的誤差。因此,評(píng)價(jià)權(quán)函數(shù)算法精度的更好的方法是針對(duì)不同裂紋長(zhǎng)度,沿著裂紋面逐點(diǎn)對(duì)權(quán)函數(shù)進(jìn)行比較,即比較格林函數(shù)[8]。這種比較方式不會(huì)受權(quán)函數(shù)與應(yīng)力乘積的積分平均效應(yīng)的影響。
格林函數(shù)G(a,x)又稱(chēng)為影響函數(shù)[9],表示裂紋面在任意位置受到一對(duì)單位集中力P作用時(shí)的無(wú)量綱應(yīng)力強(qiáng)度因子。格林函數(shù)G(a,x)與裂紋面受一對(duì)單位集中力P作用下的應(yīng)力強(qiáng)度因子K沿著集中力作用的位置逐點(diǎn)對(duì)應(yīng),其與權(quán)函數(shù)的關(guān)系如下:
(7)
對(duì)于平板結(jié)構(gòu),工程上常見(jiàn)的裂紋形式為中心裂紋和邊緣裂紋。很多學(xué)者提出了各種求解權(quán)函數(shù)的方法,本文分別采用基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法和基于多個(gè)參考載荷條件的Glinka-Shen通用權(quán)函數(shù)法計(jì)算有限矩形板中心裂紋和邊緣裂紋的應(yīng)力強(qiáng)度因子。
3.1.1 有限矩形板中心裂紋問(wèn)題
考慮有限矩形板中心裂紋問(wèn)題,板寬度為2W、長(zhǎng)度為2H,裂紋長(zhǎng)度為2a(有限矩形板中心裂紋如圖1所示)。選取板半寬W作為該裂紋幾何特征尺寸,無(wú)量綱裂紋長(zhǎng)度α和坐標(biāo)ξ分別為:α=a/W,ξ=x/W。
圖1 有限矩形板中心裂紋
選取裂紋面均布應(yīng)力作為推導(dǎo)中心裂紋權(quán)函數(shù)的參考載荷。對(duì)于裂紋面受均布應(yīng)力的中心裂紋,α=0時(shí),應(yīng)力強(qiáng)度因子和裂紋嘴位移都存在理論精確解;α=1時(shí),應(yīng)力強(qiáng)度因子存在理論極限值。
α=0時(shí):
fr(0)=1Vr(0)=2.0
α=1時(shí):
(8)
對(duì)于有限矩形板中心裂紋,本文利用有限元軟件ANSYS19.0建立模型,模擬不同裂紋長(zhǎng)度的中心裂紋,并進(jìn)行應(yīng)力分析。模擬裂紋長(zhǎng)度α包括0.01、0.05、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9共11個(gè)模型。
為了簡(jiǎn)化計(jì)算,提高運(yùn)行效率,本文建立1/4模型,并施加對(duì)稱(chēng)邊界條件。裂紋尖端處采用1/4節(jié)點(diǎn)奇異單元,網(wǎng)格類(lèi)型設(shè)置單元PLANE183,裂紋尖端網(wǎng)格如圖2所示。材料屬性:楊氏模量E=7.17×105MPa,泊松比v=0.3。
圖2 裂紋尖端網(wǎng)格
利用ANSYS進(jìn)行應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur計(jì)算,結(jié)果如表1所示。
表1 有限元計(jì)算結(jié)果
根據(jù)有限元計(jì)算得到的Kr和Ur值,采用式(9)對(duì)其進(jìn)行正則化處理,得到無(wú)量綱應(yīng)力強(qiáng)度因子fr(α)和無(wú)量綱裂紋嘴張開(kāi)位移Vr(α),然后再通過(guò)多項(xiàng)式擬合,得到fr(α)和Vr(α)的高精度擬合多項(xiàng)式。
Vr(α)=ur(α,ξ=0)E′/(σ0α)
(9)
(10)
式中的多項(xiàng)式系數(shù)λi和γn如表2所示。
表2 均布應(yīng)力下fr(α)和Vr(α)的多項(xiàng)式系數(shù)
把fr(α)和Vr(α)的高精度擬合多項(xiàng)式代入式(11),可以方便地求得?(α)的積分式。將其代入式(3),就能確定裂紋面位移函數(shù)中的Fj(α),j=1~3。
(11)
再將Fj(α)代入式(2),便能確定βi(α)的值。對(duì)βi(α)的計(jì)算結(jié)果進(jìn)行多項(xiàng)式擬合,得到有限矩形板中心裂紋βi(α)的多項(xiàng)式擬合表達(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èn)題
考慮有限矩形板邊緣裂紋問(wèn)題,板寬度為W、長(zhǎng)度為H,裂紋長(zhǎng)度為a(有限矩形板邊緣裂紋如圖3所示)。選取板寬W作為該裂紋幾何特征尺寸,無(wú)量綱裂紋長(zhǎng)度α和坐標(biāo)ξ分別為:α=a/W,ξ=x/W。
圖3 有限矩形板邊緣裂紋
選取裂紋面均布應(yīng)力作為推導(dǎo)邊緣裂紋權(quán)函數(shù)的參考載荷。對(duì)于裂紋面受均布應(yīng)力的邊緣裂紋,α=0時(shí),應(yīng)力強(qiáng)度因子和裂紋嘴位移都存在理論精確解;α=1時(shí),應(yīng)力強(qiáng)度因子存在理論極限值。
α=0時(shí):
fr(0)=1.1215
Vr(0)=2.9086
α=1時(shí):
(14)
對(duì)于有限矩形板邊緣裂紋,本文利用有限元軟件ANSYS19.0建立模型,模擬不同裂紋長(zhǎng)度的邊緣裂紋,并進(jìn)行應(yīng)力分析。模擬裂紋長(zhǎng)度α包括0.01、0.05、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9共11個(gè)模型。
為了簡(jiǎn)化計(jì)算,提高運(yùn)行效率,本文建立1/2模型,并施加對(duì)稱(chēng)邊界條件。裂紋尖端處采用1/4節(jié)點(diǎn)奇異單元,網(wǎng)格類(lèi)型設(shè)置單元PLANE183。材料屬性:楊氏模量E=7.17×105MPa,泊松比v=0.3。
利用ANSYS進(jìn)行應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur計(jì)算,結(jié)果如表3所示。
表3 有限元計(jì)算結(jié)果
根據(jù)有限元計(jì)算得到的Kr和Ur值,根據(jù)式(15)對(duì)其進(jìn)行正則化處理,得到無(wú)量綱應(yīng)力強(qiáng)度因子fr(α)和無(wú)量綱裂紋嘴張開(kāi)位移Vr(α),然后再通過(guò)多項(xiàng)式擬合,得到fr(α)和Vr(α)的高精度擬合多項(xiàng)式。
(15)
式中的多項(xiàng)式系數(shù)λi和γn如表4所示。
表4 均布應(yīng)力下fr(α)和Vr(α)的多項(xiàng)式系數(shù)
把fr(α)和Vr(α)的高精度擬合多項(xiàng)式代入公式,可以方便地求得?(α)的積分式。將其代入式(4),就能確定裂紋面位移函數(shù)中的Fj(α),j=1~3。再將Fj(α)代入式(5),便能確定βi(α)的值。對(duì)βi(α)的計(jì)算結(jié)果進(jìn)行多項(xiàng)式擬合,得到有限矩形板邊緣裂紋βi(α)的多項(xiàng)式擬合表達(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ù)。
本文求解權(quán)函數(shù)的另一種方法是基于多個(gè)參考載荷條件的Glinka-Shen通用權(quán)函數(shù)法。
3.2.1 有限矩形板中心裂紋問(wèn)題
對(duì)于有限矩形板中心裂紋,選取均布應(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)力條件下不同裂紋長(zhǎng)度的中心裂紋,并進(jìn)行應(yīng)力分析。模擬裂紋長(zhǎng)度α包括0.01、0.05、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9共11個(gè)模型。
為了簡(jiǎn)化計(jì)算,提高運(yùn)行效率,本文建立1/4模型,并施加對(duì)稱(chēng)邊界條件。在裂紋尖端處采用1/4節(jié)點(diǎn)奇異單元,其余網(wǎng)格類(lèi)型設(shè)置為單元PLANE183。材料屬性:楊氏模量E=7.17×105MPa,泊松比v=0.3。
利用ANSYS進(jìn)行應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur計(jì)算,結(jié)果如表5所示。
表5 應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur的有限元計(jì)算結(jié)果
將通過(guò)ANSYS求解得到的不同應(yīng)力分布下的應(yīng)力強(qiáng)度因子作為參考應(yīng)力強(qiáng)度因子和對(duì)應(yīng)的應(yīng)力分布函數(shù),代入式(19),并聯(lián)立方程組,可求得不同裂紋長(zhǎng)度對(duì)應(yīng)的M1、M2、M3的值,如表6所示。
(19)
表6 不同裂紋長(zhǎng)度對(duì)應(yīng)的M1、M2、M3的數(shù)值
續(xù)表6
根據(jù)計(jì)算得到的M1、M2、M3的值,對(duì)其進(jìn)行正則化處理,然后再通過(guò)多項(xiàng)式擬合,得到M1、M2、M3的高精度擬合多項(xiàng)式:
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的擬合多項(xiàng)式代入下式,就能確定有限矩形板中心裂紋的解析權(quán)函數(shù),可求解對(duì)應(yīng)載荷作用下的應(yīng)力強(qiáng)度因子。
(21)
3.2.2 有限矩形板邊緣裂紋問(wèn)題
對(duì)于有限矩形板邊緣裂紋,選取均布應(yīng)力、線性應(yīng)力、二次應(yīng)力作為推導(dǎo)邊緣裂紋權(quán)函數(shù)的參考載荷。本文利用有限元軟件ANSYS19.0建立有限矩形板的有限元模型,模擬不同應(yīng)力條件下不同裂紋長(zhǎng)度的邊緣裂紋,并進(jìn)行應(yīng)力分析。模擬的裂紋長(zhǎng)度α包括0.01、0.05、0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9共11個(gè)模型。
為了簡(jiǎn)化計(jì)算,提高運(yùn)行效率,本文建立1/2模型,并施加對(duì)稱(chēng)邊界條件。在裂紋尖端處采用1/2節(jié)點(diǎn)奇異單元,其余網(wǎng)格類(lèi)型設(shè)置為單元PLANE183。有限元模型的材料屬性:楊氏模量E=7.17×105MPa,泊松比v=0.3。
利用ANSYS進(jìn)行應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur計(jì)算,結(jié)果如表7所示。
表7 應(yīng)力強(qiáng)度因子Kr和裂紋嘴位移Ur的有限元計(jì)算結(jié)果
續(xù)表7
將通過(guò)ANSYS求解得到的不同應(yīng)力分布下的應(yīng)力強(qiáng)度因子作為參考應(yīng)力強(qiáng)度因子和對(duì)應(yīng)的應(yīng)力分布函數(shù),代入式(19),并聯(lián)立方程組,可求得不同裂紋長(zhǎng)度對(duì)應(yīng)的M1、M2、M3的值,如表8所示。
表8 不同裂紋長(zhǎng)度對(duì)應(yīng)的M1、M2、M3的數(shù)值
根據(jù)計(jì)算得到的M1、M2、M3的值,對(duì)其進(jìn)行正則化處理,然后再通過(guò)多項(xiàng)式擬合,得到M1、M2、M3的高精度擬合多項(xiàng)式。
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的擬合多項(xiàng)式代入公式,就能確定有限矩形板邊緣裂紋的解析權(quán)函數(shù),可求解對(duì)應(yīng)載荷作用下的應(yīng)力強(qiáng)度因子。
格林函數(shù)G(a,ξ/α)又稱(chēng)為影響函數(shù),表示裂紋面在任意位置受到一對(duì)集中力P作用時(shí)的無(wú)量綱應(yīng)力強(qiáng)度因子。
在有限矩形板中心裂紋的裂紋面任意位置處施加一對(duì)集中力P(單位厚度),計(jì)算格林函數(shù)并繪制函數(shù)的圖像,如圖4所示。
(23)
圖4 中心裂紋的格林函數(shù)
在有限矩形板邊緣裂紋的裂紋面任意位置處施加一對(duì)集中力P(單位厚度),計(jì)算格林函數(shù)并繪制函數(shù)的圖像,如圖5所示。
(24)
圖5 邊緣裂紋的格林函數(shù)
本文通過(guò)不同的權(quán)函數(shù)法計(jì)算有限矩形板中心裂紋和邊緣裂紋的權(quán)函數(shù),主要內(nèi)容如下:
(1)基于裂紋面位移的規(guī)范化解析權(quán)函數(shù)法,是通過(guò)參考載荷下的裂紋面位移Ur(a,x)對(duì)裂紋長(zhǎng)度a求偏導(dǎo)數(shù)推導(dǎo)確定。
(2)使用ANSYS計(jì)算參考應(yīng)力強(qiáng)度因子,確定權(quán)函數(shù)的關(guān)鍵一步是求得參考載荷情況下的高精度裂紋張開(kāi)位移Ur(a,x),然后得到權(quán)函數(shù)的解析表達(dá)式m(α,ξ)。
(3)基于Glinka-Shen通用權(quán)函數(shù)法采用三個(gè)參考應(yīng)力計(jì)算相應(yīng)的參考應(yīng)力強(qiáng)度因子,計(jì)算權(quán)函數(shù)中的待定參數(shù)M1、M2、M3,最終確定應(yīng)力強(qiáng)度因子的權(quán)函數(shù)。
(4)采用格林函數(shù)驗(yàn)證了所求權(quán)函數(shù)的正確性,為飛機(jī)結(jié)構(gòu)和材料的應(yīng)力強(qiáng)度因子確定、疲勞裂紋擴(kuò)展預(yù)測(cè)提供了一種簡(jiǎn)便、高效的方法。