基于Lyapunov特征指數(shù)的鋼制儲液罐動力失穩(wěn)概率分析

楊宏康， 高博青

(浙江大學(xué) 建筑工程學(xué)院，杭州　310058)

?

基于Lyapunov特征指數(shù)的鋼制儲液罐動力失穩(wěn)概率分析

楊宏康， 高博青

(浙江大學(xué) 建筑工程學(xué)院，杭州310058)

摘要：為有效量化儲液罐在地震激勵下的失穩(wěn)概率，參考震害報告選取536組三維地震波記錄，由壓力-位移格式的流固耦合模型建立等效動力擾動方程，然后通過計算動態(tài)Lyapunov特征指數(shù)確定儲液罐的動力失穩(wěn)概率。選取某10×104m3鋼制儲油罐作為分析對象，結(jié)果表明：水平向地震比豎向地震更易引起動力失穩(wěn)，多維地震比單維地震更加危險；動力失穩(wěn)概率隨抗風(fēng)圈的增加而減小；地震動方位、維數(shù)及抗風(fēng)圈設(shè)置對“失穩(wěn)概率-持時”曲線的影響較小。上述方法大幅降低了直接基于流固耦合模型確定儲液罐動力失穩(wěn)概率的計算成本，并能同時考慮地震動峰值與持時變化，從而全面直觀反映地震動三要素的影響。

關(guān)鍵詞：儲液罐; Lyapunov特征指數(shù); 動力失穩(wěn)概率; 多維地震; 持時

地震易損性表征結(jié)構(gòu)在特定強度地震作用下達到某一破壞狀態(tài)的條件失效概率，它是重要工程結(jié)構(gòu)地震概率安全評價的核心內(nèi)容[1]。基于美國國家標準技術(shù)研究所(NIST)的儲液罐震害報告[2]，O’Rourke[3]以及美國生命線聯(lián)盟(ALA)[4]對儲液罐的破壞現(xiàn)象進行了定量描述并劃分了儲罐損毀等級，通過完損信息的二次統(tǒng)計給出了儲液罐的地震易損性曲線及經(jīng)驗公式，但二次統(tǒng)計存在信息缺失、計量誤差等不確定性，儲罐損毀等級的確定也過于依賴主觀經(jīng)驗[5]。

