CN110489911A - 一种考虑工质储存过程的增强型地热系统的数值模拟方法及应用 - Google Patents

Abstract

本发明公开了一种考虑工质储存过程的增强型地热系统的数值模拟方法及应用，其任意时刻的求解步骤为：将连续性方程中的非稳态项作为已知值，求解连续性方程和动量方程，计算此刻热储层中工质的速度场和压力场；计算此刻注入井井筒中工质的速度场、压力场和温度场；计算此刻热储层的温度场，并估算下一时刻热储层的温度场；计算此刻产出井井筒中工质的速度场、压力场和温度场；计算此刻热储层的平均总应力场，并估算下一时刻的平均总应力场；计算此刻热储层的孔隙率和渗透率分布，并估算下一时刻的孔隙率分布；判断各物理场是否收敛，如果未收敛，则重复上述步骤，如果收敛，则停止。本发明可同时模拟热开采性能及热储层中的工质储存量。

Description

一种考虑工质储存过程的增强型地热系统的数值模拟方法及 应用

技术领域

本发明涉及地热开采技术领域，具体的说是考虑工质储存过程的增强型地热系统(EGS)的数值模拟方法及应用，用于计算EGS的热开采性能和工质储存量。

背景技术

EGS是一种从低渗透率、低孔隙率的地层(主要是干热岩)中提取热量从而获取大量热能的系统，该系统采用人工激发手段在地层中建设高渗透率、高孔隙率的热储层，进而通过工质循环置换热储层中的热量，达到开采地热能的目的。由于干热岩温度较高、资源丰富，EGS具有开采大量热能的巨大应用前景，并受到了广泛关注。参看图2，EGS系统包括注入井1、产出井2和热储层4；注入井1及产出井2在热储层4的一部分开设有孔结构3，其开采地热能的原理为：冷工质被注入到注入井1中，通过注入井1井筒竖直向下流动到达注入井开孔，然后渗透到热储层4(为多孔介质，包括岩石及其空隙中的流体工质)并与热储层4进行热交换，最后变为热工质并经过产出井2回到地面。在开采地热能的同时，一部分工质被储存在热储层中，因为随着热开采的进行，热储层温度逐渐下降，导致岩石收缩(即热储层孔隙率增大)、工质密度增大，进而导致热储层中的工质的质量增大。然而，现有研究都很少考虑由热储层孔隙率、工质密度增大引起的工质储存过程。

2016年12月，《中国石油大学学报(自然科学版)》第40卷第6期第109-117页，孙致学、徐轶、吕抒桓、徐杨、孙强、蔡明玉和姚军公开了一篇名为“增强型地热系统热流固耦合模型及数值模拟”的文章，该文章将岩体视作由离散裂隙网络和基质岩块组成的双重介质，考虑传热、渗流与应力过程中的相互作用，建立复杂裂隙岩体的THM(传热、流动、力学)耦合数学模型，基于商业有限元软件COMSOL Multiphysics进行二次开发，实现数学模型的全耦合求解；针对饱和土热弹性固结问题，将耦合模型模拟结果与解析解进行对比，以验证模型的可靠性；最后，将该模型应用于简化的增强型地热系统，研究其温度场、渗流场以及应力场的变化规律。

2016年10月，《北京工业大学学报》第42卷第10期第1560-1564页，唐志伟、米倡华、张学峰和刘爱洁公开了一篇名为“增强型地热系统热固流耦合数值模拟与分析”的文章，该文章结合北京市活动断层探测实际，运用CFD软件结构性网格，基于导热微分方程和渗流流动方程，建立裂隙-断层岩体非稳态温度场的数学模型，描述了裂隙-断层水流及岩体温度场的分布，并结合计算参数和边界条件，利用裂隙-断层模型做了一定的延伸，对增强型地热系统热固流耦合传热过程进行数值模拟与分析。

2015年2月，《工程热物理学报》第36卷第2期第388-392页，罗良、曹文炅和蒋方明公开了一篇名为“增强型地热系统采热的分形分叉网络模型”的文章，该文章基于热储裂隙开度和长度的分形特征，建立起分形分叉网络模型，推导出采热速率表达式，量化分析裂隙形貌、分布及流体驱动泵功对采热速率的影响。

