南京江北新区地铁盾构开挖引发地面沉降三维数值模拟预测

摘要：为研究地铁盾构施工引起的地面沉降问题，首先以Biot固结理论为基础，将土体本构关系拓展为黏弹塑性，考虑土体渗透性的动态变化，建立了地铁盾构引发地面沉降的三维全耦合数值模型；然后以南京地铁11号线江北新区马骡圩站至石塘公园站盾构施工区间为例，模拟预测了该盾构区间的沉降量，并完成了模型的校准；最后利用该模型模拟预测了地铁盾构开挖所引起的地面沉降规律，并模拟了盾构拱顶上方土体和盾构底部下方土体的参数变化。结果表明：实测地面沉降量与计算值吻合较好，模型可靠；以模型第二区间段为例，东线轴线处沉降量在52.41～54.52 mm范围内，西线沉降量在53.28～55.60 mm范围内；东线盾构隧道沉降以隧道轴线为中心对称分布，西线隧道开挖后，最终引起累积地面沉降量在轴线之间较大；随着开挖的进行，上方土体孔隙度、渗透系数及泊松比变小，变形模量增大，下方土体变化则相反。

关键词：地面沉降；Biot固结理论；盾构施工；三维全耦合；地铁；隧道

doi：10.13278/j.cnki.jjuese.20230178

中图分类号：P641

文献标志码：A

Supported by the Geological Science and Technology Innovation Project of Jiangsu Province （2022KY01），the Open Fund of Jiangsu Geo-Engineering Environment Intelligent Monitoring Engineering Research Centre （2021ZNJKJJ13） and the National Natural Science Foundation of China （41874014）

Three-Dimensional Numerical Simulation of Ground Settlement Caused by Metro Shield Excavation in Jiangbei New District， Nanjing， China

Xu Chenghua1，He Zhengyu2，Liu Gang1，Shi Wei1，Luo Zujiang2

1. The First Geological Brigade of the Bureau of Geology and Mineral Resources of Jiangsu，Nanjing 210041，China

2. School of Earth Sciences and Engineering，Hohai University，Nanjing 211100，China

Abstract： In order to study the ground settlement caused by subway shield construction， firstly，based on Biot consolidation theory， a visco-elastic-plastic constitutive relation of soil mass and the dynamic change of soil permeability is considered， a three-dimensional fully coupled numerical model is established. Then， taking the shield construction section from Maluowei Station to Shitang Park Station of Jiangbei New District of Nanjing Metro Line 11 as an example， the settlement of the shield section is simulated and predicted， and the model is calibrated. Finally， the model is used to simulate and predict the ground settlement law caused by shield excavation， and the parameter changes of soil above the shield arch and under the shield bottom are simulated. The results show that the measured ground settlement is in good agreement with the calculated values， and the model is reliable. Taking the second section of the model as an example， the settlement of the eastern line is in the range of 52.4154.52 mm， and the settlement of the western line is in the range of 53.2855.60 mm. The settlement of shield tunnel in the eastern line is symmetrically distributed in the center of the tunnel axis， and the cumulative ground settlement of the tunnel in the western line is larger between the axes after tunnel excavation. With the progress of excavation， the porosity， permeability coefficient and Poisson’s ratio of the soil above decreased， and the deformation modulus increased， while the change of the soil below was opposite.

Key words： ground settlement; Bio-consolidation theory; shield construction; three-dimensional full coupling；subway；tunnel

0 引言

地铁因其独特优势在我国的许多城市已成为重要交通工具之一，在常见的地铁隧道开挖修建的方法中，盾构法由于其安全、快速等优点而被广泛应用。然而盾构施工技术会影响施工区域附近地表稳定性，从而引起地面沉降等环境地质问题[1]。

在盾构施工引起地面沉降的预测研究中：Pan等[2]利用Peck公式对天津地铁施工的地面沉降进行了预测；Yu等[3]探讨了隧道开挖过程中巨石引起的地层坍塌和地表沉降规律。利用Biot固结理论对地面沉降的研究中：Adachi等[4]基于Biot的固结理论和Adachi-Oka的应变硬化、应变软化的弹塑性本构模型，对开挖引起的周围地面力学行为进行了有限元耦合分析；Verruijt[5]研究了多孔弹性中的非线性效应；Gutierrez等[6]利用Biot的固结理论对地下地层流体变形和流动的耦合进行了研究。随着计算机技术的发展：Burbey[7]基于Biot三维固结理论对三维含水层的流动和变形建模；Sabetamal等[8]在任意拉格朗日欧拉方法的框架下，提出了一个在Biot固结理论的基础上扩展而来，包含了固体和流体相的动量方程、大变形和非线性材料行为的大变形公式。近些年在对于地面沉降模型的研究中：孟世豪等[9]结合沉降过程中垂向渗透系数与水头的变化规律，建立了变渗透系数模型；袁帅等[10]基于离散化的Theis井流模型以及非线性的分层总和法，构建了孔隙率、渗透系数与压缩模量随固结情况改变的水位与沉降动态计算模型。

