南京江北新區地鐵盾構開挖引發地面沉降三維數值模擬預測

摘要：為研究地鐵盾構施工引起的地面沉降問題，首先以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.