两相流与溶液化学平衡相结合的地热井结垢模拟研究*

岑继文，蒋方明

两相流与溶液化学平衡相结合的地热井结垢模拟研究*

岑继文1,2,3，蒋方明1,2,3†

（1. 中国科学院广州能源研究所，先进能源系统研究室，广州 510640；2. 中国科学院可再生能源重点实验室，广州 510640； 3. 广东省新能源和可再生能源研究开发与应用重点实验室，广州 510640）

水热型地热系统和干热岩增强型地热系统在运行一段时间后都有可能出现井筒和管道设备结垢现象，困扰和阻碍着地热资源的高效低成本开发。由于温度和压力变化导致CO2从地热水中逸出，进而使得CaCO3析出结垢，是地热生产井堵塞的主要原因。将两相流动换热计算与水溶液物理化学模拟相结合，针对地热井内CaCO3垢进行流动换热与化学组分变化的模拟和计算方法的研究，获得了特定地热条件下井筒内流动状况与化学组分变化情况的详细分析结果，可为井筒地热水流动进行结垢分析和预测。

地热水；地热井；结垢；计算机模拟；预测

0 引 言

2020年我国首次提出了要在2030年实现碳达峰，2060年实现碳中和的目标。由于传统石化能源的不可再生性以及污染问题，大力发展新能源是实现“双碳”目标的最佳途径。地热能作为可再生能源的一种，与太阳能、风能等可再生能源相比，具有资源储量巨大、稳定性相对较高等优点。而地热开发过程中，尤其是水热型地热资源开发过程中遇到的管道设备结垢问题，困扰和阻碍地热资源的高效低成本开发。地热电站中结垢类型主要为钙垢、硅垢、金属氧化垢及地热流体中含有的微细颗粒物、污泥、微生物、硫化物等。其中钙垢主要与压力有关，通常发生在井筒内压力降低而导致的闪蒸点附近，闪蒸使得大量二氧化碳析出，卤水酸性降低，碳酸钙过饱和而形成钙垢析出[4]。另外，CaCO3结垢堵塞在地热尾水回灌过程中是一个普遍的现象[2-3]，就单个井而言，结垢现象在井底、井上及管道设备中都有可能出现，使地热能的效率大为降低。

刘明言[1]指出，目前在地热流体的腐蚀与结垢趋势预测方面的工作还相对较少，针对我国不同地区地热流体特性开展系统的腐蚀结垢趋势分析，以及三维多相传递和化学反应数值模拟研究是今后的方向。朱家玲等[5]对地热水结垢腐蚀趋势的判断方法进行了探讨，在地热水的Cl−物质的量分数大于25%时，用雷诺兹指数判断有局限性，此时，应结合拉伸指数进行水质判断，才能使判断结果符合实际。PÁTZAY等[6]建立了CaCO3-H2O-CO2体系模拟地热井结垢方法，采用了Davies和Pitzer活度计算方法，对井内流动压降和温度变化采用了简单的线性化处理方法。ZOLFAGHARROSHAN等[7]用较为详细的多相流方法模拟了井筒流动换热情况，但对结垢的模拟仅仅用地面的成分情况进行判断，并未全程对井筒进行物理化学组分变化的计算跟踪。张恒等[8]利用美国地质调查局研发的PHREEQC软件，采用水文地球化学模拟技术对康定某高温地热井结垢问题进行了分析研究。结果表明，该地热井在开采过程中随着地热流体温度、pH、压力、氧化还原环境等条件变化发生严重的结垢现象。此外，除了水热型地热井，长期运行的干热岩系统也可能会在井筒和地面设备中出现矿物结垢现象[9]。

地热水结垢的基本物理化学机理已被研究得较为清楚，很多文献已进行过详细的论述。而在实际应用过程中，要想对结垢现象进行防治，需要在了解机理的基础上建立较为精确的计算模型对其预测。本文针对地热井内CaCO3垢形成情况进行详细的两相流动换热模拟，并对地热流体中化学组分的变化进行化学平衡计算模拟，以期获得特定地热条件下结垢情况的细致而精确的预测。

1 计算模型

1.1 井筒内CaCO3结垢基本物理化学机理

