许 欣，李 宏

(大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024)

地热能作为一种清洁、可再生能源，越来越受到各国政府的重视[1]。传统的增强型地热系统是一种基于石油钻井技术的地热开采方法，技术的关键是通过水力压裂技术在钻孔的底部岩层中产生贯通裂纹。1973年，美国国家实验室和能源部在Fenton Hill开始了干热岩发电技术试验，井深3 000 m，利用花岗岩内的裂缝连接两口井，但正是由于这些裂缝的存在，注入水的回收率大大降低，严重影响水和热的回收和抽取[2]。Gringarten等[3]提出了一种从断裂热干岩石中提取热量的理论。冷水通过每个裂缝的底部进入，热通过水流从干热岩传递到裂缝岩体。这一数学模型可以大大提高干热岩地热系统的经济利用率。

在许多地热田中，有证据表明注入的流体可能会沿着“裂缝”这一优先流动路径快速迁移，且也有大量现场证据表明，在裂缝性储层中，注入井的水可以快速地迁移到生产井[4-7]。Shaw建立了多井系统中干热岩水流吸热的数学模型，研究了多井系统的井距和井径尺寸，对于地热储层的热采效率的影响[8]。不仅如此生产井和注入井的位置同样会对人工地热储层的采热程度以及效率产生影响。

以上研究可以看出，地热储层的开发具有“埋藏深、工程难、投资大”的难点，由于基于钻井技术的增强型地热模型(EGS-D)难以建立大规模稳定蓄热，水流量小，易造成污染等缺点，Tang等[9]提出了一种基于开挖技术的增强型地热模型(EGS-E)开发干热岩。它可以在一定程度上克服水力压裂技术的弊端，避免大高程抽水二次耗能，结合深部矿产资源开采开发深地热，提高综合经济效率。本文以EGS-E中钻爆致裂区为研究对象，连接并顺序执行两个计算机程序TOUGH2和FLAC3D，然后用一维固结问题的解析解验证耦合程序的正确性，最后将程序应用于厚度、孔隙度和渗透率接近均匀的多孔介质理想储层模型，模拟人工地热储层注采井中流体的快速运移及注水井和产能井的位置变化对模型"热扫掠"范围的影响。

1 TOUGH2和FLAC3D 力学耦合水热模型程序流程

1.1 水热过程数学模型

TOUGH2的控制方程除了考虑岩体变形过程中孔隙度和渗透率的变化，还考虑了孔隙压力和流体温度对变形的影响。流动和热对流的一般控制方程是[10]：

(1)

其中，Mκ是模型单元中每单位体积的质量或能量；Fκ是质量或者能量通量；qκ是控制单元中源/汇项。

对于多相流体的流动有：

Mw=φ(SlρlXwl+SgρgXwg)

(2)

Fw=Xwlρlul+Xwgρgug

(3)

其中，Mw是水的质量；φ是孔隙率；S是饱和度；ρ是密度；X是质量分数；u是达西渗流中的速度；下标w是水；l代表液相； g代表气相。

对于热对流的传导过程：

Mh=φ(SlρlUl+SgρgUg)+(1-φ)ρsUs

(4)

(5)

式中：Mh是总热量；U是内部能量；λ是导热系数；T是温度。下标h是热焓；s代表固相。

1.2 力学过程数学模型

FLAC3D是由Itasca公司开发的商业软件，可以对土壤、岩石和其他材料进行热流固耦合分析。 FLAC3D的力学模块中的控制方程是[11]：

(6)

其中，ρm是岩体的平均密度；v是指定系统的速度。在一个时间步长内的应力应变关系可写为：

(7)

FLAC3D以节点为计算对象，力和质量全都集中在节点上，然后通过运动方程使用时域有限差分法[12]和混合离散法[13]求解。节点运动方程可表示如下：

(8)

(9)

在FLAC3D中，为了确保精度和稳定性，通常在运动方程和应力应变本构方程之间设置足够小的迭代步。在每个时间步长内，首先根据初始速度和力进行运动方程的计算，并且将计算的节点速度代入到方程(9)中以获得新的应力和应变速率。

1.3 TOUGH2和FLAC3D的耦合

