基于多物理場耦合的毛細水高度研究

鄧改革, 何建國, 康寧波

(1.寧夏大學 土木與水利工程學院, 銀川 750021; 2.中國礦業大學銀川學院 信息工程學院, 銀川 750021; 3.寧夏大學 農學院, 銀川 750021)

毛細作用是指浸潤液體在細管里升高和不浸潤液體在細管里降低的現象，是液體的慣性力、黏性力、毛細力與重力等共同作用的結果[1]。土壤毛細水是靠毛細吸引力而保持于土壤毛細孔隙中的水。土壤毛細水上升高度對包括路基凍壞以及農田土壤次生鹽漬化等在內的眾多過程具有重要的影響，目前國內外學者在毛細水上升高度方面已經開展了大量的研究工作但尚不深入。

于丹等[2]進行了不同干密度重塑黃土的毛細上升速率和最大高度的研究，發現毛細最大上升高度與干密度之間非線性關系，而是呈現近似拋物線的函數關系。李先瑞等[3]針對改性黃土進行了垂直方向毛細水上升作用的研究，發現石灰能夠有效降低土體內的含水率，且隨著含量的增加，吸水作用越明顯，而水泥對于土體內部結構的改性作用更大，提升土體強度和遇水穩定性，阻礙毛細水上升作用顯著。肖紅宇等[4]進行了基于黏性土分形特征的毛細水上升高度研究，推導出土體內毛細水上升高度預測公式，并將公式預測結果與實測結果進行了對比分析，發現誤差在工程允許范圍之內，具有一定的實際應用價值。落宇杰等[5]進行的壓實黃土狀粉土毛細特性試驗發現，毛細水上升的速率隨毛細上升高度的增加而逐漸衰減。袁玉卿等[6]通過室內試驗研究了豫東黃泛區粉砂土毛細水上升規律及控制技術，研究發現隨時間的延長毛細水上升高度逐漸穩定，毛細水上升速度與壓實度呈反比，級配碎石、水泥穩定土、纖維水泥穩定土能有效阻隔毛細水上升。

胡明鑒等[1]進行了砂性土條件下，不同粒徑及潛水礦化度組合的毛細水上升試驗，研究了不同粒徑及潛水礦化度下對毛細水上升高度及上升速度的規律。王聰[7]、栗現文[8]等通過觀測毛細水上升速度和高度，研究不同濃度鹽溶液和鹽漬土對毛細水上升的影響，并對這些影響因素進行機理分析。

魏樣[9]、童玲[10]等分別研究了石油及柴油污染對土壤毛細水上升特性的影響。章求才等[11]研究了溫度和氣壓對某金屬礦山尾礦壩中毛細水上升規律的影響，發現溫度對毛細水上升的影響較氣壓大；在毛細水上升初期，溫度、氣壓對毛細水上升規律的影響不大，但在上升后期其影響逐漸明顯。

Masoodi等[12]采用達西定律模擬了液體向由纖維素和高吸水性纖維制成的紙狀膨脹多孔介質中的芯吸現象。提出了新的毛細模型，將多孔介質中達西定律與質量守恒方程相結合，通過向質量守恒方程添加附加的匯或源項，以說明基質的膨脹和液體吸收。新模型預測芯吸速率較修正的沃什伯恩方程預測芯吸速率更好。

Aghajani等[13]提出了一種改善土壤中毛細水上升高度準確性的方法，該方法重新考慮了傳統方法中被忽略的毛細區域中水分和吸力的變化，試驗結果表明，改進的方法比以前的毛細管上升解決方案更加準確和通用。

Fili等[14]進行了采用多相格子玻爾茲曼方法(LBM)檢驗土壤持水曲線中毛細現象的研究。LBM模型首先針對基準問題進行驗證，然后用于模擬靜態粒子陣列，以此來研究初始流體密度分布對系統毛細管響應的影響。結果表明，濕潤鋒的統一性直接影響土壤骨架毛細響應峰值的大小。