以往研究盾构开挖引起的地面沉降所提出的计算模型有很多，但大多模型未考虑地下水渗流场对盾构开挖引起的地面沉降的影响及土体参数的动态变化，不符合实际情况。本文以南京地铁11号线江北新区马骡圩站至石塘公园盾构施工区间为例，建立了盾构开挖和地下水渗流场与地面沉降的三维全耦合模型，将土体本构关系拓展到黏弹塑性，引入土体参数的动态变化，并考虑地下水渗流场对盾构开挖引起的地面沉降的影响，模拟预测了该盾构区间的沉降量，以研究地铁盾构开挖对地面沉降的影响。

1 区域概况

本次模型建立于南京地铁11号线马骡圩站至石塘公园站盾构施工区间，具体位置如图1所示，其纵向长度约670 m，横向宽度40.00 m，地面高程5.00～8.00 m。本区间顶板埋深为9.25 m，底板埋深约为16.25 m。

1.1 地质背景

研究区地貌类型为长江漫滩冲积平原区，第四系广泛分布，为更新世以后的长江和秦淮河沉积。其下伏前第四纪地层均为白垩系上统浦口组，岩性主要为紫红色泥质粉砂岩和泥岩。基底总体平坦，垂向上的地层结构层次分明（图2）。研究区主要为淤泥质粉质黏土和粉质黏土层，工程地质力学性质较差，且具有较高的压缩性。其中江北平均厚度约为18 m，两边厚中间薄，最厚可达35 m。

1.2 地下水

工作区的地下水可划分为松散岩类孔隙水和基岩裂隙水，其中，松散岩类孔隙水为区内主要地下水类型，根据埋藏条件、水理性质、水力特征，其又可划分为孔隙潜水和微承压水。水文地质剖面图如图3所示。

1）孔隙潜水

孔隙潜水在长江漫滩平原区广泛分布，含水层由全新统粉质黏土和淤泥质粉质黏土夹粉砂层薄层组成，厚度6～20 m。因含水层均由细颗粒土层组成，透水性和富水性差，单井出水量一般小于10 m3/d。地下水水位埋深1.0～1.5 m。区内地势平坦，地下水径流比较滞缓，水力坡度仅在

n×10-4～n×10-3之间。蒸发及向长江等地表水体排泄是其主要排泄方式。

2）微承压水

微承压水在长江漫滩平原区广泛分布。

其沉积物多呈二元或多元结构，下粗上细，含水层岩性包括粉细砂、中细砂和含砾中粗砂，砂层厚度一般为20～45 m。该含水层透水性和富水性良好，其单井涌水量为1 000～3 000 m3/d；在古河床沉积分布区单井涌水量在3 000 m3/d以上；而在接近阶地的漫滩边缘及局部含水层较薄地段，单井涌水量为100～1 000 m3/d。微承压水含水层的地下水流场基本保持天然状态，上部孔隙潜水下渗和内陆侧向渗流是其主要补给源，向长江水体侧向渗流及径流是其排泄的主要方式，而基坑降水是人工排泄的主要形式。

3）基岩裂隙水

基岩裂隙水含水层主要由白垩系上统浦口组（K2p）碎屑岩类组成，浅部以风化裂隙水为主，深部风化作用逐渐减弱，以构造裂隙水为主。岩层构造裂隙的发育程度总体较差，多为闭合状或被充填，富水性较差，一般单井涌水量小于100 m3/d；但在构造有利部位，特别是在断裂交会部位，构造裂隙发育程度较高，连通性较好，可能形成相对富水块段，单井涌水量可超过100 m3/d。由于裂隙含水层组分布位置一般较高，汇水条件差，出露面积有限，在区域上一般被视为弱含水层。工作区内该含水层均隐伏于第四系之下，主要接受孔隙潜水的垂向补给，以侧向径流的方式排泄。

2 概念模型的构建

