葛仁超

● （海軍裝備部裝沈陽局，遼寧沈陽 110031）

基于響應(yīng)面函數(shù)的攝動隨機有限元方法研究

葛仁超

● （海軍裝備部裝沈陽局，遼寧沈陽 110031）

針對隨機結(jié)構(gòu)問題，根據(jù)現(xiàn)有方法的優(yōu)缺點，將響應(yīng)面法和攝動隨機有限元相結(jié)合?？紤]結(jié)構(gòu)載荷參數(shù)、材料參數(shù)和幾何參數(shù)，結(jié)構(gòu)的響應(yīng)面近似函數(shù)采用步進回歸分析技術(shù)被確定了，其中抽樣方法為中心指數(shù)法，并利用攝動法推導(dǎo)出結(jié)構(gòu)隨機響應(yīng)的中心矩表達式。以燃氣渦輪葉片為例，計算出葉片的隨機響應(yīng)變量中心矩，及響應(yīng)變量相對于輸入變量的靈敏度，并分別與直接響應(yīng)面法、直接Monte-Carlo模擬方法的計算結(jié)果對比。結(jié)果表明，該方法適用于結(jié)構(gòu)相對復(fù)雜的模型，計算結(jié)果精度較高，計算速度較快。

響應(yīng)面函數(shù)；攝動隨機有限元；敏感度

0 引言

到目前為止，研究隨機結(jié)構(gòu)的方法主要有攝動隨機有限元法、Neumann隨機有限元法、Monte-Carlo隨機有限元法等。其中，攝動隨機有限元法、Neumann隨機有限元法為非統(tǒng)計方法[1-4]。當(dāng)研究對象為大型復(fù)雜結(jié)構(gòu)時，由于這類方法需要編程，此時，其相對應(yīng)的編程量會變得非常巨大。因此，該類方法不適用于大型復(fù)雜結(jié)構(gòu)。Monte-Carlo隨機有限元法為統(tǒng)計方法[5-8]，當(dāng)研究對象為大型復(fù)雜結(jié)構(gòu)時，由于該方法需要反復(fù)進行大量的模擬計算，因此并不適用于大型復(fù)雜結(jié)構(gòu)。

直接響應(yīng)面法[9-11]是目前較好的一種統(tǒng)計方法。該方法適用于大型的復(fù)雜結(jié)構(gòu)對象，并且不需要大量的編程。同時，與直接Monte-Carlo法相比，直接響應(yīng)面法所需要的抽樣次數(shù)也很少。但是，該方法需要進行兩輪的抽樣模擬計算，這樣就會引起前后兩次模擬計算的誤差迭加。為了改善該方法，本文對第2次的抽樣模擬計算改成利用數(shù)學(xué)分析的方法進行計算，即將響應(yīng)面法和攝動隨機有限元相結(jié)合，考慮結(jié)構(gòu)載荷參數(shù)、材料參數(shù)和幾何參數(shù)的隨機性，其中抽樣方法為中心指數(shù)法，結(jié)構(gòu)的響應(yīng)面近似函數(shù)采用步進回歸分析技術(shù)來確定，并利用攝動法推導(dǎo)出結(jié)構(gòu)隨機響應(yīng)的中心矩表達式。最后計算出結(jié)構(gòu)隨機響應(yīng)變量中心矩，及響應(yīng)變量相對于輸入變量的靈敏度。

1 攝動法

1.1 單變量下的攝動法

對于隨機變量x，概率密度函數(shù)為f (x)，則隨機變量x的m階中心矩為：

假設(shè)某結(jié)構(gòu)的單元厚度h和隨機彈性模量e，利用經(jīng)典攝動理論思想，引入小參數(shù)ε，對狀態(tài)函數(shù)進行泰勒級數(shù)展開并截斷。則其n階展開式可寫為：

為了結(jié)構(gòu)的狀態(tài)函數(shù)Y(x)的均值和m階中心矩，在此假設(shè)各隨機變量是連續(xù)的，并利用泰勒展開公式。其計算表達式為：

舉例，當(dāng)m=2時，輸入隨機變量為正態(tài)分布時，同時忽略高于6階的變量偏導(dǎo)數(shù)，經(jīng)計算簡化，得：

1.2 多變量下的攝動法

當(dāng)結(jié)構(gòu)的狀態(tài)函數(shù)為Y(x1，x2，…xn)時，表示該結(jié)構(gòu)存在n個輸入隨機變量，利用泰勒級數(shù)展開法，在均值點進行泰勒級數(shù)展開，令εi=1 (i=1， 2，…n)，同時忽略高于二階的變量偏導(dǎo)數(shù)。

