基于COMSOL-PHREEQC的砂岩型铀矿浸出性能模拟分析技术

贾明涛，金家聪，陈梅芳，苏学斌

（1.中南大学 资源与安全工程学院，湖南 长沙 410083；2.新疆中核天山铀业有限公司，新疆 伊宁 835000；3.中国铀业有限公司，北京 100010）

我国是铀矿资源储量丰富的国家之一，其中砂岩型铀矿约占三分之一［1-2］。原地浸出技术目前是砂岩型铀矿的首选开采方法，具有安全性高、开采效率高、污染小等优势［3］。然而，原地浸出过程包含复杂的物理-化学作用，如有溶浸剂的流动、铀矿物溶解与沉淀反应、浸出液元素的运移等，因此浸出率通常难以准确评价。

为量化评价铀矿浸出性能，国内外学者们做了大量研究。Liddell K.C.等［4］利用H2O2-NH4HCO3-（NH4）2CO3溶液模拟沥青铀矿和方解石碱法地浸，强调矿物反应和溶液迁移耦合的重要性。郑和秋野［5］通过现场试验模拟铀矿床的酸法地浸，认为试验单元的浸出效果与矿石品位、含矿层地球化学环境相关。李春光［6］通过对砂岩铀矿进行室内试验模拟，认为加入表面活性剂有助于提高矿石的浸出性能。随着计算机技术的发展，数值模拟可以克服现场试验和室内试验存在的缺点，逐渐成为原地浸出动力学研究和工程计算重要且有效的手段。曾晟等［7］基于编程语言开发模拟程序，建立多因素多过程耦合作用下的动力学模型，对铀矿山溶质迁移机理开展研究。黄群英等［8］基于Visual Modflow 模拟和水文地球化学模式探讨了某砂岩型铀矿地浸单元中溶质运移和溶液渗流的关系。林效宾等［9］应用PHREEQC 对铀矿床地下水系统中的水文地球化学过程进行模拟。Lagneau V 等［10］基于HYTEC代码对原地浸出铀矿山进行反应运移模拟，但没有考虑复杂地质属性的空间变异性。虽然前人在原地浸出采铀过程的数学模型方面开展了一定的研究，但是主要集中在地下水动力场模拟或地球化学模拟计算，较少考虑到铀矿浸出反应中所包含的溶解-沉淀、水相络合、氧化还原反应等一系列复杂地球化学反应与溶浸液运移的耦合效应。

COMSOL Multiphysics（简称COMSOL）是一款基于有限元理论的数值分析软件包，可模拟各种物理现象，并可实现多物理场耦合模型的仿真计算［11］。PHREEQC 是一款当前国际上通用的水文地球化学模拟软件，可实现矿物溶解沉淀、水溶物配合、氧化还原、离子交换等功能的模拟［12］，且在软件后续的迭代升级中加入了模块化开发的支持［13］，使得数值软件能够调用其强大的地球化学反应模拟功能，由此拉开了COMSOL 与PHREEQC 耦合模拟研究的序幕。相关研究包括COMSOL-IPHREEQC 模式［14］、Interface COMSOLPHREEQC（iCP）［15］等。

本研究基于MATLAB 平台开发了COMSOL Multiphysics 和PHREEQC 的耦合接口程序，构建了砂岩型铀矿原地浸出过程中浸出液的渗流场和流体颗粒运移场，实现了地浸过程中的化学动力学计算，提出了适用于砂岩型铀矿浸出性能的评价指标，对C16-1、C16-2 采区分别进行了模拟计算，数字化定量分析了砂岩型铀矿的浸出性能，解释新疆某铀矿C16-1 采区浸出性能低于C16-2 采区的问题。结果对于促进砂岩型矿床铀矿资源的综合利用、提高矿山的经济效益具有十分重要的现实意义。

1 模拟方法与评价模型

1.1 溶浸液流动和流体颗粒运动方程

在原地浸出采铀工艺过程中，假设岩石内部流体不可压缩且处于层流状态，渗流方程可采用适合低速低渗多孔介质的达西定律来描述，在COMSOL 中的达西定律表示为：

