土壤下滲問題的格子玻爾茲曼模擬

張東輝,芮孝芳,馬哲樹,孔祥雷

(1.江蘇科技大學(xué)船舶與海洋工程學(xué)院,江蘇鎮(zhèn)江 212003;2.河海大學(xué)水文水資源學(xué)院,江蘇南京 210098)

應(yīng)用連續(xù)介質(zhì)模型來分析流體流動(dòng)問題,實(shí)際上是通過分析微元體的質(zhì)量守恒、動(dòng)量定理、能量守恒推得N-S方程組,然后應(yīng)用有限差分法等方法對(duì)微分方程進(jìn)行離散,得到各結(jié)點(diǎn)的代數(shù)方程組進(jìn)行求解.這類計(jì)算方法比較直觀,目前應(yīng)用非常廣泛.但如果微分方程是強(qiáng)非線性的,數(shù)值穩(wěn)定性便成為很大的問題,往往需要通過變換和設(shè)計(jì)特殊的數(shù)值格式來分析具體問題[1].

與上述思路不同,研究流體的流動(dòng)也可從微觀離散模型出發(fā),即從更微觀的尺度——大量分子運(yùn)動(dòng)的主控方程出發(fā),向上“集成”來探討流體流動(dòng)特性[2],如圖1所示.統(tǒng)計(jì)熱力學(xué)中的玻爾茲曼方程[2]是描述系統(tǒng)中大量粒子分布函數(shù)f(x,v,t)的演化方程,它是聯(lián)系微觀與宏觀的橋梁,其簡(jiǎn)化形式為

圖1 兩種推導(dǎo)宏觀微分方程的途徑Fig.1 Two kinds of processes deriving macroscopic differential equation

式中:e——微觀粒子的速度,e=(ex,ey,ez);x——粒子的位置矢量,x=(x,y,z);F——作用在分子上的外力(也稱源項(xiàng));Ω(f*,f)——粒子與粒子之間的碰撞散射項(xiàng),表示粒子分布函數(shù)f(x,e,t)由于粒子碰撞而發(fā)生的改變率.玻爾茲曼方程的解是不同時(shí)間的粒子分布函數(shù)f(x,e,t).一旦求得f(x,e,t)的具體表達(dá)式,可通過它與質(zhì)量、速度、含水率等宏觀量之間的關(guān)系,來推得后者的變化特點(diǎn).

格子玻爾茲曼方法采用網(wǎng)格離散的方法求解玻爾茲曼方程[3].如果不考慮外力項(xiàng)的作用,式(1)在x處、α方向上可離散為[2]

式中:eα——速度矢量;τ——弛豫系數(shù);feqα(x,t)——平衡態(tài)分布函數(shù).

由于土壤含水率與土水勢(shì)(土壤水吸力)呈高度非線性關(guān)系,因而下滲方程通常是屬于強(qiáng)非線性偏微分方程.對(duì)此問題的研究具有非常重要的意義,一方面可以了解水分或溶質(zhì)在土壤中的運(yùn)移和分布情況[4];另一方面也可從中深入分析產(chǎn)流的機(jī)理.目前下滲方程的求解主要是采用解析解法、有限差分法、有限元法和基爾霍夫變換法等[1],這些方法大多都存在計(jì)算精度、穩(wěn)定性和計(jì)算效率之間的矛盾.應(yīng)用格子玻爾茲曼方法來求解下滲方程是一種新的嘗試,格子玻爾茲曼方法中蘊(yùn)涵了多尺度展開的思想,而且具有編程簡(jiǎn)潔、邊界易處理、易考慮粒子之間的相關(guān)作用等特點(diǎn),因此它為復(fù)雜系統(tǒng)的模擬提供了新的途徑,已成功應(yīng)用于很多非線性問題的研究[5].格子玻爾茲曼方法應(yīng)用于下滲問題的討論有兩類模型:第一類模型是建立在具體的土壤孔隙結(jié)構(gòu)信息基礎(chǔ)上,然后模擬其中的流動(dòng),這樣可很好描述“孔隙尺度(pore scale)”層次上的流動(dòng)[4-6].這種模型非常直觀,但由于土壤孔隙結(jié)構(gòu)的信息存儲(chǔ)問題,對(duì)所研究問題的尺度是有限制的,很難將這種方法推廣到大尺度下滲特性的研究.另一類模型是直接求解Richards方程,Ginzburg[7-8]發(fā)展出一種各向異性的格子玻爾茲曼模型(Lattice Boltzmann model,簡(jiǎn)稱LB模型),對(duì)不同類型的土壤下滲過程進(jìn)行了模擬,并考慮了水力傳導(dǎo)率沿深度的變化,其模擬結(jié)果令人滿意,但其過程非?；逎?離真正實(shí)用推廣還有一段距離;近些年,段雅麗等[9-10]提出了一種新的LB模型,成功求解了一定邊界條件下的KDV方程(淺水波方程)、Sine-Gordon方程等,本文應(yīng)用這種新的LB模型對(duì)下滲方程進(jìn)行討論分析.

