LNG储罐桩土非线性有限元模拟

赵泽钦，陈团海，张 超

中海石油气电集团有限责任公司，北京 100028

LNG储罐桩土非线性有限元模拟

赵泽钦，陈团海，张 超

中海石油气电集团有限责任公司，北京 100028

在对LNG储罐进行全罐结构分析和安全校核时，通常通过在承台底部或桩身施加等效非线性弹簧等方法来简化模拟土体对上部结构的承载作用，但这种简化的近似方法无法准确地反映桩土之间复杂的相互作用关系。采用三维非线性有限元软件ANSYS对LNG大型储罐的桩基进行桩土作用分析，分析时通过三种方法估算了土体的变形模量，同时考虑了初始地应力和土体宽度的边界效应，并把计算得到的结果与现场试桩实测结果进行比较。研究表明，三维非线性有限元模拟结果与现场实测的荷载-沉降数据拟合度较好；三种方法估算土体的变形模量，方法一（基于多孔介质理论，由压缩模量推算变形模量）的偏差相对最大，方法二（通过标贯击数估算变形模量）的准确性依赖于所选取经验公式的适用性，方法三（通过参数反分析法拟合变形模量）在缺乏原位测试数据时能得到较好的拟合结果；桩侧摩阻力总体随着深度的增加而增大，桩端处附近达到极限值，桩侧土体的总侧摩阻力占总荷载的69.7%，桩端土体的端阻力占总荷载的30.3%。

桩土作用；有限元模拟；非线性接触；变形模量；荷载-沉降曲线

LNG大型储罐的基本结构主要包括桩基、承台、预应力混凝土外罐、X7Ni9钢内罐、吊顶及顶梁框架等，上部罐体和罐内液体的重力荷载通过承台和桩基传递到下方的土体上，因而桩土在各种工况条件下的荷载传递和承载能力是LNG大型储罐设计和建造过程中主要考虑的问题之一。在对LNG储罐进行全罐结构分析和安全校核时，通常通过在承台底部或桩身施加等效非线性弹簧等方法来简化模拟土体对上部结构的承载作用。这种简化的近似方法无法准确地反映土体复杂的本构关系以及桩土之间的非线性作用关系。文章分析了桩土接触作用机理及适用的本构模型，使用国际通用有限元软件ANSYS建立了桩土非线性接触的三维实体模型，通过三种方法估算土体的变形模量，并与现场试桩实测结果进行比较和分析，为LNG大型储罐全罐应力分析中的桩基分析部分提供了一种新的思路和方法。

1 桩土作用理论

1.1 桩土作用研究方法

桩土交互作用是一种高度非线性接触问题，可以使用有限元法进行数值求解，分析时需要选择桩和土体适用的本构模型和屈服准则。对于桩体，一般软件中都有成熟的钢筋混凝土模型。土体的本构模型则比较复杂：一种为非线性弹性本构模型，如邓肯-张模型[1]、K-G模型；一种为弹塑性本构模型，主要为Dracker-Prager模型[2]（以下简称为DP模型）。LNG储罐桩基所处的土体符合弹塑性本构特征，使用DP模型进行模拟。

1.2 DP模型

ANSYS软件中目前能用于近似描述岩土材料的模型只有DP模型。DP模型是理想弹塑性模型，是对Mohr-Coulomb准则的修正。它主要包含三个参数：内聚力c（Pa或kPa）、摩擦角准（°）和剪胀角Φf（°），其屈服准则为：

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

式中：β为与摩擦角准有关的无量纲材料常量；σm为平均应力，即土力学中的静水压力，Pa或kPa；{s}为偏应力，Pa或kPa；T为矩阵转置符；[M]为Mises屈服准则中的系数项，无量纲；σy为与内聚力c和摩擦角准相关的材料屈服参数，Pa或kPa。

1.3 桩土荷载传递

桩顶在受到向下压的荷载时，桩身压缩产生向下的位移，桩土之间发生相对滑移。当施加荷载超过桩侧土体的极限侧摩阻力与桩端土体的极限端阻力之和时，桩体发生持续的大位移沉降。

