熊小鋒，施小清，吳劍鋒，吳吉春

(南京大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 210023)

彈塑性變形條件下抽水引起的地面沉降三維數(shù)值模擬

熊小鋒，施小清，吳劍鋒，吳吉春

(南京大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 210023)

基于部分耦合原理，采用TOUGH2和FLAC3D建立抽水引起的三維地面沉降彈塑性模型，模型中綜合考慮土體的彈塑性變形特征、滲流-應(yīng)力的雙向耦合作用以及參數(shù)的非線(xiàn)性，探討了持續(xù)抽水和脈沖抽水兩種抽水過(guò)程中地面沉降發(fā)展演化過(guò)程。研究結(jié)果表明：(1)集中抽水停止后地面沉降會(huì)發(fā)生回彈，抽水中心沉降量不斷減小。由于水平方向存在水力梯度，地下水繼續(xù)向地下水位漏斗中心滲流從而導(dǎo)致沉降漏斗的范圍仍繼續(xù)擴(kuò)大；(2)脈沖抽水導(dǎo)致土體的孔隙水壓力、滲透系數(shù)以及沉降量均呈周期性波動(dòng)變化，地面沉降會(huì)局部回彈，但總體仍隨著抽水的持續(xù)，沉降量不斷增加；(3)在抽水量相同前提下，對(duì)比持續(xù)抽水與脈沖抽水兩種方式引發(fā)的塑形沉降量可知，抽水速率小、脈沖式多次開(kāi)采導(dǎo)致的塑性沉降量較小，持續(xù)抽水的抽水速率越小、脈沖抽水間隔越短越有利于控制地面沉降。研究成果為地面沉降數(shù)值模擬提供了一種新方法，其中算例研究能為抽水條件下地面沉降的控制提供參考。

TOUGH2-FLAC3D；地面沉降；彈塑性變形；脈沖抽水

水-土耦合作用數(shù)值模擬是目前抽水引發(fā)地面沉降過(guò)程研究中的熱點(diǎn)[1～5]。地面沉降過(guò)程中滲流-力學(xué)場(chǎng)互相影響，滲流改變力學(xué)場(chǎng)，引起土體發(fā)生變形，土體發(fā)生變形后滲透系數(shù)和孔隙度均會(huì)減小，反過(guò)來(lái)影響滲流過(guò)程，因此地面沉降的滲流-力學(xué)作用是一個(gè)雙向耦合的過(guò)程。近幾十年來(lái)，已經(jīng)開(kāi)發(fā)了眾多的地質(zhì)力學(xué)與多相流場(chǎng)的耦合模擬器，其中Kim系統(tǒng)地驗(yàn)證了此類(lèi)耦合模擬器的優(yōu)勢(shì)[6]。TOUGH2-FLAC3D是比較成熟的THM(溫度-滲流-應(yīng)力)耦合模擬工具之一[7～10]，最先由美國(guó)勞倫斯伯克利實(shí)驗(yàn)室的Rutqvist于2002年開(kāi)發(fā)，并成功用于深部咸水CO2封存條件下的流體運(yùn)移、蓋層力學(xué)變化以及地表位移研究[7～8]。其中TOUGH2是一款國(guó)際通用非飽和多相流模擬軟件，采用積分有限差分法，通過(guò)離散求解質(zhì)量和能量平衡線(xiàn)性方程組，能夠?qū)崿F(xiàn)溫度-滲流(TH)的全耦合計(jì)算[11～12]；FLAC3D是Itasca公司開(kāi)發(fā)的一款基于快速拉格朗日方法的連續(xù)介質(zhì)力學(xué)分析的軟件。將上述兩款成熟商業(yè)軟件聯(lián)合應(yīng)用，可以發(fā)揮各自?xún)?yōu)勢(shì)。

TOUGH2-FLAC3D已在多場(chǎng)耦合計(jì)算領(lǐng)域得到廣泛的應(yīng)用，但是目前國(guó)內(nèi)外很少將TOUGH2-FLAC3D應(yīng)用于地面沉降數(shù)值模擬研究。在地面沉降數(shù)值模擬研究方面，Biot流固全耦合理論上完美，但全耦合計(jì)算效率不高，特別是當(dāng)模型復(fù)雜、大區(qū)域尺度時(shí)應(yīng)用受限[13～14]。TOUGH2-FLAC3D耦合采用TOUGH2計(jì)算滲流，將滲流場(chǎng)計(jì)算結(jié)果傳遞給FLAC3D，這樣一來(lái)FLAC3D只需要計(jì)算力學(xué)過(guò)程，從而提高計(jì)算效率，同時(shí)適當(dāng)增加耦合步數(shù)即可提高計(jì)算精度。此外，F(xiàn)LAC3D內(nèi)置多種經(jīng)典應(yīng)力應(yīng)變本構(gòu)模型，適用于各類(lèi)土體的地面沉降研究，使得TOUGH2-FLAC3D在模擬復(fù)雜土體地面沉降方面具有較大潛力。

