基于FLUENT土壤養分滲流擴散的數值模擬研究

趙崤隆,初燕芳,李云川,王 靜

(1.云南農業大學 資源與環境學院,云南 昆明 650201;2.云南省高校城鄉水安全與節水減排重點實驗室,云南 昆明 650201;3.云南農業大學 水利學院,云南 昆明 650201)

土壤作為以疏松狀態存在于陸地表面并具有肥力的動態自然體,其養分的運動行為以質流和滲流為主,養分隨下滲水流下滲融入土壤是植物根系拉伸生長和吸收營養物質的前提。在降雨量大、強降雨集中的濕潤、半濕潤地區,土壤養分傳輸、運移作用明顯,這一方面使植物更容易吸收養分,另一方面沒有被植物吸收的一部分養分會在運動過程中流失，導致土地質量和土壤肥力下降。降雨入滲是養分滲流的重要環節,從滲透機理出發研究土壤的入滲規律和特性,將所得到的滲透曲線應用于所建立的養分運移模型來模擬養分滲流特征,可以為養分有效利用、減少養分流失提供科學依據。目前對于土壤養分滲流規律的研究方法主要是試驗和數值模擬相結合的方法,例如：Viotti等建立了一維非飽和土壤中溶質遷移方程,研究了輸入不同參數情況下溶質濃度的分布特點[1];張家炳等建立了水分與溶質運移的數學模型，并用該模型模擬了入滲過程中離子的運移規律[2-3];Kumar等通過無網格法分析了二維非飽和多孔介質中溶質遷移問題[4];陶亞奇等采用Voronoi圖分裂法來構造與真實土壤相仿的虛擬土壤,并輸入到隨機行走模型中對土壤溶質遷移進行數值模擬，從而分析了土壤質地對溶質遷移的影響[5-6];孫露露等利用Richard方程建立溶質遷移模型進行數值模擬研究，得出了入滲速度增強可以促進溶質對流和彌散的結論[7-8];魏新平等研究了土壤養分隨水分運移過程中水流強度與入滲速度及溶質擴散的關系[9-11]。土壤顆粒在一般情況下處于離散狀態(尤其是表層土壤),水分溶液在通過土壤顆粒時,滲流規律的數值分析可以通過有限元軟件來進行。本研究模擬了降雨裝置試驗測定的土壤滲透系數隨土壤深度的變化曲線,將Forchhcimer公式修改后得出了適用于描述入滲速率隨土壤孔隙度變化的公式,并借助CFD-PDE耦合,采用離散元軟件FLUENT的流體模擬功能進行二次開發，探究了磷酸二氫鉀溶液在離散的多孔介質土壤中滲流速度和有效擴散系數隨土層深度的變化規律。

1 土壤滲透系數的測定

1.1 滲透系數的物理模型

Richard方程能夠模擬入滲過程中反應時間與含水率以及導水率的關系。含水量是整個測試過程中最為主要的測試變量,所以該數學模型能夠反映整個滲流過程的滲流規律[12-13]。

Richard方程如下:

(1)

式(1)中:θ為土壤體積含水率(cm3/cm3);ψm為非飽和土壤總水勢;D(θ)為非飽和土壤擴散率(cm3/cm3);K(θ)為非飽和土壤導水率(cm3/min);t為時間(min)。

1.2 滲透系數測試的機理

滲透系數測定及堵塞裝置的工作原理:在試驗裝置上安裝了2個電子水壓力傳感器和1個超聲波流速器,所有傳感器通過模數轉換器連接電腦。通過水壓力傳感器可測得試件上、下表面的水頭損失,同時通過流速傳感器測出水管內水的流速,可將其導入公式(2)。在滲透試驗過程中應注意避免土壤試件的側壁滲漏,且必須在線連續記錄滲透系數的變化過程,模擬土壤在一段時間內的滲流過程并且記錄動態滲透率變化。滲透系數的測試裝置見圖1。

圖1 滲透系數測試裝置以及各組成部件

(2)

式(2)中:k為滲透系數(mm/s);v1為土柱內部的平均流速(mm/s);i為水頭梯度;v2為出水管內出口的平均流速(mm/s);Aeff為土柱試件的有效截面積(mm2);Aou為出水管的有效截面積(mm2);L為試件長度(mm); △h為水頭損失(mm)。

1.3 滲透率試驗設計

為了研究土壤在不同水壓下土壤水分入滲的規律，設置3個壓力水頭(110、160、210 mm)。供試土壤為云南省昆明市梁王山鮮花小鎮的紅壤土。首先制作符合試驗設備的規格為50 mm×100 mm、100 mm×100 mm、150 mm×100 mm的有機玻璃磨具,并在磨具內部涂抹凡士林(起阻止側壁滲流以及增加粘性的作用，使土壤顆粒與壁面更加緊密，防止取樣時因有機玻璃內壁光滑而導致試件分離),將模具插入土壤形成試件,然后去掉周圍的土壤取出試件(見圖2和圖3)。取樣日期為2019年12月28日。

