载体桩是近年来发明并得到广泛应用的一种新桩型，该技术通过夯击能量以填料作为介质，挤密加固桩端周围土体，形成桩端扩大载体，从而大大提高单桩竖向承载力^[1]。

载体桩的施工工艺(如图 1)可以概括为：利用柱锤夯击成孔；将护筒压入地基；分批向孔内投入水泥砂拌合物填充料反复夯实挤密，并通过三击贯入度进行密实度控制，当三击贯入度满足设计要求后, 形成载体；根据需要可放置钢筋笼；灌注混凝土或放置预应力管桩而形成桩体桩^[1]。从载体桩施工工艺来看，载体桩具有以下优点：通过柱锤夯实填料和桩端土体，提高桩端周围土体的密实度和强度；桩端由水泥砂拌合物、挤密土体、影响土体形成刚度渐变的复合载体，这种变刚度的设计方式，使得载体的刚度从上往下逐渐减小，从而更加有利于桩体将竖向荷载传递到深部地基土；绿色施工，不需要采用大量泥浆护壁，同时施工不出土，保护了建筑环境；采用三击贯入度来控制载体施工质量，有利于保证载体能有效提供承载力^[1]。

图 1

图 1 载体桩施工过程 Fig. 1 Construction process of carrier pile

载体桩技术的核心在于侧限约束下的土体密实形成载体。该技术的研究对象是土体密实理论，也就是研究桩身端以下一定范围内挤密土体和影响土体的物理力学性质的变化。在设计中合理考虑挤密效应至关重要，如果高估了挤密效应，将会给工程带来安全隐患，反之，如果过低估计了载体的作用，将会造成工程造价大大提高。因此，如何从理论角度精确计算桩端挤密效应，成为制约载体桩推广和发展的一个主要因素。目前，关于载体桩的研究主要停留在试验阶段，在理论方面的研究相对较少。王建安等^[2]开展了载体桩在粉质黏土中的现场试验研究，对比了普通灌注桩和载体桩复合地基的竖向承载特性，结果表明，载体桩复合地基造价经济，同时能够提高竖向承载力。于长杰^[3]开展了软黏土中载体桩沉降控制研究，对比分析了有桩帽和无桩帽载体桩复合地基的沉降控制效果。周斌^[4]探讨了载体桩复合地基在高速铁路地基中对地基沉降变形控制的效果。罗浩^[5]开展了载体桩复合地基在高速铁路地基中的现场试验以及施工工艺的研究。李建强等^[6]开展了载体桩复合地基数值模拟研究，探讨了荷载、褥垫层、桩长、桩径、载体直径等因素对载体桩复合地基承载特性的影响，然而，该研究并未考虑桩端挤密效应的影响。张培成等^[7]开展了饱和软土地基载体桩现场试验研究，证明载体桩能够完成大承载力的要求，并可以在地基中深部有相对硬层的场地中取代预应力管桩。仇凯斌等^[8]开展了载体桩承载性状的有限元分析，通过现场取样和模型试验确定夯扩挤密区范围，利用线性插值确定挤密区土体的力学参数，获得了能够考虑挤密效应的载体桩承载特性，然而，该方法仅仅是从经验性的角度来考虑挤密效应，并未提出一种较为严格和准确的理论分析方法。

扩孔理论作为岩土力学中一种简单有效的力学模型^[9-19]，广泛用于解决岩土工程中如：圆锥贯入试验、旁压试验、扁铲侧胀试验、沉桩挤土效应、桩端承载力、隧道开挖、钻井失稳、压力注浆等问题。该模型数学求解接单，模型物理概念清楚，较为实用。笔者采用球孔扩张力学模型来模拟柱锤夯扩填料过程中的挤密效应，软黏土采用修正剑桥模型的本构关系来模拟，建立球孔扩张偏微分控制方程组，通过相似变换的求解技术将偏微分方程转化为常微分方程组，利用MATLAB中ODE45数值求解技术获得常微分方程组的数值解。基于数值解，探讨球孔扩张过程中孔周土体的挤密范围和挤密区土体力学性质的变化，从理论角度揭示载体桩载体成形过程中的加固机理，为建立载体桩竖向承载力计算方法提供理论基础。

1 载体桩载体扩孔成形力学模型

如图 2(a)所示，载体桩桩端夯扩填料的过程中，将会在端部形成一个类似球体的结构，该球体挤扩周围的土体，从而在桩端附近形成一个挤密区。这个过程与球孔扩张理论模型较为相似，因此，采用如图 2(b)所示的球孔扩张理论模型来模拟夯扩填料挤土的过程。在夯扩前，球孔的初始半径a₀ =0，随着夯扩的加剧，球孔从0不断扩张，最终扩张为半径为a的一个球体。球孔扩张理论模型是无限大土体中的力学模型，而实际问题中挤土主要在桩端以下部分，因此，球孔扩张区域为图 2(a)中桩端以下(虚线以下)部分。假定桩端的球体通过桩端左右两个角点，球体的半径为a，因此，该球体(虚线以下的部分，非完整球体)的体积可以表示为

