高温地热生产井碳酸钙结垢定量评价：两相流动

Quantitative Assessment of Calcite Scaling of a High Temperature Geothermal Production Well: Two-Phase Flow

摘要

高温地热生产井碳酸钙结垢定量评价涉及到复杂的物理和化学过程，其中井筒中的两相流研究是评价的基础.本文首先基于质量守恒、能量守恒和动量守恒方程，建立了CO₂-H₂O体系井筒两相相变稳定流动模型，提出了稳健的求解方法，并验证了其计算结果的可靠性.然后，在西藏羊八井地热田典型井开展了静止和放喷状态下的井筒中的温度和压力测试，并结合放喷试验，采用开发的模型成功评价了高温地热生产井筒两相流动过程.结果显示：气相和液相之间的速度差对井筒中温度和压力的分布有决定性的影响，不考虑气相和液体之间的速度差，会使模型计算结果远远偏离测量值.在开采速率19.10 kg/s的条件下，计算的井口温度和压力分别约为128 ℃和2.6 bar；井口的气相质量分数在6%~7%之间，对应的井口气相饱和度约为0.84；从闪蒸点往上大概20~30 m气相和液相中CO₂质量分数变化较为剧烈，也是碳酸钙结垢严重井段.

关键词

地热能 / 两相流动 / 二氧化碳 / 数值模拟 / 羊八井 / 水文地质学

Key words

geothermal energy / two-phase flow / CO₂ / numerical simulation / Yangbajing / hydrogeology

0 引言

由于低碳和储量巨大，作为可再生新能源的地热能受到越来越多的关注.理论评价显示地面以下10 km地壳内的地热能总量约为3.6✕10¹⁴ GWh（Lund， 2008），相当于2012年全球能源消耗的217万倍.目前，地热能在全球得到了广泛的开发和应用，大约82个国家在直接利用（Lund et al.， 2011）、26个国家在利用地热能发电（Bertani， 2012）.根据相关评估报告（自然资源部中国地质调查局等， 2018），美国在地热能开发和利用方面处于领先地位，拥有全球最大的地热发电装机容量，占全球的28.8%，其在地热直接利用方面排名第二；中国在地热直接利用方面全球第一，达到了总量的25.2%，但是在地热发电方面仅27.28 MWe，占全球的0.2%，排名第18.

我国地热资源量相对丰富，水热型地热能资源折合标准煤12 500亿吨，浅层地热能资源每年可开采量折合标准煤7亿吨，干热岩远景资源量折合标准煤856万亿吨（王贵玲等，2017）.地热资源分布有明显的规律性和地带性，其中高温地热资源主要分布在藏南、滇西、川西和台湾地区，中低温地热资源主要分布于大型沉积盆地和山地的断裂带上（王贵玲等，2017）.西藏羊八井是我国第一个高温地热发电站，装机容量26.18 MW，并在过去的几十年每年贡献拉萨冬季60%和夏季50%的电力需求（Zhang et al.， 2019）.但是，现在该电站处于停机状态，主要是因为发电效率低、经济效益差（Zhang et al.， 2019）、尾水污染（Guo and Wang， 2009）和井筒结垢（周大吉，2003）等问题.

井筒结垢是川藏地区高温地热利用普遍遇到的问题之一，大大影响了地热的开发效率.井筒结垢涉及到两相相变流动、气体析出、矿物沉淀和固体吸附等复杂的物理化学过程，其中两相相变流动是这些过程的基础.Barelli et al.（1982）进行了含有CO₂的两相井筒流动模拟，并与实测的温度和压力进行了对比.Pan and Oldenburg（2014）在TOUGH2框架中开发了T2Well，并可以用于模拟井筒‒储层耦合系统的多相流动.进一步，Vasini et al.（2018）耦合地热模块EWASG（Battistelli et al.， 1997）到T2Well中，扩展了T2Well对高温地热井筒中两相流动的模拟功能.Gunn and Freeston（1991）开发了地热井筒模拟器WELLSIM，可以进行放喷测试、流体组分和属性、温度和压力的历史拟合等分析.最近，Khasani et al.（2021）开展了地热井筒流动的数值模拟，并评价了CO₂对地热井筒中两相流动的影响.最近研究人员建立了地热井筒质量、能量和压降模型（卜宪标等，2021），并初步评价了河北保定岩溶地热结垢涉及的井筒流动问题（何雨江等，2021）.

地热井筒结垢问题的评价包括两相流动和水文地球化学两大部分，本文主要从两相流动进行分析.首先介绍羊八井地质概况，然后描述羊八井典型地热井的温度和压力测试，再建立两相流动模型和计算方法，最后针对典型井进行历史拟合和分析.本文是后面井筒碳酸钙评价（雷宏武等，2023）的基础.

1 井筒两相流动模型

1.1　数学模型

高温地热生产井井筒中的流体主要成分考虑为H₂O和CO₂，在井筒流动过程中会发生相变.根据质量守恒、能量守恒和动量守恒方程，可以建立CO₂-H₂O体系一维垂直井筒相变流动控制方程.考虑稳定状态，则：

