基于有限層法的承壓水開采地面沉降簡明分析

徐 進,張 豪,王斯海

(1.煙臺大學土木工程學院,山東 煙臺 264005;2.江蘇工程職業技術學院建筑工程學院,江蘇 南通 226000)

對含水層中地下水進行開采會破壞原有的水土平衡狀態,在地下水位發生降深的同時,引發土層變形,導致地面沉降。同時,人類工程活動中對潛水和承壓含水層的降水減壓也會導致工程性沉降,對周圍建筑物和市政設施造成不利影響。如何正確評估上述環境地質問題日益得到學術界的重視[1-2]。

從機理來說,地面沉降是一個復雜的流固耦合過程[3]?；谒?土耦合理論,LEWIS和SCHREFLER較早提出了地面沉降有限元計算模型,并運用于威尼斯的地面沉降問題[4]。HSI等考慮了潛水面變化,研究了抽水過程中潛水面變動對地面沉降的影響[5]。金小榮等利用二維有限元模型分析了工程降水作用下基坑周圍土體沉降性狀的影響因素[6]。駱祖江等系統建立了地下水開采引起地面沉降三維全耦合有限元模型,引入了滲透系數動態模型和鄧肯-張非線性模型[7]。鄭剛等利用有限元軟件建立三維流固耦合模型研究了承壓層長期局部降壓的土體分層沉降機理和規律[8]。ZHOU等基于半解析數值原理給出了抽水作用下三維流固耦合問題的求解格式[9]。ZENG等建立了基坑降水地面沉降的耦合模型,分析了多層地下水系統基坑降水引起的支護結構和周圍土體響應[10]。

但是,基于水-土耦合的地面沉降計算模型,由于計算工作量大、參數難以獲取,在實際應用中仍受到一定限制。鑒于此,便于實際應用的簡化計算方法得到較多關注。駱冠勇等基于一維太沙基固結方程,推導了承壓水降壓沉降的固結度計算公式[11]。龔曉南等結合地下水動力學公式與彈性半空間理論,給出了承壓水穩定流情形下的最大沉降公式[12]。為了進一步反映抽水時間的影響,王春波等將其推廣到了非穩定流情形[13]。

地下水流計算是地面沉降分析的基礎環節,如何實現三維地下水流的高效計算一直備受關注。有限層法是一種典型的混合算法,兼具解析方法和數值方法的優點,其精確高效的特點在地下水流問題中已得到充分體現[14-15]。為此,基于地下三維流高效有限層方法,針對承壓水開采引起的地面沉降問題,結合明德林解答和二維卷積公式,得到地面沉降的簡明計算公式。借助Mathematica實現求解,驗證計算公式的正確性,并對承壓水降壓下地面沉降的影響因素進行分析。

1 承壓水三維流有限層法

圖1所示為典型的下臥承壓含水層,在其中利用抽水井開采地下水,在抽水井附近將出現地下水位降深,形成地下水的滲流運動。根據達西定律和質量守恒定律,抽水作用下引起的地下水三維非穩定流控制方程如下:

圖1 承壓含水層抽水井開采示意圖

(1)

式中:s(x,y,z,t)為降深;K為含水層滲透系數;q(x,y,z,t)為匯源項;Ss為貯水率。

式(1)中q(x,y,z,t)表示單位時間從單位體積含水層中抽取的地下水量,表達式取決于抽水井的具體形式,對于圖1所示的豎向非完整井,假如井中心位置為 (x0,y0,z0),則取如下形式:

q(x,y,z,t)=Q(z)δ(x-x0)δ(y-y0),

(2)

式中:

有限層法是一種介于解析法和數值法之間的混合算法。區別于計算域的全局離散,有限層法沿z方向將含水層離散成L個層元,水位降深在x、y平面采用正交完備的解析函數系{Amn(x,y)}表示,結合z方向的有限元形函數Nj(z),得到級數形式的降深試探函數[12]:

(3)

根據標準的有限元格式,得到關于未知系數的線性方程組:

(4)

式中:Φmn={φmnj}為未知系數向量;滲透性矩陣[G]mn、貯水矩陣[B]mn和流量矢量Qmn可以由具有解析形式的公式計算得到,具體形式參見文獻[14-15]。