本文基于部分耦合原理，利用多場(chǎng)耦合模擬器TOUGH2-FLAC3D展開(kāi)抽水引起的地面沉降三維數(shù)值模擬研究工作，模型中綜合考慮了土體的彈塑性特征及參數(shù)的非線(xiàn)性變化，通過(guò)建立持續(xù)抽水及脈沖式抽水算例，探討有水位恢復(fù)時(shí)地面沉降的沉降-回彈特征，并根據(jù)塑性沉降量對(duì)地面沉降做進(jìn)一步分析。

1 基礎(chǔ)理論及方法介紹

1.1 TOUGH2滲流場(chǎng)模型

采用部分耦合方法，滲流模型和力學(xué)模型分別由TOUGH2和FLAC3D建立。TOUGH2采用積分有限差分方法進(jìn)行數(shù)值計(jì)算，其程序結(jié)構(gòu)采用模塊化設(shè)計(jì)，各子模塊采用不同的狀態(tài)方程(EOS)。模擬地面沉降時(shí)忽略溫度效應(yīng)，采用等溫模式，子模塊選用EOS9。TOUGH2中通用的質(zhì)量(或能量)平衡方程[11]為：

(1)

式中：t——時(shí)間/s；κ——不同組分；Vn——第n個(gè)單元的體積/m3；Γn——第n個(gè)單元的邊界； n——dΓn的法向量。

在只考慮質(zhì)量平衡的數(shù)學(xué)模型中，Mκ表示單位體積內(nèi)不同組分的質(zhì)量，F(xiàn)κ表示質(zhì)量通量，qκ表示單位體積的源匯項(xiàng)?；谶_(dá)西定律，Mκ和Fκ一般形式為：

(2)

(3)

在TOUGH2中，對(duì)平衡方程采用空間積分有限差分法(IFDM)和時(shí)間上的一階有限差分法[11～12]進(jìn)行離散。

1.2FLAC3D力學(xué)場(chǎng)模型

在TOUGH2中設(shè)置抽水井，抽水會(huì)引起含水層中孔隙水壓力減小，從而引起力學(xué)模型中土體有效應(yīng)力增加，土體發(fā)生垂向壓縮變形；遠(yuǎn)井處水流沿著近水平方向滲流到抽水井，在水平方向存在水力梯度，土體在滲透力作用下也會(huì)發(fā)生水平變形[15～17]。FLAC3D在力學(xué)計(jì)算方面的基本控制方程包括平衡方程、傳導(dǎo)方程和本構(gòu)方程[18～19]，由運(yùn)動(dòng)平衡方程得到節(jié)點(diǎn)的位移和速率，再根據(jù)空間導(dǎo)數(shù)得到單元的應(yīng)變率，采用本構(gòu)關(guān)系得到新的單元應(yīng)力，產(chǎn)生新的運(yùn)動(dòng)平衡，以此完成一個(gè)顯示拉格朗日計(jì)算循環(huán)[17]。本文地面沉降模型在FLAC3D中打開(kāi)流固耦合(HM)計(jì)算模式，但是在一個(gè)耦合迭代步內(nèi)將孔壓固定且不產(chǎn)生滲流，接收由TOUGH2傳遞來(lái)的孔壓作為力學(xué)擾動(dòng)。當(dāng)忽略流體流動(dòng)及溫度效應(yīng)時(shí)，流固耦合數(shù)學(xué)模型如下所示[19]：

(4)

式中：M——比奧模量/Pa； s——飽和度； φ——孔隙度； p——孔隙水壓力(不包含大氣壓)/Pa； α——比奧系數(shù)； ε——體應(yīng)變。

沉降變形采用摩爾-庫(kù)倫(Mohr-Coulomb)彈塑性本構(gòu)模型，該模型綜合了胡克定律和Coulomb破壞準(zhǔn)則，認(rèn)為土體在屈服和破壞前應(yīng)力應(yīng)變符合線(xiàn)彈性關(guān)系。FLAC3D中基于胡克線(xiàn)彈性定律的本構(gòu)函數(shù)H形式如下[19]：

(5)

式中：上標(biāo)o——初始狀態(tài)；εij——應(yīng)變；K——體積模量/Pa；G——剪切模量/Pa。

應(yīng)力達(dá)到屈服準(zhǔn)則的判別標(biāo)準(zhǔn)時(shí)，土體發(fā)生塑性變形。彈塑性本構(gòu)模型中將變形分為兩部分：一部分為可恢復(fù)的彈性變形，另外一部分為不可恢復(fù)的塑性變形。發(fā)生塑性變形時(shí)，根據(jù)流動(dòng)法則需要對(duì)塑性情況下的應(yīng)力應(yīng)變?cè)隽窟M(jìn)行修正。應(yīng)力繼續(xù)增加，達(dá)到破壞準(zhǔn)則時(shí)土體發(fā)生剪切或拉張破壞。

