鞏博瑞，陳虹旭，殷 鳴，殷國富*

(1.西北機電工程研究所，陜西 咸陽 712099；2.四川大學(xué)制造科學(xué)與工程學(xué)院， 四川 成都 610065)

乏組件水下檢測裝置是核電站中用來對乏組件進(jìn)行測量的重要輔助設(shè)備，其穩(wěn)定有效工作直接影響到核電站的安全運行。由于材料特性在統(tǒng)計上的離散性以及測量、加工、制造誤差的存在，在地震工況下，裝置薄弱部位會產(chǎn)生較大應(yīng)力，且應(yīng)力可能在較大范圍內(nèi)波動，因此有必要對其進(jìn)行可靠性分析。對裝置整機初步分析發(fā)現(xiàn)，裝置中一對齒輪副屬于薄弱部件，本文主要針對此齒輪副進(jìn)行分析。

對地震工況下結(jié)構(gòu)的響應(yīng)研究，一般通過隨機振動，采用概率與統(tǒng)計方法研究結(jié)構(gòu)系統(tǒng)的動力響應(yīng)，從能量角度分析激勵的隨機性問題[1]。目前隨機振動分析[2-4]一般是基于有限元的，如：劉士華等[2]通過有限元軟件分析了鋼架結(jié)構(gòu)在地震位移激勵下的隨機振動響應(yīng)；馬乾瑛等[4]對地震工況下的地鐵車站結(jié)構(gòu)進(jìn)行了隨機振動響應(yīng)分析。進(jìn)行可靠性分析的工程問題通常需要多次調(diào)用建立的有限元模型，而有限元模型一般為隱式函數(shù)，無法用確定的解析式表示，計算量較大，計算時間成本高，效率低。相較于蒙特卡羅法[5]、響應(yīng)面法[6]，Kriging模型[7-8]用較少的樣本點就可以擬合出符合實際工程問題的代理模型，對于給定輸入?yún)?shù)的預(yù)測也較準(zhǔn)確，效率高。陳士誠等[7]基于Kriging模型對壓力容器開孔接管區(qū)結(jié)構(gòu)進(jìn)行了可靠性分析，驗證了Kriging模型在計算效率和精度上的優(yōu)勢。為實現(xiàn)對系統(tǒng)隨機特性和風(fēng)險水平的整體把握和完整認(rèn)知，得到結(jié)構(gòu)在各輸入?yún)?shù)下的響應(yīng)特性后，需要對結(jié)構(gòu)輸出響應(yīng)進(jìn)行概率密度函數(shù)估計。核密度估計是一種從數(shù)據(jù)樣本本身出發(fā)研究數(shù)據(jù)分布特征的方法，采用平滑的峰值函數(shù)來擬合數(shù)據(jù)點，不需要數(shù)據(jù)分布形式，不附加任何假設(shè)，在統(tǒng)計領(lǐng)域有較大應(yīng)用[9]。趙淵等[10]在大電網(wǎng)可靠性樣本的基礎(chǔ)上，通過非參數(shù)核密度估計實現(xiàn)了可靠性指標(biāo)的概率密度估計，解決了傳統(tǒng)期望值指標(biāo)僅從概率平均意義角度測量系統(tǒng)的可靠性問題。

針對本文研究對象，筆者首先采用有限元軟件ANSYS對齒輪副結(jié)構(gòu)進(jìn)行隨機振動分析，得到其在地震工況下的響應(yīng)，然后對隨機振動分析中的關(guān)鍵參數(shù)通過多次有限元迭代，建立Kriging模型，最后結(jié)合不確定性變量與Kriging模型進(jìn)行不確定性分析，計算結(jié)構(gòu)失效概率并通過核密度估計得到結(jié)構(gòu)隨機振動響應(yīng)概率密度曲線。

1 隨機振動理論

