一維非飽和土壤水分運動模型的數(shù)值模擬

梁晨璟+李春光+趙文娟+王建華

摘要：分析有限體積法在一維非飽和土壤水分運動模型的數(shù)值算法上的適用性及優(yōu)劣，通過算例對有限體積法進行驗證。結(jié)果表明，有限體積法同樣適用于一維非飽和土壤水分運動，與大多數(shù)學(xué)者都采用的有限差分法計算結(jié)果比較發(fā)現(xiàn)兩種算法的數(shù)值模擬結(jié)果相近。但是改變網(wǎng)格疏密后，對有限體積法的計算結(jié)果影響不大且與實測值相近，而對有限差分法的計算結(jié)果有影響，從而說明在網(wǎng)格稀疏的情況下，有限體積法更好。

關(guān)鍵詞：有限體積法；有限差分法；非飽和土壤；水分運動

中圖分類號：S152.7； S11+9 文獻標(biāo)識碼：A 文章編號：0439-8114（2014）15-3525-04

1-D Numerical Simulation of Water Movement in Unsaturated Soil

LIANG Chen-jing1，LI Chun-guang1，2，ZHAO Wen-juan1，WANG Jian-hua3

（1.Civil Engineering and Water Conservancy， Ningxia University， Yinchuan 750021， China； 2. Beifang University of Nationalities， Yinchuan 750000， China；3. Zhejiang Construction Supervision Co.， Ltd，Hangzhou 310012，China）

Abstract： The applicability and pros or cons of finite volume method in one-dimensional model of water movement in unsaturated soil were analyzed. The finite volume method was verified by a numerical example. The results showed that the finite volume method was appliable to one-dimensional model of water movement in unsaturated soil. Compared with the finite difference method used by most scholars， the numerical simulation results of two algorithms were similar. Changing the mesh density on the finite volume method had little effect on the calculated results similiar to the measured values， while the finite difference method calculation was affected. In the case of sparse grid， the finite volume method was better.

Key words： finite volume method；finite difference method；unsaturated soil； water movement

收稿日期：2014-03-15

基金項目：國家自然科學(xué)基金重大研究計劃培養(yǎng)項目（91230111）

作者簡介：梁晨璟（1988-），女，河北邢臺人，在讀碩士研究生，研究方向為旱區(qū)河流泥沙動力學(xué)理論與數(shù)值模擬，（電話）15909510629（電子信

箱）liangmeng88-12@163.com；通訊作者，李春光，男，教授，博士，主要從事計算力學(xué)和流體力學(xué)研究，（電子信箱）

cglizd@hotmail.com。

對于一維非飽和土壤的水分運動模型研究，大多數(shù)學(xué)者都是采用有限差分法來進行數(shù)值模擬。這種方法雖然形式簡單，且任意偏微分方程都能寫出與之對應(yīng)的差分方程，但不能從差分方程中體現(xiàn)出各項的物理意義。因此，差分方程只能認(rèn)為是對微分方程的數(shù)學(xué)近似，基本上不能反映出物理特征。差分方程在計算結(jié)果上也有可能出現(xiàn)不合理的情況[1]。大量的計算實踐表明，用通常的差分方法求解土壤水鹽運移問題，只有當(dāng)彌散作用占優(yōu)時，才能取得較為滿意的結(jié)果。在求解彌散系數(shù)較小、流速相對較大的對流占優(yōu)問題時，經(jīng)常遇到數(shù)值彌散和數(shù)值振蕩現(xiàn)象[2]。而有限體積法[3]（Finite volume methed，F(xiàn)VM）是以積分型守衡方程為出發(fā)點，通過對流體運動的體積域的離散來構(gòu)造積分型離散方程。可以看出，這種方法便于用來模擬具有復(fù)雜邊界區(qū)域流體的運動。有限體積法的基本思路易于理解，并能得出直接的物理解釋。離散方程的物理意義就是因變量在有限大小的控制體積中的守恒原理，如同微分方程表示因變量在無限小的控制體積中的守恒原理一樣。有限體積法得出的離散方程對任意一組控制體積都得到滿足因變量的積分守恒，因此也滿足整個計算區(qū)域。正是因為有限體積法有這么多優(yōu)點，所以在一維非飽和土壤的水分運動模型研究中嘗試采用有限體積法。

1 非飽和土壤水分運動數(shù)學(xué)模型

