深厚饱和黏性土地基中桩筏基础沉降的时间效应是一个非常复杂的三维问题.由于黏性土渗透性较低，上部荷载引起的土体超孔隙水压力无法在短时间内完全消散，基础体系的内力和变形在施工结束后较长一段时间内仍将不断发生变化.目前，多数学者^[1-4]采用有限元数值方法对桩筏基础工作性状的变化规律以及基础内力和变形受时间的影响等进行研究，但对施工期间和施工结束后基础沉降速率随时间变化过程的研究则相对较少.Small等^[5]基于ABAQUS软件建立了考虑流固耦合作用的桩筏基础三维有限元-无限元分析模型，并利用该模型模拟基础施工和正常使用阶段的沉降发展过程.Sales等^[6]将桩筏基础用二维轴对称等代墩基模型近似表示，同样采用固结有限元理论预测基础沉降，但桩群与桩间土构成的等代墩刚度不易确定.虽然有限元方法能模拟基础变形随地基固结的发展规律，但建模过程比较复杂，计算效率和收敛性仍存在一定问题，较难在实际工程中得到推广应用.文献[7-9]在对渗流边界条件进一步简化的基础上，采用计算过程更加简便的半解析-半数值方法进行群桩基础固结沉降分析.然而，由于孔压分布和渗流边界条件十分简单，不能反映桩-土共同作用的影响，其适用性受到明显限制.

Davis等^[10]的研究表明，深厚黏性土地基中按一维固结理论计算的基础沉降在多数情况下与实际不符，二维和三维固结理论的计算结果明显优于前者，计算精度在一定程度上能满足实际工程应用要求.本文将桩-土非线性共同作用模型与二维Terzaghi-Rendulic固结理论相结合，建立了求解任意时刻基础体系固结沉降的简化且实用的分析方法.通过与Biot固结有限元分析的精确解答进行对比，验证了本文方法和假定的合理性.最后，将本文方法用于分析2个桩筏基础实际工程，模拟施工和正常使用阶段基础沉降随时间的发展过程，并将计算结果与实测值进行对比.

1 计算方法 1.1 分析步骤和计算假定

桩筏基础总沉降通常包括瞬时沉降、固结沉降和次固结沉降.3类沉降过程往往交替发生且相互影响，在时间上并不能进行严格的区分.从工程实用角度来看，对上述沉降分别进行计算有助于加深对地基变形机理的认识，同时也可使分析过程得到简化.已有研究^[11]表明，对于一般黏性土地基来讲，工作荷载作用下群桩基础因地基蠕变产生的沉降通常很小，除非荷载水平进一步增大，因此本文不考虑地基蠕变对基础沉降的影响.假定基础总沉降s_tot由施工阶段的瞬时(不排水)沉降s_u和正常使用阶段的固结沉降s_d两部分组成.主要分析过程包含以下步骤:

(1) 采用桩-土非线性共同作用模型计算桩筏基础的瞬时沉降s_u和最终沉降s_tot(固结总沉降为s_tot-s_u)，有关该模型的详细计算过程参见本文第1.2节.根据不排水分析求得基础-地基相互作用力和采用Mindlin应力解^[12]公式计算基础中心对称剖面(如图 1所示的xoz平面)内地基附加应力.

(2) 采用文献[13]的建议，假定初始超孔隙水压力完全由地基中的附加应力引起，并且在基础中心对称平面内沿竖直和水平方向同时发生消散，即满足二维Terzaghi-Rendulic固结理论要求.如图 2所示，在xoz平面内确定二维固结计算区域边界，水平计算宽度^[14]取3.5B_r(B_r为筏板宽度)，竖直方向由地表(或筏板底面)到桩端以下砂层或基岩顶面.黏土层表面(顶面边界)和右侧(远离基础一侧竖直边界)为透水边界，底部边界透水性与下卧土层渗透性有关，基础左侧边界(基础中心竖直边界)为完全不透水边界.

