第1-3章基本原理

1、水文地質(zhì)數(shù)值計算水文地質(zhì)數(shù)值計算郭巧娜郭巧娜地球科學與工程學院地球科學與工程學院第一章第一章 地下水流動定解問題概述地下水流動定解問題概述1、三維流微分方程、三維流微分方程2、數(shù)學模型、數(shù)學模型 3、微分方程、微分方程4、邊界條件、邊界條件5、模型概化、模型概化6、實例、實例tzyvxtzyxx),(| )(方向流入tzyvxtzyxxx),(| )(方向流出),(| )(tzyxxxv),(| )(tzyxxv 取右圖所示的微小六面取右圖所示的微小六面體。設與體。設與x,y,z,方向?qū)姆较驅(qū)闹鳚B透系數(shù)分別為主滲透系數(shù)分別為Kxx, Kyy,Kzz；建立均衡期；建立均衡期 t時段內(nèi)，微

2、小均衡六面體時段內(nèi)，微小均衡六面體的水量守恒方程。的水量守恒方程。1、三維流微分方程、三維流微分方程1、三維流微分方程、三維流微分方程某一點為 , yyy , 某一點為xxx某一點為 , zzzx，y，z方向流入方向流入流出分別為流出分別為:l t時段內(nèi)，六面體水量變化量為：時段內(nèi)，六面體水量變化量為：tzyxzvyvxvzyx|)(|)(|)(tzyxxvtzyvvxxxxxx|tzxyyvtzxvvyyyyyy|tyxzzvtyxvvzzzzzz|六面體內(nèi)地下水儲存量的變化為六面體內(nèi)地下水儲存量的變化為由水均衡原理得由水均衡原理得1、三維流微分方程、三維流微分方程HzyxstHwzvyvx

3、vszyx)()()(Hzyxtzyxzvyvxvszyx|)(|)(|)(方程兩端除以方程兩端除以t，并取，并取,和和，則，則zHKzyHKyxHKxzzyyxxzHKvyHKvxHKvzzzyyyxxx一般密度的空間變化率很小，故一般密度的空間變化率很小，故于是有于是有()()()yxzvvvxyz由達西定律由達西定律有有（2）水流連續(xù)性方程左端項水流連續(xù)性方程左端項zHKzyHKyxHKxzzyyxx,都很小，可以忽略。都很小，可以忽略。1、三維流微分方程、三維流微分方程()esnM 其中為比儲水系數(shù))7(tHwzHKzyHKyxHKxszzyyxx上式為非均質(zhì)各向異性承壓含水層的偏微分

4、方程。上式為非均質(zhì)各向異性承壓含水層的偏微分方程。均質(zhì)各向異性非穩(wěn)定流均質(zhì)各向異性非穩(wěn)定流)10(222222tHwzHKyHKxHKszzyyxx0222222wzHKyHKxHKzzyyxx均質(zhì)各向異性穩(wěn)定流均質(zhì)各向異性穩(wěn)定流得到地下水三維流動微分方程得到地下水三維流動微分方程1/L1、三維流微分方程、三維流微分方程微分方程微分方程定解條件定解條件邊界條件邊界條件初始條件初始條件 已知已知t=0時的因變量，時的因變量，H(x,y,z,0)=H0(x,y,z)已知水頭邊界已知水頭邊界(I(I類邊界類邊界) )H(x,y,z,t)=f(x,y,z,t) (x,y,z) B1特例：特例：定水頭邊

5、界定水頭邊界 H(x,y,z,t)=C已知流量邊界已知流量邊界特例：隔水邊界特例：隔水邊界2B),(),(zyxtzyxfnH0nH2、數(shù)學模型、數(shù)學模型 數(shù)學模型數(shù)學模型地下水運動的數(shù)學模型結(jié)構(gòu)地下水運動的數(shù)學模型結(jié)構(gòu)(1) 第一類邊界條件(簡記為Dirichlet BC)n第一類邊界條件也稱為給定水頭邊界條件，指的是滲流區(qū)某一部分邊界上，各點水頭在某一時刻是已知的，表示為：n ),(),(1tzyxtzyxHDDtzyx),(2)第二類邊界條件(簡記為Neumann BC)n第二類邊界條件也稱為定流量邊界條件，是指某一部分邊界單位面積上流入(流出時值為負)的流量已知。相應的邊界條件表示為：

