潛水一維非穩(wěn)定流方程的數(shù)值求解與工程應(yīng)用

李 茹，魯海峰

(安徽理工大學(xué) 地球與環(huán)境學(xué)院，安徽 淮南 232001)

0 引言

研究潛水非穩(wěn)定運動規(guī)律對于工程建設(shè)、地下水動態(tài)預(yù)報都有著重要的意義。潛水一維非穩(wěn)定流方程在基坑降水水位、區(qū)域地下水流水位、溝渠流水位等預(yù)測預(yù)報方面有著廣泛應(yīng)用。目前眾多學(xué)者針對潛水一維非穩(wěn)定流方程的求解問題開展了大量的研究，如任紅蕾等[1]建立了含階梯函數(shù)型源匯項的溝渠附近潛水非穩(wěn)定滲流模型，應(yīng)用Laplace變換性質(zhì)，求出了模型解析解，并利用數(shù)值法驗證了解析解的可靠性；秦志忠[2]建立基坑降水有限元模型，探究了基坑周邊水位降深變化曲線；王照亮、郭昆等[3-4]通過建立數(shù)學(xué)模型以及水文地質(zhì)概念模型，利用FEFLOW數(shù)值模擬軟件研究了地下水流場、分層的均衡性、潛水埋深、潛水水位；徐進(jìn)、劉能勝等[5-6]基于Neuman模型，利用伽遼金法與正交解析函數(shù)族推導(dǎo)了解耦形式的加權(quán)余量方程式，并在編制Fortran計算程序?qū)崿F(xiàn)數(shù)值求解的基礎(chǔ)上，提出了潛水非穩(wěn)定流的半解析數(shù)值求解格式；張年學(xué)等[7]提出了一種新的解析方法，對隔水層正傾和反傾兩種情況，選取適當(dāng)坐標(biāo)，將含水層分為上下兩部分，根據(jù)層間越流流量變化率相等原理與質(zhì)量守恒定律，導(dǎo)出全微分方程，獲得了解析解；夏強等[8]針對河渠間潛水的一維非穩(wěn)定運動，考慮了補給強度的時變性，根據(jù)Boussinesq方程的第一線性化方法，應(yīng)用Duhamel原理得到了方程的解析解。

上述研究成果在一定程度上對于潛水一維非穩(wěn)定流方程的求解及應(yīng)用起到了積極的指導(dǎo)作用，但聯(lián)合采用解析解和數(shù)值解系統(tǒng)探討此類水流方程的求解精度以便更好地指導(dǎo)實際工程的研究成果則較少。為此，本文采用分離變量法與有限差分法，以一類邊界條件為例，系統(tǒng)探討潛水一維非穩(wěn)定流方程的求解方法與精度變化規(guī)律，擬為工程實踐提供理論指導(dǎo)。

1 潛水一維非穩(wěn)定流方程模型及求解

1.1 方程建立

在描述河間區(qū)域地下水流、排水溝渠等地下水頭隨時間變化的非穩(wěn)定滲流場時，一般可建立潛水一維非穩(wěn)定流方程，在初始條件下地下水頭有一個二次函數(shù)形式的降深時，潛水一維非穩(wěn)定流方程模型初始降深曲線如圖1所示。

圖1 一維非穩(wěn)定流方程模型初始降深曲線

在一類邊界條件[9]下，潛水一維非穩(wěn)定流的偏微分方程及其定解條件為：

(1)

式中，T為含水層導(dǎo)水系數(shù)(m2/d)，μ為潛水含水層給水度，S為水位降深(m)，x為含水層寬度方向坐標(biāo)(m)，t為時間(d)，L為研究潛水含水層水頭降深變化的范圍長度。

1.2 方程求解

偏微分方程的求解方法一般有解析法和數(shù)值法等[10]，解析法通常采用分離變量法，數(shù)值法通常采用有限差分法。本文同時采用這兩種方法求解，對數(shù)值解的精度進(jìn)行系統(tǒng)討論和分析，并將顯式差分法和隱式差分法不同網(wǎng)格密度的數(shù)值解與解析解進(jìn)行對比。

1.2.1 解析法

由式(1)可知，該方程為二階偏微分方程，邊界條件為一類邊界條件，且此邊界條件為齊次[11]。所以這里采用分離變量法來求解，具體求解過程如下。

將二階齊次偏微分方程轉(zhuǎn)化為兩個常微分方程。假設(shè)給定方程的解具有以下形式：

S(x，t)=X(x)T(t)，

(2)

(3)

由式(1)中的邊界條件可得：

(4)

將方程(3)分為兩個部分進(jìn)行求解，結(jié)合給出的邊界條件，根據(jù)λ取值的正負(fù)，最后推出給定問題的特征值是：

(5)

對應(yīng)的特征函數(shù)分別具有以下形式：

(6)

這里Cn是任意常數(shù)。