(3) 群桩基础固结沉降主要由桩端以下土体压缩变形引起的，地基附加应力和超孔隙水压力具有明显的空间分布特征，基础中心和边缘位置孔压大小和消散速率并不一致^[15].考虑到筏板刚度对基础沉降和应力分布的调整作用，本文将基础半宽(B_r/2)范围内桩端下卧土层的平均固结度用于固结沉降计算，其他区域土体作为一般的排水渗流通道.

假设图示平均固结度计算区域内共包含k个计算结点，土结点i在任意时刻t的固结度

$ {U_{i, t}} = \frac{{{u_{i, 0}} - {u_{i, t}}}}{{{u_{i, 0}}}}, \;\;\;\;i = 1, 2, \cdots, k $

(1)

式中:u_{i, 0}为初始时刻结点i的超孔隙水压力; u_{i, t}为t时刻结点i的超孔隙水压力.

为考虑计算区域内各结点孔压的影响，桩端以下土层平均固结度

$ {U_z} = \frac{{\sum\limits_{i = 1}^k {\left( {{u_{i, 0}} - {u_{i, t}}} \right)} }}{{\sum\limits_{i = 1}^k {{u_{i, 0}}} }} $

(2)

施工期结束后任意时刻t桩筏基础的工后固结沉降s_tot可由下式计算:

$ {s_{{\rm{tot}}}} = {s_{\rm{u}}} + {s_{\rm{d}}}{U_z} $

(3)

1.2 桩土非线性共同作用模型

对大部分采用桩筏基础的建筑物来讲，群桩和上部结构对筏板刚度的贡献程度很大，为简化建模过程和提高计算效率，本文假定筏板为完全刚性，采用以下步骤计算桩筏基础的瞬时沉降和最终沉降:

对筏板和桩基进行计算单元剖分，建立反映桩-土-筏共同作用的土体位移方程，如下所示:

$ {\mathit{\boldsymbol{w}}_{\rm{s}}} = \mathit{\boldsymbol{F}}{\mathit{\boldsymbol{P}}_{\rm{s}}} $

(4)

式中:F为土体柔度系数矩阵，矩阵各元素反映桩、土、筏之间的相互作用; w_s和P_s分别为土体竖向位移和土体反力列向量.筏-土相互作用柔度系数按下式^[16]计算:

$ {F_{{\rm{sr}}}} = \frac{{\left( {1 - {\nu _{\rm{s}}}} \right)\left( {1 - \exp \left( { - {H_{\rm{s}}}/2{r_{{\rm{equ}}}}} \right)} \right)}}{{4{G_{\rm{s}}}{r_{{\rm{equ}}}}}} $

(5)

式中:H_s为可压缩土层深度，该深度以下为不可压缩土层; r_equ为筏板单元等效半径(按方形筏板与圆形筏板单元面积相等原则计算得到); G_s和ν_s分别为土体剪切模量和泊松比.当筏板位于非均质地基和分层地基中时，需按文献[17]的建议对G_s和ν_s进行修正.

筏-土-筏相互作用柔度系数按下式^[18]计算:

$ {F_{i, j}} = \frac{{\left( {1 - \nu _{\rm{s}}^2} \right)}}{{{\rm{ \mathsf{ π} }}{E_{\rm{s}}}{r_{{\rm{equ}}}}}}{\sin ^{ - 1}}\left( {\frac{{{r_{{\rm{equ}}}}}}{{{l_{i, j}}}}} \right), \;\;\;i \ne j $

(6)

式中:l_{i, j}为筏板单元i和j中心点之间的水平距离; E_s为土体弹性模量.

桩-土、桩-土-桩和桩-土-筏相互作用柔度系数均采用Mindlin位移解^[19]计算.对于非均质和分层地基，或可压缩土层深度有限的情况，需按文献[16, 20]建议的方法对上述土体柔度系数进行修正.

根据桩-土-筏体系的静力平衡条件和位移协调要求，可得到桩筏基础体系作用荷载与基础位移的数学关系式，即

$ {\mathit{\boldsymbol{K}}_{{\rm{tot}}}}{\mathit{\boldsymbol{w}}_{\rm{p}}} = {\mathit{\boldsymbol{Q}}_{{\rm{top}}}} $

