靳 慧

（1.東南大學(xué) 江蘇省工程力學(xué)分析重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210096；2.重慶交通大學(xué) 橋梁結(jié)構(gòu)工程重點(diǎn)實(shí)驗(yàn)室，重慶 400074）

隨機(jī)有限元方法（SFEM）是計(jì)算隨機(jī)力學(xué)的一個(gè)重要分支，是隨機(jī)結(jié)構(gòu)分析的有力工具.SFEM能夠處理工程中的結(jié)構(gòu)隨機(jī)問題，計(jì)算結(jié)構(gòu)隨機(jī)響應(yīng)量的統(tǒng)計(jì)特征.近幾十年來，國內(nèi)外學(xué)者對(duì)隨機(jī)有限元法進(jìn)行了大量深入的研究［1－3］，提出了眾多的求解方法，如Monte－Carlo隨機(jī)有限元法，Taylor展開隨機(jī)有限元法（TSFEM），攝動(dòng)隨機(jī)有限元法（PSFEM），Neumann展開隨機(jī)有限元法（NSFEM）等，成功地解決了大量工程中的隨機(jī)結(jié)構(gòu)分析問題.

除了Monte－Carlo隨機(jī)有限元法外，其他隨機(jī)有限元法都是在傳統(tǒng)的確定性有限元列式的基礎(chǔ)上，進(jìn)行推導(dǎo)演化，得出迭代式，以計(jì)算響應(yīng)量對(duì)各隨機(jī)變量的各階偏導(dǎo)數(shù)，然后根據(jù)這些偏導(dǎo)數(shù)和泰勒級(jí)數(shù)展式計(jì)算出響應(yīng)量的各階統(tǒng)計(jì)特性.可見，偏導(dǎo)數(shù)的計(jì)算是隨機(jī)有限元法的核心.

另外，隨機(jī)有限元法的重要用途是和可靠度理論相結(jié)合進(jìn)行結(jié)構(gòu)可靠性分析，稱為隨機(jī)有限元可靠度法（SFEMR）.在可靠度計(jì)算方法中，需要計(jì)算功能函數(shù)對(duì)各個(gè)隨機(jī)變量的偏導(dǎo)數(shù)，而這種偏導(dǎo)數(shù)的計(jì)算往往又通過功能函數(shù)的表達(dá)式轉(zhuǎn)化為結(jié)構(gòu)隨機(jī)響應(yīng)量的偏導(dǎo)數(shù)計(jì)算.所以，響應(yīng)量的偏導(dǎo)數(shù)計(jì)算直接影響SFEM和SFEMR的計(jì)算精度和效率.本文將復(fù)數(shù)變量求偏導(dǎo)法引入到SFEM和SFEMR中，提出一種快速、精確的求偏導(dǎo)方法，可大大簡化SFEM和SFEMR計(jì)算過程，提高計(jì)算效率.

1 隨機(jī)有限元法

一階攝動(dòng)法和一階Taylor展開法的計(jì)算過程較為簡單，適用于隨機(jī)變量變異系數(shù)較小的情況.

1.1 攝動(dòng)隨機(jī)有限元法

攝動(dòng)隨機(jī)有限元法假設(shè)結(jié)構(gòu)的某一個(gè)參數(shù)Y是隨機(jī)擾動(dòng)的，其擾動(dòng)量可以用一個(gè)隨機(jī)小參數(shù)α來表示［2］，Y＝（1＋α）為Y的均值.

有限元控制方程為

式中：K為總剛矩陣；v為位移列陣；F為等效節(jié)點(diǎn)載荷列陣.

當(dāng)有限元離散網(wǎng)格確定以后，將剛度矩陣、載荷列陣和位移列陣在隨機(jī)向量α（αi，αj，…）的均值處Taylor級(jí)數(shù)展開，并略去二階以上項(xiàng)，有

