一种考虑土体非达西固结的单桩竖向承载时变效应计算方法与流程

本发明属于土木工程领域，涉及一种考虑土体非达西固结的单桩竖向承载时变效应计算方法。

背景技术：

桩基础是深入土层的柱型构件，主要作用是将上部结构的荷载通过桩身穿越较软土层或水传递至深层相对坚硬的持力层，降低上部结构沉降并保证安全，其承载力的发挥一般以竖向为主，在某些工程中也兼具水平承载的作用。桩基础的竖向承载力来自桩侧土体提供向上的摩阻力和桩端土体提供的端阻力，当地面受到大面积堆载或浸水等环境因素作用后，桩周土体将发生明显的固结沉降，从而通过桩土界面摩擦作用引起桩基压缩量增加，桩身轴力变大甚至引起桩端沉降，桩侧原有部分正摩阻力变为负摩阻力，降低桩基础原有的竖向承载能力，并且这一复杂的桩土相互作用是一个长期过程，桩基础的负摩阻力问题在工程前期设计中应当充分考虑，因此准确计算固结土体中桩基承载时间效应显得尤为重要。

目前，诸多学者已经开展了关于固结土体中桩基础的离心模型试验、室内试验以及相关的有限元仿真计算并得出一些有价值的结论。但是这些研究成果只能反映特定情况下的桩身响应，并且没有考虑随着固结时间推移土体有效应力的动态变化特性。另外，当孔隙水的水力梯度较低时，流速与水力梯度呈非线性关系；当实际水力梯度超过流动孔隙水最大结合能所对应的水力梯度时，流速与水力梯度呈线性关系，而现有研究主要是基于达西定律的经典固结理论进行，非达西流动定律尚未得到考虑，这与实际情况存在差异。

技术实现要素：

本发明的目的在于克服现有技术中的不足，提供一种考虑土体非达西固结的单桩竖向承载时变效应计算方法，解决现有技术中大面积堆载作用下土体非达西流动固结的沉降计算问题以及桩基竖向承载力动态变化的计算问题。

为了实现上述的技术特征，本发明的目的是这样实现的：一种考虑土体非达西固结的单桩竖向承载时变效应计算方法，包括以下步骤：

步骤1：根据非达西流动定律建立一维固结平衡微分方程，通过差分格式离散方程，获得固结偏微分方程的差分形式；

步骤2：引入不同排水条件下初始边界条件和成层土交界面处渗流连续条件，求解方程获得土体的超孔隙水压力、渗透系数、固结系数、孔隙比解答；

步骤3：由孔隙比的变化计算出固结过程中每一时刻桩侧土体沉降，并通过超孔隙水压力的变化计算土体有效应力以及土体平均固结度；

步骤4：引入桩土界面荷载传递曲线，结合桩身轴力位移平衡方程和桩侧摩阻力平衡方程，加入桩端边界条件，通过假定桩顶变形迭代求解获得桩周土体固结过程中桩身变形、轴力、桩侧负摩阻力的分布情况。

所述步骤1具体步骤如下：

步骤1.1：在较低水力梯度情况下，土体孔隙水流动遵守非达西流动准则，流速与水力梯度呈非线性关系，当水力梯度超过一定界限值呈现线性关系，孔隙水的流速与水力梯度的关系具体如下式所述：

式中：为水力梯度；i0＝ic(n-1)/n，为临界水力梯度；ic为启动压力梯度；n和α为非达西定律的拟合参数，其中，α＝ic^(1-n)/n；v为流体通过土体介质的流速；当上式中α和n等于1时，非达西流动法则就退化为达西定律；

步骤1.2：当土层受到地表堆载q作用时，土体微单元中水的流入和流出沿z轴单方向流动，x和y方向单位渗流量为零，土体单元水的流出量与流入量差值等于土体单元体积变化率，存在等式：

孔隙水流量差：

土体单元体积变化率为：

式中：γw为水的重度；u为超孔隙水压力；k为土体的渗透系数；av为土体压缩系数；e为土体的孔隙比。

结合等式可以获得固结偏微分方程：

式中：固结系数

所述步骤2具体步骤如下：

步骤2.1：将土层沿深度划分为m等分，节点编号依次为0,1…m，共m+1个节点，土层微段厚度为h＝h/m，时间划分为n等分，时间编号依次为0,1…n，共n+1个节点，时间步长为δt＝t/n，采用差分格式对固结偏微分方程进行离散，近似代替各个变量的偏微分形式，把连续问题变为多个方程的离散问题，得到固结控制偏微分方程的差分格式如下：

