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

取右圖所示的微小六面體。設與x,y,z,方向對應的主滲透系數分別為Kxx,Kyy,Kzz；建立均衡期t時段內，微小均衡六面體的水量守恒方程。1、三維流微分方程1、三維流微分方程x，y，z方向流入—流出分別為:t時段內，六面體水量變化量為：六面體內地下水儲存量的變化為由水均衡原理得1、三維流微分方程方程兩端除以Δt，并取Δｘ→０,Δｙ→０,Δｚ→０和Δｔ→０，則一般密度的空間變化率很小，故于是有由達西定律有（2）水流連續性方程左端項都很小，可以忽略。1、三維流微分方程上式為非均質各向異性承壓含水層的偏微分方程。均質各向異性非穩定流均質各向異性穩定流得到地下水三維流動微分方程[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、數學模型數學模型地下水運動的數學模型結構(1)第一類邊界條件(簡記為DirichletBC)第一類邊界條件也稱為給定水頭邊界條件，指的是滲流區某一部分邊界上，各點水頭在某一時刻是已知的，表示為：

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

，單位時間抽水量為

，井壁法向

和

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

（3）

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

、滲透系數分別為

，淤泥層厚度為

，水頭、滲透系數分別為

和

，外側海水水頭為

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

（4）等值面邊界條件

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

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

，自由面方程：

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

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

數學模型4、邊界條件

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

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

數學模型5例子:河間地塊承壓水流模型

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

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

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

縱軸為水頭。（3）數學模型河間地塊承壓水流模型（續）

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

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

當非常小的時候，有

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

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

AB②由B得：

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

③由A-B可以得：

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

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

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

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

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

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

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