第二章有限差分法第一節(jié)基本概念一、離散化二、地下水流的有限差分方程

1.模型的假設條件（1）地下水流為承壓水；（2）小均衡域除側向徑流外，無其它流量交換；（3）網(wǎng)格沿行及列分別為等距的，但行距與列距可以不等，且分別記為△X與△Y；（4）離散點上的近似水頭值以u表示。2.模型的建立沿右側流入：沿左側流入：沿下部流入：

沿上部流入：單位時間流入小均衡域的總側向徑流量Qt為：

假定Δt時間內，小均衡域的水頭抬高Δu，小均衡域增加的水量Qs為：據(jù)質量守恒有：Qt=Qs

非均值各向異性承壓含水層二維非穩(wěn)定流有限差分方程。3截斷誤差任意足夠光滑的函數(shù)f(x)沿x的正向和負向分別用Taylor級數(shù)展開，有

H(x,y)是描述地下水面的函數(shù)，地下水面是光滑的曲線，所以H(x,y)是光滑的函數(shù)，可以用泰勒級數(shù)展開。

因此，對于H(x+Δx)可用泰勒級數(shù)展開：舍去高階項，得：整理得：可見，可以用差商代替微商。

同樣，對于H(x-Δx)也可用泰勒級數(shù)展開：舍去高階項，得：整理得：第1式截斷誤差：兩邊除以△x，得：其中：第2式截斷誤差：兩邊除以△x，得：其中：用差商代替導數(shù)時，O(Δx)是被舍去的，所產生的誤差是O(Δx)，叫截斷誤差。

由于Δx是一次方，所以也叫一階誤差。由于沿x方向i+1在i的前面，因此將叫前向差商。

由于沿x方向i-1在i的后面，因此將叫后向差商。中心差商兩式減可導出：由于沿x方向i在i+1、i-1的中間叫中心差商。截斷誤差為二階，精度提高了一階。同樣，二階偏導數(shù)也可用差商來近似代替X方向Y方向同理，用差商代替水位對時間的導數(shù)用有限差分方程代替偏微分方程，將產生的誤差為：這就是為什么說數(shù)值法是一種近似方法。

4.有限差分方程的建立對于微分方程沿X軸有沿Y軸有對時間軸有將三項代入微分方程，得：以上是針對（i，j）節(jié)點建立的方程，為代數(shù)方程。同理，對于其它節(jié)點同樣能夠建立相似的代數(shù)方程。求解代數(shù)方程組就可以得到各節(jié)點的水位值。當Δx、Δy、Δt取得足夠小時，可望得到足夠精確的近似解。

比較此二式可見，均衡法與用差商代替微商所得結果相同。5.解的穩(wěn)定性和收斂性以u表示差分方程的準確解，u‘表示差分方程的數(shù)值解，H表示觀測值。

|H–u|是差分值的截斷誤差；

|u-u‘|是數(shù)值解的舍入誤差。如果|H–u|0時，差分方程的數(shù)值解滿足收斂性條件。如果|u-u‘|0時，差分方程的數(shù)值解滿足穩(wěn)定性條件。有限差分法模擬地下水流系統(tǒng)，就是求整個滲流區(qū)的u‘，且要求|H–u‘|小于預先給定的誤差范圍。

H–u‘=（H–u）+（u-u‘），即截斷誤差和舍入誤差的和。三、三種主要差分格式對于上方程右端項，沒有給出時間，可取m+1時刻，也可取m時刻。

當右端項取m+1時刻的水位時，為隱式差分格式：當右端項取m時刻的水位時，為顯式差分格式：

當右端項取m和m+1時刻的平均值時，為中心差分格式或對稱差分格式：以上三種差分格式可用下式表示：

θ=0，稱為顯示差分格式；

θ=1，稱為隱式差分格式；

θ=1/2，稱為中心差分格式或對稱差分格式。第二節(jié)承壓一維穩(wěn)定流一、水文地質模型的有限差分方程組

1.模型：如圖，在0與L處有兩條長而直且相互平行的河流切穿了承壓含水層。含水層均質，頂?shù)装逅?，厚度為M（m），滲透系數(shù)為K（m/d），頂板為弱透水層，降水入滲通過弱透水層向承壓含水層補給，補給強度為W（m3/m2·d）。兩條河流的水位分別為H0與Hl，且不隨時間變化。

數(shù)學模型為：

2.解析解：

