考虑土体非线性时被动桩内力和变形特性计算

张超凡

(中铁第一勘察设计院集团有限公司， 陕西 西安 710043)

土位移法采用桩土相对位移与土体抗力指标综合反映土体侧移作用于桩身的被动荷载效应，能够从机理上较好地体现桩土相互作用对桩基性状的影响，国内外学者对此进行了一系列的研究[1～8]。张爱军等[2]采用两阶段法，基于Winkler地基模型研究不同土体侧移模式下抗滑单桩的变形和受力性状，得到了桩身响应的矩阵表达式；符金库等[3]基于竖向荷载作用下的Boussinesq水平位移解，通过对地面堆载作用效应的二重积分，得到了堆载作用下土体水平位移场的计算表达式；杨敏等[4]将p-y曲线法与Poulos弹性理论法相结合，提出了一种承受土体侧向位移的单桩计算方法。该方法依据p-y曲线确定土体的弹性模量，同时采用弹性理论法来考虑桩土之间的相互作用，计算过程更符合实际情况；周东等[5]基于室内模型试验设计了侧向位移作用下被动桩桩周土体水平位移测量装置，获得了不同条件下桩周土体自由位移场；宋修广等[6]基于双参数地基模型建立被动桩挠曲微分方程，以此研究堆载对基桩的影响；刘敦平等[7]采用ANSYS有限元软件对土体侧向位移作用下的被动桩进行分析，采用DP(Drucker-Prager)材料考虑土体的弹塑性性质，并考虑大变形的影响，研究了群桩效应和堆载作用下桩基内力和变形；梁发云等[8]采用ABAQUS软件，建立了受土体位移影响的单桩三维数值模型，讨论了不同土体侧移模式对桩身内力和变形的影响。

目前，多数学者将堆载作用下的土体自由位移场假定为土体深度的函数，或采用经验公式对其进行计算，未与堆载大小、桩周土体性质建立关系，不够准确；另一方面，计算桩周土体抗力时多采用弹性地基反力法，不能反映桩周土体的非线性。本文将位移法和p-y曲线法相结合，基于位移法计算邻近堆载对桩基产生的影响，该方法将土体水平位移与桩周软土层厚度、土体力学指标及固结度、堆载大小等因素联系起来，并采用p-y曲线法描述桩侧土抗力与桩身挠度之间的关系。同时考虑主动荷载(桩顶竖向荷载和水平荷载)及被动荷载(邻近堆载)的影响，进行被动单桩受力特性计算。通过有限差分法对桩身挠曲微分方程进行求解，并据此编制了Fortran程序，通过实例对比验证了方法的正确性与适用性。

1 模型建立

1.1 计算模型

计算模型如图1(图中：a1，a2，a3，b1，b2，b3为堆载尺寸；h为堆载高度；e为堆载边缘距桩身的距离；N0，H0，M0为桩顶主动荷载，分别为竖向荷载、水平荷载、弯矩)所示，桩顶与地面齐平，其上作用有竖向荷载、水平荷载和弯矩，桩周有邻近堆载。为便于采用有限差分法进行计算，对桩身进行分段，并在桩顶、桩端分别增加两个虚拟单元。力和位移的符号约定如下：水平力和位移沿y轴方向为正，桩身弯矩以靠近堆载侧桩身受拉为正。

图1 堆载作用下被动单桩计算模型

为了研究主动荷载(竖向荷载、水平荷载、弯矩)和邻近堆载同时作用下单桩受力特性，首先应建立其挠曲微分方程。为便于计算，本文假定桩身为弹性体，同时本文研究的重点是桩的横向内力和变形特性，故假定轴力为常数。取桩身单元如图2所示，单元上端作用有弯矩M、剪力Q和轴力N，单元下端作用有弯矩M+dM、剪力Q+dQ和轴力N，桩侧有土抗力P=b0p(z,y)和堆载产生的桩身附加荷载b0czU(E,z)。

根据桩身微单元平衡条件，可得桩身挠曲微分方程为

(1)

