基于虚位移原理的随机有限元方法研究

周宗和，尹喜庆

（海军驻武汉四三八厂军事代表室，武汉430060）

目前，研究不确定性结构的常用方法可分为两类。一类是统计的方法［1－9］，就是通过大量的随机抽样，对结构反复进行有限元计算，最后对得到的结果作统计分析，这类方法称为Monte－Carlo随机有限元法，需要反复进行大量的模拟计算，并不适用于大型复杂结构。另一类方法［10－14］是以数学、力学分析作为工具，得出结构系统的响应与输入信号之间的关系，并据此得到结构应力或位移的统计规律，这类方法有摄动随机有限元法、Neumann随机有限元法等，这类方法在大型复杂结构中的应用也非常有限。

胡海昌［15］从非统计方法出发，通过对偏微分方程弱解的积分形式——虚位移原理进行较深入的研究，王勖成、邵敏［16］和刘文廷［17］构造了八节点四边形轴对称等参单元。本文在此基础上给出轴对称静态问题的平衡方程以及力的边界条件下的等效积分“弱”形式，结合摄动技术推导出八节点四边形轴对称等参单元的随机静态有限元方程，即基于虚位移原理的随机有限元方法，最后考虑汽轮机转子的物性参数和载荷参数的随机性，利用这种方法对汽轮机转子的静态随机响应特性进行计算和分析，并将计算结果与通过直接Monte－Carlo模拟仿真法计算得出的结果相对比，以验证该方法的正确性。

1 八节点四边形轴对称等参单元

八节点四边形轴对称等参单元如图1所示。

该单元为横截面为八节点四边形的360°环形等参单元，其横截面上八个节点的位置坐标（ri，zi），i＝1，2，…，8，各个节点的位移为（ui，wi），i＝1，2，…，8。

该单元为绕z轴的环状单元，在orz平面内，八节点四边形轴对称等参单元共有16个自由度。将所有节点上的位移组成一个列阵，记作｛｝ae，同样，将所有节点上的各个力也组成一个列阵，记作｛｝pe，那么

式（1）中，采用的等参变换位移模式和应变矩阵计算式如下：

1）位移模式。

式（3）中：［N］为形函数矩阵。

写成矩阵的形式为：

2）应变矩阵。

由轴对称问题的几何方程可以推出相应的几何函数矩阵，得

式（5）中：［B］为几何矩阵。

2 线弹性小变形静态轴对称等参问题的“弱”形式

对于线弹性小变形静态轴对称问题，运用虚位移原理求位移，由平衡方程［16］得

式（6）中：δu，δw为两个方向的虚位移。

由分部积分公式将其化为“弱”形式，利用几何方程可得

式（7）中：l为单元体的边界。

式（7）表示为向量形式如下：

将本构方程的矩阵表示代入上式可得

将力和位移边界条件代入上式（对于位移边界条件，虚位移δu，δw为0，可得

由此可得出轴对称等参问题的“弱”形式如下：

式（11）中：［J］为雅可比矩阵，其中

式（12）中：S为坐标变换后求解域的面元，L为坐标变换后的单元体的边界。

将式（4）和（5）代入式（11），整理后得出轴对称等参问题的“弱”形式，具体如下：

由此可推出零阶随机积分“弱”形式和一阶随机积分“弱”形式。

3 轴对称等参问题的静力学随机有限元方程

在许多实际问题当中，影响结构静力学特性的因素，如材料特性、几何特性和载荷特性等都有一定程度随机性，因此结构的响应｛u｝和｛ε｝也具有随机性。本文为了简化书写，假设任意一个给定的空间连续函数f（｛r0｝），｛r0｝为空间的坐标向量，基本随机向量｛b｝的联合概率密度函数为g（｛b｝），εp为一个小参数，根据向量的数学期望计算公式有下：

将式（14）代入式（13），并比较小参数εp的同次幂系数，可以得到结构摄动格式的各阶随机积分“弱”形式，其中零阶和一阶随机积分“弱”形式如下：

1）零阶随机积分“弱”形式。其方程如下：

2）一阶随机积分“弱”形式。其方程如下：

