駱 勇，祝曉彬，郭 飛，吳吉春，蔣建國(guó)，曾獻(xiàn)奎，范亞民，王 棟

(1.表生地球化學(xué)教育部重點(diǎn)實(shí)驗(yàn)室/南京大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 210023；2.南京師范大學(xué)虛擬地理環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210046；3.江蘇省環(huán)境科學(xué)研究院，江蘇 南京 210036)

地下水開(kāi)發(fā)利用和疏排水過(guò)程引起的地面沉降問(wèn)題一直備受關(guān)注，抽排水引起的地面沉降問(wèn)題實(shí)際上就是一個(gè)滲流場(chǎng)與應(yīng)力場(chǎng)耦合的問(wèn)題。水流和地面沉降耦合模型按結(jié)合方式分為不耦合的兩步計(jì)算模型、部分耦合模型和全耦合模型。不耦合的兩步計(jì)算模型分為完全獨(dú)立的兩步，先計(jì)算孔隙水壓力，再計(jì)算變形，水流及沉降方程中各參數(shù)在沉降過(guò)程中不隨沉降過(guò)程發(fā)生變化，用于簡(jiǎn)單計(jì)算基坑抽排水引起的沉降計(jì)算的土力學(xué)經(jīng)驗(yàn)公式屬于此類。全耦合模型基于Biot固結(jié)理論，孔隙水壓力和變形同時(shí)算出，在理論上這種模型最符合沉降物理機(jī)制。部分耦合模型介于上述兩者之間，孔隙水壓力和變形既分步計(jì)算，兩者之間又相互影響[1]。葉淑君、施小清等[2～4]通過(guò)分析長(zhǎng)三角地區(qū)土層變形特征，采用修正的Merchant模型建立了三維變系數(shù)水流和垂向一維沉降組成的部分耦合模型，對(duì)長(zhǎng)三角區(qū)域的地面沉降進(jìn)行了模擬。熊小峰等[5]基于部分耦合原理，采用TOUGH2和FLAC3D建立抽水引起的三維地面沉降彈塑性模型。在滲流-沉降全耦合研究中，R W Lewis等[6]運(yùn)用比奧固結(jié)理論建立了威尼斯地面沉降模型。金瑋澤等[7]以比奧固結(jié)理論為基礎(chǔ)，建立了地下水疏降與地面沉降變形的水土全耦合模型，可同時(shí)求解地下水位和土體變形位移量，并將其與基于Terzaghi有效應(yīng)力原理地下水滲流與一維垂向固結(jié)模型進(jìn)行對(duì)比，取得結(jié)果更符合實(shí)際。全耦合模型相對(duì)部分耦合模型更符合實(shí)際情況、更為準(zhǔn)確，但是三維全耦合模型物理場(chǎng)控制方程更為復(fù)雜、計(jì)算量大、需要參數(shù)多，所以能進(jìn)行三維全耦合模型的建立并進(jìn)行求解的可視化軟件很少。近年來(lái)，基于偏微分方程設(shè)計(jì)專業(yè)有限元數(shù)值分析包的COMSOL Multiphysics被廣泛用于多場(chǎng)耦合問(wèn)題的模擬，只要是一個(gè)可以用偏微分方程形式數(shù)學(xué)模型描述的問(wèn)題，幾乎都可以采用COMSOL Multiphysics求解[8]。結(jié)合軟件具有的強(qiáng)大后處理功能，使其在機(jī)械制造、石油開(kāi)采等諸多領(lǐng)域應(yīng)用廣泛，但COMSOL Multiphysics應(yīng)用于求解疏排水引起的地面沉降研究相對(duì)較少[9～10]，用于計(jì)算地面沉降結(jié)果的合理性和可靠性缺乏驗(yàn)證。

基于此，本文基于不同滲透系數(shù)條件下的抽水引起地面沉降的理想算例，采用COMSOL Multiphysics進(jìn)行求解，并同時(shí)采用目前較為常用的不耦合的兩步計(jì)算方法-土力學(xué)經(jīng)驗(yàn)公式、部分耦合的地下水模擬系統(tǒng)軟件GMS中SUB模型進(jìn)行求解，通過(guò)結(jié)果對(duì)比驗(yàn)證COMSOL Multiphysics求解疏排水引起地面沉降的可行性和可靠性。

