軟土場地大型LNG儲罐考慮樁土相互作用的地震響應分析

劉 帥， 翁大根， 張瑞甫

(同濟大學 土木工程防災國家重點實驗室，上海 200092)

LNG(液化天然氣)是一種無色無味、無毒、無腐蝕性的優質高效能源，各國正逐步將其作為一種低排放量的清潔燃料加以推廣。用以加工和存儲液化天然氣的儲罐的建造量也因此不斷增多，尤其是沿海的接收站LNG儲罐項目正向著大型、超大型容量方向發展(最大可達20×104m3)。對于如此巨大而重要的結構物來說一旦在地震中發生破壞不僅會造成巨大的經濟損失，而且會導致爆炸、火災等災難性后果，所以其抗震安全性能問題引起了國內外越來越多研究人員的重視。對儲罐動力響應分析具有里程碑意義的研究是Housner[1]提出的剛性罐抗震分析理論，但實踐證明該理論未考慮罐壁柔性，使得儲罐的抗震安全性偏于非保守，所以之后的研究對罐壁的彈性進行了考慮。Edwards[2]首次采用有限元法對錨固儲罐流固耦合地震響應進行了數值模擬研究；Haroun[3-4]采用有限單元法(罐壁和液面區域)和邊界條件法(儲液區域)相結合的計算方法將儲液罐的流固耦合系統分解為相互獨立的液殼系統和剛性罐對流系統進行分析。近年來，隨著計算機技術的飛速發展，特別是一些大型通用有限元軟件的出現，更是為儲罐地震響應的數值模擬分析提供了很好的技術平臺。Ozdemir等[5]在LS-DYNA中采用ALE(任意拉格朗日歐拉法)方法實現了非錨固罐考慮幾何、材料非線性和土-罐接觸非線性條件下儲罐流固耦合地震響應的數值模擬,孫建剛等[6]在ADINA中建立了15×104m3立式儲罐的三維實體有限元模型，探討了浮頂質量、罐壁厚度等因素對儲罐晃動頻率和流固耦合頻率的影響，張營[7]在ADINA中建立了16×104m3LNG全容罐的有限元模型，對比分析了加速度和位移加載條件下儲罐的地震響應，李思[8]在ANSYS中對16×104m3LNG全容罐地震作用下的動力響應進行了分析，得出了最不利部位對儲罐結構地震響應的影響。

雖然目前儲罐的抗震研究工作已經取得了很多突破性進展，但多數都是基于剛性地基假定而進行的，對于考慮樁土相互作用對儲罐地震響應影響方面的研究卻不多見。在軟土地區，大型LNG儲罐通常采用均勻、對稱的高樁基礎形式來滿足儲罐底部散熱條件及改善地基基礎特性，高樁上端嵌固于罐體底板當中，下端埋置于半無限土體當中，當地震發生時，基巖波先是傳經樁土地基后到達上部罐體，然后返回至地基當中，最終，在地基的輻射阻尼條件下逐漸消失。在整個過程當中，樁-土-罐相互作用相互影響，形成一個整體動力系統。如果在分析過程中假設地基為剛性體，直接對儲罐輸入自由場地表波激勵，這將未能考慮樁-土-罐之間相互作用的影響，也未能進一步考慮場地類別、樁基布置形式及埋深等因素的影響。所以，為了確保大型LNG儲罐的整體地震安全性，本文將利用大型通用有限元軟件ANSYS建立起樁-土-LNG儲罐的三維整體有限元模型，對考慮樁土相互作用對LNG儲罐地震響應的影響進行分析，為工程實踐提供參考。

1 儲罐流固耦合有限單元法

1.1 基本力學方程

1.1.1 儲液域

進行流體動力分析時基本假定：液體考慮為無粘滯性的、無旋的、不可壓縮的理想流體?；炯俣l件下理想流體在柱坐標系(r,θ,z)中應滿足拉普拉斯方程[9]：

Δ2φ=0

(1)

其中φ(r,θ,z,t)是速度勢，且有：

(2)

1.1.2 罐壁域

在罐壁區域任意一點應滿足動力平衡微分方程(張量形式)：

(3)

其中σij為罐壁應力分量，ui為罐壁位移分量，fi為罐壁體積力分量，ρs為罐壁質量密度。

1.1.3 邊界條件

(1)在流體和罐壁的接觸面上流體的徑向速度和罐壁的徑向速度應相等,則有:

(4)

