基于CFD-DEM耦合的土壤滲透性數(shù)值分析

李云川 趙崤隆 石平平 李靖 南軒 王靜

摘要：基于計(jì)算流體動(dòng)力學(xué)與離散元（CFD-DEM）耦合的數(shù)值模擬方法，結(jié)合土壤入滲試驗(yàn)，進(jìn)行土壤滲透系數(shù)隨土層深度變化的研究。試驗(yàn)結(jié)果表明，隨著入滲過(guò)程的持續(xù)，3、6、9 cm土柱滲透系數(shù)分別由0.008 0、0.011 0、0.007 5 mm/s逐漸降低。在液固耦合模型中，通過(guò)計(jì)算流體動(dòng)力學(xué)（CFD）設(shè)置流體相參數(shù)，離散元（DEM）模擬土壤顆粒固體相，得出滲流速率與土層深度的定量變化關(guān)系式為y=ax2+bx+c（y為滲流速率;x為土層深度;a為常數(shù)，等于 -0.078 86;b為常數(shù)，等于-0.008 67;c為常數(shù)，等于0.092 96）。通過(guò)方差分析對(duì)模型擬合的回歸性方程進(jìn)行驗(yàn)證，得到相關(guān)系數(shù)R2=0.950 3，校正決定系數(shù)R2adj=0.999 77，都接近于1，表明擬合性好。

關(guān)鍵詞：土壤滲透系數(shù);數(shù)值分析;液固耦合;CFD;DEM;數(shù)值模擬

中圖分類(lèi)號(hào)：S152.7+2

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1002-1302（2020）16-0255-05

土壤和水是人類(lèi)社會(huì)生產(chǎn)、生活必不可少的自然資源。土壤水的運(yùn)動(dòng)行為主要分為質(zhì)流和入滲，入滲是一個(gè)復(fù)雜的過(guò)程，無(wú)論是植物根系拉力導(dǎo)致的質(zhì)流作用，還是土壤中的原始水分含量，均來(lái)源于入滲過(guò)程。降水入滲是土壤水循環(huán)過(guò)程的重要環(huán)節(jié)，開(kāi)展土壤水分入滲規(guī)律和特性的研究，可以合理判定地表徑流，進(jìn)行暴雨洪水預(yù)警，為制定因土壤侵蝕導(dǎo)致的水土流失治理措施提供依據(jù)[1-3]。有關(guān)土壤水分入滲的研究可以追溯到1856年，法國(guó)水力學(xué)家達(dá)西（Darcy，1856）提出了達(dá)西定律，可反映水在巖土空隙中的滲流規(guī)律[4-7]，該試驗(yàn)定律的提出為研究土壤水分入滲理論和入滲過(guò)程模型奠定了基礎(chǔ)。Green-Ampt入滲物理模型僅適用于在薄積水或土質(zhì)均勻條件下的土壤水分入滲研究[8]，將土壤在垂直方向上劃分為入滲濕潤(rùn)層、入滲飽和層和濕潤(rùn)鋒，該模型反映了入滲速率與土壤特性之間的關(guān)系。1931年，Richard在達(dá)西定律的研究基礎(chǔ)上提出了可用于非飽和流土壤水分入滲分析的Richard方程[8-9]，并得出導(dǎo)水率與基質(zhì)勢(shì)之間的關(guān)系，用于描述土壤水分運(yùn)動(dòng)的規(guī)律，Richard入滲模型（Richard，1931）的建立是達(dá)西定律的繼承與發(fā)揚(yáng)。Philip入滲模型（Philip，1957）以Richard入滲模型為基礎(chǔ)，描述入滲速率隨滲透時(shí)間的變化趨勢(shì)，但該模型存在模擬值與實(shí)測(cè)值吻合度不高的問(wèn)題。在眾多入滲模型中，Kostiakov入滲模型（Kostiakov，1932）的使用最為簡(jiǎn)單、便捷，因此在實(shí)際研究中被廣泛應(yīng)用，但在利用該模型進(jìn)行參數(shù)率定時(shí)，各參數(shù)的物理意義表述仍有爭(zhēng)議。以上數(shù)學(xué)模型的創(chuàng)建極大地?cái)U(kuò)展了土壤水分入滲的研究，根據(jù)長(zhǎng)期的科學(xué)研究可知，土壤入滲率隨土層深度的增加而呈現(xiàn)“高- 低-穩(wěn)定”的變化趨勢(shì)，但以往的研究很少涉及土層深度對(duì)入滲速率變化影響的定量描述[10]，本研究基于計(jì)算流體動(dòng)力學(xué)與離散元（CFD-DEM）耦合的數(shù)值模擬方法，分析土壤入滲速率隨土層深度增加的變化情況，以期獲得不同深度土層土壤的滲透系數(shù)定量變化規(guī)律。

