榫連接結構微動疲勞壽命不確定性分析

陳 實,張屹尚,張嘉誠,徐逸晗,孫 凱

(1.中國航發商用航空發動機有限責任公司,上海 200241;2上海飛機發動機工程技術中心,上海 200241)

0 引言

榫連接結構是航空發動機葉片與輪盤連接的主要形式,在航空發動機飛行剖面轉速工況下承受巨大的交變離心載荷,使得微動疲勞損傷成為其主要失效機制。而航空發動機全生命周期存在幾何、載荷[1]、材料、監測誤差等多種不確定性因素[2],這些因素導致發動機在試驗、服役過程中性能偏離設計值,帶來性能不達標甚至引發故障。因此針對榫連接結構微動疲勞壽命開展不確定性分析具有重要的意義。

蒙特卡洛法基于大數定律,適用于各種極限狀態方程形式、變量維數、變量分布特征,廣泛應用于航空發動機結構強度可靠性評估[3]。崔海濤等[4]以彈性模量、摩擦系數和壽命預測模型參數為隨機變量,通過104次蒙特卡洛模擬,得到了不同可靠度下的燕尾榫結構微動疲勞壽命;同時,這一方法被用于燕尾榫結構高周微動疲勞壽命可靠性計算,然后基于Miner疲勞累積損傷準則,建立了高低周微動疲勞壽命對數正態分布模型[5]。此外,蒙特卡洛法也應用于航空發動機其他關鍵部件的研究,例如渦輪盤持久壽命分析[6],且蒙特卡洛抽樣次數達到106。因此,對于高安全性要求的航空發動機限壽件失效風險評估問題,其失效屬于小概率事件,傳統蒙特卡洛方法在模擬時存在計算成本高、效率不足的問題。為此,工程上發展了基于抽樣的傳統蒙特卡洛模擬方法及其改進算法[7],與傳統蒙特卡洛法計算時間成本相比,重要性抽樣方法能夠將時間降低約1個數量級。

目前,國內關于航空發動機榫連接結構微動疲勞壽命不確定性的研究相對較少。為此,本文中針對航空發動機榫連接結構微動疲勞壽命不確定性分析這一工程問題,提出了改進蒙特卡洛抽樣方法,其應用貝塔分布抽樣和雙峰正態分布抽樣對樣本進行抽樣,然后基于重要抽樣思想推導得到變量原始分布特征的失效概率。結合非線性算例分析了分布參數、縮放參數和平移參數對改進蒙特卡洛抽樣方法失效概率評估性能的影響。進一步地,針對航空發動機風扇葉片榫連接結構微動疲勞壽命不確定性分析問題,采用雙峰正態分布抽樣完成了失效概率計算及全局靈敏度分析。

1 改進蒙特卡洛抽樣

蒙特卡洛法利用計算機生成服從變量空間隨機分布的抽樣樣本,然后進行隨機模擬試驗,從而根據試驗結果獲得系統失效概率。設X=(X1,X2,…,Xn)為系統變量參數,其分布概率密度函數記為fX(X),功能函數記為g(x),則系統失效概率Pf可計算為

(1)

其中If[g(x)]為功能函數g(x)的指示函數,計算如下:

(2)

蒙特卡洛法借助計算機產生服從分布fX(x)的樣本xi,然后模擬系統特性,則式(1)系統失效概率Pf的估計值為

(3)

從式(3)可知產生服從變量分布的隨機樣本為蒙特卡洛法以及改進的數字模擬方法的基礎[8]。式(3)所示的簡單隨機抽樣廣泛應用于工程可靠性分析。但是在航空發動機零部件參數不確定性分析中,設計點變量通常為對應分布均值[1,10],屬于小概率失效問題,簡單隨機抽樣存在效率低、樣本量巨大的問題。針對這一問題,研究人員提出了重要抽樣法,其通過采用重要抽樣密度函數代替變量分布密度函數,使得落入失效域的樣本概率增加,從而提高失效概率預估精度和收斂速度[9]。在重要抽樣法中,式(1)可計算為

(4)

式(4)中:hX(x)為重要抽樣密度函數,此時失效概率估計值為:

(5)

