基于損傷力學(xué)和XFEM的艦船蒸汽輪機(jī)葉片裂紋擴(kuò)展研究

李亞政，王 崢，孫亞非，張廣輝，黃洪雁，劉占生

（1.哈爾濱工業(yè)大學(xué) 能源科學(xué)與工程學(xué)院，哈爾濱 150001；2.上?？臻g推進(jìn)研究所，上海 200433）

0 引 言

葉片作為船用汽輪機(jī)的核心部件，其疲勞壽命的長(zhǎng)短直接影響整個(gè)艦船動(dòng)力裝置的穩(wěn)定可靠運(yùn)行；若能在復(fù)雜工作環(huán)境中大致預(yù)測(cè)葉片的疲勞壽命，并在危險(xiǎn)點(diǎn)出現(xiàn)前及時(shí)檢修，將能在很大程度上減少事故的發(fā)生。目前，傳統(tǒng)疲勞壽命分析方法[1]在計(jì)算疲勞壽命的過(guò)程中大多要根據(jù)實(shí)際經(jīng)驗(yàn)進(jìn)行修正，難以給出詳細(xì)的裂紋具體擴(kuò)展過(guò)程和具體發(fā)生破壞的循環(huán)次數(shù)，而連續(xù)介質(zhì)損傷力學(xué)則能夠較好地解決這一問(wèn)題。

損傷力學(xué)是研究材料在各種不同的加載條件下，其損傷的演化發(fā)展直至破壞的過(guò)程。材料疲勞破壞的實(shí)質(zhì)是材料在循環(huán)應(yīng)力下的性能發(fā)生惡化而不斷減損，其間斷裂面的有效接觸面積逐漸減小，材料的彈性模量等參數(shù)減小，間接造成應(yīng)力、應(yīng)變?cè)龃?。目前損傷力學(xué)的基本理論已經(jīng)越來(lái)越多地應(yīng)用于材料的疲勞壽命評(píng)估。陳沙古等人[3]結(jié)合低周疲勞試驗(yàn)數(shù)據(jù)，在簡(jiǎn)化的橢圓形裂紋形狀擴(kuò)展模型的基礎(chǔ)上提出了一種基于損傷演化的疲勞裂紋擴(kuò)展壽命簡(jiǎn)化計(jì)算方法。張鼎等[4]將風(fēng)暴模型和裂紋擴(kuò)展率單一曲線模型及焊趾表面裂紋應(yīng)力強(qiáng)度因子的計(jì)算方法結(jié)合起來(lái)，探討了復(fù)雜載荷作用下船舶結(jié)構(gòu)疲勞裂紋擴(kuò)展預(yù)報(bào)方法。

本文基于損傷力學(xué)的基本理論，推導(dǎo)葉片低周疲勞的壽命模型，同時(shí)利用擴(kuò)展有限元方法進(jìn)行葉片裂紋擴(kuò)展模擬和低周疲勞壽命分析，對(duì)比理論計(jì)算和數(shù)值模擬的結(jié)果，驗(yàn)證所提出的壽命模型的準(zhǔn)確性。

1 循環(huán)載荷下的損傷力學(xué)模型

1.1 損傷數(shù)學(xué)模型

Lemaitre[5]文中提出的耗散能理論認(rèn)為：疲勞損傷過(guò)程往往伴隨著能量的流動(dòng)以及損失，其間發(fā)生的材料性能劣化以及塑性應(yīng)變均為不可逆過(guò)程，所以損傷的實(shí)質(zhì)是不可逆的能量耗散過(guò)程，可以用耗散勢(shì)函數(shù)將其表示出來(lái)，最終求得可以表達(dá)疲勞損傷過(guò)程的損傷本構(gòu)方程。

在損傷演化的過(guò)程中，起主要作用的有以下幾個(gè)狀態(tài)變量，溫度、應(yīng)力、彈性應(yīng)變、塑性應(yīng)變、總應(yīng)變、累積塑性應(yīng)變以及損傷變量等。假設(shè)材料存在Helmholtz自由能ψ，且由于損傷與塑性不耦合，則其表達(dá)式為：

