1、7.1 7.1 概述概述7.2 7.2 差分的概念及分類差分的概念及分類7.3 7.3 顯式差分、隱式差分和中心差分顯式差分、隱式差分和中心差分7.4 7.4 土壤水分運(yùn)移方程的線性化方法土壤水分運(yùn)移方程的線性化方法7.5 7.5 線性化方法與土壤水分參數(shù)的取值線性化方法與土壤水分參數(shù)的取值7.6 7.6 邊界條件的處理及追趕法求解三對角方程邊界條件的處理及追趕法求解三對角方程7.7 7.7 垂直一維非飽和土壤水流計算流程圖垂直一維非飽和土壤水流計算流程圖Huyakorn, P. S. and G. F. Pinder. Computational methods in subsurface

2、flow. Academic Press, New York, 1983.薛禹群，謝春紅，水文地質(zhì)學(xué)的數(shù)值法，煤炭工業(yè)出版社，薛禹群，謝春紅，水文地質(zhì)學(xué)的數(shù)值法，煤炭工業(yè)出版社，1980.土壤水分運(yùn)動方程的求解方法：土壤水分運(yùn)動方程的求解方法： 解析解（解析解（Analytical solution） 數(shù)值解（數(shù)值解（Numerical solution）(1) 有限差分法（有限差分法（Finite Difference method）(2) 有限單元法（有限單元法（Finite Element method）7.1 7.1 概述概述 差分的由來（差分的由來（Taylor展開）展開） if =

