非飽和礦堆溶液滲流遲滯與毛細擴散行為表征

王雷鳴，李希雯，尹升華，周根茂，李 輝，劉培正，鄧博納

1) 金屬礦山高效開采與安全教育部重點實驗室，北京 100083 2) 北京科技大學土木與資源工程學院，北京 100083 3) 金屬礦山安全與健康國家重點實驗室，馬鞍山 243000 4) 武漢工程大學綠色化工過程教育部重點實驗室，武漢 430205 5) 核工業(yè)北京化工冶金研究院，北京 101149

堆浸技術憑借其基建周期短、投入成本低、作業(yè)安全性高等優(yōu)勢，被廣泛應用于低品位硫化銅礦、砂巖型金礦等戰(zhàn)略礦產資源回采作業(yè)[1–2].其中，堆浸是一種氣、固、液多相介質并存、滲流場、溫度場、化學場等多場耦合的復雜反應體系[3]；該體系中，溶液是可溶性氧、Cu2+/Cu+、Fe2+/Fe3+等有價金屬離子等溶質賦存的重要媒介，直接影響著反應傳質過程與浸礦效率[4].因此，能否有效表征堆內溶液毛細擴散、流動與分布規(guī)律是制約礦物浸出效率的重要因素.作為一種非飽和多孔介質中多相流，溶浸液的滯后現(xiàn)象是普遍存在且十分復雜的[5].已有研究發(fā)現(xiàn)堆內溶液流動主要有兩種形式[6–7]：重力驅動的優(yōu)勢流和微觀力驅動的毛細擴散流.其中，由毛細管力驅動的毛細擴散過程，對溶液滲流遲滯行為的影響是十分顯著的[8].毛細擴散過程不僅是影響礦堆滲流遲滯行為，更重要的會影響傳質過程與浸礦效率，需要深入理解和認識.

目前，針對非飽和堆內溶液滲流遲滯行為規(guī)律，特別是毛細擴散過程，國內外專家學者主要從物理實驗、無損探測與數(shù)值模擬三個層面展開.其中，英國帝國理工大學Ilankoon等[9–10]開展礦堆溶液滲流遲滯實驗，利用持液率、殘余穩(wěn)態(tài)持液率對非飽和礦石顆粒堆的持液行為進行了研究；王雷鳴等[11]立足于制粒礦堆，基于自主研發(fā)的持液行為原位監(jiān)測系統(tǒng)實現(xiàn)了對不可動液、可動液的量化表征；毛細擴散過程對非飽和堆（如破碎礦堆、制粒礦堆等）持液行為的影響不可忽略[12]，毛細入滲率與毛細速率、含水率呈正相關[13–14]；隨著多種無擾動探測設備技術的進步與引入，對溶液滲流的研究逐漸步入了微細觀層面[15–16]，英國劍橋大學Fagan等[17]利用核磁共振技術（Magnetic resonance imaging, MRI）對堆內溶液擴散和溶液優(yōu)勢流形成過程進行了研究，薛振林等[18]利用粒子圖像測速技術（Particle image velocity, PIV）對礦堆結構內溶液滲流場進行了無損探測，有效實現(xiàn)了飽和溶液滲流場的可視化等.然而，單純的實驗研究與無擾動探測存在實驗誤差大、瞬態(tài)短時、隨機性強等劣勢，難以實現(xiàn)對溶液滲流過程進行有效的預測.對此，多類數(shù)學表征模型與滲流模擬軟件逐漸被提出應用：國外 Bouffardt和 Dixon D G[19]、Liu和Hashemzadeh[20]合作研發(fā)了堆浸模擬程序HeapSIM，可以對工業(yè)堆浸噴淋、礦物浸出實現(xiàn)過程監(jiān)測與結果預測；我國Wu[21]和尹升華等[22]利用COMSOL multiphysics軟件和計算機斷層掃描（Computed tomography, CT）構建非飽和礦堆網格模型并模擬了溶液滲流過程，證實了堆內優(yōu)勢流動區(qū)與溶液停滯區(qū)；此外，Ma等[23]、Zulian等[24]聚焦裂隙巖體滲流，利用Fluent開展了數(shù)值模擬研究，實現(xiàn)了破裂砂巖內的溶液流動規(guī)律.