其中：ε為總應(yīng)變；εe為彈性應(yīng)變；εp為塑性應(yīng)變；γ表示累積塑性應(yīng)變；α為背應(yīng)變；D為損傷變量。類比斷裂力學(xué)中的能量釋放率，定義損傷變量的伴隨變量應(yīng)變能釋放率，其表達(dá)式為：

其中：σeq為Mises等效應(yīng)力，Rv為三軸應(yīng)力比，表達(dá)式為：

1.2 損傷演化方程及其修正

對(duì)于一個(gè)穩(wěn)定的循環(huán)過(guò)程，考慮在材料受損的過(guò)程中損傷變量的變化，在此式中引入損傷變量D，其中塑性應(yīng)變的微分表達(dá)式為：

選取損傷的增量形式進(jìn)行推導(dǎo)。對(duì)于單次循環(huán)而言，塑性累積增量越大，表示材料受損的程度越高。損傷增量與應(yīng)變能釋放率Y成正比。損傷增量D˙表示為：

Lemaitre在推導(dǎo)損傷演化方程時(shí)，沒(méi)有考慮材料在發(fā)生斷裂時(shí)彈性模量的變化。而現(xiàn)今大量實(shí)驗(yàn)數(shù)據(jù)表明，鋼材料在進(jìn)入塑性區(qū)時(shí)彈性模量會(huì)下降6%-10%左右，所以將彈性模量用=E·（1- ）D代替。將（2）式與（4）式代入（5）式，同時(shí)得到塑性應(yīng)變?cè)隽康谋磉_(dá)式：損傷增量定義為單次循環(huán)加載內(nèi)損傷的增加量，則將上式寫成微分形式為：

對(duì)于單次循環(huán)內(nèi)的損傷，在積分式要考慮加載（0-1積分）和卸載（1-0積分）的過(guò)程[6]，因此積分（6）式得：

（7）式中應(yīng)力部分是針對(duì)于對(duì)稱循環(huán)而言，對(duì)于非對(duì)稱循環(huán)，在求解等效應(yīng)力時(shí)考慮平均應(yīng)力的影響，在（7）式中引入平均應(yīng)力項(xiàng)加以修正，得到下式：

根據(jù)（8）式可以快速地計(jì)算出單次循環(huán)內(nèi)損傷的增量。隨著循環(huán)次數(shù)的不斷增加，損傷變量D不斷變大，危險(xiǎn)部位的應(yīng)力值不斷增大，所以單次循環(huán)內(nèi)損傷變化量也逐漸增大。與Morrow疲勞壽命公式不同的是，Morrow公式只能粗略估計(jì)材料總的疲勞壽命，而損傷變量的計(jì)算公式則體現(xiàn)出真實(shí)的材料失效過(guò)程，即損傷不是隨著循環(huán)呈線性變化，而是呈類似指數(shù)增長(zhǎng)的趨勢(shì)。

在得到單次循環(huán)內(nèi)損傷增量后，對(duì)損傷變量從0到1積分，即可得到低周疲勞壽命的表達(dá)式：

對(duì)于同一種材料，循環(huán)應(yīng)變硬化系數(shù)n′、循環(huán)強(qiáng)度系數(shù)K′、楊氏模量等不隨著構(gòu)件形狀、載荷大小而變化，而不同載荷直接影響的是應(yīng)力幅值σeq以及材料常數(shù)S。

為了得到損傷隨循環(huán)次數(shù)的變化過(guò)程，在積分的過(guò)程中，從0到D積分，可導(dǎo)出D的表達(dá)式為：

（9）（10）ND。過(guò)低周疲勞試驗(yàn)確定材料的低周循環(huán)疲勞壽命，再基于損傷力學(xué)的基本公式進(jìn)行反推[5]，即可得到材料常數(shù)S，在無(wú)法滿足實(shí)驗(yàn)條件的情況下也可以通過(guò)（6）式曲線的斜率來(lái)確定等。在已知各項(xiàng)材料參數(shù)時(shí)，可以通過(guò)有限元分析，得到裂紋位置的應(yīng)力狀態(tài)，即可得到計(jì)算材料的低周疲勞壽命。