井筒结垢的形成基本机理是由于地热水溶液上升、压力不断降低。尽管水溶液上升过程中地层温度下降使得水溶液热损失温度有所降低，但由于压力下降更为迅速对应的饱和温度下降更为快速，因此在井筒某处会开始出现闪蒸，原先溶解在地热水溶液中的酸性气体CO2逸出进入蒸汽中，使得溶液pH升高，进而使CaCO3等溶解度达到饱和，析出固体形成结垢。

1.2 计算模型建立

井筒内CaCO3结垢原理虽已获公认，但对其进行精确的建模计算具有相当的复杂性，涉及两相流、地层换热和溶液物理化学变化，尤其是闪蒸过程。不凝气体组分和含量对结垢的影响问题涉及井筒的水溶液化学反应平衡和相平衡方面，以及地热水溶液在管道中的两相流流阻压降和地层热损失换热等方面，本文参考了文献[7,10-12]，建立如下计算模型。

1.2.1 两相流计算模型

质量方程：

其中：为密度，m；为速度，m/s；为时间，s；为井筒长度，m。稳态流动情况下没有质量累积，因此有：

动量方程：

其中：为压力，Pa；为摩擦系数，kg/(m∙s2)，为井筒内径，m；为井筒内截面积，m2；为重力加速度，m/s2；为井筒与水平方向的夹角，无量纲。联合质量方程，整理可得：

如为两相流，可将流体当作气液均相混合物，井为竖直井，上式可变为：

其中：下标mix表示混合参数；为竖直方向长度，m。根据不同的流型计算出不同的摩擦因子，需要用到占空比和混合密度。气液两相流可分为bubbly、slug、churn和annular四种流型[10]。

能量方程：

其中：为焓，J/kg；为单位质量流量单位管长内的换热量，W/(kg/s)，可参考文献[11]：

式中：k为有效导热系数，W/(m∙K)，其地热井典型值为4 W/(m∙K)。

1.2.2 化学平衡、相平衡模型

拟借鉴电解质相平衡模型[12]进行闪蒸过程计算。电解质溶液相平衡模型如图1所示，含有一定量化学组分的来流溶液变化到指定的温度、压力之后，其化学组分在气、液、固三相中重新分配。所涉及的方程如下。

物料衡算式：

其中：in,i和out,i分别为电离或离解反应前后各组分的摩尔流量，mol/s；为化学计量系数，无量纲；为反应程度，mol/s；为某个平衡方程；为总的平衡方程数量。

因此，化学反应平衡后更新的各组分组成为：

总物料衡算式为：

其中：为所有相中物料的总和，为液相物料，为气相物料，单位均为mol/s；为液相组成，为气相组成，为各相中的物料组成，均为无量纲；为总的物料的种类。

能量衡算式：

其中：L、V和F分别为液相、气相和所有相总和的摩尔焓，J/mol；为换热量，W。

相平衡：

其中：为亨利系数，无量纲。上式可理解为根据亨利定律，在一定温度和平衡状态下，气体在液体中的溶解度（用摩尔分数表示）和该气体的平衡分压成正比。

化学反应平衡：

其中：为化学反应平衡常数，单位与具体反应方程式有关；为活度系数，无量纲；为逸度系数，无量纲；为逸度，Pa。液相活度系数的计算采用Pitzer电解质活度系数模型[13]，而气相逸度系数采用 PR状态方程[14]计算。

组分约束：

其中：m为气相分子种类数，m表示分子。

电中性平衡：

此时的Z表示离子种类所带电荷的真实值，而不是绝对值，无单位。

物性关联式：

1.2.3 结垢判断

其中：[Me]是二价阳离子；[An]是二价阴离子；结垢指数SI= 0表示平衡，SI > 0表示过饱和产生结垢；sp为溶解度。

对于常见的碳酸钙结垢，可采用以下关联式[7]：

式中：为温度，℃；为压力，单位为psi；离子强度为：

其中：C为各组分浓度，mol/kgw（kgw表示1 kg地热水）。

1.2.4 气泡点的判断

气泡的起始点可根据各不凝性气体分压加和是否超过了当地压力来进行判断：

1.3 模型求解算法