(7)

式中:K_tot为桩-土-筏体系整体刚度矩阵，是群桩竖向刚度矩阵K_p和地基土刚度矩阵K_s之和，K_p通过杆系有限元方法计算得到，其维数应扩充至与所有桩、筏结点数目相同(筏-土结点所在行列元素均用0代替)，K_s通过对土体柔度系数矩阵F求逆得到; w_p为所有桩、筏计算结点的竖向位移列向量; Q_top为基础外荷载列向量.

对于筏板完全刚性的桩筏基础，已知筏板和桩顶结点的竖向沉降，按乘大数法^[21]对式(7)两侧基础体系刚度矩阵K_p+K_s和荷载向量Q_top进行修正.若已知筏板和桩顶竖向沉降为w₀，将K_tot中对应筏板和桩顶结点位置的主对角线元素k_tot修改为k_totb(b可取10¹⁶~10¹⁸)，将Q_top中相同位置的元素q_top改为w₀bk_tot，其余元素均为零.此时，式(7)在基础顶部满足位移边界条件，对该式进行计算可获得基础各结点处的竖向位移和相互作用力，以及基础顶面作用的总荷载.

考虑桩-土界面可能出现滑移和土体进入塑性，按位移增量法求解不同荷载作用下的基础内力与变形.在每个增量步内，将求得的结点反力增量与之前各增量步内的结点反力增量累加，并与(桩-土、筏-土)界面反力极限值进行对比.若计算结点力超过界面极限反力，用极限反力值代替该结点力，并将土体柔度系数矩阵中相同位置主对角线元素补充一大数，非主对角线元素置零(进入极限状态的结点不再引起相邻结点的附加位移)，重复该过程直至完成所有增量步.通过上述求解步骤，可得到桩筏基础非线性荷载-沉降关系，以及各级荷载作用下的结构-地基相互作用力.

1.3 二维固结算法实现

本文采用Terzaghi-Rendulic固结理论进行平面应变条件下的地基土固结排水分析，计算过程不考虑土骨架变形与渗流的耦合作用，直接求解任意时刻地基不同位置结点的孔隙水压力.平面应变问题的二维固结方程如下所示:

$ {v_x}\frac{{{\partial ^2}u}}{{\partial {x^2}}} + {v_z}\frac{{{\partial ^2}u}}{{\partial {z^2}}} = \frac{{\partial u}}{{\partial t}} $

(8)

式中:t为固结时间; u为土体超孔隙水压力; v_x和v_z分别为x和z方向的固结系数，文献[13]证实二维条件下的固结系数与一维和三维时具有相同的数学表达式.

由于初始孔隙水压力值在整个计算域内并不具有明确的数学表达式，需采用数值算法求解式(8).显式差分格式虽然求解过程简便，但需要满足稳定性条件，当渗透系数较大或网格划分较密时，计算精度和收敛性较难保证.因此，本文采用交替方向的隐式差分算法^[22-24]进行编程求解.

由t时刻的孔隙水压力计算t+Δt时刻的孔隙水压力，在t和t+Δt时刻之间增加一个t+Δt/2时刻.对于t到t+Δt/2时刻，将$\frac{{{\partial ^2}u}}{{\partial {z^2}}} $取成隐式格式，把$ \frac{{{\partial ^2}u}}{{\partial {x^2}}}$取成显式格式.对图 3中的结点u_{i, j, t}，式(8)可以写成

$ {v_x}{\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right)_{i, j, t}} + {v_z}{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)_{i, j, t + \Delta t/2}} = {\left( {\frac{{\partial u}}{{\partial t}}} \right)_{i, j, t + \Delta t/2}} $

(9)

式(9)的差分方程形式为

