隨機(jī)有限元法

1、 第7章 隨機(jī)有限元法§7.1 緒論結(jié)構(gòu)工程中存在諸多的不確定性因素，從結(jié)構(gòu)材料性能參數(shù)到所承受的主要荷載，如車(chē)流、陣風(fēng)或地震波，無(wú)不存在隨機(jī)性。在有限單元法已成為分析復(fù)雜結(jié)構(gòu)的強(qiáng)有力的工具和廣泛使用的數(shù)值方法的今天，人們已不滿(mǎn)足精度越來(lái)越高的確定性有限元計(jì)算，而設(shè)法用這一強(qiáng)有力的工具去研究工程實(shí)踐中存在的大量不確定問(wèn)題。隨機(jī)有限元法（Stochastic FEM）,也稱(chēng)概率有限元法（Probabilistic FEM）正是隨機(jī)分析理論與有限元方法相結(jié)合的產(chǎn)物，是在傳統(tǒng)的有限元方法的基礎(chǔ)上發(fā)展起來(lái)的隨機(jī)的數(shù)值分析方法。最初是Monte-Carlo法與有限元法直接結(jié)合，形成獨(dú)特的統(tǒng)計(jì)有

2、限元方法。Astill和Shinozuka（1972）首先將Monte-Carlo法引入結(jié)構(gòu)的隨機(jī)有限元法分析。該法通過(guò)在計(jì)算機(jī)上產(chǎn)生的樣本函數(shù)來(lái)模擬系統(tǒng)的隨機(jī)輸入量的概率特征，并對(duì)于每個(gè)給定的樣本點(diǎn)，對(duì)系統(tǒng)進(jìn)行確定性的有限元分析，從而得到系統(tǒng)的隨機(jī)響應(yīng)的概率特征。由于是直接建立在大量確定性有限元計(jì)算的基礎(chǔ)上，計(jì)算量極大，不適用于大型結(jié)構(gòu)，而且最初的直接Monte-Carlo法還不是真正意義上的隨機(jī)有限元法。但與隨后的攝動(dòng)隨機(jī)有限元法（PSFEM）相比，當(dāng)樣本容量足夠大時(shí)，Monte-Carlo有限元法的結(jié)果更可靠也更精確。結(jié)構(gòu)系統(tǒng)的隨機(jī)分析一般可分為兩大類(lèi)：一類(lèi)是統(tǒng)計(jì)方法，另一類(lèi)是非統(tǒng)計(jì)方法

3、。因此，隨機(jī)有限元法同樣也有統(tǒng)計(jì)逼近和非統(tǒng)計(jì)逼近兩種類(lèi)型。前者通過(guò)樣本試驗(yàn)收集原始的數(shù)據(jù)資料，運(yùn)用概率和統(tǒng)計(jì)理論進(jìn)行分析和整理，然后作出科學(xué)推斷。這里，樣本試驗(yàn)和數(shù)據(jù)處理的工作量很大，隨著計(jì)算機(jī)的普及和發(fā)展，數(shù)值模擬法，如蒙特卡羅（Monte Carlo）模擬，已成為最常用的統(tǒng)計(jì)逼近法。后者從本質(zhì)上來(lái)說(shuō)是利用分析工具找出結(jié)構(gòu)系統(tǒng)的（確定的或隨機(jī)的）輸出隨機(jī)信號(hào)與輸入隨機(jī)信號(hào)之間的關(guān)系，采用隨機(jī)分析與求解系統(tǒng)控制方程相結(jié)合的方法得到輸出信號(hào)的各階隨機(jī)統(tǒng)計(jì)量的數(shù)字特征（如各階原點(diǎn)矩或中心矩）。在20世紀(jì)70年代初， Cambou首先采用一次二階矩方法研究線(xiàn)彈性問(wèn)題。由于這種方法將隨機(jī)變量的影響量進(jìn)

4、行Taylor級(jí)數(shù)展開(kāi)，就稱(chēng)之為T(mén)aylor展開(kāi)法隨機(jī)有限元（TSFEM）。Shinozuka和Astill（1972)分別獨(dú)立運(yùn)用攝動(dòng)技術(shù)研究了隨機(jī)系統(tǒng)的特征值問(wèn)題。隨后，Handa（1975）等人在考慮隨機(jī)變量波動(dòng)性時(shí)采用一階和二階攝動(dòng)技術(shù)，并將這種攝動(dòng)法隨機(jī)有限元成功地應(yīng)用于框架結(jié)構(gòu)分析。Vanmarcke等人（1983）提出隨機(jī)場(chǎng)的局部平均理論，并將它引入隨機(jī)有限元。局部平均理論是用隨機(jī)場(chǎng)函數(shù)在每一個(gè)離散單元上的局部平均的隨機(jī)變量來(lái)代表該單元的統(tǒng)計(jì)量的近似理論。Liu W. K.等人（1986、1988）的系列工作，提供了一種“主模態(tài)”技術(shù)，運(yùn)用隨機(jī)變量的特征正交化方法，將滿(mǎn)秩的協(xié)方差