本次研究取模型长、宽分别为670、80 m，底板高程为-42 m，模型表层为研究区地表。模型共划分为6层，隧道所在的第三层厚度约为12 m（图4a），其顶板主要穿越淤泥质粉质黏土， 底板主要穿越中粗砂。隧道的直径为6.28 m，盾构隧道西线和东线相距12.5 m（图4b）。西线盾构隧道开挖于东线开挖25 d后，此时计算的沉降量为累积地面沉降量，即西线、东线的地面沉降量之和。将模型四周边界概化为定水头边界，底板作为隔水边界，地表为自由边界，模型的初始应力取土体的自重，地下水的埋深取1.2 m。在x、y方向对模型前后两边界和左右两边界分别施加水平约束。

根据场地地层岩性划分为3个小区间段进行模拟计算，从而更准确地模拟由地铁盾构施工引发的地面沉降，分别为第1区间段（x=0～227 m）、第2区间段（x=227～450 m）、第3区间段（x=450～670 m），假定各个区间段地层岩性相近。盾构推进作为非连续的过程进行研究，并将自然充填在盾尾空隙的土体连同注浆后浆体概化为均质、等厚弹性的等代层[11]。

3 模型方法

3.1 Biot固结理论

假定土体为饱和土体，其土骨架变形为线弹性、微小变形，土体中地下水渗流符合达西定律，地下水不可压缩，其三维Biot固结方程[12]如下：

-G2wx-G1-2μxwxx+

wyy+

wzz+

ux=0，

-G2wy-G1-2μywxx+

wyy+

wzz+

uy=0，

-G2wz-G1-2μzwxx+

wyy+

wzz+

uz=-γ；

-twxx+

wyy+

wzz+1γw

［

xKxx

ux+

yKyy

uy+

zKzz

uz+γw］+W=0

。

（1）

式中：G为剪切模量；为拉普拉斯算子；μ为泊松比；wx、wy、wz分别为x、y、z方向上的位移；u为超静孔隙水压力；γ为土的重度；t为时间；Kxx、Kyy、Kzz分别为三个方向的渗透系数；γw为水的重度；W为源汇项。

3.2 土体本构模型

考虑土体的流变特性，将土体本构关系拓展到黏弹塑性，其总应变增量

dεij

分为弹性应变增量dεije、黏弹性应变增量dεijve和黏塑性应变增量dεijvp，即

dεij=

dεije+dεijve+dεijvp

。 （2）

dεije、dεijve、dεijvp分别计算如下。

1）弹性应变增量

由胡克定律可得

dεije=Cdσ′。（3）

式中：C为弹性柔度矩阵；

σ′为有效应力，

dσ′=[dσx′dσy′dτxz]T，

其中σx′=2Gμ1-2μεV+εx

，

σy′=2Gμ1-2μεV+εy

，εV

为土体固结形变，

εx=-wxx，

εy=-wyy，

τxz

为正应力在x轴方向时（即

σx

时）z轴方向的切应力，τxz=Gγxz，

γxz=-wxz+

wzx，

应力应变符号以压为正，以拉为负。

2）黏弹性应变增量

由kelvin流变模型

E1ε+Meε=σ

可得应力不变时复杂应力状态下的黏弹性应变增量为

dεijve= A{dσ′} E1（1-e-ηedt）。

（4）

式中：ε为kelvin体的应变；

σ

为kelvin体的应力；A为应力矩阵;

ηe=E1/Me，

E1为kelvin体黏弹性模量，Me为黏滞系数，根据以往工程经验，E1取7.8 MPa，Me取1×108 kPa/d。

3）黏塑性应变增量

dεijvp=2〈F〉tdtη2

2sin φ（3-sin φ）PI1+

3R

6J2

dσ。

（5）

式中：F为屈服函数；〈·〉为开关函数；

η2

为蠕变阶段的变形模量；φ为内摩擦角；I1为第1应力不变量；J2为应力偏量的第2不变量；P、R分别为系数矩阵。

3.3 定解条件

3.3.1 初始条件

1）初始孔隙水压力

u（x，y，z，t）｜t=0=u0（x，y，z）=γwh0（x，y，z，t0），（x，y，z）∈Ω。

（6）

式中：u0（x，y，z）为研究区域内已知初始孔隙水压力；

h0（x，y，z，t0）

为初始水头值； Ω为研究区；t0" 为初始时刻。

2）土体初始应力

σx=P0γh;σy=σx;σz=γh。

（7）

式中：σx、σy为土体初始水平向应力；σz为土体初始

垂向应力；h为每层深度；P0为静止侧压力系数，当土体为砂土时，