1 滲透率變化的試驗(yàn)

1.1 試件制作

供試土壤為云南省昆明市盤(pán)龍區(qū)松華壩水源區(qū)棕壤土，土壤樣品為表層土壤，采樣時(shí)剔除土壤表面植物殘?jiān)凸腆w石塊，捏碎表面大顆粒黏結(jié)土塊，使得清理后的土壤表面均勻平整，將相應(yīng)尺寸的有機(jī)玻璃圓柱體打入土體取樣。供試土樣基本理化性質(zhì)如下：容重、全氮含量、全磷含量、全鉀含量、pH值、EC（可溶性鹽濃度）值、氧化還原電位（Eh）分別為1.69 t/m3、0.14%、0.21%、0.75%、7.3、0.8 mS/cm、136 mV。

試驗(yàn)所用試件的直徑為10 cm，高度為3、6、9 cm，材質(zhì)為有機(jī)玻璃，為圓柱形筒體，滲透率測(cè)量?jī)x器采用透水性滲透設(shè)備。分別設(shè)置直徑為 10 cm，高度為3、6、9 cm的圓柱形土柱各3個(gè)，共計(jì)9個(gè)土柱試件，土壤試件如圖1所示。

1.2 滲透系數(shù)測(cè)定和試驗(yàn)裝置

滲透系數(shù)計(jì)算公式如下：

式中：K為滲透系數(shù)，mm/s;υ1為土柱試件內(nèi)部的平均流速，mm/s;i為水梯度，表示土柱試件高度與長(zhǎng)度的比值;υ2為出水管內(nèi)出口的平均流速，mm/s;S1為出水管的有效截面積，m2;S2為土柱內(nèi)的有效截面積，m2;L為土柱試件長(zhǎng)度，mm;ΔH為水頭損失，mm。

試驗(yàn)裝置（圖2）工作原理：透水性滲透測(cè)定設(shè)備上安裝了2個(gè)電子水壓力傳感器和1個(gè)超聲波流速器，壓力水頭高度為11、16、21 cm，出水口為自由出水口，傳感器、超聲波流速器與數(shù)據(jù)采集儀均連接測(cè)量設(shè)備，水壓力傳感器可測(cè)得土柱試件上下兩端水頭損失（ΔH），流速傳感器可測(cè)得土柱試件內(nèi)部流速（υ1）和出水管內(nèi)水流速（υ2），將數(shù)據(jù)采集儀收集的數(shù)據(jù)導(dǎo)入公式（1），計(jì)算得出滲透系數(shù)。

1.3 滲透系數(shù)測(cè)定試驗(yàn)結(jié)果

圖3顯示了在壓力水頭高度為11 cm條件下，土柱試件滲透系數(shù)隨時(shí)間的變化規(guī)律，3條曲線分別對(duì)應(yīng)不同高度的土柱試件，它們具有相同的線型，均反映了土壤滲透系數(shù)隨入滲時(shí)間增加而呈現(xiàn)由高到低的變化趨勢(shì)，3、6、9 cm土柱滲透系數(shù)分別由0.008 0、0.011 0、0.007 5 mm/s逐漸降低。3 cm高度的土柱試件因其較薄的土層厚度，土壤滲透系數(shù)的變化最為明顯。試驗(yàn)還設(shè)置了高度為16、21 cm 的壓力水頭，而在不同高度壓力水頭下土壤滲透系數(shù)隨時(shí)間的變化規(guī)律大致相同。

