徐 越 王躍方 李 聰 魏學(xué)敏 李盛文

（1.中國航空工業(yè)空氣動力研究院；2.低速高雷諾數(shù)氣動力航空科技重點實驗室；3.大連理工大學(xué)）

0 引言

當(dāng)前透平機械朝著高端化、節(jié)能、高可靠性的方向發(fā)展，遠程運行監(jiān)測、智能化控制及故障自愈化成為未來發(fā)展的趨勢。葉片是工業(yè)軸流式壓縮機的核心部件，承擔(dān)功能轉(zhuǎn)換的重要任務(wù)，運轉(zhuǎn)時既承受離心力引起的高平均應(yīng)力，又承受各種激振載荷引起的交變應(yīng)力。當(dāng)轉(zhuǎn)-靜子剮蹭或高階模態(tài)被激振時，葉片邊緣材料的初始缺陷易萌生微小裂紋，在轉(zhuǎn)動過程中不斷發(fā)展，形成宏觀裂紋直至使葉片斷裂。為保證旋轉(zhuǎn)葉片長周期使役安全，研究初始裂紋萌生后的裂紋尖端應(yīng)力場，判斷應(yīng)力強度因子是否超出材料斷裂韌性而導(dǎo)致葉片斷裂，對軸流葉片設(shè)計具有重要參考價值。

化工及空分裝置中使用的軸流式壓縮機葉片一般為薄壁金屬件，離心力引起的拉應(yīng)力在葉片橫截面上的分布較為均勻，葉片受力接近于薄板平面應(yīng)力狀態(tài)。因此，可參照均質(zhì)平板斷裂力學(xué)分析方法計算葉片的應(yīng)力強度因子[1]。對于含中心裂紋、處于平面應(yīng)力狀態(tài)的線彈性矩形薄板，中外學(xué)者大多采用權(quán)函數(shù)法推導(dǎo)裂紋尖端的I型應(yīng)力強度因子理論解[2-5]。近年來，樊俊鈴等[6-8]基于權(quán)函數(shù)法研究了平板含單邊裂紋、雙邊裂紋和中心裂紋的應(yīng)力強度因子，給出其與裂紋長度、位置、轉(zhuǎn)速等參數(shù)之間的關(guān)系。用權(quán)函數(shù)法計算應(yīng)力強度因子需要工程人員具備一定的斷裂力學(xué)基礎(chǔ)和解析推導(dǎo)能力。在數(shù)值解方面，尹奇志等[9]用有限元法計算了含裂紋平板的應(yīng)力強度因子，表明采用仿真技術(shù)求解此類問題具有較好的精度。姜偉等[10]將動態(tài)有限元分析和相互作用積分相結(jié)合，計算了含裂紋平板受沖擊載荷時的動態(tài)應(yīng)力強度因子。大量研究表明，擴展有限元法(XFEM)在分析裂紋等不連續(xù)性問題給出了理想的計算結(jié)果。李錄賢等[11]詳細介紹了該方法并提供了應(yīng)用實例，對以往數(shù)值結(jié)果作了修正。董玉文、余天堂等[12-13]研究表明擴展有限元法可以準確地求解應(yīng)力強度因子。盡管擴展有限元法在求解裂紋擴展問題具有較大優(yōu)勢，但對于建模復(fù)雜、工況繁多等問題分析周期較長，仿真分析很難滿足工程上對應(yīng)力強度因子預(yù)測的效率要求。因此，探索如何有效地把理論分析與數(shù)值方法結(jié)合起來，通過少量仿真分析修正理論解，從而實現(xiàn)針對不同工況下的裂紋參數(shù)，快速、準確地計算應(yīng)力強度因子。

本文開展壓縮機軸流葉片單邊貫穿型裂紋應(yīng)力強度因子的計算方法研究。將離心力作用下的薄壁金屬葉片簡化為矩形平板，在應(yīng)力強度因子權(quán)函數(shù)解法的基礎(chǔ)上，通過擴展有限元數(shù)值分析修正理論解中的形狀因子，并研究裂紋位置和板厚對理論解的影響，通過數(shù)據(jù)擬合提出計算葉片裂紋應(yīng)力強度因子的半解析公式。進一步建立實際壓縮機葉片模型，基于擴展有限元仿真分析修正基于薄板的半解析公式，建立實際葉片I 型應(yīng)力強度因子的半解析解公式。本文將應(yīng)力強度因子的經(jīng)典理論解適用范圍從中心型裂紋推廣到非中心型裂紋，綜合考慮葉片厚度和形狀的影響，為快速評估實際軸流壓縮機旋轉(zhuǎn)葉片的應(yīng)力強度因子提供新方法，也為后續(xù)葉片剩余強度及壽命分析提供分析基礎(chǔ)。