式中：E为桩身弹性模量；I为桩身惯性矩；N为桩身轴力；y为桩身水平位移；z为计算深度；b0为桩的计算宽度，柱桩及桩身直径0.8 m以下的灌注桩b0=d+1(m)，其余类型的桩参照相关资料，其中d为桩径[9]；cz为地基系数，表示单位面积土在弹性限度内产生单位变形时所需加的力，cz=p/y；U(E，z)为土体自由场水平位移。

1.2 边界条件

根据桩身实际受力特点，对边界条件作出如下假定：

(1)桩顶边界条件

本文仅考虑一种情况，即桩顶作用有竖向荷载N0、水平力H0、弯矩M0。

(2)桩底边界条件

桩底边界条件有两种情况：桩底自由或桩底固定。

1.3 堆载作用下土体位移场计算

采用两阶段法计算被动桩桩身内力和变形时，土体自由场水平位移是否准确直接关系到分析结果的合理性[10]。本文在计算时，鉴于现场或模型试验耗时费力，实施较为困难，因此通过理论公式计算堆载引起的桩周土体水平位移。Bourges等[11]针对路堤堆载作用下土体水平位移进行了长期现场监测，根据监测结果提出了路堤堆载作用下土体水平位移经验公式，其所采用的计算图图示如图3所示。

图3 路堤堆载作用下土体水平位移监测示意

Bourges等[11]定义了桩周土体侧移的三种分布曲线，均为土体深度的函数，如图4(图中：ρ为土体水平位移；ρmax为土体最大水平位移；y=ρ÷ρmax；d为软土层厚度；z为桩长范围内土体深度与软土层厚度的比值)所示。

图4 土体水平位移随深度的分布形式

三条曲线方程如下：

曲线①：整体平均曲线

y=1.83z3-4.69z2+2.13z+0.73

(2)

曲线②：可压缩层位于地表以下几米处的情况

y=3.42z3-6.37z2+2.14z+0.81

(3)

曲线③：正常固结的沉积土

y=-2z3+1.5z+0.5

(4)

土体最大水平位移ρmax由瞬时分量ρimax和固结分量ρcmax组成。ρimax与X/L(见图3)的值和边坡安全系数F的值有关。安全系数F按式(5)计算。

F=(π+2)cuav/(γH)

(5)

式中：cuav为黏土层平均不排水抗剪强度；γ为路堤填土重度；H为路堤高度。

Bourges等[11]建立了最大水平位移的无量纲量λ(λ=ρimax/d)与边坡稳定系数F的关系曲线[11]。任意时刻的最大水平位移固结分量ρcmax与路堤中心的固结沉降量Sc有关，按式(6)计算。

ρcmax,t=0.16Sc=0.16(St-Si)

(6)

式中：St为t时刻的总沉降量；Si为瞬时沉降量。

St可根据土体的固结度U和最终沉降量Sf进行计算，故任意时刻土体最大水平位移ρmax可用下式表示：

(7)

在无准确实测数据或其他确定方法时，可采用本节所述方法初步估算堆载引起的桩周土体水平位移。

1.4 p-y曲线选取

本文选取国内外五种典型的p-y曲线模型[12～16]。从模型假定、适用条件、参数计算等角度对其进行对比，最终采用章连洋等[16]提出的模型计算土体抗力与桩身水平位移的关系。因篇幅有限，对比过程不再详述。

章连洋等在Matlock的基础上进行改进，引入有效长度概念，提出了适用于我国长江中下游地区软黏土的p-y曲线计算方法。计算步骤如下：

(1)计算不同深度处土体不排水抗剪强度cu以及最大主应力差一半时的应变ε50[16]。

(2)计算静载作用下桩的有效长度Ls和影响深度zrs，如式(8)(9)所示。

(8)

zrs=Ls/4

(9)

式中：Ec为桩长范围内土体弹性模量的平均值。

(3)计算土体极限抗力Pu，如式(10)(11)所示。

Pu=FsNpcu

(10)

Np=min{2.5+6.5(z/zrs)，9}

(11)

式中：Fs为与土性和荷载密切相关的折减系数[16]。

(4)计算桩侧土反力达到极限值的50%时对应的桩身水平位移y50，其值为：

y50=4.5ε50d3/4

(12)

