胡五龍 劉國峰 晏石林 范嚴(yán)偉

?(武漢理工大學(xué)綠色智能江海直達(dá)船舶與郵輪游艇研究中心,武漢 430070)

? (武漢理工大學(xué)理學(xué)院,武漢 430070)

??(蘭州理工大學(xué)能源與動(dòng)力工程學(xué)院,蘭州 730050)

引言

土壤是人類、植物及微生物維持生命活動(dòng)所需水分及營養(yǎng)物質(zhì)的重要來源[1];同時(shí)也是生態(tài)循環(huán)的主要組成部分[2].水在土壤孔隙中的分布直接影響土壤中物質(zhì),包括水、空氣及其他有機(jī)或無機(jī)物質(zhì)的存儲(chǔ)及傳輸[3],并對(duì)土壤中微生物的活動(dòng)范圍及強(qiáng)度有著決定性的影響[1].

土壤中水分分布一直廣受農(nóng)業(yè)、環(huán)境、生態(tài)等領(lǐng)域研究人員的關(guān)注,并在不同角度不同尺度得到了大量的研究.例如周啟友和島田純[4]通過野外觀測(cè)及反演計(jì)算,探討了土壤水在空間分布的非均質(zhì)性及各向異性.張常亮等[5]通過對(duì)天然降雨入滲后黃土中不同深度土層體積含水率為期一年的監(jiān)測(cè)數(shù)據(jù)進(jìn)行分析,確定了黃土地區(qū)降雨入滲深度及浸潤線以下水分遷移特點(diǎn).杜娟等[6]采用野外觀測(cè)手段研究了沙漠不同地形下土壤中水分時(shí)空變化特征,發(fā)現(xiàn)土壤滲透性與孔隙結(jié)構(gòu)、有機(jī)質(zhì)及水分分布有密切關(guān)系.這些研究加深了人們對(duì)土壤水分時(shí)空分布變異性及復(fù)雜性的理解,但大部分研究都是對(duì)宏觀現(xiàn)象的描述,土壤中的所有物理過程及生化反應(yīng)都發(fā)生在孔隙中,從孔隙尺度模擬非飽和土壤中水及溶質(zhì)的分布,有助于對(duì)土壤中各過程和反應(yīng)的理解.

土壤一般為非飽和狀態(tài),土壤中水的運(yùn)動(dòng)實(shí)際是復(fù)雜的多相流,在相當(dāng)長時(shí)間內(nèi),對(duì)于土壤水問題一直停留在定性描述或采用各種經(jīng)驗(yàn)方法處理的層面上[7].自1856 年Darcy[8]在大量實(shí)驗(yàn)的基礎(chǔ)上提出達(dá)西定律,1907 年Buckingham[9]提出毛管勢(shì)理論,到1931 年Richard[10]在達(dá)西定律的基礎(chǔ)上推導(dǎo)出用于描述非飽和流的方程,土壤水運(yùn)動(dòng)理論模型逐步建立起來.此后土壤水運(yùn)動(dòng)研究的關(guān)注點(diǎn)轉(zhuǎn)移到了基本方程的求解,而其中的關(guān)鍵之一是土壤水分特征參數(shù)的確定,1964 年Brooks 和Corey[11]建立了飽和度、壓力梯度及水力特性之間的關(guān)系,并提出利用毛管壓力?飽和度曲線來確定水力特性的方法;van Genuchten[12]提出了估算非飽和土壤水力特征的半經(jīng)驗(yàn)公式,該公式后來成為了描述土壤水力特性的主要公式之一.還有一些其他方法,包括轉(zhuǎn)換函數(shù)和神經(jīng)網(wǎng)絡(luò)等也被用于土壤水力特征及土壤水運(yùn)動(dòng)的研究中[13-14].國內(nèi)邵明安等[15-16]也在土壤水運(yùn)動(dòng)、水力特征參數(shù)測(cè)量與計(jì)算方面做了大量的基礎(chǔ)研究工作.2019 年Jin 等[17]在前人模型基礎(chǔ)上,將薄膜流模型和毛管流模型進(jìn)行整合,提出了一種可以描述全飽和度區(qū)間土水特征曲線的分形模型.以上這些模型對(duì)于理解和模擬土壤中水的分布起到了非常關(guān)鍵的作用,但此類模型及其參數(shù)都是土壤宏觀特征的平均表征,直接從孔隙尺度研究水在每個(gè)孔隙中的具體形態(tài)和分布,有助于理解影響水分分布的因素及其作用機(jī)理,并改進(jìn)宏觀模型.早期的孔隙尺度模擬方法是將真實(shí)的土壤孔隙結(jié)構(gòu)簡(jiǎn)化為由不同形狀和尺寸的毛細(xì)管道組成的毛管束[18],或簡(jiǎn)化為由不同孔隙空間和喉道組成的孔隙網(wǎng)絡(luò)[19],然后基于這些毛管束或孔隙網(wǎng)絡(luò)模型模擬孔隙管道中的流體運(yùn)動(dòng).毛管束模型和孔隙網(wǎng)絡(luò)模型為研究土壤中的物質(zhì)輸運(yùn)和化學(xué)過程提供了很大幫助[20-22],但其無法反映土壤的真實(shí)孔隙結(jié)構(gòu)和連通特征,因此其在模擬土壤結(jié)構(gòu)和水力特征方面往往存在較大誤差[23].