式中：u—流体达西速度，m·s-1；k—渗透率，m2；μ—流体动力黏度，Pa·s；p—矿岩孔隙中流体压力，Pa；ρ—流体密度，kg·L-1；g—重力加速度，m·s-2；Hp—压力水头，m；D—沿重力方向的高程，m；∇—拉普拉斯算子。

流体流动过程中保持质量守恒，质量守恒定律可以用瞬态斯托克斯方程描述：

式中：ε—孔隙率，%；t—时间，s；Qm—质量源项，即单位体积中流体流动的速率，kg·m-3·s-1。

在COMSOL 中，抽注流量可以通过井节点设置，理论方程如下所述：

式中：dw—抽出（注入）井的直径，m；rw—抽出（注入）井的半径，m；lw—单位长度，m；S—抽出（注入）井与流体的接触面积，m2；M0—质量流率，kg·s-1，当井的类型是注入井时，符号为正，否则为负。

采用流体流动颗粒跟踪模块对浸出液流体颗粒的运动特征进行模拟，假设流体颗粒为无质量粒子，则耦合模型为：

式中：q—流体颗粒的位置矢量，m。从上式可以看出，颗粒的运动速度来源于达西定律计算出来的达西速度。

1.2 化学反应方程式

新疆某铀矿C16-1 和C16-2 采区使用CO2+O2原地浸出铀矿开采技术，主要原理是利用CO2和O2配制溶浸液，O2的氧化性使地层中的四价铀氧化成六价铀，CO2与水及地层碳酸盐反应生成HCO3-，再与矿层铀氧化物络合生成碳酸铀酰，使铀溶解浸出并提升至地表进行回收，涉及的化学反应见表1。表中总结了在PHREEQC 计算中使用的与铀酰离子相关的水解和碳酸盐络合反应及其平衡常数，列出的反应和数值来源于llnl.dat数据库和ThermoChimie 热力学数据库（https：//www.thermochimie-tdb.com/）［16］。

表1 化学反应式［16］Table 1 Chemical reaction formula[16]

1.3 COMSOL-PHREEQC 耦合模拟

本研究基于COMSOL Multiphysics 的达西定律、流体流动颗粒跟踪模块，构建可描述砂岩型铀矿原地浸出过程中浸出液的渗流场和流体颗粒运移场，并耦合PHREEQC 对地浸过程中的地球化学反应进行了化学动力学计算，基于MATLAB 平台开发了上述两款模拟软件的接口程序，通过空间插值得到的砂岩型铀矿山地质模型和品位模型将浸出反应和颗粒运移耦合起来，其基本原理见图1。

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

2 数值试验

为了加快模拟速度和减少运算量，本研究从新疆某铀矿C16-1 和C16-2 采区地质模型中截取90 m×67.5 m 大小的剖面模型为研究对象，使用有限元软件COMSOL 和地球化学模拟软件PHREEQC 进行耦合计算，将含矿介质设为连续、非均质的多孔介质，分别对两个模型进行渗流计算和化学动力学计算，并对其浸出性能进行定量评价。

研究区域的矿床地质和水文地质条件如下：矿体产状平缓，倾角为4°～7°；矿体形态为板状或似层状，平面上分布连续；承压水头高度达416.96～474.73 m；矿床渗透系数为0.08～0.23 m/d；导水系数为1.4～4.26 m2/d；含矿含水层厚度小，平均厚度为27.64 m；矿体埋深大，为441.4～603.05 m；地下水埋深小，达10.76～41.27 m。基于上述矿床地质和水文地质条件，选取模拟参数，并建立COMSOL 模型和PHREEQC 模型，实现COMSOL 与PHREEQC的耦合模拟。

2.1 COMSOL 模型建立

C16-1、C16-2 采区的几何模型示意图如图2 所示，每个模型包含7 种岩层材料属性，以颜色进行区分，岩层对应颜色及材料参数见表2。

图2 C16-1（a）和C16-2（b）采区几何模型Fig.2 Geometric model of C16-1（a）block and C16-2（b）block

