土壤水動力學(xué)SWD9土壤水運(yùn)動的數(shù)值解法與軟件應(yīng)用

1、土壤水動力學(xué) (Soil Water Dynamics)第9章 土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用 (Numerical simulation & application of Soil Water dynamics)毛 曉 敏 中國農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院Tel:mail: 土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用第9章 土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用n土壤水分運(yùn)動的數(shù)值計(jì)算方法n常用土壤水模擬軟件n應(yīng)用舉例土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用9.1 土壤水分運(yùn)動的數(shù)值計(jì)算方法n定解問題及常用數(shù)值計(jì)算方法nRichards方程的

2、差分離散n邊界條件處理n土壤水分運(yùn)動參數(shù)取值n差分方程的求解土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用9.1.1 定解問題及常用數(shù)值計(jì)算方法n定解問題：uRichards方程u邊界條件u初始條件n求解方法：u解析法、半解析法：u數(shù)值方法：土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用解析法解析法 ：求出解的函數(shù)表達(dá)式：求出解的函數(shù)表達(dá)式: 精度高，但精度高，但有很大局限性有很大局限性:簡化、特殊情況。數(shù)值法：求用數(shù)值表示的有限個離散點(diǎn)和離數(shù)值法：求用數(shù)值表示的有限個離散點(diǎn)和離散時段上的近似解散時段上的近似解: 具有廣泛適用性；容易具有廣泛適用性；容易程序化，數(shù)值模擬軟件多；隨著

3、計(jì)算機(jī)技術(shù)程序化，數(shù)值模擬軟件多；隨著計(jì)算機(jī)技術(shù)的提高和數(shù)值算法的改進(jìn)得以迅速發(fā)展。的提高和數(shù)值算法的改進(jìn)得以迅速發(fā)展。 土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用n常用數(shù)值計(jì)算方法：u全離散化方法：時間和空間進(jìn)行離散，利用微分方程推導(dǎo)出物理量在節(jié)點(diǎn)值的代數(shù)方程組，求解方程組得到微分方程在節(jié)點(diǎn)上的近似解t有限差分法(FDM)：用差商替代導(dǎo)數(shù) t控制容積法(CVM)：積分方程離散，守恒t有限元(FEM)：借助基函數(shù)，直接求函數(shù)的近似解u半離散方法：時間或空間離散t有限解析法（FAM）t線法(MOL)土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用有限差分法有限差分法利用Taylor