$ \begin{array}{l} \frac{{{v_x}}}{{{\rm{d}}{x^2}}}\left( {{u_{i, j - 1, t}} - 2{u_{i, j, t}} + {u_{i, j + 1, t}}} \right) + \frac{{{v_z}}}{{{\rm{d}}{z^2}}}\left( {{u_{i - 1, j, t + \Delta t/2}} - } \right.\\ \;\;\;\;\;\;\;\left. {2{u_{i, j, t + \Delta t/2}} + {u_{i + 1, j, t + \Delta t/2}}} \right) = \\ \;\;\;\;\;\;\;\;\frac{2}{{\Delta t}}\left( {{u_{i, j, t + \Delta t/2}} - {u_{i, j, t}}} \right) \end{array} $

(10)

式(10)可整理为

$ \begin{array}{l} \frac{{{v_z}}}{{{\rm{d}}{z^2}}}{u_{i - 1, j, t + \Delta t/2}} - 2\left( {\frac{{{v_z}}}{{{\rm{d}}{z^2}}} + \frac{1}{{\Delta t}}} \right){u_{i, j, t + \Delta t/2}} + \\ \;\;\;\;\;\;\;\;\;\frac{{{v_z}}}{{{\rm{d}}{z^2}}}{u_{i + 1, j, t + \Delta t/2}} = - \frac{{{v_x}}}{{{\rm{d}}{x^2}}}{u_{i, j - 1, t}} + \\ \;\;\;\;\;\;\;\;\;2\left( {\frac{{{v_x}}}{{{\rm{d}}{x^2}}} - \frac{1}{{\Delta t}}} \right){u_{i, j, t}} - \frac{{{v_x}}}{{{\rm{d}}{x^2}}}{u_{i, j + 1, t}} \end{array} $

(11)

将式(11)应用于xoz平面内所有内部结点(不包括所有边界结点)，所有方程联立得到反映t时刻与t+Δt/2时刻孔隙水压力的方程组，该方程组孔隙水压力未知量的个数将超过所有内部结点个数.此时，引入已知边界条件，由于各结点孔隙水压力值与上下左右4个结点孔隙水压力值均有关，当中心结点与边界结点相邻时，需按如下方式处理式(10):

当边界完全透水时，边界结点孔隙水压力u=0;当边界不透水时，边界结点孔隙水压力变化为0，即$\frac{{\partial u}}{{\partial t}} = 0 $.

将上述边界条件写成差分格式，并代入式(10)，孔隙水压力未知量的个数将缩减为与内部结点个数相等，方程组有唯一解.需要说明的是，对于五点十字型差分方程，边界上4个角点在求解过程中并不出现，可看作孔隙水压力已知的边界结点.

为保持对称性，在对t+Δt/2到t+Δt时刻进行计算时，将$ \frac{{{\partial ^2}u}}{{\partial {x^2}}}$取成隐式格式，$\frac{{{\partial ^2}u}}{{\partial {z^2}}} $取成显式格式.同样可将式(8)写成

$ {v_x}{\left( {\frac{{{\partial ^2}u}}{{\partial {x^2}}}} \right)_{i, j, t + 1}} + {v_z}{\left( {\frac{{{\partial ^2}u}}{{\partial {z^2}}}} \right)_{i, j, t + 1/2}} = {\left( {\frac{{\partial u}}{{\partial t}}} \right)_{i, j, t + 1}} $

(12)

与上述过程类似，式(12)的差分方程形式为

$ \begin{array}{l} \frac{{{v_x}}}{{{\rm{d}}{x^2}}}{u_{i, j - 1, t + 1}} - 2\left( {\frac{{{v_x}}}{{{\rm{d}}{x^2}}} - \frac{1}{{\Delta t}}} \right){u_{i, j, t + 1}} + \frac{{{v_x}}}{{{\rm{d}}{x^2}}}{u_{i, j + 1, t + 1}} = \\ \;\;\;\;\;\;\;\; - \frac{{{v_z}}}{{{\rm{d}}{z^2}}}{u_{i - 1, j, t + 1/2}} + 2\left( {\frac{{{v_z}}}{{{\rm{d}}{z^2}}} - \frac{1}{{\Delta t}}} \right){u_{i, j, t + 1/2}} - \\ \;\;\;\;\;\;\;\;\frac{{{v_z}}}{{{\rm{d}}{z^2}}}{u_{i + 1, j, t + 1/2}} \end{array} $