其中wg和w分別為罐壁隨地面一起運動(剛性脈沖)和自身彈性變形(柔性脈沖)的徑向位移分量，R為罐壁半徑。

(2) 流體自由表面由于擾動而發生晃動時，在重力的作用下總是趨向于恢復到靜止的平衡位置(形成儲液的對流晃動)，所以在液體的自由表面(z=H)應滿足自由面波條件(假設液體表面壓力為0)：

(5)

其中f為流體表面偏離平衡位置的垂直位移，g為重力加速度。

1.2 基于位移形式的流固耦合有限單元法

考慮到儲液流體運動方程和彈性罐體運動方程的相似性，同樣以流體位移為待求未知量，建立起流體基于位移的有限元格式，這樣流體與罐體均為位移計算模式，可以方便進行標準有限元格式的計算，即便是復雜的耦合界面問題也能夠較為容易地處理。采用流體和罐體的統一的基于位移的有限單元法后，儲液罐流固耦合系統在地震作用下基于變分原理的整體動力有限元方程可以表示為：

(6)

由于理想流體幾乎不存在抗剪強度，所以流體和固體僅在耦合界面的正法向發生相互作用，在切向上不存在剪力的傳遞。在采用有限單元法處理儲液罐流固耦合問題時，應將流體與罐體在交界面處的法向節點位移約束，切向位移仍然保留。

2 土動力本構關系

土動力本構關系是了解在動力荷載作用下土體動力響應的基礎，也是利用數值計算方法進行動力分析的前提條件。目前廣泛采用的能夠考慮土體非線性動力特性的計算模型大致有三種，一種是等效線性化模型[10]，該模型簡化了土的黏性和塑性能量損耗的復雜本質，將應力-應變關系用隨動應變幅值變化的等效剪切模量與阻尼比來表示，而后通過多次計算，反復迭代獲得近似的非線性解答；另一種是在非線性骨架曲線和Masing法則基礎上發展起來的經驗模型，例如Hardin-drnevich雙曲線模型[11]、Martin-seed-davidenvic 模型[12]以及陳國興等[13]提出的修正Davidenvic模型(在ABAQUS軟件平臺上成功地利用子程序開發實現了該模型)，最后一種是利用一系列具有特定物理關系的物理元件按某種組合來構造總體應力-應變關系的模型，例如Iwan模型。在這三種模型當中，由于等效線性化方法在每次迭代計算中都屬于彈性分析，不存在收斂問題，而且直觀簡單，所以在場地土層地震響應分析和土-結構相互作用分析時得到廣泛應用，本文在進行樁-土-罐相互作用地震響應分析時采用了土的等效線性化模型，并在ANSYS軟件平臺上實現了該模型的算法。

3 邊界條件

用有限單元法分析樁-土-罐動力相互作用時，需要從半無限地球介質中切取有限介質區域進行分析，在切取的邊界上需建立人工邊界以模擬連續介質的輻射阻尼，該輻射阻尼可以耗散從上部結構物反射過來的地震波能量。為減少人工邊界波動反射帶來的誤差，國內外學者[14-15]做了許多有益的工作，其中最簡單的處理方法是避開邊界條件問題，直接采用自由邊界，即取人工邊界至結構物的距離足夠遠，以使在邊界處波的衰減趨近于零。樓夢麟等[16]通過采用有限單元法分析了在水平地震波激振作用下均勻土層中有限土域的取值范圍對計算精度的影響，結果得出確定有限土層合理取值范圍的決定因素是長深比，當阻尼比為10%時，選取土層的長深比大于7之后，采用自由邊界可以將計算結果與精確解的誤差控制在5%以內。

4 樁-土-罐相互作用地震響應算例分析

4.1 罐體幾何尺寸及材料屬性

某LNG接收站項目工程采用特大型LNG全容罐，儲罐的有效容積為16×104m3，外罐內徑為82 m，外罐高49.9 m，壁厚上下均為0.8 m，外罐頂厚0.35 m，圓弧半徑為40 m；內罐內徑為80 m，內罐高35.43 m，儲液高度為33 m，壁厚上下均為0.028 6 m；內外罐之間采用保溫材料珍珠巖填充，厚度為1 m；混凝土底板的直徑為86.6 m，厚度為1.2 m(中間0.9 m)；采用摩擦型樁基礎，樁長24 m，半徑為0.6 m，共360根。特大型LNG儲罐的斷面示意圖及樁基平面布置圖如圖1、2所示。

LNG儲罐的材料屬性見表1。

表1 LNG儲罐材料屬性

