基于非均质地层模型的桩基注浆数值模拟

，， ，，，

(三峡大学 a.防灾减灾湖北省重点实验室；b.土木与建筑学院，湖北 宜昌 443002)

1 研究背景

提高桩侧和桩端阻力是工程中迫切需要解决的研究课题。桩端注浆和桩侧注浆是提高桩承载力的有力措施之一[1]，但注浆过程中桩周土体物理力学参数变化过程及注浆效果仍有待研究。

由于土具有非均质性和原位变异性[2]，土层物理性质参数在一定范围内波动，土体物理性质力学特征及空间分布呈非均匀性，导致注浆浆液扩散路径和范围具有不确定性。张继周等[3]运用随机场理论，将土层剖面模拟为随机场而非传统意义上的随机变量，提出了基于随机场理论并考虑空间趋势分量的土性参数变异性统计方法，系统地统计分析了苏中腹地湖相软土土性参数的变异特性，对比了实验室物理性质指标参数、变形参数和强度参数变异性特征，分析了土性参数的变异性及相关性；蒋水华等[4]在蒙特卡洛模拟(MCS)和极限平衡分析框架下提出了一种有效的考虑参数二维空间变异性的边坡失稳风险定量评估方法。

分形作为分析自然界极不规则和极其复杂现象的数学工具，可通过迭代计算来描述精细的任意小尺度结构等特点[5-6]。岩土介质的孔隙分布、孔隙数量和孔隙表面积都具有分形特性[7-8]。陈建飞[9]利用土性参数的相关关系公式，以随机分形方法生成的地层孔隙率为基本量，实现了对未知数据地层土体密度、渗透率、压缩模量等土性参数的空间随机分布的预测，生成了土性参数随机分布的地质剖面。

秦鹏飞等[10]进行了砂砾石土灌浆浆液扩散半径试验研究。雷进生等[11-12]基于渗流物理场和扩散物理场耦合理论，研究了浆液扩散的时变效应、土层参数的动态变化；燕乔等[13]根据变态混凝土的特点，建立了相关注浆模型，得到了可以指导变态混凝土注浆施工的公式。

本文结合桩周注浆加固工程实例，根据工程地质勘查资料和数据，基于随机分形原理进行模拟，得到了孔隙率、渗透系数、随机分布剖面图，并形成了二维随机非均质地层，并在此基础上构造相关模型，进行桩基周围土体的注浆模拟，探讨桩基与注浆加固共同作用的效果，评价桩基后注浆在工程中的利用价值。

2 注浆作用下土体参数的动态变化

2.1 孔隙率的动态变化

土体注浆加固过程中，注浆压力大小会改变土体颗粒之间的空间分布。在实际工程注浆加固时，土层的物理力学参数将随着渗透压力呈动态变化的特点，直接影响到浆液的渗透扩散过程以及注浆加固的效果。因此，为了更加贴近实际地反映出软弱地基土体注浆加固的过程，必须将土层物性参数的动态变化特点考虑进去[14]。

文献[15]中给出了渗透压力下土体孔隙率动态变化(不考虑温度效应)的关系式，即

(1)

式中：n为实时动态孔隙率；n0为初始孔隙率；εV为体积应变；Ks为体积模量，Ks=E/[3(1-2ν)],E为弹性模量,ν为泊松比；Δp为压力差,Δp=p-p0(p0为初始压力,p为当前压力)。

2.2 渗透率的动态变化

注浆过程中，土体的固结强度和固结变形受到土层多孔介质渗透率的影响，而渗透率又与土层孔隙率等物理参数相关。因此，流固耦合注浆数值模拟要考虑渗透率和孔隙率等物性参数动态变化过程。

土体渗透率的动态变化(不考虑温度变化)关系式[11]为

式中：κ为实时动态渗透率；κ0为初始渗透率。

2.3 其他参数的动态变化

注浆过程中，孔隙率的变化也会造成土体密度和弹性模量的变化，改变土层的物理指标和力学参数，在实际计算中应考虑实时动态孔隙率。考虑浆液和土体共同作用的平均密度为

(3)

式中：ρs为固相颗粒密度(kg/m3)；ρf为浆液密度(kg/m3)。

同时注浆过程会改变土体的结构，从而影响土体的力学参数，对计算结果产生影响。在应力场的计算中，弹性模量的取值直接与动态的孔隙率相关，弹性模量与孔隙率的关系式[12]为

(4)