桩土接触是一个典型的刚性面-柔性面非线性接触问题，可以引入Goodman无厚度四节点接触单元[3]。在ANSYS分析中，选用TARGE170描述桩体的刚性目标面，CONTAC173描述桩侧土的柔性接触面，组成接触单元，接触形式为“No separation”（允许滑移和禁止分离），接触摩擦形式为库伦摩擦，摩擦系数的取值：黏土为0.3、粉质黏土为0.5、粉砂为0.8[4]。

1.4 土体的变形模量

土体的变形模量E0和压缩模量Es是判断土体压缩性和计算地基压缩变形的重要指标。在计算桩基沉降时，应使用土体的变形模量，而实际地基勘察一般不直接提供土体的变形模量。

本文通过三种方法估算土体的变形模量，并比较了三种方法的沉降拟合结果：

（1） 基于多孔介质理论、由Es推导得到E0。当假定土体为理想弹性体时，根据广义虎克定律，可推导得到E0与Es的关系为[5]：

式中：μ为土体的泊松比，无量纲；Es、E0的单位为MPa。

然而，实际的土体是应力-应变关系极其复杂的弹塑性材料，其变形模量通常不满足该式。大量的统计资料表明，E0有时可能是Es的3～5倍，甚至更大。

本文根据文献[6]中基于多孔介质理论建立的地基土变形模量估算方法，计算得到工区各层土体的变形模量。

式中：Ψ为土体扰动补偿系数，取值1.49；e为地基土孔隙比。

（2）通过标贯击数N值估算E0。标准贯入试验（SPT） 是国内外应用广泛的一种现场原位测试技术，工程中可以通过相应的经验公式，拟合得到变形模量与N值的近似函数关系。通过文献调研，分别使用式（4）中的经验公式[7-11]，计算得到相应的E0，代入ANSYS中进行计算和比较。

（3）通过参数反分析法拟合E0。方法（1）、（2）只是对E0的近似估算，估算结果的准确性依赖于土体性质的相似性和所采用经验公式的适用性。由于DP模型的参数较少并且通常都能够获得，因此在缺乏地基土勘察数据时，可采用参数反分析的方法，通过拟合桩顶的最终沉降位移，来估算土体的实际E0。

2 试桩概况

试桩工区位于天津市滨海新区，桩型为钻孔灌注桩，桩径1.2 m、桩长56 m，设计单桩竖向抗压极限承载力为9 600 kN。

2.1 单桩静载试验

储罐建造的前期，需要在工程桩施工前进行单桩静载试验。根据JGJ 106-2014《建筑桩基检测技术规范》进行单桩竖向抗压静载试验，采用如图1所示的压重平台反力装置，在桩顶竖向逐级施加等量荷载（逐级荷载为设计单桩竖向抗压极限承载力9600 kN的1/10）并维持，最终加载至9 600 kN。

图1 单桩竖向抗压静载试验压重平台反力装置

2.2 桩基土体参数

工区桩基土体包含多套土层，桩身上部土层主要为冲填土、黏土和粉质黏土，下部土层主要为粉土、粉质黏土和粉砂。桩端所在持力层为⑩4粉砂层，该土层呈密实状态，工程地质性质好。

试桩工区土体各土层的参数见表1。

3 模型与计算

3.1 ANSYS模型

在ANSYS中建立如图2所示的三维桩土实体模型（1/4模型），土体为长方体，宽为桩径的15倍，厚约64 m，包含12层土层。

表1 试桩工区土体参数

图2 桩土三维实体模型

采用映射网格划分，见图3。模型的边界条件设置为：土体两侧面分别固定X方向和Y方向上的位移UX、UY，桩体和土体的对称面采用对称边界条件，土体底面固定竖直方向上的位移UZ，桩端面与土层接触面的自由度耦合，桩顶面施加竖直向下的荷载。

图3 桩土网格模型

3.2 初始地应力场

在桩浇筑之前，土体早已在自身重力的作用下产生了较大的竖向位移和初始应力，后者直接影响桩土之间的侧摩阻力。因此，在分析桩沉降时，需要在去除土体初始位移的同时保留初始地应力。