上述的低周疲勞壽命公式是在損傷變量由0到1積分得到的，即包括疲勞裂紋萌生以及疲勞裂紋擴(kuò)展的壽命，但在試件的實(shí)際工作中，尤其是汽輪機(jī)的葉片，在非外物損傷的情況下，損傷變量達(dá)到1即完全斷裂開(kāi)是非常危險(xiǎn)的，所以應(yīng)該在適當(dāng)?shù)膿p傷情況下及時(shí)更換受損葉片。大量的試驗(yàn)表征，材料的疲勞裂紋擴(kuò)展的壽命僅僅占總壽命的20%左右，損傷力學(xué)中將疲勞裂紋萌生和裂紋擴(kuò)展的界限處的損傷變量值Dc稱為損傷門檻臨界值，加載過(guò)程中當(dāng)損傷變量達(dá)到Dc時(shí)，一般可以認(rèn)為材料發(fā)生破壞。本文葉片的材料為2Cr13，劉曉燕等通過(guò)疲勞試驗(yàn)測(cè)得其損傷門檻值約為0.094[7]。

2 基于損傷力學(xué)的裂紋擴(kuò)展模擬

2.1 擴(kuò)展有限元方法

裂紋位置處的應(yīng)變不連續(xù)、幾何界面處的位移不連續(xù)性和裂紋尖端的奇異性使得疲勞裂紋成為典型的不連續(xù)問(wèn)題。求解不連續(xù)問(wèn)題最有效的數(shù)值方法就是擴(kuò)展有限元方法（XFEM），其適用于任意邊界條件、幾何形狀和幾何非線性問(wèn)題等，因而成為分析裂紋等不連續(xù)問(wèn)題的主要手段。XFEM所使用的網(wǎng)格與結(jié)構(gòu)內(nèi)部的幾何或物理界面無(wú)關(guān)，克服了在裂紋尖端等高應(yīng)力和變形集中區(qū)進(jìn)行細(xì)密的網(wǎng)格劃分所帶來(lái)的困難，同時(shí)在模擬裂紋擴(kuò)展時(shí)也無(wú)需對(duì)網(wǎng)格進(jìn)行重新劃分。

XFEM方法通過(guò)對(duì)已被裂紋分割開(kāi)的單元和裂紋尖端所處的單元上的節(jié)點(diǎn)位移函數(shù)引入富集項(xiàng)加以表征，其表達(dá)式為[8]：

式中：第一項(xiàng)代表未受裂紋影響的單元的位移形函數(shù)，其中I是網(wǎng)格單元中所有的離散點(diǎn)集，N（x）表示傳統(tǒng)有限元單元內(nèi)的形函數(shù)， ｛u0I,v0I｝T表示未被分割開(kāi)的位移函數(shù)；第二項(xiàng)代表被裂紋完全穿透的單元的位移形函數(shù)，J是該區(qū)域內(nèi)的節(jié)點(diǎn)集，如圖1中空心節(jié)點(diǎn)位置，｛a1j,a2j｝T為裂紋位置單元修正后的節(jié)點(diǎn)自由度；H（x）為單位階躍函數(shù)，模擬穿過(guò)裂紋面的位移跳躍，其表達(dá)式為：

圖1 裂紋及裂紋尖端單元節(jié)點(diǎn)圖Fig.1 Crack and crack tip unit node map

完整的XFEM方程如下：

式中：Kuq、Muq是剛度矩陣中兩種自由度的耦合項(xiàng)，F(xiàn)、V是與兩種節(jié)點(diǎn)自由度對(duì)應(yīng)的合外力，從方程中可見(jiàn)，在XFEM方法計(jì)算中，最重要的是計(jì)算出兩種節(jié)點(diǎn)自由度。

2.2 ABAQUS低周疲勞分析技術(shù)