由式(6)得到的方程通解及式(1)中的初始邊界條件S(x，0)=-0.01x2+x，可得：

(7)

由此解析解方程可表示為：

(8)

取L=100 m，T=51 m2/d，μ=0.3，根據(jù)式(8)，運用Matlab軟件求解出降深S隨時間t變化的規(guī)律曲線，如圖2所示。

圖2 解析法求出的降深隨時間變化的規(guī)律曲線

由圖2可以看出，在一定時間內(nèi)，降深S在含水層區(qū)間范圍內(nèi)呈規(guī)律性變化，降深在邊界處的水位不變，在區(qū)間中部最大。由于含水層區(qū)間兩側(cè)有水源的徑流，可補給該潛水含水層，因此降深隨時間的增大而減小，經(jīng)過第10 d后，最大降深已經(jīng)小于5 m，水位恢復(fù)較快。

1.2.2 數(shù)值法

解析法一般只能針對方程以及邊界較為簡單的情形求解，隨著計算機技術(shù)的飛速發(fā)展，數(shù)值法已成為偏微分方程的主要求解方法。本文分別采用顯式差分和隱式差分兩種方法進(jìn)行數(shù)值法求解[12]，并將求解結(jié)果與解析解進(jìn)行對比，進(jìn)而討論求解精度。

1.2.2.1 顯式差分法

將式(1)中的x方向分成n個節(jié)點，共有(n-1)個間隔；將t分成m個節(jié)點，共有(m-1)個間隔；將偏微分方程寫成差分方程，x方向步長為h，t方向步長為k?？傻玫饺缦路匠蹋?/p>

(9)

由式(1)和方程(9)聯(lián)立可得：

(10)

進(jìn)一步對式(10)中的步長h和k進(jìn)行變換，得：

(11)

Si，j+1-Si，j=γ(Si+1，j+Si-1，j-2Si，j)。

(12)

將式(12)變換，得：

Si，j+1=γSi+1，j+(1-2γ)Si，j+γSi-1，j。

(13)

式(13)為該模型的顯式差分方程式。運用Matlab軟件求解出降深S隨時間t變化的規(guī)律曲線，如圖3所示。

圖3 顯式差分法求出的降深隨時間變化的規(guī)律曲線

由于顯式差分法剖分的網(wǎng)格數(shù)會對計算結(jié)果產(chǎn)生較大的影響[13]，為探討網(wǎng)格數(shù)對精度的影響規(guī)律，本文分別設(shè)置三種不同密度的網(wǎng)格。低密度網(wǎng)格時間t的步長為0.2 d，x的步長為10 m，結(jié)果為圖3中的虛線部分，從中可以看出，在時間開始時刻即t=0 d時，函數(shù)圖像曲線較為曲折，與解析法模擬出的結(jié)果相差較大；在t=6 d時函數(shù)圖像曲線較為光滑，可近似擬合解析法計算結(jié)果。中密度網(wǎng)格時間t的步長為0.02 d，x的步長為1 m，中密度網(wǎng)格顯式差分圖像如圖4所示。由圖4可以看出，當(dāng)方程中的參數(shù)γ≥0.5時，計算結(jié)果不收斂，無法給出方程的解，這是由顯式差分法的固有缺陷造成的[14]。高密度網(wǎng)格時間t的步長為0.002 d，x的步長為1 m，結(jié)果為圖3中的實線部分，從中可以看出，函數(shù)圖像曲線整體光滑，在x=0處的同一時刻，通過此方法模擬出的結(jié)果與解析法模擬出的結(jié)果基本吻合；在接近區(qū)間中點即x=50 m時，得到的結(jié)果準(zhǔn)確度較低。

圖4 中密度網(wǎng)格下顯式差分法計算結(jié)果

1.2.2.2 隱式差分法

由于顯式差分法對網(wǎng)格剖分有著較高的要求，下面討論隱式差分法的求解過程。將式(1)寫為隱式差分格式，可得：

(14)

由式(1)和方程(14)聯(lián)立可得：

(15)

Si，j+1-Si，j=γ(Si+1，j+1+Si-1，j+1-2Si，j+1)。

(16)

將式(16)變換后為：

(1+2γ)Si，j+1-γSi+1，j+1-γSi-1，j+1=Si，j，

(17)

式中，i=2，3，…，n-1，共(n-2)個節(jié)點。由此可建立方程組為：

(18)

根據(jù)此方程可變換為以下矩陣形式：

(19)

根據(jù)式(19)，采用追趕法在Matlab軟件中進(jìn)行求解，并給出降深S隨時間t變化的規(guī)律曲線，結(jié)果如圖5所示。

圖5 隱式差分法求出的降深隨時間變化的規(guī)律曲線