Zhou等[15]研究了非飽和土的抗剪強度與毛管持水曲線的關系,提出了一種新的、考慮吸附水毛細凝結的非飽和土壤保水模型。模型中，將土壤的飽和度分為基于毛細管水的飽和度和基于吸附水的飽和度。通過對部分飽和的雙圓柱系統的分析，提出了一種新的非飽和土抗剪強度準則。研究結果將大大改進對試驗數據預測的準確性。

以上對毛細水高度的研究多側重于室內土柱試驗，研究內容多集中于不同土質結構、礦化水、污染物以及外界溫度壓力等方面，取得了一定的研究成果。而基于多物理場耦合的毛細水上升高度機理方面研究尚未見公開報道。本文使用多孔介質中多相流理論，通過多物理場耦合軟件COMSOL Multiphysics建立模型，耦合多孔介質界面中達西定律和相輸運，最后通過實驗室土柱試驗對模擬結果進行驗證。研究結果將對弄清毛細水上升機理提供新的思路和方法。

1 研究方法

1.1 數值模擬研究

1.1.1 理論模型 在經典的多孔介質模型中，構成其基體的尺寸及其之間的孔隙在毛細吸水過程中保持恒定。目前計算毛細水上升高度的模型主要有兩種，一種基于Lucas-Washburn方程，另一種基于達西定律[12]。這兩種方法在計算毛細上升高度方面均存在缺陷，致使理論計算值與實際結果之間偏差較大。下面介紹多物理場耦合多孔介質中的多相流計算模型，該多物理場耦合了多孔介質界面中的達西定律和相輸運。多孔介質界面中的相輸運遵循潤濕或非潤濕流體體積分數的獨立方程：

(1)

式中：εp為孔隙率；ρi為第i相的流體密度(kg/m3)；?為梯度算子；si為相的體積分數；κ為滲透率(m2)；κri為相對滲透率(給定流體的飽和度的函數)；μi為流體的動態黏度(Pa/s)；pi為壓力(Pa)；?pi為壓力差(Pa)；g為重力加速度(m/s2)。

由于兩相的體積分數之和為1，因此可以從以下公式計算剩余的體積分數：

s1=1-s2

(2)

式中：s1為非潤濕相體積分數；s2為潤濕相體積分數。

毛管壓力pc作為潤濕相sw(模型中為s2)的飽和度和入口毛細管壓力pec的函數進行計算。通過使用Brooks-Corey模型[16-18]，毛管壓力由下式得出：

(3)

式中:λp為孔隙分布指數；pec為入口毛細管壓力。基于Brooks-Corey模型,濕相和非濕相的相對滲透率由下式給出:

(4)

(5)

式中：κrs為濕相相對滲透率；κrsn為非濕相相對滲透率。

通過達西定律接口將達西定律與連續性方程相結合：

(6)

式中:ρ為流體密度；μ為流體動力黏度；κ為滲透率(m2)。

1.1.2 模擬參數確定 以1，2，3號管土柱為數值模擬的研究對象，建立模型水力邊界條件：地下水位線在砂土柱的最底端并保持穩定，砂土柱頂端設置為自由邊界[6]。由于豎管的對稱結構，同時為了減小網格數量以提高計算效率，建立豎管的二維模型，所建立模型寬度為41 mm，高度為1 000 mm，并對所建立模型進行網格劃分，建模相關參數見表1—2。計算時間為1 600 min，時間步長為5 min，通過COMSOL Multiphysics多物理場仿真軟件，利用有限差分法對上述方程進行離散化求解。

表1 水、空氣和多孔材料的材料參數

表2 不同編號砂土多孔介質特性參數

1.2 土柱試驗研究

本文采用豎管法研究毛細水上升高度。試驗用砂土編號為#1，#2，#3，按照相關規范進行物理指標試驗，測定相應參數，結果見表3。