綜合已有研究不難發(fā)現(xiàn)，以往研究通常利用短時、瞬態(tài)可視化表征手段、數(shù)值模擬軟件與相似物理滲流實驗等方式開展研究，多聚焦溶液優(yōu)勢流動、重力流動、飽和條件、橫向毛細擴散等方面[25–27]，對非飽和礦堆體系的上向毛細滲流擴散、滲流遲滯行為的認識仍不深刻和不清晰，嚴重制約著浸礦過程中溶液滲流調控方法與措施，相關研究有待進一步深入開展.對此，加拿大Topp等提出了時域反射方法，該法利用電磁波在不同介質中的傳播速度的差異來測定土壤含水率的一種方法，具有良好的可預測性[28–29]，因此可以被有效借鑒到本研究當中.

對此，本文聚焦非飽和礦石顆粒堆，開展毛細溶液滲流遲滯與毛細擴散行為表征模擬，構建非飽和礦堆滲流模型，利用時域反射器（Time domain reflector）實時探測并計算非飽和礦堆內持液率，并結合Design expert進行響應曲面分析，深入探究非飽和礦堆內溶液滲流過程中各影響因素之間的相互關系.運用 COMSOL multiphysics開展溶液毛細滲透過程模擬，相關研究結果對調控礦堆內溶液滲流、強化溶液毛細擴散，進而提高礦物浸出效率具有良好的指導意義.

1 非飽和礦堆溶液滲流模型構建

1.1 基本假設

為研究溶液在浸出礦堆中的滲流規(guī)律，本文構建礦石顆粒堆簡化模型，將其簡化為矩形二維平面進行分析，建立尺寸為 0.05 m × 0.2 m 的數(shù)學模型.并且，做出如下基本假設：

（1） 礦石顆粒的骨架形狀、顆粒相對位置不隨時間發(fā)生變化；

（2） 顆粒間孔隙率大小、形狀、連通性不隨時間發(fā)生變化；

（3） 滲流過程中流體不可壓縮；

（4） 礦石顆粒堆內孔隙率、滲透率各向同性；

（5） 礦石顆粒堆的孔隙空間中存在氣液兩相流動狀態(tài)；

（6） 氣液兩相流動均為低速低擾動的達西滲流；

（7） 模擬環(huán)境溫度不變.

1.2 控制方程

為探明礦堆非飽和浸出滲流規(guī)律，本文以界面作用為切入點，分析浸出液滲流毛細擴散狀態(tài)[30].在礦堆滲透過程中，基于達西滲流定律、質量守恒定律和兩相流理論[31]，引入毛細吸力和體積含水率的變量，分別推導礦堆非飽和滲流過程中液相和氣相的運移方程，如式（1）和（2）所示：

式中：C為容水度，反映了非飽和多孔介質的儲存特性，Cw和Ca分別為液相和氣相容水度；pw和pa分別 為孔隙水壓力和孔隙氣壓力；ρw和ρa分別為液相和氣相的密度，kg·m?3，κ為飽和滲透率，%；kr,w和kr,a分別為液相和氣相的相對滲透率，%；η為運動黏度，Pa·s；g為重力加速度，m·s–2；D為垂直坐標.

目前，礦堆水分特征曲線目前具有代表性的模型較多，運用較為廣泛的模型是非飽和滲流van Genuchten-Mualem（VGM）模型[32]，為求解函數(shù)毛細壓力Pc，根據(jù)持液率和滲透率關系，分別定義θ、C、Se和kr，如式（3）～（6）所示：

式中，Se為飽和度，1；Hc為毛細壓力水頭，m；θw為持液率，%；其中，θw的變化范圍從殘余持液率θr增大到飽和持液率θs.α，n，m和L是 Van Genuchten參數(shù)，由固相介質本身決定.

迭代后，可求得氣相的VGM公式：

1.3 物理模型及網格劃分

為探明礦堆非飽和浸出滲流過程中毛細上升規(guī)律，定義模型為經過對稱軸的任意截面.模型的寬為 5 cm，高為 20 cm；設置礦堆底存在高 3 cm、水平方向無窮大的恒定飽和流動區(qū).