5、矩陣變換為對(duì)角矩陣，減少計(jì)算工作量，對(duì)攝動(dòng)隨機(jī)有限元法的發(fā)展做出貢獻(xiàn)，此外，提出了一個(gè)隨機(jī)變分原理。Yamazaki和Shinozuka(1987)創(chuàng)造性地將算子的Neumann級(jí)數(shù)展開(kāi)式引入隨機(jī)有限元的列式工作。從本質(zhì)上講，Neumann級(jí)數(shù)展開(kāi)方法也是一類(lèi)正則的小參數(shù)攝動(dòng)方法，正定的隨機(jī)剛度矩陣和微小的隨機(jī)擾動(dòng)量是兩個(gè)基本要求，這兩個(gè)基本要求保證了攝動(dòng)解的正則性和收斂性，其優(yōu)點(diǎn)在于攝動(dòng)形式較簡(jiǎn)單并可以得到近似解的高階統(tǒng)計(jì)量。Shinozuka等人（1987）將隨機(jī)場(chǎng)函數(shù)的Monte-Carlo模擬與隨機(jī)剛度矩陣的Neumann級(jí)數(shù)展開(kāi)式結(jié)合,得到具有較好計(jì)算精度和效率的一類(lèi)Neumann隨

6、機(jī)有限元列式（稱(chēng)NSFEM）。Benaroya等（1988）指出，將出現(xiàn)以隨機(jī)變分原理為基礎(chǔ)的隨機(jī)有限元法來(lái)逐漸取代以攝動(dòng)法為基礎(chǔ)的隨機(jī)有限元法。Spanos和Ghanem等人（1989，1991）結(jié)合隨機(jī)場(chǎng)函數(shù)的Karhuen-Loeve展式和Galerkin（迦遼金）射影方法建立了相應(yīng)的隨機(jī)有限元列式，并撰寫(xiě)了隨機(jī)有限元法領(lǐng)域的第一本專(zhuān)著隨機(jī)有限元譜方法。國(guó)內(nèi)對(duì)隨機(jī)有限元的研究起步較晚。吳世偉等人（1988）提出隨機(jī)有限元的直接偏微分法及相應(yīng)的可靠度計(jì)算方法。陳虬、劉先斌等人（1989、1991）提出一種新的隨機(jī)場(chǎng)離散模型，建立了等參局部平均單元，并基于變分原理研究了一類(lèi)隨機(jī)有限元法的收斂

7、性和誤差界。Papadrakakis(1995)采用預(yù)處理共軛梯度法給出了空間框架的非線(xiàn)性隨機(jī)有限元列式。Schorling和Bucher(1996)基于Monte-Carlo技術(shù)，采用響應(yīng)面法研究幾何非線(xiàn)性時(shí)的可靠度隨機(jī)有限元方法。劉寧（1996）則基于偏微分法，給出了三維彈塑性隨機(jī)有限元列式。隨機(jī)有限元法的數(shù)學(xué)理論研究和非線(xiàn)性隨機(jī)問(wèn)題的有限元分析工作還有待深入。自20世紀(jì)80年代以來(lái)，隨機(jī)有限元法已在工程結(jié)構(gòu)可靠性、安全性分析領(lǐng)域以及在各種隨機(jī)激勵(lì)下結(jié)構(gòu)響應(yīng)變異研究領(lǐng)域中得到應(yīng)用，如應(yīng)用于大型水利工程的重力壩、拱壩的可靠度計(jì)算；應(yīng)用于非線(xiàn)性瞬態(tài)響應(yīng)分析；結(jié)構(gòu)振動(dòng)中隨機(jī)阻尼對(duì)響應(yīng)的影響；結(jié)構(gòu)

8、分析的隨機(jī)識(shí)別；復(fù)雜結(jié)構(gòu)地震響應(yīng)的隨機(jī)分析和兩相動(dòng)力系統(tǒng)的隨機(jī)模擬等等。隨著理論研究的深入，隨機(jī)有限元將得到更加廣泛的應(yīng)用。§7.2 隨機(jī)有限元的控制方程22從隨機(jī)有限元控制方程的獲得來(lái)看，隨機(jī)有限元可分為T(mén)aylor展開(kāi)法隨機(jī)有限元（TSFEM）、攝動(dòng)法隨機(jī)有限元(PSFEM)以及Neumann展開(kāi)Monte-Carlo法隨機(jī)有限元（NSFEM）。 Taylor展開(kāi)法隨機(jī)有限元 該隨機(jī)有限元法的基本思路是將有限元格式中的控制量在隨機(jī)變量均值點(diǎn)處進(jìn)行Taylor級(jí)數(shù)展開(kāi)（取一階或二階），經(jīng)過(guò)適當(dāng)?shù)臄?shù)學(xué)處理得出所需的計(jì)算方程式。有限元靜力分析控制方程的矩陣形式為: KU = F (7.

9、2.1)式中,U為位移矩陣，F(xiàn)為等效節(jié)點(diǎn)荷載列陣，K為整體剛度矩陣 (7.2.2)其中，B為形變矩陣,D為材料彈性矩陣。在計(jì)算出節(jié)點(diǎn)位移U后，即由下式求得應(yīng)力列陣= DBU (7.2.3)設(shè)基本隨機(jī)變量為，將位移U在均值點(diǎn)處一階Taylor級(jí)數(shù)展開(kāi)，并在兩邊同時(shí)取均值（數(shù)學(xué)期望），得 (7.2.4)式中：符號(hào)E·表示求均值，任一結(jié)點(diǎn)位移U的方差可由下式計(jì)算： (7.2.5)式中：符號(hào)Var·表示求方差；Cov(Xi,Xj)為Xi和Xj的協(xié)方差。其中 (7.2.6) (7.2.7)同樣將在均值點(diǎn)處Taylor展開(kāi)，也有與上面類(lèi)似的表達(dá)式?？梢?jiàn)，TSFEM關(guān)鍵在于對(duì)有限元方程式

