干热岩流动换热多尺度有限容积法

宇 波 李庭宇 韩东旭 孙东亮 杨福胜 魏进家

1. 北京石油化工学院机械工程学院 2. 西安交通大学化学工程与技术学院

0 引言

干热岩一般指埋藏于地下3 ～10 km，内部不存在流体（或仅有少量流体）、温度高于150 ℃且具有开发价值的热岩体，是一种储量巨大的清洁能源[1-2]。干热岩热能开采过程与致密油气田的开采有相似之处，均需要通过压裂的方式增强储层的导流性能。但是油气开采主要为“物质”的采出，而干热岩开采为“热能”的采出，二者有本质的区别。数值模拟是研究干热岩热能开发过程的重要手段之一，采用该方法研究干热岩储层中的热能提取规律，对于干热岩开发方案的制订、运行优化以及采热能力评价具有重要意义[3-6]。对干热岩储层进行数值模拟主要有两个难点：①裂缝模型的表征和流动—换热模型的构建；②高效的数值求解方法。

裂缝网络的表征模型在油藏数值模拟中就已经取得了很多进展，近年来相关模型被越来越多地用于干热岩储层的描述中，其中应用最广泛的裂缝模型为离散裂缝模型[7-8]。Sun 等[9]和Shi 等[10]分别使用COMSOL Multiphysics 基于离散裂缝模型分别建立了干热岩热能开采过程的二维、三维热水力三场耦合模型。离散裂缝模型在划分网格时，裂缝要与基岩网格相互对齐，在模拟具有复杂裂缝型储藏时必须使用大量的非结构化网格，严重影响数值计算效率。为了克服离散裂缝模型的缺点，Lee 等[11]发展了嵌入式离散裂缝模型，该模型将裂缝型计算区域分成两种介质——基岩和裂缝，前者可以采用结构化网格划分，而后者可以直接嵌入基岩结构化网格系统中，避免了传统离散裂缝模型复杂的非结构化网格剖分过程。Yan 等[12]基于该模型开发了裂缝岩体的流固耦合模型，计算结果与商业软件相差不大且计算速度较快。Karvounis 和Jenny[13]改进了嵌入式离散裂缝模型，提出采用自适应分层裂缝模型来模拟增强型地热系统的运行过程。

实际储层往往复杂多变，天然裂缝和人工裂缝的长度可以从几毫米到数十米不等。热储层内的物性参数，如渗透率、孔隙度、导热系数等均呈现出不均匀性。因此，采用传统的有限单元法[9-10]、有限容积法[13]或者二者的混合求解法[12]进行求解时，需要划分较细的网格来实现不均匀性的描述，从而导致计算耗时非常大。为此，研究人员发展了一些以多尺度求解方法为代表的模型降阶方法[14]，如多尺度有限单元法[15-16]、多尺度有限容积法[17-18]等。多尺度方法的核心思想是通过局部细尺度（细网格）的求解获得当地基函数，然后基于获得的基函数在粗尺度（粗网格）上对控制方程离散并求解得到粗尺度的解，最后采用局部基函数重构得到细尺度的解。由于不需要直接求解细尺度的离散方程，多尺度方法可以显著提高模拟效率，成为近年来的研究热点。

相比于多尺度有限单元法，多尺度有限容积法不仅可以获得局部守恒的通量且不会引入过多的求解自由度[19]。目前，多尺度有限容积法已经从简单的多孔介质单相流动问题被推广到了多相流[20]、非牛顿流[21]、流固耦合问题[22]等领域。但是，研究人员发现对于不均匀性非常强的问题或网格长宽比较大时，传统的多尺度有限容积法的求解精度明显降低，为此学者提出了迭代多尺度有限容积法[20,23-25]，其通过不断的迭代求解来改进局部问题的边界条件，从而提高多尺度方法的求解精度。Hajibeygi 等[20]提出了类似于几何多重网格的迭代多尺度有限容积法，并将其应用于不可压缩多孔介质单相和多相流动问题，其中细尺度光顺采用线松弛技术；Lunati 等[24]基于多尺度算子提出了一种健壮性较好的多尺度有限容积法，使局部问题的边界条件更为精确。Møyner 等[25]提出了模拟多孔介质流动的多尺度投影限定光滑基方法。该方法预先定义基岩和裂缝粗网格内多尺度基函数的支撑区域，并通过迭代光顺法使得多尺度基函数的求解更加便捷。随后Shah 等[26]将该方法推广到了求解裂缝型多孔介质的流动问题。

