基于Cosserat連續(xù)體理論的疲勞壽命概率可靠性分析

任遠(yuǎn)++尹進(jìn)++張成成++侯乃先++張盛++陳飆松

摘要： 基于Cosserat連續(xù)體理論建立考慮應(yīng)變梯度效應(yīng)的粉末高溫合金彈塑性模型.該模型考慮微尺度長度效應(yīng)，并能消除缺陷局部化問題求解時(shí)的網(wǎng)格依賴性.基于參變量變分原理提出一種高效的求解方法，該方法將原非線性問題轉(zhuǎn)化為互補(bǔ)問題求解，可大大提高求解效率和收斂性.針對一種假想渦輪盤模型，將Cosserat尺度參數(shù)和Basquin斜率作為隨機(jī)變量，分別采用蒙特卡洛法與二階矩法進(jìn)行疲勞壽命概率可靠度計(jì)算.結(jié)果表明二階矩法在保證精度的情況下具有較高的效率.

關(guān)鍵詞： 應(yīng)變梯度理論； 粉末合金； 尺度效應(yīng)； 彈塑性； 參變量變分原理； 渦輪盤； 不確定性； 可靠度方法

中圖分類號： TB33文獻(xiàn)標(biāo)志碼： B

Probabilistic reliability analysis of fatigue life based on

Cosserat continuum theory

REN Yuan1， YIN Jin2， ZHANG Chengcheng1， HOU Naixian1，

ZHANG Sheng2， CHEN Biaosong2

（1. AVIC Commercial Aircraft Engine Co.， Ltd.， Shanghai 201108； 2. State Key Laboratory of Structural Analysis for

Industrial Equipment， Dalian University of Technology， Dalian 116024， Liaoning， China）

Abstract： Considering the strain gradient effect， an elastoplastic model of powder metallurgy superalloy is proposed on the basis of Cosserat continuum theory. The model takes into account the effect of microscale length and can eliminate the mesh dependence when solving the problem of localization of defects. Based on parametric variational principle， an efficient solution method is proposed. The original nonlinear problem is transformed into the complementary problem to solve by the method， which can greatly improve the efficiency and convergence of the solution. As to a hypothetical turbine disk model， the Cosserat scale parameter and the Basquin slope are considered as random variables， and Monte Carlo method and the secondmoment method are used to calculate the probabilistic reliability of fatigue life. The results show that the efficiency of secondmoment method is high in the case of ensuring accuracy.

Key words： strain gradient theory； powder alloy； scale effect； elastoplasticity； parametric variational principle； turbine disk； uncertainty； reliability method