式中：K0，F(xiàn)0，v0分別表示剛度矩陣、載荷列陣和位移列陣的均值矩陣，它們是確定性量；下標(biāo)i表示對(duì)αi求偏導(dǎo)數(shù).將式（2），（3），（4）代入控制方程（1），運(yùn)用中心攝動(dòng)方法可以得到如下的遞歸方程組：

由式（5）可以得出位移的均值E和協(xié)方差Cov

1.2 Taylor展開隨機(jī)有限元法

該法的基本思路是將有限元格式中的控制變量在隨機(jī)變量均值點(diǎn)處進(jìn)行Taylor展開，經(jīng)適當(dāng)數(shù)學(xué)處理得出所需計(jì)算式［3］.

有限元控制方程為Kv＝F，應(yīng)力向量σ由下式計(jì)算：

式中：D為彈性矩陣；B為應(yīng)變矩陣.

X＝ （X1，X2，…，Xn）為表示結(jié)構(gòu)材料、幾何、載荷等隨機(jī)性質(zhì)的隨機(jī)向量.將應(yīng)力σ在均值點(diǎn)處一階Taylor展開，并在兩邊同時(shí)取均值，可得

式中：D0，B0分別表示彈性矩陣和應(yīng)變矩陣的均值，而

以上兩式與常規(guī)確定性有限元的計(jì)算完全相同.應(yīng)力的方差Var可由下式計(jì)算：

1.3 隨機(jī)有限元法中的偏導(dǎo)數(shù)計(jì)算

可以看出，SFEM中結(jié)構(gòu)位移和應(yīng)力的均值計(jì)算等同于確定性有限元；方差計(jì)算式（7）和（11）的關(guān)鍵在于求解位移和應(yīng)力對(duì)隨機(jī)變量的偏導(dǎo).對(duì)于線彈性結(jié)構(gòu)，直接對(duì)有限元方程求偏導(dǎo)可得

式中：ve表示單元位移陣列.

這種直接求偏導(dǎo)法可得到解析解，但需要對(duì)有限元矩陣進(jìn)行求導(dǎo)、求逆等多次運(yùn)算，過程復(fù)雜.

另外一種簡便求解偏導(dǎo)數(shù)的方法是有限差分法，函數(shù)f（x）對(duì)變量x的偏導(dǎo)數(shù)由下式計(jì)算

式中：h為微小擾動(dòng)量.有限差分法雖然避免了對(duì)矩陣進(jìn)行計(jì)算，只需求解函數(shù)值f，但精度差，且對(duì)擾動(dòng)量h的依賴性很強(qiáng).

復(fù)數(shù)變量法是一種精度和效率都極高的偏導(dǎo)數(shù)求解方法，本文將此法運(yùn)用到SFEM和SFEMR中.

2 復(fù)數(shù)變量求導(dǎo)法

復(fù)數(shù)變量求導(dǎo)法由Lyness和Moler于1967年提出［4］，Vatsa V N用于解 Navier－Stokes流動(dòng)方程的靈敏度偏導(dǎo)［5］，Daniel Contreras Mundstock等用在二維彈性的靈敏度估算問題中［6］，Gao X W用在非線性邊界元的求解中［7］.

將函數(shù)f（x）的變量x替換為復(fù)數(shù)變量x＋i h，h為微小擾動(dòng)量，f（x＋i h）可以展開為Taylor級(jí)數(shù)展式：

一階偏導(dǎo)數(shù)可以表示為

Im表示復(fù)數(shù)的虛部，忽略高階微小量O（h2），則有

式（17）即為復(fù)數(shù)變量法求一階導(dǎo)數(shù)的計(jì)算式，和其他方法相比，復(fù)數(shù)變量求導(dǎo)法有以下幾個(gè)特點(diǎn)：

（1）精度高.由式（17）可見，在復(fù)數(shù)空間進(jìn)行計(jì)算，一次導(dǎo)數(shù)計(jì)算時(shí)只需求一個(gè)復(fù)數(shù)函數(shù)值f（x＋i h），不存在多個(gè)函數(shù)之間的運(yùn)算誤差，算例表明具有很高的精度.