(5)建立p-y曲线模型：

曲线段：

p/pu=0.5(y/y50)1/3

(13)

直线段：

(14)

p/pu=1，z>zrs

(15)

2 模型求解

采用有限差分法求解桩身内力和变形，将桩身沿长度方向分为n个单元，分段长度为λ，为便于计算，避免出现零下标和负下标，将桩顶和桩底各增加两个虚拟节点，节点编号从1开始，共n+5个节点，如图 1所示。

对于分段点i，由一维中心差分公式可得i节点的差分方程为：

(16)

式中：yi为i节点处对应的桩身水平位移；λ为桩身分段长度(见图 1)；czi为i节点深度处单位面积土在弹性限度内产生单位变形时所需加的力，czi=pi/yi；Ui为深度z处土体自由场水平位移。

yi+2-ryi+1+(2r-2+Ai)yi-ryi-1+
yi-2-AiUi=0

(17)

根据微分关系，可以求得各截面的转角θi、弯矩Mi、剪力Qi为：

(18)

(19)

(20)

由图 1可知，桩身范围内共有n+1个节点，对每个节点i(i=3，4，…，n+2，n+3)均可写出如式(16)所示的平衡方程，共n+1个方程，但有n+5个未知的桩身水平位移yi，还需补充4个方程才能求解，这4个方程根据1.2节边界条件确定。至此，n+5个方程对应有n+5个未知的yi，方程数与未知数个数相等，方程可解。

采用追赶法求解方程组，为便于表达，引入Bi，Ci，Di作为计算i节点处桩身变形特性的三个无量纲中间变量。此处以桩底自由为例，计算结果如下：

结合桩底边界条件可得：

yn+5=(r-2)yn+4-(r-2)yn+2+yn+1

(21)

yn+4=2yn+3-yn+2

(22)

yn+3=Bn+3yn+2+Cn+3yn+1+Dn+3

(23)

对于中间节点i=n+2至i=3，有

yi=Biyi-1+Ciyi-2+Di

(24)

(25)

(26)

式中：zj1=B3(B4-2)+C4+1；zj2=C3(B4-2)；zj3=D3(B4-2)+D4；zj4=(B5-k+2)(B3B4+C4)+C5B3+k-2；zj5=C3B4(B5-k+2)+C3C5-1；zj6=D5+(B5-k+2)(D4+D3B4)+D3C5。

3 程序编制及验证

3.1 DSN程序功能与流程

根据上述求解思路，采用Fortran语言编制了DSN程序，该程序可用于求解考虑土体非线性的被动单桩桩身内力和变形。求解思路如下：首先输入关于土层、桩及堆载的计算参数，然后根据实际情况选取边界条件，采用Bourges等[11]提出的桩周土体自由场水平位移，将所得的土体水平位移施加于桩身，通过有限差分法计算出桩身内力和变形。程序流程如图5所示。

图5 考虑土体非线性时被动单桩受力特性计算流程

3.2 注意事项

(1)关于地基系数的说明

本文在推导理论公式时，地基系数czi=pi/yi，对于桩身位移非常小甚至接近于0时，按此公式计算所得czi为无穷大，程序不能继续运行。故程序对于桩身水平位移yi<0.1 mm的yi值，取yi=0.1 mm，此精度可满足工程要求。

(2)关于程序迭代终止条件的说明

程序迭代终止判别条件采用相邻两次计算所得的桩身水平位移差值ε进行衡量，两次计算所得的ε值小于0.01 mm时终止迭代。

3.3 程序验证

为验证DSN程序的正确性，本文选取上海宝山钢铁厂群桩试验案例[17]以及浙江某地区不平衡堆载试验[18]分别进行计算，并与现场实测数据进行对比。

3.3.1 工程案例一

上海宝山钢铁厂群桩堆载试验位于长江入海口附近，平面布置如图 6所示。试验桩为空心钢管桩，桩长60 m，桩径609 mm，壁厚11 mm，桩身弹性模量E=2.1×106kPa。承台尺寸为5.4 m×5.4 m。试验场地土层自上而下分别为：0～2.5 m为黄褐色亚黏土；2.5～7.5 m为淤泥质亚黏土；7.5～22.5 m为淤泥质黏土；22.5～60 m为亚黏土；60 m以下为细砂。土的主要物理力学指标如表1所示。