式(5)中:xi不再服從變量分布概率密度函數fX(x),而是根據重要抽樣密度函數hX(x)隨機生成。

(6)

1.1 貝塔分布抽樣

貝塔分布是一組定義在(0,1)區間的連續概率分布,其概率密度函數計算如下

(7)

式中:α、β為貝塔分布的2個參數,B(α,β)計算為

(8)

α、β均小于1時呈現中心概率密度小、邊緣概率密度大的特點,均大于1時正相反,而均為1時則轉化為均勻分布。

針對簡單隨機抽樣在失效區域樣本分布不足問題,將貝塔分布函數構造為重要抽樣密度函數。引入縮放參數kB,將貝塔分布區間其擴增到目標變量區間,使抽樣樣本xi能夠囊括失效區域,計算如下

(9)

(10)

1.2 雙峰正態分布抽樣

雙峰正態分布由2個正態分布疊加得到,其表達式為

(11)

式中:μ1、μ2、σ1、σ2分別為2個正態分布的均值和標準差,r為第一個正態分布所占的比重,本文取σ1=σ2=σ,r=0.5。抽樣密度如圖1所示。為了增加分布樣本在失效區域的抽樣概率,引入平移參數kD,使得

圖1 雙峰正態分布曲線

(12)

式中:μ為設計點,則式(11)為重要抽樣密度函數hX(x)。

1.3 失效概率評估流程

失效概率評估流程包括固定抽樣次數和依據收斂準則2種流程。前者按照抽樣次數N直接進行隨機模擬試驗,然后根據試驗結果評估失效概率。相比于簡單隨機抽樣,重要抽樣法在計算失效概率時考慮了抽樣密度函數的影響。

按照收斂準則的失效概率評估流程如圖2所示,收斂判據為前后兩次預測失效概率的相對誤差連續5次小于收斂指標。初始抽樣次數為10,樣本池樣本為空集,按照抽樣次數進行蒙特卡洛模擬,計算樣本池的預測失效概率,當預測失效概率滿足收斂指標時收斂次數加1;當預測失效概率不滿足收斂指標時,收斂次數歸零,然后按照抽樣次數進行蒙特卡洛模擬,并將模擬樣本放入樣本池中,在樣本池數目達到抽樣次數的10倍時,則應更改抽樣次數為樣本池數目,如此往復,直至連續5次滿足收斂指標后輸出預測失效概率。

圖2 依據收斂準則的失效概率評估流程

2 抽樣性能分析

2.1 算例1

計算案例采用西北工業大學呂震宙等[8]提出的非線性功能函數,如下所示:

(13)

式中:X1、X2相互獨立,且X1服從正態分布N(10,25),X2服從正態分布N(9.9,25),屬于小概率失效。

針對上述功能函數,采用簡單隨機抽樣、貝塔分布抽樣和雙峰正態分布抽樣3種抽樣方法,設置蒙特卡洛抽樣次數分別為[10,102,103,104,105,106,107,108],開展100輪蒙特卡洛試驗。

2.1.1簡單隨機抽樣

圖3 簡單隨機抽樣預測結果

2.1.2貝塔分布抽樣

影響貝塔抽樣密度函數及其失效概率評估的因素包括分布參數α、β以及縮放參數kB,通過隨機模擬,分析這些參數對失效概率預估性能的影響。

1) 貝塔分布參數。設置縮放參數kB為8,對應樣本分布范圍為(μ-4σ,μ+4σ)。分布參數α=β且分別取值為0.75、1.00、1.25、1.50、1.75。如圖4(a)所示,隨著抽樣次數的增加,不同分布參數的失效概率預估值均趨于穩定。在α=β=1.50時波動迅速減小并穩定在收斂值上,此時樣本能夠囊括絕大多數失效區域,且中心區域仍占有較多的樣本,使得失效概率評估收斂較快。預估失效概率標準差如圖4(b)所示,α=β=1.50情況下預估失效概率標準差在不同抽樣次數下均處于較低水平。

圖4 不同分布參數下貝塔分布抽樣預測結果

