超固结黏土中不排水柱孔扩张弹塑性解及应用

吕桂阳，和西良，梁汝鸣，付艳斌，刘 洋，葛序尧，祝人杰，张 波

(1.中建八局第二建设有限公司设计研究院，山东，济南 250014；2.深圳大学土木与交通工程学院，深圳 518060；3.深圳大学滨海城市韧性基础设施教育部重点实验室，深圳 518060)

圆柱形孔穴扩张问题的解决方案由于形式简单，力学原理明确，广泛应用于解释一些现场实测数据以及实际工程现象，如孔压静力触探技术(CPTU)试验贯入过程和旁压试验[1-2]，分析桩基工程中静压沉桩引起的应力和孔压变化[3-4]以及桩基承载[5]等问题。

过去半个世纪，数值求解能力有限，本构模型较为简单，使得柱孔扩张问题大多把土体看作理想弹塑性材料，求解柱孔周围土体单元的响应规律。例如采用Mohr-Coulomb 准则[6-7]和Tresca准则[8]。VESIĆ[6]考虑了土体压缩特性，采用Mohr-Coulomb 准则，求解了柱孔扩张问题解析解；CARTER 等[8]运用小应变理论，结合Mohr-Coulomb准则以及非相关联流动法则，求得柱孔扩张问题弹塑性解。

后来，一些学者开始应用各向同性硬化本构模型，例如原始剑桥模型和修正剑桥模型推导孔穴扩张问题(例如文献[9 - 14])。周航等[9]结合保角变换与复变函数理论推导了非圆形孔的弹性解，并分析了正方形孔的力学特性。CAO 等[10]采用修正剑桥模型分析了不排水孔穴扩张理论，给出了孔穴扩张周围有效应力的闭合解，表达了超静孔隙水压力和总压力。张亚国等[12]运用坐标变换的方法，结合复变函数理论解答了半无限体中考虑斜坡体自重圆孔扩张问题的解析解。

最近，在求解柱孔扩张问题中，一些较为复杂的本构模型得到了应用[15-19]。例如，SU 和YANG[16]采用剪胀性与状态相关的边界面模型来研究各参数对砂土中圆柱形扩张影响。武孝天和徐永福[17]采用边界面模型概念，假定加载一旦开始土体立即产生塑性变形，模拟了砂/黏土柱孔扩张过程中土体应力场和应变场的变化规律。周茗如等[18]以塑性力学与大变形理论为基础，运用结构性黄土修正剑桥模型推导了圆孔扩张问题的理论解并应用到劈裂注浆压力理论计算中。

上述研究对孔穴扩张理论做出了重要贡献，对孔穴扩张周围土体响应有了更好的认识。由于土体材料存在不确定性，静止土压力参数K0对柱孔扩张的影响研究较少、超固结黏土中预测桩基工程中静压沉桩和CPTU 实测孔压数据方法不多。

为进一步完善柱孔扩张理论，本文采用CASM屈服面以及非相关流动法则，结合空间应力转化方法，建立超固结黏土中考虑K0影响的柱孔不排水扩张模型，分析超固结黏土中考虑K0影响的柱孔周围土体单元的响应规律，并结合桩基工程中静压沉桩以及CPTU 现场孔压实测数据验证模型的实用性以及正确性，同时为类似实际工程应用提供理论支撑。

1 柱孔扩张模型分析

图1 展示了在瞬时的柱孔扩张模型中，土体受到初始应力 σ0和初始孔隙压力u0作用。通过施加均匀分布的扩张压力，当压力足够大时，孔壁上的土将首先屈服。随着扩张压力继续增大，围绕柱孔就会形成一个塑性区，其中临界状态区将在半径rc中存在。塑性区以外的土体被认为处于弹性状态区。当扩张压力从 σ0增加到扩孔压力σa时，初始柱孔从半径a0扩张到当前半径ax；对应弹塑性边界半径从rp0移动到rp。用柱形极坐标系统(r, θ,ν)来描述柱孔扩张特性，总应力分量分别为 σr、 σθ、 σν，有效应力分量为。用总应力和相应有效应力之差计算孔隙水压力

总平均应力和偏应力分别用p和q表示，根据柱面问题的分析定义如下：