ABAQUS在低周疲勞分析方面有獨(dú)特之處，它通過(guò)Direct cyclic分析直接計(jì)算結(jié)構(gòu)穩(wěn)態(tài)的循環(huán)響應(yīng)，利用傅里葉級(jí)數(shù)及材料非線性行為的時(shí)間積分獲得穩(wěn)定的循環(huán)求解。該方法的基礎(chǔ)是構(gòu)造一個(gè)位移函數(shù)描述結(jié)構(gòu)在整個(gè)循環(huán)載荷周期內(nèi)所有時(shí)間點(diǎn)的響應(yīng)。低周疲勞分析基礎(chǔ)是損傷起始標(biāo)準(zhǔn)的判斷。當(dāng)循環(huán)中的結(jié)構(gòu)響應(yīng)穩(wěn)定時(shí)，使用每個(gè)循環(huán)的累積非彈性滯后能△ω和材料常數(shù)來(lái)確定開(kāi)始損傷的循環(huán)次數(shù)。一旦在材料滿足損傷起始標(biāo)準(zhǔn)時(shí)，則需要基于穩(wěn)定循環(huán)的非彈性滯后能來(lái)計(jì)算和更新?lián)p傷狀態(tài)。

由于裂紋尖端不可能停留在單元內(nèi)，只能在邊界上，所以對(duì)于穩(wěn)定循環(huán)N結(jié)束時(shí)如果單元的任一裂紋尖端均滿足裂紋擴(kuò)展標(biāo)準(zhǔn)，則將從當(dāng)前循環(huán)向前延伸△N個(gè)周期，通過(guò)壓裂使裂紋長(zhǎng)度向前延伸至少一個(gè)單元，當(dāng)前循環(huán)的損傷變量DN也增大，新的損傷變量DN+△N由下式給出：

式中：L為單元節(jié)點(diǎn)相關(guān)的特征長(zhǎng)度，c3和c4為材料常數(shù)。

通過(guò)指定損傷外推公差△Dtol，使用自適應(yīng)算法基于下式來(lái)確定損傷外推的周期數(shù)：

如果最大損傷增量大于指定的損傷外推公差，則向前推出損傷的循環(huán)次數(shù)相應(yīng)減少，以確保最大損傷增量小于損傷外推公差。反之如果所有材料點(diǎn)的最大損傷增量小于指定的損傷外推公差的一半，則循環(huán)次數(shù)相應(yīng)增加，以確保最大損傷增量等于損傷外推公差。

圖2 預(yù)置裂紋的葉片模型Fig.2 The blade model with pre-crack

2.3 艦船蒸汽輪機(jī)葉片裂紋擴(kuò)展數(shù)值模擬

在損傷起始標(biāo)準(zhǔn)中，選用POWER準(zhǔn)則來(lái)計(jì)算臨界應(yīng)變能釋放率，其表達(dá)式為：

其中：Gequiv為等效應(yīng)變能釋放率。 GIC、GIIC、GIIIC為沿三個(gè)方向的剪切模態(tài)斷裂能量。arn、ars和art為各冪指數(shù)。

本文選取實(shí)際艦船蒸汽輪機(jī)葉片來(lái)進(jìn)行裂紋擴(kuò)展模擬，并在葉片上預(yù)置裂紋，如圖2所示。

葉片材料選擇2Cr13，其主要參數(shù)如表1所示。

參考艦船蒸汽輪機(jī)葉片實(shí)際的工作過(guò)程，不考慮溫度引起的熱應(yīng)力的影響，對(duì)葉片施加旋轉(zhuǎn)載荷條件，旋轉(zhuǎn)半徑414 mm，大小為2 500弧度/分析步，步長(zhǎng)設(shè)置為3 000步。計(jì)算結(jié)果如圖3所示。

表1 2Cr13材料參數(shù)Tab.1 2Cr13 material parameters

圖3 始末狀態(tài)的應(yīng)力云圖Fig.3 Stress cloud charts of beginning and end

裂紋擴(kuò)展過(guò)程中的應(yīng)力云圖和塑性應(yīng)變累積過(guò)程如圖4所示。