6、),(1tzyxqnHKNNtzyx),(注水井或抽水井也可以作為內(nèi)邊界來處理，此時的承壓含水層可忽略其水頭的垂向變化，注水井或抽水井開口于整個厚度，如圖所示，設井的半徑為 ，單位時間抽水量為 ，井壁法向 和 軸方向相反，如圖中所示。R)(tQ承壓抽水井示意圖承壓抽水井示意圖則由Darcy 定律，有：)()(2tQtKMRvrHKnHKt)(vRMtQrHKRr2)( （3） 第三類邊界條件(Rolin BC)n第三類邊界條件也稱為混合邊界條件，是指知道某段邊界上 和 的線性組合：HnHn臨海含水層與海水之間有一層薄淤泥層，則淤泥層兩側(cè)同一位置上的點存在水頭差，內(nèi)側(cè)含水層水頭 、滲透系數(shù)分別為

7、 ，淤泥層厚度為 ，水頭、滲透系數(shù)分別為 和 ，外側(cè)海水水頭為 ，忽略淤泥層內(nèi)孔隙水的彈性貯存的變化，在淤泥層與含水層交界面上根據(jù)流量的相等關(guān)系，有：HKmhKseaHn淤泥層與含水層交界面上根據(jù)流量的相等關(guān)系，mHHKnhKnHKseaseaHKmKHKmKnHKmKseaHKmK其中其中 （4）等值面邊界條件n 等值面邊界條件適用于不完全井的情況，即井壁和底部都有含水層中的水流滲入。等值面邊界條件的表達形式是： 已知，并有正、負之分，分別對應于單位時間內(nèi)的注入水量和抽出水量。 是井壁或渠壁與含水層的所有接觸面。)(tQdsnHK )(tQ（5） 自由面邊界條件n自由面邊界條件即潛水面邊界條

8、件，若潛水含水層中水頭為 ，自由面方程： ，那么，在自由面上，有：),(tzyxH),(tyxz),(),(,(tyxttyxyxH（5） 自由面邊界條件0),()tyxzezzzyyyyxxxxtnHKHKHK3、微分方程、微分方程tHwxHTxtHwyHTyxHTxtHwzHKzxHKxtHwzHKzyHKyxHKxxxyyxxszzxxszzyyxx一維流：平面二維流：剖面二維流：三維流： 數(shù)學模型數(shù)學模型4、邊界條件、邊界條件 邊界條件邊界條件：滲流區(qū)邊界上水力特征，即邊界上的水頭分布和變：滲流區(qū)邊界上水力特征，即邊界上的水頭分布和變化特征或流入流出含水層的水量分布和變化情況。主要有兩

9、類：化特征或流入流出含水層的水量分布和變化情況。主要有兩類： 1 已知邊界上的水頭分布規(guī)律已知邊界上的水頭分布規(guī)律 HB1=(x,y,t),其中其中(x,y,t)為已知函數(shù)。主要常見的是滲流為已知函數(shù)。主要常見的是滲流區(qū)與地表水體相接觸。區(qū)與地表水體相接觸。2已知邊界上的單位寬度流量已知邊界上的單位寬度流量q隨時間的變化規(guī)律隨時間的變化規(guī)律),(2tyxnHTB如已知流量為如已知流量為Q的承壓含水層中完整的抽水井，其井壁可以看作的承壓含水層中完整的抽水井，其井壁可以看作此類邊界。此類邊界。rQnHTB2 數(shù)學模型數(shù)學模型5例子例子:河間地塊承壓水流模型河間地塊承壓水流模型 設兩條河流平行、完全

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

11、軸原點位于左端河，右側(cè)為正向，設兩河流間距為軸原點位于左端河，右側(cè)為正向，設兩河流間距為L. L. 縱軸為水頭。縱軸為水頭。（3 3）數(shù)學模型）數(shù)學模型22eHHKMtx)(00 xHHt)(),(210tHtHLXxLx 00;0tLx0t河間地塊承壓水流模型（續(xù)）河間地塊承壓水流模型（續(xù)） n1 有限差分法原理有限差分法原理n2 導數(shù)的有限差分近似表示導數(shù)的有限差分近似表示n3 承壓一維流動有限差分法承壓一維流動有限差分法n4 承壓二維不穩(wěn)定流有限差分法承壓二維不穩(wěn)定流有限差分法n5 源匯項的處理及井孔水頭校正源匯項的處理及井孔水頭校正n6不規(guī)則邊界問題不規(guī)則邊界問題n7 矩形變格距網(wǎng)格差

12、分矩形變格距網(wǎng)格差分第二章第二章 有限差分法有限差分法1 有限差分法的基本原理有限差分法的基本原理將連續(xù)的問題離散后求解：將連續(xù)的問題離散后求解： 方法一方法一以地下水流基本微分以地下水流基本微分方程及其定解條件為基礎，方程及其定解條件為基礎， 在在滲流區(qū)剖分基礎上，用差商代滲流區(qū)剖分基礎上，用差商代替微商，將地下水流微分方程替微商，將地下水流微分方程的求解轉(zhuǎn)化為差分方程（代數(shù)的求解轉(zhuǎn)化為差分方程（代數(shù)方程）求解。方程）求解。 方法二方法二在滲流區(qū)剖分的基礎在滲流區(qū)剖分的基礎上，直接由達西定律和水均衡上，直接由達西定律和水均衡原理，建立各個均衡區(qū)的水均原理，建立各個均衡區(qū)的水均衡方程，即差分方

13、程。衡方程，即差分方程。矩形網(wǎng)格矩形網(wǎng)格多邊形網(wǎng)格多邊形網(wǎng)格1 有限差分法的基本原理有限差分法的基本原理網(wǎng)格劃分的基本類型網(wǎng)格劃分的基本類型n（1）先劃格線，）先劃格線，格點位于網(wǎng)格中心格點位于網(wǎng)格中心均衡網(wǎng)格均衡網(wǎng)格節(jié)點網(wǎng)格節(jié)點網(wǎng)格o（2）先規(guī)定格點位）先規(guī)定格點位置，再垂直平分兩相置，再垂直平分兩相鄰結(jié)點的連線作格線，鄰結(jié)點的連線作格線，形成的網(wǎng)格即為水均形成的網(wǎng)格即為水均衡區(qū)衡區(qū)基本概念基本概念為時間步長。通常稱為等或不相等的時段將連續(xù)的時間分割成相稱為空間步長。方向上的格距和方向上的格距nnttx,yyx將時間離散點和空間離散點聯(lián)合組成的網(wǎng)格稱為時空網(wǎng)格。將時間離散點和空間離散點聯(lián)合組

14、成的網(wǎng)格稱為時空網(wǎng)格。有限差分的基本原理：某點處水頭函數(shù)的導數(shù)用該點和幾個相鄰點處水頭值及其間有限差分的基本原理：某點處水頭函數(shù)的導數(shù)用該點和幾個相鄰點處水頭值及其間距近似表示。距近似表示。MODFLOW網(wǎng)格系統(tǒng)網(wǎng)格系統(tǒng)x000()()()f xxf xfxx 導導數(shù)數(shù)的的有有限限差差商商近近似似導數(shù)的定義導數(shù)的定義 當當非常小的時候，有非常小的時候，有 上式右端項即為上式右端項即為f(x)(x)在在x0 0處的差商。處的差商。 這樣定義的差商很容易理解，但不知道用差商代替微商所產(chǎn)生的誤差。下面利用泰勒公式導出差商及其誤差。xxfxxfxfx)()(lim)(0000方法一：方法一：差商代替微

15、商差商代替微商2 導數(shù)的有限差分近似表示導數(shù)的有限差分近似表示2000( )()()()()2!ff xxf xfxxx 2000( )()()()()2!ff xxf xfxxx 000()()()()f xxf xfxOxx 已知泰勒公式 由A得： AB 由B 得： 000()()()()f xf xxfxOxx 稱稱 為為f(x)在在x0處的處的一階前向差商，一階前向差商， 為為截斷誤差截斷誤差。xxfxxf)()(00)( xO 稱稱 為為f(x)在在x0處的處的一階后向差商，一階后向差商， 為為截斷誤差截斷誤差。xxxfxf)()(00)( xO 方法一方法一200002()2 ()

16、()()()()f xxf xf xxfxOxx 由A-B可以得： 由A+B可以得：2000()()()()2f xxf xxfxOxx (4)23400000()()( )()()()()()()2!3!4!fxfxff xxf xfxxxxx (4)23400000()()( )()()()()()()2!3!4!fxfxff xxf xfxxxxx AB稱稱 為為f(x)在在x0處的處的一階中心差商，一階中心差商， 為為截斷誤差截斷誤差。xxxfxxf2)()(002)( xO 稱稱 為為f(x)在在x0處的處的二階二階中心中心差商，差商， 為為截斷誤差截斷誤差。2000)()()(2)