$ {V_{\rm{s}}} = \frac{2}{3}\pi \left[ {{a^3} + \left( {{a^2} - r_0^2} \right)3/2} \right] $

(1)

图 2

图 2 载体桩载体扩孔成形力学模型 Fig. 2 Mechanics model of the cavity expansion formation of carrier pile

式中：r₀为桩体的半径。

球体的体积V_s可以根据实际的填料体积来确定，因此，通过式(1)，可以确定球体扩张后的半径a。式(1)是关于半径a的隐函数，无法求得闭合解析解，因此，可以利用MATLAB计算出球体体积V_s和球体半径的关系，为了方便，可以将式(1)进行归一化，得到

$ \frac{{{V_{\rm{s}}}}}{{r_0^3}} = \frac{2}{3}\pi \left\{ {{{\left( {\frac{a}{{{r_0}}}} \right)}^3} + \left[ {{{\left( {\frac{a}{{{r_0}}}} \right)}^2} - 1} \right]3/2} \right\} $

(2)

通过数值计算，可以直接获得归一化的球体体积V_s/r₀³与归一化的球体半径a/r₀的关系，结果如图 3所示。通过对数据做拟合，可以得到一个关于归一化的球体体积V_s/r₀³与归一化的球体半径a/r₀的简单闭合方程

$ \frac{a}{{{r_0}}} = 0.665{\left( {\frac{{{V_{\rm{s}}}}}{{r_0^3}}} \right)^{0.325}} $

(3)

图 3

图 3 载体桩球体半径与体积的关系 Fig. 3 Relation between the radius of the sphere of carrier pile and the volume

如果已知填料的体积用量，便可以通过式(3)计算出载体扩大头球体的半径，从而利用图 2(b)所示的球孔扩张理论来计算挤密效应。

图 2(b)给出了球孔扩张理论模型，初始的球孔半径为a₀(对于本文问题a₀=0)，位于无限大的土体中，在无穷远处受到的初始均匀应力场应力为σ_h0。初始圆孔在均匀内压σ_a的作用下扩张到半径为a的圆孔。孔洞周围的土体服从修正剑桥模型。随着球孔内压的不断增大，孔周围土体要发生屈服，孔洞附近土体进入塑性状态，而离孔洞比较远的地方土体处于弹性状态，继续扩张会使得孔周土体进入临界状态。这样，在孔洞周围会形成3块区域：孔洞附近的临界状态区和塑性区以及远离孔洞的弹性区。此外，推导解析解基于两个基本假设：扩孔过程，孔周土体假定是完全排水的，因为在柱锤夯扩的过程中，土中的水要被排出，土体被挤密；土体弹性关系服从各向同性的胡克定律，塑性关系服从修正剑桥模型。

2 球孔扩张基本控制方程及求解

对于基于修正剑桥模型的弹塑性球孔扩张问题，孔洞周围将会出现3个区域：临界状态区域、塑性区域以及弹性区。在临界状态区域和塑形区中，应力必须满足平衡方程和修正剑桥模型屈服准则，而在弹性区中，应力需要满足平衡方程和胡克定律。此外，孔扩张的运动学方程则由土体的排水条件、应变-位移关系以及边界条件来控制。在弹塑性分界面(弹性区和塑形区的分界面)上，应力和位移必须保持连续。基于这样的基本原理有以下的推导：

2.1 临界状态区和塑性区

1) 平衡方程

根据Collins等^[19]的研究，平均应力p′和偏应力q可以写成

$ {{p^\prime } = \frac{{\sigma _{\rm{r}}^\prime + \sigma _\theta ^\prime }}{2}} $

(4)

$ {q = (\sigma _{\rm{r}}^\prime - \sigma _\theta ^\prime )} $

(5)

式中：σ′_r，σ′_θ为球坐标下的径向和环向有效应力分量。对于球孔扩张问题另外一个应力分量σ′_φ=σ′_θ。因此，采用应力分量p′和q表达的应力平衡方程可以写成

$ \frac{{\partial {p^\prime }}}{{\partial r}} + \frac{2}{3}\frac{{\partial q}}{{\partial r}} + 2\frac{q}{r} + \frac{{\partial {u_{\rm{w}}}}}{{\partial r}} = 0 $

(6)

式中:u_w为土体的孔隙水压。

2) 本构方程

采用Collins等^[19]的方法来定义两个应变分量

$ {{\varepsilon _{\rm{p}}} = {\varepsilon _{\rm{r}}} + 2{\varepsilon _\theta }} $

(7)

$ {{\varepsilon _{\rm{q}}} = \frac{2}{3}({\varepsilon _{\rm{r}}} - {\varepsilon _\theta })} $

(8)

根据应变-位移关系，可以知道球坐标下的应变分量可以写成

$ {{\varepsilon _{\rm{r}}} = - \frac{{\partial {u_{\rm{r}}}}}{{\partial r}}} $

