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

权利要求书

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——初始时刻热储层的平均总应力；

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

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

T_res——热储层的温度；

T_sur——地表温度；

a——地温梯度；

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

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

ρ_ro——岩石的密度；

g——重力加速度；

p_f——工质的压力；

ρ_w——水的密度；

θ——积分变量；

t——时间；

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

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

——热储层的孔隙率；

ρ_f——工质的密度；

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

——微分算子；

将方程(2)中的非稳态项作为已知值：若n＝1，则忽略方程中的非稳态项；若n1，则由(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，则将上一时刻注入井井筒底部的工质压力赋值给该边界压力；若i1，则将上一次迭代中步骤2计算的注入井井筒底部的工质压力赋值给该边界压力；在热储层与产出井井筒的交界处，给定边界压力，即若i＝1，则将上一时刻产出井井筒底部的工质压力赋值给该边界压力；若i1，则将上一次迭代中步骤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)时刻的平均总应力场；

步骤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为注入井井口的工质速度。