将井筒内地热水流动简化为一维流动，在长度方向上将井筒分为若干个单元，从底部向上依次逐步推进计算各单元参数，包括换热、流动、组分变化、结垢情况等。单元计算流程如图2所示。

图2 单元计算流程图

闪蒸算法是最为复杂的部分，可参考化工计算模拟的闪蒸分离过程算法[12]，该方法是对于包含个组分（所有离子和分子）并且发生个电离或水解等化学反应的两相体系，根据电离平衡和相平衡关系得到一个+维方程组，并通过Marquardt迭代方法求解非线性方程组实现。

每个单元的计算最后都做一次饱和指数的计算，判断是否开始结垢，如果结垢发生，则析出一定的量，使得溶液正好恢复当时条件下的饱和状态，从而获得各组分平衡关系。

1.4 化学反应方程和相变平衡方程

由于各地地热流体成分的复杂性，不可能做到全部覆盖，本文目的是阐述CaCO3结垢的计算模拟方法，计算程序中目前仅考虑CaCO3、CO2、H2O、NaCl、CH4这几种物质的相互作用，在该物质体系下，涉及的化学反应方程和相变平衡见表1[15]。

表1 计算程序所涉及的化学反应方程和相变平衡

注：下角标aq表示在水溶液中，g表示气相，s表示固相。

2 结果与分析

计算程序采用Python计算机语言，在Visual studio 2019集成开发环境中编写了上述计算模型程序。使用Tkinter库制作简单的参数输入界面，可方便地输入各个算例参数。

目前常用的结垢模拟软件有Wellsim、Phreeqc等，均为国外软件，无法掌握内部计算模型细节，国内尚未有专门的地热井结垢计算软件。本文模拟计算工作的特点是可以考虑各种不凝性气体的影响，两相流计算可细化至不同流型的影响等。而且计算过程中使用到的各个计算模型，如两相流压降计算、活度系数、逸度系数模型、迭代模型等各方面均可随时替换或优化，所有参数均可根据需要做出调整。

2.1 算例分析

算例输入参数：井筒长度为1 000 m、直径为200 mm，管道表面粗糙度0.015，井底温度70℃，压力11.2 MPa，流量2.5 kg/s。输入流体组分Ca2+浓度为4 mg/kgw，HCO3−浓度为12.18 mg/kgw，CO2浓度为50 mg/kgw，CH4浓度为200 mg/kgw。地温梯度按照地面温度25℃至井底线性变化计算。

将井筒划分成500段（可任意设置），经模拟软件计算后获得各单元的气−液组分组成、温度、压力、流速、pH、结垢饱和指数等参数，如表2所示。各气体组分分压为总压乘以其气相摩尔分数及逸度系数；pH可从H+浓度求出；CaCO3(s)累积结垢量为7.0 × 10−5mol/kgw；泡点/开始结垢点位置为118 m。因此，表2中的4个代表深度，分别为井底、井口、泡点前（150 m）和泡点后（60 m）各选取一点作为代表井内各部位参数的变化。

表2 井内地热水不同深度下各组分浓度

图3 （a）温度压力随井深变化曲线；（b）泡点之后CO2分压和pH值随井深变化曲线；（c）泡点之后气相组成随井深变化曲线；（d）泡点后溶液中各组分变化情况

图3为算例的输出数据结果绘制的温度、压力变化，泡点之后的地热水pH，气相、液相组分以及CO2分压变化曲线图。在2.5 kg/s流量的地热条件下，温度下降较小，压力呈线性下降。泡点之后由于CO2逸出，地热水pH逐渐加速升高，CO2分压呈现加速下降的模式，曲线变化趋势符合碳酸钙垢形成的机理。图3c中气体的组分变化表明CH4气体比CO2有更快的增长速度，是由于CH4在水中具有较低的溶解度。图3d为液相中各组分的变化曲线，可知溶液中除了游离CO2分子、CH4分子、HCO3−、Ca2+外，其他组分浓度较低，小于1 × 10−5mol/kgw。总之，泡点后相应数据逐步变化明显，而泡点之前除总压外其他参数几乎不变。

2.2 不凝性气体含量对结垢量和泡点位置的影响分析