3.有限差分法將0-L剖分成n份，空間步長為Δx=L/n，0與L分別位于0號小均衡域與n號小均衡域的中心，0節(jié)點的水頭為H0，L節(jié)點的水頭為HL。如圖，以i點為例，建立有限均衡方程。由i-1號小均衡域流進第i號小均衡域的水量為：由i號小均衡域流進第i+1號小均衡域的水量為：由于是穩(wěn)定流，流入的水量等于流出的水量：化簡后，對其它點，可建立同樣的方程，組成的方程組如下：方程組為線性，解該方程組可得到各節(jié)點的水位。該方程組的特點：（1）系數(shù)矩陣中的元素僅位于三條對角線上；（2）系數(shù)矩陣不僅對稱而且正定。這樣的方程有一種比較簡單的解法。二、解三對角型線性方程組的追趕法上方程可表示為寫成矩陣形式為：或系數(shù)矩陣用A表示：根據(jù)正定對稱性，A可被分解為兩個矩陣之積：或于是，方程組可寫成

方程組可寫成如下形式由此，可通過求的{y}。在通過求的{H}。求的{H}關鍵的問題變成了求[L]、[U]。

根據(jù)正定對稱性，A可被分解為兩個矩陣之積：矩陣對應元素必須相等。由此得：對于有代入上模型，設L=1000m，空間步長取100m，H0=10m，HL=20m，W=0.005m/d，T=100m2/d，n=10。有11個節(jié)點，0、10號節(jié)點水頭是已知的，1--9號節(jié)點未知水頭方程組為因為，W=0.005m/d，Δx=100m，T=100m2/d，所以，W（Δx）2/T=0.5ai=-1，bi=2，ci=-1將W=0.005m/d，Δx=100m，T=100m2/d，L=1000m，HL=20m，H0=10m

差分解解析解

上例n=10，請計算n=4，n=6，n=8時，節(jié)點處的水位值。同時計算相應的解析解，并進行比較。三、流量邊界的處理流量邊界是給定流量的邊界。特例為零流量邊界。零流量邊界有兩種：一種是隔水邊界；一種是水力坡度等于0。隔水邊界，邊界上的導水系數(shù)為零，這種情況不需要知道邊界上的水頭，處理的方法是將小均衡域的邊界落在含水層的邊界上，圖a，令Hn=Hn-1。

水力坡度等于0，這種情況需要知道邊界上的水頭，處理的方法是將小均衡域的節(jié)點落在含水層的邊界上，圖b，令Hn+1=Hn-1，以形成水力坡度等于0。一般情況下，不需要知道邊界上的水位，圖a即可。

對于側向流量不為零，可先按零流量邊界處理，然后將給定的流量再按垂向交換水量的形式加入邊界上相應小均衡域內，按源、匯處理。例，如圖模型，x=0處水位為H0，x=L處是單寬流量為q的排水溝。數(shù)學模型為：解析解為：積分，得當x=L時，，代入上方程，得解得C1：積分，得當x=0時，得C2=H0，將C1、C2代入方程，得滲流區(qū)域剖分如圖n號節(jié)點方程為：整理得：

將W=0.005m/d，Δx=100m，T=100m2/d，L=1050m，H0=10m，q=3m2/d代入方程組，得分解后得第三節(jié)潛水一維穩(wěn)定流一、潛水一維穩(wěn)定流水文地質模型如圖微分方程為：相應的數(shù)學模型為：解析解為：有限差分法：將L剖分成n份，空間步長Δx=L/n，h0與hL為0號與n號節(jié)點的潛水面高度值。對于i號小均衡域通過斷面1進入i單元的水量為通過斷面2進入i單元的水量為均衡方程為：化簡進一步化簡，得此為對i點建立的有限差分方程。在1~n-1的每一個節(jié)點上均可建立這樣的方程。形成一個代數(shù)方程組。這個代數(shù)方程組是非線性的。二、潛水有限差分方程組的迭代解法解非線性方程組的方法有兩種：（1）將非線性方程用線性方程組來近似，然后用迭代法逐漸逼近所求的解；（2）求極小值的方法。常采用線性化方法求解。該方程可化為：可利用迭代法求解上方程組，計算步驟如下：（1）先給（hi+hi-1）與（hi+hi+1）一個近似值（初值）與。這時將上方程線性化了。在1~n-1個節(jié)點上得到n-1階三對角型線性方程組。

（2）利用追趕法解此方程組，得1~n-1個節(jié)點上潛水位近似值，。（3）在用與代替（hi+hi-1）與（hi+hi+1），再次形成n-1階三對角型線性方程組。利用追趕法解此方程組，得1~n-1個節(jié)點上潛水位近似值，。（4）重復（3）的做法，直到m次。當時，終止計算。him就作為所要求的結果，i=1，2，…，n-1。

ε叫迭代收斂的誤差限，通常取10-5m。

設L=1000m，W=0.005m/d，H0=10m，HL=20m，K=100m/d，Δx=100m，代入上方程，并取H00=10m，

H10=11m，H20=12m，H30=13m，H40=14m，H50=15m，

