增强型地热系统关键技术研究现状及发展趋势

巩 亮 韩东旭 陈 峥 汪道兵 焦开拓 张 旭 宇 波

1.中国石油大学（华东）新能源学院 2. 北京石油化工学院机械工程学院 3. 西安交通大学动力工程多相流国家重点实验室

0 引言

地热能因具有清洁性高、运行稳定和空间分布广泛等优点，逐渐成为世界各国重点研究和开发的新型清洁能源。国家能源局《关于促进地热能开发利用的若干意见》[1]明确指出，到2025年全国地热能供暖（制冷）面积比2020年增加50%，地热能发电装机容量比2020年翻一番；到2035年，地热能供暖（制冷）面积及地热能发电装机容量力争比2025年翻一番。地热资源根据地质构造特征、热流传输方式、温度范围以及开发利用方式等因素可分为浅层地热能、水热型地热（地下热水）和干热岩3种类型[2]。其中干热岩是指不含或仅含少量流体，温度高于180℃，其热能在当前技术经济条件下可以利用的岩体[3]。经初步测算，我国地下3～10 km范围内干热岩资源折合标准煤860×1012t[4]，利用其中的2%，相当于我国2020年能源总消耗量的约3 400倍。我国青海共和地区已钻探发现温度高达236 ℃的高品质干热岩[5]，干热岩开发在我国具有广阔的应用前景。

1 增强型地热系统研究简况

增强型地热系统（Enhanced Geothermal System，EGS）是开发干热岩型地热资源的主要手段。EGS形成过程为：钻注入井，经注入井向储层内注入高压水，迫使岩石开裂从而形成具有渗透性的多尺度人工缝网结构（图1-a），之后依据人工缝网走向确定采出井位置，形成一注一采、一注多采或多注多采等开发模式。EGS采热过程为：经注入井将水、CO2等低温工质注入具有人工缝网的干热岩储层，与高温岩体充分换热后，经采出井采出，由地面发电系统将热能转化为电能（图1-b）。据统计，全球在建及投运的EGS工程已达30个，其中14个实现了运行发电，目前尚在运行的有5处，总装机容量为12.2 MW。整体来看EGS尚未实现规模化、商业化运行[6]，究其原因为一系列瓶颈有待突破，主要有：①EGS开发过程中，缺乏有效的干热岩人工造缝调控技术，导致人工裂隙单一且尺寸较大、流体短路循环，造成过早热突破、采热效率低；②EGS形成和采热过程受渗流、传热、介质变形、水岩反应等多种因素的影响，地热储层内多尺度多场耦合发展规律和机理仍不明确；③EGS采出井举升过程中高压降引起的流体闪蒸问题，影响井内的流动换热特性，制约井内热流体的高效提取；④EGS地上发电模式众多，但具体适用场景并不明晰，仍存在热电转换效率低的问题。

图1 增强型地热系统示意图

因此，围绕上述重大需求和瓶颈问题，笔者针对EGS开采干热岩型地热资源涉及的4项关键技术进行了综述，包括干热岩储层人工压裂技术、干热岩开采数值模拟技术、井筒热流体高效提取技术以及干热岩地热发电技术，以期为我国干热岩高效开发相关研究的发展提供一定的指导。

2 干热岩储层人工压裂技术

干热岩储层较为致密，渗透率极低，为保证注入流体能够与储层岩石进行充分热交换并顺利经采出井采出，一般需要通过储层改造形成连通良好、导流能力高的流动通道，提高储层渗透性。但在储层改造过程中，易形成单一的高渗裂缝，造成流体流动短路、过早产生热突破，影响EGS的可持续开发利用。因此，明晰人工缝网形成机制、精准预测压裂裂缝以及有效调控缝网结构，是形成复杂人工缝网和提高EGS采热效率的一大关键。下面将从干热岩储层改造方法、人工缝网形成机制和裂缝扩展预测模型3个部分总结干热岩储层人工压裂技术方面的研究进展和发展趋势。

2.1 干热岩储层改造方法

干热岩储层改造方法主要包括[7]：水力改造、化学改造和热改造，表1对3种方法进行了简要总结。水力改造又称水力压裂，指将低温高压的流体注入地层，在井底憋压产生人工裂缝，从而提高地热储层渗透性的方法[8-10]，是最主要的储层改造方法。注入流体一般为清水，也可以是高黏压裂液，一般垂直于最小应力方向，岩石最易发生开裂和扩展。油气开发中压裂易造成岩石拉伸破坏，而EGS主要为水力剪切破坏，能够依靠裂缝粗糙度使裂缝维持张开状态，裂缝较宽但开度较小，有利于流体与储层的充分换热[11]，但这种自支撑机制并不总是有效，依然需要注入支撑剂，保证裂缝张开[12]。

表1 干热岩储层改造的主要方法表

化学改造，包括酸化和碱性化学刺激等手段，通过向井中注入酸性液体/碱性液体，溶解储层的矿物质[13]，达到改善裂缝连通性和渗透率的目的。酸化过程一般需加缓蚀剂，减缓酸液对完井管柱的腐蚀。酸化过程通常包括前置酸、主体酸和后置酸3个阶段[14]，前置酸一般用HCl溶解碳酸盐类岩石矿物，主体酸用土酸（HCl/HF混合物）来溶解硅酸盐类岩石矿物，后置酸一般用清水将管柱中的酸液顶替进入地层中。室内实验结果表明[15]，12%HCl+5%HF的土酸可显著改善花岗岩的渗透性，在最优反应时间下，渗透率可提升4个数量级。微观结构和酸岩反应离子探测结果表明[15]，土酸与长石和云母矿物反应较快，而与石英反应相对不活跃。因此，石英颗粒可充当压裂支撑剂，使裂缝保持张开。Watanabe等[16]研制了一种新的化学改造方法，主要使用螯合剂来持续溶解特定的矿物，进而形成流动通道，研究结果表明该方法可使得储层渗透率增大近6倍。Xu等[17]对比研究了高温高压条件下酸性液体和碱性液体对干热岩样品的溶蚀情况，实验结果发现，酸蚀后岩石渗透率提高了1.62倍，而与碱性溶性反应后，岩心渗透率可提高2.45倍。环境扫描电镜微结构分析结果表明，干热岩与酸液和碱性液体反应后，岩心内部均产生了二次矿物，包括绿泥石、球形二氧化硅和蒙脱石等。但现阶段化学改造主要用于地层近井附近，溶解岩石或裂缝里的充填物[18-19]，作用范围有限，是进行储层改造的次要手段。

热改造，指在低于岩石破裂压力时将冷流体注入高温地层，随着注入时间延长，地层的注水能力逐步得到提高[20-21]。热改造机理包括以下2个方面：①注入的冷流体会引起岩石冷却收缩，使裂缝张开；②地热储层冷却引起的热应力会促进裂缝扩展[22-23]，缝面剥落碎片可起到支撑裂缝的作用，使裂缝保持永久性张开，进而提升注水能力。武治盛等[24]利用超声波检测技术，实验研究了裂隙充填干热岩体热冲击后纵波波速演化规律，0 ℃水冷却冲击对岩石的损伤比室温冷却冲击更强；对比了不同热冲击实验条件下的波速降幅变化规律，由于胶结面附近的矿物颗粒尺度差异最大，热破裂程度较强，含胶结面充填花岗岩的波速降幅比母岩、充填体大。Baujard等[25]研究了法国Rittershoffen EGS工程中花岗岩的热性质和水力性质，由于初始生产指数较低，通过低排量注入冷水、局部化学注入和水力刺激改造，可极大地改善井间水力连通性，从而提高EGS取热效果。目前，热改造方法的研究与应用较少，一般作为储层改造的辅助手段，与水力压裂和化学改造方法联合使用。