总的质量守恒： \frac{\mathrm{d}}{\mathrm{d}z} \left[{\rho }_{\mathrm{l}}{u}_{\mathrm{l}}\left(\mathrm{1}-\alpha \right)+{\rho }_{\mathrm{g}}{u}_{\mathrm{g}}\alpha \right] = \mathrm{0}，(1)

CO₂质量守恒：\frac{\mathrm{d}}{\mathrm{d}z} \left[{\rho }_{\mathrm{l}}{u}_{\mathrm{l}}\left(\mathrm{1}-\alpha \right)x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}+{\rho }_{\mathrm{g}}{u}_{\mathrm{g}}\alpha y_{\mathrm{g}}^{\mathrm{C}\mathrm{O}\mathrm{2}}\right] = \mathrm{0}，(2)

动量守恒： \frac{\mathrm{d}}{\mathrm{d}z} \left({\dot{m}}_{\mathrm{l}}{u}_{\mathrm{l}}+{\dot{m}}_{\mathrm{g}}{u}_{\mathrm{g}}\right) + A \frac{\mathrm{d}P}{\mathrm{d}z} + {\rho }_{\mathrm{m}} g A + \frac{{\rho }_{\mathrm{m}}fA}{\mathrm{2}d} u_{\mathrm{m}}^{\mathrm{2}} = \mathrm{0}，(3)

能量守恒： \frac{\mathrm{d}}{\mathrm{d}z} [ {\dot{m}}_{\mathrm{l}} \left(\frac{u_{\mathrm{l}}^{\mathrm{2}}}{\mathrm{2}}+gz+h_{\mathrm{l}}\right) + {\dot{m}}_{\mathrm{g}} \left(\frac{u_{\mathrm{g}}^{\mathrm{2}}}{\mathrm{2}}+gz+h_{\mathrm{g}}\right) ] + q = \mathrm{0}，(4)

其中： {\rho }_{\mathrm{l}}和 {\rho }_{\mathrm{g}}分别为液相和气相密度，单位kg/m³； u_{\mathrm{l}}和 u_{\mathrm{g}}分别为液相和气相速度，单位m/s； \alpha为气相体积分数； x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}和 y_{\mathrm{g}}^{\mathrm{C}\mathrm{O}\mathrm{2}}分别为液相和气相中CO₂的质量分数； {\dot{m}}_{\mathrm{l}}和 {\dot{m}}_{\mathrm{g}}为液相和气相质量速度，单位kg/s； A为面积，m²； P = P_{\mathrm{H}\mathrm{2}\mathrm{O}} + P_{\mathrm{C}\mathrm{O}\mathrm{2}}为总压力， P_{\mathrm{H}\mathrm{2}\mathrm{O}}和 P_{\mathrm{C}\mathrm{O}\mathrm{2}}分别为H₂O和CO₂的分压，单位Pa； f为井壁摩擦系数，当雷洛数 ， f = \frac{\mathrm{16}}{Re}，否则 \frac{\mathrm{1}}{\sqrt[]{f}} = - \mathrm{4} \mathrm{l} \mathrm{g} [ \frac{\frac{\mathrm{2}\epsilon }{d}}{\mathrm{3.7}} - \frac{\mathrm{5.02}}{Re} \mathrm{l} \mathrm{g} ( \frac{\frac{\mathrm{2}\epsilon }{d}}{\mathrm{3.7}} + \frac{\mathrm{13}}{Re} ) ]， {\mu }_m为混合物粘度，单位Pa·s， \epsilon为井壁粗糙度，单位m； d为井筒直径，单位m； h_{\mathrm{l}}和 h_{\mathrm{g}}分别为液相和气相热焓，单位J/kg； {\rho }_{\mathrm{m}}为混合密度， {\rho }_{\mathrm{m}} = \left(\mathrm{1}-\alpha \right) {\rho }_{\mathrm{l}} + \alpha {\rho }_{\mathrm{g}}，单位kg/m³； u_{\mathrm{m}}为混合物的平均速度， u_{\mathrm{m}} = \frac{\left(\mathrm{1}-\alpha \right){\rho }_{\mathrm{l}}{u}_{\mathrm{l}}+\alpha {\rho }_{\mathrm{g}}{u}_{\mathrm{g}}}{{\rho }_{\mathrm{m}}}，单位m/s； q为井筒与围岩单位长度的热交换量，单位J/（s·m）， q = U ( T_{\mathrm{f}} - T_{\mathrm{w}\mathrm{b}} )， U为集中参数，受井径、以及围岩、套管和固井水泥等热传导系数和比热容综合影响，一般需要采用实际数据校正，单位为W/（m·℃）， T_{\mathrm{f}}和 T_{\mathrm{w}\mathrm{b}}分别为井筒流体和井壁的温度，单位为℃； z为井筒高程，单位为m； g = \mathrm{9.8}m²/s为重力加速度.

忽略动量随高程的变化，即方程（3）中的第一项，以及忽略速度对能量守恒的贡献，考虑井底一类边界条件，方程（1）~（4）可以

{\rho }_{\mathrm{l}} u_{\mathrm{l}} \left(\mathrm{1}-\alpha \right) + {\rho }_{\mathrm{g}} u_{\mathrm{g}} \alpha = \dot{m}，(5)

