基于响应面函数的摄动随机有限元方法研究

葛仁超

● （海军装备部装沈阳局，辽宁沈阳 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-），男，工程师。主要从事蒸汽、燃气动力装置建造和质量监督工作。