呂 遠(yuǎn)， 孫建剛， 孫宗光， 崔利富， 王 振

(1.大連海事大學(xué) 道路與橋梁工程研究所，遼寧 大連 116026；2.大連民族大學(xué) 土木工程學(xué)院，遼寧 大連 116650)

球形儲(chǔ)罐是石油化工領(lǐng)域非常重要的存儲(chǔ)設(shè)備，其存儲(chǔ)介質(zhì)往往為生產(chǎn)過(guò)程中易燃易爆的原料、過(guò)程料和成品。此類結(jié)構(gòu)一旦遭遇強(qiáng)震，可能發(fā)生儲(chǔ)罐容器墜落倒塌引起爆炸，進(jìn)而引發(fā)火災(zāi)等。因此研究球形儲(chǔ)罐在地震作用下的動(dòng)力響應(yīng)及抗震、減震措施是很有必要的?！稑?gòu)筑物抗震設(shè)計(jì)規(guī)范》[1]對(duì)球形儲(chǔ)罐的抗震設(shè)計(jì)做出了相應(yīng)規(guī)定，其忽略儲(chǔ)液晃動(dòng)效應(yīng)，采用等效附加質(zhì)量法將球罐簡(jiǎn)化為一個(gè)單質(zhì)點(diǎn)體系進(jìn)行抗震設(shè)計(jì)，但在地震動(dòng)作用下球罐的振動(dòng)為罐內(nèi)儲(chǔ)液與罐壁及支承結(jié)構(gòu)的液固藕聯(lián)振動(dòng)，其受力過(guò)程較為復(fù)雜。事實(shí)上，有關(guān)儲(chǔ)罐類的液固耦合振動(dòng)國(guó)內(nèi)外學(xué)者已經(jīng)進(jìn)行了大量研究，但主要針對(duì)立式圓筒形儲(chǔ)罐，而對(duì)舉架式球形儲(chǔ)罐的研究相對(duì)較少。20世紀(jì)，學(xué)者們對(duì)球形儲(chǔ)罐的晃動(dòng)問(wèn)題進(jìn)行了理論和實(shí)驗(yàn)研究[1-7]。Ramaneyulu等[8]采用有限元數(shù)值仿真手段，對(duì)LPG球形儲(chǔ)罐進(jìn)行了地震荷載作用下的可靠性評(píng)估。Papaspyrou等[9]采用速度勢(shì)理論針對(duì)50%儲(chǔ)液球形儲(chǔ)罐，推導(dǎo)出了一種用于計(jì)算儲(chǔ)液動(dòng)液壓力等儲(chǔ)液晃動(dòng)效應(yīng)的數(shù)學(xué)解析模型，并研究了不同外部激勵(lì)條件時(shí)的儲(chǔ)液晃動(dòng)效應(yīng)。Lazaros等[10-11]采用速度勢(shì)理論推導(dǎo)出了球形儲(chǔ)罐和水平圓柱形儲(chǔ)罐線性晃動(dòng)效應(yīng)的數(shù)學(xué)模型，并以50%儲(chǔ)液時(shí)為例進(jìn)行了地震動(dòng)響應(yīng)研究。肖志剛[12]利用有限元軟件對(duì)球形儲(chǔ)罐進(jìn)行了考慮儲(chǔ)液晃動(dòng)的地震動(dòng)響應(yīng)研究，認(rèn)為罐內(nèi)儲(chǔ)液的晃動(dòng)吸收了部分地震能量，有利于球罐結(jié)構(gòu)安全。戴鴻哲等[13]采用有限元數(shù)值模擬手段，對(duì)液化石油球罐進(jìn)行了地震響應(yīng)研究，并分析了儲(chǔ)液晃動(dòng)對(duì)球罐受力的影響，認(rèn)為球罐罐體在地震作用下是偏于安全的，其支承體系才是薄弱點(diǎn)。Seyyed等[14]采用線性勢(shì)流體理論及Mehler-Fock變換求解球罐儲(chǔ)液晃動(dòng)問(wèn)題，并研究了不同儲(chǔ)液高度對(duì)動(dòng)液壓力的影響。