1 疏排水引起地面沉降計(jì)算理論與方法

1.1 土力學(xué)經(jīng)驗(yàn)公式法

土力學(xué)經(jīng)驗(yàn)公式為典型的滲流-地面沉降兩步計(jì)算模型。第一步計(jì)算疏排水引起的水位變化ΔH；第二步根據(jù)ΔH求出土層承載壓力變化ΔP，進(jìn)而求出地層壓縮量ΔB。計(jì)算含水層和弱透水黏性土層的土力學(xué)經(jīng)驗(yàn)公式[11]:

(1)

(2)

(3)

式中:ΔB1——含水層壓縮量/m；

ΔB2——黏性土層壓縮量/m；

ΔP——水位變化施加于土層的平均荷載/MPa；

M——計(jì)算土層厚度/m；

E——含水層彈性模量/MPa；

αv——壓縮系數(shù)；

e0——孔隙比；

ΔH——水位降深/m；

γw——水的重度/(N·m-3)。

1.2 GMS中SUB模型計(jì)算理論

地下水模擬系統(tǒng)軟件GMS中SUB模型為滲流-沉降部分耦合模型，水流和沉降模型控制方程：

(4)

b1=-ΔHS′skeΔb

(5)

b2=-ΔHS′skvΔb

(6)

式中：Ss——含水層儲(chǔ)水率/m-1；

K——滲透系數(shù)/(m·d-1)；

H——水頭/m；

Qm——源匯項(xiàng)/d-1；

b1——土體的彈性變形量/m；

b2——土體的非彈性變形量/m；

Δb——壓縮土層的厚度/m；

ΔH——水頭變化值/m。

兩場(chǎng)控制方程通過(guò)水頭項(xiàng)進(jìn)行耦合。數(shù)值計(jì)算中，先求解水流控制方程，然后將計(jì)算得到的水頭變化代入沉降控制方程計(jì)算土層變形。上述水流模型與沉降模型的耦合方式屬于分步耦合[12]，也可稱為部分耦合。

SUB模塊可以用垂向一維擴(kuò)散方程模擬較厚弱透水夾層隨時(shí)間變化的滯后排水過(guò)程：

(7)

(8)

由式(8)計(jì)算的夾層釋水時(shí)間參數(shù)遠(yuǎn)超過(guò)模擬時(shí)間步長(zhǎng)時(shí)，就必須考慮夾層的滯后釋水。

1.3 COMSOL Multiphysics模型計(jì)算理論

COMSOL Multiphysics采用滲流-沉降完全耦合模型，考慮滲流場(chǎng)變化對(duì)沉降影響，而沉降的變化又會(huì)影響滲流場(chǎng)，兩者之間相互耦合、相互影響。滲流場(chǎng)控制方程：

(9)

(10)

式中：K——滲透系數(shù)/(m·d-1)；

ux,uy,uz——x、y、z方向上的位移分量；

α——比奧固結(jié)系數(shù)；

γw——水的重度/(N·m-3)；

Ss——儲(chǔ)水率/m-1；

H——水頭/m；

εv——體應(yīng)變。

應(yīng)力場(chǎng)控制方程：

G=E/[2(1+ν)]

(11)

式中：uij——位移變量；

G——土體的剪切模量/MPa；

v——土體的泊松比；

E——楊氏模量/MPa；

ρf——水的密度/(kg·m-3)；

Hij——水頭/m。

(12)

Ss=nχf+(1-n)χp

(13)

式中：n0——初始孔隙度；

K0——初始滲透系數(shù)/(m·d-1)；

χf——多孔介質(zhì)骨架壓縮系數(shù)；

χp——孔隙壓縮系數(shù)。

2 疏排水引起地面沉降算例模型

算例模型為一單井抽水引起地面沉降模型。含水系統(tǒng)為二層結(jié)構(gòu)，表層為弱透水的黏土層(層厚20 m)，下層為含水砂層(層厚30 m)，水平范圍是以抽水井為中心的500 m×500 m區(qū)域。含水系統(tǒng)結(jié)構(gòu)見(jiàn)圖1。