當(dāng)系統(tǒng)受到不能用時間的確定函數(shù)描述的激勵作用時，產(chǎn)生的不確定性的振動過程稱為隨機振動，其分析一般是通過輸入功率譜密度函數(shù)進(jìn)行的。功率譜密度是隨機動態(tài)載荷激勵下系統(tǒng)響應(yīng)的統(tǒng)計結(jié)果，是功率譜密度值與頻率值的關(guān)系曲線。功率譜密度可以是位移功率譜密度、速度功率譜密度、加速度功率譜密度、力功率譜密度等[11-12]，其定義為

(1)

式中：Sx(x)為激勵信號的功率譜密度函數(shù)；Rxx(τ)為激勵信號自相關(guān)函數(shù)；信號在時間間隔為τ的兩數(shù)值之間的相互關(guān)系，用Rx(τ)表示。

(2)

當(dāng)系統(tǒng)受到一平穩(wěn)隨機激勵時，響應(yīng)的功率譜密度可以表示為

(3)

式中H(ω)為傳遞函數(shù)。模態(tài)分析是隨機振動分析的基礎(chǔ)，通過模態(tài)分析，將機構(gòu)固有頻率、振型結(jié)果等作為隨機振動分析輸入，再通過隨機振動分析，可得到分布在正態(tài)區(qū)間上的應(yīng)力、應(yīng)變值，從而為設(shè)計提供參考。

2 Kriging模型

Kriging模型是線性回歸分析的一種改進(jìn)技術(shù)，包含了線性回歸部分和非參數(shù)部分[8]，隨機過程的實現(xiàn)取決于非參數(shù)部分，對于預(yù)近似的關(guān)于x多項式函數(shù)y(x)可以表示為

y(x)=F(β，x)+Z(x)=fTβ+Z

(4)

式中：f為x的多項式函數(shù)；β為回歸系數(shù)；Z(x)為一均值為0、方差為σ2、協(xié)方差非零的統(tǒng)計過程，協(xié)方差可表示為

Cov[Z(xi)，Z(xj)]=σ2R[(xi，xj)]

(5)

R(xi，xj)為樣本點xi和xj的相關(guān)函數(shù)，通常采用Gaussian相關(guān)方程，形式為

(6)

式中θk和Pk為待定參數(shù)。對于任意點x，Kriging模型對y的預(yù)測為

(7)

其中r為待測點x和樣本點間的相關(guān)向量，

r(x)=[R(x，x1)，R(x，x2)，…，R(x，xn)]

(8)

β*為極大似然估計因子，

β*=(FTR-1F)-1FTR-1Y

(9)

由此，可得Kriging模型對梯度dy/dx的預(yù)測值為

(10)

式中Jf(x)和Jr(x)分別為f和r的Jacobian矩陣。

3 核密度估計

分布密度能以曲線的形式直觀表示可靠性指標(biāo)圍繞其平均值的變化趨勢和變化程度，尾部特征可給出系統(tǒng)在嚴(yán)重風(fēng)險下的信息，在實踐中有重要應(yīng)用價值。概率統(tǒng)計學(xué)基本問題之一是由給定樣本點集合求解隨機變量的分布密度，通常有參數(shù)法和非參數(shù)法2種。參數(shù)法需要提前確定樣本的分布，由于分布種類較多，各分布特性相對單一，理論及實際應(yīng)用表明采用該方法參數(shù)模型與實際物理模型差異較大，不能獲得理想的分布函數(shù)[14]。核密度估計(kernel density estimation，KDE)是一種估計概率密度函數(shù)的非參數(shù)方法，由Rosenblatt和Parzen提出，其不利用數(shù)據(jù)分布的先驗知識及概率分布形式的假設(shè)，對模型概率密度有很好的估計[10]。

設(shè)X1，X2，…，Xn，是滿足概率密度為f(x)的總體X的樣本，則其核概率密度估計定義為

(11)

式中：hn為窗寬；K(*)為核概率密度函數(shù)，通常選用以0為中心的對稱單峰概率密度函數(shù)。為保證合理性，K(*)需滿足：

(12)

(13)

其中σ為樣本標(biāo)準(zhǔn)差。

4 齒輪副結(jié)構(gòu)的有限元分析