(13)

通过求解式(12)和式(13)，可获得t+Δt时刻的u值.该差分格式解法在数值上是无条件稳定的，故可以采用较大的时间步长以提高算法的计算效率.

2 计算程序验证 2.1 均布条形荷载地基固结计算(二维问题)

文献[25]采用Biot固结有限元方法对均布条形荷载下的地基固结变形问题进行计算.饱和地基土层厚度为12 m，顶面透水，底面不透水，其上作用20 kPa竖向均布荷载，荷载宽度为6 m，土体弹性模量和泊松比分别为4 MPa和0.3，渗透系数为2.5×10^-6 cm·s^-1.图 4~6给出了本文方法和文献[25]有限元方法计算的不同位置土体孔隙水压力与时间的关系曲线.不难看出，虽然本文方法未考虑固结过程中总应力随时间的变化，理论上不如Biot固结理论严密，但对于基础竖向变形问题的计算值与严格的有限元解答之间的差异并不十分明显，误差大小完全能够满足工程应用要求.

2.2 桩筏基础固结沉降计算(三维问题)

文献[5]对桩筏基础在均布竖向荷载作用下的固结沉降问题进行了数值计算分析.如图 7所示，方形筏板宽度为5.0S(S为桩间距)，厚为0.25S，桩长和桩径分别为6S和0.2S，筏板底面与地表面齐平且完全接触，计算土层厚度为15S，底部假定为刚性不可压缩土层.桩与筏板的弹性模量为土体弹性模量的1 000倍，桩筏基础和土体的泊松比分别为0.2和0.3.

图 8给出了按不同方法计算的基础中心沉降固结度，时间因子T_v的表达式为

$ {T_{\rm{v}}} = \frac{{vt}}{{{B^2}}} $

(14)

式中:v为固结系数，由$v = \frac{{kE'\left( {1-\nu '} \right)}}{{{\gamma _{\rm{w}}}\left( {1 + \nu '} \right)\left( {1-2\nu '} \right)}} $计算得到，其中γ_w为土中水的容重，k为土体渗透系数，E′、ν′分别为排水条件下的土体弹性模量和泊松比; B为筏板宽度的一半(2.5S).

可以看出，本文方法与有限元数值计算结果在总体上具有较好的一致性，最大偏差不超过20%，表明本文计算方法与假定具有一定合理性.

3 工程实例分析 3.1 案例1

如图 9所示，法兰克福Messe-Torhaus大厦桩筏基础平面尺寸24.5 m×17.5 m，厚2.5 m，弹性模量34 MPa，筏板埋深5.5 m，筏板以下为42根钻孔灌注桩，桩长和桩径分别为20.0 m和0.9 m，桩身弹性模量为23.5 GPa，沿筏板长边和短边的桩间距为3.5 m和2.9 m.

该场地黏土为典型的Gibson地基土，土体弹性模量与深度之间具有较明确的数学关系.文献[5]采用阶梯状分布的土体弹性模量假定进行有限元建模分析(见图 10)，本文用线性关系对其进行拟合，有