式中a，b，c是常数，根据本次工程勘察条件分别取a=45.50，b=0.056 78，c=0.365 4。

3 二维随机非均质地层模型的构建

3.1 构建非均质地质模型的基本方法

土层整体参数与土层局部参数之间存在着因果关系。基于勘察数据采用随机分形算法构建非均质地质模型的基本方法如下：

(1)分析和预处理地质勘探数据。

(2)以已知点坐标为控制点，该坐标点对应的土性参数地质信息作为数据源，推断控制点间地质特征的分布情况，并绘制勘探区域地质剖面。

(3)基于地勘控制点的平面或者空间分布，分别从上下左右方向逐步构建平面四边形，且平面四边形有序相连，在数据不足的情况下将四边形改成三角形，体内无交叉线。

(4)土层物理力学参数分形分布场按随机分形插值算法生成，并采用二维随机分形算法生成地质剖面的土性参数分形分布场。

(5)依据其他参数与基本量之间的函数关系，采用以上方法产生其他参数分形分布场。

3.2 土层参数基本参量的选取

随机分形二维算法是将地勘数据进行分析划分之后，形成许多个四边形，再运用随机分形算法可以构建成非均匀地质模型。

对于正方形的4个点1，2，3，4，第1步获取这4个点的参数值；第2步分别计算正方形中点0和边中点5，6，7，8的参数值。中点的参数值为

I0=(I1+I2+I3+I4)/4+g0。

(5)

边中点5,6,7,8的参数值分别为：

I5=(I1+I2+I0)/3+g5；

(6)

I6=(I1+I3+I0)/3+g6；

(7)

I7=(I2+I4+I0)/3+g7；

(8)

I8=(I3+I4+I0)/3+g8。

(9)

(10)

式中：H是赫斯特指数；σ是原始数据的标准方差；G是服从N(0，1)分布的高斯随机数；d为网格间距;g5，g6，g7，g8的计算与g0类似。

如图1所示，Diamond-Square算法以正方形ABDC的4个端点为初始点，将相邻点之间的线段均分，每均分一次便可得到一部分点，其中各个点均是以周围点的参数值平均为基准值，然后再加上一个高斯随机数。

图1 Diamond-Square算法实现原理Fig.1 Principle of Diamond-Square algorithm

Diamond-Square算法的递归细分过程如下[16]。

(1)Diamond步：获取正方形4个点的参数值，在中心生成一个随机值，中点的参数值等于中心点周边的4点参数值加上一个随机增量，形成棱形网格。

(2)Square步：取4个点组成一个棱形，在棱形的中心处生成一个随机值，取平均角点参数值再加上与Diamond步相同的增量，计算出每边中点值。这样又得到一个正方形格网。

按以上方法每划分一次得到1个种子正方形并经过一次细分过程得到4个方形，第2次得到16个，第3次64个，数目等于22i，其中i为细分次数，重复该过程，直到达到要求为止。

按以上随机分形算法计算随机增量g0，g5，g6，g7，g8等，从式(10)中可以看出这些量与赫斯特指数H、标准方差σ和随机增量G之间相关，其中H能够很直观地反映出事物的复杂性和非均匀性。

3.3 非均质地质模型的构建

以宜昌某实际工程为例，基于地勘成果对研究区域内地层的孔隙率和弹性模量进行随机分形数值模拟。对研究区域进行地基加固处理设计,根据已知的压缩模量可计算得到粉质黏土的弹性模量指标,见表1，粉质黏土层厚度为0～8.1 m。在研究区域中选取2个地质剖面Ⅰ-Ⅰ(孔点编号A1，A2，A3,A4)、Ⅱ-Ⅱ(孔点编号B1，B2，B3,B4)，各孔点坐标及对应孔隙率指标见表2。

表1 粉质黏土层弹性模量指标Table 1 Elastic modulus of silty clay layer

表2 各剖面模拟计算区域孔隙率指标值

注:A1,A2,A3,A4点的平面相对坐标(x,y)分别为(0，0),(0，1.5),(0，3.5),(0，5.5),单位为m；B1,B2,B3,B4点的平面相对坐标(x,y)分别为(4，0),(4，1.5),(4，3.5),(4，5.5),单位为m

土体初始弹性模量和渗透系数可以根据土力学公式计算得出，如式(11)和式(12)。

(11)

k=κρg/μg。

(12)

式中：k为渗透系数；ρg为流体比重；μg为流体黏度。