4、 級數(shù)用差商代替導(dǎo)數(shù)：用差商代替導(dǎo)數(shù)：前向差分：前向差分：后向差分：后向差分：中心差分：中心差分：土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用一階導(dǎo)數(shù)的近似表達(dá)式：一階導(dǎo)數(shù)的近似表達(dá)式： xffxxfxxfdxdfjiji, 1)()( 前差前差后差后差有限差分法的離散有限差分法的離散 ：i1 2 。 I j1 2 。 。 。 。 。 。 。 J f i,jxxffxxxfxfdxdfjiji, 1,)()( 中心差中心差xffxxxfxxfdxdfjiji,21,21)2()2(土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用？22dxfd若采用中心差的形式，二階若采用中心差的形

5、式，二階導(dǎo)數(shù)如何表示？導(dǎo)數(shù)如何表示？dxdfdxdD土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用有限元法用簡單的插值函數(shù)代替每個單元上未有限元法用簡單的插值函數(shù)代替每個單元上未知函數(shù)分布，形成一系列代數(shù)方程組的方法知函數(shù)分布，形成一系列代數(shù)方程組的方法 有限元法的離散有限元法的離散 ：有限元法變化形式較多，如：有限元法變化形式較多，如：1）從權(quán)函數(shù)的選擇來說，有配置法、矩量法、最小二乘法和伽遼金法； 2）從計(jì)算單元網(wǎng)格的形狀來劃分，有三角形網(wǎng)格、四邊形網(wǎng)格和多邊形網(wǎng)格； 3）從插值函數(shù)的精度來劃分，又分為線性插值函數(shù)和高次插值函數(shù)等。 f j土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與

6、軟件應(yīng)用Galerkin法引入基函數(shù)法引入基函數(shù)i 并采用并采用i作權(quán)函數(shù)，使得作權(quán)函數(shù)，使得： 舉例舉例Galerkin法：法：2222hxhySTht對于：對于：02222dthTSyhxh),()(),(),(1yxthtyxhtyxhnNjjj02222dxdythTSyhxhi土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用可用矩陣表示為：可用矩陣表示為：式中式中H、dH/dt分別為由分別為由h1，h2, hn和和dh1/dt, dh2/dt, dhn/dt組成的列矢量；組成的列矢量；D為導(dǎo)水矩陣，為導(dǎo)水矩陣，P為貯水矩陣。為貯水矩陣。若采用隱式方法：若采用隱式方法：通過求解大型

7、稀疏矩陣方程，得到通過求解大型稀疏矩陣方程，得到Hk+1FPDHdtdH FPDHtHHkkk 11BkAH 1土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用n數(shù)值計(jì)算的收斂性與穩(wěn)定性u收斂性：t、x 0時，差分方程解微分方程的解u穩(wěn)定性：計(jì)算誤差有界n數(shù)值彌散問題：u對流占優(yōu)情況下容易出現(xiàn)數(shù)值彌散，利用特征線法等來克服土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用9.1.2 Richards方程的差分離散n定解問題：u初始：u地表瞬時濕潤zKKzDDtkikikikikikikikikiki2/12/1212/112/111)()()(LzztzKzDztisi00)()(zKK

8、zDDtkikikikikikikikikiki2/12/1212/112/11)()()( 0 1 . k-1 k k+1 . M t 0 1 . i-1 i i+1 . N z zi zi+1 差分網(wǎng)格 tk+1 土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用n時間差分：u顯式差分格式：=0不穩(wěn)定u中心差分(Crank-Nicholson )格式：=0.5可能不穩(wěn)定u隱式差分格式：=1穩(wěn)定n參數(shù)取值：u時段初(=0)u時段中(=0.5)u時段末(=1)效果較好，需要迭代zKKzDDtkikikikikikikikikiki2/12/1212/112/11)()()(ikiikiiki

9、ihcba11111土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用9.1.3 邊界條件處理n一類邊界：0=sn二類邊界（通量邊界）：u直接差分u半節(jié)點(diǎn)水量平衡u虛擬節(jié)點(diǎn)土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用9.1.4 土壤水分運(yùn)動參數(shù)取值n時間：u顯式：用時段初含水率剖面計(jì)算參數(shù)u預(yù)報(bào)校正：以時段初值求解方程，用新的含水率校正參數(shù)，再次求解u迭代：n半節(jié)點(diǎn)土壤水分運(yùn)動參數(shù)：u算術(shù)平均u幾何平均u調(diào)和平均u土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用9.1.5 差分方程的求解n對于非顯式差分，得到方程：n代入邊界條件得到三對角方程組如下，用追趕法（Thomas算法）求

10、解ikiikiikiihcba11111 JkIkkJJdddbacbacb.211121122211()ALUAdL Ud土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用求解步驟求解步驟求解土壤水分運(yùn)動方程時，首先要給出有關(guān)土壤水分運(yùn)動參數(shù)的表達(dá)式、初始條件與邊界條件，然后按一定的時間、空間步長進(jìn)行離散。在每一離散時段內(nèi)，進(jìn)行以下迭代計(jì)算：(1) 以時段初的土壤含水率作為時段末含水率的預(yù)報(bào)值。(2) 計(jì)算相關(guān)的土壤水分運(yùn)動參數(shù)，利用追趕法求解土壤水分運(yùn)動的差分方程組（三對角方程組），得到時段末含水率。(3) 如果時段末含水率的預(yù)報(bào)值和計(jì)算值相差不大，可以認(rèn)為迭代收斂，完成本時段計(jì)算，進(jìn)入

11、下一時段繼續(xù)求解；否則，以時段末含水率的計(jì)算值作為新的預(yù)報(bào)值，返回步驟(2)繼續(xù)進(jìn)行迭代。在求解過程中，也可以根據(jù)各時段迭代情況對時間步長進(jìn)行調(diào)整。一般情況下，如果迭代次數(shù)較少，表明時段內(nèi)含水率變化不大，可以增大時間步長；否則，如果迭代次數(shù)較多，需要減小時間步長。這種自動調(diào)節(jié)步長的方法既可以保證模擬計(jì)算的精度，也可以有效地提高計(jì)算速度。土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用 不同差分處理方法的模擬結(jié)果比較(a) CN格式 (b) 隱式格式 (c) 特征差分格式 降水入滲過程：降水強(qiáng)度0.2427mm/min，歷時180min 土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用0

12、5101520253000.10.20.30.4含水率深度 z (cm)試驗(yàn)點(diǎn)t=1mint=5min05101520253000.10.20.30.4含水率深度 z (cm)試驗(yàn)點(diǎn)全隱式格式指數(shù)加權(quán)隱式格式 (a) 中心差分格式 (b) 全隱式及指數(shù)加權(quán)隱式格式 降雨入滲試驗(yàn)?zāi)M結(jié)果（降雨強(qiáng)度0.05374cm/min，歷時140min） 土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用9.2 常用土壤水模擬軟件n土壤水模擬軟件：一般基于飽和-非飽和土壤水分運(yùn)動方程，采用有限差或有限元進(jìn)行數(shù)值求解。n常用軟件：（一般不僅僅局限于土壤水分運(yùn)動模擬，有的也包括地下水運(yùn)動、溶質(zhì)運(yùn)移、熱遷移、甚

