第1-3章 基本原理

水文地質(zhì)數(shù)值計算郭巧娜地球科學(xué)與工程學(xué)院第一章地下水流動定解問題概述1、三維流微分方程2、數(shù)學(xué)模型3、微分方程4、邊界條件5、模型概化6、實例

取右圖所示的微小六面體。設(shè)與x,y,z,方向?qū)?yīng)的主滲透系數(shù)分別為Kxx,Kyy,Kzz；建立均衡期t時段內(nèi)，微小均衡六面體的水量守恒方程。1、三維流微分方程1、三維流微分方程x，y，z方向流入—流出分別為:t時段內(nèi)，六面體水量變化量為：六面體內(nèi)地下水儲存量的變化為由水均衡原理得1、三維流微分方程方程兩端除以Δt，并取Δｘ→０,Δｙ→０,Δｚ→０和Δｔ→０，則一般密度的空間變化率很小，故于是有由達西定律有（2）水流連續(xù)性方程左端項都很小，可以忽略。1、三維流微分方程上式為非均質(zhì)各向異性承壓含水層的偏微分方程。均質(zhì)各向異性非穩(wěn)定流均質(zhì)各向異性穩(wěn)定流得到地下水三維流動微分方程[1/L]1、三維流微分方程微分方程定解條件邊界條件初始條件已知t=0時的因變量，H(x,y,z,0)=H0(x,y,z)已知水頭邊界(I類邊界)H(x,y,z,t)=f(x,y,z,t)(x,y,z)B1特例：定水頭邊界H(x,y,z,t)=C已知流量邊界特例：隔水邊界2、數(shù)學(xué)模型數(shù)學(xué)模型地下水運動的數(shù)學(xué)模型結(jié)構(gòu)(1)第一類邊界條件(簡記為DirichletBC)第一類邊界條件也稱為給定水頭邊界條件，指的是滲流區(qū)某一部分邊界上，各點水頭在某一時刻是已知的，表示為：

(2)第二類邊界條件(簡記為NeumannBC)第二類邊界條件也稱為定流量邊界條件，是指某一部分邊界單位面積上流入(流出時值為負)的流量已知。相應(yīng)的邊界條件表示為：注水井或抽水井也可以作為內(nèi)邊界來處理，此時的承壓含水層可忽略其水頭的垂向變化，注水井或抽水井開口于整個厚度，如圖所示，設(shè)井的半徑為

，單位時間抽水量為

，井壁法向

和

軸方向相反，如圖中所示。承壓抽水井示意圖則由Darcy定律，有：

（3）

第三類邊界條件(RolinBC)第三類邊界條件也稱為混合邊界條件，是指知道某段邊界上和的線性組合：臨海含水層與海水之間有一層薄淤泥層，則淤泥層兩側(cè)同一位置上的點存在水頭差，內(nèi)側(cè)含水層水頭

、滲透系數(shù)分別為

，淤泥層厚度為

，水頭、滲透系數(shù)分別為

和

，外側(cè)海水水頭為

，忽略淤泥層內(nèi)孔隙水的彈性貯存的變化，在淤泥層與含水層交界面上根據(jù)流量的相等關(guān)系，有：淤泥層與含水層交界面上根據(jù)流量的相等關(guān)系，其中

（4）等值面邊界條件

等值面邊界條件適用于不完全井的情況，即井壁和底部都有含水層中的水流滲入。等值面邊界條件的表達形式是：已知，并有正、負之分，分別對應(yīng)于單位時間內(nèi)的注入水量和抽出水量。是井壁或渠壁與含水層的所有接觸面。（5）

自由面邊界條件自由面邊界條件即潛水面邊界條件，若潛水含水層中水頭為

，自由面方程：

，那么，在自由面上，有：（5）

自由面邊界條件3、微分方程

數(shù)學(xué)模型4、邊界條件

邊界條件：滲流區(qū)邊界上水力特征，即邊界上的水頭分布和變化特征或流入流出含水層的水量分布和變化情況。主要有兩類：1．已知邊界上的水頭分布規(guī)律

HB1=ψ(x,y,t),其中ψ(x,y,t)為已知函數(shù)。主要常見的是滲流區(qū)與地表水體相接觸。2．已知邊界上的單位寬度流量q隨時間的變化規(guī)律如已知流量為Q的承壓含水層中完整的抽水井，其井壁可以看作此類邊界。

數(shù)學(xué)模型5例子:河間地塊承壓水流模型

設(shè)兩條河流平行、完全切割乘壓含水層，含水層等厚、均質(zhì)各向同性,無垂向補給或排泄，對于如圖所示的坐標系，已知某時刻的含水層各處的水頭為20米，自該時刻后，河水位分別為如圖所示的函數(shù)。試根據(jù)條件作合理簡化建立其數(shù)學(xué)模型。（1）模型概化

由所述水文地質(zhì)條件，可以概化為一維承壓水流問題。（2）建立坐標系（如圖）

