饱和结构性黄土不排水柱孔扩张问题弹塑性解

李镜培，周 攀，李 亮，谢 峰，崔纪飞

（1.同济大学岩土及地下工程教育部重点实验室，上海 200092；2.同济大学土木工程学院，上海 200092）

西部大开发战略的实施使黄土分布最广的西北地区岩土工程问题研究变得十分紧迫［1-2］。尤其是在流域低阶地区，因地质环境改变、地下水位上升而形成了饱和黄土，此类黄土虽然有着与软土相近的某些特征，但是具有独特的结构性，因此对饱和黄土地区的岩土工程问题研究十分重要。对于许多复杂岩土工程问题，如原位土工试验（旁压测试、静力触探测试等）、桩基工程、隧道和地下工程以及井壁失稳等，圆孔扩张理论是既简单又适用的分析工具［3-7］。

在岩土体中的圆孔扩张或收缩实质上是一个一维边值问题，其求解往往需要一个合理描述岩土材料应力-应变关系的本构模型。理想弹塑性本构模型（Mohr-Coulomb 理论、Tresca 准则和统一强度理论等）和临界状态应变硬化或软化本构模型（剑桥模型、修正剑桥模型以及众多用于特殊土体的改进修正剑桥模型）是目前广泛应用于反映土体力学特征的2 类本构模型。圆孔扩张理论应用初期，许多学者采用摩尔-库伦准则求解柱孔和球孔扩张问题。Vesic［8］基于土体塑性屈服后的体积改变推导出了柱孔扩张和球孔扩张问题的近似解；Carter等［9］采用非关联流动法则描述土体塑性流动，求得了小应变下孔压与孔扩张关系的显式解，同时求得了有限变形下孔压与圆孔扩张关系的数值解；邹金锋等［10］同样采用非关联流动法则并考虑体变守恒和大变形理论推导了砂性土中极限扩孔压力表达式。同时，其他理想弹塑性模型适用性也较为广泛，丰富了岩土工程中圆孔扩张问题的解答。饶平平等［11］假定土体服从Lade-Duncan 屈服准则，采用大变形和非关联流动法则对圆柱孔扩张问题进行了解答；曹黎娟等［12］、綦春明等［13］在弹性区采用小变形理论，在塑性区采用统一强度准则与大变形理论，推导出了扩孔问题的塑性区半径、极限扩孔压力以及超静孔压的解析式。

理想弹塑性模型无法考虑土体应力历史与超固结比的影响，因此临界状态模型逐渐受到广大学者的青睐而广泛应用于柱孔扩张问题的研究。由于临界状态模型自身的复杂性，因此在使用该类模型时或多或少会进行适当简化，其中包括假定塑性区内剪切模量恒定与偏应力线性分布［14］、应力路径线性简化［15-17］以及平均塑性应变值估算［18］等，以便求解柱孔和球孔扩张问题解析解。为了反映土体各向异性、强度三维力学特性对扩张问题计算结果的影响，李镜培等［19-20］采用各向异性修正剑桥模型和基于空间滑动面（SMP）准则改进的修正剑桥模型分析了不排水和排水条件下柱孔扩张机制。为了得到柱孔扩张问题的精确解，Chen等［21-22］采用修正剑桥模型，在未对平均有效应力和偏应力做任何简化的前提下求得了柱孔不排水与排水扩张半解析解。随后，郑金辉等［23］在此求解框架基础上采用考虑砂土颗粒破碎的本构模型，分析了砂土圆孔扩张问题。

对于天然的岩土材料，特别是广泛分布于我国西北地区的黄土，其结构性对土体强度和变形特性的影响是土力学最本质的一个因素。然而，以上研究大部分是针对于求解饱和黏土、软土和砂土等的柱孔扩张问题，现有解析解中对考虑土体结构性的扩孔机理分析仍存在不足，尤其是黄土地区的圆孔扩张问题解答较少。因此，基于Chen等［21］的分析框架，在遵循2个有效应力不变的严格定义前提下，通过采用邓国华等［24-25］提出的黄土结构性修正剑桥模型来描述土体塑性屈服后的材料强度和变形特性，为结构性黄土的圆孔排水扩张边值问题求解推导出精确的半解析解，最终得到孔周围土体精确的有效径向、切向和竖向应力以及超静孔压的分布规律。