H60=16m，H70=17m，H80=18m，H90=19m，H100=20m，得：以Δx=200m為例，設L=1000m，W=0.005m/d，H0=10m，HL=20m，K=100m/d，代入上方程，并取H00=10m，H10=12m，H20=14m，H30=16m，H40=18m，H50=20m，ε取10-3m。得：第四節(jié)承壓水一維非穩(wěn)定流一、承壓一維非穩(wěn)定流的有限差分方程組

1.模型及解析解假設：（1）兩河流平行，且切穿含水層；（2）含水層的頂、底板水平、均值各向同性、沿河流方向無限延伸；（3）兩河水位長期在同一高度，為H0（基準面為O—O'；（4）含水層中的水頭標高與河流水位一致為H0。（5）兩河水位于某時刻突然下降到O—O'。其數(shù)學模型為：解析解為：

2.有限差分法模擬將L分成n份，空間步長Δx=L/n，對i節(jié)點所在小均衡域建立差分方程。由水均衡原理得：由于是非穩(wěn)定流，地下水位是隨時間變化，設地下水位的變化率第i號小均衡域由于水頭變化引起地下水流的貯存量在單位時間的變化量為：

對于1斷面有兩個時階，tk時階流入單元的水量為，tk+1時階流入單元的水量為從tk時刻到tk+1時刻流入單元的總水量應取加權值：同理流出2斷面的水量為代入得：整理后得：對1—n-1號節(jié)點均可以建立如上的方程，這樣已知tk時刻的水位就可求得tk+1時刻的水位。二、三種有限差分格式的應用

1、θ=0時，上式變?yōu)椋夯虼藶轱@式差分格式。

2、θ=1時，上式變?yōu)椋捍藶殡[式差分格式。

3、θ=1/2時，上式變?yōu)椋捍藶閷ΨQ差分格式或中心隱式差分格式。令

4、顯式差分格式算例模型如上，T=1000m2/d，Δt=0.01d，L=1000m，Δx=100m，S=0.01，將參數(shù)代入下顯式差分方程得：結論：顯式差分格式，只有時，才是穩(wěn)定。由上式知，要使r減少只能是縮小Δt或增大Δx，縮短時間步長就要增加計算工作量，增大空間步長就要降低計算精度。顯式差分格式的優(yōu)點：不需要解聯(lián)立方程組。缺點：要保證穩(wěn)定，就會增加計算工作量，和降低計算精度。

5、隱式差分格式

因為：上式可變?yōu)椋喝=1，可得：對于5、6、7、8、9號節(jié)點可得方程組：6、中心差分格式

r=1，上式可變?yōu)椋簩τ?、6、7、8、9號節(jié)點可得方程組：結論：（1）當時，，其差分格式是穩(wěn)定的。（2）當時，r取任何值，差分格式都是穩(wěn)定的。（3）穩(wěn)定就收斂。三、有限差分方程的迭代解法常用的迭代法有：簡單迭代法、高斯-塞德爾迭代法、超松馳迭代法。

1、簡單迭代法（Jacobi法）以隱式差分格式為例。對于i點差分方程有：如果各項系數(shù)用A、B、C表示，常數(shù)項用F表示。則上方程可寫為：寫出各個節(jié)點上的方程，組成方程為：其中：簡單迭代法迭代終止準則：和

2、高斯-塞德爾迭代法簡單迭代法也叫同步迭代法。高斯-塞德爾迭代法也叫異步迭代法。

3、超松弛迭代法將高斯-塞德爾迭代法的方程寫成如下形式：進行上述迭代前的值為我們可以取上二值的加權平均值，作為第m+1次的結果：ω叫作松弛因子。一般取1—2，也叫作超松弛。方程也可寫為：這三種迭代法的特點：簡單迭代法收斂最慢；高斯-塞德爾迭代法收斂較快；超松弛迭代法，ω選1-2時，收斂最快。第五節(jié)潛水一維非穩(wěn)定流一、潛水一維非穩(wěn)定流的有限差分方程依據(jù)均衡原理，可推導出無垂向水量交換的差分方程，對i點有：因此，在所有1—n-1個節(jié)點上均可建立這樣的方程。形成一大型方程組。對于方程如果令可寫成這是一個三對角方程組，可用追趕法求解。其迭代追趕求解過程為：也可寫成迭代的形式：迭代求解過程為：對于方程如果令可寫成這是一個三對角方程組，可用追趕法求解。也叫雙重迭代法。其迭代追趕求解過程為：也可寫成迭代的形式：迭代求解過程為：第六節(jié)承壓二維非穩(wěn)定流一、差分方程下式為一維非穩(wěn)定流的差分方程：利用水量均衡原理，可得二維流的差分方程：沿x方向流入量和流出量之差為：沿y方向流入量和流出量之差為：貯存量為：對于i，j單元，差分方程為：整理得：二、三種差分格式

