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

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

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

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

中圖分類號：S152.7； S11+9 文獻標識碼：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

基金項目：國家自然科學基金重大研究計劃培養項目（91230111）

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

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

cglizd@hotmail.com。

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

1 非飽和土壤水分運動數學模型

田間水分運動是非常復雜的，受多種因素影響諸如降雨、畦灌、噴灌及蒸發蒸騰等。因此，有多種數學模型。這里只研究最簡單的地表濕潤條件下的土壤水分運動問題的數學模型，假定土壤均質，各向同性，根據達西定律和質量守恒方程，可得到如下數學模型[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（θ）為導水率（cm/min）。z為位置坐標，地表取為原點，向下為正方向。

2 數值模型模擬

2.1 建立數值模型

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

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

（θik+1-θik）=

-

-

K（

θ）-K（

θ）

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

Aθ+Bθ+Cθ=b

式中，

A=-

B=++

Ci=-

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

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

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

（θ-θ）=

-

-

[K（θ）-K（θ）]

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

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

（θ-θ）=

-

-

[K（θ）-K（θ）]

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

從而第一個節點到第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（θ）這項采用算數平均法，即：

D（θ）=

D（θ）=

對于對流項中，K（θ）這項采用迎風格式，即：

K（θ）=K（θ）

K（θ）=K（θ）

2.2模型的實現

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

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（θ）為導水率（cm/min）。設置時間步長為1 min，空間步長1 cm，用Matlab編程序進行數值模擬，其計算結果如圖1。

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

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

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

數學模型如下：

=

[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（θ）導水率（cm/min）由K（θ）=C·D（θ）導出，其中C為比水容量，通過土壤水分特征曲線得出。

當時間步長為30 min，空間步長為2 cm，網格數為30時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結果如圖2。圖2中，可明顯看出有限體積法模擬效果優于有限差分法，此時的網格數為30。

當時間步長不變仍為30 min，空間步長為1cm，網格數為60時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結果如圖3。圖3中，有限體積法和有限差分法兩者的模擬的效果基本相同，并沒有太大差別。此時的網格數為60，網格加密了一倍。與圖2比較，有限體積法并沒太大的改變，而有限差分法模擬效果明顯變好。

當時間步長不變仍為30 min，空間步長為2/3cm，網格數為90時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結果如圖4。圖4中，有限體積法和有限差分法模擬結果基本重合在一起。此時的網格數為90，較之圖2網格加密了兩倍，有限差分法模擬效果明顯優于圖2。

4 結論

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

參考文獻：

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

[2] 呂歲菊，李春光.土壤水鹽運移規律數值模擬研究綜述[J].農業科學研究，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] 雷志棟，楊師秀，謝森傳.土壤水動力學[M].北京：清華大學出版社，1988.

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

當時間步長為30 min，空間步長為2 cm，網格數為30時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結果如圖2。圖2中，可明顯看出有限體積法模擬效果優于有限差分法，此時的網格數為30。

當時間步長不變仍為30 min，空間步長為1cm，網格數為60時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結果如圖3。圖3中，有限體積法和有限差分法兩者的模擬的效果基本相同，并沒有太大差別。此時的網格數為60，網格加密了一倍。與圖2比較，有限體積法并沒太大的改變，而有限差分法模擬效果明顯變好。

當時間步長不變仍為30 min，空間步長為2/3cm，網格數為90時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結果如圖4。圖4中，有限體積法和有限差分法模擬結果基本重合在一起。此時的網格數為90，較之圖2網格加密了兩倍，有限差分法模擬效果明顯優于圖2。

4 結論

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

參考文獻：

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

[2] 呂歲菊，李春光.土壤水鹽運移規律數值模擬研究綜述[J].農業科學研究，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] 雷志棟，楊師秀，謝森傳.土壤水動力學[M].北京：清華大學出版社，1988.

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

當時間步長為30 min，空間步長為2 cm，網格數為30時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結果如圖2。圖2中，可明顯看出有限體積法模擬效果優于有限差分法，此時的網格數為30。

當時間步長不變仍為30 min，空間步長為1cm，網格數為60時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結果如圖3。圖3中，有限體積法和有限差分法兩者的模擬的效果基本相同，并沒有太大差別。此時的網格數為60，網格加密了一倍。與圖2比較，有限體積法并沒太大的改變，而有限差分法模擬效果明顯變好。

當時間步長不變仍為30 min，空間步長為2/3cm，網格數為90時，計算1 230 min后的入滲情況。有限體積法和有限差分法的計算結果如圖4。圖4中，有限體積法和有限差分法模擬結果基本重合在一起。此時的網格數為90，較之圖2網格加密了兩倍，有限差分法模擬效果明顯優于圖2。

4 結論

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

參考文獻：

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

[2] 呂歲菊，李春光.土壤水鹽運移規律數值模擬研究綜述[J].農業科學研究，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] 雷志棟，楊師秀，謝森傳.土壤水動力學[M].北京：清華大學出版社，1988.

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