0 引言

高压旋喷桩技术是公路工程中常用的地基加固手段^[1]，具有施工便捷、桩型灵活等优点。高压旋喷桩施工过程中，随着高压浆液持续注入土层，会引起较大的超孔隙水压力，挤压周边一定区域土体发生变形^[2]，在一些地质环境脆弱的城市地区，可能对附近的建(构)筑物和地下管线产生不利影响，甚至引发工程事故，具有重要的研究意义。

目前，相关研究主要基于圆孔扩张理论展开，采用多种半理论半经验分析方法，求解圆柱或球形空腔在无限/半无限空间中膨胀引起的应变场^[3-4]。也有学者结合有限元等数值模拟手段，在孔壁设置均布应力^[5]或对塑性区施加体积膨胀^[6]，间接反映高压浆液对周围土体挤压作用。这些方法一般可获得较好的分析结果，但其主要计算参数如塑性区半径、体积膨胀系数、孔壁应力等无法预先确定，需要依靠工程实测数据结果反演，或基于经验公式近似取值。另外，由于存在网格畸变问题等，难以模拟浆液冲击土体的大变形过程。为此，亟需一种可直接考虑浆液与土体相互作用的计算手段。

光滑粒子流体动力学方法(Smoothed Particle Hydrodynamics，SPH)是一种基于拉格朗日框架的无网格粒子法，将问题域离散为携带质量、密度、应力等场变量的“粒子”^[7]。由于SPH方法的无网格特性，在处理大应变、动边界和自由表面的问题中具有显著优势，已被应用于模拟滑坡、泥石流等岩土材料的大变形问题^[8]。同时，SPH方法可用于求解多相耦合问题，已发展均相模型、分相模型等求解方法^[9-10]。近年来，土体冲击问题成为SPH方法的研究热点，相关学者在SPH框架下模拟了高压射流开挖沟槽^[11]、冲蚀土体^[12]的过程，探究了土体应力分布和变形情况，在材料本构选择、水-土两相流模拟等方面取得了积极进展。

在此基础上，本研究尝试进一步在构建工程尺度的SPH模型，模拟高压旋喷施工引起的饱和土体变形，为该类问题的研究提供一种新的思路。通过构建基本控制方程和材料本构方程的SPH形式，采用混合物理论解决两相介质耦合问题，并选择合适的数值技术，编制了相应的分析程序。作为初步探究，本研究主要关注高压浆液对饱和地基的挤压效应，在后续研究中将进一步考虑对土体的扰动破坏作用。

1 SPH数值方法 1.1 SPH基本原理

SPH本质上是一种核密度估计方法，其基本思想为“核近似”和“粒子近似”，如图 1所示^[13]。核近似是将宏观连续介质动力学方程表示为适合粒子计算的SPH形式，将场函数f(x)的微分运算转化为光滑函数W的微分运算：

(1)

式中，▽为位置矢量x处的梯度；Ω为x的积分域；l为光滑长度，代表光滑函数W支持域的大小。

粒子近似是将连续域离散为有限数量的粒子，每个粒子上的场变量可由影响半径内所有近邻粒子的场变量加权求和得到，由此进一步将式(1)变换为离散形式：

(2)

式中，m和ρ分别为粒子j的质量和密度；m_j/ρ_j为与粒子j关联的体积ΔV_j；N为影响域内粒子总数。

1.2 SPH形式的基本控制方程和材料本构方程

本研究中涉及孔隙水、水泥浆和土体骨架3种材料，流体、固体两相介质。

采用弱可压缩形式的Navier-Stokes方程作为流体基本控制方程。假设流体连续且与流体有关的物理场可微，并引入Monaghan型人工黏度项Π_ij防止流体粒子发生非物理性穿透，得到质量守恒和动量守恒方程的SPH形式：

(3)

式中，α，β为坐标方向；ρ，v^α，x^α，τ^αβ，t分别代表密度、速度、空间坐标、总应力张量和时间；f^α为作用在单位质量流体上的体力。为避免额外求解压力泊松方程，采用Monaghan状态方程将不可压缩流体转化为弱可压流体^[9]。

土体同样遵守质量守恒和动量定理，因此在SPH中常采用相似的基本控制方程。同时，为解决粒子之间吸引力过大而聚集成块的问题，可在动量守恒方程中增加Monaghan提出的人工应力项f_abⁿ(R_a^αβ+ R_b^αβ)为粒子提供短程排斥力。

对于孔隙水，可以近似地采用牛顿流体本构^[14]；而对于水泥浆，其流变特征与水灰比等因素相关，一般将水灰比W/C为0.8~1.0的水泥浆考虑为宾汉流体^[11]。

针对土体粒子采用D-P屈服准则的弹塑性本构模型。在大变形问题中，可引入Jaumann应力变化率反映刚性转动对力学行为的影响^[15]，最终得到土体粒子应力率方程的SPH形式：

(4)

式中，为旋转率张量；为偏应变张量；s_a为偏应力张量。

1.3 两相介质模拟方法

将饱和土体视为孔隙水与土颗粒均匀混合的两相介质，基于混合物理论对两相耦合流进行描述，其原理如图 2所示。在生成初始模型时，分别生成土粒子和水粒子，并对其各自的材料属性进行赋值。在每个时间步内，水、土粒子先按照各自的控制方程计算，进行粒子搜索、密度求解、应力求解等。然后，求解两相之间的作用力，包括渗流力和孔隙水压力，并施加于近邻粒子。对于浆-土关系，在本模型中视为与水-土关系相似，仅考虑水泥浆和孔隙水在材料参数方面的区别。

以上标波浪线符号“~”代表材料在单相流中的取值，根据混合物理论可作如下假设^[16]：