干热岩一般埋藏较深且地应力各向异性较强，仅依靠上述常规储层改造技术一般难以形成复杂的裂缝网络，注水采热时易产生快速热突破现象，导致采热效率较低。近年来，有学者提出暂堵转向压裂技术，实现人工缝网准靶向调控，是未来提高干热岩采热能力的一种新途径[26]。该技术一般使用可降解纤维或变形颗粒作为暂堵材料，通过暂堵增加缝内净压力，从而逼迫裂缝转向形成多条裂缝，以提高人工缝网的形成效率（图2）。初步研究结果表明[26]，高温条件下，干热岩人工裂缝暂堵过程中会出现明显的压力波动，温度越高，压力波动次数越多。同时，高压作用下缝内突破封堵的次数与压力波动存在明显的关系，可以通过提升排量、增加缝宽来提高暂堵剂对裂缝的封堵效果。

图2 暂堵转向压裂示意图

2.2 人工缝网形成机制

合理的人工缝网形成不仅依赖于相关调控技术，还需要明晰干热岩所处环境下的岩石力学特性以及人工缝网形成机制。

干热岩由于高温高压作用，岩石塑性增强，压裂呈现剪切破坏模式，实验研究结果表明[26]：温度增加40 ℃，可使干热岩的断裂韧性减小30%。高温将使得干热岩的峰值强度、弹性模量等岩石力学特征参数出现显著的劣化，从而加剧干热岩内部微裂隙的产生。此外，干热岩的破裂特征一定程度上受热流过程的影响[27]，在干热岩热采过程中，注入冷水可使裂缝面冷却，导致缝宽增加，这将影响人工裂缝扩展，或者诱导二次裂缝的产生（图3）。热冲击形成的裂缝具有自相似特征，且热冲击形成新裂缝面所需的能量与特征长度、地应力等因素有关。实验研究结果表明，热破裂的阈值温度为260 ℃，此温度下发育大量裂隙或形成贯通网络。热应力诱导的裂缝开度是增大近井区域导流能力的主要原因。在高温条件下，干热岩将出现轴向劈裂、双缝共轭剪切、多缝平行剪切、单缝剪切4种破坏模式[28]。Liu等[29]通过热流固耦合三轴压缩实验系统，研究了干热岩在不同破坏模式下水力传导能力的演化特征，干热岩缝面粗糙度和特征应力随着围压增加而增加，复合破坏模式下干热岩渗透率最大（11.4 μm2），“Y”形剪切破坏模式下渗透率次之（9.7 μm2），最小渗透率情况为单纯剪切破坏裂缝（7.2 μm2），增加围压大小可减小渗透率。Liu等[30]利用真三轴水力压裂实验系统，研究了循环注入条件下干热岩的裂缝起裂与扩展机制，水力裂缝的起裂可分为2个阶段：仅2～4条简单裂缝阶段和大于4条的复杂裂缝网络阶段，水力裂缝的起裂与扩展主要受流体循环、循环注入时间和注入速率等因素影响。在低排量注入条件下，高频循环下（循环时间为10 s、20 s）水力裂缝的起裂与扩展主要受注入压力影响。干热岩可压性，即储层形成复杂缝网的能力，Wang等[31]通过高温高压三轴试验压机，以声发射为监测手段，开展不同温度条件下干热岩可压性评价的实验研究，基于应力应变曲线获得不同应力加载水平下的体积应变，据此建立了考虑热效应的脆性指数计算模型，该脆性指数不仅可以描述初始加载阶段的热损伤微裂纹密度，还可以描述不同应力水平下的裂纹起裂、扩展和成核特征。在上述脆性指数模型基础上，发展了考虑热效应的干热岩可压性评价新方法，该方法可适用于单轴和三轴应力状态情形。

图3 注冷水诱导热应力产生的二次裂缝图

此外，干热岩水力压裂过程中需注入成千上万立方米的压裂液以形成复杂的裂缝网络，流体与干热岩间相互作用对干热岩缝网扩展的影响研究较为缺乏。为揭示干热岩在不同饱和流体、不同加载应力水平下的损伤断裂特征，国内学者利用高温高压岩石力学三轴压机，通过监测干热岩的超声波波速和声发射响应特征，开展了蒸馏水、纳流乳液等不同流体饱和干热岩的岩石力学性质实验研究。干热岩处于饱和流体状态后，纵波波速比干燥干热岩升高40%，并且裂缝损伤应力比、断裂韧性和黏聚力等力学参数明显降低。饱和流体干热岩具有更高的声发射率，累积声发射曲线呈现多个台阶状上升特征。实验结果表明[32]：饱和流体可以增加干热岩的孔隙压力，促进干热岩水力压裂过程中复杂缝网的形成，且纳米乳液比蒸馏水形成复杂缝网的能力更强。

近年来，国内学者借鉴物体“热胀冷缩”物理思想，对干热岩实施交变温载处理，增强压裂改造过程中干热岩储层人工缝网复杂度，并研究了交变温载作用下干热岩人工缝网扩展规律[33]。如图4所示，交变温载后形成缝网的分形维数较单纯热处理后明显提升，随温度升高，裂缝分形维数也逐渐增大，在500 ℃时分形维数最大提高了20%。交变温载后，水力压裂过程中干热岩的破裂压力和裂缝扩展时间显著降低，在500～600 ℃时，破裂压力降幅最大为51%，说明交变温载作用能够降低岩心破裂强度，并加速岩心破裂，形成宏观贯通裂纹。

图4 交变温载与单纯加热处理干热岩样品水力压裂裂缝形态对比图[33]

2.3 裂缝扩展预测模型

干热岩中人工缝网的形成机制复杂，要实现对储层缝网的准确预测，需要构建综合考虑热应力、裂缝内流体流动以及岩石应力等因素影响的裂缝扩展预测模型。目前，干热岩水力压裂模型主要包括二维模型、拟三维模型和全三维模型，各种模型进展情况详述如下：

二维模型假定裂缝高度恒定，一般可得到缝长、缝宽、压力等参数的解析解，较为著名的模型有PKN模型、KGD模型和Penny模型[34-39]。PKN模型[9-10]适用于缝长远大于缝高的情形，它假设每个垂向截面是相互独立的，相对于缝长方向的压力来说，缝高方向上的压力起主导作用，而KGD模型[36-37]适用于缝高远大于缝长的情况，模型假设平面应变在水平方向上，即所有水平截面相互独立，裂缝宽度在垂向上比水平方向变化的更加缓慢。实际上，由于岩石具有渗透性，水力压裂过程中流体沿着缝面滤失，因此PKN和KGD模型需要考虑缝内液体滤失和流体存储效应。1957年，Carter提出了流体滤失方程，随后Nordgren将滤失效应和流体存储效应引入经典的PKN和KGD模型中，弥补了原有模型的不足。Penny模型[38-39]也称径向裂缝模型，适用于水力压裂形成水平缝的情形。