2 響應(yīng)表面方法

用數(shù)學(xué)函數(shù)來表達隨機輸入變量對于隨機輸出變量的影響。如果該函數(shù)是一個二次多項式，那么，擬合函數(shù)Y可以表達成：

式中：c0為常數(shù)項；ci(i=1， 2，…N)為線性項系數(shù)；cij(j=1， 2，…N)為二次項系數(shù)。

采用中心指數(shù)設(shè)計進行抽樣，其抽樣包括一個中心點、2N個軸線點和位于2N-f階乘個N維超立方體的頂點。其中，f為中心指數(shù)設(shè)計階乘因子表達式中的一個參數(shù)，N為隨機輸入變量的數(shù)目。

響應(yīng)表面方法[10]分為兩個步驟：1）計算對應(yīng)于隨機輸入變量空間樣本點的隨機響應(yīng)變量的數(shù)據(jù)：采用有限元方法進行循環(huán)計算；2）確定近似函數(shù)：采用步進回歸分析技術(shù)。

為了始終確保仿真循環(huán)數(shù)目總是合理的，當(dāng)隨機輸入變量數(shù)目逐步增加時，可逐步增大階乘因子參數(shù) f。圖 1是1個有3個隨機輸入變量的樣本點位置示意圖。表1是中心指數(shù)設(shè)計抽樣中需要的樣本點數(shù)目(即仿真循環(huán)次數(shù))與隨機輸入變量數(shù)目的關(guān)系。

圖1 有3個隨機輸入變量的樣本點位置示意圖

表1 中心指數(shù)設(shè)計抽樣中需要的樣本點數(shù)目與隨機輸入變量數(shù)目的關(guān)系

確定該近似函數(shù)后，利用上文所推導(dǎo)出的公式求出結(jié)構(gòu)隨機響應(yīng)變量的均值和各階中心矩，并計算響應(yīng)變量相對于輸入變量的靈敏度。

其中隨機響應(yīng)變量相對于各隨機輸入變量的靈敏度計算公式為：

Sxi(i=1， 2，…n)為靈敏系數(shù)，其反映了各隨機輸入變量變化對結(jié)構(gòu)隨機響應(yīng)變量的影響相對強弱。

3 實例分析

本文以某燃氣渦輪第5級葉片為研究對象，葉片為斜齒長葉根變截面帶冠葉片。根據(jù)結(jié)構(gòu)圖紙，首先利用ANSYS軟件對燃氣渦輪葉片進行有限元參數(shù)化建模，并采用結(jié)構(gòu)分析單元solid185號單元對葉片進行自由劃分網(wǎng)格。圖2為燃氣渦輪葉片的幾何模型。

圖2 燃氣渦輪葉片的幾何模型

利用本文提出的方法，把燃氣渦輪葉片角速度 ω作為隨機載荷參數(shù)，葉片密度ρ、彈性模量E作為隨機材料參數(shù)，分布類型均為正態(tài)；其他各參數(shù)為確定性參數(shù)。為了節(jié)約篇幅，只將燃氣渦輪葉片的隨機變量參數(shù)列出，如表2所示。

表2 葉片隨機變量參數(shù)

將E、ρ、ω等3個變量作為隨機輸入變量，葉片的δmax和σmax作為隨機響應(yīng)變量。利用中心指數(shù)設(shè)計抽樣，4個隨機輸入變量需要25個樣本點，即調(diào)用確定性有限元25次。然后，通過采用步進回歸技術(shù)過濾掉不重要的項，利用帶有交叉項的二次近似函數(shù)回歸模型擬合樣本，計算得出最大變形和最大應(yīng)力的近似函數(shù)。

根據(jù)式(9)、(11)和(15)計算得出最大變形和最大應(yīng)力的均值、方差和各變量的靈敏度。同時，采用Monte-Carlo方法計算葉片的均值和標(biāo)準(zhǔn)方差，樣本點為1000個，對照結(jié)果如表3所示。

表3 葉片最大變形、最大應(yīng)力的統(tǒng)計參數(shù)

從表3可知，與直接響應(yīng)面法、Monte-Carlo方法相比，驗證了本文方法的正確性，并且該方法費時較少。同時，通過計算，我們可以得知：1）對葉片的最大應(yīng)力σmax影響較大的有葉片密度ρ、轉(zhuǎn)速ω，其中葉片密度對其影響最大，影響比例為61.03%，而彈性模量E對其沒有影響；2）對葉片的最大變形σmax影響較大的有彈性模量E、葉片密度ρ和轉(zhuǎn)速ω，其中彈性模量E對其影響最大，影響比例為38.3%。

