热补偿条件下双井EGS产能和寿命预测方法研究

张 霄，孙玉学，张庆松，李相辉，尹占超，李 壮，周 政

(山东大学 岩土与结构工程研究中心，山东 济南 250061)

近年来，化石资源临近枯竭，环境问题日益突出[1]。干热岩资源因储量巨大、清洁环保，决定了其必然会成为未来重要能源之一[2]。增强型地热系统（EGS）作为开发干热岩资源的有效手段，受到了广泛关注[3-5]。

研究EGS取热过程时，若利用建造干热岩电站的方法获得实际工程的地质资料、工程经验和数据，对研究干热岩储层热提取极为有利，但目前的干热岩储层改造技术尚不完善，投资巨大、周期长。因此，现阶段研究大多基于数值模拟方法建立数学模型，利用多场耦合过程对EGS取热过程进行评估，指导EGS的设计与优化[6-8]。在热储改造和EGS运行过程中，热储岩体的物理化学性质都会随时间和空间不断变化[9-10]，改造后热储尺寸多以km为量级，内部裂隙通常以mm为量级，在利用精细化网格直接模拟取热过程时计算量巨大[11]，需要对热储内的裂隙进行简化。目前，常用的简化模型有单孔隙模型[12]、双孔隙模型[13-14]和多孔隙模型。其中，双孔隙模型和多孔隙模型因精度要求较高，需要大量试验数据支撑，在实际中应用困难[15]。

干热岩热量资源预测方面，目前主要有两种预测方式：一种是利用目标热储体积及热储温度和某一确定温度（如地表温度）的差值计算，按百分比的可开采量预测[16-19]，主要有体积法、蒙特卡罗法及地表流量法等[20-23]，但这些方法在应用于EGS产能预测时，热储开采前后温差的设定值偏大，导致预测不准确；另一种是数值法[24]，根据详细的热储参数计算和评价勘查程度较高、具有一定开采历史的地热田[23]，此种方法较为繁琐，计算周期长，需要详细的热储参数，不适合未开采的热储产能的预测。在EGS运行系统可控因素的研究方面，目前的研究多集中于以不同的热媒介质注入温度、不同的热媒介质注入速度及不同井间距等人为可控因素控制EGS的运行[8,25]。以注水速度控制时，随着EGS的运行，注入井的压力不断产生变化，大大增加了EGS运行时的不可控性[26]。EGS取热效果方面，大多数研究将热储视为单孔隙模型，针对热储改造后的孔隙比、渗透率、注入井温度、布井方式等多种因素分析EGS取热效果[27-28]。如：王昌龙[29]考虑了11种因素对EGS取热效果的影响。在这些研究中，大多局限于单因素的定性研究，较少考虑多因素组合的定量研究，且鲜有将大地热流的热补偿作为初始条件加入到模型之中[30-33]。有研究发现，热储周围的大地热流对热储层长期积累的热补偿热量会对EGS运行寿命产生不可忽略的影响[34]。

以吸放热公式[35]和Dupuit公式为基础，建立EGS产能和寿命控制方程；依据产能和热储参数，通过调节热媒介质注入温度和注采压差两个可控因素取值预测EGS寿命，增加EGS运行时的可控性。同时，根据地下水流动方程和热传导方程分析出热储初始温度、热媒介质注入温度、EGS寿命、热储体积和热储比表面积5个因素对发生热突破（生产井产出温度低于423.16 K[36]）时对控制方程中热储减少热量、大地热补偿热量、生产井平均产出温度和热储形状系数4个未知参数的影响。在考虑大地热补偿的情况下，利用双井EGS的3维数值模型，对不同工况的EGS在发生热突破时的4个未知参数的变化规律进行研究，得到4个未知参数随5个影响因素变化的具体预测公式，实现EGS产能预测和多因素定量化分析。

1 流热耦合下EGS换热过程分析

如图1所示，假设热储周围围岩渗透率为0，热媒介质在流动过程中无热量损失，人工储层原有热量Qm和大地补偿的热量Qt完全由热媒介质携带并提取，其中，热媒介质携带热量由Qw表示。

图1 EGS换热过程的流热耦合机制Fig.1 Flow-heat coupling mechanism of EGS heat transfer process

根据热量平衡，由吸放热公式（式（1））结合Dupuit公式（式（2））可得热媒介质携带热量（EGS产出热量），如式（3）所示：