然而，上述三篇文献都未分析由热储层孔隙率、工质密度增大引起的工质储存过程，而工质储存过程可能对系统热开采性能有影响，并且工质储存量也是系统的重要性能参数：对于以水为工质的EGS(H₂O-EGS)，工质储存量意味着被消耗的水量，而为了使系统正常运行，必须要补充相应的水量；对于以CO₂为工质的EGS(CO₂-EGS)，工质储存可以实现CO₂的封存，工质储存量即为CO₂封存量。综上所述，现有的EGS数值模拟方法都有一定的局限性。因此，有必要探究同时模拟EGS热开采性能和工质储存量的新方法。

发明内容

1.要解决的技术问题

EGS系统在实现热开采的同时，还伴随着一定量的工质储存，并且工质储存与系统中的THM(传热、流动、力学)多场耦合密切相关，鉴于此，为了解决如何模拟EGS系统的热开采性能和工质储存量的问题，本发明提供了一种考虑工质储存过程的增强型地热系统的数值模拟方法及应用；本发明综合分析了注入井、产出井及热储层中的THM多场耦合过程，同时考虑了由热储层孔隙率、工质密度增大引起的工质储存过程，能够较为全面地模拟EGS系统。

2.技术方案

为了解决上述问题，本发明所采用的技术方案如下：

本发明的一种考虑工质储存过程的增强型地热系统的数值模拟方法，包括以下步骤：

步骤1：将热储层中工质的连续性方程中的非稳态项作为已知值，求解热储层中工质的连续性方程和动量方程，计算n时刻热储层中工质的速度场和压力场；

步骤2：基于注入压力、注入温度及步骤1计算的热储层与注入井井筒的交界处的工质速度，计算n时刻注入井井筒中工质的速度场、压力场和温度场；

步骤3：基于步骤1计算的热储层中工质的速度场及步骤2计算的热储层与注入井井筒的交界处的工质温度，计算n时刻热储层的温度场，并根据n时刻热储层的温度场估算(n+1)时刻热储层的温度场；

步骤4：基于产出压力、步骤1计算的热储层与产出井井筒的交界处的工质速度及步骤3计算的热储层与产出井井筒的交界处的工质温度，计算n时刻产出井井筒中工质的速度场、压力场和温度场；

步骤5：求解热储层的平均总应力方程，计算n时刻的平均总应力场，并根据所估算的(n+1)时刻热储层的温度场估算(n+1)时刻的平均总应力场；

步骤6：基于平均总应力场，计算n时刻热储层的孔隙率和渗透率分布，并根据所估算的(n+1)时刻的平均总应力场估算(n+1)时刻热储层的孔隙率分布；

步骤7：判断n时刻各物理场在经过i次迭代后是否收敛：如果各物理场还未收敛，则将(i+1)赋值给i，进入下一次迭代，并返回到步骤1；如果各物理场已经收敛，则计算n时刻的热开采速率及工质储存量，然后将(n+1)赋值给n，进入下一个时刻的计算，并令i＝1，返回到步骤1，直至完成所有时刻的计算。

本发明的一种考虑工质储存过程的增强型地热系统数值模拟方法的应用，该方法适用于以水为工质的EGS，也适用于以CO₂为工质的EGS。

3.有益效果

相比于现有技术，本发明的有益效果为：

(1)本发明的一种考虑工质储存过程的增强型地热系统的数值模拟方法，综合分析了注入井、产出井及热储层中的THM(传热、流动、力学)多场耦合过程，同时考虑了由热储层孔隙率、工质密度增大引起的工质储存过程，能够较为全面地模拟EGS系统；

(2)本发明的一种考虑工质储存过程的增强型地热系统的数值模拟方法，不仅可以模拟EGS系统的热开采性能，而且可以模拟系统的工质储存量(对于H₂O-EGS为被消耗的水量，而对于CO₂-EGS为CO₂封存量)。

附图说明

图1为本发明所提出的考虑工质储存过程的EGS系统的数值模拟方法的流程图；

图2为本发明EGS系统的示意图；

图3为本发明2个实施例中的EGS系统热储层的尺寸图；