顯微成像技術(shù),尤其是X 射線斷層掃描技術(shù)(XCT)的發(fā)展,使直接觀測(cè)土壤中微納米尺度物質(zhì)及孔隙分布成為了可能[24].但目前的X-CT 和圖像處理技術(shù)仍很難辨別土壤中的水和空氣,更難直接觀測(cè)土壤中的動(dòng)態(tài)過程[25-26],因此仍需要結(jié)合孔隙尺度模擬方法來研究土壤中的流體動(dòng)力學(xué)問題.利用X-CT 獲取土壤的真實(shí)三維孔隙結(jié)構(gòu),結(jié)合計(jì)算流體力學(xué)(CFD)方法可直接模擬流體物質(zhì)在孔隙中的運(yùn)動(dòng).直接孔隙尺度模擬中常用的CFD 方法主要包括基于粒子的格子玻爾茲曼方法[27-28](LBM)和光滑粒子動(dòng)力學(xué)方法[29](SPH),以及基于網(wǎng)格的水平集方法[30](Level Set)、流體體積法[31](VOF)以及有限體積法[32](FVM)等.這些CFD 方法在模擬孔道中流動(dòng)問題方面各有優(yōu)劣勢(shì)[33],其中LBM 因其在處理復(fù)雜邊界及多相流方面具有天然優(yōu)勢(shì),并且方便并行化計(jì)算,近年來得到了迅速發(fā)展[34-37].目前應(yīng)用的LBM多相流基本模型主要有三種：顏色函數(shù)模型[38-39]、自由能模型[40-41]和Shan-Chen 偽勢(shì)模型[42-44],本文選用偽勢(shì)模型進(jìn)行模擬研究.偽勢(shì)模型通過引入偽勢(shì)來描述流體粒子之間的相互作用,具有較好的物理意義,并且易于程序?qū)崿F(xiàn),因此得到了較廣泛的應(yīng)用.但該模型缺少熱力學(xué)一致性,在模擬大密度比的多相流時(shí)容易造成數(shù)值不穩(wěn)定,雖然通過改進(jìn)狀態(tài)方程和力項(xiàng),使模型適用的密度比有了很大的提升[43],但是在模擬真實(shí)的復(fù)雜孔隙結(jié)構(gòu)中多相流時(shí)仍非常困難,尤其是對(duì)于孔隙直徑分布廣且微孔較多的孔隙結(jié)構(gòu).

