張 衡，項(xiàng) 煦，劉宇浩，黃 斌，曾 磊

(1.長江大學(xué)城市建設(shè)學(xué)院，荊州 434023；2.武漢理工大學(xué)土木工程與建筑學(xué)院，武漢 430070)

結(jié)構(gòu)穩(wěn)定性是關(guān)系到結(jié)構(gòu)安全的主要問題之一[1?2]。傳統(tǒng)的結(jié)構(gòu)穩(wěn)定性分析，如復(fù)合板材[3]、混凝土灌注樁[4?5]、殼體結(jié)構(gòu)[6? 8]和管道[9 ? 10]等的彈性屈曲荷載求解，均是在確定性結(jié)構(gòu)系統(tǒng)的框架下進(jìn)行。然而現(xiàn)實(shí)工程中，材料性能、幾何尺寸以及邊界約束總是存在一定程度的不確定性，一些研究已注意到這些不確定性將對(duì)結(jié)構(gòu)的穩(wěn)定性產(chǎn)生不可忽視的影響[4,11?12]。因此，在結(jié)構(gòu)建模中考慮這些不確定性將更加符合工程實(shí)際，對(duì)結(jié)構(gòu)的安全評(píng)估具有重要意義。

為研究材料不確定性(如彈性模量的隨機(jī)性)對(duì)結(jié)構(gòu)屈曲荷載的影響，學(xué)者提出了不同的分析方法，主要包括蒙特卡洛模擬法[13?15]和基于攝動(dòng)技術(shù)的隨機(jī)有限元法[16?23]。

SHINOZUKA 和ASTILL [13]利用蒙特卡洛模擬法研究了支撐于彈性地基上梁柱構(gòu)件屈曲荷載的均值和方差；SONG 等[14]采用逆迭代法結(jié)合蒙特卡洛模擬實(shí)現(xiàn)了高斯隨機(jī)場(chǎng)下第一階屈曲特征值的求解，并指出當(dāng)?shù)谝浑A、二階屈曲特征值為密特征值時(shí)，由于屈曲模態(tài)在結(jié)構(gòu)參數(shù)變化后可能互換，需將兩個(gè)屈曲模態(tài)同時(shí)作為候選向量才能保證結(jié)果的準(zhǔn)確性。蒙特卡洛模擬法的應(yīng)用不受結(jié)構(gòu)規(guī)模、隨機(jī)變量變異性大小以及變量維度的影響，但其求解過程十分耗時(shí)。一些學(xué)者在攝動(dòng)法的基礎(chǔ)上提出了計(jì)算效率優(yōu)異的隨機(jī)有限元法，應(yīng)用于結(jié)構(gòu)的彈性穩(wěn)定性分析。

對(duì)于考慮材料不確定性的軸心受壓柱，JEONG[16]用攝動(dòng)隨機(jī)有限元法得到了屈曲荷載的均值和方差；ZHANG 和ELLINGWOOD [17]用低階攝動(dòng)法計(jì)算彈性地基上梁屈曲荷載的均值和方差，并指出隨著隨機(jī)場(chǎng)變異性的增大，二階攝動(dòng)求得屈曲荷載統(tǒng)計(jì)特征值精度將降低；RAMU 和GANESAN[18]考慮彈性模量和質(zhì)量密度的隨機(jī)性，研究了置于Winkler 地基模型上梁結(jié)構(gòu)的穩(wěn)定性；GRAHAM和 SIRAGY [19]對(duì)加筋板進(jìn)行了均方差為0.1 的小變異參數(shù)下的彈性穩(wěn)定性分析；ALTUS 和TOTRY[20]利用函數(shù)式攝動(dòng)法得到了考慮材料不確定時(shí)形式較簡單結(jié)構(gòu)屈曲荷載表達(dá)式；近期，GUPTA 和ARUN[21]將鋼材的彈性模量考慮為高斯分布或?qū)?shù)分布隨機(jī)場(chǎng)，結(jié)合無網(wǎng)格伽遼金法和攝動(dòng)技術(shù)求解了隨機(jī)參數(shù)變異性較小時(shí)軸心受壓柱隨機(jī)屈曲荷載的均值和方差?？紤]到低階攝動(dòng)法的計(jì)算精度只在隨機(jī)參數(shù)變異性較小時(shí)有保證，KAMI?SKI和?WITA[22]提出用廣義隨機(jī)有限元法(GSFEM)分析彈性穩(wěn)定性問題，該方法最高采用了十階攝動(dòng)技術(shù)求解臨界應(yīng)力的前四階統(tǒng)計(jì)矩。但該方法中，為便于推導(dǎo)隨機(jī)響應(yīng)級(jí)數(shù)表達(dá)式的高階項(xiàng)，要求隨機(jī)變量服從高斯分布且為單一變量。