图4为本发明实施例1所模拟的H₂O-EGS热开采速率与忽略工质储存过程情况下的H₂O-EGS热开采速率的对比图；

图5为本发明实施例1所模拟的H₂O-EGS工质储存量；

图6为本发明实施例2所模拟的CO₂-EGS热开采速率与忽略工质储存过程情况下的CO₂-EGS热开采速率的对比图；

图7为本发明实施例2所模拟的CO₂-EGS工质储存量。

示意图中的标号说明：

1、注入井；2、产出井；3、孔结构；4、热储层。

具体实施方式

下面结合具体实施例对本发明进一步进行描述。

本发明的一种考虑工质储存过程的增强型地热系统的数值模拟方法，考虑一个增强型地热系统(EGS)，其初始条件为：

T_res——热储层的温度；

T_sur——地表温度；

a——地温梯度；

z——直角坐标，在数值上等于深度；

σ_m——热储层的平均总应力；

ρ_ro——岩石的密度；

g——重力加速度；

p_f——工质的压力；

ρ_w——水的密度；

θ——积分变量；

t——时间。

对EGS系统进行网格划分，其中，将热储层划分为一定数量的尺寸不等的长方体网格，将注入井井筒及产出井井筒在深度方向上分别等分为一定数量的圆柱体网格；并设置所模拟的总时长及时间步长。

以n和i分别表示第n个时刻和第i次迭代，首先令n＝1，i＝1，则模拟该EGS系统的步骤如下：

步骤1：热储层中工质的连续性方程如下：

——热储层的孔隙率；

ρ_f——工质的密度；

——工质速度的矢量形式；

——微分算子；

将方程(2)中的非稳态项(即方程左边的第一项)作为已知值：若n＝1，则忽略方程中的非稳态项；若n>1，则由(n-1)时刻估算的n时刻的孔隙率分布和温度场计算得到非稳态项：

——(n-1)时刻估算的n时刻的热储层孔隙率；

——(n-1)时刻估算的n时刻的热储层温度所对应的工质密度；

——(n-1)时刻的热储层孔隙率；

——(n-1)时刻的工质密度；

tⁿ——n时刻对应的时间；

t^n-1——(n-1)时刻对应的时间。

热储层中工质的动量方程(也即上文所述工质速度的矢量形式)如下：

k——热储层的渗透率；

——重力加速度的矢量形式；

μ_f——工质的粘性系数。

热储层中工质的流动边界条件为：在热储层与注入井井筒的交界处，给定边界压力(若i＝1，则将上一时刻注入井井筒底部的工质压力赋值给该边界压力；若i>1，则将上一次迭代中步骤2计算的注入井井筒底部的工质压力赋值给该边界压力)；在热储层与产出井井筒的交界处，给定边界压力(若i＝1，则将上一时刻产出井井筒底部的工质压力赋值给该边界压力；若i>1，则将上一次迭代中步骤4计算的产出井井筒底部的工质压力赋值给该边界压力)；在其余边界，采用无滑移边界条件。

基于方程(2)、(3)和(4)及流动边界条件，采用SIMPLE算法计算n时刻热储层中工质的速度场和压力场。

步骤2：注入井井筒中工质的动量方程如下：

u——工质的速度；

τ——摩擦阻力，其计算式如下：

τ＝fρ_fu²/(4R_wel) (6)

R_wel——井筒的外半径；

f——摩擦系数，其满足的关系式如下：

f^-0.5＝-1.8log₁₀{[K/(7.4R_wel)]^1.11+6.9/Re} (7)

K——井筒的等效绝对粗糙度；

Re——雷诺数，Re＝2R_welρ_fu/μ_f；

注入井井筒中工质的能量方程为：

A_wel——井筒的流动面积，A_wel＝π(R_wel)²；

c_f——工质的比热容；

T_f——工质的温度；

dQ/dz——单位深度上井筒与周围地层之间的热流，其满足如下的方程(9)和(10)：

h_f——对流换热系数；

T_wel——井筒与周围地层交界处的温度；

λ_for——周围地层的热导率；

T_inf——无限远处的周围地层温度，T_inf＝T_sur+az；

f(t)——瞬态导热函数，其表达式如下：

ρ_wel——井筒的密度；

