LNG儲罐樁土非線性有限元模擬

趙澤欽，陳團(tuán)海，張 超

中海石油氣電集團(tuán)有限責(zé)任公司，北京 100028

LNG儲罐樁土非線性有限元模擬

趙澤欽，陳團(tuán)海，張 超

中海石油氣電集團(tuán)有限責(zé)任公司，北京 100028

在對LNG儲罐進(jìn)行全罐結(jié)構(gòu)分析和安全校核時，通常通過在承臺底部或樁身施加等效非線性彈簧等方法來簡化模擬土體對上部結(jié)構(gòu)的承載作用，但這種簡化的近似方法無法準(zhǔn)確地反映樁土之間復(fù)雜的相互作用關(guān)系。采用三維非線性有限元軟件ANSYS對LNG大型儲罐的樁基進(jìn)行樁土作用分析，分析時通過三種方法估算了土體的變形模量，同時考慮了初始地應(yīng)力和土體寬度的邊界效應(yīng)，并把計算得到的結(jié)果與現(xiàn)場試樁實測結(jié)果進(jìn)行比較。研究表明，三維非線性有限元模擬結(jié)果與現(xiàn)場實測的荷載-沉降數(shù)據(jù)擬合度較好；三種方法估算土體的變形模量，方法一（基于多孔介質(zhì)理論，由壓縮模量推算變形模量）的偏差相對最大，方法二（通過標(biāo)貫擊數(shù)估算變形模量）的準(zhǔn)確性依賴于所選取經(jīng)驗公式的適用性，方法三（通過參數(shù)反分析法擬合變形模量）在缺乏原位測試數(shù)據(jù)時能得到較好的擬合結(jié)果；樁側(cè)摩阻力總體隨著深度的增加而增大，樁端處附近達(dá)到極限值，樁側(cè)土體的總側(cè)摩阻力占總荷載的69.7%，樁端土體的端阻力占總荷載的30.3%。

樁土作用；有限元模擬；非線性接觸；變形模量；荷載-沉降曲線

LNG大型儲罐的基本結(jié)構(gòu)主要包括樁基、承臺、預(yù)應(yīng)力混凝土外罐、X7Ni9鋼內(nèi)罐、吊頂及頂梁框架等，上部罐體和罐內(nèi)液體的重力荷載通過承臺和樁基傳遞到下方的土體上，因而樁土在各種工況條件下的荷載傳遞和承載能力是LNG大型儲罐設(shè)計和建造過程中主要考慮的問題之一。在對LNG儲罐進(jìn)行全罐結(jié)構(gòu)分析和安全校核時，通常通過在承臺底部或樁身施加等效非線性彈簧等方法來簡化模擬土體對上部結(jié)構(gòu)的承載作用。這種簡化的近似方法無法準(zhǔn)確地反映土體復(fù)雜的本構(gòu)關(guān)系以及樁土之間的非線性作用關(guān)系。文章分析了樁土接觸作用機(jī)理及適用的本構(gòu)模型，使用國際通用有限元軟件ANSYS建立了樁土非線性接觸的三維實體模型，通過三種方法估算土體的變形模量，并與現(xiàn)場試樁實測結(jié)果進(jìn)行比較和分析，為LNG大型儲罐全罐應(yīng)力分析中的樁基分析部分提供了一種新的思路和方法。

1 樁土作用理論

1.1 樁土作用研究方法

樁土交互作用是一種高度非線性接觸問題，可以使用有限元法進(jìn)行數(shù)值求解，分析時需要選擇樁和土體適用的本構(gòu)模型和屈服準(zhǔn)則。對于樁體，一般軟件中都有成熟的鋼筋混凝土模型。土體的本構(gòu)模型則比較復(fù)雜：一種為非線性彈性本構(gòu)模型，如鄧肯-張模型[1]、K-G模型；一種為彈塑性本構(gòu)模型，主要為Dracker-Prager模型[2]（以下簡稱為DP模型）。LNG儲罐樁基所處的土體符合彈塑性本構(gòu)特征，使用DP模型進(jìn)行模擬。

1.2 DP模型

ANSYS軟件中目前能用于近似描述巖土材料的模型只有DP模型。DP模型是理想彈塑性模型，是對Mohr-Coulomb準(zhǔn)則的修正。它主要包含三個參數(shù)：內(nèi)聚力c（Pa或kPa）、摩擦角準(zhǔn)（°）和剪脹角Φf（°），其屈服準(zhǔn)則為：

3βσm+[0.5{s}T[M]{s}]0.5-σy=0 （1）