ANSYS中可以使用inistate命令模拟初始地应力的自平衡，平衡前后的初始位移场见图4。在自重与初始地应力平衡前，土体在自身重力作用下发生了较大的竖向位移，顶部土层位移最大，达209 mm，远大于试桩的最终沉降量。平衡后，土体的最大位移为7.86×10-3mm，原始位移场基本清零。图5说明，在消除土体原始位移场的同时，也很好地保留了原始地应力场。

图4 初始位移场/m

3.3 土体宽度的选择

在建立三维桩土模型时，需要确定合适的土体宽度。过小的土体宽度会产生较大的边界效应，影响计算结果的准确性，过大的土体宽度则会加大求解的计算量。分别采用5、10、15、20、30倍桩径的土体宽度，计算桩端的最终沉降量。如图6所示，随着土体宽度的增大，边界效应减小，当土体宽度为桩径的15倍及以上时，最终沉降量收敛于同一值，说明此时边界效应已经小到可以忽略不计。因此，在后文的算例中，均采用15倍桩径的土体模型。

图5 初始地应力场/Pa

图6 桩端最终沉降量-土体宽度与桩径之比关系曲线

3.4 三种变形模量估算方法的比较

在1.4节中介绍了三种估算土体变形模量的方法，在保持其他参数不变的条件下，比较和分析了这三种方法的计算结果。

（1）方法一，基于多孔介质理论，由ES推导得到E0。前期的勘察工作已经获得了各土层的天然孔隙比、泊松比、压缩模量等参数，根据式（3）估算得到了各土层的变形模量，代入模型计算得到图7所示的荷载-沉降曲线。

图7 由方法一得到的沉降-荷载曲线

（2）方法二，通过标贯击数N值估算E0。根据式（4）中的各经验分公式，分别估算得到各土层的变形模量（见图8），并进一步计算得到沉降-荷载曲线（见图9）。

图8 由方法二各经验分公式计算得到的土变形模量

图9 由方法2得到的沉降-荷载曲线

可以看出，当施加荷载较低时，荷载主要通过桩身上部的侧摩阻力传递到上部桩周的土体，而上部土层的变形模量差别不大，因此各经验公式计算得到的沉降量差别不大；当施加荷载＞4 000 kN时，下部土体侧摩阻力逐渐发挥作用，荷载也传到桩端土体，由于下部土层尤其是持力层土层的变形模量差别较大，因此计算的沉降量也产生较大偏差。估算的变形模量越大，计算的沉降量越接近于真实值，因此在后续的方法比较中，选择式（4）中的分公式（4-5）的计算结果作为方法二的典型算例。

（3）方法三，通过参数反分析法拟合E0。工程上一般难以获得土体的真实变形模量，因此当仅有室内测得的压缩模量时，可以通过倍数放大压缩模量的方法拟合沉降-荷载曲线，从而估算变形模量。分别假定1、1.5、2、3倍压缩模量作为变形模量，计算得到如图10所示的荷载-沉降曲线。当取1倍压缩模量作为变形模型模量时，计算沉降偏大，说明工区土体的实际变形模量应大于压缩模量；当取2倍压缩模量时，计算结果拟合最好，认为工区土体的实际变形模量约为压缩模量的2倍。

图10 由方法三得到的荷载-沉降曲线

（4）三种估算方法的比较。如图11所示，选择由这三种方法拟合得最好的算例进行比较：方法一（基于多孔介质理论、由Es推导E0）在荷载与沉降的早期与实测数据的吻合度较高，但随着荷载和沉降量的增大，开始发生偏移；方法二（通过标贯击数N值估算E0）、方法三（通过参数反分析法拟合E0）都是通过放大压缩模量来估算变形模量的，拟合结果相近，在荷载与沉降的早期不如方法一，但后期的吻合度更高。

3.5 结果分析