3、 (x) exists, thenkkkxkxxxxxx!)()(kkkkxkxxxxxx!) 1()()(1)(2)(0)()(xxxxxx)(0)()(xxxxxxwhen x 0, thenxxxxxxxxxxx)()(lim)()(lim00 向前差分向前差分xxxxx)()( 向后差分向后差分xxxxx)()( 中心差分（中心差分（Eqn.(1) - Eqn. (2)）xxxxxxxx)(02)()(2xxxxxx2)()(i.e.截斷誤差截斷誤差0 ( x)截斷誤差截斷誤差0 ( x)截斷誤差截斷誤差0 ( x2)0 x ABCEFxi-1xixi+1如圖所示：如圖所示： BC表示

4、向前差分；表示向前差分； AB表示向后差分；表示向后差分； AC表示中心差分；表示中心差分；In addition, Eqn.(1) + Eqn.(2) )(0)(2)()(22222xxxxxxxxx)(0)()(2)(2222xxxxxxxxi.e.:其截斷誤差亦為與其截斷誤差亦為與 x2同階的無窮小量。同階的無窮小量。i=1, z(1)=0i=2i=3.i=n, z(n)=Lz考虛垂直一維問題（如右圖所示），考虛垂直一維問題（如右圖所示），泛定方程可寫為（以泛定方程可寫為（以 方程為例）：方程為例）：剖面深度為剖面深度為L，共剖分出，共剖分出n個節(jié)點(diǎn)（個節(jié)點(diǎn)（i = 1, 2, 3, ,

5、 n），），z (1) = 0, z (n) = L，已，已知知t = tk時刻各節(jié)點(diǎn)含水率時刻各節(jié)點(diǎn)含水率 的分布，求的分布，求t = t k+1時刻各節(jié)點(diǎn)含水率時刻各節(jié)點(diǎn)含水率 的分布，的分布， kkttt1zKzDzt 7.3.17.3.1顯式差分格式顯式差分格式kiizKzDztttkikii121212121iikikikizzzDzDzDz iikikikikizzDzD112121 112121iikikikikizzDzD7.3.17.3.1顯式差分格式顯式差分格式 2112111211121 zzzzDzzDzDziiikikikiiikikikiki21212121iiki

6、kikizzKKzK21212121112111211iiiiiikikiiiikikiikikizzKKzzDzzDt)(21112121iiiizzzz) 1, 3, 2( 2212111211121111niKKzzDzzDzztiiiikikiiiikikiiiikiki7.3.17.3.1顯式差分格式顯式差分格式 若取等步長若取等步長 z，即，即 z =zi+1 - zi = zi - zi-1，則有：，則有：) 1, 3, 2(212112112121niKKzDDztiikikiikikiikiki7.3.17.3.1顯式差分格式顯式差分格式 顯式差分格式是有條件收斂的顯式差分格

7、式是有條件收斂的，一般應(yīng)滿足，一般應(yīng)滿足max2Dzrt其中：其中：r多取多取1/2。由于土壤接近飽和時，。由于土壤接近飽和時，Dmax很很大，故大，故 t一般要求取值很小，耗費(fèi)機(jī)時。一般要求取值很小，耗費(fèi)機(jī)時。7.3.17.3.1顯式差分格式顯式差分格式 1kiiKzDzt) 1, 3, 2( 22/12/1111121111121111niKKzzDzzDzztiiiikikiiiikikiiiikiki7.3.27.3.2隱式差分格式隱式差分格式 )(1111iiiizzzztr)(1112iiiizzzztrLetikiikiikiiHGFE11111then ( i = 2, 3,

8、, n - 1)三對角方程三對角方程7.3.27.3.2隱式差分格式隱式差分格式 where;211212iiiiDrGDrE11211212iiiiiDrDrGEF2121112iiiikiiKKzztH隱式差分格式是無條件收斂的。隱式差分格式是無條件收斂的。7.3.27.3.2隱式差分格式隱式差分格式 kikikiiKzDzKzDzKzDzt12/1217.3.37.3.3中心差分格式中心差分格式 ) 1, 3, 2()(2)(22212111111211111121111niKKzzDzzDzztiiiikikikikiiiikikikikiiiikiki7.3.37.3.3中心差分格式

9、中心差分格式 )(1111iiiizzzztr)(1112iiiizzzztrLetikiikiikiiHGFE11111then ( i = 2, 3, , n - 1)where;211212iiiiDrGDrE7.3.37.3.3中心差分格式中心差分格式 11211212iiiiiDrDrGEF中心差分格式也是無條件收斂的。中心差分格式也是無條件收斂的。212111112) 1(iiiikiikiiikiiiKKzztGGEEH7.3.37.3.3中心差分格式中心差分格式 7.4 7.4 土壤水分運(yùn)移方程的線性化方法土壤水分運(yùn)移方程的線性化方法 以基質(zhì)勢水頭以基質(zhì)勢水頭h為因變量的一維垂

10、直向土壤為因變量的一維垂直向土壤水分運(yùn)動基本方程為：水分運(yùn)動基本方程為： K hhhC hK htzzz7.4 7.4 土壤水分運(yùn)移方程的線性化方法土壤水分運(yùn)移方程的線性化方法7.4.17.4.1顯式差分格式顯式差分格式相應(yīng)h方程的差分方程為：11/211/21211kkkkkkkkkiiiiiiiiiChhKhhKhhtz1/21/2kkiiKKz7.4 7.4 土壤水分運(yùn)移方程的線性化方法土壤水分運(yùn)移方程的線性化方法7.4.27.4.2隱式差分格式隱式差分格式相應(yīng)h方程的差分方程為1111111111/211/2111/21/22kkkkkkkkkkiiiiiiiikiiihhKhhKhh

11、KKCtzz7.4 7.4 土壤水分運(yùn)移方程的線性化方法土壤水分運(yùn)移方程的線性化方法7.4.37.4.3中心差分格式中心差分格式相應(yīng)h方程的差分方程為1111111221221122111111122212222kkkkkkkkiiiiiiiikkkkkkiikiiiiihhhhhhCKtzKKhhhhKzz 土壤水分運(yùn)移方程中的各參數(shù)均依賴于變量土壤水分運(yùn)移方程中的各參數(shù)均依賴于變量 或或 m，從而使得方程呈現(xiàn)出較強(qiáng)的非線性性，求解前必須將其線從而使得方程呈現(xiàn)出較強(qiáng)的非線性性，求解前必須將其線性化，得到性化，得到n元一次線性代數(shù)方程組，以便于求解。元一次線性代數(shù)方程組，以便于求解。（1）顯式

12、線性化）顯式線性化 計算過程中，方程中的各參數(shù)如：計算過程中，方程中的各參數(shù)如：C ( m)，K( m) or D( )，K( )等均以時段初的值（即前一時段的值）代入。等均以時段初的值（即前一時段的值）代入。 適用條件：此法只能用于適用條件：此法只能用于 m或或 變化緩慢的情況；當(dāng)變化緩慢的情況；當(dāng)變化劇烈時，此法會導(dǎo)致較大的偏差。變化劇烈時，此法會導(dǎo)致較大的偏差。7.5.17.5.1線性化方法線性化方法例如：已知例如：已知k - 1與與k時刻各節(jié)點(diǎn)的時刻各節(jié)點(diǎn)的 值，需求值，需求k + 1時刻各時刻各節(jié)點(diǎn)的節(jié)點(diǎn)的 值，應(yīng)用顯式外推線性化參數(shù)的方法如下（以值，應(yīng)用顯式外推線性化參數(shù)的方法如下

13、（以導(dǎo)水率導(dǎo)水率K的求法為例，其他各參數(shù)的求法與此類似）：的求法為例，其他各參數(shù)的求法與此類似）：1111)(kkkkkikikikitttt(*) 求求 ： 方法（方法（a）：）：21kiK)(21121kikiki2121kikiKK方法方法(b)：算術(shù)平均算術(shù)平均:幾何平均幾何平均:調(diào)和平均調(diào)和平均:121)(21kikikiKKK)()(121kikikiKKK)()()()(21121kikikikikiKKKKK 求求 ： 方法（方法（a）：）：kiK21)(21121kikikikikiKK2121方法方法(b)：算術(shù)平均算術(shù)平均幾何平均幾何平均調(diào)和平均調(diào)和平均kikikiKKK

14、121)(21)()(121kikikiKKK)()()()(21121kikikikikiKKKKK 計算時，先用顯式線性化的方法，采用隱式差分格式，計算時，先用顯式線性化的方法，采用隱式差分格式，求出求出k + 1/2時刻各節(jié)點(diǎn)的時刻各節(jié)點(diǎn)的 m或或 （即（即 or ）；根據(jù)）；根據(jù)這些值可求得相應(yīng)的這些值可求得相應(yīng)的 or 等參數(shù)，以此作為計算時始末（等參數(shù)，以此作為計算時始末（tk to tk+1）的平均值，代）的平均值，代入方程求解入方程求解 ，從而完成方程的線性化工作。，從而完成方程的線性化工作。21kmi21ki2121kmikiKK21kiK1ki7.5.17.5.1線性化方法

15、線性化方法7.5.1.27.5.1.2預(yù)報校正法預(yù)報校正法（1）顯式外推法（先求顯式差分方程，得）顯式外推法（先求顯式差分方程，得 作為預(yù)報作為預(yù)報值）值） 計算參數(shù)時，據(jù)前一時段始末的函數(shù)值計算參數(shù)時，據(jù)前一時段始末的函數(shù)值 m或或 ，用，用線性外推近似求出本時段末的函數(shù)值，然后線性外推近似求出本時段末的函數(shù)值，然后 方法方法(a)：求出計算時段始末的平均值，再通過各參數(shù)與：求出計算時段始末的平均值，再通過各參數(shù)與 m或或 的關(guān)系，求出各參數(shù)；或的關(guān)系，求出各參數(shù)；或 方法方法(b)：先求出計算時段始末的參數(shù)值，然后對參數(shù)取：先求出計算時段始末的參數(shù)值，然后對參數(shù)取平均。平均。 平均值的計算

16、可采用算術(shù)平均、幾何平均、調(diào)和平平均值的計算可采用算術(shù)平均、幾何平均、調(diào)和平均等方法。均等方法。7.5.17.5.1線性化方法線性化方法7.5.1.27.5.1.2預(yù)報校正法預(yù)報校正法1ki（2）線性外推法（已知線性外推法（已知 時刻，外推時刻，外推 時刻）時刻） 假定含水率在相鄰時段內(nèi)是線性變化的，根據(jù)前假定含水率在相鄰時段內(nèi)是線性變化的，根據(jù)前一個時段末的含水率，用直線外推近似求出計算時段一個時段末的含水率，用直線外推近似求出計算時段末的含水率：末的含水率：7.5.17.5.1線性化方法線性化方法7.5.1.27.5.1.2預(yù)報校正法預(yù)報校正法,1k k 1k 112kkkiii11231

17、22kkkiii 對于第一時段，無法應(yīng)用線性外推法，可用其他對于第一時段，無法應(yīng)用線性外推法，可用其他線性化方法。線性化方法。 7.5.17.5.1線性化方法線性化方法7.5.1.27.5.1.2預(yù)報校正法預(yù)報校正法（3 3）顯式預(yù)測法（兩次求解隱式：第一次預(yù)報值；）顯式預(yù)測法（兩次求解隱式：第一次預(yù)報值；第二次得校正值）第二次得校正值） 首先利用顯式線性化方法，由時段初的含水率首先利用顯式線性化方法，由時段初的含水率及相應(yīng)的參數(shù)值，按隱式差分格式求解方程，求及相應(yīng)的參數(shù)值，按隱式差分格式求解方程，求得本時段末（或時段中間）的含水率及相應(yīng)的參得本時段末（或時段中間）的含水率及相應(yīng)的參數(shù)值（預(yù)報

18、值），再求解隱式差分方程，得到時數(shù)值（預(yù)報值），再求解隱式差分方程，得到時段末的含水率的值（校正值）。段末的含水率的值（校正值）。計算時，先假定本時段末（計算時，先假定本時段末（tk+1）的）的 m或或 （一般做（一般做法是先取時段初的函數(shù)值法是先取時段初的函數(shù)值 or ） ， or 等參數(shù)的迭代初值等參數(shù)的迭代初值 求解線性代數(shù)方程組得求解線性代數(shù)方程組得 or 直至前后兩次迭代值的誤差滿足精度要求為止。直至前后兩次迭代值的誤差滿足精度要求為止。kmiki)1(1kiK)1(1kiD)1(1kiC)1(1kmi)1 (1ki)2(1)2(1)2(1or ,kikikiCDK)2(1)2(1o

19、r kikmi)(1)(1or pkipkmi7.5.1.3迭代法迭代法（a）若取絕對誤差，判斷條件為：）若取絕對誤差，判斷條件為：), 2, 1(max1) 1( 1)( 1nipkipkii), 2, 1(max1)1(1)(1nipkmipkmiior7.5.17.5.1線性化方法線性化方法7.5.1.3迭代法迭代法（b）若取相對誤差，判斷條件為：）若取相對誤差，判斷條件為：), 2, 1(max2)(1)1(1)(1nipkipkipkii), 2, 1(max2)(1)1(1)(1nipkmipkmipkmiior7.5.17.5.1線性化方法線性化方法7.5.1.3迭代法 1與與

20、2 的取值可根據(jù)實際問題的精度要求而定。的取值可根據(jù)實際問題的精度要求而定。 此法計算稍復(fù)雜，但精度較高，數(shù)值計算時此法計算稍復(fù)雜，但精度較高，數(shù)值計算時經(jīng)常采用。為避免迭代次數(shù)太多，一般可設(shè)置一經(jīng)常采用。為避免迭代次數(shù)太多，一般可設(shè)置一個判斷語句：當(dāng)?shù)螖?shù)個判斷語句：當(dāng)?shù)螖?shù)m5時，將時間步長時，將時間步長 減半，再重新計算。減半，再重新計算。kkkttt117.5.17.5.1線性化方法線性化方法7.5.1.3迭代法7.5 線性化方法與土壤水分參數(shù)的取值線性化方法與土壤水分參數(shù)的取值7.5.2.1直接法直接法先求出兩結(jié)點(diǎn)的參數(shù)，然后先求出兩結(jié)點(diǎn)的參數(shù)，然后(1)取兩結(jié)點(diǎn)參數(shù)的算術(shù)平均值

