王 芳，楊虎山

（忻州師范學(xué)院數(shù)學(xué)系，山西忻州 034000）

復(fù)雜邊界非飽和模型有限元數(shù)值模擬

王 芳，楊虎山

（忻州師范學(xué)院數(shù)學(xué)系，山西忻州 034000）

利用有限元集中質(zhì)量法對一維復(fù)雜邊界土壤水分運動數(shù)學(xué)模型（非飽和Richards模型）進(jìn)行了數(shù)值模擬，并進(jìn)行了模型驗證．模型驗證結(jié)果表明，該方法計算精度很高，符合實際工作要求．

土壤水分運動；有限元；數(shù)值模擬

非飽和土壤水分運動模型可用一類非線性的對流擴(kuò)散方程表示．研究此類問題，對干旱半干旱地區(qū)的土壤水分預(yù)測尤為重要，對發(fā)展節(jié)水型農(nóng)業(yè)有重要意義．

非飽和土壤水分運動問題可歸結(jié)為下面的模型方程：

求h（z ,t ）使得對于任意的T＞0，滿足：

其中，λ為指標(biāo)變量；h（L）為負(fù)壓水頭（cm）；C（h）為比水容量（1cm）；K（h）為非飽和導(dǎo)水率（cmm in）；Z為空間坐標(biāo)（cm）；h0（z ）為土壤的初始負(fù)壓水頭（cm）；為變量，是隨時間t變化的水頭值（cm），在實際計算時，采用插值函數(shù)進(jìn)行處理；L和T分別是空間和時間變化的上界．

對有關(guān)滲透問題的一維非飽和流問題的研究[1-3]，用有限差分法[4]對方程進(jìn)行離散時，結(jié)果對邊界條件以及土壤參數(shù)極為敏感．本文通過變分的辦法將其直接轉(zhuǎn)化為已知通量的計算，較好地解決了這類邊界條件問題．

1 數(shù)學(xué)模型的建立

1.1 基本假定

1）各層土壤均質(zhì)且各向同性；

2）入滲水流為連續(xù)介質(zhì)且不可壓縮，在土壤水分運動過程中，土壤骨架不變形．

1.2 方程的離散

方程（1）中，由于比水容量C（h）和非飽和導(dǎo)水率K（h）為基質(zhì)勢h的函數(shù)，所以方程（1）屬于二階非線性偏微分方程，除少量問題外，一般情況下求解析解是困難的，須使用數(shù)值求解．下面考慮方程（1）在初始條件（2）及上、下邊界條件（3）、（4）的定解問題，以此說明本模型的數(shù)值模擬過程．

這里的（·,·）表示區(qū)域Ω上的L2內(nèi)積，以下類同．

1.2.1 半離散有限元逼近

先在空間方向上用有限元逼近[6]，對Ω=（0,L）進(jìn)行有限元剖分：

其中，

方程（6）定義了一類非線性常微分方程組，要數(shù)值求解，仍須對時間變量離散化．

1.2.2 全離散格式

在時間方向上引入有限差分格式逼近矩陣方程的時間導(dǎo)數(shù)．

對時間t進(jìn)行離散，設(shè)時間步長為Δt，tm=mΔt ．在時間層m與m+1之間，tm≤t≤tm+1，將變量t變換到變量τ，t=（m +τ）Δt，其中0≤τ≤1．

在tm≤t≤tm+1，對Xj（t）進(jìn)行線性插值，令

其中，{Xm}=（X1（mΔt）,X2（mΔt）,…,Xn（mΔt ））T，{Xm+1}也類似．

同理也可對{F}進(jìn)行插值．這樣方程（6）可化為：

若（7）兩邊乘上一個非負(fù)的權(quán)函數(shù)，記為w，再令τ從0到1積分，就得到：

2 數(shù)值模擬

用本文所建立的有限元集中質(zhì)量法的數(shù)值模型對山西省太谷縣三臺村的擾動黃土土樣[7]進(jìn)行數(shù)值模擬．解決非飽和土壤水分運動數(shù)值模擬的求解問題，土壤非飽和導(dǎo)水率、水分?jǐn)U散率、土壤水分特征曲線3個土壤水分運動參數(shù)是必不可少的重要數(shù)據(jù)．其中，土壤水分特征曲線是定量描述和預(yù)測土壤水分運動的重要基礎(chǔ)．一般地，水分特征曲線的經(jīng)驗?zāi)Ｐ涂煞譃閮绾瘮?shù)、指數(shù)函數(shù)、雙曲線余弦函數(shù)和誤差函數(shù)4類[8]，其中，冪函數(shù)模型是描述土壤水分特征的最普遍模型[8]．本文采用較常用的指數(shù)形式：