c_wel——井筒的比热容；

t_D——无量纲时间，其表达式为：

ρ_for——周围地层的密度；

c_for——周围地层的比热容；

注入井井筒中工质的流动边界条件及热边界条件为：在注入井井筒的顶部，给定注入压力及注入温度；在注入井井筒的底部，给定工质的速度(即步骤1计算的热储层与注入井井筒的交界处的工质速度)。

基于方程(5)、(6)、(7)、(8)、(9)、(10)、(11)和(12)及边界条件，采用迭代法计算n时刻注入井井筒中工质的速度场、压力场和温度场。

步骤3：在热储层中，假定岩石和工质之间达到瞬态热平衡，则热储层的能量方程为：

ρ_res——热储层的密度，

c_res——热储层的比热容，其中c_ro为岩石的比热容；

λ_res——热储层的热导率，其中λ_f和λ_ro分别为工质和岩石的热导率；

热储层的热边界条件为：在热储层与注入井井筒的交界处，给定边界温度(即步骤2计算的注入井井筒底部的工质温度)；在其余边界，采用绝热边界条件。

将步骤1计算的n时刻热储层中工质的速度场作为已知条件，基于方程(13)及热边界条件计算n时刻热储层的温度场，并根据n时刻热储层的温度场和热储层中工质的速度场估算(n+1)时刻热储层的温度场，其中，温度场的计算和估算均采用隐式迭代法。

步骤4：产出井井筒中工质的动量方程如下：

产出井井筒中工质的流动边界条件及热边界条件为：在产出井井筒的顶部，给定产出压力；在产出井井筒的底部，给定工质的速度(即步骤1计算的热储层与产出井井筒的交界处的工质速度)和工质的温度(即步骤3计算的热储层与产出井井筒的交界处的热储层温度)；

基于方程(6)、(7)、(8)、(9)、(10)、(11)、(12)和(14)及边界条件，采用迭代法计算n时刻产出井井筒中工质的速度场、压力场和温度场。

步骤5：热储层的平均总应力方程如下：

ν——岩石的泊松比；

α——Biot系数；

β——岩石的线膨胀系数；

E——岩石的杨氏模量；

热储层的应力边界条件为：在所有的边界上，给定平均总应力(等于初始值)。

将步骤1计算的n时刻热储层中工质的压力场及步骤3计算的n时刻热储层的温度场作为已知条件，并基于方程(15)及应力边界条件，采用迭代法计算n时刻的平均总应力场。

将步骤1计算的n时刻热储层中工质的压力场及步骤3估算的(n+1)时刻热储层的温度场作为已知条件，并基于方程(15)及应力边界条件，采用迭代法估算(n+1)时刻的平均总应力场。

步骤6：热储层的孔隙率关于平均总应力的关系式如下：

——初始时刻热储层的孔隙率；

C——热储层的孔隙压缩系数；

σ_m0——初始时刻热储层的平均总应力；

p_f0——初始时刻工质的压力；

热储层的渗透率关于孔隙率的关系式如下：

k₀——初始时刻热储层的渗透率。

基于步骤5计算的n时刻的平均总应力场，由方程(16)和(17)计算n时刻热储层的孔隙率和渗透率分布。

基于步骤5估算的(n+1)时刻的平均总应力场，由方程(16)估算(n+1)时刻热储层的孔隙率分布。

步骤7：判断n时刻各物理场在经过i次迭代后是否收敛：如果各物理场还未收敛，则将(i+1)赋值给i，进入下一次迭代，并返回到步骤1；如果各物理场已经收敛，则计算n时刻的热开采速率及工质储存量，然后将(n+1)赋值给n，进入下一个时刻的计算，并令i＝1，然后返回到步骤1，直至完成所有时刻的计算。

其中，热开采速率及工质储存量的计算式分别为：

Q_pro＝q_pro(h_pro-h_inj) (18)

Q_pro——热开采速率；

q_pro——产出流量，q_pro＝u_proA_welρ_f，其中u_pro为产出井井口的工质速度；

h_pro——产出井井口的工质比焓；

h_inj——注入井井口的工质比焓；

M_s——工质储存量；

