隨機有限元法

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

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

3、分為兩大類：一類是統(tǒng)計方法，另一類是非統(tǒng)計方法。因 此，隨機有限元法同樣也有統(tǒng)計逼近和非統(tǒng)計逼近兩種類型。前者通過樣本試驗收集原始的數(shù)據(jù)資料， 運用概率和統(tǒng)計理論進行分析和整理， 然后作出科學推斷。 這里，樣本試驗和數(shù) 據(jù)處理的工作量很大，隨著計算機的普及和開展，數(shù)值模擬法，如蒙特卡羅Monte Carlo 模擬， 已成為最常用的統(tǒng)計逼近法。 后者從本質(zhì)上來說是利用分析工具找出結(jié)構(gòu)系統(tǒng)的 確 定的或隨機的 輸出隨機信號與輸入隨機信號之間的關(guān)系， 采用隨機分析與求解系統(tǒng)控制方 程相結(jié)合的方法得到輸出信號的各階隨機統(tǒng)計量的數(shù)字特征如各階原點矩或中心矩。在20世紀70年代初，Cambout先采用一次

4、二階矩方法研究線彈性問題。由于這種方法將隨機變量的影響量進行 Taylor 級數(shù)展開，就稱之為 Taylor 展開法隨機有限元 TSFEM。 Shinozuka 和 Astill1972 分別獨立運用攝動技術(shù)研究了隨機系統(tǒng)的特征值問題。隨后，Handa 1975等人在考慮隨機變量波動性時采用一階和二階攝動技術(shù)，并將這種攝動法隨 機有限元成功地應(yīng)用于框架結(jié)構(gòu)分析。Vanmarcke 等人 1983提出隨機場的局部平均理論，并將它引入隨機有限元。 局部平均理論是用隨機場函數(shù)在每一個離散單元上的局部平均的隨 機變量來代表該單元的統(tǒng)計量的近似理論。Liu W. K.等人1986、1988的系列工作，提

5、供了一種“主模態(tài) 技術(shù)， 運用隨機變量的特征正交化方法，將滿秩的協(xié)方差矩陣變換為對 角矩陣，減少計算工作量，對攝動隨機有限元法的開展做出奉獻，此外，提出了一個隨機變分原理。Yamazaki和Shi no zuka1987創(chuàng)造性地將算子的 Neuma nn級數(shù)展開式引入隨機有限元 的列式工作。從本質(zhì)上講，Neuma nn級數(shù)展開方法也是一類正那么的小參數(shù)攝動方法，正定的隨機剛度矩陣和微小的隨機擾動量是兩個根本要求， 這兩個根本要求保證了攝動解的正那么性 和收斂性，其優(yōu)點在于攝動形式較簡單并可以得到近似解的高階統(tǒng)計量。Shinozuka 等人1987將隨機場函數(shù)的 Monte-Carlo模擬與隨機

6、剛度矩陣的Neumann級數(shù)展開式結(jié)合,得到具有較好計算精度和效率的一類Neumann隨機有限元列式稱NSFEMoBenaroya等1988指出，將出現(xiàn)以隨機變分原理為根底的隨機有限元法來逐漸取代以攝動法為根底的隨機有限 元法。Spa nos和Gha nem等人1989,1991結(jié)合隨機場函數(shù)的 Karhue n-Loeve展式和 Galerki n 迦遼金 射影方法建立了相應(yīng)的隨機有限元列式， 并撰寫了隨機有限元法領(lǐng)域的第一本專 著?隨機有限元譜方法? 。國內(nèi)對隨機有限元的研究起步較晚。吳世偉等人 1988提出隨機有限元的直接偏微分法及相應(yīng)的可靠度計算方法。陳虬、劉先斌等人(1989、199

7、1)提出一種新的隨機場離散模型，建立了等參局部平均單元，并基于變分原理研究了一類隨機有限元法的收斂性和誤差 界。Papadrakakis(1995)采用預(yù)處理共軛梯度法給出了空間框架的非線性隨機有限元列式。Schorling 和Bucher(1996)基于Monte-Carlo 技術(shù)，采用響應(yīng)面法研究幾何非線性時的 可靠度隨機有限元方法。 劉寧(1996)那么基于偏微分法，給出了三維彈塑性隨機有限元列式。 隨機有限元法的數(shù)學理論研究和非線性隨機問題的有限元分析工作還有待深入。自20世紀80年代以來，隨機有限元法已在工程結(jié)構(gòu)可靠性、平安性分析領(lǐng)域以及在 各種隨機鼓勵下結(jié)構(gòu)響應(yīng)變異研究領(lǐng)域中得到應(yīng)