收稿日期： 2017[KG*9〗01[KG*9〗03修回日期： 2017[KG*9〗02[KG*9〗01

基金項(xiàng)目： 國家自然科學(xué)基金（11232003， 91315302）；國家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（“九七三”計(jì)劃）（2010CB832704）；高等學(xué)校學(xué)科創(chuàng)新引智計(jì)劃（B08014）；教育部博士點(diǎn)基金（20130041110050）；上海市科學(xué)技術(shù)委員會(huì)科研計(jì)劃（13DJ1400200）

作者簡介： 任遠(yuǎn)（1982—），男，四川南充人，高級工程師，博士，研究方向?yàn)楹娇瞻l(fā)動(dòng)機(jī)結(jié)構(gòu)可靠性，（Email）renyuan@acae.com.cn0引言

航空發(fā)動(dòng)機(jī)是飛機(jī)的心臟，其中渦輪盤等熱端部件，由于處在高溫的惡劣工作條件下，其故障最多，是影響發(fā)動(dòng)機(jī)可靠性的主要因素.實(shí)際的不確定性因素，一是來自于外部因素的影響，如載荷和環(huán)境等的不確定性；二是來自于內(nèi)部因素的影響，如材料性能的分散性.在起飛—巡航—著陸的多次循環(huán)中，部件材料實(shí)際上經(jīng)受著疲勞作用.隨著航空發(fā)動(dòng)機(jī)的發(fā)展，機(jī)械應(yīng)力越來越高，使用條件越來越惡劣，疲勞問題也越來越突出.因此，疲勞壽命評估和可靠性分析是構(gòu)件設(shè)計(jì)、防止失效、解除對發(fā)動(dòng)機(jī)使用破壞的憂慮所必須進(jìn)行的一項(xiàng)重要工作.相關(guān)學(xué)者[14]進(jìn)行大量的研究工作，提出多種疲勞可靠性分析方法，取得一定的成果.

渦輪盤常用的材料是粉末合金材料，該材料由粉末冶金工藝生產(chǎn)，常存在缺陷.制造過程中產(chǎn)生的缺陷與傳統(tǒng)的鑄鍛合金的缺陷有所不同[5]，如熱誘導(dǎo)孔洞、原始顆粒邊界、夾雜物與異金屬等，其對疲勞性能具有很大的影響.熱誘導(dǎo)孔洞由殘留氬氣或氦氣不溶于合金產(chǎn)生.原始顆粒邊界作為裂紋源是由粉末冶金工藝過程產(chǎn)生的，并且在熱處理過程中很難消除，嚴(yán)重影響疲勞壽命.夾雜物與異金屬嚴(yán)重降低粉末合金的疲勞壽命和力學(xué)性能，包括塑性，其中主要有合金中的夾雜物以及粉末中的夾雜物.

傳統(tǒng)連續(xù)介質(zhì)力學(xué)的研究對象的共同特點(diǎn)是其部件或變形所涉及的尺度較大，常大于若干毫米.近年來連續(xù)介質(zhì)力學(xué)開始在其他工程中獲得新的應(yīng)用，例如材料微結(jié)構(gòu)、微電子、微納米薄膜、微傳感器、顆?；蚶w維夾雜復(fù)合材料微加工等.這些領(lǐng)域應(yīng)用的共同特點(diǎn)是部件或變形所涉及的特征尺寸很小，約0.1～10.0 μm，即所謂的細(xì)觀尺度.在細(xì)觀尺度下，材料的力學(xué)性質(zhì)呈現(xiàn)出很強(qiáng)的尺寸效應(yīng)，而應(yīng)變梯度理論是能夠很好地解釋尺度效應(yīng)的有效方法之一.常用的商業(yè)有限元軟件不適合模擬細(xì)觀尺度下的力學(xué)現(xiàn)象，是因?yàn)闆]有考慮特征長度尺度.近年來，研究者通過引進(jìn)不同長度量綱的參數(shù)，提出多種應(yīng)變梯度理論.[6]

1基于Cosserat理論的彈塑性模型

采用Cosserat模型[7]，其材料點(diǎn)除具有平動(dòng)自由度ui外，還引入獨(dú)立的旋轉(zhuǎn)自由度ωi及與其對應(yīng)的微曲率和偶應(yīng)力.應(yīng)變向量ε依賴于位移向量和旋轉(zhuǎn)向量，含有18個(gè)應(yīng)變分量，即εij=uj，i+ejikωk

χij=ωj，i （1）式中：ejik為置換向量；當(dāng)i=j時(shí)χij為扭轉(zhuǎn)微曲率分量，當(dāng)i≠j時(shí)χij為彎曲微曲率分量.合并應(yīng)變和曲率向量，得到廣義應(yīng)變向量εc=εχT，同樣，廣義應(yīng)力向量包含應(yīng)力和偶應(yīng)力，即σc=σmT.Cosserat連續(xù)體模型平衡方程表達(dá)為σji，j+bi=0

mij，i+ejikσik+ψj=0 （2）式中：bi為體力；ψj為體力偶.本構(gòu)關(guān)系表示為σ=Duuε

m=Dωωχ

σc=Dεc （3）彈性剛度陣D=D1000

0D200

00D20

000D20

0D3 （4）其中：D1=λ+2Gλλ

λλ+2Gλ

λλλ+2G （5）

D2=G+GcG-Gc

G-GcG+Gc （6）

D3=4G·diag l2tl2tl2tl2cl2cl2cl2cl2cl2c （7）式中：λ為拉梅常數(shù)；G為剪切模量；Gc為Cosserat剪切模量；lt和lc分別為Cosserat材料扭轉(zhuǎn)和彎曲的內(nèi)部特征長度尺度，可通過實(shí)驗(yàn)測定，本文取lt=lc.邊界條件可以表示為σijni=j，on Γp

ui=i，on Γu

mijni=j，on Γm

ωi=i，on Γω （8）采用壓力無關(guān)的von Mises屈服準(zhǔn)則模擬不可逆的變形.應(yīng)力應(yīng)變關(guān)系表示為σ=D：（ε-εp） （9）式中：σ和ε分別為應(yīng)力和應(yīng)變張量；εp為塑性應(yīng)變；D為彈性材料矩陣.塑性流動(dòng)方程為dεpij=λgσij （10）式中：g為塑性勢；σij為流動(dòng)因子.屈服函數(shù)f滿足一致性條件f≤0，λ≥0，f·λ=0 （11）2參變量變分原理

采用von Mises屈服準(zhǔn)則及關(guān)聯(lián)的流動(dòng)法則f（σij，εij，p，κ）=σeq-σy（εij，p，κ） （12）式中：κ為材料硬化/軟化內(nèi)變量；σeq為有效應(yīng)力，表示為σeq=32sijsij1/2

sij=σij-13σkkδij

σm=13σkk （13）σy為屈服應(yīng)力，對各向同性線性硬化材料，有σy=σ0+ph （14）式中：σ0為初始屈服應(yīng)力；p為等效塑性應(yīng)變；h為硬化模量.

對本構(gòu)關(guān)系線性化，屈服函數(shù)進(jìn)行一階展開得f=f0+fσijdσij+fεij，pdεij，p+fκdκ=

f0+fσijDijkldεkl+-fσijDijklgσkl+

fεij，pgσij+fκhλ（15）利用關(guān)系dκ=hλ，引入非負(fù)松弛因子v得c1kldεkl+u1λ+w1+v=0，

vλ=0，λ≥0，v≥0 （16）其中： c1kl=fσijDijkl，w1=f0，u1=-fσijDijklgσkl+fεpijgσij+fκh （17）定義勢能泛函Π=∫Ω12Dijkldεijdεkl-λDijklgσijdεkldΩ-

∫ΩdbiduidΩ+∫ΓdiduidΓ （18）則參變量變分原理為min Π（dεij，dui，λ）

s.t. c1kldεkl+u1λ+w1+v=0

vλ=0，λ≥0，v≥0 （19）由于 dσkl=Dijkldεij-Dijklλgσij （20）

對泛函變分得δΠ=∫Ωdσijδ（dεij）dΩ-

∫Ωdbiδ（dui）dΩ+∫Γdiδ（dui）dΓ （21）分部積分得δΠ=-∫Ω（dσij，j+dbi）δ（dui）dΩ+

∫Γ（dσijnj-di）δ（dui）dΓ （22）令變分等于0，并考慮δ（dui）的任意性，可得平衡方程和應(yīng)力邊界條件.進(jìn)一步，由于δ2Π≥0，則泛函取總體最小值.由于泛函是在狀態(tài)控制方程控制下取極值的，因此整個(gè)過程中也滿足塑性流動(dòng)定律.

對方程進(jìn)行離散，du=Nudu⌒，dε=Bd （21）式中：u為單元內(nèi)位移向量；u⌒為節(jié)點(diǎn)位移向量；Nu為位移插值形函數(shù)；B為應(yīng)變陣.

假定狀態(tài)方程在單元上積分滿足，得U1λ+C1d-d1+v=0

Φ1λ+Kd-dP=0

vλ=0，λ≥0，v≥0 （22）其中：

U1=∫Ω-fσTDgσ+fεpTgσ+fκhdΩ，

CT1=∫ΩBTDfσdΩ，d1=-∫Ω f0dΩ，

Φ1=-∫ΩBTDgσdΩ，K=∫ΩBTDBdΩ，

dP=∫ΩNTudbdΩ+∫ΓNTudpdΓ.

上述線性互補(bǔ)問題可采用自由變量二次規(guī)劃算法進(jìn)行求解.[78]

3疲勞壽命分析

3.1假想渦輪盤結(jié)構(gòu)應(yīng)力分析

假想的渦輪盤結(jié)構(gòu)見圖1.材料的彈性模量為200 GPa，泊松比為0.3，密度為7 850 kg/m3，屈服極限為500 MPa，塑性硬化切線模量為10 GPa，Cosserat連續(xù)體剪切模量為50 GPa，受到均勻轉(zhuǎn)速所導(dǎo)致的離心力載荷，結(jié)構(gòu)內(nèi)側(cè)受到約束作用，轉(zhuǎn)速工況為2 100 r/min.

考慮不同尺度參數(shù)對彈塑性分析應(yīng)力結(jié)果的影響，尺度參數(shù)取為0.05～0.25 m，不同尺度參數(shù)的最大von Mises應(yīng)力見表1.a）渦輪盤結(jié)構(gòu)網(wǎng)格

b）渦輪盤結(jié)構(gòu)約束

Tab.1Maximum von Mises stress of different scale parameters尺度參數(shù)/m0.0500.0750.1000.1250.1500.1750.2000.2250.250最大von Mises應(yīng)力/MPa588.745587.355585.463583.187580.649578.16575.566572.939570.357

由表1可以發(fā)現(xiàn)：尺度參數(shù)與最大von Mises應(yīng)力呈近似線性關(guān)系，尺度參數(shù)越大，最大von Mises應(yīng)力越小，其中，尺度參數(shù)為0.150 m的位移、應(yīng)力、應(yīng)變結(jié)果分別見圖2～4.

4.2疲勞壽命計(jì)算

考慮不同尺度參數(shù)對疲勞特性的影響，針對特定工況下Cosserat連續(xù)體理論彈塑性分析計(jì)算出的應(yīng)力結(jié)果計(jì)算疲勞壽命.在循環(huán)載荷中，考慮完全加載的循環(huán)載荷形式.采用基于應(yīng)力的BrownMiller算法[910]，并采用Goodman平均應(yīng)力修正公式.

BrownMiller算法認(rèn)為疲勞損傷是剪應(yīng)變和正應(yīng)變的耦合造成的，計(jì)算公式為γmax2+εN2=1.65σ′fE（2Nf）b+1.75ε′f（2Nf）c（24）式中：γmax/2為最大剪切應(yīng)變幅值；εN/2為正應(yīng)變幅值；σ′f為拉伸疲勞強(qiáng)度系數(shù)；E為彈性模量；b為Basquin斜率；ε′f為疲勞韌性系數(shù)；c為疲勞韌性指數(shù)；Nf為疲勞壽命循環(huán)次數(shù).

基于應(yīng)力的BrownMiller算法的損傷參數(shù)用應(yīng)力表示為τmax2+σN2=E1.65σ′fE（2Nf）b+1.75ε′f（2Nf）c （25） 假設(shè)疲勞極限Sf取為500 MPa，Basquin斜率b取為-0.138 3，SN曲線見圖5.

同樣考慮尺度參數(shù)取為0.050～0.250 m，不同尺度參數(shù)的臨界循環(huán)次數(shù)見表2.

考慮不同的Basquin斜率對疲勞壽命的影響.疲勞極限Sf取500 MPa，Cosserat連續(xù)體尺度參數(shù)l取0.150 m，Basquin斜率b取為-0.16～-0.12，不同斜率b的臨界循環(huán)次數(shù)見表3.

表 2不同尺度參數(shù)的臨界循環(huán)次數(shù)

Tab.2Critical cycle number under different scale parameters尺度參數(shù)/m0.0500.0750.1000.1250.1500.1750.2000.2250.250壽命/（106 次）1.601.641.691.761.851.932.032.132.23

表 3不同斜率b的臨界循環(huán)次數(shù)

Tab.3Critical cycle number of different slope b斜率b-0.16-0.15-0.14-0.13-0.12壽命/（106 次）2.322.111.881.661.43

由表3可以發(fā)現(xiàn)：尺度參數(shù)越大，壽命越大；b越大，壽命越小.需針對這2個(gè)參數(shù)討論概率可靠性問題.

4概率可靠性分析

針對給定載荷工況的假想渦輪盤算例，考慮尺度參數(shù)為隨機(jī)變量，采用二次多項(xiàng)式回歸方法，對不同尺度參數(shù)得到的壽命數(shù)據(jù)進(jìn)行擬合，可以得到擬合后的方程，擬合曲線與原始數(shù)據(jù)見圖6.

考慮安全的循環(huán)壽命為1.7×106，則可得到極限狀態(tài)方程Z，Z>0為安全.Z=N-1.700 0×106=7.168 8×106×l2+

1.076 0×106×l-0.179 2×106 （26）考慮l為隨機(jī)變量，假設(shè)l的均值數(shù)學(xué)期望值為0.150，標(biāo)準(zhǔn)差為0.033 3，變異因數(shù)為0.222.圖 6二次多項(xiàng)式擬合曲線與原始數(shù)據(jù)對比

Fig.6Comparison of quadratic polynomial fitting curve and original data

采用蒙特卡洛方法得到的Z>0的安全概率為93.38%，采用二階矩的可靠度方法[11]得到的Z>0的安全概率為93.35%，可以發(fā)現(xiàn)2種方法相差不大，說明二階矩的可靠度方法的合理性.

針對給定載荷工況的渦輪盤算例，考慮b為隨機(jī)變量，采用一次多項(xiàng)式回歸方法，尺度參數(shù)取0.150 m時(shí)的應(yīng)力結(jié)果，對不同斜率b得到的壽命數(shù)據(jù)進(jìn)行擬合，可以得到擬合后的方程，擬合曲線與原始數(shù)據(jù)見圖7.

Fig.7Comparison of linear polynomial

fitting curve and original data

考慮安全的循環(huán)壽命為1.7×106，則可得到極限狀態(tài)方程Z，Z>0為安全.Z=N-1.700×106=

-2.230×107×b-2.942×106 （27）考慮b為隨機(jī)變量，假設(shè)b的均值數(shù)學(xué)期望為-0.140，標(biāo)準(zhǔn)差為0.003 33，變異因數(shù)為0.023 8.

采用蒙特卡洛方法得到的Z>0的安全概率為99.24%，采用二階矩的可靠度方法得到的Z>0的安全概率為99.23%，可以發(fā)現(xiàn)2種方法相差不大，說明二階矩的可靠度方法合理.

同時(shí)考慮l和b為隨機(jī)變量，根據(jù)壽命數(shù)據(jù)，采用二元二次多項(xiàng)式進(jìn)行回歸分析，可以得到擬合后的方程，擬合曲線與原始數(shù)據(jù)見圖8.

考慮安全的循環(huán)壽命為1.7×106，則可得到極限狀態(tài)方程Z，Z>0為安全.

Z=N-1.7×106=

p00-1.7×106+p10l+p01b+p20l2+p11lb（28）

其中，

p00=-1.291×106，p10=-9.839×105，

p01=2.040×107，p20=7.457×106，

p11=1.400×107.

考慮l和b為隨機(jī)變量，假設(shè)l的均值數(shù)學(xué)期望值為0.150，標(biāo)準(zhǔn)差為0.033 3，變異系數(shù)為0.222，假設(shè)b的均值數(shù)學(xué)期望值為-0.140，標(biāo)準(zhǔn)差為0.003 33，變異因數(shù)為0.023 8.

采用蒙特卡洛方法得到的Z>0的安全概率為93.25%，采用二階矩的可靠度方法得到的Z>0的安全概率為93.58%，可以發(fā)現(xiàn)2種方法相差不大，但二階矩方法的計(jì)算量要遠(yuǎn)小于蒙特卡洛方法，說明二階矩的可靠度方法合理.在實(shí)際問題中，需要考慮各種參數(shù)對疲勞壽命的敏感性，如果考慮的隨機(jī)變量不只2個(gè)甚至很多，可以采用代理模型的方法[12]構(gòu)建極限狀態(tài)方程，并以此模型進(jìn)行可靠度計(jì)算.

6結(jié)束語

基于Cosserat連續(xù)體理論可以通過尺度參數(shù)描述粉末合金的細(xì)觀微結(jié)構(gòu)的特征，因此建立Cosserat連續(xù)體理論的彈塑性模型.采用參變量變分原理構(gòu)建求解算法，這種算法要比傳統(tǒng)的求解算法具有更好的收斂性和穩(wěn)定性，算例表明材料的尺度參數(shù)會(huì)影響結(jié)構(gòu)的應(yīng)力，進(jìn)而影響結(jié)構(gòu)的疲勞壽命，可以同時(shí)考慮材料的不確定性和強(qiáng)度的不確定性來進(jìn)行疲勞壽命概率可靠性分析.本文將尺度參數(shù)和Basquin斜率作為隨機(jī)變量，并用代理模型和回歸分析的方法構(gòu)建極限狀態(tài)方程，采用蒙特卡洛方法和二階矩方法進(jìn)行可靠度計(jì)算，可以說明這種求解策略的合理性.

參考文獻(xiàn)：

[1]王應(yīng)珍， 呂文林. 用攝動(dòng)法研究幾何尺寸為隨機(jī)量時(shí)輪盤的可靠性[J]. 航空動(dòng)力學(xué)報(bào)， 1994， 9（4）： 371374.

WANG Y Z， LYU W L. Perturbation technique for reliability analysis of a rotating disk with uncertain geometry[J]. Journal of Aerospace Power， 1994， 9（4）： 371374.

[2]陸山， 呂文林. 輪盤低周疲勞壽命可靠性分析新方法[J]. 航空學(xué)報(bào)， 1997， 18（2）： 135138.

LU S， LYU W L. New method for reliability analysis of discs LCF life[J]. Acta Aeronautica et Astronautica Sinica， 1997， 18（2）： 135138.

[3]樂曉斌. 用剩余疲勞損傷強(qiáng)度預(yù)測機(jī)械零件疲勞可靠度及壽命的方法[J]. 機(jī)械強(qiáng)度， 1996， 18（1）： 3436.

LE X B. The method for calculating fatigue reliability and life of mechanical parts using residual fatigue damage strength[J]. Journal of Mechanical Strength， 1996， 18（1）： 3436.

[4]高陽， 白廣忱， 張瑛莉. 渦輪盤多軸低循環(huán)疲勞壽命可靠性分析[J]. 航空學(xué)報(bào)， 2009， 30（9）： 16781682.

GAO Y， BAI G C， ZHANG Y L. Reliability analysis of multiaxial low cycle fatigue life for turbine disk[J]. Acta Aeronautica et Astronautica Sinica， 2009， 30（9）： 16781682.

[5]張義文， 上官永恒. 粉末高溫合金的研究與發(fā)展[J]. 粉末冶金工業(yè)， 2004， 14（6）： 3043.

ZHANG Y W， SHANGGUAN Y H. Research and development in P/M superalloy[J]. Powder Metallurgy Industry， 2004， 14（6）： 3043.

[6]陳少華， 王自強(qiáng). 應(yīng)變梯度理論進(jìn)展[J]. 力學(xué)進(jìn)展， 2003， 33（2）： 207216.

CHEN S H， WANG Z Q. Advances in strain gradient theory[J]. Advances in Mechanics， 2003， 33（2）： 207216.

[7]張洪武. 參變量變分原理與材料和結(jié)構(gòu)力學(xué)分析[M]. 北京： 科學(xué)出版社， 2010： 32105.

[8]鐘萬勰， 張洪武， 吳承偉. 參變量變分原理及其在工程中的應(yīng)用[M]. 北京： 科學(xué)出版社， 1997： 138205.

[9]MILLER K J. Fatigue under complex stress[J]. Metal Science， 1977， 11（89）： 432438. DOI： 10.1179/msc.1977.11.89.432.

[10]BROWN M W， MILLER K J. A theory for fatigue failure under multiaxial stressstrain conditions[J]. Proceedings of Institution of Mechanical Engineers， 1973， 187： 745755. DOI： 10.1243/PIME_PROC_1973_187_069_02.

[11]姚衛(wèi)星. 結(jié)構(gòu)疲勞壽命分析[M]. 北京： 國防工業(yè)出版社， 2003： 187220.

[12]李堅(jiān). 代理模型近似技術(shù)研究及其在結(jié)構(gòu)可靠度分析中的應(yīng)用[D]. 上海： 上海交通大學(xué)， 2013： 2240.（編輯武曉英）