0.   引言

群桩基础已成为软土地区广泛采用的基础形式之一。目前，诸多学者已对群桩基础进行深入的研究，例如程泽海等^[1]、崔春义等^[2]、梁发云等^[3]。但已有的研究^[1-3]仅针对饱和地基固结效应，无法考虑其长期流变性质。实际上，软土骨架在应力作用下，颗粒的重新排列和骨架错动具有明显的时间效应，土体变形和内部应力重新调整，表现为土体的流变现象^[4-5]。饱和软黏土中群桩基础的固结流变效应十分显著，为了避免过大的长期沉降和不均匀沉降，研究饱和软黏土中群桩基础的固结流变效应具有重要意义。针对群桩基础与黏弹性地基共同作用问题，陆建飞等^[6]利用边界积分方程法求解了黏弹性均质饱和半空间中的竖向受荷群桩。曾庆有等^[7]采用Mesri蠕变模型模拟土体的流变性质，提出了一种计算桩基长期沉降的方法。Ai等^[8]采用有限元-边界元耦合的方法，对层状黏弹性饱和土中竖向受荷群桩的时变行为进行了研究。然而，这些研究^[6-8]均基于经典整数阶导数模型来描述土的流变特性。近来，一些试验结果^[9-10]表明分数阶导数模型能以少量的参数准确描述软土蠕变曲线；为此凭借分数阶导数模型的优势，Ai等^[11-12]采用精细积分法求得了多层分数阶黏弹性地基的三维固结解答。

在作者已有工作的基础^[11-12]上，本文引入分数阶导数Merchant模型描述饱和软黏土的固结流变效应，再耦合边界元法与有限元法，对多层饱和软土与群桩相互作用进行研究，并讨论分数阶次对群桩时效行为的影响。

1.   竖向受荷群桩的有限元求解方程

对于本文所研究的桩基，桩长要远大于截面尺寸，因此本文将桩基础简化为一维有限元模型，以避免将桩基假设为三维模型时繁杂的运算。将轴向受荷单桩模拟为一维单元杆件时，2节点轴力杆中单元e的刚度矩阵为

式中，E_p和A_p分别代表单桩的弹性模量和横截面积，${L^{(e)}}$为桩单元e的长度。

考虑N根群桩中轴向受荷单桩的受力平衡条件，采用以下总集刚度矩阵可求得第j根桩的节点位移^[8]：

式中${\mathit{\boldsymbol{U}}_{{\rm{p}}j}} = \sum {_e\mathit{\boldsymbol{U}}_{{\rm{p}}j}^{(e)}}$，${\mathit{\boldsymbol{Q}}_{{\rm{p}}j}} = \sum {_e\mathit{\boldsymbol{Q}}_{{\rm{p}}j}^{(e)}}$和${\mathit{\boldsymbol{K}}_{{\rm{p}}j}} = \sum {_e\mathit{\boldsymbol{K}}_{{\rm{p}}j}^{(e)}}$分别代表了单桩总的节点位移向量、总的等效边界力向量和总集刚度矩阵；${\mathit{\boldsymbol{F}}_{{\rm{p}}j}}$为外荷载等效节点向量；$\mathit{\boldsymbol{U}}_{{\rm{p}}j}^{(e)} = {[\begin{array}{*{20}{c}} {{u_{z1}}}&{{u_{z2}}} \end{array}]^{\rm{T}}}$和$\mathit{\boldsymbol{Q}}_{{\rm{p}}j}^{(e)} = {[\begin{array}{*{20}{c}} {{Q_{z1}}}&{{Q_{z2}}} \end{array}]^{\rm{T}}}$分别为单元节点位移向量和等效节点向量，下标z1和z2分别对应单元e上下节点坐标；t为时间。

Laplace变换及其逆变换定义为^[13]

式中，$\tilde f(s)$是$f(t)$经过Laplace变换的函数，$s$为Laplace变换的参数，${\text{i}}$为单位虚数。

对式（2）进行Laplace变换，并总集N根桩的平衡控制方程可得

将式（5）写为矩阵形式，即得竖向受荷群桩的有限元求解方程：