8、用，如應(yīng)用于大型水利工程的重力壩、拱壩的可靠度計算；應(yīng)用于非線性瞬態(tài)響應(yīng)分析；結(jié)構(gòu)振動中隨機阻尼對響應(yīng)的影響；結(jié)構(gòu)分析的隨機識別；復雜結(jié)構(gòu)地震響應(yīng)的隨機分析和兩相動力系統(tǒng)的隨機模擬等等。隨著理論研究的深入，隨機有限元將得到更加廣泛的應(yīng)用。 7.2隨機有限元的控制方程22從隨機有限元控制方程的獲得來看，隨機有限元可分為Taylor展開法隨機有限元(TSFEM、攝動法隨機有限元(PSFEM)以及Neumann展開Monte-Carlo法隨機有限元(NSFEM。 Taylor展開法隨機有限元該隨機有限元法的根本思路是將有限元格式中的控制量在隨機變量均值點處進行Taylor級數(shù)展開(取一階或二階)，經(jīng)

9、過適當?shù)臄?shù)學處理得出所需的計算方程式。有限元靜力分析控制方程的矩陣形式為：KU = F(7.2.1)式中，U為位移矩陣，F(xiàn)為等效節(jié)點荷載列陣，K為整體剛度矩陣K 八.BTDBdv(7.2.2)其中，B為形變矩陣，D為材料彈性矩陣。在計算出節(jié)點位移 U后，即由下式求得應(yīng)力列陣bd = DBU(7.2.3)設(shè)根本隨機變量為 X =(X1,X2，,Xn)T，將位移U在均值點X =(X1,X2 , Xn)T處一階Taylor級數(shù)展開，并在兩邊同時取均值(數(shù)學期望)，得(7.2.4)式中：符號E 表示求均值，任一結(jié)點位移U的方差可由下式計算:n n z-j -VarUJ : Ui =1 j 1 Xi：U

10、X 文：X jCov(Xi , Xj)xM(7.2.5)式中：符號Var 表示求方差；Cov(Xi,Xj)為X和X的協(xié)方差。其中.:Uk4(u)Xi cXi(7.2.6)同樣將d BU + DB空X:Xi:Xi(7.2.7)d在均值點處Taylor展開，也有與上面類似的表達式。TSFEM關(guān)鍵在于故該法也稱可見，對有限元方程式直接進行偏微分計算，計算出有限元輸出量對隨機變量的梯度，直接偏微分法或梯度分析法。由于一階TSFEM只需一次形成剛度矩陣，也只需一次求剛度矩陣的逆，因此效率較高。但由于忽略了二階以上的高次項，使 TSFEM對隨機變量的變異性有所限制。一般要求一階 TSFEM隨機變量的變異系

