泥浆侵入土层泥膜状态细观分析与堆积单向性验证

刘 成，杨 平，王海波，潘晓明

(1. 南京林业大学土木工程学院，江苏南京 210037； 2. 深圳市地铁集团有限公司，深圳 518026)

泥水盾构法施工依靠在开挖面上形成的微透水致密泥膜，通过平衡开挖面前方地层的水土压力保证开挖面的稳定。泥膜形成受到多种因素影响，形成机理复杂，泥膜的最终状态是多种类型的组合。

目前，泥膜形成机理和状态演变规律尚无完整的理论，相关理论和试验研究仍滞后于工程实践。泥膜(泥饼)的研究起源于化学工业和矿产工业中过滤处理、废物处理等方面的研究。B.F.Ruth和P.C.Carman等最早进行了过滤介质固-液分离和泥膜形成与增长规律研究，提出了泥膜过滤方程的经典公式。该式引入了较多的假设，无法描述泥膜结构固结变形特征，但该式较为简单，应用较广。C.Tien等[1]根据多相流理论求解分析，给出了考虑泥膜内部颗粒运动解析解，后被白云等[2]引入到泥水盾构泥膜形成模型研究，取得了较好的分析效果。刘成等[3]假设泥膜为饱和重塑土，根据应力状态和土的体积状态的唯一性关系，在泥膜动态形成过程中同时考虑了其固结变形。上述泥膜的过滤过程不考虑泥浆颗粒进入过滤介质，可统称为泥皮型泥膜。另有学者根据细颗粒可进入多孔介质并在多孔介质及表面形成泥浆颗粒填充和表面泥膜堆积，提出了渗透性泥膜模型，分析了泥浆颗粒侵入土层深度的影响因素，提出了适用于高透水性地层中开挖面稳定的分析模型[4-6]。

泥膜形成过程是泥浆固-液分离的过程，也是多孔介质堵塞、填充堆积的过程，泥膜内部结构、渗透系数及压降随着时间增长不断变化。因此，常规的连续介质分析手段无法模拟这个过程，颗粒流离散元法则可以从细观尺度将颗粒与周围介质的作用描述出来，包括颗粒间接触、非接触相互作用，颗粒-流体相互作用，特别是粒径介于10～100 μm之间的颗粒，粒间范德华力不可忽略。K.J.Dong等[7-9]采用颗粒离散元法分析了流速恒定和压力恒定下沉降和滤失过程中泥膜的形成与增长规律，考虑了范德华力对其的影响。C.Eichholz等[10]采用三维离散元法分析了磁场对泥膜形成的影响。这些分析多局限于泥皮型泥膜，与细颗粒进入粗粒土间隙形成其他状态的泥膜有很大的差异。土层堆积和土层间隙是三维问题，传统的过滤理论和一维、二维的离散元无法分析细颗粒进入粗颗粒土间隙的问题。本文基于三维颗粒流离散元程序YADE，编制颗粒-流体相互作用和范德华力作用模型，模拟初始恒速阶段泥浆侵入土层过程和泥膜形成过程，根据在局部区域形成稳定堆积的颗粒数来量化泥膜的稳定堆积状态，进行泥膜状态的判别，并与泥膜状态的相关理论进行比较分析。通过砂-膨润土混合物变水头渗透试验验证了泥浆颗粒在土层中的堆积具有一定的单向性特点。

1 泥膜形成判定与状态划分

1.1 泥膜形成判定

泥膜是在一定的压力作用下由颗粒逐渐堆积形成，随着时间延长泥膜厚度增加，在厚度方向表现出不同程度的压缩性。泥膜针入度特征曲线表明泥膜具有层状结构物理形态，从里到外分别是致密层、密实层、可压缩层以及虚泥膜层[11]。表层的虚泥膜层是一层疏松的呈胶凝状态的泥浆，其强度接近于0，容易被表层流动的泥浆冲蚀，只有在泥浆静止时才能形成，因此不是真正的实泥膜。可压缩层表现一定的强度和可压缩性，可作为分析的泥膜外侧。判别泥膜形成常用两种方法：① 泥浆平均浓度大于泥膜形成临界浓度[12-13]；② 单个颗粒接触数应大于等于3，顶层颗粒接触数可能等于3，颗粒速度接近于0[14]。文中采用后者计算泥膜稳定堆积颗粒的数量，并根据各层堆积比率确定泥膜状态。