除隨機(jī)有限元法外，學(xué)者們也提出了不確定彈性穩(wěn)定性分析的其它方法。如WU 等 [23 ? 24]將信息不完備的不確定性參數(shù)描述為區(qū)間量，給出了屈曲荷載的區(qū)間值，但各區(qū)間參數(shù)的波動(dòng)范圍不超過 ±24%；LI 等[25]指出對(duì)數(shù)正態(tài)分布比高斯分布更適合于描述彈性模量的隨機(jī)性，并提出用等幾何譜隨機(jī)有限元法來研究考慮材料不確定時(shí)板的彈性穩(wěn)定性，但隨機(jī)變量的變異系數(shù)小于0.1；在響應(yīng)面法的基礎(chǔ)上，ALIBRANDI 等[26]對(duì)隨機(jī)框架結(jié)構(gòu)的屈曲極限狀態(tài)進(jìn)行了可靠性分析；徐亞洲和白國良 [27]采用概率密度演化法研究了大型冷卻塔在自重、內(nèi)吸力及外部風(fēng)壓綜合作用下的屈曲承載能力，并進(jìn)行可靠性分析；SU 等[28]利用格林函數(shù)法得到了材料和幾何參數(shù)小變異時(shí)結(jié)構(gòu)的屈曲荷載?？偨Y(jié)來說，這些方法都是試圖對(duì)隨機(jī)屈曲荷載和屈曲模態(tài)進(jìn)行準(zhǔn)確而高效的求解。遺憾的是，除概率密度演化法之外，這些方法在求解隨機(jī)屈曲荷載時(shí)通常要求結(jié)構(gòu)參數(shù)變異性較小且為高斯分布，對(duì)于屈曲荷載的分析也局限于低階統(tǒng)計(jì)矩(均值和方差)。眾所周知，現(xiàn)實(shí)工程中結(jié)構(gòu)參數(shù)常為非高斯分布且具有較大變異性；同時(shí)，僅參考低階統(tǒng)計(jì)矩對(duì)于精確描述隨機(jī)響應(yīng)的概率分布特征是不夠的。特別是當(dāng)結(jié)構(gòu)失穩(wěn)為小概率事件時(shí)，屈曲荷載概率分布的尾部能否準(zhǔn)確擬合就顯得尤為重要。而利用高階統(tǒng)計(jì)矩(如偏度和峰度)可進(jìn)一步保證對(duì)隨機(jī)響應(yīng)概率分布(特別是尾部)擬合的準(zhǔn)確性。因此，建立一種不局限于隨機(jī)變量分布類型、能處理較大變異性并獲得高精度高階統(tǒng)計(jì)矩的隨機(jī)屈曲荷載高效求解方法具有重要意義。

近期文獻(xiàn)[29 ? 31]基于同倫分析法提出了同倫隨機(jī)有限元法(HSFEM)來處理線性和幾何非線性靜力響應(yīng)問題[29?30]及隨機(jī)特征值問題[31]。該方法繼承了同倫分析法[32?33]不受小參數(shù)限制的特點(diǎn)，當(dāng)結(jié)構(gòu)參數(shù)有較大變異性時(shí)仍具有良好的求解精度，且對(duì)隨機(jī)變量的分布類型沒有限制。與現(xiàn)有常用隨機(jī)有限元法如攝動(dòng)隨機(jī)有限元法和正交多項(xiàng)式展開法(PC 法)相比，在展開階數(shù)相同的情況下，HSFEM 比攝動(dòng)隨機(jī)有限元法具有更大的收斂域和更優(yōu)的求解精度，且有效避免了高階攝動(dòng)法可能出現(xiàn)的發(fā)散問題。同時(shí)，其計(jì)算工作量和求解難度顯著小于PC 法：設(shè)結(jié)構(gòu)自由度為N，含n個(gè)隨機(jī)變量，級(jí)數(shù)取r階展開，對(duì)于隨機(jī)特征值問題，利用伽遼金法確定PC 法的投影系數(shù)，將把原N×N維特征方程擴(kuò)階為初始值未知的N(n+r)!/n!r!維非線性方程組迭代求解問題，不易保證所求解的穩(wěn)定性。相比之下，HSFEM 中的系數(shù)則是通過1+rn(rn?r?n+5)/4 次N維矩陣的加、減和乘法運(yùn)算求得。

在同倫隨機(jī)有限元法中，結(jié)構(gòu)的隨機(jī)響應(yīng)首先以同倫級(jí)數(shù)表示，然后，通過在隨機(jī)樣本空間內(nèi)選取一定數(shù)量的樣本來確定同倫級(jí)數(shù)的最優(yōu)形式，因此，該方法的計(jì)算精度會(huì)受所選樣本點(diǎn)的影響，且對(duì)不同問題樣本點(diǎn)的選擇需依賴一定的人工經(jīng)驗(yàn)。故尋求一種最優(yōu)同倫級(jí)數(shù)的自適應(yīng)確定方法有利于大大提高該方法的穩(wěn)定性。