初始條件：設(shè)定算例研究區(qū)初始水位標(biāo)高為45 m。

邊界條件：設(shè)置研究區(qū)上、下邊界為零通量邊界，研究區(qū)四周邊界為給定水頭(45 m)邊界。

基本參數(shù)：含水砂層水平滲透系數(shù)取Kx、Ky設(shè)為4 m/d、垂向滲透系數(shù)Kz取0.4 m/d，黏土層Kx、Ky取0.004 m/d、Kz取0.000 4 m/d；含水層孔隙度n0取0.3，黏土層孔隙度n0取0.6。

模擬情景：基于模型，考慮單井抽水的情形(抽水量為1 000 m3/d)，設(shè)置兩個(gè)時(shí)間應(yīng)力期，抽水應(yīng)力期時(shí)長(zhǎng)50 d(第0～50天)，停抽恢復(fù)應(yīng)力期時(shí)長(zhǎng)200 d(第51～250天)。

計(jì)算結(jié)果：對(duì)含水系統(tǒng)的水位變化及引起的沉降量進(jìn)行計(jì)算。

圖1 理想算例含水系統(tǒng)結(jié)構(gòu)Fig.1 A schematic diagram of the structure of an ideal example water-bearing system

3 不同計(jì)算方法求解疏排水引起的地面沉降

使用土力學(xué)經(jīng)驗(yàn)公式法、GMS中SUB模型以及COMSOL Multiphysics分別計(jì)算上述單井抽水模擬情景下引起的地面沉降。

3.1 經(jīng)驗(yàn)公式法求解

采用經(jīng)驗(yàn)公式計(jì)算單井抽水引起不同水位變幅時(shí)的地面沉降值。降水面以下的土層通常不產(chǎn)生明顯固結(jié)壓縮量，抽水產(chǎn)生的地面沉降主要由最終降水面至原始地下水面之間土層變形產(chǎn)生[11]。本文算例模型中黏土層厚20 m，下層承壓含水層厚30 m，初始水位45 m，當(dāng)水位降低幅度在15 m以內(nèi)時(shí)只計(jì)算黏土層壓縮量。采用經(jīng)驗(yàn)公式計(jì)算不同水位變幅下的土層沉降量涉及有關(guān)參數(shù)取值見(jiàn)表1，水位下降引起的地面沉降量計(jì)算結(jié)果見(jiàn)圖2，沉降和水位下降同步發(fā)生，沉降量隨水位降深增大而逐漸增大。

圖2 經(jīng)驗(yàn)公式計(jì)算不同水位變幅下的土層沉降量Fig.2 Analytical solution under different water levels under the soil settlement

土層孔隙比e0壓縮系數(shù)α0/kPa-1計(jì)算土層厚度M/m粘土層1.50.1與水位變動(dòng)幅度一致

3.2 GMS中SUB模型求解

采用GMS中SUB模型求解算例模型抽水引起地面沉降的有關(guān)參數(shù)取值[16～17]見(jiàn)表2。在模型中，承壓含水層未分配滯后夾層，將弱透水層進(jìn)行多層劃分。計(jì)算得到距抽水井不同距離處含水層水位變化與地面沉降量變化對(duì)比見(jiàn)圖3。距抽水井不同距離處(10 m、25 m、50 m)，地面沉降量隨水位變化規(guī)律基本一致，水位下降導(dǎo)致沉降發(fā)生，水位降深越大、沉降越明顯，水位回升導(dǎo)致地面沉降逐步恢復(fù)，水位回升至初始狀態(tài)時(shí)，沉降并未完全恢復(fù)，存在一個(gè)不可恢復(fù)的永久沉降量。在0～50 d時(shí)抽水井抽水，各處水位都在下降，各處沉降量隨著水位的下降而逐漸增大。距抽水井越近水位下降幅度越大，沉降量越大。距抽水井越遠(yuǎn)水位下降幅度越小，沉降量越小。51～250 d時(shí)抽水井停止抽水，距抽水井不同距離處水位立即同時(shí)回升，在第120 d左右水位恢復(fù)到初始水位。根據(jù)前人相關(guān)研究[17～19]，SUB沉降模型可以模擬疏水時(shí)弱透水夾層的沉降滯后效應(yīng)，但在本文模型設(shè)置情景下，沉降滯后現(xiàn)象不明顯。