土壤作為強(qiáng)非均勻性的多孔三相系統(tǒng),具有從納米尺度至毫米尺度豐富的孔隙結(jié)構(gòu),而土壤中的物理及生化反應(yīng)都發(fā)生在孔隙壁面,即液?固或氣?固界面處.接觸角反映了土壤固相的潤濕性,控制著土壤孔隙中水的狀態(tài)及分布.最近研究表明植物根系分泌物[45]、微生物活動(dòng)產(chǎn)物[46]及有機(jī)物[47]會(huì)改變土壤顆粒表面的潤濕性,并進(jìn)而影響土壤孔隙中的水分分布.流體與固體壁面的接觸角通常都采用實(shí)驗(yàn)方法測(cè)得,但是對(duì)于結(jié)構(gòu)復(fù)雜的多孔介質(zhì),通過實(shí)驗(yàn)觀測(cè)流體在孔隙內(nèi)的流動(dòng)非常困難.利用LBM 孔隙尺度模擬流體接觸角是一種簡(jiǎn)便可行的辦法,自Briant[48]于2004 年提出利用LBM 來模擬多相流體的接觸線運(yùn)動(dòng),此后LBM 被經(jīng)常用于流體接觸角研究,尤其是模擬液滴在材料表面的行為.張博[49]將LBM 應(yīng)用于材料領(lǐng)域,分析了液滴潤濕行為與材料表面微納結(jié)構(gòu)的關(guān)系.方可寧等[50]利用LBM 對(duì)加熱基板上的液滴進(jìn)行模擬,研究了不同接觸角液滴的鋪展及蒸發(fā)過程.胡夢(mèng)丹等[51]采用LBM 對(duì)結(jié)構(gòu)表面液滴的冷凝行為進(jìn)行模擬研究,分析了材料表面的幾何尺寸和接觸角的局部不均勻性對(duì)冷凝液滴形核位置及最終潤濕狀態(tài)的影響規(guī)律.這些研究都是通過LBM 方法模擬液滴在材料表面或簡(jiǎn)單管道中的運(yùn)動(dòng)狀態(tài),鮮見針對(duì)真實(shí)三維孔隙結(jié)構(gòu)中多相流的模擬,更無對(duì)土壤中孔隙水受接觸角影響的模擬研究.

本文采用改進(jìn)的Shan-Chen 格子玻爾茲曼模型,直接模擬了不同接觸角情況下水分在土壤孔隙中的分布,探討了接觸角對(duì)孔隙水分分布狀態(tài)的影響,以期為進(jìn)一步模擬土壤中流體運(yùn)動(dòng)及物質(zhì)遷移提供基礎(chǔ),并為土壤水資源研究及優(yōu)化農(nóng)業(yè)生產(chǎn)提供理論指導(dǎo).

1 格子玻爾茲曼方法

在格子玻爾茲曼模型中,流體的運(yùn)動(dòng)被描述成一系列離散的單密度分布函數(shù).標(biāo)準(zhǔn)的格子玻爾茲曼方程可表示為[52]

其中fi(x,t)為流體粒子密度分布函數(shù),i表示速度離散方向;為平衡態(tài)分布函數(shù);τ 為無量綱松弛時(shí)間;c=?x/?t,為流體粒子的格子速度;?x和?t分別為格子尺寸和格子時(shí)間步長,一般取1.F(x,t)為流體粒子受到的作用力,一般包括流流作用力Fint(x,t)、流固作用力Fs(x,t)以及其他外體力Fb(x,t).由于土壤中水受到的基質(zhì)吸力遠(yuǎn)大于重力,因此本文研究中忽略了重力影響.

本文采用D3Q19 格式,即將粒子速度在空間上按照19 個(gè)方向進(jìn)行分解,如圖1 所示.離散后i方向的速度ei可表示為[53]

ρ 和u分別表示格點(diǎn)處的宏觀密度和速度;wi為權(quán)重系數(shù),在D3Q19 格子系統(tǒng)中取[53]

流體密度ρ 和速度u可通過流體粒子分布函數(shù)的一階矩和二階矩求得

圖1 D3Q19 LBM 速度離散模型示意圖Fig.1 Schematic sketch of LBM discrete velocities in the D3Q19 scheme

1.1 流流作用力

在Shan-Chen 模型中,流體粒子之間的相互作用力可表示為[55]

其中ψ 為偽勢(shì)函數(shù),或稱有效質(zhì)量,是流體密度的函數(shù);G為格林函數(shù),表征流流粒子間相互作用力強(qiáng)度,控制液、氣相的分離及液氣界面的表面張力,G<0表示流體粒子相互吸引.對(duì)于D3Q19 模型,格林函數(shù)為[56]

g為流流作用力強(qiáng)度系數(shù).x位置處的流體粒子受到周邊流體粒子作用力的合力為[55]

1.2 流固作用力

流體粒子與固相的流固作用力用下式表示[55]