（2）擾動(dòng)量h可取微小值.在復(fù)數(shù)空間進(jìn)行計(jì)算，高階微小量O（h2）隨微小擾動(dòng)量h的減小而高次減小，h取微小值時(shí)，偏導(dǎo)計(jì)算精度高.有限差分法需要計(jì)算兩個(gè)函數(shù)值f（x＋h）和f（x－h(huán)）的差，當(dāng)h取微小值，這兩個(gè)函數(shù)值在一個(gè)數(shù)量級(jí)時(shí)，相減計(jì)算可能導(dǎo)致誤差.而復(fù)數(shù)變量求一次導(dǎo)數(shù)計(jì)算時(shí)只需求一個(gè)函數(shù)值，對(duì)擾動(dòng)量h的依賴性低.

（3）應(yīng)用簡單方便.

圖1為平面桁架結(jié)構(gòu)，載荷P1＝35500N，P2＝36100N，彈性模量E＝2.0×1011Pa，桿件橫截面積A＝0.00320m2.求桿件最大應(yīng)力值σmax對(duì)隨機(jī)變量A，P1的偏導(dǎo)數(shù).有限元程序采用平面桿單元編制，由直接求偏導(dǎo)法求得解析解.

列出一段復(fù)數(shù)變量求偏導(dǎo)的Matlab程序，這段程序是算例1中計(jì)算桁架結(jié)構(gòu)中的最大桿應(yīng)力響應(yīng)量Ymax對(duì)桿件橫截面積A的偏導(dǎo)數(shù)?Ymax／?A.

A＝0.00320m2；%A為桁架中所有桿件的橫截面積

P1＝35500N；%P1，P2為載荷值

P2＝36100N；

E＝2×1011Pa；%E為彈性模量

h＝1×10－6；%h為擾動(dòng)量

A＝A＋hi；%設(shè)置A為復(fù)數(shù)變量

FEM程序；%有限元程序計(jì)算Ymax

?Ymax／?A＝Im（Ymax）／h；%計(jì)算?Ymax／?A

可以看出，復(fù)數(shù)變量法只需在FEM程序前通過語句A＝A＋hi設(shè)置A為復(fù)數(shù)變量，在FEM程序后只用一條語句?Ymax／?A＝Im（Ymax）／h取出Ymax的虛部，并根據(jù)式（17）計(jì)算偏導(dǎo)，其余的程序完全與實(shí)數(shù)空間的FEM程序相同，計(jì)算所耗機(jī)時(shí)與普通有限元程序運(yùn)行機(jī)時(shí)并無差異.與直接求偏導(dǎo)法相比，省去了矩陣的求逆和求導(dǎo)等運(yùn)算，計(jì)算過程大大簡化，并且程序編制極其簡單，只需對(duì)實(shí)數(shù)空間的FEM程序進(jìn)行簡單的設(shè)置即可.

復(fù)數(shù)變量法和有限差分法的計(jì)算結(jié)果見表1.

圖1 平面桁架結(jié)構(gòu)Fig.1 Plane truss structure

表1 偏導(dǎo)數(shù)計(jì)算結(jié)果Tab.1 Results of derivative

由表1可見，計(jì)算?σmax／?A 時(shí)，只有h＝1.0×10－8時(shí)，有限差分法計(jì)算結(jié)果與解析解比較接近，而當(dāng)h取其他值時(shí)，都會(huì)產(chǎn)生較大的誤差；復(fù)數(shù)變量法在h取值足夠小，小于1.0×10－6時(shí)，都會(huì)得出很精確的結(jié)果.?σmax／?P1的計(jì)算也有類似的結(jié)論.

3 隨機(jī)有限元可靠度法

隨機(jī)有限元法的重要用途是和可靠度理論相結(jié)合進(jìn)行結(jié)構(gòu)可靠性分析，目前使用最多的是隨機(jī)有限元一次二階矩法.

3.1 改進(jìn)的一次二階矩法

