张太平,刘 洋,丁若伟,凡倩莹

(1.山东省地质科学研究院，山东 济南 250013；2.自然资源部金矿成矿过程与资源利用重点实验室，山东 济南 250013；3.山东省金属矿产成矿地质过程与资源利用重点实验室，山东 济南 250013；4.中国地质大学(武汉)环境学院，湖北 武汉 430078)

地热能是一种清洁、可持续和可再生的能源，由于其具有储量巨大、产量稳定且基本不受季节、气候影响等优势，其开发利用与存储问题越来越受到人们的广泛关注。实际上，更多的热能贮存于地下的岩石中，也被称为“干热岩型地热资源”。在干热岩型地热资源的开发利用过程中，裂隙-基岩渗流传热问题的基础研究是认识和揭示相关地热能开发利用的热学理论基础，也是解决能源与环境领域诸多问题的关键环节之一，具有较大的工程应用价值。注抽试验广泛用于地热能的开发利用，即先将冷水注入干热岩中获取热量，然后将热水抽出来使用，有时完整的干热岩的储水能力有限，必须采用水压致裂方法提高地热能的开采效率。通常注抽试验分为多井注抽试验和单井注抽(Single-Well Push-Pull,SWPP)试验，其中多井注抽试验至少需要两口井，开展多井注抽试验时，首先将冷水通过一口或者多口井注入到裂隙中，然后通过观测井将热水抽出；而单井注抽试验只需要一口井，其操作是首先将冷水通过一口井注入裂隙介质中，经过一段时间后再将热水从该井中抽出。相比之下，单井注抽试验具有成本低，耗时短和操作简单等优点。

到目前为止，国内外学者已经建立了大量的模型用于研究热在裂隙-基岩含水层介质中物质和能量的迁移转化机理，如Neretnieks、Tang等、Mahmoudzadeh等、Zhu等、Zhou等、Zhu等、Trinchero等。在裂隙-基岩含水层系统中的SWPP热质运移解析解方面，Kocabas采用双重拉普拉斯变换法构建了两阶段的SWPP热质运移试验模型；Jung等采用拉普拉斯变换法和傅立叶变换法构建了三阶段的SWPP热质运移试验解析解;Wang等考虑混合效应物理过程，采用拉普拉斯变换法和格林函数法构建了三阶段的SWPP热质运移试验解析解。混合效应是指注入的冷水与井筒中的热水混合的过程。现有的SWPP试验模型包含很多的假定条件，实际很难满足，例如现有的模型忽略了裂隙-基岩含水层介质中的热弥散和横向热扩散的影响，导致现有的模型难以推广。在数值解方面，目前有很多商业的数值模拟软件(FEFLOW，GMS，COMSOL，TOUGH2等)可以构建SWPP热质运移试验模型，如魏永霞等采用GMS构建了韶关温泉三维地热数值模型，开展了地热资源量评价；刘东林等采用TOUGH2软件开展了东丽湖地区雾迷山组地热资源可开采潜力评价研究；单丹丹等采用COMSOL软件研究了地热深井生产过程中井筒热能损失。然而，Wang等指出现有的商业模拟软件在刻画混合效应过程时假定井筒的水位等于含水层的厚度，因此现有的商业模拟软件难以准确地开展地热资源量评价、地热相关参数的敏感性分析和参数反演方面的研究。

为了更好地理解裂隙-基岩含水层系统中的热质运移机理，本文采用对流-扩散方程来构建裂隙-基岩含水层系统中SWPP热质运移试验的数学模型，但由于数学模型考虑了裂隙和基岩中的横向和纵向热传导和热扩散以及井筒中的混合效应等多个过程，很难导出该数学模型的解析解，因此本文采用有限差分法构建SWPP热质运移试验的数值模型来研究这些过程对试验结果的影响，数值模型采用MATLAB软件中的ODE15s求解器求解。

1 SWPP热质运移试验数学模型的建立

SWPP试验概念模型示意图如图1所示，该概念模型已经广泛用于研究裂隙-基岩含水层系统中地下水的流动、溶质运移和热质运移问题。类似于Wang等的研究，假定井筒完整穿透裂隙介质，裂隙宽度是均一的，基岩中只考虑热传导作用，裂隙中的渗流服从达西定律，本研究采用直角坐标系建立裂隙-基岩含水层中考虑混合效应的SWPP热质运移试验模型，坐标原点位于井壁处的基岩上，以裂隙方向为

x

方向，垂直裂隙方向为

z

方向。SWPP热质运移试验模型的控制方程中包括热对流、热扩散和热传导过程，SWPP热质运移试验模型的控制方程为

图1 SWPP试验概念模型示意图Fig.1 Schematic diagram of the SWPP test conceptual model

0≤

z

≤

b

,

t

>0)

(1a)

(1b)

(2a)

(2b)

(2c)

式中：

k

和

k

