第7章 隨機有限元法.doc

第7章 隨機有限元法7.1 緒論結(jié)構工程中存在諸多的不確定性因素，從結(jié)構材料性能參數(shù)到所承受的主要荷載，如車流、陣風或地震波，無不存在隨機性。在有限單元法已成為分析復雜結(jié)構的強有力的工具和廣泛使用的數(shù)值方法的今天，人們已不滿足精度越來越高的確定性有限元計算，而設法用這一強有力的工具去研究工程實踐中存在的大量不確定問題。隨機有限元法（Stochastic FEM）,也稱概率有限元法（Probabilistic FEM）正是隨機分析理論與有限元方法相結(jié)合的產(chǎn)物，是在傳統(tǒng)的有限元方法的基礎上發(fā)展起來的隨機的數(shù)值分析方法。最初是Monte-Carlo法與有限元法直接結(jié)合，形成獨特的統(tǒng)計有限元方法。Astill和Shinozuka（1972）首先將Monte-Carlo法引入結(jié)構的隨機有限元法分析。該法通過在計算機上產(chǎn)生的樣本函數(shù)來模擬系統(tǒng)的隨機輸入量的概率特征，并對于每個給定的樣本點，對系統(tǒng)進行確定性的有限元分析，從而得到系統(tǒng)的隨機響應的概率特征。由于是直接建立在大量確定性有限元計算的基礎上，計算量極大，不適用于大型結(jié)構，而且最初的直接Monte-Carlo法還不是真正意義上的隨機有限元法。但與隨后的攝動隨機有限元法（PSFEM）相比，當樣本容量足夠大時，Monte-Carlo有限元法的結(jié)果更可靠也更精確。結(jié)構系統(tǒng)的隨機分析一般可分為兩大類：一類是統(tǒng)計方法，另一類是非統(tǒng)計方法。因此，隨機有限元法同樣也有統(tǒng)計逼近和非統(tǒng)計逼近兩種類型。前者通過樣本試驗收集原始的數(shù)據(jù)資料，運用概率和統(tǒng)計理論進行分析和整理，然后作出科學推斷。這里，樣本試驗和數(shù)據(jù)處理的工作量很大，隨著計算機的普及和發(fā)展，數(shù)值模擬法，如蒙特卡羅（Monte Carlo）模擬，已成為最常用的統(tǒng)計逼近法。后者從本質(zhì)上來說是利用分析工具找出結(jié)構系統(tǒng)的（確定的或隨機的）輸出隨機信號與輸入隨機信號之間的關系，采用隨機分析與求解系統(tǒng)控制方程相結(jié)合的方法得到輸出信號的各階隨機統(tǒng)計量的數(shù)字特征（如各階原點矩或中心矩）。在20世紀70年代初， Cambou首先采用一次二階矩方法研究線彈性問題。由于這種方法將隨機變量的影響量進行Taylor級數(shù)展開，就稱之為Taylor展開法隨機有限元（TSFEM）。Shinozuka和Astill（1972)分別獨立運用攝動技術研究了隨機系統(tǒng)的特征值問題。隨后，Handa（1975）等人在考慮隨機變量波動性時采用一階和二階攝動技術，并將這種攝動法隨機有限元成功地應用于框架結(jié)構分析。Vanmarcke等人（1983）提出隨機場的局部平均理論，并將它引入隨機有限元。局部平均理論是用隨機場函數(shù)在每一個離散單元上的局部平均的隨機變量來代表該單元的統(tǒng)計量的近似理論。Liu W. K.等人（1986、1988）的系列工作，提供了一種“主模態(tài)”技術，運用隨機變量的特征正交化方法，將滿秩的協(xié)方差矩陣變換為對角矩陣，減少計算工作量，對攝動隨機有限元法的發(fā)展做出貢獻，此外，提出了一個隨機變分原理。Yamazaki和Shinozuka(1987)創(chuàng)造性地將算子的Neumann級數(shù)展開式引入隨機有限元的列式工作。從本質(zhì)上講，Neumann級數(shù)展開方法也是一類正則的小參數(shù)攝動方法，正定的隨機剛度矩陣和微小的隨機擾動量是兩個基本要求，這兩個基本要求保證了攝動解的正則性和收斂性，其優(yōu)點在于攝動形式較簡單并可以得到近似解的高階統(tǒng)計量。Shinozuka等人（1987）將隨機場函數(shù)的Monte-Carlo模擬與隨機剛度矩陣的Neumann級數(shù)展開式結(jié)合,得到具有較好計算精度和效率的一類Neumann隨機有限元列式（稱NSFEM）。Benaroya等（1988）指出，將出現(xiàn)以隨機變分原理為基礎的隨機有限元法來逐漸取代以攝動法為基礎的隨機有限元法。Spanos和Ghanem等人（1989，1991）結(jié)合隨機場函數(shù)的Karhuen-Loeve展式和Galerkin（迦遼金）射影方法建立了相應的隨機有限元列式，并撰寫了隨機有限元法領域的第一本專著隨機有限元譜方法。國內(nèi)對隨機有限元的研究起步較晚。吳世偉等人（1988）提出隨機有限元的直接偏微分法及相應的可靠度計算方法。陳虬、劉先斌等人（1989、1991）提出一種新的隨機場離散模型，建立了等參局部平均單元，并基于變分原理研究了一類隨機有限元法的收斂性和誤差界。Papadrakakis(1995)采用預處理共軛梯度法給出了空間框架的非線性隨機有限元列式。Schorling和Bucher(1996)基于Monte-Carlo技術，采用響應面法研究幾何非線性時的可靠度隨機有限元方法。劉寧（1996）則基于偏微分法，給出了三維彈塑性隨機有限元列式。隨機有限元法的數(shù)學理論研究和非線性隨機問題的有限元分析工作還有待深入。自20世紀80年代以來，隨機有限元法已在工程結(jié)構可靠性、安全性分析領域以及在各種隨機激勵下結(jié)構響應變異研究領域中得到應用，如應用于大型水利工程的重力壩、拱壩的可靠度計算；應用于非線性瞬態(tài)響應分析；結(jié)構振動中隨機阻尼對響應的影響；結(jié)構分析的隨機識別；復雜結(jié)構地震響應的隨機分析和兩相動力系統(tǒng)的隨機模擬等等。隨著理論研究的深入，隨機有限元將得到更加廣泛的應用。7.2 隨機有限元的控制方程22從隨機有限元控制方程的獲得來看，隨機有限元可分為Taylor展開法隨機有限元（TSFEM）、攝動法隨機有限元(PSFEM)以及Neumann展開Monte-Carlo法隨機有限元（NSFEM）。 Taylor展開法隨機有限元 該隨機有限元法的基本思路是將有限元格式中的控制量在隨機變量均值點處進行Taylor級數(shù)展開（取一階或二階），經(jīng)過適當?shù)臄?shù)學處理得出所需的計算方程式。有限元靜力分析控制方程的矩陣形式為: KU = F (7.2.1)式中,U為位移矩陣，F(xiàn)為等效節(jié)點荷載列陣，K為整體剛度矩陣 (7.2.2)其中，B為形變矩陣,D為材料彈性矩陣。在計算出節(jié)點位移U后，即由下式求得應力列陣= DBU (7.2.3)設基本隨機變量為，將位移U在均值點處一階Taylor級數(shù)展開，并在兩邊同時取均值（數(shù)學期望），得 (7.2.4)式中：符號E表示求均值，任一結(jié)點位移U的方差可由下式計算： (7.2.5)式中：符號Var表示求方差；Cov(Xi,Xj)為Xi和Xj的協(xié)方差。其中 (7.2.6) (7.2.7)同樣將在均值點處Taylor展開，也有與上面類似的表達式。可見，TSFEM關鍵在于對有限元方程式直接進行偏微分計算，計算出有限元輸出量對隨機變量的梯度，故該法也稱直接偏微分法或梯度分析法。由于一階TSFEM只需一次形成剛度矩陣，也只需一次求剛度矩陣的逆，因此效率較高。但由于忽略了二階以上的高次項，使TSFEM對隨機變量的變異性有所限制。一般要求一階TSFEM隨機變量的變異系數(shù)小于0.3。如果隨機變量的變異系數(shù)較大，可以采用有限元控制方程的二階Taylor展開： (7.2.8) (7.2.9)上式可見，二階TSFEM可以放寬隨機變量變異性大小的限制，但隨機變量數(shù)目較多時，計算量將十分龐大，而且一階或二階TSFEM均無法計算響應量三階以上的統(tǒng)計特性。由于TSFEM簡單明了、效率高，為我國許多學者所采用。 攝動法隨機有限元攝動技術最初被用于非線性力學分析。Handa等人成功地將一階、二階攝動技術用于隨機問題，給出攝動法有限元列式。該法假定基本隨機變量在均值點處產(chǎn)生微小攝動，利用Taylor級數(shù)把隨機變量表示為確定部分和由攝動引起的隨機部分，從而將有限元控制方程（非線性的）轉(zhuǎn)化為一組線性的遞推方程，求解得出位移的統(tǒng)計特性，進而求出應力的統(tǒng)計特性。假設為隨機變量在均值點處的微小攝動量，即。于是 (7.2.10)對于U、F，也有類似上式K的表達式，式中：K0、U0、F0分別為K、U、F在隨機變量均值點的值。根據(jù)二階攝動法，可得 (7.2.11) (7.2.12) (7.2.13)由上式可得位移的均值和協(xié)方差： (7.2.14) (7.2.15)由于任何量的隨機性都可以引入攝動量，而且更易于考慮非線性問題，因此PSEFEM適用范圍較廣，對于結(jié)構幾何特性的隨機性（包括隨機邊界問題）易得出隨機有限元控制方程。一階PSFEM和一階TSFEM一樣，只需一次形成剛度矩陣、一次對剛度矩陣求逆，計算效率較高。但PSFEM需以微小的攝動量為條件，一般應小于均值的20%或30%。 Neumann展開Monte-Carlo隨機有限元20世紀80年代后期，Shinozuka等人提出基于Neumann展開式的隨機有限元法，使Monte-Carlo法與有限元法較完美地結(jié)合起來。Monte-Carlo法是最直觀、最精確、獲取信息最多、對非線性問題最有效的計算統(tǒng)計方法。Neumann展開式的引入是為了解決矩陣求逆的效率問題。如果對每一次隨機抽樣，只需形成剛度矩陣，進行前代、回代以及矩陣乘和矩陣加減，而無需矩陣分解，則可大大減少工作量。在一般有限元控制方程KU = F中，假定荷載F為確定值，在隨機變量波動值的影響下剛度矩陣K分解為K = K0+K，根據(jù)Neumann級數(shù)展開，有 K-1=(K0+K)-1=(I-P+P2-P3+)K0-1 (7.2.16)式中：K0為隨機變量均值處的剛度矩陣；K為剛度矩陣的波動量；I為單位矩陣。對于Monte-Carlo隨機抽樣，剛度矩陣只改變K項，而P = K0-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)用確定性有限元的計算程序，而TSFEM在編程上較為復雜，PSFEM則更為復雜。由于采用Monte-Carlo隨機模擬技術，NSFEM不受隨機變量波動范圍的限制，當變異系數(shù)小于0.2時，NSFEM與一階TSFEM或一階PSFEM精度相當；當變異系數(shù)大于0.2時，后兩者已不能滿足精度要求，但NSFEM仍能得出滿意的結(jié)果。7.3 隨機場的離散模型許多物理現(xiàn)象和物體系統(tǒng)具有隨機分布特性，包括系統(tǒng)本身的不確定或系統(tǒng)的激勵和響應的不確定，都可以模型化為隨機空間分布的隨機場或隨時間分布的隨機過程。隨機有限元法除了必須考慮材料參數(shù)等的空間變異性，需要獲得隨機有限元方程列式以及解決隨機算子和隨機矩陣的求逆問題外，還須包含對隨機場的離散處理。均勻各向同性隨機場的特征量1. 隨機場S(t)的均值E(S(t)為常數(shù)m.2. 隨機場S(t)的標準相關函數(shù) (7.3.1)式中，隨機場的協(xié)方差函數(shù) R()= CovS(t+),S(t) (7.3.2)對于一切t，隨機場的方差為 VarS(t)= R(0)= 2 (7.3.3)相關函數(shù)也可以用譜分解表示（即Wiener-Khintchine變換對）： (7.3.4) 常用的標準相關函數(shù)有：非協(xié)調(diào)階躍型、協(xié)調(diào)階躍型、三角型、指數(shù)型、二階AR型、高斯型等多種形式。3. 方差折減系數(shù)2(h)設Sh(t)是隨機場S(t)的局部平均隨機過程，即 (7.3.5)則方差 VarSh(t)= 22(h) (7.3.6)式中，方差折減系數(shù) (7.3.7)可見，方差折減系數(shù)起著使Sh(t)的方差比原來S(t)的方差縮小的作用。4. 相關距離SS可以看成是任意兩個相隔距離為的隨機變量不相關的最小距離（也稱相關偏度）。，則不相關，否則完全相關。利用相關距離S便于對隨機場作近似處理，其計算式為： (7.3.8) 以上公式均對一維問題列出，二維、三維問題也可以類似得出。 隨機場的中心離散中心離散是隨機場最簡單的一種離散方法。該法用隨機場在每個單元中心點的值來表征該隨機場在每個單元的屬性，因而隨機場在每個單元內(nèi)部都是常量，且等于它在各個單元中心的值。該法程序簡單，但精度欠佳，現(xiàn)較少采用。 隨機場的局部平均 該法將隨機場在每個單元的屬性用隨機場在單元上的局部平均來表征。Baecher、Vanmarcke首先提出隨機場局部平均的離散方法，我國學者曾對該法進行了深入研究2324。1. 一維隨機場的局部平均設一個一維連續(xù)平穩(wěn)隨機場S(x)，其均值為m，方差為2，則隨機場在一個離散單元t-T/2,t+T/2上的局部平均定義為 (7.3.9)式中：T是局部平均單元的長度，ST(t)稱為局部平均隨機場，其均值也為m，方差按（7.3.6）式計算，h=T ，S(t)的標準相關函數(shù)為()=R()/2，無量綱方差折減系數(shù)2（h）有如下性質(zhì)：2（T）0；2（0）=1；2（-T）=2（T）?？紤]隨機場在兩個長度分別為T和T的單元上的局部平均。如果局部平均隨機場分別為ST和ST，則其協(xié)方差為 (7.3.10)2. 二維隨機場的局部平均設S(x1，x2)為二維連續(xù)參數(shù)連續(xù)狀態(tài)隨機場，Ai=L1iL2i為一矩形單元中心點（x1i，x2i）為中心，邊平行于坐標軸x1與x2，且長為L1i與L2i的矩形面積，隨機場在Ai內(nèi)的局部平均定義為 (7.3.11)如果S(x1,x2)是一個均勻隨機場，則可用均值m、方差2及歸一化協(xié)方差函數(shù)（r1,r2）近似描述，r1和r2分別為沿兩個方向的距離。對應的局部平均隨機場可用E(Si)、Var(Si)及互協(xié)方差Cov(Si,Sj)近似描述。E(Yi)、Var(Yi)及Cov(Yi,Yj)可由m、2及(r1，r2)計算獲得。第i和第j單元局部平均Si、Sj的互協(xié)方差可表示為 (7.3.12)式中：的約定如圖7.3.1所示。圖7.3.1 二維局部平均單元 如果有限元的網(wǎng)格已劃分，且單元總數(shù)為n，隨機場實際上被離散成n個隨機變量，這n個隨機變量的統(tǒng)計特性可由ESi、VarSi及互協(xié)方差CovSi，Sj反映。由于有限元網(wǎng)格的疏密是由應力梯度決定的，與隨機場無關，當單元數(shù)很多時，隨機場可另劃網(wǎng)格，網(wǎng)格的疏密可由相關距離S決定。當然，隨機場網(wǎng)格越密精度越高。隨機場的局部平均法由于對原始數(shù)據(jù)的要求低、收斂快、精度高，是隨機有限元計算中最常用的方法。 隨機場的插值Liu W. K.提出隨機場的插值法。該法將隨機場在單元內(nèi)的值用單元結(jié)點處值的插值函數(shù)來表示，于是隨機場的統(tǒng)計特性可由各單元結(jié)點處隨機變量間的統(tǒng)計特性近似反映。利用形函數(shù)Ni(X)，隨機場b(X)離散式表示為 (7.3.13)式中：X表示空間位置；bi為隨機場在結(jié)點i處的值（i=1，2，q）；q為單元結(jié)點數(shù)。隨機場b(X)在單元內(nèi)的均值和方差可表示為 (7.3.14) (7.3.15) 隨機場的插值法將原連續(xù)狀態(tài)的隨機場仍離散成一個連續(xù)函數(shù)，未直接計算隨機場引起的單元間的相關性，只需給定隨機場在各結(jié)點上的值，計算相對簡單，易于考慮非線性和非均勻隨機場問題。但需要已知相關函數(shù)，并且要求隨機場對空間參數(shù)具有較高的連續(xù)性。 隨機場的加權積分法Takada、Shinozuka及Deodatis提出隨機場加權積分方法。該法在單元剛度矩陣的推導過程中采用隨機場在單元高斯點上的加權積分，以表征單元上的隨機場。假設單元e的彈性模量為 (7.3.16)式中，E0(e)(X)為彈性模量的均值；f(e)(X)為一維零均值均勻隨機場，其值域為 -1+f(e)(X)-1- (01)兩端鉸接桿單元的剛度矩陣可近似表示為 K(e)K0(e)+X(e)K0(e) (7.3.17)式中，K0(e)為彈性模量取均值時的單元剛度矩陣，而 （L為桿長） (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）實質(zhì)上是在考慮彈性模量隨機場的情況下，將單元剛度矩陣分解成確定性部分和隨機部分。式（7.3.18）等表征隨機場在各單元的均值和方差，而單元間彈性模量的相關性則由下式表示 (7.3.20) 不難看出，局部平均法是加權積分法的特例（即權系數(shù)全部相同）。由于采用加權積分，其計算精度相對較高，而且該法積分只需一次進行，剛度矩陣的波動性也由此得出，因此計算效率也較高。 隨機場的正交展開Spanos和Ghanem提出的隨機場正交展開法，將材料特性參數(shù)隨機場進行Karhumen-Loeve正交展開，并由此推導出剛度矩陣的級數(shù)展開式，從而獲得位移、應力的統(tǒng)計特性。設隨機場為 (7.3.21)其中 式中：n、n(x)分別為隨機場S(x)相關函數(shù)的特征值和特征函數(shù)。n(x)具有如下正交性： （Kronecker函數(shù)） (7.3.22)式（7.3.21）對于任何分布的S(x)均收斂，可取至第r階以滿足精度要求。對于一維桿單元，單元剛度矩陣可表示為 其中 式中：Pe(x)=De(x)/S(x), De(x)為單元彈性矩陣。集成整體剛度矩陣，得 (7.3.23)從上式，可得位移： (7.3.24)利用Neumann展開式，可進一步得出位移的統(tǒng)計特性。該法關鍵在于獲得特征值和特征函數(shù)。7.4 隨機有限元與結(jié)構可靠度 結(jié)構可靠性分析的一次二階矩法結(jié)構的安全性、適用性和耐久性統(tǒng)稱為結(jié)構的可靠性?？煽啃缘臄?shù)學量度為可靠度，其定義為：在規(guī)定的條件下和規(guī)定的時間內(nèi)完成預定功能的概率?！耙?guī)定條件”是指正常設計、正常施工、正常使用的條件；“規(guī)定時間”是指結(jié)構的設計基準期。結(jié)構的使用時間超過基準期后，其失效的概率將增大。結(jié)構的一系列基本變量都具有不確定性，因此結(jié)構的可靠性分析屬于隨機性分析的范疇，隨機有限元法將是十分有效的工具。在結(jié)構可靠性分析中，結(jié)構的極限狀態(tài)（包括承載能力極限狀態(tài)、正常使用極限狀態(tài)和條件極限狀態(tài)）是通過功能函數(shù)來描述的。當有n個隨機變量影響結(jié)構可靠度時，結(jié)構的功能函數(shù)為: z=g(X1,X2,Xn)，基本變量Xi(i=1,2,n)是結(jié)構上的各種外因作用、材料性能和幾何參數(shù)等。z0時，結(jié)構處于可靠狀態(tài)；z0,則處于失效狀態(tài)；z=0稱結(jié)構極限狀態(tài)方程（一般難以用顯式表示），結(jié)構處于極限狀態(tài)。如果z的概率密度函數(shù)或概率分布函數(shù)可求得，則結(jié)構可靠度的數(shù)量指標便可基于各種狀態(tài)出現(xiàn)的概率而確定。若功能函數(shù)僅與兩個隨機變量有關（如結(jié)構抵抗破壞或變形的能力R和荷載引起結(jié)構內(nèi)力、應力、位移等效應S），即z=g(R,S)。假設R、S均為正態(tài)分布，其均值和標準差分別為和，此時z=R-S也是正態(tài)分布的隨機變量，具有均值和標準差。其概率密度函數(shù)為 (7.4.1)其分布圖示于圖7.4.1，陰影部分是結(jié)構失效概率Pf，非陰影部分面積即結(jié)構的可靠度Pr。 圖7.4.1 正態(tài)分布概率密度函數(shù) 圖7.4.2 失效邊界在工程實踐中R和S不一定是正態(tài)分布，但可以變換成標準正態(tài)分布。他們的統(tǒng)計特征量是均值、標準差、相關偏度或變異系數(shù)等，變異系數(shù)=/，表示隨機變量相對于均值的變異。目前工程上一般用無量綱的可靠指標來反映結(jié)構的可靠度，越大，失效概率Pf越小，其互補的可靠度Pr就越大。一次二階矩法采用只需已知均值和標準差的數(shù)學模型去求解結(jié)構的可靠度。此法將功能函數(shù)z=g(X1,X2,Xn)在某點用Taylor級數(shù)展開，并近似地取一次項使極限狀態(tài)方程線性化，然后求得可靠指標。改進的一次二階矩法將線性化點選在失效邊界上，而且選在結(jié)構最大可能失效點P*上（如圖7.4.2）。選擇設計驗算點P*(Xi*i=1,2,n)作為線性化點X0i時，線性化的極限狀態(tài)方程為 (7.4.2)則z的均值為 由于設計驗算點P*選在失效邊界上，有g(X1*,X2*,Xn*)=0，因此z成為 但隨機變量Xi互不相關時，z的標準差z為 其中 稱i為靈敏度系數(shù)，表示第i個隨機變量對標準差的相對影響。于是可靠指標為 (7.4.3)變換上式為 即 (i=1，2，n) (7.4.4)上式中，Xi，Xi為已知的各隨機變量的均值和標準差，待求的量為Xi*和，可迭代求解。 隨機有限元的一次二階矩法設可靠性分析中的一組基本變量X=(X1,X2,Xn)T為相互獨立的正態(tài)變量（不滿足時可作變換），并進一步變換為標準正態(tài)變量Y=(Y1,Y2,Yn)T，其中Yi=(Xi-Xi)/Xi。于是功能函數(shù)也可轉(zhuǎn)換到標準正態(tài)空間，即 g(X)=g(R(X),S(X)=g(R(T-1(Y),S(T-1(Y)=G(Y)失效概率Pf=1-(),由的幾何含義可知，值為在標準正態(tài)空間中從原點到失效面的最短距離。在標準正態(tài)空間中概率密度是關于原點（均值點）旋轉(zhuǎn)對稱，并且隨著到原點的距離的平方呈指數(shù)下降。在一次二階矩法中，標準正態(tài)空間中的極限狀態(tài)面（失效面）被一個到原點最小距離點處的切平面代替，也即按一階Taylor展開式把功能函數(shù)在設計點處展開，使之線性化。采用迭代法可以近似確定極限狀態(tài)面G(Y)=0上距原點最近的點Y*,然后按距離公式確定結(jié)構的可靠指標。具體的迭代格式如下： (7.4.5) (7.4.6)以Mises屈服準則為例，功能函數(shù) g(S)=s02-(s2xx-sxxsyy+s2yy+3s2xy)于是極限狀態(tài)方程 g(S(X(Y)=G(Y)=0 (7.4.7)計算可靠指標時要用到功能函數(shù)G(Y)的梯度向量，由求導鏈式法則 (7.4.8)由于S(X)難于用顯式表達，求存在困難。有效方法是利用隨機有限元一次二階矩法計算結(jié)構的可靠指標和設計點Y*。步驟如下：（1） 確定隨機變量轉(zhuǎn)換關系Y=Y(X)；（2） 給定初值Y(0)=0，計算式（7.4.8）右邊的三個偏導數(shù)；（3） 按式（7.4.5）和（7.4.6）計算設計驗算點坐標Y*；（4） 用隨機有限元法計算S(i)和；（5） 重復（3）、（4）兩步驟，直至收斂（G(Y(i)=0）。（6） 再按式（7.4.5）和（7.4.6）計算可靠指標值。 隨機有限元的最大熵法一次二階矩法的計算量隨迭代次數(shù)成倍增加，使該法的使用受到限制。而最大熵法用于結(jié)構的可靠性分析時，可根據(jù)隨機變量的二階矩來擬合概率密度曲線。因此隨機有限元結(jié)合最大熵法可用于求結(jié)構響應的概率密度曲線，從而計算結(jié)構的失效概率。熵被定義為信息的均值，信息是對個別X值不確定性的度量。不確定性越大，熵也越大。對于一個連續(xù)隨機變量，熵為 ，而對于離散隨機變量，熵為，這里，f(X)是隨機變量X的概率密度函數(shù)；f(xi)是離散點概率函數(shù)。最大熵法通過調(diào)整概率密度函數(shù)f(X)使熵S取得最大值，并滿足約束條件： ，利用最大熵法可得近似的概率密度函數(shù)f(X)。為了求得最大值，可引入拉格朗日乘子法，把問題轉(zhuǎn)化為確定拉格朗日乘子。選用基于Neumann級數(shù)的隨機有限元法，寫出遞推方程?？蓪⒕植科骄碚撘隢eumann隨機有限元法，以進一步提高計算效率。隨機有限元最大熵法計算結(jié)構失效概率的步驟如下：(1) 隨機場離散，隨機變量化為隨機向量，并計算協(xié)方差矩陣；(2) 利用特征正交化技術，用一組統(tǒng)計獨立的隨機變量描述隨機場；(3) 利用Neumann隨機有限元法計算響應量的前若干階矩；(4) 擬合所有節(jié)點或部分節(jié)點的最大熵密度函數(shù)；(5) 利用數(shù)值積分計算結(jié)構的失效概率。 算例例1 如圖7.4.3所示桁架下端受一集中力，假設各設計變量都是均勻分布的隨機變量，三根桿彈性模量E，E=2.0106N/cm2，Cov(E)=2.5%；截面積A，A=0.8cm2，Cov(A)=2%；桿2長度L2,L2=20cm，Cov(L2)=1.5%；桿1和桿3的長度均為L1，L1=28.284cm，Cov(L1)=2.3%；外力P，P=4000N；Cov(P)=10%；材料屈服應力y，y=3200N/cm2，Cov(y)=7%。當極限狀態(tài)方程取為 y-S=0（S為構件中最大應力值）時，試求該結(jié)構的失效概率。圖7.4.3 簡單超靜定桁架解： 用隨機有限元一次二階矩法計算，假設由各結(jié)點坐標值變異引起的方向余弦的變異很小，形成確定矩陣，則由隨機有限元遞推方程組可得各結(jié)點位移分布特征值。計算得結(jié)點4的位移均值和標準差為4=2.2928cm和4=0.0037cm；各桿件中最大應力為桿元2的應力，其均值和標準差為2=2928N/cm2和2=377N/cm2。由一次二階矩法可計算得該結(jié)構的失效概率為0.23（由于結(jié)構超靜定，桿2失效并不意味著整個結(jié)構失效），與蒙特卡羅有限元法計算結(jié)果比較接近（蒙特卡羅取樣數(shù)為40時，失效概率達0.224）。隨機有限元法在分析混凝土重力壩可靠度方面已得到較多應用。研究者采用抽樣方法對重力壩施行隨機自動剖分，研究其可靠指標的統(tǒng)計規(guī)律和收斂于真解的規(guī)律，并進行結(jié)構優(yōu)化25。在對碾壓混凝土重力壩抗滑穩(wěn)定的體系可靠度研究中，還采用剛體極限平衡法、線彈性隨機有限元法和三維彈塑性隨機有限元法進行計算比較。7.5 隨機有限元動力分析方法 不確定結(jié)構的自振特性假設某一梁-柱結(jié)構承受隨機分布的軸向靜荷載N(x)，沿長度有微小擾動，同時在端部有軸向推力c為隨機變量，其方差為，軸向分布荷載可表示為式中p(x)是一維、均值為零的均勻隨機場，方差是P2，相關偏度是P，互相關函數(shù)是RPP()。結(jié)構彈性模量和質(zhì)量密度也隨機變化，分別表示為，其中分別是彈性模量E和質(zhì)量密度m的均值，方差分別為2E和2m；隨機變量a(x)和b(x)是兩個一維互不相關的、均勻的、均值為零的隨機場。自相關函數(shù)分別為Raa()和Rbb()或等價的功率譜密度函數(shù)為Saa()和Sbb()，相關偏度分別是E和m。任一點橫向位移的有限元模式為 w(x,t)=Neqe，其中，單元位移列陣qe=w1,1,w2,2T,形函數(shù)Ne取三次多項式插值。單元的動能表達為： (7.5.1)式中：，Ae為單元橫截面面積。單元的應變能為 (7.5.2)式中：，I是慣性矩。單元的外力功為 (7.5.3)式中：，p為任一隨機變量據(jù)此集成結(jié)構總能量和總外力功，利用Hamilton原理，即對總位移向量q取變分，考慮到q的任意性，可得系統(tǒng)運動方程，進而得到如下特征問題方程： （7.5.4）式中，剛度系數(shù)，其中是確定性分量，而隨機擾動項：。于是可以寫出剛度系數(shù)、質(zhì)量系數(shù)以及幾何剛度系數(shù)的均值和方差的表達式26。在總剛度矩陣中兩個剛度系數(shù)kij和krs的互協(xié)方差，當單元尺寸相等時，利用方差函數(shù)可作進一步簡化。利用局部平均理論，在每個單元上E、m、Q的均值為零，方差為： (7.5.5)式中，和分別為隨機變量a(x)、b(x)和p(x)的方差函數(shù)和相關偏度。利用協(xié)方差函數(shù)寫出兩剛度系數(shù)或兩質(zhì)量系數(shù)的互協(xié)方差。最后，特征值問題的方程可表達為 (7.5.6)或 K*x=M*x 。相應平均問題的方程為 (7.5.7) 利用單元剛度、質(zhì)量和幾何剛度之間的協(xié)方差的表達式，等效剛度矩陣K*的協(xié)方差矩陣能用獨立的a(x)、b(x)和p(x)來表示。解未擾動的特征值問題（平均值問題）可得特征值的均值。對于滿足小擾動情況，可推導出求兩個特征值之間的