1 矩形板I型裂紋應(yīng)力強度因子

軸流壓縮機葉片示意圖如圖1所示，由于葉片厚度遠遠小于長度和寬度，故可將其簡化為xy面內(nèi)、以角速度w旋轉(zhuǎn)的矩形平板，簡化模型如圖2所示。設(shè)板的長度為h，寬度為b，厚度為W，板的底邊坐標y=-h/2。考慮在平面應(yīng)力狀態(tài)下產(chǎn)生的裂紋以I型（張開型）為主，在左側(cè)板邊處(x=-b)設(shè)置了貫穿板厚度的長度為a的裂紋。與其他研究不同的是：本文設(shè)置的裂紋不限于板的中心位置（y=0），即起裂位置可以是任意值y=l-h/2-r，且旋轉(zhuǎn)軸不是板底邊，而是平行于x軸的任意AB(y=-h/2-r)軸。

圖1 軸流壓縮機葉片示意圖Fig.1 Schematic diagram of axial flow compressor blade

圖2 含單邊貫穿型裂紋矩形板Fig.2 Rectangular plate of single-edged penetrating crack

依據(jù)斷裂力學(xué)理論，裂紋尖端的I型應(yīng)力強度因子的一般形式為[14]：

其中，σy為裂紋尖端名義應(yīng)力，F(xiàn)I為形狀因子。

當(dāng)葉片旋轉(zhuǎn)時，裂紋所在橫截面處的離心力為：

式中，ρ為材料質(zhì)量密度。于是，裂紋面上任意點x 處的拉應(yīng)力為：

由此得到旋轉(zhuǎn)葉片上裂紋應(yīng)力強度因子：

對于含中心裂紋、平面應(yīng)力狀態(tài)的線彈性矩形薄板，大多采用權(quán)函數(shù)法推導(dǎo)I型裂紋尖端的應(yīng)力強度因子理論解。基于此方法，在板中面(z=0)上由名義應(yīng)力σy(x)產(chǎn)生的I型裂紋應(yīng)力強度因子表達式為[14]：

其中，m(x，a)為權(quán)函數(shù)，其與載荷條件無關(guān)，由裂紋體結(jié)構(gòu)自身的性質(zhì)決定，且具有唯一性。E為材料彈性模量；vy(x，a)為外載荷引起的裂紋面張開位移。根據(jù)文獻[3]，權(quán)函數(shù)取為x∈(0,a)內(nèi)一致收斂的級數(shù)

式中，Mn為待定系數(shù)。根據(jù)文獻[3]，級數(shù)展開式一般取到n=3即可達到精度要求。

在理論上，若已知互相獨立的3個應(yīng)力場對應(yīng)的應(yīng)力強度因子，則可根據(jù)式(5)、式(7)確定M1；M2；M3。對于中心型裂紋，由于位移場的對稱性，裂紋面的位移場曲率為零，即

可知M2為常數(shù)，不影響應(yīng)力強度因子解，M1；M3表達式為：

其中，F(xiàn)1和F2是矩形板在遠場均布拉力和純彎曲兩個獨立載荷作用下的形狀因子，通過解析推導(dǎo)或數(shù)值計算得到。

對于受離心力作用的細長軸流葉片，可使用Brown等[15]提出的適于有限寬、無限長板條的公式計算形狀因子F1。即

其中，a/b＜0.6 的中心型裂紋，上式計算誤差不超過0.5%。也可采用Tada等[16]提出的理論解計算形狀因子F1，即

對于a/b＜0.2 的中心型短裂紋，上式的計算誤差不超過1%；對更長的裂紋，其誤差不超過0.5%。

若純彎加載情形，可使用文獻[15]提出的公式計算形狀因子F2，表達式為

對于a/b＜0.6 的中心型裂紋，上式計算誤差不超過0.5%。將式(10)、(12)依次代入式(9)、(7)，再積分式，即可求出應(yīng)力強度因子KI的理論解。