圖3 GMS SUB模型計(jì)算下距抽水井不同距離處水位及沉降量變化對(duì)比圖Fig.3 Change in the water level and settlement at different distances from the pumping well calculated with the GMS SUB module

土層彈性骨架儲(chǔ)水率S′ske/m-1非彈性骨架儲(chǔ)水率S′skv/m-1彈性骨架儲(chǔ)水系數(shù)非彈性骨架儲(chǔ)水系數(shù)垂向滲透系數(shù)K′v/(m·d-1)含水砂層——3.4×10-43.4×10-2—粘土層4.0×10-54.2×10-3——4.0×10-6

3.3 COMSOL Multiphysics求解

采用COMSOL Multiphysics求解算例模型抽水引起地面沉降的有關(guān)參數(shù)取值[20]見(jiàn)表3，含水層中距抽水井不同距離處(10 m、25 m、50 m)水位變化與沉降量變化對(duì)比見(jiàn)圖4，不同距離處水位和沉降量總體變化規(guī)律一致。在0～50 d時(shí)抽水井抽水，各處水位都在下降，各處沉降量隨著水位的下降逐漸增大。距抽水井越近水位下降幅度越大，產(chǎn)生的沉降量越大；距抽水井越遠(yuǎn)水位下降幅度越小，沉降量越小。51～150 d時(shí)抽水井停止抽水，距抽水井不同距離處(10 m、25 m、50 m)水位立即產(chǎn)生較大回升然后逐漸恢復(fù)到初始狀態(tài)，但距抽水井不同距離處(10 m、25 m、50 m)沉降量變化出現(xiàn)滯后現(xiàn)象。距抽水井越遠(yuǎn)處，沉降滯后時(shí)間越長(zhǎng)，距抽水井10 m處發(fā)生沉降的趨勢(shì)在第60 d開(kāi)始停止，隨著水位的恢復(fù)地面發(fā)生回彈，沉降量逐漸變小，在模擬期結(jié)束后仍有約0.63 mm的沉降量。距抽水井25 m處發(fā)生沉降的趨勢(shì)在第64 d開(kāi)始停止，隨著水位的恢復(fù)地面發(fā)生回彈，在模擬期結(jié)束后仍有一個(gè)0.5 mm左右的沉降量。距抽水井50 m處發(fā)生沉降的趨勢(shì)在第70 d開(kāi)始停止，隨著水位的恢復(fù)地面發(fā)生回彈，沉降量逐漸變小，在模擬期結(jié)束后仍有約0.24 mm的沉降量。距離抽水井越遠(yuǎn)，產(chǎn)生的最大沉降量越小、模擬期結(jié)束時(shí)的沉降量越小，但滯后于水位變化的時(shí)間卻越長(zhǎng)。

圖4 COMSOL Multiphysics計(jì)算下距抽水井不同距離處水位及沉降量變化對(duì)比圖Fig.4 Comparison of water level and settlement variation at different distances from the pumping well calculated with the COMSOL Multiphysics

土層泊松比ν0密度ρ/(kg·m-3)楊氏模量E/MPa儲(chǔ)水率S/m-1含水層0.4892 150400.008黏土層0.4982 000500.000 04

4 地面沉降計(jì)算結(jié)果分析與討論

4.1 計(jì)算結(jié)果對(duì)比分析

4.1.1COMSOL Multiphysics與經(jīng)驗(yàn)公式法計(jì)算結(jié)果對(duì)比

計(jì)算應(yīng)力期內(nèi)，由COMSOL Multiphysics計(jì)算得到距抽水井10 m、25 m、50 m處最大水位變幅分別為4.3 m、3.1 m、2.2 m。在相同水位變幅條件下，用經(jīng)驗(yàn)公式和COMSOL Multiphysics計(jì)算出距抽水井不同距離處對(duì)應(yīng)最大沉降量對(duì)比見(jiàn)表4。不同水位變幅下，COMSOL Multiphysics與經(jīng)驗(yàn)公式最大沉降量計(jì)算值雖有偏差但較為接近，說(shuō)明本文建立的COMSOL Multiphysics模型可靠。