1.2 泥膜状态的划分

泥水在掘削面上的渗透形态也叫泥膜状态，可分为3种[15-17]：①泥皮型，即“表面泥膜”，泥膜在开挖面表面完全形成，对应“薄膜模型”。认为泥浆压力作为外力完全有效地作用于不透水的薄膜上，这种情况多发生在黏粒土、粉粒土及细砂土等土层，对应地层的有效间隙L3dmax(泥水最大粒径)，其解决措施是增大泥水的粒径，即在泥水中添加砂粒；③泥皮+渗透带型，相当于“中间状态”，这种情形多发生于砂性地层(中、细)，对应地层的有效间隙L满足条件dmin

(a) 泥皮型 (b) 渗透带型 (c) 泥皮+渗透带型图1 泥膜状态3种类型Fig.1 Three types of slurry membrane state

这种泥膜状态的划分标准简单实用，可以较好地应用于工程实际。但是，泥膜形成受到多种因素的影响，根据泥浆粒径与土层有效间隙的关系确定泥膜状态时忽略了土层条件、泥浆配比和操作压力等重要因素。另外，土层是非均匀体，粒径范围较广，孔隙率或者孔洞的存在导致土层有效间隙难以确定。采用离散元法确定特定土层粒径范围和土体颗粒堆积状态下的泥膜状态可以作为实际复杂影响因素分析的基础。

2 泥浆侵入地层细观分析

2.1 接触模型

泥浆侵入地层过程中，泥浆颗粒与周围土层颗粒和泥浆颗粒间相互作用包括接触法向力、接触摩擦力及非接触范德华力。接触本构模型采用Cundall-Strack本构模型，模型由线弹性压缩阶段和不含黏聚力的摩尔库伦屈服面组成，颗粒间不承受拉应力，模型仅定义了切向接触摩擦，未考虑滚动摩擦，对于粒径介于10～100 μm之间的泥浆颗粒受力分析是合适的。

粒间接触受力由YADE根据接触位移和接触刚度确定。YADE定义了两个接触刚度：法向刚度kn和切向刚度ks。切向刚度ks根据法向刚度kn乘上接触刚度比ks/kn确定，法向刚度kn则用两个长度等于球体半径的弹簧链表示，根据下式计算：

(1)

两个球体之间的法向力F可以表示为：

(2)

2.2 颗粒受力分析

在离散元分析中，一个颗粒的运动包括平动和转动两个部分，由于分析中未考虑滚动摩擦，仅列出颗粒平动运动模型公式如下：

(3)

式中：vi,mi和Fi分别表示颗粒i的运动速度、质量和作用在其上的合力。

泥浆在操作压力作用下进入前方土层，颗粒受到各种力的作用，为简化分析，进行了如下假设：(1)假设流体为一维流动，只考虑流体对颗粒的拖曳力和浮力，不考虑颗粒对流体的拖曳力作用。(2)泥膜形成过程中泥膜内部压降不断增加，渗透系数随时间降低。泥膜形成过程分为恒定流速阶段和恒定压力阶段，这里仅分析第一阶段即初始恒定流速阶段。(3)假设土层含水率较大，忽略毛细力作用。

(4)

图2 颗粒i上的作用力示意图Fig.2 Forces acting on particle i

(7)

范德华力是非接触力，使得颗粒处于团聚的状态。在离散元计算中，仅考虑范德华力引力部分，不考虑斥力作用。范德华力与Hamaker常数Ha成正比，与粒间间隔h平方成反比。当粒间间隔无限接近于0或者已发生接触，范德华力将接近于无穷大，这将导致离散元计算过程发生奇异。为了克服数值奇异，可以引入一个“截止距离hmin”[7]，通常取hmin为1 nm，h小于1 nm按1 nm计算。

2.3 离散元模型介绍

(a) 立体图 (b) 俯视图图3 最密堆积(六方密堆积)示意图Fig.3 Schematic diagrams of the closest packing (hexagonal close packing)

泥浆初始浓度为1-ε0，孔隙率为ε0，对应空间的泥浆颗粒数量为N0。计算中每隔0.1 s在顶部0.1 s×vf0高度空间内加入N0/10个颗粒，如此循环计算至预定时间，分析各阶段泥浆颗粒的堆填状态。