11、數(shù)小于 0.3。如果隨機變量的變異系數(shù)較大，可以采用有限元控制 方程的二階Taylor展開：2 / 2cU d F=K (X i iX ji CX ptX j:K.Xi ；Xj2:K?U -K U-X i -X jX i .X j(728)Xi ：XjdXi ?XjBU；:D+：Xj.XiDB；：2UrAZ lZXplXj(7.2.9)上式可見，二階TSFEM可以放寬隨機變量變異性大小的限制，但隨機變量數(shù)目較多時， 計算量將十分龐大，而且一階或二階TSFEM均無法計算響應(yīng)量三階以上的統(tǒng)計特性。由于TSFEM簡單明了、效率高，為我國許多學者所采用。攝動法隨機有限元攝動技術(shù)最初被用于非線性力學分析

12、。Handa等人成功地將一階、二階攝動技術(shù)用于隨機問題，給出攝動法有限元列式。該法假定根本隨機變量在均值點處產(chǎn)生微小攝動，利用 Taylor級數(shù)把隨機變量表示為確定局部和由攝動引起的隨機局部，從而將有限元控制方程(非線性的)轉(zhuǎn)化為一組線性的遞推方程，求解得出位移的統(tǒng)計特性，進而求出應(yīng)力的統(tǒng)計特性。假設(shè):i為隨機變量Xi在均值點Xi處的微小攝動量，即r =Xi - Xi。于是n:K 1 n；：2kK。；二：r j? i2 i,j 二：? i j(7.2.10)對于U F,也有類似上式 K的表達式，式中：K0、U0、F0分別為K、U F在隨機變量均值點 的值。根據(jù)二階攝動法，可得(7.2.11)K

13、o1U。二 K。F。(7.2.12)蘭U J cF , cKcU cKcK cU(7.2.13)=K。_U。_網(wǎng)/j 宮刊必/j丿由上式可得位移的均值和協(xié)方差:1 n n(7.2.14)eU h UUijCovC i / j)n nn n nCovU：二二 UiUjCovGi, ： j)亠二二二 UiUjkEGi： j k)i d j di =i j =i k=i(7.2.15)n n n nUijUkiE(：i：j)E(*i) E(：i：k)E(ji)li=1 j 1 k=1 I 由于任何量的隨機性都可以引入攝動量，而且更易于考慮非線性問題，因此PSEFEM適用范圍較廣，對于結(jié)構(gòu)幾何特性的隨

14、機性(包括隨機邊界問題)易得出隨機有限元控制方程。 一階PSFEM和一階TSFEMH樣，只需一次形成剛度矩陣、一次對剛度矩陣求逆，計算效率較 高。但PSFEM需以微小的攝動量為條件，一般應(yīng)小于均值的20%或 30% Neumann展開Monte-Carlo 隨機有限元20世紀80年代后期，Shinozuka等人提出基于 Neumann展開式的隨機有限元法，使 Mon te-Carlo法與有限元法較完美地結(jié)合起來。Mon te-Carlo法是最直觀、最精確、獲取信息最多、對非線性問題最有效的計算統(tǒng)計方法。Neuma nn展開式的引入是為了解決矩陣求逆的效率問題。如果對每一次隨機抽樣，只需形成剛度

15、矩陣，進行前代、回代以及矩陣乘和矩陣加減，而無需矩陣分解，那么可大大減少工作量。在一般有限元控制方程KU = F中，假定荷載F為確定值，在隨機變量波動值的影響下剛度矩陣K分解為K = K 0+A K,根據(jù)Neumann級數(shù)展開，有-1- 12 3- 1K =( Kg+ A K) =(I-P+P-P+, ) Kg(7.2.16)式中：Kg為隨機變量均值處的剛度矩陣；A K為剛度矩陣的波動量；I為單位矩陣。對于Monte-Carlo隨機抽樣，剛度矩陣只改變A K項，而P = Ko-1 A K(7.2.17)U0 =K-1F(7.2.18)將式(7.2.16 )代入式(7.2.1 ),并利用式(7.

16、2.17 :)和式(7.2.18 ),得23U = Uo- PUo+P U0- P Uo+,(7.2.19)令U=PU0,那么得如下的遞推公式：U = U0 + U1+U2+,(7.2.20)1Ui =- Ko- A KU1(i=1,2,3,)(7.2.21)由式(7.2.18 )求出 U0后，可由式(7.2.21 )求出 U(i=1,2,3,)。上述三種方法中，NSFEM可以方便地調(diào)用確定性有限元的計算程序，而 TSFEM在編程上 較為復雜，PSFEM那么更為復雜。由于采用 Mon te-Carlo隨機模擬技術(shù)，NSFEM不受隨機變量 波動范圍的限制，當變異系數(shù)小于 0.2時，NSFEM與一