設(shè)結(jié)構(gòu)構(gòu)件功能函數(shù)為［8］

式中：Xi（i＝1，2，…，n）為統(tǒng)計(jì)獨(dú)立正態(tài)隨機(jī)變量.

將功能函數(shù)在設(shè)計(jì)驗(yàn)算點(diǎn)x＊（，…）作Taylor級(jí)數(shù)展開，僅保留線性項(xiàng)，有

可得Z的均值與方差為

式中：uXi表示Xi的均值；σXi表示Xi的方差.

可靠度系數(shù)β為

設(shè)計(jì)驗(yàn)算點(diǎn)坐標(biāo)為

β的求取只能采用迭代法，計(jì)算步驟為：

（2）求cosθXi，采用式（23）.

（3）求β，采用式（21）.

3.2 隨機(jī)有限元一次二階矩法

利用式（21），（23）計(jì)算可靠度指標(biāo)β時(shí)，要用到功能函數(shù)的值g （，…）和功能函數(shù)對(duì)隨機(jī)變量的偏導(dǎo)數(shù)值.功能函數(shù)g往往是結(jié)構(gòu)響應(yīng)量s（位移，應(yīng)力等）的函數(shù)g＝f（s），所以求g的值和偏導(dǎo)數(shù)轉(zhuǎn)化為求s的值和偏導(dǎo)數(shù)，即

響應(yīng)量s的值和偏導(dǎo)可通過SFEM法求出，也就是將SFEM嵌入到一次二階矩的迭代中，形成隨機(jī)有限元一次二階矩法，其迭代步驟為：

（2）由SFEM 求s和?s／?X.

（3）求g和?g／?X，采用式（24），（25）.

（5）求β，采用式（21）.

每迭代一次求解β，就要解一次隨機(jī)有限元，求解的關(guān)鍵是計(jì)算?s／?X，可采用復(fù)數(shù)偏導(dǎo)法.

4 復(fù)數(shù)變量求偏導(dǎo)SFEMR法

由式（15）所示的f（x＋i h）的Taylor級(jí)數(shù)展式可得

忽略高階微小量O（h2），則有

對(duì)于FEM，在復(fù)數(shù)空間進(jìn)行運(yùn)算，可取復(fù)數(shù)計(jì)算結(jié)果的實(shí)部作為響應(yīng)量的值，取復(fù)數(shù)計(jì)算結(jié)果的虛部計(jì)算響應(yīng)量的偏導(dǎo)，無需再在實(shí)數(shù)空間進(jìn)行計(jì)算.

在3.2節(jié)的隨機(jī)有限元一次二階矩法的迭代格式中，用復(fù)數(shù)空間的FEM取代SFEM，迭代格式為：

（2）進(jìn)行復(fù)數(shù)空間FEM計(jì)算.取復(fù)數(shù)響應(yīng)量的實(shí)部作為s的值，采用式（27）；取復(fù)數(shù)響應(yīng)量的虛部求解?s／?X 的值，采用式（17）.

（3）求g和?g／?X，采用式（24），（25）.

（5）求β，采用式（21）.

和隨機(jī)有限元一次二階矩法相比，本節(jié)方法用復(fù)數(shù)空間的FEM計(jì)算代替了SFEM計(jì)算，其程序的編制只需普通的FEM程序，無需再編制SFEM程序，運(yùn)算過程大大簡化.

以算例1中的桁架結(jié)構(gòu)為例，設(shè)A，P1，P2分別為復(fù)數(shù)變量且隨機(jī)擾動(dòng)，根據(jù)式（27）計(jì)算得到的最大桿件應(yīng)力值σmax和實(shí)數(shù)空間的有限元結(jié)果比較見表2.

設(shè)桿件強(qiáng)度為R＝8.0×107MPa，A，P1，P2的變異系數(shù)為0.2，復(fù)數(shù)擾動(dòng)量取h＝10－30，計(jì)算危險(xiǎn)點(diǎn)的強(qiáng)度可靠度.功能函數(shù)為g＝R－σmax.根據(jù)第4節(jié)的迭代步驟計(jì)算可靠度系數(shù)β和驗(yàn)算點(diǎn)的值，迭代過程如表3所示.