本文基于隨機(jī)殘差最小化思想改進(jìn)同倫隨機(jī)有限元法，并將其應(yīng)用于隨機(jī)結(jié)構(gòu)的彈性穩(wěn)定性分析。首先，推導(dǎo)了彈性結(jié)構(gòu)隨機(jī)屈曲荷載和屈曲模態(tài)的同倫級(jí)數(shù)表達(dá)式，并給出了同倫級(jí)數(shù)中任意階系數(shù)的顯式遞推表達(dá)。然后，定義了彈性穩(wěn)定方程解的隨機(jī)殘余誤差，通過最小化該隨機(jī)殘余誤差實(shí)現(xiàn)了同倫級(jí)數(shù)的優(yōu)化求解。本文方法解決了目前同倫隨機(jī)有限元法計(jì)算精度會(huì)受所選樣本點(diǎn)影響的問題，具有很好的穩(wěn)定性。同時(shí)，本文方法不受參數(shù)分布類型限制，適用于結(jié)構(gòu)參數(shù)變異性較大的情形，有效避免了高階攝動(dòng)法計(jì)算大變異隨機(jī)參數(shù)結(jié)構(gòu)屈曲荷載時(shí)容易出現(xiàn)發(fā)散的現(xiàn)象，并能獲得屈曲荷載及模態(tài)的高精度高階統(tǒng)計(jì)矩。通過函數(shù)算例以及變截面軸心受壓柱和平面框架結(jié)構(gòu)的彈性穩(wěn)定性分析，驗(yàn)證了本文方法的有效性。

1 結(jié)構(gòu)彈性穩(wěn)定性分析的有限元方程

不失一般性，以承受外荷載P的確定性梁結(jié)構(gòu)為例，設(shè)每個(gè)梁單元的內(nèi)力為Fe，則每個(gè)有限單元的勢(shì)能可寫為：

式中：E、I和le分別為彈性模量、截面慣性矩以及單元長度；we為單元位移場(chǎng)，是單元節(jié)點(diǎn)位移De與形狀函數(shù)Ne的乘積，可表示為we=NeDe，將它代入式(1)可得：

式中，λ 和D分別為屈曲特征值和特征向量。其中，λ 也稱為荷載比例因子，λ 中的最小特征值λcr對(duì)應(yīng)無缺陷結(jié)構(gòu)的屈曲荷載Pcr=λcrP。

2 求解結(jié)構(gòu)屈曲特征值的改進(jìn)同倫隨機(jī)有限元法

2.1 隨機(jī)屈曲特征對(duì)的多變量同倫級(jí)數(shù)表達(dá)

本文考慮材料彈性模量的不確定性，將其描述為隨機(jī)場(chǎng)或用相互獨(dú)立的隨機(jī)變量來表達(dá)，這樣式(3)中的彈性剛度矩陣K可用式(4)表示：

式中：K0 為彈性模量取均值時(shí)的整體剛度矩陣；ΔK 為彈性剛度矩陣的隨機(jī)部分，它包含n 個(gè)相互

獨(dú) 立 的 隨 機(jī) 變 量，其 向 量 形 式 為ξ={ξ1, ξ2, …,ξn}；Ki 為對(duì)應(yīng)第i 個(gè)隨機(jī)變量的確定性系數(shù)矩陣，可通過K-L 展開、中心點(diǎn)法、局部平均法等隨機(jī)場(chǎng)離散方法得到。那么，屈曲特征值λ 和屈曲特征向量D 將是隨機(jī)向量ξ 的函數(shù)，此時(shí)隨機(jī)結(jié)構(gòu)彈性穩(wěn)定性方程為：

通常采用攝動(dòng)法和正交多項(xiàng)式展開法對(duì)形如式(5)的含隨機(jī)變量的廣義特征方程進(jìn)行求解。攝動(dòng)法一般采用低階攝動(dòng)技術(shù)，因此，在隨機(jī)變量變異性較大時(shí)計(jì)算精度難以保證；正交多項(xiàng)式展開法采用低階展開表達(dá)式時(shí)的計(jì)算精度明顯優(yōu)于攝動(dòng)法，但是在隨機(jī)變量維度和結(jié)構(gòu)自由度增加時(shí)，方程的求解工作量會(huì)急劇增大，同時(shí)，因方程為迭代初始值未知的高維非線性方程組，求解十分困難[34]。

廖世俊[32]最早在博士論文中提出了同倫分析法用以處理強(qiáng)非線性問題，該方法較傳統(tǒng)的攝動(dòng)法有著顯著的優(yōu)越性，已在流體力學(xué)、固體力學(xué)、應(yīng)用數(shù)學(xué)和金融等領(lǐng)域得到廣泛的應(yīng)用[35]。本文基于同倫分析方法的思想，對(duì)式(5)進(jìn)行求解。下面給出具體求解過程。

對(duì)于式(5)，設(shè)隨機(jī)變量取均值時(shí)，對(duì)應(yīng)的解為λ0和D0，同倫分析方法通過重新構(gòu)造式(5)，建立起真實(shí)解λ(ξ)和D(ξ)與均值系統(tǒng)下λ0和D0之間的關(guān)系。首先將式(5)重新構(gòu)造如下：

式 中：p∈[0, 1]、h≠0 為 輔 助 參 數(shù)；Γ(ξ;h,p)和Θ(ξ;h,p)為隨機(jī)向量ξ 和輔助參數(shù)的函數(shù)。為簡化表達(dá)，引入運(yùn)算符號(hào)L和N，將式(6)寫為：