式（1）～（3）中：t为EGS寿命；k为热储渗透率；mw为热媒介质质量； ρw为热媒介质密度；L为热媒介质路径；A1为热媒介质在热储中流过的截面面积；cw为热媒介质的热容；vw为热媒介质平均流速；P为注采压差；ΔT1为注采温差，ΔT1=Tout-Tin。

若热储在tn时刻发生热突破，tm时刻的产出温度Tout(tm)可作为Tout(t)在t0至tn时刻的平均值，即：

整理式（3）和（4），可得EGS产能和寿命控制方程：

利用达西定律和多孔介质传热方程描述流热耦合过程，如式（6）～（9）所示：

连续性方程：

流体动量方程：

流体能量方程：

周围围岩热补偿能量方程：

牛顿冷却公式：

式（5）～（10）中，t为热媒介质流动时间，xi为热媒介质在i方向上的位移，ρw、ρm分别为流体、围岩的热储密度，vi为热储中流体的流速，Fi为体积力， ∇2为拉普拉斯算子，cw、cm为流体和围岩比热容，A2为热储与围岩接触面积，η为动力黏度，T为热储温度场，T∞为围岩初始温度，λw、λp分别为流体和围岩的导热系数，h为表面传热系数，Φ为热流量，A3为流体与固体接触面积，Ts、Tw分别为固体和流体温度。

利用式（6）～（9），结合牛顿冷却公式（10），对式（5）中的4个参数进行分析，可确定出影响4个参数的主要因素，如表1所示。

表1 4个未知变量的影响因素汇总Tab. 1 Summary of influencing factors of four unknown parameters

Qm由热储终温Tml和热储初始温度Tm控制，所以在预测时仅需考虑Tml即可，整理后如式（11）所示。

若利用式（5）达到讨论Tin和P与t之间关系的目的，可利用数值模拟预测4个参数的具体表达式，并代入式（5）中得到Tin、P、t三者关系。

2 EGS取热模型概述

2.1 地质模型构建

如图2所示，EGS地下部分由人工热储、热储周围的围岩、注入井和生产井组成。为了更准确地模拟EGS地下热交换过程，在模型中作出主要假设：1）热储裂隙小而均匀，可看做单孔隙率的多孔介质；2）热媒介质流动为单相流；3）热储中热媒介质流动是饱和层流流动（0.003 m/s以下），初始时刻热储层中充满了与周围岩石温度相同的热媒介质；4）周围围岩渗透率为0，热媒介质在流动过程中无损失；5）热储及周围围岩不会产生形变和各种化学反应。

图2 EGS地下部分概念图Fig.2 Concept picture of EGS underground part

图3 双井EGS模型3维图Fig.3 Three-dimensional pictures of the double well EGS model

如图3所示，热储周围包覆200 m厚的花岗岩围岩，在外表面施加以3 K/100 m随深度线性增加的温度边界，模拟大地热流给予热储的热量补给；分别建立不同边长的正六面体，研究体积对4个参数的影响，如图3（a）所示；分别建立同体积下不同形状的热储，研究比表面积对4个参数的影响，如图3（b）所示；其余模型边界条件如表2所示。模型采用“下注上采”的开采方式，最大模拟时间50 a，采用瞬态求解方式，计算结果每年采集一次，且最终结果取符合EGS运行25 a，生产井出水速率达到80 kg/s以上，保证EGS商业化利用[35]。

表2 模型边界条件取值Tab. 2 Values of model boundary conditions

2.2 模型热物性参数取值

模型其他热物性参数如表3所示[30-31]。

表3 模型热物性参数取值Tab. 3 Values of model thermal property parameters

2.3 网格独立性验证

本文主要对计算结果中的生产井平均产出温度和平均速度、热储平均温度和热储中流体平均速度进行独立性检验。

进行网格无关性检验时，按时间步长0.01 s，对不同自由四面体的网格数量进行无关性检验。采用热储边长为500 m正六面体，热储初始温度-热媒介质注入温度-注采压差为565.66 K-353.16 K-20 MPa的工况，分别选取第15年和第30年的计算结果进行比较（图4）。由图4可知：网格数量对计算结果有着较大的影响，当网格数量小于25万个时，计算结果存在较大差异；当网格数量大于30万个时，模拟结果变化不大。

图4 网格无关性检验结果Fig.4 Results of the grid independence verification

进行时间独立性检验时，以第12种（36.76万个）网格数量为网格模型，时间步长为0.007 5、0.010 0、0.012 5和0.015 0 s，其余条件与网格独立性检验时相同，结果如表4所示。由表4可知，当网格数量达到30万以后，时间步长在0.012 5 s以下时对模拟结果影响不大。