式中，${\mathit{\boldsymbol{\tilde U}}_{\rm{g}}}(s)$和${\mathit{\boldsymbol{K}}_{\rm{g}}}$分别为Laplace变换域内群桩总集位移向量和群桩总集刚度矩阵，${\mathit{\boldsymbol{\tilde F}}_{\rm{g}}}(s)$和${\mathit{\boldsymbol{\tilde Q}}_{\rm{g}}}(s)$分别为Laplace变换域内群桩所受外荷载和等效边界力向量。

2.   分数阶黏弹性地基的边界元解

2.1   分数阶黏弹性模型

为描述土的流变固结特性，常采用元件模型以模拟土的松弛蠕变现象，而Merchant模型是岩土工程中应用最广泛的黏弹性模型之一^[11-12]。如图 1（a）所示，经典的整数阶Merchant模型由弹簧元件和Kelvin元件串联而成。该模型含有两个胡克弹性体H₀和H₁以及一个牛顿黏壶N。其中，E₀和E₁分别为弹簧元件和Kelvin元件对应的弹性模量，$\eta$为黏滞系数。

图  1  黏弹性本构模型

Figure  1.  Viscoelastic constitutive model

下载: 全尺寸图片 幻灯片

经典Merchant模型的应力–应变关系为^[8]

式中，$\sigma '$为有效应力量，$\varepsilon$为应变量。

如图 1（b）所示，将牛顿黏壶N替换为Abel黏壶A后，改进的分数阶Merchant模型应力–应变关系为^[11]

式中，${{\text{d}}^\alpha }/{\text{d}}{t^\alpha }$为$\alpha$阶分数阶导数（$0$＜$\alpha \leqslant 1$）。

此处分数阶导数采用Riemann-Liouville定义^[14]：当$t \in (0, \;T)$时，函数$f(t)$的分数阶导数为

式中，$\Gamma (x)$为Gamma函数，其定义如下：

根据分数阶导数的定义（9），Laplace变换域内的式（8）可写为

2.2   饱和软土地基的黏弹性解答

竖向荷载V作用下群桩–土相互作用问题如图 2（a）所示。假定刚性承台连接各桩，且未与地基相接触；地基划分为n层横观各向同性黏弹性土层，底部为固定不透水边界。其中，$\alpha$，${E_0}$，${E_1}$和$\eta$为式（8）中的分数阶黏弹性模型参数；${\nu _{\text{h}}}$为水平向应力引起竖直向应变的泊松比；${\nu _{{\text{vh}}}}$为竖直向应力引起水平向应变的泊松比；${k_{\text{v}}}$和${k_{\text{h}}}$分别为竖直向和水平向的渗透系数；$L$和$D$分别为桩长和桩直径。

图  2  层状地基与群桩模型示意图

Figure  2.  Diagram of layered soils and pile groups

下载: 全尺寸图片 幻灯片

根据弹性-黏弹性对应原理^[15]，式（12）所获得的应力–应变关系参数替换弹性地基中相关弹性参数，即可获得横观各向同性分数阶黏弹性饱和软土地基的解答^[11-12]。基于该解答，可导出地基内Laplace变换域内的柔度系数，即地基中$\varphi$处作用一z向单位环向荷载引起节点$\beta$处Laplace变换域内的竖向位移$\tilde G(\varphi , \beta , s)$和超孔压$\tilde R(\varphi , \beta , s)$。

如图 2（b）所示，本文基于有限单元法划分竖向受荷群桩，桩–土界面也被划分为相应的空间环状曲面单元。对于单元e，利用拉格朗日插值，环状表面上$\varphi$点的侧摩阻力${\tilde q^{(e)}}$可表示为^[8]

式中，$\boldsymbol{N}^{(e)}=\left[\begin{array}{ll}N_{1} & N_{2}\end{array}\right]$为两节点拉格朗日单元插值函数构成的向量，$\boldsymbol{q}^{(e)}=\left[\begin{array}{ll}q_{1} & q_{2}\end{array}\right]^{\mathrm{T}}$为单元e环状表面上下边缘所受到的边界力所构成的向量。