式中g(shù)s為流固作用系數(shù),表征流體粒子與固相之間的作用力強(qiáng)度.在有固相的LBM 中,g與gs共同決定液固的本征接觸角,但在本文,由于對(duì)模型進(jìn)行了改進(jìn),采用的Peng-Robinson(P-R)狀態(tài)方程,最終接觸角僅依賴于gs.s(x+ei)為標(biāo)示函數(shù),當(dāng)(x+ei)為固相時(shí),s(x+ei)=1;當(dāng)(x+ei)為流體時(shí),s(x+ei)=0;ψ(ρw)取固定值1.

1.3 模型改進(jìn)

在Shan-Chen 模型中,Fint(x)??c0ψ(x)g?ψ(x)對(duì)應(yīng)的狀態(tài)方程為,偽勢(shì)函數(shù)可表示為

其中c0為常數(shù),對(duì)于D3Q19 模型,c0=6[53].

很多研究已經(jīng)表明,原始的Shan-Chen 模型不適用于密度比大的多相流體.而實(shí)際上,無論是水和水蒸氣還是油與氣,其密度比都達(dá)到好幾百甚至上千.采用接近真實(shí)流體的狀態(tài)方程可以大幅度提高LBM適用的密度比,本文采用P-R 狀態(tài)方程,改進(jìn)后模型實(shí)現(xiàn)的最大密度比超過7000.

P-R 狀態(tài)方程為

α(T)={1+(0.374 64+1.542 26ω ?0.269 92ω2)×[1 ?(T/Tc)0.5]}2,ω 為偏心因子,對(duì)于水,ω=0.344;T為溫度,Tc為臨界溫度;a=,b=0.077 8RTc/pc.本文計(jì)算中取a=1/49,b=2/21[43].

為減小模型計(jì)算中產(chǎn)生的偽流,提高穩(wěn)定性,將流體粒子間作用力改寫為

在引入作用力時(shí)采用精確差分法(EDM),此時(shí),外力項(xiàng)為

其中u通過式(6)求得,而u′=u+F?t/ρ.EDM 方法可以提高模型的熱穩(wěn)定性,結(jié)合以上方法,該LBM可模擬的最低溫度達(dá)0.4Tc.

2 材料及土樣