TOUGH2能够实现流体的流动和热量的传输，而不能反映变形与孔隙度和渗透率之间的关系，所以本文建立连接的TOUGH-FLAC3D模拟器来进行热流固耦合的计算。当连接两个代码时，耦合方程不能同时被求解，而是顺序求解。本研究通过连接外部耦合模块，实现两个代码的顺序执行，并在实验室估计的经验函数上，对各物理场的参数校准，采用数学插值方法实现相互转化。耦合函数通过连接两个代码，实现了最重要的耦合，包括岩石变形对孔隙度和渗透率的影响，以及孔隙压力和温度对岩石变形的影响。

在TOUGH2中，网格坐标由单元的中心点定义。元素大小由相邻单元格的中心点到连接平面的距离确定，而FLAC3D中的网格坐标由角点坐标来定义。由于TOUGH2和FLAC3D分别是以单元中心点和节点为计算单位的，为了实现TOUGH2和FLAC3D的耦合计算，需要通过插值对主要变量进行由单元到单元节点转换。耦合流程如图1所示，可归为以下四步：(1)在TOUGH2中进行多相渗流计算；(2)将单元中心点的温度和孔隙压力通过插值的方法传至节点；(3)在FLAC3D中设置集合收敛和偏差，进行应力应变分析；(4)用得到的平均应力修正孔隙度和渗透率。然后将得到的新的孔隙度和新的渗透率放在TOUGH2的输入文件中进行多相渗流计算，循环直至到达指定时间。

2 一维固结沉降模型解析解和数值解对比

为了验证 TOUGH2-FLAC3D程序的准确性，本文选取存在解析解的Terzaghi一维固结沉降问题进行验证。该模型描述了地面施加载荷后土体排水固结过程(见图2)，究其本质为渗流与力的耦合。固结模型的具体参数如表1所示。

图1 一维固结沉降的弹簧活塞模型[14]

表1 一维固结沉降参数设置

2.1 太沙基单向固结微分方程及其解析解

太沙基单向固结微分方程可表示为如下形式:

(10)

其中,Cv称为土的固结系数，cm2/s，其值为：

(11)

上述固结微分方程表示了超静孔隙水压力u与位置z及时间t的关系，然后根据土层渗流固结的起始条件与边界条件，求出其解析解，当附加应力P沿土层均匀分布时孔隙水压力u(z,t)的解答如下：

(12)

其中，m为奇正整数(1，3，5，……)；

(13)

其中，Es=κ(1+v)/[3(1-v)]为压缩模量。

2.2 结果分析

图2(a)为在t=25 s、409 s、1 638 s、6 553 s和26 214 s时，孔隙压力随深度变化的分布。

在外荷载P的作用下，随着固结时间的增大，土中孔隙水逐渐排出，超静孔隙水压力逐步消散，土骨架的有效应力逐步增大，直至超静孔隙水压力为0，水完全排出。在固结过程中，随着孔隙水的排出，土体产生压缩，使得土体更加密实，强度也逐渐增大。

从图2可以看出，随着时间的推移，孔隙水被排出，土体产生压缩；顶部单元的竖向位移随着孔隙水的排出逐渐增大，在孔隙水全部排出之后，竖向位移逐渐趋于稳定。图2(a)通过对数值解与解析解进行对比，发现整体趋势基本一致，误差较小，这个误差可能由模型网格稀疏以及渗流和力学计算过程中数据插值误差引起的；图2(b)为模型顶部单元的位移与太沙基一维固结问题的解析解对比分析，可以看出，TOUGH2和FLAC3D对一维固结问题模拟的数值模拟结果，与其解析解比较基本吻合，印证了上述力学水热耦合程序的准确性。

图2 一维固结沉降解析解和TOUGH2-FLAC3D数值解对比结果

3 人工地热储层的模拟

3.1 模型建立

以EGS-D难以建立大规模稳定蓄热，水流量小，易造成污染等缺点，唐春安等人提出了一种新的EGS-E开发深部干热岩储层：该方法主要由竖井、上下水管道、换热池、冷- 热水流通道和爆破致裂区组成(如图3所示)。其工作原理：通过开凿竖井、换热池、横巷，形成地下空间。在横巷中实施凿岩爆破致裂，形成水-岩换热接触面积。地下空间施工完成后，一次性注水，并通过爆破致裂区换热，形成高温高压地下热湖。

图3 基于开挖技术的增强型地热系统(EGS-E)[15]

