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

徐 進(jìn),張 豪,王斯海

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

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

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

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

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

1 承壓水三維流有限層法

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

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

(1)

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

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

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

(2)

式中:

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

(3)

根據(jù)標(biāo)準(zhǔn)的有限元格式,得到關(guān)于未知系數(shù)的線性方程組:

(4)

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

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

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

2.1 水位降深引起的附加應(yīng)力

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

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

(5)

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

2.2 地面沉降計(jì)算公式

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

(6)

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

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

w(x,y,t)=

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

(7)

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

3 對比驗(yàn)證與影響因素分析

3.1 對比驗(yàn)證

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

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

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

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

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

圖4 地面沉降隨抽水時(shí)間變化

3.2 影響因素分析

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

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

圖5 滲透系數(shù)對地面沉降的影響

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

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

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

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

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

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

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

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

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

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

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

表1 群井問題計(jì)算參數(shù)

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

4 結(jié) 論

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