周宗和，尹喜慶

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

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

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

1 八節(jié)點(diǎn)四邊形軸對(duì)稱等參單元

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

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

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

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

1）位移模式。

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

寫(xiě)成矩陣的形式為：

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

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

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

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

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

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

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

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

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

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

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

由此可得出軸對(duì)稱等參問(wèn)題的“弱”形式如下：

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

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

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

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

3 軸對(duì)稱等參問(wèn)題的靜力學(xué)隨機(jī)有限元方程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

忽略位移的二階及二階以上攝動(dòng)展開(kāi)式，進(jìn)行協(xié)方差運(yùn)算，可以得到位移的協(xié)方差矩陣

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

將單元應(yīng)力在基本隨機(jī)變量的均值點(diǎn)處進(jìn)行二階Taylor級(jí)數(shù)展開(kāi)，得

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

將位移的一階變異量寫(xiě)成矩陣形式

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

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

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

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

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

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

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

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

隨機(jī)場(chǎng)的相關(guān)長(zhǎng)度越長(zhǎng)，表示隨機(jī)場(chǎng)的空間兩點(diǎn)之間的相關(guān)性越強(qiáng)，當(dāng)隨機(jī)場(chǎng)的相關(guān)長(zhǎng)度等于無(wú)窮大時(shí)，隨機(jī)場(chǎng)退化為一個(gè)隨機(jī)變量。

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

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

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

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

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

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

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

5 小結(jié)

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

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

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

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

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

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

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

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

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

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

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

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

［13］薛小鋒，馮蘊(yùn)雯，馮元生．基于隨機(jī)有限元法的平面多裂紋結(jié)構(gòu)可靠性研究［J］．西北工業(yè)大學(xué)學(xué)報(bào)，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è)計(jì)手冊(cè)［M］．北京：國(guó)防工業(yè)出版社，2008．

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

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