21、：取兩結(jié)點(diǎn)參數(shù)的算術(shù)平均值：7.5.2參數(shù)取值參數(shù)取值11212iiiDDD(2)取兩結(jié)點(diǎn)參數(shù)的調(diào)和平均值：取兩結(jié)點(diǎn)參數(shù)的調(diào)和平均值：11122/iiiiiDDDDD(3) 取兩結(jié)點(diǎn)參數(shù)的幾何平均值：取兩結(jié)點(diǎn)參數(shù)的幾何平均值：112iiiDDD7.5 線性化方法與土壤水分參數(shù)的取值線性化方法與土壤水分參數(shù)的取值7.5.2 2間接法間接法先求出兩結(jié)點(diǎn)處含水率的算術(shù)、調(diào)和、幾何平均值，先求出兩結(jié)點(diǎn)處含水率的算術(shù)、調(diào)和、幾何平均值，然后再求出相應(yīng)的參數(shù)值。然后再求出相應(yīng)的參數(shù)值。(1)算術(shù)平均值：算術(shù)平均值：11212iii1122iiDD(2)調(diào)和平均值：調(diào)和平均值：11122/i iiii112

22、2iiDD(3)幾何平均值：幾何平均值：112i ii1122iiDD7.5 線性化方法與土壤水分參數(shù)的取值線性化方法與土壤水分參數(shù)的取值7.5.2 3其他方法其他方法（1）求出的參數(shù)和兩結(jié)點(diǎn)處的參數(shù)和，三點(diǎn)參）求出的參數(shù)和兩結(jié)點(diǎn)處的參數(shù)和，三點(diǎn)參數(shù)取平均，故可稱為數(shù)取平均，故可稱為“三點(diǎn)式三點(diǎn)式”：11122124iiiiDDDD（2） 11121iiiiiDDd （3）由鄰近結(jié)點(diǎn)的參數(shù)代替，如）由鄰近結(jié)點(diǎn)的參數(shù)代替，如11122,iiiiDD DD前已敘及，采用顯式差分格式可直接求解，采前已敘及，采用顯式差分格式可直接求解，采用隱式或中心差分格式時總可以得到如下形式的三對用隱式或中心差分格