圖4（a）為葉片裂紋擴(kuò)展過(guò)程中的應(yīng)力變化云圖，從圖中可以看出應(yīng)力集中基本發(fā)生在已完全裂開(kāi)單元的上下表面和裂紋尖端，當(dāng)循環(huán)達(dá)到1 075次時(shí)，由于此時(shí)裂紋長(zhǎng)度較長(zhǎng)，葉片有效承載面積較小，整個(gè)葉身的有效應(yīng)力處于較高水平。圖4（b）為累積塑性應(yīng)變?cè)茍D，同應(yīng)力云圖一樣，塑性應(yīng)變區(qū)集中在裂紋上下表面以及裂紋尖端，隨著葉片損傷程度的增加，在保持循環(huán)內(nèi)外載荷不變的情況下，塑性應(yīng)變量也大幅增加。

圖4 （a）應(yīng)力云圖；（b）累計(jì)塑形應(yīng)變?cè)茍DFig.4 (a)Stress cloud charts;(b)Cumulative plastic strain cloud

3 結(jié)果分析

ABAQUS模擬裂紋擴(kuò)展時(shí)每次向前推進(jìn)至少一個(gè)單元，且裂紋尖端一般都在網(wǎng)格界面處，因此，對(duì)于整個(gè)葉片的損傷變量可以通過(guò)已斷裂的網(wǎng)格面積和葉片橫截面積的比值來(lái)確定。表2和圖5給出了裂紋擴(kuò)展模擬時(shí)損傷變量與循環(huán)次數(shù)的對(duì)應(yīng)關(guān)系。

表2 N次循環(huán)后的損傷值Tab.2 The damage value after N times cycle

圖5 損傷隨循環(huán)次數(shù)的變化曲線Fig.5 Damage curve with the number of cycles

由于載荷較大，當(dāng)材料進(jìn)入塑性區(qū)時(shí)，隨著塑性的逐漸積累，ABAQUS在每一個(gè)循環(huán)加載后會(huì)出現(xiàn)由于塑性引起的殘余應(yīng)力場(chǎng)，這使得循環(huán)特性的值不為0，間接影響平均應(yīng)力的大小。

選擇初始循環(huán)中裂紋尖端位置的前緣單元的應(yīng)力狀態(tài)作為理論計(jì)算的基礎(chǔ)。峰值應(yīng)力出現(xiàn)在載荷達(dá)到最大時(shí)，大小為533 MPa，谷值應(yīng)力出現(xiàn)在載荷卸載后的殘余應(yīng)力場(chǎng)，大小為403 MPa，同時(shí)可計(jì)算出三軸應(yīng)力Rv為1.79。

根據(jù)損傷增量的表達(dá)式來(lái)確定材料常數(shù)S，即：

通過(guò)ABAQUS輸出材料隨循環(huán)加載得到的累積塑性應(yīng)變，得出S≈89 MPa。

將上述值代入（9）式中，可得葉片的低周疲勞壽命約1 132次。將計(jì)算結(jié)果與ABAQUS模擬結(jié)果對(duì)比如圖6所示。

圖6 損傷過(guò)程變化曲線Fig.6 Damage process curve

從上圖可以看出，數(shù)值模擬曲線與理論計(jì)算曲線結(jié)果基本相近，誤差約為5%左右。同時(shí)曲線變化的趨勢(shì)大致相同，在損傷變量小于0.1之前，材料的損傷程度隨著循環(huán)次數(shù)呈現(xiàn)近似比例增長(zhǎng)，且斜率較小。達(dá)到損傷門檻值前后時(shí)，損傷的增長(zhǎng)速率逐漸增加，直到損傷變量達(dá)到0.3左右時(shí)，較高的應(yīng)力集中程度使得材料進(jìn)入嚴(yán)重塑性區(qū)，葉片失去作用而迅速斷裂，此時(shí)模型曲線也接近垂直于橫軸。

4 結(jié) 論

本文基于損傷力學(xué)的基本理論給出了葉片低周疲勞壽命的計(jì)算方法，考慮實(shí)際葉片的受力情況以及材料塑性應(yīng)變的影響，在公式中加入平均應(yīng)力進(jìn)行修正。同時(shí)基于XFEM方法模擬葉片裂紋擴(kuò)展過(guò)程，分析了低周疲勞壽命，得出以下結(jié)論。

