基于虛位移原理的隨機有限元方法研究

周宗和，尹喜慶

（海軍駐武漢四三八廠軍事代表室，武漢430060）

目前，研究不確定性結(jié)構(gòu)的常用方法可分為兩類。一類是統(tǒng)計的方法［1－9］，就是通過大量的隨機抽樣，對結(jié)構(gòu)反復(fù)進行有限元計算，最后對得到的結(jié)果作統(tǒng)計分析，這類方法稱為Monte－Carlo隨機有限元法，需要反復(fù)進行大量的模擬計算，并不適用于大型復(fù)雜結(jié)構(gòu)。另一類方法［10－14］是以數(shù)學(xué)、力學(xué)分析作為工具，得出結(jié)構(gòu)系統(tǒng)的響應(yīng)與輸入信號之間的關(guān)系，并據(jù)此得到結(jié)構(gòu)應(yīng)力或位移的統(tǒng)計規(guī)律，這類方法有攝動隨機有限元法、Neumann隨機有限元法等，這類方法在大型復(fù)雜結(jié)構(gòu)中的應(yīng)用也非常有限。

胡海昌［15］從非統(tǒng)計方法出發(fā)，通過對偏微分方程弱解的積分形式——虛位移原理進行較深入的研究，王勖成、邵敏［16］和劉文廷［17］構(gòu)造了八節(jié)點四邊形軸對稱等參單元。本文在此基礎(chǔ)上給出軸對稱靜態(tài)問題的平衡方程以及力的邊界條件下的等效積分“弱”形式，結(jié)合攝動技術(shù)推導(dǎo)出八節(jié)點四邊形軸對稱等參單元的隨機靜態(tài)有限元方程，即基于虛位移原理的隨機有限元方法，最后考慮汽輪機轉(zhuǎn)子的物性參數(shù)和載荷參數(shù)的隨機性，利用這種方法對汽輪機轉(zhuǎn)子的靜態(tài)隨機響應(yīng)特性進行計算和分析，并將計算結(jié)果與通過直接Monte－Carlo模擬仿真法計算得出的結(jié)果相對比，以驗證該方法的正確性。

1 八節(jié)點四邊形軸對稱等參單元

八節(jié)點四邊形軸對稱等參單元如圖1所示。

該單元為橫截面為八節(jié)點四邊形的360°環(huán)形等參單元，其橫截面上八個節(jié)點的位置坐標(biāo)（ri，zi），i＝1，2，…，8，各個節(jié)點的位移為（ui，wi），i＝1，2，…，8。

該單元為繞z軸的環(huán)狀單元，在orz平面內(nèi)，八節(jié)點四邊形軸對稱等參單元共有16個自由度。將所有節(jié)點上的位移組成一個列陣，記作｛｝ae，同樣，將所有節(jié)點上的各個力也組成一個列陣，記作｛｝pe，那么

式（1）中，采用的等參變換位移模式和應(yīng)變矩陣計算式如下：

1）位移模式。

式（3）中：［N］為形函數(shù)矩陣。

寫成矩陣的形式為：

2）應(yīng)變矩陣。

由軸對稱問題的幾何方程可以推出相應(yīng)的幾何函數(shù)矩陣，得

式（5）中：［B］為幾何矩陣。

2 線彈性小變形靜態(tài)軸對稱等參問題的“弱”形式

對于線彈性小變形靜態(tài)軸對稱問題，運用虛位移原理求位移，由平衡方程［16］得

式（6）中：δu，δw為兩個方向的虛位移。

由分部積分公式將其化為“弱”形式，利用幾何方程可得

式（7）中：l為單元體的邊界。

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

將本構(gòu)方程的矩陣表示代入上式可得

將力和位移邊界條件代入上式（對于位移邊界條件，虛位移δu，δw為0，可得

由此可得出軸對稱等參問題的“弱”形式如下：

式（11）中：［J］為雅可比矩陣，其中