基于Diamond-Square算法，利用MatLab进行编程，可以得到基于随机分形的土层孔隙率剖面如图2和图3。

图2 基于随机分形的孔隙率二维云图Fig.2 Two-dimensional nephogram of porosity based on random fractal

图3 基于随机分形的孔隙率三维云图Fig.3 Three-dimensional nephogram of porosity based on random fractal

注浆加固前，在基础传递荷载作用下，泊松比为0.3，弹性模量最大为13.7 MPa，计算得到土层密度最大为1 827 kg/m3、最大位移为113 mm，Mises应力为293 kPa。

4 桩基与注浆加固共同作用数值模拟及效果分析

4.1 桩基后注浆加固模型

根据研究区域地层资料取剖面范围为8 m×8 m，计算采用三角形网格，单元数6 903个，自由度41 418个，计算时长600 s。桩基与注浆加固模型网格划分如图4。

图4 地基网格划分Fig.4 Mesh generation of foundation

边界条件：加固区土体底部固定约束，两侧滚轴约束，顶部自由约束。

荷载条件：考虑注浆施工工艺条件和上部结构对地基承载力的要求，桩基为混凝土灌注桩，直径1.2 m，深度为5 m，距地面1.0 m。基础传递荷载的局部压力为300 kPa。桩基两侧注浆孔直径为0.08 m，注浆压力0.3 MPa，水灰比0.7。

4.2 桩基与注浆加固共同作用效果分析

通过土层孔隙率动态模型、渗透率动态模型及其他参数的动态模型在有限元模拟软件中进行调用[11]，结合流体的运动方程，得到注浆过程中土层的实时动态变化的物理参数。

利用MatLab软件对有限元计算数据提取和分析，可以计算注浆结束后有限元模型中每个网格结点在注浆结束后土体的弹性模量、泊松比以及密度等参数。

考虑浆液流动特性和动态扩散过程，假定注浆压力Pf>0.1 MPa时，浆液可渗透土层，形成注浆结石体。不同时刻的注浆加固浆液扩散过程见图5,浆液渗透的范围为红色区域。

图5 注浆加固作用扩散过程Fig.5 Diffusion of grouting slurry

由于混凝土灌注桩不可渗透，故浆液在注浆压力的作用下，沿着桩侧和桩端渗透。随着注浆浆液的渗透、劈裂和挤密作用，桩侧和桩端周围形成浆液和土体的结石体，对桩周泥土进行置换，空隙进行充填。

由图5中可以看出，受地层非均质土层物理和力学参数空间变化的影响，不同部位和不同方向的扩散距离并不相同，注浆范围随时间也呈现不均匀和非等速变化。t=40 s前变化明显，40 s之后逐渐趋于稳定。t=5 s时，浆液开始沿与注浆孔相连接的低密度、高孔隙率和高渗透率的区域向前扩展运动。t=40 s时，浆液以指形扩散模式，从注浆孔壁以优势路径(低密度、高孔隙率和高渗透率区域)为主向四周扩散。从浆液已到达的某一位置开始，继续沿着优势路径运动继续扩散，形成指形扩散锋面。t=600 s时，整体浆液扩散范围已基本保持不变。

桩基周边土体注浆加固过程中，上部恒荷载作用下的土层物理力学参数发生变化，注浆扩散范围稳定及注浆过程终止后Ⅰ-Ⅰ剖面土层物理力学参数如图6所示。

图6 荷载作用下桩基与注浆加固后土层物理参数变化Fig.6 Changes in the physical parameters of pile foundation and soil layer after grouting reinforcement under load

由图6可以看出，混凝土灌注桩被注浆结石体所包围形成一个扩大的柱体，使桩端受力面积增大，桩土接触面的条件得到改善，进而使桩侧土单位面积的侧阻力提高。

将图6数据与土层数据进行对比可以看出：注浆加固前，在基础传递荷载作用下，I-I剖面土层密度最大值为1 827 kg/m3;注浆加固后，注浆复合地基区域结石体密度最大值为2 314 kg/m3，相较未注浆加固状态增加了约26.7%，周边未注浆区域受注浆挤压作用影响，密度也发生明显变化。