23、式時總可以得到如下形式的三對角方程組（以角方程組（以 方程為例）：方程為例）：共共n - 2個方程，再加上上、下邊界條件，即可組成個方程，再加上上、下邊界條件，即可組成n n 階的階的n元一次方程組，線性化后求解該方程組可獲得元一次方程組，線性化后求解該方程組可獲得原定解問題的解。原定解問題的解。) 1, 3 , 2( 11111niHGFEikiikiikii7.6.17.6.1邊界條件的處理邊界條件的處理 上、下邊界不論為第一、第二或第三類邊界，上、下邊界不論為第一、第二或第三類邊界，總可以化為：總可以化為：1121111HGFkknknnknnHFE111(2)(1)7.6.17.6.1

24、邊界條件的處理邊界條件的處理 驗驗 證：證： 當(dāng)上、下邊界為第一類邊界時當(dāng)上、下邊界為第一類邊界時)()()(),()(), 0(121111121kknkkztftftftLtft)(, 0, 1)2.(Eqn)(, 0, 1) 1.(Eqn12n11111knnktfHEFtfHGF7.6.17.6.1邊界條件的處理邊界條件的處理1121111HGFkknknnknnHFE111 當(dāng)上、下邊界為第二、三類邊界時當(dāng)上、下邊界為第二、三類邊界時)4()()3(), 0(210tqKzDttqKzDzLzz)(),(122/111112111112312111223knnnknknnkkkktq

25、KzzDtqKzzD1121111HGFkknknnknnHFE111211212123111111112231)(,),(,nknnnnnnnkkKtqHEFzzDEKtqHFGzzDFnknnknnnknnknnknnkkkkkkkkHFEHGFEHGFEHGFEHGF111111111121314313312321321221121121111經(jīng)線性化后，原定解問題可化為經(jīng)線性化后，原定解問題可化為n n階線性代數(shù)方程組，即：階線性代數(shù)方程組，即：以矩陣的形式表示即為：以矩陣的形式表示即為：nnknknkkknnnnnHHHHHFEGFEGFEGFEGF13211111312111113

26、33222117.6.27.6.2 追趕法求解三對角方程追趕法求解三對角方程1原理原理1121111HGFkk(1)(1a) ,( 111111111111HBFAAGABkk7.6.27.6.2 追趕法求解三對角方程追趕法求解三對角方程1原理原理2132122112HGFEkkk1122132121212ABEHGAEGFkk13222212kkAGAB(2)(2a)1122212122,ABEHBAEGFA.1111111iiiikiikiiiiiABEHGAEGF1111,iiiiiiiiiiABEHBAEGFA.111kiiiiikiAGAB(ia)(i)ikiikiikiiHGFE1

27、11117.6.27.6.2 追趕法求解三對角方程追趕法求解三對角方程111111121nknnknnknnHGFE221111112121nnnnknnknnnnnABEHGAEGF1111111knnnnnknAGAB2211121211,nnnnnnnnnnABEHBAEGFA(n-1)a)(n-1)7.6.27.6.2 追趕法求解三對角方程追趕法求解三對角方程nknnknnHFE11111111nnnnknnnnnABEHAEGF11111nnnnnnnnknAEGFABEH(na)(n)7.6.27.6.2 追趕法求解三對角方程追趕法求解三對角方程 解出解出 后，利用式（后，利用式（

28、(n-1)a）,，(ia),，(2a), (1a)往前遞推，依次即可解出往前遞推，依次即可解出 ，于，于是得到是得到k+1時刻各節(jié)點(diǎn)的時刻各節(jié)點(diǎn)的 值（值（i = 1, 2, ，n）。）。 若若 均為已知（第一類邊界條件），則可直均為已知（第一類邊界條件），則可直接從式（接從式（2）開始求解，依次遞推至式（）開始求解，依次遞推至式（n - 1）即可。）即可。1kn11121211 , , , , kkknkn111 , knk7.6.27.6.2 追趕法求解三對角方程追趕法求解三對角方程(1) 輸入輸入D、K含水量含水量SITA關(guān)系（建議采用關(guān)系（建議采用函數(shù)函數(shù)子程序調(diào)用）；子程序調(diào)用）；節(jié)

29、點(diǎn)剖分總數(shù)節(jié)點(diǎn)剖分總數(shù)n及對應(yīng)的空間坐標(biāo)及對應(yīng)的空間坐標(biāo)z(1:n)（若為等步長，則（若為等步長，則輸入空間步長輸入空間步長DZ）；）；初始含水量剖面分布初始含水量剖面分布SITA0（1:n）；）；上、下邊界條件上、下邊界條件UBC、LBC（=1，2，or3）及相應(yīng)的邊界）及相應(yīng)的邊界值值QSUBC、QSLBCT（建議采用子程序調(diào)用）；（建議采用子程序調(diào)用）；初始時間步長初始時間步長DT0、最大時間步長、最大時間步長DTMAX，計算結(jié)束時，計算結(jié)束時間間TEND，迭代控制標(biāo)準(zhǔn)，迭代控制標(biāo)準(zhǔn)DET。(2) 賦初值賦初值T(0)=0 初始時間初始時間T(1)=DT0 第一次計算時間第一次計算時間*（以下賦含水量計算中間值）（以下賦含水量計算中間值） do i=1, nSITAM(i)=SITA0(i) enddo時段計數(shù)時段計數(shù)KT=1(3) 判斷計算時間是否已超過需模擬時間判斷計算時間是否已超過需模擬時間if (T(KT)=DTMAX) then DT=DTMAX T(KT)=T(KT-1)+DT endif(5) 開始計算系數(shù)矩陣，求解線性代數(shù)方程組開始計算系數(shù)矩陣，求解線性代數(shù)方程組 賦迭代計次數(shù)賦迭代計次數(shù) MP1=0(6) if (MP1=5)then (else goto (12)(7) 計算本時