1)当0≤ih＜ic时，孔隙水的流速与水力梯度呈现非线性关系，固结控制偏微分方程的差分格式为：

2)当ih≥ic时，孔隙水的流速与水力梯度呈线性关系，固结控制偏微分方程的差分格式为：

式中：ui,j指固结土层深度方向第i点在时间j时的超孔隙水压力；ki,j指固结土层深度方向第i点在时间j时的渗透系数；cvi,j指固结土层深度方向第i点在时间j时的固结系数；

步骤2.2：在求解偏微分方程时初始边界条件分为单面排水和双面排水两种情况，当上表面排水底面不排水时边界条件为：

当上表面排水底面排水时边界条件为：

步骤2.3：成层土计算时两层土之间具有边界连续条件和渗流连续条件；

在土层交界面处超孔压相等：u1[m][j]＝u2[0][j]；

考虑渗流连续条件时相邻土层交界面流速相等：

式中：k1指第一层土体的渗透系数；k2指相邻第二层土体的渗透系数；

所述步骤3具体步骤如下：

步骤3.1：由孔隙比的变化计算出固结过程中每一时刻桩侧土体沉降：

步骤3.2：在固结过程中随着孔隙水的排出，土体内部孔隙逐渐减小，有效重度随之发生改变，根据有效应力原理来计算单层土某一深度实时有效应力为：

对于成层土，第i层土的有效应力可以按下式计算：

式中：γ′为某一时刻土体的有效重度；σ′(z,t)为固结t时刻深度z处的有效应力；u(z,t)为固结t时刻深度z处的超孔隙水压力；gs为土粒比重；e(z,t)为固结t时刻深度z处的孔隙比；

步骤3.3：某一特定时间t＝j时刻土体的平均固结度u可以按照下式计算：

式中：ui,j指土层i处在j时间的超孔隙水压力，ui,0是在土层i处初始超孔隙水压力。

对于任意给定固结时间t，能够通过下述迭代求解计算获得所有桩节点的摩阻力、沉降、轴力，步骤四的具体步骤如下：

步骤4.1：根据土层的分层情况将嵌入土体的桩划分为k段，每个桩段再划分为与周围土层相同的m微段，同时桩顶施加上部结构传递的竖向荷载；

步骤4.2：假定一个微小量作为第一段桩桩头变形初始值，其节点编号i＝0，此处对应的剪应力按照如下公式计算：

式中：为第一段桩第i节点在第一轮迭代获得的桩侧剪应力；为第一段桩第i节点在第一轮迭代的桩身变形；st1(i)为第一段桩第i节点相邻土体沉降；ks1(i)为第一段桩第i节点桩侧剪切刚度，按照ks1(i)＝gsoil1(i)/(r0ln(rm/r0))计算，gsoil(i)为桩周土的剪切模量；r0为桩的半径；rm为最大影响半径，对于摩擦桩rm＝2.5ρ(1-υ)l，对于端承桩rm＝{0.25+[2.5ρ(1-υ)-0.25]ξ}l；式中：υ为按深度加权平均的桩侧土泊松比；l为桩嵌入土体的深度；ρ为非均匀系数，体现桩身中部与桩尖处桩侧土剪切模量的比值，按照公式计算；ξ为桩尖处桩侧土体剪切模量与桩端以下土体剪切模量的比值，按ξ＝gl/gb计算；为桩嵌入深度范围内剪切模量平均值，按计算；gl为z＝l处桩周土的剪切模量；gb为桩端土的初始切线剪切模量,按照公式计算；es为桩端以下土体压缩模量；vs为桩端以下土体泊松比；

桩土界面剪应力需小于等于桩土界面抗剪强度τf(i)，τf1(i)＝k01(i)tanδ1·σ′1(i)；式中：k01为第一层土的静止土压力系数；σ′1(i)为第一层土i节点的竖向有效应力；δ1为第一层土桩土界面摩擦角，按照δ1＝κ·φ′1计算；κ为强度折减因子；φ′1为第一层土的内摩擦角；

步骤4.3：相邻节点桩身轴力p1⁽¹⁾(1)按照公式算得，编号为1的桩节点沉降按照公式计算，该节点摩阻力按公式计算；

步骤4.4：然后将第一桩段由桩顶向下按照公式依次计算第2节点至第n节点的桩身轴力、桩身变形以及桩侧剪应力；

步骤4.5：随后计算第二桩段的内力和变形，第二桩段第0节点沉降和轴力与第一桩段最后一节点相等，将第二个桩段由上向下按照公式

依次计算其余节点的桩身轴力、桩身变形以及桩侧剪应力；运用同样的方法将第三至第k所有桩段计算完毕，最后可获得第一轮迭代后的桩端轴力和桩端沉降

