基于COMSOL-PHREEQC的砂巖型鈾礦浸出性能模擬分析技術

賈明濤，金家聰，陳梅芳，蘇學斌

（1.中南大學 資源與安全工程學院，湖南 長沙 410083；2.新疆中核天山鈾業(yè)有限公司，新疆 伊寧 835000；3.中國鈾業(yè)有限公司，北京 100010）

我國是鈾礦資源儲量豐富的國家之一，其中砂巖型鈾礦約占三分之一［1-2］。原地浸出技術目前是砂巖型鈾礦的首選開采方法，具有安全性高、開采效率高、污染小等優(yōu)勢［3］。然而，原地浸出過程包含復雜的物理-化學作用，如有溶浸劑的流動、鈾礦物溶解與沉淀反應、浸出液元素的運移等，因此浸出率通常難以準確評價。

為量化評價鈾礦浸出性能，國內外學者們做了大量研究。Liddell K.C.等［4］利用H2O2-NH4HCO3-（NH4）2CO3溶液模擬瀝青鈾礦和方解石堿法地浸，強調礦物反應和溶液遷移耦合的重要性。鄭和秋野［5］通過現(xiàn)場試驗模擬鈾礦床的酸法地浸，認為試驗單元的浸出效果與礦石品位、含礦層地球化學環(huán)境相關。李春光［6］通過對砂巖鈾礦進行室內試驗模擬，認為加入表面活性劑有助于提高礦石的浸出性能。隨著計算機技術的發(fā)展，數(shù)值模擬可以克服現(xiàn)場試驗和室內試驗存在的缺點，逐漸成為原地浸出動力學研究和工程計算重要且有效的手段。曾晟等［7］基于編程語言開發(fā)模擬程序，建立多因素多過程耦合作用下的動力學模型，對鈾礦山溶質遷移機理開展研究。黃群英等［8］基于Visual Modflow 模擬和水文地球化學模式探討了某砂巖型鈾礦地浸單元中溶質運移和溶液滲流的關系。林效賓等［9］應用PHREEQC 對鈾礦床地下水系統(tǒng)中的水文地球化學過程進行模擬。Lagneau V 等［10］基于HYTEC代碼對原地浸出鈾礦山進行反應運移模擬，但沒有考慮復雜地質屬性的空間變異性。雖然前人在原地浸出采鈾過程的數(shù)學模型方面開展了一定的研究，但是主要集中在地下水動力場模擬或地球化學模擬計算，較少考慮到鈾礦浸出反應中所包含的溶解-沉淀、水相絡合、氧化還原反應等一系列復雜地球化學反應與溶浸液運移的耦合效應。

COMSOL Multiphysics（簡稱COMSOL）是一款基于有限元理論的數(shù)值分析軟件包，可模擬各種物理現(xiàn)象，并可實現(xiàn)多物理場耦合模型的仿真計算［11］。PHREEQC 是一款當前國際上通用的水文地球化學模擬軟件，可實現(xiàn)礦物溶解沉淀、水溶物配合、氧化還原、離子交換等功能的模擬［12］，且在軟件后續(xù)的迭代升級中加入了模塊化開發(fā)的支持［13］，使得數(shù)值軟件能夠調用其強大的地球化學反應模擬功能，由此拉開了COMSOL 與PHREEQC 耦合模擬研究的序幕。相關研究包括COMSOL-IPHREEQC 模式［14］、Interface COMSOLPHREEQC（iCP）［15］等。

本研究基于MATLAB 平臺開發(fā)了COMSOL Multiphysics 和PHREEQC 的耦合接口程序，構建了砂巖型鈾礦原地浸出過程中浸出液的滲流場和流體顆粒運移場，實現(xiàn)了地浸過程中的化學動力學計算，提出了適用于砂巖型鈾礦浸出性能的評價指標，對C16-1、C16-2 采區(qū)分別進行了模擬計算，數(shù)字化定量分析了砂巖型鈾礦的浸出性能，解釋新疆某鈾礦C16-1 采區(qū)浸出性能低于C16-2 采區(qū)的問題。結果對于促進砂巖型礦床鈾礦資源的綜合利用、提高礦山的經濟效益具有十分重要的現(xiàn)實意義。