基于主應(yīng)力的摩爾庫(kù)倫破壞準(zhǔn)則表達(dá)式為：

(6)

式中：σ1、σ3——最大、最小有效主應(yīng)力/Pa； c——黏聚力/Pa； φ——摩擦角/(°)。 當(dāng)fs=0或ft=0時(shí)，土體將發(fā)生破壞。

1.3 參數(shù)的非線(xiàn)性關(guān)系

滲流-應(yīng)力相互作用導(dǎo)致土體應(yīng)力應(yīng)變狀態(tài)發(fā)生改變，從而使得土體的參數(shù)發(fā)生變化。目前已有的研究中針對(duì)參數(shù)的非線(xiàn)性關(guān)系并沒(méi)有統(tǒng)一標(biāo)準(zhǔn)，大多采取經(jīng)驗(yàn)公式方法[2,4]。本文考慮土體單元孔隙度和滲透率隨有效應(yīng)力變化，其表達(dá)式如下[7,20]：

(7)

(8)

平均有效應(yīng)力在FLAC3D中不是內(nèi)置變量，可以通過(guò)fish代碼為每個(gè)單元增加一個(gè)額外變量實(shí)現(xiàn)同步計(jì)算。使用下式計(jì)算平均有效應(yīng)力：

(9)

P——孔隙水壓力/Pa。

1.4TOUGH2-FLAC3D參數(shù)傳遞

為順次執(zhí)行TOUGH2和FLAC3D，需要在原基礎(chǔ)上對(duì)TOUGH2代碼修改，并經(jīng)由FLAC3D的fish語(yǔ)言控制迭代及參數(shù)傳遞過(guò)程。修改TOUGH2代碼后，每個(gè)耦合時(shí)步中TOUGH2計(jì)算達(dá)到設(shè)定時(shí)間即停止，將計(jì)算得到的孔隙水壓力、飽和度、溫度等變量傳遞給FLAC3D，由FLAC3D計(jì)算力學(xué)變形。其中TOUGH2傳遞變量給FLAC3D時(shí)需要插值，因?yàn)榭紫端畨毫Α柡投?、溫度等變量賦值在TOUGH2網(wǎng)格中心，而在FLAC3D模型中則賦值在網(wǎng)格節(jié)點(diǎn)；FLAC3D傳遞有效應(yīng)力、應(yīng)變給TOUGH2時(shí)則不需要插值，這兩個(gè)變量本身就賦值在FLAC3D單元中心。TOUGH2-FLAC3D變量傳遞比較耗時(shí)，但是在高版本的fish語(yǔ)言中新增加了節(jié)點(diǎn)與單元鏈接函數(shù)，不再需要重復(fù)執(zhí)行節(jié)點(diǎn)和單元的遍歷循環(huán)，使得TOUGH2-FLAC3D的變量傳遞效率大大提高[9]。

TOUGH2到FLAC3D的插值過(guò)程中，采用距離反比計(jì)算權(quán)重，將單元變量轉(zhuǎn)化到節(jié)點(diǎn)上。采用規(guī)則六面體網(wǎng)格剖分時(shí)，1≤n≤8。網(wǎng)格節(jié)點(diǎn)與單元連接距離示意如圖1所示。

圖1 TOUGH2-FLAC3D節(jié)點(diǎn)到單元中心距離示意圖Fig.1 The distance of grid points to zone node for TOUGH2-FLAC3D

2 三維地面沉降耦合模型建立

2.1 網(wǎng)格剖分

建立抽水條件下的三維地面沉降網(wǎng)格模型。模型徑向尺寸為1×105m，XY面上剖分900個(gè)網(wǎng)格，抽水井處網(wǎng)格尺寸加密到20 m，遠(yuǎn)離抽水井網(wǎng)格尺寸適當(dāng)增加，靠近邊界處網(wǎng)格尺寸較小以方便設(shè)置邊界條件；Z方向?yàn)?00 m，靠近模型頂部和底部厚度為5 m，其他每層厚度15 m，總共剖分8層。坐標(biāo)原點(diǎn)位于含水層頂部正中間。抽水井為點(diǎn)源式，位于模型中間底部。模型網(wǎng)格剖分如圖2所示。

圖2 地面沉降模型三維網(wǎng)格剖分(Z軸放大200倍)Fig.2 The three-dimensional grid for land subsidence model (Z axis amplification equals 200)

2.2 模型概化

為方便研究，將模型概化為上下兩層含水層，上層為弱透水層，下層為含水層。弱透水層土體材料為砂質(zhì)黏土，范圍為自頂部0 ～-50 m；下層含水層為細(xì)砂，范圍自-50～-100 m。力學(xué)模型中，左右兩側(cè)為輥軸邊界，底部為固定邊界，頂部為自由邊界。滲流模型中，四周為給定水頭邊界，孔壓保持不變，上下邊界均為不透水邊界，含水層為承壓含水層，模型上頂板初始孔隙水壓力為3.22×105Pa。