分别为裂隙和基岩的热传导率[W/(L·℃)];

β

和

β

分别为裂隙的纵向和横向热弥散度。

在许多前人的研究中，裂隙和基岩中的初始温度分布都假定为常数。实际上，在很多含水层中，温度的分布在空间上并不是均匀的。类似于Chen等对于初始温度的处理方式，本研究假设裂缝和基岩中的初始温度分布是一个任意函数：

T

(

x

,

z

,

t

)|=0=

T

0(

x

,

z

)

(3a)

T

(

x

,

z

,

t

)|=0=

T

0(

x

,

z

)

(3b)

式中：

T

0(

x

,

z

)和

T

0(

x

,

z

)为分布在裂隙和基岩上的温度的任意函数，在交界面上

T

0(

x

,

z

=

b

)=

T

0(

x

,

z

=

b

)。

与Wang等的工作类似，本研究的数学模型的内边界条件考虑混合效应，在SWPP热质运移试验的注入阶段，考虑混合效应的内边界条件为

(4a)

T

(

x

,

t

)|=0=

T

0(

r

)

(4b)

在SWPP热质运移试验的抽取阶段，考虑混合效应的内边界条件为

(5)

其外边界条件为

(6)

裂隙-基岩交界面的温度和温度通量都是连续的，因此交界面的边界条件为

T

(

x

,

z

,

t

)|==

T

(

x

,

z

,

t

)|=

(7a)

(7b)

式中：

φ

为基岩的孔隙度。

2 数值模型的构建

本文采用有限差分法构建SWPP热质运移试验数值模型。为了减小数值计算的误差，本研究采用对数网格系统来对时间和空间进行离散。首先将裂隙空间[

x

,

x

]离散为

N

-1个区域，即

N

个节点，这里采用

x

表示

x

方向上外边界的长度，则

x

方向上的空间离散为

(8)

其中，

x

+12的计算公式如下：

(9)

裂隙介质区域[-

b

,0]和基岩区域[0,

Z

]在

z

方向上分别离散成

M

-1(

M

个节点)个子区域和

M

-1(

M

个节点)个子区域，这里采用

Z

表示

z

方向上外边界的长度,则

z

方向上的空间离散为

(10a)

(10b)

其中，

z

+12和

z

+12的计算公式分别如下：

(11a)

(11b)

因此，本研究构建的SWPP热质运移试验数学模型的空间离散为

(

i

=0,1,…,

N

;

j

=0,1,…,

M

)

(12a)

(12b)

(12c)

(

t

>

t

)

(12d)

(12e)

T

,(

x

,

z

,

t

)|==

T

(

x

,

z

,

t

)|=

(12f)

(12g)

上述方程(12a)～(12g)即为考虑混合效应、横向热传导和热弥散作用的二维SWPP热质运移试验数值模型，数值模型采用MATLAB软件中的ODE15s求解器求解。

3 结果与讨论

3.1 模型的校正

为了检验上述构建的二维SWPP热质运移数值模型和计算结果的精度，本研究将构建的数值模型在特定条件下与前人的模型进行了对比，当

D

=

D

=0和

β

=0以及

t

=0时，本研究构建的数值模型将还原为Wang等的模型。本研究所选取的参数源于Wang等和Jung等研究中，参数设置为：

k

=2.1 W/(m·℃),

k

=2.1 W/(m·℃),

ρ

=2 650 kg/m,

c

=1 000 J/(kg·℃),

ρ

=1 000 kg/m,

c

=4 200 J/(kg·℃),

φ

=0.5,

φ

=0.5,

b

=0.1 m,

B

=5 m,

Q

=-

Q

=1.44 m/h,

r

=0.5 m,

t

=250 min,

t

=400 min,

h

,inj=50 m,

h

,res=47.5 m,

h

,ext=45 m,

D

=

D

=0,

β

=

β

=0，

T

=1℃，初始条件设置为0。图2为本研究构建的数值解在特定条件下与Wang等的数值解的对比。

图2 两模型在井壁上温度穿透曲线的对比Fig.2 Comparison of BTCs at the wellbore computed by two models

由图2可见，修改后两个模型在井壁上的温度穿透曲线能够很好地拟合，从而验证了本研究构建的数值模型的正确性。

3.2 横向热扩散作用对试验结果的影响

相比前人的工作，本研究考虑了裂隙介质中的横向热扩散作用，为了研究横向热扩散作用对SWPP热质运移试验结果的影响，采用数值模型计算了不同有效热扩散系数

D

值下井壁上的温度穿透曲线，其计算结果见图3。其中，

D

=0表示不考虑横向热扩散作用，其他参数值的设置同图2。

图3 不同有效热扩散系数Dfz值下井壁上温度穿透 曲线的对比Fig.3 Comparison of BTCs at the wellbore for the different values of Dfz

由图3可见，随着

D