觀察式(7)可知，當(dāng)p=0 時(shí)，有：

當(dāng)p=1 時(shí)，有：

那么，伴隨著參數(shù)p從0 逐漸增加至1，Γ(ξ;h,p) 和 Θ(ξ;h,p)由λ0和D0變化到式(5)的真實(shí)解λ(ξ)和D(ξ)，這樣就建立了關(guān)于Γ (ξ;h,p) 和Θ(ξ;h,p)的同倫。將 Γ(ξ;h,p) 和 Θ(ξ;h,p)做關(guān)于參數(shù)p的泰勒展開，將有：

關(guān)于式(10)的收斂性，文獻(xiàn)[31]討論了一個(gè)變量的情況，并說明只要適當(dāng)選擇非線性方程中的線性算子，式(10)總會(huì)在p=1 處收斂；對(duì)于多個(gè)變量的情形，這個(gè)性質(zhì)仍然成立，具體證明可參考文獻(xiàn)[31]。那么，當(dāng)p=1 時(shí)有：

式中，λm和Dm(m≥1)可通過以下步驟得到：

1) 將式(7)對(duì)參數(shù)p求導(dǎo)m次，然后令p=0可得：

具體如下：

3) 將式(14)代入式(13)可得系數(shù)Dm的顯式遞推表達(dá)式：

觀察式(14)和式(15)可知，λm和Dm的遞推表達(dá)式是由已知矩陣K0、ΔK、Kg、λ0和D0等簡單的乘法與求和得到，且只涉及一次矩陣求逆(即η的計(jì)算)，故而可方便地通過MAPLE 等符號(hào)運(yùn)算軟件推出。將λm和Dm代入式(11)整理后，便可得到隨機(jī)屈曲特征值λ(ξ)和屈曲模態(tài)D(ξ)的表達(dá)式如下：

式中：Z(ξ;h) 為 屈曲特征值或屈曲模態(tài)；Zi1,Zi1i2,···等為確定性系數(shù)，其與λm和Dm之間有如下對(duì)應(yīng)關(guān)系：

式中：m=0, 1, …, ∞為式(16)的展開階數(shù)；βm,l(h)(l=1, 2, …,m)為只含輔助參數(shù)h的函數(shù)，稱為趨近函數(shù)，其中h∈(?2, 0)。

式(16)為同倫級(jí)數(shù)，該級(jí)數(shù)隨參數(shù)h∈(?2, 0)取值的不同而變化，相應(yīng)有不同的收斂范圍。可以證明，傳統(tǒng)的泰勒級(jí)數(shù)只是同倫級(jí)數(shù)的特解(即h=?1 時(shí))。同時(shí)，注意到推導(dǎo)過程無需指定隨機(jī)變量的分布類型，同倫級(jí)數(shù)可適用于任意分布類型的隨機(jī)變量。

2.2 同倫級(jí)數(shù)的降維與降階

理論上，當(dāng)同倫級(jí)數(shù)式(16)的展開階數(shù)m趨近無窮時(shí)，就是式(5)的精確解。在實(shí)際應(yīng)用中，為減少計(jì)算工作量，可只取有限階展開式；另外，考慮到當(dāng)變量增多時(shí)，高維交叉項(xiàng)對(duì)級(jí)數(shù)精度的貢獻(xiàn)將減小[29?30]，可忽略同倫級(jí)數(shù)中的部分交叉項(xiàng)以進(jìn)一步提高計(jì)算效率。

如果同倫級(jí)數(shù)式(16)中保留至多含q個(gè)不同隨機(jī)變量的交叉項(xiàng)，則稱其為q維同倫級(jí)數(shù)展開，通常一維或二維同倫級(jí)數(shù)展開即可獲得較滿意的計(jì)算精度[30]。表1 列出了一維、二維和完整同倫級(jí)數(shù)在不同隨機(jī)變量個(gè)數(shù)、不同展開階數(shù)時(shí)，需要計(jì)算系數(shù)的個(gè)數(shù)。

表1 同倫級(jí)數(shù)中需要計(jì)算系數(shù)的個(gè)數(shù)Table 1 Number of terms in homotopy series

從表1 可見，當(dāng)隨機(jī)變量個(gè)數(shù)或級(jí)數(shù)展開階數(shù)增加時(shí)，完整的同倫級(jí)數(shù)中，所需計(jì)算的系數(shù)個(gè)數(shù)將急劇增加；而采用一維或者二維同倫級(jí)數(shù)時(shí)，相應(yīng)的系數(shù)計(jì)算工作量將大幅減少。具體地，當(dāng)隨機(jī)變量數(shù)為20 時(shí)，3 階展開同倫級(jí)數(shù)中，一維系數(shù)個(gè)數(shù)約為完整級(jí)數(shù)的1/30，二維系數(shù)個(gè)數(shù)約為完整級(jí)數(shù)的1/3；若取4 階展開，則一維和二維系數(shù)個(gè)數(shù)將分別降為完整級(jí)數(shù)展開的1/130 和1/9。

當(dāng)通過級(jí)數(shù)截?cái)嗪徒稻S后，設(shè)同倫級(jí)數(shù)的近似解為Z?(ξ,h) ，該近似解與精確解Ze(ξ)之間的相對(duì)誤差為：