式（12）中：S為坐標(biāo)變換后求解域的面元，L為坐標(biāo)變換后的單元體的邊界。

將式（4）和（5）代入式（11），整理后得出軸對稱等參問題的“弱”形式，具體如下：

由此可推出零階隨機積分“弱”形式和一階隨機積分“弱”形式。

3 軸對稱等參問題的靜力學(xué)隨機有限元方程

在許多實際問題當(dāng)中，影響結(jié)構(gòu)靜力學(xué)特性的因素，如材料特性、幾何特性和載荷特性等都有一定程度隨機性，因此結(jié)構(gòu)的響應(yīng)｛u｝和｛ε｝也具有隨機性。本文為了簡化書寫，假設(shè)任意一個給定的空間連續(xù)函數(shù)f（｛r0｝），｛r0｝為空間的坐標(biāo)向量，基本隨機向量｛b｝的聯(lián)合概率密度函數(shù)為g（｛b｝），εp為一個小參數(shù)，根據(jù)向量的數(shù)學(xué)期望計算公式有下：

將式（14）代入式（13），并比較小參數(shù)εp的同次冪系數(shù)，可以得到結(jié)構(gòu)攝動格式的各階隨機積分“弱”形式，其中零階和一階隨機積分“弱”形式如下：

1）零階隨機積分“弱”形式。其方程如下：

2）一階隨機積分“弱”形式。其方程如下：

基于攝動理論［18］將影響結(jié)構(gòu)的隨機參數(shù)在基本隨機變量｛b｝的均值｛b0｝處作小參數(shù)攝動展開，如

式（15）、（16）中，δ｛a｝eT是δ｛a｝e的轉(zhuǎn)置向量。

由于δ｛a｝e是任意的，故由式（16）可得到單元的零階靜態(tài)隨機有限元方程：

式（17）中：［K0］e為單元的線彈性剛度矩陣均值，｛F0｝e為單元的體積力均值，｛T0F｝e為單元的面力向量均值，其相應(yīng)的單元矩陣如下：

求解區(qū)域離散為n個單元，經(jīng)組集后可以得到結(jié)構(gòu)的零階靜態(tài)隨機有限元方程：

式（19）中：［K0］為結(jié)構(gòu)的線彈性剛度矩陣均值，｛F0｝為結(jié)構(gòu)的體積力均值，｛T0F｝為結(jié)構(gòu)的邊界面力向量均值。

同理，由式（17）可得到單元的一階靜態(tài)隨機有限元方程：

式（20）中：［K］e′j為單元的線彈性剛度矩陣的變異矩陣，｛F｝e′j為單元的體積力向量的變異向量，｛TF｝e′j為單元的邊界面力向量的變異向量，分別如下表示：

求解區(qū)域離散為n個單元，經(jīng)組集后可以得到結(jié)構(gòu)的一階靜態(tài)隨機有限元方程：

4 汽輪機轉(zhuǎn)子靜態(tài)隨機有限元分析

由于汽輪機轉(zhuǎn)子的材料參數(shù)、幾何參數(shù)和載荷參數(shù)都存在隨機性，所以轉(zhuǎn)子的變形和應(yīng)力也存在隨機性。因此，本文采用基于虛位移原理的隨機有限元方法和直接 Monte－Carlo模擬仿真法［19］2種方法分別對汽輪機轉(zhuǎn)子的靜態(tài)隨機響應(yīng)特性進行研究，計算汽輪機轉(zhuǎn)子變形和應(yīng)力的均值和方差，對這2種方法的計算結(jié)果進行比較，以驗證基于虛位移原理的隨機有限元方法的正確性。

4．1 汽輪機轉(zhuǎn)子的隨機有限元方程

對推導(dǎo)出的軸對稱等參問題的靜態(tài)隨機有限元方程求解，可以得到位移的統(tǒng)計特性。

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