在傳統的數值方法中,數值積分是進行單元分析的必要環節,占到相當大的計算工作量。相比來說,有限層法不涉及數值積分,在地下水三維非穩定流分析時,計算工作量更低。同時,由于試探函數中解析函數的正交性,求解格式具有天然的解耦性,通過并行計算可以進一步提高計算效率[16]。

2 承壓水開采引起的地面沉降

2.1 水位降深引起的附加應力

在抽水作用下,承壓含水層中形成的降落漏斗將打破原有的平衡狀態,導致上覆土層浮托力的減少,上覆土層底部(承壓含水層頂板)形成豎向的附加應力, 這種由水位下降導致的豎向附加分布力為

f(x,y,t)=γw·s(x,y,0,t),

(5)

式中,γw是水的單位重度。

2.2 地面沉降計算公式

承壓含水層頂板為隔水層,不用考慮含水層水位下降引起的潛水層越流和上覆土層的固結,地面沉降主要由含水層頂板處的附加作用力(5)引起。根據彈性力學的明德林解答可知,無限半空間體內,單位豎向集中力作用下,地表的豎向位移響應為

(6)

式中:E和ν分別為彈性模量和泊松比;h是上覆土層的厚度。

基于二維卷積公式,結合單位集中力位移響應(6)和級數降深函數(3),可以得到分布力(5)作用下的地表沉降計算公式

w(x,y,t)=

w′(x-η,y-θ)dηdθ。

(7)

利用地下水流的有限層求解格式(4),將得到的降深級數計算解代入式(7)進行積分,可以分析抽水作用下的地表沉降過程。如果需要對土層內部的位移、應力等情況進行分析,可以用相應物理量的明德林解在式(7)中對地表沉降w′進行替換。式(7)難以進一步展開成解析形式,本文利用Mathematica進行求解。

3 對比驗證與影響因素分析

3.1 對比驗證

為了驗證該方法的正確性,選取典型算例進行了分析,其中下臥含水層中抽水井為完整井,地下水流為穩定流。圖2給出了不同彈性模量取值(E=1 MPa, 2 MPa, 5 MPa)下,地表中心沉降w0隨覆蓋層厚度變化的計算結果,并將其與已有解進行了對比?？梢钥闯?三種情況下本文給出的計算結果都與已有解[12]吻合良好,驗證了本文方法的正確性和合理性。同時,圖2計算結果也反映了上覆土層厚度h對地面沉降的影響,在特定彈性模量取值下,隨著h的增加,含水層埋深越大,地下水開采的影響越小,所引起的地表沉降越小。當土層厚度h=30 m,泊松比ν=0.4,E=1 MPa時,在承壓水降深s=15 m的條件下,地面中心沉降接近0.015 m,這與文獻[12]中江南工作井的計算和實測結果相吻合。

圖2 地表中心沉降隨覆蓋層厚度的變化

地下水開采是地面沉降發生的誘因,本文研究了開采量對地面沉降的影響,見圖3?？梢钥闯?開采量對地面沉降具有直接影響,隨著開采量的增加(Qw從500 m3/d增加至1000 m3/d),地面沉降明顯加劇,中心沉降從0.013 m增大到接近0.27 m,兩者呈正相關,保持良好的線性關系。因此,結合當地水文地質條件,特別是開采含水層的補給條件,制定合理地下水開采量是地面沉降災害預防的關鍵。

圖3 抽水量對地面沉降的影響

文獻[12]的地面沉降計算公式可以分析中心沉降,不能給出地面沉降的空間分布。利用本文方法給出了抽水井附近的地面沉降分布,如圖4所示?？梢钥闯?豎向單井抽水作用引起的地面沉降呈漏斗狀分布,與水位降深分布保持一致,最大沉降出現在抽水井中心,隨距井徑向距離r的增加逐漸減小。同時,由于本文方法中考慮了地下水的非穩定流,可以給出地面沉降的發展過程,圖4給出了不同時刻(t=0.5, 1, …, 10, 20 d)的沉降曲線,可以看出,隨著抽水時間的延續,井附近的地面沉降不斷加劇。

圖4 地面沉降隨抽水時間變化

3.2 影響因素分析