泥膜形成一般分2个阶段，即初始恒定流速阶段和后期的恒压过滤阶段。当土层压降和土层表面泥皮压降之和等于施加的最大压降时，进入恒压过滤阶段。这里仅研究第1阶段即初始恒定流速阶段。

模型中箱体由光滑刚性墙体组成，摩擦系数为0，刚度为球体的1 000倍；土层和泥浆主要细观参数见表1。由于泥浆小颗粒受流体作用在土层孔隙间流动和填充，土体变形量相对较小，在参数标定的基础上降低颗粒弹性模量约100～1 000倍，增加计算效率。

表1 泥浆侵入土层细观参数Tab.1 Microscopic parameters for slurry invading stratum

模型增加了长程引力非接触判断以及流体对每个球体的浮力与拖曳力计算，需遍历颗粒周围非接触极限距离范围内的所有颗粒来确定范德华力，耗时较长，分析时间随颗粒数量增加近似呈指数增长。

3 泥膜状态模拟结果分析

泥浆颗粒受流体拖曳力作用在土层中迁移，与土层产生接触和继续运动具有一定随机性。与重力沉降堆积和不考虑泥浆侵入地层的过滤过程不同，颗粒堆积厚度并非均匀增长，而是局部形成堆积并向外扩散最终堵塞土层间隙。堵塞也并非首先发生在土层表面，堆积的增长模式可能是横向扩展，也可能是纵向延伸。

根据泥膜形成判定要求，选择各阶段具有稳定堆积的颗粒数量来分析泥膜的增长过程，根据各层堆积比率来判断泥膜的状态。

3.1 范德华力对泥膜状态影响分析

模型中土层颗粒的位移受到边界条件约束，位移几乎为0。范德华力是长程引力，而土层粒径较大，对泥浆颗粒的引力较大，粒间接触数和配位数Z(各颗粒配位数之和除以颗粒总数)随着范德华力增大而增大，如图4所示。粒间接触数和配位数增大表明堆积变得更紧密。因此，与K.J.Dong等[8]的研究相比较可知，自重沉降堆积形成的泥膜孔隙率受范德华力影响较大，泥膜结构变疏松的程度较大；不考虑泥浆侵入地层的过滤过程形成的泥膜孔隙率受范德华力影响较小，泥膜结构变疏松的程度较小；而考虑泥浆侵入地层的过滤过程中颗粒堆积变得更紧密。产生这种变化的根本原因是3种过程中作用在泥浆颗粒上的颗粒-流体间作用力和不同颗粒-颗粒间作用力的差异。K.J.Dong等指出颗粒-流体间向下的作用力加强了重力的挤密作用而降低了范德华力的疏松作用。这里仅分析了泥浆侵入土层的第1阶段即恒速阶段，拖曳力的作用相对较大；而土层颗粒相对位置固定，增加了吸附在上面的尺寸较小的泥浆颗粒的范德华引力作用，因此粒间接触数和配位数增大表明堆积变得更紧密。

图5给出了土层表面泥皮(图例top，线型虚线)和整个分布区域(图例all，线型实线)的颗粒数量随时间延长的变化规律。根据颗粒运动速度和与周围颗粒接触数量判定是否形成稳定堆积。图5为分布在各层已形成稳定堆积的颗粒的统计结果。

图4 范德华力与配位数的关系 图5 范德华力与泥浆颗粒在土层中分布规律的关系 Fig.4 Relationships between Van der Waals force and coordination number Fig.5 Relationships between Van der Waals force and distribution regularities of slurry particles in soil

泥浆颗粒经过土层颗粒附近时，在土层颗粒的范德华力作用下降低了颗粒移动速度，局部区域泥浆颗粒聚集，形成堆积；随着堆积结构的形成，透水性降低，流体对泥浆颗粒的托曳力降低，这也是造成堆积速度增加的原因。因此，范德华力对颗粒受力的贡献比例随时间变化，其影响也发生变化。图5表明范德华力对堆积速度增加的影响呈非线性，范德华力增加至一定值时，虽对土层表面泥膜(top)的增长幅度影响较大，但对颗粒堆积总数(all)的增长幅度的影响却在降低。

形成稳定堆积的颗粒数量占颗粒总数比例较大，大约在80%至95%之间，随着范德华力的增大而增大，但各层分布的比例变化规律不甚明显。