以1 km3的花岗岩储层地热能开发为例，以崩落区的十分之一为目标储层，此工程结构包含25个分支“热湖”构件，则每个热储的容积为0.04 m3，若硐室崩落长1 m，则正方形截面边长为200 m，模拟该平面问题。假设此200 m2的截面内，绝大部分是崩落的岩体，孔隙度为2%。人工地热储层模型以不可渗透岩石的半无限半空间为界，整个初始温度为300℃。在100℃温度下的水以4 kg/s的恒定速率注入储层的一侧，而另一侧则在特定的井筒压力下进行生产，供热产能。在表2中给出了人工地热储层注采模型问题的参数，并且低温水由左侧点I处注入，在右侧点P处产能，模型的详细信息如图4所示。

图4 人工地热储层注采模型示意图

表2 人工地热储层参数设置

3.2 优先流动路径分析

在经过37个时间步后达到规定的最大时间1.577 88×108s，此时生产速度为4.991 4 kg/s，热焓为872.461 kJ/kg； 生产元件的温度为204.73℃。用软件Tec plot 360对模拟结果进行整理，如图5所示为人工地热储层模型在1.5 a、9 a、10 a、13 a、13.6 a和16 a的温度分布云图。

图5 不同时刻人工地热储层温度分布云图(℃)

从图5中我们可以看出，在模型左侧用100 ℃的水以7.5 kg/s的速率注入后，且在右侧以7.5 kg/s的速率抽取。在1.5 a时，在注入井周围由于低温水的注入，岩体与流体之间产生强烈的热交换，此时生产井周围单元的温度急剧下降；在产能井周围，由于液体被抽离，即温度也会有一定程度的下降。随着生产井模拟时间的增大，整个模型的温度逐渐下降。注入模型的流体沿注入井和生产井之间的优先流动路径运移，在换热效率足够高的情况下，岩石与流体之间的换热效果最好，即在这个连线上的温度梯度变化最快。在9 a之后，模型注入井周围的温度都散落在注入流体的温度附近，由于生产井中抽出流体的温度在下降，因此，模型的产热效率逐渐的下降。

由生产的流体温度随时间变化的关系曲线(见图6)，可以看到随着抽取时间的增大，生产井单元的温度在逐渐下降，产生流体的温度变化的斜率也在逐渐减小，最后趋于恒定，说明随着时间的增大，模型中的热量被逐渐抽离，且整个人工地热储层产热的效率在逐渐下降。

图6 生产的流体温度随时间变化的关系

3.3 井位置变化对产能的影响

影响人工地热储层开发的效率除了储层的渗透率、孔隙度、注入井中流体的速率和温度、以及产能井抽取流体的速率和井底压力外，注入井和产能井的位置也会在很大程度上对其产生影响。为此，本小节通过对同一理想的人工地热储层模型，仅改变注入井和生产井的位置来探究其对人工地热储层产热的效能的影响。通过改变注入井和生产井的位置，如图7所示，对称布置的①、②和③处，分别在左侧设置注入井，在右侧设置产能井。

图7 井位置分布

图8显示了注入井和生产井分别布置于①、②和③处时，人工地热储层在13.6 a时的温度云图。从图可以看出，未设置井的两角的区域面积(图中标注)，随着注入井和生产井在相对位置距离的变近在逐渐的增大，即在注入的流体在沿优先流动路径迁移时，扫掠过的面积在逐渐减小；随着注入井和生产井在相对位置距离的变近，各模型中产能井单元温度的变化在逐渐的增大；此外，在同样时间刻度，按对角布置井的位置，储层的热抽取更快。结果表明，按对角布置井的位置可以将储层的热在最大程度上且最快的时间内产出。

图8 在13.6 a不同模型的温度云图(℃)

4 结 论

本文在TOUGH2软件框架的基础上，引入了FLAC3D的应力应变分析，结合TOUGH2共同建立通用的三维力学水热耦合模型，并用解析解验证程序的耦合功能的准确性，得出结论如下：

(1) 建立人工地热储层注采模型，验证了人工地热储层中流体由注入井到生产井这一“优先流动路径”的快速迁移，与现场地热田证据相吻合[4-7]。同时也详细监测了在非等温注入下的生产温度随时间的变化曲线。

(2) 通过改变注入井和生产井的位置，研究了人工地热储层注采模型随井位置的变化所引起采热效率的不同，得出井沿对角这一布置方式，注入流体在优先流动路径下扫掠的面积最大，即储层的生产效率高，开发也较为完全。