圖3 模型概化及邊界設(shè)計(jì)Fig.3 Conceptual model and boundary design

摩爾-庫(kù)倫公式中土體的強(qiáng)度主要取決于黏聚力和摩擦角。土體破壞形式主要為剪切破壞和拉張破壞，土體一般不承受拉力作用[21]，在FLAC3D模型中人為將抗拉強(qiáng)度設(shè)定為一個(gè)較大的值，使得本文模型不產(chǎn)生拉張破壞。模型中流體密度取1×103kg/m3，土體密度2.0×103kg/m3，重力加速度取9.81 m/s2。殘留孔隙度根據(jù)相關(guān)文獻(xiàn)[7]取初始孔隙度的90%，其他與土體有關(guān)的參數(shù)取值如表1。

表1 模型參數(shù)取值

2.3 模型分析

模型為單井抽水三維地面沉降彈塑性模型，模型XY方向尺寸足夠大，以減少人為設(shè)定邊界的影響。模型的四周邊界設(shè)定為給定水頭邊界，停止抽水后水位能夠自由恢復(fù)。以研究點(diǎn)距離抽水井的水平距離r作為井距。當(dāng)以定流量400 m3/d連續(xù)抽水150 d時(shí)抽水井上方井距r=10 m和r=29 633 m處孔隙水壓力變化如圖4。從圖中可以看到，一旦開(kāi)始抽水，抽水井附近的孔壓會(huì)“瞬時(shí)”下降，隨后緩慢下降；而遠(yuǎn)離抽水井位置，抽水一段時(shí)間后才影響到，孔壓下降為增速趨勢(shì)。單井抽水在模型中產(chǎn)生了非穩(wěn)定流場(chǎng)，其孔壓(或水頭)下降規(guī)律符合Theis理論?？紤]到沉降與水位的對(duì)應(yīng)關(guān)系，抽水引起的地面沉降變化也是非穩(wěn)定的，沉降漏斗的中心沉降量會(huì)迅速趨于穩(wěn)定，抽水后期主要增加沉降漏斗的擴(kuò)展范圍。

圖4 孔隙水壓力隨時(shí)間變化(Z=0)Fig.4 The pore pressure evolution with time(Z=0)

3 彈塑性地面沉降算例結(jié)果分析

3.1 持續(xù)抽水沉降規(guī)律

在保證抽水量相同的前提下，根據(jù)持續(xù)抽水時(shí)間長(zhǎng)短及抽水速率的不同建立如表2所示三種情景算例?？偰M時(shí)間均為150 d，其中持續(xù)抽水60 d和100 d之后停止的情況下，水位存在恢復(fù)過(guò)程，地面沉降也會(huì)跟著回彈。抽水期間地面沉降量不斷增加，停止抽水后近井處地面沉降立即回彈。水位通過(guò)滲流自由恢復(fù)，水流在水力梯度下從遠(yuǎn)井處流向近井處，逐漸補(bǔ)充近井空缺，在模擬期內(nèi)即使彈性變形量也未完全恢復(fù)。

表2 持續(xù)抽水情景設(shè)計(jì)

(1)沉降規(guī)律分析

一旦開(kāi)始抽水，抽水井附近的孔壓會(huì)迅速降低，地面產(chǎn)生“瞬時(shí)”沉降，這種瞬時(shí)沉降現(xiàn)象與該處孔壓變化是對(duì)應(yīng)的。隨后沉降量隨抽水時(shí)間逐漸增加，沉降增速減緩。沉降最大的位置位于抽水井附近，隨距離增加地面沉降量呈現(xiàn)梯度減小。圖5～圖7顯示不同開(kāi)采情景地面沉降隨時(shí)間發(fā)展規(guī)律，其中情景一抽水速率較大，因此抽水60 d時(shí)最大沉降量達(dá)到了24.1 cm，情景二、三的歷史最大沉降量分別為17.5 cm和13.9 cm。

情景一、二均存在一個(gè)水位恢復(fù)過(guò)程，情景一抽水60 d后停止，后90 d水位自由恢復(fù)。停止抽水后遠(yuǎn)井和近井處存在一定的水力梯度，遠(yuǎn)井處的水會(huì)繼續(xù)流向近井處，沉降漏斗的范圍會(huì)繼續(xù)擴(kuò)展，而沉降漏斗中心沉降量則不斷減小，見(jiàn)圖8～圖9沉降云圖。

圖5 持續(xù)抽水60 d沉降隨時(shí)間變化曲線(xiàn)(情景一)Fig.5 Displacement evolution for pumping lasting 60 days (the first scenario)

圖6 持續(xù)抽水100 d沉降隨時(shí)間變化曲線(xiàn)(情景二)Fig.6 Displacement evolution for pumping lasting 100 days (the second scenario)