取x-軸原點位于左端河，右側(cè)為正向，設(shè)兩河流間距為L.

縱軸為水頭。（3）數(shù)學(xué)模型河間地塊承壓水流模型（續(xù)）

1有限差分法原理2導(dǎo)數(shù)的有限差分近似表示3承壓一維流動有限差分法4承壓二維不穩(wěn)定流有限差分法5源匯項的處理及井孔水頭校正6不規(guī)則邊界問題7矩形變格距網(wǎng)格差分第二章有限差分法1有限差分法的基本原理將連續(xù)的問題離散后求解：方法一．以地下水流基本微分方程及其定解條件為基礎(chǔ)，在滲流區(qū)剖分基礎(chǔ)上，用差商代替微商，將地下水流微分方程的求解轉(zhuǎn)化為差分方程（代數(shù)方程）求解。方法二．在滲流區(qū)剖分的基礎(chǔ)上，直接由達西定律和水均衡原理，建立各個均衡區(qū)的水均衡方程，即差分方程。矩形網(wǎng)格多邊形網(wǎng)格1有限差分法的基本原理

——網(wǎng)格劃分的基本類型（1）先劃格線，格點位于網(wǎng)格中心均衡網(wǎng)格節(jié)點網(wǎng)格（2）先規(guī)定格點位置，再垂直平分兩相鄰結(jié)點的連線作格線，形成的網(wǎng)格即為水均衡區(qū)基本概念將時間離散點和空間離散點聯(lián)合組成的網(wǎng)格稱為時空網(wǎng)格。有限差分的基本原理：某點處水頭函數(shù)的導(dǎo)數(shù)用該點和幾個相鄰點處水頭值及其間距近似表示。MODFLOW網(wǎng)格系統(tǒng)導(dǎo)數(shù)的有限差商近似導(dǎo)數(shù)的定義

當非常小的時候，有

上式右端項即為f(x)在x0處的差商。

這樣定義的差商很容易理解，但不知道用差商代替微商所產(chǎn)生的誤差。下面利用泰勒公式導(dǎo)出差商及其誤差。方法一：差商代替微商2導(dǎo)數(shù)的有限差分近似表示已知泰勒公式①由A得：

AB②由B得：

稱為f(x)在x0處的一階前向差商，為截斷誤差。稱為f(x)在x0處的一階后向差商，為截斷誤差。方法一

③由A-B可以得：

④由A+B可以得：AB稱為f(x)在x0處的一階中心差商，為截斷誤差。稱為f(x)在x0處的二階中心差商，為截斷誤差。方法一對于偏導(dǎo)數(shù)（偏微商），類似可以得到相應(yīng)的差商：方法一

取右圖所示得微小六面體。設(shè)與x,y,z,方向?qū)?yīng)得主滲透系數(shù)分別為Kx,Ky,Kz；建立均衡期t時段內(nèi)，微小均衡六面體的水量守恒方程。方法二：達西定律和水均衡原理基于達西定律，x，y，z方向流入—流出分別為:t時段內(nèi)，側(cè)向流入與源匯項導(dǎo)致六面體水量變化量為：ABC方法二：達西定律和水均衡原理DA+B+C+D源匯項六面體內(nèi)地下水儲存量的變化為由水均衡原理得三維地下水流動方程的有限差分格式方法二：達西定律和水均衡原理有限差分法:三維（MODFLOW）差商代替微商3承壓一維流動有限差分法顯式差分格式隱式差分格式方法一控制方程網(wǎng)格剖分nx個算例：顯式有限差格式

設(shè)兩條河流平行、完全切割含水層，含水層等厚、均質(zhì)各向同性。應(yīng)用實例：河間地塊承壓水流模型3承壓一維流動有限差分法步驟：（1）基礎(chǔ)資料的分析（2）概念模型（3）數(shù)學(xué)模型（4）數(shù)值方法及計算機程序（5）參數(shù)（6）結(jié)果分析3承壓一維流動有限差分法建立數(shù)學(xué)模型（1）模型概化由所述水文地質(zhì)條件，可以概化為一維承壓水流問題。（2）建立坐標系（如圖），將地下水流動系統(tǒng)空間結(jié)構(gòu)放在坐標系內(nèi)，從而量化各變量的取值范圍。本例，取x-軸原點位于左端河，右側(cè)為正向，設(shè)兩河流間距為L.（3）數(shù)學(xué)模型3承壓一維流動有限差分法差分方程及其解法—顯式格式①將（0—L）分成N等份，1)網(wǎng)格剖分：②取時間步長，記（n=0、1、2、3、4……）記，(i=0,1,2,3，4……N)2）建立差分方程：在網(wǎng)格系統(tǒng)中任意取一點設(shè)是問題的解，則在處有記為（i，n）3承壓一維流動有限差分法用差商代替微商：