17、(xxxfxfxxf2)( xO 方法一方法一n對于偏導數(shù)（偏微商），類似可以得到相應的差商對于偏導數(shù)（偏微商），類似可以得到相應的差商：ttxHttxHttxH),(),(),(000000 xtxHtxxHxtxH),(),(),(00000020000002002)(),(),(2),(),(xtxxHtxHtxxHxtxH方法一方法一),(|tzyxxxv),(|tzyxxvtzyvxtzyxx),(| )(方向流入tzyvxtzyxxx),(| )(方向流出 取右圖所示得微小六面體。設與取右圖所示得微小六面體。設與x,y,z,方向?qū)弥鳚B透系數(shù)分別為方向?qū)弥鳚B透系數(shù)分別為Kx,

18、 Ky,Kz；建立均衡期；建立均衡期 t時段內(nèi)，時段內(nèi)，微小均衡六面體的水量守恒方程。微小均衡六面體的水量守恒方程。方法二：達西定律和水均衡原理基于基于達西定律達西定律，x，y，z方向流入方向流入流出分別為流出分別為:l t時段內(nèi)，側(cè)向流入與源匯項導致六面體水量變化量為：時段內(nèi)，側(cè)向流入與源匯項導致六面體水量變化量為：tyxzyHHyHHtzxvvkjikjikjikjiyyyyy|, 1, 1,tyxzxHHxHHtzyvvkjikjikjikjixxxxx|, 1, 1tyxzzHHzHHtyxvvkjikjikjikjizzzzz|,1,1,ABC方法二：達西定律和水均衡原理DtyxzA