表4 相同水位變幅下經(jīng)驗(yàn)公式與COMSOL Multiphysics沉降量計(jì)算結(jié)果對(duì)比Table 4 Comparison of the settlement calculation results under the same water level variation calculated with the empirical formula and COMSOL Multiphysics

4.1.2COMSOL Multiphysics與GMS中SUB計(jì)算結(jié)果對(duì)比

距抽水井25 m處GMS的SUB模型與COMSOL Multiphysics水位變化、沉降量變化計(jì)算值對(duì)比見(jiàn)圖5，結(jié)果表明2種方法計(jì)算沉降量值總體接近。在本文模型設(shè)置情景下，GMS中SUB模型計(jì)算結(jié)果的水位與沉降幾乎同步變化，沉降變化滯后現(xiàn)象不明顯；而COMSOL Multiphysics中沉降滯后現(xiàn)象相比前者較為明顯，且圖4表明距抽水井越遠(yuǎn)沉降滯后越明顯。在抽水初期，不同模型計(jì)算水位呈現(xiàn)大致相同幅度的降低，但GMS中SUB模型相比COMSOL Multiphysics沉降量增大速率和幅度明顯，表明GMS中SUB模型沉降量對(duì)水位變化響應(yīng)快。沉降回彈時(shí)，COMSOL Multiphysics計(jì)算的沉降量回彈過(guò)程較為緩慢，而GMS中SUB模型計(jì)算的沉降在抽水停止后瞬間產(chǎn)生一個(gè)較大幅度的回彈。對(duì)比實(shí)際沉降對(duì)水位變化的響應(yīng)規(guī)律，COMSOL Multiphysics計(jì)算結(jié)果更符合實(shí)際沉降特征。這是因?yàn)镚MS中SUB模型是滲流-沉降部分耦合模型，通過(guò)先計(jì)算地下水滲流方程，求出水頭變化后，即由孔隙水壓力變化計(jì)算地面沉降，未將孔隙水壓力變化與沉降同時(shí)考慮。COMSOL Multiphysics中充分考慮了滲流-沉降場(chǎng)的耦合，不僅考慮了孔隙度、滲透率等土體參數(shù)的動(dòng)態(tài)變化，還考慮了土體的黏彈性、黏塑性，抽水引起滲流場(chǎng)改變，滲流場(chǎng)改變使得孔隙水壓力發(fā)生變化，土體將產(chǎn)生沉降，而土體的沉降又會(huì)改變土體的孔隙度、滲透率從而影響滲流場(chǎng)，上述過(guò)程即實(shí)現(xiàn)了流固全耦合。

圖5 距抽水井25 m處GMS SUB模塊及COMSOL Multiphysics計(jì)算下水位及沉降量變化對(duì)比圖Fig.5 Comparison of water level and settlement change at the distance of well 25 m from the pumping well calculated with the GMS SUB module and COMSOL Multiphysics

4.2 滲透系數(shù)K對(duì)地面沉降量計(jì)算結(jié)果的影響討論

為進(jìn)一步研究不同模型計(jì)算疏排水引起地面沉降的差異，考慮不同滲透系數(shù)對(duì)地面沉降計(jì)算結(jié)果的影響。因?yàn)闈B透系數(shù)對(duì)流場(chǎng)影響較大，而流場(chǎng)中水頭的變化是導(dǎo)致沉降的直接因素，因此滲透系數(shù)大小對(duì)沉降量影響同樣較大。采用不同方法計(jì)算不同滲透系數(shù)條件下抽水引起的地面沉降，可以進(jìn)一步對(duì)方法的有效性進(jìn)行比較。GMS中SUB模型的滲透系數(shù)K為一定值，雖然也有學(xué)者對(duì)SUB模塊這一不足進(jìn)行了研究，如施小清等[21]針對(duì)黏性土壓縮過(guò)程中各土體參數(shù)都發(fā)生變化的問(wèn)題，對(duì)MODFLOW-2000中的SUB模塊進(jìn)行了修正，較好刻畫(huà)了地面沉降過(guò)程中黏性土釋水的特點(diǎn)。但相比COMSOL Multiphysics從機(jī)理上實(shí)現(xiàn)完全的流固耦合，無(wú)疑COMSOL Multiphysics具有更好應(yīng)用條件。