目前球形儲(chǔ)罐考慮儲(chǔ)液晃動(dòng)的地震動(dòng)響應(yīng)研究已經(jīng)取得了一些成果，但仍處于理論與實(shí)驗(yàn)的研究階段，且理論求解過(guò)程十分復(fù)雜，不便于工程應(yīng)用。因此為了更加真實(shí)地反應(yīng)球形儲(chǔ)罐在地震動(dòng)作用時(shí)的液固藕聯(lián)振動(dòng)動(dòng)力響應(yīng)，本文采用速度勢(shì)剛性理論，根據(jù)邊界條件推導(dǎo)出合理的勢(shì)函數(shù)，并進(jìn)一步推導(dǎo)出球形儲(chǔ)罐在地震動(dòng)作用時(shí)的動(dòng)液壓力、儲(chǔ)液晃動(dòng)波高、支承底部剪力及傾覆彎矩表達(dá)式，構(gòu)建了便于工程應(yīng)用的球形儲(chǔ)罐考慮儲(chǔ)液晃動(dòng)簡(jiǎn)化動(dòng)力學(xué)模型，進(jìn)行了地震動(dòng)響應(yīng)研究并與有限元數(shù)值仿真模擬計(jì)算結(jié)果進(jìn)行對(duì)比，輔證了理論推導(dǎo)的準(zhǔn)確性。

1 基本理論

1.1 基本假定

由于球形儲(chǔ)罐相對(duì)于大型立式儲(chǔ)罐來(lái)說(shuō)體積較小，但壁厚較厚，且由于其球形形狀及主要存儲(chǔ)輕質(zhì)油，所以在液體的晃動(dòng)過(guò)程中罐壁發(fā)生的變形極其小，對(duì)工程應(yīng)用來(lái)說(shuō)達(dá)到了可忽略不計(jì)的程度。故將其假設(shè)為剛性罐壁，設(shè)x0是球罐罐體在地震動(dòng)作用時(shí)的相對(duì)位移，其是時(shí)間相關(guān)變量，同時(shí)我們也知道在由于球罐支承的存在，罐體各部位的位移未必相同，因此x0也是空間相關(guān)變量，記作x0(R,θ,β,t)。球形儲(chǔ)罐在地震動(dòng)作用中主要激發(fā)第一階振型，故主要以第一階振型作為研究對(duì)象，現(xiàn)利用有限元軟件，對(duì)球形儲(chǔ)罐進(jìn)行模態(tài)分析，研究其第一階振型的振動(dòng)形態(tài)。取某一工程實(shí)例，建立有限元數(shù)值仿真模型，并進(jìn)行模態(tài)分析。

由圖1可知，對(duì)球罐罐體來(lái)說(shuō)，第一階振型可近似分解為球罐罐體水平平動(dòng)及罐體繞赤道面的旋轉(zhuǎn)，由于球罐支承系統(tǒng)的抗彎剛度比抗剪剛度大得多，屬于剪切型結(jié)構(gòu)，故第一階振型以水平平動(dòng)為主，轉(zhuǎn)動(dòng)十分微弱，在工程計(jì)算時(shí)為簡(jiǎn)化計(jì)算可忽略轉(zhuǎn)動(dòng)影響。Lazaros等在進(jìn)行儲(chǔ)液晃動(dòng)效應(yīng)對(duì)工業(yè)球罐抗震設(shè)計(jì)影響的研究時(shí)便做過(guò)類似假定，其忽略了罐體轉(zhuǎn)動(dòng)影響，假設(shè)球罐罐體整體平動(dòng)，在地震動(dòng)作用時(shí)罐體各部位位移僅為時(shí)間相關(guān)變量，記作x0(t)。