13、至作物生長等。)n一維軟件：SWAP、HYDRUS系列、SHAW 等n二維/三維軟件：土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用nSWAPn由瓦赫寧根大學(xué)集成的SWAP(SoilWaterAtmospherePlant) 模型是由詳盡完善的土壤水運(yùn)移、溶質(zhì)遷移、熱量傳輸、土壤蒸發(fā)、植物蒸騰和作物生長子模塊組成。 土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用nHYDRUS1D、HYDRUS2D/3D系列軟件是美國農(nóng)業(yè)部鹽土實(shí)驗(yàn)室開發(fā)的模擬飽和-非飽和介質(zhì)中水分、熱、溶質(zhì)運(yùn)移的模型。其中HYDRUS1D在該系列模型中主要進(jìn)行垂向一維條件下的模擬。模型采用水頭形式的Richards

14、方程，因此，可以對飽和-非飽和區(qū)進(jìn)行統(tǒng)一模擬，可以處理諸如有壓入滲、地下水位線波動情況下的區(qū)域水流運(yùn)動問題。土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用nSHAW 是由美國農(nóng)業(yè)部Flerchinger 等人建立的，該模式是一個最初發(fā)展模擬土壤凍結(jié)和融化的一維模式 ， 通過對土壤凍結(jié)、融雪、土壤溫度、水量、蒸騰、水平衡影響的處理， 可以準(zhǔn)確預(yù)測土壤的凍結(jié)融化過程土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用二維/三維軟件n 1) 專注于地表產(chǎn)匯流過程的模型，如SWAT；n2）專注于干旱區(qū)區(qū)域水均衡的模型，如清華大學(xué)開發(fā)的干旱區(qū)散耗型水文模型；n3）專注于區(qū)域地下水運(yùn)動的模型，如FE