式中，|·|為絕對(duì)值運(yùn)算符?？梢杂^察到，該誤差會(huì)隨著h值的調(diào)整而變化。

2.3 確定輔助參數(shù)h 值的隨機(jī)殘差最小化法

為使近似解Z?(ξ,h)的誤差盡可能小，需要確定參數(shù)h的合適取值。在確定性問題中，同倫分析方法給出了三種確定h值的方法：h曲線法、殘差法和比值法[33,35]，這些方法均不便于直接用于隨機(jī)問題的求解。已有的HSFEM 方法采用抽樣技術(shù)來確定h值[29?31]，但該方法計(jì)算精度依賴于樣本值的選取，且對(duì)不同問題需要有一定的選點(diǎn)經(jīng)驗(yàn)。本文提出利用隨機(jī)殘差最小化法來實(shí)現(xiàn)輔助參數(shù)h值的自動(dòng)求解。

由于精確解Ze(ξ)未知，無法直接通過優(yōu)化式(19)來確定h值。注意到隨機(jī)彈性穩(wěn)定性方程式(5)是已知的，屈曲特征值和屈曲模態(tài)的精確解λe(ξ)與De(ξ)對(duì)隨機(jī)變量樣本空間內(nèi)任意一點(diǎn)，有式(20)恒成立：

式中，N為結(jié)構(gòu)自由度。第j(j=1, 2, …,N)階屈曲特征值與特征向量和在整個(gè)樣本區(qū)間( ξ ∈?)上滿足：

式中，gξ(ξ)為聯(lián)合概率密度函數(shù)。將經(jīng)過2.2 節(jié)降維與降階后的屈曲特征值和特征向量同倫級(jí)數(shù)代入式(21)，則有：

式中：fξi(ξi) (i=1, 2, …,n)為 ξi的概率密度函數(shù)；E{?}為求期望運(yùn)算符；wi(h)(i=1, 2, …,N)為關(guān)于參數(shù)h的標(biāo)量函數(shù)。

式(22)中的W(h)描述了在隨機(jī)變量全域上第j階屈曲特征對(duì)近似解基于彈性穩(wěn)定性方程產(chǎn)生的隨機(jī)殘差，是關(guān)于參數(shù)h的向量函數(shù)。注意到W(h)越小意味著近似效果越好，于是有如下優(yōu)化問題：

這里推薦最小化W(h)的2 范數(shù)來確定h值，即：

因式(24)僅含一個(gè)未知量h，且任一wi(h)均為多項(xiàng)式形式，故式(24)的求解工作量不大。

通過以上隨機(jī)殘差最小化法得到參數(shù)h的優(yōu)化值h*后，便可確定隨機(jī)屈曲特征值和屈曲模態(tài)的優(yōu)化同倫級(jí)數(shù)展開，稱此方法為改進(jìn)的同倫隨機(jī)有限元法。同時(shí)約定，經(jīng)由2.2 節(jié)降維的同倫級(jí)數(shù)中，采用只含單個(gè)隨機(jī)變量時(shí)的方法稱為本文方法-1，至多含兩個(gè)隨機(jī)變量交叉項(xiàng)時(shí)的方法稱為本文方法-2。

3 算例分析

算例分析中將采用不同方法擬合非線性函數(shù)，求解軸心受壓柱和平面框架結(jié)構(gòu)的屈曲特征值及屈曲模態(tài)，以說明本文方法(包括本文方法-1 和本文方法-2)的有效性。其它方法包括泰勒級(jí)數(shù)展開法、正交多項(xiàng)式展開法(PC 法)，文獻(xiàn)[17]提出的低階攝動(dòng)隨機(jī)有限元法(以下簡記為文獻(xiàn)[17]方法)，文獻(xiàn)[22]提出的廣義隨機(jī)有限元法(以下簡記為文獻(xiàn)[22]方法)以及基于樣本點(diǎn)的同倫隨機(jī)有限元法HSFEM，并以蒙特卡洛模擬法的結(jié)果作為參考值。所有方法均采用MATLAB 編程計(jì)算。

3.1 算例1.強(qiáng)非線性函數(shù)

函數(shù)Y如式(25)所示，該函數(shù)具有強(qiáng)非線性。分別采用泰勒級(jí)數(shù)展開法、正交多項(xiàng)式展開法(PC 法)和同倫級(jí)數(shù)展開法逼近。各方法均取4 階展開，得到的表達(dá)式見式(26)～式(28)：

分別采用HSFEM 方法和本文方法確定同倫級(jí)數(shù)展開中參數(shù)h的值。在HSFEM 方法中的樣本點(diǎn)坐標(biāo)分別取為1σ、2σ 和3σ(σ=1 為隨機(jī)變量的均方差)，所得h值分別為?0.96、?0.53 和?0.90，相應(yīng)計(jì)算結(jié)果以HSFEM(X)表示，其中X為樣本點(diǎn)坐標(biāo)；本文方法利用隨機(jī)殘差最小化法得到h值為?0.86。另以蒙特卡洛模擬法20 萬樣本的計(jì)算結(jié)果作為標(biāo)準(zhǔn)參考值。不同方法計(jì)算Y的前4 階統(tǒng)計(jì)矩結(jié)果見表2。