(9)

$ {{\varepsilon _\theta } = - \frac{{{u_{\rm{r}}}}}{r}} $

(10)

式中：u_r为土体的径向位移；r为土颗粒的径向坐标。

定义符号(°)表示土体颗粒的物质导数，可以表达为

$ {(^ \circ }) = ( \cdot ) + v\partial ({\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} )/\partial r $

(11)

式中：v为土颗粒的径向速度；(·)为局部的时间导数。

土颗粒的径向和环向应变量可以表示为

$ {{{\mathop \varepsilon \limits^{\ \ 。} }_{\rm{r}}} = - \frac{{\partial v}}{{\partial r}}} $

(12)

$ {{{\mathop \varepsilon \limits^{{\ \ 。}} }_\theta } = - \frac{v}{r}} $

(13)

将式(12)和式(13)带入到式(7)和式(8)中，可以得到

$ {{{\mathop \varepsilon \limits^{{\ \ 。}} }_{\rm{p}}} = - \frac{{\partial v}}{{\partial r}} - k\frac{v}{r}} $

(14)

$ {{{\mathop \varepsilon \limits^{{\ \ 。}} }_{\rm{q}}} = \frac{2}{3}\left( { - \frac{{\partial v}}{{\partial r}} + \frac{v}{r}} \right)} $

(15)

进一步地，弹性本构关系可以写成