(5)

式中，I为二阶单位张量；p_f为水的压强，也即孔隙水压力，SPH表达式为：

(6)

渗流力可由Den Adel半经验公式计算。假设受到高压浆液冲击时土中的渗流为紊流状态，采用Forchheimer方程计算水力梯度。

1.4 SPH方法的程序架构

本研究在开源程序PersianSPH基础上，对前述Wendland光滑函数、密度变化率、应力变化率、两相耦合方程等进行编译，选择合适的数值技术，基于Linux平台，采用C++语言编写SPH分析程序，计算流程基本框架如图 3所示。

2 高压浆液冲击土体模拟示例 2.1 模型建立及参数设立

SPH模型初始粒子分布如图 4所示。土体部分总尺寸为5 m×10 m(宽×高)。考虑地下水位在地面以下0.5 m，设置干土和饱和土两区域。其中干土位于H=9.5~10 m范围内，为图中黑色粒子指示区域，任一位置只有土粒子；饱和土位于H=0~9.5 m范围内，任一位置同时存在土粒子及水粒子，如左侧示意图所示。

模型右侧浅色粒子区域指示水泥浆，总尺寸为0.5 m×4.5 m，作用范围为H=5.0~9.5 m。计算开始后，水泥浆以较高速度冲击土体，模拟高压浆液引起土体变形的过程。可根据伯努利方程建立浆液粒子初始速度与浆液压力之间的关系：

(7)

式中，v_c为浆液粒子初始速度；P_i为浆液压力；P_c为浆液自重产生的压力；ρ_c为水泥浆密度；h为埋深。例如考虑浆液压力30 MPa时，对应浆液粒子的初始速度为199 m/s。

为阻止计算过程中水、土粒子的穿透，边界粒子层设置为4层，采用六方堆积方式布置。经试算，在水、土粒子间距设为0.1 m，浆液粒子间距设为0.05 m时，计算可较为稳定。模型中最终共包含11 260颗粒子，其中土粒子数量4 980颗，水粒子数量4 650颗，浆液粒子950颗，固壁边界粒子数量680颗。

以上海地区②₃灰色粉性土为背景，参照上海市岩土工程勘察规范^[17] (DGJ08-37—2012)，模型主要参数取值如表 1和表2所示。

2.2 土体变形过程分析

图 5展示了高压浆液挤压饱和土体发生变形的过程。在水泥浆高速冲击、侵蚀和挤压作用下，土体侧表面形成一个“匚”字形凹槽，其高度约等于浆液深度，宽度随着模拟时间的推进逐渐增大，土体在较大范围内出现明显的侧向变形。同时，总位移矢量图显示，地表以下一定范围内还存在向上移动的趋势，随着距右边界距离的增加，地表隆起量值先增大后减小，大致在距离右边界1 m位置处达到最大值。

分别提取距离右边界2 m和地表以下埋深约3 m两个剖面的土体水平位移曲线，如图 6所示。土体侧向变形沿深度大致呈现“抛物线”形分布，水泥浆冲击区域对应土体变形较大，最大值出现在地表以下约1.5 m位置；冲击区域以下土体也会产生一定的水平位移；土体变形随着距右边界距离的增大而逐渐减小，在距离浆液3 m范围内，土体变形量收敛较快，至模型左边界位置处(水平距离5 m，约10倍浆液宽度)，土体变形仍未收敛为零。

图 7和图 8绘制了有效应力增量和超孔隙水压力的模拟结果。可以看出，土体有效应力、超孔隙水压力在浆液冲击土体的瞬间激增，影响范围较大。在模拟的终止时刻，超孔隙水压力最大值达到165.2 kPa。较高超孔隙水压力可能对周围环境产生不良影响。

图 9绘制了不同浆液压力下，距离右边界2 m剖面的土体变形情况。平均来看，浆液压力每增加5 MPa，土体侧向变形和地表变形最大值分别将增加22.4%和20.2%。可见，为控制对周围土体的扰动作用，应选择与施工地层土体特性相匹配的施工参数和浆液材料。

根据工程现场报告，在高压浆液作用下，周边土体侧向受扰动范围可达15~20倍的浆液宽度^[18]，地表最大隆起变形大致在直接挤压区域临近位置产生^[5]，并且挤土效应随水平距离增加而降低^[2]。本研究模拟结果与现场实测的土体变形情况在定性上取得较好的一致性，说明采用SPH方法模拟高压浆液挤压饱和土体过程是可行的。

3 结论

本研究引入光滑粒子流体动力学(SPH)方法，建立高压浆液挤压饱和土体变形问题的数值分析模型，为该类工程问题的模拟提供了新途径。

(1) 选择了适当的控制方程和数值技术，并针对大变形问题需要，引入人工黏度项、应力变化率等。采用了基于混合物理论的水-土耦合方法，将渗流力和孔隙水压力作为水-土相互作用力。采用C++语言编译了相应的SPH分析程序。

(2) 基于SPH模型，分析了饱和土体侧向位移和地表隆起发展过程和规律，研究了地层中孔隙水压力及土体有效应力的分布情况，探讨了浆液压力对土体变形的影响。模拟结果与实际工程在定性趋势上有较好的一致性。

事实上，高压旋喷桩施工过程中，包含高压浆液对周围土体的冲切、劈裂、渗透和挤压的复杂作用。作为初步探究，本研究首先探讨了高压浆液对饱和土体的挤压效应。下一步可考虑在土体的控制方程中加入断裂本构，并引入更加复杂的浆-土耦合关系，分析高压水泥浆冲切、破坏土体的过程，实现饱和土体中高压旋喷桩施工更真实的模拟效果。