由于二维模型假设缝高恒定，难以适用于多层水力压裂缝高垂向扩展的情况，因此拟三维模型应运而生。该模型假设缝长方向的延伸大于缝高方向的延伸，并将缝内流体流动简化为一维或二维流动，水力裂缝可以三维延伸，各个截面上变形为二维平面应变状态且相互独立，因此计算量要比全三维模型低[37-40]。基于平衡高度假设，Simonson等[40]推导了垂向含有3个小层、对称分布的拟三维压裂模型；随后，Fung等[41]得到了垂向含有多个非对称小层的拟三维压裂模型；Zhang等[42]考虑了层状岩石材料属性和地应力的差异性，构建了一种新的拟三维压裂模型，垂向平面裂缝沿着横向划分为一定数目的位移不连续单元，每个单元内横截面变形为平面应变，流体压力可沿着垂向和水平向2个方向发生变化。

二维模型和拟三维水力压裂模型仅能考虑平面裂缝扩展，无法考虑非平面裂缝转向扩展。因此，国内外学者发展了全三维水力压裂模型，该模型可以真实反映水力压裂在三维空间的裂缝扩展情况，裂缝扩展可以同时在3个方向上进行且可发生空间扭曲，流体可以在2个方向流动，因此全三维模型的计算量也最大[43-44]。Settari等[45]发展了一种全三维水力压裂模型，该模型耦合了垂直平面内三维水力裂缝内的两相流、支撑剂运移和热传递过程，此外还可模拟停泵后的裂缝闭合过程，能预测压力的变化趋势。Yan等[46]建立了一种新的全三维水力压裂模型，该模型中裂缝内流体流动通过裂隙渗流来表征，基质中的流体流动通过孔隙渗流来表征，可以模拟任意复杂裂缝网络的扩展。Wang等[47]发展了基于生死弹簧单元的暂堵本构模型，构建了干热岩变形、基岩渗流、裂隙流完全耦合的三维计算模型，该模型可自适应建立复杂人工缝网结构。弹簧元在加入暂堵剂时可以自动激活，停止注入暂堵剂时弹簧元自动删除，可模拟暂堵剂加入过程中缝内净压力的增加过程。

3 干热岩开采数值模拟技术

干热岩储层含有多储渗结构，包括孔隙和裂缝，其中裂缝尺寸从厘米级跨越到米级（甚至千米级），具有多尺度、强非均质特征；另外，低温流体注入高温储层，破坏储层内温度场、压力场、应力场和化学场，产生强烈扰动，热流固化多场耦合效应显著。因此，干热岩开采属于典型的多尺度多场耦合问题。笔者将从岩心尺度多场耦合模型、储层尺度多场耦合模型和尺度升级方法3个部分总结干热岩开采数值模拟技术方面的研究进展和发展趋势。

3.1 岩心尺度多场耦合模型

在岩心尺度下可利用计算机图像处理技术或者分形建模技术对岩石的孔隙结构与矿物组成特征进行精细刻画，进而对岩石内部的流动传热和受力变形特征进行精确描述，方便定量研究各种微观因素对岩石物理特性的影响。此外，岩心尺度数值模拟也可从微观角度揭示宏观现象，弥补传统岩石物理实验的诸多不足，揭示物理实验无法展现的机理[48-49]，常见的岩心尺度多场耦合模型如表2所示。

表2 岩心尺度多场耦合模型表

岩心尺度中的力场求解涉及多孔介质内部应力较小时的连续变形过程，和内部应力较大时由完整材料到起裂以及裂缝扩展的非连续变形过程，大致可分为连续介质力学模型和离散力学模型。连续介质力学模型下，除裂缝以外其他区域中孔隙和固体骨架被假设为连续变形的块体，通过有限元方法（Finite Element Method，FEM）对基岩中连续变形区域的静力学平衡方程进行离散求解。对存在裂缝的非连续变形区域，其周围应变梯度往往较大，通常需结合塑性变形理论和损伤力学理论求解[50-51]。但当岩心内存在复杂裂缝网络时，由于裂缝周围区域具有较强的非连续性，连续介质力学模型的模拟效果往往不佳[52]。相较连续介质力学模型而言，离散力学模型将材料视为离散粒子或块体的集合（图5），采用拉格朗日方法追踪每个粒子的运动，更适于计算非连续性大应变。常见的离散力学模型有离散元法（Discrete Element Method，DEM）[53]和非连续变形分析法（Discontinuous Deformation Analysis，DDA）[54]。离散力学模型中时间步长较小，计算负荷较大，适于岩心尺度下的力学分析，随着计算机性能的提升也可进行一些大尺度的破裂模拟[55-56]。随着模拟技术的发展，连续介质力学模型和离散力学模型之间的界限已不再清晰，连续介质力学模型在非连续变形处与离散力学模型结合能够形成优势互补。例如有限元/离散元耦合方法（Finite-Discrete Element Method，FDEM）[57-58]，即在FEM中的黏结单元失效后，两侧单元由黏结关系转变为接触关系，接触力的计算采用DEM中的方法。

根据孔隙的描述方法，岩心尺度下的流动与热质传递模型可分为以下2种类型：①孔隙重构模型，即在详细刻画岩石孔隙结构的基础上，求解岩石内部裂缝、孔隙中的流动，再根据流体速度场求解孔隙内部的传热和传质方程（图6-a）；②孔隙网络模型[61]（Pore Network Model，PNM），将岩石内部裂缝、孔隙的流动通道简化为球形孔隙和喉道的流通关系，流体存储在孔隙中，孔隙之间流体的质量交换由一维的喉道描述，流体与固体之间的相互作用如对流换热、流固作用力等也需经喉道计算得到（图6-b）。

图6 岩心尺度下裂缝及孔隙的描述方法图

在第一种类型中，可采用传统的CFD方法如有限体积法（Finite Volume Method，FVM）、有限差分方法（Finite Difference Method，FDM）对孔隙空间中的Naiver-Stokes方程、传热和传质的对流扩散方程进行求解，具有较高的计算效率和数值稳定性，也适用于岩心多相流动中大密度比、黏度比的模拟。另一方面，近年来快速发展的格子玻尔兹曼方法（Lattice Boltzmann Method，LBM），由于在网格上近似求解时间、空间、速度离散化后的玻尔兹曼方程而具有粒子动力学内在特征，在处理复杂形状的流固界面时具有显著优势，LBM中多分布函数模型的发展也使其适合多场耦合方面的模拟[62]。在无网格方法中，光滑粒子流体力学（Smoothed Particle Hydrodynamics，SPH）方法近年来也广泛应用于多场耦合模拟[63]，其核心是根据粒子之间的距离计算粒子的未知量场[64]。值得注意的是，LBM、SPH在多场耦合中的数值稳定性较差，容易产生数值振荡。对于孔隙网络模型而言，真实裂缝和孔隙的复杂几何结构被简化，因此具有很高的计算效率，是研究多孔介质中热质传输规律的有效办法。因为孔隙和喉道的尺度大小和空间分布存在多种选择方式，且与多孔介质的结构特征存在密切关系，如何从岩心的CT数字扫描中提取出合理的孔隙网络模型是该类问题研究的难点和热点之一[65]。具体到流动模拟时，一维喉道中的流动通常采用达西或泊肃叶稳态流动进行近似计算，孔隙之间的流动总体满足质量守恒，并且孔隙压力波动由孔隙中流体体积变化引起[66-67]。另外，在计算喉道的传热和传质方程时，为描述流固界面的信息交换，不可避免地需要引入简化后的经验公式。