表2 不同方法計(jì)算Y 的前4 階統(tǒng)計(jì)矩Table 2 First four order statistical moments of Y calculated by different methods

從表2 中可以看出，對(duì)于均值和均方差，四種方法均能得到滿意的計(jì)算結(jié)果；對(duì)于偏度和峰度，泰勒級(jí)數(shù)結(jié)果已經(jīng)發(fā)散，PC 法精度較好；HSFEM 方法的結(jié)果較泰勒級(jí)數(shù)展開法有明顯改進(jìn)，但是精度受所選樣本點(diǎn)影響較大；本文方法的結(jié)果精度最高，與蒙特卡洛模擬法相比峰度誤差僅為8%。

為進(jìn)一步說明本文方法較HSFEM 方法的優(yōu)勢(shì)，繪制了隨機(jī)殘余誤差式(24)隨h值變化的關(guān)系曲線如圖1 所示。圖1 中，h=?1 處即為泰勒級(jí)數(shù)展開對(duì)應(yīng)的隨機(jī)殘余誤差?？捎^察到，HSFEM方法和本文方法的誤差較泰勒級(jí)數(shù)展開法有明顯減小，但HSFEM 方法所選樣本點(diǎn)并不能保證找到h值的最優(yōu)點(diǎn)，且樣本點(diǎn)坐標(biāo)的選取依賴于經(jīng)驗(yàn)，本文方法則能自動(dòng)找到全域上的誤差最小值。

圖1 隨機(jī)殘余誤差關(guān)于輔助參數(shù)h 值變化曲線Fig.1 Stochastic residual error with respect to the auxiliary parameter h

3.2 算例2.變截面軸心受壓柱

如圖2 所示變截面軸心受壓柱，AB和BC兩段長均為3 m，在柱的端部有一集中力荷載P=1 kN。柱的有限元模型被劃分為12 個(gè)單元，每個(gè)單元節(jié)點(diǎn)含橫向側(cè)移和轉(zhuǎn)角兩個(gè)自由度。下面分兩種工況討論AB段和BC段抗彎剛度隨機(jī)性對(duì)該柱屈曲特征值的影響。

圖2 變截面柱 /mFig.2 Variable cross-section column

工況1：假定AB段抗彎剛度為確定值EIAB=16 kN·m2，考慮BC段抗彎剛度的隨機(jī)性并以式(29)表達(dá)：

圖3 變截面柱屈曲特征值前4 階統(tǒng)計(jì)矩(工況1)Fig.3 First four order statistical moments of the buckling eigenvalue for Case 1

從圖3 可觀察到，由于本工況中抗彎剛度變異性較大，隨著展開階數(shù)的增加，文獻(xiàn)[22]方法只是對(duì)于均值存在收斂的趨勢(shì)，對(duì)于均方差、偏度和峰度，隨著展開階數(shù)的增加計(jì)算結(jié)果反而會(huì)發(fā)散；相比之下，基于同倫級(jí)數(shù)展開的HSFEM(7σ)和本文方法得到的各階統(tǒng)計(jì)矩值在4 階展開后均趨于穩(wěn)定；對(duì)于高階統(tǒng)計(jì)矩，如峰度，HSFEM(7σ)在8 階展開后可收斂，而本文方法只需要4 階展開即可，說明本文方法在計(jì)算隨機(jī)彈性屈曲特征值時(shí)具有良好的收斂性。

工況2：考慮多維隨機(jī)變量，AB段的抗彎剛度為確定值EIAB=16 kN·m2。BC段每個(gè)單元抗彎剛度為對(duì)數(shù)正態(tài)分布并以式(30)表示：

圖4 變截面柱屈曲特征值前4 階統(tǒng)計(jì)矩(工況2)Fig.4 First four order statistical moments of the buckling eigenvalue in Case 2

圖4 表明：與蒙特卡洛模擬法的結(jié)果相比，文獻(xiàn)[17]方法在計(jì)算均值和均方差時(shí)有較好的精度，但是對(duì)于偏度和峰度在 δ大于0.3 后會(huì)有顯著的誤差；如果采用高階攝動(dòng)(4 階攝動(dòng))，其求解精度不僅沒有改進(jìn)，反而在 δ大于0.5 后迅速發(fā)散?；谕瑐惣?jí)數(shù)展開的HSFEM(7σ)、本文方法-1 和本文方法-2 均能得到更加準(zhǔn)確的結(jié)果，特別是對(duì)偏度和峰度這些高階矩信息的求解沒有出現(xiàn)發(fā)散現(xiàn)象。同時(shí)，本文方法計(jì)算的偏度值比HSFEM(7σ)更好。

