LNG儲罐在地震作用下彈性罐壁Mises應力理論解

侯鋼領，張春龍，2，賈曉飛，宋天舒

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

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

1 我國規范與歐洲規范計算動水壓力

1.1 我國規范的動水壓力計算公式

我國規范是基于居榮初和曾心傳等學者的研究成果［1］，并采用反應譜方法。我國規范給出的圓柱形容器的動水壓力計算公式為

1.2 基于歐洲規范動水壓力計算公式

歐洲LNG儲罐的現行規范《歐洲鋼制常壓立式圓筒形儲罐抗震鑒定標準》(BS EN 1998-4:2006)，該規范將彈性壁儲罐動水壓力分解為:脈沖壓力、對流壓力和柔性壓力3個部分，分別為:

儲罐內任意一點所受到的脈沖壓力為

儲罐內任意一點所受到的對流壓力為

儲罐內任意一點所受到的柔性壓力為

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

脈沖壓力與對流壓力合力即是剛性壁動水壓力，而以上3個壓力合力即是彈性壁動水壓力。

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

2 彈性罐壁Mises應力理論解

2.1 圓柱殼的基本方程

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

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

3個位移分量u、v、w來表示一般情況下圓柱殼的基本方程，設其半徑為a，則R=a=const，并引入無量綱量ξ=x/a，β=h2/12a2，分別用角標1、2來代替x、φ方向，3方向為垂直于罐壁的方向，則基本方程如下:

與其對應的內力表達式為

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

如果已知動水壓力的三維空間分布，且在滿足給定的邊界條件下，求解方程組(14)～(16)可得到位移分量的解。將位移分量代入到式(17)～(22)，可得到各內力素。但由于動水壓力公式比較復雜，如果將動水壓力公式直接代入到方程組中進行求解，計算量將非常大，并且難以解出。

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

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

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

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

分別將靜水壓力與動水壓力作用下罐壁的應力進行疊加，最后便可求解出彈性壁儲罐在地震作用下罐壁的Mises應力。

3 LNG儲罐算例

3.1 LNG儲罐彈性罐壁Mises理論解

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

假設儲罐建造地區為II類場地，設計地震分組為第2組，抗震設防基本烈度為8度。選用EL Centro波和Manual波2組地震波進行時程分析，峰值加速度為0.7 m/s2;EL波的峰值加速度出現時間為2.12 s;Manual的峰值加速度出現時間為 6.26 s，如圖3所示。

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

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

圖4 反應譜法與歐洲規范對比圖Fig.4 Result comparison between response spectrum method and European standard

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

圖5 彈性壁儲罐動水壓力值Fig.5 Hydrodynamic pressure on elastic wall

圖6 彈性壁儲罐靜動合壓力值Fig.6 Hydrodynamic and hydrostatic pressure on elastic wall

圖7為進行曲線擬合精度的對比。可以看出，擬合曲線除在最低點的誤差為4%以外，其余點的誤差均小于3%，最大為2.89%，平均誤差為1%，說明曲線擬合已經達到精度要求。

圖7 計算數值與曲線擬合對比Fig.7 Contrast between numerical result and curve fitting

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

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

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

取單位長度，b=1，由上面動水壓力作用下的彈性罐壁內力公式得到

式中:下標wd表示動水壓力作用下的彎矩，md為動水壓力作用下的膜應力，d代表動水作用。

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

式中:下標wj表示動水壓力作用下的彎矩，mj為動水壓力作用下的膜應力，j代表靜水作用。

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

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

由于自重影響相對于靜水壓力和動水壓力的影響很小，自重影響可以忽略不計。上述彈性罐壁的Mises應力計算未考慮彈性罐壁自重的影響。

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

ABAQUS中的CEL算法是完全的流固耦合算法，分析類型為動態顯示，具有較高的計算精度，并且已經成為檢驗理論簡化解的重要工具。流體構型通過計算流體在歐拉單元中所占的體積分數確定。流體材料可以在歐拉單元之間流動，并與固體單元相互作用。計算模型時需要一個足夠大的歐拉網格，能夠將固體單元可能移動到的所有位置全部包裹在內，由于CEL計算量較大，所以模型采用對稱建模方式如圖8所示，取一半模型進行分析。罐壁單元尺寸為0.8 m，歐拉體單元尺寸為 0.8 m，如圖 9 所示。

圖8 LNG儲罐整體模型Fig.8 Model of LNG tank

圖9 LNG儲罐網格劃分Fig.9 Meshing model of LNG tank

4 計算結果及分析

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

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

從圖中還可以看出，對罐壁Mises應力起主導作用的是靜水壓力，而非動水壓力;彈性儲罐在EI波作用時罐壁最大Mises應力在底部，理論計算數值為 461.14 MPa，仿真數值為 483.5 MPa，最大 Mises應力出現時間為3.5 s滯后于峰值加速度1.38 s。Manual波作用時，罐壁底部最大Mises應力理論計算數值為 461.14 MPa，仿真數值為 483.175 MPa，最大Mises應力出現時間為4.7 s，提前于峰值加速度1.56 s。

在Manual波作用下，彈性儲罐罐壁最大Mises應力出現時間提前于峰值加速度時間其原因如下。分別對EI波和Manual波加速度時程曲線進行傅里葉變換，可以得到加速度頻譜曲線如圖14～15。

圖10 EI波對比0°Fig.10 Comparison results under EI wave at 0°

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

圖12 Manual波對比180°Fig.12 Comparison results under Manual wave at 180°

圖 13 Manual結果 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

通過對頻譜特性［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。

根據文獻［12-14］，儲罐滿液時的自振頻率為

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

5 結論

本文通過對LNG罐壁Mises應力理論計算與有限元仿真結果相對比，得出結論:

1)基于我國反應譜法計算出的彈性壁儲罐的動水壓力在儲液上部的對流壓力值和儲液下部的脈沖壓力值都偏小。

2)通過彈性罐壁Mises應力隨高度變化的趨勢圖中可以看出，對罐壁Mises應力起到主導作用的是靜水壓力對罐壁的作用，而非動水壓力。

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

4)當地震加速度的卓越頻率與儲液罐的自振頻率比較接近時，導致儲液罐出現共振或亞共振狀態時，會造成罐壁最大的Mises應力提前于地震峰值加速度出現。

［1］居榮初，曾心傳.彈性結構與液體的耦聯振動理論［M］.北京:地震出版社，1983:40-48.

［2］SH/T 3026-2005鋼制常壓立式圓筒形儲罐抗震鑒定標準［S］.上海:中國石化出版社，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］胡盈輝，莊茁，由小川.大型儲液罐在地震作用下的附加質量法研究［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］杜顯赫，沈新普，劉應華.預應力LNG儲罐在地震作用下的流固耦合數值模擬［C］//力學與工程應用.鄭州，中國，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］孫建剛，崔利富，張營.全容式LNG儲罐地震響應數值模擬研究［C］//中國國際建筑科技大會論文集.北京，中國，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］.北京:清華大學出版社，1987:137-146.

［9］史潔玉，孔玲軍.MATLAB R2012a超級學習手冊［M］.北京:人民郵電出版社，2013:72-85.

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

［11］翁智遠.圓柱薄殼容器的振動與屈曲［M］.上海:上海科學技術出版社，2007:3-16.

［12］翁智遠，楊建中.圓柱形彈性貯液容器基頻近似計算法［J］.力學季刊，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.