10、直接進(jìn)行偏微分計(jì)算，計(jì)算出有限元輸出量對(duì)隨機(jī)變量的梯度，故該法也稱(chēng)直接偏微分法或梯度分析法。由于一階TSFEM只需一次形成剛度矩陣，也只需一次求剛度矩陣的逆，因此效率較高。但由于忽略了二階以上的高次項(xiàng)，使TSFEM對(duì)隨機(jī)變量的變異性有所限制。一般要求一階TSFEM隨機(jī)變量的變異系數(shù)小于0.3。如果隨機(jī)變量的變異系數(shù)較大，可以采用有限元控制方程的二階Taylor展開(kāi)： (7.2.8) (7.2.9)上式可見(jiàn)，二階TSFEM可以放寬隨機(jī)變量變異性大小的限制，但隨機(jī)變量數(shù)目較多時(shí)，計(jì)算量將十分龐大，而且一階或二階TSFEM均無(wú)法計(jì)算響應(yīng)量三階以上的統(tǒng)計(jì)特性。由于TSFEM簡(jiǎn)單明了、效率高，為我國(guó)許多

11、學(xué)者所采用。 攝動(dòng)法隨機(jī)有限元攝動(dòng)技術(shù)最初被用于非線(xiàn)性力學(xué)分析。Handa等人成功地將一階、二階攝動(dòng)技術(shù)用于隨機(jī)問(wèn)題，給出攝動(dòng)法有限元列式。該法假定基本隨機(jī)變量在均值點(diǎn)處產(chǎn)生微小攝動(dòng)，利用Taylor級(jí)數(shù)把隨機(jī)變量表示為確定部分和由攝動(dòng)引起的隨機(jī)部分，從而將有限元控制方程（非線(xiàn)性的）轉(zhuǎn)化為一組線(xiàn)性的遞推方程，求解得出位移的統(tǒng)計(jì)特性，進(jìn)而求出應(yīng)力的統(tǒng)計(jì)特性。假設(shè)為隨機(jī)變量在均值點(diǎn)處的微小攝動(dòng)量，即。于是 (7.2.10)對(duì)于U、F，也有類(lèi)似上式K的表達(dá)式，式中：K0、U0、F0分別為K、U、F在隨機(jī)變量均值點(diǎn)的值。根據(jù)二階攝動(dòng)法，可得 (7.2.11) (7.2.12) (7.2.13)由上式可

12、得位移的均值和協(xié)方差： (7.2.14) (7.2.15)由于任何量的隨機(jī)性都可以引入攝動(dòng)量，而且更易于考慮非線(xiàn)性問(wèn)題，因此PSEFEM適用范圍較廣，對(duì)于結(jié)構(gòu)幾何特性的隨機(jī)性（包括隨機(jī)邊界問(wèn)題）易得出隨機(jī)有限元控制方程。一階PSFEM和一階TSFEM一樣，只需一次形成剛度矩陣、一次對(duì)剛度矩陣求逆，計(jì)算效率較高。但PSFEM需以微小的攝動(dòng)量為條件，一般應(yīng)小于均值的20%或30%。 Neumann展開(kāi)Monte-Carlo隨機(jī)有限元20世紀(jì)80年代后期，Shinozuka等人提出基于Neumann展開(kāi)式的隨機(jī)有限元法，使Monte-Carlo法與有限元法較完美地結(jié)合起來(lái)。Monte-Carlo法是

13、最直觀(guān)、最精確、獲取信息最多、對(duì)非線(xiàn)性問(wèn)題最有效的計(jì)算統(tǒng)計(jì)方法。Neumann展開(kāi)式的引入是為了解決矩陣求逆的效率問(wèn)題。如果對(duì)每一次隨機(jī)抽樣，只需形成剛度矩陣，進(jìn)行前代、回代以及矩陣乘和矩陣加減，而無(wú)需矩陣分解，則可大大減少工作量。在一般有限元控制方程KU = F中，假定荷載F為確定值，在隨機(jī)變量波動(dòng)值的影響下剛度矩陣K分解為K = K0+K，根據(jù)Neumann級(jí)數(shù)展開(kāi)，有 K-1=(K0+K)-1=(I-P+P2-P3+)K0-1 (7.2.16)式中：K0為隨機(jī)變量均值處的剛度矩陣；K為剛度矩陣的波動(dòng)量；I為單位矩陣。對(duì)于Monte-Carlo隨機(jī)抽樣，剛度矩陣只改變K項(xiàng)，而P = K0-

14、1K (7.2.17)U0 = K0-1F (7.2.18)將式（7.2.16）代入式（7.2.1），并利用式（7.2.17）和式（7.2.18），得 U = U0-PU0+P2U0-P3U0+ (7.2.19)令Ui=PiU0，則得如下的遞推公式： U = U0+U1+U2+ (7.2.20) Ui=-K0-1KUi-1 (i=1,2,3,) (7.2.21)由式（7.2.18）求出U0后，可由式（7.2.21）求出Ui(i=1,2,3,)。 上述三種方法中，NSFEM可以方便地調(diào)用確定性有限元的計(jì)算程序，而TSFEM在編程上較為復(fù)雜，PSFEM則更為復(fù)雜。由于采用Monte-Carlo隨機(jī)