為探討網(wǎng)格數(shù)對精度的影響規(guī)律，同顯式差分法一樣采用三種不同密度的網(wǎng)格。低密度網(wǎng)格時間t的步長為0.2 d，x的步長為10 m，結(jié)果為圖5中的虛線部分，從中可以看出，在網(wǎng)格密度過低時，計算結(jié)果存在較大的偏差。高密度網(wǎng)格時間t的步長為0.002 d，x的步長為1 m，計算結(jié)果為圖5中的實線部分，從中可以看出，隨著網(wǎng)格數(shù)的增加，曲線更為光滑，計算結(jié)果更為精確，當(dāng)t=0 d時，所呈現(xiàn)的曲線為初始降深，初始最大降深為25 m，隨著時間的增加，在地層同一點處的降深不斷減小。

為了與顯式差分法進(jìn)行對比，同樣采用顯式差分法中的中密度網(wǎng)格，時間t的步長為0.02 d，x的步長為1 m，結(jié)果如圖6所示。對比圖6和圖4可以看出，隱式差分的網(wǎng)格剖分不受γ取值的限制，在不滿足γ<0.5時，計算結(jié)果依然收斂，由此顯示出隱式差分法的優(yōu)越性。

圖6 中密度網(wǎng)格下隱式差分法計算結(jié)果

1.3 結(jié)果對照

為反映各種計算方法的精確度，選取其中一組數(shù)據(jù)進(jìn)行誤差對比。這里選取時間為6 d時的數(shù)據(jù)，對照結(jié)果如表1所示。

表1 部分節(jié)點降深取值對照表 m

2 實例應(yīng)用

下面對潛水一維非穩(wěn)定流方程在一類邊界下的各種求解方法進(jìn)行實例運用。某地鐵線路在建設(shè)施工過程中，在基坑開挖時遇到了潛水含水層，這將會對工程建設(shè)以及基坑開挖產(chǎn)生重大影響[15]。圖7為某地鐵線路建設(shè)過程中基坑開挖示意圖?；娱_挖寬度為20 m，深度為15 m，長度為100 m，可視為潛水一維非穩(wěn)定滲流。根據(jù)巖土勘察結(jié)果，場地潛水含水層水位埋深為3 m，含水層厚度為17 m，巖性主要為粉砂，含水層滲透系數(shù)K=5.1 m/d。工程設(shè)計時采用井點抽水方式對基坑進(jìn)行降水，為準(zhǔn)確掌握抽水停止后水位恢復(fù)情況，采用潛水一維非穩(wěn)定流方程進(jìn)行求解。

2.1 建立模型

假設(shè)抽水一定時間后，地下水位的降深曲線滿足二次函數(shù)形式(如圖7所示)，以下將討論抽水停止后基坑水位恢復(fù)情況，其水位降深方程為：

(20)

其定解條件可描述為：

(21)

含水層滲透系數(shù)K=5.1 m/d，含水層厚度M=20 m，給水度μ=0.3，將潛水含水層線性化[16]，可得到含水層導(dǎo)水系數(shù)T=K(M/2)=51 m2/d。

圖7 某地鐵線路基坑開挖示意圖

2.2 求解方程

根據(jù)前面的研究結(jié)果，采用解析法和隱式差分法進(jìn)行求解，圖8為隱式差分法和解析法模擬抽水停止后水位恢復(fù)的動態(tài)演變過程。

圖8 某地鐵線路基坑降深恢復(fù)曲線

通過圖8可知，地鐵工程的不同時間段對地下水位的影響程度不同。停止抽水時，該基坑地下含水層中部即x=50 m處的水頭降深最大，為25 m。隨著時間的變化，地下水位慢慢恢復(fù)，在停止抽水后第4 d水位恢復(fù)速率減慢，到第10 d水位基本恢復(fù)，直至最后地下水降落漏斗消失，水位恢復(fù)到初始狀態(tài)，且地下水環(huán)境趨于平穩(wěn)。根據(jù)結(jié)果分析可知，基坑停止抽水1 d后基坑中部水位已上升到坑底。

3 結(jié)論

(1)對潛水一維非穩(wěn)定流水位降深規(guī)律進(jìn)行研究，建立降深隨時間變化的方程模型，借助分離變量法和有限差分法求解方程的解析解和數(shù)值解，對數(shù)值解剖分不同密度的網(wǎng)格，采用顯式差分法和隱式差分法進(jìn)行求解，以此討論計算精度變化規(guī)律。

(2)顯式差分和隱式差分求解方法的模擬結(jié)果相差不大，在相同密度網(wǎng)格剖分下，顯式差分法較隱式差分法的求解精度要高，但顯式差分法存在網(wǎng)格約束條件，在實際應(yīng)用時需注意網(wǎng)格剖分密度。

(3)實例應(yīng)用表明，潛水一維非穩(wěn)定流方程可很好地用于實際工程，為工程建設(shè)過程中地下水滲流計算提供理論基礎(chǔ)。