0 引言
上述问题主要涉及井筒流动及结垢过程模拟。地热开发利用较早的国家有美国、土耳其、菲律宾、新西兰和意大利等,对于井筒流动和结垢模拟,积累了比较丰富的经验[9],包括二氧化碳在地热水中的溶解特性,碳酸钙在地热水中的溶解沉淀规律,井筒模拟,井筒两相流压降确定,结垢机理分析,闪蒸位置预测,防垢措施等[10,11,12]。我国地热开发利用较晚,针对上述问题积累较少。井口流型和气液含量,牵涉到设备容量和选型;闪蒸位置决定了下泵深度以及除垢措施的选择;热储特性的变化决定了开发年限和经济效益。为解决上述问题,数值模拟由于简单快捷、成本低,是一种可选的方法,而且通过多个现场验证的数学模型,可以用于指导生产。本文以较常见的碳酸钙结垢为例,开展井筒流动和结垢过程模拟,为解决上述问题提供一种简便易行的方法。
1 模拟计算依据
$\mathrm{Ca}^{2+}+2 \mathrm{HCO}_{3}^{-} \leftrightarrow \mathrm{CaCO}_{3}+\mathrm{CO}_{2} \uparrow+\mathrm{H}_{2} \mathrm{O}$ (1)
地热流体一般是从地下几千米深处采出,在井筒内从下往上的流动过程中,由于压力变化比较大,最容易出现闪蒸现象[16]。因此,开展井筒的流动模拟,既可以确定不同位置处的温度、压力、干度、流型,又可以确定闪蒸位置和井底参数甚至热储参数。
地热流体中不凝气体种类包括二氧化碳、氮气、氨气和硫化氢等,以二氧化碳为主。在计算时,假设不凝气体全部为二氧化碳。
1.1 质量和能量守恒
$\frac{\mathrm{d} M_{\mathrm{t}}}{\mathrm{d} L}=0$ (2)
液相中二氧化碳的溶解浓度按亨利定律给出,气相中二氧化碳的含量由分压确定[18]：
$n_{1}=\alpha(T) P_{\mathrm{c}}$ (3)
$n_{\mathrm{g}}=\frac{44 P_{\mathrm{c}}}{44 P_{\mathrm{c}}+18 P_{\mathrm{w}}}$ (4)
式（3）中的系数根据文献[19]计算确定。
在闪蒸点,流体的总压等于二氧化碳的分压与流体温度对应的饱和水蒸气压力之和[17]：
$P_{\mathrm{t}}=P_{\mathrm{c}}+P_{\mathrm{w}}$ (5)
$P_{\mathrm{wtb}}=P_{\mathrm{w}}+\frac{n_{1}}{\alpha(T)}$ (6)
流体的干度以及各项流量由下式给出[17]：
$x_{\mathrm{s}}=\frac{H_{\mathrm{wt}}-H_{\mathrm{w}}}{H_{\mathrm{s}}-H_{\mathrm{w}}}$ (7)
$M_{\mathrm{wt}}=M_{\mathrm{t}}(1-\gamma)$ (8)
$M_{\mathrm{s}}=M_{\mathrm{wt}} x_{\mathrm{s}}$ (9)
$M_{\mathrm{w}}=M_{\mathrm{wt}}\left(1-x_{\mathrm{s}}\right)$ (10)
$M_{\mathrm{ct}}=M_{\mathrm{t}} \gamma$ (11)
$M_{\mathrm{gc}}=\frac{n_{\mathrm{g}} M_{\mathrm{s}}}{1-n_{\mathrm{g}}}$ (12)
$M_{\mathrm{lc}}=\frac{n_{1} M_{\mathrm{w}}}{1-n_{1}}$ (13)
地热井运行一段时间后,流体向井筒周围岩石的传热基本达到平衡,相对于井口采出的热量,通过井筒向围岩散热造成的热损失很小,在计算时忽略不计。对于二氧化碳含量一定的地热流体,在各个不同深度的截面上气液相的总焓值相等[15]。
$\begin{aligned}G(\gamma, P, T)=& \frac{n_{1}-\gamma}{n_{\mathrm{g}}-n_{1}}\left[n_{\mathrm{g}} H_{\mathrm{gc}}+\left(1-n_{\mathrm{g}}\right) H_{\mathrm{gw}}\right]+\\& \frac{\gamma-n_{\mathrm{g}}}{n_{\mathrm{g}}-n_{1}}\left[n_{1} H_{1 \mathrm{c}}+\left(1-n_{1}\right) H_{1 \mathrm{w}}\right]\end{aligned}$ (14)
$G(\gamma, P, T)=\text { constant }=G\left(\gamma, P_{0}, T_{0}\right)$ (15)
1.2 单相流压降
$\Delta P_{\mathrm{t}}-\left(\Delta P_{\mathrm{a}}+\Delta P_{\mathrm{g}}+\Delta P_{\mathrm{f}}\right)=0$ (16)
$\Delta P_{\mathrm{g}}=\rho g \Delta L$ (17)
$\Delta P_{\mathrm{a}}=\left(\frac{M_{\mathrm{t}}}{A}\right)^{2}\left(v_{1}-v_{2}\right)$ (18)
$\Delta P_{\mathrm{f}}=\frac{\lambda \Delta L}{2 D \rho}\left(\frac{M_{\mathrm{t}}}{A}\right)^{2}$ (19)
1.3 两相流压降
总的体积流量：
$Q_{\mathrm{t}}=Q_{\mathrm{g}}+Q_{1}$ (20)
$V_{\mathrm{t}}=\frac{Q_{\mathrm{t}}}{A}$ (21)
$V_{\mathrm{gd}}=\frac{Q_{\mathrm{g}}}{A \sqrt[4]{\frac{\rho_{1}}{g \sigma}}}$ (22)
$L_{\mathrm{B}}=1.071-0.7277 \frac{V_{\mathrm{t}}^{2}}{D}$ (23)
式中：LB—bubble-slug转变边界,无量纲,LB ≥ 0.13。
$L_{\mathrm{S}}=50+36 V_{\mathrm{gd}} \frac{Q_{1}}{Q_{\mathrm{g}}}$ (24)
式中：LS—slug-transition转变边界,无量纲。
$L_{\mathrm{M}}=75+84\left(V_{\mathrm{gd}} \frac{Q_{1}}{Q_{\mathrm{g}}}\right)^{0.75}$ (25)
式中：LM—transition-mist转变边界,无量纲。
流型判断准则：
$ \frac{Q_{\mathrm{g}}}{Q_{\mathrm{t}}}<L_{\mathrm{B}} bubble $ (26)
$ \left.\begin{array}{l}\frac{Q_{\mathrm{g}}}{Q_{\mathrm{t}}}>L_{\mathrm{B}} \\V_{\mathrm{gd}}<L_{\mathrm{S}}\end{array}\right\} slug $ (27)
$ L_{\mathrm{S}}<V_{\mathrm{gd}}<L_{\mathrm{M}} transition $ (28)
$ V_{\mathrm{gd}}>L_{\mathrm{M}} mist $ (29)
假如流型是bubble,按如下公式计算
$F_{\mathrm{g}}=\frac{1}{2}\left[1+\frac{Q_{\mathrm{t}}}{V_{\mathrm{s}} A}-\sqrt{\left(1+\frac{Q_{\mathrm{t}}}{V_{\mathrm{s}} A}\right)^{2}-\frac{4 Q_{\mathrm{g}}}{V_{\mathrm{s}} A}}\right]$ (30)
$\bar{\rho}=\left(1-F_{\mathrm{g}}\right) \rho_{1}+F_{\mathrm{g}} \rho_{\mathrm{g}}$ (31)
$\tau_{\mathrm{f}}=\lambda \rho_{1} \frac{u_{1}^{2}}{2 D}$ (32)
$u_{1}=\frac{Q_{1}}{A_{1}}=\frac{Q_{1}}{A\left(1-F_{\mathrm{g}}\right)}$ (33)
假如流型是slug,按如下公式计算：
$R e=\frac{\rho_{1} V_{\mathrm{t}} D}{\mu_{1}}$ (34)
$R e_{\mathrm{b}}=\frac{\rho_{1} V_{\mathrm{b}} D}{\mu_{1}}$ (35)
假如$R e_{\mathrm{b}} \leq 3000$
$V_{\mathrm{b}}=\left(0.546+8.74 \times 10^{-6} \operatorname{Re}\right) \sqrt{g D}$ (36)
假如$3000<R e_{\mathrm{b}}<8000$
$V_{\mathrm{bi}}=\left(0.251+8.74 \times 10^{-6} \mathrm{Re}\right) \sqrt{g D}$ (37)
$V_{\mathrm{b}}=\frac{1}{2}\left(V_{\mathrm{b} \mathrm{i}}+\sqrt{V_{\mathrm{bi}}^{2}+\frac{11170 \mu_{1}}{\rho_{1} \sqrt{D}}}\right)$ (38)
假如$N_{\mathrm{b}} \geq 8000$
$V_{\mathrm{b}}=\left(0.35+8.74 \times 10^{-6} \operatorname{Re}\right) \sqrt{g D}$ (39)
$\bar{\rho}=\frac{M_{\mathrm{t}}+\rho_{1} A V_{\mathrm{b}}}{Q_{\mathrm{t}}+A V_{\mathrm{b}}}+\Gamma \rho_{1}$ (40)
$\tau_{\mathrm{f}}=\frac{\lambda \rho_{1} u_{\mathrm{t}}^{2}}{2 D}\left(\frac{Q_{1}+A V_{\mathrm{b}}}{Q_{\mathrm{t}}+A V_{\mathrm{b}}}+\Gamma\right)$ (41)
两个截面之间的压力损失为：
$\Delta P=\left(\frac{\bar{\rho} g+\tau_{\mathrm{f}}}{1-\frac{M_{\mathrm{t}} Q_{\mathrm{g}}}{A^{2} \bar{P}}}\right) \Delta L$ (42)
1.4 求解方法
上述方程有两种求解方法,一是已知井底参数（压力、温度、总质量流量、二氧化碳含量）,由井底逐步向井口方向求解,确定各个截面的参数,称为Down-Top算法[20];二是已知井口参数（压力、温度、干度、总流量、二氧化碳含量）,由井口逐步向井底和热储方向求解,称为Top-Down。本文以Down-Top算法为例,简要介绍求解方法。
（1）根据井筒长度L,将井筒从上到下等距离分成若干份,每份长度ΔL;
（2）根据流体温度计算水蒸气饱和蒸汽压,再根据式（3）,已知二氧化碳含量计算其分压;
（3）根据式（6）计算闪蒸开始的压力;
（4）根据式（16）~ 式（19）,计算从井底到闪蒸点的距离,确定闪蒸位置;
（5）对于两相流压降的计算,首先假设下一点的压力,根据假设的压力,通过式（14）和式（15）求解该点的温度;
（6）待温度求解平衡后,根据温度和假设的压力,利用式（7）~ 式（13）,更新干度和各相的流量以及物性;
（7）根据两相流的判定准则以及压降计算公式,计算假设点的压力;
（8）比较假设压力和计算的压力,如果偏差小于设定值,进入下一点继续计算;如果大于设定值,继续计算,直至满足要求;
（9）输出各个截面的压力、温度、干度、流量等参数。
如果是已知井口参数,计算步骤与上述类似,先从井口判断两相流流型,确定压降,逐步向下计算到单相流和闪蒸点。
1.5 模型验证
Table 1 Parameter list used for calculation表1 计算参数表 |
| 参数 | 单位 | 数值 |
|---|---|---|
| 井底压力 | bar | 75 |
| 井底温度 | ℃ | 260 |
| 质量流量 | kg/s | 50 |
| 井筒直径 | m | 0.2 |
| 井深 | m | 1 500 |
| 井管表面粗糙度 | m | 4.6×10-5 |
| 二氧化碳质量含量 | % | 0.5 |
1.6 热储特性模拟
热储的下降因子由下式计算：
$ C_{\mathrm{D}}=\frac{P_{\mathrm{R}}-P_{\mathrm{wb}}}{M}$ (43)
下降因子是热储的固有特性,可以通过两个稳态流动工况进行确定：调节井口压力分别为Pwh1和Pwh2,待系统流动平衡后,对应的流量和井底压力分别为M1、M2和Pwb1、Pwb2。根据井口成分和上述两相流理论,通过井口压力可以计算出对应的井底压力。则有：
$ C_{\mathrm{D}}=\frac{P_{\mathrm{R}}-P_{\mathrm{wb} 1}}{M_{1}}$ (44)
$ C_{\mathrm{D}}=\frac{P_{\mathrm{R}}-P_{\mathrm{wb} 2}}{M_{2}}$ (45)
整理上述方程可得：
$ C_{\mathrm{D}}=\frac{P_{\mathrm{wb} 2}-P_{\mathrm{wbl}}}{M_{1}-M_{2}}$ (46)
$ P_{\mathrm{R}}=\frac{M_{2} P_{\mathrm{wb} 1}-M_{1} P_{\mathrm{wb} 2}}{M_{2}-M_{1}}$ (47)
式（46）和式（47）中井底压力可由井口压力计算获得,因此,可以通过井口的压力和流量计算下降因子和热储压力。另外,下降因子由下述参数确定：
$ C_{\mathrm{D}}=\frac{\mu_{\mathrm{R}} \ln \left(\frac{r_{\mathrm{R}}}{r_{\mathrm{W}}}\right)}{2 \pi K L_{\mathrm{R}} \rho_{\mathrm{R}}}$ (48)
上述下降因子的表达式中,仅有热储渗透率、热储厚度和半径未知。因此,测出三个参数中的两个,可以求解另一个,或者测出一个,可以求出另两个的组合值。在不同的年份计算多个稳态工况下的下降因子和热储压力,可以画出二者随时间的变化规律,以此判断热储的衰减程度和后续开采能力。
2 结果及分析
影响井筒流动的因素众多,本文仅分析质量流量、二氧化碳含量、热储温度、热储压力、井管表面粗糙度、结垢厚度、结垢长度这七个因素。
基本条件包括热储压力25 MPa、温度120℃、井管直径177.8 mm、壁厚8.05 mm、管壁粗糙度0.045 7 mm、井深3 000 m、二氧化碳质量含量0.01 kg/kg、热储下降因子0.022 8 MPa/(kg•s)、质量流量200 t/h。
2.1 质量流量的影响
Fig. 1 (a) Relationship between mass flow rate and wellhead pressure; (b) variation of flashing depth with mass flow rate; (c) variation of wellhead temperature, CO2 partial pressure and dryness with mass flow rate图1 （a）质量流量与井口压力的关系;（b）质量流量与闪蒸深度的关系;（c）井口温度、二氧化碳分压和干度与质量流量的关系 |
图2反映了在质量流量为200 t/h时,流体温度、干度和压力随深度的变化。由于热储压力和下降因子不变,当流量不变时,井底压力也保持恒定。在流体从井底向井口流动的过程中,在未闪蒸前,流体温度不变,压力损失主要是重力和摩擦引起的。由图2a可知,从井底3 000 m深度,到闪蒸点深度794 m,流体压力线性降低,流体温度保持不变。在闪蒸点,流体总压小于二氧化碳分压与流体温度对应的饱和水蒸气压力之和,出现闪蒸汽化现象。在闪蒸点附近,流体温度和压力均出现较大变化：由于闪蒸汽化吸热,流体温度突然降低;由于气相密度降低,流体压降梯度较单相流变小,如图2a所示。自闪蒸点开始,流体从单相流进入两相流,流体的压力损失按照两相流计算。图2b中,流体的干度从闪蒸点到井口逐渐增加,越向上,干度随深度的变化梯度越大,这说明越向上,液体闪蒸为气体的量越大。液体闪蒸为气体的量越大,流体温度降低也越大,从图2a和图2c可知,干度变化和温度变化的斜率是一致的。
Fig. 2 (a) Variation of pressure and temperature with depth; (b) variation of dryness with depth; (c) variation of total pressure, partial pressure and temperature with depth图2 （a）流体压力和温度随深度的变化;（b）流体干度随深度的变化;（c）总压、分压和温度随深度的变化 |
2.2 二氧化碳含量
地热水因溶解了CO2而具有侵蚀性,能在热储中和固态碳酸钙岩发生溶蚀作用。同样,在地热流体运移过程中,由于压力减小（地热流体从热储开采到地面,温度是逐渐降低的,而温度的降低能促进碳酸钙的溶解,在本文中暂不考虑温度降低对碳酸钙结垢的影响）,CO2从液相逸出,造成了碳酸钙结垢。因此,研究CO2含量对流动和结垢的影响具有重要的意义。
Fig. 3 (a) Effect of CO2 content on wellhead dryness and flashing depth; (b) effect of CO2 content on wellhead temperature, pressure and CO2 partial pressure; (c) variation of temperature with depth and CO2 content图3 （a）CO2质量含量对井口干度和闪蒸深度的影响;（b）CO2质量含量对井口温度、压力和二氧化碳分压的影响;（c）不同CO2含量条件下,流体温度随深度的变化 |
2.3 热储温度和压力
由图4a可知,热储温度对井口压力和井口温度有较大影响,二者均随热储温度的升高而增大。闪蒸深度受热储温度的影响较小,当热储温度在90 ~ 150℃之间变化时,闪蒸深度在790 m附近变化不大。
Fig. 4 (a) Variation of wellhead temperature, pressure and flashing depth with reservoir temperature; (b) relationship of wellhead pressure and CO2 partial pressure with reservoir pressure; (c) effect of reservoir pressure on wellhead dryness, temperature and flashing depth图4 （a）井口温度、压力和闪蒸深度与热储温度的关系;（b）井口压力和CO2分压与热储压力的关系;（c）热储压力对井口干度、温度和闪蒸深度的影响 |
2.4 井管表面粗糙度
Fig. 5 (a) Relationship of wellhead pressure and flashing depth with the surface roughness of well tube; (b) relationship of wellhead temperature and dryness with the surface roughness of well tube图5 （a）井口压力和闪蒸深度与井管表面粗糙度的关系;（b）井口温度和干度与井管表面粗糙度的关系 |
2.5 结垢厚度和结垢长度的影响
Fig. 6 (a) Effect of scale thickness on wellhead pressure; (b) effect of scale length on wellhead pressure图6 （a）结垢厚度对井口压力的影响;（b）结垢长度对井口压力的影响 |
3 结论
在地热开发利用过程中需要探索的问题包括结垢位置和井口干度的确定,井口压力和流量的关系,结垢厚度和长度对生产的影响,以及如何通过井口参数反推热储特性等。如果能通过简单易行的方法解决上述问题,将对地热产业的发展起到巨大的促进作用。数学模拟是一种被证明能够模拟地热井筒流动和碳酸钙结垢过程的有效方法。本文利用地热流体从井底到井口的流动模型分析了井底、井口以及不同深度处各个参数的变化,得到如下结论：
（1）所建立的数学模型,既可以通过热储和井底参数推导井筒和井口的参数,又可以通过井口参数反推井筒和热储的参数。
（2）通过数学模拟,可以确定井底、井口以及不同深度处的温度、压力、二氧化碳分压和气液各组分的质量流量,并可以确定闪蒸压力和闪蒸深度。这些参数的确定,既可以指导设备的选型,又可以指导实际运行。
（3）对闪蒸和结垢深度有较大影响的参数主要是热储压力、流量和二氧化碳含量;所得到的井口压力和质量流量的关系,可以指导采暖季非稳态工况下的质量调节;通过两个稳态流量,可以反推热储压力和下降因子。
本文的数学模型仅是通过文献中的模拟结果进行了验证,下一步,本课题组将通过文献报道的实验数据以及国内地热井的现场运行数据验证数学模型。另外,根据井口参数反推井筒以及热储的参数,对生产具有重要的指导意义,下一步本课题组将开展相关的研究工作。
符号表：
Mt 总质量流量,kg/s
L 深度坐标,m
nl 液相中二氧化碳的质量含量
ng 气相中二氧化碳的质量含量
α(T) 亨利系数
Pc 二氧化碳分压,Pa
Pt 总压力,Pa
Pw 水蒸气分压,Pa
Pwtb 闪蒸开始压力,Pa
Hwt 未闪蒸前水的焓值,J/kg
Hs 两相区水蒸气的焓值,J/kg
Hw 两相区液体水的焓值,J/kg
xs 两相区的干度
$\gamma$ 单相区二氧化碳的质量含量
Mwt 单相区水的质量流量,kg/s
Ms 两相区水蒸气的质量流量,kg/s
Mw 两相区液态水的质量流量,kg/s
Mct 单相区二氧化碳的总质量流量,kg/s
Mgc 两相区气相二氧化碳质量流量,kg/s
Mlc 两相区液相二氧化碳质量流量,kg/s
G 气液相总的焓值J/kg
Hlc 液相中二氧化碳的焓值,J/kg
Hlw 液态水的焓值,J/kg
Hgc 气相中二氧化碳的焓值,J/kg
Hgw 气态水的焓值,J/kg
ΔPt 总压降,Pa
ΔPg 加速压降,Pa
ΔPh 重力压降,Pa
ΔPf 摩擦压降,Pa
v1, v2 流体的比容,m3/kg
g 重力加速度,m/s2
λ 摩擦因子
D 井管内径,m
A 井管截面积,m2
ρ 流体的密度,kg/ m3
Vt 总的流速,m/s
Vgd 无量纲气体流速
σ 液体的表面张力,N/m
Fg 气泡率,无量纲
Vs 气液速度差,m/s
τf 摩擦损失梯度,Pa/m
$\bar{\rho}$ 平均密度,kg/m3
Re 雷诺数
Reb 气泡雷诺数
Vb 气泡上升速度,m/s
μl 黏度,Pa∙s
Γ 液体分布系数,无量纲
ΔL 垂直方向的计算距离,m
ΔP 在ΔL距离上的压降,Pa
$\bar{P}$ 在ΔL距离上的平均压力,Pa
CD 下降因子,Pa/(kg/s)
PR 热储压力,Pa
Pwb 在流量为M时的井底压力,Pa
Pwh 井口压力,Pa
μR 热储流体黏度,Pa∙s
rR 热储半径,m
rW 井管半径,m
K 热储渗透率,m2
LR 热储厚度,m
ρR 热储流体密度,kg/m3