1 物理模型

1.1 Richards下滲方程

土壤水流下滲過程可由Richards方程來描述:

這里所取的邊界條件為

式中:K——水力傳導(dǎo)率;D——水力擴(kuò)散率;θ——土壤含水率;θ0——初始土壤含水率;θn——土壤表層含水率;t——下滲時(shí)間;z——土壤深度;L——土壤總深度.

如果水力傳導(dǎo)率K、水力擴(kuò)散率D與土壤含水率θ之間關(guān)系滿足常用的Brooks-Corey公式時(shí),即有:

式中 α和β是常數(shù).方程(3)化為

其等價(jià)表達(dá)式為

當(dāng)不考慮重力下滲的影響時(shí),即水力傳導(dǎo)率K=0時(shí),上述方程即變?yōu)閿U(kuò)散方程:

Richards方程(3)只有在一些簡(jiǎn)化條件下才能獲得分析解[1]:

a.當(dāng)水力擴(kuò)散率D為常數(shù)和水力傳導(dǎo)率K=0時(shí),方程(5)演變?yōu)榫€性擴(kuò)散方程,當(dāng)邊界條件滿足式(4)時(shí),一定時(shí)間后不同深度的土壤含水率分布滿足[1]:

erfc在數(shù)學(xué)上稱為余誤差函數(shù).

b.當(dāng)水力擴(kuò)散率D為常數(shù)以及水力傳導(dǎo)率K與土壤含水率θ呈線性關(guān)系時(shí),方程(5)演變?yōu)榫€性Richards方程,相同的邊界條件下,一定時(shí)間后不同深度的土壤含水率分布滿足[1]:

1.2 Richards方程的LB模型

本文所采用的格子玻爾茲曼模型對(duì)空間坐標(biāo)并不進(jìn)行多重尺度化處理,而對(duì)時(shí)間坐標(biāo)引入3個(gè)時(shí)間尺度.對(duì)于一維LB模型網(wǎng)格,常用的是3速模型(α=0～2),如圖2所示,每一網(wǎng)格結(jié)點(diǎn)上,存在3類不同速度的粒子,速度大小分別是-v,0,v.由于LB模型中的運(yùn)動(dòng)粒子在每一時(shí)間步長(zhǎng)下都是移動(dòng)單個(gè)的網(wǎng)格步長(zhǎng),C0定義為網(wǎng)格步長(zhǎng)Δx與時(shí)間步長(zhǎng)Δt之比,因而C0=v=Δx/Δt.

那么時(shí)間和空間坐標(biāo)的多重尺度表達(dá)為

圖2 3速模型示意圖Fig.2 3-bit velocity model

空間坐標(biāo)的微分形式不變,時(shí)間坐標(biāo)的微分其多重尺度的形式為

粒子分布函數(shù)同樣可展開為