利 用 COMSOL multiphysics 5.0 數(shù) 值 模 擬 軟件，構建二維網格模型，如圖1所示.由圖1可見，較礦堆中部而言，礦堆頂部和底部的邊界網格結構更為細密，其可更好的表征溶液毛細擴散規(guī)律.

圖1 柱浸顆粒堆物理模型構建及網格劃分Fig.1 Physical model and mesh of the packed heap in column leaching

借鑒工業(yè)礦堆滲流的實際情況及相關文獻數(shù)據(jù)，選取模型主要參數(shù)，如表1所示.

表1 數(shù)學模型的關鍵參數(shù)Table 1 Key parameters of the mathematical model

1.4 邊界條件

建模前對模型邊界條件進行設定.初始狀態(tài)，礦堆中的氣液是符合靜壓分布的.而液相從底部向上運動，氣相是在礦堆頂部出現(xiàn)，且不會消失.假定氣液兩相流動均為低速低擾動的達西滲流，在入口處、表面、邊壁，氣液兩相滿足式（11）.

在底部，液相和氣相的壓力為恒壓，如式（12）所示：

式中，p為液體t時的總壓強，p0液面上方的壓強，t為時間.

1.5 數(shù)值模擬方案與賦值

為進一步探究溶液在浸出礦堆內的溶液滲流遲滯規(guī)律，將物理實驗結果和數(shù)值模擬相結合，利用孔隙率對滲透系數(shù)進行賦值，考察不同孔隙率、噴淋強度對非飽和礦堆滲流過程；探索滲流速度、毛細吸力、持液率、非飽和系數(shù)等相關特性.具體模擬和實驗方案，如表2和表3所示.

表2 不同噴淋強度和孔隙率條件下溶液毛細滲流模擬方案Table 2 Experimental scheme of the liquid capillarity seepage under different irrigation rate and porosity condition

表3 不同孔隙率條件下溶液毛細滲流模擬方案Table 3 Experimental scheme of the liquid capillarity seepage under different porosity condition

依據(jù)初始實驗條件差異，各組互為對照，對實驗方案進行簡要闡釋，如下：（1）不同噴淋強度與孔隙率條件下，共設置6個實驗組，噴淋強度分3 個階梯，即：0、10、50 L·m–2·h–1，礦石顆粒堆孔隙結構用孔隙比進行表征，選擇2種，即：1.04、1.17.（2）不同孔隙率條件下，不進行噴淋作業(yè)，即僅依賴毛細吸力和下部飽和溶液區(qū)供給顆粒堆，孔隙比介于1.04～1.17，孔隙率介于51%～54%，相應地，礦石顆粒堆的水力傳導系數(shù)介于0.02～0.08 cm·s–1.

2 非飽和堆溶液毛細滲流物理實驗

流體在毛細管中流動時，由于各相的張力不同，相界面會產生毛細管吸力，毛細上升過程的實質是氣液兩相流的過程.為進一步揭示非飽和堆溶液毛細滲流規(guī)律，利用時域反射器（Time domain reflector， TDR），基于式（13），對非飽和堆內持液率（或稱含水率）進行原位實時探測和計算，揭示不同噴淋強度、孔隙比等條件下非飽和堆溶液毛細滲流實驗.

式中，θTDR為校正的持液率，%；為平均持液率，%；θ1～4為4個測量點位的持液率，%.溶液滲流物理實驗裝置，如圖2所示.該毛細滲流實驗研究裝置是由噴淋系統(tǒng)（1）、柱體（2）、非飽和礦堆（3）、恒液位槽（4）、多孔隔篩（5）、支撐結構（6）、TDR（7）、數(shù)顯系統(tǒng)（8）、溢流口（9）、集液罐（10）共10部分組成.實驗開始前，將非飽和礦石堆的底部浸泡在恒定液面的溶液區(qū)域，使之在礦堆最底部形成一個飽和溶液區(qū).實驗開始后，溶液在毛細吸力作用下，溶液自下而上發(fā)生毛細擴散，形成浸潤面，并不斷向上遷移擴散，隨著堆內溶液毛細擴散過程趨于停滯，非飽和礦石顆粒堆逐步達到穩(wěn)態(tài)持液.該過程中，利用毛細吸力等關鍵參數(shù)，實現(xiàn)溶液毛細擴散過程的有效監(jiān)測.