(a) 第一階振型

(b) 罐體水平平動(dòng)

(c) 罐體轉(zhuǎn)動(dòng)圖1 球罐第一階振型Fig.1 Spherical tank first order mode

同時(shí)假定儲(chǔ)液為無(wú)旋、無(wú)粘、不可壓縮的理想流體。在地震作用下，根據(jù)速度勢(shì)剛性理論[15-16]，儲(chǔ)液地震動(dòng)響應(yīng)可分為對(duì)流運(yùn)動(dòng)和剛性運(yùn)動(dòng)，儲(chǔ)液總速度勢(shì)可分為剛性速度勢(shì)和對(duì)流晃動(dòng)速度勢(shì)，記為Φ(r,θ,β,t)=φr(r,θ,β,t)+φs(r,θ,β,t)，滿足球坐標(biāo)系下的Laplace方程：

(1)

球罐計(jì)算簡(jiǎn)圖如圖2所示，地面運(yùn)動(dòng)xg(t)作用下，罐體會(huì)產(chǎn)生相對(duì)位移為x0(t)的振動(dòng)。因?yàn)棣誶、φs均滿足方程(1)，故分別討論剛性速度勢(shì)和晃動(dòng)速度勢(shì)。

圖2 球罐簡(jiǎn)圖Fig.2 Spherical tank diagram

1.2 剛性速度勢(shì)

根據(jù)上述假定，剛性速度勢(shì)φr(r,θ,β,t)應(yīng)滿足Laplace方程(1)和邊界條件：

φr(β+2π)=φr(β)

(2)

(3)

(4)

(5)

1.3 晃動(dòng)速度勢(shì)

同剛性速度勢(shì)，晃動(dòng)速度勢(shì)φs(r,θ,β,t)也滿足Laplace方程(1)和邊界條件：

φs(β+2π)=φs(β)

(6)

(7)

(8)

m=0,1,2,3…

(9)

其中hv為表面波動(dòng)方程，因此可得：

(10)

而又因?yàn)閤=rsinθcosβ，所以存在

式(10)可以寫為：

(11)

將速度勢(shì)方程式(5)、(9)代入式(11)可得：

(12)

由于式(12)在自由液面S2上均滿足，根據(jù)疊加原理，可將兩邊同時(shí)乘以調(diào)和函數(shù)φ*(r,θ)，并在圓域S2內(nèi)積分，所以可以得到如下方程：

sinθφ*(r,θ)rdrdβ=

(13)

(14)

同時(shí)由邊界條件(8)可得：

(15)

式(15)在儲(chǔ)液與罐壁耦合面S1(r=R)上均滿足，同樣可根據(jù)疊加原理，將兩邊同時(shí)乘以調(diào)和函數(shù)φ*(r,θ)，并在圓域S1內(nèi)積分，所以可以得到如下方程：

cosθR2sinθdθdβ=0

(16)

將式(14)、(16)相加可得：

(17)

式(17)可以寫成矩陣的形式：

(18)

其中：

(19)

(20)

(21)

(22)

[f]=[fn(t)]m×1

(23)

(24)

其中

[M]m×m=[P][N1]T

(25)

[K]m×m=[P]([N2]T-[N3]T+[N4]T)

(26)

(27)

(28)

(29)

(30)

參照結(jié)構(gòu)動(dòng)力學(xué)振型疊加原理[18]，可以將式(24)看作是耦合在一起的m階線性無(wú)阻尼運(yùn)動(dòng)控制方程，現(xiàn)在需要對(duì)其進(jìn)行解耦，首先求出各階振型及對(duì)應(yīng)晃動(dòng)頻率：

(31)

式中：ωn為n階晃動(dòng)頻率，{ψn}為對(duì)應(yīng)的n階晃動(dòng)振型；根據(jù)上述內(nèi)容可知，式(24)為m階矩陣方程，其中m=1,2,3,4，為截?cái)鄶?shù)，其取值不同時(shí)結(jié)果精度會(huì)產(chǎn)生差異；