式中，s為土壤水吸力（cm）．非飽和導(dǎo)水率的測定有直接法和間接法兩種[4]．由于直接測定法耗時多、昂貴、繁瑣，且測定范圍較窄，可獲得的實驗數(shù)據(jù)不能代表土壤水分特征的完整關(guān)系，因此許，多研究提出了間接推求導(dǎo)水率的方法[4]．本文采用較常用的指數(shù)形式，推求得土壤非飽和導(dǎo)水率K與負(fù)壓水頭h的關(guān)系為：

利用土壤水分特征曲線，推求得土壤比水容量C與負(fù)壓水頭h的關(guān)系為：

C（h）=3.42×10-6e9.3865（8.045-ln（-h））-9.3865

土壤水分運動的初始和邊界條件為：

將區(qū)域Ω=（0,100），每個單元的長度為1cm，時間步長Δt=5s，模擬1小時，將這些數(shù)據(jù)輸入該程序，得到入滲率分布圖，如圖1所示．

圖1 垂直入滲率與時間的關(guān)系曲線

從圖1可知，本文所得到數(shù)值模擬結(jié)果與文獻(xiàn)[7]提到的試驗結(jié)果是相吻合的．從圖中還可以看出入滲率的變化規(guī)律：從入滲開始到0.1h這一段時間，入滲率從0.2迅速下降到0.04以下，此后，地表層z = 0處的入滲率變化不大，然后地表逐步接近飽和含水率．

3 結(jié)語

本文利用非飽和流方程，建立了變水頭非恒定水位復(fù)雜邊界條件下一維變水頭均質(zhì)水入滲及土壤水分運動的數(shù)學(xué)模型．對于入滲邊界條件，通過變分的辦法將其直接轉(zhuǎn)化為已知通量的計算，從而較好地處理了邊界條件問題．在每一個新的時間層，利用插值函數(shù)近似代替通過試驗實測的負(fù)壓水頭．?dāng)?shù)值模擬結(jié)果表明，采用非飽和流方程可以很好地模擬均質(zhì)土壤一維變水頭的土壤水分運動，該方法對田間土壤水分的預(yù)測預(yù)報有重要意義．

[1] Haverkamp R, Vauclin M, Touma J, et al. A comparsion of numerical simulation models for one-dimensional infiltration [J]. Soil Science Society of America Journal, 1977, 41：285-294.

[2] Ciarlet P C. The finite element methods for elliptic problems [M]. Amsterdam：NorthHolland, 1978：1-50.

[3] Vachaud G, Thony J L. Hyteresis during infiltration and redistribution in a soil column at different initial water contents [J]. Water Resources Research, 1971, 7：111-127.

[4] 雷志棟, 楊詩秀, 謝森傳. 土壤水動力學(xué)[M]. 北京：清華大學(xué)出版社, 1988：270-319.

[5] Adams R A. Sobolev spaces [M]. New York：Academ ic Press, 1975：26-27.

[6] O. C. Zienkiew icz. 有限元法[M]. 尹澤勇, 柴家振, 譯. 科學(xué)出版社, 1985：445-473.

[7] 王芳, 孫西歡, 馬娟娟, 等. 不同水頭作用下土壤垂直入滲特性試驗研究[J]. 山西水利, 2005（1）：75-76.

[8] Van Genuchten M Th, Leij F J. Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils [M]. California：University of California Press, 1992：1-14

A Finite Element Numerical Simulation for Unsaturated Richards Model w ith Complex Boundary

WANG Fang, YANG Hushan
（Department of Mathematics, Xinzhou Teacher’s College, Xinzhou, China 034000）

Through the method of the lumped mass of the finite element, the paper implements the numerical simulation and the model validation toward the mathematical model of the complex one-dimensional boundary soil water movement （unsaturated Richards model）, the results of the model validation show that this method possesses a very high accuracy of calculation which is in line w ith the requirements of the practical work.

Soil Water Movement；Finite Element；Numerical Simulation

O241.82

A

1674-3563（2013）01-0031-05

10.3875/j.issn.1674-3563.2013.01.006 本文的PDF文件可以從xuebao.wzu.edu.cn獲得

（編輯：王一芳）

2012-03-06

王芳（1978- ），女，河南焦作人，講師，碩士，研究方向：偏微分方程數(shù)值模擬，運籌學(xué)