由图4和表4可知，在计算时，最大网格取40 m，最小网格取0.3 m，网格增长率取2.4，单元格数量41.21万个，时间步长选择0.01 s较为合理。

表4 时间独立性检验结果Tab. 4 Results of the time independence verification

3 结果分析

图5 Tm、Tin、V、SSa、t对Tml、Qt、、 α等参数的影响结果Fig.5 Influences of five factors (Tm、Tin、V、SSa、t) to four parameters (Tml、Qt、 Tout、 α)

由图5（a）、（c）、（d）、（e）可知：当Tm和Tin增加时，Tml逐渐上升，且Tm对Tml的影响要低于Tin对Tml的影响；随着SSa的增加，Tml也会逐渐升高；随着V的增大，Tml表现出逐渐减小的趋势；以Tm=565.66 K、Tin=293.16 K为例，EGS不同寿命所对应的Tml相差不大，近似服从高斯分布，与时间项无关。由于热突破条件并非与热媒介质注入温度相同，所以在注入井与生产井之间必定存在温度梯度，在均匀介质渗流场中，两井之间必存在流速最快的最优通路。在最优通路间，温度梯度介于注入温度和热突破温度之间，在最优通路外其温度梯度上限与热储初始温度有关。因此，随着Tm的升高，最优通路外的热储温度升高；随着Tin的升高，整个热储温度梯度下限升高；随着SSa的增加，热储与围岩接触面积增加，热补偿的热量增多，都会导致Tml的升高。此外，其他条件相同时，V的增加使得热突破时低温区域增大，Tml逐渐减小。

由图5（f）、（g）、（h） 可知：受围岩热补偿的影响，25年后 ，Qt随着EGS运行时间逐渐增加；随着Tm、V和SSa的增大，Qt也逐渐增加；对于Tin对Qt的影响，由图5（f）分析可以发现，除热储初始温度处于485.66 K时，Qt随热媒介质注入温度逐渐增加，其余工况均无明显规律。分别取第27年和第29年Qt的值，利用极差分析法将Tm与Tin对Qt的影响进行比较，如表5和6所示，发现热储初始温度Tm对Qt的影响要远远大于Tin对Qt的影响，可近似将Tin对Qt的影响忽略。因此，可以确定，Qt主要受到Tm、Tin、V和SSa的影响。

表5 第27年Tm与Tin极差分析Tab. 5 Range analysis of Tm and Tin in the 27th year

表6 第29年Tm与Tin极差分析Tab. 6 Range analysis of Tm and Tin in the 29th year

由图5（g）、（h）可知：随着V的增大，α表现出先减小后增大的趋势，其最低点大致出现在边长为400 m的正六面体热储附近；随着SSa的增大，α表现出先增大后减小的趋势，其最高点出现在正六面体热储附近。

经上述分析，将表1中的影响因素予以修正，如表7所示，同时将式（11）修正为式（12）：

表7 修正后的Tml、Qt、、α影响因素汇总Tab. 7 Summary of revised influencing factors of Tml, Qt,,α

表7 修正后的Tml、Qt、、α影响因素汇总Tab. 7 Summary of revised influencing factors of Tml, Qt,,α

参数 影响因素Tml Tm、Tin、V、SSa Qt Tm、V、SSa、t Tout Tm、Tin、V、SSa α V、SSa

为得到式（12）的具体表现形式，需要对上述计算结果继续分析。

3.1 热储终温Tml 变化规律

为探究表7中4个因素与Tml的定量关系，以Tin为353.16 K时为基准，确定出Tml353.16随Tm的变化；然后，基于此分别讨论Tin、SSa与V对Tml的增量的变化，如图6所示；最后，将上述4项进行线性叠加，如式（13）所示：

如图6（a）所示，极易确定出Tml353.16随Tm变化公式。如图6（b）所示，不同Tm的ΔTmlTin随着Tin呈一致的线性负相关关系，因此，ΔTmlTin不受Tm的影响。取6种不同Tm的ΔTmlTin的平均值，得到ΔTmlTin随Tin的变化公式，如图6（c）所示。由于比表面积在不同体积时难以一致表示，因此，预测比表面积引起的热储终温增量ΔTmlS时，应以在同体积下与比表面积变化规律一致的量—各表面积与相应体积下正六面体表面积的比值S/SR表示，以消除热储体积对ΔTmlS的影响。如图6（d）所示，得到了ΔTmlS随S/SR的变化规律。考虑ΔTmlV随不同体积热储与边长为500 m的正六面体热储体积比值V/V500的变化规律，用以量化体积对Tml的影响，如图6（e）所示。所以得到式（13）中各分项的具体表达式为：