(a) Ha=1×10-21J (b) Ha=20×10-21J图6 不同范德华力作用下空间力链结构(t=18 s)Fig.6 Spatial force chain structure under the action of different Van der Waals forces (t=18 s)

图6给出了范德华力对空间力链结构的影响，图中球体被缩小至30%以便观察力链结构，这里仅考虑作用在泥浆颗粒上的力链，未考虑墙体与颗粒的力链。图6表明：(1)形成的力链基本分成4组：① 泥浆颗粒受到2个以上土层颗粒作用；② 泥浆颗粒受到1个土层颗粒和其他泥浆颗粒作用；③ 泥浆颗粒受到其他泥浆颗粒作用；④ 泥浆颗粒与其他颗粒不发生作用。其中①组的力链最粗，粒间作用力较大，其他组依次递减。(2)力链的形成往往是关联的，力链上下传递有一定的连续性，堆积之间往往有泥浆颗粒作为媒介形成空间结构，证明了稳定堆积的形成往往从某个堆积位置向周围扩散。(3)颗粒堆积具有一定的黏结单向性，导致渗透系数和强度等表现出一定的单向性，如渗流从上至下作用堆积趋于稳定，而从下至上作用堆积易于破坏。因为向下渗流粒间的作用力最终传递给稳定的土层颗粒，而向上渗流粒间的作用力最终传递给连接力较低的泥浆颗粒。由于渗透系数随着堆积结构的改变而变化，后面将利用砂-膨润土混合物的变水头渗透试验进行验证。因此颗粒堆积单向性是泥浆在土层中填充堆积的一个重要特点。

上述分析表明，其他条件相同时范德华力的变化使得泥膜形成堆积的速度和各层颗粒数量占总数的比例产生量的变化；范德华力增大使得泥膜状态向形成泥皮方向发展。不同流体介质中，Hamaker常数Ha不同，因此，改变泥浆配比增大Hamaker常数，将是在土层表面更快形成泥皮的一种有效手段。目前这方面的理论和试验研究还有待继续发展和深入。

3.2 粒径比对泥膜状态影响分析

根据上述划分选择5种情形(case 1～case 5)进行比较分析，case 1～case 4土层粒径为1 mm，粒径比分别为0.051 5，0.075 0，0.100 0，0.154 7，case 5土层粒径为0.5 mm，粒径比为0.154 7。粒径比为0.051 5时，泥水最大粒径dpmax=0.051 5 mm<0.051 5ds，对应逸泥状态；粒径比为0.051 5，0.075 0和0.100 0时，满足dpmin<0.154 7ds<3dpmax，对应泥浆渗入土层一定深度的状态；粒径比为0.154 7时，泥水最小粒径dpmin≥0.154 7ds，对应泥皮状态。

图7分别给出了不同粒径比下泥浆颗粒在土层中形成稳定堆积的分布规律曲线，图示表明粒径比对泥浆颗粒分布起关键作用。case 1粒径比为0.051 5时，泥浆颗粒易于穿过土层分析区域，仅在两个土颗粒接触的位置形成堆积，最终未形成有效的堵塞，此种泥膜状态属于逸泥状态。粒径比为0.075 0和0.100 0时，泥浆颗粒渗入土层一定深度后开始形成堆积，最终形成有效的堵塞。堆积向土层表面发展在土层表面形成泥皮，此种泥膜状态属于泥浆渗入土层一定深度的状态。粒径比为0.154 7时，泥浆颗粒相对较大容易形成堵塞，经过数秒，泥浆颗粒无法进入土层区域，且进入的泥浆颗粒未被流体带走，此种泥膜状态属于泥皮状态，其中t=6.5 s时不同粒径比的泥浆颗粒分布和堆积状态如图8所示。

上述分析表明，离散元分析的结果与张凤祥等对泥膜状态的划分标准基本一致，不过这里为了控制土层有效间隙均匀性，土层颗粒采用了六方密堆积，孔隙率较低。泥浆浓度为0.001 5，较实际工程采用的泥浆配比浓度低，而土层颗粒和泥浆颗粒采用等值粒径，对分析结果产生的影响将在后续分析中作进一步研究。

(a) case 1 (b) case 2 (c) case 4 (d) case 5 图7 粒径比与泥浆颗粒在土层中分布规律的关系 图8 泥浆颗粒分布和堆积状态(t=6.5 s)Fig.7 Relationships between partical diameter ratio and distribution regularities of slurry particles in soil Fig.8 Distribution and packing state of slurry particles(t=6.5 s)