圖7 持續(xù)抽水150 d沉降隨時(shí)間變化曲線(xiàn)(情景三)Fig.7 Displacement evolution for pumping lasting 150 days (the third scenario)

圖8 持續(xù)抽水60d沉降云圖(Z放大200倍，最小顯示值人為設(shè)定-0.01 m)Fig.8 Displacement contour evolution for pumping last 60 days (Z axis amplification equals 200, the minimize value is -0.01 m)

圖9 持續(xù)抽水60 d沉降云圖切片(Z放大200倍，最小顯示值人為設(shè)定-0.01 m)Fig.9 Displacement contour cut for pumping last 60 days (Z axis amplification equals 200, the minimize value is-0.01 m)

(2)塑性沉降量分析

抽水過(guò)程中土體有效應(yīng)力會(huì)增加，當(dāng)有效應(yīng)力達(dá)到彈塑性模型的屈服準(zhǔn)則時(shí)產(chǎn)生塑性變形。停止抽水后水位會(huì)逐漸恢復(fù)，土體的彈性變形也隨即回彈，而已經(jīng)發(fā)生的塑性變形則無(wú)法回彈。圖10顯示與彈性模型計(jì)算結(jié)果相比，彈塑性模型的沉降量總保持一個(gè)差值，兩者的差值情況如圖11，曲線(xiàn)在抽水前期波動(dòng)較大，后期趨于穩(wěn)定。由于本文三維地面沉降模型中綜合考慮了參數(shù)的非線(xiàn)性關(guān)系，因此差值不僅僅代表塑性應(yīng)變量，還包括參數(shù)非線(xiàn)性的影響，但是在總體上可以認(rèn)為彈性與彈塑性模型沉降量差值能夠代表塑性變形狀態(tài)。三種開(kāi)采情景對(duì)比，持續(xù)抽水60 d、100 d和150 d的塑性變形量分別為4.0 cm、3.7 cm和3.4 cm。由此可見(jiàn)塑性沉降量主要與抽水速率有關(guān)，抽水速率越大，歷史沉降量最大值越大，最后導(dǎo)致的塑性沉降量也越大。從地面沉降控沉管理角度來(lái)分析，由于彈性變形量可恢復(fù)，而塑性變形量不可恢復(fù)，因此實(shí)際工程中人們更關(guān)注的是抽水引起的地面沉降塑性沉降量。采用長(zhǎng)時(shí)間慢速抽水時(shí)導(dǎo)致的塑性變形量更小，在保證日開(kāi)采量滿(mǎn)足需求的情況下，可適當(dāng)減小抽水速率以控制地面沉降的塑性沉降量。

圖10 r=10 m處彈性與彈塑性沉降量對(duì)比圖Fig.10 Compare displacement for elastic and elasto-plastic model (r=10 m)

圖11 持續(xù)抽水r=10 m彈性、彈塑性沉降量差值Fig.11 The displacement difference for elastic and elasto-plastic model (r=10 m)

3.2 脈沖抽水沉降規(guī)律

地下水開(kāi)采過(guò)程中，為了控制地面沉降，往往采用多次抽取地下水的開(kāi)采方式。設(shè)計(jì)如表3所示脈沖抽水情景方案，不同方案均從零時(shí)刻開(kāi)始抽水，抽水總時(shí)間均為80 d，各方案保持最終開(kāi)采量一致，均為6×104m3。

表3 脈沖抽水情景設(shè)計(jì)

(1)孔壓及參數(shù)變化分析

采用脈沖式抽水，近井處水位恢復(fù)較快，遠(yuǎn)井處水位恢復(fù)相對(duì)延遲。由于存在脈沖作用，脈沖間隔期遠(yuǎn)井處的水量能夠補(bǔ)充近井，使得近井處的孔壓下降有所減小。以脈沖抽水情景方案四為例，近井和遠(yuǎn)井處孔隙水壓力隨時(shí)間的變化規(guī)律如圖12。

圖12 脈沖抽水含水層頂板壓強(qiáng)變化(情景四)Fig.12 The pore pressure evolution with time at the top of the aquifer for pulse pumping (the fourth scenario)

根據(jù)理論及實(shí)踐經(jīng)驗(yàn)得出[4, 22～23]，研究地面沉降時(shí)考慮參數(shù)的非線(xiàn)性十分必要，實(shí)際工程中地下水位或地表高程并非持續(xù)下降，而是在時(shí)空上呈現(xiàn)波動(dòng)狀態(tài)。本文在TOUGH2-FLAC3D耦合地面沉降模型中綜合考慮了參數(shù)的非線(xiàn)性。抽水過(guò)程中土體被壓縮產(chǎn)生沉降，土體單元的孔隙度和滲透系數(shù)也會(huì)相應(yīng)減小，脈沖式抽水過(guò)程中土體的壓縮與回彈以波動(dòng)式進(jìn)行，參數(shù)的變化規(guī)律也呈相同變化規(guī)律。