有限元仿真分析可分為有限元模型建立、模態(tài)分析和隨機振動分析3個過程。 本文采用UG建立三維模型，利用ANSYS軟件Workbench模塊進(jìn)行有限元分析。

4.1 有限元模型建立

通過三維建模軟件UG對聯(lián)軸器、軸承、鍵等零部件進(jìn)行刪除、簡化及合并，齒輪利用GC工具箱通過參數(shù)化方法建立，得到的三維模型由一對嚙合的圓柱齒輪副及與齒輪配合的軸組成，倒角、細(xì)小孔、鍵槽特征被簡化刪除。

將簡化模型導(dǎo)入Workbench，進(jìn)行網(wǎng)格劃分。材料屬性如表1所示。

表1 材料屬性

根據(jù)實際工作狀況，將齒輪與軸約束設(shè)定為綁定接觸，兩齒輪之間約束設(shè)定為摩擦約束，檢驗并修改軟件自動生成接觸。網(wǎng)格劃分對分析結(jié)果有較大影響，本次分析在軟件自動劃分的基礎(chǔ)上，通過添加體尺寸策略、接觸尺寸策略等對網(wǎng)格進(jìn)行優(yōu)化，網(wǎng)格劃分結(jié)果如圖1所示，共有6萬6 622個節(jié)點，3萬6 846個單元。

圖1 齒輪副有限元模型

4.2 模態(tài)分析

模態(tài)分析是結(jié)構(gòu)動力學(xué)分析的基礎(chǔ)，通過其可得到結(jié)構(gòu)的各階固有頻率及振型。通過功率譜密度函數(shù)可得到各模態(tài)下的模態(tài)響應(yīng)幅值，各模態(tài)響應(yīng)幅值疊加的結(jié)果就是隨機振動分析的結(jié)構(gòu)響應(yīng)。根據(jù)齒輪副實際工況，在齒輪軸兩端添加固定約束，分析得到結(jié)構(gòu)各階固有頻率如表2所示，部分振型如圖2所示。

表2 模態(tài)分析結(jié)構(gòu)各階固有頻率

4.3 隨機振動分析

地震激勵的加速度功率譜密度(PSD)如表3所示，將其施加在2個齒輪軸兩端的固定約束處。選用模態(tài)為全部10階頻率，設(shè)定系統(tǒng)的阻尼比為2%。得到隨機振動分析應(yīng)力云圖，如圖3所示，其3σ應(yīng)力為233.38 MPa，在齒輪輪齒嚙合處最大，即在嚙合位置處應(yīng)力有99.73%的概率為233.38 MPa。

表3 加速度功率譜

(a)1階振型

(b)4階振型

(c)7階振型

(d)9階振型

(a)齒輪副整體應(yīng)力云圖 (b)齒輪局部應(yīng)力云圖

5 基于代理模型的不確定性分析

由于代理模型與實際模型有一定的誤差，傳統(tǒng)基于代理模型的響應(yīng)分析結(jié)果必然存在誤差。有限元分析顯示齒輪副嚙合過程中應(yīng)力較大，而且在實際工作過程中材料屬性、零件尺寸參數(shù)、零件安裝配合參數(shù)等變量都會產(chǎn)生波動，當(dāng)各變量在特定組合時，有可能導(dǎo)致齒輪副應(yīng)力過大，從而使齒輪結(jié)構(gòu)發(fā)生破壞。

本文將代理模型隨機振動下的3σ應(yīng)力和參數(shù)不確定性相結(jié)合，通過Kriging建立代理模型，以有限元分析計算各不確定狀態(tài)下響應(yīng)，進(jìn)而得到基于代理模型隨機振動分析的可靠度，最后通過核密度函數(shù)得到結(jié)構(gòu)輸出響應(yīng)的應(yīng)力最大值的概率分布函數(shù)。不確定性分析流程如圖4所示。

圖4 不確定性分析流程

5.1 Kriging模型

本文以隨機振動3σ應(yīng)力為目標(biāo)變量，以齒輪的軸半徑r、峰值功率譜a(即100～1 000 Hz段功率譜)及軸彈性模量E為不確定性變量。假設(shè)各不確定性變量滿足正態(tài)分布，為：