15、模擬技術(shù)，NSFEM不受隨機(jī)變量波動(dòng)范圍的限制，當(dāng)變異系數(shù)小于0.2時(shí)，NSFEM與一階TSFEM或一階PSFEM精度相當(dāng)；當(dāng)變異系數(shù)大于0.2時(shí)，后兩者已不能滿(mǎn)足精度要求，但NSFEM仍能得出滿(mǎn)意的結(jié)果。§7.3 隨機(jī)場(chǎng)的離散模型許多物理現(xiàn)象和物體系統(tǒng)具有隨機(jī)分布特性，包括系統(tǒng)本身的不確定或系統(tǒng)的激勵(lì)和響應(yīng)的不確定，都可以模型化為隨機(jī)空間分布的隨機(jī)場(chǎng)或隨時(shí)間分布的隨機(jī)過(guò)程。隨機(jī)有限元法除了必須考慮材料參數(shù)等的空間變異性，需要獲得隨機(jī)有限元方程列式以及解決隨機(jī)算子和隨機(jī)矩陣的求逆問(wèn)題外，還須包含對(duì)隨機(jī)場(chǎng)的離散處理。均勻各向同性隨機(jī)場(chǎng)的特征量1. 隨機(jī)場(chǎng)S(t)的均值E(S(t)為常數(shù)

16、m.2. 隨機(jī)場(chǎng)S(t)的標(biāo)準(zhǔn)相關(guān)函數(shù) (7.3.1)式中，隨機(jī)場(chǎng)的協(xié)方差函數(shù) R()= CovS(t+),S(t) (7.3.2)對(duì)于一切t，隨機(jī)場(chǎng)的方差為 VarS(t)= R(0)= 2 (7.3.3)相關(guān)函數(shù)也可以用譜分解表示（即Wiener-Khintchine變換對(duì)）： (7.3.4) 常用的標(biāo)準(zhǔn)相關(guān)函數(shù)有：非協(xié)調(diào)階躍型、協(xié)調(diào)階躍型、三角型、指數(shù)型、二階AR型、高斯型等多種形式。3. 方差折減系數(shù)2(h)設(shè)Sh(t)是隨機(jī)場(chǎng)S(t)的局部平均隨機(jī)過(guò)程，即 (7.3.5)則方差 VarSh(t)= 22(h) (7.3.6)式中，方差折減系數(shù) (7.3.7)可見(jiàn)，方差折減系數(shù)起著使S

17、h(t)的方差比原來(lái)S(t)的方差縮小的作用。4. 相關(guān)距離SS可以看成是任意兩個(gè)相隔距離為的隨機(jī)變量不相關(guān)的最小距離（也稱(chēng)相關(guān)偏度）。，則不相關(guān)，否則完全相關(guān)。利用相關(guān)距離S便于對(duì)隨機(jī)場(chǎng)作近似處理，其計(jì)算式為： (7.3.8) 以上公式均對(duì)一維問(wèn)題列出，二維、三維問(wèn)題也可以類(lèi)似得出。 隨機(jī)場(chǎng)的中心離散中心離散是隨機(jī)場(chǎng)最簡(jiǎn)單的一種離散方法。該法用隨機(jī)場(chǎng)在每個(gè)單元中心點(diǎn)的值來(lái)表征該隨機(jī)場(chǎng)在每個(gè)單元的屬性，因而隨機(jī)場(chǎng)在每個(gè)單元內(nèi)部都是常量，且等于它在各個(gè)單元中心的值。該法程序簡(jiǎn)單，但精度欠佳，現(xiàn)較少采用。 隨機(jī)場(chǎng)的局部平均 該法將隨機(jī)場(chǎng)在每個(gè)單元的屬性用隨機(jī)場(chǎng)在單元上的局部平均來(lái)表征。Baeche

18、r、Vanmarcke首先提出隨機(jī)場(chǎng)局部平均的離散方法，我國(guó)學(xué)者曾對(duì)該法進(jìn)行了深入研究2324。1. 一維隨機(jī)場(chǎng)的局部平均設(shè)一個(gè)一維連續(xù)平穩(wěn)隨機(jī)場(chǎng)S(x)，其均值為m，方差為2，則隨機(jī)場(chǎng)在一個(gè)離散單元t-T/2,t+T/2上的局部平均定義為 (7.3.9)式中：T是局部平均單元的長(zhǎng)度，ST(t)稱(chēng)為局部平均隨機(jī)場(chǎng)，其均值也為m，方差按（7.3.6）式計(jì)算，h=T ，S(t)的標(biāo)準(zhǔn)相關(guān)函數(shù)為()=R()/2，無(wú)量綱方差折減系數(shù)2（h）有如下性質(zhì)：2（T）0；2（0）=1；2（-T）=2（T）。考慮隨機(jī)場(chǎng)在兩個(gè)長(zhǎng)度分別為T(mén)和T的單元上的局部平均。如果局部平均隨機(jī)場(chǎng)分別為ST和ST，則其協(xié)方差為 (

19、7.3.10)2. 二維隨機(jī)場(chǎng)的局部平均設(shè)S(x1，x2)為二維連續(xù)參數(shù)連續(xù)狀態(tài)隨機(jī)場(chǎng)，Ai=L1iL2i為一矩形單元中心點(diǎn)（x1i，x2i）為中心，邊平行于坐標(biāo)軸x1與x2，且長(zhǎng)為L(zhǎng)1i與L2i的矩形面積，隨機(jī)場(chǎng)在A(yíng)i內(nèi)的局部平均定義為 (7.3.11)如果S(x1,x2)是一個(gè)均勻隨機(jī)場(chǎng)，則可用均值m、方差2及歸一化協(xié)方差函數(shù)（r1,r2）近似描述，r1和r2分別為沿兩個(gè)方向的距離。對(duì)應(yīng)的局部平均隨機(jī)場(chǎng)可用E(Si)、Var(Si)及互協(xié)方差Cov(Si,Sj)近似描述。E(Yi)、Var(Yi)及Cov(Yi,Yj)可由m、2及(r1，r2)計(jì)算獲得。第i和第j單元局部平均Si、Sj的互

20、協(xié)方差可表示為 (7.3.12)式中：的約定如圖7.3.1所示。 二維局部平均單元 如果有限元的網(wǎng)格已劃分，且單元總數(shù)為n，隨機(jī)場(chǎng)實(shí)際上被離散成n個(gè)隨機(jī)變量，這n個(gè)隨機(jī)變量的統(tǒng)計(jì)特性可由ESi、VarSi及互協(xié)方差CovSi，Sj反映。由于有限元網(wǎng)格的疏密是由應(yīng)力梯度決定的，與隨機(jī)場(chǎng)無(wú)關(guān)，當(dāng)單元數(shù)很多時(shí)，隨機(jī)場(chǎng)可另劃網(wǎng)格，網(wǎng)格的疏密可由相關(guān)距離S決定。當(dāng)然，隨機(jī)場(chǎng)網(wǎng)格越密精度越高。隨機(jī)場(chǎng)的局部平均法由于對(duì)原始數(shù)據(jù)的要求低、收斂快、精度高，是隨機(jī)有限元計(jì)算中最常用的方法。 隨機(jī)場(chǎng)的插值Liu W. K.提出隨機(jī)場(chǎng)的插值法。該法將隨機(jī)場(chǎng)在單元內(nèi)的值用單元結(jié)點(diǎn)處值的插值函數(shù)來(lái)表示，于是隨機(jī)場(chǎng)的統(tǒng)計(jì)特

21、性可由各單元結(jié)點(diǎn)處隨機(jī)變量間的統(tǒng)計(jì)特性近似反映。利用形函數(shù)Ni(X)，隨機(jī)場(chǎng)b(X)離散式表示為 (7.3.13)式中：X表示空間位置；bi為隨機(jī)場(chǎng)在結(jié)點(diǎn)i處的值（i=1，2，q）；q為單元結(jié)點(diǎn)數(shù)。隨機(jī)場(chǎng)b(X)在單元內(nèi)的均值和方差可表示為 (7.3.14) (7.3.15) 隨機(jī)場(chǎng)的插值法將原連續(xù)狀態(tài)的隨機(jī)場(chǎng)仍離散成一個(gè)連續(xù)函數(shù)，未直接計(jì)算隨機(jī)場(chǎng)引起的單元間的相關(guān)性，只需給定隨機(jī)場(chǎng)在各結(jié)點(diǎn)上的值，計(jì)算相對(duì)簡(jiǎn)單，易于考慮非線(xiàn)性和非均勻隨機(jī)場(chǎng)問(wèn)題。但需要已知相關(guān)函數(shù)，并且要求隨機(jī)場(chǎng)對(duì)空間參數(shù)具有較高的連續(xù)性。 隨機(jī)場(chǎng)的加權(quán)積分法Takada、Shinozuka及Deodatis提出隨機(jī)場(chǎng)加權(quán)積分

22、方法。該法在單元?jiǎng)偠染仃嚨耐茖?dǎo)過(guò)程中采用隨機(jī)場(chǎng)在單元高斯點(diǎn)上的加權(quán)積分，以表征單元上的隨機(jī)場(chǎng)。假設(shè)單元e的彈性模量為 (7.3.16)式中，E0(e)(X)為彈性模量的均值；f(e)(X)為一維零均值均勻隨機(jī)場(chǎng)，其值域?yàn)?-1+f(e)(X)-1- (01)兩端鉸接桿單元的剛度矩陣可近似表示為 K(e)K0(e)+X(e)K0(e) (7.3.17)式中，K0(e)為彈性模量取均值時(shí)的單元?jiǎng)偠染仃?，?（L為桿長(zhǎng)） (7.3.18)固結(jié)桿單元的剛度矩陣可表示為 K(e)K0(e)+X0(e)K0(e)+X1(e)K1(e)+X2(e)K2(e) (7.3.19)其中： 式（7.3.17）實(shí)質(zhì)上