本文選取兩個(gè)多孔材料的X-CT 圖像作為模擬對(duì)象,如圖2 所示,其中模擬土樣是由微小玻璃珠粘結(jié)而成,孔隙結(jié)構(gòu)相對(duì)比較均勻,尺寸為150×150×150 像素(voxel),像素分辨率為40 μm,孔隙率為0.368 1,后文用土樣1 表示;真實(shí)土樣為2015 年取自英國洛桑實(shí)驗(yàn)站的Highfield 草地試驗(yàn)田(51.810 3N,?0.374 8E),該試驗(yàn)田為粉質(zhì)黏壤土,其基本理化性質(zhì)見文獻(xiàn)[58]表1 中的Grassland,后文用土樣2 表示.土樣2 的提取過程如下：將田間提取的土壤樣芯振碎,然后分別過孔徑為4 mm,2 mm 和0.71 mm 的篩子,隨機(jī)挑選0.71～2 mm 的團(tuán)聚體,再將挑選的團(tuán)聚體通過Phoenix Nanotom 掃描系統(tǒng)進(jìn)行X-CT 掃描,掃描電壓和電流分別為90 kV 和65μA,分辨率為1.51μm.拍攝曝光時(shí)間為500 ms,整個(gè)土樣掃描時(shí)間為69 min,共獲取1440 張投影圖片.利用掃描系統(tǒng)配套的軟件Phoenix datos|x2 將投影圖片重構(gòu)為三維圖像,并導(dǎo)出三維圖像的原始連續(xù)切片圖,然后利用開源軟件ImageJ (https://imagej.nih.gov/ij/)將圖片進(jìn)行二值化處理,土樣2 的詳細(xì)獲取及處理過程見文獻(xiàn)[58].由于獲得的土樣初始圖像較大,為減少計(jì)算量,截取200×200×200 像素的立方體作為本文研究的土樣2(如圖2),其孔隙率為0.232 5.為便于分析,本文中所有物理量包括長度、密度等均直接采用格子單位,其都可以轉(zhuǎn)換成實(shí)際的物理單位,1 lu=1 voxel(lu 為格子單位).

圖2 兩個(gè)土樣的三維結(jié)構(gòu):左邊是模擬土樣(土樣1),右邊是真實(shí)土樣[58] (土樣2)Fig.2 Three dimensional structure of two soils used in the simulations,the artificial soil shows in left column and real soil shows in right

假定水只在連通孔隙中流動(dòng),為減少模擬中的計(jì)算量和所需內(nèi)存,計(jì)算前先提取連通孔隙,并采用形態(tài)模型[59]計(jì)算每個(gè)孔隙的直徑.兩個(gè)土樣的連通孔隙結(jié)構(gòu)分別如圖2 所示,土樣1 的平均孔徑為6.12像素,土樣2 的平均孔徑為5.68 像素.由圖可以看出兩個(gè)土樣的孔隙結(jié)構(gòu)有明顯差異,其中土樣1 連通孔隙率為0.368,孔徑分布為正態(tài)分布,而土樣2 連通孔隙率為0.221,孔徑分布近似泊松分布.

3 計(jì)算結(jié)果及討論

3.1 接觸角計(jì)算

研究表明大部分土壤既不是完全親水的,也不是疏水的,即接觸角在0?～90?范圍內(nèi)變化[47,60].在微生物或植物根系的作用下,土壤與水的接觸角會(huì)發(fā)生變化,但仍屬于親水物質(zhì),因此本文只研究接觸角小于等于90?的情況.

假定式(10)中ψ(ρ)=1,改變gs的值,即可獲得不同的接觸角.為獲取接觸角與gs之間的關(guān)系,先模擬一個(gè)水滴在固體平面上的穩(wěn)定狀態(tài).如圖3(a)所示,構(gòu)建一長寬高(XYZ)分別為200 lu×200 lu×100 lu 的空間,下邊界設(shè)置2 lu 厚的固體平面,其他邊界均設(shè)置為周期邊界,固體壁面采用反彈邊界條件.系統(tǒng)恒溫,溫度設(shè)置為T=0.7Tc,Tc為臨界溫度.初始時(shí)刻在固體平面中央放置一長寬相等有一定高度的方體液滴,密度為ρl0;其余空間為水蒸氣,密度為ρv0,如圖3(a).每1000 次迭代計(jì)算輸出一次結(jié)果,當(dāng)兩次輸出的液態(tài)水體積分?jǐn)?shù)變化小于10?7時(shí),則認(rèn)為水的空間分布已經(jīng)達(dá)到穩(wěn)定狀態(tài).設(shè)液態(tài)水與水蒸氣的臨界密度ρc=(ρl0+ρv0)/2,圖3(b)和圖3(c)為穩(wěn)定后的液態(tài)水在平面上的形態(tài).

液滴達(dá)到穩(wěn)定狀態(tài)后,測(cè)量液滴的高度h、液滴與固體平面接觸面直徑d,則接觸角θ 可通過下式計(jì)算得出

圖3 水滴在固體平面上的狀態(tài)變化Fig.3 The state of liquid water on a solid plane

計(jì)算得到的接觸角θ 與gs值如圖4 所示,其關(guān)系可用直線θ=?114gs+90?近似表示,假設(shè)土壤中水的接觸角θ 與gs也服從此規(guī)律,則后文中接觸角θ均可由gs依據(jù)此公式求得.土壤孔隙中的流體接觸角不僅與固液氣相互作用及液體表面張力有關(guān),還受孔隙直徑、壁面粗糙度的影響,但對(duì)于同一土樣,其固體組分及內(nèi)部孔隙結(jié)構(gòu)不變,流體接觸角僅隨gs而改變,因此以上假設(shè)不會(huì)影響本文的研究.

圖4 T=0.7Tc時(shí)液滴與光滑固體平面接觸角θ 與gs的關(guān)系Fig.4 Relationship between contact angle θ of droplet on smooth solid plane and gsunder temperature T=0.7Tc

3.2 接觸角對(duì)土壤水分布的影響

