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

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

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

中圖分類號(hào)：S152.7； S11+9 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：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

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

作者簡(jiǎn)介：梁晨璟（1988-），女，河北邢臺(tái)人，在讀碩士研究生，研究方向?yàn)楹祬^(qū)河流泥沙動(dòng)力學(xué)理論與數(shù)值模擬，（電話）15909510629（電子信

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

cglizd@hotmail.com。

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

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

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

=

[D（θ）

]-

θ=θa t=0，z=0

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

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

2 數(shù)值模型模擬

2.1 建立數(shù)值模型

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

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

（θik+1-θik）=

-

-

K（

θ）-K（

θ）

上述形式可化簡(jiǎn)為如下形式：

Aθ+Bθ+Cθ=b

式中，

A=-

B=++

Ci=-

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

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

對(duì)于上邊界第一個(gè)節(jié)點(diǎn)，原方程可以寫成如下形式：

（θ-θ）=

-

-

[K（θ）-K（θ）]

可化簡(jiǎn)為Bθ+Cθ=b這樣的形式

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

（θ-θ）=

-

-

[K（θ）-K（θ）]

同樣也可以可化簡(jiǎn)為Aθ+Bθ=b

從而第一個(gè)節(jié)點(diǎn)到第n個(gè)節(jié)點(diǎn)都可把原方程轉(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

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

其中，對(duì)于擴(kuò)散項(xiàng)中D（θ）這項(xiàng)采用算數(shù)平均法，即：

D（θ）=

D（θ）=

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

K（θ）=K（θ）

K（θ）=K（θ）

2.2模型的實(shí)現(xiàn)

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

2.3算例驗(yàn)證

對(duì)于一個(gè)程序的有效性是需要算例驗(yàn)證的。該程序的驗(yàn)證算例為地表濕潤(rùn)條件下的土壤水分運(yùn)動(dòng)問題[5]。土壤的垂向深度是100 cm，地表含水率θa=0.38 cm3/cm3，深度在100 cm時(shí)含水率θ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（θ）為擴(kuò)散率（cm2/min），K（θ）為導(dǎo)水率（cm/min）。設(shè)置時(shí)間步長(zhǎng)為1 min，空間步長(zhǎng)1 cm，用Matlab編程序進(jìn)行數(shù)值模擬，其計(jì)算結(jié)果如圖1。

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

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

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

數(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時(shí)含水率為0.09 cm3/cm3。初始含水率為0.09 cm3/cm3。D（θ）擴(kuò)散率（cm2/min）是通過水平土柱吸滲法測(cè)定得出。K（θ）導(dǎo)水率（cm/min）由K（θ）=C·D（θ）導(dǎo)出，其中C為比水容量，通過土壤水分特征曲線得出。

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

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

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

4 結(jié)論

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

參考文獻(xiàn)：

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

[2] 呂歲菊，李春光.土壤水鹽運(yùn)移規(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] 雷志棟，楊師秀，謝森傳.土壤水動(dòng)力學(xué)[M].北京：清華大學(xué)出版社，1988.

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

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

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

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

4 結(jié)論

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

參考文獻(xiàn)：

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

[2] 呂歲菊，李春光.土壤水鹽運(yùn)移規(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] 雷志棟，楊師秀，謝森傳.土壤水動(dòng)力學(xué)[M].北京：清華大學(xué)出版社，1988.

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

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

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

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

4 結(jié)論

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

參考文獻(xiàn)：

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

[2] 呂歲菊，李春光.土壤水鹽運(yùn)移規(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] 雷志棟，楊師秀，謝森傳.土壤水動(dòng)力學(xué)[M].北京：清華大學(xué)出版社，1988.

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