式中：β為與摩擦角準(zhǔn)有關(guān)的無量綱材料常量；σm為平均應(yīng)力，即土力學(xué)中的靜水壓力，Pa或kPa；{s}為偏應(yīng)力，Pa或kPa；T為矩陣轉(zhuǎn)置符；[M]為Mises屈服準(zhǔn)則中的系數(shù)項，無量綱；σy為與內(nèi)聚力c和摩擦角準(zhǔn)相關(guān)的材料屈服參數(shù)，Pa或kPa。

1.3 樁土荷載傳遞

樁頂在受到向下壓的荷載時，樁身壓縮產(chǎn)生向下的位移，樁土之間發(fā)生相對滑移。當(dāng)施加荷載超過樁側(cè)土體的極限側(cè)摩阻力與樁端土體的極限端阻力之和時，樁體發(fā)生持續(xù)的大位移沉降。

樁土接觸是一個典型的剛性面-柔性面非線性接觸問題，可以引入Goodman無厚度四節(jié)點接觸單元[3]。在ANSYS分析中，選用TARGE170描述樁體的剛性目標(biāo)面，CONTAC173描述樁側(cè)土的柔性接觸面，組成接觸單元，接觸形式為“No separation”（允許滑移和禁止分離），接觸摩擦形式為庫倫摩擦，摩擦系數(shù)的取值：黏土為0.3、粉質(zhì)黏土為0.5、粉砂為0.8[4]。

1.4 土體的變形模量

土體的變形模量E0和壓縮模量Es是判斷土體壓縮性和計算地基壓縮變形的重要指標(biāo)。在計算樁基沉降時，應(yīng)使用土體的變形模量，而實際地基勘察一般不直接提供土體的變形模量。

本文通過三種方法估算土體的變形模量，并比較了三種方法的沉降擬合結(jié)果：

（1） 基于多孔介質(zhì)理論、由Es推導(dǎo)得到E0。當(dāng)假定土體為理想彈性體時，根據(jù)廣義虎克定律，可推導(dǎo)得到E0與Es的關(guān)系為[5]：

式中：μ為土體的泊松比，無量綱；Es、E0的單位為MPa。

然而，實際的土體是應(yīng)力-應(yīng)變關(guān)系極其復(fù)雜的彈塑性材料，其變形模量通常不滿足該式。大量的統(tǒng)計資料表明，E0有時可能是Es的3～5倍，甚至更大。

本文根據(jù)文獻(xiàn)[6]中基于多孔介質(zhì)理論建立的地基土變形模量估算方法，計算得到工區(qū)各層土體的變形模量。

式中：Ψ為土體擾動補(bǔ)償系數(shù)，取值1.49；e為地基土孔隙比。

（2）通過標(biāo)貫擊數(shù)N值估算E0。標(biāo)準(zhǔn)貫入試驗（SPT） 是國內(nèi)外應(yīng)用廣泛的一種現(xiàn)場原位測試技術(shù)，工程中可以通過相應(yīng)的經(jīng)驗公式，擬合得到變形模量與N值的近似函數(shù)關(guān)系。通過文獻(xiàn)調(diào)研，分別使用式（4）中的經(jīng)驗公式[7-11]，計算得到相應(yīng)的E0，代入ANSYS中進(jìn)行計算和比較。

（3）通過參數(shù)反分析法擬合E0。方法（1）、（2）只是對E0的近似估算，估算結(jié)果的準(zhǔn)確性依賴于土體性質(zhì)的相似性和所采用經(jīng)驗公式的適用性。由于DP模型的參數(shù)較少并且通常都能夠獲得，因此在缺乏地基土勘察數(shù)據(jù)時，可采用參數(shù)反分析的方法，通過擬合樁頂?shù)淖罱K沉降位移，來估算土體的實際E0。

2 試樁概況

試樁工區(qū)位于天津市濱海新區(qū)，樁型為鉆孔灌注樁，樁徑1.2 m、樁長56 m，設(shè)計單樁豎向抗壓極限承載力為9 600 kN。

2.1 單樁靜載試驗

儲罐建造的前期，需要在工程樁施工前進(jìn)行單樁靜載試驗。根據(jù)JGJ 106-2014《建筑樁基檢測技術(shù)規(guī)范》進(jìn)行單樁豎向抗壓靜載試驗，采用如圖1所示的壓重平臺反力裝置，在樁頂豎向逐級施加等量荷載（逐級荷載為設(shè)計單樁豎向抗壓極限承載力9600 kN的1/10）并維持，最終加載至9 600 kN。