{\rho }_{\mathrm{l}} u_{\mathrm{l}} \left(\mathrm{1}-\alpha \right) x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}} + {\rho }_{\mathrm{g}} u_{\mathrm{g}} \alpha y_{\mathrm{g}}^{\mathrm{C}\mathrm{O}\mathrm{2}} = \gamma \dot{m},(6)

\frac{dP}{dz} = - {\rho }_{\mathrm{m}} g - \frac{{\rho }_{\mathrm{m}}f}{\mathrm{2}d} u_{\mathrm{m}}^{\mathrm{2}},(7)

\frac{d{h}_{\mathrm{m}}}{dz} = - g - \frac{q}{\dot{m}},(8)

其中： \dot{m}为总的质量速度，单位kg/s，为已知量； \gamma为CO₂的质量分数，也为已知量； h_{\mathrm{m}} = \left(\mathrm{1}-x\right) h_{\mathrm{l}} + x h_{\mathrm{g}}为混合物的热焓，单位J/kg， x为气相质量分数.

单相时方程（5）~（8）中有效方程为（5）、（7）和（8）三个，待求变量选择P、 h和 u三个，其他参数均可为这三个变量的函数；两相时未知数有 {\rho }_{\mathrm{l}}、 {\rho }_{\mathrm{g}}、 u_{\mathrm{l}}、 u_{\mathrm{g}}、 \alpha、 x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}、 y_{\mathrm{g}}^{\mathrm{C}\mathrm{O}\mathrm{2}}、 P、 h_{\mathrm{l}}、 h_{\mathrm{g}}和 x，还需要状态辅助方程.两相时，水的分压是温度的函数，即：

P_{\mathrm{H}\mathrm{2}\mathrm{O}} = P^{\mathrm{s}\mathrm{a}\mathrm{t}} ( T ).(9)

CO₂在液相中的质量分数 x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}是CO₂气相中分压和温度的函数，即：

P_{\mathrm{C}\mathrm{O}\mathrm{2}} = H ( T ) x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}.(10)

液相密度 {\rho }_{\mathrm{l}}是总压 P、温度 T和液相中CO₂质量分数 x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}的函数；气相密度 {\rho }_{\mathrm{g}}和气相CO₂质量分数 y_{\mathrm{g}}^{\mathrm{C}\mathrm{O}\mathrm{2}}均为温度 T、H₂O的分压 P_{\mathrm{H}\mathrm{2}\mathrm{O}}和CO₂的分压 P_{\mathrm{C}\mathrm{O}\mathrm{2}}的函数；液相混合物的热焓 h_{\mathrm{l}}为：

\begin{array}{l}{h}_{\mathrm{l}}=\left(u_{\mathrm{l}\_\mathrm{H}\mathrm{2}\mathrm{O}}+\frac{P}{{\rho }_{\mathrm{l}\_\mathrm{H}\mathrm{2}\mathrm{O}}}\right)\left(\mathrm{1}-x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}\right)+\\ (h_{\mathrm{C}\mathrm{O}\mathrm{2}}+h_{\mathrm{S}\mathrm{O}\mathrm{L}})x_{\mathrm{l}}^{\mathrm{C}\mathrm{O}\mathrm{2}}\end{array},(11)

其中： u_{\mathrm{l}\_\mathrm{H}\mathrm{2}\mathrm{O}}和 {\rho }_{\mathrm{l}\_\mathrm{H}\mathrm{2}\mathrm{O}}分别为压力 P和温度 T下的液相H₂O的内能和密度； h_{\mathrm{C}\mathrm{O}\mathrm{2}}为分压 P_{\mathrm{C}\mathrm{O}\mathrm{2}}和温度 T下的CO₂的热焓； h_{\mathrm{S}\mathrm{O}\mathrm{L}}为CO₂溶解于水放热量.气相混合物的热焓 h_{\mathrm{g}}为：

h_{\mathrm{g}} = \left(u_{\mathrm{g}\_\mathrm{H}\mathrm{2}\mathrm{O}}+\frac{P_{\mathrm{H}\mathrm{2}\mathrm{O}}}{{\rho }_{\mathrm{g}\_\mathrm{H}\mathrm{2}\mathrm{O}}}\right) \left(\mathrm{1}-y_{\mathrm{g}}^{\mathrm{C}\mathrm{O}\mathrm{2}}\right) + h_{\mathrm{C}\mathrm{O}\mathrm{2}} y_{\mathrm{g}}^{\mathrm{C}\mathrm{O}\mathrm{2}},(12)

其中： u_{\mathrm{g}\_\mathrm{H}\mathrm{2}\mathrm{O}}和 {\rho }_{\mathrm{g}\_\mathrm{H}\mathrm{2}\mathrm{O}}分别为压力 P_{\mathrm{H}\mathrm{2}\mathrm{O}}和温度 T下的气相H₂O的内能和密度.

气相体积分数 \alpha在两相计算中是一个至关重要的参数，根据定义，其为：

\alpha = [ \mathrm{1} + ( \frac{\mathrm{1}-x}{x} ) ( \frac{{\rho }_{\mathrm{g}}}{{\rho }_{\mathrm{l}}} ) {(\frac{u_{\mathrm{g}}}{u_{\mathrm{l}}})]}^{-\mathrm{1}},(13)