注浆加固对荷载作用变形抑制作用明显。以Ⅰ-Ⅰ剖面为例，未加固前桩基受上部恒荷载作用时，计算桩周土体最大位移为113 mm，而注浆加固后的注浆固结增强体最大位移为43.7 mm。随着注浆压力增加和注浆过程持续，注浆结束后注浆扩散区域的结石体及周边未注浆土体的弹性模量、泊松比、Mises应力都明显提高，加固后桩周结石增强体最大弹性模量超过600 MPa，约为加固前的44倍，土层Mises应力增加到4 000～6 900 kPa，为加固前的13～23倍，应力状况也得到改善，说明注浆行为大幅减小了地基的变形。

5 结 语

(1) 土性参数是影响注浆扩散的重要指标，考虑土性参数，基于随机分形方法构建地层剖面随机模型，可以实现桩周注浆模拟计算及加固效果评估。

(2) 考虑浆液流动特性和动态扩散过程，通过土层孔隙率动态模型、渗透率动态模型及其他参数的动态模型在有限元模拟软件中进行调用，结合流体的运动方程，得到桩周浆土结石体形成过程和不同时刻的分布范围。

(3) 受地层非均质和土性参数空间变化的影响，注浆范围呈现空间上的不均匀和时间上的非等速变化。注浆过程中土体的弹性模量、泊松比以及密度等参数发生动态变化。

(4) 注浆过程终止后，混凝土灌注桩被注浆结石体所包围，结石体强度增大，土层物理和力学参数明显提高，结石体形成及桩周土性参数的改善明显提高了桩基承载力，桩基后注浆能有效控制地基沉降,提高桩端承载力。合理利用桩基后注浆技术，能有效控制建筑地基的沉降量，适当减少桩长,还可达到节约工程造价的目的。

参考文献：

[1] 张忠苗, 吴世明, 包 风. 钻孔灌注桩桩底后注浆机理与应用研究[J]. 岩土工程学报，1999，21(6)：681-686.

[2] HUANG J, LYAMIN A V, GRIFFITHS D V,etal. Quantitative Risk Assessment of Landslide by Limit Analysis and Random Fields[J]. Computers and Geotechnics, 2013, 53: 60-67.

[3] 张继周, 缪林昌, 陈俊波. 苏中腹地湖相软土土性参数变异性统计描述[J]. 岩土力学, 2010, 31(2): 471-477.

[4] 蒋水华，杨建华，姚 池，等. 考虑土体参数空间变异性边坡失稳风险定量评估[J]. 工程力学, 2018, 35(1): 136-147.

[5] FALCONER K. Fractal Geometry: Mathematical Foundations and Applications [M]. Chichester: Wiley, 2003.

[6] 肖高逾, 周源华. 基于分形插值的地貌生成技术[J]. 上海交通大学学报, 2000, 34(5): 705-707.

[7] 宫英振, 牛海霞, 董正茂, 等. 分形多孔介质孔道网络模型的构建[J]. 农业机械学报, 2009, 40(11): 109-114，108.

[8] BARTOLI F, PHILIPPY R, DOIRISSE M,etal. Structure and Self-similarity in Silty and Sandy Soils: The Fractal Approach[J]. Journal of Soil Science, 1991, 42(2): 167-185.

[9] 陈建飞. 非均匀地层物性参数随机模型及其应用研究[D]. 宜昌:三峡大学, 2014.

[10] 秦鹏飞, 符 平, 王 春, 等. 砂砾石土灌浆浆液扩散半径试验研究[J]. 长江科学院院报，2014，31(5):84-86, 101.

[11] 雷进生, 刘 非, 王乾峰, 等. 非均质土层的注浆扩散特性与加固力学行为研究[J]. 岩土工程学报, 2015, 37(12): 2245-2253.

[12] 雷进生, 刘 非, 彭 刚, 等. 考虑参数动态变化和相互关联的浆液扩散范围研究[J]. 长江科学院院报, 2016, 33(2): 57-61.

[13] 燕 乔, 徐 枫, 莫世远, 等. 变态混凝土成孔注浆浆液扩散研究[J]. 长江科学院院报, 2016, 33(2): 91-93, 99.

[14] 冉启全, 李士伦. 流固耦合油藏数值模拟中物性参数动态模型研究[J]. 石油勘探与开发, 1997, 24(3): 61-65.

[15] CHENG P, LI L, TANG J,etal. Application of Time-varying Viscous Grout in Gravel-foundation Anti-seepage Treatment[J]. Journal of Hydrodynamics, 2011, 23(3): 391-397.

[16] 吴 稳, 周云轩, 李 行, 等. 分形布朗运动在水下地形模拟中的应用[J]. 计算机工程与应用, 2010, 46(25): 242-245.