在弹性和塑性区，距离柱形扩张中心为r位置处的径向和环向应力满足平衡微分方程：

假设：受压应力和应变为正。土体弹性阶段服从胡克定律，土体塑性阶段服从CASM 模型准则。

2 建立超固结黏土柱孔扩张模型

2.1 弹性区分析

假定在弹性区域r≥rp为小变形，柱孔径向和环向应变增量 dεr和 dεθ在r处可以写成：

式中， dur为位移增量。

弹性阶段服从胡克定律表示为：

式中， dεν为体应变增量，与比体积的关系为：

在根据不排水条件总体积应变率为0：

在弹性阶段，比体积和平均有效应力不变。

瞬时弹性模量E和剪切模量G也应该不变且等于对应的初始值E0和G0。

现结合以上公式，径向和环向有效应力可以表示为：

通过上述方程可求得弹性区任意一点处土体单元的平均有效应力p′。在弹性区，总平均有效应力p没有变化，平均有效应力p′是恒定的常数，可以相应的推出没有超静孔隙水压力的变化，u=u0。

2.2 本构模型简介

在临界状态下，土体表现为具有恒定体积和应力的摩擦性流体材料；对于给定的土体来说，无论其应力路径如何变化，临界状态线CSL 是唯一的，定义为：

式中：q和p′分别为偏应力和平均有效应力；M为p′-q应力空间中临界状态线的斜率；μ=1+e0为比体积，e0为孔隙比；λ、κ和Γ 为临界状态常量(见图2)。

图2 ln p′-μ空间状态参数图Fig.2 ln p′-μ Space state parameter diagram

为了解决柱孔扩张过程中屈服面高估“干面”破坏应力，本文采用了YU[20]基于临界状态理论，引入应力状态系数n和参考状态参数r∗(ξR=(λ-κ)lnr∗)的边界状态曲面，称为CASM。能够同时在排水和不排水的加载条件下，较好模拟黏土和砂土整体受力特性。描述CASM 的状态边界面为：

状态参数ξ：在相同的平均有效应力下，当前状态与临界状态之间比体积差值，ξ=μ+λlnp′-Γ，也被BEEN 和JEFFERIES[21]表达为颗粒材料在一定围压和密度范围内的力学性能。q/p′为应力比；n为应力状态系数，是一种新的材料常数，通常在1.0～5.0；ξR=(λ-κ)lnr∗为参考状态参数；r*为间距比，定义为。同时，原始剑桥模型屈服面对应式(14)中n=1和r∗=2.7183的情况。修正剑桥模型屈服面对应选择r*=2和一个合适的n=1.5 ∼2.0可以很好模拟“湿面”屈服面以及克服“干面”高估破坏应力的情况。根据YU 和COLLINS[22]在用临界状态模型考虑不同超固结比分析时取n=2.0和r∗=2.718，本文取其值进行分析。此时CASM 和MCC 本构模型的屈服面形状如图3 所示。

图3 CASM 与修正剑桥模型屈服面的比较Fig.3 Comparison of yield surface between CASM and modified Cambridge model

对于平均有效应力和偏应力以及状态参数下，屈服函数的形式可以表示为：

同时，为了能够反映强超固结黏土的应变软化现象，本文采用ROWE[23]的应力-剪胀关系表示为：

塑性势能可以通过应力-剪胀关系(式(16))的积分得到：

对于任意给定的应力状态p′、q，通过求解式(17)可以很容易地确定尺寸参数 β。

塑性流动法则采用的是非相关联的，硬化法则采用了经典的各向同性体积塑性应变硬化，表示为：

2.3 弹塑性分析

在p′-q空间中，当q增加到屈服面时，土体开始进入屈服，此时总应变增量为：

式中，上标e 和 p分别为弹性和塑性应变。弹性应变分量的增量，满足胡克定律与应力增量的关系。

扩张过程中，孔周土体刚度E是弹性模量，定义为2G(1+ν′)；G是剪切模量，在WOOD[24]中用平均有效应力p′和比体积μ表示：

式中：ν′为有效泊松比；κ为ν-lnp′平面上卸载再加载曲线的斜率。在扩张过程中ν′为常数，不排水条件下体积不变，比体积μ增量为0。因此，剪切模量G与平均有效应力p′成正比。

