王 芳，楊虎山

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

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

王 芳，楊虎山

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

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

土壤水分運(yùn)動(dòng)；有限元；數(shù)值模擬

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

非飽和土壤水分運(yùn)動(dòng)問(wèn)題可歸結(jié)為下面的模型方程：

求h（z ,t ）使得對(duì)于任意的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）；為變量，是隨時(shí)間t變化的水頭值（cm），在實(shí)際計(jì)算時(shí)，采用插值函數(shù)進(jìn)行處理；L和T分別是空間和時(shí)間變化的上界．

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

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

1.1 基本假定

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

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

1.2 方程的離散

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

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

1.2.1 半離散有限元逼近

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

其中，

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

1.2.2 全離散格式

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

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

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

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

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

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

2 數(shù)值模擬

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

式中，s為土壤水吸力（cm）．非飽和導(dǎo)水率的測(cè)定有直接法和間接法兩種[4]．由于直接測(cè)定法耗時(shí)多、昂貴、繁瑣，且測(cè)定范圍較窄，可獲得的實(shí)驗(yàn)數(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

土壤水分運(yùn)動(dòng)的初始和邊界條件為：

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

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

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

3 結(jié)語(yǔ)

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

[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] 雷志棟, 楊詩(shī)秀, 謝森傳. 土壤水動(dòng)力學(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] 王芳, 孫西歡, 馬娟娟, 等. 不同水頭作用下土壤垂直入滲特性試驗(yàn)研究[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ù)值模擬，運(yùn)籌學(xué)