以脈沖抽水情景方案四為例，參數(shù)的動(dòng)態(tài)變化規(guī)律如圖13。近井處的土體除了塑性應(yīng)變量之外能快速回彈，因此孔隙度和滲透率也基本能夠恢復(fù)到初始水平，而遠(yuǎn)井處由于抽水的延遲影響作用，參數(shù)在波動(dòng)變化的同時(shí)，總體上呈現(xiàn)減小的規(guī)律。

圖13 脈沖抽水含水層頂板參數(shù)動(dòng)態(tài)變化曲線(xiàn)(情景四)Fig.13 The parameters evolution with time at the top of the aquifer for pulse pumping (the fourth scenario)

(2)沉降規(guī)律分析

三種脈沖抽水方案的區(qū)別主要在于脈沖間隔不同，而抽水速率及總抽水時(shí)間相同。情景四、五、六的脈沖周期分別為8、4和2個(gè)。

脈沖抽水條件下彈塑性模型沉降變化曲線(xiàn)如圖14。與持續(xù)抽水沉降規(guī)律類(lèi)似，圖中顯示在抽水井附近(r=10 m)，一旦開(kāi)始抽水即會(huì)產(chǎn)生一個(gè)瞬時(shí)沉降量，而在遠(yuǎn)井處沉降從零開(kāi)始逐漸發(fā)展。隨后的每次脈沖時(shí)段內(nèi)，沉降會(huì)在某個(gè)范圍內(nèi)循環(huán)回彈，脈沖間隔較長(zhǎng)則水位有足夠的恢復(fù)時(shí)間，沉降回彈也越多。情景四中，沉降中心(r=10 m)最大沉降量在時(shí)間上總體上不斷增加，最大沉降量為15.8 cm，情景五和情景六最大歷史沉降量為17.1 cm和20.9 cm。

(3)塑性沉降量分析

以土體應(yīng)力應(yīng)變符合線(xiàn)彈性本構(gòu)關(guān)系建立彈性沉降模型，基于相同抽水條件計(jì)算彈性沉降量，與彈塑性模型的沉降量之間的差值如圖15。與持續(xù)抽水變化規(guī)律類(lèi)似，塑性沉降前期仍然處于發(fā)展期，后期穩(wěn)定。根據(jù)計(jì)算結(jié)果，以差值代表塑性沉降量，在模擬期內(nèi)情景四、五、六產(chǎn)生的塑性沉降量分別為3.7 cm、3.8 cm和4.1 cm。

從150 d以后脈沖抽水也完全停止，水位自由恢復(fù)，由于邊界為給定水頭邊界，最終土體彈性部分完全回彈而僅保留塑性沉降，可以塑性沉降量作為控沉指標(biāo)。與持續(xù)抽水情景方案對(duì)比，脈沖抽水速率均為750 m3/d，比情景一小，而比情景二、三大，塑性沉降量在很大程度上與抽水速率有關(guān)，因此持續(xù)抽水與脈沖抽水不能僅從結(jié)果上進(jìn)行對(duì)比。在持續(xù)抽水算例中，認(rèn)為抽水速率小引起的塑性沉降量也較??；在脈沖抽水算例中，認(rèn)為脈沖間隔短，引起的塑性沉降量小。綜合各類(lèi)影響因素，塑性沉降量取決于抽水過(guò)程中歷史最大沉降量，在工程降水或開(kāi)采地下水過(guò)程中結(jié)合慢速、多次開(kāi)采方式可以有效控制地面沉降的塑性沉降量。

圖15 脈沖抽水r=10 m彈性、彈塑性沉降量差值Fig.15 The displacement difference for elastic and elasto-plastic model (r=10 m)

4 結(jié)論

本文通過(guò)TOUGH2-FLAC3D建立抽水引起的三維地面沉降彈塑性模型，分析了集中抽水一段時(shí)間和脈沖抽水的地面沉降彈塑性變形規(guī)律，并根據(jù)塑性沉降量提出了控制沉降的方法。得到以下結(jié)論：

(1)在模擬區(qū)范圍較大的情況下，抽水過(guò)程產(chǎn)生了非穩(wěn)定的流場(chǎng)。近井處水位會(huì)迅速降低，遠(yuǎn)井處水位下降相對(duì)延遲；抽水停止后，滲流驅(qū)動(dòng)力為前期抽水形成的水力梯度，近井處水位恢復(fù)較快，而遠(yuǎn)井處水位恢復(fù)遲緩。地面沉降與水位變化規(guī)律具有一致性。