需要指出的是：上述理論解的前提條件是中心型裂紋假設(shè)，而實際裂紋可能萌生在葉片的任意位置，此時的權(quán)函數(shù)應(yīng)與式(7)有所不同。此外，實際葉片具有一定厚度，裂紋尖端是一條曲線而不是平面情形下的單個點。因此，上述應(yīng)力強度因子理論解在計算實際軸流葉片時存在一定誤差，需對其進行驗證并做必要的修正。

2 應(yīng)力強度因子理論解修正

本節(jié)基于ABAQUS軟件的擴展有限元法(XFEM)模塊，選取不同的裂紋參數(shù)分析單邊貫穿型裂紋的I型應(yīng)力強度因子，驗證形狀因子并對其修正，討論裂紋位置對理論解精度的影響。

擴展有限元法是當(dāng)前數(shù)值上解決不連續(xù)力學(xué)問題的最佳選擇之一，已廣泛應(yīng)用于含裂紋結(jié)構(gòu)的應(yīng)力強度因子、J積分等斷裂力學(xué)性能分析。它在常規(guī)位移場插值的基礎(chǔ)上引入了Heaviside函數(shù)描述裂紋兩側(cè)的間斷位移，以及裂尖函數(shù)反映裂紋尖端的應(yīng)力奇異[11]。擴展有限元方法優(yōu)于傳統(tǒng)有限元方法的一個重要方面是裂紋面無需創(chuàng)建在單元的邊界，即裂紋獨立于網(wǎng)格，并且無須隨著裂紋的擴展而重新劃分網(wǎng)格。在用擴展有限元計算應(yīng)力強度因子時，積分路徑原則上由裂紋下表面開始，沿封閉圍道環(huán)繞裂紋尖端，終止于裂紋上表面。本文在實際計算時嘗試了不同的圍道、富集半徑、裂尖網(wǎng)格密度，最終選出精度滿意的參數(shù)組合，保證了應(yīng)力強度因子解的收斂性。

2.1 形狀因子驗證

采用ABAQUS 軟件，建立帶裂紋的矩形平板仿真模型如圖3(a)所示，模型包括50898 個節(jié)點，16765 個單元。如圖3(b)所示，在裂紋尖端附近布置長度為1mm的奇應(yīng)變?nèi)切螁卧?，采用環(huán)形區(qū)域積分，設(shè)置積分路徑為4條，選擇最大主應(yīng)力準則為裂紋初始化準則。

圖3 帶裂紋的矩形板仿真模型Fig.3 Simulation model of a rectangular plate with crack

令板寬b=100mm，彈性模量E=200GPa，板頂邊和底邊受100MPa 均布拉力作用。給定50 組不同的板長度h和裂紋長度a，根據(jù)ABAQUS 仿真分析得到應(yīng)力強度因子KI，再通過式反推出裂紋的形狀因子。依據(jù)數(shù)值分析得到的形狀因子與理論解對比，結(jié)果如圖4所示。

圖4 形狀因子隨裂紋長度的變化Fig.4 Variation of crack shape factor with crack length

從圖4可以看出，形狀因子數(shù)值解與式預(yù)測趨勢一致，均隨裂紋延伸而單調(diào)增大。其中，長寬比h/b越小，形狀因子仿真值與理論值的偏差越大，對于長度大于0.05b 的裂紋，當(dāng)h/b＜0.8 時，仿真值與理論值的偏差已不可忽略。這種情形多見于軸流壓縮機的高壓級短葉片，在以往文獻中對此問題研究不夠充分。

基于上述分析結(jié)果，對形狀因子的理論解進行修正。本文采用LM法和通用全局優(yōu)化算法（UGO），擬合形狀因子的多項式解如下：其中，α=a/b，β=h/b，c0-c9為系數(shù)。對于本例，它們依次為：1.9472、9.1361、-4.4591、-1.9440、-14.3208、6.3277、9.7379、1.3786、6.5427、-2.6999。擬合結(jié)果與仿真結(jié)果的相關(guān)系數(shù)為0.998，說明式(13)具有較高精度。形狀因子的函數(shù)曲面如圖5所示。

圖5 形狀因子關(guān)于裂紋長度及板長寬比的曲面圖Fig.5 Surface diagram of the crack shape factor with crack length and plate’s length to width ratio

2.2 裂紋位置對應(yīng)力強度因子解的影響