1 不排水柱孔扩张问题力学模型

1.1 基本假设与模型建立

图1为无限空间内的饱和黄土柱孔扩张力学模型。土体受到的初始水平应力与初始竖向应力分别为和，且=，其中K0为静止侧压力系数。假设柱孔初始半径为a0，初始扩孔压力为σ0，扩孔稳定后柱孔半径为a，同时扩孔压力增加到σa。当扩张完成后，孔周任意一点土体从原始位置rx0移动到当前位置rx，其位移用Ux表示，孔周塑性区半径从原始位置rp0移动到当前位置rp，其位移用Up表示。当孔压增大到一定值后，孔周土体应力处于3 种状态：孔壁附近为临界状态区（a≤r≤rf，其中rf为破坏区半径），应力不再改变，应变不断发展；向孔外延伸一定厚度的土体为塑性区（rf＜r ＜rp），随着孔压增大应力应变仍在改变；距离孔壁较远处为弹性区（r≥rp），随着距孔壁距离越大，土体应力状态越接近初始应力值。

图1 结构性黄土中的柱孔扩张力学模型Fig.1 Mechanical model of cylindrical cavity expansion in structured loess

柱孔扩张问题的求解前提是基于以下假设：

（1）土体为均质、连续且各向不等压的饱和黄土。

（2）孔隙水不可压缩，应力和应变以压为正。

（3）平面应变和不排水条件下柱孔扩张，因此扩孔过程中土体体积不变。

（4）在弹性区土体服从线弹性和小变形假定，在塑性区采用大变形理论和结构性黄土修正剑桥模型描述土体的力学特性。

孔周任意土体微单元总应力满足以下平衡微分方程：

式中：σr和σθ分别为孔周土体总径向应力和切向应力。式（1）还可表示为有效应力形式，如下所示：

1.2 结构性黄土本构模型

在孔周弹性区采用胡克定律描述土体应力状态与变形特征，应力-应变关系的增量形式如下所示：

式中：E为弹性模量，E=2G（1+ν），其中ν为泊松比，G为剪切模量；εij，e为应变分量；为应力分量；为主应力分量；δij为克罗内克函数。G的计算式如下所示：

式中：υ为土体比容；p′为有效平均应力；κ为回弹曲线斜率。

在孔周塑性区采用邓国华等［24-26］、邵生俊等［27］提出的结构性黄土修正剑桥模型，如下所示：

式中：a、b、A为试验拟合参数；φ为土体内摩擦角。mη的表达式如下所示［24］：

式中：下标i、r、s分别表示原状土、重塑土和饱和土；ηi、ηr、ηs分别为原状土、重塑土和饱和土在一般应力路径下的应力比；pi、pr、ps分别为原状土、重塑土和饱和土在一定体应变条件下的等向压应力或球应力；qi、qr、qs分别为原状土、重塑土和饱和土在某一剪切变形条件下的剪应力。

考虑应变硬化规律后，式（5）可进一步表示为

式中：λ为土体压缩曲线斜率；为初始有效平均应力；e0为土体初始孔隙比；εv，p为塑性体应变。当综合结构性参数C=0 时，式（12）退化为修正剑桥模型屈服方程表达式。

2 柱孔扩张问题半解析解

2.1 弹性区分析

由Chen等［21-22］分析可知，根据小变形理论、胡克定律与孔周土体微单元的应力平衡方程可推得孔周弹性区应力分量、位移与超静孔隙水压力的解析式，如下所示：

式中：为弹塑性界面上的有效径向应力；G0为初始剪切模量。G0的计算式如下所示：

式中：υ0为土体初始比容。

2.2 弹塑性界面分析

对于初始应力各向异性的结构性黄土，孔周初始有效径向应力、切向应力和竖向应力之间的关系如下所示：

