利用FEFLOW求解非飽和帶土壤水運(yùn)動問題.doc

精品論文推薦利用 feflow 求解非飽和帶土壤水運(yùn)動問題曹英杰1 河海大學(xué) 水文水資源學(xué)院e-mail：yingjie_摘要：土壤水是地表水和地下水聯(lián)系的紐帶，在水循環(huán)過程中起著重要作用，地下水?dāng)?shù)值模擬中的許多問題都與土壤，即包氣帶密切相關(guān)，如地下水?dāng)?shù)值模擬中的補(bǔ)給和 蒸發(fā)問題。feflow 是德國 wasy 水資源規(guī)劃和系統(tǒng)研究所基于有限元開發(fā)的地下水 數(shù)值模擬軟件，功能強(qiáng)大齊全。本文用其對一長 1m，寬 0.01m 土柱中的土壤水進(jìn)行數(shù) 值模擬，得到了基膜勢隨深度變化的曲線。模擬中利用 van genuchten-mualem model 計(jì)算非飽和導(dǎo)水率，由于資料有限，沒有能夠?qū)?shù)值模擬的結(jié)果與解析解進(jìn)行對比，驗(yàn) 證計(jì)算精度，這方面的工作需要進(jìn)一步研究探討。關(guān)鍵詞：feflow；van genuchten-mualem model；土壤水；數(shù)值模擬1 概 述土壤是一種由固、液、汽三相組成的多孔系統(tǒng)。固相基質(zhì)由大小、形狀和排列不同的土 粒組成，土粒的排列和組合，決定了土壤孔隙的特征、水及空氣在土壤孔隙中傳輸和存在的 形式。當(dāng)土壤孔隙全部被水充滿，則土壤中水處于飽和狀態(tài)，存在于地面以下的飽和土壤水 即地下水。土壤孔隙未被水充滿，土壤水分處于非飽和狀態(tài)，一般稱為土壤水。由此可見，土壤水和地下水、地表水有著密切的聯(lián)系，地表水入滲進(jìn)入土壤，可以補(bǔ)給 土壤水，土壤水進(jìn)一步下滲，進(jìn)入地下飽和帶以下，又補(bǔ)給地下水；同時，地下水通過毛細(xì) 作用，不斷向非飽和帶的土壤中輸運(yùn)地下水，補(bǔ)充土壤水。在一定條件下地表水，土壤水和 地下水可以相互轉(zhuǎn)化。土壤水是地表水和地下水聯(lián)系的紐帶，在水循環(huán)過程中起著重要作用， 它和農(nóng)業(yè)、水文、環(huán)境等領(lǐng)域有著密切關(guān)系。2 feflow 簡介feflow是德國wasy水資源規(guī)劃和系統(tǒng)研究所于 20 世紀(jì) 70 年代末開發(fā)的基于有限元 法的地下水?dāng)?shù)值模擬軟件，它是迄今為止功能最為齊全的地下水模擬軟件包之一。該軟件包 具有圖形人機(jī)對話、地理信息系統(tǒng)數(shù)據(jù)接口、自動產(chǎn)生空間各種有限單元網(wǎng)格、空間參數(shù)區(qū) 域化及快速精確的數(shù)值算法和先進(jìn)的圖形視覺化技術(shù)等特點(diǎn)。在feflow系統(tǒng)中，用戶可以 很方便迅速地產(chǎn)生空間有限單元網(wǎng)格，設(shè)置模型的參數(shù)和定義邊界條件，運(yùn)行數(shù)值模擬以及 實(shí)時圖形顯示結(jié)果與成圖3。feflow的應(yīng)用領(lǐng)域非常廣泛，例如：水量模擬可以模擬水源地開采或者油田注水 對區(qū)域地下水流場的影響、模擬水庫放水或者河流斷流時，河道沿線地下水流場的變化以及 非飽和帶土壤水運(yùn)動問題等；水質(zhì)模擬模擬污染物在地下水中的遷移過程及其時間空間 分布模式、模擬沿海地區(qū)抽取地下水引起的海水入侵等；溫度模擬模擬非飽和帶以及飽 和帶溫度場的分布3。- 5 -3 土壤水運(yùn)動的基本方程和求解方法土壤水分運(yùn)動遵從達(dá)西定律。在達(dá)西定律和質(zhì)量守恒的基礎(chǔ)上，richard 于 1931 年利用 數(shù)學(xué)物理方法導(dǎo)出了非飽和流基本方程： = k() t式中： 為土壤含水率； 為土水勢； k ( ) 為非飽和導(dǎo)水率。（3.1）基本方程的一維、二維、三維以及在柱坐標(biāo)、球坐標(biāo)的表達(dá)式均可通過一定的途徑導(dǎo)出。 由于土壤水多垂向運(yùn)動，故忽略土壤水水平運(yùn)動，得到土壤水運(yùn)動的一維 richard 方程： = k() （3.2）tz z 式中符號與（3.1）式相同。（3.2）式配合相應(yīng)的初始條件和邊界條件即可對土壤水的垂向分布進(jìn)行求解。實(shí)際應(yīng) 用中多采用數(shù)值法進(jìn)行求解，常用的數(shù)值求解方法為有限差分法和有限單元法。有限差分法(finite difference method)就是以差商近似代替微商，將土壤水分運(yùn)動的偏微 分方程變?yōu)椴罘址匠?，組成可以直接求解的代數(shù)方程組。有限單元法(finite element method)即是用簡單的插值函數(shù)來代替每個單元上的未知函 數(shù)分布，然后，集合起來形成可以直接求解的代數(shù)方程組。由于建立代數(shù)方程組的出發(fā)點(diǎn)不 同，有限單元法又有各種變化：如將求解土壤水分運(yùn)動的偏微分方程轉(zhuǎn)化為等價的泛函求極 值問題里茲(ritz)法；建立在余項(xiàng)的加權(quán)和趨于零的基礎(chǔ)上的伽遼金(galerkin)法；以及 根據(jù)一般的水均衡概念建立方程組的均衡法等。大量的計(jì)算實(shí)踐證明，在各種求解土壤水分運(yùn)動問題的數(shù)值方法中，都不同程度地包含 著數(shù)值彌散和數(shù)值振蕩這兩類計(jì)算誤差，并且這兩類誤差之間還存在密切的聯(lián)系。為了克服 上述數(shù)值困難，許多研究者對這兩種方法進(jìn)行了改進(jìn)。改進(jìn)方法有如下兩種： 有限元集中質(zhì)量法，是對有限元通常離散方法的一種改進(jìn)，就是將質(zhì)量矩陣每行的非對角元集中到主對角線上做代數(shù)和，這些矩陣的元素依賴于負(fù)壓水頭，因而必須估計(jì)在半時間 層上的水頭分布。采取這種方法，可以避免振蕩非物理解的產(chǎn)生，并且也可以較好地處理邊 界條件2。有限體積法，又稱為控制體積法。其基本思路：將計(jì)算區(qū)域劃分為一系列不重復(fù)的控制 體積，并使每個網(wǎng)格點(diǎn)周圍有一個控制體積;將待解的微分方程對每一個控制體積積分，便 得出一組離散方程。其中的未知數(shù)是網(wǎng)格點(diǎn)上的因變量的數(shù)值。為了求出控制體積的積分， 必須假定值在網(wǎng)格點(diǎn)之間的變化規(guī)律，即假設(shè)值的分段的分布剖面。簡言之，子區(qū)域法加離 散就是有限體積發(fā)的基本方法2。4 非飽和導(dǎo)水率的確定與地下水?dāng)?shù)值模擬相比，非飽和帶的土壤水?dāng)?shù)值模擬更為復(fù)雜，原因在于土壤水的非飽 和導(dǎo)水率隨土壤的含水率會發(fā)生顯著的變化，在脫濕和吸濕的過程中呈復(fù)雜的繩套曲線關(guān) 系。因此，往往需要大量的實(shí)驗(yàn)來獲得非飽和導(dǎo)水率的值及土壤水分特征曲線。本文采用經(jīng)典的 van genuchten-mualem model 來確定非飽和導(dǎo)水率和飽和率及土水勢（如果簡化問題，可以忽略溶質(zhì)勢和重力項(xiàng)，僅僅考慮基膜勢）。van genuchten-mualem model的形式如下1：1 0sf = 1 + (a )n m（4.1）e 1 0kr = (se )1 1 (s )1/ m2（4.2）4.2 式中：f 1/ 2 fesf sff fsf = r = r esf sff f（4.3）s r s r（4.1）至（4.3）中：a、n、m為水分特征曲線的適線參數(shù)；kr為非飽和導(dǎo)水率；sef為土壤相對飽和率； 為土壤含水率；其它符號同前。5 非飽和土柱含水量數(shù)值模擬以 van genuchten-mualem model 為基礎(chǔ)，利用 feflow 對如下問題進(jìn)行求解： = k() 1m z 0m tz z t = 0 = 1.0m = 0.75m z = 0 z =1 = 10.0m（5.1）利用 van genuchten-mualem model 確定非飽和含水率，模型中的各個適線參數(shù)選取如 下：n=2，m 取 n 的 1/4，即 m=0.5，a=3.351/m。r其他土壤水動力參數(shù)取值如下：孔隙度 =0.368，殘留飽和率 sf =0.277，飽和含水率ssf =1.0，飽和水力傳導(dǎo)度（滲透系數(shù)）k=0.922 10-4m/s。設(shè)土柱高度 1m，寬 0.01m，計(jì)算時將。模擬期為一天。本次模擬采用feflow中的矩形有限元計(jì)算，垂向剖分為高 0.005m，寬 0.01m的網(wǎng)格，網(wǎng)格剖分?jǐn)?shù) 200。時間步長方案選用ab/tr自動迭代方法，初始時間步長為 1 10-5天。計(jì)算結(jié)果如表 1 和圖 1 所示：表1 不同深度土壤中基膜勢計(jì)算值深度(m)1.000.900.800.700.60基膜勢(kpa)-17.65-23.32-104.08-105.42-104.44深度(m)0.500.400.300.200.10基膜勢(kpa)-103.46-102.48-101.50-100.52-99.540.00-20.00基膜勢（kpa）-40.00-60.00-80.00-100.00-120.000.00 0.20 0.40 0.60 0.801.00深度（m）圖1 不同深度上土壤中水分基膜勢計(jì)算結(jié)果6 結(jié)論及建議本文采用 van genuchten-mualem model 對非飽和導(dǎo)水率進(jìn)行計(jì)算，以此為基礎(chǔ)利用矩形 有限元法對長 1m 土柱中土壤水的基膜勢進(jìn)行了數(shù)值模擬，得到了定解條件下土柱中土壤水 基膜勢隨深度的變化情況。為進(jìn)一步研究土壤水變化情況提供了基礎(chǔ)。參考文獻(xiàn)1 o. ippisch, h.-j. vogel， p. bastian validity limits for the van genuchtenmualem model andimplications for parameter estimation and numerical simulationjadvances in water resources，2006，volume 29，issue 12： pages 1780-17892 候憲東，汪志榮 非飽和土壤水分運(yùn)動數(shù)值模擬研究綜述j 水資源與水工程學(xué)報(bào)，2006 年，17 卷4 期：41-45，493 賀國平，邵景力，崔亞莉等 feflow 在地下水流模擬方面的應(yīng)用j 成都理工大學(xué)學(xué)報(bào)（自然科 學(xué)版），2003 年，30 卷 4 期：356-361numerical simulation of soil water in unsaturated zone byfeflowcao yingjie11 college of hydrology and water resource hohai university, nanjing, prcabstractsoil water is conjunction of groundwater and surface water, and it plays a very important role in the hydrology cycle. there are many problem correlated with the movement of soil water in groundwater numerical simulation, for example, how to deal with the recharge and evaportanspiration more exactly. feflow is developed by wasy of german. it is a famous groundwater numerical simulation software and has powerful abilities to simulate lots of groundwater flow and transport problems. in this paper, numerical simulation of soil watermov