1、顯式差分格式

θ=0，二維差分格式的穩(wěn)定條件：方程組的求解方法：迭代法

2、隱式差分格式

θ=1，這是一個5對角型方程組。

3、對稱差分格式

θ=1/2，

產生的差分方程形式上同上。方程組的求解方法：迭代法、第七節(jié)不規(guī)則網(wǎng)格有限差分法

規(guī)則差分網(wǎng)格的特點：（1）建立有限差分方程時物理概念清楚，計算簡單編制計算程序容易。（2）有現(xiàn)成的軟件，如ProcessingModflow（美國環(huán)保署）和VisualModflow（美國地質調查所開發(fā)）。（3）擬合自然邊界及非均值界線時顯的不夠靈活。（4）在需要進行網(wǎng)格加密的地區(qū)進行網(wǎng)格加密時，會使那些不需要加密的地區(qū)也進行了加密，增加了計算工作量。不規(guī)則網(wǎng)格的特點：與規(guī)則網(wǎng)格相反。

一、地下水流區(qū)域的剖分不規(guī)則網(wǎng)格一般采用三角網(wǎng)格。三角形叫面元，三角形的邊叫線元，三角形的頂點叫點元或節(jié)（結）點，區(qū)域內部的點叫內點元，邊界上的點元叫邊界點元。剖分時注意：（1）點元必須是具有公共線元和面元的頂點。（2）面元的任一內角不應大于90°及小于30°。（3）剖分時點元的密度可根據(jù)在某些地區(qū)的需要靈活選擇，如在漏斗中心及地下水流場變化比較復雜地段要密。

（4）點元和面元均要編號。（5）點元編號時，不必分內節(jié)點與第二類邊界點，自小到大依次編號，編完之后再編第一類邊界點元的號。（6）編號時，第一類邊界與第二類邊界的交界點要取第一類邊界點元。（7）滲流區(qū)域內的各種界線及自然邊界均可用線元構成的折線逼近，如非均值界線、承壓水與潛水分界線等。（8）滲流區(qū)域內的抽水孔不一定要求與點元重合，但滲流區(qū)域內的觀測孔應盡量與點元重合。二、水頭模式的建立與面元水力坡度在一個滲流區(qū)，水頭函數(shù)是一個非常復雜的連續(xù)函數(shù)，用H（x，y，t）表示，如圖，難以找到合適的解析式來表示。解析解是坐標的函數(shù)，當要計算某點的水位時，只要將坐標代入方程就可以了。但是對于復雜的滲流區(qū)域是給不出解析解的。數(shù)值解的處理方法：將滲流區(qū)域剖分成n個小的單元，每一個單元的水面是一個平面，在整個滲流區(qū)將變成為由n個小三角平面組成的一個折曲面。用這個折曲面來代替地下水面?？芍?，單元越密折曲面越逼近地下水面，但是會造成計算工作量增大。

在三角面元上，水面是一個平面，可以用平面方程表示，以β面元為例：這就可以用來逼近

H（x，y，t）對面元水頭函數(shù)對x，y求導于是，面元上沿x，y方向的水力坡度為因此，求得a，b，c后，就可得面元上的水頭函數(shù)和水力坡度。

設β面元三節(jié)點的編號分別為i，j，m。

水頭H（x，y，t）在i，j，m節(jié)點上的水位分別為Hi（t），

Hj（t），

Hm（t）。

點元i，j，m的坐標為（xi，yi），（xj，yj），（xm，ym）。

現(xiàn)在的問題是：如果能確定出系數(shù)a，b，c，面元的水頭函數(shù)就確定了，已知條件是三點的坐標和三點的水位，就是，將a，b，c轉化成三點的坐標和三點的水位。按照克萊姆法則，有

面元三節(jié)點上的水位為：由知。或由克萊姆法則和行列式的計算知引入下列符號，也叫坐標輪減數(shù)。下面確定Da，Db，Dc，D。令則

將其代入水頭方程：得：令水頭方程變?yōu)椋涸诿恳粋€單元上，都定義著一個水頭模式，將這些水頭模式“拼接”起來，可近似表示滲流區(qū)的水頭函數(shù)H（x，y，t）因此，水頭方程為：其中，面元水力坡度：三、均衡單元的基本均衡方程

1、均衡單元的建立（1）首先進行三角網(wǎng)格剖分。（2）i節(jié)點所在均衡單元Ri的建立：Ri的圍線由以點i為公共頂點的面元外心的連線組成。

2、均衡單元的基本均衡方程

依據(jù)均衡原理：（