圖2 非飽和礦堆溶液毛細上升實驗裝置構成.(a)毛細擴散實驗實物圖; (b)時域反射儀; (c)實驗裝置結構構成Fig.2 Composition of the experimental device for the capillary rise of the unsaturated ore pile solution: (a) macroscale image of the capillary diffusion experiment; (b) time-domain reflector; (c) detailed structure of the experimental device

3 結果與討論

3.1 持液率隨浸潤時間變化規(guī)律

持液率是非飽和滲流中最具代表性的參數(shù)之一，探明持液率變化是揭示非飽和礦堆內部滲流遲滯行為規(guī)律的重要前提.其中，噴淋強度和礦石粒徑均對礦堆持液率有重要影響[33–34].對此，本文取圖1模型頂部中點為研究對象，坐標為（0.025,0.2），基于前述構建的毛細擴散數(shù)學模型，開展對液體流動過程中持液率變化過程數(shù)值模擬，生成不同孔隙比和噴淋強度條件下非飽和堆內持液率曲線，如圖3所示.

圖3 不同噴淋強度和孔隙比條件下持液率隨時間變化Fig.3 Changes of the liquid holdup with time under different irrigation and porosity ratio conditions

結果表明，隨著溶液噴淋作業(yè)的進行，礦堆持液率的收斂性逐漸增加.礦堆持液率與破碎礦石顆粒的幾何平均粒徑呈負相關、與孔隙率呈正相關、與噴淋強度呈正相關.由圖3可知：（1）對比A1-A2、A3-A4和A5-A6可知，不考慮溶液噴淋強度影響時，礦堆持液率與孔隙比、水力傳導系數(shù)呈正相關；以t= 20 s為例，具有更高孔隙比（1.174）和水力傳導系數(shù)（0.08 cm·s–1）的 A2 組持液率（約為 52%）明顯高于 A1組持液率（約為 30%）.（2）對比A1-A3-A5、A2-A4-A6可知，不考慮水力傳導系數(shù)和孔隙比的影響時，礦堆持液率與噴淋強度呈正相關；以t= 10 s為例，具有更高噴淋強度的 A1組持液率（45%）要高于A3組持液率（44%）和A5組持液率（41%）.并且，在噴淋初期（0～20 s），上述因素（噴淋強度、水力傳導系數(shù)和孔隙比）對礦堆持液率的影響更為顯著.

分析認為：礦石顆粒堆持液率對礦石粒徑（孔隙率）的敏感度要明顯高于對噴淋強度的敏感度；礦石粒徑越大，溶液噴淋強度對礦堆含水率的影響越小[35].在相同孔隙比條件下，噴淋強度小的模型持液率在初期增長緩慢、后期增速較快.結果表明：在溶液噴淋初期顆粒堆內毛細吸力極高，噴淋溶液可以自下部飽和溶液區(qū)發(fā)生上向位移，形成浸潤鋒，快速進入顆粒堆內部，此時，固–液接觸主要以潤濕作用為主.當噴淋強度越小時，模型內部的溶液量較少，宏觀表現(xiàn)為持液率較低；反之，噴淋強度較大時，宏觀礦堆持液率較高.在溶液噴淋后期，堆內原有顆粒間、顆粒內孔隙被溶液大量填充，重力場和滲流場共同影響溶液毛細滲流過程[36]，堆內溶液量趨于穩(wěn)態(tài)，宏觀表現(xiàn)為持液率增速減慢.

3.2 堆內不同位置處毛細吸力與持液率關聯(lián)規(guī)律

忽略頂部噴淋作用的影響，本文探究了不同孔隙率條件下礦石顆粒堆毛細滲透行為規(guī)律，獲取毛細吸力–持液率特征曲線，結果如圖4所示.由圖4可見，在相同孔隙結構（孔隙率）條件下，礦堆頂部的毛細吸力要明顯高于礦堆底部的毛細吸力，這種毛細吸力的差異在物理模型的中部位置更為顯著.這是由于礦石顆粒堆底部與飽和溶液接觸，該處的持液率最大且毛細吸力最小.以孔隙率52.0%為例，礦堆底部穩(wěn)態(tài)毛細吸力約為1669.5 Pa，明顯高于礦堆頂部的穩(wěn)態(tài)毛細吸力（139.0 Pa）.