假定土壤孔隙中的水分布由毛管壓力p∝cos θ·σ/d控制,其中σ 為氣水表面張力,θ 為接觸角,d為孔隙直徑.對(duì)于給定的接觸角和表面張力,每個(gè)孔都對(duì)應(yīng)一毛管壓力.假定土樣的底面與水液壓相連,并施加吸力p,在毛管壓力與吸力p的共同作用下,水就會(huì)進(jìn)入土樣中毛管壓力小于p且連通的孔隙[59].依據(jù)上述方法對(duì)土樣1 設(shè)置初始水分分布,然后模擬不同θ 情況下水分在土樣中的運(yùn)動(dòng)及重分布,土樣1 的初始飽和度約為0.35.模擬中土樣的所有邊界均設(shè)置為周期邊界,固體壁面采用反彈邊界條件.溫度T=0.7Tc,每迭代計(jì)算1000 次輸出一次結(jié)果,當(dāng)兩次輸出的液態(tài)水體積分?jǐn)?shù)變化小于10?7時(shí),則認(rèn)為此時(shí)刻土樣孔隙中的水分分布已經(jīng)達(dá)到穩(wěn)定狀態(tài).土樣1 在不同θ 情況下的水分布如圖5 所示.在不考慮重力的情況下,θ=90.0?時(shí),液態(tài)水和水蒸氣分散分布于土壤的所有孔隙,與孔隙直徑無關(guān).隨著θ 減小,液態(tài)水的分布逐漸轉(zhuǎn)移到固體土顆粒附近,緊密圍繞固體土顆粒,形成水膜.

圖5 不同θ 情況下土樣1 孔隙中液態(tài)水及水蒸氣的分布:紅色為固體土顆粒,綠色為液態(tài)水,紫色為水蒸氣Fig.5 Distribution of liquid and vapor water in soil 1 under different contact angles θ,with solid in red,liquid in green and vapor in purple

圖6 展示了不同θ,水蒸氣在土樣1 中的分布.當(dāng)θ=90.0?時(shí),水蒸氣的分布圖像棱角尖銳,即水蒸氣在大孔隙與小孔隙中均有分布,說明接觸角接近90.0?時(shí),土壤固體壁面的表面自由能較低,對(duì)液態(tài)水的吸引力較弱;此時(shí)土壤中的毛管吸力幾乎為零,孔隙中水的基質(zhì)勢(shì)為零.隨著θ 減小,水蒸氣的分布圖像棱角逐漸變得圓滑,即水蒸氣逐漸被擠出小孔隙,說明隨著接觸角減小,固體壁面對(duì)液態(tài)水會(huì)產(chǎn)生較大的吸引力,導(dǎo)致液態(tài)水會(huì)迅速填充孔徑較小的孔隙,并排擠出氣態(tài)水,使水蒸氣在大孔隙中匯集形成氣泡;土壤中的毛管吸力隨接觸角減小而增大,此時(shí)孔隙水由于土壤固相的吸附產(chǎn)生負(fù)的基質(zhì)勢(shì),并且其值隨接觸角減小而減小.

圖6 不同θ 時(shí)水蒸氣在土樣1 中的分布Fig.6 Spatial distribution of vapor(red)in simulated soil under different contact angles θ;the voxels of solid and liquid were made transparent

土壤中的物理及生化過程都發(fā)生在界面處,因此在研究土壤中水分布時(shí)需要進(jìn)一步分析流體與固體以及流體與流體之間的界面面積.接觸角會(huì)影響液態(tài)水與空氣的分布,也必然會(huì)影響液固、氣固及液氣界面面積.為此對(duì)土樣1 中不同接觸角情況下的液固(LS)、氣固(VS)及液氣(LV)界面面積(A)進(jìn)行統(tǒng)計(jì),為研究方便,對(duì)所有界面面積A均除以土樣的固相總表面積A0進(jìn)行無量綱化.結(jié)果如圖7 所示,LS和LV 的面積均隨θ 增大而減小,而VS 的面積隨θ增大而增大.

圖7 土樣1 中液固(LS)、氣固(VS)及液氣(LV)界面面積(A/A0)隨θ 的變化Fig.7 Changes of interface area(A/A0)in soil 1 with contact angle θ

圖8 土樣2 的驗(yàn)證結(jié)果Fig.8 Verification results of soil 2