對(duì)上式兩邊求和,根據(jù)質(zhì)量守恒可得:

由此推得:

對(duì)離散方程(2)進(jìn)行多重尺度處理,即可得到3個(gè)不同時(shí)間尺度的玻爾茲曼方程:

選擇3速模型平衡態(tài)分布函數(shù)的各階矩的形式為

為還原得到Richards方程,首先需對(duì)不同時(shí)間尺度的玻爾茲曼方程在不同方向上求和,然后將不同時(shí)間尺度對(duì)應(yīng)的方程按不同數(shù)量級(jí)合并可得:

只要令B=1/(τ-0.5)ε,即可還原得到二階精度的Richards方程(6),還原過程主要是保證宏觀量選擇和多尺度方案的正確性.

由于在計(jì)算中必須要知道平衡態(tài)分布函數(shù),對(duì)于三速模型,3個(gè)方向的平衡態(tài)分布函數(shù)的求取可由式(15)確定.利用方程的對(duì)稱性,可求得:

本文應(yīng)用格子玻爾茲曼方法討論下滲問題時(shí),先從線性擴(kuò)散方程和線性Richards方程開始,由于它們都存在分析解,這樣一方面可驗(yàn)證LB模型的正確性,另一方面也可了解其誤差特性.

2 格子玻爾茲曼方法的分析結(jié)果

如果僅考慮水力擴(kuò)散效應(yīng),那么水分在土壤中的下滲過程由擴(kuò)散方程(7)描述.如果水力擴(kuò)散率D為常數(shù),給定上邊界含水率和初始土壤含水率,不同時(shí)間下不同深度的土壤含水率分布則呈余誤差函數(shù)形式,即由式(8)決定.采用格子玻爾茲曼方法模擬擴(kuò)散方程時(shí),計(jì)算所取的參數(shù)為:初始土壤含水率 θ0=0.028,土壤表層含水率 θn=0.45,水力擴(kuò)散率D分別為42.4cm2/min和424cm2/min.土壤總深度L=10m,網(wǎng)格結(jié)點(diǎn)數(shù)N=201,網(wǎng)格步長(zhǎng)Δx=L/N=0.05m,時(shí)間步長(zhǎng)Δt=0.01s,C0=Δx/Δt=5m/s,弛豫系數(shù) τ=1.5.不同時(shí)間的土壤含水率分布如圖3所示,LB模型的計(jì)算結(jié)果與分析解符合得非常好.水力擴(kuò)散率越大,相同時(shí)間下深層土壤的含水率就越大.

圖4是線性Richards方程的LB模型的計(jì)算結(jié)果,其水力擴(kuò)散率D和水力傳導(dǎo)率K與土壤含水率θ的關(guān)系是:D=42.4cm2/min,K=6θ cm/min.

圖3 不同水力擴(kuò)散率時(shí)水分在土壤中不同時(shí)間的分布Fig.3 Distribution of water in soils with different diffusion coefficients at different time

圖4 線性 Richards方程的LB模型計(jì)算結(jié)果與分析解的比較Fig.4 Comparison between simulated results by LB model and analytical solutions of linear Richards equation

線性Richards方程考慮了重力下滲的影響,含水率在土壤中呈類似“S”形分布.其計(jì)算結(jié)果和真實(shí)土壤的含水率分布已經(jīng)比較接近,在土壤上部存在“飽和帶”,但濕潤(rùn)鋒面并沒有出現(xiàn),這是線性Richards方程美中不足之處.

為了解各種因素的影響,可將格子玻爾茲曼方法的模擬結(jié)果與分析解進(jìn)行比較,這里應(yīng)用相對(duì)總體誤差指標(biāo)來衡量,其定義是:

