第7章土壤水分運動的數值解

1、7.1 7.1 概述概述7.2 7.2 差分的概念及分類差分的概念及分類7.3 7.3 顯式差分、隱式差分和中心差分顯式差分、隱式差分和中心差分7.4 7.4 土壤水分運移方程的線性化方法土壤水分運移方程的線性化方法7.5 7.5 線性化方法與土壤水分參數的取值線性化方法與土壤水分參數的取值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.薛禹群，謝春紅，水文地質學的數值法，煤炭工業(yè)出版社，薛禹群，謝春紅，水文地質學的數值法，煤炭工業(yè)出版社，1980.土壤水分運動方程的求解方法：土壤水分運動方程的求解方法： 解析解（解析解（Analytical solution） 數值解（數值解（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é)點（個節(jié)點（i = 1, 2, 3, ,

5、 n），），z (1) = 0, z (n) = L，已，已知知t = tk時刻各節(jié)點含水率時刻各節(jié)點含水率 的分布，求的分布，求t = t k+1時刻各節(jié)點含水率時刻各節(jié)點含水率 的分布，的分布， 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、式是有條件收斂的，一般應滿足，一般應滿足max2Dzrt其中：其中：r多取多取1/2。由于土壤接近飽和時，。由于土壤接近飽和時，Dmax很很大，故大，故 t一般要求取值很小，耗費機時。一般要求取值很小，耗費機時。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 土壤水分運移方程的線性化方法土壤水分運移方程的線性化方法 以基質勢水頭以基質勢水頭h為因變量的一維垂

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

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

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

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

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

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

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

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

18、值），再求解隱式差分方程，得到時數值（預報值），再求解隱式差分方程，得到時段末的含水率的值（校正值）。段末的含水率的值（校正值）。計算時，先假定本時段末（計算時，先假定本時段末（tk+1）的）的 m或或 （一般做（一般做法是先取時段初的函數值法是先取時段初的函數值 or ） ， or 等參數的迭代初值等參數的迭代初值 求解線性代數方程組得求解線性代數方程組得 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 的取值可根據實際問題的精度要求而定。的取值可根據實際問題的精度要求而定。 此法計算稍復雜，但精度較高，數值計算時此法計算稍復雜，但精度較高，數值計算時經常采用。為避免迭代次數太多，一般可設置一經常采用。為避免迭代次數太多，一般可設置一個判斷語句：當迭代次數個判斷語句：當迭代次數m5時，將時間步長時，將時間步長 減半，再重新計算。減半，再重新計算。kkkttt117.5.17.5.1線性化方法線性化方法7.5.1.3迭代法7.5 線性化方法與土壤水分參數的取值線性化方法與土壤水分參數的取值7.5.2.1直接法直接法先求出兩結點的參數，然后先求出兩結點的參數，然后(1)取兩結點參數的算術平均值

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

22、2iiDD(3)幾何平均值：幾何平均值：112i ii1122iiDD7.5 線性化方法與土壤水分參數的取值線性化方法與土壤水分參數的取值7.5.2 3其他方法其他方法（1）求出的參數和兩結點處的參數和，三點參）求出的參數和兩結點處的參數和，三點參數取平均，故可稱為數取平均，故可稱為“三點式三點式”：11122124iiiiDDDD（2） 11121iiiiiDDd （3）由鄰近結點的參數代替，如）由鄰近結點的參數代替，如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、邊界條件的處理邊界條件的處理 驗驗 證：證： 當上、下邊界為第一類邊界時當上、下邊界為第一類邊界時)()()(),()(), 0(121111121kknkkztftftftLtft)(, 0, 1)2.(Eqn)(, 0, 1) 1.(Eqn12n11111knnktfHEFtfHGF7.6.17.6.1邊界條件的處理邊界條件的處理1121111HGFkknknnknnHFE111 當上、下邊界為第二、三類邊界時當上、下邊界為第二、三類邊界時)4()()3(), 0(210tqKzDttqKzDzLzz)(),(122/111112111112312111223knnnknknnkkkktq

25、KzzDtqKzzD1121111HGFkknknnknnHFE111211212123111111112231)(,),(,nknnnnnnnkkKtqHEFzzDEKtqHFGzzDFnknnknnnknnknnknnkkkkkkkkHFEHGFEHGFEHGFEHGF111111111121314313312321321221121121111經線性化后，原定解問題可化為經線性化后，原定解問題可化為n n階線性代數方程組，即：階線性代數方程組，即：以矩陣的形式表示即為：以矩陣的形式表示即為：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é)點的時刻各節(jié)點的 值（值（i = 1, 2, ，n）。）。 若若 均為已知（第一類邊界條件），則可直均為已知（第一類邊界條件），則可直接從式（接從式（2）開始求解，依次遞推至式（）開始求解，依次遞推至式（n - 1）即可。）即可。1kn11121211 , , , , kkknkn111 , knk7.6.27.6.2 追趕法求解三對角方程追趕法求解三對角方程(1) 輸入輸入D、K含水量含水量SITA關系（建議采用關系（建議采用函數函數子程序調用）；子程序調用）；節(jié)

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