侯鋼領(lǐng)，張春龍，2，賈曉飛，宋天舒

(1.哈爾濱工程大學(xué) 航天與建筑工程學(xué)院，黑龍江哈爾濱150001;2.中國(guó)核電工程有限公司，北京100000)

LNG儲(chǔ)罐的地震響應(yīng)由于涉及流體與固體的耦合，導(dǎo)致該類結(jié)構(gòu)的動(dòng)水壓力計(jì)算復(fù)雜［1-2］，難以直接應(yīng)用到實(shí)際工程中。國(guó)內(nèi)外廣大學(xué)者［3-4］的研究大多數(shù)是基于有限元法研究動(dòng)水壓力，比如清華大學(xué)的胡盈輝和莊茁［5］提出的附加質(zhì)量法，把動(dòng)水壓力轉(zhuǎn)化成附加質(zhì)量作用到罐壁上，實(shí)現(xiàn)流固耦合解耦;杜顯赫［6］和孫建剛［7］等人通過(guò)液體勢(shì)能函數(shù)的建立，進(jìn)而考慮其流固耦合作用，但其數(shù)值模擬結(jié)果是建立在經(jīng)驗(yàn)判斷上，缺乏理論依據(jù)。本文提出了一種從理論上計(jì)算彈性壁儲(chǔ)罐罐壁Mises應(yīng)力的計(jì)算方法，從而為有限元求解及工程實(shí)際應(yīng)用提供定量評(píng)估。

1 我國(guó)規(guī)范與歐洲規(guī)范計(jì)算動(dòng)水壓力

1.1 我國(guó)規(guī)范的動(dòng)水壓力計(jì)算公式

我國(guó)規(guī)范是基于居榮初和曾心傳等學(xué)者的研究成果［1］，并采用反應(yīng)譜方法。我國(guó)規(guī)范給出的圓柱形容器的動(dòng)水壓力計(jì)算公式為

1.2 基于歐洲規(guī)范動(dòng)水壓力計(jì)算公式

歐洲LNG儲(chǔ)罐的現(xiàn)行規(guī)范《歐洲鋼制常壓立式圓筒形儲(chǔ)罐抗震鑒定標(biāo)準(zhǔn)》(BS EN 1998-4:2006)，該規(guī)范將彈性壁儲(chǔ)罐動(dòng)水壓力分解為:脈沖壓力、對(duì)流壓力和柔性壓力3個(gè)部分，分別為:

儲(chǔ)罐內(nèi)任意一點(diǎn)所受到的脈沖壓力為

儲(chǔ)罐內(nèi)任意一點(diǎn)所受到的對(duì)流壓力為

儲(chǔ)罐內(nèi)任意一點(diǎn)所受到的柔性壓力為

式中:H儲(chǔ)罐內(nèi)液高，R為儲(chǔ)罐半徑，ρ為儲(chǔ)液密度，ρs為儲(chǔ)罐密度，θ為沿環(huán)向任意角度，s(ζ)為罐壁厚度，Ag(t)為地震加速度，Acn(t)為自由振蕩的單自由度加速時(shí)程響應(yīng)，Afn(t)為n階阻尼加速度響應(yīng)，f(ζ)為儲(chǔ)罐裝有液體時(shí)的一階振型函數(shù)(設(shè)f(ζ)=sinπz/2H)，J1為第一類貝塞爾函數(shù)，如圖1所示。

脈沖壓力與對(duì)流壓力合力即是剛性壁動(dòng)水壓力，而以上3個(gè)壓力合力即是彈性壁動(dòng)水壓力。

圖1 儲(chǔ)液罐幾何模型Fig.1 Geometric model of the tank

2 彈性罐壁Mises應(yīng)力理論解

2.1 圓柱殼的基本方程

圓柱殼中曲面的一部分［8］，如圖2所示，其中BC為母線方向，DE為準(zhǔn)線方向，它們是中曲面上的一組曲率線族。取母線作為x線，準(zhǔn)線作為φ線，以任意一曲率作為參考軸，于是中曲面上任一點(diǎn)A的位置可以用(x，φ)，其中x為沿母線方向的線性長(zhǎng)度，而φ為兩條法線之間的夾角。

圖2 圓柱殼示意圖Fig.2 Diagram of cylindrical shell