图6 堆载试验平面布置/m

表1 土层参数

堆载采用分级加载模式，分别为60，90，120，150 kPa。第一级持续32 d，第二级持续98 d，第三级持续126 d，第四级持续161 d，共计417 d，测得各级荷载作用下基桩挠度以及桩周土体沉降和侧移值。

如图6所示，以第四级荷载施加完毕后B桩为研究对象，采用本文研制的DSN程序对上述工程进行计算。所取的基本计算参数如表2所示。

表2 计算参数取值

将堆载作用下土体自由场水平位移计算值与实测值进行对比，结果如图7所示。

图7 桩周土体水平位移计算值与实测值对比

参照1.3节土体水平位移经验公式，本工程在曲线②适用范围内。故本文将根据曲线②计算所得土体水平位移施加于桩身，采用DSN程序计算B桩水平位移，计算值与实测值对比如图8所示。

图8 B桩桩身水平位移计算值与实测值对比

由图8可知，桩身水平位移计算值与实测值相比，二者趋势相同，沿深度方向均先增大后减小。桩顶水平位移实测值为17.9 mm，计算值为13.9 mm，计算值比实测值小22%；桩身水平位移实测最大值为36.82 mm，发生在地下深度约9.0 m处，计算最大值为38.8 mm，发生在地下深度约8.5 m处。

3.2.2 工程案例二

不平衡堆载试验[18]简介如下：基桩共三根，桩长71.5 m，桩径1.5 m，桩间距5.6 m，持力层为黏性土圆砾。地表下存在流塑状高压缩性淤泥质黏土，厚度约24.5 m。桩周堆载高度3 m，形状为梯形，上表面尺寸14.2 m×5 m，下表面尺寸22.6 m×13.4 m，长边方向与基桩中轴线平行。

由于桩周存在深厚淤泥质黏土，满足1.3节中曲线②适用条件，故采用曲线②计算桩周土体水平位移，以文献[18]中SZ2为研究对象，采用DSN程序对桩身水平位移进行计算，计算值与实测值对比如图9所示。

图9 桩身水平位移计算值与实测值对比

由图9可知，桩身水平位移计算值与实测值沿深度方向均逐渐减小。桩顶水平位移实测值为105.4 mm，计算值为121.87 mm，计算值比实测值大15.6%，在误差允许范围以内；桩身水平位移实测值在地下23 m左右基本为0，桩身水平位移计算值在地下10 m左右基本为0。

对以上两个工程案例计算结果进行分析，误差产生的原因有以下几点：(1)现场实测数据不准确，现场试验时受人为及外部因素影响，测试结果会有一定偏差；(2)工程案例一中现场试验为2×2根的群桩试验，本文在计算时，通过将地面以下2 m深度范围内的土体变形模量适当提高考虑承台及群桩效应影响，存在一定的误差；(3)本文所述土体自由场水平位移公式有一定的局限性。

综上所述，在误差允许范围内，计算结果与现场实测值吻合良好。在无准确实测数据或其他确定方法时，可采用本文所述方法计算考虑土体非线性时被动桩内力和变形。

4 结 论

本文将位移法与p-y曲线法相结合，对考虑土体非线性的被动桩桩身内力和变形展开研究。主要结论如下：

(1)采用典型p-y曲线法描述桩侧土抗力与桩身水平位移的关系，并考虑桩顶竖向荷载N对桩身内力和变位的影响(P-Δ效应)，能综合反映桩周土体的非线性，桩的刚度和外荷载作用性质等特点，更符合实际情况；

(2)计算土体位移场时，将土体位移场分布情况同桩周土体性质、地面堆载等因素联系起来，可计算具体工况下被动桩桩身内力和变形，为工程实际提供参考；

(3)基于上述理论，采用Fortran语言编制了程序DSN，该程序可用于计算被动桩桩身内力和变形，将程序DSN应用于典型案例，通过计算结果与实测结果的对比，验证了计算方法和程序的正确性。