当力场模型和流动与热质传递模型耦合计算时，上述方法两两之间结合可以产生多种组合模型，如图7所示，即流固（Hydro-Mechanical）[55-56]、热固（Thermal-Mechanical）[52,68]、热流固（Thermal-Hydro-Mechanical）[50,59,69-70]等岩心尺度模型。化学反应会造成裂缝通道、孔隙结构的变化，其对岩石力学特征影响较大，但缺乏有效数值手段表征，因此目前有关岩心尺度下力场与化学场耦合的数值模型如流固化（Hydro-Mechanical-Chemical）、热流固化（Thermal-Hydro-Mechanical-Chemical）模型的相关研究较少。

图7 岩心尺度热—流—固—化多场耦合图

3.2 储层尺度多场耦合模型

储层尺度的模拟区域更大、时间更长、更符合工程需求，对于整个采热系统的取热性能预测和评价具有重要意义。储层尺度模拟主要涉及裂缝型储层中裂缝的描述、多场耦合控制方程与求解方法2个方面。

裂缝型储层的表征方法主要分为3类：等效介质方法、双重孔隙介质方法、离散裂缝类方法（图8[71]）。等效介质方法不直接描述岩石中的裂缝和孔隙，而是将裂缝和孔隙中的孔隙度和渗透率等效为岩石的原生孔隙度和渗透率，之后采用连续介质对应的一系列理论进行建模[72-74]。对于岩石总体近似表现为具有对称渗透率张量的均质多孔介质，该方法可以得到符合实际的结果，然而对于裂缝主导的流动，该方法忽略了岩石基质和裂缝介质间的流量交换，模拟结果往往误差较大。双重孔隙介质方法认为基质原生孔隙与天然裂缝相互连通，同时基质和裂缝即是流体的储存场所又是流动通道，但模型中基质被分割为正方体、圆柱体、球形、层叠等规则形状[75]，不能充分体现裂缝的各向异性和不连续特征。相较上述2种方法而言，离散裂缝类模型以显式表征裂缝为核心，再结合连续介质力学和裂缝渗流力学等理论进行求解，可以提高储层中裂缝的计算精度。根据物理情形不同，离散裂缝类模型又可分为离散裂缝模型（Discrete Fracture Model，DFM）[76-77]、离散裂缝网络模型（Discrete Fracture Network，DFN）[78]、嵌入式离散裂缝模型（Embedded Discrete Fracture Model，EDFM）[79-80]。但离散裂缝模型满足高精度要求的同时，存在网格划分困难、计算量大等问题，难以用于具有复杂裂缝结构的三维实际干热岩多场耦合数值模拟。

图8 裂缝岩体表征示意图（据李庭宇[71]修改）

基于上述介绍的裂缝描述方法，已开发出很多储层尺度下热流[81-82]、流固[83-84]、热流固[85-86]、热流化[87-88]、热流固化[89-90]耦合模型，此处对热流固化四场中的经典控制方程和其之间的部分耦合关系进行介绍。

岩石中力场求解多基于线弹性假设的静力学平衡方程，在孔隙压力、热应力的影响下应力应变本构方程为：

式中，σc表示柯西应力张量，Pa；εc表示应变张量；λ和G表示拉梅系数，Pa；α表示Biot系数；β表示热膨胀因子，Pa/℃；I表示单位张量；p0表示参考压力，Pa；T0表示参考温度，℃。

基岩和裂缝中的质量守恒方程分别如下式所示：

式中，上标m和cr分别表示基岩和裂缝系统；p表示孔隙压力，Pa；ρf表示流体密度，kg/m3；φ表示孔隙度；ct表示压缩系数，Pa-1；u表示达西速度，m/s；Qm-cri表示裂缝i到基岩的网格体积平均后的流体质量交换，1/s；Qcri-m表示基岩到裂缝i的网格体积平均后的流体质量交换，1/s；Qcri-crj表示裂缝j到裂缝i的网格体积平均后的流体质量交换，1/s；QB表示注入井或抽出井造成的源项，1/s。

基岩和裂缝中的压力场和渗流场之间满足达西定律：

式中，μ表示孔隙流体的黏度，(N·s)/m2；K表示渗透率；g表示重力加速度，m2/s。

基岩和裂缝中的能量守恒方程分别如式（5）和式（6）所示，其中基岩内部流体和固体骨架之间被假设为局部热平衡状态。

式中，T表示温度，℃；cp表示比热容，J/(kg·℃)；Γ表示热扩散系数，m2/s；下标eff表示孔隙水和固体骨架综合考虑后的等效参数；E表示能量交换，J/(s·m3)；EB表示注入井或抽出井造成的源项，J/(s·m3)。

基岩和裂缝中的溶质运移方程分别为：

式中，C表示溶质浓度，mol/m3；D表示溶质扩散系数，m2/s；HB表示注入井或抽出井造成的源项，mol/(m3·s)；γ表示溶质浓度交换，mol/(m3·s)；HC表示化学反应产生的源项，mol/(m3·s)，其与水岩反应速率有关，水岩反应速率可以表示为：

式中，kC表示化学反应速率常数，mol/(m2·s)；KC表示反应平衡常数；AS表示反应比表面积，m2/kg；γ表示离子活度积；指数θ和η与化学反应级数有关。

化学反应速率常数采用Arrhenius公式描述：

水岩反应会对基岩和裂缝中的孔隙度产生影响：

式中，φ0表示初始孔隙度；ρs表示固体矿物密度，kg/m3；Ms表示固体矿物摩尔质量，kg/mol。

孔隙度变化造成渗透率变化，采用Kozeny-Carman方程描述基质中孔隙度（φ）和渗透率（K）之间的关系：

式中，K0表示初始渗透率，m2。

孔隙度φ的变化进一步造成反应比表面积AS的变化，它们之间的关系为：

式中，A0表示初始反应比表面积，m2/kg。

上述控制方程的求解方法众多，常见的传统方法有FEM、FVM、FDM等，不同方法之间的原理和所得到的离散方程差别较大。近年来，扩展有限元方法（Extended Finite Element Method， XFEM）和扩展有限体积法（Extended Finite Volume Method，XFVM）[91]被提出，其核心思想是在裂缝所在单元的节点处引入增强函数来描述单元内部裂缝面处的非连续变形，同时这些增强函数不会影响传统FEM中形函数的特性。XFEM相较于传统FEM的主要优点是无需在非连续变形区域加密网格或者采用大量的非结构化网格。XFEM一经提出便被广泛关注且快速发展[92-93]，著名商业软件ABAQUS也植入了XFEM方法，可用于进行非交叉裂缝的模拟。XFEM缺点是增强函数引入后整体离散方程容易接近线性相关从而增加了收敛难度，同时裂缝分布较多时大量裂缝贯穿单元增加了求解自由度使得计算耗时较大。借鉴XFEM的部分思想，Deb等[94-95]提出了基于有限容积法离散非连续变形问题的XFVM方法，其可以将XFEM中裂缝贯穿单元增加的8个自由度降低至4个，从而适于快速求解。XFVM求解力场并结合FVM求解THC控制方程可以使得流动和力学变量处于相同网格节点，易于通用式编程和网格量较大时的并行计算[71]。