式（25）中：［K0］為汽輪機轉(zhuǎn)子的線彈性剛度矩陣均值，｛F0｝為汽輪機轉(zhuǎn)子離心力向量均值，｛T0F｝為汽輪機轉(zhuǎn)子的邊界面力向量均值。

求解式（25），可以得到位移的一階變異量

位移一階變異量的矩陣形式為：

式（26）稱為各節(jié)點位移相對于各隨機變量的靈敏度矩陣，其中m為基本隨機變量的個數(shù)，S｛a｝ij表示第i個位移分量對第j個基本隨機變量的靈敏度。

忽略位移的二階及二階以上攝動展開式，進行協(xié)方差運算，可以得到位移的協(xié)方差矩陣

式（28）中：［C｛b｝］為基本隨機變量的協(xié)方差矩陣。求出轉(zhuǎn)子的節(jié)點位移后，則轉(zhuǎn)子的單元應(yīng)力為：

將單元應(yīng)力在基本隨機變量的均值點處進行二階Taylor級數(shù)展開，得

比較小參數(shù)εp的同次冪系數(shù)，分別得出單元應(yīng)力的均值和一階變異量：

將位移的一階變異量寫成矩陣形式

式（33）稱為單元應(yīng)力相對于各隨機變量的靈敏度矩陣。同樣，可以得到應(yīng)力的協(xié)方差矩陣

式（34）中：［C｛b｝］為基本隨機變量的協(xié)方差矩陣；［S｛σ｝］ij表示第i個應(yīng)力分量對第j個基本隨機變量的靈敏度。

4．2 汽輪機轉(zhuǎn)子的二維幾何模型

由于汽輪機轉(zhuǎn)子本身是一個復(fù)雜的部件，所以在不影響轉(zhuǎn)子計算準(zhǔn)確性的前提下，應(yīng)盡量使計算趨于簡單化，因此，將其簡化為軸對稱問題，同時，為了減少邊界條件的設(shè)定誤差對計算結(jié)果可能造成的影響，建模時取轉(zhuǎn)子在汽缸內(nèi)的整段剖面圖作為研究對象。

轉(zhuǎn)子的結(jié)構(gòu)、形狀、尺寸及材料取自某型船設(shè)計資料中的轉(zhuǎn)子加工用圖以及總裝配圖，利用網(wǎng)格劃分軟件Grid 9．0對汽輪機轉(zhuǎn)子進行二維建模，如圖2所示。

4．3 汽輪機轉(zhuǎn)子的隨機有限元程序流程

汽輪機轉(zhuǎn)子的隨機靜態(tài)有限元分析程序流程圖如圖3所示。

4．4 汽輪機轉(zhuǎn)子的隨機靜態(tài)響應(yīng)特性分析

隨機場的相關(guān)長度越長，表示隨機場的空間兩點之間的相關(guān)性越強，當(dāng)隨機場的相關(guān)長度等于無窮大時，隨機場退化為一個隨機變量。

本文把轉(zhuǎn)子的彈性模量、密度和轉(zhuǎn)速簡化為3個隨機變量，利用基于虛位移原理的隨機有限元方法，對汽輪機轉(zhuǎn)子的隨機靜態(tài)響應(yīng)特性進行分析。

轉(zhuǎn)子的參數(shù)特性如表1所示。

首先，讀入節(jié)點、單元、材料、載荷比值等信息，利用Fortran語言進行編程，再采用Grid 9．0處理技術(shù)顯示運行結(jié)果。

圖4和圖5為3個隨機變量同時作用時汽輪機轉(zhuǎn)子的應(yīng)力、變形的均值和方差云圖。

由圖4和圖5可知：汽輪機轉(zhuǎn)子的應(yīng)力、變形的均值和方差在轉(zhuǎn)子上的分布類型近乎一致。

利用直接Monte－Carlo模擬仿真法進行驗證計算，抽樣數(shù)為1000組，通過計算，比較結(jié)果如表2和表3所示。

