基于數值模擬的地下水開采引起的地面沉降機理和對策研究
----以定襄縣賈家莊為例

盛登寶，呂義清

(太原理工大學 地球科學與工程系，太原 030024)

0 引 言

過度抽取地下水引起的地表不均勻沉降，影響地表建筑物、道路、管網等設施的安全，所以，開展地面沉降的研究，分析其產生沉降的原因以及動態降水的時間與空間效應具有重大的現實意義。抽取地下水過程中引起的固結變形是一個典型的三維滲流固結耦合問題，國內外許多學者開展了大量的理論和實踐研究。

研究表明，地下水在含水層孔隙內流動時，導致孔隙水壓力的變化，會引起與土骨架變形有關的應力變化，同時改變巖土體的物理力學性質[1,2]。前人在巖土體流固耦合分析方面已經做了大量研究[3-6]，本文借助Comsol Multiphysics有限元數值模擬軟件，基于biot固結理論，結合彈塑性力學理論和滲透力學理論，建立數值模型進行求解，得出了大量抽取地下水情況下地表變形位移情況以及含水層水頭、土體孔隙水壓力分布情況，為地質災害預防提供一定的幫助和技術指導。

1 固結數學模型

1.1 基本假設

(1)土骨架線彈性，變形微??；

(2)土顆粒和孔隙水均不可壓縮；

(3)滲流符合達西定律；

滲流場變化過程則主要基于達西定律形式的連續方程來描述，即：

▽·(-K▽H)=0

(1)

式中：Sh為存儲系數，m-1；K為滲透系數，m/s；H為水頭，m。

將方程(1)中水頭用孔壓表示進而轉化為方程(2):

(2)

式中：S=Sh/(ρg)；H=P/(ρg)+D；K=κρg/μ。

1.2 基本方程

(1)平衡方程。假設一均質，各向同性的飽和土單元體dxdydz，若體力只考慮重力，z坐標向上為正，以土體為隔離體則三維平衡微分方程為：

(3)

(2)本構方程。Biot理論最初假設土骨架為線彈性體，服從廣義胡克定律，根據彈性力學本構方程，應力用應變來表示：

(4)

式中：G、υ分別為剪切模量和泊松比；ευ為體應變ευ=εx+εy+εz。

(3)幾何方程。用幾何方程將應變表示成位移，設x、y、z方向的位移為us、vs、ws在小變形的假定下，6個應變分量為：

(5)

式中：εx,εy,εz為x、y、z方向的正應變。

(4)固結微分方程。將(4)和(5)帶入(3)中就可以得到以位移和空隙壓力表示的平衡微分方程：

(6)

上面方程中有4個未知量us,vs,ws,u，求解還需要一個方程，由達西定律得：

▽2u

(7)

展開用位移表示得：

▽2u=0

(8)

式中：K為滲流系數；γw為水的容重。

公式(6)、(8)便是Biot固結方程。

2 應用實例

2.1 工程概況

賈家莊村位于山西省定襄縣西北部，村莊東側有1眼抽水井，井深130 m，井口直徑30 cm。地下水位埋深較淺，含水層為泥河灣組、午城組、離石組黃土層中的礫石層、流沙層及鈣質結核層、中—粗砂層，地下水主要接受河水及大氣降水補給。不良地質現象主要為由于地下水的過量抽取引起了沉積層固結壓縮，引起地表沉降變形，最后造成居民房屋嚴重破壞，見圖1。

圖1 賈家莊地面沉降房屋變形破壞圖

2.2 模型的建立

根據該礦區的地層資料，建立長為1 000 m、寬為1 000 m、高為130 m的三維模型，共三層，厚度從上至下依次為40、10、80 m。在抽水井的設定中，本文采用桿單元來描述抽水井，它能精細的控制實際流量，很好的反映動態降水的時間與空間效應，與采用定水頭的點模擬抽水井的方法相比，此處理方法考慮了尺寸效應對模型的影響。

2.3 邊界、初始條件、參數

結合定襄縣水文地質資料，確定了巖土體物理力學參數(表1所示)。模型中流體的密度為1 000 kg/m3；黏度為1×10-3Pa·s；流體壓縮率為4×10-101/Pa。通過從最低含水層抽取地下水來降低含水層，抽水井日出水量約8 200 m3/d，所以本模型直接將抽水邊界的水頭定為H=-3 (m/a)*t，即假設水以平均3 m/a的滲透速度流動，且滲流場從初始即進入穩定狀態，10年時間保持不變。而影響模型分析計算的主要因素是水頭的變化而不是水頭值自身，所以本模型直接將表層的水頭定為H0=0 m，這樣處理可避免隨時間變化要不斷變化表層水頭邊界的麻煩。這種處理方法便于建立初始邊界條件。結合實際監測情況，在計算模型中距抽水井0、250、500 m位置分別添加3個域點探針來模擬實際監測。

表1 巖層物理力學參數

2.4 計算結果及分析

地面沉降是土層中孔隙水承擔的孔隙水壓力和土骨架承擔的有效應力發生變化的結果。處于平衡狀態的含水系統，當地下水被抽出后，孔隙水壓力減小，原先的土、水平衡狀態被破壞，有效應力發生變化土體產生變形，且土體的力學性質、貯水性和透水性都將隨之變化。地面沉降是土和水相互作用、內部應力發生變化的外在表現，它與土的變形特性和水的滲流情況密切相關。