3.3 尺度升级方法

由于储层中裂缝和孔隙存在明显的多尺度特征并且不同尺度的裂缝和孔隙在储层模拟中表现出的规律和影响不同，尺度升级方法成为储层描述的新概念和新发展方向。截止目前，针对裂缝型储层尺度升级方法的研究正处于发展阶段且多为HM耦合的水力压裂模型，按照粗、细尺度之间信息交换的方法主要分为计算均匀化方法、连接尺度方法、多尺度有限元法、多尺度有限体积法等。

计算均匀化方法是建立不同尺度联系时最常用的方法，在尺度连接时多基于平均场理论，假设在粗尺度和细尺度下由微元变形造成的虚功是相等的来进行求解[96]。粗尺度下采用FEM方法求解，在有限单元的每个高斯积分点处设置一个可以反映该位置处微观结构下受力、流动等表征属性的细尺度单元，通过对细尺度单元界面处的边界值问题的计算求解获得单元内部响应，之后进一步获得积分点处的应力—应变张量、渗透率等参数用于粗尺度求解。在细尺度单元边界值处理时，前人研究表明当采用均匀变形、零转动以及流动传热中的Dirichlet边界条件时会高估表征属性，而均匀受力、自由转动以及流动传热中的Neumann边界条件会低估表征属性，采用周期性边界条件可以得到处于两者中间较合理的结果[97-99]。计算均匀化方法对于裂缝、岩性分布均匀的储层，能很好地融合宏观和微观的特征，摒弃宏观唯象本构假定，根据微观结构特性估计其宏观的等效特性[100-101]。此外，细尺度单元可以采用多种方法离散求解，常见的有FEM方法、DEM方法[102]等（图9）。

图9 采用计算均匀化方法的多尺度建模技术图

连接尺度方法最早是由Wagner等[103-104]在将分子动力学模型和宏观连续体模型连接计算时所提出，其将分子动力学的计算区域设定在局部小部分区域，而在其他绝大部分区域中依然采用宏观连续体模型计算。该方法将节点位移未知量分解为相互解耦且正交的粗尺度和细尺度两个部分，粗尺度的解由FEM形函数插值表示，而细尺度的解被表示为偏离插值的一部分，由离散力学模型如DEM求解[105]。连接尺度之间的投影算子是该方法的关键，其准则是使总动能的粗细尺度部分解耦，即总动能或总动力虚功可分解为粗、细两种尺度的动能或动力虚功之和，进而解耦求解表示颗粒材料的运动行为。万柯等[106-107]提出了一种基于DEM-FEM的连接尺度法模型，其在多尺度交界面处提出了高效的准静态界面条件，并基于此模型计算了边坡稳定问题[108]，其中在边坡上方大变形的梯形区域内采用了DEM粒子计算（图10）。此外，连接尺度方法也已成功应用于材料动态断裂[109]、流体流动[110]的多尺度模拟中。在储层人工压裂和干热岩开采过程中，可在裂缝周围采用DEM求解来模拟非连续变形和裂纹扩展，即进行细尺度表征，在热响应较小的连续变形区域采用FEM方法求解，即进行粗尺度表征，在两者界面处采用连接尺度方法，进而实现高效准确的储层模拟。

图10 采用连接尺度方法的多尺度建模技术图（据Wan等[108]修改）

多尺度有限元法、多尺度有限体积法的核心思想在于通过对单元内细尺度局部子问题的数值模拟，得到粗尺度下当地单元内的基函数，进一步在粗尺度下基于获得的基函数对控制方程离散求解[111]。这些基于细尺度获得的基函数不同于传统方法中采用多项式表示的基函数，其可以准确描述材料粗单元内部的非均质性，并可以自动地将细尺度下解的信息代入到粗尺度范围内，得到粗尺度下较为准确的结果。在多尺度有限元法中，若粗尺度单元尺寸与非均质尺寸相近会产生共振效应降低计算精度，此时需采用超样本技术[112]，即在大于粗尺度单元尺寸的超样本中求解临时基函数。多尺度有限元法可以较好地捕捉小尺度下大变形问题，进而解决传统FEM模拟裂纹扩展时网格重构和加密问题，适用于研究金属、岩石等细观破裂损伤评估[113]。相较而言，多尺度有限体积法通过全局守恒的通量进行计算，求解自由度低，因此被广泛应用于多孔介质流动换热、流固耦合问题中[114]。此外，若细观、宏观尺度求解控制方程相同时，多尺度有限元法、有限体积法较直接细尺度网格求解均大大降低了计算负荷。

4 井筒热流体高效提取技术

地热开采过程中，注入的冷流体在储层裂隙内流动并与高温岩石进行充分换热后，经过采出井传输至地面，进行后续温泉、供暖、发电等应用。EGS中，流体与岩石换热后进入采出井时的温度较高，同时由于井筒深度可达3～10 km，流体经井筒传输至地面的过程中将产生几十MPa的压降。在采出井口附近，若压力降低至流体饱和压力之下，将引发流体的闪蒸相变现象，极大地影响井内热流体的高效提取。因此，明晰地热井内流体闪蒸过程的流动换热特性，实现对井内流体闪蒸的精确预测，并进行有效预防是实现地热井内热流体高效提取的关键。笔者将从地热井内流体闪蒸的原理，以及井内闪蒸特性研究的实验和数值模拟方法3个部分，总结井筒热流体高效提取技术方面的研究进展和发展趋势。

4.1 地热井内的闪蒸

闪蒸现象主要可以分为以下2种：一种是将过冷流体的压力迅速降低到初始温度所对应的饱和压力之下，使其由过冷状态瞬间转变为过热状态，从而产生剧烈相变的闪蒸现象[115]；另一种闪蒸现象主要发生在循环系统中，由于流体在较长的上升段内举升时，其饱和温度随着流体静压的减小不断减小，当流体温度大于当地饱和温度时将发生闪蒸[116]，产生闪蒸点以下为单相流体，闪蒸点以上为气、液两相流体的流动状态。地热井筒内流体的闪蒸现象与后者类似。

流体中的热不平衡程度和汽化核心是诱发闪蒸汽化的充分条件和必要因素[117]，汽化核心可以是流体内部的杂质、不凝性气体或者从流动壁面产生。在流体举升过程中，一旦发生闪蒸，大量的流体显热转化为蒸汽潜热释放[118]，汽化核心处的流体产生汽泡形成泡状流，随着流体向下游流动过程中的静压逐渐减小，闪蒸稳定流动状态下的流型由上游的泡状流逐步向下游的弹（帽）状流和搅混流发展[117]（图11）。由于两相流动状态的换热系数要大于单相流动状态[119]，若井筒没有良好的保温措施[120]，井内流体闪蒸的发生将增大井筒与地层之间的换热损失，降低热流体出口温度，从而导致地热井的开采热量减产。此外，闪蒸的发生会导致井筒内流体的压力发生剧烈波动，诱发系统产生流动震荡[121]，甚至诱发水击现象[122]，威胁系统的安全性。