4 堆积单向性试验验证

4.1 试验概况

室内试验直接验证泥浆颗粒在砂土中形成的堆积具有一定的单向性比较困难，一方面是室内试验(如变水头渗透试验)土样内部的颗粒堆积状态无法直接观测，而对土样切片并进行图像处理的常规方法由于对土样的扰动量过大造成土样结构破坏也不适用；另一方面是难以给出定量的数据对单向性进行验证。泥浆颗粒在土样间隙中的填充堆积必然引起渗透系数随着时间变化。这里采用砂-膨润土变水头渗透试验对渗流方向改变引起的渗透系数改变进行定量分析，从而间接验证了泥浆颗粒在土样中堆积的发展过程和堆积单向性特点。

采用TST-55渗透仪进行变水头渗透试验，砂-膨润土混合物由粒径为0.25～0.50 mm砂和南京汤山膨润土有限公司生产的GCL专用膨润土按一定比例拌合而成。掺砂率分别为80%和85%，干法制样，控制干土总质量为197.83 g，进行3组平行试验。装样后，抽真空饱和，连接变水头管，记录渗透系数随时间的变化数据，4 d后数据已基本稳定，第11天将仪器倒置，调整进水口和出水口，保证进水口在下出水口在上，记录数据。

4.2 试验结果分析

图9为砂-膨润土渗透系数与试验时间的关系曲线，表明随着膨润土含量的增大，混合物的渗透系数降低。随着渗透次数的增加和试验时间的延长，膨润土在砂土中迁移并逐渐形成类似于泥膜的结构，导致渗透系数下降，稳定后渗透系数约为4.4×10-8cm/s。仪器倒置后，泥膜结构受到扰动产生不完全的破坏，因此渗透系数上升，但不超过仪器正放时的初始值。随着渗流的进行，新的泥膜结构逐渐形成，渗透系数逐渐降低。从渗透系数的变化规律可知，膨润土在混合物中形成的堆积具有一定的单向性，间接验证了泥浆颗粒在砂土中形成的堆积具有一定的单向性。

图9 砂-膨润土渗透系数与试验时间的关系Fig.9 Variation of permeability of sand-bentonite mixtures with testing time

由于渗透试验中土样和泥浆参数、水力条件和YADE模型尚存在一定的差异，仅能从变化规律上进行间接对比分析，因此泥浆侵入土层堆积单向性和空间结构特点有待进一步试验验证。

5 结 语

在三维颗粒流离散元程序YADE基础上编制了流体与颗粒相互作用以及范德华力作用模型，通过离散元模拟，重点分析了泥浆侵入土层初始恒速阶段泥膜的范德华力和粒径比对泥膜状态的影响；通过砂-膨润土混合物变水头渗透试验对泥浆颗粒在土层中形成的堆积单向性进行了间接试验验证，得到如下结论：

(1)与自重沉降堆积和不考虑泥浆侵入的过滤过程受力特点不同，范德华力增大使得泥浆颗粒堆积更为紧密。

(2)范德华力是影响泥膜形成速度和泥浆颗粒堆积状态的重要因素；范德华力增大使得泥膜状态向形成泥皮方向发展。通过理论和试验手段改变泥浆配比增大Hamaker常数，将是在土层表面更快形成泥皮的一种有效手段。

(3)泥浆颗粒在土层中的堆积具有一定的单向性和连续性规律，这和泥浆颗粒与周围土层颗粒和泥浆颗粒形成的空间力链结构有关。

(4)粒径比是泥膜状态的决定因素，粒径比分析验证了泥膜状态划分标准相关理论。

(5)变水头试验中渗流方向改变了膨润土在砂土中形成的泥膜结构；渗透系数的定量分析表明泥浆颗粒在土层中形成的堆积具有一定的单向性。

泥膜形成机理复杂，影响因素较多，篇幅所限文中仅对部分因素进行分析，尚无法揭示泥膜形成过程的全貌，而从泥膜形成的恒速阶段向恒压阶段转变将是今后继续研究的重点。

参 考 文 献：

[1]TIEN C, BAI R. An assessment of the conventional cake filtration theory[J]. Chemical Engineering Science, 2003, 58(7)： 1323-1336.