r～N(35，0.12)
a～N(0.02，0.22)
E～N(2×1011，0.12)

(14)

為滿足精度要求、模型收斂條件并減小仿真計算量，根據(jù)齒輪的軸半徑r、峰值功率譜a和軸材料的彈性模量E的分布參數(shù)，產(chǎn)生65組隨機樣本點(每組樣本點包含1個軸半徑r值、1個峰值功率譜a值和1個軸材料的彈性模量E的值)，將每組樣本點分別作為Workbench中有限元分析的模型參數(shù)，得到齒輪副結(jié)構(gòu)隨機振動的有限元計算3σ應(yīng)力結(jié)果值。將65組樣本點和對應(yīng)的應(yīng)力結(jié)果值作為Kriging模型的訓(xùn)練樣本點，從而建立關(guān)于齒輪的軸半徑、峰值功率譜密度、軸彈性模量與隨機振動下3σ應(yīng)力的Kriging模型。

5.2 失效概率計算

為計算結(jié)構(gòu)的可靠性，需要預(yù)測各變量在隨機波動下目標(biāo)變量3σ應(yīng)力。根據(jù)結(jié)構(gòu)參數(shù)變量的分布特性，采用蒙特卡洛法生成106組樣本點，代入建立的Kriging模型，預(yù)測各樣本點下隨機振動3σ應(yīng)力值。

3σ應(yīng)力可能失效功能函數(shù)g(x)為

g(x)=σs-σKriging

(15)

式中：σs為極限應(yīng)力值，根據(jù)強度要求取380 MPa；σKriging為基于Kriging模型預(yù)測的3σ應(yīng)力。齒輪副結(jié)構(gòu)失效概率P的計算公式為

P=P(g(x)<0)=P(σs-σKriging<0)

(16)

5.1節(jié)建立了齒輪的軸半徑、峰值功率譜密度、軸彈性模量與隨機振動下 3σ應(yīng)力的Kriging模型，基于此模型，可以采用傳統(tǒng)的蒙特卡洛方法來求解3σ應(yīng)力齒輪副結(jié)構(gòu)的失效概率。具體的求解流程如下：

1)根據(jù)齒輪的軸半徑r、峰值功率譜a和軸材料的彈性模量E的分布參數(shù)，隨機產(chǎn)生N個輸入樣本點A；

2)將A代入5.2節(jié)建立的Kriging模型中，可得輸出響應(yīng)B(即3σ應(yīng)力值，為N×1數(shù)組)；

3)根據(jù)式(15)可得失效功能函數(shù)值C(為N×1數(shù)組)；

4)統(tǒng)計數(shù)組C中小于零的項，總數(shù)記為M；

5)按式(17)計算失效概率。

(17)

計算得到齒輪副結(jié)構(gòu)失效概率為 0.000 573 9，齒輪副結(jié)構(gòu)在隨機振動下失效概率低，可靠。

5.3 3σ應(yīng)力概率密度估計

由Kriging模型預(yù)測得到106個隨機振動3σ應(yīng)力值，采用核密度估計，得到3σ應(yīng)力值概率密度曲線，如圖5所示?？梢钥闯?，3σ應(yīng)力在[200，370]MPa區(qū)間內(nèi)概率密度值較高，在區(qū)間兩側(cè)概率密度較低，分布情況與實際工程經(jīng)驗一致。

圖5 概率密度曲線

6 結(jié)束語

1)對水下檢測裝置中一對齒輪副結(jié)構(gòu)進(jìn)行了多軟件聯(lián)合隨機振動分析，仿真得到了該結(jié)構(gòu)在地震激勵功率譜下的響應(yīng)。

2)考慮到各不確定性變量，對齒輪副結(jié)構(gòu)進(jìn)行基于代理模型的不確定性分析。通過建立的Kriging模型預(yù)測各隨機狀態(tài)下的隨機振動3σ應(yīng)力。計算得到的失效概率顯示，該結(jié)構(gòu)可靠度較高，可滿足核電工作條件。