23、是在考慮彈性模量隨機(jī)場(chǎng)的情況下，將單元?jiǎng)偠染仃嚪纸獬纱_定性部分和隨機(jī)部分。式（7.3.18）等表征隨機(jī)場(chǎng)在各單元的均值和方差，而單元間彈性模量的相關(guān)性則由下式表示 (7.3.20) 不難看出，局部平均法是加權(quán)積分法的特例（即權(quán)系數(shù)全部相同）。由于采用加權(quán)積分，其計(jì)算精度相對(duì)較高，而且該法積分只需一次進(jìn)行，剛度矩陣的波動(dòng)性也由此得出，因此計(jì)算效率也較高。 隨機(jī)場(chǎng)的正交展開(kāi)Spanos和Ghanem提出的隨機(jī)場(chǎng)正交展開(kāi)法，將材料特性參數(shù)隨機(jī)場(chǎng)進(jìn)行Karhumen-Loeve正交展開(kāi)，并由此推導(dǎo)出剛度矩陣的級(jí)數(shù)展開(kāi)式，從而獲得位移、應(yīng)力的統(tǒng)計(jì)特性。設(shè)隨機(jī)場(chǎng)為 (7.3.21)其中 式中：n、n(x)

24、分別為隨機(jī)場(chǎng)S(x)相關(guān)函數(shù)的特征值和特征函數(shù)。n(x)具有如下正交性： （Kronecker函數(shù)） (7.3.22)式（7.3.21）對(duì)于任何分布的S(x)均收斂，可取至第r階以滿(mǎn)足精度要求。對(duì)于一維桿單元，單元?jiǎng)偠染仃嚳杀硎緸?其中 式中：Pe(x)=De(x)/S(x), De(x)為單元彈性矩陣。集成整體剛度矩陣，得 (7.3.23)從上式，可得位移： (7.3.24)利用Neumann展開(kāi)式，可進(jìn)一步得出位移的統(tǒng)計(jì)特性。該法關(guān)鍵在于獲得特征值和特征函數(shù)。§7.4 隨機(jī)有限元與結(jié)構(gòu)可靠度 結(jié)構(gòu)可靠性分析的一次二階矩法結(jié)構(gòu)的安全性、適用性和耐久性統(tǒng)稱(chēng)為結(jié)構(gòu)的可靠性?？煽啃缘臄?shù)學(xué)