圖1 LNG儲罐斷面圖

圖2 樁基平面布置圖

4.2 場地條件及土的動力本構模型

由本工程《場地地震安全性評價報告》得知，該場地的等效剪切波速為56 m/s，小于150 m/s，且覆蓋層厚度為100 m,大于80 m，判定為IV類軟土場地。土體共分為7層， 各層土的密度和現場測得的土層剪切波速度如表2所示；由共振柱試驗測得的動剪切模量比、阻尼比與動剪應變幅值的關系如圖3所示。

表2 土層物理力學參數

圖3 土剪切模量比和阻尼比與動剪應變幅值的關系

土的動力本構關系采用等效線性化模型，并利用有限元軟件ANSYS 程序參數化設計語言(APDL)的功能將該模型在程序中實現。采用參數化編程實現該模型循環迭代的邏輯流程如圖4所示。

圖4 參數化語言實現等效線性化模型流程圖

4.3 單元類型及有限元模型

有限元軟件ANSYS擁有非常完善的單元庫，針對不同的材料和分析類型，單元庫提供了豐富的單元類型。依據樁土和儲罐結構本身的特點，本文中樁-土-罐整體有限元模型所選用的單元類型如表3所示。

儲液和內罐壁之間的相互作用是通過適當的耦合兩個單元之間接觸面上的節點來實現的，這就要求液體單元劃分出來的節點應和內罐壁殼單元劃分出來的節點相對應，而且此時還應將笛卡爾直角坐標系變換至柱坐標系下，變換后的X軸為圓周法向方向，Y軸為圓周切向方向，Z軸方向保持不變，然后在法向相同位置的節點進行耦合，保證儲液不能與內罐壁分離，但在切向可以相對運動，并對內罐壁施加液動壓力。液體單元和內罐壁殼單元之間的耦合作用如圖5所示。

表3 有限元模型單元類型

圖5 流固耦合作用示意圖

用有限單元法分析樁-土-罐動力相互作用時所切取的有限土層的范圍為半徑750m的圓柱體，該范圍滿足自由邊界長深比大于7的條件。鑒于樁土動力相互作用的復雜性，假設樁土之間不出現滑移和分離等相對位移現象，變形過程中樁土接觸點的位移始終滿足位移協調條件。由此建立起的樁-土-罐整體有限元模型如圖6所示。

圖6 樁-土-罐整體有限元模型

4.4 地震波選取

為了能夠分析IV軟土場地條件下考慮樁土相互作用對LNG儲罐地震響應的影響，本文對以下兩種模型狀態進行了地震響應分析：

狀態1(ST1)：樁-土-罐相互作用分析模型(基巖地震動波輸入)

在進行樁-土-罐相互作用地震響應分析時，依據《建筑抗震設計規范》[17]中I類場地(特征周期0.2 s)標準反應譜并通過傅里葉變換擬合了三條人工地震動時程，將其作為基巖地震動輸入。原理上基巖地震動是從基巖頂面輸入，而在有限元軟件ANSYS中則是通過對樁-土-罐結構體系施加等價的外部激勵來實現。為了考慮LNG儲罐在SSE(安全停運地震)工況下的地震響應，將加速度峰值調整至0.3 g。其中一條歸一化加速度時程曲線和其傅里葉譜如圖7、8所示，3條地震波反應譜曲線如圖9所示。

圖7 PC1加速度時程曲線

圖8 PC1加速度傅里葉譜

圖9 歸一化加速度反應譜(阻尼比5%)

狀態2(ST2)：假設地基為剛性的儲罐分析模型ST2(自由場地表波輸入，時程曲線如圖10)。

地震波從基巖傳至地表時，經過多次反射和折射，其頻譜特征發生了很大的變化。圖11給出了本文所選基巖地震動波在傳經軟土場地到達地表后其加速度反應譜變化情況，從圖中可知，自由場地表波的長周期成分較基巖地震動波獲得顯著的增強。自由場地表波的特征周期在1s左右，這與Ⅳ類軟土場地條件相符。

圖10 PC1地表波加速度時程曲線

圖11 地表波歸一化加速度反應譜

4.5 計算結果與分析

在進行樁-土-罐相互作用地震響應分析之前，首先針對計算結構進行了模態分析，兩種結構狀態下模態分析結果的對比如表4所示。

表4 模態分析結果對比

4.5.1 儲罐基底剪力、傾覆力矩及晃動波高

