基于Lyapunov特征指数的钢制储液罐动力失稳概率分析

杨宏康， 高博青

(浙江大学 建筑工程学院，杭州　310058)



基于Lyapunov特征指数的钢制储液罐动力失稳概率分析

杨宏康， 高博青

(浙江大学 建筑工程学院，杭州310058)

摘要：为有效量化储液罐在地震激励下的失稳概率，参考震害报告选取536组三维地震波记录，由压力-位移格式的流固耦合模型建立等效动力扰动方程，然后通过计算动态Lyapunov特征指数确定储液罐的动力失稳概率。选取某10×104m3钢制储油罐作为分析对象，结果表明：水平向地震比竖向地震更易引起动力失稳，多维地震比单维地震更加危险；动力失稳概率随抗风圈的增加而减小；地震动方位、维数及抗风圈设置对“失稳概率-持时”曲线的影响较小。上述方法大幅降低了直接基于流固耦合模型确定储液罐动力失稳概率的计算成本，并能同时考虑地震动峰值与持时变化，从而全面直观反映地震动三要素的影响。

关键词：储液罐; Lyapunov特征指数; 动力失稳概率; 多维地震; 持时

地震易损性表征结构在特定强度地震作用下达到某一破坏状态的条件失效概率，它是重要工程结构地震概率安全评价的核心内容[1]。基于美国国家标准技术研究所(NIST)的储液罐震害报告[2]，O’Rourke[3]以及美国生命线联盟(ALA)[4]对储液罐的破坏现象进行了定量描述并划分了储罐损毁等级，通过完损信息的二次统计给出了储液罐的地震易损性曲线及经验公式，但二次统计存在信息缺失、计量误差等不确定性，储罐损毁等级的确定也过于依赖主观经验[5]。