q_inj——注入流量，q_inj＝u_injA_welρ_f，其中u_inj为注入井井口的工质速度。

实施例1

本实施例针对一个虚拟的H₂O-EGS，模拟其热开采速率与工质储存量，并将计算结果与忽略工质储存过程(即忽略方程(2)中的非稳态项)情况下的结果进行对比。

图1为本发明所提出的考虑工质储存过程的EGS系统的数值模拟方法的流程图，图2为EGS系统的示意图(为了便于分析，本实施例将注入井1及产出井2的孔结构3也当作热储层4的一部分)，图3为热储层的尺寸图，而本实施例中H₂O-EGS系统的相关参数如表1所示。工质的热物性参数为工质压力及温度的函数，可由PROPATH软件分析得到。

表1 H₂O-EGS系统的相关参数

对该H₂O-EGS系统进行网格划分，其中，将热储层划分为35×47×20个尺寸不等的长方体网格，将注入井井筒在深度方向上等分为3700个圆柱体网格，并将产出井井筒在深度方向上也等分为3700个圆柱体网格；模拟的总时长设为20年，并设置时间步长，其中前3个时间步长分别为10天、10天和80天，其余72个时间步长均为100天。

首先令时刻n＝1，迭代次数i＝1，则本实施例的具体计算步骤如下：

步骤1：若n＝1，则忽略方程(2)中的非稳态项；若n>1，则由方程(3)计算方程(2)中的非稳态项。

设置热储层中工质的流动边界条件：在热储层与注入井井筒的交界处，给定边界压力(若i＝1，则将上一时刻注入井井筒底部的工质压力赋值给该边界压力；若i>1，则将上一次迭代中步骤2计算的注入井井筒底部的工质压力赋值给该边界压力)；在热储层与产出井井筒的交界处，给定边界压力(若i＝1，则将上一时刻产出井井筒底部的工质压力赋值给该边界压力；若i>1，则将上一次迭代中步骤4计算的产出井井筒底部的工质压力赋值给该边界压力)；将其余边界设置为无滑移边界条件。

根据方程(2)、(3)和(4)及流动边界条件，采用SIMPLE算法计算n时刻热储层中各个网格的工质速度和压力。

步骤2：设置注入井井筒中工质的流动边界条件及热边界条件：在注入井井筒的顶部，给定注入压力及注入温度；在注入井井筒的底部，给定工质的速度(即步骤1计算的热储层与注入井井筒的交界处的工质速度)。

根据方程(5)、(6)、(7)、(8)、(9)、(10)、(11)和(12)及边界条件，计算n时刻注入井井筒中各个网格的工质速度、压力和温度。

步骤3：设置热储层的热边界条件：在热储层与注入井井筒的交界处，给定边界温度(即步骤2计算的注入井井筒底部的工质温度)；将其余边界设置为绝热边界条件。

将步骤1计算的n时刻热储层中工质的速度场作为已知条件，基于方程(13)及边界条件计算n时刻热储层中各个网格的温度，并根据n时刻热储层中各个网格的温度和速度估算(n+1)时刻热储层中各个网格的温度。

步骤4：设置产出井井筒中工质的流动边界条件及热边界条件：在产出井井筒的顶部，给定产出压力；在产出井井筒的底部，给定工质的速度(即步骤1计算的热储层与产出井井筒的交界处的工质速度)和工质的温度(即步骤3计算的热储层与产出井井筒的交界处的热储层温度)。

根据方程(6)、(7)、(8)、(9)、(10)、(11)、(12)和(14)及边界条件，计算n时刻产出井井筒中各个网格的工质速度、压力和温度。

步骤5：设置热储层的应力边界条件：在所有的边界上，给定平均总应力(等于初始值)。

将步骤1计算的n时刻热储层中工质的压力场及步骤3计算的n时刻热储层的温度场作为已知条件，基于方程(15)及应力边界条件计算n时刻热储层中各个网格的平均总应力。

将步骤1计算的n时刻热储层中工质的压力场及步骤3估算的(n+1)时刻热储层的温度场作为已知条件，基于方程(15)及应力边界条件估算(n+1)时刻热储层中各个网格的平均总应力。