17、階TSFEM或一階PSFEMW度相當；當 變異系數(shù)大于0.2時，后兩者已不能滿足精度要求，但 NSFEM仍能得出滿意的結(jié)果。 7.3隨機場的離散模型許多物理現(xiàn)象和物體系統(tǒng)具有隨機分布特性，包括系統(tǒng)本身的不確定或系統(tǒng)的鼓勵和響應(yīng)的不確定，都可以模型化為隨機空間分布的隨機場或隨時間分布的隨機過程。隨機有限元法除了必須考慮材料參數(shù)等的空間變異性，需要獲得隨機有限元方程列式以及解決隨機算子和隨機矩陣的求逆問題外，還須包含對隨機場的離散處理。均勻各向同性隨機場的特征量1.隨機場S(t)的均值E(S(t)為常數(shù)m.2.隨機場S(t)的標準相關(guān)函數(shù):：()R()R(0)(7.3.1)式中，隨機場的協(xié)方差函數(shù)

18、R( T)=CovS(t+t ),S(t)(7.3.2)對于一切t，隨機場的方差為VarS(t)= R(0)=/(7.3.3)相關(guān)函數(shù)也可以用譜分解表示(即Wie ner-Khi ntchi ne變換對)：1 -beG(G0)=f RG)coa dT2兀 q(7.3.4)-boRO) = LG)coa di常用的標準相關(guān)函數(shù)有：非協(xié)調(diào)階躍型、協(xié)調(diào)階躍型、三角型、指數(shù)型、二階AR型、高斯型等多種形式。3.方差折減系數(shù)r 2(h)設(shè)Sh(t)是隨機場S(t)的局部平均隨機過程，即114hSh(t- t S(t)dth 4(7.3.5)那么方差VarSh(t)= d 2r 2(h)(7.3.6)式中

19、，方差折減系數(shù)22 h-(h)0(1)( )dhh(7.3.7)可見，方差折減系數(shù)起著使Sh(t)的方差比原來S(t)的方差縮小的作用。4.相關(guān)距離S SS S可以看成是任意兩個相隔距離為T的隨機變量不相關(guān)的最小距離(也稱相關(guān)偏度)t| 6s，那么不相關(guān)，否那么完全相關(guān)。利用相關(guān)距離S S便于對隨機場作近似處理，其計算式為：hV2 (h) , h t、s =三2(7.3.8)2G( ) / 二，二0以上公式均對一維問題列出，二維、三維問題也可以類似得出。 隨機場的中心離散中心離散是隨機場最簡單的一種離散方法。該法用隨機場在每個單元中心點的值來表征該隨機場在每個單元的屬性，因而隨機場在每個單元內(nèi)

20、部都是常量，且等于它在各個單元中心的值。該法程序簡單，但精度欠佳，現(xiàn)較少采用。隨機場的局部平均Baecher、2324。該法將隨機場在每個單元的屬性用隨機場在單元上的局部平均來表征。Vanmarcke首先提出隨機場局部平均的離散方法，我國學者曾對該法進行了深入研究1. 一維隨機場的局部平均設(shè)一個一維連續(xù)平穩(wěn)隨機場S(x)，其均值為 m,方差為d 2,那么隨機場在一個離散單元t-T/2,t+T/2上的局部平均定義為ST(t)t (T / 2)T t_(T/2)S(x)dx(7.3.9)式中：T是局部平均單元的長度，ST(t)稱為局部平均隨機場，其均值也為 m方差按(7.3.6 )2 2式計算，h

21、=T , S(t)的標準相關(guān)函數(shù)為P ( T )=R( T )/ d ,無量綱方差折減系數(shù)r (h)有如下性質(zhì)： r2(T) o； r2(0)=1 : r2(-t)=r2(T)?？紤]隨機場在兩個長度分別為T和T的單元上的局部平均。如果局部平均隨機場分別為St和St那么其協(xié)方差為COV0 , St )二2TT3八(-1)kTj2(Tk)k =0(7.3.10)2.二維隨機場的局部平均設(shè)S(x1, X2)為二維連續(xù)參數(shù)連續(xù)狀態(tài)隨機場，A=L1iL2i為一矩形單元中心點(X1, X2 )為中心，邊平行于坐標軸X1與X2，且長為L1i與L2i的矩形面積，隨機場在A內(nèi)的局部平均、1xii*Lii/2)