塑性应变分量的增量，假设CASM 模型中的塑性流动是不相关联的与应力增量的关系：

式中： dL为比例因子； η为应力比，定义为q/p′。根据YU[20]的研究，扩张过程中塑性体积应变为：

通过将式(20)、式(23)代入式(22)， dL求得：

其中：

将式(6)、式(28)代入到式(19)弹塑性本构方程中得：

或者求逆矩阵：

柱孔扩张过程中塑性区变形由大应变理论来描述，因此将环向应变增量定义为：

对于不排水的条件，体积变化为0，即dεv=0，因此可以给出径向应变增量。

将式(33)代入式(32)，建立一阶常微分方程，求解塑性区柱孔扩张有效径向、环向应力与扩张半径之间的关系：

式(35)可以计算塑性区位置r对应的有效径向和环向应力。

2.4 边界条件求解

随着柱孔扩张压力不断增大，在扩张过程中土体单元从初始位置r0移动到一个新的位置r。在CHEN 和ABOUSLEIMAN[14]方法中，土体单元弹性部分的总位移计算如下：

由于不排水的变形必然是等体积的，所以径向半径位移与扩张半径位移关系:

通过设置r和rxp均等于rp，这代表了当前弹塑性边界上土体单元的位置，推导结果：

在弹塑性边界上由应力连续性条件可得：

联立屈服面方程式(15)求得屈服时的偏应力qy：

考虑到静止土压力参数K0的影响：

联立式(48)～式(50)可以得到弹塑性边界上的有效应力：

2.5 超静孔隙水压力

本文超静孔压计算方法为：通过有效应力原理，将总应力等于有效应力与孔隙水压力之和的关系，即σr=+u和σθ=+u。代入平衡方程式(4)，从扩孔周围土体单元的弹塑性边界半径rp处，积分到扩孔后土体单元对应半径rx处，便可计算出超静孔隙水压力：

3 理论模型验证与分析

通过与CHEN 和ABOUSLEIMAN[14]、李镜培等[15]及武孝天和徐永福[17]三种方法进行对比，验证和分析本文理论模型方法的正确性与实用性，黏土基本物理力学参数选取CHEN 和ABOUSLEIMAN[14]方法进行参数分析，见表1。

表1 黏土基本物理力学参数Table 1 Basic physical and mechanical parameters of clay

3.1 应力路径

对比图4(a)、图4(b)和图4(c)可见，MCC 模型在p-q平面内初始和最终屈服面为椭圆形，而CASM 模型通过应力状态系数n=2.0和参考状态参数r∗=2.718旋转屈服面，来模拟“湿面”屈服面以及克服“干面”屈服面高估屈服应力以及破坏应力的现象。

在图4(a)中可知，当土体OCR≤2 时，本文方法破坏时偏应力略高于MCC 模型，但在Roscoe物态边界面以内，且本文方法与物态边界面较为接近。有效应力路径在进入屈服面时表现为土体的应变硬化现象。在图4(b)和图4(c)中对比Hvorslev面可以得出，当土体OCR>2 时本文方法更接近土体在强超固结时发生屈服应力和破坏应力。有效应力路径在进入屈服面时表现为土体的应变软化现象。

当静止土压力参数K0≠1 时，土体单元初始应力状态是在屈服面内；当静止土压力参数K0=1时，初始应力状态为各项同性，不存在偏应力。

3.2 应力分布

为进一步验证本文方法的正确性与适用性，通过对比CHEN 和ABOUSLEIMAN[14]、李镜培等[15]及武孝天和徐永福[17]方法对应力分布分析如下：图5 表示不同超固结比OCR 和不同静止土压力参数K0下，黏土柱孔扩张周围归一化量纲的径向应力、环向应力、和竖向应力随径向位移r/a的分布规律。图5 中横坐标采用对数坐标。

图5 应力分布Fig.5 Stress distribution

对比图5 可知，本文方法与修正剑桥模方法[14-15]和CSUH 模型方法[17]弹塑性边界的位置不同，且随着超固结比的增大，弹塑性边界的位置相差越大，这是由于本文方法采用了CASM 模型，修正了重度超固结黏土中屈服应力被高估的现象；在CSUH 模型方法[17]中采用了加载一旦开始，土体立即产生塑性应变，因此扩孔周围无弹塑性边界位置。在塑性区，有效应力变化的快慢与硬化法则有关。