表3 試驗用砂土物理參數

試驗采用自制的毛細管水上升高度裝置(圖1)。裝置由儀器支架、3支有機玻璃管、有機玻璃水箱、標尺、微型循環水泵及相連接的橡膠管組成。有機玻璃管長度為100 cm，內徑41 mm，外徑45 mm，頂端敞口可以與大氣連通。將編號為#1，#2，#3的砂土分別裝入編號為1，2，3的有機玻璃管中。三根有豎管放置在有機玻璃水箱托盤中，托盤高出水箱液面一定距離，通過微型循環水泵將水箱中水補充至托盤中，使托盤中液面保持恒定。

圖1 毛細管水上升高度裝置

毛細水上升高度主要試驗步驟如下：將篩分配制好的一定質量的砂樣，通過漏斗將土樣裝入有機玻璃管中，使其具有相同的壓實度。打開微型循環水泵并開始計時，觀測濕潤鋒的位置，記錄毛細水上升高度與時間。前期毛細水上升較快，每隔10 min記錄濕潤鋒所在位置，之后隨著毛細水上升速度變慢，讀數時間間隔相應加大。

1.3 現場試驗研究

在實訓基地內部選擇合適的場地，在場地上方提前30 d布置遮擋物防止外界雨水滲入，然后開挖長寬高分別為1 m×1 m×0.8 m的土坑，通過水泵向坑底部供水，不斷調整水泵的流量使水面保持恒定，通過鋼卷尺讀取土壤濕潤鋒位置的變化見圖2。

圖2 毛細水上升高度現場試驗

2 結果與分析

2.1 數值模擬結果及分析

從圖3可知在模擬1 600 min后1，2，3號管濕潤鋒位置高度分別為16.8，23.1，37.5 cm，該模擬結果與后面實驗室室內土柱試驗結果較為接近。由圖4可知，在模擬初期毛細水高度隨時間增加較快，在模擬后期毛細水高度增加緩慢并逐漸趨于穩定。在毛細水上升高度不高的情況下，由于基質吸力做功遠大于水分重力勢，所以毛細水上升速度較快；但隨著高度的上升，重力勢增加明顯，基質吸力做的功相對減少，上升的速度就逐漸變慢[6]。

圖3 不同編號毛細管1 600 min毛細水高度

圖4 3號管7個不同時刻毛細水高度情況

2.2 土柱試驗結果及分析

由圖5可知，毛細水上升試驗開始的前120 min毛細水上升速度較快。試驗開始供水的第2小時，1，2，3號管毛細水分別上升至14.2，18，28 cm，這表明在相同壓實度條件下不同粒徑的土柱毛細水上升高度不同。3號管(粒徑范圍：0.15～0.28 mm)毛細水上升高度高于1號管(粒徑范圍：0.36～0.45 mm)及2號管(粒徑范圍：0.28～0.36 mm)，說明砂土粒徑越小毛細水上升的高度越大。隨后毛細水高度繼續增大，但是其上升速度逐漸減小。至510 min，1，2，3號管毛細水分別上升至14.6，19.6，34 cm。

圖5 毛細水前510 min上升高度

由圖6可知，在隨后的試驗時間內，試驗土柱內毛細水高度持續增加，3號管毛細水高度最大，2號管次之，1號管最小。至1 600 min 1，2，3號管毛細水分別上升至15.6，22.0，37.5 cm。毛細高度隨時間變化曲線斜率代表毛細水上升速度(圖7)，前120 min曲線斜率較大，表明毛細水上升速度較大，隨后曲線斜率逐漸減小，毛細水上升速度也隨之逐漸減小但并未停止，毛細水高度增加緩慢并逐漸趨于穩定。

圖6 毛細水1 600 min上升高度

圖7 毛細水上升速度隨時間變化(0～1 600 min)