25、量度為可靠度，其定義為：在規(guī)定的條件下和規(guī)定的時(shí)間內(nèi)完成預(yù)定功能的概率。“規(guī)定條件”是指正常設(shè)計(jì)、正常施工、正常使用的條件；“規(guī)定時(shí)間”是指結(jié)構(gòu)的設(shè)計(jì)基準(zhǔn)期。結(jié)構(gòu)的使用時(shí)間超過(guò)基準(zhǔn)期后，其失效的概率將增大。結(jié)構(gòu)的一系列基本變量都具有不確定性，因此結(jié)構(gòu)的可靠性分析屬于隨機(jī)性分析的范疇，隨機(jī)有限元法將是十分有效的工具。在結(jié)構(gòu)可靠性分析中，結(jié)構(gòu)的極限狀態(tài)（包括承載能力極限狀態(tài)、正常使用極限狀態(tài)和條件極限狀態(tài)）是通過(guò)功能函數(shù)來(lái)描述的。當(dāng)有n個(gè)隨機(jī)變量影響結(jié)構(gòu)可靠度時(shí)，結(jié)構(gòu)的功能函數(shù)為: z=g(X1,X2,Xn)，基本變量Xi(i=1,2,n)是結(jié)構(gòu)上的各種外因作用、材料性能和幾何參數(shù)等。z0時(shí)，結(jié)構(gòu)

26、處于可靠狀態(tài)；z0,則處于失效狀態(tài)；z=0稱(chēng)結(jié)構(gòu)極限狀態(tài)方程（一般難以用顯式表示），結(jié)構(gòu)處于極限狀態(tài)。如果z的概率密度函數(shù)或概率分布函數(shù)可求得，則結(jié)構(gòu)可靠度的數(shù)量指標(biāo)便可基于各種狀態(tài)出現(xiàn)的概率而確定。若功能函數(shù)僅與兩個(gè)隨機(jī)變量有關(guān)（如結(jié)構(gòu)抵抗破壞或變形的能力R和荷載引起結(jié)構(gòu)內(nèi)力、應(yīng)力、位移等效應(yīng)S），即z=g(R,S)。假設(shè)R、S均為正態(tài)分布，其均值和標(biāo)準(zhǔn)差分別為和，此時(shí)z=R-S也是正態(tài)分布的隨機(jī)變量，具有均值和標(biāo)準(zhǔn)差。其概率密度函數(shù)為 (7.4.1)其分布圖示于圖7.4.1，陰影部分是結(jié)構(gòu)失效概率Pf，非陰影部分面積即結(jié)構(gòu)的可靠度Pr。 圖7.4.1 正態(tài)分布概率密度函數(shù) 圖7.4.2 失

27、效邊界在工程實(shí)踐中R和S不一定是正態(tài)分布，但可以變換成標(biāo)準(zhǔn)正態(tài)分布。他們的統(tǒng)計(jì)特征量是均值、標(biāo)準(zhǔn)差、相關(guān)偏度或變異系數(shù)等，變異系數(shù)=/，表示隨機(jī)變量相對(duì)于均值的變異。目前工程上一般用無(wú)量綱的可靠指標(biāo)來(lái)反映結(jié)構(gòu)的可靠度，越大，失效概率Pf越小，其互補(bǔ)的可靠度Pr就越大。一次二階矩法采用只需已知均值和標(biāo)準(zhǔn)差的數(shù)學(xué)模型去求解結(jié)構(gòu)的可靠度。此法將功能函數(shù)z=g(X1,X2,Xn)在某點(diǎn)用Taylor級(jí)數(shù)展開(kāi)，并近似地取一次項(xiàng)使極限狀態(tài)方程線(xiàn)性化，然后求得可靠指標(biāo)。改進(jìn)的一次二階矩法將線(xiàn)性化點(diǎn)選在失效邊界上，而且選在結(jié)構(gòu)最大可能失效點(diǎn)P*上（如圖7.4.2）。選擇設(shè)計(jì)驗(yàn)算點(diǎn)P*(Xi*i=1,2,n)作

28、為線(xiàn)性化點(diǎn)X0i時(shí)，線(xiàn)性化的極限狀態(tài)方程為 (7.4.2)則z的均值為 由于設(shè)計(jì)驗(yàn)算點(diǎn)P*選在失效邊界上，有g(shù)(X1*,X2*,Xn*)=0，因此z成為 但隨機(jī)變量Xi互不相關(guān)時(shí)，z的標(biāo)準(zhǔn)差z為 其中 稱(chēng)i為靈敏度系數(shù)，表示第i個(gè)隨機(jī)變量對(duì)標(biāo)準(zhǔn)差的相對(duì)影響。于是可靠指標(biāo)為 (7.4.3)變換上式為 即 (i=1，2，n) (7.4.4)上式中，Xi，Xi為已知的各隨機(jī)變量的均值和標(biāo)準(zhǔn)差，待求的量為Xi*和，可迭代求解。 隨機(jī)有限元的一次二階矩法設(shè)可靠性分析中的一組基本變量X=(X1,X2,Xn)T為相互獨(dú)立的正態(tài)變量（不滿(mǎn)足時(shí)可作變換），并進(jìn)一步變換為標(biāo)準(zhǔn)正態(tài)變量Y=(Y1,Y2,Yn)T，其