受投影限定光滑基模型的启发，笔者提出了模拟干热岩热储流动和换热过程一体化多尺度有限容积求解方法（Thermal-Hydraulic Coupling Multi-scale Finite Volume Method, HT-MsFVM）。与传统多尺度有限容积求解方法不同，该方法在求解裂缝型流动换热问题时主要有两点优势：①在求解温度基函数时无需单独计算各个粗网格上的局部瞬态对流换热问题，而是通过细尺度上的离散矩阵及预先定义的延拓算子（基函数），通过雅克比迭代即可求得局部温度基函数；②通过两步迭代求解策略，可逐步消除细网格上的高频误差和粗网格上的低频误差，对于复杂裂缝流动换热问题，精度更高。HT-MsFVM 方法主要实施过程如下：①基于嵌入式离散裂缝模型，分别采用两套控制方程描述基岩和裂缝介质的流动和换热过程，其中基岩和裂缝介质之间通过质量和能量交换项进行耦合；②采用传统有限容积法将两套方程在细尺度上分别离散为稀疏矩阵形式；③通过定义多尺度网格系统及限定和延拓算子，将细尺度上的离散系统分别在细网格上光顺和粗网格上迭代求解，达到最终所需要的求解精度。

1 基于嵌入式离散裂缝模型的储层流动换热描述

1.1 嵌入式离散裂缝模型

在干热岩开发过程中，储层中的裂缝是流体的主要流动通道和换热场所，因此裂缝的刻画至关重要。嵌入式离散裂缝模型是目前较为流行的裂缝表征方式。该模型有两个关键点：①采用降维的方法描述裂缝网络，在二维问题中，裂缝被看作是一维线段（图1），而在三维问题里裂缝则是一个二维平面；②将基岩和裂缝看成两个独立的介质，根据局部守恒定律构建裂缝、基岩和裂缝与裂缝之间的质量和能量交换。嵌入式离散裂缝模型最大的优点是基岩和裂缝可以通过结构化网格单独划分，避免了局部网格的加密（图1）。

图1 二维嵌入式离散裂缝模型示意图

1.2 流动方程

如图1 所示，嵌入式离散裂缝模型将基岩和裂缝看成是两套独立的介质。基岩与裂缝、裂缝与裂缝之间通过窜流函数进行耦合。基岩系统质量守恒方程为：

式中TI 表示相交裂缝间的交换系数；bi表示裂缝段的开度；Li表示裂缝段的长度。为了提高模拟效率，式（5）的推导过程中忽略了相交裂缝处的网格尺寸。

1.3 换热方程

基于嵌入式离散裂缝模型的思想，仍然采用两套能量方程来分别描述基岩系统和裂缝系统，同时需要考虑二者之间的能量交换，并假设流体和岩石之间处于局部热平衡状态。对于基岩系统有

对于第i 条裂缝介质，采用以下能量方程进行描述：

采用有限容积法可将上述流动和换热方程分别离散为如下矩阵形式：

或矢量形式：

式中A 表示离散的系数矩阵；B 表示方程右端项；x表示待求变量，x=p 或x=T。

式（10）或（11）可以采用传统的迭代法进行求解。但当计算区域较大，离散裂缝较多时，数值模拟较耗时，因此笔者采用多尺度有限容积法进行求解。

2 多尺度有限容积法

2.1 多尺度网格系统

图2 多尺度网格系统图

2.2 多尺度求解方法

根据图2 所示的网格系统，需要定义粗细网格之间的映射关系：