储液罐的震害事故有很大比例属于屈曲破坏[2,4]。Iervolino等[6]采用响应面法考虑设计参数变异性，通过质量弹簧模型对储液罐进行增量动力分析(简称IDA)，针对象足屈曲问题建立了易损性曲线。孙建刚等[7采用随机变量概率模型，基于质量弹簧模型就储液罐的失稳和提离等失效模式进行了易损性分析。文献[6-7]采用的质量弹簧模型需结合许用临界压应力[σcr] =CEt/D进行罐壳稳定核算，[σcr]主要由均布轴压下的理想柱壳弹性稳定理论确定[8]。大型钢制储罐常需设置环向抗风圈来提高稳定性能，该重要特征无法利用质量弹簧模型来反映。储液罐在地震下受到静动液压的复合作用，受力状态与均布轴压区别较大，GB50341、API 650、JIS B8501各取C= 0.15、0.413、0.33[8]，GB50341的新版征求意见稿已将C值调整为0.22。Berahman等[9]基于动力有限元分析确定储罐倾覆弯矩与环向膜力，采用Bayesian方法针对象足屈曲破坏建立了易损性曲线。Buratti等[10通过附加质量模型的IDA分析确定储罐的临界地震动峰值，据此定义破坏状态并求解易损性曲线。IDA分析通过多次全过程时程分析来构造伪平衡路径，直接应用至流固耦合模型的计算成本过高。受计算成本影响，文献[9-10分别选用了12、14条地震波，地震动的复杂频谱特性没有得到全面反映。文献[6-10]给出了储罐易损性与水平地震动强度的关系，没有考虑地震动多维特性和持时的影响。

本文参考ALA报告建立实际地震动数据库，基于压力-位移格式的流固耦合模型建立等效动力扰动方程并引入地震动多维特性，通过求解动态Lyapunov特征指数确定储液罐的动力失稳概率。本文同时考虑地震动峰值与持时变化，研究地震动多维特性与环向抗风圈对大型钢制储液罐动力失稳概率的影响。

1动力失稳概率的求解方法

1.1动态Lyapunov特征指数

储液罐为非线性定常流固耦合系统，基于势流理论与标准波动方程[11，其动力扰动方程如下：

(1)

(2)

(3)

(4)

式中：Ma=Ms+Mu，Mu=ρfQKf-1QT；ügi(t)为地面加速度，i=1,2,3分别表示x,y,z方向；G为静力作用的结构域节点力列阵；当i= 1、2时，si=Maιsi，ιsi为地震影响向量；耦合作用沿两域接触面法向传递，当i= 3时，仅以-Maιs3üg3(t)引入竖向地震作用并不完整，取si=Maιsi-G/g，本文统一取g= 10 m/s2。

当大位移效应不显著且材料保持弹性时，采用线弹性刚度矩阵Ke替代Ks，并以几何刚度矩阵K0、Kg(t)分别考虑G的预应力效应与时变应力软(刚)化效应：

(5)

(6)

基于Newmark法离散化式(6)，记yk=y(tk)，tk=t0+kτ，t0为初始扰动时刻，τ为时间步长，则有：

(7)

(8)

Le表示微小扰动的平均指数变化率：当Le>0时，系统发生动力失稳；当Le=0时，系统处于临界状态；当Le<0时，系统保持稳定。结构动力失稳的地震动临界峰值加速度acr及初始动力失稳时刻tins可定义为[13]：

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

(9)

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

1.2动力失稳概率

借鉴地震易损性概念，动力失稳概率可定义为[1]：

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

(10)

式中:“INS”表示动力失稳的定量标准；EQ(td,a)表示地震动记录样本，td为持时变量。选取足量地震动记录样本进行多次动力稳定性分析，Pins可表达为[14]：

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

式中：EQ(·)表示满足特定条件的地震动记录样本，num(·)表示样本数量。地震动记录的实际持时长短不一，离散性较大，式(11)是按储液罐动力失稳的充分条件来确定的，即同时满足td≥tins，a≥acr。

地震动的频谱特性无法作为控制变量，Pins可理解为以地震动参数td、a作为变量的二元函数，Pins-td-a曲面可称为动力失稳概率曲面。当td=tins(或a=acr)时，Pins-a(或Pins-td)曲线即为动力失稳概率曲线。

2储液罐的动力失稳概率分析

2.1模型参数及模态分析

本文选用某10×104m3储油罐(见图1)作为分析对象[13]：罐壁共9层，t1～t9为壁厚，钢材弹性模量为206 GPa；储液密度为860 kg/m3，声速取1 400 m/s，液面高度与第9层壁板中部平齐；包边角钢位于罐壁顶部，1 ～ 4号环形抗风圈均由槽钢与钢板焊接而成，依次位于第9、8、7层壁板中部与第5层壁板顶部。

采用有限元软件ANSYS建立压力-位移格式的流固耦合模型，通过科学计算软件MATLAB完成压力变量凝聚及模态分析，而后组集动力扰动方程并编制基于Lyapunov特征指数的动力失稳概率分析模块。罐壁和储液分别采用Shell 181和Fluid 30单元模拟，罐壁网格沿周向均分为120段，每层罐壁沿竖向均分为2段，流体网格与之协调；抗风圈钢板采用Shell 181单元模拟，槽钢与角钢以Beam 188单元模拟，网格与罐壁协调。图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

储液罐的模态分析结果如表1所示，fn为模态频率，m(=Ⅰ,Ⅱ,Ⅲ)按Mni*升序排列。在相同i下，表1所列模态的Mni*所占比重均在95%以上，因此无需取全部模态来确定Kg(t)。由表1可知，在三维地震动作用下，考虑阶次n的重复情况，仅需取7阶不同模态(n= 19, 20, 21, 60, 97, 98, 99)即可确定Kg(t)。

2.2 Lyapunov特征指数方法的验证

由下至上依次递减储罐的抗风圈，考虑水平x向地震动作用，基于Lyapunov特征指数方法(简称LCE方法)对图1所示模型进行动力稳定性分析，并采用基于B-R准则的增量动力分析方法(IDA)进行验证。其中，模态阻尼比取0.02，模态截断数r= 100。Cape Mendocino与Kobe地震波的详细信息见文献[13]，其卓越频率均在2 Hz左右。由图3可知，LCE方法与IDA方法确定的acr值基本一致，前者相对后者的最大误差仅在10%左右，详细数据和论述见文献[13]。

图3　临界地震动峰值加速度αcr的对比Fig.3 Comparison of critical peak ground acceleration αcr

2.3地震波特性分析

根据美国生命线联盟(ALA)的储液罐震害调查报告，结合地震动数据库PEER、COSMOS的记录情况，综合确定表2所示的536组地震动信息，并进行傅里叶分析与反应谱分析(阻尼比取0.02)，地震动卓越频率fp与标准动力放大系数β谱如图4所示，fp表示地震波最大幅值分量所对应的频率，β表示谱加速度与地震动峰值加速度的比值。地震动方位统一采用笛卡尔坐标系标记，x向与东西向(或较小记录角)一致，y向与南北向(或较大记录角)一致，z向对应竖直分量。由图4可知：表2的地震动记录频谱覆盖范围广，反应谱特性丰富；x,y,z向地震动的fp分别集中于1～6 Hz、1～6 Hz、1～10 Hz范围内，最大β值分别为3.20、3.18、3.26；x,y向地震动的fp分布及β均值谱均十分相似。考虑到x,y向地震动的fp相近，z向fp较高，后文以x向的fp值作为地震动记录的主频代表值。

表2　储液罐的地震动数据库

图4　地震动记录的频谱特性Fig.4 Frequency spectrum characteristic of earthquake records

2.4地震动维数及方位的影响

地震动存在很大不确定性，具有多维、频谱复杂等特点，本节按x,y,z,xy,xyz5类方位组合处理表2地震波，然后基于动态Lyapunov特征指数对储液罐进行动力稳定性分析，acr、tins以及动力失稳概率曲面如图5所示。当采用多维地震动时，acr对应x向的地震动峰值加速度，y,z向的a值按原比例调整。

图5　储液罐的动力失稳概率曲面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向地震动下的动力失稳概率曲面非常相似，z向失稳概率曲面最为平缓，xy,xyz向的失稳概率曲面最为陡峭，5类方位组合下的动力失稳概率曲面在原点位置均存在明显底角。

如图6所示，当td=tins(或a=acr)时，图5所示的动力失稳概率曲面退化为曲线形式。由图6(a)可知，当地震峰值加速度a相同时，Pins满足xyz>xy>x=y>z的大小关系，且z向地震动下的Pins-a曲线“坡度”较小，变化趋势明显有异于其余方位；从概率角度而言，5类方位下的Pins-td曲线基本一致，地震动方位及维数对Pins-td曲线的影响不大，相互差异不超过10%。

图6　储液罐的动力失稳概率曲线Fig.6 Dynamic instability probability curves of liquid storage tank

2.5环向抗风圈的影响

环向抗风圈主要用于加强储液罐在风荷载作用下的稳定性，具体设置道数和位置会因周边气候环境而异[15]。由下至上递减储液罐的抗风圈，依次编号为M4、M3、M2、M1、M0，数字代表抗风圈数量。模型M0～M4的动力失稳概率曲线如图7所示。我国战略石油储备基地的最大设防烈度为独山子区的8度(0.20 g)，规范建议此时用于时程分析的罕遇地震加速度峰值为0.40 g，储液罐的动力失稳概率P0.2g、P0.4g分别如表3、4所示，下标表示相应峰值水平。

表3　动力失稳概率P0.2g

表4　动力失稳概率P0.4g

图7　抗风圈对动力失稳概率曲线的影响Fig.7 Effect of wind girders on dynamic instability probability curves

由图7可知：

(1) 从Pins-a曲线来看，Pins随抗风圈道数的增加而显著降低，而且Pins-a曲线在不同地震动方位组合下的收缩程度满足z>x=y>xyz≈xy的递进关系。上述结果表明，从缓解储液罐地震动危险性的相对效果来看，抗风圈在z向地震动下最为有效。

(2) 不同抗风圈设置下的Pins-td曲线非常相近，最大差异不超过20%。按规范建议，时程分析的地震有效持时应满足5～10倍的结构基本周期(约2.5～5 s)，但td= 5 s时的Pins仅在8%左右，故为保证储液罐的动力稳定性，建议持时td取10 s以上。

由表3、4可知：①Pins随着抗风圈的增加而减小，大体符合xyz>xy>x=y>z的大小关系；② 当a= 0.20 g时，储罐M4在x,y向地震动下甚至未发生动力失稳，但其在z,xy,xyz向地震动下仍有可能失稳；③ 当a= 0.40 g时，储罐M4在5类地震动下均有可能发生动力失稳，其在xy,xyz向地震动下的失稳概率可达10%。

本文采用4G内存、2.33 GHz主频的4核计算机，若直接基于压力-位移格式的流固耦合模型进行增量动力分析，完成单次失稳概率分析的时间将数以年计。本文经凝聚压力自由度并参数化结构域抗力，藉由计算动态Lyapunov特征指数极大提高了动力失稳概率的求解效率，完成单次失稳概率分析在8小时左右。

3结论

(1) 本文建立的地震动数据库包含536组3向地震波信息，数量充足，频谱与反应谱特性丰富，能合理反映地震动的不确定性及多维特点。

(2) 经凝聚压力变量并采用时变几何刚度矩阵参数化结构域抗力，Lyapunov特征指数法大幅提升了储液罐动力失稳概率分析的计算效率。相比直接基于流固耦合模型的增量动力分析方法，Lyapunov特征指数法可节约时间成本2 000倍左右。

(3) 竖向地震相对不易引起动力失稳，多维地震下的失稳概率比单维地震提高50%～75%，动力失稳概率与地震动峰值与持时成正相关，地震动方位及维数导致的Pins-td曲线的最大差异不超过10%。

(4) 动力失稳概率随抗风圈的增加而减小，抗风圈对动力失稳的抑制作用在竖向地震下最为显著，Pins-td曲线在不同抗风圈下的最大差别不到20%。

针对储液罐的动力失稳现象，设计分析时应结合规范量化不同抗震设防水准下的动力失稳概率容许值，从而有效评价储液罐选址和设计选型的合理性。

参 考 文 献

[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] 孙建刚, 张荣花, 蒋峰. 储罐地震易损性数值仿真分析[J]. 哈尔滨工业大学学报, 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特征指数及弹性动力失稳[J]. 上海交通大学学报,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. 立式圆筒形钢制焊接油罐设计规范[S]. 北京: 中国计划出版社, 2003.

基金项目：国家自然科学基金资助项目(51178414)

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

通信作者高博青 男，教授，博士生导师，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年生