步骤6：基于步骤5计算的n时刻的平均总应力场，由方程(16)和(17)计算n时刻热储层中各个网格的孔隙率和渗透率。

基于步骤5估算的(n+1)时刻的平均总应力场，由方程(16)估算(n+1)时刻热储层中各个网格的孔隙率。

步骤7：判断n时刻各物理场在经过i次迭代后是否收敛，收敛判据为：EGS系统中各个网格的温度在i次迭代与(i-1)迭代的差值均小于0.05K，EGS系统中各个网格的工质压力在i次迭代与(i-1)迭代的差值均小于10Pa，EGS系统中各个网格的工质速度在i次迭代与(i-1)迭代的差值均小于10^-8m/s，热储层中各个网格的平均总应力在i次迭代与(i-1)迭代的差值均小于20Pa。

如果各物理场还未收敛，则将(i+1)赋值给i，进入下一次迭代，并返回到步骤1；如果各物理场已经收敛，则计算n时刻的热开采速率及工质储存量，然后将(n+1)赋值给n，进入下一个时刻的计算，并令i＝1，然后返回到步骤1，直至完成所有时刻的计算。

本实施例的结果如图4和图5所示。图4为本实施例所模拟的H₂O-EGS热开采速率与忽略工质储存过程情况下的H₂O-EGS热开采速率的对比图，图5为本实施例所模拟的H₂O-EGS工质储存量。可以看出忽略工质储存过程将过高估计热开采速率，并且20年后的工质储存量约为0.6×10⁸kg。

实施例2

本实施例针对一个虚拟的CO₂-EGS，模拟其热开采速率及工质储存量，并将计算结果与忽略工质储存过程情况下的结果进行对比。

本实施例的具体步骤及EGS参数与实施例1相同，而不同之处在于CO₂与水的热物性参数不同，并且注入压力、注入温度及产出压力也不同。本实施例的注入压力、注入温度及产出压力分别为13×10⁶Pa、293K及13×10⁶Pa。

本实施例的结果如图6和图7所示。图6为本实施例所模拟的CO₂-EGS热开采速率与忽略工质储存过程情况下的CO₂-EGS热开采速率的对比图，图7为本实施例所模拟的CO₂-EGS工质储存量。结果也显示忽略工质储存过程将过高估计热开采速率，并且20年后的工质储存量约为2.2×10⁸kg。

Claims (10)

1.一种考虑工质储存过程的增强型地热系统的数值模拟方法，其特征在于，包括以下步骤：

步骤1：将热储层中工质的连续性方程中的非稳态项作为已知值，求解热储层中工质的连续性方程和动量方程，计算n时刻热储层中工质的速度场和压力场；

步骤2：基于注入压力、注入温度及步骤1计算的热储层与注入井井筒的交界处的工质速度，计算n时刻注入井井筒中工质的速度场、压力场和温度场；

步骤3：基于步骤1计算的热储层中工质的速度场及步骤2计算的热储层与注入井井筒的交界处的工质温度，计算n时刻热储层的温度场，并根据n时刻热储层的温度场估算(n+1)时刻热储层的温度场；

步骤4：基于产出压力、步骤1计算的热储层与产出井井筒的交界处的工质速度及步骤3计算的热储层与产出井井筒的交界处的工质温度，计算n时刻产出井井筒中工质的速度场、压力场和温度场；

步骤5：求解热储层的平均总应力方程，计算n时刻的平均总应力场，并根据所估算的(n+1)时刻热储层的温度场估算(n+1)时刻的平均总应力场；

步骤6：基于平均总应力场，计算n时刻热储层的孔隙率和渗透率分布，并根据所估算的(n+1)时刻的平均总应力场估算(n+1)时刻热储层的孔隙率分布；

步骤7：判断n时刻各物理场在经过i次迭代后是否收敛：如果各物理场还未收敛，则将(i+1)赋值给i，进入下一次迭代，并返回到步骤1；如果各物理场已经收敛，则计算n时刻的热开采速率及工质储存量，然后将(n+1)赋值给n，进入下一个时刻的计算，并令i＝1，返回到步骤1，直至完成所有时刻的计算。

2.根据权利要求1所述的一种考虑工质储存过程的增强型地热系统的数值模拟方法，其特征在于：增强型地热系统，其初始条件为：