圖1 單樁豎向抗壓靜載試驗壓重平臺反力裝置

2.2 樁基土體參數(shù)

工區(qū)樁基土體包含多套土層，樁身上部土層主要為沖填土、黏土和粉質(zhì)黏土，下部土層主要為粉土、粉質(zhì)黏土和粉砂。樁端所在持力層為⑩4粉砂層，該土層呈密實狀態(tài)，工程地質(zhì)性質(zhì)好。

試樁工區(qū)土體各土層的參數(shù)見表1。

3 模型與計算

3.1 ANSYS模型

在ANSYS中建立如圖2所示的三維樁土實體模型（1/4模型），土體為長方體，寬為樁徑的15倍，厚約64 m，包含12層土層。

表1 試樁工區(qū)土體參數(shù)

圖2 樁土三維實體模型

采用映射網(wǎng)格劃分，見圖3。模型的邊界條件設(shè)置為：土體兩側(cè)面分別固定X方向和Y方向上的位移UX、UY，樁體和土體的對稱面采用對稱邊界條件，土體底面固定豎直方向上的位移UZ，樁端面與土層接觸面的自由度耦合，樁頂面施加豎直向下的荷載。

圖3 樁土網(wǎng)格模型

3.2 初始地應(yīng)力場

在樁澆筑之前，土體早已在自身重力的作用下產(chǎn)生了較大的豎向位移和初始應(yīng)力，后者直接影響樁土之間的側(cè)摩阻力。因此，在分析樁沉降時，需要在去除土體初始位移的同時保留初始地應(yīng)力。

ANSYS中可以使用inistate命令模擬初始地應(yīng)力的自平衡，平衡前后的初始位移場見圖4。在自重與初始地應(yīng)力平衡前，土體在自身重力作用下發(fā)生了較大的豎向位移，頂部土層位移最大，達(dá)209 mm，遠(yuǎn)大于試樁的最終沉降量。平衡后，土體的最大位移為7.86×10-3mm，原始位移場基本清零。圖5說明，在消除土體原始位移場的同時，也很好地保留了原始地應(yīng)力場。

圖4 初始位移場/m

3.3 土體寬度的選擇

在建立三維樁土模型時，需要確定合適的土體寬度。過小的土體寬度會產(chǎn)生較大的邊界效應(yīng)，影響計算結(jié)果的準(zhǔn)確性，過大的土體寬度則會加大求解的計算量。分別采用5、10、15、20、30倍樁徑的土體寬度，計算樁端的最終沉降量。如圖6所示，隨著土體寬度的增大，邊界效應(yīng)減小，當(dāng)土體寬度為樁徑的15倍及以上時，最終沉降量收斂于同一值，說明此時邊界效應(yīng)已經(jīng)小到可以忽略不計。因此，在后文的算例中，均采用15倍樁徑的土體模型。

圖5 初始地應(yīng)力場/Pa

圖6 樁端最終沉降量-土體寬度與樁徑之比關(guān)系曲線

3.4 三種變形模量估算方法的比較

在1.4節(jié)中介紹了三種估算土體變形模量的方法，在保持其他參數(shù)不變的條件下，比較和分析了這三種方法的計算結(jié)果。

（1）方法一，基于多孔介質(zhì)理論，由ES推導(dǎo)得到E0。前期的勘察工作已經(jīng)獲得了各土層的天然孔隙比、泊松比、壓縮模量等參數(shù)，根據(jù)式（3）估算得到了各土層的變形模量，代入模型計算得到圖7所示的荷載-沉降曲線。

圖7 由方法一得到的沉降-荷載曲線

（2）方法二，通過標(biāo)貫擊數(shù)N值估算E0。根據(jù)式（4）中的各經(jīng)驗分公式，分別估算得到各土層的變形模量（見圖8），并進(jìn)一步計算得到沉降-荷載曲線（見圖9）。

圖8 由方法二各經(jīng)驗分公式計算得到的土變形模量

圖9 由方法2得到的沉降-荷載曲線

可以看出，當(dāng)施加荷載較低時，荷載主要通過樁身上部的側(cè)摩阻力傳遞到上部樁周的土體，而上部土層的變形模量差別不大，因此各經(jīng)驗公式計算得到的沉降量差別不大；當(dāng)施加荷載＞4 000 kN時，下部土體側(cè)摩阻力逐漸發(fā)揮作用，荷載也傳到樁端土體，由于下部土層尤其是持力層土層的變形模量差別較大，因此計算的沉降量也產(chǎn)生較大偏差。估算的變形模量越大，計算的沉降量越接近于真實值，因此在后續(xù)的方法比較中，選擇式（4）中的分公式（4-5）的計算結(jié)果作為方法二的典型算例。