19、+B+C+D源匯項源匯項六面體內(nèi)地下水儲存量的變化為六面體內(nèi)地下水儲存量的變化為由水均衡原理得三維地下水流動方程的有限差分格式由水均衡原理得三維地下水流動方程的有限差分格式),(),(tzyxHttzyxHzyxHzyxssnkjinkjiskjikjikjikjikjikjikjikjikjikjikjikjiHHzyxtyxztyxzzHHzHHtyxzyHHyHHtyxzxHHxHH,1,1,1, 1, 1, 1, 1方法二：達西定律和水均衡原理有限差分法:三維（MODFLOW）差商代替微商差商代替微商3 承壓一維流動有限差分法ttxhxtxhT),(),(22ttxHttxHxtxxH

20、txHtxxH),(),()(),(),(2),(00002000000ttxHttxHxttxxHttxHttxxH),(),()(),(),(2),(00002000000顯式差分格式顯式差分格式隱式差分格式隱式差分格式方法一方法一1, 3 , 2 )(21211 nxitHHxHHHTninininini控制方程控制方程1, 3 , 2 )(21211111 nxitHHxHHHTninininini網(wǎng)格剖分網(wǎng)格剖分nx個個算例：顯式有限差格式 設兩條河流平行、完全切割含水層，含水層等厚、均質(zhì)設兩條河流平行、完全切割含水層，含水層等厚、均質(zhì)各向同性。各向同性。應用實例：河間地塊承壓水流模

21、型應用實例：河間地塊承壓水流模型3 承壓一維流動有限差分法承壓一維流動有限差分法步驟：（1 1）基礎資料的分析）基礎資料的分析（2 2）概念模型）概念模型（3 3）數(shù)學模型）數(shù)學模型（4 4）數(shù)值方法及計算機程序）數(shù)值方法及計算機程序（5 5）參數(shù)）參數(shù)（6 6）結(jié)果分析）結(jié)果分析 3 承壓一維流動有限差分法承壓一維流動有限差分法建立數(shù)學模型0,0tLxLx 022eHHKMtx)(00 xHHt)(),(210tHtHLXx（1 1）模型概化）模型概化 由所述水文地質(zhì)條件，可以概化為一維承壓水流問題。由所述水文地質(zhì)條件，可以概化為一維承壓水流問題。（2 2）建立坐標系）建立坐標系（如圖），將

22、地下水流動系統(tǒng)空間結(jié)構(gòu)放在坐標系內(nèi)，（如圖），將地下水流動系統(tǒng)空間結(jié)構(gòu)放在坐標系內(nèi)，從而量化各變量的取值范圍。本例，取從而量化各變量的取值范圍。本例，取x-x-軸原點位于左端河，右側(cè)軸原點位于左端河，右側(cè)為正向，設兩河流間距為為正向，設兩河流間距為L. L. （3 3）數(shù)學模型）數(shù)學模型0t3 承壓一維流動有限差分法承壓一維流動有限差分法差分方程及其解法顯式格式將（0L）分成 N 等份，Nlx 1)1)網(wǎng)格剖分：網(wǎng)格剖分：xixittntn取時間步長取時間步長 ，記，記 （n n=0=0、1 1、2 2、3 3、44）記 ，(i=0,1,2,3，4N)2 2）建立差分方程）建立差分方程：在網(wǎng)格

