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)时刻热储层的孔隙率分布；其中，热储层的孔隙率关于平均总应力的关系式如下：——初始时刻热储层的孔隙率；

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

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

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

热储层的渗透率关于孔隙率的关系式如下：k0——初始时刻热储层的渗透率；

基于步骤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，直至完成所有时刻的计算；

增强型地热系统，其初始条件为：

Tres——热储层的温度；

Tsur——地表温度；

a——地温梯度；

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

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

ρro——岩石的密度；

g——重力加速度；

pf——工质的压力；

ρw——水的密度；

θ——积分变量；

t——时间；

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

步骤1所述热储层中工质的连续性方程如下：——热储层的孔隙率；

ρf——工质的密度；

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

——微分算子；

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

*n

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

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

n‑1

ρf ——(n‑1)时刻的工质密度；

n

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

n‑1

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

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

k——热储层的渗透率；

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

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

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

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

步骤2中注入井井筒中工质的动量方程如下：u——工质的速度；

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

2

τ＝fρfu/(4Rwel)    (6)Rwel——井筒的外半径；

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

f ＝‑1.8log10{[K/(7.4Rwel)] +6.9/Re}    (7)K——井筒的等效绝对粗糙度；

Re——雷诺数，Re＝2Rwelρfu/μf；

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

2

Awel——井筒的流动面积，Awel＝π(Rwel) ；

cf——工质的比热容；

Tf——工质的温度；

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

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

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

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

f(t)——瞬态导热函数，其表达式如下：ρwel——井筒的密度；

cwel——井筒的比热容；

tD——无量纲时间，其表达式为：ρfor——周围地层的密度；

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

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

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

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

cres——热储层的比热容， 其中cro为岩石的比热容；

λ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)时刻的平均总应力场；

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

其中，热开采速率及工质储存量的计算式分别为：Qpro＝qpro(hpro‑hinj)    (18)Qpro——热开采速率；

qpro——产出流量，qpro＝uproAwelρf，其中upro为产出井井口的工质速度；

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

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

Ms——工质储存量；

qinj——注入流量，qinj＝uinjAwelρf，其中uinj为注入井井口的工质速度。