结合变换域内的柔度系数$\tilde G(\varphi , \beta , s)$和$\tilde R(\varphi , \beta , s)$，本文利用边界元法的基本思想，将桩–土界面上侧摩阻力的作用表达为边界积分方程，再将积分区间离散为桩单元e的长度${L^{(e)}}$，以便与群桩有限元解答（6）相互耦合。每一根桩含有M个${L^{(e)}}$段单元，则对于N根桩中所有单元e上的侧摩阻力${\tilde q^{(e)}}$引起$\beta$处的位移${\tilde u_z}(\beta , s)$和超孔压$\tilde \sigma (\beta , s)$，可通过如下积分求和得到^[8]

对式（13）进行Laplace变换，代入式（14）可得

式（15）按照式（1）~（6）的组装方式转化为矩阵表达：

式中，$\tilde{\boldsymbol{W}}_{\mathrm{g}}(s), \tilde{\boldsymbol{\sigma}}_{\mathrm{g}}(s)$和$\tilde{\boldsymbol{R}}_{\mathrm{g}}$分别为群桩节点竖向位移、超孔压和边界力向量，$\tilde{\boldsymbol{G}}_{\mathrm{g}}$和$\tilde{\boldsymbol{R}}_{\mathrm{g}}$分别是群桩地基对应的总位移柔度矩阵和总孔压柔度矩阵。

3.   有限元与边界元耦合

边界元方程（16）中的边界力向量$\tilde{\boldsymbol{q}}_{\mathrm{g}}(s)$代表了节点处的侧摩阻力值，而为了将连续分布的侧摩阻力等效至节点处，单元e上的等效节点力向量$\tilde{\boldsymbol{Q}}^{(e)}$可通过单元转换矩阵得到^[8]

式中，$\boldsymbol{T}^{(e)}=\int_{L^{(e)}} \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N }{\text{d}} L^{(e)}$代表单元e的转换矩阵。

借鉴组装有限元总刚度矩阵的方法，群桩有限元方程（6）中$\tilde{\boldsymbol{Q}}_{\mathrm{g}}(s)$可根据式（17）建立与边界元方程（16）中$\tilde{\boldsymbol{q}}^{(e)}(s)$的联系，即

式中，$\tilde{\boldsymbol{T}_{\mathrm{g}}}$为群桩总体转换矩阵。

由桩-土界面处的位移连续条件可知

于是联立式（6），式（16a），式（18）和Laplace变换后的式（19），可将有限元方程（6）与边界元方程（16a）耦合成群桩–土相互作用方程：

4.   群桩承台中的位移协调

求解群桩–土相互作用方程（20）得

式中，逆矩阵$\tilde{\boldsymbol{D}}_{\mathrm{g}}=\left(\boldsymbol{K}_{\mathrm{g}}+\boldsymbol{T}_{\mathrm{g}} \tilde{\boldsymbol{G}}_{\mathrm{g}}^{-1}\right)^{-1}$。

群桩的荷载V直接作用于承台，再由承台再分配给各基桩桩顶，因此只考虑$\tilde{\boldsymbol{U}}_{\mathrm{g}}(s)$中桩顶节点的位移；消去$\tilde{\boldsymbol{F}}_{\mathrm{g}}(s)$中的零值，同时缩减矩阵$\tilde{\boldsymbol{D}}_{\mathrm{g}}(s)$，那么式（21）可以改写为各基桩桩顶位移与桩顶力之间的关系：

式中，$\tilde{\boldsymbol{U}}_{\mathrm{d}}(s)$和$\tilde{\boldsymbol{F}}_{\mathrm{d}}(s)$分别为各基桩桩顶位移和桩顶力向量，$\tilde{\boldsymbol{D}}_{\mathrm{d}}$中保留了矩阵$\tilde{\boldsymbol{D}}_{\mathrm{g}}$中桩顶节点所对应行列的元素。

如图 3所示，刚性承台受到竖向集中荷载V时，各基桩的桩顶力有如下关系：

图  3  群桩承台受力变形示意图

Figure  3.  Stresses and deformations of cap for pile group

下载: 全尺寸图片 幻灯片

式中，$\boldsymbol{B} = [\begin{array}{*{20}{c}} 1&1& \cdots &1 \end{array}]$为N维行向量。

由承台底部与各基桩桩顶的位移协调条件可得