2 入滲研究的數(shù)學(xué)模型

2.1 Green-Ampt入滲物理模型

改進(jìn)前的Green-Ampt入滲模型適用條件范圍單一，僅適用于土質(zhì)均勻、入滲初期土壤水分含量低的土壤入滲過(guò)程研究，導(dǎo)致該模型的使用有很大的局限性。國(guó)內(nèi)外學(xué)者出于拓寬Green-Ampt入滲模型適用范圍的目的，對(duì)其進(jìn)行了改進(jìn)，修正后的Green-Ampt入滲模型可表示為

式中：K為滲透系數(shù)，mm/s;Ks為有效飽和導(dǎo)水率，mm/s;θ1為初始含水率，mm3/mm3;θ2為飽和含水率，mm3/mm3;I為累計(jì)入滲量，mm;SK為入滲土層濕潤(rùn)鋒處的固定吸力。

2.2 土壤水分運(yùn)移方程

水分在土壤中的運(yùn)動(dòng)軌跡可分為水平側(cè)移和縱向垂直運(yùn)移，Richard模型可以描述土壤水在一維垂直方向上的運(yùn)動(dòng)規(guī)律，數(shù)學(xué)方程式可表示為

式中：θ為土壤體積含水率，%;H為壓力水頭高度，mm;t為滲透時(shí)間，s;α為流向與垂直方向夾角，rad，該模型為一維垂直方向的一階偏導(dǎo)，因此α=0;Kf為非飽和滲透系數(shù)，mm/s;S為源匯項(xiàng)。

3 數(shù)值模擬

3.1 設(shè)置網(wǎng)格和時(shí)間步長(zhǎng)

根據(jù)土柱透水性滲透試驗(yàn)，選取9 cm高度土柱為模擬對(duì)象，將土柱均勻劃分，以1 cm為尺度標(biāo)準(zhǔn)將土柱等距剖分，每隔1 cm設(shè)置1個(gè)觀測(cè)點(diǎn)（圖4）。初始時(shí)間步長(zhǎng)為0.001 s，最小時(shí)間步長(zhǎng)為0.001 s，最大時(shí)間步長(zhǎng)為10 s，迭代控制參數(shù)采用默認(rèn)值。

3.2 初始條件和邊界條件

入滲過(guò)程中土壤水的運(yùn)移軌跡分為水平和垂直方向，本研究建立模型時(shí)只考慮水分在垂直面內(nèi)的運(yùn)動(dòng)。土壤滲透試驗(yàn)從土柱上端連續(xù)進(jìn)水，壓力水頭設(shè)定為11 cm高度定水頭，模擬模型上邊界選擇恒定壓力水頭邊界（constant pressure head），下邊界為自由出水邊界（free drainage）。

3.3 參數(shù)設(shè)置

CFD設(shè)置：液相為試驗(yàn)用水，選用RNG k-ε湍流模型進(jìn)行模擬，采用壓力進(jìn)口，自由出口，壁面條件選擇增強(qiáng)壁面函數(shù)和無(wú)滑移邊界條件，動(dòng)量和湍流動(dòng)能采用二階迎風(fēng)格式，采用壓力耦合方程組的半隱式（SIMPLE）算法進(jìn)行求解。

DEM設(shè)置：固相為土壤，進(jìn)出口條件與液相相同，模擬過(guò)程為瞬時(shí)模擬，時(shí)間步長(zhǎng)為0.001 s，模擬時(shí)間總時(shí)長(zhǎng)為5 s。

4 模擬結(jié)果與分析

4.1 模擬場(chǎng)內(nèi)流體受到的動(dòng)態(tài)壓強(qiáng)和顆粒總能量

根據(jù)實(shí)際試驗(yàn)條件，設(shè)置高度為11 cm的壓力水頭，選取直徑為10 cm，高度為9 cm的土柱試件作為仿真模擬對(duì)象，進(jìn)行CFD-DEM耦合數(shù)值模擬分析，通過(guò)流體力學(xué)軟件FLUENT得到CFD-DEM耦合場(chǎng)下的流體動(dòng)態(tài)壓強(qiáng)和固相土壤顆粒總能量。