2) 縮放參數。設置分布參數α=β=1.50,縮放參數kB分別取值6、8、10、12、14,不同抽樣次數下的單輪預測失效概率如圖5(a)所示。在kB取值為6時,由于模擬樣本分布范圍為(μ-3σ,μ+3σ),使得超出區域并沒有樣本,導致預估失效概率值低于準確值。如圖5(b)所示,由于縮放參數增加導致其在概率密度較大的失效區域分布的樣本量反而減少,使得預測失效概率標準差隨縮放參數增加而增大,在取值為10時可以達到較好的收斂效果。

圖5 不同縮放參數下貝塔分布抽樣預測結果

2.1.3 雙峰正態分布抽樣

影響雙峰正態分布抽樣密度函數的因素為平移參數kD。如圖6(a)所示,不同平移參數單輪抽樣結果均在106時趨于穩定,圖6(b)所示的標準差表明kD在取值為2.0時可以得到較好的收斂速度,此時2個正態分布的均值分別為μ-2σ和μ+2σ。預估失效概率的標準差并非隨著平移參數kD的增加而單調遞增或單調遞減,這是因為當平移參數kD較小時其抽樣落入失效區域的樣本偏少,其預測失效概率標準差與簡單隨機抽樣相似;而當平移參數kD較大時其抽樣樣本落入失效區域的概率增加,但是落入概率密度較大的失效區域的樣本量減少,導致其標準差也會增加。

圖6 不同平移參數下雙峰正態分布抽樣預測結果

2.1.4失效概率評估性能分析

設置貝塔分布參數α=β=1.50、kB=10,雙峰正態分布參數kD=2.0,進一步對比三種抽樣方法對失效概率預測性能的影響。1 000次隨機抽樣下不同抽樣方式樣本點的分布如圖7所示,簡單隨機抽樣只有極少數樣本分布于非失效區域(黑色曲線左側),貝塔分布抽樣和雙峰正態分布抽樣樣本分布較為均勻,從而能夠提高失效概率預測精度。

圖7 不同抽樣方式樣本分布

如圖8(a)所示,在抽樣次數為10時,由于雙峰正態分布和貝塔分布抽樣樣本點數目在失效和非失效區域分布量級一致,兩者給出了與理論值量級一致的預測失效概率。另一方面,圖8(b)為預測失效概率標準差相對于簡單隨機抽樣的值,可以看出雙峰正態分布失效概率預測標準差大約為簡單隨機抽樣的一半甚至更少,證明其在收斂性能上要優于另外2種抽樣方式。

圖8 特定參數下不同抽樣方法預測性能對比

按照圖2流程,設置收斂指標δ為5%,進行10輪蒙特卡洛試驗,不同抽樣方式達到收斂所需的抽樣次數如表1所示。取簡單隨機抽樣108次預測失效概率期望值0.581 2%作為理論值,雙峰正態分布抽樣最為接近理論計算結果,其平均抽樣次數不到簡單隨機抽樣的1/10;而且簡單隨機抽樣標準差最大,雙峰正態分布抽樣最小。

表1 不同抽樣方法算例1預測結果對比

2.2 算例2

九盒段結構由64個桿元件和42個板元件構成,材料為鋁合金,其2個主要失效模式的功能函數為[8]:

(14)

式中:R68、R77、R68為單元強度,外載荷P均值為120 kg,所有參數均滿足正態分布,如表2所示。

表2 算例2計算參數

類似地,設置收斂指標δ為5%,進行10輪蒙特卡洛試驗,其中貝塔分布參數α、β設置為1.5,縮放參數kB=10,雙峰正態分布平移參數kD=2.0,3種抽樣方法計算得到的結果如表3所示,雙峰正態分布抽樣達到收斂的平均次數為7.19×105,小于簡單隨機抽樣和貝塔分布抽樣,驗證了改進蒙特卡洛抽樣的準確性和高效性。

表3 不同抽樣方法算例2預測結果對比

3 微動疲勞壽命不確定性分析

3.1 有限元模型

參考某型號風扇葉片、輪盤提取榫連接結構,簡化為平面應力模型,如圖9所示,榫頭頂端施加風扇葉片離心載荷產生的壓強P,榫槽兩側施加周向約束,榫槽底部施加全約束,設置接觸區域摩擦系數為0.3。榫頭、榫槽均采用TC4鈦合金材料,如表4所示,參數均值來源于《中國航空材料手冊》[11]TC4棒材室溫下的低周疲勞數據,考慮隨機因素為壓強P和應變疲勞參數σ′f,假設兩者均服從正態分布。