設(shè)板長度h=180mm，寬度為b=100mm，裂紋位置l分別取10mm、36mm、60mm 和90mm，以5mm 為間隔計算不同裂紋長度對應(yīng)的應(yīng)力強度因子理論解和有限元解，結(jié)果如圖6所示。

圖6 不同裂紋位置的應(yīng)力強度因子數(shù)值解和理論解Fig.6 Numerical and theoretical solutions of stress intensity factor at different crack locations

由圖6(b)、(c)和(d)中可見，裂紋位置分別為l=36mm、60mm 和90mm 時，應(yīng)力強度因子理論解與擴展有限元解吻合得較好。圖6(a)中l(wèi)=10mm 時，即裂紋足夠接近板底邊時，理論解與有限元解偏差較大，原因在于理論解建立在假設(shè)裂紋處于平板中心線上，計算精度隨裂紋位置距離中心線越遠而逐漸下降。

根據(jù)以上分析，采用擬合算法對式的應(yīng)力強度因子理論解進行修正。以裂紋長度和位置為參數(shù)，用LM法和通用全局優(yōu)化算法（UGO）擬合出應(yīng)力強度因子的半解析解如下：

對本算例，系數(shù)d0～d9依次為189.2635、18.8499、6.3932、-0.1711、0.4155、-0.2272、0.01228、-0.00349、-0.00403、0.00166。上述擬合的相關(guān)系數(shù)為0.999，說明精度滿足要求。

2.3 板厚對應(yīng)力強度因子解的影響

分別取板厚W=6mm、10mm、20mm 和40mm，并分別設(shè)置長度a=10mm、30mm 的裂紋。由于裂紋尖端應(yīng)力場關(guān)于板中面對稱分布，故圖7只展示從上表面到中面的I型應(yīng)力強度因子XFEM解。其中橫坐標為沿厚度方向的無量綱尺寸2z/W，圖中水平線代表應(yīng)力強度因子的理論解。

圖7 應(yīng)力強度因子沿板厚度的變化Fig.7 Variation of the stress intensity factor along plate thickness

由圖7可見，應(yīng)力強度因子隨板厚的變化量大致不超過1%。在厚度方向上，中面(z=0)的應(yīng)力強度因子值最大，上、下表面的因子值最小，與實際裂紋尖端接近于圓弧的規(guī)律相符。

考慮板厚度影響后，可對式的應(yīng)力強度因子進行修正，得到受離心力作用的矩形平板I型應(yīng)力強度因子半解析解：

依據(jù)仿真分析結(jié)果，修正系數(shù)dw在1.06至1.09之間波動，為保守起見，對實際軸流葉片建議取dw=1.1。

3 軸流葉片應(yīng)力強度因子的半解析解

實際壓縮機軸流葉片雖是薄壁結(jié)構(gòu)，但其安裝角沿葉片高度方向不斷變化，從葉根至葉頂呈現(xiàn)出一定的空間扭曲，且橫截面的厚度中間大、兩邊小。因此，基于平板公式計算的應(yīng)力強度因子必然存在誤差。以下針對含單邊貫穿裂紋的實際壓縮機軸流葉片，基于擴展有限元法對半解析公式進行進一步修正以提高其預(yù)測精度。

仿真分析對象為圖1 所示的某軸流壓縮機第6 級葉片，寬度b=100mm，高寬比為1.8。有限元模型網(wǎng)格尺寸為5mm，在葉片側(cè)邊及葉根處，網(wǎng)格尺寸加密為1mm，模型共劃分127867 個單元。計算時分別在葉片上部和底邊兩個不同位置（l=36mm和l=90mm）設(shè)置不同長度的貫穿裂紋。圖8給出了I型應(yīng)力強度因子隨裂紋擴展的變化趨勢，同時給出了按式(15)算出的半解析解結(jié)果。

圖8 實際葉片應(yīng)力強度因子隨裂紋長度的變化Fig.8 Variation of actual blade stress intensity factor with crack length

從圖中看到，半解析解與擴展有限元解變化趨勢基本一致，但數(shù)值上存在一定偏差，原因是實際葉片橫截面積沿葉高方向逐漸減小，裂紋面實際承受的拉應(yīng)力小于式計算值，對應(yīng)的應(yīng)力強度因子值必然小于等截面薄板情況。此外，由于葉片形狀扭曲，受力后實際處于拉伸、剪切和扭轉(zhuǎn)共存的組合應(yīng)力狀態(tài)，裂紋前緣變形復(fù)雜，同時存在I型、II型和III型應(yīng)力強度因子，與推導(dǎo)理論解時假設(shè)的平面應(yīng)力狀態(tài)差別較大，故必然產(chǎn)生一定偏差。