式中：σh0、σv0分别为总初始水平、竖向应力；u0为初始孔隙水压力。式（19）可转化为

由弹性区应力解析式可知，扩孔过程中孔周土体的有效平均应力p′保持不变，则根据应力连续条件弹塑性界面处有效平均应力可表示为

将（，q0）代入式（5）屈服方程中，可推得初始加载点所对应的当前屈服面与平均正应力p轴的交点p′A，p′A的计算式如下所示：

式中：q0为初始有效偏应力。

将弹塑性界面上点（，qp）代入式（5）屈服方程中，并联立式（21）和式（22），可解得弹塑性界面上偏应力

式中：qp为弹塑性界面上的偏应力；R为结构性黄土超固结比。R的计算式如下所示：

弹塑性界面上的有效径向应力、切向应力、竖向应力和偏应力qp分别表示为如下关系式：

结合式（25）～（27），可解得弹塑性界面上应力分量

设刚进入塑性状态的微单元位置为rxp，其初始位置为rx0，由式（16）可得该土体单元的径向位移

对于不排水柱孔扩张，土体单元体积不变，因此孔周任意一点的当前径向坐标rx和初始径向坐标rx0与圆孔当前半径a和初始半径a0之间的关系可表示为

结合式（31）和式（32）可推得柱孔不排水过程中孔周任意一点土体进入塑性状态时径向位置

令rx=rxp=rp，则由式（33）可得柱孔不排水扩张过程中孔周塑性区半径rp与当前扩张半径a之间的关系为

2.3 塑性区分析

孔周土体塑性屈服后总应变增量由弹性体应变增量和塑性体应变增量两部分组成，如下所示：

式中：ε为应变。土体屈服后弹性体应变增量任然采用胡克定律计算，由结构性黄土剑桥模型和相关联流动法则，可通过下式计算塑性体应变增量：

式中：Λ为塑性因子，表示塑性应变增量大小的一个标量；f为屈服函数。

塑性体应变增量dεv，p与主应力方向上的塑性应变增量dεr，p、dεθ，p和dεz，p满足

联立式（36）和式（37），由塑性一致性条件可将结构性黄土修正剑桥模型的塑性因子表示为

有效平均应力增量dp′与偏应力增量dq可表示为

联立式（27）～（35），孔周塑性区应力-应变关系式采用对称矩阵形式表示，如下所示：

结合孔周弹性区应力-应变关系式（3）和塑性区应力-应变关系式（47），能够得到用矩阵表示的弹塑性本构方程，如下所示：

式中：εr、εθ、εz分别为径向、切向和竖向应变。

式（52）转化为如下形式：

基于大变形理论，孔周塑性区r处的径向应变增量dεr和切向应变增量dεθ的拉格朗日描述形式为

柱孔扩张属于典型的平面应变问题，因此竖向应变增量亦为零，即

根据不排水条件可知，柱孔扩张过程中土体体应变增量为零，即

联立式（64）～（67），可得

将式（68）代入式（53），可将柱孔扩张过程的本构方程转化为如下形式的一阶非线性常微分方程组：

注意到b11、b12、b21、b22、b31和b32均为关于3 个应力分量、和的显式函数，因此能够通过3 个应力分量在弹塑性界面上的初始值对式（69）～（71）进行数值求解。

孔周塑性区内任意位置处超静孔隙水压力Δurx可通过对式（2）积分得出，如下所示：

式中：为当前计算点处的有效径向应力。

3 解析解验证与讨论

上述解析式和求解程序为结构性黄土不排水柱孔扩张问题提供了严格的理论框架，然而对于控制方程式（69）～（71）无法直接得到显式解析解，必须通过数值方法求解。首先，通过Matlab 数值软件对Chen等［21］的算例参数进行编程求解，用于验证本方法的正确性；然后，研究初始结构性与超固结比对柱孔扩张的影响与孔周土体应力和超孔隙压力分布。为了与已有计算结果进行对比以研究土体结构性在扩孔过程中对土体受力和变形的影响规律，采用文献［21］中算例参数，并且分别对超固结比为1.2、3.0和10.0 的土体假定综合结构性参数C为0、50、100 kPa，计算参数如表1 所示，其中υcs表示土体正应力为1 kPa时的土体比容。