步骤4.6：再通过桩端边界条件计算新的桩端沉降根据上一步计算获得的桩端轴力按照公式计算；式中：gb为桩端土剪切模量；vs为桩端土泊松比；wcb为桩端相邻节点桩周土体沉降；αb为桩端承载力水平系数，按αb＝rfpb/pfb计算；rf是用来描述应力应变双曲线模型关系的破坏比，按rf＝τf/τult计算；τult为桩土界面极限剪切应力；pfb为桩端极限承载力，计算；cs为桩尖以下土体粘聚力；ab为桩端横截面面积；为桩端土的内摩擦角；为桩端平面处桩侧平均竖向压应力；为桩端以上土体的有效重度，对于多层土而言按照深度加权平均法计算，

步骤4.7：若基于假定桩顶变形计算的桩端沉降和按照桩端边界条件计算的桩端沉降两者之间误差小于容许误差，即说明在第一步中假定的桩顶变形是满足要求的；否则，需假定新的桩顶变形并重复多次迭代直到第i轮迭代计算的沉降误差即停止迭代，输出土体固结过程中土体沉降、桩身各节点侧摩阻力、桩身变形、轴力等数据。

迭代收敛标准建议值取容许误差ζ＝1×10^-6m。

为保证数值解的稳定性和精度要求，在采用差分方法求解固结偏微分方程时，常量时间步δt和土层厚度h小于1/4。

在进行桩土微单元划分时，应遵循同一土层对应的桩土微单元厚度相等，不同土层对应的桩土微单元厚度接近的原则。

本发明有如下有益效果：

1、本发明从软土固结定律的角度出发，考虑固结过程孔隙比和渗透系数等重要参数的动态变化，计算桩侧土体固结沉降，并结合桩土界面荷载传递模型，计入了土体有效应力增加对桩土界面强度的影响，最后计算桩基在堆载作用下的长期响应。

2、本发明所阐述的计算方法既可以考虑非达西流动也可以考虑达西流动，同时固结方程的解答非常简洁，易于编程计算，可计算地表大面积堆载与桩顶竖向荷载的组合效应，能够更准确地解决堆载作用下桩周土体在固结过程中桩基长期承载特性问题。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1为基于非达西定律固结土体中单桩竖向时变效应计算流程图。

图2为流速与水力梯度关系曲线。

图3为土体剪切模量比，其中图3(a)为摩擦桩，图3(b)为端承桩。

图4为桩土单元离散示意图。

图5为桩土荷载传递模型示意图，其中图5(a)为桩微段，图5(b)为桩单元，图5(c)为轴力分布。

图6为桩单元受力示意图，其中图6(a)为桩土受荷示意图，图6(b)为侧摩阻力，图6(c)力学简化模型。

图7为实测沉降与计算沉降对比曲线。

图8为摩擦桩桩身下拉荷载实测结果与本发明计算结果对比曲线。

图9为端承桩桩身下拉荷载实测结果与本发明计算结果对比曲线。

具体实施方式

下面结合附图对本发明的实施方式做进一步的说明。

实施例1：

参见图1-9，本发明基于非达西流动定律推导出土体一维固结非线性方程，将固结过程中土体内部孔隙比、渗透系数、土体重度等重要参数的动态变化计入计算过程，通过有限差分法来求解控制微分方程，获得超孔隙水压力的数值解；基于孔隙比变化与沉降的关系建立桩侧土体固结沉降计算模型，再结合桩土荷载传递模型，计算固结土体中桩基竖向承载的长期响应。

如图1所示，为本发明的计算方法流程示意图，具体计算步骤如下：

(1)明确固结计算过程中使用的非达西流动定律计算公式，即土体内部孔隙水的流速与水力梯度按照计算，如图2所示为非达西定律下流速与水力梯度关系曲线；式中：为水力梯度；i0＝ic(n-1)/n，为临界水力梯度；ic为启动压力梯度；n和α为非达西定律的拟合参数，其中，α＝ic^(1-n)/n；v为流体通过土体介质的流速；当上式中α和n等于1时，非达西流动法则就退化为达西定律；

(2)建立一维固结平衡偏微分方程。假设地表堆载大小为q，土体微单元中水的流入和流出沿z轴单方向流动，x和y方向单位渗流量为零，土体单元水的流出量与流入量差值等于土体单元体积变化率，即：

等式左边孔隙水流量差为：

等式右边土体单元体积变化率为：

式中：γw为水的重度；u为超孔隙水压力；k为土体的渗透系数；av为土体压缩系数；e为土体的孔隙比；