IV類軟土場地條件下考慮樁土相互作用時(ST1)和假設地基為剛性時(ST2)，該LNG儲罐在SSE工況下的基底剪力、傾覆力矩及晃動波高的時程曲線如圖12所示，峰值和最大值如表5所示。

由表5可知，在IV類軟土場地條件下考慮樁土相互作用以后，該LNG儲罐的基底剪力和傾覆力矩的峰值分別減小了24%和7%；最大晃動波高增大了15%。

兩種狀態下儲液的對流晃動形式相同，圖13為在ST1狀態下不同時刻儲液對流晃動的云圖。

圖12 基底剪力、傾覆力矩及晃動波高時程曲線對比(PC1)

表5 基底剪力和傾覆力矩峰值及最大晃動波高對比

圖13 儲液對流晃動云圖(5倍真實尺寸)

4.5.2 內罐壁等效應力

通過對比得知，兩種狀態下內罐壁Mises等效應力值的分布相似，即在地震作用方向均由一邊向另一邊逐漸增大，應力較大值也均出現在內罐壁與地震作用方向成零度角處，這和儲液對流晃動的形式是一致的。圖14為在ST1狀態下不同時刻內罐壁Mises等效應力云圖。

圖15 內罐Mises應力峰值沿罐高分布

圖15為在兩種狀態下內罐壁Mises等效應力峰值沿罐高的分布曲線(內罐壁與地震作用方向成零度角處)。從圖中可知，內罐壁等效應力峰值先是在頂部較小，然后沿罐高向底部逐漸增大，直至接近罐底部位時出現內罐壁等效應力峰值最大值，然后又開始減小至罐底，峰值最大值均出現在離罐底約3 m處，未考慮樁土相互作用時為329 MPa，考慮樁土相互作用后有稍微減小，為323 MPa。由美國石油學會(API)標準可知，在低溫狀態下9%Ni鋼的屈服應力達到515 MPa，所以在0.3 g地震作用下內罐均未出現屈服，不會發生破壞。

4.5.3 保溫層對儲罐地震響應的影響

內外罐之間采用珍珠巖保溫材料填充，在數次循環荷載作用下保溫材料會壓縮變得密實，剛度也因此會增大，這將對內外罐之間的動力相互作用產生影響。為了分析這一影響，本文分別對保溫材料彈性模量放大為5倍和10倍狀態下罐體的動力響應進行了分析，分析結果如表6所示。

狀態1(ST1)：原結構;

狀態3(ST3)：5倍保溫材料彈性模量;

狀態4(ST4)：10倍保溫材料彈性模量。

由表6可知，相比于狀態1，狀態3條件下的基底剪力峰值減小0.6%，傾覆力矩峰值增大3%，最大晃動波高減小0.8%；狀態4條件下的基底剪力峰值減小1.2%，傾覆力矩峰值增大4%，最大晃動波高減小0.8%，這說明隨著保溫層剛度在一定范圍內的增加，基底剪力峰值、傾覆力矩峰值和最大晃動波高變化較小。

圖16 內罐Mises應力峰值沿罐高分布

表6 基底剪力和傾覆力矩峰值及最大晃動波高對比

圖16為三種狀態下內罐壁Mises等效應力峰值沿罐高的分布曲線，由圖可知，內罐壁等效應力峰值隨著保溫層剛度的增大有明顯的減小，相比于狀態1，狀態3條件下內罐壁等效應力峰值最大值減小41%(189 MPa)，狀態4條件下內罐壁等效應力峰值最大值減小60%(129 MPa)。這說明隨著保溫層剛度在一定范圍內的增加，內外罐之間的動力相互作用越來越明顯，部分液動壓力可以通過保溫材料傳遞至外罐，由混凝土外罐來承擔，對內罐起到了保護作用。

5 結 論

(1) 該LNG儲罐的基底剪力和傾覆力矩峰值在IV類軟土場地條件下考慮樁土相互作用以后分別有24%和7%的減小，故通常在以罐底基底剪力和力矩為控制條件進行儲罐抗震設計時將地基視為剛體是偏于保守的。

(2) 雖然考慮樁土相互作用以后該LNG儲罐基底剪力和力矩峰值均有減小，但其液面最大晃動波高卻有15%的增大，這是儲液的對流分量在IV類軟土場地條件下考慮樁土相互作用以后被放大的結果。所以，對于波高比較敏感的LNG儲罐在以基底剪力和力矩為控制條件進行抗震設計時應重視考慮樁土相互作用對波高的影響，否則將有可能造成儲液的溢出或灌頂的破壞。