P0=1-sin φ′，

当土体为黏土时，

P0=0.95-sin φ′，

其中φ′为有效内摩擦角。

3）初始位移条件

w（x，y，z，t）｜t=0=0。

（8）

式中，

w（x，y，z，t0）

为初始位移。

3.3.2 边界条件

1）孔隙水压力边界条件Γ1：

u（x，y，z）Γ1=us。

（9）

式中，us为水头边界已知孔隙水压力。

2）自由面边界条件Γ2：

u=γwZ;q=-Sutcos θ。

（10）

式中：Z为自由面所在高程；S为给水度；θ为自由面外法线方向与垂线的交角；q为边界Γ２上的已知单位面积流量。

3）位移边界条件Γ3：

wx｜Γ3=w-x;

wy｜Γ3=w-y;

wz｜Γ3=w-z。

（11）

式中，w-x、

w-y、

w-z

分别为位移边界Γ3上三个方向的已知位移。

3.4 参数的动态变化

本文在建立数学模型的过程中考虑了盾构施工过程中土体参数，如土体固结形变εV、有效孔隙度n和渗透率k等参数的动态变化。基于固体骨架的连续性方程和Kozeny-Carman 方程最终得到关系方程[13 14]如式（12）所示。

n=1-（1-n0）（1-Δp/Es）1+εV

;

εV=-wxx+wyy+wzz

;

k=k01+εV·1+εVn03。

（12）

式中：Es为多孔介质骨架固体颗粒的体积弹性压缩模量；Δp为孔隙水压力；n0、k0分别为初始有效孔隙度和初始渗透率。

另外采用Duncan-Chang非线性模型[15]，得到切线弹性模量Et和切线泊松比μt的公式如式（13）所示。

Et=E0paσ3pam

1-Rf

（1-sin φ）（σ1-σ3）

2ccos φ+2σ3sin φ

2，

Rf=（σ1-σ3）f

（σ1-σ3）ult

，

（σ1-σ3）f=

2ccos φ+2σ3sin φ

1-sin φ;

μt=μ0-Flg（σ3/pa）

（1-Dε1）2

，

ε1=

σ1-σ3

kpa

σ3pa

n

1-Rf

（1-sin φ）（σ1-σ3）

2ccos φ+2σ3sin φ

2。

（13）

式中：pa为大气压强；E0为

σ3

等于100 kPa时初始弹性模量；m为弹性模量与固结压力曲线的斜率，一般在0.2～1.0之间；c为土体的黏聚力；Rf为破坏比；σ1、σ3

分别为最大主应力和最小主应力；

（σ1-σ3）f为摩尔库伦准则中的抗剪强度；

（σ1-σ3）ult

为主应力差极限值，是应变趋于无限大的主应力差；μ0为围压等于100 kPa时初始泊松比；ε1

为轴应变；

F、D为试验常数。

4 模型识别验证

根据上述数学模型，运用伽辽金有限元方法求解，采用Fortran95语言编制了地铁盾构引发地面沉降的三维全耦合数值分析程序模型，见图5。

模型在校准过程中分别完成了地下水位的校准和地面沉降的校准。其中，地下水位的校准选取的是位于计算区内且距离隧道最近的观测点，地面沉降的校准利用位于计算区内且距离隧道最近的监测点的数据。校准的结果以第1区间段为例（表1）。

对比地面沉降的实测值和计算值（图6）可知，本次模型所计算的沉降量与实测值趋势一致，该模型的地面沉降量的计算值和实测值在距离隧道-10～10 m之间的误差满足规范[16]要求，误差在合理范围内。

研究区地下水位拟合结果见图7。

选取计算区的地下水位观测点，并使用2016年12月—2021年12月的观测资料作为其实测值，其中每个月作为一个应力期，每个应力期为1个时间步长。对比地下水位的实测及计算值（图7），可见拟合度较好。综上所述，本次研究所构建的模型可靠。

5 地面沉降预测及影响

5.1 地面沉降预测

本次研究区区间全长670 m，将模型共划分为

224个开挖步，每个开挖步长度3 m，设定推进速率为9 m/d。随着每推进一个开挖步，将开挖单元模量参数设为最小值。图8为模型计算的东线盾构开挖至第2区间第19步（即东线开挖至284 m，西线开挖至57 m）时地面沉降平面等值线。东、西线轴线的位置分别为y=-10、10 m。