结合等式左右两边建立固结偏微分方程：

式中：固结系数

(3)将方程离散获得固结偏微分方程的差分解答。将土层沿深度划分为m等分，节点编号依次为0,1…m，共m+1个节点，土层微段厚度为h＝h/m，时间划分为n等分，时间编号依次为0,1…n，共n+1个节点，时间步长为δt＝t/n，采用差分格式对固结偏微分方程进行离散，近似代替各个变量的偏微分形式，把连续问题变为多个方程的离散问题，为保证数值解的稳定性和精度要求，常量时间步δt和土层厚度h建议小于1/4，得到固结控制偏微分方程的差分解答如下：

当0≤ih＜ic时，孔隙水的流速与水力梯度呈现非线性关系，固结控制偏微分方程的差分解答为：

当ih≥ic时，孔隙水的流速与水力梯度呈线性关系，固结控制偏微分方程的差分解答为：

式中：ui,j指固结土层深度方向第i点在时间j时的超孔隙水压力；ki,j指固结土层深度方向第i点在时间j时的渗透系数；cvi,j指固结土层深度方向第i点在时间j时的固结系数；

(4)选择固结计算流动法则输入非达西固结指数n的数值，根据土体单面排水或双面排水边界条件计算所有时刻排水边界土体的固结参数；

当上表面排水底面不排水时边界条件为：

当上表面排水底面排水时边界条件为：

(5)计算出t＝0时刻土体所有节点的孔隙比、有效应力、渗透系数、固结系数，按照第(3)步由地表节点开始向下依次计算各土层节点的超孔隙水压力、孔隙比、有效应力、渗透系数以及固结系数；计算至土层交界面时通过边界连续条件和渗流连续条件公式进行计算，然后向下继续计算直到土体最后一个节点即完成超孔隙水压力计算过程，获得所有节点的超孔隙水压力、渗透系数、固结系数、孔隙比；

在土层交界面处边界连续条件为相邻土层交界面超孔隙水压力相等：

u1[m][j]＝u2[0][j]；

在土层交界面处渗流连续条件为相邻土层交界面流速相等：

式中：k1指第一层土体的渗透系数；k2指相邻第二层土体的渗透系数；

(6)由孔隙比的变化通过公式计算出固结过程中每一时刻桩侧土体沉降；固结过程中孔隙水排出，土体内部孔隙减小，有效重度随之发生改变，根据有效应力原理按照公式计算某一深度实时有效应力为；式中：γ′为某一时刻土体的有效重度；σ′(z,t)为固结t时刻深度z处的有效应力；u(z,t)为固结t时刻深度z处的超孔隙水压力；gs为土粒比重；e(z,t)为固结t时刻深度z处的孔隙比。

按照公式计算t＝j时刻土体的平均固结度；式中：ui,j指土层i处在j时间的超孔隙水压力，ui,0是在土层i处初始超孔隙水压力；

(7)根据土层的分层情况将嵌入土体的桩划分为k段，每个桩段再划分为与周围土层相同的m微段，同时桩顶施加上部结构传递的竖向荷载q，图4表示桩土单元离散示意图；

(8)假定一个微小量作为第一段桩桩头变形初始值，其节点编号i＝0，0号节点轴力p1⁽¹⁾(0)＝q，对应剪应力按照公式计算，图5表示桩土荷载传递模型；式中：为第一段桩第i节点在第一轮迭代获得的桩侧剪应力；为第一段桩第i节点在第一轮迭代的桩身变形；st1(i)为第一段桩第i节点相邻土体沉降；ks1(i)为第一段桩第i节点桩侧剪切刚度，按照ks1(i)＝gsoil1(i)/(r0ln(rm/r0))计算，gsoil(i)为桩周土的剪切模量；r0为桩的半径；rm为最大影响半径，对于摩擦桩rm＝2.5ρ(1-υ)l，对于端承桩rm＝{0.25+[2.5ρ(1-υ)-0.25]ξ}l；式中：υ为按深度加权平均的桩侧土泊松比；l为桩嵌入土体的深度；ρ为非均匀系数，体现桩身中部与桩尖处桩侧土剪切模量的比值，按照公式计算；ξ为桩尖处桩侧土体剪切模量与桩端以下土体剪切模量的比值，按ξ＝gl/gb计算，图3(a)表示摩擦桩的土体剪切模量比，图3(b)表示端承桩的土体剪切模量比；为桩嵌入深度范围内剪切模量平均值，按计算；gl为z＝l处桩周土的剪切模量；gb为桩端土的初始切线剪切模量,按照公式计算；es为桩端以下土体压缩模量；vs为桩端以下土体泊松比；