2）粗网格修正：消除粗网格上的低频误差，同时保证粗网格上的质量和能量守恒。

值得注意的是，利用多尺度压力解p 计算得到的速度场v 只能保证在粗网格上的质量守恒。这是由于在计算局部基函数时忽略了粗网格内边界上的通量传递[19]。为了得到细网格上守恒的速度场v'，还需要计算新的压力场p'：

根据新的压力场p'可得到守恒的速度场v'：

利用守恒的速度场v'即可通过方程（17）和（18）计算多尺度温度解。压力和温度的多尺度循环求解过程见图3：首先，根据方程（17）和（18）计算多尺度压力解；其次，根据方程（19）和（20）重构温度方程所需要的速度场；最后，根据方程（17）和（18）计算多尺度温度解。重复上述步骤直至达到所规定的多尺度循环次数（IC）。

图3 压力和温度多尺度循环求解过程伪代码图

3 数值算例

3.1 交叉裂缝算例

本文基于嵌入式离散裂缝模型建立的热开采模型的有效性已经经过实践验证[30-31]。为了验证多尺度有限容积法的精度，设计了图4-a 所示的“十”字交叉裂缝算例。模型边长为100 m，水平裂缝和垂直裂缝的长度均为40 m。模型四周为非渗透绝热边界条件。左下角蓝点为注入井的位置，注入流量为500 m3/d，注入工质的温度为20 ℃。图4 右上角红点为开采井的位置，开采压力为1 MPa，其他物性参数如表1 所示。图4-b 为多尺度计算网格。细网格总数为10 000，裂缝细网格总数为82。基岩粗网格个数为25，裂缝粗网格个数为4。

图4 交叉裂缝算例图

表1 模型物性参数表

为了定量比较多尺度解和细尺度解之间的差距，定义了如下二范数相对偏差：

式中xfs表示细尺度压力解或温度解；xms表示多尺度压力解或温度解。

图5 为压力和温度基函数的空间分布情况。图5-a中①和②及图5-b 中①和②分别为第13 个基岩粗网格和第3 个裂缝粗网格内的压力和温度基函数空间分布。从基岩的基函数可以看出，相交裂缝间的导流作用比导热作用更明显。裂缝支撑域内的压力和温度基函数的空间分布差别不大。图5-a 中③和④及图5-b 中③和④分别为基岩和裂缝粗网格压力和温度基函数各自单独求和后的云图。可以看出，无论是压力还是温度基函数，单独将基岩或裂缝基函数求和并不满足单位分解特性，这是由于根据式（15）迭代计算得到的基函数包含了基岩和裂缝的耦合影响，但将图5-a 中③和④及图5-b 中③和④分别求和后即可使压力或温度基函数满足单位分解。

图5 压力和温度基函数空间分布图

图6-a 为采用传统有限容积法在细尺度上计算开采10 年后得到的压力解。图6-b 为采用本文提出的多尺度有限容积法计算开采10 年后得到的压力场，图6-c为采用传统多尺度有限容积法计算开采10 年后得到的压力场。从图中可以看出图6-a 和图6-b 几乎呈现出了一致的压力分布，根据式（21）算得的相对偏差为2.41×10－3。图6-a与图6-c之间的相对偏差为0.758。图7 为采用两种多尺度方法计算得到的温度解和细尺度温度解间的比较，从图中可以看出无论是基岩区域还是裂缝区域，图7-a 和图7-b 中的两组温度结果均非常的接近，相对偏差为2.30×10－3，而采用传统多尺度有限体积法计算得到的温度解（图7-c）则与细尺度温度解（图7-a）差距较大，相对偏差为0.186。

3.2 复杂裂缝算例

为了进一步验证多尺度有限容积法的计算性能，设计了图8-a 所示的复杂裂缝算例。该算例中，除裂缝条数为101，其余计算参数和前述交叉裂缝算例相同。如图8-b 所示，基岩和裂缝的细尺度网格数分别为10 000 和2 226。基岩和裂缝的粗网格数分别为25 和59。为了表征储层的非均质性，笔者将渗透率、孔隙度以及热导率设定为高斯随机分布（图9）。