图4为以上述算例为基准，改变井底温度和CO2含量两个变量，获得的不同地温条件下CO2含量与最终结垢量的变化曲线。最终结垢量是根据CaCO3饱和指数，将其维持小于零来计算析出量。从图中可以看出，在70℃地温条件下，结垢量随着CO2含量的增加而减少，增加到一定浓度后，结垢量趋近于0。而相同CO2含量下，随着地温的升高，结垢量更多，但由于Ca2+数量的限制，结垢量也有一定的最大值。

图4 不同温度下结垢量与CO2含量关系

图5为上述算例条件下，分别改变每千克地热水中CH4和CO2的含量对泡点位置的影响。从图中可以看出泡点位置与这两种不凝性气体的含量呈现线性增加关系。一般情况下，泡点的位置即为结垢开始出现的位置。随着不凝性气体含量的增加，气体更容易从水溶液中逸出，闪蒸点的位置会更深。

此外，需要说明的是本文算例是根据各类情况综合考虑，通过改变相应参数条件来考察结垢模拟方法和结垢形成特性问题，并非某个特定井的具体参数。如需要对某个地热井实际数据进行验证对比，还需要根据该实际井添加相应的化学组分，进而在程序中关联相应的化学反应方程、相平衡方程，结垢也不仅限于CaCO3垢，可能是硅酸盐、硫酸盐或多种成分都存在，气体也不仅有CH4、CO2，还可能有H2S、NH3、N2等，这些都需要在程序中重新做针对性的修改和录入相关参数。但模型方程和算法与上述一致，是通用的。上述算例中计算出的泡点位置与参考文献[6]大致数量级相同，因此计算结果在合理范围内吻合。

3 结 论

将两相流计算模型与化学反应平衡模型以及化工闪蒸计算模型相结合，对地热井内的流动进行精细模拟，获得温度、压力、溶液及气相中各组分的变化情况，分析其pH变化、气相分压、沸腾起始点，并进行结垢判断等。经过试运算，模拟结果的数值较为合理，能够测算出CO2分压的指数性减小，pH的变化亦呈指数增大趋势。同时，经过对比分析还得出以下两个结论：（1）随着CO2含量的增加，结垢量逐渐减少直至不结垢，而地热水温度越高越容易发生结垢；（2）随着不凝性气体的增加，泡点位置即结垢位置更深。

本文的地热水化学成分组成仅考虑了Ca2+、HCO3−以及CO2和CH4不凝性气体组分，真实的地热井卤水含有更多更复杂的组分组成，可在此模型框架的基础上，增加更多物质组分，并尽可能利用更精准的化学反应平衡常数、亨利系数等化工常数，即可准确地模拟真实的地热井流动和化学变化情况。

[1] 刘明言. 地热流体的腐蚀与结垢控制现状[J]. 新能源进展, 2015, 3(1): 38-46. DOI: 10.3969/j.issn.2095-560X. 2015.01.007.

[2] 云智汉, 马致远, 周鑫, 等. 碳酸盐结垢对中低温地热流体回灌的影响——以咸阳地热田为例[J]. 地下水, 2014, 36(2): 31-33. DOI: 10.3969/j.issn.1004-1184.2014. 02.012.

[3] 田涛, 陈玉林, 姚杰. 地热回灌系统水结垢预测[J]. 地下水, 2011, 33(6): 27-28, 91. DOI: 10.3969/j.issn.1004- 1184.2011.06.012.

[4] 蔡正敏, 李刚, 李源, 等. 肯尼亚地热电站结垢问题的日常维护[J]. 科技视界, 2018(25): 41-43. DOI: 10.19694/ j.cnki.issn2095-2457.2018.25.018.

[5] 朱家玲, 姚涛. 地热水腐蚀结垢趋势的判断和计算[J]. 工业用水与废水, 2004, 35(2): 23-25. DOI: 10.3969/ j.issn.1009-2455.2004.02.007.

[6] PÁTZAYG, STÁHLG, KÁRMÁNF H, et al. Modeling of scale formation and corrosion from geothermal water[J]. Electrochimica acta, 1998, 43(1/2): 137-147. DOI: 10.1016/S0013-4686(97)00242-9.