(2)在保證抽水總量一致的前提下，抽水時(shí)間短則抽水速率相對(duì)較大，在地面沉降歷史上引起的最大沉降量較大，最終塑性沉降量也相對(duì)較大；抽水次數(shù)多，則每次持續(xù)的時(shí)間短，同樣減小了歷史最大沉降量。三種持續(xù)抽水情景中，持續(xù)抽水60 d、100 d、150 d的歷史最大沉降量分別為24.1 cm、17.5 cm和13.9 cm，產(chǎn)生的塑性沉降量分別為4.0 cm、3.7 cm和3.4 cm。三種脈沖抽水情景中，脈沖間隔分別為10 d、20 d、40 d是，歷史最大沉降量分別為15.8 cm、17.1 cm和20.9 cm，塑性沉降量分別為3.7 cm、3.8 cm和4.1 cm。

(3)從控制地面沉降塑性沉降量的角度來(lái)看，持續(xù)抽水速率小，則含水層歷史最大沉降量小，最終塑性沉降量較小。脈沖抽水中，保持抽水速率一致，則脈沖間隔小引起的塑性沉降量較小?？傮w來(lái)看，地下水有補(bǔ)充源時(shí)慢速率、多次開(kāi)采能有效控制地面沉降的塑性沉降量。

[1] 薛禹群.論地下水超采與地面沉降[J].地下水,2012, 34(6):1-5.[XUE Y Q. Discussion on Groundwater overexploitation and ground settlement[J]. Ground Water, 2012, 34(6):1-5. (in Chinese)]

[2] 李培超. 地面沉降變形非線(xiàn)性完全耦合數(shù)學(xué)模型[J]. 河海大學(xué)學(xué)報(bào)(自然科學(xué)版), 2011(6):665-670.[LI P C. A nonlinear fully coupled mathematical model for land subsidence[J]. Journal of Hohai University(Natural Sciences), 2011, 39(6):665-670. (in Chinese)]

[3] 張弘懷,鄭銑鑫,唐仲華,等. 寧波平原地面沉降全耦合數(shù)值模擬研究[J]. 水文地質(zhì)工程地質(zhì),2013, 40(4):77-82 [ZHANG H H, ZHENG X X, TANG Z H,etal. A study of full coupling numerical simulation of land subsidencein the Ningbo plain[J]. Hydrogeology & Engineering Geology, 2013, 40(4):77-82. (in Chinese)]

[4] 陳興賢,駱祖江,金瑋澤,等. 高層建筑荷載、地下水開(kāi)采與地面沉降耦合研究[J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào),2015, 23(2):285-298.[CHEN X X, LUO Z J, JIN W Z,etal. Study of high-rise building load，groundwater seepage and land subsidence[J]. Journal of Basic Science and Engineering,2015, 23(2):285-298. (in Chinese)]

[5] Zhou X L, Huang K Y, Wang J H. Numerical simulation of groundwater flow and land deformation due to groundwater pumping in cross-anisotropic layered aquifer system[J]. Journal of Hydro-environment Research, 2017,14:19-33.

[6] Kim J. Sequential methods for coupled geomechanics and multiphase flow[J]. Dissertations & Theses-Gradworks. 2010.

[7] Rutqvist J, Wu Y S, Tsang C F,etal. A modeling approach for analysis of coupled multiphase fluid flow, heat transfer, and deformation in fractured porous rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2002,39(4):429-442.

[8] Rutqvist J, Tsang C-F. TOUGH-FLAC: a numerical simulator for analysis of coupled thermal-hydrologic-mechanical processes in fractured and porous geological media under multi-phase flow conditions[C]// Proceedings of the TOUGH Symposium. 2003.

[9] Rutqvist J. Status of the TOUGH-FLAC simulator and recent applications related to coupled fluid flow and crustal deformations[J]. Computers & Geosciences,2011,37(6): 739-750.

[10] Rutqvist J, Blanco Martin L, Mukhopadhyay S,etal. Modeling coupled THMC processes and brine migration in salt at high temperatures[R]. Berkeley:Ernest Orlando Lawrence Berkeley National Laboratory, 2014.

[11] Pruess K, Oldenburg C, Moridis G. TOUGH2 user’s guide version 2[M]. Berkeley :Lawrence Berkeley National Laboratorym,1999.

[12] 施小清, 張可霓, 吳吉春. TOUGH2軟件的發(fā)展及應(yīng)用[J]. 工程勘察,2009(10):29-34.[SHI X Q, ZHANG K N, WU J C. The history and application of TOUGH2 code[J]. Geotechnical Investigation &Surveying,2009(10):29-34.(in Chinese)]

[13] 蘇晨, 崔亞莉, 邵景力. 基于理論法的地面沉降數(shù)值模型進(jìn)展[J]. 水文地質(zhì)工程地質(zhì),2014, 41(6):147-152.[SU C, CUI Y L, SHAO J L. Advances in numerical models of land subsidence based on comprehensive theories[J]. Hydrogeology & Engineering Geology,2014, 41(6): 147-152. (in Chinese)]