表4 榫連接結構計算參數

圖9 榫連接結構有限元模型

取表4參數均值,計算結果如圖10所示,接觸區域應力分布與渦輪葉片榫連接結構[12]、樅樹形榫連接結構[13]有限元計算結果相近,最大Mises應力均位于接觸區邊緣。

圖10 榫連接結構應力分布

3.2 微動疲勞壽命計算

SWT損傷參量在鈦合金燕尾榫連接結構微動疲勞預測方面具有較好的效果,該參量綜合考慮了臨界平面上的循環正應變幅εn,a和最大正應力σn,max為

SWT=σn,maxεn,a

(15)

其中臨界平面為SWT值最大的平面,對應微動疲勞壽命為:

(16)

提取榫連接結構接觸區域的應力和應變信息,計算各個節點的SWT損傷參量,如圖11所示,在滑動區域兩側邊緣存在2個SWT極值,與榫連接結構試驗件計算結果一致[12]。將最大SWT值代入式(16)得到微動疲勞壽命為91 720,滿足10 000周循環的壽命要求。

圖11 接觸區域損傷參量

3.3 微動疲勞壽命失效概率

考慮拉伸壓強P和應變疲勞參數σ′f兩個隨機因素(表4),基于雙峰正態分布抽樣進行蒙特卡洛法模擬,計算榫頭微動疲勞壽命失效概率。設定功能函數為

g(X)=Nf(X)-104

(17)

雙峰正態分布抽樣結果如圖12所示,拉伸壓強P的增加和材料應變疲勞參數σ′f的減少,都將導致榫連接結構微動疲勞壽命下降,甚至失效。對于同一個失效的樣本點,相比于增加材料應變疲勞參數σ′f,改變更少的拉伸壓強P就能滿足壽命要求,如圖12兩個箭頭所示。因此,降低轉速波動更能有效降低結構失效概率。將抽樣結果代入式(5)計算得到失效概率為0.046 4%。

3.4 失效概率影響因素分析

考慮拉伸壓強P和應變疲勞參數σ′f的不確定性對榫連接結構失效概率的影響時,采用Cui等[14]提出的失效概率全局靈敏度指標,其表達式為

(18)

根據式(18)計算得到榫連接結構微動疲勞壽命失效概率影響因素的全局靈敏度指標如表5所示。拉伸壓強P的全局靈敏度指標最大,說明其對榫連接結構微動疲勞壽命失效概率影響最為顯著。這是因為拉伸壓強P增大時將同時增加接觸面的摩擦力和正壓力,使得循環正應變幅εn,a和最大正應力σn,max增大,導致微動疲勞損傷參量SWT值迅速增加,且在同等標準差下拉伸壓強P相對增加值更大,使得榫連接結構微動疲勞壽命迅速下降。因此,拉伸壓強P,即轉速波動對榫連接結構微動疲勞壽命影響不可忽略。

4 結論

針對蒙特卡洛法簡單隨機抽樣在評估失效概率效率不足這一問題,本文中提出了基于貝塔抽樣和雙峰正態分布抽樣的蒙特卡洛模擬方法,在不引入大量計算復雜度的前提下提高了不確定性評估效率。得到如下結論:

1) 算例分析結果表明:影響貝塔分布抽樣預測失效概率性能的因素包括分布參數α、β和縮放參數kB,縮放參數過小時將會導致失效區域分布樣本不足。

2) 算例分析結果表明:影響雙峰正態分布抽樣預測失效概率性能的因素為平移參數kD,平移參數過小時其失效概率預測性能將會趨于簡單隨機抽樣,過大將會導致落在概率密度較大的失效區域的樣本量減少。給定算例下雙峰正態分布抽樣在預測失效概率精度和收斂性方面均表現更好。

3) 基于雙峰正態分布抽樣完成了風扇榫連接結構微動疲勞壽命不確定性分析,計算得到失效概率為0.046 4%。基于全局靈敏度指標分析表明拉伸壓強P對榫連接結構微動疲勞壽命失效概率影響最為顯著,在設計中應給予考慮。