3個(gè)位移分量u、v、w來(lái)表示一般情況下圓柱殼的基本方程，設(shè)其半徑為a，則R=a=const，并引入無(wú)量綱量ξ=x/a，β=h2/12a2，分別用角標(biāo)1、2來(lái)代替x、φ方向，3方向?yàn)榇怪庇诠薇诘姆较?，則基本方程如下:

與其對(duì)應(yīng)的內(nèi)力表達(dá)式為

式中:E為彈性模量，μ為泊松比，h為罐壁厚度。

如果已知?jiǎng)铀畨毫Φ娜S空間分布，且在滿足給定的邊界條件下，求解方程組(14)～(16)可得到位移分量的解。將位移分量代入到式(17)～(22)，可得到各內(nèi)力素。但由于動(dòng)水壓力公式比較復(fù)雜，如果將動(dòng)水壓力公式直接代入到方程組中進(jìn)行求解，計(jì)算量將非常大，并且難以解出。

在此，研究歐洲規(guī)范中的動(dòng)水壓力規(guī)律，即分析研究在不同情況下，同一水平面和隨高度變化2種情況下動(dòng)水壓力的分布。發(fā)現(xiàn)其規(guī)律為:在同一水平面上，彈性壁儲(chǔ)罐的動(dòng)水壓力可近似為正弦分布;而隨高度變化，最大動(dòng)水壓力的分布也可近似正弦變化。應(yīng)用對(duì)比分析和擬合分析，函數(shù)形式Y(jié)=Psin(Ax+B)sin(πφ/φ0)可以與彈性壁儲(chǔ)罐動(dòng)水壓力值進(jìn)行很好的擬合［9］，并將此代替動(dòng)水壓力計(jì)算公式代入上述方程中進(jìn)行求解。

2.2 靜水壓力作用下圓柱殼有矩理論

圓柱殼在靜水壓力作用下，薄膜內(nèi)力起主導(dǎo)作用，但是圓柱殼底部的彎矩也非常大，不容忽視。現(xiàn)考慮薄膜內(nèi)力與彎曲共同作用情況下圓柱殼內(nèi)力表達(dá)式如下:

式中:H為液面高度，γ為液體重度，μ為泊松比，a為儲(chǔ)罐半徑，h為儲(chǔ)罐壁厚度。

分別將靜水壓力與動(dòng)水壓力作用下罐壁的應(yīng)力進(jìn)行疊加，最后便可求解出彈性壁儲(chǔ)罐在地震作用下罐壁的Mises應(yīng)力。

3 LNG儲(chǔ)罐算例

3.1 LNG儲(chǔ)罐彈性罐壁Mises理論解

將LNG儲(chǔ)罐簡(jiǎn)化計(jì)算取半徑a=40 m，罐高L=30 m，液高H=30 m，彈性罐壁厚h=0.02 m。內(nèi)罐材料為鋼材，密度為 ρs=7.85 × 103kg/m3，彈性模量為E=2.1 × 1011Pa，泊松比 μ =0.3，屈服強(qiáng)度為500 MPa。液體密度 ρ=0.48 × 103kg/m3。

假設(shè)儲(chǔ)罐建造地區(qū)為II類場(chǎng)地，設(shè)計(jì)地震分組為第2組，抗震設(shè)防基本烈度為8度。選用EL Centro波和Manual波2組地震波進(jìn)行時(shí)程分析，峰值加速度為0.7 m/s2;EL波的峰值加速度出現(xiàn)時(shí)間為2.12 s;Manual的峰值加速度出現(xiàn)時(shí)間為 6.26 s，如圖3所示。

圖3 地震波加速度時(shí)程曲線Fig.3 Time history curve of seismic wave acceleration

比較我國(guó)的反應(yīng)譜法與歐洲規(guī)范的動(dòng)水壓力。在水平夾角為0°，在峰值加速度時(shí)刻彈性壁儲(chǔ)罐動(dòng)水壓力沿高度方向的分布，如圖4所示。從圖4可以看出:我國(guó)的規(guī)范方法在計(jì)算儲(chǔ)液上部起主導(dǎo)作用的對(duì)流壓力部分以及儲(chǔ)液下部起主導(dǎo)作用的脈沖壓力部分明顯偏小。在此，采用歐洲規(guī)范為依據(jù)，計(jì)算靜、動(dòng)合壓力作用下彈性罐壁Mises應(yīng)力。