1 模擬方法與評價模型

1.1 溶浸液流動和流體顆粒運動方程

在原地浸出采鈾工藝過程中，假設巖石內部流體不可壓縮且處于層流狀態(tài)，滲流方程可采用適合低速低滲多孔介質的達西定律來描述，在COMSOL 中的達西定律表示為：

式中：u—流體達西速度，m·s-1；k—滲透率，m2；μ—流體動力黏度，Pa·s；p—礦巖孔隙中流體壓力，Pa；ρ—流體密度，kg·L-1；g—重力加速度，m·s-2；Hp—壓力水頭，m；D—沿重力方向的高程，m；?—拉普拉斯算子。

流體流動過程中保持質量守恒，質量守恒定律可以用瞬態(tài)斯托克斯方程描述：

式中：ε—孔隙率，%；t—時間，s；Qm—質量源項，即單位體積中流體流動的速率，kg·m-3·s-1。

在COMSOL 中，抽注流量可以通過井節(jié)點設置，理論方程如下所述：

式中：dw—抽出（注入）井的直徑，m；rw—抽出（注入）井的半徑，m；lw—單位長度，m；S—抽出（注入）井與流體的接觸面積，m2；M0—質量流率，kg·s-1，當井的類型是注入井時，符號為正，否則為負。

采用流體流動顆粒跟蹤模塊對浸出液流體顆粒的運動特征進行模擬，假設流體顆粒為無質量粒子，則耦合模型為：

式中：q—流體顆粒的位置矢量，m。從上式可以看出，顆粒的運動速度來源于達西定律計算出來的達西速度。

1.2 化學反應方程式

新疆某鈾礦C16-1 和C16-2 采區(qū)使用CO2+O2原地浸出鈾礦開采技術，主要原理是利用CO2和O2配制溶浸液，O2的氧化性使地層中的四價鈾氧化成六價鈾，CO2與水及地層碳酸鹽反應生成HCO3-，再與礦層鈾氧化物絡合生成碳酸鈾酰，使鈾溶解浸出并提升至地表進行回收，涉及的化學反應見表1。表中總結了在PHREEQC 計算中使用的與鈾酰離子相關的水解和碳酸鹽絡合反應及其平衡常數(shù)，列出的反應和數(shù)值來源于llnl.dat數(shù)據庫和ThermoChimie 熱力學數(shù)據庫（https：//www.thermochimie-tdb.com/）［16］。

表1 化學反應式［16］Table 1 Chemical reaction formula[16]

1.3 COMSOL-PHREEQC 耦合模擬

本研究基于COMSOL Multiphysics 的達西定律、流體流動顆粒跟蹤模塊，構建可描述砂巖型鈾礦原地浸出過程中浸出液的滲流場和流體顆粒運移場，并耦合PHREEQC 對地浸過程中的地球化學反應進行了化學動力學計算，基于MATLAB 平臺開發(fā)了上述兩款模擬軟件的接口程序，通過空間插值得到的砂巖型鈾礦山地質模型和品位模型將浸出反應和顆粒運移耦合起來，其基本原理見圖1。

圖1 基于MATLAB 耦合COMSOL-PHREEQC 原理圖Fig.1 Schematic diagram of COMSOL-PHREEQC coupling based on MATLAB

2 數(shù)值試驗

為了加快模擬速度和減少運算量，本研究從新疆某鈾礦C16-1 和C16-2 采區(qū)地質模型中截取90 m×67.5 m 大小的剖面模型為研究對象，使用有限元軟件COMSOL 和地球化學模擬軟件PHREEQC 進行耦合計算，將含礦介質設為連續(xù)、非均質的多孔介質，分別對兩個模型進行滲流計算和化學動力學計算，并對其浸出性能進行定量評價。