（3）方法三，通過參數(shù)反分析法擬合E0。工程上一般難以獲得土體的真實變形模量，因此當(dāng)僅有室內(nèi)測得的壓縮模量時，可以通過倍數(shù)放大壓縮模量的方法擬合沉降-荷載曲線，從而估算變形模量。分別假定1、1.5、2、3倍壓縮模量作為變形模量，計算得到如圖10所示的荷載-沉降曲線。當(dāng)取1倍壓縮模量作為變形模型模量時，計算沉降偏大，說明工區(qū)土體的實際變形模量應(yīng)大于壓縮模量；當(dāng)取2倍壓縮模量時，計算結(jié)果擬合最好，認(rèn)為工區(qū)土體的實際變形模量約為壓縮模量的2倍。

圖10 由方法三得到的荷載-沉降曲線

（4）三種估算方法的比較。如圖11所示，選擇由這三種方法擬合得最好的算例進(jìn)行比較：方法一（基于多孔介質(zhì)理論、由Es推導(dǎo)E0）在荷載與沉降的早期與實測數(shù)據(jù)的吻合度較高，但隨著荷載和沉降量的增大，開始發(fā)生偏移；方法二（通過標(biāo)貫擊數(shù)N值估算E0）、方法三（通過參數(shù)反分析法擬合E0）都是通過放大壓縮模量來估算變形模量的，擬合結(jié)果相近，在荷載與沉降的早期不如方法一，但后期的吻合度更高。

3.5 結(jié)果分析

從三種估算方法的比較中可以看出，方法二和方法三的擬合度更高，擬合結(jié)果差別不大，而方法二依賴于所選取經(jīng)驗公式的適用性，因此選取方法三的擬合結(jié)果做進(jìn)一步的分析。圖12為土體在最大荷載下（9 600 kN） 的最終沉降云圖，在樁頂?shù)呢Q向荷載下，樁周土體產(chǎn)生以樁為中心的沉降盆：橫向上，離樁越近，土體沉降越大，土體側(cè)邊界的沉降為0；縱向上，第2、3層土體的變形模量最小，壓縮性高，因而沉降量較大。樁端處的土體提供端阻力，受到向下的壓應(yīng)力，因而相比附近土體也有明顯的下沉。

圖11 三種估算方法計算結(jié)果對比

圖12 土體最終沉降云圖/m

圖13為樁體的最終位移云圖，樁身各點的位移隨著深度的增加逐漸減小，樁頂處沉降量最大，為14.85 mm，樁端處沉降最小，為5.17 mm。可以看出，5.17 mm為樁端處土體受壓產(chǎn)生的變形，而樁頂?shù)目偝两禐?5.87 mm，其差值為樁土接觸面的相對滑移和樁體自身的變形。這說明對于這種超長樁，其在上部載荷作用下的變形沉降，不僅來自樁端持力土層的變形，還包括樁土的相對滑移和樁體自身的變形，其分擔(dān)比例取決于樁側(cè)土和樁端土的變形模量和壓縮性。

圖14為最大荷載時樁頂?shù)綐抖说膫?cè)摩阻力分布曲線，可以看出，側(cè)摩阻力總體上隨著深度（樁側(cè)接觸壓力，見圖15）的增加而增大；樁頂處土體受到向下的剪應(yīng)力，發(fā)生明顯的相對滑動（見圖16）；樁身下部側(cè)摩阻力隨著深度的增加而增大，樁端處達(dá)到極限值。計算得到樁側(cè)土體的總側(cè)摩阻力為6 785 kN，占總荷載的69.7%；樁端土體的端阻力為2 914 kN，占總荷載的30.3%。

圖13 樁最終位移云圖/m

圖14 樁側(cè)摩阻力-距樁頂距離曲線（9 600 kN）

圖15 樁側(cè)接觸壓力-距樁頂距離曲線（9 600 kN）

4 結(jié)論

（1）使用ANSYS建立了樁土非線性接觸的三維實體模型，考慮了初始地應(yīng)力和土層寬度的邊界效應(yīng)，模擬結(jié)果與現(xiàn)場實測的荷載-沉降數(shù)據(jù)擬合度較好，說明了該模型的準(zhǔn)確性。