因此式(24)可以轉(zhuǎn)化為m個(gè)非耦合的方程：

(32)

其中

(33)

可將式(32)寫為如下形式：

(34)

其中

(35)

式(34)為晃動(dòng)分量無(wú)阻尼運(yùn)動(dòng)方程，有阻尼運(yùn)動(dòng)方程如下所示：

(36)

由于儲(chǔ)液晃動(dòng)主要以第一階振型為主，所以只以i=1時(shí)為主要研究對(duì)象，記xc1(t)=xc(t)。

根據(jù)式(32)～(35)可得儲(chǔ)液晃動(dòng)速度勢(shì)為：

(37)

所以總速度勢(shì)為：

(38)

有了總速度勢(shì)，可得得出重力場(chǎng)作用下自由液面的振動(dòng)形態(tài)，及儲(chǔ)液作用在罐壁上的動(dòng)液壓力：

(39)

(40)

將式(38)分別代入式(39)、(40)中可得：

(41)

(42)

則由儲(chǔ)液運(yùn)動(dòng)產(chǎn)生的水平方向基底剪力為：

(43)

可以寫為：

(44)

由儲(chǔ)液慣性力而產(chǎn)生的對(duì)支柱底部的傾覆彎矩也是球罐地震動(dòng)響應(yīng)研究時(shí)的主要研究對(duì)象：

M1(t)=

(45)

其中：

(46)

(47)

(a) e值

(b) 晃動(dòng)分量系數(shù)

(c) 剛性分量等效高度

(d) 晃動(dòng)分量等效高度

(e) 晃動(dòng)頻率圖3 不同截?cái)鄶?shù)m時(shí)各參數(shù)曲線Fig.3 The parameters curve of different truncated

由圖3可知，隨著截?cái)鄶?shù)m的增大，各參數(shù)曲線趨近于某一值。當(dāng)m≥8時(shí)各參數(shù)數(shù)值變化已十分微弱，基本滿足計(jì)算精度要求，因此取截?cái)鄶?shù)m=9。由于上述各參數(shù)也受儲(chǔ)罐半徑R的影響，因此研究不同儲(chǔ)罐半徑時(shí)，各參數(shù)的變化情況，計(jì)算結(jié)果如圖4所示。

(a) e值

(b) 晃動(dòng)分量系數(shù)

(c) 剛性分量等效高度

(d) 晃動(dòng)分量等效高度

(e) 晃動(dòng)頻率圖4 不同儲(chǔ)液高度時(shí)各參數(shù)曲線Fig.4 The parameters curve of different height of the fluid

(48)

s=0.517 9+0.538 1cos(1.377x)-0.072 68sin(1.377x)-0.059 57cos(2.754x)-0.064 12sin(2.754x)

(49)

(50)

hc=(0.035 01x5-0.149 3x4+0.232x3-0.147x2+0.372x-0.002 425)2R+h

(51)

h0=(-0.013 76x5+0.050 26x4-0.066 86x3+0.030 77x2+0.274 8x-0.000 725 7)2R+h

(52)

表1 曲線擬合優(yōu)度R2Tab.1 Goodness of fit of curves R2

將本文推導(dǎo)的晃動(dòng)分量系數(shù)與晃動(dòng)頻率近似解析解與Lazaros[11]值作對(duì)比，如圖5所示，發(fā)現(xiàn)兩者十分接近，由此也輔證了本文推導(dǎo)過(guò)程的正確性。

(a) 晃動(dòng)分量系數(shù)對(duì)比

(b) 晃動(dòng)頻率對(duì)比圖5 數(shù)值對(duì)比Fig.5 Numerical comparison

2 簡(jiǎn)化力學(xué)模型