[7] ZOLFAGHARROSHAN M, KHAMEHCHI E. A rigorous approach to scale formation and deposition modelling in geothermal wellbores[J]. Geothermics, 2020, 87: 101841. DOI: 10.1016/j.geothermics.2020.101841.

[8] 张恒, 胡亚召, 云智汉, 等. 水文地球化学模拟技术在康定某高温地热井结垢研究中的应用[J]. 新能源进展, 2016, 4(2): 111-117. DOI: 10.3969/j.issn.2095-560X.2016. 02.006.

[9] YANAGISAWA N, MATSUNAGA I, SUGITA H, et al. Temperature-dependent scale precipitation in the Hijiori Hot Dry Rock system, Japan[J]. Geothermics, 2008, 37: 1-18. DOI: 10.1016/j.geothermics.2007.08.003.

[10] HASAN A R, KABIR C S. Modeling two-phase fluid and heat flows in geothermal wells[J]. Journal of petroleum science and engineering, 2010, 71(1/2): 77-86. DOI: 10.1016/j.petrol.2010.01.008.

[11] GARG S K, PRITCHETT J W, Alexander J H. A new liquid hold-up correlation for geothermal wells[J]. Geothermics, 2004, 33(6): 795-817. DOI: 10.1016/j.geothermics.2004. 07.002.

[12] 韩莎莎. 含电解质过程模拟研究与开发[D]. 青岛: 青岛科技大学, 2019.

[13] 韩莎莎, 郑俊强, 孙晓岩, 等. Pitzer活度系数模型研究与开发[J]. 当代化工, 2020, 49(1): 221-227. DOI: 10.3969/j.issn.1671-0460.2020.01.053.

[14] LI J, WEI L L, LI X C. Modeling of CO2-CH4-H2S-brine based on cubic EOS and fugacity-activity approach and their comparisons[J]. Energy procedia, 2014, 63: 3598- 3607. DOI: 10.1016/j.egypro.2014.11.390.

[15] DUAN Z H, LI D D. Coupled phase and aqueous species equilibrium of the H2O–CO2–NaCl–CaCO3system from 0 to 250oC, 1 to 1000 bar with NaCl concentrations up to saturation of halite[J]. Geochimica et cosmochimica acta, 2008, 72(20): 5128-5145. DOI:10.1016/j.gca.2008.07.025.

Geothermal Well Scaling Simulation by Combining Two Phase Flow with Chemical Equilibrium

CEN Ji-wen1,2,3, JIANG Fang-ming1,2,3

(1. Laboratory of Advanced Energy Systems, Guangzhou Institute of Energy Conversion, Chinese Academy of Sciences,Guangzhou 510640, China; 2. CAS Key Laboratory of Renewable Energy, Guangzhou 510640, China;3. Guangdong Provincial Key Laboratory of New and Renewable Energy Research and Development, Guangzhou 510640, China)

Both hydrothermal geothermal system and hot dry rock enhancement geothermal system may experience scaling problem in the wells or facility on the ground after an operation period. Scaling problems inhibit the efficient exploitation of geothermal resources. For carbonate scale in production wells, the temperature and pressure changes in the fluid cause the release of CO2from the geothermal solution, and lead to supersaturation of the solubility of carbonate calcium. In this paper, two-phase flow model and chemical reaction simulation were integrated together to simulate the scaling problems in production wells to reveal detail temperature and pressure of the flow states and the changes of chemical components simultaneously. The model can be used to be a scaling prediction and analytical tool for geothermal production wells.

geothermal water; geothermal well; scaling; simulation; prediction

2095-560X（2022）01-0020-07

TK521

A

10.3969/j.issn.2095-560X.2022.01.004

2021-10-20

2021-11-09

国家重点研究发展计划项目（2019YFB1504104）；中国科学院 A 类战略性先导科技专项项目（XDA21060700）

蒋方明，E-mail：jiangfm@ms.giec.ac.cn

岑继文（1979-），男，博士，副研究员，硕士生导师，主要从事地热能开发利用、制冷热泵、电子散热、电动汽车热管理等方面的研究。

蒋方明（1973-），男，博士，研究员，博士生导师，主要从事电化学能量/动力系统、增强型地热系统、微热流体系统、燃料电池水、热管理及高效节能技术/产品等研发工作。