图6 热储终温Tml计算图Fig.6 Calculation of final storage temperature (Tml)

3.2 大地热补偿热量 Qt变化规律

首先，得到前25年的补偿热量随热储初始温度Tm的变化量Qt25，基于此分别讨论Qt随t和SSa的增量变化ΔQtt和ΔQtS，将上述3个量进行线性叠加。如图5（g）所示，当热储体积为0时，Qt也一定为0。令不同体积得到的补偿热量与边长为500 m的正六面体热储体积得到的补偿热量的比值为ηV，将ηV作为体积系数与Qt25、ΔQtt和ΔQtS叠加得到的补偿热量相乘，以考虑体积因素对Qt的影响，如式（15）所示：

Qt25与Tm的关系、ΔQtS与S/SR的关系及ηV与V/V500的关系如图7（a）、（b）和（d）所示。由图7（e）可知：ΔQtt随着时间的延长不断增加，近似呈线性变化；但不同Tm对应的ΔQtt的增量不同，且随着Tm的上升，ΔQtt随着时间变化的线性曲线斜率不断增加，截距不断减小。故在考虑ΔQtt随着时间的变化时，还要考虑Tm的影响。分别计算出斜率c和截距d随Tm的变化，如图7（c）所示。最终得到式（15）中各分项的具体表达式为：

图7 大地热补偿热量Qt计算Fig.7 Calculation of geothermal compensation heat (Qt)

3.3 平均产出温度 Tout变化规律

式中，Tout353.16、ΔToutS和ΔToutV极易通过数值计算结果通过拟合曲线得到，如图8（a）、（c）和（d）。如图8（e）所示，随着Tin的减小，ΔToutTin逐渐增加，可近似看做线性曲线，且不同Tm的ΔToutTin在同一注入温度的增量不同，随着注入温度的上升而不断增加。因此，在计算ΔToutTin时，需要考虑Tm的影响。与ΔQtt计算方法类似，如图8（b）所示，分别计算出ΔToutTin随Tin变化线性曲线的斜率a和截距b随Tm的变化，即可得到ΔToutTin的计算值，则式（17）中各分项的具体表达式为：

图8 平均产出温度 计算Fig.8 Calculation of average output temperature

3.4 系数 α的变化规律

由表7可知，系数α只与V和SSa有关。为得到三者间的定量关系，将α分为αV和ΔαS。首先，得到αV随V的变化，如图9（a）所示；然后，基于αV计算由SSa引起的增量ΔαS，如图9（b）所示；最后，将αV与ΔαS进行线性叠加得到系数α的计算式，如式（19）所示：

如图9所示，最终得到式（19）中各分项的具体表达式为：

图9 系数 α的计算图Fig.9 Calculation of coefficientsα

式（21）中各分式可由式（14）、（16）、（18）、（20）表示。将式（21）代入式（5）中，结合勘测到的热储渗透率即可得到Tin、P、t三者关系，指导EGS建设。

3.5 计算方程的适用性验证

目前，大多数研究主要利用工质循环流量表示热媒介质注入的量及热媒介质产出量，利用注采压差P与已有研究比较来验证计算公式适用性的方法难以直接实现。根据Dupuit公式可知，通过热储尺寸和渗透率可将注采压差P换算为热媒介质在生产井的平均流速，由生产井的尺寸即可得到热媒介质产出量，即工质循环流量。为验证计算方程的可靠性，将其应用至文献[15]的模型中。将文献[15]中的热储尺寸、热储初始温度、热媒介质注入温度和不同EGS运行时间代入式（5）中，得到不同EGS运行时间对应的注采压差P；然后，仅使用 COMSOL Multiphysics中达西定律模块计算得到工质循环流量，所得结果如图10所示。研究发现，利用式（5）计算结果与文献[15]所得曲线结果一致，说明式（5）具有良好的适用性。

图10 计算公式适用性验证Fig.10 Applicability verification of calculation formula

4 结 论

利用理论分析与数值模拟相结合，描述EGS地下取热过程，主要结论如下：

1）将吸放热公式与Dupuit公式相结合，得到了双井EGS产能与寿命控制方程，揭示了EGS产能随热储渗透率、EGS寿命、热媒介质注入温度、注采压差和热储形状系数的变化关系。

3）将控制方程和预测公式与现有研究进行对比，控制方程和预测公式具有良好的适用性，可有效指导EGS建设，加强EGS运行期间的实时调控。