上述推導(dǎo)中所求出的Q1,M1分別為儲(chǔ)液在地震動(dòng)作用下而產(chǎn)生的作用于支承底部的剪力和彎矩。地震作用時(shí)，支柱底部總剪力和總彎矩還應(yīng)包含由球罐罐體、支承質(zhì)量、球罐配件等產(chǎn)生的剪力和彎矩?！稑?gòu)筑物抗震設(shè)計(jì)規(guī)范》(GB 50191—2012)[1]中將球罐簡(jiǎn)化為一個(gè)單質(zhì)量體系:

meq=m1+m2+0.5m3+m4+m5

(53)

式中:m1為球殼質(zhì)量；m2為液體有效質(zhì)量；m3為支柱和拉桿質(zhì)量；m4為保溫層質(zhì)量；m5為球罐其他附件質(zhì)量，包括各開(kāi)口，噴淋裝置，梯子和平臺(tái)等。

由于m4、m5一般附加于罐體，1/2支柱和拉桿質(zhì)量集中于支承頂部附加在罐體上，為便于工程設(shè)計(jì)，此時(shí)也可將m1、0.5m3、m4、m5整合為一個(gè)集中質(zhì)量，記作

ms=m1+0.5m3+m4+m5

(54)

可設(shè)地震動(dòng)作用時(shí)集中質(zhì)量ms對(duì)應(yīng)的相對(duì)位移為罐體相對(duì)位移x0(t)。

進(jìn)而可得由于其慣性作用而產(chǎn)生的基底剪力為：

(55)

總基底剪力為：

(56)

由罐體，支承體系及附件等產(chǎn)生的基底彎矩表達(dá)式為：

(57)

式中：h1為球殼集中質(zhì)量等效高度；h3為支承質(zhì)量等效高度；h4為保溫層質(zhì)量等效高度；h5…h(huán)n為各附件等效高度，根據(jù)不同球罐結(jié)構(gòu)形式，表示多種附件對(duì)應(yīng)的等效高度。

對(duì)球殼來(lái)說(shuō)其對(duì)支柱底部產(chǎn)生的彎矩可表達(dá)為

m1h1=

(58)

求解積分，可得：

(59)

式中：R1為球罐內(nèi)壁半徑；R2為球罐外壁半徑；ρ1為罐體密度；R1-R2為球罐罐壁厚度，相對(duì)于R1、R2來(lái)說(shuō)十分小，因此式(59)可以寫為：

(60)

(61)

總基底彎矩為：

(62)

則根據(jù)式(56)、(62)可以構(gòu)建球形儲(chǔ)罐考慮儲(chǔ)液晃動(dòng)時(shí)的簡(jiǎn)化動(dòng)力學(xué)模型

圖6中kc,cc,k0,c0,ms分別為：

kc=mcω2,cc=2ξωmc

(63)

式中:k0,c0分別為支承結(jié)構(gòu)的剛度和阻尼，可根據(jù)規(guī)范[1]相關(guān)公式可得出；ms為儲(chǔ)罐和支承結(jié)構(gòu)質(zhì)量，其等效高度為h+R。由Hamilton原理，可推導(dǎo)相應(yīng)的運(yùn)動(dòng)方程：

圖6 球罐簡(jiǎn)化力學(xué)模型Fig.6 Spherical tank simplified mechanical model δ(T-V)dt+δWdt=0

(64)

式中：T,V分別為系統(tǒng)的動(dòng)能和勢(shì)能；W為非保守力做的功。

由式(64)可得運(yùn)動(dòng)方程為：

(65)

3 算例分析

選取某一1 000 m3液化石油氣罐作為算例，儲(chǔ)液高度為H=1.5R，忽略其內(nèi)壓影響，儲(chǔ)液密度為480 kg/m3，球罐直徑為12.3 m，球心距地面8 m，拉桿上部連接處距地面6 m，球殼與支承體系均采用雙線性強(qiáng)化模型，算例參數(shù)如表2所示。分別選擇四類場(chǎng)地中滿足規(guī)范[1]的五條天然波和兩條人工波對(duì)球形儲(chǔ)罐簡(jiǎn)化力學(xué)模型和有限元數(shù)值仿真模型進(jìn)行地震動(dòng)響應(yīng)對(duì)比研究，四種場(chǎng)地波加速度反應(yīng)譜如圖7所示，調(diào)整加速度時(shí)程曲線峰值為0.2 g，計(jì)算結(jié)果,如表3～6所示。