由图8可知：东线盾构隧道开挖的x方向地面沉降影响范围在y=-25～5 m之间，基本轴线两边对称分布；东、西线在y方向的沉降影响范围在y=-29～29 m之间，且其轴线之间的沉降量较大。

图9为东线盾构隧道开挖至第2区间第19步时，分别选取x=40、60、80 m所统计的横向断面地面沉降量计算值。当x=40 m时，东、西线盾构开挖的最大沉降量分别位于y=-8.9 m和y=9.3 m

处，相应的沉降量分别为54.20 mm和53.47 mm；x=60 m时，东、西线最大沉降量分别位于y=-10.3 m和y=9.5 m处，相应的沉降量分别为53.27 mm和43.05 mm；x=80 m时，东线于y=-10.9 m处沉降最多，沉降了51.65 mm，西线于9.0 m处沉降了9.75 mm，""" 该处为沉降最大处。此外，东线轴线处沉降量在52.41～54.52 mm范围内，西线轴线处沉降量在53.28～55.60 mm范围内。

东线盾构隧道先行开挖引起的地面沉降基本以隧道轴线为中心对称分布，其施工引起的地面沉降影响半径约为15.2 m。西线隧道开始开挖后，沉降中心向东偏移，这可能是由于西线开挖导致的土体应力重分布，东侧土体强度降低，故而沉降中心向东偏移。但纵观整体，其沉降规律保持不变。

5.2 参数变化

土体的压缩势必会造成土体参数的变化。盾构施工引发土体的变形，会造成相应的土体力学及水力学参数的变化，本文将探讨地铁盾构施工作用下这些参数动态的变化关系。

整个模型划分3个区间段进行开挖，以第2区间段上模型第2层一个拱顶上方25811单元和模型第4层一个拱底下方20358单元为例，分别对25811单元和20358单元的孔隙度、泊松比、垂向和水平渗透系数及变形模量的变化进行分析，结果见图10。图10表明：拱顶上方单元的孔隙度、泊松比、垂向渗透系数和水平渗透系数均随着开挖时间的增加而减小，变形模量则随着开挖时间的增加而增大；拱顶下方单元参数的变化和上方的相反，即孔隙度、泊松比、垂向渗透系数和水平渗透系数均随着开挖时间的增加而增大，变形模量随着开挖时间的增加而减小。

6 结论

1）基于Biot固结理论，考虑土体参数的动态变化，构建了考虑地下水渗流场影响的南京市11号线地铁盾构开挖所引起地面沉降的三维全耦合模型。利用数值模拟对软土地区地铁盾构引起的土体变形进行研究，并描绘了土体的变形。

2）将场地的地层岩性划分为3个小区间段，以第2区间段（x=227～450 m）为例，模型计算第2区间东线隧道轴线处地面沉降在52.41～54.52 mm范围内，西线隧道轴线处地面沉降在53.28～55.60 mm范围内。可见研究区的沉降较大，开挖时应注意支护，防止产生大变形。

3）随着开挖的进行，渗透系数随着有效孔隙度的变化而发生变化。盾构拱顶上方土体固结压缩使其出水能力减弱，孔隙度、渗透系数及泊松比变小，变形模量则缓慢增大，这些变化趋势在开挖的后期逐渐变缓乃至不变；盾构底部下方土体变化刚好相反。

参考文献（References）：

[1] Yan Y G， Yang Q， Jia Z H， et al. Integration of Multiband InSAR and Leveling Measurements for Analyzing the Surface Subsidence of Shield Tunneling at Beijing-Zhangzhou High-Speed Railway[J]. Journal of Sensors， 2021， 2021（1）：124.

[2] Pan H Z， Jiang B Y， Huang T. An Adaptability Study of Gaussian Equation Applied to Predicting Ground Settlements Induced by Shield Tunneling in Tianjin[J]. Science of Surveying and Mapping， 2010， 35（3）： 5355.

[3] Yu G M， Xia P J， Zhang Y Y， et al. Land Lapse Phenomenon and Surface Subsidence Law of Shield Tunnel Passing Through Rock-Bearing Formation[J/OL]. Frontiers in Earth Science，2022，10：113[20230604]. doi：10.3389/feart.2022.934360.

[4] Adachi T， Jun L. Finite Element Analysis of Biot’s Consolidation in Tunnel Excavation Based on a Constitutive Model with Strain Softening[C]//Yuan J X. Proceedings of the 9th International Conference on Computer Methods and Advances in Geomechanics. Wuhan：9th International Conference on Computer Methods and Advances in Geomechanics， 1997： 12691272.

[5] Verruijt A. Non-Linear Effects in Poroelasticity[C]//Thimus J F. Proceedings of the Biot Conference on Poromechanics." Louvain La Neuve：[s.n.]， 1998： 525530.

[6] Gutierrez M S， Lewis R W. Coupling of Fluid Flow and Deformation in Underground Formations[J]. Journal of Engineering Mechanics， 2002， 128（7）： 779787.

[7] Burbey T J. Three Dimensional Aquifer Deformation and Flow Modeling[C]//Hassanizadeh S M. Proceedings of the 14th International Conference on Computational Methods in Water Resources. Delft：[s.n.]， 2002： 129136.

[8] Sabetamal H， Nazem M， Sloan S W， et al. Numerical Simulation of Dynamic Pore Fluid-Solid Interaction in Fully Saturated Non-Linear Porous Media[C]//Onate E. Proceedings of the 11th International Conference on Computational Plasticity （COMPLAS）-Fundamentals and Applications. Barcelona：[s.n.]， 2011： 12631273.

[9] 孟世豪，崔亚莉，田芳，等.基于MODFLOWSUB建立变渗透系数的地下水流地面沉降模型[J].吉林大学学报（地球科学版）， 2022， 52（2）：550559.

Meng Shihao，Cui Yali，Tian Fang，et al. Modeling of Groundwater Flow-Land Subsidence with Variable Hydraulic Conductivity Based on MODFLOW-SUB[J]. Journal of Jilin University （Earth Science Edition）， 2022， 52（2）：550559.

[10] 袁帅，王君，吴朝峰，等.虹吸排水法处理软土地基的水位与沉降计算模型[J].吉林大学学报（地球科学版）， 2024， 54（1）：208218.

Yuan Shuai，Wang Jun，Wu Zhaofeng，et al. Calculation Model for Water Level and Settlement of Soft Foundation Treated by Siphon Drainage[J]. Journal of Jilin University （Earth Science Edition） ， 2024， 54（1）：208218.

[11] 贾善坡，陈卫忠，于洪丹，等.渗流应力耦合作用下深埋黏土岩隧道盾构施工特性及其动态行为研究[J].岩石力学与工程学报， 2012，" 31（增刊1）：26812691.

Jia Shanpo， Chen Weizhong， Yu Hongdan， et al. Study of Construction Characteristics and Dynamicbehavior of Deep Clay Stone During Shield Tunneling Under Seepage-Stress Coupling Effect[J]. Chinese Journal of Rock Mechanics and Engineering， 2012，" 31（Sup.1）：26812691.

[12] 骆祖江，刘金宝，李朗.第四纪松散沉积层地下水疏降与地面沉降三维全耦合数值模拟[J].岩土工程学报， 2008， 30（2）：193198.

Luo Zujiang， Liu Jinbao， Li Lang. Three-Dimensional Full Coupling Numerical Simulation of Groundwater Dewatering and Land-Subsidence in Quaternary Loose Sediments[J]. Chinese Journal of Geotechnical Engineering， 2008， 30（2）：193198.

[13] 冉启全，李士伦.流固耦合油藏数值模拟中物性参数动态模型研究[J].石油勘探与开发， 1997， 24（3）： 6165，100.

Ran Qiquan， Li Shilun. Study on Dynamic Models of Reservoir Parameters in the Coupled Simulation of Multiphase Flow and Reservoir Deformation[J]. Petroleum Exploration and Development， 1997， 24（3）： 6165，100.

[14] 田杰，刘先贵，尚根华.基于流固耦合理论的套损力学机理分析[J].水动力学研究与进展：A辑， 2005， 20（2）： 221225.

Tian Jie， Liu Xiangui， Shang Genhua. Casing Damage Mechanism Based on Theory of Fluid-Solid Coupling Flow Through Underground Rock[J]. Journal of Hydrodynamics：A， 2005， 20（2）： 221225.

[15] 罗刚，张建民.邓肯张模型和沈珠江双屈服面模型的改进[J].岩土力学， 2004， 25（6）： 887890.

Luo Gang， Zhang Jianmin. Improvement of Duncan-Chang Nonlinear Model and Shen Zhujiang’s Elastoplastic Model for Granular Soils[J]. Rock and Soil Mechanics， 2004， 25（6）： 887890.

[16] 地下水流数值模拟技术要求： GW1D1[Z].北京：中国地质调查局， 2004.

Technical Requirements for Numerical Simulation of Groundwater Flow： GW1D1[Z]. Beijing：China Geological Survey， 2004.