表2 采区模型各岩层材料参数Table 2 Material parameters of each rock layer in the mining area model

渗流模拟的边界条件设置如图3 所示。上下边界选用无流动边界，左边界设置0 m 的压力水头，右边界设置45 m 的压力水头。通过井节点在注入井和生产井的过滤器设置抽注流量，井直径均设为0.3 m，间距设为30 m。

图3 渗流模拟条件设置Fig.3 Setting of seepage simulation conditions

采用自由三角形网格，在井的过滤器周围将网格加密，计算模型网格如图4 所示。其中C16-1 采区模型共剖分44 920 个单元，22 523个节点，共计257 725 个自由度；C16-2 采区模型共剖分69 753 个单元，34 944 个节点，共计286 280个自由度。

图4 计算模型网格图Fig.4 Grid diagram of computational model

由于顶板和底板的孔隙率和渗透率较低，隔水性良好，模拟期间只考虑在含矿含水层中释放流体颗粒。在COMSOL 中，流体颗粒可以指定从数据文件释放，通过数据文件的形式将流体颗粒的位置矢量信息传入软件中进行模拟计算。品位模型中含有每个品位单元（2.5m×1.25m）的中心坐标和品位信息，由此可以实现品位模型和流体颗粒示踪的耦合。

2.2 PHREEQC 模型建立

研究区域中含矿含水层地下水水质类型为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，模拟期间对含矿含水层持续注入适量的CO2和O2，使溶液中CO2和O2浓度维持550 mg/L。在PHREEQC 中，设反应容积为一个品位单元区域体积与有效孔隙率的乘积（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—气体常数，8.31446 J·mol-1·K-1；T—温度，K；ai—物种i的活性度；Ω—矿物饱和指数，Q·K-1；p，q—经验指数。

由于PHREEQC 支持模块化开发，本研究基于MATLAB 平台进行地球化学模拟：首先读取品位模型的中心坐标信息，利用COMSOL 后处理函数计算对应位置的孔隙率，进而求得PHREEQC 每个计算单元的反应容积；然后读取品位模型中的品位信息，求得PHREEQC 中每个计算单元内的铀矿物摩尔量。

2.3 浸出性能的评价指标

由于砂岩型铀矿床赋存条件较为复杂，半定量或是定量评价铀元素浸出性能是一项较为困难的任务。可考虑以下几个重要技术经济指标以对浸出性能进行合理评价：铀的浸出率、浸出液中铀的平均浓度、从矿石中提取铀的生产率等。

2.3.1 铀的浸出率

本研究定义铀的浸出率为铀元素在浸出液中的含量与铀元素在含矿层中的含量之比，如式（9）：

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

2.3.2 浸出液的平均铀浓度

本研究定义模拟区域内浸出液的平均铀浓度如式（10）：

式中：Ca—浸出液的平均铀浓度，μmol·kg-1；Ci—PHREEQC 中每个计算单元中铀元素的浓度，μmol·kg-1；n—PHREEQC 的计算单元数量。

2.3.3 铀的生产率

本研究定义从矿石中提取铀的生产率为浸出液中流体颗粒经由抽出井的过滤器被泵送至地表的铀元素总量与含矿层中铀元素总量之比，如式（11）：

式中：RP—从矿石中提取铀的生产率，%；Uj—在含矿含水层中运移至抽出井过滤器的浸出液中的铀含量，g；m—被泵送至地表的浸出液单元数量。

3 结果与讨论

图5 为在抽注循环作用及初始水力条件下计算区域速度场分布图。由图可知，在抽出井和注入井的过滤器之间的流体速度较于其他区域的流体速度更高，溶浸液在含矿含水层中的流动主要集中在粗砂岩和中砂岩区域，即存在溶浸液在具有复杂地质属性的多孔介质中沿一些渗透性能较好、流体流动阻力小的区域集中运移的优势流现象。