田間水分運動是非常復(fù)雜的，受多種因素影響諸如降雨、畦灌、噴灌及蒸發(fā)蒸騰等。因此，有多種數(shù)學(xué)模型。這里只研究最簡單的地表濕潤條件下的土壤水分運動問題的數(shù)學(xué)模型，假定土壤均質(zhì)，各向同性，根據(jù)達西定律和質(zhì)量守恒方程，可得到如下數(shù)學(xué)模型[4]：

=

[D（θ）

]-

θ=θa t=0，z=0

θ=θb t=0，z→∞或z=L

式中θ為含水率（cm3/cm3），θa為均勻分布的初始含水率（cm3/cm3），θb為地表因濕潤條件而維持不變的含水率（cm3/cm3）， D（θ）為擴散率（cm2/min）， K（θ）為導(dǎo)水率（cm/min）。z為位置坐標(biāo)，地表取為原點，向下為正方向。

2 數(shù)值模型模擬

2.1 建立數(shù)值模型

將全部求解區(qū)域z=0-L離散化為n個單元，共為n個節(jié)點，但上邊界點和下邊界點并沒有計入。距離步長為Δz，時間步長為Δt。

對2到n-1的節(jié)點，按有限體積法的完全隱格式可把原方程寫出如下形式：

（θik+1-θik）=

-

-

K（

θ）-K（

θ）

上述形式可化簡為如下形式：

Aθ+Bθ+Cθ=b

式中，

A=-

B=++

Ci=-

bi=θ-[K（θ）-K（θ）]

i=2，3，4，……， n

對于上邊界第一個節(jié)點，原方程可以寫成如下形式：

（θ-θ）=

-

-

[K（θ）-K（θ）]

可化簡為Bθ+Cθ=b這樣的形式

對于下邊界第n個節(jié)點，原方程可以寫成如下形式：

（θ-θ）=

-

-

[K（θ）-K（θ）]

同樣也可以可化簡為Aθ+Bθ=b

從而第一個節(jié)點到第n個節(jié)點都可把原方程轉(zhuǎn)化成a+aθ+aθ=b這樣的形式：

B1 C1 0 0 0

A2 B2 C2 0 0

0 ? ? ?

0 0 An-1 Bn-1 Cn-1

0 0 0 An Bn ·θ

θ

[…]

θ

+

θ

+=b1

b2

[…]

bn-1

bn-2

即三對角矩陣[A][θ]k+1=[b]，此時只需利用追趕法便可進行求解。

其中，對于擴散項中D（θ）這項采用算數(shù)平均法，即：

D（θ）=

D（θ）=

對于對流項中，K（θ）這項采用迎風(fēng)格式，即：

K（θ）=K（θ）

K（θ）=K（θ）

2.2模型的實現(xiàn)

上述模型的運算是基于Matlab數(shù)學(xué)軟件。程序采用有限體積法進行計算，經(jīng)過方程的離散后，方程便化寫成一個三對角矩陣。對于三對角矩陣采用的是TDMA法（追趕法）進行求解。TAMA算法是一種直接解法，是由高斯消元法應(yīng)用于這種特殊方程而得到的一種解法。對于數(shù)學(xué)模型中的土壤擴散率D（θ）和土壤導(dǎo)水率K（θ）均是由試驗測得。在程序中一般是定步長，根據(jù)需要改變步長。

2.3算例驗證

對于一個程序的有效性是需要算例驗證的。該程序的驗證算例為地表濕潤條件下的土壤水分運動問題[5]。土壤的垂向深度是100 cm，地表含水率θa=0.38 cm3/cm3，深度在100 cm時含水率θb=0.28 cm3/cm3，初始含水率θ0=0.28 cm3/cm3，飽和含水率θs=0.4 cm3/cm3，則可以表示成如下形式：

=

[D（θ）

]-

θa=0.38 t=0，z≥0

θb=0.28 t=0，z=100 cm

K（θ）=10.25（θ/θs）3.6

D（θ）=0.0379（θ/θs）3.487

式中，D（θ）為擴散率（cm2/min），K（θ）為導(dǎo)水率（cm/min）。設(shè)置時間步長為1 min，空間步長1 cm，用Matlab編程序進行數(shù)值模擬，其計算結(jié)果如圖1。

從圖1可以看出，土壤含水率隨著時間的推移，逐步向下入滲。當(dāng)t=5 min時入滲到20 cm處，在t=25 min時入滲到35 cm左右，最后到t=55 min時，入滲到了65 cm左右處。從整個模擬的趨勢上看是符合客觀規(guī)律的。從而也能說明，有限體積法也可以對一維非飽和土壤水分運動進行數(shù)值模擬。

3 有限體積法和有限差分法比較