对于速度非均质模型（即 u_{\mathrm{g}} \ne u_{\mathrm{l}}），还需要建立气相体积分数 \alpha与速度的关系.最为常用的模型是最早由Zuber and Findlay（1965）建立的漂移流模型（Drift flux model，简写DFM）：

u_{\mathrm{g}} = C_{\mathrm{0}} j + u_{\mathrm{d}},(14)

其中：体积平均速度 j = \alpha u_{\mathrm{g}} + ( \mathrm{1} - \alpha ) u_{\mathrm{l}}； u_{\mathrm{d}}为气相漂移速度， C_{\mathrm{0}}为分布参数，受局部气相饱和度、管道的横向速度分布、流态、井斜和流动方向的影响.联合（13）和（14）最终可以得到：

\alpha = \frac{{\rho }_{\mathrm{g}}}{x} [ C_{\mathrm{0}} \left(\frac{\mathrm{1}-x}{{\rho }_{\mathrm{l}}}+\frac{x}{{\rho }_{\mathrm{g}}}\right) + \frac{u_{\mathrm{d}}}{\dot{m}/A} ].(15)

分布参数和气相漂移速度对气相体积分数有决定性的作用，并最终影响压力和温度的分布.虽然DFM模型有很多的形式，但是本文仅考虑了Shi et al.（2005）模型，主要原因有：（1）该模型最大的优势是把不同流态样式（比如泡状流、弹状流、环状流和滴状流）统一到一个模型中，避免了流态转化时因采用不同模型引起的不连续性；（2）通过对比各种模型（为了节约篇幅，本文未讨论），结果显示采用该模型计算的结果与实测的误差最小.Shi et al.（2005）根据实验数据，优化了参数 C_{\mathrm{0}}和 u_{\mathrm{d}}，其模型为：

C_{\mathrm{0}} = \frac{C_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\mathrm{1}+(C_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{1}){\eta }^{\mathrm{2}}},(16)

u_{\mathrm{d}} = \frac{(\mathrm{1}-\alpha C_{\mathrm{0}})K{u}_{\mathrm{c}}}{\alpha C_{\mathrm{0}}\sqrt[]{\frac{{\rho }_{\mathrm{g}}}{{\rho }_{\mathrm{l}}}}+\mathrm{1}-\alpha C_{\mathrm{0}}},(17)

其中： \mathrm{参数} \eta = \frac{\beta -B}{\mathrm{1}-B}， \beta = \mathrm{m} \mathrm{a} \mathrm{x} ( \alpha , \frac{\alpha |u_{\mathrm{m}}|}{u_{\mathrm{s}\mathrm{g}\mathrm{f}}} )， u_{\mathrm{s}\mathrm{g}\mathrm{f}} = K_{\mathrm{u}} {(\frac{{\rho }_{\mathrm{l}}}{{\rho }_{\mathrm{g}}})}^{\mathrm{0.5}} u_{\mathrm{c}}， B = \frac{\mathrm{2}}{C_{\mathrm{m}\mathrm{a}\mathrm{x}}} - \mathrm{1.066} \mathrm{7}， u_{\mathrm{c}} = {[\frac{g\sigma ({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}})}{{{\rho }_{\mathrm{l}}}^{\mathrm{2}}}]}^{\mathrm{0.25}}， K_{\mathrm{u}} = [ \frac{\mathrm{142}}{\sqrt[]{N_{\mathrm{B}}}} {(\sqrt[]{\mathrm{1}+\mathrm{0.006}\mathrm{2}{N}_{\mathrm{B}}}-\mathrm{1})]}^{\mathrm{0.5}}， N_{\mathrm{B}} = d^{\mathrm{2}} \left[\frac{g\left({\rho }_{\mathrm{l}}-{\rho }_{\mathrm{g}}\right)}{\sigma }\right]；参数 K是气相体积分数 \alpha的函数，即：

