賈明濤，金家聰，陳梅芳，蘇學(xué)斌

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

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

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

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

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

1 模擬方法與評(píng)價(jià)模型

1.1 溶浸液流動(dòng)和流體顆粒運(yùn)動(dòng)方程

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

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

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

式中：ε—孔隙率，%；t—時(shí)間，s；Qm—質(zhì)量源項(xiàng)，即單位體積中流體流動(dòng)的速率，kg·m-3·s-1。

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

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

采用流體流動(dòng)顆粒跟蹤模塊對(duì)浸出液流體顆粒的運(yùn)動(dòng)特征進(jìn)行模擬，假設(shè)流體顆粒為無質(zhì)量粒子，則耦合模型為：

式中：q—流體顆粒的位置矢量，m。從上式可以看出，顆粒的運(yùn)動(dòng)速度來源于達(dá)西定律計(jì)算出來的達(dá)西速度。

1.2 化學(xué)反應(yīng)方程式

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

表1 化學(xué)反應(yīng)式［16］Table 1 Chemical reaction formula[16]

1.3 COMSOL-PHREEQC 耦合模擬

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

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

2 數(shù)值試驗(yàn)

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

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

2.1 COMSOL 模型建立

C16-1、C16-2 采區(qū)的幾何模型示意圖如圖2 所示，每個(gè)模型包含7 種巖層材料屬性，以顏色進(jìn)行區(qū)分，巖層對(duì)應(yīng)顏色及材料參數(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

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

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

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

圖4 計(jì)算模型網(wǎng)格圖Fig.4 Grid diagram of computational model

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

2.2 PHREEQC 模型建立

研究區(qū)域中含礦含水層地下水水質(zhì)類型為HCO3·SO4·Cl-Ca·Na·Mg 型，設(shè)模擬溫度為20 ℃，pH 值為7.3。查閱地質(zhì)報(bào)告和相關(guān)文獻(xiàn)［9，17-18］，動(dòng)力學(xué)和平衡反應(yīng)的初始條件離子濃度見表3。

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

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

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

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

2.3 浸出性能的評(píng)價(jià)指標(biāo)

由于砂巖型鈾礦床賦存條件較為復(fù)雜，半定量或是定量評(píng)價(jià)鈾元素浸出性能是一項(xiàng)較為困難的任務(wù)。可考慮以下幾個(gè)重要技術(shù)經(jīng)濟(jì)指標(biāo)以對(duì)浸出性能進(jìn)行合理評(píng)價(jià)：鈾的浸出率、浸出液中鈾的平均濃度、從礦石中提取鈾的生產(chǎn)率等。

2.3.1 鈾的浸出率

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

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

2.3.2 浸出液的平均鈾濃度

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

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

2.3.3 鈾的生產(chǎn)率

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

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

3 結(jié)果與討論

圖5 為在抽注循環(huán)作用及初始水力條件下計(jì)算區(qū)域速度場分布圖。由圖可知，在抽出井和注入井的過濾器之間的流體速度較于其他區(qū)域的流體速度更高，溶浸液在含礦含水層中的流動(dòng)主要集中在粗砂巖和中砂巖區(qū)域，即存在溶浸液在具有復(fù)雜地質(zhì)屬性的多孔介質(zhì)中沿一些滲透性能較好、流體流動(dòng)阻力小的區(qū)域集中運(yùn)移的優(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ū)模型中溶浸液流動(dòng)顆粒的初始分布圖（a、c）及某時(shí)刻下兩采區(qū)中流動(dòng)顆粒的分布圖（b、d）。由圖可知，隨著時(shí)間的推移，含礦含水層中的流體顆粒在天然壓力水頭和抽注流量補(bǔ)給排泄的擾動(dòng)下逐漸向抽出井過濾器或左右邊界運(yùn)移。

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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