[2]白云， 孔祥鹏， 廖少明. 泥水盾构泥膜动态形成机制研究[J]. 岩土力学， 2010, 31(增2)： 19-24. (BAI Yun, KONG Xiang-peng, LIAO Shao-ming. Research on dynamic formation mechanism of slurry membrane for slurry shield[J]. Rock and Soil Mechanics, 2010, 31(Suppl2)： 19-24. (in Chinese))

[3]刘成， 孙钧， 赵志峰， 等. 泥水盾构泥膜形成二维理论分析[J]. 岩土力学, 2013, 34(6)： 1593-1597. (LIU Cheng, SUN Jun, ZHAO Zhi-feng, et al. Two-dimensional theoretical analysis of slurry membrane formation process in slurry shield[J]. Rock and Soil Mechanics, 2013, 34(6)： 1593-1597. (in Chinese))

[4]BROERE W, VAN TOL A F. Influence of infiltration and groundwater flow on tunnel face stability[C]∥Geotechnical Aspects of Underground Construction in Soft Ground, Japan: Kusakabe, 2000： 339-344.

[5]KIRSCH A. Experimental investigation of the face stability of shallow tunnels in sand[J]. Acta Geotechnica， 2010, 5(1)： 43-62.

[6]AHMED M. Investigation of tunnel face stability and ground movements using transparent soil models[D]. New York： Polytechnic Institute of New York University, 2011.

[7]DONG K J, YANG R Y, ZOU R P, et al. Simulation of the cake formation and growth in sedimentation and filtration[C]∥3rd International Conference on CFD in Minerals and Process Industries, Melbourne, Australia. 2003： 189-194.

[8]DONG K J, ZOU R P, YANG R Y, et al. DEM simulation of cake formation in sedimentation and filtration[J]. Minerals Engineering, 2009, 22(11)： 921-930.

[9]DONG K J, YANG R Y, ZOU R P, et al. Settling of particles in liquids： Effects of material properties[J]. AIChE Journal, 2012, 58(5)： 1409-1421.

[10]EICHHOLZ C, NIRSCHL H, CHEN F, et al. DEM-simulation of the magnetic field enhanced cake filtration[J]. AIChE Journal, 2012, 58(12)： 3633-3644.

[11]周风山， 王世虎， 李继勇， 等. 泥饼结构物理模型与数学模型研究[J]. 钻井液与完井液, 2003(3)： 7-11. (ZHOU Feng-shan, WANG Shi-hu, LI Ji-yong, et al. Mathematical & physical models for drilling fluid filter cake[J]. Drilling Fluid & Completion Fluid, 2003(3)： 7-11. (in Chinese))

[12]TIEN C, RAMARAO B V. Revisiting the laws of filtration： An assessment of their use in identifying particle retention mechanisms in filtration[J]. Journal of Membrane Science, 2011, 383(1-2)： 17-25.

[13]BAI R, TIEN C. Further work on cake filtration analysis[J]. Chemical Engineering Science, 2005, 60(2)： 301-313.

[14]NI L A, YU A B, LU G Q, et al. Simulation of the cake formation and growth in cake filtration[J]. Minerals Engineering, 2006, 19(10)： 1084-1097.

[15]张凤祥， 傅德明， 杨国祥， 等. 盾构隧道施工手册[M]. 北京： 人民交通出版社, 2005. (ZHANG Feng-xiang, FU De-ming, YANG Guo-xiang, et al. Construction manual of shield tunnel[M]. Beijing： China Communications Press, 2005. (in Chinese))

[16]韩晓瑞， 朱伟， 刘泉维， 等. 泥浆性质对泥水盾构开挖面泥膜形成质量影响[J]. 岩土力学, 2008, 29(增1)： 288-292. (HAN Xiao-rui, ZHU Wei, LIU Quan-wei, et al. Influence of slurry property on filter-cake quality on working face of slurry shield[J]. Rock and Soil Mechanics, 2008, 29(Suppl1)： 288-292. (in Chinese))

[17]邓宗伟， 伍振志， 曹浩， 等. 基于流固耦合的泥水盾构隧道施工引发地表变形[J]. 中南大学学报： 自然科学版, 2013(2)： 785-791. (DENG Zong-wei, WU Zhen-zhi, CAO Hao, et al. Surface deformation of slurry shield tunneling using fluid-solid coupling theory[J]. Journal of Central South University (Science and Technology), 2013(2)： 785-791. (in Chinese))