当土体OCR≤2 时，本文方法应力分布与CHEN 和ABOUSLEIMAN[14]方法较为吻合，表明本文方法能够较好地模拟土体单元“湿面”屈服应力和破坏应力。

当土体OCR>2 时，本文方法应力分布与李镜培等[15]方法更接近，这是由于SMP-MCC 模型引入应力罗德角，其屈服面大小和形状随应力状态变化而改变，与本文采用的CASM 模型较为吻合，克服了“干面”高估屈服应力和破坏应力的情况。同时发现，当OCR=3 时，本文方法与李镜培等[15]、武孝天和徐永福[17]两种方法最终的破坏应力值基本相同，进一步验证了本文方法的正确性与适用性。

3.3 超静孔隙水压力与径向位移的关系

对比图6(a)、图6(b)和图6(c)可知，本文方法和修正剑桥模型方法[14-15]和CSUH 模型方法[17]的超静孔隙水压力之间的异同。

图6 超静孔隙水压力分布Fig.6 Excess pore water pressure distribution

当土体OCR≤2 时，本文方法超静孔隙水压力分布与CHEN 和ABOUSLEIMAN[14]方法较为吻合；当土体OCR>2 时，本文方法超静孔隙水压力分布与李镜培等[15]方法更接近；且本文方法比CSUH 模型方法[17]超静孔隙水压力分布较低，这是由于CSUH 模型方法假定加载一旦开始，土体立即产生塑性变形，土体无弹性变形。同时在OCR=10时，超静孔隙水压力分布出现负孔压，这是由于重度超固结黏土中径向和环向有效应力急剧变化引起，反映了重度超固结黏土的剪胀特性。

3.4 剪切模量与径向位移的关系

图7 表示柱孔扩张半径a/a0=2 时，不同超固结比OCR 和不同静止土压力参数K0下扩孔周围土体单元的剪切模量变化规律，在弹性区和临界状态区土体的剪切模量基本为一个定值，在塑性区会产生急剧的变化。当土体为轻微超固结和重度超固结时，本文方法剪切模量分布规律与CHEN和ABOUSLEIMAN[14]、李镜培等[15]两种方法基本一致；当土体为重度超固结时，本文方法剪切模量分布低于CHEN 和ABOUSLEIMAN[14]和李镜培等[15]两种方法，这是由于本文方法采用了CASM模型，修正了屈服应力和破坏应力被高估的现象。

图7 剪切模量G/G0 与径向位移r/a 的关系Fig.7 Relationship between shear modulus G/G0 and radial displacement r/a

3.5 分析参数K0 与扩孔关系的影响

图8 和图9 展示了不同参数K0下归一化量纲扩孔压力[σa-(1+2K0)σ1/3]/su和超静孔隙水压力ua/su与归一化量纲扩孔半径a/a0之间的关系。当K0=1 时与CHEN 等[19]方法基本相同。

图8 不同参数K0 下扩孔压力与扩孔半径的关系Fig.8 Variation of normalised internal cavity pressure with cavity radius with different parameters K0

图9 不同参数K0 下扩孔超静孔隙水压力与扩孔半径的关系Fig.9 Variation of normalised excess pore pressure with cavity radius with different parameters K0

当0＜K0＜1 时，归一化量纲极限扩孔压力和极限超静孔隙水压力随着K0的增加而增加；当K0＞1 时，归一化量纲极限扩孔压力和极限超静孔隙水压力随着K0的增加而减小；当K0=1 时，归一化极限扩孔压力和极限超静孔隙水压力达到最大值，这是由于处于土体处于三向等压状态，受力稳定，一旦土体处于三向不等压状态，就会有偏应力存在，先于三向等压状态发生破坏。

同时，当K0由1.0 变化到2.5 时，归一化极限扩孔压力由5.4 变化到3.2；当K0由1.0 变化到2.5 时，归一化极限超静孔隙水压力由2.4 变化到1.0。当K0较大时(如K0=2.5)在扩孔孔壁周围会产生负的超静孔隙水压力，这是由于竖向应力与水平应力比值较大的影响。如果忽略K0的影响，得到的极限扩孔压力和极限超静孔隙水压力偏大，高估了土体实际应力状态，对实际工程不利。