23、系統(tǒng)中任意取一點：在網(wǎng)格系統(tǒng)中任意取一點),(nitx),(txH),(nitx設設是問題的解，則在是問題的解，則在處有處有記為（記為（i i，n n）22( , )nneiniiHHKMx ttx3 承壓一維流動有限差分法承壓一維流動有限差分法22( , )nneiniiHHKMx ttx1()nnniiiHHHOttt2211222()()nnnniiiiHHHHOxxx用差商代替微商用差商代替微商： 將上述兩式舍去余項，代入方程并記將上述兩式舍去余項，代入方程并記niHnih為為11122(,)()nnnnniiiiiinehhhhhKMx txt2)()(xOtO顯然該式具有截斷誤差得

24、到得到顯式格式（續(xù)1）111(12 )( , )nnnniiiiinethhhhx t引入無量綱變量引入無量綱變量：將該式子代入得到：將該式子代入得到： （i=1，2，3，.N-1），（n=1，2，3，.） 顯式格式（續(xù)2）2)( xtKMe11122(,)()nnnnniiiiiinehhhhhKMx txt3）顯示差分方程的求解n計算各結(jié)點初始時刻水頭值計算各結(jié)點初始時刻水頭值n利用差分方程計算各結(jié)點利用差分方程計算各結(jié)點t1t1時刻水頭值時刻水頭值n利用邊界條件計算邊界結(jié)點水頭值利用邊界條件計算邊界結(jié)點水頭值n重復重復2 2、3 3步，直到計算出擬計算的各個時刻的水頭步，直到計算出擬計算

25、的各個時刻的水頭值值)(00 xHHt顯式格式（續(xù)顯式格式（續(xù)3）22eHHKMtx)(00 xHHt)(),(210tHtHLXx算例（續(xù)4）004. 0em10)(m,20)( m,20 m/d,821ttMK在上述模型中，設在上述模型中，設L=1000米米取空間步長為取空間步長為200200米，時間步長為米，時間步長為0.250.25天，分別計算各節(jié)點天，分別計算各節(jié)點各各時刻的時刻的水頭值。水頭值。0 x20m,0 x10m,)(0 xH4/1)200(*004. 025. 0*20*8)(22xtKMeTime/dayx=0 mx=200 mx=400 mx=600 mx=800 m

26、x=1000 m02010101010100.252012.5101010100.502013.7510.6251010100.752014.531 11.250 10.156 10101.002015.078 11.797 10.391 10.039 101.252015.488 12.266 10.654 10.117 10nieninininithhhh,111/)21 ()2(4141)4121 (4111111nininininininihhhhhhh算例（續(xù)5）Time/dayx=0 mx=200 mx=400 mx=600 mx=800 mx=1000 m020101010101

27、00.252020101010100.502010201010100.7520300 2010101.0020-10 101.252010算例（續(xù)6）1)200(*004. 01*20*8)(22xtKMenieninininithhhh,111/)21 (nininininininihhhhhhh11111) 121 (如果如果 t=1,t=1,則則算例：隱式格式22eHHKMtx)(00 xHHt)(),(210tHtHLXx004. 0e在上述模型中，設在上述模型中，設L=1000米米取空間步長為取空間步長為200米，時間步長為米，時間步長為0.25天，用隱式差分格式計算各節(jié)天，用隱式差

28、分格式計算各節(jié)點個時刻的水頭值。點個時刻的水頭值。m10)(m,20)( 10m,)( m,20 m/d,821ttxHMK11122332210111112131211/./.2121.212121nNnNenNnNenNnennennnennNnNnnnhththththhthhhhhh4/1)200(004. 025. 0208)(22xtKMe在這個例子中，在這個例子中，.0,1,2,3,4, m,10,m2050nhhnnm10,m10,m10,m1004030201hhhh154321011413121125. 025. 0225. 00025. 0225. 00025. 0225

29、. 00025. 02nnnnnnnnnnhhhhhhhhhh解：隱式格式一般方程為解：隱式格式一般方程為于是有于是有根據(jù)初始條件得根據(jù)初始條件得根據(jù)邊界條件得根據(jù)邊界條件得.0,1,2,3,4, m,10,m2050nhhnnm10,m10,m10,m1004030201hhhh1504030210011413121125. 025. 0225. 00025. 0225. 00025. 0225. 00025. 02hhhhhhhhhh由初始條件和邊界條件由初始條件和邊界條件5 .121010151025. 01010102025. 010225. 00025. 0225. 00025. 0

30、225. 00025. 0214131211hhhh由此解得由此解得t1t1時刻的水頭值為時刻的水頭值為m10,m10,m10,m1014131211hhhh在上述方程中取在上述方程中取 n=0,n=0,可以得到計算可以得到計算t1t1時刻水頭值的方程時刻水頭值的方程所以上述方程變成所以上述方程變成同理，可計算同理，可計算t2t2時刻的時刻的水頭值水頭值差分方程求解n一維顯式差分格式 )(2)(2 )(2)(2111212212212331324322122321tHHxHHHtHHxHHHtHHxHHHtHHxHHHnninninninninninninninninninninnnnnnnn

31、nn網(wǎng)格個數(shù)為網(wǎng)格個數(shù)為ni直接求解直接求解差分方程求解n一維隱式差分格式一維隱式差分格式 )(2)(2 )(2)(2111211112212211121331321413122122131211tHHxHHHtHHxHHHtHHxHHHtHHxHHHnninninninninninninninninninninnnnnnnnnn網(wǎng)格個數(shù)為網(wǎng)格個數(shù)為ni迭代求解迭代求解方方程程組組PCGSIPSORWHSSAMGGMGMODFLOW差分方程的收斂性和穩(wěn)定性n截斷誤差截斷誤差：用差商代替微商時，地下水流動方程產(chǎn)生用差商代替微商時，地下水流動方程產(chǎn)生的誤差為截斷誤差。的誤差為截斷誤差。n收斂性收斂

32、性：當空間步長和時間步長趨于當空間步長和時間步長趨于0 0時，有限差分方時，有限差分方程的精確解趨于地下水流動問題微分方程定解問題的程的精確解趨于地下水流動問題微分方程定解問題的精確解。則稱該差分格式是收斂的。精確解。則稱該差分格式是收斂的。n穩(wěn)定性穩(wěn)定性：如果在求解差分方程過程中，某時間步引入如果在求解差分方程過程中，某時間步引入某個誤差，而在以后的各時段計算中，該誤差不再擴某個誤差，而在以后的各時段計算中，該誤差不再擴大，則稱該差分格式是穩(wěn)定的。大，則稱該差分格式是穩(wěn)定的。102一維一維顯示格式顯示格式的收斂條件和穩(wěn)定條件是：的收斂條件和穩(wěn)定條件是：2()()OtOx第二類邊界條件的處理n

33、前面介紹的幾種差分格式，都是以第一類邊界條件為例，說明其具體計算方法的。如果是第二類邊界條件，怎么使用這些格式進行計算呢？為此，我們以下面的定解問題來討論這個問題。設問題數(shù)學模型為出量。入量；取負值時，為流當流量取正值時，為流處的單寬流量。分別表示左、右兩邊界其中)(),(0)()(0)(0,0212100022tqtqttqxHTtqxHTLxxHHtLxtHxHTLxxt)()(),.,2 , 1 , 0()(,)(),.,2 , 1 , 0()()(0,211102121010100nnNnNnnnnnNnNnNnnNnnnnnnLxxtqTxhhtqTxhhMntqxhhTxHtqxH

34、TLxMntqxhhTxHtqxHTxxHxH得用后向差分公式代替，有對于右邊界，得用前向差分公式代替處，有。對于左邊界用后向差分公式代替用前向差分公式代替界條件中的偏導數(shù)，即用單側(cè)差分公式代替邊用一階中心差分公式代替邊界條件中的偏導數(shù),.)1 , 0()(2)(2)(2)(2-,21111121111110111ntqTxhhtqTxhhtqxhhTtqxhhTxxxxxxxxnnNnNnnnnnNnNnnnNNN于是由邊界條件得他們滿足和為的邊界點到虛線位置，記新擴充為此，首先將邊界擴充4 承壓二維不穩(wěn)定流有限差分法n我們考慮無垂向補給、排泄、均質(zhì)、各向同性、等厚的承壓流動問題，邊界條件第

35、一類，滲流區(qū)為矩形，則可用下述定解問題描述),(),(),(),()0 ;0(),()0 ;0(),(210210002222txHtxHtyHtyHbyaxyxHHbyaxtyxyHxHKMtHayyaxxte4 承壓二維不穩(wěn)定流有限差分法（1）顯式格式-網(wǎng)格剖分.)3 , 2 , 1 , 0(,),.2 , 1 , 0()b, 0(),.2 , 1 , 0(), 0(ntnttNjyjyNbyNNixixNaxNanyjyyxixx記取時間步長等份，記分成將等份，記分成將（1）顯式格式2）建立差分方程：22,1,21,1,2, 1, 1,22,22,)()()(),()(2)(2),(),