基于摄动理论［18］将影响结构的随机参数在基本随机变量｛b｝的均值｛b0｝处作小参数摄动展开，如

式（15）、（16）中，δ｛a｝eT是δ｛a｝e的转置向量。

由于δ｛a｝e是任意的，故由式（16）可得到单元的零阶静态随机有限元方程：

式（17）中：［K0］e为单元的线弹性刚度矩阵均值，｛F0｝e为单元的体积力均值，｛T0F｝e为单元的面力向量均值，其相应的单元矩阵如下：

求解区域离散为n个单元，经组集后可以得到结构的零阶静态随机有限元方程：

式（19）中：［K0］为结构的线弹性刚度矩阵均值，｛F0｝为结构的体积力均值，｛T0F｝为结构的边界面力向量均值。

同理，由式（17）可得到单元的一阶静态随机有限元方程：

式（20）中：［K］e′j为单元的线弹性刚度矩阵的变异矩阵，｛F｝e′j为单元的体积力向量的变异向量，｛TF｝e′j为单元的边界面力向量的变异向量，分别如下表示：

求解区域离散为n个单元，经组集后可以得到结构的一阶静态随机有限元方程：

4 汽轮机转子静态随机有限元分析

由于汽轮机转子的材料参数、几何参数和载荷参数都存在随机性，所以转子的变形和应力也存在随机性。因此，本文采用基于虚位移原理的随机有限元方法和直接 Monte－Carlo模拟仿真法［19］2种方法分别对汽轮机转子的静态随机响应特性进行研究，计算汽轮机转子变形和应力的均值和方差，对这2种方法的计算结果进行比较，以验证基于虚位移原理的随机有限元方法的正确性。

4．1 汽轮机转子的随机有限元方程

对推导出的轴对称等参问题的静态随机有限元方程求解，可以得到位移的统计特性。

求解式（19），可以得到位移｛a｝的均值

式（25）中：［K0］为汽轮机转子的线弹性刚度矩阵均值，｛F0｝为汽轮机转子离心力向量均值，｛T0F｝为汽轮机转子的边界面力向量均值。

求解式（25），可以得到位移的一阶变异量

位移一阶变异量的矩阵形式为：

式（26）称为各节点位移相对于各随机变量的灵敏度矩阵，其中m为基本随机变量的个数，S｛a｝ij表示第i个位移分量对第j个基本随机变量的灵敏度。

忽略位移的二阶及二阶以上摄动展开式，进行协方差运算，可以得到位移的协方差矩阵

式（28）中：［C｛b｝］为基本随机变量的协方差矩阵。求出转子的节点位移后，则转子的单元应力为：

将单元应力在基本随机变量的均值点处进行二阶Taylor级数展开，得

比较小参数εp的同次幂系数，分别得出单元应力的均值和一阶变异量：

将位移的一阶变异量写成矩阵形式

式（33）称为单元应力相对于各随机变量的灵敏度矩阵。同样，可以得到应力的协方差矩阵

式（34）中：［C｛b｝］为基本随机变量的协方差矩阵；［S｛σ｝］ij表示第i个应力分量对第j个基本随机变量的灵敏度。

4．2 汽轮机转子的二维几何模型

由于汽轮机转子本身是一个复杂的部件，所以在不影响转子计算准确性的前提下，应尽量使计算趋于简单化，因此，将其简化为轴对称问题，同时，为了减少边界条件的设定误差对计算结果可能造成的影响，建模时取转子在汽缸内的整段剖面图作为研究对象。

转子的结构、形状、尺寸及材料取自某型船设计资料中的转子加工用图以及总装配图，利用网格划分软件Grid 9．0对汽轮机转子进行二维建模，如图2所示。

4．3 汽轮机转子的随机有限元程序流程

汽轮机转子的随机静态有限元分析程序流程图如图3所示。

4．4 汽轮机转子的随机静态响应特性分析

随机场的相关长度越长，表示随机场的空间两点之间的相关性越强，当随机场的相关长度等于无穷大时，随机场退化为一个随机变量。

本文把转子的弹性模量、密度和转速简化为3个随机变量，利用基于虚位移原理的随机有限元方法，对汽轮机转子的随机静态响应特性进行分析。