（1）本文給出了計(jì)算葉片低周疲勞壽命的數(shù)學(xué)模型，并在考慮葉片實(shí)際非對(duì)稱循環(huán)的工作狀態(tài)下，引入平均應(yīng)力的概念對(duì)模型進(jìn)行了修正。

（2）使用該模型計(jì)算低周疲勞壽命的結(jié)果與數(shù)值模擬的結(jié)果相比基本接近，誤差約在5%。

（3）材料損傷程度達(dá)到損傷門檻值0.094前，單次循環(huán)內(nèi)損傷增量較小，損傷變量隨循環(huán)次數(shù)變化幅值較小，此階段大致占總壽命的80%；損傷達(dá)到臨界值以后，單次循環(huán)內(nèi)增量迅速提高，材料進(jìn)入嚴(yán)重的塑性變形區(qū)，進(jìn)而迅速失效，表現(xiàn)為葉片斷裂。

參 考 文獻(xiàn)：

[1]朱寶田,劉東遠(yuǎn).計(jì)算汽輪機(jī)葉片低周和高周疲勞壽命的改進(jìn)的局部應(yīng)力應(yīng)變方法—第一部分 模型與方法[J].熱力發(fā)電,1997(04)：9-14.Zhu Baotian,Liu Dongyuan.Improved local stress-strain method for calculating low cycle and high cycle fatigue life of steam turbine blades-part I models and methods[J].Thermal Power Generation,1997(04)：9-14.

[2]朱寶田,劉東遠(yuǎn).計(jì)算汽輪機(jī)葉片低周和高周疲勞壽命的改進(jìn)的局部應(yīng)力應(yīng)變方法—第二部分 計(jì)算和分析[J].熱力發(fā)電,1997(05)：3-8.Zhu Baotian,Liu Dongyuan.Improved local stress-strain method for calculating low cycle and high cycle fatigue life of steam turbine blades-part II calculation and analysis[J].Thermal Power Generation,1997(05)：3-8.

[3]陳沙古,黃進(jìn)浩,萬(wàn)正權(quán),卞如岡.基于損傷演化的裂紋擴(kuò)展壽命簡(jiǎn)化計(jì)算方法[J].船舶力學(xué),2014,18(07)：815-820.Chen Shagu,Huang Jinhao,Wan Zhengquan,Bian Rugang.Research on calculation method of crack propagation life based on damage evolution[J].Journal of Ship Mechanics,2014,18(07)：815-820.

[4]張 鼎,黃小平,崔維成.一種復(fù)雜載荷作用下船舶結(jié)構(gòu)疲勞裂紋擴(kuò)展預(yù)報(bào)方法[J].船舶力學(xué),2015,19(05)：541-552.Zhang Ding,Huang Xiaoping,Cui Weicheng.A procedure to predict fatigue crack growth of ship structures under complex loading condition[J].Journal of Ship Mechanics,2015,19(05)：541-552.

[5]Lemaitre J.損傷力學(xué)教程[M].北京：科學(xué)技術(shù)出版社,1996.Lemaitre J.Damage mechanics tutorial[M].Beijing：Science and Technology Press,1996.

[6]周勝田.航空發(fā)動(dòng)機(jī)葉片疲勞的損傷力學(xué)研究及外物損傷影響[D].沈陽(yáng)：東北大學(xué),2008.Zhou Shengtian.Study on fatigue of air-engine blades by damage mechanics and infulence of object damage[D].Shenyang：Northeastern University,2008.

[7]劉曉燕,何曉梅,董 潔.2Cr13鋼的疲勞裂紋萌生與擴(kuò)展行為[J].熱加工工藝,2012(02)：51-54.Liu Xiaoyan,He Xiaomei,Dong Jie.Fatigue crack initiation and propagation behavior of 2Cr13 steel[J].Hot Working Technology,2012(02)：51-54.

[8]孫 勇.基于Monte-Carlo法的燃?xì)廨啓C(jī)壓氣機(jī)葉片可靠性分析[D].哈爾濱：哈爾工業(yè)大學(xué),2015.Sun Yong.Reliability analysis of gas turbine compressor blade based on Monte-Carlo method[D].Harbin：Harbin University of Technology,2015.