本文方法中h值和趨近函數(shù) βm,l(h) 隨 δ變化的關(guān)系如圖5 所示。從圖5(a)可見，隨著隨機(jī)變量變異性的增加，h值由接近?1 逐漸變化至?0.6，HSFEM(7σ)中確定的h值在隨機(jī)殘差法得到h值附近波動(dòng)，說明對(duì)于變異性不同的隨機(jī)變量，需采用不同的樣本點(diǎn)才能取得最優(yōu)h值；圖5(b)中，βm,l(h)、sigβm,l(h)和dou βm,l(h) (m=4,l=1, 2, 3, 4)分別為HSFEM(7σ)、本文方法-1 和本文方法-2 同倫級(jí)數(shù)中的趨近函數(shù)系數(shù)。可見，隨著隨機(jī)變量變異性的增大，3 階、4 階趨近函數(shù)值逐漸接近于0，可以顯著避免高階展開系數(shù)帶來的發(fā)散現(xiàn)象。

圖5 輔助參數(shù)h 值和趨近函數(shù)隨δ 變化關(guān)系圖(工況2)Fig.5 Values of the auxiliary parameter h and the approaching function in the homotopy series expansion for Case 2

3.3 算例3.平面框架結(jié)構(gòu)

如圖6 所示7 層平面框架結(jié)構(gòu)，每根框架柱受P=150 kN 軸向荷載。該框架有限元建模中，每隔1 m 將梁柱劃分為一個(gè)梁單元，每單元3 個(gè)自由度，模型共有383 個(gè)單元，自由度為1059 個(gè)。首先進(jìn)行內(nèi)力分析計(jì)算圖示荷載作用下各單元應(yīng)力，進(jìn)而確定各單元幾何剛度矩陣。

圖6 7 層框架結(jié)構(gòu) /mFig.6 A 7-story frame structure

1)考慮6 根框架柱彈性模量的隨機(jī)性，形成6個(gè)相互獨(dú)立的對(duì)數(shù)分布隨機(jī)場(chǎng)，每個(gè)隨機(jī)場(chǎng)表示為：

式中，ξi為第i個(gè)標(biāo)準(zhǔn)高斯分布隨機(jī)變量，指數(shù)部分是協(xié)方差核滿足式(32)的高斯隨機(jī)場(chǎng)：

式中：σ為高斯隨機(jī)場(chǎng)的均方差；y1和y2為沿Y方向的坐標(biāo)值；相關(guān)長度l=10 m。6 根柱子彈性模量隨機(jī)場(chǎng)中的μc均為 5.327。為將該對(duì)數(shù)分布隨機(jī)場(chǎng)轉(zhuǎn)化為級(jí)數(shù)展開的形式，可用正交多項(xiàng)式展開法(PC 法)來逼近式(31)，具體表達(dá)式如下：

式 中：ξ={ξ1,ξ2,···,ξ4}，gk(ξ)為 關(guān) 于 隨 機(jī) 變 量ξi(i=1,2,···,4 )的正交基；lk(x)為各正交基相應(yīng)的確定性系數(shù)，其求解方法可參見文獻(xiàn)[36]。

框架柱隨機(jī)場(chǎng)的變異系數(shù) δEc分別假設(shè)為 0.30和0.60。表3 和表4 分別列出了用攝動(dòng)法、HSFEM和本文方法采用2 階～4 階展開時(shí)求出的屈曲特征值前4 階統(tǒng)計(jì)矩值，其中HSFEM 選用樣本點(diǎn)坐標(biāo)為7 σ，并以蒙特卡洛模擬法(10 萬樣本)計(jì)算結(jié)果作為參考值。

表3 框架結(jié)構(gòu)屈曲特征值的前4 階統(tǒng)計(jì)矩( δEc=0.3)Table 3 First four order statistical moments of the buckling eigenvalue ( δEc=0.3)

表4 框架結(jié)構(gòu)屈曲特征值的前4 階統(tǒng)計(jì)矩( δEc=0.6)Table 4 First four order statistical moments of the buckling eigenvalue ( δEc=0.6)

從表3 可看出，在隨機(jī)場(chǎng)變異性較小的情況下，攝動(dòng)法、HSFEM 和本文方法求得的4 階矩值都與蒙特卡洛模擬方法很吻合，隨著展開階數(shù)的提高，計(jì)算精度也越高，且HSFEM 和本文方法計(jì)算的偏度和峰度略優(yōu)于攝動(dòng)法。

當(dāng)隨機(jī)場(chǎng)變異性增大時(shí)，參見表4，可發(fā)現(xiàn)2 階展開的攝動(dòng)法計(jì)算的峰度和偏度相對(duì)誤差超過20%。若提高展開階數(shù)，攝動(dòng)法的計(jì)算精度并未提高，反而出現(xiàn)了發(fā)散現(xiàn)象，說明高階攝動(dòng)并不一定適用于變異性較大時(shí)的非高斯分布隨機(jī)變量。HSFEM 的結(jié)果較攝動(dòng)法雖然有了明顯的改進(jìn)，但偏度和峰度的結(jié)果仍有較為顯著的誤差(特別是3 階展開時(shí))，原因是本算例中樣本點(diǎn)坐標(biāo)在隨機(jī)變量變異性不同時(shí)、或展開階數(shù)不同時(shí)均統(tǒng)一選取為7 σ。換言之，HSFEM 需針對(duì)不同變異性大小和展開階數(shù)提出更為具體的選點(diǎn)原則才能保證該方法計(jì)算精度的穩(wěn)定性。同時(shí)，可發(fā)現(xiàn)在隨機(jī)變量變異性較大時(shí)，HSFEM 精度對(duì)所選樣本將更加敏感。