在實驗室模擬土壤水分運動，只考慮地表濕潤條件，不考慮其他情況。選取的土壤為中壤土，初始含水率為0.09 cm3/cm3，土樣容重為1.4 g/cm3。為保證土柱初始含水率均勻和密度均一，先對土樣進行了風(fēng)干和過篩，然后土樣按控制的容重裝入直徑為10 cm，高為80 cm的有機玻璃管，土柱長60 cm。在土柱的頂部加水，形成一個薄水層（1 cm），由自動加水器供水，使水位保持不變，并計時。然后用γ射線測得1 230 min時刻的土柱剖面含水率。

數(shù)學(xué)模型如下：

=

[D（θ）

]-

θa=0.42 t=0，z≥0 （1）

θb=0.09 t=0，z=60 cm （2）

D（θ）=0.0002e19.902θ cm2/min

K（θ）=555.22θ13.181 cm/min

（1）式為上邊界條件，即在表層含水率0.42 cm3/cm3。（2）式為下邊界條件，垂直深度60 cm時含水率為0.09 cm3/cm3。初始含水率為0.09 cm3/cm3。D（θ）擴散率（cm2/min）是通過水平土柱吸滲法測定得出。K（θ）導(dǎo)水率（cm/min）由K（θ）=C·D（θ）導(dǎo)出，其中C為比水容量，通過土壤水分特征曲線得出。

當(dāng)時間步長為30 min，空間步長為2 cm，網(wǎng)格數(shù)為30時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結(jié)果如圖2。圖2中，可明顯看出有限體積法模擬效果優(yōu)于有限差分法，此時的網(wǎng)格數(shù)為30。

當(dāng)時間步長不變?nèi)詾?0 min，空間步長為1cm，網(wǎng)格數(shù)為60時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結(jié)果如圖3。圖3中，有限體積法和有限差分法兩者的模擬的效果基本相同，并沒有太大差別。此時的網(wǎng)格數(shù)為60，網(wǎng)格加密了一倍。與圖2比較，有限體積法并沒太大的改變，而有限差分法模擬效果明顯變好。

當(dāng)時間步長不變?nèi)詾?0 min，空間步長為2/3cm，網(wǎng)格數(shù)為90時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結(jié)果如圖4。圖4中，有限體積法和有限差分法模擬結(jié)果基本重合在一起。此時的網(wǎng)格數(shù)為90，較之圖2網(wǎng)格加密了兩倍，有限差分法模擬效果明顯優(yōu)于圖2。

4 結(jié)論

經(jīng)過算例的驗證，有限體積法同樣適用于一維非飽和土壤水分運動。與有限差分法比對發(fā)現(xiàn)，兩種方法的數(shù)值模擬結(jié)果相近。但是在網(wǎng)格稀疏的時候，有限體積法模擬的效果更好。在數(shù)值計算中，網(wǎng)格的疏密會影響計算精度和計算時間。網(wǎng)格越密計算精度越高，但是運算時間也會隨之增長。通過改變網(wǎng)格的疏密后有限體積法的計算值與實測值并無太大誤差，表明可選擇有限體積法在網(wǎng)格相對較疏時進行數(shù)值運算，效果較好。

參考文獻：

[1] 李人憲.有限體積法基礎(chǔ)[M].北京：國防工業(yè)出版社，2008.

[2] 呂歲菊，李春光.土壤水鹽運移規(guī)律數(shù)值模擬研究綜述[J].農(nóng)業(yè)科學(xué)研究，2005，26（1）：80-84.

[3] VERSTEEG H K， MALALASEKERA W. An Introduction to Computational Fluid Dynamics： The Finite V Method[M]. Edinburgh， England： Pearson Education Limited，1995.

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

[5] 趙慧麗，張 彌.降雨入滲對基坑非飽和土瞬態(tài)含水率影響的數(shù)值計算[J].石家莊鐵道學(xué)院學(xué)報，2000，12（13）：43-46.

當(dāng)時間步長為30 min，空間步長為2 cm，網(wǎng)格數(shù)為30時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結(jié)果如圖2。圖2中，可明顯看出有限體積法模擬效果優(yōu)于有限差分法，此時的網(wǎng)格數(shù)為30。

當(dāng)時間步長不變?nèi)詾?0 min，空間步長為1cm，網(wǎng)格數(shù)為60時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結(jié)果如圖3。圖3中，有限體積法和有限差分法兩者的模擬的效果基本相同，并沒有太大差別。此時的網(wǎng)格數(shù)為60，網(wǎng)格加密了一倍。與圖2比較，有限體積法并沒太大的改變，而有限差分法模擬效果明顯變好。

