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

盛登寶，呂義清

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

0 引 言

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

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

1 固結(jié)數(shù)學(xué)模型

1.1 基本假設(shè)

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

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

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

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

▽·(-K▽H)=0

(1)

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

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

(2)

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

1.2 基本方程

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

(3)

(2)本構(gòu)方程。Biot理論最初假設(shè)土骨架為線彈性體，服從廣義胡克定律，根據(jù)彈性力學(xué)本構(gòu)方程，應(yīng)力用應(yīng)變來表示：

(4)

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

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

(5)

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

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

(6)

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

▽2u

(7)

展開用位移表示得：

▽2u=0

(8)

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

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

2 應(yīng)用實(shí)例

2.1 工程概況

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

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

2.2 模型的建立

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

2.3 邊界、初始條件、參數(shù)

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

表1 巖層物理力學(xué)參數(shù)

2.4 計(jì)算結(jié)果及分析

地面沉降是土層中孔隙水承擔(dān)的孔隙水壓力和土骨架承擔(dān)的有效應(yīng)力發(fā)生變化的結(jié)果。處于平衡狀態(tài)的含水系統(tǒng)，當(dāng)?shù)叵滤怀槌龊螅紫端畨毫p小，原先的土、水平衡狀態(tài)被破壞，有效應(yīng)力發(fā)生變化土體產(chǎn)生變形，且土體的力學(xué)性質(zhì)、貯水性和透水性都將隨之變化。地面沉降是土和水相互作用、內(nèi)部應(yīng)力發(fā)生變化的外在表現(xiàn)，它與土的變形特性和水的滲流情況密切相關(guān)。

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

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

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

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

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

圖4 t=0～10 a，流線分布與Von Mises應(yīng)力云圖

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

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

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

3 結(jié) 論

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

圖6 研究區(qū)各監(jiān)測點(diǎn)沉降擬合圖

壓力分布情況。

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

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

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

□

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

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

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

[4] 張 云,薛禹群.抽水地面沉降數(shù)學(xué)模型的研究現(xiàn)狀與展望[J].中國地質(zhì)災(zāi)害與防治學(xué)報(bào),2002,13(2):1-6.

[5] 李元雄.李元雄.抽汲地下水引起地面沉降三維數(shù)值模擬及工程應(yīng)用研究[D]. 武漢：中國地質(zhì)大學(xué),2007.

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

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