$ \left[ {\begin{array}{*{20}{c}} {\mathop \varepsilon \limits^{{\ \ 。}} {\kern 1pt} _{\rm{p}}^{\rm{e}}}\\ {\mathop \varepsilon \limits^{{\ \ 。}} {\kern 1pt} _{\rm{d}}^{\rm{e}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\frac{1}{K}}&0\\ 0&{\frac{k}{{2(k + 1)G}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathop p\limits^{{\ \ 。}} {{\kern 1pt} ^\prime }}\\ {\mathop q\limits^{{\ \ 。}} } \end{array}} \right] $

(16)

式中：$ G = \frac{{3\left( {1 - 2v'} \right)\upsilon p'}}{{2\left( {1 + v'} \right)\kappa }} $，是土体的剪切模量；$ K = \frac{{2\left( {1 + v'} \right)}}{{3\left( {1 - 2v'} \right)}}G $，是土体的体积模量，v′是土体的有效泊松比。

根据Wood^[20]的研究，剑桥模型弹塑性的本构关系为

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{{\mathop \varepsilon \limits^{{\ \ 。}} }_{\rm{p}}}}\\ {{{\mathop \varepsilon \limits^{{\ \ 。}} }_{\rm{q}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\frac{1}{K}}&0\\ 0&{\frac{1}{{3G}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\mathop p\limits^{{\ \ 。}} {{\kern 1pt} ^\prime }}\\ {\mathop q\limits^{{\ \ 。}} } \end{array}} \right] + \frac{{2\frac{q}{{{p^\prime }}}(\lambda - \kappa )}}{{(1 + e)\left( {{M^2} + \frac{{{q^2}}}{{{p^{\prime 2}}}}} \right){p^\prime }}} \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ \begin{array}{l} \frac{{\left( {{M^2} - \frac{{{q^2}}}{{{p^{\prime 2}}}}} \right)}}{{2\frac{q}{{{p^\prime }}}}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 1\\ 1\ \ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{{2\frac{q}{{{p^\prime }}}}}{{{M^2} - \frac{{{q^2}}}{{{p^{\prime 2}}}}}} \end{array} \right]\left[ {\begin{array}{*{20}{l}} {\mathop p\limits^{{\ \ 。}} {{\kern 1pt} ^\prime }}\\ {\mathop q\limits^{{\ \ 。}} } \end{array}} \right] \end{array} $

(17)

将式(17)中的两个应变分量采用式(14)和式(15)中的表达式代替，可以得到

$ {{A_{\rm{p}}}\mathop p\limits^{\ \ 。} {{\kern 1pt} ^\prime } + {A_{\rm{q}}}\mathop q\limits^{\ \ 。} + \frac{{\partial v}}{{\partial r}} + 2\frac{v}{r} = 0} $

(18)

$ {{B_{\rm{p}}}\mathop p\limits^{\ \ 。} {{\kern 1pt} ^\prime } + {B_{\rm{q}}}\mathop q\limits^{\ \ 。} + \frac{2}{3}\left( {\frac{{\partial v}}{{\partial r}} - \frac{v}{r}} \right) = 0} $

(19)

其中：

$ {A_{\rm{p}}} = \frac{1}{K} + \frac{{2\frac{q}{{}}(\lambda - \kappa )}}{{(1 + e)\left( {{M^2} + \frac{{{q^2}}}{{{p^{\prime 2}}}}} \right){p^\prime }}}\left( {\frac{{{M^2} - \frac{{{q^2}}}{{{p^{\prime 2}}}}}}{{2\frac{q}{{{p^\prime }}}}}} \right) $

(20)

$ {{A_{\rm{q}}} = \frac{{2\frac{q}{{{p^\prime }}}(\lambda - \kappa )}}{{(1 + e)\left( {{M^2} + \frac{{{q^2}}}{{{p^{\prime 2}}}}} \right){p^\prime }}}} $

(21)

$ {{B_{\rm{p}}} = {A_{\rm{q}}}} $

(22)

$ {{B_{\rm{q}}} = \frac{1}{{3G}} + \frac{{2\frac{q}{{{p^\prime }}}(\lambda - \kappa )}}{{(1 + e)\left( {{M^2} + \frac{{{q^2}}}{{{p^{\prime 2}}}}} \right){p^\prime }}}\left( {\frac{{2\frac{q}{{{p^\prime }}}}}{{{M^2} - \frac{{{q^2}}}{{{p^{\prime 2}}}}}}} \right)} $

(23)

3) 一致性条件

此外，屈服面需要满足如下的一致性条件

$ {{\rm{d}}f = \frac{{\partial f}}{{\partial {p^\prime }}}\mathop p\limits^{} {{\kern 1pt} ^\prime } + \frac{{\partial f}}{{\partial q}}\mathop q\limits^{\ \ 。} + \frac{{\partial f}}{{\partial {p_{\rm{c}}}^\prime }}\mathop p\limits^{\ \ 。} {{\kern 1pt} _{\rm{c}}}^\prime = 0} $

(24)

$ {f = {q^2} - \left[ {{M^2}{p^\prime }\left( {{p^\prime }_{\rm{c}} - {p^\prime }} \right)} \right]} $

(25)

$ {\frac{{\partial f}}{{\partial {p^\prime }}} = {M^2}\left( {{p^\prime }_{\rm{c}} - 2{p^\prime }} \right)} $

(26)

$ {\frac{{\partial f}}{{\partial q}} = 2q} $

(27)

$ {\frac{{\partial f}}{{\partial {p^\prime }_{\rm{c}}}} = {M^2}{p^\prime }} $

(28)

其中：p′_c表示卡各向同性压缩下的屈服压力。

4) 连续性条件

在球孔扩张过程中，土体需要满足质量守恒条件

$ \frac{{\mathop \upsilon \limits^{。} }}{\upsilon } = \frac{{\partial v}}{{\partial r}} + 2\frac{v}{r} $

(29)

其中：υ为体积比。

5) 排水条件

由于采用了完全排水的条件，也即在扩孔过程中土体孔隙水压保持不变，因此有

$ {\mathop u\limits^{\ \ 。}{ _{\rm{w}}}} = 0 $

(30)

上面5个条件即构成了球孔扩张的控制方程组，采用相似变换技术求解该方程组。采用如下变换方式

$ \eta = \frac{r}{{ct}} $

(31)

$ \left( {\widetilde {{p^\prime }},\tilde q,{\rm{ }}\widetilde {{p_{\rm{c}}}^\prime },\widetilde {{u_{\rm{w}}}},\tilde v} \right) = \left( {\frac{{{p^\prime }}}{{{p^\prime }_0}},\frac{q}{{{p^\prime }}},\frac{{{p_{\rm{c}}}^\prime }}{{{p^\prime }_0}},\frac{{{u_{\rm{w}}}}}{{{p^\prime }_0}},\frac{v}{c}} \right) $

(32)

式中：c为弹塑性边界的扩张速度；t为时间。将上述两个变换方程带入到前面5个条件构成的方程组中，可以得到方程组

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} 1&{\frac{2}{3}}&0&1&0&0\\ {{A_{\rm{p}}}(\tilde v - \eta )}&{{A_{\rm{q}}}(\tilde v - \eta )}&0&0&0&1\\ {{B_{\rm{p}}}(\tilde v - \eta )}&{{B_{\rm{q}}}(\tilde v - \eta )}&0&0&0&{\frac{2}{3}}\\ {{M^2}({{\widetilde {{p^\prime }}}_{\rm{c}}} - 2{{\tilde p}^\prime })}&{2\tilde q}&{{M^2}{{\tilde p}^\prime }}&0&0&0\\ 0&0&0&0&{\frac{1}{\upsilon }(\tilde v - \eta )}&{ - 1}\\ 0&0&0&1&0&0 \end{array}} \right] \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} {{\rm{d}}\widetilde {{p^\prime }}/{\rm{d}}\eta }\\ {{\rm{d}}\tilde q/{\rm{d}}\eta }\\ {{\rm{d}}{{\widetilde {{p^\prime }}}_{\rm{c}}}/{\rm{d}}\eta }\\ {{\rm{d}}{{\tilde u}_{\rm{w}}}/{\rm{d}}\eta }\\ {{\rm{d}}\upsilon /{\rm{d}}\eta }\\ {{\rm{d}}\tilde v/{\rm{d}}\eta } \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - 2\tilde q/\eta }\\ { - 2\tilde v/\eta }\\ {\frac{2}{3}\tilde v/\eta }\\ 0\\ {2\tilde v/\eta }\\ 0 \end{array}} \right] \end{array} $