當(dāng)時間步長不變?nèi)詾?0 min，空間步長為2/3cm，網(wǎng)格數(shù)為90時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結(jié)果如圖4。圖4中，有限體積法和有限差分法模擬結(jié)果基本重合在一起。此時的網(wǎng)格數(shù)為90，較之圖2網(wǎng)格加密了兩倍，有限差分法模擬效果明顯優(yōu)于圖2。

4 結(jié)論

經(jīng)過算例的驗證，有限體積法同樣適用于一維非飽和土壤水分運動。與有限差分法比對發(fā)現(xiàn)，兩種方法的數(shù)值模擬結(jié)果相近。但是在網(wǎng)格稀疏的時候，有限體積法模擬的效果更好。在數(shù)值計算中，網(wǎng)格的疏密會影響計算精度和計算時間。網(wǎng)格越密計算精度越高，但是運算時間也會隨之增長。通過改變網(wǎng)格的疏密后有限體積法的計算值與實測值并無太大誤差，表明可選擇有限體積法在網(wǎng)格相對較疏時進行數(shù)值運算，效果較好。

參考文獻：

[1] 李人憲.有限體積法基礎(chǔ)[M].北京：國防工業(yè)出版社，2008.

[2] 呂歲菊，李春光.土壤水鹽運移規(guī)律數(shù)值模擬研究綜述[J].農(nóng)業(yè)科學(xué)研究，2005，26（1）：80-84.

[3] VERSTEEG H K， MALALASEKERA W. An Introduction to Computational Fluid Dynamics： The Finite V Method[M]. Edinburgh， England： Pearson Education Limited，1995.

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

[5] 趙慧麗，張 彌.降雨入滲對基坑非飽和土瞬態(tài)含水率影響的數(shù)值計算[J].石家莊鐵道學(xué)院學(xué)報，2000，12（13）：43-46.

當(dāng)時間步長為30 min，空間步長為2 cm，網(wǎng)格數(shù)為30時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結(jié)果如圖2。圖2中，可明顯看出有限體積法模擬效果優(yōu)于有限差分法，此時的網(wǎng)格數(shù)為30。

當(dāng)時間步長不變?nèi)詾?0 min，空間步長為1cm，網(wǎng)格數(shù)為60時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結(jié)果如圖3。圖3中，有限體積法和有限差分法兩者的模擬的效果基本相同，并沒有太大差別。此時的網(wǎng)格數(shù)為60，網(wǎng)格加密了一倍。與圖2比較，有限體積法并沒太大的改變，而有限差分法模擬效果明顯變好。

當(dāng)時間步長不變?nèi)詾?0 min，空間步長為2/3cm，網(wǎng)格數(shù)為90時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結(jié)果如圖4。圖4中，有限體積法和有限差分法模擬結(jié)果基本重合在一起。此時的網(wǎng)格數(shù)為90，較之圖2網(wǎng)格加密了兩倍，有限差分法模擬效果明顯優(yōu)于圖2。

4 結(jié)論

經(jīng)過算例的驗證，有限體積法同樣適用于一維非飽和土壤水分運動。與有限差分法比對發(fā)現(xiàn)，兩種方法的數(shù)值模擬結(jié)果相近。但是在網(wǎng)格稀疏的時候，有限體積法模擬的效果更好。在數(shù)值計算中，網(wǎng)格的疏密會影響計算精度和計算時間。網(wǎng)格越密計算精度越高，但是運算時間也會隨之增長。通過改變網(wǎng)格的疏密后有限體積法的計算值與實測值并無太大誤差，表明可選擇有限體積法在網(wǎng)格相對較疏時進行數(shù)值運算，效果較好。

參考文獻：

[1] 李人憲.有限體積法基礎(chǔ)[M].北京：國防工業(yè)出版社，2008.

[2] 呂歲菊，李春光.土壤水鹽運移規(guī)律數(shù)值模擬研究綜述[J].農(nóng)業(yè)科學(xué)研究，2005，26（1）：80-84.

[3] VERSTEEG H K， MALALASEKERA W. An Introduction to Computational Fluid Dynamics： The Finite V Method[M]. Edinburgh， England： Pearson Education Limited，1995.

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

[5] 趙慧麗，張 彌.降雨入滲對基坑非飽和土瞬態(tài)含水率影響的數(shù)值計算[J].石家莊鐵道學(xué)院學(xué)報，2000，12（13）：43-46.