圖5(a)是線性擴(kuò)散方程的LB模型的誤差特性分析.對(duì)一定的擴(kuò)散時(shí)間,存在一個(gè)最佳的弛豫系數(shù)τ,大約為1.0,此時(shí)計(jì)算誤差最小,大約是0.1%.當(dāng)弛豫系數(shù)增大或減小時(shí)(偏離最佳弛豫系數(shù)),計(jì)算誤差幾乎呈直線上升;在圖5(a)中,上下兩簇曲線分別代表不同的網(wǎng)格步長(zhǎng)(時(shí)間步長(zhǎng)均相同),這一結(jié)果表明:對(duì)應(yīng)一定的時(shí)間步長(zhǎng),存在最佳的網(wǎng)格步長(zhǎng),即C0=Δx/Δt存在最優(yōu)值.圖5(b)是應(yīng)用格子玻爾茲曼方法求解線性Richards方程的誤差特性分析.這里仍然采用相對(duì)總體誤差,它主要受到速度模型、邊界格式、弛豫系數(shù)和格子長(zhǎng)度等因素的制約.由于上下邊界的含水率已經(jīng)設(shè)定,邊界節(jié)點(diǎn)各時(shí)段的分布函數(shù)用平衡態(tài)分布函數(shù)取定,所以沒有討論邊界格式對(duì)計(jì)算誤差的影響.

從圖5(a)和圖5(b)的比較可以發(fā)現(xiàn):線性Richards方程的誤差特性要更為復(fù)雜一些,其最佳弛豫系數(shù)隨著下滲時(shí)間的延長(zhǎng),從1.0變?yōu)?.0,而且在本文計(jì)算的下滲時(shí)間范圍內(nèi),隨下滲時(shí)間的增大,最小計(jì)算誤差是先增大后減小.對(duì)一定的下滲時(shí)間,存在一個(gè)最佳的弛豫系數(shù),最小計(jì)算誤差可達(dá)到0.1%.當(dāng)下滲時(shí)間較短時(shí),最佳弛豫系數(shù)在1.0附近;當(dāng)下滲時(shí)間較長(zhǎng)時(shí),最佳弛豫系數(shù)是在3.0左右.這和大多數(shù)研究的結(jié)論不同,以往很多研究一般都認(rèn)為弛豫系數(shù)范圍在0.5～1.0,但本文通過對(duì)計(jì)算誤差的分析,發(fā)現(xiàn)其最佳值有可能大于這一范圍.當(dāng)偏離最佳弛豫系數(shù),計(jì)算誤差幾乎呈直線上升.

圖5 LB模型誤差特性分析Fig.5 Error analysis of LB model

為了了解真實(shí)土壤中的下滲情況,需直接求解非線性的Richards方程.這里采用了文獻(xiàn)[1]中的算例,其中水力傳導(dǎo)率K、水力擴(kuò)散率D與土壤含水率θ及土壤飽和含水率θs之間的關(guān)系滿足Brooks-Corey公式[1]:

Richards方程就變換為

LB模型宏觀量取法如式(15)所示,其中m,n,α,β分別取:

土壤飽和含水率 θs=0.485,邊界條件和初始條件仍然是:

Richards方程的分析解可由Philip提出的方法求出[1],Philips解的缺點(diǎn)是當(dāng)下滲時(shí)間足夠長(zhǎng)時(shí)誤差比較大.圖6對(duì)LB模型的計(jì)算結(jié)果與Philips解[1]進(jìn)行了比較.在下滲歷時(shí)較短時(shí),兩者比較吻合;在一定的下滲時(shí)間后,則表現(xiàn)出一定的差別,由LB模型計(jì)算所得的濕潤(rùn)鋒下滲速度要略大于Philip解[1]的結(jié)果,兩者比較,最大計(jì)算誤差大約是15%,土壤水分分布形態(tài)則基本一致.圖7是下滲歷時(shí)較長(zhǎng)后,由LB模型計(jì)算所得的土壤水分分布特性.由圖7所見,“飽和區(qū)”和“濕潤(rùn)鋒”等典型的土壤下滲特征都可從中得到反映.這說明格子玻爾茲曼方法可以應(yīng)用于實(shí)際土壤的下滲模擬研究.