(9)按照公式τf1(i)＝k01(i)tanδ1·σ′1(i)计算桩土界面抗剪强度τf(i)，判断是否小于等于τf(i)，若大于则取τf1(i)＝τf(i)；式中：k01为第一层土的静止土压力系数；σ′1(i)为第一层土i节点的竖向有效应力；δ1为第一层土桩土界面摩擦角，按照δ1＝κ·φ′1计算；κ为强度折减因子；φ′1为第一层土的内摩擦角；

(10)第一桩段第1节点至最后一节点按照公式计算桩身轴力，按照公式计算桩节点沉降，按公式计算桩侧剪应力，直至第一桩段所有节点桩身轴力、桩身变形以及桩侧剪应力全部获得，图6表示桩单元受力示意图；

(11)随后计算第二桩段的内力和变形，第二桩段第0节点沉降和轴力与第一桩段最后一个节点相等，将第二个桩段由上向下按照公式依次计算其余节点的桩身轴力、桩身变形以及桩侧剪应力；运用同样的方法将第三至第k所有桩段计算完毕，最后获得第一轮迭代后的桩端轴力和桩端沉降

(12)再通过桩端边界条件计算新的桩端沉降根据上一步计算获得的桩端轴力按照公式计算；式中：gb为桩端土剪切模量；vs为桩端土泊松比；wcb为桩端相邻节点桩周土体沉降；αb为桩端承载力水平系数，按αb＝rfpb/pfb计算；rf是用来描述应力应变双曲线模型关系的破坏比，按rf＝τf/τult计算；τult为桩土界面极限剪切应力；pfb为桩端极限承载力，计算；cs为桩尖以下土体粘聚力；ab为桩端横截面面积；为桩端土的内摩擦角；为桩端平面处桩侧平均竖向压应力；为桩端以上土体的有效重度，对于多层土而言按照深度加权平均法计算，

(13)若基于假定桩顶变形计算的桩端沉降和按照桩端边界条件计算的桩端沉降两者之间误差小于容许误差，即说明在第一步中假定的桩顶变形是满足要求的；否则，需假定新的桩顶变形并重复多次迭代直到第i轮迭代计算的沉降误差即停止迭代，输出土体固结过程中土体沉降、桩身各节点侧摩阻力、桩身变形、轴力等数据；

上述实施例仅是本发明的优选实施方式，但本发明的具体保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易得到改进和变换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围以权利要求的保护范围为准。

实施例2：

为验证本发明方法的可行性，将本发明方法与其他学者所做堆载作用下单桩离心机试验进行对比，试验中离心加速度为60g，g为重力加速度，土体和桩的几何尺寸都放大为模型尺寸的60倍。试验1为端承桩，将一根铝制单桩打入18m的饱和软黏土中，桩尖嵌入底端砂层0.25d，d为桩径，试验2为摩擦桩，将一根17.7m的桩打入饱和软黏土，桩尖位于底部密砂层以上3m处，桩外径为1140mm，壁厚为90mm，内径为960mm，桩身弹性模量e为70gpa，在桩体打入完成后在地表回填3m厚砂子作为表面堆载并作排水面，回填砂重度为15kn/m³，黏土的初始孔隙比取1.6，剪切模量取2.4gpa，泊松比为0.35，持力层砂土的剪切模量取38gpa，泊松比为0.3；黏土初始渗透系数为0.13m/year，非达西固结指数n＝1.25，启动压力梯度ic＝10，试验监测了土体固结过程中单桩下拉荷载的变化。

如图7所示为土体固结过程中实测地表沉降随时间变化曲线与预测结果对比，本文理论解与试验发展趋势非常接近，在固结80个月后本文计算所得地表沉降为613mm，实测结果为615mm。如图8所示为摩擦桩实测下拉荷载与预测结果对比，摩擦桩中性面位于桩尖以上大约3.4m位置，在固结1年、2年及5年后桩身中性面处本发明计算的最大下拉荷载分别为641kn，870kn和1065kn，而对应实测值分别为604kn，874kn和1070kn，计算值与实测结果接近。如图9所示为总长19.2m的端承桩实测结果与本发明计算结果对比图，端承桩中性面位于桩尖以上大约1.2m平面处，在固结2，4及6.7年后中性面处下拉荷载计算值分别为921kn，1530kn和1914kn，对应实测值分别为908kn，1449kn和1707kn，计算值与实测数据比较接近。对比结果证明了本发明方法推导的正确性，也验证了本发明可以正确计算软土固结沉降并预测单桩基础的竖向长期承载响应。