基于響應面函數的攝動隨機有限元方法研究

葛仁超

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

基于響應面函數的攝動隨機有限元方法研究

葛仁超

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

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

響應面函數；攝動隨機有限元；敏感度

0 引言

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

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

1 攝動法

1.1 單變量下的攝動法

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

假設某結構的單元厚度h和隨機彈性模量e，利用經典攝動理論思想，引入小參數ε，對狀態函數進行泰勒級數展開并截斷。則其n階展開式可寫為：

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

舉例，當m=2時，輸入隨機變量為正態分布時，同時忽略高于6階的變量偏導數，經計算簡化，得：

1.2 多變量下的攝動法

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

2 響應表面方法

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

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

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

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

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

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

表1 中心指數設計抽樣中需要的樣本點數目與隨機輸入變量數目的關系

確定該近似函數后，利用上文所推導出的公式求出結構隨機響應變量的均值和各階中心矩，并計算響應變量相對于輸入變量的靈敏度。

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

Sxi(i=1， 2，…n)為靈敏系數，其反映了各隨機輸入變量變化對結構隨機響應變量的影響相對強弱。

3 實例分析

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

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

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

表2 葉片隨機變量參數

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

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

表3 葉片最大變形、最大應力的統計參數

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

4 結語

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

通過對燃氣輪機葉片結構隨機響應的計算，驗證了該方法的計算結果。其精度較高，計算速度較快，并且適用于結構相對復雜的模型。

[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隨機有限元法的復合材料壓力容器的應力分析[J]. 西北工業大學學報,2002, 20(3): 363-367.

[4]安偉光, 朱衛兵. 隨機有限元法在不確定性分析中的應用[J]. 哈爾濱工程大學學報, 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]. 西安電子科技大學學報, 2007, 34(3): 453-457.

[7]王東, 陳建康. Monte-Carlo隨機有限元結構可靠度分析新方法[J]. 四川大學學報, 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]段巍, 王璋奇. 基于響應面方法的汽輪機葉片概率強度設計及敏感性分析[J]. 中國電機工程學報, 2007,27(5): 99-104.

[11]田四朋. 固體火箭發動機藥柱三維粘彈性響應面隨機有限元分析[J]. 固體火箭技術, 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-），男，工程師。主要從事蒸汽、燃氣動力裝置建造和質量監督工作。