(3) 考慮樁土相互作用以后內罐壁Mises等效應力峰值最大值有稍微減小，且在兩種狀態下最大值均出現在離罐底約3 m處。

(4) 保溫層對儲罐的地震響應具有重要的影響，隨著保溫層剛度在一定范圍內的增加，LNG儲罐的基底剪力峰值、傾覆力矩峰值和最大晃動波高變化較小，但內罐壁等效應力峰值最大值減小明顯，這對于防止罕遇地震作用下內罐底部發生屈曲破壞是有利的。

參 考 文 獻

[1]Housner G W. Dynamic pressures on accelerated fluid containers [J]. Bulletin of the Seismological Society of America, 1957, 47(1):15-35.

[2]Edwards N W. A Procedure for dynamic analysis of thin walled cylindrical liquid storage tanks subjected to lateral ground motions[D].Michigan: University of Michigan-Ann Arbor, 1969.

[3]Haroun M A. Dynamic analysis of liquid storage tanks[R]. 1981010009051, Washington: National Science Foundation, 1980.

[4]Haroun M A. Vibration studies and tests of liquid storage tanks[J]. Earthquake Engineering and Structure Dynamics, 1983,11(2): 179-206.

[5]Ozdemir Z, Souli M, Fahjan Y M. Application of nonlinear fluid structure interaction methods to seismic analysis of anchored and unanchored tanks [J]. Engineering Structures, 2010, 32(2): 409-423.

[6]孫建剛,崔利富,趙長軍,等. 15×104m3立式儲罐隔震設計分析[J]. 地震工程與工程振動,2010, 30(4): 153-158.

SUN Jian-gang, CUI Li-fu, ZHAO Chang-jun, et al. Design and analysis of seismic isolation for 15×104m3vertical storage tank[J]. Journal of Earthquake Engineering and Engineering Vibration,2010, 30(4):153-158.

[7]張營. 大型全容式LNG儲罐地震響應數值模擬研究[D]. 大慶:東北石油大學，2011.

[8]李思. 全容式LNG儲罐的地震響應分析[D].天津:天津大學，2010.

[9]居榮初,曾心傳. 彈性結構與液體的耦聯振動理論[M]. 北京: 地震出版社, 1983.

[10]齊文浩, 薄景山. 土層地震反應等效線性化方法綜述[J].世界地震工程, 2007, 23(4): 221-226.

QI Wen-hao, BO Jing-shan. Summarization on equivalent linear method of seismic responses for soil layers[J].World Earthquake Engineering, 2007, 23(4): 221-226.

[11]Hardin B O, Drnevich V P. Shear modulus and damping in soil: design equations and curves[J]. Journal of the Soil Mechanics and Foundation Engineering Division, ASCE, 1972, 98(7): 667-692.

[12]Martin P P, Seed H B. One dimensional dynamic ground response analysis[J]. Journal of Geotechnical Engineering, ASCE, 1982, 108(7): 935-952.

[13]陳國興, 莊海洋. 基于Davidenkov 骨架曲線的土體動力本構關系及其參數研究[J]. 巖土工程學報, 2005, 27(8): 860-864.

CHEN Guo-xing, ZHUANG Hai-yang. Developed nonlinear dynamic constitutive relations of soils based on Davidenkov skeleton curve[J]. Chinese Journal of Geotechnical Engineering, 2005, 27(8): 860-864.

[14]劉晶波, 谷音, 杜義欣. 一致粘彈性人工邊界及粘彈性邊界單元[J]. 巖土工程學報, 2006,28(9) :1070-1075.

LIU Jing-bo,GU Yin,DU Yi-xin. Consistent viscous-spring artificial boundaries and viscous-spring boundary elements[J]. Chinese Journal of Geotechnical Engineering, 2006, 28 (9) :1070-1075.

[15]Lysmer J, Kuhlemeyer R L. Finite dynamic model for infinite media[J]. Journal of the Engineering Mechanics, ASCE,1969,95(4):859-877.

[16]樓夢麟, 潘旦光, 范立礎. 土層地震反應分析中側向人工邊界的影響[J]. 同濟大學學報, 2003, 31(7):757-761.

LOU Meng-lin,PAN Dan-guang,FAN Li-chu. Effect of vertical artificial boundary on seismic response of soil layer[J].Journal of Tongji University, 2003, 31(7):757-761.

[17]GB50011-2010.建筑抗震設計規范[S]. 北京: 中國建筑工業出版社,2010.