從圖5流體受到的動(dòng)態(tài)壓強(qiáng)云圖可以看出，當(dāng)水流與土壤接觸時(shí)，流體動(dòng)態(tài)壓強(qiáng)基本保持不變，在入滲過(guò)程持續(xù)一段時(shí)間后，土壤顆粒中的空隙壓力會(huì)驟然上升，這是由于流體的持續(xù)入滲填充了原本土壤顆粒間的空隙，當(dāng)水分入滲量達(dá)到飽和時(shí)，流體所受壓強(qiáng)不再發(fā)生明顯的變化，基本保持不變。

從圖6可以看出，在滲流初期土壤顆粒總能量最大，這是因?yàn)榱黧w從高處落下時(shí)攜帶一定的勢(shì)能，在沖刷到土壤時(shí)產(chǎn)生曳力作用，導(dǎo)致土壤顆粒受力，從而增加土壤顆粒的動(dòng)能，土壤顆粒總能量也隨之增大。隨著入滲過(guò)程的持續(xù)，土壤水分含量達(dá)到飽和，顆粒總能量減少，少部分區(qū)域的土壤顆粒因土壤吸水產(chǎn)生的黏結(jié)作用而呈現(xiàn)總能量增大的變化，但總的來(lái)說(shuō)，顆粒總能量基本保持不變。

4.2 數(shù)值模擬耦合場(chǎng)中土壤顆粒相與流體相的分布

如圖7、圖8所示，結(jié)合顆粒模型可以發(fā)現(xiàn)，流體在多孔顆粒域整體所占體積極小，且大部分位置的體積分?jǐn)?shù)都小于0.187，同時(shí)可以發(fā)現(xiàn)，上部的壓實(shí)程度較小，流體相體積分?jǐn)?shù)較大，而下部壓實(shí)較嚴(yán)重，體積分?jǐn)?shù)相對(duì)較小，這與滲透系數(shù)隨深度增加而減少一致。

4.3 數(shù)值模擬耦合場(chǎng)內(nèi)入滲速率的變化

進(jìn)行數(shù)值模擬時(shí)選擇直徑為10 cm，高度為 9 cm 的土柱試件為研究對(duì)象。每隔1 cm土層深度設(shè)定1個(gè)監(jiān)測(cè)面，通過(guò)FLUENT軟件進(jìn)行液固2相耦合仿真，得出每個(gè)監(jiān)測(cè)面流體的平均流速。結(jié)合圖9、圖10可以看出，土壤模型上層孔隙度較大，監(jiān)測(cè)面所監(jiān)測(cè)的平均入滲速率大于下層。

5 模型驗(yàn)證

為了確保模型的準(zhǔn)確性和適應(yīng)性，需要對(duì)模型的預(yù)測(cè)能力進(jìn)行評(píng)估，根據(jù)需要對(duì)響應(yīng)面方程進(jìn)行顯著性檢驗(yàn)，一般采用相關(guān)系數(shù)R2和校正決定系數(shù)R2adj來(lái)了解其逼近程度。

式中：ST=∑ni=1（y1-y），為總平方和;SA為回歸平方和;SE為殘差平方和，相關(guān)系數(shù)為完全擬合的度量值，反映了響應(yīng)面與試驗(yàn)數(shù)據(jù)的符合程度，其值要在0.9以上才能保證擬合效果。

式中：校正決定系數(shù)R2adj為自變量與因變量的相關(guān)程度，其值越接近1擬合效果越好。對(duì)擬合回歸性方程進(jìn)行方差分析可以得到，滲透速率回歸方程的相關(guān)系數(shù)R2=0.950 3，校正決定系數(shù)R2adj=0.999 77，兩者都接近于1，表示方程擬合性良好，可認(rèn)為擬合得到的回歸方程能夠反映參數(shù)之間的關(guān)系，可以用來(lái)擬合滲流速度與土壤深度的變化關(guān)系。得到如下擬合方程：