式中，${\tilde U_{\text{a}}}(s)$为Laplace变换域内承台的竖向位移。

联立式（22）~（24）得到

式（25）求得承台位移${\tilde U_{\text{a}}}(s)$，然后利用式（24）求得各基桩的桩顶位移$\tilde{\boldsymbol{U}}_{\mathrm{d}}(s)$；进而通过式（22）即可获得各基桩的桩顶反力$\tilde{\boldsymbol{F}}_{\mathrm{d}}(s)$。以上均为Laplace变换域内的解，通过数值逆变换^[16]，可求得其在物理域内的解。

5.   数值计算与分析

目前，分数阶黏弹性饱和地基与群桩共同作用研究鲜有报道。陆建飞等^[6]采用边界积分方程法推导出整数阶Merchant模型黏弹性饱和地基与群桩共同作用解答，显然该解答可视为本文分数阶黏弹性地基的特殊情况（$\alpha = 1.0$）。鉴于Ai等^[11-12]已验证了分数阶黏弹性饱和地基解答的正确性，故本文仅对整数阶黏弹性地基中的群桩解答进行验证，以表明本文理论的正确性。验证算例中，泊松比$\nu = 0.49$，渗透系数$k =$ $1.0 \times {10^{ - 7}}{\text{m}}/{\text{s}}$，桩间距与桩径比${s_{\text{p}}}/D = 5$，长径比$L/D = 15$，${R_{\text{p}}} = 0.5D$，${E_{\text{p}}}/{E_0}{\text{ = 250}}$，${E_{\text{p}}}/{E_1} = 1000$，${R_{\text{p}}}{E_{\text{p}}}/k\eta = 5.0 \times {10^4}$。无量纲时间$\tau = 2{\zeta _{\text{v}}}{E_1}kt/{\gamma _{\text{w}}}{R_{\text{p}}}^2$，${\zeta _{\text{v}}} = (1 - \nu )/(1 - 2\nu )$，水重度${\gamma _{\text{w}}} = 0.01{\text{ MN}}/{{\text{m}}^3}$；第i根桩的无量纲轴力为$V_i^* = 9{V_i}/V$；承台无量纲位移$U_{\text{a}}^* = 9{E_1}{R_{\text{p}}}{U_{\text{a}}}/V$。

表 1列举了不同时刻群桩刚性承台的位移，而各桩沿深度变化的轴力见图 4。由图 4可知，本文基于边界元-有限元耦合法求解群桩时效行为（桩单元离散数量M = 10），与文献[6]积分方程法所得结果吻合得较好，这也证明本文将桩模拟为一维杆件的合理性。

表  1  不同时刻下刚性承台无量纲位移

Table  1.  Non-dimensional displacements at different time

文献	${\tau _1} = 5.1 \times {10^{ - 3}}$	${\tau _2} = 5.1 \times {10^4}$
陆建飞等[6]	0.0384	0.1737
本文	0.0402	0.1749

下载: 导出CSV 

| 显示表格

图  4  不同时刻下黏弹性地基中桩身轴力

Figure  4.  Axial forces of each pile in viscoelastic soils

下载: 全尺寸图片 幻灯片

Ai等^[8]研究了基于整数阶黏弹性地基模型中群桩时变行为，并分析了地基黏滞系数、横观各向同性参数、分层性、群桩长径比和桩间距的影响，因此，本文仅对分数阶次的影响进行研究。关于分数阶黏弹性参数如何选取，可参考文献[9，10]中的试验方法：利用三轴固结试验获得土的蠕变曲线，再根据应力松弛公式拟合分数阶黏弹性模型参数。为了研究分数阶次的影响，本文采用3个不同阶次的分数阶黏弹性地基模型（$\alpha =$0.7，0.8，1.0）模拟饱和软土的固结流变特性。算例考虑竖向荷载$V$作用在群桩承台上，群桩布置如图 4所示。三层饱和软土地基参数如表 2所示，厚度分别为${h_1} = {h_2} = 15D$，${h_3} = 200D$，桩直径$D =$ $0.5{\text{ m}}$，桩间距${s_{\text{p}}}/D = 4$，桩身模量${E_{\text{p}}}/{E_{{\text{0v1}}}} = 2000$，长径比$L/D = 30$，泊松比分别为${\nu _1} = 0.25$，${\nu _2} =$ $0.3889$和${\nu _3} = 0.3158$，无量纲时间${\tau ^{\text{*}}} = {E_{{\text{0v1}}}}{k_1}t/{\gamma _{\text{w}}}{D^2}$，无量纲孔压$\sigma (z){\text{π }}DL/V$，无量纲承台位移${w^*} =$ $9D{E_{{\text{0v1}}}}{U_{\text{a}}}/V$。其它参数为：${k_1} = 2.0 \times {10^{ - 7}}{\text{m}}/{\text{s}}$，${k_2} =$ $6.0 \times {10^{ - 8}}{\text{m}}/{\text{s}}$，${k_3} = 1.0 \times {10^{ - 7}}{\text{m}}/{\text{s}}$，${\zeta _{\text{G}}} = {G_{\text{v}}}/{E_{\text{v}}}$，${\zeta _{\text{E}}} =$${E_{\text{h}}}/{E_{\text{v}}}$。