將上述兩式舍去余項，代入方程并記為顯然該式具有截斷誤差得到顯式格式（續(xù)1）引入無量綱變量：將該式子代入得到：

（i=1，2，3，.....N-1），（n=1，2，3，.....）

顯式格式（續(xù)2）3）顯示差分方程的求解計算各結(jié)點初始時刻水頭值利用差分方程計算各結(jié)點t1時刻水頭值利用邊界條件計算邊界結(jié)點水頭值重復(fù)2、3步，直到計算出擬計算的各個時刻的水頭值顯式格式（續(xù)3）算例（續(xù)4）在上述模型中，設(shè)L=1000米取空間步長為200米，時間步長為0.25天，分別計算各節(jié)點各時刻的水頭值。Time/dayx=0mx=200mx=400mx=600mx=800mx=1000m02010101010100.252012.5101010100.502013.7510.6251010100.752014.53111.25010.15610101.002015.07811.79710.39110.039101.252015.48812.26610.65410.11710算例（續(xù)5）Time/dayx=0mx=200mx=400mx=600mx=800mx=1000m02010101010100.252020101010100.502010201010100.75203002010101.0020-10101.252010算例（續(xù)6）如果⊿t=1,則算例：隱式格式在上述模型中，設(shè)L=1000米取空間步長為200米，時間步長為0.25天，用隱式差分格式計算各節(jié)點個時刻的水頭值。在這個例子中，解：隱式格式一般方程為于是有根據(jù)初始條件得根據(jù)邊界條件得由初始條件和邊界條件由此解得t1時刻的水頭值為在上述方程中取n=0,可以得到計算t1時刻水頭值的方程所以上述方程變成同理，可計算t2時刻的水頭值差分方程求解一維顯式差分格式網(wǎng)格個數(shù)為ni直接求解差分方程求解一維隱式差分格式網(wǎng)格個數(shù)為ni迭代求解方程組PCGSIPSORWHSSAMGGMGMODFLOW差分方程的收斂性和穩(wěn)定性截斷誤差：用差商代替微商時，地下水流動方程產(chǎn)生的誤差為截斷誤差。收斂性：當空間步長和時間步長趨于0時，有限差分方程的精確解趨于地下水流動問題微分方程定解問題的精確解。則稱該差分格式是收斂的。穩(wěn)定性：如果在求解差分方程過程中，某時間步引入某個誤差，而在以后的各時段計算中，該誤差不再擴大，則稱該差分格式是穩(wěn)定的。一維顯示格式的收斂條件和穩(wěn)定條件是：第二類邊界條件的處理前面介紹的幾種差分格式，都是以第一類邊界條件為例，說明其具體計算方法的。如果是第二類邊界條件，怎么使用這些格式進行計算呢？為此，我們以下面的定解問題來討論這個問題。設(shè)問題數(shù)學(xué)模型為用一階中心差分公式代替邊界條件中的偏導(dǎo)數(shù)4承壓二維不穩(wěn)定流有限差分法我們考慮無垂向補給、排泄、均質(zhì)、各向同性、等厚的承壓流動問題，邊界條件第一類，滲流區(qū)為矩形，則可用下述定解問題描述4承壓二維不穩(wěn)定流有限差分法（1）顯式格式--網(wǎng)格剖分（1）顯式格式2）建立差分方程：（1）顯式格式（1）顯式格式3）二維顯示差分方程的求解計算各結(jié)點初始時刻水頭值計算t1時刻水頭值計算邊界結(jié)點水頭值重復(fù)2、3步，直到計算出擬計算的各個時刻水頭值（1）顯式格式差分方程收斂性（2）隱式格式（2）隱式格式4）差分方程的收斂性和穩(wěn)定性條件5源匯項的處理及井孔水頭校正越流、入滲和抽水井等問題的處理如果考慮垂直滲流項（即源匯項），則二維承壓流動微分方程可寫成建立差分方程時，在結(jié)點處應(yīng)加上這一項，它可具體表示為：井水位校正對比圖

圖3-5（a）初始網(wǎng)格有限差分法計算井水位的校正圖3-5（b）加密網(wǎng)格6不規(guī)則邊界問題當研究區(qū)的幾何形狀屬于簡單形式（如矩形滲流區(qū)）時，差分網(wǎng)格的劃分往往將結(jié)點設(shè)在邊界上。然而，對于實際問題來說，邊界通常不是那么規(guī)則，邊界的某些部分，甚至大部分不能與結(jié)點重合，我們稱這種邊界是不規(guī)則的。關(guān)于不規(guī)則邊界問題，直接取最靠近邊界的曲折格線為近似邊界河流、水庫、湖泊等問題7矩形變格距網(wǎng)格差分本節(jié)直接由達西定律和水均衡原理建立差分方程達西定律：水均衡原理：對某一研究對象，流入-流出=體系內(nèi)質(zhì)量變化量研究對象可以是大區(qū)域，也可以是微分單元體。大區(qū)域的水均