對(duì)不同初始滲透系數(shù)條件下，距抽水井25 m處，GMS中SUB模型和COMSOL Multiphysics求解的沉降結(jié)果進(jìn)行比較。滲透系數(shù)逐漸增大(Kx、Ky由4 m/d增大到10 m/d，再增大到30 m/d；Kz由0.4 m/d增大到1 m/d，再增大到3 m/d)，地面沉降量計(jì)算結(jié)果見(jiàn)圖6。結(jié)果表明，算例條件下，兩種方法計(jì)算結(jié)果均反應(yīng)出滲透系數(shù)越大沉降量越小的現(xiàn)象，滲透系數(shù)越大，相同抽水條件下水位變化越小，和實(shí)際含水系統(tǒng)中對(duì)應(yīng)的水位和沉降特征吻合。COMSOL Multiphysics計(jì)算的沉降量滯后時(shí)間隨著滲透系數(shù)的變大而減小，GMS中 SUB模型計(jì)算的沉降量在本文模擬情景下滯后效應(yīng)不明顯，在抽水停止時(shí)會(huì)發(fā)生較大回彈，并且隨著滲透系數(shù)的變大趨于穩(wěn)定的時(shí)間變快；在抽水初期，GMS中SUB模型沉降量計(jì)算值增幅較COMSOL Multiphysics更大。兩種方法計(jì)算結(jié)果的差異性是由方法理論控制方程的差異導(dǎo)致的：完全耦合模型以比奧固結(jié)理論為基礎(chǔ)，對(duì)土體變形與孔隙水壓力同時(shí)進(jìn)行求解，將沉降模型與水流模型統(tǒng)一于同一物理空間，并且不作固結(jié)過(guò)程中總應(yīng)力為常量的假設(shè)；部分耦合模型考慮了土體僅垂向上發(fā)生變形的情況，假設(shè)總應(yīng)力不隨時(shí)間變化，且滲透系數(shù)等土體參數(shù)不變。完全耦合的COMSOL Multiphysics模型計(jì)算結(jié)果更為合理。滲透系數(shù)越小、沉降越易發(fā)生的地層，COMSOL Multiphysics模型計(jì)算的沉降過(guò)程滯后越明顯，水位回升后沉降越難以恢復(fù)，結(jié)果比GMS中SUB模型計(jì)算結(jié)果越合理。

圖6 不同滲透系數(shù)下距抽水井25 m處GMS SUB模塊與COMSOL Multiphysics沉降量計(jì)算值對(duì)比圖Fig.6 Comparison of settlement at the distance of 25 m from the pumping well at different coefficients of permeability with the GMS SUB module and COMSOL Multiphysics

5 結(jié)論

分析比較不同方法計(jì)算抽水引起的沉降量計(jì)算結(jié)果和反應(yīng)的沉降過(guò)程，得到如下結(jié)論：

(1)COMSOL Multiphysics與經(jīng)驗(yàn)公式、GMS中SUB模型計(jì)算算例條件下的沉降量計(jì)算結(jié)果總體較為接近，表明COMSOL Multiphysics求解疏排水引起的地面沉降可行。

(2)在本文模擬情景下，COMSOL Multiphysics計(jì)算的沉降量出現(xiàn)明顯滯后于水位變化的規(guī)律，且距抽水井越遠(yuǎn)沉降滯后時(shí)間越長(zhǎng)；而在相同模擬情景下GMS SUB模塊沉降滯后現(xiàn)象不明顯?；跐B流-應(yīng)力場(chǎng)全耦合的COMSOL Multiphysics模型計(jì)算結(jié)果更合理。

(3)與GMS中SUB模型計(jì)算結(jié)果相比，滲透系數(shù)越小，COMSOL Multiphysics模型計(jì)算結(jié)果更能反應(yīng)實(shí)際沉降過(guò)程特征。