表  2  多层黏弹性地基参数

Table  2.  Parameters of layered viscoelastic soils

土层	${\zeta _E}$	${\zeta _G}$	${E_{{\text{0v}}i}}/{\text{MPa}}$	${E_{{\text{1v}}i}}/{\text{MPa}}$	${\eta _{1i}}/({\text{MPa}} \cdot {{\text{s}}^{ - 1}})$
1	2.0	0.40	10	8	2×107
2	1.5	0.36	20	5	5×107
3	1.7	0.38	25	5	1×108

下载: 导出CSV 

| 显示表格

图 5，6分别给出了承台竖向位移和群桩桩顶轴力分配随时间的变化。相较于整数阶黏弹性地基（$\alpha = 1.0$），本文理论可根据软土蠕变曲线确定分数阶地基参数，从而更好地预测饱和软黏土的蠕变对群桩基础时效行为的影响。随着时间的变化，承台竖向位移逐渐增大，而边桩与中心桩的轴力逐渐转移由角桩承担；同时，随着分数阶次减小，固结流变所需时间呈指数级增长，而且软土蠕变过程的时间越长，承台竖向位移和轴力分配达到稳定的时间就越长；此外，不同的分数阶次不影响初始阶段和最终阶段承台竖向位移和各基桩桩顶的轴力。现有的建筑桩基技术规范^[17]对桩基础的设计计算进行了简化，只考虑了最终阶段的沉降和桩顶轴力，而本文方法能够考虑沉降和桩顶轴力随时间的变化，在各个时间段对群桩基础进行沉降和受力分析。显然，这对工程精细化设计以及全过程质量监测具有重要的指导意义。图 7表明，超孔隙水压力的消散不受分数阶次的影响，这是因为饱和软黏土中超孔隙水压力的消散集中在固结过程；而蠕变过程由分数阶次描述，故图 7中不同分数阶次同一时刻的超孔隙水压力曲线重合。

图  5  承台竖向位移随时间的变化

Figure  5.  Variation of vertical displacement of the pile cap with time

下载: 全尺寸图片 幻灯片

图  6  群桩桩顶轴力分配随时间的变化

Figure  6.  Variation of distribution of axial force at pile top with time

下载: 全尺寸图片 幻灯片

图  7  超孔隙水压力沿中心桩轴线变化情况

Figure  7.  Distribution of pore water pressure along axis of central pile

下载: 全尺寸图片 幻灯片

6.   结论

为了更好地模拟饱和软黏土的蠕变对群桩基础时效的影响，本文引入分数阶导数Merchant模型研究层状饱和软黏土中群桩基础的时效行为，数值算例的结果表明：

（1）分数阶次的大小显著影响了土体固结流变的进程。分数阶次越小，软土流变过程的时间越长，群桩轴力分配达到稳定所需要的时间越长。

（2）对于分数阶Merchant模型而言，分数阶次的改变，不影响初始阶段和最终阶段的桩顶轴力值。

（3）相较于现有规范^[17]，本文方法能够综合考虑软土固结流变过程对群桩基础沉降和各基桩轴力变化的影响，这有助于工程精细化设计以及全过程的质量监测。