将表1中参数代入第2节公式中，通过Matlab数值方法，可得到孔周土体应力和超静孔压分布、孔壁处扩孔压力、超静孔压以及应力路径，计算结果如图2～4所示。图2～4中，虚线表示Chen等［21］基于修正剑桥模型（未考虑结构性影响）的计算值，实线为本方法在相同土体参数中考虑结构性影响的解。当C=0 时，本方法的解退化为Chen 等［21］的解，间接验证本方法的正确性。当C≠0 时，由本方法解与Chen 等［21］解的差异性来分析土体结构性对沉桩扩孔力学机理的影响。

表1 计算参数Tab.1 Calculated parameters（M=1.2, λ=0.15, κ=0.03,ν=0.278, υcs=2.740）

图2为不同先期固结压力与结构性参数条件下孔壁总应力σa、孔壁超静孔压Δu（a）（通过除以不排水强度su进行归一化处理）与扩孔半径之间的关系。

图2 归一化的扩孔压力与超静孔压随扩孔半径的变化Fig.2 Variation of normalized internal cavity pressure and excess pore pressure with cavity radius

如图2 所示，在扩孔过程中随着扩孔半径的增加扩孔压力与孔壁超静孔压均随之增大，当a/a0介于1 和2 之间时，σa/su与Δu（a）/su均变化剧烈，直到a/a0超过2时两者逐渐趋于稳定，表明土体已达到临界状态，此时应力不再变化而剪应变持续发展。还可以看出，归一化的扩孔压力和孔壁超静孔压随超固结比增加而减小。由此可知，无论是否考虑土体结构性影响，扩孔压力、孔壁超静孔压随扩孔半径的变化规律都与Chen等［21］的分析结果一致，当结构性参数C=0 时，两者的计算结果完全相同，因此Chen 等［21］的计算结果为本解析解的一种特殊情况（C=0）。还可知，当考虑土体结构性时，扩孔压力σa/su、孔壁超静孔压Δu（a）/su均大于未考虑结构性时的计算结果，但是当土体受到的先期固结压力较大时，扩孔压力σa/su、孔壁超静孔压Δu（a）/su的变化幅值受结构性影响减弱。

图3为不同超固结比和结构性参数下孔周有效径向、切向和竖向应力分布规律。如图3所示，无论是否考虑结构性影响，扩孔稳定后孔周土体分为弹性、塑性和破坏3种状态，靠近孔壁区域的土体达到了临界状态而破坏，应力不再变化，塑性区土体应力变化急剧，弹性区土体随着距孔壁越远应力状态越接近初始应力状态。由图3a 可知，当忽略土体结构性影响时，本方法计算结果与Chen等［21］的计算结果一致，表明采用修正剑桥模型分析柱孔扩张是本方法的特例（C=0）。由图3b、3c可知，当考虑土体结构性影响时，相比于Chen 等［21］的计算结果，本方法计算的塑性区半径与破坏区半径均减小，并且有效径向、切向和竖向应力也均减小，但有效径向应力相比于有效切向应力减小的幅度要小，即（σ′r-σ′θ）增大，表明扩孔过程中考虑结构性参数影响提高了土体抵抗剪切变形的能力，还可知修正剑桥模型的计算结果低估了扩孔过程结构性土体的强度。

图3 不同超固结比结构性黄土的有效径向、切向与竖向应力分布（a/a0=2）Fig.3 Effective radial,tangential and vertical stress distributions for different over-consolidated structured loesses（a/a0=2）

图4为不同超固结比和结构性参数下孔周土体有效平均应力p′、偏应力q与超静孔压Δu的分布规律。

图4 不同超固结比结构性黄土的有效平均应力、偏应力与超静孔压分布（a/a0=2）Fig.4 Effective mean stress,deviator stress and excess pore pressure distributions for different overconsolidated structured loesses（a/a0=2）