(33)

为了求解方程式(33)，需要给出初始条件和边界条件，可以利用孔洞和弹塑性边界处的信息来确定。

在弹塑性边界处

$ {\eta = 1} $

(34)

$ {\widetilde {{p^\prime }}(\eta = 1) = 1} $

(35)

$ {\tilde q(\eta = 1) = M\sqrt {\left( {{p^\prime }_{{\rm{c0}}}/{p^\prime }_0 - 1} \right)} } $

(36)

$ {\widetilde {{u_{\rm{w}}}}(\eta = 1) = {u_{{\rm{w}}0}}/{p^\prime }_0} $

(37)

$ {u(\eta = 1) = {v_0}} $

(38)

$ {G(\eta = 1) = {G_0} = \frac{{3\left( {1 - 2{v^\prime }} \right){v_0}{p^\prime }_0}}{{2\left( {1 + {v^\prime }} \right)\kappa }}} $

(39)

$ {\tilde v(\eta = 1) = \frac{{q(\eta = 1)}}{{2G(\eta = 1)}}} $

(40)

在孔口处

$ \eta = \tilde v $

(41)

结合控制方程和上述边界条件编写Matlab微分方程组数值求解程序，可以获得临界状态区和塑性区的半解析解。

2.2 弹性区

在弹性区，本构关系服从胡克定律，因此有

$ \Delta {\sigma ^\prime }_{\rm{r}} = - \frac{E}{{(1 + \nu )(1 - 2\nu )}}\left[ {(1 - \nu )\frac{{\partial {u_{\rm{r}}}}}{{\partial r}} + 2\nu \frac{{{u_{\rm{r}}}}}{r}} \right] $

(42)

$ \Delta {\sigma ^\prime }_\theta = - \frac{E}{{(1 + \nu )(1 - 2\nu )}}\left[ {\nu \frac{{\partial {u_{\rm{r}}}}}{{\partial r}} + \frac{{{u_{\rm{r}}}}}{r}} \right] $

(43)

同时，在弹性区域，应力需要满足平衡方程，将平衡方程写成关于σ_r和σ_θ的形式

$ \frac{{\partial {\sigma ^\prime }_{\rm{r}}}}{{\partial r}} + \frac{2}{r}({\sigma ^\prime }_{\rm{r}} - {\sigma ^\prime }_\theta ) = 0 $

(44)

将式(42)和式(43)代入式(44)，可以得到

$ \frac{{{\partial ^2}{u_{\rm{r}}}}}{{\partial {r^2}}} + \frac{2}{r}\frac{{\partial {u_{\rm{r}}}}}{{\partial r}} - \frac{{2{u_{\rm{r}}}}}{r} = 0 $

(45)

采用如下的变换

$ \eta = \frac{r}{{ct}},\widetilde {{u_{\rm{r}}}} = \frac{{{u_{\rm{r}}}}}{{ct}} $

(46)

式(45)可以简化为常微分方程

$ \frac{{{{\rm{d}}^2}\widetilde {{u_r}}}}{{{\rm{d}}{\eta ^2}}} + \frac{2}{\eta }\frac{{{\rm{d}}\widetilde {{u_r}}}}{{{\rm{d}}\eta }} - \frac{2}{{{\eta ^2}}}\widetilde {{u_r}} = 0 $

(47)

式(47)存在闭合解析解

$ {\bar u_{\rm{r}}} = \frac{A}{{{\eta ^2}}} + \frac{B}{3}\eta $

(48)

式中：A和B为常数，当η→∞时，u_r必须衰减为0，因此，常数B必须为0。式(48)可以简化为

$ \widetilde {{u_{\rm{r}}}} = \frac{A}{{{\eta ^2}}} $

(49)

此外

$ \begin{array}{*{20}{l}} {\tilde q(\eta = 1) = [\widetilde {{\sigma ^\prime }_{\rm{r}}}(\eta = 1) - \widetilde {{\sigma ^\prime }_\theta }(\eta = 1)] = }\\ {M\sqrt {{p^\prime }_{{\rm{c0}}}/{p^\prime }_0 - 1} } \end{array} $

(50)

因此，球坐标下的应力分量便可以通过式(42)和式(43)获得

$ \widetilde {{\sigma ^\prime }_{\rm{r}}} = \frac{{{\sigma ^\prime }_{\rm{r}}}}{{{p^\prime }_0}} + \frac{{{\sigma ^\prime }_0}}{{{p^\prime }_0}} = \frac{{2M\sqrt {{p^\prime }_{{\rm{c0}}}/{p^\prime }_0 - 1} }}{3}\frac{1}{{{\eta ^3}}} + \widetilde {{\sigma ^\prime }_0} $

(51)