研究區(qū)域的礦床地質和水文地質條件如下：礦體產狀平緩，傾角為4°～7°；礦體形態(tài)為板狀或似層狀，平面上分布連續(xù)；承壓水頭高度達416.96～474.73 m；礦床滲透系數(shù)為0.08～0.23 m/d；導水系數(shù)為1.4～4.26 m2/d；含礦含水層厚度小，平均厚度為27.64 m；礦體埋深大，為441.4～603.05 m；地下水埋深小，達10.76～41.27 m。基于上述礦床地質和水文地質條件，選取模擬參數(shù)，并建立COMSOL 模型和PHREEQC 模型，實現(xiàn)COMSOL 與PHREEQC的耦合模擬。

2.1 COMSOL 模型建立

C16-1、C16-2 采區(qū)的幾何模型示意圖如圖2 所示，每個模型包含7 種巖層材料屬性，以顏色進行區(qū)分，巖層對應顏色及材料參數(shù)見表2。

圖2 C16-1（a）和C16-2（b）采區(qū)幾何模型Fig.2 Geometric model of C16-1（a）block and C16-2（b）block

表2 采區(qū)模型各巖層材料參數(shù)Table 2 Material parameters of each rock layer in the mining area model

滲流模擬的邊界條件設置如圖3 所示。上下邊界選用無流動邊界，左邊界設置0 m 的壓力水頭，右邊界設置45 m 的壓力水頭。通過井節(jié)點在注入井和生產井的過濾器設置抽注流量，井直徑均設為0.3 m，間距設為30 m。

圖3 滲流模擬條件設置Fig.3 Setting of seepage simulation conditions

采用自由三角形網格，在井的過濾器周圍將網格加密，計算模型網格如圖4 所示。其中C16-1 采區(qū)模型共剖分44 920 個單元，22 523個節(jié)點，共計257 725 個自由度；C16-2 采區(qū)模型共剖分69 753 個單元，34 944 個節(jié)點，共計286 280個自由度。

圖4 計算模型網格圖Fig.4 Grid diagram of computational model

由于頂板和底板的孔隙率和滲透率較低，隔水性良好，模擬期間只考慮在含礦含水層中釋放流體顆粒。在COMSOL 中，流體顆粒可以指定從數(shù)據文件釋放，通過數(shù)據文件的形式將流體顆粒的位置矢量信息傳入軟件中進行模擬計算。品位模型中含有每個品位單元（2.5m×1.25m）的中心坐標和品位信息，由此可以實現(xiàn)品位模型和流體顆粒示蹤的耦合。

2.2 PHREEQC 模型建立

研究區(qū)域中含礦含水層地下水水質類型為HCO3·SO4·Cl-Ca·Na·Mg 型，設模擬溫度為20 ℃，pH 值為7.3。查閱地質報告和相關文獻［9，17-18］，動力學和平衡反應的初始條件離子濃度見表3。

表3 含礦含水層中地下水化學組分/（mg·L-1）Table 3 Chemical composition of groundwater in orebearing aquifer/ (mg·L-1)

在原地浸出過程中，設流體密度和粘滯性保持不變，密度為1 000 kg/m3，動力黏度為0.001 Pa·s，模擬期間對含礦含水層持續(xù)注入適量的CO2和O2，使溶液中CO2和O2濃度維持550 mg/L。在PHREEQC 中，設反應容積為一個品位單元區(qū)域體積與有效孔隙率的乘積（2.5m×1.25m×1m×ε），平衡模擬設置CO2、O2、碳酸鈣的溶解沉淀，動力學設置鈾礦的初始摩爾量、反應表面積，動力學反應速率表達式［19］如下所述：

式中：r—礦物的反應速率，mol·L-1·s-1；SA—每單位水中礦巖表面積，m-2·L-1；A—阿倫尼烏斯指前因子，mol·m-2·s-1；Ea—反應活化能，J·mol-1；R—氣體常數(shù)，8.31446 J·mol-1·K-1；T—溫度，K；ai—物種i的活性度；Ω—礦物飽和指數(shù)，Q·K-1；p，q—經驗指數(shù)。

由于PHREEQC 支持模塊化開發(fā)，本研究基于MATLAB 平臺進行地球化學模擬：首先讀取品位模型的中心坐標信息，利用COMSOL 后處理函數(shù)計算對應位置的孔隙率，進而求得PHREEQC 每個計算單元的反應容積；然后讀取品位模型中的品位信息，求得PHREEQC 中每個計算單元內的鈾礦物摩爾量。