值的增加，注入阶段井壁上的温度降低，但是在抽取阶段井壁上的温度升高，拖尾现象严重。以上分析表明：裂隙介质中的横向热扩散作用会导致温度穿透曲线发生明显的变化，说明横向热扩散作用不容忽视。

3.3 裂隙宽度对试验结果的影响

横向热扩散作用可能会受到裂隙介质宽度(

b

)的影响，通常来说，当裂隙介质的宽度较小时，横向热扩散作用较小。因此，为了研究裂隙宽度对SWPP热质运移试验结果的影响，采用本研究建立的数值模型，计算了

t

=240 min时不同裂隙宽度下在

x

=0.5 m位置处裂隙-基岩中的温度分布曲线，其计算结果见图4。其中，

D

=0.063 1 m/d，其他参数设置同图2。

图4 不同裂隙宽度b下在x=0.5 m处裂隙-基岩中的 温度分布曲线Fig.4 Distribution curves of temperature at x=0.5 m in fracture-matrix for the different values of fracture width (b)

由图4可见，在横向扩散系数一定的条件下，裂隙宽度为1.00 m裂隙和基岩介质中的温度值最大。但值得注意的是，当

z

<

b

时，裂隙介质中的温度在横向热扩散作用下缓慢降低。以上分析表明：裂隙介质的宽度会对SWPP热质运移试验造成明显的影响，裂隙介质的宽度越大，裂隙和基岩介质中温度值越高。

3.4 裂隙-基岩交界面上的热交换量

在裂隙-基岩交界面上的热交换量能够直接反映裂隙和基岩介质的性质对热质运移的影响，因此本研究采用建立的数值模型，计算了不同时间下裂隙-基岩交界面在

x

方向上的热交换量，其计算结果见图5。其中，

t

=0.5 d和

t

=1.0 d表示注入阶段裂隙-基岩交界面上热交换量的分布，

t

=20 d和

t

=25 d表示抽取阶段裂隙-基岩交界面上热交换量的分布。参数设置如下：

k

=2.1 W/(m·℃),

k

=2.1 W/(m·℃),

ρ

=2 650 kg/m,

c

=1 800 J/(kg·℃)，

c

=1 000 J/(kg·℃),

ρ

=1 000 kg/m,

c

=4 200 J/(kg·℃),

φ

=0.5,

b

=0.1 m,

B

=5 m,

Q

=-

Q

=1.44 m/h,

r

=0.5 m,

t

=1 d,

t

=35 d,

h

,inj=50 m,

h

,res=47.5 m，

h

,ext=45 m。

图5 不同时间下裂隙-基岩交界面在x方向上的热交换量Fig.5 Distribution of the heat flux across the interface of fracture-matrix along the x-direction at different times

由图5可见：在注入阶段，

t

=0.5 d和

t

=1.0 d时裂隙-基岩交界面上的热交换量为负值，这说明在注入阶段热量由裂隙介质往基岩介质中运移，但值得注意的是，

t

=0.5 d时裂隙-基岩交界面上的热交换量大于

t

=1.0 d时的热交换量，这主要是因为在注入阶段的初期，裂隙-基岩中的温度梯度较大，随着注入时间的增加，裂隙-基岩中的温度梯度逐渐减小，进而造成裂隙-基岩交界面上的热交换量减小；相反地，在抽取阶段，裂隙-基岩交界面上的热交换量为正值，这说明在抽取阶段热量由基岩介质往裂隙介质中运移，随着抽取时间的增加，裂隙-基岩交界面上的热交换量增加，这主要是因为随着抽取时间的增加，裂隙中的热量逐渐被抽出，造成裂隙-基岩中的温度梯度逐渐增加，从而造成基岩介质中的热量往裂隙介质中运移。

4 结 论

本文采用考虑横向热扩散作用下的对流扩散方程建立了SWPP热质运移试验数学模型，数学模型的内边界条件考虑混合效应，数学模型的求解采用有限差分法构建数值模型并用MATLAB软件中的ODE15s求解器求解，通过与现有的解析解模型进行对比，验证了数值模型的稳定性和精度，并分析了横向热扩散作用对SWPP热质运移试验结果的影响，得到以下结论：

(1) 在SWPP热质运移试验中，横向热扩散作用对热质运移的影响不能忽视，裂隙介质中的横向热扩散作用会导致注入阶段井壁上的温度穿透曲线值降低，而抽取阶段井壁上的温度穿透曲线值则升高，拖尾现象明显。

(2) 裂隙介质的宽度对SWPP热质运移试验结果的影响明显，裂隙宽度越大，裂隙和基岩介质中温度值越高。

(3) 在SWPP热质运移试验的注入阶段，随着注入时间的增加，裂隙介质中的热量往基岩介质中的运移量逐渐降低，在抽取阶段，随着抽取时间的增加，基岩介质中的热量往裂隙介质中的运移量逐渐增加。