图11 闪蒸稳定流动状态下的两相流型图[117]

EGS和水热型地热系统均为开式系统，即从地层内部抽采出热流体，井内流体含有大量Ca2+等矿物质离子及二氧化碳等不凝性气体。当井内流体发生闪蒸汽化后，气相水蒸汽增多，导致二氧化碳分压降低，溶解在液相的二氧化碳随之逸出，造成最常见的碳酸钙结垢[123]（图12）。目前，美国[124]、匈牙利[125]以及中国西藏羊八井[126]、甘孜[127]、那曲[128]等地的中高温地热发电站地热井均有井筒结垢问题的相关报道。严重的井筒结垢问题将导致地热井的减产和关闭，制约地热能的可持续开发利用。墨西哥Los Humeros地热井由于井筒碳酸钙结垢问题，蒸汽产量不断下降，由投产时的38 t/h降至4 t/h[129]，羊八井、那曲、川西地区的地热井因为严重的结垢问题，需要长期依靠机械除垢来维持正常运行，甚至直接停产[128]。目前，国内外学者已经针对地热采出井筒内流体的闪蒸问题展开了大量研究，但由于采出井高井深的结构特点以及大压降、闪蒸相变、两相流动等复杂的物理特征，相关研究存在着实验研究难以有效展开、数值模型不够精确等难点。

图12 地热井筒及阀门的碳酸钙结垢图[123,127]

4.2 流动系统中的可视化闪蒸实验

针对3～10 km的地热井筒开展全尺度的实验研究显然是不切实际的，目前以地热井筒为实验对象，开展井筒内流体流动换热特性以及闪蒸现象研究的实验鲜有报道。值得关注的是，地热井筒内流体的闪蒸过程与自然循环系统中的类似，根本区别在于循环的驱动力不同，自然循环依靠流体自身密度差驱动循环，而地热井筒内流体的流动主要依靠泵功。

关于自然循环系统中流体闪蒸实验的研究目前广泛存在于核反应堆安全系统的设计中[130]，主要的实验思路为[117]：冷流体首先进入加热段，通过电加热等方式加热至设定温度；流体经充分加热后进入上升管路模拟自然循环中的流体举升过程，管路一般采用透明材料，结合高速摄像机可以开展可视化研究，管路沿程布置测温点、测压点测量流体的温度、压力参数；流体经上升段进入上部水箱，后经下降段进入加热段完成循环，沿程可布置电磁流量计等测量循环流量（图13）。

图13 自然循环系统中的可视化流体闪蒸实验平台图[117]

国内外学者基于上述实验过程开展了大量研究，包括对流体闪蒸流型演变的可视化研究、闪蒸的稳定性研究、闪蒸的流动换热特性以及对闪蒸进行控制等。Manera等[131]首次搭建出了垂直管道的流体闪蒸实验平台，并对管内静止状态与流动状态的流体瞬态闪蒸过程进行了可视化实验，采用丝网传感器实现了对气泡尺寸分布的测量，探究了不同流速下流体闪蒸流型的变化。闪蒸的发生将引起系统内部强烈的不稳定性流动，大大影响系统内部的阻力特性，郭雪晴等[132]通过可视化实验分析了自然循环系统中的瞬态运行特性和不稳定性机理，发现流体的降压闪蒸是影响系统稳定的流动驱动压头的一大根本原因。在流动特性研究的基础上，通过测量系统的流量和温度参数，可以对流动过程中的闪蒸换热特性进行研究，李鲁宁等[133]探究了不同过热度、压力、流量等工况下闪蒸换热系数的变化规律，张友森等[118]通过搭建循环闪蒸实验台研究了循环闪蒸显热转化率受流量、压力、液膜高度的影响规律。此外，在了解闪蒸流动换热特性的同时，预防闪蒸的发生也是备受关注的方向，研究表明，通过提升系统压力或减小阻力可以有效减少闪蒸现象的发生，使系统流动趋于稳定[134-135]。程俊等[136]还通过在上升段内设置内插物使系统流型实现了由流动不稳定状态下的间歇流向稳定流型的转变，有效地减小了系统的震荡幅度。上述研究可作为地热井内流体闪蒸流动换热特性研究以及预防井内流体闪蒸的参考，若要准确把握地热井内流体闪蒸的特征，仍需进一步开展地热井内的流体流动闪蒸实验。

参考自然循环系统中的闪蒸实验，要开展地热井筒内闪蒸过程的实验需要解决以下2个问题：①改变流体循环的驱动方式；②实现对高井深、大压降特点下地热井筒内流体流动状态的等效实验还原。Ishii等[137]根据小扰动原理，在单相模化方法基础上，应用漂移流模型推导出两相稳态自然循环模化准则，意图构建合理的小比例模化准则，在准确反应原型系统流动和传热特性的前提下，减小实验系统的几何尺寸和运行参数。郭雪晴等[138]对大尺度开式自然循环系统内的闪蒸过程进行了模化分析，所得模型系统的运行特性与原型实验结果吻合较好。由此可以推断，通过合理的模化准则来进行地热井筒闪蒸实验的设备尺寸设计及运行参数设计，即通过小比例的实验来实现对大尺度地热井筒内流体闪蒸过程的模拟，在设计上具有一定的可行性。

4.3 地热井内流动传热的数值模型

相较于实验研究，数值模拟方法在地热井筒方面的应用更为广泛，主要可以分为稳态下井内流动的模拟和瞬态下井内流动的模拟2个方面。