综上，本文建立的超固结黏土柱孔扩张模型，克服了“干面”高估屈服应力和破坏应力的情况，较好地模拟土体单元“湿面”屈服应力和破坏应力，考虑了参数K0的影响，真实反映了柱孔周围土体单元在扩孔过程中应力位移响应规律。

4 工程应用

应用1：加拿大魁北克市西部，ROY 等[25]进行了一系列静压桩试验并测得了土体初始参数以及桩身周围孔隙水压力的变化。试验场地的工程地质情况为：粉质黏土具有中等至低塑性，并受到中度的先期固结作用，OCR 在2.0～3.0 范围内。原位有效应力从2.0 m 深处的16 kPa 线性增加到6.0 m 深处的36 kPa。现场试验测量的不排水抗剪强度su的从2.0 m 深处的12 kPa 线性增加到6.0 m深处的32 kPa。土壤刚度指数Ir，定义为剪切模量G与不排水抗剪强度su之比约为150，有效内摩擦角φ′为30°。

通过对比CAO 等[10]基于修正剑桥模型方法和RANDOLPH 等[26]经验方法分析如下：从图10中可知，经验方法不能考虑黏土超固结、静止土压力参数K0以及打桩和贯入过程中土体参数变化的影响，高估了实测值。

应用2：深圳地铁12 号线海上田园东站基坑试验区进行了一系列CPTU 试验如图11、图12和图13 所示，并通过室内土工试验测得相应的土体初始参数见表2。

表2 黏土基本物理力学参数Table 2 Basic physical and mechanical parameters of clay

图11 CPTU 试验设备Fig.11 CPTU test equipment

图12 CPTU 探头室内饱和Fig.12 CPTU probe indoor saturation test

图13 现场试验Fig.13 Field test

通过试验得知场地附近工程地质情况：地表以下土层分布情况如图14 所示。

图14 土层分布Fig.14 Soil layer distribution

本文选取③3黏土层CPTU 现场实测数据分析孔隙水压力分布规律。该黏土层受到先期固结作用，超固结比OCR 约在2.3～3.3，原位竖向应力从224 kPa 线性增加至310 kPa，根据CPTU 现场实测数据，确定黏土不排水抗剪强度su从70 m深处的12 kPa 线性增加至19 m 深处的90 kPa，土体刚度指数Ir约为50。

通过对比CAO 等[10]方法和李镜培等方法[15]分析如下：从图15 中可知，基于SMP 准侧改进的修正剑桥模型方法考虑了打桩过程中土体参数变化的影响，但在超固结黏土中还存在高估的现象。

图15 CPTU 孔压预测值与实测值对比Fig.15 Comparison of CPTU pore pressure predicted valueand measured value

综上，本文方法充分考虑了土体的应力历史、土体参数变化及三维特性，比较接近地预测了现场实测数据。

5 结论

本文考虑了土体静止土压力参数K0和重度超固结比对柱孔扩张问题的影响，基于CASM 屈服准则，采用经典的应力剪胀关系，运用非相关联流动法则推导了黏土中考虑超固结与K0影响的柱孔扩张问题的弹塑性刚度矩阵。并通过数值方法分析了孔壁周围土体应力和超静孔隙水压力分布规律，应用到打桩和CPTU 贯入过程中扩孔力学机制研究中，得出以下结论：

(1) 建立了黏土中考虑超固结与K0影响的柱孔不排水扩张问题求解模型。

(2) 通过对比物态边界面中Hvorslev 线表明，本文方法在中度与重度超固结黏土中，克服了屈服应力和破坏应力被高估的现象。

(3) 静止土压力参数K0≠1 时，柱孔不排水扩张周围土体处于三向不等压状态，会有偏应力存在，先于三向等压状态发生破坏。如果忽略K0的影响，高估了土体实际应力状态，对实际工程不利。

(4) 对比现场实测数据得出，本文方法能够预测中度与重度超固结土中打桩与CPTU 贯入过程中扩孔周围孔隙水压力的变化规律。