圖2 供試土壤樣品試件

圖3 試驗整體流程

1.4 滲透系數測定實驗結果

如圖4所示，3種顏色的曲線分別代表50 mm×100 mm、100 mm×100 mm、150 mm×100 mm試件的滲透系數,滲透系數分別由0.008、0.011、0.0075 mm/s逐漸降低；這些滲透系數的變化具有相同的變化特征,反映了土壤滲透系數隨入滲時間增加而變慢的趨勢,而50 mm高度的土柱試件因其土層厚度較薄,其土壤滲透系數的變化最為顯著。實驗還設置了不同高度的壓力水頭,而在不同高度壓力水頭下土壤滲透系數隨時間的變化大致相同。

入滲速率隨土層深度的不同而發生變化,對直徑為100 mm,高度為150 mm的試件進行固相和液相的耦合仿真模擬,結果如圖5所示，最大入滲速率為0.09 cm/s,最小入滲速率為0.35 cm/s,土壤的入滲速率隨土層深度增加而降低。

1.5 分析與校正模型

將孔隙平均滲流速度乘以孔隙度可計算出土壤的滲流速度；通過Forchhcimer公式可以將滲流速度與滲透系數聯系起來。由于在初始階段土壤中含有一定量的空氣,所以對Forchhcimer公式中的相關參數進行校正時無法做到準確定位,并且其滲透系數處于一個變化過程,因此很難準確地分析具體的滲透系數的轉換,因此可以引進一個關于非飽和度的參數S,使得在整個滲流過程中存在一個非飽和度S，以此對滲透系數進行校正，從而得到了以下基于Darcy-Forchhcimer多孔滲流的滲透系數變化的公式，亦即改進的Forchhcimer公式:

(3)

(4)

式(3)中:J為壓力梯度(Pa/m);S為非飽和度(%); μ為流體動力粘滯系數1.005;V為滲流速度(m/s);β為非達西系數(m-1)。

式(4)中：K為滲透率(10-10m2);k為水力傳導系數(m/s)；μ為流體動力粘滯系數1.005。

修正后Forchhcimer公式成為一個與孔隙度有關的函數,從而使之更加符合土壤模型的適用性。將滲透曲線與改進后的Forchhcimer公式結合,能反映整體滲流過程中孔隙體積的變化規律。

圖4 土壤滲透系數的變化曲線

圖5 入滲速率隨土層深度的變化曲線

2 養分運移的離散元土壤模型設置

土壤是典型的松散介質,利用離散元模型模擬土壤水化后聚團的形成及破裂是目前較為合理的方法,目前常用離散元法模擬土壤的模型有整合延遲彈性模型(hysteretic spring contact model, HSCM)和線性內聚力模型(liner cohesion model, LCM)。最近提出的EEPA建立云南紅土農田土壤模型最為合適,其彈性塑性形變模擬能反映顆粒之間在受力作用下更為緊密的現象,能夠表現顆粒之間孔隙壓縮的過程和反映不同壓實程度下的養分滲透機理。本實驗選取的梁王山鮮花小鎮的表層原狀土在長時間形成過程中受到壓實固化以及沉積作用，其顆粒粘結較為緊密,因此在進行實驗設計時需要采用DEM中模擬抗壓的設計進行一定的壓實，進而確保模擬的試件與實驗試件保持相同的密度[14-17]。

2.1 顆粒參數的設置

由于在模擬時不好掌握原狀土層顆粒的分布情況,因此本實驗先將取到的原狀土層壓碎，然后篩分土壤顆粒粒徑,選取2.0、1.0、0.5 mm三種孔徑的篩子進行篩分,其質量占比分別為20%、65%、15%。設置紅黏土顆粒的相關離散元參數如下:恢復系數0.55,靜摩擦系數0.20,動摩擦系數0.10,密度2600 kg/m3,剪切模量1107 Pa,泊松比0.25。

2.2 數學模型

組分運移方程如下:

(5)

式(5)中:U為流速;ρ為密度;Q為某種流體的屬性,可以是標量(如濃度C)，也可以是向量(如流速U)。

2.3 物理模型的邊界條件和網格劃分

在滲流過程中,養分滲流擴散的運動分為水平運動和垂直運動。本文以多孔介質土壤中的養分滲流擴散為研究對象來建立模型,沿X軸正負方向長度為250 mm,沿Z軸負方向的深度為250 mm,在模型中距離Z軸0 mm和80 mm處各設置1條監測線。設計邊界條件為:上界面為速度進口,在以坐標原點為中心處設40 mm×40 mm區域,該區域含有10%磷酸二氫鉀養分釋放點；下界面為速度出口,設為自由出入邊界,與Z軸對稱的左右側邊界設為壓力遠場,結構如圖6所示。