36、(),(),(),(yxtthhtyxyhhhxhhhKMtyxyHxHKMtHtxtxHnjityxnjinjienjinjinjinjinjinjinjinjinjinjinjieninji顯然該式具有截斷誤差用差商代替微商，得到處有是問題的解，則在設記為點在網(wǎng)格系統(tǒng)中任意取一（1）顯式格式.)3 , 2 , 1 , 0() 1N,.3 , 2 , 1, 1.N, 3 , 2 , 1(),()221 ()(,)(yx,1,1, 1, 11,22njityxthhhhhhytKMxtKMnjienjiyxnjiynjiynjixnjixnjieyex將該式子代入得到：引入無量綱變量：（1）顯

37、式格式3）二維顯示差分方程的求解計算各結(jié)點初始時刻水頭值計算t1時刻水頭值計算邊界結(jié)點水頭值重復2、3步，直到計算出擬計算的各個時刻水頭值（1）顯式格式n差分方程收斂性21)()(2yxxoto穩(wěn)定條件是隱式格式的收斂條件和截斷誤差（2）隱式格式.)3 , 2 , 1; 1,.3 , 2 , 1; 1,.2 , 1(),()221 (-)(,)()()()(),()(2)(2),(1,1,1, 1, 122221 -,21,1,2, 1, 1,22,22,nNyjNxityxthhhhhhytKMxtKMyxtthhtyxyhhhxhhhKMtyxyHxHKMtHnjienjinjiyxnji

38、ynjiynjixnjixeyexnjinjienjinjinjinjinjinjinjinjinjinjinjie將該式子代入得到：引入無量綱變量：顯然該式具有截斷誤差用差商代替微商，得到對于（2）隱式格式4）差分方程的收斂性和穩(wěn)定性條件0, 0)()(2yxxoto穩(wěn)定條件是隱式格式的收斂條件和截斷誤差5 源匯項的處理及井孔水頭校正越流、入滲和抽水井等問題的處理越流、入滲和抽水井等問題的處理n如果考慮垂直滲流項（即源匯項），則二維承壓流動如果考慮垂直滲流項（即源匯項），則二維承壓流動微分方程可寫成微分方程可寫成n建立差分方程時，在結(jié)點處應加上這一項，它可具體建立差分方程時，在結(jié)點處應加上這

39、一項，它可具體表示為：表示為：2222HHHTTxy,()i juzi ji ji ji jzQKhhWMx y 為抽水或注水流量。；為入滲強度或蒸發(fā)強度為主含水層水頭；為相鄰含水層水頭；為相鄰弱透水層厚度；透系數(shù)；為相鄰弱透水層垂向滲流等）；越流、入滲或蒸發(fā)、井為垂向滲流強度（包括jijijiujizzjiQWhhMK,井水位校正對比圖井水位校正對比圖 圖圖3-53-5（a a）初始網(wǎng)格）初始網(wǎng)格, ,()22wi j kwQxhhInTr有限差分法計算井水位的校正有限差分法計算井水位的校正 圖圖3-53-5（b b）加密網(wǎng)格）加密網(wǎng)格6不規(guī)則邊界問題不規(guī)則邊界問題l當研究區(qū)的幾何形狀屬于簡

40、當研究區(qū)的幾何形狀屬于簡單形式（如矩形滲流區(qū)）時，單形式（如矩形滲流區(qū)）時，差分網(wǎng)格的劃分往往將結(jié)點差分網(wǎng)格的劃分往往將結(jié)點設在邊界上。設在邊界上。l然而，對于實際問題來說，然而，對于實際問題來說，邊界通常不是那么規(guī)則，邊邊界通常不是那么規(guī)則，邊界的某些部分，甚至大部分界的某些部分，甚至大部分不能與結(jié)點重合，我們稱這不能與結(jié)點重合，我們稱這種邊界是不規(guī)則的。種邊界是不規(guī)則的。l關(guān)于不規(guī)則邊界問題，直接關(guān)于不規(guī)則邊界問題，直接取最靠近邊界的曲折格線為取最靠近邊界的曲折格線為近似邊界近似邊界河流、水庫、湖泊等問題河流、水庫、湖泊等問題7 矩形變格距網(wǎng)格差分n本節(jié)直接由達西定律和水均衡原理建立差分方程n達西定律：n水均衡原理： 對某一研究對象，流入-流出=體系內(nèi)質(zhì)量變化量 研究對象可以是大區(qū)域，也可以是微分單元體。 大區(qū)域的水均衡計算經(jīng)常用于區(qū)域的水資源評價 表示當水頭下降一個單位時，從單位體積飽水空隙介質(zhì)中釋放的水量或體積-KdH/dxq nn)(sss參考教材參考教材n周志芳主編，地下水運動數(shù)值模擬，自編教材，2003n薛禹群，謝春紅，水文地質(zhì)學的數(shù)值法，煤炭工業(yè)出版社，1980n陳崇希，唐仲華