考虑到地热井筒存在高度方向尺寸与径向尺寸比例极大的结构特征以及井筒内可能存在的闪蒸相变及复杂两相流等问题，在进行数值模拟计算过程中往往需要引入一些假设：①将井筒简化为沿高度方向的一维物理模型，忽略径向的参数变化或将三维方程沿径向截面积分简化为一维的守恒方程[139]；②假设井筒内的流动状态为稳定状态[140]，忽略时间项导数；③假设井内为单相流体，对于两相的模拟常采用均相模型[141]，将流体看作是一种伪单相流体。对于稳态流动的模型，流量参数可以通过边界条件获得，在整个模拟计算过程中，只需要考虑两个主要变量即可求解整个模型。求解方法可以是从井口边界到井底的自上而下的求解（TOP-DOWN），也可以是从井底边界到井口的自下而上的求解（DOWN-TOP）[142]。压力、温度、流动质量分数或含气率常作为主要的变量用于稳态模型的求解[143]，其中质量分数无法用于含逆流的模型求解，因为当质量流量为零时流体质量分数是无法定义的，此时可选择采用含气率来求解相关模型[144]。压力和流体的焓也是常用的主要变量，当流体不存在相变时，流体的焓可以看作是温度、压力的函数，但地热井内存在闪蒸相变时，这种方法就不可行了，需要对方程进行离散，根据边界条件逐步求解各节点的压力、焓、流体质量分数等参数[140]。基于上述方法，学者们采用一维稳态模拟实现了井筒内流体沿程速度、压力、温度以及含气率等参数的快速计算，完成了包括采出井出口参数和井底储层环境参数的反演、井内流体闪蒸位置的预测及预防结垢等问题的模拟研究。卜宪标等[145]针对地热井内的结垢问题，采用DOWN-TOP的求解方法，根据已知的压力、温度、总质量流量、二氧化碳含量等井底参数对沿程压力、温度分布进行了计算，确定了井内流体的闪蒸位置，明确了影响井内闪蒸和结垢的主要参数。梁海军等[146]采用DOWN-TOP算法，计算获得了井筒沿程压力、温度分布，反演得到井底流体的温度参数，预测了流体闪蒸的起始位置，并提出通过调控井口压力和流量，可以调节闪蒸位置，优化流量和防垢成本之间的关系。表3给出了主要的地热井筒模拟器，Bjornsson等[147]开发出了第一个可以处理多进料层的稳态地热井筒模拟器HOLA，自此之后，学者们经过对模型的不断修正和开发，又相继提 出 了 GWELL[148]、WELLSIM[149]、MULFEWS[150]、FLOWELL[151]、SWELFLO[152]、PROFILI[153]等 一 系列求解器，可以满足多井、蒸汽井等不同地热生产条件下井筒流动传热过程的模拟计算。但闪蒸作为一种井内流体相态突变的现象，采用稳态模型计算不具备足够的说服力，因此对于井内流体闪蒸流动的模拟需进一步结合瞬态模型进行。

表3 主要的地热井筒流动模拟器统计表

对于瞬态下井筒内的流动模拟，Miller[154]在1980年开发出了第一个瞬态井筒模拟器WELBORE，后续相继有学者发展出用于求解地热井内两相流动的三方程模型，Akbar等[157]采用相关模型实现了高焓值地热井内闪蒸过程中流体性质的预测，同时考虑了浮力、相变、压缩性、热相互作用、壁面摩擦以及相间滑移等因素的影响。近来，Tonkin等[142]总结归纳了大量关于地热井筒流动的稳态和瞬态数学模型，指出了Akbar等一些学者所应用的动量/能量守恒方程中所存在的错误，并提出了一个完整的、物理上精确的三方程模型用于模拟地热井中的不稳定流动。此外，石油行业内井筒内的多相流动与地热井内的流动具有几何相似性和物理模型的相似性，目前在石油行业的油气井多相流模拟中，已经发展出瞬态的五方程模型和六方程模型[158-159]，因此对比石油行业中的前沿研究成果，发展用于求解地热井筒内复杂流动的数学模型是值得关注的研究方向。

5 干热岩地热发电技术

5.1 地热发电基本原理

通常来说，地热资源的利用方式取决于地热流体温度的高低，低温的地热资源（小于90 ℃）一般是对其进行直接利用，例如种植业、温泉等；对于中温的地热资源（90～150 ℃），分布式的制冷和供热是重要的利用方式；对于温度较高的地热资源（大于150 ℃），地热发电是实现能量利用的有效途径[160-161]。

地热发电原理是基于朗肯循环，利用高温蒸汽驱动膨胀机作功，实现热电转换。依据地热资源特性，地热发电厂的循环利用方式分为直接发电（干蒸汽发电系统、闪蒸发电系统）和间接发电（双工质发电系统）[162]（图14）。

针对高温水热型、干热岩型地热资源中含气率较高的地热流体的利用，一般通过降压闪蒸的方式获取高温蒸汽，直接驱动膨胀机作功，即干蒸汽发电（图14-a）、单级闪蒸（图14-b）、双级闪蒸（图14-c）。

干蒸汽发电系统结构简单、操作维修方便，从地下开采的高焓地热蒸汽经过必要的净化和脱酸处理后即可驱动膨胀机作功。但是，由于干蒸汽发电对地热条件要求较高，只能在少数地方使用而未能得到推广。相较于干蒸汽发电，闪蒸发电则具有更强的灵活性，通过降压闪蒸的方式，将地热水转换为饱和蒸汽，进而实现热电转换。闪蒸发电多适用于170 ℃以上的中、高温地热流体，目前在全球地热装机容量发电占据半数以上[163-164]。我国广东丰顺的邓屋地热发电厂即采用闪蒸发电模式，装机容量为300 kW[165]。此外，为了提高闪蒸循环的效率，在单闪蒸循环的基础上进行了双级闪蒸循环改进，即对一次闪蒸液相部分进行第二次闪蒸，获取更高的净输出效益，例如冰岛东北部的Krafla地热田[166]。针对闪蒸发电模式，Jalilinasrabady等[167]对伊朗的Sabalan地热发电厂单级闪蒸和双级闪蒸循环系统进行热力学和经济性的比较，研究表明双级闪蒸循环系统的净输出功率高于单闪循环系统，但双级闪蒸模式的投资较高，回收期较长，综合考虑热经济效益，双级闪蒸发电模式具有更加广阔的应用前景。

闪蒸发电的循环方式有结构简单，热效率较高的优点，但是存在一定的局限性：当地热流体温度低于130 ℃时，受到水的物性条件限制，闪蒸模式无法获取满足地热发电所需的高焓蒸汽，采用引入中间介质的双工质发电模式（有机朗肯循环[168]、卡琳娜循环[169]）是实现对中低温地热资源利用的有效途径。由于有机朗肯循环（Organic Rankine Cycle，ORC）属于闭口循环，纯净的低沸点工质通过换热器吸收地热能量，避免了地热流体对管路的腐蚀，因此ORC在中低温发电模式中展现出明显的优越性。在ORC循环中，根据地热流体的特性引入低沸点的工质，使其能够在较低地热温度条件下蒸发，进入膨胀机推动作功，乏汽进入冷凝器被冷凝，再由泵升压进入蒸发器，从而完成如图14-d所示的循环过程。随着深层地热资源的开发，EGS技术与ORC相结合的EGS-ORC技术逐渐展现优势。2009年法国苏尔士地热田率先建成了1.5 MW地热电站，采用异丁烷作为工质驱动膨胀机作功，并成功实现了并网运行[170]。与闪蒸发电模式相似，ORC循环衍生出的双级ORC循环发电具有更高的作功能力[171-172]，双级ORC通过对地热能的分段分品利用，有效降低了循环不可逆损失，提升了系统性能。在ORC循环中工质的选择非常关键，基于不同类型的地热资源，选择不同沸点的工质与热源进行匹配，能够有效提高地热能的利用率[173]。此外，基于闪蒸和双工质两种循环方式建立的闪蒸—双工质发电厂[174]，能够充分发挥两种循环模式的优势，实现对地热的温度分区利用：高温区间的地热流体用于闪蒸发电，闪蒸后的中温地热流体为ORC循环供热，从而完成对地热能的梯级利用，提高地热发电的效率。

根据地热资源特性选择合理的循环发电模式是实现提高热电转换效率的关键：干蒸汽、闪蒸等直接发电模式多适用于高温地热资源，ORC等间接发电模式则是中低温地热资源利用的重要方式。然而，目前针对各种发电技术的具体适用温度、压力及地质条件尚未明晰，仅以温度为技术选择标准存在局限性，且地热采出温度明显低于传统火力电厂，这将导致地热电效率低下。因此，改善地热发电技术、探明地热发电技术与地热资源条件的适配性是未来地热发电技术的发展趋势。