T_res——热储层的温度；

T_sur——地表温度；

a——地温梯度；

z——直角坐标，在数值上等于深度；

σ_m——热储层的平均总应力；

ρ_ro——岩石的密度；

g——重力加速度；

p_f——工质的压力；

ρ_w——水的密度；

θ——积分变量；

t——时间；

对EGS系统进行网格划分，并设置所模拟的总时长及时间步长。

3.根据权利要求1所述的一种考虑工质储存过程的增强型地热系统的数值模拟方法，其特征在于：步骤1所述热储层中工质的连续性方程如下：

——热储层的孔隙率；

ρ_f——工质的密度；

——工质速度的矢量形式；

——微分算子；

将方程(2)中的非稳态项作为已知值：若n＝1，则忽略方程中的非稳态项；若n>1，则由(n-1)时刻估算的n时刻的孔隙率分布和温度场计算得到非稳态项：

——(n-1)时刻估算的n时刻的热储层孔隙率；

ρ_f ^*n——(n-1)时刻估算的n时刻的热储层温度所对应的工质密度；

——(n-1)时刻的热储层孔隙率；

ρ_f ^n-1——(n-1)时刻的工质密度；

tⁿ——n时刻对应的时间；

t^n-1——(n-1)时刻对应的时间；

热储层中工质的动量方程如下：

k——热储层的渗透率；

——重力加速度的矢量形式；

μ_f——工质的粘性系数；

热储层中工质的流动边界条件为：在热储层与注入井井筒的交界处，给定边界压力，即若i＝1，则将上一时刻注入井井筒底部的工质压力赋值给该边界压力；若i>1，则将上一次迭代中步骤2计算的注入井井筒底部的工质压力赋值给该边界压力；在热储层与产出井井筒的交界处，给定边界压力，即若i＝1，则将上一时刻产出井井筒底部的工质压力赋值给该边界压力；若i>1，则将上一次迭代中步骤4计算的产出井井筒底部的工质压力赋值给该边界压力；在其余边界，采用无滑移边界条件；

基于方程(2)、(3)和(4)及流动边界条件，采用SIMPLE算法计算n时刻热储层中工质的速度场和压力场。

4.根据权利要求1所述的一种考虑工质储存过程的增强型地热系统的数值模拟方法，其特征在于：步骤2中注入井井筒中工质的动量方程如下：

u——工质的速度；

τ——摩擦阻力，其计算式如下：

τ＝fρ_fu²/(4R_wel) (6)

R_wel——井筒的外半径；

f——摩擦系数，其满足的关系式如下：

f^-0.5＝-1.8log₁₀{[K/(7.4R_wel)]^1.11+6.9/Re} (7)

K——井筒的等效绝对粗糙度；

Re——雷诺数，Re＝2R_welρ_fu/μ_f；

注入井井筒中工质的能量方程为：

A_wel——井筒的流动面积，A_wel＝π(R_wel)²；

c_f——工质的比热容；

T_f——工质的温度；

dQ/dz——单位深度上井筒与周围地层之间的热流，其满足如下的方程(9)和(10)：

h_f——对流换热系数；

T_wel——井筒与周围地层交界处的温度；

λ_for——周围地层的热导率；

T_inf——无限远处的周围地层温度，T_inf＝T_sur+az；

f(t)——瞬态导热函数，其表达式如下：

ρ_wel——井筒的密度；

c_wel——井筒的比热容；

t_D——无量纲时间，其表达式为：

ρ_for——周围地层的密度；

c_for——周围地层的比热容；

注入井井筒中工质的流动边界条件及热边界条件为：在注入井井筒的顶部，给定注入压力及注入温度；在注入井井筒的底部，给定工质的速度，即步骤1计算的热储层与注入井井筒的交界处的工质速度；

基于方程(5)、(6)、(7)、(8)、(9)、(10)、(11)和(12)及边界条件，采用迭代法计算n时刻注入井井筒中工质的速度场、压力场和温度场。

5.根据权利要求1所述的一种考虑工质储存过程的增强型地热系统的数值模拟方法，其特征在于：步骤3假定在热储层中，岩石和工质之间达到瞬态热平衡，则热储层的能量方程为：