转子的参数特性如表1所示。

首先，读入节点、单元、材料、载荷比值等信息，利用Fortran语言进行编程，再采用Grid 9．0处理技术显示运行结果。

图4和图5为3个随机变量同时作用时汽轮机转子的应力、变形的均值和方差云图。

由图4和图5可知：汽轮机转子的应力、变形的均值和方差在转子上的分布类型近乎一致。

利用直接Monte－Carlo模拟仿真法进行验证计算，抽样数为1000组，通过计算，比较结果如表2和表3所示。

从表2和表3可知：2种方法计算的结果吻合度较高，从而验证了本文方法的正确性。另外，弹性模量对转子的最大应力没有影响，对转子最大应力影响最大的是密度，其次是角速度；对转子最大变形影响最大的是密度，其次是弹性模量和角速度。

5 小结

1）在对八节点四边形等参单元和轴对称问题研究的基础上，本文提出了采用八节点四边形轴对称等参单元解决轴对称问题。根据偏微分方程弱解的积分形式推导出线弹性小变形静态轴对称等参问题的“弱”形式，结合摄动技术推导出了轴对称等参问题的静力学随机有限元方程。

2）以汽轮机转子为研究对象，根据推导的随机有限元方程分别计算汽轮机转子变形和应力的均值和方差，并与直接Monte－Carlo模拟仿真法的计算结果进行比较，验证了本文方法的正确性。

［1］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．

［2］Pellissetti M F，Schueller G L．Scalable uncertainty and reliability analysis by integration of advanced Monte Carlo simulation and generic finite element solvers［J］．Computer and Structures，2009，87：937－947．

［3］王东，陈建康，王启智，等．Monte－Carlo随机有限元结构可靠度分析新方法［J］．四川大学学报工程科学报，2008，40（3）：20－26．

［4］段巍，王璋奇．基于响应面方法的汽轮机叶片概率强度设计及敏感性分析［J］．中国电机工程学报，2007，27（5）：99－104．

［5］张旭方，PANDEY Mahesh D．张义民．结构随机响应计算的一种数值方法［J］．中国科学：E 辑，2012，42（1）：103－114．

［6］黄会荣，朱怡婕，郭家元，等．基于APDL的钢板剪力墙可靠性研究［J］．应用力学学报，2012，29（2）：210－215．

［7］刘效勇，卢佩，黄旭初，等．基于有限元方法光学材料热损伤的研究［J］．石河子大学学报：自然科学版，2011，29（1）：125－128．

［8］朱健文，王卫兵．基于Solidworks与ANSYS的机器鱼壳体有限元分析［J］．石河子大学学报：自然科学版，2011，29（4）：509－513．

［9］卢勇鹏．基于有限元法对乌拉斯台水库浆砌石拱坝的除险加固效应分析研究［D］．新疆：石河子大学，2011．

［10］安伟光，朱卫兵，等．随机有限元法在不确定性分析中的应用［J］．哈尔滨工程大学学报，2002，23（1）：132－135．

［11］安立强，王璋奇．随机约束汽轮机叶片频率的有限元分析［J］．中国电机工程学报，2009，29（2）：95－100．

［12］戎志祥，林少芬．基于随机有限元法的连杆可靠性分析［J］．舰船科学技术，2011，33（9）：68－71．

［13］薛小锋，冯蕴雯，冯元生．基于随机有限元法的平面多裂纹结构可靠性研究［J］．西北工业大学学报，2012，30（4）：508－513．

［14］ Marcin Kaminski．Perturbation－based stochastic finite element method using polynomial response function for the elastic beams［J］．Mechanics Research Communications，2009，36：381－390．

［15］胡海昌．弹性力学的变分原理及其应用［M］．北京：科学出版社，1981．

［16］王勖成，邵敏．有限单元法基本原理和数值方法［D］．北京：清华大学出版社，1997．

［17］刘文廷．结构可靠性设计手册［M］．北京：国防工业出版社，2008．

［18］安利强．汽轮机叶片静动特性的随机有限元方法研究［D］．保定：华北电力大学，2005．

［19］曾攀．有限元方法基本原理［M］．北京：清华大学出版社，2008．