從表2和表3可知：2種方法計算的結(jié)果吻合度較高，從而驗證了本文方法的正確性。另外，彈性模量對轉(zhuǎn)子的最大應(yīng)力沒有影響，對轉(zhuǎn)子最大應(yīng)力影響最大的是密度，其次是角速度；對轉(zhuǎn)子最大變形影響最大的是密度，其次是彈性模量和角速度。

5 小結(jié)

1）在對八節(jié)點四邊形等參單元和軸對稱問題研究的基礎(chǔ)上，本文提出了采用八節(jié)點四邊形軸對稱等參單元解決軸對稱問題。根據(jù)偏微分方程弱解的積分形式推導(dǎo)出線彈性小變形靜態(tài)軸對稱等參問題的“弱”形式，結(jié)合攝動技術(shù)推導(dǎo)出了軸對稱等參問題的靜力學(xué)隨機有限元方程。

2）以汽輪機轉(zhuǎn)子為研究對象，根據(jù)推導(dǎo)的隨機有限元方程分別計算汽輪機轉(zhuǎn)子變形和應(yīng)力的均值和方差，并與直接Monte－Carlo模擬仿真法的計算結(jié)果進行比較，驗證了本文方法的正確性。

［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隨機有限元結(jié)構(gòu)可靠度分析新方法［J］．四川大學(xué)學(xué)報工程科學(xué)報，2008，40（3）：20－26．

［4］段巍，王璋奇．基于響應(yīng)面方法的汽輪機葉片概率強度設(shè)計及敏感性分析［J］．中國電機工程學(xué)報，2007，27（5）：99－104．

［5］張旭方，PANDEY Mahesh D．張義民．結(jié)構(gòu)隨機響應(yīng)計算的一種數(shù)值方法［J］．中國科學(xué)：E 輯，2012，42（1）：103－114．

［6］黃會榮，朱怡婕，郭家元，等．基于APDL的鋼板剪力墻可靠性研究［J］．應(yīng)用力學(xué)學(xué)報，2012，29（2）：210－215．

［7］劉效勇，盧佩，黃旭初，等．基于有限元方法光學(xué)材料熱損傷的研究［J］．石河子大學(xué)學(xué)報：自然科學(xué)版，2011，29（1）：125－128．

［8］朱健文，王衛(wèi)兵．基于Solidworks與ANSYS的機器魚殼體有限元分析［J］．石河子大學(xué)學(xué)報：自然科學(xué)版，2011，29（4）：509－513．

［9］盧勇鵬．基于有限元法對烏拉斯臺水庫漿砌石拱壩的除險加固效應(yīng)分析研究［D］．新疆：石河子大學(xué)，2011．

［10］安偉光，朱衛(wèi)兵，等．隨機有限元法在不確定性分析中的應(yīng)用［J］．哈爾濱工程大學(xué)學(xué)報，2002，23（1）：132－135．

［11］安立強，王璋奇．隨機約束汽輪機葉片頻率的有限元分析［J］．中國電機工程學(xué)報，2009，29（2）：95－100．

［12］戎志祥，林少芬．基于隨機有限元法的連桿可靠性分析［J］．艦船科學(xué)技術(shù)，2011，33（9）：68－71．

［13］薛小鋒，馮蘊雯，馮元生．基于隨機有限元法的平面多裂紋結(jié)構(gòu)可靠性研究［J］．西北工業(yè)大學(xué)學(xué)報，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］胡海昌．彈性力學(xué)的變分原理及其應(yīng)用［M］．北京：科學(xué)出版社，1981．

［16］王勖成，邵敏．有限單元法基本原理和數(shù)值方法［D］．北京：清華大學(xué)出版社，1997．

［17］劉文廷．結(jié)構(gòu)可靠性設(shè)計手冊［M］．北京：國防工業(yè)出版社，2008．

［18］安利強．汽輪機葉片靜動特性的隨機有限元方法研究［D］．保定：華北電力大學(xué)，2005．

［19］曾攀．有限元方法基本原理［M］．北京：清華大學(xué)出版社，2008．