為驗(yàn)證以上規(guī)律對(duì)其他土樣是否也成立,選用一個(gè)孔隙結(jié)構(gòu)與土樣1 相差較大且非均勻性強(qiáng)的土樣2 (如圖2),按照前面的方法模擬飽和度約為0.60時(shí),不同θ 情況下土樣中液態(tài)水及水蒸氣的分布.結(jié)果與土樣1 一致,如圖8 所示,在土樣2 中,隨著接觸角減小,液態(tài)水進(jìn)入貼近土壤固相表面的小孔隙,并沿著固相表面擴(kuò)展形成水膜;氣態(tài)水逐漸被擠入大孔隙形成氣泡;液態(tài)水的表面面積隨θ 減小而增大,而水蒸氣與固體的界面面積隨θ 減小而減小.由圖7與圖8(b)對(duì)比發(fā)現(xiàn),雖然兩個(gè)土樣中LS,LV 及VS 隨θ 變化規(guī)律相同,但土樣1 中的曲線變化明顯比土樣2 的更顯著,這是由于土樣1 飽和度較低,在θ 較大時(shí),大部分小孔隙仍未被液態(tài)水填充,LS 界面面積小于VS 界面面積;隨著θ 減小,液態(tài)水進(jìn)入小孔隙,導(dǎo)致LS 界面面積顯著增大,而VS 界面面積減小.土樣2 中孔隙結(jié)構(gòu)非均勻性很強(qiáng),有非常多的小孔隙,但由于土樣2 的飽和度較高,大部分小孔隙在θ 較大時(shí)就已被液態(tài)水填充,因此LS 及VS 界面面積隨θ 變化較小.

3.3 接觸角對(duì)物質(zhì)輸運(yùn)通道直徑的影響

由式(17)可知接觸角由液滴在固體表面的厚度(h)和接觸面直徑(d)決定,其中d影響液固和氣固界面面積;而h決定液態(tài)水和水蒸氣的厚度,在三維孔隙結(jié)構(gòu)中即為液態(tài)水和水蒸氣占據(jù)的孔隙直徑,后文中直接簡(jiǎn)稱液態(tài)水或水蒸氣厚度,用hl和hv表示.若將土壤中液態(tài)水看作水及溶質(zhì)運(yùn)移的通道,而將水蒸氣看作氣體輸運(yùn)的通道,則hl和hv實(shí)際分別是土壤中液體和氣體的輸運(yùn)通道直徑.采用前面提到的計(jì)算孔隙直徑的方法,分別計(jì)算hv和hl.在不同θ,土樣1 中液態(tài)水的平均厚度ˉhl均隨θ 減小而減小,而ˉhv則均隨θ 減小先增大后減小,如圖9(a)所示.這是由于隨著θ 減小,水優(yōu)先填充小孔隙,積聚在大孔隙中的液態(tài)水沿著土壤固相表面擴(kuò)展,形成水膜,液態(tài)水的厚度hl減小,hv則相應(yīng)增大;但是當(dāng)θ 減小到一定程度后,液態(tài)水充滿所有小孔隙,并且逐漸浸潤了大孔隙的固相表面,使原來在大孔隙中與土壤固相接觸的水蒸氣脫離固相表面,形成游離在液態(tài)水中的氣泡,hv減小.對(duì)于土樣1,ˉhl和ˉhv均可以用公式ˉh=αθ2+βθ+A進(jìn)行擬合,兩個(gè)土樣的擬合結(jié)果如圖9(a)中曲線,擬合參數(shù)值如圖9 所示,參數(shù)α,β 值與土壤的孔隙結(jié)構(gòu)及飽和度有關(guān).為驗(yàn)證以上規(guī)律與土壤孔隙結(jié)構(gòu)無關(guān),又計(jì)算了土樣2 中不同gs的hl和hv,結(jié)果如圖9(b)所示.在土樣2 中,hl和hv隨θ 的變化規(guī)律與土樣1 相同,ˉhl和ˉhv也可以擬合為二次多項(xiàng)式.

3.4 水分密度分布

圖9 液態(tài)水及水蒸氣通道平均直徑隨θ 的變化Fig.9 Changes of channel diameter for liquid water and vapor with θ