$ \widetilde {{\sigma ^\prime }_\theta } = \frac{{{\sigma ^\prime }_\theta }}{{{p^\prime }_0}} + \frac{{{\sigma ^\prime }_0}}{{{p^\prime }_0}} = - \frac{{M\sqrt {{p^\prime }_{{\rm{c0}}}/{p^\prime }_0 - 1} }}{3}\frac{1}{{{\eta ^3}}} + \widetilde {{\sigma ^\prime }_0} $

(52)

$ \widetilde {{u_{\rm{r}}}} = \frac{{M\sqrt {{p^\prime }_{{\rm{c0}}}/{p^\prime }_{\rm{0}} - 1} }}{{3\sqrt 3 G/{p^\prime }_{\rm{0}}}}\frac{1}{{{\eta ^2}}} $

(53)

式中：η≥1。

3 球孔扩张引起的土体挤密效应讨论

载体桩夯扩过程中土体被挤密，土体的强度和剪切模量会发生变化。根据修正剑桥模型，土体的强度和剪切模量可以表示为

$ {{s_{\rm{u}}} = \frac{{\sqrt 3 }}{3}M{p^\prime }\left( {\frac{R}{2}} \right)1 - \kappa /\lambda } $

(54)

$ {G = [3(1 - 2{v^\prime })\upsilon {p^\prime }]/[2(1 + {v^\prime })\kappa ]} $

(55)

对于排水孔扩张问题，土体的体积比υ=1+e随着扩孔的过程发生变化，此外，土体的平均有效应力p′也发生变化，这两个参数的变化导致了土体强度和剪切模量的变化。而在扩孔过程中，其余的本构参数κ、λ、M、ν′、R可以认为保持一个常量，不发生变化。因此，可以定义

$ {\delta = \frac{{{s_{\rm{u}}}}}{{{s_{{\rm{u0}}}}}} = \frac{{{p^\prime }}}{{{p^\prime }_0}}} $

(56)

$ {\tilde \omega = \frac{G}{{{G_0}}} = \frac{{{p^\prime }\upsilon }}{{{p^\prime }_0{\upsilon _0}}}} $

(57)

式中：δ表示球孔扩张后土体强度和初始土体强度之比；s_u0表示土体的初始强度；ω表示球孔扩张后土体剪切模量和初始土体剪切模量之比；G₀表示土体的初始剪切模量。

4 参数分析

图 4和图 5分别给出了不同M条件下球孔扩张后孔周土体强度和剪切模量变化系数随着孔周径向距离的变化规律。M的取值范围从0.5到2，其余参数如图 4、图 5所示。可以看出，随着M的不断增大，孔周附近(r/a＜2)的土体强度和剪切模量系数不断增加，挤密效应不断增强。然而，当径向距离r/a > 2时，随着M的不断增大，土体强度和刚度系数不断减小，挤密效应不断减弱。在2＜r/a＜3的位置处，土体强度和剪切模量系数小于1，说明该处土体的强度和剪切模量在扩孔后发生退化，土体产生剪胀，而并没有被挤密，但该部分区域范围很小。此外，在2＜r/a＜3的位置处，还可以发现，当M > 1时，土体的强度和剪切模量系数随着径向距离的增加，出现先减小后增加的现象。

图 4

图 4 不同M条件下强度变化系数随着径向距离的变化 Fig. 4 Variation of strength coefficient with the radial distance for different M

图 5

图 5 不同M条件下剪切模量变化系数随着径向距离的变化 Fig. 5 Variation of shear modulus coefficient with the radial distance for different M

图 6和图 7分别给出了不同超固结比R条件下球孔扩张后孔周土体强度和剪切模量变化系数随着孔周径向距离的变化规律。R的取值范围从1到10，涵盖了正常固结土，轻微超固结土和重度超固结土，其余参数如图 6、图 7所示。可以看出，随着R的不断增大，孔周的土体强度和剪切模量系数不断增加，挤密效应不断增强，而且在M取值为1的条件下，孔周土体的强度和剪切模量没有出现软化的现象。对于重度超固结土，球孔扩张后土体强度可以提高7倍左右，剪切模量可以提高6.5倍左右。此外，土体超固结比R较大时，土体强度和剪切模量系数随着径向距离的变化率会出现先减小，再增加的转折现象。而对于正常固结土或轻微超固结土，变化率随着径向距离的增加基本逐渐减小。

图 6

图 6 不同R条件下强度变化系数随着径向距离的变化 Fig. 6 Variation of strength coefficient with the radial distance for different R

图 7

图 7 不同R条件下剪切模量变化系数随着径向距离的变化 Fig. 7 Variation of shear modulus coefficient with the radial distance for different R

图 8和图 9给出了不同λ条件下强度和刚度变化系数随着径向距离的变化规律。λ的取值范围从0.05到0.25，反映了土体不同的可压缩特性。可以发现，随着λ的不断增加，孔周的土体强度和剪切模量系数不断减小，挤密效应不断减弱。说明土体可压缩性越小，土体强度和刚度的提高就越大。此外，不同λ条件下，球孔扩张后土体的强度和剪切模量也并未出现软化现象。

图 8