图10 为开采20 a 时的细尺度压力解和分别采用两种多尺度方法求得的压力解间的比较。首先从图10-a 中可以看出，由于裂缝的渗透率较大，导流作用明显，致使裂缝附近的压力空间分布与远场基岩的压力分布几乎保持一致，同时可以看出：由本文提出的多尺度方法计算得到的压力解（图10-b）与细尺度压力解（图10-a）吻合较好，而由传统多尺度方法得到的压力解（图10-c）误差较大。图11 为开采20 a 时的采用不同求解方法得到的温度解对比。从图11-a 可以看出，由于储层的非均质性和裂缝的存在，使得温度空间分布呈现出各向异性。从图11-b 可以看出，本文方法得到的多尺度温度解能够非常细致地捕捉到温度变化的各个细节，而传统多尺度计算方法则失效（图11-c）。

不同开采时间下HT-MsFVM 多尺度解和细尺度解间的二范数相对偏差如表2 所示。从表2 中可以看出，压力和温度的相对偏差随着开采时间的增加略有增多，这主要是前期数值计算累积误差的影响，但总体精度完全可以满足实际需要。表3 为在开采10 a 时，不同多尺度循环次数（即图3 中Ic 值不同）下多尺度解和细尺度解间的相对偏差。从表中可以看出，随着循环次数的增加多尺度求解精度增加。这主要是循环次数较少时，细尺度上的高频误差残留较多，从而导致累积误差也会增多。

图6 交叉裂缝算例压力解比较图

多尺度有限容积法所需要的整体计算时间主要消耗在粗网格划分、基函数获取和粗细网格间的迭代求解上。针对特定的干热岩储层，进行大量工况的数值模拟时，网格划分只需进行一次。因此，为了最大程度的利用多尺度方法的计算性能，可以采用离线和在线相结合的方式。多尺度网格系统的划分可以采用离线的方式一次计算完成，而基函数计算和多尺度迭代求解则采用在线的方式根据不同的工况的条件多次求解。图12 为两个算例分别采用传统有限容积法和本文提出的HT-MsFVM 方法的在线耗时比较。可以看出，多尺度方法可以明显提高数值计算效率，加速比可达4 ～5 倍。

图7 交叉裂缝算例温度解比较图

图8 复杂裂缝算例图

图9 物性参数随机分布图

图10 复杂裂缝算例压力解比较图

4 结论

1）针对采用传统有限容积法模拟干热岩储层采热过程效率低的问题，基于嵌入式离散裂隙模型，提出了一体化的流动与换热多尺度有限容积法求解技术（HT-MsFVM）。相比于传统的粗化降阶方法，该方法采用延拓算子和限定算子建立粗细尺度间的交互关系，在提高数值模拟计算效率的同时，保证了计算精度。

图11 复杂裂缝算例温度解比较图

2）多尺度有限容积法的压力场和温度场求解精度会随着模拟时间的增长而略有降低，主要是由于误差累积导致的，但总体得到的结果完全可以满足工程需求。增加多尺度循环次数可以提高计算精度，但循环次数过多会降低多尺度的加速效果。对于本文的算例而言，多尺度循环次数介于50 ～100 时，相对误差的量级介于10－3～10－4，数值计算效率可提高4 ～5 倍。

表2 不同开采时间下HT-MsFVM 多尺度解和细尺度解间的相对偏差表

表3 不同多尺度循环次数下HT-MsFVM 多尺度解和细尺度解间的相对偏差表

图12 多尺度求解方法和细尺度计算方法时间比较图

3）在干热岩的热开采模拟过程中，需要数值计算大量的算例，从而实现热开采参数的优化和不确定性评估。对于通过改变边界条件或物性参数等进行大量数值实验的问题，采用本文提出的多尺度计算方法可有效提高计算效率。

4）为了方便比较，本文暂未考虑物性参数随温度变化带来的物理场空间分布的影响，因此得到的细尺度压力和温度解会与实际情况有一定的偏差。