多物理場耦合分析軟件Comsol Multiphysics中的多孔彈性物理場接口建立的地下水抽取—地面沉降流固耦合模型，實現了達西定律場和固體力學場的真耦合，能夠精確反映動態降水的時間與空間效應；通過多孔彈性物理場接口，建立流固耦合模型，模型的求解時間為10 a，1 a為一個時間步。

經過計算分析得出地面及地下水時空效應變化。如圖2所示，該結果為10 a(現如今)后地面沉降結果，抽水后，地下水水位下降，在抽水井周圍產生較大的水力坡度，隨著水力坡度的增大，相應地滲透壓力也增大。在地下水自上而下的滲透過程中，地表變形，形成降落漏斗，表現為在地面呈以降水井為圓心的同心圓分布，影響半徑約300 m。圖3為t=10 a地面沉降二維剖面圖，礦井抽水區域漏斗中心沉降最大為Smax=1.19 m。

圖2 t=10 a，Z方向地面沉降位移圖

圖3 t=0～10 a，地面沉降總位移二維剖圖

由于土體是由固體顆粒(固相)、水(液相)和空氣(氣相)三部分組成，作用在土體上的應力是通過顆粒間的接觸和孔隙水來傳遞的。由顆粒間接觸點傳遞的應力是對土體變形和強度有效的粒間力，稱為有效應力。由孔隙水傳遞的應力稱為孔隙水壓力，它僅對土顆粒產生壓縮，由于固體顆粒的壓縮模量非常大，可以認為是不可壓縮的，工程上經常忽略不計，故認為不能直接引起土體的變形和強度變化，所以又稱為中性應力[7]。因此，土體在垂直方向所受的總應力σ為有效應力σ′與孔隙水壓力u之和。這就是有名的有效應力原理，當水位下降到一定深度后，土中的孔隙水壓力降低，由于總壓力基本不變，則有效應力相應增大，如圖3、圖4所示，t=10 a時，抽水井附近有效應力較大，最大為4.14×105Pa，水頭最大為Hmax=-30 m。

圖4 t=0～10 a，流線分布與Von Mises應力云圖

通過計算得到模型中各監測點近10年的地表沉降量(圖5)，根據計算結果，在降水固結過程的初期，由于有效應力的增加，引起土體固結壓縮，且隨著時間的推移而發生，呈線性變化；靠近抽水井附近的監測點1的10年累計下沉值最大，達到1.12 m，監測點2的累計下沉值為0.72 m，監測點3的累計下沉值最小，為0.41 m；并將其分別與實際的監測結果進行擬合對比，整體效果較好，則可以認定所建的數值模型正確，參數較合理，能很好地反映研究區的沉降特征。

圖5 t=0～10 a，水頭和流速分布圖

通過計算不同開采工況得到t=10 a時地表沉降量(圖6)，根據計算結果，在H1=-0.5t(日出水量Q=1 400 m3)的工況條件下，地表變形較小，根據《建筑地基基礎設計規范》，不會引起墻體及房屋的破壞，所以，將日開采量控制在1 400 m3以下能有效避免地面沉降造成的破壞。

3 結 論

(1)采用多物理場耦合分析軟件Comsol Multiphysics中的多孔彈性物理場接口建立的地下水抽取-地面沉降流固耦合模型，實現了達西定律場和固體力學場的真耦合，能夠精確反映地表變形位移以及含水層水頭、土體孔隙水

圖6 研究區各監測點沉降擬合圖

壓力分布情況。

(2)通過分析計算，確定了地下水抽取引起的地面變形的影響范圍；根據有效應力原理，當水位下降到一定深度后，土中的孔隙水壓力降低，由于總壓力基本不變，則有效應力相應增大，引起土體固結壓縮并得出t=10 a時，抽水井附近最大效應力σ′max=4.14×105Pa，水頭最大值為Hmax=-30 m。

(3)土體固結壓縮隨著時間的推移而發生，呈線性變化；并對位移的計算結果與監測結果進行擬合對比分析，認定所建的數值模型正確，參數較合理，能很好地反映研究區的沉降特征。

(4)通過計算不同開采工況，根據計算結果，給出最佳的沉降控制措施：將日開采量控制在1 400 m3以下，有效避免地面沉降造成的破壞。

□

[1] 紀佑軍,劉建軍,程林松.考慮流-固耦合的隧道開挖數值模擬[J].巖土力學,2011,32(4):1 229-1 233.

[2] 陳仕闊,楊天鴻,王 航,等.高富水砂層地鐵豎井動態降水優化及固結變形預測[J].東北大學學報(自然科學版),2011,32(9):1 328-1 331.

[3] 陳 杰,朱國榮,顧阿明,等.Biot固結理論在地面沉降計算中的應用[J].水文地質工程地質,2003,30(2):28-31.

[4] 張 云,薛禹群.抽水地面沉降數學模型的研究現狀與展望[J].中國地質災害與防治學報,2002,13(2):1-6.

[5] 李元雄.李元雄.抽汲地下水引起地面沉降三維數值模擬及工程應用研究[D]. 武漢：中國地質大學,2007.

[6] 唐金穎.地下水開采引起地面沉降的數值模擬分析[D]. 沈陽：東北大學, 2008.

[7] 李志平.基坑降水引起的地面沉降分析[D]. 長沙：中南大學,2008.