由图4可以看出，无论是否考虑结构性影响，在破坏区，有效平均应力与偏应力均不变，但随着距孔壁越近超静孔压呈线性急剧增大。在塑性区，随着距孔壁越近，有效平均应力变化较小，对于轻度超固结土（R=1.2），p′逐渐减小且q逐渐增加到临界值，超静孔压急剧增大，但对于中度超固结土（R=3.0）和重度超固结土（R=10.0），p′逐渐增大且q先增大到最大值再逐渐减小至临界值，并且超固结比越大，此种变化趋势越明显。在重度超固结土中超静孔压先减小至最小值再增大，表明在重度超固结土体中出现了应变软化特性伴随着剪胀效应，但是考虑结构性影响能够减小土体的剪胀效应，如图4c所示。在弹性区超静孔压与有效平均应力保持不变，随着距孔壁越近，偏应力逐渐增加。由图4a 可知，当忽略土体结构性影响时，本方法计算结果与Chen等［21］计算结果一致，表明采用修正剑桥模型分析柱孔扩张是本研究的特例（C=0），这与上述分析一致。由图4b、c可知，当考虑土体结构性影响时，相比于Chen等［21］的计算结果，本方法计算的超静孔压与偏应力均增大，有效平均应力减小，但是超固结比增大，三者的变化幅值均减小，表明在重度超固结下的结构性土体中，结构性参数对土体的应力与变形影响减弱，这与图3中分析结果一致。

图5为不同超固结比与结构性参数影响下孔壁处黄土微单元在p′-q平面上的有效应力路径（ESP），图中E、F和G分别表示初始加载点、初始屈服点和临界状态点。

图5 孔壁处不同超固结比结构性黄土在p′-q平面上的应力路径（a/a0=2）Fig.5 Stress path in p′-q plane around cavity wall for different over-consolidated structured loesses（a/a0=2）

由图5a可知，对于轻度超固结土（R=1.2），采用的黄土结构性修正剑桥模型的初始屈服面与Chen等［21］采用的修正剑桥模型的初始屈服面相同，并且土体应力路径也一致，进一步说明当不考虑结构性影响时本方法解析解退化为Chen 等［21］的解。对于中度超固结土（R=3.0）和重度超固结土（R=10.0），从图5b、c 可以看出，考虑土体结构性参数时，两者初始屈服面均向左平移，并且屈服面形状发生变化，临界状态线也向左平移相同距离且斜率不变。无论是否考虑土体结构性影响，在3 种超固结土体扩孔时孔壁土体均先经历弹性阶段（p′不变，q增加），在应力点运动到初始屈服面并产生塑性体应变后则进入塑性阶段（p′增大或减小，q增加），直到应力点达到临界状态线后进入破坏阶段（p′不变，q不变）。

4 结论

（1）考虑了土体结构性，推导出的半解析解能够更加真实地反映扩孔过程中土体应力状态。与采用修正剑桥模型的计算结果相比，本方法中结构性黄土的初始屈服面向左平移并且椭圆尺寸发生改变，平移量与尺寸改变量取决于结构性参数的取值。当忽略结构性影响时，本方法计算结果与Chen等［21］提出方法的计算结果一致。

（2）孔周土体分为3种状态：外侧土体处于弹性状态，并且距孔壁越远土体应力状态越接近扩孔前的初始状态；中间土体处于塑性状态，有效应力分量与超静孔压急剧变化；内侧土体处于临界状态，土体有效应力保持不变，但超静孔压仍变化较大。

（3）土体结构性与超固结比对孔周应力与变形影响显著，弹塑性半径与临界区半径均随着超固结比和综合结构性参数增大而减小。对于轻度与中度超固结土，当产生塑性体应变后土体出现应变硬化现象。当超固结比较大时出现负超静孔压，表明出现塑性软化现象，但考虑结构性能够减小土体软化程度。

作者贡献声明:

李镜培：给出论文研究方法，稿件的审核和监督。

周 攀：给出论文研究方法，软件实施，研究结果验证，稿件初稿写作和返修。

李 亮：稿件审阅和修改。

谢 峰：稿件审阅和修改。

崔纪飞：稿件审阅和修改。