22、X2i*L2i/2)(7.3.11)定義為Si =SAi(Xii,X2i)二/2)S(Xi,X2)dxidx2A x1i _(L1i / 2) x2i (L2i 1 2)r1,r 2)如果S(X1,x 2)是一個均勻隨機場，那么可用均值 m方差b 2及歸一化協(xié)方差函數(shù)p近似描述，n和“分別為沿兩個方向的距離。對應(yīng)的局部平均隨機場可用E(Si)、Var(Si)及互協(xié)方差 Cov(Si,Sj)近似描述。E(Yi)、Var(Yi)及Cov(Yi,Yj)可由m 異及p (r 1,2)計算 獲得。第i和第j單元局部平均S、S的互協(xié)方差可表示為Cov(S , Sj)=33 (-1)5訂21)2】fl)k

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

24、值用單元結(jié)點處值的插值函 數(shù)來表示，于是隨機場的統(tǒng)計特性可由各單元結(jié)點處隨機變量間的統(tǒng)計特性近似反映。利用形函數(shù)N(X)，隨機場b(X)離散式表示為qb(X)八 Ni(X)bi(7.3.13)式中：X表示空間位置；bi為隨機場在結(jié)點i處的值(i=1，2，,， q) ; q為單元結(jié)點數(shù)。 隨機場b(X)在單元內(nèi)的均值和方差可表示為qE(b(X)Ni(X)E(bi)(7.3.14)i 二iqVar(b(X)=遲 Ni(X)Nj(X)Cov(b ,bj)(7.3.15)i,U隨機場的插值法將原連續(xù)狀態(tài)的隨機場仍離散成一個連續(xù)函數(shù)，未直接計算隨機場引起的單元間的相關(guān)性，只需給定隨機場在各結(jié)點上的值，

25、計算相對簡單，易于考慮非線性和非 均勻隨機場問題。但需要相關(guān)函數(shù)，并且要求隨機場對空間參數(shù)具有較高的連續(xù)性。隨機場的加權(quán)積分法Takada、Shi no zuka及Deodatis提出隨機場加權(quán)積分方法。該法在單元剛度矩陣的推 導過程中采用隨機場在單元高斯點上的加權(quán)積分，以表征單元上的隨機場。假設(shè)單元e的彈性模量為E(X)二E0e)f (e)(X)(7316)式中，曰(e)(x)為彈性模量的均值；f(e)(x)為一維零均值均勻隨機場，其值域為-1+n w f(e) (X) w -1- n (0n(7.3.20)不難看出，局部平均法是加權(quán)積分法的特例(即權(quán)系數(shù)全部相同)。由于采用加權(quán)積分，其計算

26、精度相對較高，而且該法積分只需一次進行，剛度矩陣的波動性也由此得出，因此計算效率也較高。隨機場的正交展開Spa nos和Gha nem提出的隨機場正交展開法，將材料特性參數(shù)隨機場進行 Karhumen-Loeve正交展開，并由此推導出剛度矩陣的級數(shù)展開式，從而獲得位移、應(yīng)力的 統(tǒng)計特性。_QO設(shè)隨機場為S(x) = S(x)亠二 bn - n ：n (x)(7.3.21)n =01其中* ： .LS(X)n(X)dX式中：入n、U n(x)分別為隨機場S(x)相關(guān)函數(shù)的特征值和特征函數(shù)。U n(x)具有如下正交性：m(x) n(x)dx =、mn ( Kronecker 函數(shù))(7.3.22)

27、式(7.3.21 )對于任何分布的S(x)均收斂，可取至第r階以滿足精度要求。對于一維桿單e -e re元，單元剛度矩陣可表示為ke = k a bnk：n =0其中k： = kn 護n(x)B(x)TPe(x)Be(x)dx式中：Pe(x)=De(x)/S(x), De(x)為單元彈性矩陣。集成整體剛度矩陣，得(7323)從上式，可得位移:U =(1、bnKKjKFn衛(wèi)(7.3.24)該法關(guān)鍵在于獲得特征值和特征利用Neumann展開式，可進一步得出位移的統(tǒng)計特性。 函數(shù)。 7.4隨機有限元與結(jié)構(gòu)可靠度結(jié)構(gòu)可靠性分析的一次二階矩法結(jié)構(gòu)的平安性、適用性和耐久性統(tǒng)稱為結(jié)構(gòu)的可靠性。可靠性的數(shù)學量

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

29、z0時，結(jié)構(gòu)處于可靠狀態(tài)；zv0,那么處于失效狀態(tài)；z=0稱結(jié)構(gòu)極限狀態(tài)方程(一般難以用顯式表示)，結(jié)構(gòu)處于極限狀態(tài)。如果 z的概率密度函數(shù)或概率分布函 數(shù)可求得，那么結(jié)構(gòu)可靠度的數(shù)量指標便可基于各種狀態(tài)出現(xiàn)的概率而確定。假設(shè)功能函數(shù)僅與兩個隨機變量有關(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)分布，其均值和標準差分別為z = R-S和標準差R、S和二r、s，此時z=R-S也是正態(tài)分布的隨機變量，具有均值二z. ；R飛S。其概率密度函數(shù)為(7.4.1)其分布圖示于圖7.4.1 ,陰影局部是結(jié)構(gòu)失效概率Pf，非陰影局部面積即結(jié)

30、構(gòu)的可靠度P。圖正態(tài)分布概率密度函數(shù)圖7.4.2失效邊界在工程實踐中R和S不一定是正態(tài)分布，但可以變換成標準正態(tài)分布。他們的統(tǒng)計特征 量是均值、標準差6、相關(guān)偏度或變異系數(shù)等，變異系數(shù)V =6 /卩，表示隨機變量相對于均值的變異。目前工程上一般用無量綱的可靠指標B來反映結(jié)構(gòu)的可靠度，一：=z/；z, 3越大，失效概率 Pf越小，其互補的可靠度Pr就越大。一次二階矩法采用只需均值和標準差的數(shù)學模型去求解結(jié)構(gòu)的可靠度。此法將功能函數(shù)z=g(Xi,X2, ,Xn)在某點用Taylor級數(shù)展開，并近似地取一次項使極限狀態(tài)方程線性 化，然后求得可靠指標3。改良的一次二階矩法將線性化點選在失效邊界上，而且

31、選在結(jié)構(gòu)最大可能失效點P*上* *(如圖742 )。選擇設(shè)計驗算點 P(Xi | i=1,2, , ,n)作為線性化點Xoi時，線性化的極限狀n態(tài)方程為z g(X；,X；,，X；廠(Xiid-Xi二 0X*(742)那么z的均值為n)z =g(x；,x；,x；)、-X；)磊由于設(shè)計驗算點P*選在失效邊界上，有g(shù)(X1* ,X 2*,X)=0 ,因此3成為但隨機變量X互不相關(guān)時，z的標準差6 z為2=LIi 土X；n=11i 土其中-i_nfh n 1z 嚴JcXiX*丿X2CXi ：g:Xi稱a i為靈敏度系數(shù)，表示第i個隨機變量對標準差的相對影響。于是可靠指標3為ZCJZ(7.4 3)變換上

32、式為“Xi -X； - l&xi) =0X*即Xi；-xi 一 - Xi (i=1,2，,，n) (7.4.4)上式中，3 Xi, 6 Xi為的各隨機變量的均值和標準差，待求的量為X*和3,可迭代求解。隨機有限元的一次二階矩法設(shè)可靠性分析中的一組根本變量X=(X1,X2, , ,Xn)T為相互獨立的正態(tài)變量(不滿足時可作變換)，并進一步變換為標準正態(tài)變量丫=(丫1,丫2, , ,Yn)T,其中Yi=(Xi-卩Xi)/ gi。于是功能函數(shù)也可轉(zhuǎn)換到標準正態(tài)空間，即-1-1g(X)=g(R(X),S(X)=g(R(T(Y),S(T(Y)=G(Y)失效概率Pf=1- U ( 3 ),由B的幾何含義可

33、知，3值為在標準正態(tài)空間中從原點到失效面 的最短距離。在標準正態(tài)空間中概率密度是關(guān)于原點(均值點)旋轉(zhuǎn)對稱，并且隨著到原點 的距離的平方呈指數(shù)下降。在一次二階矩法中，標準正態(tài)空間中的極限狀態(tài)面(失效面)被 一個到原點最小距離點處的切平面代替，也即按一階Taylor展開式把功能函數(shù)在設(shè)計點處展開，使之線性化。采用迭代法可以近似確定極限狀態(tài)面G(Y)=0上距原點最近的點 Y,然后按距離公式確定結(jié)構(gòu)的可靠指標。具體的迭代格式如下：Y(i)GN)芒砂丫)丫j 土CY j丫2nzj 4；：G(Y)(1,2,n)心(746)以Mises屈服準那么為例， 于是極限狀態(tài)方程計算可靠指標3時要用到功能函數(shù) G(

34、Y)的梯度向量，由求導鏈式法那么功能函數(shù) g(S)=sg(S(X(Y)=G(Y)=02 20 -(s xx-s2 2:xSyy+S yy+3s xy)(7.4.7)X：(748),G (Y)；:g (S) ：S(X) ：X(Y)浮-:S:X:Y由于S(X)難于用顯式表達，求；S(X)/；X存在困難。有效方法是利用隨機有限元一次二階矩法計算結(jié)構(gòu)的可靠指標3和設(shè)計點Y*。步驟如下：(1) 確定隨機變量轉(zhuǎn)換關(guān)系 Y=Y(X);(2) 給定初值Y(0) =0，計算式(7.4.8 )右邊的三個偏導數(shù)；(3) 按式(7.4.5 )和(7.4.6 )計算設(shè)計驗算點坐標丫；(4) 用隨機有限元法計算 S(i)

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

36、nf(Xi)丨，這里，f(X)是隨機變量 X的概率密度函數(shù)；f(x i)是離散點概率i 函數(shù)。最大熵法通過調(diào)整概率密度函數(shù)f(X)使熵S取得最大值，并滿足約束條件：Rf(X)dX =1和JRXi f (X)dm ,利用最大熵法可得近似的概率密度函數(shù)f(X)。為了求得最大值，可引入拉格朗日乘子法，把問題轉(zhuǎn)化為確定拉格朗日乘子。選用基于Neumann級數(shù)的隨機有限元法，寫出遞推方程。可將局部平均理論引入 Neumann 隨機有限元法，以進一步提高計算效率。隨機有限元最大熵法計算結(jié)構(gòu)失效概率的步驟如下：(1) 隨機場離散，隨機變量化為隨機向量，并計算協(xié)方差矩陣；(2) 利用特征正交化技術(shù)，用一組統(tǒng)計

37、獨立的隨機變量描述隨機場；(3) 利用Neumann隨機有限元法計算響應(yīng)量的前假設(shè)干階矩；(4) 擬合所有節(jié)點或局部節(jié)點的最大熵密度函數(shù)；(5) 利用數(shù)值積分計算結(jié)構(gòu)的失效概率。算例例1如圖743所示桁架下端受一集中力，假設(shè)各設(shè)計變量都是均勻分布的隨機變量，三根桿彈性模量 E,卩 e=2.0 X 106N/cm2, Cov(E)=2.5% ;截面積 A,卩 A=0.8cm2, Cov(A)=2%; 桿 2 長度 L2,卩 L2=20cm, Cov(L2)=1.5% ;桿 1 和桿 3 的長度均為 L1,卩 “=28.284cm , Cov(L=2.3% ;外力 P,卩 p=4000N; Cov

38、(P)=10%;材料屈服應(yīng)力 b( E )或等價的功率譜密度函數(shù)為Saa( 3 )和Sbb( 3 )，相關(guān)偏度分別是 E和E任一點橫向位移的有限元模式為w(x,t)=Neqe,其中，單元位移列陣qe=wi, B i,W2, 0 2T,形函數(shù)Nf取三次多項式插值。單元的動能表達為：e 1 le盯 p.m(x)Ae巴 dxqeTMeqe?t2M M(7.5.1)e式中：Me=MMe()=jem(Ne)T AeNedx + mf b(x)(Ne)TAeNedx, A為單元橫截面面積。單元的應(yīng)變能為e 1 leUeE(x)_ 2 W-2rx2I eT e e dx q K q12(7.5.2)e式中：

39、Ke=KKe():El(N；dx+El fa(x)(N； f N；dx, I 是慣性矩。單元的外力功為we；FieqeTKGcqe ；QqeTK；Dqe(7.5.3)式中：KGcN：N：dx KGdT(L2le/l)N；dx, Qx/I =(Qx/l) Qxp(x)/l , p為任一隨機變量據(jù)此集成結(jié)構(gòu)總能量和總外力功,利用t212Hamilton 原理 、(T - U)dt: WCdt = 0 ,t1飛即對總位移向量q取變分，考慮到q的任意性，可得系統(tǒng)運動方程， 進而得到如下特征問題(7.5.4 )方程：-s2M (-K) pKGc QKGDq=0le式中，剛度系數(shù)kj二kj - kj C1

40、)，其中kj是確定性分量，kj二El Nj(x)N j(x)dx，而隨機擾動項：匕(門)=El ； a(x)N(x)Nj(x)dx。于是可以寫出剛度系數(shù)、質(zhì)量系數(shù)以及 幾何剛度系數(shù)的均值和方差的表達式。在總剛度矩陣中兩個剛度系數(shù)kij和krs的互協(xié)方差,當單元尺寸相等時，利用方差函數(shù)可作進一步簡化。利用局部平均理論，在每個單元上E、m Q的均值為零，方差為：Var(Ee) *E E(le) *Ee/leVar(me) *m m(le)(7.5.5) Var(Qe) V p(le-2p/le式中，E(le)、m(le)、p (le)和GE、刁m、 p分別為隨機變量a(x)、b(x)和p(x)的方

41、差函數(shù)和相關(guān)偏度。利用協(xié)方差函數(shù)寫出兩剛度系數(shù)或兩質(zhì)量系數(shù)的互協(xié)方差。最后，特征值問題的方程可表達為K K()PKgg PKgcD QKgd QKgd)x= M M)x(7.5.6)* *或 K x=入M x。相應(yīng)平均問題的方程為(7.5.7)K PKgc QKgd x = - M x利用單元剛度、質(zhì)量和幾何剛度之間的協(xié)方差的表達式，等效剛度矩陣K*的協(xié)方差矩陣能用獨立的a(x)、b(x)和p(x)來表示。解未擾動的特征值問題(平均值問題)可得特征值 的均值。對于滿足小擾動情況，可推導出求兩個特征值之間的協(xié)方差公式。在建立了材料特性和幾何尺寸有隨機擾動的梁-柱的隨機特征值問題的隨機有限元公式后

42、，利用局部平均理論，隨機過程由均值、方差和相關(guān)偏度來定義。上述公式可推廣到二維或三維24。 線性結(jié)構(gòu)系統(tǒng)的瞬態(tài)響應(yīng)設(shè)結(jié)構(gòu)系統(tǒng)的質(zhì)量矩陣M是確定性的，阻尼、剛度和外力的概率分布由廣義協(xié)方差矩陣Cov(bi,bj) ( i,j=1,2, ,q )表示，線性結(jié)構(gòu)系統(tǒng)的運動方程為Mx(b,t) C(b)x(b,t) K(b)x(b,t) = F(b,t)(7.5.8)x和荷載向量F為n式中，質(zhì)量陣 M阻尼陣C(b)和剛度陣K(b)都是n階方陣，位移向量階列陣。用Taylor級數(shù)對隨機向量b展開，并保存二階項，那么位移向量 x(b,t)關(guān)于b的二階攝動式為-q -x(b,t) =x(t) 、Xbi(t) bi 二Xbibj (t)心b 心bj(7.5.9)i,j呂其中，丫是小參數(shù)；b是b的均值；X(t)、Xbi (t)和Xbibi (t)分別表示位移向量的均值、位移向量相對于bi的