5.2 地热发电市场应用

随着石化行业钻井、压裂等技术的突破，对高温深层地热资源例如高温水热型、干热岩型地热资源的开发促进了地热发电的发展。1904年世界上第一座地热发电厂在意大利投建并成功运行[175]，地热发电的研究自此蓬勃展开，并因其清洁可持续的特点在世界范围广泛应用。除了意大利之外，美国、菲律宾、墨西哥、意大利、冰岛等地热资源丰富的国家的地热发电装机容量近年来迅速增加[163,176-177]，表4、图15所示为全球典型地热发电市场应用。截至2020年底，全球地热装机容量已达15 608 MW，其中，美国地热发电装机3 714 MW，占据世界首位[178]。在我国，高温的地热资源集中分布在西南地区，基于此条件，我国先后在西藏建立了羊八井、羊易、那曲、郎久等地热发电站，截至2021年西藏地区地热电站装机达42.18 MW，有效缓解了市区供电压力。80 ℃以上中低温地热资源的地热电站则是在广东丰顺、河北后郝窑等地展开，由于资源和技术限制，现仅剩广东丰顺的300 kW机组在间断运行[179-181]。

表4 全球典型地热厂统计表[182-186]

图15 典型地热发电厂现场图

针对干热岩的地热发电利用，目前尚未实现大规模的商业应用。20世纪70年代，美国洛斯阿拉莫斯实验室在新墨西哥州Fenton Hill地区[191]实施了全球首例干热岩项目，自此之后英国、日本、德国、澳大利亚、法国等国也相继开展了干热岩的项目建设[192]。值得关注的是，法国的Soultz[193]、德国的Landau[194]、澳大利亚的Habanero[195]、美国的Desert Peak[196]和The Geysers[197]等试验项目目前已初步具备1～5 MW的发电能力，实现了商业运营。上述案例也将为后续干热岩项目的规模化发展提供宝贵经验，表5对具有代表性的干热岩发电项目进行了总结。我国在干热岩开发利用方面仍处于资源勘探及EGS技术研究的初步发展阶段，目前已在青海、西藏、福建等地进行了地热资源勘探工作，其中在青海共和盆地发现高品质干热岩体，GR1井底最高温度236 ℃，孔内3 366 m以下深度平均地温梯度可达8.8 ℃/100 m。此外，2012年我国启动了“863”计划项目“干热岩热能开发与综合利用关键技术研究”，相关研究成果为我国干热岩资源开发利用提供了理论依据和关键技术支撑。但距离我国真正实现干热岩地热资源的发电利用，仍需进一步地解决本文所述EGS开发过程中的一系列关键技术问题。

6 结论与展望

干热岩开采中存在人工压裂缝网不合理、多尺度多场耦合规律不明、井筒流体闪蒸造成采热效率低下、地热资源热电转换效率低等问题，需要解决4个关键技术问题。

1）干热岩储层人工压裂技术：干热岩储层较为致密，渗透率极低，一般需要通过人工压裂的手段为流体创造流动换热通道，但不合理的压裂缝网易引发流体短路，造成过早的热突破，制约地热系统的可持续开发利用。应用暂堵转向压裂、液氮压裂等技术可帮助形成复杂的人工缝网，但仍需明晰缝网的形成机制，开发精确的裂缝扩展预测模型来指导压裂工作。对于干热岩水力压裂模型方向，将向能考虑干热岩热弹塑性变形与破坏、压裂液与干热岩相互作用、应力腐蚀引起的亚临界裂纹扩展等现象的全三维压裂模型发展；同时对于干热岩水力压裂算法，由于全三维压裂模型计算成本较高，水力压裂算法向最佳正交分解（POD）、广义最佳分解（PGD）、减基（Reduced Basis）等降阶算法发展，从而提高水力压裂数值模拟的运算效率。

2）干热岩开采数值模拟技术：干热岩储层具有多尺度、强非均质特征，开采过程受到热—流—固—化多场耦合作用的影响，属于典型的多尺度多场耦合问题，但相关耦合作用机理尚不明确。对于干热岩岩心尺度多场耦合模型，将向更精细刻画孔隙结构、更精确计算孔隙内部流动换热、更精确计算固体骨架局部变形及断裂扩展的方向发展，离散力学模型与新型流体力学求解方法如LBM、SPH将更深入的结合和应用，此外岩心尺度THMC耦合数值模型待需开发。对于储层尺度多场耦合模型，精确描述裂缝分布和走向的离散裂缝类模型是开发热点，各场控制方程之间会引入更全面的关联项和耦合方程进而符合实际物理过程，对于方程求解将大力应用XFEM、XFVM等先进方法以实现高效模拟。另一方面，尺度升级方法可以综合考虑储层中裂缝和孔隙多尺度特征，将与多场耦合模型结合，实现宏观和微观尺度耦合。

3）井筒热流体高效提取技术：地热井筒由于高井深、大压降的特点，内部流体在举升过程中易发生闪蒸，极大地影响井内的流动换热特性，引发流动震荡、井内碳酸钙结垢等问题，威胁地热系统的安全性和可持续开采。通过实验可以实现对地热井内流体闪蒸过程的直观描述，目前相关研究已在核反应堆设计的自然循环方面广泛展开，未来可通过建立合理的模化准则开展地热井筒内流体闪蒸的等效实验。对于地热井筒内的数值模拟，稳态、单相的模型已受到广泛应用，瞬态、两相模型可更准确地描述井内的闪蒸现象，已初步应用于地热井筒领域。此外油气开发领域具有更加成熟和前沿的井内两相流动模型，可扩展应用于地热井筒。同时，将地热井筒全尺度的数值模拟与局部区域的可视化实验研究相结合，可更好地实现对地热井内流体闪蒸机理的深入探究，为预防井内流体的闪蒸以及提高井筒热流体提取效率提供理论支持。

4）干热岩地热发电技术：地热发电的主要方式包括直接发电模式和间接发电模式，适配高温和中低温的地热资源条件。近年来，地热发电技术以其清洁可持续的特点在全球范围内得到广泛应用。然而，目前在地热发电技术的选择存在较大局限性，地热压力、流体成分、地质条件等也应作为技术选择的重要参考。针对地热资源利用效率低的问题，多能互补和能量梯级利用的分布式能源利用将是提高能源利用率的有效途径。此外，对于地热发电的评估不应仅限于电厂的投建运行，以全生命周期的角度将地热勘探、开采、利用等环节的投入与电力输出目标效益结合，明晰系统运行周期内的经济收益、“碳”收益，为地热发电的产业化发展提供指导。

5）EGS是一个技术密集型系统，对其工程开发需先从机理分析和可行性上着手。然而，由于地下储层工况的复杂性和地面设备的不稳定性，机理研究往往与实际脱节，需与先导试验项目密切结合、互相指导，开展产—研结合模式，在不断往复中提高认知、突破关键点，提升EGS的应用价值，为我国实现“双碳”目标助力。