圖4 非飽和堆毛細吸力特征–孔隙率關聯(lián)關系.(a) 堆頂毛細吸力與孔隙率間的關系; (b) 毛細吸力差值與孔隙率間的關系Fig.4 Relationship of capillarity suction features and porosity of ore heap: (a) relationship between capillarity force (top of the column) and porosity;(b) relationship between capillarity force differences and porosity

無論是顆粒堆頂部還是底部，毛細吸力均與孔隙率的變化呈負相關，即：礦堆該處的毛細吸力隨孔隙率的增加而減少.以礦堆頂部為例，當?shù)V堆孔隙率由51%增加至52.5%時，堆頂毛細吸力由1670.1 Pa 下降至 1669.4 Pa.此外，堆頂?shù)拿毼ψ兓潭让黠@大于堆底.由圖4（a）可見，對比孔隙率為51%和54%的礦堆，可見二者堆底毛細吸力的凈差值為59.189 Pa，明顯高于堆頂毛細吸力的凈差值（0.725 Pa）.分析認為：液體總是從吸力小往吸力大的方向運動，模型頂部和底部的毛細吸力之差的物理意義為液相自下而上的運動的驅動力[37]，毛細吸力克服重力和摩擦阻力做功并趨于平衡.圖4（b）表明，溶液上升過程中，由于孔隙率較大、持液率差較大，液體自下而上的運動的驅動力也就越大.雖然孔隙率越大，礦堆中的毛細驅動力越大，但由于液體自重和摩擦阻力的共同作用，溶液上升速度反而越小.結合已有研究與工業(yè)現(xiàn)場運行情況可知，這種溶液的滲流遲滯行為對于礦石顆粒潤濕和礦堆持液率的影響是十分明顯的，礦堆礦石內孔隙結構的發(fā)育程度越高，礦堆內的溶液滲流遲滯行為越顯著，礦石顆粒的潤濕效果更好.

3.3 礦石–溶液特征曲線

假定礦石顆粒堆內孔隙與大氣聯(lián)通，且不考慮摩擦阻力和蒸發(fā)作用，那么可知：礦石顆粒堆內壓力水頭與毛細水頭的值大小相等但方向相反.基于此假設，本文取模型頂部中點為研究對象，繪制礦石顆粒堆持液率–水頭壓力（藍線）和毛細吸力（紅線）關系曲線，如圖5所示.

由圖5可知，在相同持液率條件下，壓力水頭與毛細水頭的對稱特征被有效的證實.隨著噴淋時間的增長，非飽和礦石顆粒堆持液率的收斂性增加.對于孔隙比小的礦堆，需要更長的時間達到穩(wěn)態(tài)持液.對于孔隙比相同條件時，噴淋強度主要影響初期持液率的值；噴淋強度越大，初期持液率增長越快.當持液率為0%時，毛細吸力水頭達到峰值，約為2×104Pa；隨著礦堆內溶液量的逐漸增加，持液率明顯上升，相應地，堆內毛細吸力逐漸減小.當?shù)V石顆粒堆持液率達到穩(wěn)態(tài)持液率時，堆內毛細吸力下降至0 Pa.對比不同噴淋強度和孔隙比條件下持液行為差異可知，當持液率相同條件下，礦堆孔隙較為發(fā)育時毛細吸力較小；換言之，毛細吸力的變化對孔隙率較小礦石顆粒堆更敏感.

圖5 持液率與水頭壓力（藍線）和毛細吸力（紅線）關系.(a)A1 組；(b)A2 組；(c)A3 組；(b)A4 組；(c)A5 組；(b)A6 組Fig.5 Relationship of the liquid holdup and the pressure heap (blue) and capillarity suction (red): (a) A1 group; (b) A2 group; (c) A3 group; (d) A4 group: (e) A5 group; (f) A6 group

3.4 非飽和堆內毛細吸力、孔隙比與噴淋強度響應規(guī)律

