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

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

摘要：基于計算流體動力學(xué)與離散元（CFD-DEM）耦合的數(shù)值模擬方法，結(jié)合土壤入滲試驗，進(jìn)行土壤滲透系數(shù)隨土層深度變化的研究。試驗結(jié)果表明，隨著入滲過程的持續(xù)，3、6、9 cm土柱滲透系數(shù)分別由0.008 0、0.011 0、0.007 5 mm/s逐漸降低。在液固耦合模型中，通過計算流體動力學(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）。通過方差分析對模型擬合的回歸性方程進(jìn)行驗證，得到相關(guān)系數(shù)R2=0.950 3，校正決定系數(shù)R2adj=0.999 77，都接近于1，表明擬合性好。

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

中圖分類號：S152.7+2

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

文章編號：1002-1302（2020）16-0255-05

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

1 滲透率變化的試驗

1.1 試件制作

供試土壤為云南省昆明市盤龍區(qū)松華壩水源區(qū)棕壤土，土壤樣品為表層土壤，采樣時剔除土壤表面植物殘渣和固體石塊，捏碎表面大顆粒黏結(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。

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

1.2 滲透系數(shù)測定和試驗裝置

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

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

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

1.3 滲透系數(shù)測定試驗結(jié)果

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

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

2.1 Green-Ampt入滲物理模型

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

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

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

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

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

3 數(shù)值模擬

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

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

3.2 初始條件和邊界條件

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

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

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

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

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

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

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

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

從圖6可以看出，在滲流初期土壤顆?？偰芰孔畲?，這是因為流體從高處落下時攜帶一定的勢能，在沖刷到土壤時產(chǎn)生曳力作用，導(dǎo)致土壤顆粒受力，從而增加土壤顆粒的動能，土壤顆粒總能量也隨之增大。隨著入滲過程的持續(xù)，土壤水分含量達(dá)到飽和，顆?？偰芰繙p少，少部分區(qū)域的土壤顆粒因土壤吸水產(chǎn)生的黏結(jié)作用而呈現(xiàn)總能量增大的變化，但總的來說，顆?？偰芰炕颈３植蛔?。

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

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

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

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

5 模型驗證

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

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

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

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

6 結(jié)論

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

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

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

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

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

參考文獻(xiàn)：

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

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

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

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

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

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

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

[8]孫丹焱，鄭 濤，徐竟成，等. 城市綠地土壤滲透性改良對雨水徑流污染的削減效果及去除規(guī)律[J]. 環(huán)境工程學(xué)報，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]蘭簡琪，謝世友. 有機(jī)復(fù)合肥對土壤水分入滲特性的影響[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.