15、FLOW，HYDRUS2D/3D，MODFLOW，SUTRA等。土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用地表產(chǎn)匯流模型nSWAT模型是近年來發(fā)展起來的分布式水文模型，可用于模擬預(yù)測各種管理措施及氣候變化對水資源供給的影響，評價(jià)流域非點(diǎn)源污染現(xiàn)象。模型陸面模擬組成包括水文、氣候、侵蝕、土壤溫度、植物生長、營養(yǎng)物質(zhì)、殺蟲劑、以及土地利用和水資源管理。土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用干旱區(qū)散耗型水文模型n由于描述流域水文循環(huán)的分布式水文模型多考慮了產(chǎn)匯流過程，并把產(chǎn)匯流模擬結(jié)果作為檢驗(yàn)?zāi)M準(zhǔn)確性的重要標(biāo)準(zhǔn)。這與干旱區(qū)綠洲的情形是很不相同的。綠洲水文循環(huán)的特點(diǎn)是降水稀少

16、、蒸發(fā)強(qiáng)烈，水資源多依靠上游河渠輸水并經(jīng)過各級渠道分配到田間，最終消耗于田間蒸散發(fā)。對于這種情況，清華大學(xué)開發(fā)了相應(yīng)的干旱區(qū)散耗型水文模型。土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用區(qū)域地下水運(yùn)動模型nMODFLOW是由美國地質(zhì)調(diào)查局 (USGS) 開發(fā)的孔隙介質(zhì)中三維地下水流數(shù)值模擬的模型，其源代碼、執(zhí)行文件和使用說明書都可以在其官方網(wǎng)站上免費(fèi)下載。MODFLOW自從開發(fā)以來，已經(jīng)經(jīng)過數(shù)次更新，從起初的MODFLOW88更新為MODFLOW96、MODFLOW2000、MODFLOW2005等。以MODFLO為基礎(chǔ)與地表水、溶質(zhì)運(yùn)移軟件等耦合產(chǎn)生一些新的模擬軟件，MODFLOW也被

17、嵌入一些用戶界面中（VP，GMS，PMwin etc.)，使其具有方便的前后處理界面。土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用nFEFLOW 是由德國WASY水資源規(guī)劃系統(tǒng)研究所研制開發(fā)的地下水模型軟件包，可以進(jìn)行復(fù)雜二維和三維飽和、非飽和帶水流（包括密度流）、熱量運(yùn)移和溶質(zhì)對流彌散的模擬，同時可以模擬溶質(zhì)遷移過程中發(fā)生的線性非線性吸附反應(yīng)和一階動力學(xué)化學(xué)反應(yīng)等。土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用nSUTRA由USGS開發(fā)，采用有限元法進(jìn)行二維/三維、飽和/非飽和、變密度情況下地下水流動和溶質(zhì)/能量運(yùn)移的模擬。我國已有報(bào)道介紹該軟件，但應(yīng)用尚不多見。國外多采用該模

18、型進(jìn)行海水入侵方面的模擬研究。 土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用土壤水模擬軟件HYDRUS nHYDRUS1D、HYDRUS2D/3D系列軟件 n模型采用水頭形式的Richards方程，可以對飽和-非飽和區(qū)進(jìn)行統(tǒng)一模擬，可以處理諸如有壓入滲、地下水位線波動情況下的區(qū)域水流運(yùn)動問題。n在給定充分供水條件下的地面蒸發(fā)和作物蒸騰量條件下，該模型可以模擬土壤水鹽脅迫下的實(shí)際根系吸水（作物蒸騰）和地面實(shí)際蒸發(fā)過程。n可以處理多種靈活的邊界條件（大氣邊界、自由排水、定水位、動水位等），具有良好的用戶界面和靈活的輸入輸出功能。模型采用有限元方法進(jìn)行數(shù)值離散求解。土壤水動力學(xué) 第9章土壤水動力學(xué)的數(shù)值求解與軟件應(yīng)用9.3 HYDRUS應(yīng)用舉例 （1） 應(yīng)用軟件模擬前要考慮：應(yīng)用軟件模擬前要考慮：n計(jì)算區(qū)域和計(jì)算時段？n土壤水分運(yùn)動參數(shù)？n如何給定初始條件？u設(shè)起初為零通量穩(wěn)定剖面，下邊界含水率為0.3，各處