利用Origin軟件，對圖7數據的曲線進行函數擬合，所得關系式見表4。上升高度關系式兩邊對時間求導即可得毛細水上升速度表達式(表4)。

表4 毛細水上升高度擬合關系式

由表4可知，毛細水上升高度符合冪函數性質，為增函數，其變化特征是初期高度增加快，到一定程度后增加緩慢。而毛細水上升速度屬于指數減函數，毛細水上升速度隨時間減小。

將數值模擬結果與室內試驗結果進行對比(表5)，在誤差范圍內數值模擬與室內試驗結果吻合較好。

表5 模擬高度與試驗高度對比

2.3 現場試驗結果及分析

由于現場土坑中砂土參數與3號管中砂土參數較為接近，因此圖8試驗數據與3號管試驗數據較為接近。通過分析數據可知毛細高度在30 cm以內，現場試驗毛細高度小于實驗室內土柱試驗高度，其試驗誤差小于10%，當毛細高度超過30 cm后，試驗誤差逐漸增大至15%左右。現場觀測發現，土坑中30 cm以上高度區域砂土粒徑增大，且砂土中石塊等雜物增多，對毛細水上升高度產生了影響。將試驗誤差考慮在內，可以認為現場毛細水上升規律與室內土柱試驗以及多物理場仿真結果相吻合。

圖8 毛細水1 600 min上升高度現場試驗

3 討論與結論

相同壓實度不同粒徑砂土柱內，毛細水上升高度與時長近似呈冪函數關系，為增函數，其變化特征是初期高度增加快，到一定程度后增加緩慢；毛細水上升速度與時長近似呈指數函數關系，為減函數，毛細水上升速度表現為初始階段較大但隨后迅速減小以致各組間差異不明顯[8]。

在本試驗中，相同壓實度不同粒徑砂土柱毛細水上升高度由大到小依次為：3號管(粒徑范圍0.36～0.45 mm)，2號管(粒徑范圍0.28～0.36 mm)，1號管(粒徑范圍0.15～0.28 mm)，毛細水上升高度與砂土粒徑呈負相關，總體表現為粒徑越小，毛細水上升高度越大。在較小粒徑條件下，黏性力的作用較強，在小粒徑范圍內起主要作用。相反，在大粒徑范圍下，重力作用表現更加明顯，尤其是到毛細運移的后期，重力的影響增大，作用不可忽略[19]。砂柱粒徑越小其內部毛細孔隙越小，毛細管力越大，毛細潤濕峰上移迅速，同樣時間間隔內，上升距離較高，且在試驗結束時的高度也遠遠高于大粒徑的高度[19]。這一結論和胡明鑒[1]、董斌[19]等的研究結論具有相同之處。

本文通過耦合多物理場多孔介質中達西定律和相輸運，較好地模擬了毛細水上升高度，并且在試驗誤差范圍內理論值與室內試驗結果吻合較好，該研究方法在其他參考文獻中尚未見公開報道。

壓實度作為土柱試驗的參數之一對毛細水高度具有一定的影響，本文僅研究了相同壓實度條件下不同粒徑的毛細水高度，而基于多物理場耦合條件下，不同壓實度對毛細水高度的影響尚有待今后進一步研究。

試驗用砂土作為一種多孔介質，其內部毛細流動現象是由慣性力、黏性力、毛細力與重力共同作用產生的，但在不同毛細高度各項作用力對毛細上升的貢獻并不相同，這也是進一步理解毛細流動現象機理的關鍵所在[19]。

由于現場砂土粒徑及土中雜物的影響使得現場試驗與實驗室試驗及數值模擬結果存在一定的誤差，但是將試驗誤差考慮在內，可以認為現場毛細水上升規律與室內土柱試驗以及多物理場仿真結果相吻合。該結論應用于工程領域時，應該充分考慮現場因素對毛細水高度的影響，以現場試驗作為補充，盡量減小試驗誤差。