為進一步探討堆內毛細吸力、孔隙比和噴淋強度間的關聯(lián)響應規(guī)律，利用Design expert進行響應曲面分析.當t= 400 s時，堆內毛細滲透過程達到穩(wěn)態(tài)，用模擬軟件計算不同噴淋強度和孔隙率條件下的毛細吸力，分別獲得頂部毛細吸力、底部毛細吸力、毛細吸力差與孔隙比和噴淋強度的響應規(guī)律，如圖6所示.

圖6 毛細吸力對噴淋強度和孔隙率的等值面特征.(a) 堆頂; (b) 堆底; (c) 毛細吸力差值Fig.6 Equivalent surface characterization of the capillarity suction to the irrigation intensity and porosity: (a) top of column; (b) bottom of column; (c)differences of capillarity force

結果表明：噴淋強度、孔隙比和毛細力間存在明顯的響應關聯(lián)關系，噴淋強度較小時，頂部毛細吸力的值基本上不會隨孔隙率發(fā)生變化.數(shù)值模型底部的毛細吸力受噴淋強度和孔隙比的影響較為明顯.當噴淋強度越大（50 L·m–2·h–1）時，孔隙比越?。?.04），底部毛細吸力越大；反之，當噴淋強度越?。?0 L·m–2·h–1），孔隙比越大（1.17）時，自下而上的毛細吸力越大，也越不利于內部液體的滲透擴散.

為進一步考察堆內滲透特征與毛細吸力間的內在關聯(lián)，本文利用構建的溶液毛細滲流模型與COMSOL multiphysics數(shù)值模擬軟件，獲取了沿模型垂直方向的相對滲透率分布規(guī)律，由圖7所示.

圖7 基于 COMSOL multiphysics 的礦堆內毛細吸力穩(wěn)態(tài)分布特征.(a) 毛細吸力; (b) 相對滲透率Fig.7 Steady distribution characterization of the capillarity suction in the ore heap relied on COMSOL Multiphysics: (a) capillarity forces; (b) relative permeability

結果表明：相對滲透率kr隨著高度的增加呈現(xiàn)非線性遞減特征，隨高度增加，kr加速減小，在模型頂部達到最小值.表明礦堆達到穩(wěn)態(tài)持液時，礦堆上部的滲透率低，下部的滲透率高.礦石顆粒堆的滲透特性主要由頂部滲透率決定的.結合礦堆穩(wěn)態(tài)毛細吸力、相對滲透系數(shù)的計算結果發(fā)現(xiàn)：在孔隙率越小和噴淋強度較大時，礦堆更易達到穩(wěn)態(tài)持液狀態(tài)，并且，底部毛細吸力越大、相對滲透率越大.在毛細吸力驅動下，礦石顆粒堆持液率逐步上升達到穩(wěn)態(tài)持液率.礦堆內溶液滲流擴散過程受重力和毛細吸力共同作用，兩者方向相反.孔隙率較大時，毛細吸力較大但滲透較好，總體表現(xiàn)為浸潤面擴散過程較慢.堆內溶液毛細上升速度分布不均勻，具體表現(xiàn)為：礦堆頂部毛細上升速度較慢，礦堆下部毛細上升較快.

4 結論

本次建模主要基于礦堆達西滲流和兩相流模型，借助多相多場模擬軟件COMSOL multiphysics的多物理場功能，模擬驗證礦堆非飽和滲流過程，獲得非飽和礦堆內持液率、毛細吸力、孔隙率、噴淋強度等之間的相互聯(lián)系.通過模擬分析，得出以下結論：

（1） 相較于噴淋強度，礦堆持液率對堆孔隙率更為敏感，礦堆持液率隨噴淋時間的增長收斂性增加，且孔隙率小礦堆需要更長時間達到穩(wěn)態(tài)持液；

（2） 毛細吸力的變化對孔隙率較小的礦石顆粒堆更敏感.礦堆底部毛細吸力最小、頂部最大，且孔隙率越大兩端的壓力差越大、內部溶液毛細上升速度越慢，利于溶浸液滲透；

（3） 噴淋強度較小時，孔隙比越大，毛細吸力越大，但此時重力對溶液滲流影響較大，不利于堆內溶液的滲透擴散；在噴淋強度較大時，孔隙比越小，底部毛細吸力越大，毛細吸力的影響提升，礦堆更易達到穩(wěn)態(tài)持液狀態(tài).