圖6 LB模型計(jì)算結(jié)果與Philips解[1]的比較Fig.6 Comparison between simulated results by LB model and Philips solutions[1]

圖7 非線性Richards方程LB模型的計(jì)算結(jié)果Fig.7 Simulated results by LB model for nonlinear Richards equation

由于Richards方程是屬于流體力學(xué)的Burgers方程,而Burgers方程常被用于激波問題的模擬,激波現(xiàn)象最明顯的特征是具有物理量梯度變化很大的前鋒面.在土壤下滲問題中,也存在類似現(xiàn)象,比如說濕潤(rùn)峰(即干濕界面的前移)的存在.目前有些文獻(xiàn)[11-12]就是討論如何將激波傳播的解的特性應(yīng)用于Richards方程的求解.兩者的區(qū)別在于:激波問題主要適用于超音速的可壓縮流動(dòng),而Richards方程適用于下滲速度很慢的孔隙性介質(zhì)流動(dòng).

3 結(jié) 語

與其他數(shù)值計(jì)算方法相比,本文采用的LB模型很好地解決了計(jì)算精度、穩(wěn)定性和計(jì)算效率之間的矛盾,在下滲問題的應(yīng)用上展現(xiàn)出可期待的前景.但研究中也發(fā)現(xiàn):此LB模型需要水力擴(kuò)散率與土壤含水率的函數(shù)關(guān)系可積,當(dāng)土壤水分特征曲線滿足Mualem-Van-Genuchten關(guān)系式時(shí),水力擴(kuò)散率與土壤含水率的函數(shù)關(guān)系表達(dá)式很難給出可積形式,對(duì)這一困難還需在未來研究中加以克服.

[1]雷志棟.土壤物理學(xué)[M].北京:科學(xué)出版社,1986:180-200.

[2]林宗涵.熱力學(xué)與統(tǒng)計(jì)物理學(xué)[M].北京:北京大學(xué)出版社,2007:265-280.

[3]CHOPARD B.物理系統(tǒng)的元胞自動(dòng)機(jī)模擬[M].北京:清華大學(xué)出版社,2003:20-78.

[4]VER MA N,SALEM K,MEWES D.Simulation of micro-and macro-transport in a packed bed of porous adsorbents by lattice Boltzmann method[J].Chemical Engineering Science,2007,62:3685-3698.

[5]YEOMANS J M,Mesoscale simulations:Lattice Boltzmann and particle algorithms[J].Physica A:Statistical and Theoretical Physics,2006,369:159-184.

[6]PAPAFOTIOU A,HELMIG R.From the pore scale to the lab scale:3-D lab experiment and numerical simulation of drainage in heterogeneous porousmedia[J].Advances in Water Resources,2008,31:1253-1268..

[7]GINZBURG I.Variably saturated flow described with the anisotropic Lattice Boltzmann methods[J].Computer&Fluids,2006,35:831-848.

[8]GINZBURG I.Lattice Boltzmann and analytical modeling of flow processes in anisotropic and heterogeneous stratified aquifers[J].Advances in Water Resources,2007,30:2202-2234.

[9]段雅麗.格子Boltzmann方法及其在流體動(dòng)力學(xué)上的一些應(yīng)用[D].合肥:中國(guó)科學(xué)技術(shù)大學(xué),2007.

[10]ZHANG Jian-ying,YAN Guang-wu.Lattice Boltzmann method for one and two-dimensional Burgers Equation[J].Physica A:Statistical Mechanics and its Applications,2008,387(19):4771-4786

[11]CAPUTO J,STEPANYANTS Y.Front Solutions of Richards'Equation[J].Transport in PorousMedia,2008,74:1-20.

[12]MILLER C,ABHISHEK C,Farthing M.A Spatially and temporally adaptive solution of Richards'equation[J].Advances in Water Resources,2006,29:525-545.