Numerical Analysis of Multiphysics Coupling of Grout Penetration
-
摘要:
非饱和地层中的注浆渗透扩散是一个复杂的多物理场过程，为更加精确分析浆液在饱和与非饱和地层中的扩散特性，并估算注浆渗透扩散范围和注浆挤密区域，以混合理论为基础，建立一个非饱和多孔介质中的多物理场耦合模型. 通过ABAQUS二次开发，构建一种新型八节点五自由度四边形轴对称Serendipity单元，实现对注浆过程中土体变形、土体孔隙率、孔隙压力和浆液浓度分布的数值求解，以及对土体饱和度、渗透系数等状态变量的实时更新；结合一个三维轴对称注浆算例，分析浆液水灰比、注浆压力、土体初始干密度以及土体初始含水率对粉砂地层注浆效果的影响，并得到浆液水平和竖向扩散距离随不同因素变化的拟合曲线. 研究结果表明：浆液扩散范围受水灰比影响最显著，受注浆压力影响次之，受含水率和干密度影响最小；浆液扩散范围随水灰比增加而增长，水灰比大于1.0时增长显著；注浆管壁周围会形成挤密区域，浆液扩散区域内土体同时受到注浆压力的挤压和孔隙压力的支撑作用；随着远离注浆管壁，土体孔隙率在挤密区域内逐渐减小，在挤密区域外逐渐恢复，且挤密区域随注浆压力的增加而增大；研究成果可为土体注浆加固范围计算提供理论指导.
Abstract:Grout penetration in the unsaturated formation is a complex multiphysics process. In order to analyze more accurately the penetration characteristics of grout in saturated and unsaturated formations and estimate the range of grout penetration and grout compaction area, a multiphysics coupling model for unsaturated porous media was developed based on mixture theory. A novel five-degree-of-freedom eight-node quadrilateral axisymmetric Serendipity element was constructed through the secondary development of ABAQUS, so as to realize the numerical solution of soil deformation, soil porosity, pore pressure, and grout concentration distribution during grouting. The state variables, such as soil saturation and permeability coefficient, were updated in real time. The effects of grout water-cement ratio, grouting pressure, initial dry density, and water content of soil on the grouting in silty sand were analyzed with a three-dimensional axisymmetric grouting example, and the fitting curves of horizontal and vertical diffusion distance of grout with the above factors were obtained. The results show that the range of grout penetration is most significantly influenced by the water-cement ratio, followed by the grouting pressure, and it is least influenced by the water content and dry density of soil; the range of grout penetration increases with the increase in the water-cement ratio, especially when the water-cement ratio is greater than 1.0; a compacted area will be formed around the grouting pipe, where the soil is simultaneously subject to the grouting pressure and pore pressure; with the increase in the distance from the grouting pipe, the porosity of the soil gradually decreases in the compacted area and gradually recovers outside the compacted area, while the compacted area increases with the increase in the grouting pressure. The research results can provide theoretical guidance for calculating the grouting reinforcement in the soil.
-
Key words:
- grout penetration /
- multiphysics /
- finite element method /
- unsaturated soil
-
随着我国经济建设的蓬勃发展，人们对交通基础设施建设的需求与日俱增，同时，我国隧道与地下工程建设面临地质构造复杂、工程灾害频发的严峻考验[1]. 大量工程实践证明，注浆是解决地下工程灾害的重要方法之一，可有效实现破碎地层加固、减小沉降、突水突泥防控等一系列治理目标，因此，在隧道、矿山等各类地下工程建设中均具有广泛的应用[2].
然而，注浆是一个复杂的水力耦合过程，当注浆材料在压力作用下进入注浆区域后，浆液会在注浆压力的作用下对周围土体进行挤压，同时渗透扩散进入周围土层[3]. 加之注浆过程的隐蔽性以及地质条件的复杂性，导致注浆扩散及加固机理的研究十分困难，注浆理论的发展远远落后于工程实践[4]. 总体而言，目前的注浆工程大多仍依靠实际工程经验，工程设计缺乏严密的理论依据. 鉴于此，有必要提出一种针对浆液在地层中迁移扩散的计算方法来分析浆液在注浆地层中的扩散特性，以达到优化注浆参数的目的.
针对地层中的注浆研究，现有方法多根据浆液扩散压力和地下水压力之间的平衡确定浆液扩散范围[5]. 事实上，非饱和地层中注浆扩散是一个复杂的多物理场水力耦合过程. 在注浆过程中，浆液会在注浆压力的作用下驱替孔隙水和空气，与此同时，浆液在地层中扩散改变了浆液的浓度，进而改变了浆液在地层中的流动特性（如黏度等）以及土体特性（如孔隙水压力分布、饱和度、土体孔隙率等）[6-9]. 因此，为更加精确地分析浆液在饱和与非饱和地层中的扩散特性，并估算注浆渗透扩散范围，亟需一个固-液-气多物理场耦合模型.
本文以连续介质力学框架下混合物理论为基础[10]，建立一个准静态条件下的多孔介质固-液-气多物理场耦合模型，用以分析和预测浆液在饱和与非饱和地层中的扩散. 该方法以基于混合物理论的质量平衡方程描述浆液的扩散及其对土体饱和度、渗透性等的影响. 同时，该模型考虑了浆液浓度和黏度随时间的变化、孔隙压力及基质吸力变化对土体的作用. 然后，通过ABAQUS的二次开发，建立一种新型的八节点五自由度四边形轴对称Serendipity单元，实现对注浆过程中的土体变形、土体孔隙率、孔隙压力和浆液浓度分布的数值求解，以及对土体饱和度、渗透系数等状态变量的实时更新. 最后，结合一个三维轴对称注浆算例，分析了浆液水灰比、注浆压力、土体初始干密度以及土体初始含水率对粉砂地层注浆效果的影响，并结合数值计算结果，分别得到不同水灰比、注浆压力、干密度和含水率下的浆液水平扩散距离和竖向扩散距离经验公式.
1. 非饱和土注浆渗流模型
1.1 质量平衡
浆液在非饱和土中的渗流可视为由固-液-气三相组分构成的混合物系统. 其中：固相组分为土体颗粒（ss）；液相组分即孔隙液体，包括土体中初始存在的孔隙水（ff）及注入土体的浆液悬浮液（fp，包括浆液中的水和颗粒材料）；气相组分为孔隙中的空气（a）. 分别用nss、nfp、nff、na表示4个组分的体积分数[10-11]. 对于给定体积（dV）的特征单元，单个组分的体积分数nα为
$$ {n_\alpha }\left( {x,t} \right) = \frac{{{\text{d}}{V_\alpha }\left( {x,t} \right)}}{{{\text{d}}V}} \text{，} $$ (1) 式中：${\text{d}}{V_\alpha } $为单个组分α所占体积，α$\in ${ss,ff,fp,a}；x为空间坐标；t为时间.
因此，单位混合物体积中各组分的质量${\rho _\alpha }$可表示为
$$ {\rho _\alpha }\left( {x,t} \right) = {n_\alpha }\rho _{{\mathrm{R}}\alpha} \text{，} $$ (2) 式中：$ \rho _{{\mathrm{R}}\alpha} $为混合物中各组分的固有密度.
土体中的孔隙率ϕ(x,t)、孔隙液体饱和度sw(x,t)和浆液相对浓度c(x,t)分别为
$$ \phi\left(x,t\right)=1-n_{\mathrm{ss}}=n_{\mathrm{ff}}+n_{\mathrm{fp}}+n_{\mathrm{a}}\text{，} $$ (3) $$ s_{\mathrm{w}}\left(x,t\right)=\frac{n_{\mathrm{ff}}+n_{\mathrm{fp}}}{n_{\mathrm{ff}}+n_{\mathrm{fp}}+n_{\mathrm{a}}}=\frac{n_{\mathrm{ff}}+n_{\mathrm{fp}}}{\phi\left(x,t\right)}\text{，} $$ (4) $$ c\left( {x,t} \right) = \frac{{{n_{{\mathrm{fp}}}}}}{{{n_{{\mathrm{ff}}}} + {n_{{\mathrm{fp}}}}}} = \frac{{{n_{{\mathrm{fp}}}}}}{{{s_{\mathrm{w}}}\phi (x,t)}} . $$ (5) 当sw=1时，土体中的孔隙完全由液体填充；当c=0时，孔隙液体中不含有注入的浆液悬浮液；当c=1时，孔隙液体仅包含浆液悬浮液.
在连续介质力学的框架下，混合物系统特征单元体积中的各组分质量平衡通式[10-11]可表示为
$$ \frac{{\partial {{\rho _\alpha }} }}{{\partial t}} + {\text{div}}\left( {{\rho_\alpha }{{\boldsymbol{v}}_\alpha }} \right) = {\rho _{{\mathrm{ex}},\alpha }} \text{，} $$ (6) 式中：${\rho _{{\mathrm{ex}},\alpha }} $为组分$\alpha $的质量交换项，vα为相应组分的速度矢量；div(·)为矢量散度算子.
假定：1） 本文模型是以连续介质力学框架下的混合物理论为基础的宏观模型，不考虑非饱和土渗流的细微观尺度特性，而是通过水土特征关系以及饱和度与渗透系数折减关系考虑非饱和土的渗流特性；2） 浆液悬浮液中的材料颗粒随着孔隙水以相同速度运动；3） 土体中空气与大气层完全联通，且土体中的空气可快速自由流动，在准静态条件下，超孔隙气体压力可以快速地扩散，从而使孔隙气体和大气压力时刻保持平衡，因此，在平衡方程中，可不考虑气体压力的影响；4） 不考虑土颗粒的脱落或浆液材料在土体孔隙中的淤堵，即ρex,α=0；5） 不考虑土体颗粒及孔隙液体的任何受力变形.
将式（2）代入式（6），可得
$$ \frac{{\partial {{n_\alpha }} }}{{\partial t}} + {\text{div}}\left( {{n_\alpha }{{\boldsymbol{v}}_\alpha }} \right) = 0 . $$ (7) 结合式（3）～（5）对多孔介质中各个组分含量的定义，利用式（7）将各组分进行展开，可得
$$ \frac{{\partial \phi }}{{\partial t}} - {\text{div}}\; {{{\boldsymbol{v}}_{\mathrm{s}}}} + {\text{div}}\left( {\phi {{\boldsymbol{v}}_{\mathrm{s}}}} \right) = 0\text{，} $$ (8) $$ \frac{{\partial \left( {{s_{\mathrm{w}}}\phi } \right)}}{{\partial t}} + {\text{div}} \;{{{\boldsymbol{q}}_{\mathrm{w}}}} + {\text{div}}\left( {{s_{\mathrm{w}}}\phi {{\boldsymbol{v}}_{\mathrm{s}}}} \right) = 0\text{，} $$ (9) $$ \frac{{\partial \left( {{s_{\mathrm{w}}}c\phi } \right)}}{{\partial t}} + {\text{div}}\left( {c{{\boldsymbol{q}}_{\rm{w}}}} \right) + {\text{div}}\left( {{s_{\mathrm{w}}}c\phi {{\boldsymbol{v}}_{\mathrm{s}}}} \right) = 0\text{，} $$ (10) 式中：qw为液相组分的渗流速度（式11）；vs为固相组分的位移速度（式12）.
$$ {{\boldsymbol{q}}_{\mathrm{w}}} = \phi {s_{\mathrm{w}}}\left( {{{\boldsymbol{v}}_{\mathrm{f}}} - {{\boldsymbol{v}}_{\mathrm{s}}}} \right)\text{，} $$ (11) $$ \left\{\begin{array}{l} {{\boldsymbol{v}}_{\mathrm{s}}} = \dfrac{{\partial {\boldsymbol{u}}\left( {x,t} \right)}}{{\partial t}} \text{，}\\ \dfrac{{\partial {{\varepsilon _{\mathrm{v}}}} }}{{\partial t}} = {\text{div}} \;{{{\boldsymbol{v}}_{\mathrm{s}}}} \text{，} \end{array}\right. $$ (12) 式中：u(x,t)为固相组分的位移矢量，vf为液相组分的位移速度矢量，εv为土体孔隙体积变化引起的土体体积应变.
式（8）描述了固相组分（土体颗粒ss）的质量平衡；式（9）描述了液相组分（包括土体中初始存在的孔隙水ff及注入土体的浆液悬浮液fp）的质量平衡；式（10）描述了注入土体的浆液悬浮液（fp）的质量平衡. 气相组分的含量变化由sw的变化表征. sw通过土体材料的水土特征曲线（SWCC）由土体基质吸力确定.
1.2 浆液渗流速度与渗压
根据Uzuoka等[12]对非饱和土体孔隙水的动量平衡推导可知，静力作用下的非饱和土体孔隙液体渗流速度可表示为
$$ {{\boldsymbol{q}}_{\mathrm{w}}} = - \frac{k}{\mu }\left( {{\text{grad}}\; {{p_{\mathrm{w}}}} - \rho {{g}}} \right)\text{，} $$ (13) 式中：pw为孔隙液体压力；k为土体固有渗透系数（m2），随孔隙液体饱和度及孔隙率而变化，如式（14）；$ \rho $为孔隙液体密度，如式（15）；$\mu $为孔隙液体的动态黏度（Pa·s），随浆液相对浓度呈线性变化[13-14]，如式（16）.
$$ k = {k_{\mathrm{i}}}{k_{\mathrm{s}}}{k_{\text{ϕ}} } = {k_{\mathrm{i}}}{ {{s_{\mathrm{w}} ^3}}}\frac{{{\phi ^3}}}{{{{\left( {1 - \phi } \right)}^2}}}\text{，} $$ (14) $$ \rho = c{\rho _{\mathrm{g}}} + \left( {1 - c} \right){\rho _{\mathrm{w}}}\text{，} $$ (15) $$ \mu = {\mu _{\mathrm{g}}}c + {\mu _{\mathrm{w}}}\left( {1 - c} \right)\text{，} $$ (16) 式中： ki为初始固有渗透系数，ks为非饱和土中与饱和度相关的渗透修正系数[13]，kϕ为Kozeny-Carman系数[15]，μg和μw分别为注入的浆液悬浮液和纯水的动态黏度，ρg为注入的浆液悬浮液密度，ρw为纯水密度.
浆液的黏度也会随着浆液的胶凝而逐渐增加. 根据文献[16]针对不同水胶比（胶凝材料为水泥）条件下的水泥基浆液黏度时变性研究结果，μg随时间呈指数变化[17]，如式（17）.
$$ \mu_{\mathrm{g}}=\mu_{\mathrm{g}0}\mathrm{e}^{\xi t}\text{，} $$ (17) 式中：μg0为浆液初始黏度，ξ为需要通过试验确定的材料常数.
1.3 多孔介质的力学控制方程
多孔固体的力学行为由式（18）控制.
$$ {\left( {{\sigma _{ij}} - {s_{\mathrm{w}}}{p_{\mathrm{w}}}} \right),_{j}} = {w_i}\text{，} $$ (18) $$ {\text{d}}{\sigma _{ij}} = {\text{d}}{\sigma _{ij1}} + B{\delta _{ij}}{\text{d}}{p_w} = D_{ijkl}{\text{d}}{\varepsilon _{kl}} + B{\delta _{ij}}{\text{d}}{p_w}\text{，} $$ (19) 式中：$ {\sigma _{ij}} $为总应力张量的分量，$ {w_i} $为单位体积的体力，$ {\sigma _{ij1}} $为有效应力张量的分量，$ {\varepsilon _{ij}} $为总应变张量的分量，$ D_{ijkl} $为弹塑性刚度矩阵的分量，$ B $为Biot有效应力系数（本研究中取1），$ {\delta _{ij}} $为克罗内克符号，下标i、j、k、l表示空间张量分量，(·), j表示对j方向求偏导致.
模型中的渗流-力耦合可以分为2个方面：1） 根据式（8）～（10），浆液在注浆压力的作用下扩散，导致土体孔隙率、渗透系数、浆液浓度发生时空变化，从而进一步引起土体孔隙压力及饱和度的改变. 因此，结合式（18）、（19），在边界注浆压力及孔隙压力、饱和度及应力平衡的共同作用下，土体的有效应力及位移响应不断发展. 2） 土体的位移及变形响应也会影响浆液的扩散过程，通过式（8）～（10）中vs实现. 因此，本模型可以实现对注浆扩散过程中，渗流、应力及变形相互作用的耦合分析. 需要注意的是，本模型假设土体中的空气与大气层完全联通，忽略孔隙气压力的变化. 气体含量的变化以及其对渗流的影响是通过式（8）～（10）中的饱和度结合SWCC的基质吸力，以及饱和度对渗透系数的折减来考虑. 同时，本文模型同样适用于对饱和土体中的渗流与注浆问题的分析，以及对非饱和土体在渗流或注浆作用下逐渐饱和过程的分析. 对于饱和介质，控制方程中的sw始终等于1，模型即退化为饱和土体中的渗流及浆液扩散问题.
1.4 模型参数
本文中所提出的非饱和土注浆渗流模型参数主要包括：
1） 初始状态变量（初始孔隙率ϕ0、土体中的初始浆液浓度c0、土体初始饱和度sw0、ρg、ρw）.
2） 渗流相关参数（ki、ks、浆液动态黏度参数μg0、ξ）. ks由非饱和土特性根据饱和度确定，本次研究取${k_{\mathrm{s}}} = {{s_{\mathrm{w}}^3}}$[13]，结合式（14）和初始孔隙率ϕ0，ki可换算为饱和土的初始渗透系数k0.
3） 力学参数. 取决于土体本构模型. 本次分析的算例中，土体本构模型为线弹性模型. 若实验数据充足，在后续的研究中，可采用基于临界状态和颗粒级配的和考虑注浆加固作用的弹塑性模型来更准确地模拟浆液扩散过程中的土体力学行为.
1.5 数值模拟求解
式（8）～（10）、（18）构成了非饱和多孔介质中的固-液-气多物理场耦合模型控制方程，其主要的未知量包括： $ {\boldsymbol{u}}\left( {x,t} \right) $、$ {p_{\mathrm{w}}} $、$ \phi $和$ c $. 其他未知量，如$ {v_{\mathrm{s}}} $、$ {q_{\mathrm{w}}} $、$ \bar \rho $、$ k $及$ {s_{\mathrm{w}}} $可由式（11）～（17）确定.
式（8）～（10）、（18）的数值求解是通过将其弱形式植入有限元软件ABAQUS进行二次开发实现[18]. ABAQUS允许用户通过用户自定义单元子程序（UEL）定义新的自由度和新的单元类型，因而，可利用ABAQUS求解器求解自定义偏微分方程组. 图1展示了ABAQUS调用UEL子程序的计算流程. 图中：$\Delta \varepsilon、 \Delta \sigma $分别为应变增量和应力增量. ABAQUS/UEL采用牛顿法对非线性方程组进行求解（式（20））.
图 1 ABAQUS调用UEL子程序计算流程Figure 1. Calculation flow chart of ABAQUS calling UEL subroutine$$ {\boldsymbol{A}}{{\boldsymbol{d}}_{\mathrm{N}}} = {\boldsymbol{R}} \text{，} $$ (20) 式中：${\boldsymbol{A}} $为待求解偏微分方程组弱形式的Jacobian矩阵，R为待求解偏微分方程组弱形式的节点不平衡力矩阵，dN为待求解节点自由度.
在UEL子程序中，用户需要根据传入UEL中的位移、位移增量、荷载及状态变量，计算应变、应变增量和节点荷载，并更新应力和状态变量，最终更新矩阵A和矩阵R. ABAQUS/UEL根据用户更新的A和R进行迭代求解.
本次研究通过ABAQUS/UEL开发一种新型的八节点四边形轴对称Serendipity单元，包含5个自由度（dN=[ux uy pw ϕ c]）和9个积分点，ux、uy分别为x、y方向的位移自由度，如图2所示.
图 2 新型五自由度轴对称单元示意Figure 2. Illustration of novel five-degree-of-freedom axisymmetric element为满足Inf-Sup或LBB （Ladyzhenskaya-Babuška-Brezzi）条件[19-20]，土体的位移自由度（ux，uy）为二阶插值离散，而其他3个变量（pw，ϕ，c）为线性插值离散. 因此，待求解变量可表示为
$$ \left\{\begin{array}{l}\boldsymbol{u}=(u_x,u_y)=\boldsymbol{N}_1\boldsymbol{u}_{\mathrm{N}}, \\ p_{\mathrm{w}}=\boldsymbol{N}_2\boldsymbol{P}_{\mathrm{wN}}, \\ \boldsymbol{\phi}=\boldsymbol{N}_2\mathbf{\boldsymbol{\boldsymbol{\phi}}}_{\mathrm{N}}, \\ c=\boldsymbol{N}_2\boldsymbol{c}_{\mathrm{N}},\end{array}\right. $$ (21) 式中：uN、PwN、ϕN、cN分别为对应待求解自由度的节点值组成的矩阵， N1和N2分别为二阶和一阶插值的单元形函数，其具体形式可参见文献[21].
式（8）～（10）、（18）的弱形式可表示为
$$ \begin{split} & {R_1} = \int {{\boldsymbol{B}}_1^{\mathrm{T}}\left( {\sigma_1 - {p_{\mathrm{w}}}} \right){\mathrm{d}}V} - \\ &\quad \oint {{\boldsymbol{N}}_1^{\mathrm{T}}{\tau _n}{\mathrm{d}}S} - \int {{\boldsymbol{N}}_1^{\mathrm{T}}w{\mathrm{d}}V} \end{split} \text{，} $$ (22) $$ \begin{split} & {R_2} = \frac{1}{{\Delta t}}\int { {{\boldsymbol{N}}_2^{\mathrm{T}}{{\boldsymbol{N}}_2}\Delta {{\boldsymbol{\phi}}_{\mathrm{N}} }} {\mathrm{d}}V} + \frac{{\phi - 1}}{{\Delta t}} \times \\ &\quad \int { {{\boldsymbol{N}}_2^{\mathrm{T}}{{\boldsymbol{m}}^{\mathrm{T}}}{{\boldsymbol{B}}_1}\Delta {\boldsymbol{u}}_{\mathrm{N}}} {\mathrm{d}}V} + \frac{1}{{\Delta t}}\int { {{\boldsymbol{N}}_2^{\mathrm{T}}{{{\boldsymbol{\phi}} }_{\mathrm{N}}^{\mathrm{T}}}{\boldsymbol{B}}_2^{\mathrm{T}}{{\boldsymbol{N}}_1}\Delta {\boldsymbol{u}}_{\mathrm{N}}} } {\mathrm{d}}V \end{split} \text{，} $$ (23) $$ \begin{split} & {R_3} = - \int { {{\boldsymbol{B}}_2^{\mathrm{T}}{q_{\mathrm{w}}}} {\mathrm{d}}V} {\text{ + }}\int {\frac{{{s_{\mathrm{w}}}}}{{\Delta t}} {{\boldsymbol{N}}_2^{\mathrm{T}}{{\boldsymbol{m}}^{\mathrm{T}}}{{\boldsymbol{B}}_1}\Delta {\boldsymbol{u}}_{\mathrm{N}}} {\mathrm{d}}V}+ \\ &\quad \int {\phi \frac{{\partial {s_{\mathrm{w}}}}}{{\partial {p_{\mathrm{w}}}}}\frac{1}{{\Delta t}} {{\boldsymbol{N}}_2^{\mathrm{T}}{\boldsymbol{N}}_2^{}\Delta {{\boldsymbol{p}}_{\rm{wN}}}} {\mathrm{d}}V} {\text{ + }}\oint {{\boldsymbol{N}}_2^{\mathrm{T}}q_{{\mathrm{w,n}}}^{}{\mathrm{d}}S}+ \\ &\quad \int {\phi \frac{{\partial {s_{\mathrm{w}}}}}{{\partial {p_{\mathrm{w}}}}}\frac{1}{{\Delta t}} {{\boldsymbol{N}}_2^{\mathrm{T}}{{\boldsymbol{p}}_{\rm{wN}}^{\mathrm{T}}}{\boldsymbol{B}}_2^{\mathrm{T}}{{\boldsymbol{N}}_1}\Delta {\boldsymbol{u}}_{\mathrm{N}}} {\mathrm{d}}V} ， \end{split} $$ (24) $$ \begin{split} & {R_4} = \frac{1}{{\Delta t}}\phi {s_{\mathrm{w}}}\int { {{\boldsymbol{N}}_2^{\mathrm{T}}{{\boldsymbol{N}}_2}\Delta {\boldsymbol{c}}_{\mathrm{N}}} {\mathrm{d}}V} + \int { {{\boldsymbol{N}}_2^{\mathrm{T}}{{\boldsymbol{c}}_{\mathrm{N}}^{\mathrm{T}}}{\boldsymbol{B}}_2^{\mathrm{T}}{q_{\mathrm{w}}}} {\mathrm{d}}V}+ \\ &\quad \phi {s_{\mathrm{w}}}\frac{1}{{\Delta t}}\int { {{\boldsymbol{N}}_2^{\mathrm{T}}{{\boldsymbol{c}}_{\mathrm{N}}^{\mathrm{T}}}{\boldsymbol{B}}_2^{\mathrm{T}}{{\boldsymbol{N}}_1}\Delta {\boldsymbol{u}}_{\mathrm{N}}} } {\mathrm{d}}V , \end{split} $$ (25) 式（22）～（25）中：σ1为有效应力；∆t为时间增量；∆ϕN为孔隙率增量矩阵；∆pwN为孔隙水压力增量矩阵；∆uN为位移增量矩阵；∆cN为浆液相对浓度增量矩阵；τn为沿法线方向的边界应力；m为单位矩阵，二维情况下是[1 1 0]，三维情况下是[1 1 1 0 0 0]；B1和B2分别是形函数N1和N2对局部坐标的空间导数[21]；W=[Wi]，qw,n为qw沿法线方向的分量.
上述非线性控制方程组可统一写成
$$ {\boldsymbol{R}}\left( {{{{\boldsymbol{u}}}},{{{p}}_{\mathrm{w}}},{\phi },{{c}}} \right) = 0 \text{，} $$ (26) 采用牛顿法求解非线性方程组的第m次迭代形式可表示为
$$ \begin{split} & \boldsymbol{A}\mathbf{\mathit{\boldsymbol{d}}}_{\mathrm{N}}=\boldsymbol{R}= \\ & \quad\left[ \begin{array}{*{20}{c}}\dfrac{\partial R_1}{\partial\boldsymbol{u}_{\mathrm{N}}} & \dfrac{\partial R_1}{\partial\boldsymbol{p}_{\mathrm{wN}}} & \dfrac{\partial R_1}{\partial\boldsymbol{\phi}_{\mathrm{N}}} & \dfrac{\partial R_1}{\partial\boldsymbol{c}_{\mathrm{N}}} \\ \dfrac{\partial R_2}{\partial\boldsymbol{u}_{\mathrm{N}}} & \dfrac{\partial R_2}{\partial\boldsymbol{p}_{\mathrm{wN}}} & \dfrac{\partial R_2}{\partial\boldsymbol{\phi}_{\mathrm{N}}} & \dfrac{\partial R_2}{\partial\boldsymbol{c}_{\mathrm{N}}} \\ \dfrac{\partial R_3}{\partial\boldsymbol{u}_{\mathrm{N}}} & \dfrac{\partial R_3}{\partial\boldsymbol{p}_{\mathrm{wN}}} & \dfrac{\partial R_3}{\partial\boldsymbol{\phi}_{\mathrm{N}}} & \dfrac{\partial R_3}{\partial\boldsymbol{c}_{\mathrm{N}}} \\ \dfrac{\partial R_4}{\partial\boldsymbol{u}_{\mathrm{N}}} & \dfrac{\partial R_4}{\partial\boldsymbol{p}_{\mathrm{wN}}} & \dfrac{\partial R_4}{\partial\boldsymbol{\phi}_{\mathrm{N}}} & \dfrac{\partial R_4}{\partial\boldsymbol{c}_{\mathrm{N}}}\end{array} \right]_{(m)}\left[ \begin{array}{*{20}{c}}\Delta\boldsymbol{u}_{\mathrm{N}} \\ \Delta\boldsymbol{p}_{\mathrm{wN}} \\ \Delta\boldsymbol{\phi}_{\mathrm{N}} \\ \Delta\boldsymbol{c}_{\mathrm{N}}\end{array} \right]_{(m+1)}=\\ & \quad -\left[ R_1 \quad R_2 \quad R_3 \quad R_4 \right]_{(m)}^{\mathrm{T}}.\end{split} $$ (27) Jacobian矩阵可依据待求解偏微分方程组的弱形式对待求解自由度的偏导数确定，节点不平衡力由待求解方程组弱形式的牛顿迭代残余值确定.
通过ABAQUS/UEL求解非线性偏微分方程组，可以直接利用ABAQUS进行前处理建模和后处理结果输出. 结合本文开发的新八节点五自由度四边形轴对称单元，可以模拟浆液在不同初边界值条件下（注浆压力、复合地层等）的运移以及地层的变形. 需要注意的是，实际注浆过程中浆液在土体孔隙中的淤堵需要结合大量试验以标定淤堵速率. 由于缺乏试验数据，本文不对浆液淤堵进行探讨，仅考虑浆液在注浆压力的扩散，即浆液的最大扩散范围.
由于篇幅所限，针对本文所开发的八节点五自由度四边形Serendipity单元的验证可参见文献[22]，其中包括利用本文所开发的模型对室内轴对称冲刷实验，以及一维注浆试验的模拟对比与验证.
2. 非饱和地层中的注浆渗透扩散模拟
2.1 算例描述
本节以一典型的非饱和地层注浆为例，着重展示多物理场方法对注浆过程的模拟，主要分析浆液水灰比、注浆压力、土体初始干密度以及土体初始含水率对粉砂地层注浆效果的影响.
2.1.1 数值模型
注浆土层的轴对称截面如图3（a）所示，截面为3 m （长） × 5 m （高）的矩形，注浆孔深3 m，直径0.24 m，注浆区域高度1.5 m，注浆时间30 min. 图3（b）为轴对称注浆问题的三维示意. 本次分析采用5 062个自开发的八节点五自由度四边形轴对称单元模拟土体，并对注浆孔附近的网格进行加密. 模型的竖向边界面为平面对称边界（滚轴边界），模型的底部边界为全固定边界.
2.1.2 材料性质
本次分析中的注浆参数及计算工况汇总如表1所示. 其中，不同水灰比的水泥浆液黏度时变公式根据文献[16]确定.
表 1 数值模拟主要试验变量Table 1. Main test variables for numerical simulation工况 注浆压力/MPa A 黏度时变[16]/
（MPa·s）初始干密
度/（g·cm−3）初始含水率/% 1 1.0 0.5 $ 407.41{{\mathrm{e}}^{0.0005t}} $ 1.6 10 2 1.0 0.7 $ 43.864{{\mathrm{e}}^{0.0088t}} $ 1.6 10 3 1.0 0.8 $ 28.183{{\mathrm{e}}^{0.0114t}} $ 1.6 10 4 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 5 1.0 1.5 $ 10.867{{\mathrm{e}}^{0.0017t}} $ 1.6 10 6 1.0 3.0 $ 6.8812{{\mathrm{e}}^{0.0009t}} $ 1.6 10 7 0.3 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 8 0.6 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 9 0.9 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 10 1.2 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 11 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.5 10 12 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.4 10 13 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.3 10 14 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 5 15 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 15 16 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 20 水泥浆的密度可按式（28）计算.
$$ \rho_{\mathrm{g}}=\frac{1+A}{1/\rho_{\mathrm{c}}+A/\rho_{\mathrm{w}}}\text{，} $$ (28) 式中：$ {\rho _{\mathrm{c}}} $为水泥密度，取3.1 g/cm3；A为水灰比.
根据式（28）确定的不同水灰比浆液密度如表2所示.
表 2 不同水灰比下的水泥浆液密度Table 2. Density of grout with different water-cement ratiosA 0.5 0.7 0.8 1.0 1.5 3.0 水泥浆密度/
（g·cm−3）1.823 1.662 1.603 1.512 1.372 1.204 本次着重分析浆液的渗透扩散特性，因缺乏土体力学特性数据，因此，土体采用线弹性模型. 土体力学参数根据文献[23]中的粉质砂土确定：弹性模量为7 MPa，泊松比为0.3，渗透系数为7.27 × 10−6 m/s，初始孔隙率及初始饱和度可分别根据初始干密度及初始含水率确定，不同孔隙率和饱和度对土体渗透系数的影响根据式（14）确定，粉质砂土的水土特征见表3所示[24].
表 3 粉质砂土的水土特征[23]Table 3. Soil and water characteristic of silty sand[23]饱和度/% 100 94.9647 78.3324 28.0630 23.5915 20.0685 18.6214 17.3607 基质吸力/kPa 0 15.2868 19.6129 34.2297 53.6164 91.5849 156.0620 219.9930 2.2 计算结果
现分别从浆液渗流扩散模式及影响因素展开分析各工况下的注浆模拟结果.
2.2.1 典型浆液扩散云图
典型的浆液扩散模式如图4所示（工况6）. 图中：r为水平扩散半径，h竖向扩散深度. 可以看出：浆液在土体中呈椭球形分布，浆液的扩散区域可分为填充区域和过渡区域；填充区域靠近注浆管壁，且孔隙内的浆液浓度为1.0，即该区域内的孔隙液体全部为浆液；随着扩散距离的增加，浆液浓度在过渡区域内逐渐降低至0. 图5所示为工况6中土体孔隙的液体饱和度分布. 可以看出：饱和区域靠近注浆管壁，随着扩散距离的增加，土体的饱和度逐渐减小至初始饱和度. 通过对比图4和图5可以发现：注浆后的饱和区域范围大于浆液填充区域，这是由于浆液扩散对初始孔隙水及孔隙气体产生驱替作用，使浆液填充区域内的孔隙水被驱至浆液填充区域外，因此，浆液填充区域外土体的饱和度也相应增大，同时，饱和度增大区域的范围也大于浆液的扩散范围. 图6为工况6中土体孔隙压力分布，随着远离注浆管壁，孔隙压力逐渐减小，其分布范围与图4浆液填充区域相当.
2.2.2 浆液扩散范围的影响因素
1） 水灰比
图7为工况1～6注浆30 min后，不同A下，浆液浓度沿水平和竖向测量路径的分布曲线. 图8为不同水灰比浆液的最大扩散范围. 结果发现：在A<1.0时，浆液的最大扩散半径（r）和最大扩散深度（h）基本相同，浆液集中在注浆面附近0.1～0.2 m内；当A≥1.0后，浆液流动性显著增加，浆液的扩散范围随水灰比增大有显著提升；在A=3.0时，浆液的扩散距离达到0.5～0.7 m，且r较h增长更显著.
图 7 不同水灰比浆液的浓度分布Figure 7. Grout concentration distribution under different water cement ratios
图 8 不同水灰比浆液的最大扩散范围Figure 8. The maximum diffusion range of grout with different water cement ratio2） 注浆时间
以注浆压力1 MPa，水灰比1.0 （即工况4）为例，其沿水平测量路径的浆液浓度分布随注浆时间的变化曲线如图9所示. 由图9可知：随着注浆时间的增加，浆液扩散范围逐渐增大，注浆填充区域向土体内部扩展；但由于浆液黏度随时间增加，流动能力逐渐减弱，因此，扩散范围趋于稳定. 在本次分析中，注浆约15 min后浆液达到最大扩散范围.
图 9 水平路径上的浆液浓度分布随时间变化Figure 9. Grout concentration distribution along horizontal path at different time不同水灰比下，浆液的r随注浆时间的变化规律如图10所示. 可以看出：浆液的扩散半径随注浆时间的增长而增加，当A<1.0时，浆液初始黏度增长更为迅速，在注浆15 min即达到最大扩散半径，而后保持恒定；A=1.0时，浆液扩散稳定时间略有增长，约在15～20 min时达到最大扩散范围；当A>1.0后，在注浆25 min后仍有较大的扩散范围增长；高水灰比的浆液更容易渗流扩散，但浆液达到最大扩散范围的时间也会随之增加.
2.2.3 浆液扩散范围估算
为分析不同变量（水灰比、注浆压力、干密度、含水率）的增减幅度对浆液扩散距离的影响，对4种影响因素进行归一化处理：
$$ \left(x_m-\min x_m\right)/\left(\max x_m-\min x_m\right)\text{，} $$ (29) 式中：xm为4种因素归一化数据组中的第m个数据.
不同水灰比（工况1～6）、注浆压力（工况7～10）、土体干密度（工况4、11～13）和含水率（工况4、14～16）下的r和h如图11所示. 针对各因素进行单一变量拟合，所得表达式见表4.
图 11 浆液扩散距离随水灰比、注浆压力、干密度和含水率变化的拟合曲线Figure 11. Fitting curves of grout penetration distance with water-cement ratio, grouting pressure, dry density, and water content表 4 各因素对浆液扩散距离影响的拟合曲线方程Table 4. Fitting equations of grout penetration distance with various factors因素 最大扩散半径 最大扩散深度 水灰比 $ \begin{array}{l} y_1 = {{ - }}0.16 + 0.55x_1 -\\\quad 0.09{x_1^2}\end{array} $ $ \;y_1 = {{ - }}0.09 + 0.44x_1 - 0.08{x_1^2} $ 注浆
压力$ \begin{array}{l} y_1 = 0.12 + 0.55x_1 -\\\quad 0.22{x_1^2}\end{array} $ $ y_1 = 0.16 + 0.38x_1 - 0.14{x_1^2} $ 干密度 $ y_1 = 1.21{{ - }}0.54x_1 $ $ y_1 = 0.93{{ - }}0.37x_1 $ 含水率 $ y_1 = 0.35{{ - 0}}{\text{.01}}x_1 $ $ y_1 = 0.34{{ - 0}}{\text{.15}}x_1 $ 从图11可以看出，各工况的扩散半径增长幅度略小于扩散深度. 其中：浆液扩散范围受水灰比影响最明显，极限扩散范围最大可相差约0.51 m，在水灰比从1.0变化到3.0时增长幅度最大；受注浆压力和干密度对扩散范围的影响其次，其各自极限扩散范围最大相差约0.18、0.16 m；受土体含水率的影响最小，浆液扩散距离几乎没有变化.
数值计算结果表明浆液扩散范围受干密度和含水率的影响相对较小，并近似呈线性变化. 根据式（14）可知，由于在含水率相同时，干密度的增加一方面减小了孔隙率，对渗透系数有降低作用；另一方面，孔隙率的减小同时提高了土体饱和度，会导致渗透系数的增大，两者的共同影响导致浆液的渗流扩散范围变化较小. 同理，在干密度相同的情况下，含水率改变引起的孔隙率和土体饱和度变化的共同作用导致浆液的渗流扩散范围变化不大.
2.3 注浆压力对土层的挤密作用
灌浆土体受注浆压力和浆液渗流的作用会发生挤压变形，在本模型中可通过孔隙率来衡量土体固体骨架的体积变形程度. 图12（a）为不同注浆压力下（工况7～10），注浆段浆管壁监测点处（如图3（b）所示）土体孔隙率随时间的变化. 不同注浆压力下，孔隙率随注浆时间的变化大致可分为2个阶段：初始1～2 min内，浆液尚未扩散至监测点所处土体位置，土体主要在注浆压力的力学作用下产生压缩变形，孔隙率在较短时间内迅速减小；此后，监测点所在土体受浆液渗流的影响，土体孔隙因孔隙压力作用而扩张，故而孔隙率回升，在大约1 200 s后，浆液扩散达到稳定，孔隙率基本不再变化. 与此同时，孔隙率的变化幅度及变化率随注浆压力的增大而增大.
图 12 注浆土体的孔隙率变化及分布特征Figure 12. Porosity variation and distribution characteristics of soil during grouting of soil during grouting图12（b）为不同注浆压力下注浆30 min后，孔隙率沿水平监测路径（如图3（b）所示）的分布. 可以看出：随距管壁距离的增大，土体孔隙率呈现出先减小后增大的规律，这是因为靠近注浆管壁的土体同时受到注浆压力的挤压和孔隙压力的支撑作用；随着远离注浆管壁，孔隙压力逐渐减小，注浆压力的挤密作用占据优势，土体产生压缩变形，孔隙率降低；压密区域外的土体孔隙率随着距离的增加而增大；与此同时，土体孔隙率的降低幅度及最小孔隙率处距管壁的距离随着注浆压力的增大而增大，即注浆压力越大，挤密区域越大.
基于上述结果，本文所提出的多物理场耦合模型和数值方法可以很好地重现注浆过程中浆液的渗透、孔隙水的驱替、孔隙压力的传递与消散以及注浆压力对周围土体的挤密等特性.
3. 结 论
1） 提出了一种计算非饱和地层注浆渗透扩散的固-液-气多物理场耦合模拟方法，通过ABAQUS的二次开发，建立一种新型的八节点五自由度四边形轴对称Serendipity单元，实现对注浆过程中的浆液浓度分布、孔隙压力及土体孔隙率的数值求解.
2） 分析浆液水灰比、注浆压力、土体初始干密度以及土体初始含水率对粉砂地层注浆效果的影响. 浆液扩散范围受水灰比影响最明显，受注浆压力和干密度影响其次，受含水率影响最小. 改变干密度或含水率，将引起孔隙率和土体饱和度变化，其共同作用导致浆液的渗流扩散范围变化不显著.
3） A<1.0时，浆液扩散范围较小，不同水灰比下浆液最大扩散范围相近，且达到最大范围的扩散时间相对较短；A>1.0的浆液更容易渗流扩散，但浆液达到最大扩散范围的时间也会随之增加.
4） 注浆管壁周围会形成挤密区域. 在时间尺度上，浆液扩散区域内的土体孔隙率随注浆时间的增加，先减小后增大；在空间尺度上，由于靠近注浆管壁的土体单元同时受到注浆压力的挤压和孔隙压力的支撑作用，随着远离注浆管壁，土体孔隙率在挤密区域内逐渐减小，在挤密区域外逐渐恢复，且挤密区域随着注浆压力的增加而增大.
需要注意的是，本文所提出的数值方法是一个模块化的数值模拟框架，可根据实际问题以及分析精度要求，继续植入不同的本构模型、水土特征曲线以及渗透系数矩阵等. 本文算例着重展示文中所提数值方法的计算效果，并对注浆问题进行定性分析. 同时，由于缺乏土体力学特性数据，土体采用了线弹性模型. 这并不影响本文中所获得的规律性结论.
-
图 1 ABAQUS调用UEL子程序计算流程
Figure 1. Calculation flow chart of ABAQUS calling UEL subroutine
图 2 新型五自由度轴对称单元示意
Figure 2. Illustration of novel five-degree-of-freedom axisymmetric element
图 7 不同水灰比浆液的浓度分布
Figure 7. Grout concentration distribution under different water cement ratios
图 8 不同水灰比浆液的最大扩散范围
Figure 8. The maximum diffusion range of grout with different water cement ratio
图 9 水平路径上的浆液浓度分布随时间变化
Figure 9. Grout concentration distribution along horizontal path at different time
图 11 浆液扩散距离随水灰比、注浆压力、干密度和含水率变化的拟合曲线
Figure 11. Fitting curves of grout penetration distance with water-cement ratio, grouting pressure, dry density, and water content
图 12 注浆土体的孔隙率变化及分布特征
Figure 12. Porosity variation and distribution characteristics of soil during grouting of soil during grouting
表 1 数值模拟主要试验变量
Table 1. Main test variables for numerical simulation
工况 注浆压力/MPa A 黏度时变[16]/
（MPa·s）初始干密
度/（g·cm−3）初始含水率/% 1 1.0 0.5 $ 407.41{{\mathrm{e}}^{0.0005t}} $ 1.6 10 2 1.0 0.7 $ 43.864{{\mathrm{e}}^{0.0088t}} $ 1.6 10 3 1.0 0.8 $ 28.183{{\mathrm{e}}^{0.0114t}} $ 1.6 10 4 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 5 1.0 1.5 $ 10.867{{\mathrm{e}}^{0.0017t}} $ 1.6 10 6 1.0 3.0 $ 6.8812{{\mathrm{e}}^{0.0009t}} $ 1.6 10 7 0.3 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 8 0.6 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 9 0.9 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 10 1.2 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 10 11 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.5 10 12 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.4 10 13 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.3 10 14 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 5 15 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 15 16 1.0 1.0 $ 15.632{{\mathrm{e}}^{0.0021t}} $ 1.6 20 
下载: 导出CSV
表 2 不同水灰比下的水泥浆液密度
Table 2. Density of grout with different water-cement ratios
A 0.5 0.7 0.8 1.0 1.5 3.0 水泥浆密度/
（g·cm−3）1.823 1.662 1.603 1.512 1.372 1.204
下载: 导出CSV
表 3 粉质砂土的水土特征[23]
Table 3. Soil and water characteristic of silty sand[23]
饱和度/% 100 94.9647 78.3324 28.0630 23.5915 20.0685 18.6214 17.3607 基质吸力/kPa 0 15.2868 19.6129 34.2297 53.6164 91.5849 156.0620 219.9930
下载: 导出CSV
表 4 各因素对浆液扩散距离影响的拟合曲线方程
Table 4. Fitting equations of grout penetration distance with various factors
因素 最大扩散半径 最大扩散深度 水灰比 $ \begin{array}{l} y_1 = {{ - }}0.16 + 0.55x_1 -\\\quad 0.09{x_1^2}\end{array} $ $ \;y_1 = {{ - }}0.09 + 0.44x_1 - 0.08{x_1^2} $ 注浆
压力$ \begin{array}{l} y_1 = 0.12 + 0.55x_1 -\\\quad 0.22{x_1^2}\end{array} $ $ y_1 = 0.16 + 0.38x_1 - 0.14{x_1^2} $ 干密度 $ y_1 = 1.21{{ - }}0.54x_1 $ $ y_1 = 0.93{{ - }}0.37x_1 $ 含水率 $ y_1 = 0.35{{ - 0}}{\text{.01}}x_1 $ $ y_1 = 0.34{{ - 0}}{\text{.15}}x_1 $
下载: 导出CSV
-
[1] 李天胜,何川,方砚兵,等. 基于围岩变形失效的隧道结构可靠度设计方法[J]. 西南交通大学学报,2023,58(3): 613-621. doi: 10.3969/j.issn.1000-7598.2009.05.020LI Tiansheng, HE Chuan, FANG Yanbing, et al. Reliability-based design method of tunnel structures based on deformation failure of surrounding rock[J]. Journal of Southwest Jiaotong University, 2023, 58(3): 613-621. doi: 10.3969/j.issn.1000-7598.2009.05.020 [2] 潘钦锋,张丙强,黄志斌. 隧道下穿诱发既有管道-土体非协调变形解析研究[J]. 西南交通大学学报,2024,59(3): 637-645.PAN Qinfeng, ZHANG Bingqiang, HUANG Zhibin. Analytical study for uncoordinated deformation of existing pipeline and soil induced by tunnel undercrossing[J]. Journal of Southwest Jiaotong University, 2024, 59(3): 637-645. [3] 秦鹏飞,钟宏伟,刘坚,等. 考虑浆土应力耦合作用的劈裂注浆机理分析[J]. 西南交通大学学报,2023,58(3):584-591.QIN Pengfei, ZHONG Hongwei, LIU Jian, et al. Analysis of split grouting mechanism considering coupling effect of slurry and soil stress[J]. Journal of Southwest Jiaotong University, 2023, 58(3): 584-591. [4] 刘奇,陈卫忠,袁敬强,等. 基于渗流-侵蚀理论的岩溶充填介质注浆加固效果评价[J]. 岩石力学与工程学报,2020,39(3):572-580.LIU Qi, CHEN Weizhong, YUAN Jingqiang, et al. Evaluation of grouting reinforcement effect for karst filling medium based on seepage-erosion theory[J]. Chinese Journal of Rock Mechanics and Engineering, 2020, 39(3): 572-580. [5] 马明杰,杨新安,周建,等. 非达西渗流作用下饱和黏土压密注浆扩孔理论分析[J]. 哈尔滨工业大学学报,2022,54(3):32-40.MA Mingjie, YANG Xin’an, ZHOU Jian, et al. Theoretical analysis of cavity expansion for compaction grouting in saturated clay under non-Darcy seepage[J]. Journal of Harbin Institute of Technology, 2022, 54(3): 32-40. [6] YIN Z Y, YANG J, LAOUAFA F, et al. A framework for coupled hydro-mechanical continuous modelling of gap-graded granular soils subjected to suffusion[J]. European Journal of Environmental and Civil Engineering, 2020, 27(1): 1-22. [7] YANG J, YIN Z Y, LAOUAFA F, et al. Analysis of suffusion in cohesionless soils with randomly distributed porosity and fines content[J]. Computers and Geotechnics, 2019, 111: 157-171. doi: 10.1016/j.compgeo.2019.03.011 [8] YANG J, YIN Z Y, LAOUAFA F, et al. Hydromechanical modeling of granular soils considering internal erosion[J]. Canadian Geotechnical Journal, 2020, 57(2): 157-172. doi: 10.1139/cgj-2018-0653 [9] YANG J, YIN Z Y, LAOUAFA F, et al. Three-dimensional hydromechanical modeling of internal erosion in dike-on-foundation[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2020, 44(8): 1200-1218. doi: 10.1002/nag.3057 [10] RAATS P A C. Applications of the theory of mixtures in soil science[J]. Mathematical Modelling, 1987, 9(12): 849-856. doi: 10.1016/0270-0255(87)90003-0 [11] DE BOER R. Contemporary progress in porous media theory[J]. Applied Mechanics Reviews, 2000, 53(12): 323-370 doi: 10.1115/1.3097333 [12] UZUOKA R, BORJA R I. Dynamics of unsaturated poroelastic solids at finite strain[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2012, 36(13): 1535-1573. doi: 10.1002/nag.1061 [13] ZHOU Z L, ZANG H Z, WANG S Y, et al. Filtration behaviour of cement-based grout in porous media[J]. Transport in Porous Media, 2018, 125(3):435-463. doi: 10.1007/s11242-018-1127-x [14] LIU X X, SHEN S L, XU Y S, et al. A diffusion model for backfill grout behind shield tunnel lining[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2021, 45(4): 457-477. doi: 10.1002/nag.3164 [15] REVIL A, CATHLES L. Permeability of shaly sands[J]. Water Resources Research, 1999, 35(3): 651-662. doi: 10.1029/98WR02700 [16] 雷进生. 碎石土地基注浆加固力学行为研究[D]. 武汉:中国地质大学,2013. [17] LIANG Y, ZHANG J, LAI Z S, et al. Temporal and spatial distribution of the grout pressure and its effects on lining segments during synchronous grouting in shield tunnelling[J]. European Journal of Environmental and Civil Engineering, 2020, 24(1): 79-96. doi: 10.1080/19648189.2017.1364299 [18] SMITH M. ABAQUS/Standard: user’s manual[M]. Providence: [s.n.], 1997. [19] LADYZHENSKAYA O A. The mathematical theory of viscous incompressible flow[M]. New York:Gordon & Breach, 1969. [20] BREZZI F. On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers[J]. Revue Française D’automatique, Informatique, Recherche Opérationnelle Analyse Numérique, 1974, 8: 129-151. [21] 王勖成,邵敏. 有限单元法基本原理和数值方法[M]. 2版. 北京: 清华大学出版社,1997. [22] YANG J. Numerical analyses of the multi-physics problem of sinkholes in the vicinity of a dike or a linear geo-structure[D]. Nantes: École centrale de Nantes, 2019. [23] 林鸿州,李广信,于玉贞,等. 基质吸力对非饱和土抗剪强度的影响[J]. 岩土力学,2007,28(9):1931-1936. doi: 10.3969/j.issn.1000-7598.2007.09.031LIN Hongzhou, LI Guangxin, YU Yuzhen, et al. Influence of matric suction on shear strength behavior of unsaturated soils[J]. Rock and Soil Mechanics, 2007, 28(9): 1931-1936. doi: 10.3969/j.issn.1000-7598.2007.09.031 [24] 李文. 路基沉陷劈裂注浆处治试验研究及力学参数计算[D]. 长沙: 长沙理工大学,2013. 期刊类型引用(1)
1. 勾奇庆，冯皓龙，李文杰，梁斌. 双层采空区垮落带破碎岩注浆充填浆液扩散特性研究. 金属矿山. 2025(02): 189-197 . 
百度学术其他类型引用(0)
-
下载:
下载:
百度学术
点击查看大图
计量
- 文章访问数: 360
- HTML全文浏览量: 73
- PDF下载量: 51
- 被引次数: 1