图 8 不同λ条件下强度变化系数随着径向距离的变化 Fig. 8 Variation of strength coefficient with the radial distance for different λ

图 9

图 9 不同λ条件下剪切模量变化系数随着径向距离的变化 Fig. 9 Variation of shear modulus coefficient with the radial distance for different λ

图 10和图 11给出了不同κ条件下强度和刚度变化系数随着径向距离的变化规律。κ的取值范围从0.01到0.05，反映了土体不同的回弹特性。可以发现，随着κ的不断增加，孔周的土体强度和剪切模量系数轻微增加，说明参数κ对孔周的土体强度和剪切模量系数的影响很小。另外，在3＜r/a＜4的位置处，球孔扩张后土体的强度和剪切模量出现软化现象。

图 10

图 10 不同κ条件下强度变化系数随着径向距离的变化 Fig. 10 Variation of strength coefficient with the radial distance for different κ

图 11

图 11 不同κ条件下剪切模量变化系数随着径向距离的变化 Fig. 11 Variation of shear modulus coefficient with the radial distance for different κ

图 12和图 13给出了不同泊松比ν′条件下强度和刚度变化系数随着径向距离的变化规律。泊松比ν′的取值范围从0.15到0.49。可以看出，泊松比ν′对孔周的土体强度和剪切模量系数的影响也很小，随着泊松比ν′的增加，土体强度和剪切模量系数只是略微减小。在3＜r/a＜4的位置处，球孔扩张后土体的强度和剪切模量出现软化现象。

图 12

图 12 不同ν′条件下强度变化系数随着径向距离的变化 Fig. 12 Variation of strength coefficient with the radialdistance for different v′

图 13

图 13 不同v′条件下剪切模量变化系数随着径向距离的变化 Fig. 13 Variation of shear modulus coefficient with the radial distance for different v′

图 14给出了不同本构参数条件下挤密扰动区半径大小的变化规律。可以看出，随着参数M从0.5增加到2，挤密扰动区的半径从5倍的球孔半径减小到4.2倍的球孔半径。对于超固结比R，当土体属于正常固结土(R=1)或者轻微超固结土(R＜1.5)时，挤密区的半径随着超固结比的不断增加，迅速减小；当超固结比1.5＜R＜5时，超固结比对挤密区的半径影响较小，挤密区的半径随着超固结比的增加轻微减小。当超固结比大于5时，可以发现，挤密区的半径近似等于球孔半径，说明对于重度超固结土，球孔扩张产生的挤密区可以忽略不计。此外，参数λ和κ的增加会导致挤密区不断减小。对于泊松比ν′，可以看出, 当泊松比ν′＜0.4时，随着泊松比ν′的不断增加，挤密区轻微减小，当泊松比ν′>0.4时，挤密区的半径急剧减小。

图 14

图 14 挤密区半径随着不同参数的变化规律 Fig. 14 Variation of the radius of the compaction zone with different parameters

5 结论

探讨了软黏土中载体桩桩端挤密效应，建立了基于球孔扩张的载体桩桩端夯扩挤密力学模型，分析了桩端挤密效应，可以得到如下的结论：

1) 给出了载体桩桩端夯扩填料体积用量与球孔扩张理论模型中球孔半径关系的闭合解析表达式。

2) 利用相似变换技术求解球孔扩张偏微分控制方程，获得了修正剑桥模型中球孔扩张力学响应的数值解。

3) 采用球孔扩张前后土体强度比和剪切模量比两个物理力学指标来定量分析载体桩桩端夯扩挤密后土体力学性质变化，获得了土体强度比和剪切模量比随着土体参数的变化规律。采用塑性区半径来分析挤密扰动区的大小，获得了挤密扰动区的大小随着土体参数的变化关系。

4) 提出了考虑载体桩桩端挤密效应的理论计算方法，该方法可以较为准确地计算载体桩夯扩过程中周围土体强度和剪切模量的变化以及挤密扰动区的大小，可为建立考虑载体桩桩端挤密效应的承载力计算方法提供理论基础。

参考文献