[14] 楊林權(quán). 荷載和地下水開(kāi)采作用下的地面沉降數(shù)值模擬及應(yīng)用研究[D]. 武漢: 中國(guó)地質(zhì)大學(xué), 2014.[YANG L Q. Numerical simulation of land subsidence considering both effects of load and groundwater exploitation and applications[D]. Wuhan: China University of Geoscience, 2014. (in Chinese)]

[15] Helm D C. Horizontal aquifer movement in a Theis‐Thiem confined system[J]. Water Resources Research, 1994,30(4):953-964.

[16] 劉紅云. 抽水引起含水層水平運(yùn)動(dòng)及與地裂縫的關(guān)系研究[D]. 西安: 長(zhǎng)安大學(xué), 2007.[Horizontal aquifer movement induced from groundwater Pumping and its relation with ground fissure[D]. Xi’an: Chang’an University, 2007. (in Chinese)]

[17] 王哲成, 張?jiān)? 地下水超采引起的地裂縫災(zāi)害的研究進(jìn)展[J]. 水文地質(zhì)工程地質(zhì),2012, 39(2):88-93.[WANG Z C, ZHANG Y. Advances of research on the earth fissure induced bygroundwater over-exploitation[J]. Hydrogeology & Engineering Geology, 2012, 39(2):88-93. (in Chinese)]

[18] 陳育民, 徐鼎平. FLAC/FLAC3D基礎(chǔ)與工程實(shí)例[M]. 2版. 北京: 中國(guó)水利水電出版社; 2011.[CHEN Y M, XU D P. Theory and Engineering Application for FLAC/FLAC3D[M]. 2nd ed. Beijing: China Water & Power Press, 2011. (in Chinese)]

[19] ITASCA FD. Version 3.1 User’s Manual[M]. Minneapolis: Itasca Consulting Group, 2005.

[20] Davies J, Davies D. Stress-dependent permeability: characterization and modeling[C]//SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 1999.

[21] 湯連生, 桑海濤, 羅珍貴, 等. 土體抗拉張力學(xué)特性研究進(jìn)展[J]. 地球科學(xué)進(jìn)展,2015, 30(3):297-309.[TANG L S, SANG H T, LUO Z G. Advances in research on the mechanical behavior of the tensile strength of soils[J]. Advances in Earth Science,2015, 30(3):297-309. (in Chinese)]

[22] 陳崇希, 裴順平. 地下水開(kāi)采一地面沉降模型研究[J]. 水文地質(zhì)工程地質(zhì),2001,28(2):5-8.[CHEN C X, PEI S P. Groundwater extract and land subsidence model research[J]. Hydrogeology & Engineering Geology, 2001,28(2):5-8. (in Chinese)]

[23] 房浩,何慶成,徐斌,等.滄州地區(qū)地面沉降災(zāi)害風(fēng)險(xiǎn)評(píng)價(jià)研究[J].水文地質(zhì)工程地質(zhì), 2016,43(4):159-164. [FANG H, HE Q C, XU B,etal. A study of risk assessment of the land subsidence in Cangzhou[J].Hydrogeology & Engineering Geology,2016,43(4):159-164. (in Chinese)]

3D numerical simulation of elasto-plastic land subsidence induced by groundwater pumping

XIONG Xiaofeng, SHI Xiaoqing, WU Jianfeng, WU Jichun

(SchoolofEarthSciencesandEngineering,NanjingUniversity,Nanjing,Jiangsu210023,China)

In this paper, TOUGH2 and FLAC3Dare sequentially executed and linked to solve partial coupled three-dimensional land subsidence model, which considers the elastic-plastic deformation behavior, coupled hydro-mechanical effect and the varying nonlinear parameters. The evolution of land subsidence due to constant pumping and pulse pumping is compared to study the strategy of land subsidence control. Results show that: (1) Settlements rebound near the pumping well when stopping pumping, but the land subsidence cone continues enlarge. The plastic deformation can’t recover; (2) Pulse pumping results in fluctuation of the pore pressure, parameters and deformation, but they still dwindle globally; (3) The constant pumping strategy was compared with the pulse pumping with the same pumpage. Under this condition, the plastic deformation decrease with lower pumping rate and small intervals. This is mainly due to the fact that the groundwater level can recover rapidly with lower pumping rate and small intervals. This research provides a new method to simulate land subsidence, and the case study would contribute to the control of land subsidence.

TOUGH2-FLAC3D；land subsidence；elastic-plastic model；pulse pumping

10.16030/j.cnki.issn.1000-3665.2017.02.23

2016-10-08;

2017-01-11

國(guó)家自然科學(xué)基金項(xiàng)目(U1503282, 41672229，41172206)

熊小鋒(1991-)，男，碩士研究生，主要研究方向?yàn)榈孛娉两禂?shù)值模擬。E-mail：xiongxf_2016@163.com

施小清(1979-)，男，副教授，主要從事地下水?dāng)?shù)值模擬的教學(xué)和科研工作。E-mail：shixq@nju.edu.cn

P642.26

A

1000-3665(2017)02-0151-09