考慮到完好葉片的寬、高、厚度均已確定，因此影響斷裂力學(xué)性能的主要參數(shù)僅為裂紋長度及起裂位置。引入修正系數(shù)ζ(α，β)，其中β=l/h為無量綱裂紋位置，將式（16）的半解析解改寫為

根據(jù)圖7 的擴展有限元仿真結(jié)果，采用LM 法和UGO 法擬合出修正系數(shù)ζ(α，β)關(guān)于裂紋長度和裂紋位置的響應(yīng)面。本例設(shè)置了75個(α，β)組合，得到修正系數(shù)ζ(α，β)的響應(yīng)面如圖9 所示，其中黑點為原始數(shù)據(jù)點。對應(yīng)的函數(shù)表達式為:

圖9 修正系數(shù)隨裂紋長度和裂紋位置的變化Fig.9 The variation of modified parameter with crack lengths and crack locations

在本文算例中，表達式(17)的系數(shù)g0～g9依次是：0.19988,-0.05989,1.52846,2.26096,-0.22049,-4.93579,-2.30420,-0.34714,0.34261,5.30313。由此可得實際壓縮機軸流葉片應(yīng)力強度因子的半解析解，式（16）可用于計算不同位置、不同長度的貫穿型裂紋尖端應(yīng)力強度因子。

分析圖9 中曲面的變化特征可知，當(dāng)裂紋長度a不變時，若改變裂紋位置l，則修正系數(shù)ζ(α，β)隨裂紋位置l的增大而增大；當(dāng)裂紋位置l不變時，改變裂紋長度a，則修正系數(shù)ζ(α，β)隨裂紋長度a的增加而迅速增大，表明了基于數(shù)值仿真來修正半解析解的必要性。

盡管本例葉片的裂紋長度和裂紋位置的數(shù)據(jù)離散度均較大，但應(yīng)力強度因子半解析解與擴展有限元法解的相對誤差仍在5%至10%之間，表明精度滿足工程要求。

4 結(jié)論

本文針對軸流壓縮機行業(yè)快速評估含裂紋葉片應(yīng)力強度因子的實際需求，根據(jù)結(jié)構(gòu)和載荷特點建立了含單邊貫穿性裂紋的矩形平板模型，結(jié)合斷裂力學(xué)理論解公式，采用擴展有限元方法研究了形狀因子、板厚及裂紋位置對I 型裂紋應(yīng)力強度因子理論解精度的影響，通過擬合提出了改進后的半解析解公式。進一步開展實際軸流葉片的仿真分析，修正基于平板模型的半解析解公式，構(gòu)建了用于計算包含不同位置、不同長度的貫穿型裂紋的實際軸流葉片應(yīng)力強度因子的半解析解。主要結(jié)論如下：

1）依據(jù)斷裂力學(xué)平板理論預(yù)測軸流葉片裂紋尖端的I 型應(yīng)力強度因子存在一定誤差，對于高壓級短葉片，形狀因子影響不可忽略；裂紋位置距離平板中心線越遠，計算精度越低；應(yīng)力強度因子隨板厚的變化量大致不超過1%。

2）實際軸流葉片裂紋尖端的I 型應(yīng)力強度因子與基于平板模型的計算結(jié)果存在差異，裂紋越長、裂紋距離板中心線越遠，計算誤差越大。

3）本文提出的實際軸流葉片I 型應(yīng)力強度因子的半解析解，考慮了實際葉片厚度及葉片形狀變化、不同裂紋位置和裂紋長度的影響，能夠滿足工程對計算精度的要求，與理論解和仿真分析相比，可以滿足一般的工程使用需求,有效提升軸流葉片應(yīng)力強度因子的評估精度及效率，節(jié)省斷裂力學(xué)分析所需的工作量。

4）采用代理模型技術(shù)及優(yōu)化設(shè)計算法有望進一步改善解的預(yù)測精度，此外，未來還可探索代理模型、深度學(xué)習(xí)等技術(shù)在本研究方向的工程應(yīng)用。

致謝：本文得到國家自然科學(xué)基金及國家科技重大專項項目資助，關(guān)萍同志參加了有關(guān)分析計算以及審稿人對論文寫作提出了寶貴意見，在此一并致謝。