[1]	王继忠. 载体桩技术的诞生与发展[J]. 建筑结构, 2008, 38(4): 118-119. WANG J Z. Birth and development of carrier pile technology[J]. Building Structure, 2008, 38(4): 118-119. (in Chinese)
[2]	王建安, 张慧海, 陈琳, 等. 载体桩在陕南粉质黏土地区的适用性研究[J]. 岩土工程技术, 2019, 33(2): 109-114. WANG J A, ZHANG H H, CHEN L, et al. Study on bearing capacity of cast-in-place pile in silty clay area of southern Shaanxi Province[J]. Geotechnical Engineering Technique, 2019, 33(2): 109-114. (in Chinese)
[3]	于长杰. 载体桩对软土地基沉降量控制效果的试验分析[J]. 四川建筑, 2018, 38(3): 101-103. YU C J. Experimental analysis of control effect of carrier pile on settlement of soft soil foundation[J]. Sichuan Architecture, 2018, 38(3): 101-103. (in Chinese)
[4]	周斌.高速铁路载体桩复合地基沉降变形规律研究[D].成都: 西南交通大学, 2011. ZHOU B. Study on settlement regularity of pile with bearing base composite foundation in high-speed railway[D].Chengdu: Southwest Jiaotong University, 2011. http://cdmd.cnki.com.cn/Article/CDMD-10613-1011233723.htm
[5]	罗浩.高速铁路载体桩复合地基沉降特性现场试验及施工工艺研究[D].成都: 西南交通大学, 2010. LUO H. Filed research on settlement characters and construction technology of ram-compaction piles with composite bearing base in high-speed railway[D]. Chengdu: Southwest Jiaotong University, 2010. http://cdmd.cnki.com.cn/Article/CDMD-10613-1011233570.htm
[6]	李建强, 刘洪滨, 张家尊, 等. 载体桩复合地基性状数值模拟研究[J]. 工业建筑, 2018, 48(6): 107-110. LI J Q, LIU H B, ZHANG J Z, et al. Research on the numerical simulation of the carrier pile composite foundation[J]. Industrial Construction, 2018, 48(6): 107-110. (in Chinese)
[7]	张培成, 孙玉文, 张殿树, 等. 基于饱和软土地基大承载力载体桩试验研究[J]. 河北水利电力学院学报, 2020(1): 19-23. ZHANG P C, SUN Y W, ZHANG D S, et al. Experimental study on piles with ram-compacted bearing sphere for large saturated soft bearing capacity based on soil foundation[J]. Journal of Hebei University of Water Resources and Electric Engineering, 2020(1): 19-23.
[8]	仇凯斌, 介玉新, 李广信, 等. 载体桩承载力性状有限元分析[J]. 工程勘察, 2009, 37(Sup1): 15-20. QIU K B, JIE Y X, LI G X, et al. Finite element analysis of bearing capacity of carrier pile[J]. Geotechnical Investigation & Surveying, 2009, 37(Sup1): 15-20. (in Chinese)
[9]	张续萱, 魏殿兴. 锚定板挡土结构的实测土压力及其分析研究[J]. 岩土工程学报, 1994, 16(2): 73-79. ZHANG X X, WEI D X. Measurement and analysis of earth pressure in anchor slab retaining structures[J]. Chinese Journal of Geotechnical Engineering, 1994, 16(2): 73-79. (in Chinese)
[10]	ZHOU H, KONG G Q, LIU H L, et al. Similarity solution for cavity expansion in thermoplastic soil[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2018, 42(2): 274-294. DOI:10.1002/nag.2724
[11]	ZHOU H, KONG G Q, LIU H L. A semi-analytical solution for cylindrical cavity expansion in elastic-perfectly plastic soil under biaxialin situ stress field[J]. Géotechnique, 2016, 66(7): 584-595. DOI:10.1680/jgeot.15.P.115
[12]	ZHOU H, LIU H L, ZHA Y H, et al. A general semi-analytical solution for consolidation around an expanded cylindrical and spherical cavity in modified Cam Clay[J]. Computers and Geotechnics, 2017, 91: 71-81. DOI:10.1016/j.compgeo.2017.07.005
[13]	ZHOU H, LIU H L, KONG G Q, et al. Analytical solution of undrained cylindrical cavity expansion in saturated soil under anisotropic initial stress[J]. Computers and Geotechnics, 2014, 55: 232-239. DOI:10.1016/j.compgeo.2013.09.011
[14]	ZHOU H, KONG G Q, LI P, et al. Flat cavity expansion:theoretical model and application to the interpretation of the flat dilatometer test[J]. Journal of Engineering Mechanics, 2016, 142(1): 04015058. DOI:10.1061/(ASCE)EM.1943-7889.0000957
[15]	YU H S. Cavity expansion methods in geomechanics[M]. Springer Science & Business Media, 2013.
[16]	CHEN S L, ABOUSLEIMAN Y N. Exact undrained elasto-plastic solution for cylindrical cavity expansion in modified Cam Clay soil[J]. Geotechnique, 2012, 62(5): 447-456. DOI:10.1680/geot.11.P.027
[17]	CHEN S, ABOUSLEIMAN Y N. Exact drained solution for cylindrical cavity expansion in modified Cam Clay soil[J]. Geotechnique, 2013, 63(6): 510-517. DOI:10.1680/geot.11.P.088
[18]	LI L, LI J P, SUN D A. Anisotropically elasto-plastic solution to undrained cylindrical cavity expansion in K0-consolidated clay[J]. Computers and Geotechnics, 2016, 73: 83-90. DOI:10.1016/j.compgeo.2015.11.022
[19]	COLLINS I F, STIMPSON J R. Similarity solutions for drained a undrained cavity expansions in soils[J]. Geotechnique, 1994, 44(1): 21-34.
[20]	WOOD D M. Soil behaviour and critical state soil mechanics[M]. Cambridge: Cambridge University Press, 1991.

X