29、中Yi=(Xi-Xi)/Xi。于是功能函數(shù)也可轉(zhuǎn)換到標(biāo)準(zhǔn)正態(tài)空間，即 g(X)=g(R(X),S(X)=g(R(T-1(Y),S(T-1(Y)=G(Y)失效概率Pf=1-(),由的幾何含義可知，值為在標(biāo)準(zhǔn)正態(tài)空間中從原點(diǎn)到失效面的最短距離。在標(biāo)準(zhǔn)正態(tài)空間中概率密度是關(guān)于原點(diǎn)（均值點(diǎn)）旋轉(zhuǎn)對(duì)稱(chēng)，并且隨著到原點(diǎn)的距離的平方呈指數(shù)下降。在一次二階矩法中，標(biāo)準(zhǔn)正態(tài)空間中的極限狀態(tài)面（失效面）被一個(gè)到原點(diǎn)最小距離點(diǎn)處的切平面代替，也即按一階Taylor展開(kāi)式把功能函數(shù)在設(shè)計(jì)點(diǎn)處展開(kāi)，使之線(xiàn)性化。采用迭代法可以近似確定極限狀態(tài)面G(Y)=0上距原點(diǎn)最近的點(diǎn)Y*,然后按距離公式確定結(jié)構(gòu)的可靠指標(biāo)。具體的迭代

30、格式如下： (7.4.5) (7.4.6)以Mises屈服準(zhǔn)則為例，功能函數(shù) g(S)=s02-(s2xx-sxxsyy+s2yy+3s2xy)于是極限狀態(tài)方程 g(S(X(Y)=G(Y)=0 (7.4.7)計(jì)算可靠指標(biāo)時(shí)要用到功能函數(shù)G(Y)的梯度向量，由求導(dǎo)鏈?zhǔn)椒▌t (7.4.8)由于S(X)難于用顯式表達(dá)，求存在困難。有效方法是利用隨機(jī)有限元一次二階矩法計(jì)算結(jié)構(gòu)的可靠指標(biāo)和設(shè)計(jì)點(diǎn)Y*。步驟如下：（1） 確定隨機(jī)變量轉(zhuǎn)換關(guān)系Y=Y(X)；（2） 給定初值Y(0)=0，計(jì)算式（7.4.8）右邊的三個(gè)偏導(dǎo)數(shù)；（3） 按式（7.4.5）和（7.4.6）計(jì)算設(shè)計(jì)驗(yàn)算點(diǎn)坐標(biāo)Y*；（4） 用隨機(jī)有限元

31、法計(jì)算S(i)和；（5） 重復(fù)（3）、（4）兩步驟，直至收斂（G(Y(i)=0）。（6） 再按式（7.4.5）和（7.4.6）計(jì)算可靠指標(biāo)值。 隨機(jī)有限元的最大熵法一次二階矩法的計(jì)算量隨迭代次數(shù)成倍增加，使該法的使用受到限制。而最大熵法用于結(jié)構(gòu)的可靠性分析時(shí)，可根據(jù)隨機(jī)變量的二階矩來(lái)擬合概率密度曲線(xiàn)。因此隨機(jī)有限元結(jié)合最大熵法可用于求結(jié)構(gòu)響應(yīng)的概率密度曲線(xiàn)，從而計(jì)算結(jié)構(gòu)的失效概率。熵被定義為信息的均值，信息是對(duì)個(gè)別X值不確定性的度量。不確定性越大，熵也越大。對(duì)于一個(gè)連續(xù)隨機(jī)變量，熵為 ，而對(duì)于離散隨機(jī)變量，熵為，這里，f(X)是隨機(jī)變量X的概率密度函數(shù)；f(xi)是離散點(diǎn)概率函數(shù)。最大熵法通過(guò)

32、調(diào)整概率密度函數(shù)f(X)使熵S取得最大值，并滿(mǎn)足約束條件： ，利用最大熵法可得近似的概率密度函數(shù)f(X)。為了求得最大值，可引入拉格朗日乘子法，把問(wèn)題轉(zhuǎn)化為確定拉格朗日乘子。選用基于Neumann級(jí)數(shù)的隨機(jī)有限元法，寫(xiě)出遞推方程?？蓪⒕植科骄碚撘隢eumann隨機(jī)有限元法，以進(jìn)一步提高計(jì)算效率。隨機(jī)有限元最大熵法計(jì)算結(jié)構(gòu)失效概率的步驟如下：(1) 隨機(jī)場(chǎng)離散，隨機(jī)變量化為隨機(jī)向量，并計(jì)算協(xié)方差矩陣；(2) 利用特征正交化技術(shù)，用一組統(tǒng)計(jì)獨(dú)立的隨機(jī)變量描述隨機(jī)場(chǎng)；(3) 利用Neumann隨機(jī)有限元法計(jì)算響應(yīng)量的前若干階矩；(4) 擬合所有節(jié)點(diǎn)或部分節(jié)點(diǎn)的最大熵密度函數(shù)；(5) 利用數(shù)值積分

33、計(jì)算結(jié)構(gòu)的失效概率。 算例例1 如圖7.4.3所示桁架下端受一集中力，假設(shè)各設(shè)計(jì)變量都是均勻分布的隨機(jī)變量，三根桿彈性模量E，E=2.0×106N/cm2，Cov(E)=2.5%；截面積A，A=0.8cm2，Cov(A)=2%；桿2長(zhǎng)度L2,L2=20cm，Cov(L2)=1.5%；桿1和桿3的長(zhǎng)度均為L(zhǎng)1，L1=28.284cm，Cov(L1)=2.3%；外力P，P=4000N；Cov(P)=10%；材料屈服應(yīng)力y，y=3200N/cm2，Cov(y)=7%。當(dāng)極限狀態(tài)方程取為 y-S=0（S為構(gòu)件中最大應(yīng)力值）時(shí)，試求該結(jié)構(gòu)的失效概率。圖7.4.3 簡(jiǎn)單超靜定桁架解： 用隨機(jī)有限

34、元一次二階矩法計(jì)算，假設(shè)由各結(jié)點(diǎn)坐標(biāo)值變異引起的方向余弦的變異很小，形成確定矩陣，則由隨機(jī)有限元遞推方程組可得各結(jié)點(diǎn)位移分布特征值。計(jì)算得結(jié)點(diǎn)4的位移均值和標(biāo)準(zhǔn)差為4=2.2928cm和4=0.0037cm；各桿件中最大應(yīng)力為桿元2的應(yīng)力，其均值和標(biāo)準(zhǔn)差為2=2928N/cm2和2=377N/cm2。由一次二階矩法可計(jì)算得該結(jié)構(gòu)的失效概率為0.23（由于結(jié)構(gòu)超靜定，桿2失效并不意味著整個(gè)結(jié)構(gòu)失效），與蒙特卡羅有限元法計(jì)算結(jié)果比較接近（蒙特卡羅取樣數(shù)為40時(shí)，失效概率達(dá)0.224）。隨機(jī)有限元法在分析混凝土重力壩可靠度方面已得到較多應(yīng)用。研究者采用抽樣方法對(duì)重力壩施行隨機(jī)自動(dòng)剖分，研究其可靠指標(biāo)

35、的統(tǒng)計(jì)規(guī)律和收斂于真解的規(guī)律，并進(jìn)行結(jié)構(gòu)優(yōu)化25。在對(duì)碾壓混凝土重力壩抗滑穩(wěn)定的體系可靠度研究中，還采用剛體極限平衡法、線(xiàn)彈性隨機(jī)有限元法和三維彈塑性隨機(jī)有限元法進(jìn)行計(jì)算比較。§7.5 隨機(jī)有限元?jiǎng)恿Ψ治龇椒?不確定結(jié)構(gòu)的自振特性假設(shè)某一梁-柱結(jié)構(gòu)承受隨機(jī)分布的軸向靜荷載N(x)，沿長(zhǎng)度有微小擾動(dòng)，同時(shí)在端部有軸向推力c為隨機(jī)變量，其方差為，軸向分布荷載可表示為式中p(x)是一維、均值為零的均勻隨機(jī)場(chǎng)，方差是P2，相關(guān)偏度是P，互相關(guān)函數(shù)是RPP()。結(jié)構(gòu)彈性模量和質(zhì)量密度也隨機(jī)變化，分別表示為，其中分別是彈性模量E和質(zhì)量密度m的均值，方差分別為2E和2m；隨機(jī)變量a(x)和b(x)

36、是兩個(gè)一維互不相關(guān)的、均勻的、均值為零的隨機(jī)場(chǎng)。自相關(guān)函數(shù)分別為Raa()和Rbb()或等價(jià)的功率譜密度函數(shù)為Saa()和Sbb()，相關(guān)偏度分別是E和m。任一點(diǎn)橫向位移的有限元模式為 w(x,t)=Neqe，其中，單元位移列陣qe=w1,1,w2,2T,形函數(shù)Ne取三次多項(xiàng)式插值。單元的動(dòng)能表達(dá)為： (7.5.1)式中：，Ae為單元橫截面面積。單元的應(yīng)變能為 (7.5.2)式中：，I是慣性矩。單元的外力功為 (7.5.3)式中：，p為任一隨機(jī)變量據(jù)此集成結(jié)構(gòu)總能量和總外力功，利用Hamilton原理，即對(duì)總位移向量q取變分，考慮到q的任意性，可得系統(tǒng)運(yùn)動(dòng)方程，進(jìn)而得到如下特征問(wèn)題方程： ）式

37、中，剛度系數(shù)，其中是確定性分量，而隨機(jī)擾動(dòng)項(xiàng)：。于是可以寫(xiě)出剛度系數(shù)、質(zhì)量系數(shù)以及幾何剛度系數(shù)的均值和方差的表達(dá)式26。在總剛度矩陣中兩個(gè)剛度系數(shù)kij和krs的互協(xié)方差，當(dāng)單元尺寸相等時(shí)，利用方差函數(shù)可作進(jìn)一步簡(jiǎn)化。利用局部平均理論，在每個(gè)單元上E、m、Q的均值為零，方差為： (7.5.5)式中，和分別為隨機(jī)變量a(x)、b(x)和p(x)的方差函數(shù)和相關(guān)偏度。利用協(xié)方差函數(shù)寫(xiě)出兩剛度系數(shù)或兩質(zhì)量系數(shù)的互協(xié)方差。最后，特征值問(wèn)題的方程可表達(dá)為 (7.5.6)或 K*x=M*x 。相應(yīng)平均問(wèn)題的方程為 (7.5.7) 利用單元?jiǎng)偠?、質(zhì)量和幾何剛度之間的協(xié)方差的表達(dá)式，等效剛度矩陣K*的協(xié)方差矩陣能用獨(dú)立的a(x)、b(x)和p(x)來(lái)表示。解未擾動(dòng)的特征值問(wèn)題（平均值問(wèn)題）可得特征值的均值。對(duì)于滿(mǎn)足小擾動(dòng)情況，可推導(dǎo)出求兩個(gè)