2.3 浸出性能的評價指標

由于砂巖型鈾礦床賦存條件較為復雜，半定量或是定量評價鈾元素浸出性能是一項較為困難的任務?？煽紤]以下幾個重要技術經濟指標以對浸出性能進行合理評價：鈾的浸出率、浸出液中鈾的平均濃度、從礦石中提取鈾的生產率等。

2.3.1 鈾的浸出率

本研究定義鈾的浸出率為鈾元素在浸出液中的含量與鈾元素在含礦層中的含量之比，如式（9）：

式中：RL—鈾的浸出率，%；Uw—鈾元素在浸出液中的含量，g；Uo—鈾元素在礦石中的含量，g。

2.3.2 浸出液的平均鈾濃度

本研究定義模擬區(qū)域內浸出液的平均鈾濃度如式（10）：

式中：Ca—浸出液的平均鈾濃度，μmol·kg-1；Ci—PHREEQC 中每個計算單元中鈾元素的濃度，μmol·kg-1；n—PHREEQC 的計算單元數(shù)量。

2.3.3 鈾的生產率

本研究定義從礦石中提取鈾的生產率為浸出液中流體顆粒經由抽出井的過濾器被泵送至地表的鈾元素總量與含礦層中鈾元素總量之比，如式（11）：

式中：RP—從礦石中提取鈾的生產率，%；Uj—在含礦含水層中運移至抽出井過濾器的浸出液中的鈾含量，g；m—被泵送至地表的浸出液單元數(shù)量。

3 結果與討論

圖5 為在抽注循環(huán)作用及初始水力條件下計算區(qū)域速度場分布圖。由圖可知，在抽出井和注入井的過濾器之間的流體速度較于其他區(qū)域的流體速度更高，溶浸液在含礦含水層中的流動主要集中在粗砂巖和中砂巖區(qū)域，即存在溶浸液在具有復雜地質屬性的多孔介質中沿一些滲透性能較好、流體流動阻力小的區(qū)域集中運移的優(yōu)勢流現(xiàn)象。

圖5 C16-1（a）和C16-2（b）模型速度場分布云圖Fig.5 Cloud diagram of velocity field distribution of block model C16-1（a）and C16-2（b）

圖6 為C16-1、C16-2 采區(qū)模型中溶浸液流動顆粒的初始分布圖（a、c）及某時刻下兩采區(qū)中流動顆粒的分布圖（b、d）。由圖可知，隨著時間的推移，含礦含水層中的流體顆粒在天然壓力水頭和抽注流量補給排泄的擾動下逐漸向抽出井過濾器或左右邊界運移。

礦石密度為1.90 g/cm3，基于品位模型計算C16-1 和C16-2 采區(qū)模型的鈾礦儲量，分別為184.34 kg、476.59 kg，即C16-1 采區(qū)的鈾礦儲量比C16-2 采區(qū)高2～3 倍。設模擬時間為60 天，則在模擬條件下基于COMSOL-PHREEQC 耦合的鈾浸出率計算結果如圖7 所示。

由圖7 可知，隨著反應運移時間的延長，鈾浸出率呈近似線性增長，C16-2 的增長速率遠大于C16-1。前人研究發(fā)現(xiàn)，液固比對鈾的浸出率影響顯著［20］，而在PHREEQC 和COMSOL耦合的計算單元內，C16-1 采區(qū)的液固比要遠小于C16-2 采區(qū)，具體的表現(xiàn)是液固比明顯與地層的滲透性能有關，而C16-1 模型的滲透性能較C16-2 更弱。HE K 等［21］提出了基于凸包搜索策略的空間信息熵方法來評價砂巖鈾礦（品位、巖性）非均質性，結果表明：C16-2 采區(qū)各地質屬性空間集中成片、連續(xù)性較好，而C16-1 采區(qū)則相反。

圖7 C16-1 和C16-2 中鈾浸出率隨時間的變化情況Fig.7 Changes of uranium leaching rate over time in Block C16-1 and Block C16-2