而對(duì)于本文提出的兩種降維方法，隨著展開階數(shù)的增加，仍能保證計(jì)算精度的提高，說明本文方法具有良好的收斂性和很好的穩(wěn)定性。

圖7 給出了框架柱隨機(jī)場(chǎng)變異系數(shù)0.6 時(shí)，由4 階展開本文方法-1 計(jì)算的框架結(jié)構(gòu)第1 階屈曲模態(tài)的均值、均方差、偏度和峰度。從圖7 中可觀察到，對(duì)屈曲模態(tài)的前4 階統(tǒng)計(jì)矩，本文方法的計(jì)算精度都很高。

圖7 一階屈曲模態(tài)前4 階統(tǒng)計(jì)矩Fig.7 First four order statistical moments of the 1st-order buckling mode

表5 是攝動(dòng)法和本文方法采用不同展開階數(shù)時(shí)計(jì)算屈曲特征值所需要的時(shí)間。本算例中共考慮了24 個(gè)隨機(jī)變量，由于本文方法采取降維策略舍棄了交叉項(xiàng)，在高階展開時(shí)其計(jì)算量將少于同階展開的攝動(dòng)隨機(jī)有限元法，同時(shí)遠(yuǎn)少于蒙特卡洛模擬法的計(jì)算量。而HSFEM 方法需要額外進(jìn)行24 次確定性有限元分析以確定h值，此部分工作量需耗時(shí)約170 s。

表5 各方法計(jì)算時(shí)間 /sTable 5 Computation time for different methods

2)除考慮6 根框架柱彈性模量的隨機(jī)性外，同時(shí)考慮橫梁彈性模量的隨機(jī)性。設(shè)該7 層框架每層橫梁的彈性模量可表示為式(34)的形式：

表6 框架結(jié)構(gòu)屈曲特征值的前4 階統(tǒng)計(jì)矩( δEc=0.6, δb=0.3)Table 6 First four order statistical moments of the buckling eigenvalue ( δEc=0.6, δb=0.3)

4 結(jié)論

本文基于同倫分析法推導(dǎo)了屈曲特征值和屈曲模態(tài)的多變量同倫級(jí)數(shù)表達(dá)式，并給出級(jí)數(shù)中任意階系數(shù)的顯式遞推關(guān)系式；借助關(guān)于彈性穩(wěn)定性方程的隨機(jī)殘差最小化法，確定了隨機(jī)屈曲特征值和屈曲模態(tài)同倫級(jí)數(shù)解的最優(yōu)形式。

函數(shù)數(shù)值算例表明：本文方法能夠?qū)崿F(xiàn)隨機(jī)響應(yīng)同倫級(jí)數(shù)表達(dá)式中輔助參數(shù)h值的自動(dòng)尋優(yōu)，解決了已有同倫隨機(jī)有限元法HSFEM 結(jié)果易受所選樣本的影響和不穩(wěn)定的問題。并且本文方法結(jié)果比同階泰勒級(jí)數(shù)展開結(jié)果精度更高。

變截面軸心受壓桿和框架結(jié)構(gòu)的穩(wěn)定性分析數(shù)值算例表明：相比于現(xiàn)有的各類低階或高階攝動(dòng)法，本文方法求解的隨機(jī)屈曲特征值和屈曲模態(tài)均展現(xiàn)出優(yōu)異的計(jì)算精度，相比于蒙特卡洛模擬法計(jì)算效率大大提高。同時(shí)發(fā)現(xiàn)，當(dāng)隨機(jī)參數(shù)變異性較大時(shí)，利用高階攝動(dòng)法得到的高階統(tǒng)計(jì)矩(如偏度和峰度)可能比低階攝動(dòng)更差，而本文方法十分有效地避免了高階展開級(jí)數(shù)解可能出現(xiàn)的發(fā)散現(xiàn)象，體現(xiàn)了本文方法的穩(wěn)定性。此外，對(duì)于平面框架結(jié)構(gòu)，當(dāng)同時(shí)考慮柱和橫梁彈性模量的隨機(jī)性時(shí)，本文方法中交叉項(xiàng)對(duì)屈曲特征值計(jì)算精度的貢獻(xiàn)已不可忽略，建議采用至少含兩個(gè)隨機(jī)變量交叉項(xiàng)的同倫級(jí)數(shù)求解屈曲特征值的高階統(tǒng)計(jì)矩。

本文方法的提出為隨機(jī)結(jié)構(gòu)的彈性穩(wěn)定性分析提供了一種新的高效求解途徑。考慮到現(xiàn)有同倫隨機(jī)有限元法已被應(yīng)用于幾何非線性問題的求解，本文方法也可方便拓展到處理隨機(jī)結(jié)構(gòu)非線性靜力及幾何非線性穩(wěn)定問題。