土壤水運動數(shù)值模擬

土壤水運動數(shù)值模擬

土壤水運動是多孔介質(zhì)流運動的重要形式。土壤和地下水的各種文件過程，如大氣降水、地表水入滲、地表徑流、采礦蒸發(fā)和植物蒸發(fā)，與土壤水運動密切相關(guān)。對土壤水運動的數(shù)值模擬研究在大氣科學(xué)、土壤學(xué)、農(nóng)業(yè)工程、環(huán)境工程和地下水動力學(xué)等方面均具有重要意義。本文通過有限元法建立一維土壤水的數(shù)值模擬模型,針對土壤水運動數(shù)值模擬中對含水量的時間項比較敏感的問題,利用質(zhì)量集中法進行處理,并采用特殊的迭代格式以減少模擬過程中的水量平衡誤差。1土壤水價值1的模型1.1土壤含水率[]以垂直向上為正,一維土壤水的運動方程(Richard方程)可表達如下?θ?t=??z[Κ(h)(?h?z+1)]-S(1)?θ?t=??z[K(h)(?h?z+1)]?S(1)式中:θ為土壤體積含水率[L3L-3];h為土壤水負壓[L];S為根系吸水率或其它源匯項[T-1];K(h)為飽和/非飽和水力傳導(dǎo)度[LT-1]。本文模型采用ven-Genuchten的修正模型來描述三者之間的聯(lián)系。1.2有限元方程的求解用伽遼金有限元法對式(1)進行空間離散,則對于每一個計算節(jié)點∫Ω{?θ?t-??z[Κ?h?z+Κ]+S}?ndΩ=0(2)∫Ω{?θ?t???z[K?h?z+K]+S}?ndΩ=0(2)式中:n為計算節(jié)點編號;?n為線性插值基函數(shù);Ω為一維模擬空間。將式(2)展開并利用分部積分法,有∫Ω?θ?t?ndΩ+∫Ω(Κ?h?z+Κ)??n?zdΩ-(Κ?h?z+Κ)?n|ΖΤopΖBot+∫ΩS?ndΩ=0(3)∫Ω?θ?t?ndΩ+∫Ω(K?h?z+K)??n?zdΩ?(K?h?z+K)?n|ZTopZBot+∫ΩS?ndΩ=0(3)式中:ZTop和ZBot分別為模擬區(qū)域上、下邊界的垂向位置。將模擬區(qū)域內(nèi)的土壤水負壓h(z,t)用下式進行近似h′(z,t)=Ν∑n=1?n(z)hn(t)(4)h′(z,t)=∑n=1N?n(z)hn(t)(4)式中:hn(t)為待求系數(shù),其解為各離散節(jié)點處的土壤水負壓值;N為區(qū)域內(nèi)的計算節(jié)點總數(shù)。在式(3)中用h′(z,t)替換h(z,t),并將積分化在各計算單元上進行可得∑e∫Ωe?θ?t?ndΩ+∑e∫ΩeΚ?h′?z??n?zdΩ=Κ(?h′?z+1)?n|ΖΤopΖBot-∑e∫ΩeΚ??n?zdΩ-∑e∫ΩeS?ndΩ(5)∑e∫Ωe?θ?t?ndΩ+∑e∫ΩeK?h′?z??n?zdΩ=K(?h′?z+1)?n|ZTopZBot?∑e∫ΩeK??n?zdΩ?∑e∫ΩeS?ndΩ(5)對上式的時間項用質(zhì)量集中法進行處理,即:dθndt=∑e∫Ωe?θ?t?ndΩ/∑e∫Ωe?ndΩ(6)dθndt=∑e∫Ωe?θ?t?ndΩ/∑e∫Ωe?ndΩ(6)對計算區(qū)域內(nèi)的N個計算節(jié)點,應(yīng)用上式(5)和式(6)進行離散處理可以得到矩陣方程組[F]d{θ}dt+[A]{h}={Q}-{B}-{D}(7)式中:Fnm=δnm∑e∫Ωe?ndΩ;Anm=∑e∫ΩeΚ??m?z??n?zdΩ;Bn=∑e∫ΩeΚ??n?zdΩ;Dn=∑e∫ΩeS?ndΩ;Ωn為節(jié)點邊界流量項,Qn=Κ(?h′?z+1)?nΖΤopΖBot,對于非邊界節(jié)點,由于插值基函數(shù)?n在上下邊界處的值為0,因此Qn=0;n=m時δnm=1,n≠m時δnm=0。對于時間的離散,為保持計算的穩(wěn)定性,模型中采用隱式差分的格式,即:[F]{θ}j+1-{θ}jΔtj+[A]j+1{h}j+1={Q}j+1-{B}j+1-{D}j+1(8)式中:j+1代表當(dāng)前的時間層;j代表前一時間層;Δtj=tj+1-tj。模型的求解通過h進行,由于方程系數(shù)θ,A,B,D都是h的函數(shù),因此以上方程組具有高度的非線性,在每個Δt計算時段都需要通過迭代法進行求解。迭代過程對如何處理上述矩陣方程中的含水率項十分敏感。模型中采用了一種稱為“質(zhì)量守恒”的方法,以減少計算過程中的水量平衡誤差。這種方法在迭代過程中將式(8)中的第一項分成兩部分,即:[F]{θ}j+1-{θ}jΔtj=[F]{θ}k+1j+1-{θ}kj+1Δtj+[F]{θ}kj+1-{θ}jΔtj(9)式中:k+1表示當(dāng)前迭代;k表示上一次迭代。式(9)中右端第二項對于當(dāng)前迭代過程來說是已知的,將式(9)中右端第一項轉(zhuǎn)化為用水頭表示得[F]{θ}j+1-{θ}jΔtj=[F][C]kj+1{h}k+1j+1-{h}kj+1Δtj+[F]{θ}kj+1-{θ}jΔtj(10)式中:[C]、[F]為對角矩陣,其矩陣元素為Cnm=δnmCn,Cn為計算節(jié)點處的容水度。當(dāng)?shù)Y(jié)果滿足精度要求時,式(10)中右端第一項基本上等于0,這種特性可以有效減少求解過程中的水量平衡誤差。經(jīng)過上述處理之后,暫不考慮邊界條件的處理,迭代求解過程中矩陣方程可以表達如下([F][C]kj+1Δtj+[A]kj+1){h}k+1j+1=[F][C]kj+1Δtj{h}kj+1-[F]{θ}kj+1-{θ}jΔtj-{B}kj+1-{D}kj+1(11)當(dāng)前第k+1次迭代時,式(11)右端及左端的各項系數(shù)矩陣均已知,經(jīng)整理可得編程求解形式[G]{h}={g}(12)式中:G為系數(shù)矩陣;{g}為右端已知項。如果節(jié)點編號由上自下順序進行,則G為主對角線占優(yōu)的三對角矩陣,采用追趕法求解。1.3地表蒸發(fā)狀態(tài)土壤水方程的定解條件有兩個:一是初始條件h(z,0)=h0(z),二是邊界條件。即使是在一維情況下,土壤水的上邊界條件也具有相當(dāng)?shù)膹?fù)雜性。土壤水的上邊界主要接受降雨/灌溉,以及土表蒸發(fā),受限于土表的入滲能力,在降雨/灌溉過程中由于強度的大小不同,土表可能積水或不積水,也可能超出積水滯蓄深度產(chǎn)流。在地表處于蒸發(fā)狀態(tài)時,根據(jù)土壤輸水能力,也可分為達到或未達到蒸發(fā)極限的情況。目前大多數(shù)土壤水模型對于土表積水(SurfacePonding)過程以及產(chǎn)流過程尚無嚴(yán)格的數(shù)值方法,通常積水和非積水過程人為分離模擬,本文試圖從邊界處理的角度解決這一問題。1.3.1土壤水負壓h當(dāng)灌溉/降雨強度,或蒸發(fā)強度并未超出土壤入滲/蒸發(fā)能力時,上邊界條件為流量邊界條件(二類邊界條件)-Κ(?h/?z+1)=Ι(t)?hA<h<0(13)式中:hA為對應(yīng)土壤風(fēng)干含水率時的土壤水負壓,由土質(zhì)而定;I(t)為作用在土表之上的潛在凈流入/流出強度。由于模型取z軸向上為正,I(t)在數(shù)值上等于t時刻的外界蒸發(fā)強度減去降雨/灌溉強度。二類邊界條件可以直接應(yīng)用,在上邊界節(jié)點處,有Q=Κ(?h′?z+1)?ΖΤopΖΤopΖBot=-ˉΙΔtj,其中?ZTop為上邊界節(jié)點處的線性插值基函數(shù),ˉΙΔtj為計算時段Δtj內(nèi)的平均潛在凈流入/流出強度(下同)。作為己知項,-ˉΙΔtj可以并入式(12)的右端項[g]中。1.3.2tj邊界節(jié)點的流量強度計算當(dāng)降雨/灌溉強度超出土壤入滲能力時,地表將會出現(xiàn)積水現(xiàn)象,積水深度在沒有超出最大深度(對于田間灌溉可以理解為田埂高度,對于降雨產(chǎn)流則可以理解為地表洼地滯蓄深度)之前,超滲的水量將會在地表隨時間積累。對于一維問題,該邊界條件可用下式描述-?h/?t=Ι(t)+Κ(?h/?z+1)?0≤h≤hS(14)式中:hS為地表最大積水滯蓄深度。這種邊界情況下h既是地表的負壓水頭(此時為靜水壓力),又代表了地表的積水深度。由于此時邊界表達式中含有時間項,該邊界形式不同于一般的流量或水頭邊界,本文稱之為積水邊界條件。應(yīng)用該邊界條件時需要對時間進行離散。為保持與土壤水控制方程一致,離散時同樣采用隱式格式。對于計算時段Δtj,有Q=Κ(?h?z+1)?ΖΤop|ΖΤopΖBot=-[Ι(t)+?h?t)?ΖΤop]|ΖΤopΖBot=-ˉΙΔtj-hj+1ΖΤop-hjΖΤopΔtj(15)上式中含有未知項hj+1ΖΤop,即上邊界節(jié)點處當(dāng)前待求的水頭值,其系數(shù)1/Δtj可以并入式(12)的系數(shù)矩陣[G]的主對角線中,而已知項-(ˉΙΔtj+hjΖΤop/Δtj)則可以并入式(12)的右端項{g}中。式(15)具有較明顯的物理意義,右端第一項為Δtj時間內(nèi)上邊界作用的平均潛在凈流入/流出強度,第二項為Δtj時間內(nèi)上邊界積水深度從hjΖΤop變化到hj+1ΖΤop引起的積水量變化強度,兩者的綜合疊加結(jié)果即為上邊界通過的實際流量強度。當(dāng)Δtj時段計算收斂后,即已經(jīng)求出tj+1時刻各計算節(jié)點處的土壤水負壓{h}j+1,則可以通過上邊界處的節(jié)點方程確定上邊界的實際流量強度Q。如果在模型中n代表上邊界節(jié)點的編號,由式(8)可知Δtj時段內(nèi)上邊界節(jié)點處的實際流量強度為Qn=Ν∑m=1Fnmθj+1m-θjmΔtj+Ν∑m=1Aj+1nmhj+1m+Bj+1n+Dj+1n(16)其中:N為模擬區(qū)域內(nèi)的計算節(jié)點個數(shù)。式(15)適用于地表已經(jīng)積水時的情況,即在tj和tj+1時刻,上邊界節(jié)點處的土壤水負壓hZTop均大于0(靜水壓力)。但在邊界積水過程中有兩個特殊的時段需要單獨處理,一是上邊界從非積水狀態(tài)過渡到積水狀態(tài)時(hjΖΤop<0,hj+1ΖΤop≥0),二是從積水狀態(tài)過渡到非積水狀態(tài)時(hjΖΤop≥0,hj+1ΖΤop<0)。從邊界上的水量平衡出發(fā),結(jié)合式(15)的物理意義。對于第一種情況,在該Δtj時段內(nèi)上邊界實際的流量強度為Q=-ˉΙΔtj-(hj+1ΖΤop-0)/Δtj(17)上式右端第二項為Δtj時間內(nèi)邊界上的積水厚度從0變化到hj+1ΖΤop引起的積水量變化強度。與積水狀態(tài)時的處理相同,系數(shù)1Δtj并入式(12)的系數(shù)矩陣[G],而已知項-ˉΙΔtj則可以并入式(12)的右端項{g}。在該計算時段求解完畢后,邊界上實際通過的流量大小仍可由式(16)計算。對于第二種情況,從上邊界的水量平衡觀點出發(fā),此時邊界上的實際流量強度可以表達為Q=-ˉΙΔtj-(0-hjΖΤop)/Δtj(18)上式右端第二項為Δtj時間內(nèi)邊界上的積水厚度從hjΖΤop變化到0引起的積水量變化強度。注意到對于當(dāng)前時段Δtj,邊界節(jié)點在tj時刻的積水深度hjΖΤop已知,因此在該時段內(nèi)Q為已知的,將-ˉΙΔtj+hjΖΤop/Δtj并入式(12)中的右端項{g}即可。1.3.3建立邊界模型當(dāng)積水深度超出地表最大滯蓄深度hS時,積水深度不再增加,超滲的水量即形成地表徑流。此時的上邊界條件可以用水頭邊界條件(一類邊界)描述。h=hS,h>hS(19)為了模擬上邊界積水深度達到極限深度后不再增加的情況,在數(shù)值模擬過程中可如下處理:計算時段Δtj內(nèi)每次應(yīng)用上述的積水邊界條件進行計算時,迭代計算收斂之后都判斷一次積水深度hj+1ΖΤop的大小,如果超過極限深度hS,則將該時段的計算結(jié)果作廢,上邊界處的節(jié)點方程用h=hS代替,并對該時段進行重新計算。該時段計算完畢后上邊界的實際入滲強度Qn可由式(16)確定。根據(jù)上邊界處的水量平衡關(guān)系,可得Δt時段內(nèi)土表的產(chǎn)流量為RΔtj=-[-ˉΙΔtj-Qn]Δt。1.3.4土壤水負壓h在外界的蒸發(fā)力超過土壤的輸水能力時,土表處于風(fēng)干狀態(tài),土表的實際蒸發(fā)強度與外界蒸發(fā)力的大小將不再存在直接關(guān)系,此時的邊界條件為水頭邊界,可描述為h=hA?h≥hA(20)式中:hA為對應(yīng)土壤風(fēng)干含水率時的土壤水負壓。在滿足式(20)的判斷條件時,迭代過程中上邊界處的節(jié)點方程將由h=hA代替,上邊界的實際蒸發(fā)強度可由式(16)確定。2模型驗證2.1土壤水運動特性分析模型的正確性首先用試驗室一維入滲試驗驗證,二維有限元土壤水模型SWMS—2D曾利用該入滲試驗作過模型模擬測試。本文模型對該試驗進行同樣的數(shù)值模擬,并將模擬結(jié)果與SWMS—2D的模擬結(jié)果進行對比。該入滲試驗過程如下(圖1):一維砂柱的高度為61cm,初始土壤水負壓為-150cm;砂柱中的土質(zhì)均一,各向同性,飽和滲透系數(shù)為0.043cm/min;試驗過程中土柱頂端維持2cm厚的水層,底端封閉不排水;試驗持續(xù)的時間為90min。本次數(shù)值模擬時的計算節(jié)點總數(shù)為56個,計算單元55個,節(jié)點編號和單元編號的分布情況見圖1,所取土壤水運動參數(shù)見表1(與SWMS—2D所取參數(shù)一致)。模擬過程中上邊界應(yīng)用定水頭邊界(h=2cm),下邊界應(yīng)用隔水邊界(Q=0)。數(shù)值模擬過程中土柱剖面土壤含水率變化情況如圖2所示。從本模型90min時模擬的土壤含水率分布與SWMS—2D模擬的結(jié)果對比來看,兩模型之間模擬的結(jié)果十分接近,說明本文模型在基本計算方面是正確的。2.2砂柱剖面含水率結(jié)果分析以上用簡單的一維入滲模擬驗證了模型的基本正確性,為了驗證在降雨/灌溉、蒸發(fā)強度給定的情況下,上邊界對積水過程的模擬能力,維持以上模擬環(huán)境中的土壤水運動參數(shù)、初始條件和單元剖分情況不變,僅對上、下邊界的情況進行修改。模擬的總時長為480min,對于上邊界,假設(shè)在0～60min時間內(nèi)作用0.5cm/min的降雨強度,60～240min時間內(nèi)無降雨、蒸發(fā),240～480min時間內(nèi)作用0.0625cm/min的蒸發(fā)強度;對于下邊界,假設(shè)0～240min時間內(nèi)封閉不排水(Q=0),240～480min時間內(nèi)自由排水(排水期間下邊界節(jié)點處h=0)。實際情況下不可能出現(xiàn)如此大強度的降雨和蒸發(fā),這里只是作為檢驗?zāi)Ｐ退谩ＤM過程中暫不限制上邊界積水深度(hS取很大的值)。首先觀察上、下邊界節(jié)點在模擬過程中的變化情況。由圖3可知,在0～60min的降雨過程中,由于降雨強度很大,上邊界節(jié)點在計算過程中迅速飽和并開始積水,在降雨結(jié)束時刻(60min)邊界積水深度達到最大值21.27cm。隨后60～240min之內(nèi)上邊界無降雨、蒸發(fā),且下邊界無排水,然而在160min之前積水深度減少,這是因為砂柱內(nèi)部尚未完全飽和,積水可繼續(xù)入滲。160min之后砂柱完全飽和,積水深度維持13.33cm不變直到第240min。240～480min之內(nèi)上邊界開始有蒸發(fā)作用,同時下邊界開始排水,因此上邊界積水深度開始下降。在360min以后上邊界已無積水,節(jié)點處的土壤水負壓在高強度的蒸發(fā)條件下迅速增加。由圖4可知,由于在240min之前底邊界隔水,在上邊界的入滲影響尚未到達砂柱底端之前,下邊界節(jié)點處的土壤水負壓變化很小,115min之后入滲水量到達下邊界,該邊界節(jié)點的負壓迅速減小并向正值過渡,160min時砂柱完全飽和,下邊界節(jié)點處的土壤水負壓達到正的最大值74.33cm,并維持該值一直到240min之后下邊界開始自由排水,此時模型將下邊界切換為定水頭邊界(h=0)并維持到模擬結(jié)束。模擬過程中不同時刻的砂柱剖面含水率分布見圖5。從圖中分析可知邊界水頭的變化與剖面含水率之間存在良好對應(yīng)關(guān)系。模擬過程中上、下邊界處的累計流入/流出水量見圖6,以流入為正,流出為負。由圖5可知,對于上邊界,在160min之前上邊界持續(xù)有水量入滲,水量持續(xù)正累積。160min之后砂柱完全飽和,上邊界雖然有積水(圖3)但水量不能入滲,因此在160～240min之間累積水量維持16.7cm3無變化。240min之后下邊界開始排水,上邊界的積水又開始入滲。在第360min時積水在入滲和上邊界蒸發(fā)的雙重作用下消耗完畢,此時上邊界已經(jīng)正累積入滲22.4cm3的水量,此后上邊界在蒸發(fā)作用之下開始流量負累積直到模擬結(jié)束。對于下邊界,由于0～240min之內(nèi)為隔水邊界,因此累積流入/流出水量維持為0;240min之后邊界開始自由排水,邊界上持續(xù)有水量負累積。從以上模擬結(jié)果分析中可知,模擬過程中邊界處的土壤水負壓變化曲線、邊界的流入/流出水量累積曲線及各時刻砂柱剖面含水率分布之間均存在良好對應(yīng)關(guān)系,從一方面證明了模型的正確性,然而模擬過程中各個時刻的水量是否平衡、有無明顯計算誤差尚需分析確定。根據(jù)水量平衡關(guān)系,對于本次模擬應(yīng)該有:(1)在上邊界積水過程中任意時刻(0～360min):累計降雨量=累計蒸發(fā)量+上邊界累計入滲水量+上邊界積水量;(2)在數(shù)值模擬過程中任意時刻(0～480min):砂柱含水量=砂柱初始含水量+上邊界累積流入/流出水量+下邊界累積流入/流出水量。在模擬過程中選取不同的時刻進行分析,其結(jié)果見表2和表3。分析結(jié)果表明模擬計算過程中以上兩種水量平衡關(guān)系確實成立,而且計算誤差很小。以上積水過程驗證并未考慮積水滯蓄深度限