儲液罐的震害事故有很大比例屬于屈曲破壞[2,4]。Iervolino等[6]采用響應(yīng)面法考慮設(shè)計參數(shù)變異性，通過質(zhì)量彈簧模型對儲液罐進行增量動力分析(簡稱IDA)，針對象足屈曲問題建立了易損性曲線。孫建剛等[7采用隨機變量概率模型，基于質(zhì)量彈簧模型就儲液罐的失穩(wěn)和提離等失效模式進行了易損性分析。文獻[6-7]采用的質(zhì)量彈簧模型需結(jié)合許用臨界壓應(yīng)力[σcr] =CEt/D進行罐殼穩(wěn)定核算，[σcr]主要由均布軸壓下的理想柱殼彈性穩(wěn)定理論確定[8]。大型鋼制儲罐常需設(shè)置環(huán)向抗風(fēng)圈來提高穩(wěn)定性能，該重要特征無法利用質(zhì)量彈簧模型來反映。儲液罐在地震下受到靜動液壓的復(fù)合作用，受力狀態(tài)與均布軸壓區(qū)別較大，GB50341、API 650、JIS B8501各取C= 0.15、0.413、0.33[8]，GB50341的新版征求意見稿已將C值調(diào)整為0.22。Berahman等[9]基于動力有限元分析確定儲罐傾覆彎矩與環(huán)向膜力，采用Bayesian方法針對象足屈曲破壞建立了易損性曲線。Buratti等[10通過附加質(zhì)量模型的IDA分析確定儲罐的臨界地震動峰值，據(jù)此定義破壞狀態(tài)并求解易損性曲線。IDA分析通過多次全過程時程分析來構(gòu)造偽平衡路徑，直接應(yīng)用至流固耦合模型的計算成本過高。受計算成本影響，文獻[9-10分別選用了12、14條地震波，地震動的復(fù)雜頻譜特性沒有得到全面反映。文獻[6-10]給出了儲罐易損性與水平地震動強度的關(guān)系，沒有考慮地震動多維特性和持時的影響。

本文參考ALA報告建立實際地震動數(shù)據(jù)庫，基于壓力-位移格式的流固耦合模型建立等效動力擾動方程并引入地震動多維特性，通過求解動態(tài)Lyapunov特征指數(shù)確定儲液罐的動力失穩(wěn)概率。本文同時考慮地震動峰值與持時變化，研究地震動多維特性與環(huán)向抗風(fēng)圈對大型鋼制儲液罐動力失穩(wěn)概率的影響。

1動力失穩(wěn)概率的求解方法

1.1動態(tài)Lyapunov特征指數(shù)

儲液罐為非線性定常流固耦合系統(tǒng)，基于勢流理論與標準波動方程[11，其動力擾動方程如下：

(1)

(2)

(3)

(4)

式中：Ma=Ms+Mu，Mu=ρfQKf-1QT；ügi(t)為地面加速度，i=1,2,3分別表示x,y,z方向；G為靜力作用的結(jié)構(gòu)域節(jié)點力列陣；當i= 1、2時，si=Maιsi，ιsi為地震影響向量；耦合作用沿兩域接觸面法向傳遞，當i= 3時，僅以-Maιs3üg3(t)引入豎向地震作用并不完整，取si=Maιsi-G/g，本文統(tǒng)一取g= 10 m/s2。

當大位移效應(yīng)不顯著且材料保持彈性時，采用線彈性剛度矩陣Ke替代Ks，并以幾何剛度矩陣K0、Kg(t)分別考慮G的預(yù)應(yīng)力效應(yīng)與時變應(yīng)力軟(剛)化效應(yīng)：

(5)

(6)

基于Newmark法離散化式(6)，記yk=y(tk)，tk=t0+kτ，t0為初始擾動時刻，τ為時間步長，則有：

(7)

(8)

Le表示微小擾動的平均指數(shù)變化率：當Le>0時，系統(tǒng)發(fā)生動力失穩(wěn)；當Le=0時，系統(tǒng)處于臨界狀態(tài)；當Le<0時，系統(tǒng)保持穩(wěn)定。結(jié)構(gòu)動力失穩(wěn)的地震動臨界峰值加速度acr及初始動力失穩(wěn)時刻tins可定義為[13]：

acr=a(Lm=0),tins=t(Le=Lm=0)

(9)

式中：a表示地震動峰值加速度變量，Lm= max(Le)。

1.2動力失穩(wěn)概率

借鑒地震易損性概念，動力失穩(wěn)概率可定義為[1]：

Pins=P[ INS | EQ(td,a)]

(10)

式中:“INS”表示動力失穩(wěn)的定量標準；EQ(td,a)表示地震動記錄樣本，td為持時變量。選取足量地震動記錄樣本進行多次動力穩(wěn)定性分析，Pins可表達為[14]：

Pins= num[EQ(td≥tins,a≥acr)] / num(EQ) (11)

式中：EQ(·)表示滿足特定條件的地震動記錄樣本，num(·)表示樣本數(shù)量。地震動記錄的實際持時長短不一，離散性較大，式(11)是按儲液罐動力失穩(wěn)的充分條件來確定的，即同時滿足td≥tins，a≥acr。

地震動的頻譜特性無法作為控制變量，Pins可理解為以地震動參數(shù)td、a作為變量的二元函數(shù)，Pins-td-a曲面可稱為動力失穩(wěn)概率曲面。當td=tins(或a=acr)時，Pins-a(或Pins-td)曲線即為動力失穩(wěn)概率曲線。

2儲液罐的動力失穩(wěn)概率分析

2.1模型參數(shù)及模態(tài)分析

本文選用某10×104m3儲油罐(見圖1)作為分析對象[13]：罐壁共9層，t1～t9為壁厚，鋼材彈性模量為206 GPa；儲液密度為860 kg/m3，聲速取1 400 m/s，液面高度與第9層壁板中部平齊；包邊角鋼位于罐壁頂部，1 ～ 4號環(huán)形抗風(fēng)圈均由槽鋼與鋼板焊接而成，依次位于第9、8、7層壁板中部與第5層壁板頂部。

采用有限元軟件ANSYS建立壓力-位移格式的流固耦合模型，通過科學(xué)計算軟件MATLAB完成壓力變量凝聚及模態(tài)分析，而后組集動力擾動方程并編制基于Lyapunov特征指數(shù)的動力失穩(wěn)概率分析模塊。罐壁和儲液分別采用Shell 181和Fluid 30單元模擬，罐壁網(wǎng)格沿周向均分為120段，每層罐壁沿豎向均分為2段，流體網(wǎng)格與之協(xié)調(diào)；抗風(fēng)圈鋼板采用Shell 181單元模擬，槽鋼與角鋼以Beam 188單元模擬，網(wǎng)格與罐壁協(xié)調(diào)。圖2所示模型滿足自錨固約束[13]，本文通過約束底部罐壁的所有位移自由度來模擬邊界條件。

圖1　10×104m3儲油罐的幾何模型Fig.1 Geometric model of 10×104 m3 oil storage tank

圖2　10×104 m3儲油罐的有限元模型Fig.2 Finite element model of 10×104 m3 oil storage tank

mnfn/HzM*ni×106/kgi=1i=2i=3i=1i=2i=3i=1i=2i=3Ⅰ2019212.22.22.222.922.942.4Ⅱ1920602.22.23.92.52.56.0Ⅲ9798995.35.35.30.30.32.0

儲液罐的模態(tài)分析結(jié)果如表1所示，fn為模態(tài)頻率，m(=Ⅰ,Ⅱ,Ⅲ)按Mni*升序排列。在相同i下，表1所列模態(tài)的Mni*所占比重均在95%以上，因此無需取全部模態(tài)來確定Kg(t)。由表1可知，在三維地震動作用下，考慮階次n的重復(fù)情況，僅需取7階不同模態(tài)(n= 19, 20, 21, 60, 97, 98, 99)即可確定Kg(t)。

2.2 Lyapunov特征指數(shù)方法的驗證

由下至上依次遞減儲罐的抗風(fēng)圈，考慮水平x向地震動作用，基于Lyapunov特征指數(shù)方法(簡稱LCE方法)對圖1所示模型進行動力穩(wěn)定性分析，并采用基于B-R準則的增量動力分析方法(IDA)進行驗證。其中，模態(tài)阻尼比取0.02，模態(tài)截斷數(shù)r= 100。Cape Mendocino與Kobe地震波的詳細信息見文獻[13]，其卓越頻率均在2 Hz左右。由圖3可知，LCE方法與IDA方法確定的acr值基本一致，前者相對后者的最大誤差僅在10%左右，詳細數(shù)據(jù)和論述見文獻[13]。

圖3　臨界地震動峰值加速度αcr的對比Fig.3 Comparison of critical peak ground acceleration αcr

2.3地震波特性分析

根據(jù)美國生命線聯(lián)盟(ALA)的儲液罐震害調(diào)查報告，結(jié)合地震動數(shù)據(jù)庫PEER、COSMOS的記錄情況，綜合確定表2所示的536組地震動信息，并進行傅里葉分析與反應(yīng)譜分析(阻尼比取0.02)，地震動卓越頻率fp與標準動力放大系數(shù)β譜如圖4所示，fp表示地震波最大幅值分量所對應(yīng)的頻率，β表示譜加速度與地震動峰值加速度的比值。地震動方位統(tǒng)一采用笛卡爾坐標系標記，x向與東西向(或較小記錄角)一致，y向與南北向(或較大記錄角)一致，z向?qū)?yīng)豎直分量。由圖4可知：表2的地震動記錄頻譜覆蓋范圍廣，反應(yīng)譜特性豐富；x,y,z向地震動的fp分別集中于1～6 Hz、1～6 Hz、1～10 Hz范圍內(nèi)，最大β值分別為3.20、3.18、3.26；x,y向地震動的fp分布及β均值譜均十分相似。考慮到x,y向地震動的fp相近，z向fp較高，后文以x向的fp值作為地震動記錄的主頻代表值。

表2　儲液罐的地震動數(shù)據(jù)庫

圖4　地震動記錄的頻譜特性Fig.4 Frequency spectrum characteristic of earthquake records

2.4地震動維數(shù)及方位的影響

地震動存在很大不確定性，具有多維、頻譜復(fù)雜等特點，本節(jié)按x,y,z,xy,xyz5類方位組合處理表2地震波，然后基于動態(tài)Lyapunov特征指數(shù)對儲液罐進行動力穩(wěn)定性分析，acr、tins以及動力失穩(wěn)概率曲面如圖5所示。當采用多維地震動時，acr對應(yīng)x向的地震動峰值加速度，y,z向的a值按原比例調(diào)整。

圖5　儲液罐的動力失穩(wěn)概率曲面Fig.5 Dynamic instability probability surfaces of liquid storage tank

由圖5可知：

(1)x,y向地震動下的acr值分布范圍較一致，多分布于30 m/s2以下，z向多分布于60 m/s2以下，xy,xyz向多分布于15 m/s2以下，xyz向的acr值更低一些。多維地震下的acr值比單維地震動下降低約50%～75%，z向地震動對x,y向地震動的影響有限。

(2)x,y向地震動下的動力失穩(wěn)概率曲面非常相似，z向失穩(wěn)概率曲面最為平緩，xy,xyz向的失穩(wěn)概率曲面最為陡峭，5類方位組合下的動力失穩(wěn)概率曲面在原點位置均存在明顯底角。

如圖6所示，當td=tins(或a=acr)時，圖5所示的動力失穩(wěn)概率曲面退化為曲線形式。由圖6(a)可知，當?shù)卣鸱逯导铀俣萢相同時，Pins滿足xyz>xy>x=y>z的大小關(guān)系，且z向地震動下的Pins-a曲線“坡度”較小，變化趨勢明顯有異于其余方位；從概率角度而言，5類方位下的Pins-td曲線基本一致，地震動方位及維數(shù)對Pins-td曲線的影響不大，相互差異不超過10%。

圖6　儲液罐的動力失穩(wěn)概率曲線Fig.6 Dynamic instability probability curves of liquid storage tank

2.5環(huán)向抗風(fēng)圈的影響

環(huán)向抗風(fēng)圈主要用于加強儲液罐在風(fēng)荷載作用下的穩(wěn)定性，具體設(shè)置道數(shù)和位置會因周邊氣候環(huán)境而異[15]。由下至上遞減儲液罐的抗風(fēng)圈，依次編號為M4、M3、M2、M1、M0，數(shù)字代表抗風(fēng)圈數(shù)量。模型M0～M4的動力失穩(wěn)概率曲線如圖7所示。我國戰(zhàn)略石油儲備基地的最大設(shè)防烈度為獨山子區(qū)的8度(0.20 g)，規(guī)范建議此時用于時程分析的罕遇地震加速度峰值為0.40 g，儲液罐的動力失穩(wěn)概率P0.2g、P0.4g分別如表3、4所示，下標表示相應(yīng)峰值水平。

表3　動力失穩(wěn)概率P0.2g

表4　動力失穩(wěn)概率P0.4g

圖7　抗風(fēng)圈對動力失穩(wěn)概率曲線的影響Fig.7 Effect of wind girders on dynamic instability probability curves

由圖7可知：

(1) 從Pins-a曲線來看，Pins隨抗風(fēng)圈道數(shù)的增加而顯著降低，而且Pins-a曲線在不同地震動方位組合下的收縮程度滿足z>x=y>xyz≈xy的遞進關(guān)系。上述結(jié)果表明，從緩解儲液罐地震動危險性的相對效果來看，抗風(fēng)圈在z向地震動下最為有效。

(2) 不同抗風(fēng)圈設(shè)置下的Pins-td曲線非常相近，最大差異不超過20%。按規(guī)范建議，時程分析的地震有效持時應(yīng)滿足5～10倍的結(jié)構(gòu)基本周期(約2.5～5 s)，但td= 5 s時的Pins僅在8%左右，故為保證儲液罐的動力穩(wěn)定性，建議持時td取10 s以上。

由表3、4可知：①Pins隨著抗風(fēng)圈的增加而減小，大體符合xyz>xy>x=y>z的大小關(guān)系；② 當a= 0.20 g時，儲罐M4在x,y向地震動下甚至未發(fā)生動力失穩(wěn)，但其在z,xy,xyz向地震動下仍有可能失穩(wěn)；③ 當a= 0.40 g時，儲罐M4在5類地震動下均有可能發(fā)生動力失穩(wěn)，其在xy,xyz向地震動下的失穩(wěn)概率可達10%。

本文采用4G內(nèi)存、2.33 GHz主頻的4核計算機，若直接基于壓力-位移格式的流固耦合模型進行增量動力分析，完成單次失穩(wěn)概率分析的時間將數(shù)以年計。本文經(jīng)凝聚壓力自由度并參數(shù)化結(jié)構(gòu)域抗力，藉由計算動態(tài)Lyapunov特征指數(shù)極大提高了動力失穩(wěn)概率的求解效率，完成單次失穩(wěn)概率分析在8小時左右。

3結(jié)論

(1) 本文建立的地震動數(shù)據(jù)庫包含536組3向地震波信息，數(shù)量充足，頻譜與反應(yīng)譜特性豐富，能合理反映地震動的不確定性及多維特點。

(2) 經(jīng)凝聚壓力變量并采用時變幾何剛度矩陣參數(shù)化結(jié)構(gòu)域抗力，Lyapunov特征指數(shù)法大幅提升了儲液罐動力失穩(wěn)概率分析的計算效率。相比直接基于流固耦合模型的增量動力分析方法，Lyapunov特征指數(shù)法可節(jié)約時間成本2 000倍左右。

(3) 豎向地震相對不易引起動力失穩(wěn)，多維地震下的失穩(wěn)概率比單維地震提高50%～75%，動力失穩(wěn)概率與地震動峰值與持時成正相關(guān)，地震動方位及維數(shù)導(dǎo)致的Pins-td曲線的最大差異不超過10%。

(4) 動力失穩(wěn)概率隨抗風(fēng)圈的增加而減小，抗風(fēng)圈對動力失穩(wěn)的抑制作用在豎向地震下最為顯著，Pins-td曲線在不同抗風(fēng)圈下的最大差別不到20%。

針對儲液罐的動力失穩(wěn)現(xiàn)象，設(shè)計分析時應(yīng)結(jié)合規(guī)范量化不同抗震設(shè)防水準下的動力失穩(wěn)概率容許值，從而有效評價儲液罐選址和設(shè)計選型的合理性。

參 考 文 獻

[1] Ellingwood B R. Earthquake risk assessment of building structures[J]. Reliability Engineering & System Safety, 2001, 74(3): 251-262.

[2] Cooper T W. A study of the performance of petroleum storage tanks during earthquakes, 1933-1995[R]. Gaithersburg, MD: US National Institute of Standards and Technology, 1997.

[3] O’Rourke M J. Seismic fragility curves for on-grade steel tanks[J]. Earthquake Spectra, 2000, 16(4): 801-815.

[4] Eidinger J M, Avila E A, Ballantyne D B, et al. Seismic fragility formulations for water systems[R]. Washington, DC: American Lifelines Alliance, 2001.

[5] Berahman F, Behnamfar F. Seismic fragility curves for un-anchored on-grade steel storage tanks: Bayesian approach[J]. Journal of Earthquake Engineering, 2007, 11(2): 166-192.

[6] Iervolino I, Fabbrocino G, Manfredi G. Fragility of standard industrial structures by a response surface based method[J]. Journal of Earthquake Engineering, 2004, 8(6): 927-945.

[7] 孫建剛, 張榮花, 蔣峰. 儲罐地震易損性數(shù)值仿真分析[J]. 哈爾濱工業(yè)大學(xué)學(xué)報, 2009, 41(12): 138-142.

SUN Jian-gang, ZHANG Rong-hua, JIANG Feng. Numerical simulation analysis on the seismic fragility of storage-tank[J]. Journal of Harbin Institute of Technology, 2009,41(12):138-142.

[8] Chen Z, Sun B, Yu C, et al. Comparison of the strength design and prevention method of elephant foot buckling among countries’ standards of oil tanks[C] // Proceedings of the ASME 2009 pressure vessels and piping division conference. Prague, Czech Republic: ASME, 2009: 461-466.

[9] Berahman F, Behnamfar F. Probabilistic seismic demand model and fragility estimates for critical failure modes of un-anchored steel storage tanks in petroleum complexes[J]. Probabilistic Engineering Mechanics, 2009, 24(4): 527-536.

[10] Buratti N, Tavano M. Dynamic buckling and seismic fragility of anchored steel tanks by the added mass method[J]. Earthquake Engineering & Structural Dynamics, 2014,43(1): 1-21.

[11] Moslemi M, Kianoush M R. Parametric study on dynamic behavior of cylindrical ground-supported tanks[J]. Engineering Structures, 2012, 42: 214-230.

[12] Wilson E L. Three-dimensional static and dynamic analysis of tructures[M]. 3rd ed. Berkeley, California: Computers and Structures, Inc, 2002: 160-180.

[13] 楊宏康, 高博青. 鋼制儲液罐的Lyapunov特征指數(shù)及彈性動力失穩(wěn)[J]. 上海交通大學(xué)學(xué)報,2014,48(11):1660-1666.

YANG Hong-kang, GAO Bo-qing. Lyapunov characteristic exponents and dynamic elastic instability of steel liquid storage tanks[J]. Journal of Shanghai Jiaotong University,2014,48(11):1660-1666.

[14] Zareian F, Krawinkler H. Assessment of probability of collapse and design for collapse safety[J]. Earthquake Engineering & Structural Dynamics, 2007, 36(13): 1901-1914.

[15] GB50341-2003. 立式圓筒形鋼制焊接油罐設(shè)計規(guī)范[S]. 北京: 中國計劃出版社, 2003.

基金項目：國家自然科學(xué)基金資助項目(51178414)

收稿日期：2014-04-18修改稿收到日期：2014-10-09

通信作者高博青 男，教授，博士生導(dǎo)師，1963年生

中圖分類號:TU 33; TE 972

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.01.019

Dynamic instability probability analysis for liquid storage steel tanks subjected to earthquake excitations

YANG Hong-kang, GAO Bo-qing

(College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China)

Abstract:To effectively quantify the instability probability of liquid storage steel tanks under earthquake excitations, 536 sets of three-dimensional seismic wave records were selected according to the earthquake damage reports, and the equivalent dynamic perturbation equations were then established based on a fluid-solid coupled model with pressure-displacement form, the dynamic instability probabilities of liquid storage tanks were determined by calculating Lyapunov characteristic exponents. A 10×104m3 oil storage steel tank was chosen as an analysis object, the results showed that horizontal earthquakes are more likely to cause dynamic instability than vertical ones be; multidimensional earthquakes are more dangerous than unidimensional ones be; the dynamic instability probability decreases with increase in wind girder; earthquake direction, dimension and wind girder have a slight influence on “instability probability - time duration” curves; the proposed method greatly reduces the computing cost of dynamic risk analysis of liquid storage tanks based on the fluid-solid coupled model, it can simultaneously consider changes of peak acceleration and time duration of ground motion, so it can comprehensively and intuitively reflect the effects of three essential factors of ground motions.

Key words:liquid storage tanks; Lyapunov characteristic exponents; dynamic instability probability; multidimensional earthquakes; time duration

第一作者 楊宏康 男，博士生，1986年生