从三种估算方法的比较中可以看出，方法二和方法三的拟合度更高，拟合结果差别不大，而方法二依赖于所选取经验公式的适用性，因此选取方法三的拟合结果做进一步的分析。图12为土体在最大荷载下（9 600 kN） 的最终沉降云图，在桩顶的竖向荷载下，桩周土体产生以桩为中心的沉降盆：横向上，离桩越近，土体沉降越大，土体侧边界的沉降为0；纵向上，第2、3层土体的变形模量最小，压缩性高，因而沉降量较大。桩端处的土体提供端阻力，受到向下的压应力，因而相比附近土体也有明显的下沉。

图11 三种估算方法计算结果对比

图12 土体最终沉降云图/m

图13为桩体的最终位移云图，桩身各点的位移随着深度的增加逐渐减小，桩顶处沉降量最大，为14.85 mm，桩端处沉降最小，为5.17 mm。可以看出，5.17 mm为桩端处土体受压产生的变形，而桩顶的总沉降为15.87 mm，其差值为桩土接触面的相对滑移和桩体自身的变形。这说明对于这种超长桩，其在上部载荷作用下的变形沉降，不仅来自桩端持力土层的变形，还包括桩土的相对滑移和桩体自身的变形，其分担比例取决于桩侧土和桩端土的变形模量和压缩性。

图14为最大荷载时桩顶到桩端的侧摩阻力分布曲线，可以看出，侧摩阻力总体上随着深度（桩侧接触压力，见图15）的增加而增大；桩顶处土体受到向下的剪应力，发生明显的相对滑动（见图16）；桩身下部侧摩阻力随着深度的增加而增大，桩端处达到极限值。计算得到桩侧土体的总侧摩阻力为6 785 kN，占总荷载的69.7%；桩端土体的端阻力为2 914 kN，占总荷载的30.3%。

图13 桩最终位移云图/m

图14 桩侧摩阻力-距桩顶距离曲线（9 600 kN）

图15 桩侧接触压力-距桩顶距离曲线（9 600 kN）

4 结论

（1）使用ANSYS建立了桩土非线性接触的三维实体模型，考虑了初始地应力和土层宽度的边界效应，模拟结果与现场实测的荷载-沉降数据拟合度较好，说明了该模型的准确性。

图16 桩土相对滑移距离-距桩顶距离曲线（9 600 kN）

（2）使用三种方法估算土体的变形模量，并进行比较：方法一（基于多孔介质理论，由Es推导E0）的偏差相对最大；方法二（通过标贯击数N值估算E0）依赖于所选取经验公式的适用性；方法3（通过参数反分析法拟合E0）在缺乏原位测试数据时，能得到较好的拟合结果。

（3）侧摩阻力总体随着深度的增加而增大，桩端处附近达到极限值。桩侧土体的总侧摩阻力占总荷载的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]许宏发，吴华杰.桩土接触面单元参数分析[J].探矿工程，2002（5）：10-12.

[5]陈希哲.土力学地基基础[M].第四版.北京：清华大学出版社，2004：210-215.

[6]梁发云.基于多孔介质理论的地基土变形模量估算方法[J].岩土力学，2004，25（7）：1 147-1 150.

[7]郭淋，王春艳，张飞，等.标贯试验N值与土体物理力学参数的相关性分析[J].安全与环境工程，2012，19（4）：148-152.

[8]徐光大，徐光黎，李俊杰.日本标准贯入试验方法及其N值在岩土工程中的应用[J].安全与环境工程，2011，18（4）：33-38.

[9]李修建，李松.标准贯入与载荷试验数据的对比分析[J].山西建筑，2011，37（10）：58-59.

[10]祝龙根，刘利民，耿乃兴.地基基础测试新技术[M].北京：机械工业出版社，2002：175-176.

[11]杨文卫，岳中琦.世界各地标准贯入试验比较和共同问题[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年毕业于北京大学力学专业，硕士，现主要从事LNG储罐结构分析方向的研究工作。

2017-01-03

中国海洋石油总公司科技课题“大型LNG储罐结构设计软件开发”（CNOOC-KJ 135 KJXM QD 2016-002）

Email：zhaozq4@cnooc.com.cn