式中：y為滲流速率，m/s;x為土層深度，m;a為常數(shù)項(xiàng)，等于-0.078 86;b為常數(shù)項(xiàng)，等于 -0.008 67，c為常數(shù)項(xiàng)，等于0.092 96。

6 結(jié)論

通過(guò)設(shè)置11、16、21 cm高度的壓力水頭，對(duì)直徑為10 cm，高度為3、6、9 cm土柱分別進(jìn)行滲流試驗(yàn)，分析土壤滲透系數(shù)隨土壤深度的變化規(guī)律，并基于CFD-DEM耦合進(jìn)行土壤滲透性數(shù)值模擬研究，得出結(jié)論如下：

（1）通過(guò)相分布云圖和滲透系數(shù)的曲線可知，入滲過(guò)程中土壤深度越低的區(qū)域，土壤含水量越高，土壤深度越高的區(qū)域，土壤含水量越低，這是由滲透系數(shù)隨土壤深度的變化決定的。

（2）滲透系數(shù)受壓力水頭高度的影響較小，在11、16、21 cm高度的壓力水頭下，滲透系數(shù)無(wú)明顯變化。

（3）通過(guò)數(shù)值模擬分析得出入滲速率隨土層深度增加而降低的定量變化關(guān)系，可表示為y=ax2+bx+c，其中y為入滲速率，x為土層深度，a為常數(shù)項(xiàng)，等于-0.078 86;b為常數(shù)項(xiàng)，等于-0.008 67，c為常數(shù)項(xiàng)，等于0.092 96。

（4）對(duì)擬合回歸性方程進(jìn)行方差分析，得到相關(guān)系數(shù)R2=0.950 3，校正決定系數(shù)R2adj=0.999 77，二者均接近于1，表明擬合度好。

參考文獻(xiàn)：

[1]高 巖，胡曉蕾，劉 然. 雨水滲流過(guò)程中土壤水分動(dòng)態(tài)變化規(guī)律研究[J]. 建筑技術(shù)，2019，50（1）：4-7.

[2]楊 帆，武桂芝，馮增帥，等. 大沽河河床非飽和入滲模型研究[J]. 水土保持通報(bào)，2017，37（6）：152-156.

[3]王 俊，黃歲樑. 土壤水分特征曲線模型對(duì)數(shù)值模擬非飽和滲流的影響[J]. 水動(dòng)力學(xué)研究與進(jìn)展，2010，25（1）：16-22.

[4]李俊燁，蘇寧寧，胡敬磊. 基于CFD-DEM耦合的磨粒流微小孔加工數(shù)值分析與試驗(yàn)[J]. 農(nóng)業(yè)工程學(xué)報(bào)，2018，34（16）：80-88，299.

[5]楊慶璐，李子涵，李洪文. 基于CFD-DEM耦合的集排式分肥裝置顆粒運(yùn)動(dòng)數(shù)值分析[J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào)，2019，50（8）：81-89.

[6]祁文燕. 紙坊溝流域坡面尺度土壤水分動(dòng)態(tài)變化及其運(yùn)移數(shù)值模擬[D]. 蘭州：蘭州大學(xué)，2018.

[7]丁海晶，姜 姜，張金池，等. 土壤滲透性的區(qū)域變化規(guī)律及因子分析[J]. 水土保持學(xué)報(bào)，2019，33（1）：51-56.

[8]孫丹焱，鄭 濤，徐竟成，等. 城市綠地土壤滲透性改良對(duì)雨水徑流污染的削減效果及去除規(guī)律[J]. 環(huán)境工程學(xué)報(bào)，2019，13（2）：372-380.

[9]Moret-Fernández D，Latorre B，Giner M L，et al. Estimation of the soil hydraulic properties from the transient infiltration curve measured on soils affected by water repellency[J]. Catena，2019，178：298-306.

[10]蘭簡(jiǎn)琪，謝世友. 有機(jī)復(fù)合肥對(duì)土壤水分入滲特性的影響[J]. 江蘇農(nóng)業(yè)科學(xué)，2020，48（5）：280-286.

[11]Solekhudin I. Steady infiltration in heterogeneous soil[J]. Journal of Physics（Conference Series），2018，1108（1）：012030.