表2 模型參數(shù)Tab.2 Model parameters

從表3～表6中數(shù)據(jù)可以看出，對(duì)不同場(chǎng)地地震動(dòng)輸入，球罐地震動(dòng)響應(yīng)差異較大，從(類場(chǎng)地到Ⅳ類場(chǎng)地，晃動(dòng)波高、基底剪力及傾覆彎矩均值逐漸增大。變異系數(shù)反應(yīng)了各工況計(jì)算數(shù)值在同場(chǎng)地不同地震波條件下的離散程度，從表中數(shù)據(jù)可知，變異系數(shù)均較小，說(shuō)明離散程度較小，選波較為合理。而均值差異率則反應(yīng)了本文推導(dǎo)的理論解與有限元解的差異大小，從表中可以看出，均值差異率均在10%以下，有限元解與理論解十分接近，且理論解均值均大于有限元解均值，計(jì)算結(jié)果偏于安全，由此則進(jìn)一步佐證了本文所構(gòu)建的簡(jiǎn)化力學(xué)模型的準(zhǔn)確性及可靠性。

(a) Ⅰ類場(chǎng)地

(b) Ⅱ類場(chǎng)地

(c) Ⅲ類場(chǎng)地

(d) Ⅳ類場(chǎng)地圖7 四類場(chǎng)地加速度反應(yīng)譜Fig.7 Four types of site acceleration response spectrum

表3 Ⅰ類場(chǎng)地地震動(dòng)響應(yīng)對(duì)比Tab.3 ClassⅠsite ground motion response comparison

注：差異率=(理論解-有限元解)/理論解

表4 Ⅱ類場(chǎng)地地震動(dòng)響應(yīng)對(duì)比Tab.4 ClassⅡsite ground motion response comparison

表5 Ⅲ類場(chǎng)地地震動(dòng)響應(yīng)對(duì)比Tab.5 Class Ⅲ site ground motion response comparison

表6 Ⅳ類場(chǎng)地地震動(dòng)響應(yīng)對(duì)比Tab.6 Class Ⅳ site ground motion response comparison

4 結(jié) 論

(1) 采用速度勢(shì)剛性理論，根據(jù)邊界條件推導(dǎo)出合理的勢(shì)函數(shù)，并進(jìn)一步推導(dǎo)出球形儲(chǔ)罐在地震動(dòng)作用時(shí)的動(dòng)液壓力、儲(chǔ)液晃動(dòng)波高、支承底部剪力及傾覆彎矩表達(dá)式，并分析了不同截?cái)鄶?shù)m及求球罐半徑R對(duì)各參量的影響，將本文推導(dǎo)的晃動(dòng)分量系數(shù)和晃動(dòng)頻率近似解析解與與Lazaros值作對(duì)比，兩者十分接近，由此也輔證了本文推導(dǎo)過(guò)程的正確性。

(2) 構(gòu)建了便于工程應(yīng)用的球形儲(chǔ)罐考慮儲(chǔ)液晃動(dòng)簡(jiǎn)化動(dòng)力學(xué)模型，進(jìn)行了地震動(dòng)響應(yīng)研究并與有限元數(shù)值仿真模擬計(jì)算結(jié)果進(jìn)行對(duì)比，有限元解與理論解十分接近，且理論解均值均大于有限元解均值，計(jì)算結(jié)果偏于安全，由此則進(jìn)一步佐證了本文所構(gòu)建的簡(jiǎn)化力學(xué)模型的準(zhǔn)確性及可靠性。