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

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

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

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

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

1 隨機(jī)振動(dòng)理論

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

(1)

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

(2)

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

(3)

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

2 Kriging模型

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

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

(4)

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

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

(5)

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

(6)

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

(7)

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

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

(8)

β*為極大似然估計(jì)因子，

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

(9)

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

(10)

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

3 核密度估計(jì)

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

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

(11)

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

(12)

(13)

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

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

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

4.1 有限元模型建立

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

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

表1 材料屬性

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

圖1 齒輪副有限元模型

4.2 模態(tài)分析

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

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

4.3 隨機(jī)振動(dòng)分析

地震激勵(lì)的加速度功率譜密度(PSD)如表3所示，將其施加在2個(gè)齒輪軸兩端的固定約束處。選用模態(tài)為全部10階頻率，設(shè)定系統(tǒng)的阻尼比為2%。得到隨機(jī)振動(dòng)分析應(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 基于代理模型的不確定性分析

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

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

圖4 不確定性分析流程

5.1 Kriging模型

本文以隨機(jī)振動(dòng)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ì)算量，根據(jù)齒輪的軸半徑r、峰值功率譜a和軸材料的彈性模量E的分布參數(shù)，產(chǎn)生65組隨機(jī)樣本點(diǎn)(每組樣本點(diǎn)包含1個(gè)軸半徑r值、1個(gè)峰值功率譜a值和1個(gè)軸材料的彈性模量E的值)，將每組樣本點(diǎn)分別作為Workbench中有限元分析的模型參數(shù)，得到齒輪副結(jié)構(gòu)隨機(jī)振動(dòng)的有限元計(jì)算3σ應(yīng)力結(jié)果值。將65組樣本點(diǎn)和對(duì)應(yīng)的應(yīng)力結(jié)果值作為Kriging模型的訓(xùn)練樣本點(diǎn)，從而建立關(guān)于齒輪的軸半徑、峰值功率譜密度、軸彈性模量與隨機(jī)振動(dòng)下3σ應(yīng)力的Kriging模型。

5.2 失效概率計(jì)算

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

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

g(x)=σs-σKriging

(15)

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

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

(16)

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

1)根據(jù)齒輪的軸半徑r、峰值功率譜a和軸材料的彈性模量E的分布參數(shù)，隨機(jī)產(chǎn)生N個(gè)輸入樣本點(diǎ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)計(jì)數(shù)組C中小于零的項(xiàng)，總數(shù)記為M；

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

(17)

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

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

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

圖5 概率密度曲線

6 結(jié)束語

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

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