ρ_res——热储层的密度，

c_res——热储层的比热容，其中c_ro为岩石的比热容；

λ_res——热储层的热导率，其中λ_f和λ_ro分别为工质和岩石的热导率；

热储层的热边界条件为：在热储层与注入井井筒的交界处，给定边界温度，即步骤2计算的注入井井筒底部的工质温度；在其余边界，采用绝热边界条件；

将步骤1计算的n时刻热储层中工质的速度场作为已知条件，基于方程(13)及热边界条件计算n时刻热储层的温度场，并根据n时刻热储层的温度场和热储层中工质的速度场估算(n+1)时刻热储层的温度场，其中，温度场的计算和估算均采用隐式迭代法。

6.根据权利要求1所述的一种考虑工质储存过程的增强型地热系统的数值模拟方法，其特征在于：步骤4中产出井井筒中工质的动量方程如下：

产出井井筒中工质的流动边界条件及热边界条件为：在产出井井筒的顶部，给定产出压力；在产出井井筒的底部，给定工质的速度，即步骤1计算的热储层与产出井井筒的交界处的工质速度；和工质的温度，即步骤3计算的热储层与产出井井筒的交界处的热储层温度；

基于方程(6)、(7)、(8)、(9)、(10)、(11)、(12)和(14)及边界条件，采用迭代法计算n时刻产出井井筒中工质的速度场、压力场和温度场。

7.根据权利要求1所述的一种考虑工质储存过程的增强型地热系统的数值模拟方法，其特征在于：步骤5中热储层的平均总应力方程如下：

ν——岩石的泊松比；

α——Biot系数；

β——岩石的线膨胀系数；

E——岩石的杨氏模量；

热储层的应力边界条件为：在所有的边界上，给定平均总应力；

将步骤1计算的n时刻热储层中工质的压力场及步骤3计算的n时刻热储层的温度场作为已知条件，并基于方程(15)及应力边界条件，采用迭代法计算n时刻的平均总应力场；

将步骤1计算的n时刻热储层中工质的压力场及步骤3估算的(n+1)时刻热储层的温度场作为已知条件，并基于方程(15)及应力边界条件，采用迭代法估算(n+1)时刻的平均总应力场。

8.根据权利要求1所述的一种考虑工质储存过程的增强型地热系统的数值模拟方法，其特征在于：步骤6中热储层的孔隙率关于平均总应力的关系式如下：

——初始时刻热储层的孔隙率；

C——热储层的孔隙压缩系数；

σ_m0——初始时刻热储层的平均总应力；

p_f0——初始时刻工质的压力；

热储层的渗透率关于孔隙率的关系式如下：

k₀——初始时刻热储层的渗透率；

基于步骤5计算的n时刻的平均总应力场，由方程(16)和(17)计算n时刻热储层的孔隙率和渗透率分布；

基于步骤5估算的(n+1)时刻的平均总应力场，由方程(16)估算(n+1)时刻热储层的孔隙率分布。

9.根据权利要求1所述的一种考虑工质储存过程的增强型地热系统的数值模拟方法，其特征在于：步骤7中判断n时刻各物理场在经过i次迭代后是否收敛：如果各物理场还未收敛，则将(i+1)赋值给i，进入下一次迭代，并返回到步骤1；如果各物理场已经收敛，则计算n时刻的热开采速率及工质储存量，然后将(n+1)赋值给n，进入下一个时刻的计算，并令i＝1，然后返回到步骤1，直至完成所有时刻的计算；

其中，热开采速率及工质储存量的计算式分别为：

Q_pro＝q_pro(h_pro-h_inj) (18)

Q_pro——热开采速率；

q_pro——产出流量，q_pro＝u_proA_welρ_f，其中u_pro为产出井井口的工质速度；

h_pro——产出井井口的工质比焓；

h_inj——注入井井口的工质比焓；

M_s——工质储存量；

q_inj——注入流量，q_inj＝u_injA_welρ_f，其中u_inj为注入井井口的工质速度。

10.一种考虑工质储存过程的增强型地热系统数值模拟方法的应用，其特征在于：该方法适用于以水为工质的EGS，也适用于以CO₂为工质的EGS。