駱 勇，祝曉彬，郭 飛，吳吉春，蔣建國，曾獻奎，范亞民，王 棟

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

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

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

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

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

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

(1)

(2)

(3)

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

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

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

M——計算土層厚度/m；

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

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

e0——孔隙比；

ΔH——水位降深/m；

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

1.2 GMS中SUB模型計算理論

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

(4)

b1=-ΔHS′skeΔb

(5)

b2=-ΔHS′skvΔb

(6)

式中：Ss——含水層儲水率/m-1；

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

H——水頭/m；

Qm——源匯項/d-1；

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

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

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

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

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

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

(7)

(8)

由式(8)計算的夾層釋水時間參數(shù)遠超過模擬時間步長時，就必須考慮夾層的滯后釋水。

1.3 COMSOL Multiphysics模型計算理論

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

(9)

(10)

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

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

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

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

Ss——儲水率/m-1；

H——水頭/m；

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

應(yī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)見圖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è)置兩個時間應(yīng)力期，抽水應(yīng)力期時長50 d(第0～50天)，停抽恢復(fù)應(yīng)力期時長200 d(第51～250天)。

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

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

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

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

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

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

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

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

3.2 GMS中SUB模型求解

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

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

土層彈性骨架儲水率S′ske/m-1非彈性骨架儲水率S′skv/m-1彈性骨架儲水系數(shù)非彈性骨架儲水系數(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]見表3，含水層中距抽水井不同距離處(10 m、25 m、50 m)水位變化與沉降量變化對比見圖4，不同距離處水位和沉降量總體變化規(guī)律一致。在0～50 d時抽水井抽水，各處水位都在下降，各處沉降量隨著水位的下降逐漸增大。距抽水井越近水位下降幅度越大，產(chǎn)生的沉降量越大；距抽水井越遠水位下降幅度越小，沉降量越小。51～150 d時抽水井停止抽水，距抽水井不同距離處(10 m、25 m、50 m)水位立即產(chǎn)生較大回升然后逐漸恢復(fù)到初始狀態(tài)，但距抽水井不同距離處(10 m、25 m、50 m)沉降量變化出現(xiàn)滯后現(xiàn)象。距抽水井越遠處，沉降滯后時間越長，距抽水井10 m處發(fā)生沉降的趨勢在第60 d開始停止，隨著水位的恢復(fù)地面發(fā)生回彈，沉降量逐漸變小，在模擬期結(jié)束后仍有約0.63 mm的沉降量。距抽水井25 m處發(fā)生沉降的趨勢在第64 d開始停止，隨著水位的恢復(fù)地面發(fā)生回彈，在模擬期結(jié)束后仍有一個0.5 mm左右的沉降量。距抽水井50 m處發(fā)生沉降的趨勢在第70 d開始停止，隨著水位的恢復(fù)地面發(fā)生回彈，沉降量逐漸變小，在模擬期結(jié)束后仍有約0.24 mm的沉降量。距離抽水井越遠，產(chǎn)生的最大沉降量越小、模擬期結(jié)束時的沉降量越小，但滯后于水位變化的時間卻越長。

圖4 COMSOL Multiphysics計算下距抽水井不同距離處水位及沉降量變化對比圖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儲水率S/m-1含水層0.4892 150400.008黏土層0.4982 000500.000 04

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

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

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

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

表4 相同水位變幅下經(jīng)驗公式與COMSOL Multiphysics沉降量計算結(jié)果對比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計算結(jié)果對比

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

圖5 距抽水井25 m處GMS SUB模塊及COMSOL Multiphysics計算下水位及沉降量變化對比圖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對地面沉降量計算結(jié)果的影響討論

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

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

圖6 不同滲透系數(shù)下距抽水井25 m處GMS SUB模塊與COMSOL Multiphysics沉降量計算值對比圖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é)論

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

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

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

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