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

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

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

關鍵詞：土壤滲透系數;數值分析;液固耦合;CFD;DEM;數值模擬

中圖分類號：S152.7+2

文獻標志碼：A

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

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

1 滲透率變化的試驗

1.1 試件制作

供試土壤為云南省昆明市盤龍區松華壩水源區棕壤土，土壤樣品為表層土壤，采樣時剔除土壤表面植物殘渣和固體石塊，捏碎表面大顆粒黏結土塊，使得清理后的土壤表面均勻平整，將相應尺寸的有機玻璃圓柱體打入土體取樣。供試土樣基本理化性質如下：容重、全氮含量、全磷含量、全鉀含量、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，材質為有機玻璃，為圓柱形筒體，滲透率測量儀器采用透水性滲透設備。分別設置直徑為 10 cm，高度為3、6、9 cm的圓柱形土柱各3個，共計9個土柱試件，土壤試件如圖1所示。

1.2 滲透系數測定和試驗裝置

滲透系數計算公式如下：

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

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

1.3 滲透系數測定試驗結果

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

2 入滲研究的數學模型

2.1 Green-Ampt入滲物理模型

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

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

2.2 土壤水分運移方程

水分在土壤中的運動軌跡可分為水平側移和縱向垂直運移，Richard模型可以描述土壤水在一維垂直方向上的運動規律，數學方程式可表示為

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

3 數值模擬

3.1 設置網格和時間步長

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

3.2 初始條件和邊界條件

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

3.3 參數設置

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

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

4 模擬結果與分析

4.1 模擬場內流體受到的動態壓強和顆粒總能量

根據實際試驗條件，設置高度為11 cm的壓力水頭，選取直徑為10 cm，高度為9 cm的土柱試件作為仿真模擬對象，進行CFD-DEM耦合數值模擬分析，通過流體力學軟件FLUENT得到CFD-DEM耦合場下的流體動態壓強和固相土壤顆粒總能量。

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

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

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

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

4.3 數值模擬耦合場內入滲速率的變化

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

5 模型驗證

為了確保模型的準確性和適應性，需要對模型的預測能力進行評估，根據需要對響應面方程進行顯著性檢驗，一般采用相關系數R2和校正決定系數R2adj來了解其逼近程度。

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

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

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

6 結論

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

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

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

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

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

參考文獻：

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

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

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

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

[5]楊慶璐，李子涵，李洪文. 基于CFD-DEM耦合的集排式分肥裝置顆粒運動數值分析[J]. 農業機械學報，2019，50（8）：81-89.

[6]祁文燕. 紙坊溝流域坡面尺度土壤水分動態變化及其運移數值模擬[D]. 蘭州：蘭州大學，2018.

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

[8]孫丹焱，鄭 濤，徐竟成，等. 城市綠地土壤滲透性改良對雨水徑流污染的削減效果及去除規律[J]. 環境工程學報，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]. 江蘇農業科學，2020，48（5）：280-286.

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