4 結(jié)語

該方法的主要思想是將響應(yīng)面法和攝動隨機有限元相結(jié)合，考慮結(jié)構(gòu)載荷參數(shù)、材料參數(shù)和幾何參數(shù)的隨機性，其中抽樣方法為中心指數(shù)法，結(jié)構(gòu)的響應(yīng)面近似函數(shù)采用步進回歸分析來確定，并利用攝動法推導(dǎo)出結(jié)構(gòu)隨機響應(yīng)的中心矩表達式。計算出隨機結(jié)構(gòu)的隨機響應(yīng)變量中心矩及響應(yīng)變量相對于輸入變量的靈敏度。

通過對燃氣輪機葉片結(jié)構(gòu)隨機響應(yīng)的計算，驗證了該方法的計算結(jié)果。其精度較高，計算速度較快，并且適用于結(jié)構(gòu)相對復(fù)雜的模型。

[1]Marcin Kaminski. Generalized perturbation-based stochastic finite element method in elastostatics[J].Computers and Structures, 2007, 85: 586-594.

[2]Xue Xiaofeng, Feng Yunwen. Research on the plane multiple cracks stress intensity factors based on stochastic finite element method [J]. Chinese journal of Aeronautics, 2009, 22: 257-261.

[3]雷開貴，鄧子辰. 基于Neumann隨機有限元法的復(fù)合材料壓力容器的應(yīng)力分析[J]. 西北工業(yè)大學(xué)學(xué)報,2002, 20(3): 363-367.

[4]安偉光, 朱衛(wèi)兵. 隨機有限元法在不確定性分析中的應(yīng)用[J]. 哈爾濱工程大學(xué)學(xué)報, 2002, 23(1): 132-135.

[5]Yang Z J, Su X T. Monte Carlo simulation of complex cohesive fracture in random heterogeneous quasi-brittle materials [J]. International Journal of Solids and Structures, 2009, 46: 3222-3234.

[6]李金平, 陳建軍. 基于Neumann展開Monte-Carlo有限元法的隨機溫度場分析[J]. 西安電子科技大學(xué)學(xué)報, 2007, 34(3): 453-457.

[7]王東, 陳建康. Monte-Carlo隨機有限元結(jié)構(gòu)可靠度分析新方法[J]. 四川大學(xué)學(xué)報, 2008, 40(3): 20-26.

[8]Marcin Kaminski. Perturbation-based stochastic finite element method using polynomial response function for the elastic beams [J]. Mechanics Research Communications, 2009, 36: 381-390.

[9]Afeefa Shaker, Wael Abdelrahman. Stochastic finite element analysis of the free vibration of functionally graded material plates [J]. Comput Mehthods Appl Mech Engrg, 2008, 41: 707-714.

[10]段巍, 王璋奇. 基于響應(yīng)面方法的汽輪機葉片概率強度設(shè)計及敏感性分析[J]. 中國電機工程學(xué)報, 2007,27(5): 99-104.

[11]田四朋. 固體火箭發(fā)動機藥柱三維粘彈性響應(yīng)面隨機有限元分析[J]. 固體火箭技術(shù), 2010, 33(1): 17-20.

Study on Perturbation Stochastic Finite Element Method Based on Response Surface Function

GE Ren-chao
(Naval Armaments Department, Bureau of Shenyang Military Representative, Shenyang 110031, China)

Aiming at problem of the stochastic structure, according to the merits and demerits of the existing methods, the response surface method and perturbation-based stochastic finite element method are combined. Considering structure load parameters,material parameters and geometrical parameters, the approximate function of response surface structure is determined by using step regression analysis techniques, the sample of which is taken by the center index design method. And the central moment expression of the structure random response is derived by using perturbation method. Taking the gas turbine blade as an example, the central moment of the random response variables and the sensitivity of response variable corresponding to the input variables for the blade are calculated. And it is contrasted with the calculation results of direct response surface method and direct Monte-Carlo simulation method respectively. The results show that the method is suitable for a relatively complex structure model. The precision of computing result is higher, and the calculation speed is faster.

response surface function; perturbation stochastic finite element; sensitivity

O242

A

葛仁超（1981-），男，工程師。主要從事蒸汽、燃氣動力裝置建造和質(zhì)量監(jiān)督工作。