圖4 反應(yīng)譜法與歐洲規(guī)范對(duì)比圖Fig.4 Result comparison between response spectrum method and European standard

根據(jù)式(4)～(13)，并應(yīng)用MATLAB軟件實(shí)現(xiàn)了彈性壁儲(chǔ)罐分別在EI波和Manual波峰值加速度作用時(shí)的動(dòng)水壓力空間分布，見(jiàn)圖5和圖6所示。圖6為靜水和動(dòng)水壓力合力的三維空間圖。應(yīng)用擬合函數(shù)形式Y(jié)=Psin(Ax+B)sin(πφ/φ0)和最小二乘法，得到P=18 586.7，A=0.038，B=7.39，φ0=π 。

圖5 彈性壁儲(chǔ)罐動(dòng)水壓力值Fig.5 Hydrodynamic pressure on elastic wall

圖6 彈性壁儲(chǔ)罐靜動(dòng)合壓力值Fig.6 Hydrodynamic and hydrostatic pressure on elastic wall

圖7為進(jìn)行曲線擬合精度的對(duì)比。可以看出，擬合曲線除在最低點(diǎn)的誤差為4%以外，其余點(diǎn)的誤差均小于3%，最大為2.89%，平均誤差為1%，說(shuō)明曲線擬合已經(jīng)達(dá)到精度要求。

圖7 計(jì)算數(shù)值與曲線擬合對(duì)比Fig.7 Contrast between numerical result and curve fitting

根據(jù)方程組(13)～(16)，可知P1=0，P2=0，P3=Psin(Ax+B)sin(πφ/φ0)，設(shè) 3 個(gè)位移分量如下:

它們滿足所有的位移方程，同時(shí)也滿足了四邊簡(jiǎn)支的邊界條件，其中U、V、W是需要計(jì)算的未知參數(shù)，再將P1、P2、P3和 3個(gè)位移分量代入式(13)～(16)，化簡(jiǎn)后得方程組:

求解上述方程組，得到位移分量U=0.000 6，V=-0.008 3，W=0.014 5。將上面 3 個(gè)位移分量代入到式(17)～(22)，可得到簡(jiǎn)后得彈性壁儲(chǔ)罐在動(dòng)水壓力作用下的內(nèi)力公式:

取單位長(zhǎng)度，b=1，由上面動(dòng)水壓力作用下的彈性罐壁內(nèi)力公式得到

式中:下標(biāo)wd表示動(dòng)水壓力作用下的彎矩，md為動(dòng)水壓力作用下的膜應(yīng)力，d代表動(dòng)水作用。

同理，由式(23)～(26)靜水壓力作用下的彈性罐壁內(nèi)力公式得到

式中:下標(biāo)wj表示動(dòng)水壓力作用下的彎矩，mj為動(dòng)水壓力作用下的膜應(yīng)力，j代表靜水作用。

則由上面公式得到彈性罐壁在地震峰值加速度作用下的主應(yīng)力分別為

最后得到彈性罐壁在地震峰值加速度作用下的Mises應(yīng)力為

由于自重影響相對(duì)于靜水壓力和動(dòng)水壓力的影響很小，自重影響可以忽略不計(jì)。上述彈性罐壁的Mises應(yīng)力計(jì)算未考慮彈性罐壁自重的影響。

3.2 耦合歐拉-朗格朗日(CEL)算法

ABAQUS中的CEL算法是完全的流固耦合算法，分析類型為動(dòng)態(tài)顯示，具有較高的計(jì)算精度，并且已經(jīng)成為檢驗(yàn)理論簡(jiǎn)化解的重要工具。流體構(gòu)型通過(guò)計(jì)算流體在歐拉單元中所占的體積分?jǐn)?shù)確定。流體材料可以在歐拉單元之間流動(dòng)，并與固體單元相互作用。計(jì)算模型時(shí)需要一個(gè)足夠大的歐拉網(wǎng)格，能夠?qū)⒐腆w單元可能移動(dòng)到的所有位置全部包裹在內(nèi)，由于CEL計(jì)算量較大，所以模型采用對(duì)稱建模方式如圖8所示，取一半模型進(jìn)行分析。罐壁單元尺寸為0.8 m，歐拉體單元尺寸為 0.8 m，如圖 9 所示。

圖8 LNG儲(chǔ)罐整體模型Fig.8 Model of LNG tank

