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

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

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

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

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

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

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

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

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

1 物理模型

1.1 Richards下滲方程

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

這里所取的邊界條件為

式中:K——水力傳導率;D——水力擴散率;θ——土壤含水率;θ0——初始土壤含水率;θn——土壤表層含水率;t——下滲時間;z——土壤深度;L——土壤總深度.

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

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

其等價表達式為

當不考慮重力下滲的影響時,即水力傳導率K=0時,上述方程即變?yōu)閿U散方程:

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

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

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

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

1.2 Richards方程的LB模型

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

那么時間和空間坐標的多重尺度表達為

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

空間坐標的微分形式不變,時間坐標的微分其多重尺度的形式為

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

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

由此推得:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Richards方程就變換為

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

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

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

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

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

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

3 結(jié) 語

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

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

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

[3]CHOPARD B.物理系統(tǒng)的元胞自動機模擬[M].北京:清華大學出版社,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].合肥:中國科學技術(shù)大學,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.