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

摘要：基于計算流體動力學(xué)與離散元（CFD-DEM）耦合的數(shù)值模擬方法，結(jié)合土壤入滲試驗(yàn)，進(jìn)行土壤滲透系數(shù)隨土層深度變化的研究。試驗(yàn)結(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)行驗(yà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]，該試驗(yàn)定律的提出為研究土壤水分入滲理論和入滲過程模型奠定了基礎(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ǔ)，描述入滲速率隨滲透時間的變化趨勢，但該模型存在模擬值與實(shí)測值吻合度不高的問題。在眾多入滲模型中，Kostiakov入滲模型（Kostiakov，1932）的使用最為簡單、便捷，因此在實(shí)際研究中被廣泛應(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 滲透率變化的試驗(yàn)

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。

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

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

滲透系數(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。

試驗(yàn)裝置（圖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ù)測定試驗(yàn)結(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ù)的變化最為明顯。試驗(yàn)還設(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ù)土柱透水性滲透試驗(yàn)，選取9 cm高度土柱為模擬對象，將土柱均勻劃分，以1 cm為尺度標(biāo)準(zhǔn)將土柱等距剖分，每隔1 cm設(shè)置1個觀測點(diǎn)（圖4）。初始時間步長為0.001 s，最小時間步長為0.001 s，最大時間步長為10 s，迭代控制參數(shù)采用默認(rèn)值。

3.2 初始條件和邊界條件

入滲過程中土壤水的運(yùn)移軌跡分為水平和垂直方向，本研究建立模型時只考慮水分在垂直面內(nèi)的運(yùn)動。土壤滲透試驗(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ù)和無滑移邊界條件，動量和湍流動能采用二階迎風(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í)際試驗(yàn)條件，設(shè)置高度為11 cm的壓力水頭，選取直徑為10 cm，高度為9 cm的土柱試件作為仿真模擬對象，進(jìn)行CFD-DEM耦合數(shù)值模擬分析，通過流體力學(xué)軟件FLUENT得到CFD-DEM耦合場下的流體動態(tài)壓強(qiáng)和固相土壤顆?？偰芰俊?/p>

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

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

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

如圖7、圖8所示，結(jié)合顆粒模型可以發(fā)現(xiàn)，流體在多孔顆粒域整體所占體積極小，且大部分位置的體積分?jǐn)?shù)都小于0.187，同時可以發(fā)現(xiàn)，上部的壓實(shí)程度較小，流體相體積分?jǐn)?shù)較大，而下部壓實(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 模型驗(yàn)證

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

式中：ST=∑ni=1（y1-y），為總平方和;SA為回歸平方和;SE為殘差平方和，相關(guān)系數(shù)為完全擬合的度量值，反映了響應(yīng)面與試驗(yàn)數(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)行滲流試驗(yà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ù)值分析與試驗(yàn)[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.