計算C16-1、C16-2 采區(qū)中各巖層的占比，見圖8。從圖中可以發(fā)現(xiàn)：C16-1 中含礦含水層含有大片區(qū)域的泥層和砂泥混合，極大地影響了含礦含水層中溶浸液的流動和鈾礦石的充分溶解；相對地，C16-2 中含礦含水層泥層和砂泥混合分布較少，且表現(xiàn)出良好滲透性能的粗砂和中砂層分布較多。綜上，可得出以下結論：在控制反應浸出時間、溶浸劑濃度相同的條件下，C16-2 的浸出反應比C16-1 更充分。

圖8 C16-1 和C16-2 模型巖層分布圖Fig.8 Rock layer distribution dirgram of block model C16-1 and C16-2

圖9 也說明了這一點，可以看到：在反應浸出時間介于0～20 天時，隨著時間的延長，浸出液的平均鈾濃度均快速增長，此時C16-1 和C16-2 的增長速率相當。然而在20～60 天的范圍內，C16-1、C16-2 采區(qū)鈾礦石中鈾元素逐漸減少，但由于兩采區(qū)地質屬性的差異，增長速率有不同程度的減緩。據統(tǒng)計，現(xiàn)階段蒙其古爾鈾礦床C16-1、C16-2 采區(qū)浸出液的平均鈾濃度分別為2.09×103μmol/L 和5.83×103μmol/L，即：隨著原地浸出鈾礦山的持續(xù)運行，C16-2 采區(qū)浸出液的平均鈾濃度優(yōu)于C16-1 采區(qū)，這一點與模擬結果相符合。C16-1 中礦石鈾品位明顯高于C16-2 但滲透性能低于C16-2，溶浸液與鈾礦石表面積之間具有較低的反應接觸面積，使得其浸出反應速率較慢，鈾礦石不能充分溶解，最終導致平均鈾濃度逐漸低于C16-2。

圖9 浸出液中鈾的平均濃度的比較Fig.9 Comparison of the average uranium concentration in the leachate

基于上述分析，可以預見C16-1 模型中鈾生產率低于C16-2 模型，見圖10。據統(tǒng)計，當前新疆某鈾礦C16-1、C16-2 采區(qū)鈾的生產率分別為5.81 %和42.67 %，即：實際生產中C16-2 采區(qū)鈾的生產率遠遠高于C16-1 采區(qū)，模擬結果與實際情況相符。工程上通常以礦層位置和礦石品位作為布置抽/注鉆孔和過濾器的重要依據，正如C16-1 采區(qū)鉆孔和過濾器布置沒有很好地考慮到地層滲透性能對生產率的影響。因此，礦山在采礦設計中還應當考慮礦層深度、厚度和地層滲透性能等地質屬性。

圖10 C16-1 和C16-2 中鈾生產率隨時間的變化情況Fig.10 Changes of uranium productivity over time in Block C16-1 and C16-2

4 結論

采用基于有限元理論和地球化學模擬理論的COMSOL-PHREEQC 耦合模型對砂巖型鈾礦原地浸出過程進行數(shù)值模擬，對新疆某鈾礦的兩個采區(qū)進行了浸出性能評價，探討了地層滲透性能對浸出性能的影響，得到以下結論：

1）溶浸液在含礦含水層中的流動主要集中在粗砂巖和中砂巖區(qū)域，存在溶浸液沿一些滲透性能較好、流體流動阻力小的區(qū)域集中運移的優(yōu)勢流現(xiàn)象。

2）在控制反應浸出時間、溶浸劑濃度相同的條件下，礦石品位較低但滲透性能較好的C16-2采區(qū)的浸出反應比高品位、低滲透的C16-1 更充分，表現(xiàn)為：隨著時間的延長，C16-2 的鈾浸出率遠大于C16-1，C16-2 的浸出液的平均鈾濃度高于C16-1。

3）C16-1模型中鈾生產率遠低于C16-2模型，建議礦山在采礦設計中綜合考慮礦層位置、深度、厚度和礦石品位以及地層滲透性能等地質屬性因素對砂巖型鈾礦浸出性能的影響，以提高可地浸礦山的浸出性能，進而提高礦山的經濟效益。