圖9 LNG儲(chǔ)罐網(wǎng)格劃分Fig.9 Meshing model of LNG tank

4 計(jì)算結(jié)果及分析

針對(duì)上節(jié)算例，分別計(jì)算在EI波和Manual波峰值加速度作用下，彈性罐壁在0°和180°位置罐壁的Mises應(yīng)力隨高度的變化，如圖10～13所示。

因?yàn)楣薇贛ises應(yīng)力起主導(dǎo)作用的位置主要在中下部［10］，所以本文提取了罐壁0～15 m的Mises應(yīng)力隨高度變化的數(shù)值，從圖10和圖12可以看出，罐壁Mises應(yīng)力隨高度的變化為底部大致呈正弦變化趨勢(shì)，再往上則逐漸成線性變化趨勢(shì)，罐壁底部的正弦變化趨勢(shì)主要是由于靜水壓力作用時(shí)罐壁彎曲應(yīng)力引起的，而再往上的線性變化趨勢(shì)則是由罐壁的膜應(yīng)力所引起的。

從圖中還可以看出，對(duì)罐壁Mises應(yīng)力起主導(dǎo)作用的是靜水壓力，而非動(dòng)水壓力;彈性儲(chǔ)罐在EI波作用時(shí)罐壁最大Mises應(yīng)力在底部，理論計(jì)算數(shù)值為 461.14 MPa，仿真數(shù)值為 483.5 MPa，最大 Mises應(yīng)力出現(xiàn)時(shí)間為3.5 s滯后于峰值加速度1.38 s。Manual波作用時(shí)，罐壁底部最大Mises應(yīng)力理論計(jì)算數(shù)值為 461.14 MPa，仿真數(shù)值為 483.175 MPa，最大Mises應(yīng)力出現(xiàn)時(shí)間為4.7 s，提前于峰值加速度1.56 s。

在Manual波作用下，彈性儲(chǔ)罐罐壁最大Mises應(yīng)力出現(xiàn)時(shí)間提前于峰值加速度時(shí)間其原因如下。分別對(duì)EI波和Manual波加速度時(shí)程曲線進(jìn)行傅里葉變換，可以得到加速度頻譜曲線如圖14～15。

圖10 EI波對(duì)比0°Fig.10 Comparison results under EI wave at 0°

圖 11 EI波結(jié)果 t=3.5 sFig.11 Computing results under EI wave on 3.5 s

圖12 Manual波對(duì)比180°Fig.12 Comparison results under Manual wave at 180°

圖 13 Manual結(jié)果 t=4.7 sFig.13 Computing results under Manual wave on 4.7 s

圖14 EI波頻譜曲線Fig.14 EI spectrum

圖15 Manual波頻譜曲線Fig.15 Manual spectrum

通過(guò)對(duì)頻譜特性［11］分析可知，EI波前20 s的卓越頻率為 1.785 7 Hz，前 3.5 s 的卓越頻率為 1.785 7 Hz。而 Manual波前20 s的卓越頻率為3.083 Hz，前4.7 s的卓越頻率為 3.756 Hz。

根據(jù)文獻(xiàn)［12-14］，儲(chǔ)罐滿液時(shí)的自振頻率為

式中:fT為空罐子時(shí)的自振頻率;R為儲(chǔ)罐半徑;ρ為儲(chǔ)液密度;h為罐壁厚;ρs為儲(chǔ)罐密度;β值查表，本文中 β 取 0.23。

5 結(jié)論

本文通過(guò)對(duì)LNG罐壁Mises應(yīng)力理論計(jì)算與有限元仿真結(jié)果相對(duì)比，得出結(jié)論:

1)基于我國(guó)反應(yīng)譜法計(jì)算出的彈性壁儲(chǔ)罐的動(dòng)水壓力在儲(chǔ)液上部的對(duì)流壓力值和儲(chǔ)液下部的脈沖壓力值都偏小。

2)通過(guò)彈性罐壁Mises應(yīng)力隨高度變化的趨勢(shì)圖中可以看出，對(duì)罐壁Mises應(yīng)力起到主導(dǎo)作用的是靜水壓力對(duì)罐壁的作用，而非動(dòng)水壓力。