根據土體試件模型內部的特性,將各計算域劃分為四面體結構化和非結構化網格并進行分塊。為了能夠獲得更加精確的數據和觀察流場的特性,對質量分數網格進行細化,將試件以2.5 cm為標準尺度進行等距離剖分,初試的迭代控制參數設為默認值,計算區域的網格的質量百分數為95%以上，其中流體網格和介質網格劃分的結果如圖7所示。

圖6 養分滲流擴散物理模型

圖7 擴散模型的網格劃分

3 數值模擬及結果分析

多孔介質擴散是指氣體或液體進入固態物質孔隙的擴散。液態分子在多孔介質土壤中擴散是一種常見的現象,是由于濃度梯度的存在而發生的質量轉移過程。數值模擬以長500 mm、寬500 mm、高250 mm的土體為研究對象,在中心點處和距離中心點80 mm處設置兩個監測點。結合公式(3)～公式(5)通過計算流體力學CFD軟件中的FLUENT進行液固兩相的耦合仿真,以模擬磷酸二氫鉀溶液滲流過程速度、質量分數和有效擴散系數的變化特征。

3.1 模擬流場內流體速度的變化特征

在溶液下滲擴散過程中兩個監測點的速度變化如圖8所示,兩個監測點的擴散速度隨土層深度增加總體上呈下降趨勢,監測點1的擴散速度總體上處于0.00500～0.00506 m/s,在Z軸方向-0.02 m處達到峰值，這是因為從釋放點自由落體到達上界面處存在始速度，其比擴散介質的阻力大;監測點2距離中心擴散源0.08 m,由于多孔介質的阻力項的存在,其質量分數擴散速度整體小于監測點1且呈緩慢下降趨勢;模擬場中磷酸二氫鉀溶液的滲流速度在中心點最大,在土壤0.02 m深處達到峰值,并隨土層深度增加向下呈輻散式遞減。

圖8 兩個監測點滲流速度隨土層深度的變化

3.2 模擬流場內磷酸二氫鉀溶液相對質量分數的變化

相對質量分數的變化可以反映溶液濃度在滲流過程中的變化狀態。兩個監測點磷酸二氫鉀溶液相對質量分數的變化如圖9、圖10所示,其中監測點1的相對質量分數隨土壤深度增加呈下降趨勢,在Z軸方向0～0.03 m處下降幅度明顯,之后下降幅度趨于平緩;監測點2的相對質量分數隨土壤深度增加呈上升趨勢,上升幅度在Z軸方向0.05～0.15 m處最明顯。結合圖11磷酸二氫鉀溶液相對質量分數云圖得其濃度在中部高,并沿土層的縱深方向和橫向遞減。

圖9 監測點1相對質量分數隨土層深度的變化

圖10 監測點2相對質量分數隨土層深度的變化

圖11 模擬場內不同時間相對質量分數擴散云圖

3.3 模擬流場內磷酸二氫鉀溶液有效擴散系數的變化

有效擴散系數表示流體在多孔固體中的擴散狀態。從圖12可以看出，磷酸二氫鉀溶液在監測點1的有效擴散系數先增加后減少,總體處于0.001～0.0012 m2/s，在土壤0.04 m深處達到峰值,之后有效擴散系數逐漸減少。由圖13可見,監測點2的有效擴散系數整體上隨土壤深度的增加呈上升趨勢，在Z軸方向10～20 cm處上升幅度大。結合圖14的云圖得模擬場內中心點土層下0.02～0.03 m處擴散系數最大,向下和向兩側逐漸減小,最后呈彌散狀態。

4 結論與討論

同高度壓力水頭下土壤滲透系數隨時間的變化大致相同,隨時間增加呈下降趨勢;在滲流過程中入滲速率隨土層深度增加而降低,土壤滲透系數會出現動態下降的過程。擬合出與土壤滲透系數變化趨勢一致的滲流曲線,并引入非飽和度參數S,得到修正后的能反映整體滲流過程中孔隙體積變化規律的Forchhcimer公式。

圖12 監測點1有效擴散系數隨土層深度的變化

修正的Forchhcimer公式與FLUENT結合模擬液相的養分溶液填充土壤孔隙的滲流過程,通過滲流速度、相對質量分數以及有效擴散系數的變化曲線和液相分布的數值模擬流場云圖發現,多孔介質土壤中的養分在中心擴散源0.2～0.3 m處的滲流擴散能力最強,且隨著土層深度和距中心點橫向距離的增加呈輻散性減弱。多孔介質土壤作為阻力項的存在導致養分擴散系數在兩個監測點的變化不同。

圖13 監測點2有效擴散系數隨土層深度的變化

圖14 模擬場內不同時間有效擴散系數變化云圖

本實驗并沒有考慮到氣相對土壤顆粒中養分滲流狀態的影響,且開發的求解器只適用于固液兩相物體,因此本課題組今后將繼續開發相關的適用于固液氣三項數值分析的求解器和模型，以便探明氣相對土壤顆粒中養分滲流狀態的影響，從而更加深入地了解養分滲流過程的變化狀態。