圖16 樁土相對滑移距離-距樁頂距離曲線（9 600 kN）

（2）使用三種方法估算土體的變形模量，并進(jìn)行比較：方法一（基于多孔介質(zhì)理論，由Es推導(dǎo)E0）的偏差相對最大；方法二（通過標(biāo)貫擊數(shù)N值估算E0）依賴于所選取經(jīng)驗公式的適用性；方法3（通過參數(shù)反分析法擬合E0）在缺乏原位測試數(shù)據(jù)時，能得到較好的擬合結(jié)果。

（3）側(cè)摩阻力總體隨著深度的增加而增大，樁端處附近達(dá)到極限值。樁側(cè)土體的總側(cè)摩阻力占總荷載的69.7%，樁端土體的端阻力占總荷載的30.3%。

[1]DUNCAN J M，CHANG C Y.Nonlinear analysis of stress and strain in soils[J].Asce Soil Mechanics&Foundation Division Journal，1970，96（5）：1 629-1 653.

[2]DRUCKER D C，PRAGER W J.Soil mechanics and plastic analysis or limit design[J].Q.appl.math，1952，10（2）：157-165.

[3]GOODMAN R E，TAYLOR R L，BREKKE T L.A model for the mechanics of jointed rock[J].Journal of Soil Mechanics& Foundations Div，1968，94：637-660.

[4]許宏發(fā)，吳華杰.樁土接觸面單元參數(shù)分析[J].探礦工程，2002（5）：10-12.

[5]陳希哲.土力學(xué)地基基礎(chǔ)[M].第四版.北京：清華大學(xué)出版社，2004：210-215.

[6]梁發(fā)云.基于多孔介質(zhì)理論的地基土變形模量估算方法[J].巖土力學(xué)，2004，25（7）：1 147-1 150.

[7]郭淋，王春艷，張飛，等.標(biāo)貫試驗N值與土體物理力學(xué)參數(shù)的相關(guān)性分析[J].安全與環(huán)境工程，2012，19（4）：148-152.

[8]徐光大，徐光黎，李俊杰.日本標(biāo)準(zhǔn)貫入試驗方法及其N值在巖土工程中的應(yīng)用[J].安全與環(huán)境工程，2011，18（4）：33-38.

[9]李修建，李松.標(biāo)準(zhǔn)貫入與載荷試驗數(shù)據(jù)的對比分析[J].山西建筑，2011，37（10）：58-59.

[10]祝龍根，劉利民，耿乃興.地基基礎(chǔ)測試新技術(shù)[M].北京：機(jī)械工業(yè)出版社，2002：175-176.

[11]楊文衛(wèi)，岳中琦.世界各地標(biāo)準(zhǔn)貫入試驗比較和共同問題[J].工程勘察，2008（1）：5-15.

Finite element simulation of non-linear interaction between pile foundation of L NG tank and soil

ZHAO Zeqin，CHEN Tuanhai，ZHANG Chao
CNOOC Gas&Power Group Co.，Ltd.，Beijing 100028，China

This paper presents a three dimensional finite element stress analysis method to simulate the non-linear interactions between pile foundation of LNG storage tank and soil.A three dimensional solid model is built by ANSYS and the deformation modulus is calculated by three different methods.The initial ground stress and the boundary effect of soil width are also considered.Compared with pile testing data，the results show that the simulated load-settlement curve is in agreement with the actualdata；The first method evaluating the deformation modulus（based on porous medium theory to calculate deformation modulus from compression modulus） has the maximum deviation，the accuracy of the second method（based on SPT to calculation deformation modulus） relies on the applicability of selected empirical formula，the third method（based on parameter back analysis to fit deformation modulus） has good fitting result in lack of in-situ test data.The side friction of the pile increases generally with the depth along the pile and reaches the maximum value at the pile end；the totalside friction is 69.7%of totalload and the pile end resistance is 30.3%of totalload.

interactions between pile foundation and soil；finite element simulation；non-linear contact；deformation modulus；load-settlement curve

10.3969/j.issn.1001-2206.2017.03.004

趙澤欽（1989-），男，廣東汕頭人，2015年畢業(yè)于北京大學(xué)力學(xué)專業(yè)，碩士，現(xiàn)主要從事LNG儲罐結(jié)構(gòu)分析方向的研究工作。

2017-01-03

中國海洋石油總公司科技課題“大型LNG儲罐結(jié)構(gòu)設(shè)計軟件開發(fā)”（CNOOC-KJ 135 KJXM QD 2016-002）

Email：zhaozq4@cnooc.com.cn