图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 采区模型中溶浸液流动颗粒的初始分布图（a、c）及某时刻下两采区中流动颗粒的分布图（b、d）。由图可知，随着时间的推移，含矿含水层中的流体颗粒在天然压力水头和抽注流量补给排泄的扰动下逐渐向抽出井过滤器或左右边界运移。

矿石密度为1.90 g/cm3，基于品位模型计算C16-1 和C16-2 采区模型的铀矿储量，分别为184.34 kg、476.59 kg，即C16-1 采区的铀矿储量比C16-2 采区高2～3 倍。设模拟时间为60 天，则在模拟条件下基于COMSOL-PHREEQC 耦合的铀浸出率计算结果如图7 所示。

由图7 可知，随着反应运移时间的延长，铀浸出率呈近似线性增长，C16-2 的增长速率远大于C16-1。前人研究发现，液固比对铀的浸出率影响显著［20］，而在PHREEQC 和COMSOL耦合的计算单元内，C16-1 采区的液固比要远小于C16-2 采区，具体的表现是液固比明显与地层的渗透性能有关，而C16-1 模型的渗透性能较C16-2 更弱。HE K 等［21］提出了基于凸包搜索策略的空间信息熵方法来评价砂岩铀矿（品位、岩性）非均质性，结果表明：C16-2 采区各地质属性空间集中成片、连续性较好，而C16-1 采区则相反。

图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 采区中各岩层的占比，见图8。从图中可以发现：C16-1 中含矿含水层含有大片区域的泥层和砂泥混合，极大地影响了含矿含水层中溶浸液的流动和铀矿石的充分溶解；相对地，C16-2 中含矿含水层泥层和砂泥混合分布较少，且表现出良好渗透性能的粗砂和中砂层分布较多。综上，可得出以下结论：在控制反应浸出时间、溶浸剂浓度相同的条件下，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 采区铀矿石中铀元素逐渐减少，但由于两采区地质属性的差异，增长速率有不同程度的减缓。据统计，现阶段蒙其古尔铀矿床C16-1、C16-2 采区浸出液的平均铀浓度分别为2.09×103μmol/L 和5.83×103μmol/L，即：随着原地浸出铀矿山的持续运行，C16-2 采区浸出液的平均铀浓度优于C16-1 采区，这一点与模拟结果相符合。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。据统计，当前新疆某铀矿C16-1、C16-2 采区铀的生产率分别为5.81 %和42.67 %，即：实际生产中C16-2 采区铀的生产率远远高于C16-1 采区，模拟结果与实际情况相符。工程上通常以矿层位置和矿石品位作为布置抽/注钻孔和过滤器的重要依据，正如C16-1 采区钻孔和过滤器布置没有很好地考虑到地层渗透性能对生产率的影响。因此，矿山在采矿设计中还应当考虑矿层深度、厚度和地层渗透性能等地质属性。

图10 C16-1 和C16-2 中铀生产率随时间的变化情况Fig.10 Changes of uranium productivity over time in Block C16-1 and C16-2

4 结论

采用基于有限元理论和地球化学模拟理论的COMSOL-PHREEQC 耦合模型对砂岩型铀矿原地浸出过程进行数值模拟，对新疆某铀矿的两个采区进行了浸出性能评价，探讨了地层渗透性能对浸出性能的影响，得到以下结论：

1）溶浸液在含矿含水层中的流动主要集中在粗砂岩和中砂岩区域，存在溶浸液沿一些渗透性能较好、流体流动阻力小的区域集中运移的优势流现象。

2）在控制反应浸出时间、溶浸剂浓度相同的条件下，矿石品位较低但渗透性能较好的C16-2采区的浸出反应比高品位、低渗透的C16-1 更充分，表现为：随着时间的延长，C16-2 的铀浸出率远大于C16-1，C16-2 的浸出液的平均铀浓度高于C16-1。

3）C16-1模型中铀生产率远低于C16-2模型，建议矿山在采矿设计中综合考虑矿层位置、深度、厚度和矿石品位以及地层渗透性能等地质属性因素对砂岩型铀矿浸出性能的影响，以提高可地浸矿山的浸出性能，进而提高矿山的经济效益。