$ \left\{ \begin{array}{l} {E_{\rm{u}}} = 13.7 + 1.8z\\ E' = 10.5 + 1.4z \end{array} \right. $

(15)

式中:E_u和E′分别为土体不排水弹性模量和有效弹性模量，MPa; z为深度，m.

类似地，土体不排水抗剪强度c_u(单位:kPa)与深度z(单位:m)具有如下关系:

$ {c_{\rm{u}}} = 127 + 3.93z $

(16)

按照文献[26]的建议，基础瞬时沉降计算取土体泊松比为0.5，桩侧摩阻力系数α为0.56，桩端和筏板底面地基极限压应力分别为9c_u和6c_u，基础最终沉降计算时土体泊松比为0.15.由此可得，基础瞬时沉降和最终沉降分别为53 mm和110 mm，在施工期结束后基础产生的平均固结沉降为57 mm.

地基二维固结计算模型水平宽度为60 m，竖向宽度为90 m(往下至基岩顶面).土体渗透系数k为0.007 7 m·月^-1，体积压缩系数m_v为0.022 m²·MN^-1.差分法网格间距设置为1 m，时间步长为1月.表 1给出了本文计算的不同时刻土体平均固结度和相应的桩筏基础沉降.不难看出，基础底面地基土固结速率相对较快，在施工期结束后约30个月时固结度可达90%以上，5年后基础固结沉降过程基本完成.

图 11为计算和实测的基础平均沉降与时间的关系曲线.可以看出，本文与有限元方法^[5-6]计算结果总体上较为一致，在基础沉降初期(30个月以前)与实测值相差也很小.随后，基础实测沉降大于计算值，这可能与该桩筏基础受到相邻高层建筑(间隔10 m)基础荷载影响，基底竖向附加应力水平进一步提高，进而导致地基土发生蠕变变形有关，故3种理论方法(仅考虑主固结变形)的计算沉降量均比实测值偏小.

3.2 案例2

文献[27]提供了伦敦某高层建筑桩筏基础的现场实测资料.筏板总面积约为621 m²，埋深为8.8 m，筏板弹性模量为34 GPa.筏板下布置51根扩底灌注桩，桩身直径为0.91 m，桩长为24.80 m，桩端扩径约为2.44 m，弹性模量为23.5 GPa.

如图 12所示，地表 2.0 m深度为填土，填土以下为厚度约89 m的伦敦黏土.由文献[27]可知，该场地土体不排水弹性模量E_u与深度z具有以下关系^[27]:

$ {E_{\rm{u}}} = 10 + 5.2z $

(17)

土体不排水抗剪强度c_u=E_u/600.在排水条件下，土体有效弹性模量E′=0.75E_u，泊松比为0.12.基础施工过程荷载工况列于表 2.

由于开挖深度较大，地基土在不同施工阶段所处的应力状态不同，需考虑基坑开挖卸荷对土体变形参数的影响.本文基础施工阶段和结构施工阶段A过程采用土体卸荷-再加荷模量^[6]进行沉降计算; 结构施工阶段B采用不排水弹性模量进行计算，桩侧极限摩阻力和桩端极限阻力等参数的确定方法同案例1;最终沉降采用土体有效参数进行计算.由此可以得到不同阶段桩筏基础的平均沉降分别为2.8 mm(基础施工阶段和结构施工阶段A产生的沉降)、10.5 mm(结构施工阶段B产生的沉降)、19.8 mm(基础最终长期沉降)、9.3 mm(工后基础固结沉降).

地基二维固结计算模型水平宽度为84.0 m，竖向宽度为90.7 m(自筏板底面到基岩顶面).土体渗透系数k为4.09×10^-4 m·月^-1，体积压缩系数m_v为0.007 2 m²·MN^-1.表 3给出了不同时刻的地基固结度和基础平均沉降.可以看出，施工结束后约4.6年(55个月)桩端下卧土层的平均固结度还不足60%，距离基础沉降完成还需较长一段时间，固结速率明显不如案例1，这与该场地土体渗透系数偏低有关.图 13为计算和实测基础平均沉降与时间的关系曲线.可以看出，本文方法由于考虑了基坑开挖卸荷对土体变形参数的影响，施工阶段计算沉降值与实测吻合较好，工后沉降速率比实测结果略微偏小，但与三维有限元方法比较接近.

4 结语

将桩-土非线性共同作用模型与二维固结理论结合，建立求解饱和深厚黏土地基中桩筏基础固结沉降的实用计算方法.算例分析结果表明，本文模型的建立与求解过程比固结有限元法大为简化，计算的孔隙水压力和基础沉降固结度与有限元法总体较为接近，初始超孔隙水压力分布和渗流边界假定具有一定合理性.将本文方法用于模拟实际桩筏基础沉降量随时间的发展过程，得到了与实测值较为接近的预测结果，表明本文方法具有一定的工程应用价值.需要指出的是，本文方法不能考虑软土地基蠕变变形的影响.