\begin{array}{l}K=\\ \left\{\begin{array}{c}\mathrm{1.53}\alpha \le a_{\mathrm{1}}\\ \mathrm{1.53}+\frac{C_{\mathrm{0}}{K}_{\mathrm{u}}-\mathrm{1.53}}{\mathrm{2}}\left[\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\left(\pi \frac{\alpha -a_{\mathrm{1}}}{a_{\mathrm{2}}-a_{\mathrm{1}}}\right)\right]a_{\mathrm{1}}\le \alpha \le a_{\mathrm{2}}\\ C_{\mathrm{0}}{K}_{\mathrm{u}}\alpha \ge a_{\mathrm{2}}\end{array}\right.\end{array}，(18)

其中：A、a ₁和a ₂均为经验参数，本论文中根据Shi et al.（2005）的建议，取C _max=1.2，a ₁=0.06，a ₂=0.12.

1.2　计算流程

对于方程（1）~（4）的求解，最严格的方法是把井筒进行网格剖分，然后建立离散方程，采用Newton-Raphson迭代方法求解，该方法已被成功应用到井筒‒储层耦合模拟器T2Well（Pan and Oldenburg， 2014）中.但是，采用该方法将遇到两个难点：（1）对于井筒两相，初始条件非常难于设定，设置不同的初始条件可能得到不同的结果；（2）迭代收敛性方面，特别是相变区域，经常会计算时收敛失败，无法得到计算结果.因此，参考文献（Barelli et al.， 1982）给出的方法，给出了以下改进的计算流程（图1）.

1.3　可靠性验证

为了验证程序的可靠性，利用本程序计算的结果与文献公开报道的结果进行对比.第一个对比实例是Khasani et al.（2021）给出的模型和计算结果，具体模型和参数见文献，模型中考虑的CO₂质量分数是1%.从图2a~2c可以看到，整体对比结果一致，但本文计算的压力比文献中的稍微偏低，温度在井筒上部有些偏差，气相中的CO₂质量分数在底部刚开始闪蒸的位置附近有些偏差.这些偏差的原因主要有：（1）本文与文献描述CO₂在水中的溶解度模型有差异，（2）水相和气相的速度差模型不同.第二个对比实例是与T2Well-EWASG（Vasini et al.， 2018）计算结果的对比.从图2d和2e可以看到，对比结果非常一致，仅温度有些偏差，主要原因是模型中本文忽略了围岩换热过程.

2 羊八井ZK304井两相流动评价

2.1　地热田地质概况

我国高温地热田主要集中在川藏地区（汪新伟等， 2022），其中的羊八井地热田是我国最早开发的高温地热田.羊八井地热田位于西藏拉萨市西北约90 km处，念青唐古拉山山前断陷盆地的西南端，海拔高程在4 290~4 500 m，地势上具有西北高、东南低的特征.热田内构造活动强烈，发育着北东、北西和近南北向三组断裂.以中尼公路为界，浅部热田可分为南、北两区：南区是第四系孔隙型热储，北区由第四系孔隙型和喜山期花岗岩裂隙型两类热储组成.

浅部热储埋深在地表以下180~280 m，岩性由第四系冲洪积砂砾石层、冰碛砂砾石、基岩顶部花岗岩风化壳组成.热储顶部由厚度不等的泥砾层或粉砂质粘土层构成盖层.热储底部基岩为喜马拉雅山早期花岗岩和凝灰岩，热田北部局部见有糜棱岩花岗岩.浅层热储流体温度一般在140~160 ℃，流体主要为液态，气体含量较少，水质以Cl·HCO₃-Na型为主，属深部流体与浅表冷水混合产物.

大地电磁探测结果揭示羊八井地热田北区深部约5 km以上存在电阻率为5 Ω·m的低阻层，推断为正在冷却的岩浆房，其为地热田的热源.北区山前发育有断层破碎带，局部破碎带宽达1 km，浅部正断层和深部走滑断层有效沟通了深部热源和浅部的环境.热田地下水氢、氧同位素分析结果表明，本区地下热水是念青唐古拉山的大量冰雪融水沿断裂带渗入地下经深循环，并被深部岩浆热源加热对流的结果（图3）.

地热田地质勘查工作开始于1976年，在勘查阶段开钻45口勘探井，圈定了42 km²的热田面积，并于20世纪80年代开始陆续开钻25口生产井和13口回灌井.自1977年9月羊八井第一台1 MW试验机组发电成功，到1991年装机容量已扩大到25.18 MW.截至2017年底，地热田南区有4口生产井，北区有12口生产井，5口回灌井，装机规模24.18 MW，累计发电超过33×10⁸ kW·h.水汽开采速率高峰期约为 12 000 t/d，大部分生产井的初始工作温度为125~140 ℃，工作压力为1.76~3.72 bar，生产井压力和温度2012年以前呈现逐年降低的趋势，2012年回灌常态化后压力和温度都逐渐恢复（许多龙等，2018）.

2.2　地热井测试

为了分析井筒的流动情况，2020年4月-5月在北区典型井ZK304中进行了放喷测试.ZK304井深度为206.54 m，埋深0~142 m为第四系孔隙型热储层，埋深142 m以下为基岩裂隙热储层.该井机械捅井除垢深度为117.17 m，结垢约在80~ 130 m之间.测试前探棒探井，下至124.00 m遇阻.

测试主要采用詹姆斯端压法，测试工具为Kuster K10 PT Geothermal Instrument.首先进行静态PT测试，然后引喷，在持续稳定喷发状态下由井底往上进行PT测试，同时在井口进行流量测试.测试结果见图4.从图中可以看到，前后两次测试条件下，井底的温度和压力差异不大，但是往上差别逐渐增大.在静态状态下，水位低于井口高程，水位在埋深50~60 m之间.动态测试条件下，压力由井底的7.3 bar逐渐降低到井口的2.7 bar.受地表环境、热储层和围岩的影响，温度有较大的变化.静态状态下，水位埋深以下温度主要受周围储层热交换的影响表现出由井底往上温度逐渐降低；在水位埋深位置到埋深约10 m的位置，温度稳定在86 ℃左右，与水在当地大气压的蒸发温度一致，结合压力可以看到这一段主要为蒸发的水蒸汽；受地表气温的影响，井口温度急剧下降到约20 ℃.在动态状态下，温度由井底的154 ℃逐渐下降到埋深10 m的130 ℃；

为了进一步分析动态和静态条件下由井底到井口的相态特征，把温度和压力投影到水的相图上.从图5可以看到，在动态测试下，流体在井底处于单相液体状态，在上升过程中，由于克服重力压力逐渐降低，降低到对应温度的饱和压力时开始蒸发，一部分液态水变成气态，并伴随着温度的下降，整个体系处于两相状态.在两相状态下，温度和压力将始终沿着饱和线移动.由于地表环境的影响，测试过程中流体在井口发生强烈的冷凝作用，导致流体状态回到单相液态.静态测试下，由井底到水位面均为单相液体状态，然后往上一段处于水蒸汽状态，到井口受地表环境影响冷凝为液体.

虽然相图显示动态测试条件下流体的相变过程与单组分水的饱和线有很高的一致性，推断流体中CO₂的含量并不高，但是在井口的蒸汽组分测试中显示仍然有一定量的CO₂，其对碳酸钙结垢起到不可忽略的影响（表1），同时多年的实际运行情况显示井筒中出现了大量的碳酸钙结垢物（周大吉，2003），而其出现与CO₂含量有很大的关系.因此，在两相流动评价中，需要考虑CO₂的存在.

2.3　模拟评价

2.3.1　模型建立

用于模型计算的井筒几何参数和井底边界条件见表2.由于压力和温度测试阶段没有获取井底保真样品，套管底部最为关键的参数之一的CO₂质量分数未知.但是，根据以前机械除垢情况，该井的闪蒸深度在套管底部以上，也就是说明套管底部位置均为单相状态，由此可以计算得到井套管底部位置CO₂的质量分数应不大于0.08%.为了分析CO₂的不确定性，笔者进行了3种情况（情况1：饱和CO₂含量（0.08%），情况2：饱和CO₂含量的一半（0.04%），情况3：没有CO₂）的模拟.由于流量相对较大和井筒长度相对较短，围岩换热的影响较小，在模型中忽略围岩换热过程.

2.3.2　压力和温度拟合结果

井筒中能够调整的参数是用于描述两相速度差的漂移流模型和相关参数.速度均质模型（即 u_{\mathrm{g}} = u_{\mathrm{l}}）常用于两相流动评价（Lu and Connell， 2014a， 2014b； Lei et al.， 2020）.为了验证速度均质模型是否有效，本文也采用该模型（即令 C_{\mathrm{0}} = \mathrm{1}， u_{\mathrm{d}} = \mathrm{0}）计算的结果与实际测量对比.从图6a和6b可以看到，采用速度均质模型计算的温度和压力与实际测量的偏差较大.而采用速度非均质模型（DFM）整体上拟合度非常好，仅井口温度有较大的差异，主要原因是井口流体与大气环境发生了较强的热交换，使得测量温度明显偏低.同时也可以看到，流体中CO₂含量越高，温度和压力也越高.但是，因为井底的CO₂含量较低，其对温度和压力整体分布的影响基本可以忽略不计.

从温度的变化（图6d）可以看到，在埋深120 m到110 m左右，温度变化非常小，往上开始减小，有明显的转折点.转折点的出现，预示着闪蒸的开始.这也进一步说明井底位置处没有CO₂气相的存在，CO₂全部溶解在水中.

2.3.3　其他关键结果

从图7a可以看到，不考虑CO₂情况下，闪蒸深度约为97 m，而考虑质量分数0.04%的情况下，其约为108 m，可见CO₂的存在会增加闪蒸深度.由闪蒸深度往上，流体往上运动过程中克服重力会伴随压力的降低，导致液态水蒸发为气态，液态中CO₂也会解析出来进入到气态中，并同时伴随着温度的降低（因为蒸发吸热）.气相质量分数由闪蒸深度的0逐渐增加到井口的0.064.气相体积分数从闪蒸面开始快速增加到埋深约75 m的0.6左右，然后逐渐增加到井口的0.84（图7b）.根据式（15）看到，这个变化主要由气相质量分数、密度和两相速递非均质模型的参数共同控制.从图7c可以看到，在井底CO₂含量为0.04%的情况下，闪蒸深度位置CO₂初始析出占气相质量的0.27，往上随着水的蒸发，气相中CO₂的质量分数迅速降低，变化最为剧烈的是从闪蒸深度到埋深80 m左右位置，这也预示该井段可能是碳酸钙结垢最为严重的位置.计算的井口气相中CO₂质量分数为0.6%，大于测量的0.12%.主要原因是测量是在1个大气压下进行的，明显小于计算的井口压力，随着压力的降低，气相中水蒸气会增加，CO₂浓度会稀释降低.液相中CO₂的质量分数与气相中的变化趋势相同.从气相和液相速度上看（图7e和7f），井底单相时约为0.2 m/s，闪蒸面开始往上气相和液相速度都基本逐渐线性增加，到井口时气相达到约11 m/s，液相为约1.3 m/s，相差了近10倍.

3 结论

高温地热生产井碳酸钙结垢是全世界高温地热开发普遍遇到的一个问题，其评价涉及到井筒两相流动过程.为了深入理解这个问题，本文首先建立了相关的评价模型，然后在典型地热井——羊八井开展了放喷测试，并应用开发的模型评价了井筒两相流动，得到以下结论：

（1）根据质量守恒、能量守恒和动量守恒方程，建立了CO₂-H₂O体系一维垂直井筒相变稳定流动模型，提出了稳健的求解方法，并与相关文献对比验证了其计算结果的可靠性.

（2） 通过典型井的温度和压力测试，得到羊八井地热田浅层热储温度一般在150 ℃左右，压力在埋深120 m左右大概7~8 bar之间，这些井大部分不能自喷.在开采状态下，井口压力2~ 3 bar之间，井下10 m处的温度在130 ℃左右.

（3）压力和温度的历史拟合表明模型中液体和气体之间的速度差异模型（漂移流模型）对井筒中温度和压力的分布有决定性的影响，不考虑气相和液体之间的速度差，会使模型计算结果远远偏离测量值.

（4）模拟结果显示在不考虑CO₂的情况下，闪蒸深度为97 m，考虑套管底部位置对应温度下饱和CO₂的质量分数的一半时，闪蒸深度增加到107 m；井口计算的温度和压力分别约为128 ℃和2.6 bar；井口的气相质量分数在6%~7%之间，对应的井口气相饱和度约为0.84；从闪蒸点往上大概20~30 m气相和液相中CO₂质量分数变化较为剧烈，也是碳酸钙结垢严重井段.

本论文模型仅考虑了稳定流状态，未来将考虑非稳定状态，同时也需要耦合储层因素考虑储层开采衰竭对井筒流动动力学的改变.另外，最为重要的是要进行井底保真取样定量获取CO₂的含量.本文开发的模型可以耦合水文地球化学评价井筒结垢问题.

参考文献

[1] Barelli, A., Corsi, R., Del Pizzo, G., et al., 1982. A Two-Phase Flow Model for Geothermal Wells in the Presence of Non-Condensable Gas. Geothermics, 11(3): 175-191. https://doi.org/10.1016/0375-6505(82)90026-8

[2] Battistelli, A., Calore, C., Pruess, K., 1997. The Simulator TOUGH2/EWASG for Modelling Geothermal Reservoirs with Brines and Non-Condensible Gas. Geothermics, 26(4): 437-464. https://doi.org/10.1016/S0375-6505(97)00007-2

[3] Bertani, R., 2012. Geothermal Power Generation in the World 2005-2010 Update Report. Geothermics, 41: 1-29. https://doi.org/10.1016/j.geothermics.2011.10.001

[4] Bu, X. B., Guo, Z. P., Wang, L. B., 2021. Process Simulation of Geothermal Fluid Flow in Wellbore and Calcium Carbonate Scaling. Advances in New and Renewable Enengy, 9(5): 434-442 (in Chinese with English abstract).

[5] China Geological Survey of Natural Resources Ministry, Department of New and Renewable Energy of National Energy Administration, Chinese Academy of Science and Technolgy for Development, et al., 2018. China Geothermal Energy Development Report (2018). China Petrochemical Press, Beijing (in Chinese).

[6] Duo, J., 2003. The Basic Characteristics of the Yangbajing Geothermal Field: A Typical High Temperature Geothermal System. Engineering Science, 5(1): 42-47 (in Chinese with English abstract).

[7] Gunn, C., Freeston, D., 1991. An Integrated Steady-State Wellbore Simulation and Analysis Package. The 13th New Zealand Geothermal Workshop. Auckland.

[8] Guo, Q. H., Wang, Y. X., 2009. Trace Element Hydrochemistry Indicating Water Contamination in and around the Yangbajing Geothermal Field, Tibet, China. Bulletin of Environmental Contamination and Toxicology, 83(4): 608-613. https://doi.org/10.1007/s00128-009-9812-7

[9] He, Y. J., Liu, X., Xing, L. X., et al., 2021. Scaling Process Simulation and Anti-Scaling Measures of Karst Geothermal Field in Baoding of Hebei. Earth Science Frontiers, 29(4): 430-437 (in Chinese with English abstract).

[10] Khasani, Deendarlianto, Itoi, R., 2021. Numerical Study of the Effects of CO₂ Gas in Geothermal Water on the Fluid-Flow Characteristics in Production Wells. Engineering Applications of Computational Fluid Mechanics, 15(1): 111-129. https://doi.org/10.1080/19942060.2020.1862709

[11] Lei, H. W., Bai, B., Cui, Y. X., et al., 2023. Quantitative Assessment of Calcite Scaling of a High Temperature Geothermal Production Well: Hydrogeochemistry—Application to the Yangbajing Geothermal Fields, Tibet. Earth Science, 48(3): 935-945 (in Chinese with English abstract).

[12] Lei, H. W., Cai, Y. N., Lu, M., et al., 2020. A Study on the Thermal-Hydrodynamical-Coupled CO₂ Flow Process in the Ordos CCS-Geological-Formation. International Journal of Greenhouse Gas Control, 95: 102999. https://doi.org/10.1016/j.ijggc.2020.102999

[13] Lu, M., Connell, L. D., 2014a. The Transient Behaviour of CO₂ Flow with Phase Transition in Injection Wells during Geological Storage-Application to a Case Study. Journal of Petroleum Science and Engineering, 124: 7-18. https://doi.org/10.1016/j.petrol.2014.09.024

[14] Lu, M., Connell, L. D., 2014b. Transient, Thermal Wellbore Flow of Multispecies Carbon Dioxide Mixtures with Phase Transition during Geological Storage. International Journal of Multiphase Flow, 63: 82-92. https://doi.org/10.1016/j.ijmultiphaseflow.2014.04.002

[15] Lund, J. W., 2008. Development and Utilization of Geothermal Resources. In: Goswami, D. Y., Zhao, Y., eds., Proceedings of ISES World Congress 2007 (Vol. I-Vol. V). Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-540-75997-3_13

[16] Lund, J. W., Freeston, D. H., Boyd, T. L., 2011. Direct Utilization of Geothermal Energy 2010 Worldwide Review. Geothermics, 40(3): 159-180. https://doi.org/10.1016/j.geothermics.2011.07.004

[17] Pan, L. H., Oldenburg, C. M., 2014. T2Well—An Integrated Wellbore-Reservoir Simulator. Computers & Geosciences, 65: 46-55. https://doi.org/10.1016/j.cageo.2013.06.005

[18] Shi, H., Holmes, J. A., Durlofsky, L. J., et al., 2005. Drift-Flux Modeling of Two-Phase Flow in Wellbores. SPE Journal, 10(1): 24-33. https://doi.org/10.2118/84228-pa

[19] Vasini, E. M., Battistelli, A., Berry, P., et al., 2018. Interpretation of Production Tests in Geothermal Wells with T2Well-EWASG. Geothermics, 73: 158-167. https://doi.org/10.1016/j.geothermics.2017.06.005

[20] Wang, G. L., Zhang, W., Liang, J. Y., 2017. Evaluation of Geothermal Resources Potential in China. Acta Geoscientica Sinica, 38(4): 449-459 (in Chinese with English abstract).

[21] Wang, X. W., Wang, T. H., Gao, N. A., et al., 2022. Formation Mechanism and Development Potential of Geothermal Resource along the Sichuan-Tibet Railway. Earth Science, 47(3): 995-1011 (in Chinese with English abstract).

[22] Xu, D. L., Jing, T. Y., Tan, J. Q., 2018. Analysis of Production Wells Monitoring and Its Impacts in Yangbajing Geothermal Field. Sino-Global Energy, 23(12): 22-28 (In Chinese with English abstract).

[23] Zhang, L., Chen, S., Zhang, C., 2019. Geothermal Power Generation in China: Status and Prospects. Energy Science & Engineering, 7(5): 1428-1450. https://doi.org/10.1002/ese3.365

[24] Zhao, P., Duo, J., Liang, T. L., et al., 1998. Characteristics of Gas Geochemistry in Yangbajing Geothermal Field, Tibet. Chinese Science Bulletin, 43(21): 1770-1777. https://doi.org/10.1007/BF02883369

[25] Zhou, D. J., 2003. Operation, Problems and Countermeasures of Yangbajing Geothermal Power Station in Tibet. Electric Power Construction, 24(10): 1-3, 9 (In Chinese with English abstract).

[26] Zuber, N., Findlay, J. A., 1965. Average Volumetric Concentration in Two-Phase Flow Systems. Journal of Heat Transfer, 87(4): 453-468. https://doi.org/10.1115/1.3689137

[27] 卜宪标, 郭志鹏, 王令宝, 2021. 地热流体在井筒中的流动及碳酸钙结垢过程模拟. 新能源进展, 9(5): 434-442.

[28] 自然资源部中国地质调查局, 国家能源局新能源和可再生能源司, 中国科学院科技战略咨询研究院, 等, 2018. 中国地热能发展报告(2018). 北京: 中国石化出版社.

[29] 多吉, 2003. 典型高温地热系统: 羊八井热田基本特征. 中国工程科学, 5(1): 42-47.

[30] 何雨江, 刘肖, 邢林啸, 等, 2021. 河北保定岩溶地热结垢过程模拟及防垢对策. 地学前缘, 29(4): 430-437.

[31] 雷宏武, 白冰, 崔银祥, 等, 2023. 高温地热生产井碳酸钙结垢定量评价：水文地球化学——以西藏羊八井为例. 地球科学, 48(3): 935-945.

[32] 王贵玲, 张薇, 梁继运, 等, 2017. 中国地热资源潜力评价. 地球学报, 38(4): 449-459.

[33] 汪新伟, 王婷灏, 高楠安, 等, 2022. 川藏铁路沿线地热资源形成机理与开发潜力, 地球科学, 47(3): 995-1011.

[34] 许多龙, 荆铁亚, 谭金群, 2018. 羊八井热田生产井监测及变化影响分析. 中外能源, 23(12): 22-28.

[35] 周大吉, 2003. 西藏羊八井地热发电站的运行、问题及对策. 电力建设, 24(10): 1-3, 9.

基金资助

国家自然科学基金项目(41972316)

中国核工业集团有限公司集中研发项目“地热勘查开发利用关键技术研究”

四川省科技创新人才项目(2022JDRC0027)