3)彈性罐壁Mises應(yīng)力隨高度變化的趨勢(shì)為底部大致呈正弦趨勢(shì)，此主要是由于靜水壓力作用部分對(duì)罐壁的彎曲應(yīng)力引起的。往上則逐漸呈線性變化趨勢(shì)，此主要是由于罐壁的膜應(yīng)力所引起的。所以彈性罐壁Mises應(yīng)力在底部起主導(dǎo)作用的是罐壁的彎曲應(yīng)力，而再往上則罐壁的膜應(yīng)力逐漸起到主導(dǎo)作用;

4)當(dāng)?shù)卣鸺铀俣鹊淖吭筋l率與儲(chǔ)液罐的自振頻率比較接近時(shí)，導(dǎo)致儲(chǔ)液罐出現(xiàn)共振或亞共振狀態(tài)時(shí)，會(huì)造成罐壁最大的Mises應(yīng)力提前于地震峰值加速度出現(xiàn)。

［1］居榮初，曾心傳.彈性結(jié)構(gòu)與液體的耦聯(lián)振動(dòng)理論［M］.北京:地震出版社，1983:40-48.

［2］SH/T 3026-2005鋼制常壓立式圓筒形儲(chǔ)罐抗震鑒定標(biāo)準(zhǔn)［S］.上海:中國(guó)石化出版社，2006.

［3］MOSLEMI M，KIANOUSH M R.Parameter study on dynamic behavior of cylindrical ground-supported tanks［J］.Engineering Structures，2012，42:214-230.

［4］MALEKI A，ZIYAEIFAR M.Damping enhancement of seismic isolated cylindrical liquid storage tanks using baffles［J］.Engineering Structures，2007，29(12):3227-3240.

［5］胡盈輝，莊茁，由小川.大型儲(chǔ)液罐在地震作用下的附加質(zhì)量法研究［J］.壓力容器，2009，26(8):1-6.HU Yinghui，ZHUANG Zhuo，YOU Xiaochuan.Added mass approach to a large-scale liquid-storage tank under seismic excitations［J］.Pressure Vessel Technology，2009，26(8):1-6.

［6］杜顯赫，沈新普，劉應(yīng)華.預(yù)應(yīng)力LNG儲(chǔ)罐在地震作用下的流固耦合數(shù)值模擬［C］//力學(xué)與工程應(yīng)用.鄭州，中國(guó)，2012:178-182.DU Xianhe，SHEN Xinpu，LIU Yinghua.Numerical simulation of fluid and structure coupling of prestressed LNG storage tank under seismic action［C］//Mechanics and Engineering Application.Zhengzhou，China，2012:178-182.

［7］孫建剛，崔利富，張營(yíng).全容式LNG儲(chǔ)罐地震響應(yīng)數(shù)值模擬研究［C］//中國(guó)國(guó)際建筑科技大會(huì)論文集.北京，中國(guó)，2010:225-230.SUN Jian’gang，CUI Lifu，ZHANG Ying.Numerical simulation of seismic response of the full capacity LNG storage tank［C］//Chinese International Building Science and Technology Conference.Beijing，China，2010:225-230.

［8］黃克智.板殼理論［M］.北京:清華大學(xué)出版社，1987:137-146.

［9］史潔玉，孔玲軍.MATLAB R2012a超級(jí)學(xué)習(xí)手冊(cè)［M］.北京:人民郵電出版社，2013:72-85.

［10］大崎順彥.地震動(dòng)的譜分析入門［M］.田琪，譯.第2版.北京:地震出版社，2008:49-58.

［11］翁智遠(yuǎn).圓柱薄殼容器的振動(dòng)與屈曲［M］.上海:上?？茖W(xué)技術(shù)出版社，2007:3-16.

［12］翁智遠(yuǎn)，楊建中.圓柱形彈性貯液容器基頻近似計(jì)算法［J］.力學(xué)季刊，2000，21(3):354-356.WENG Zhiyuan，YANG Jianzhong.An approximate calculational method for fundamental frequency of elastic cylindrical liquid storage tank［J］.Chinese Quarterly of Mechanics，2000，21(3):354-356.

［13］FISCHER F D，RAMMERSTORFER F G，SCHARF K.Earthquake resistant design of anchored and unanchored liquid storage tanks under three-dimensional earthquake excitation［M］//SCHU?LLER G I.Structural Dynamics.Berlin:Springer Verlag，1991:317-371.

［14］FISCHER F D，RAMMERSTORFER F G，SCHARF K，et al.Collapse of earthquake excited tanks［J］.Mechanica，1988，25:129-143.