土壤中同時(shí)存在液態(tài)水和水蒸氣,在干旱半干旱地區(qū),土壤中水蒸氣的作用尤為重要,水蒸氣的密度直接關(guān)系著土壤中水分含量.為分析接觸角對(duì)水蒸氣密度的影響,將不同θ 下土樣1 中水分密度分布進(jìn)行統(tǒng)計(jì),結(jié)果如表1.在土樣1 中,緊貼土壤固相表面的液態(tài)水密度最大.在飽和度一定的情況下,隨著接觸角減小,固體對(duì)流體的吸引力增強(qiáng),使得土粒表面的液態(tài)水由自由水向結(jié)合水轉(zhuǎn)變,流體勢(shì)能降低.最靠近土粒表面的為強(qiáng)結(jié)合水,其密度ρlmax隨θ減小而增大,但增大幅度較小;由于吸引力增強(qiáng),靠近土粒表面的一部分水蒸氣會(huì)被土粒吸附變成液態(tài)水,導(dǎo)致水蒸氣密度ρv隨θ 減小顯著減小,因此兩相流體的密度比ρlmax/ρv隨θ 減小而迅速增大.這說明接觸角對(duì)土壤孔隙中水蒸氣的密度有很大影響.側(cè)面反映土壤中微生物和植物根系可以通過分泌物提高附近土壤固相的潤濕性,減小水在土壤中的接觸角,吸附孔隙中的水蒸氣,增大液態(tài)水的含量,導(dǎo)致土壤中水蒸氣含量比周邊土壤低,從而達(dá)到吸取周邊水分的目的.為驗(yàn)證以上分析的正確性,另外計(jì)算了不同θ 下土樣2 中液態(tài)水與氣態(tài)水的密度,結(jié)果如表1所示,其規(guī)律與土樣1 一致.

為進(jìn)一步分析土壤中水分密度分布隨飽和度的變化,對(duì)兩個(gè)土樣在不同飽和度的水分密度分布進(jìn)行了統(tǒng)計(jì),表2 列出了θ=90.0?和θ=10.2?時(shí)土樣2的水分密度分布情況.由表2 可以看出,飽和度對(duì)土壤中液態(tài)水的最大密度并無影響.θ=90.0?時(shí),飽和度對(duì)水蒸氣密度無明顯影響,但θ=10.2?時(shí),飽和度的增大會(huì)顯著提高水蒸氣的密度,降低土壤中液相與氣相密度之比.土樣1 中水分密度隨飽和度的變化規(guī)律與土樣2 相同(文中未列出詳細(xì)數(shù)據(jù)).這一現(xiàn)象有重要的現(xiàn)實(shí)意義,在干旱半干旱地區(qū),土壤的潤濕性較弱,表層土壤中的飽和度一般都非常低,水主要以水蒸氣的形態(tài)存在,此時(shí)深層土中的水蒸氣密度與表層土壤基本相同.當(dāng)表層土壤中的水蒸氣進(jìn)入大氣或者被植物、微生物吸收時(shí),深層土壤的水蒸氣又會(huì)補(bǔ)充上來,這為表層植物和微生物的生存提供了水源,但也造成大量水資源的蒸發(fā)浪費(fèi).提高表層土壤的潤濕性,減小接觸角,一方面可以降低表層土壤的水蒸氣含量,減少地表水蒸發(fā),同時(shí)可為植物和微生物提供更多的液態(tài)水.但目前尚缺乏這方面的研究,將在后續(xù)研究中進(jìn)一步論證.

表1 不同θ 兩個(gè)土樣中水分密度分布Table 1 Density of liquid and vapor water in two soils under different θ

表2 不同飽和度下土樣2 中水分密度分布Table 2 Density of liquid and vapor water in soil 2 under different saturations

4 結(jié)論

本文采用改進(jìn)的Shan-Chen 格子玻爾茲曼模型結(jié)合X-CT,模擬了兩個(gè)不同孔隙結(jié)構(gòu)土樣中不同潤濕性條件下水分的分布,分析了接觸角對(duì)孔隙水分分布狀態(tài)的影響,初步得到了以下結(jié)論：

(1)接觸角較大時(shí),孔隙直徑對(duì)液態(tài)水和水蒸氣分布影響較小;隨著接觸角減小,液態(tài)水進(jìn)入靠近土壤固相的小孔隙,并沿著固相表面擴(kuò)展形成水膜,氣態(tài)水逐漸被擠入大孔隙形成氣泡;接觸角影響液態(tài)水與空氣分布的同時(shí),也會(huì)影響液、氣、固之間的界面面積.

(2)隨著接觸角減小,液體輸運(yùn)通道的平均直徑減小;而氣體輸運(yùn)通道平均直徑先增大,后減小.

(3)土壤中液態(tài)水的密度隨接觸角變化很小,水蒸氣密度隨接觸角減小顯著降低.接觸角較大時(shí),飽和度對(duì)水蒸氣密度無明顯影響;接觸角較小時(shí),飽和度的增大會(huì)顯著提高水蒸氣密度.