地面沉降涉及地下水運動和土體變形的相互作用,影響因素很多。本文計算公式建立過程中,雖然采用了一定簡化(忽略上覆土層的固結變形),卻考慮了一些重要水文地質參數(如滲透系數K和貯水率Ss)與土力學參數(如彈性模量E和泊松比ν)的影響。

圖5給出了滲透系數K對沉降的影響,可以看出,K數值越大,含水層內水力聯系越強,井附近的水位降深越均勻,中心沉降和差異沉降越小。圖6表明含水層貯水率Ss對沉降也有顯著影響,隨著Ss的增加,含水層中的地下水儲量越大,抽水作用下含水層內的水位降深越小,所引起的沉降值也越小。

圖5 滲透系數對地面沉降的影響

圖6 貯水率對地面沉降的影響

圖7和圖8給出了上覆土層力學參數的影響,可以看出,其他條件保持不變的情況下,覆蓋土層剛度越大,土體變形越小,沉降量越小,當彈性模量E從1增加到10 MPa時,最大沉降量下降了87.5%;相對而言,泊松比ν對沉降的影響較小,當ν在從0.15增加到0.55時,最大沉降值減少了25%。

圖7 彈性模量對地面沉降的影響

圖8 泊松比對地面沉降的影響

3.3 多種抽水形式下的地面沉降

近年來,水平井憑借比常規豎向井更高的抽水效率,在環境和水文方面得到越來越多的應用[17-19]。與普通的井網結構不同,水平井的井過濾器水平放置(圖9)。在地下水運動特征上,一般來說,水平井引起的水流模式比傳統豎向井更加復雜,具有明顯的三維流,并缺乏軸對稱性[20]。

圖9 含水層中水平抽水井示意圖

在有限層求解格式中,水平井過濾器可以完全嵌入剖分節面,從而只需對求解格式(4)中{Qmn}進行適當修改,無需離散就可以實現求解[15]?；诒疚姆椒?水平井附近的地表沉降曲線如圖10,為了便于對比,圖中也給出了相同抽水量下豎向井的地表沉降分布?？梢钥闯?相比于豎向井,水平井開采引起的地表沉降不存在軸對稱性,沿過濾器方向的沉降和垂直方向的沉降并不一致,沿井軸方向沉降值更大。此外,總體上,水平井引起的沉降分布比豎向井更為平緩,這與水平井的水力學特征相吻合[15, 20]。

圖10 水平井和豎向井引起的地面沉降分布曲線對比

為了優化開采方案,在含水層中通常會設置群井系統。在傳統有限元數值方法中,為了獲得精確的近井水流模擬結果,每一個抽水井附近,都需要更精細的網格劃分,這會帶來額外的工作量,甚至造成網格不規則和離散誤差。有限層法中,由于采用具有解析表達式的流量矢量Qmn代替網格劃分,不存在上述問題,因此在處理多井問題時更加方便。

為了說明本文方法處理多井問題的能力,設計了由三口豎向井和一口水平井組成的群井系統,具體計算參數見表1。圖11給出了群井開采下地面沉降的等值線圖,可以看出,在抽水井附近形成了四個局部沉降錐。三口豎向井附近的等高線近似于同心圓,水平井附近的等高線為橢圓,正確反映了豎向井和水平井的水力特性,說明了本文方法應用于群井問題的適用性。

表1 群井問題計算參數

圖11 多井抽水地面沉降等值線

4 結 論

針對開采承壓水誘發的地面沉降問題,基于地下水三維非穩定流的有限層方法,利用二維卷積公式和彈性半空間明德林解,提出了一種簡明的地面沉降計算方法。通過與已有解的對比,驗證了本文方法的正確性。本文簡化計算公式可以反映地面沉降的空間分布及其時間演化過程。通過算例分析,研究了多種因素對抽水地面沉降的影響規律,結果表明:(1)上覆土層厚度越大,地面沉降量越大;彈性模量對地面沉降影響較大,覆蓋土層剛度越大,地面沉降越小。相對而言,泊松比的影響較小;(2) 含水層滲透系數和貯水率對地面沉降也有一定影響,滲透系數越大,含水層內水力聯系越強,差異沉降越小;含水層的儲水量越大,地面沉降越小;(3) 對水平井、群井問題進行了分析,所得結果符合相應的水力學特征,說明該方法在處理復雜抽水問題時也具有一定的適用性。