表2表明，在復(fù)數(shù)空間運(yùn)用式（27）計(jì)算有限元響應(yīng)量的值具有很高的精度.表3表明，在求解β的迭代過程中，只需在復(fù)數(shù)空間進(jìn)行普通FEM計(jì)算，過程大大簡化.

表2 最大桿件應(yīng)力值σmax計(jì)算結(jié)果Tab.2 Results of the max stressσmaxin plane truss

表3 求解β迭代過程Tab.3 Iterative process forβ

5 結(jié)論

復(fù)數(shù)變量求導(dǎo)法效率高，精度高，應(yīng)用簡單方便，只需在復(fù)數(shù)空間進(jìn)行有限元計(jì)算，便可求出響應(yīng)量的偏導(dǎo)數(shù)，用在隨機(jī)有限元中，可求得響應(yīng)量的方差.與直接求偏導(dǎo)法相比，省去了矩陣的求逆和求導(dǎo)等大量運(yùn)算，計(jì)算過程大大簡化.

在隨機(jī)有限元一次二階矩的迭代格式中，取復(fù)數(shù)空間有限元的計(jì)算結(jié)果的實(shí)部作為響應(yīng)量的值，虛部用來求解響應(yīng)量的偏導(dǎo).這樣在求β的迭代過程中，只需進(jìn)行復(fù)數(shù)空間的有限元計(jì)算，無需再在實(shí)數(shù)空間進(jìn)行計(jì)算.

隨機(jī)有限元法可解決工程中的隨機(jī)結(jié)構(gòu)分析問題，是隨機(jī)結(jié)構(gòu)分析的有力工具.但其實(shí)現(xiàn)需要專門的程序，目前還沒有商業(yè)軟件可供利用.本文提出的復(fù)數(shù)變量法大大簡化了SFEM和SFEMR的實(shí)現(xiàn)過程，為工程應(yīng)用提供了一種現(xiàn)實(shí)可行的途徑.

［1］王建軍，于長波，李其漢.工程中的隨機(jī)有限元方法［J］.應(yīng)用力學(xué)學(xué)報(bào)，2009，26（2）：297.WANG Jianjun，YU Changbo，LI Qihan.Stochastic finite element methods in engineering［J］.Chinese Journal of Applied Mechanics，2009，26（2）：297.

［2］陳虬，劉先斌.隨機(jī)有限元法及其工程應(yīng)用［M］.成都：西南交通大學(xué)出版社，1993.CHEN Qiu，LIU Xianbin.Stochastic finite element methods and its application［M］.Chengdu：Press of Southwest Jiaotong University，1993.

［3］劉寧.可靠度隨機(jī)有限元法及其工程應(yīng)用［M］.北京：水利水電出版社，2001.LIU Ning.Reliability stochastic finite element methods and its application to engineering［M］.China Water Power Press，2001.

［4］Lyness J N，Moler C B.Numerical differentiation of analytic functions［J］.SIAM Journal on Numerical Analysis，1967（4）：202.

［5］Vatsa V N.Computation of sensitivity derivatives of Navier－Stokes equations using complex variables［J］.Advances in Engineering Software，2000，31：655.

［6］Daniel Contreras Mundstock，Rogério JoséMarczak.Boundary element sensitivity evaluation for elasticity problems using complex variable method［J］.Struct Multidisc Optim，2009，38：423.

［7］Gao X W，Liu D D，Chen P C.Internal stresses in inelastic BEM using complex－variable differentiation ［J］. Computational Mechanics，2002，28：40.

［8］張新培.建筑結(jié)構(gòu)可靠度分析與設(shè)計(jì)［M］.北京：科學(xué)出版社，2001.ZHANG Xinpei.Reliability analysis and design of building strictures［M］.Beijing：Science Press，2001.