封闭地热系统水平井采热储层内渗流与传热的有限元模拟

2023-01-10张连平王世民

地球物理学报 2023年1期

张连平，王世民*

1 中国科学院大学地球与行星科学学院，北京 100049 2 中国科学院计算地球动力学重点实验室，北京 100049

0 引言

地热能是来源稳定、储量巨大的可再生清洁能源.开发利用地热能既能有效应对能源短缺的挑战，也能大幅降低碳排放，因而愈来愈得到重视(Zeng et al., 2016; 汪集暘等, 2017; Chen et al., 2019; Li et al., 2019; Liu et al., 2019; Shi et al., 2019).地热能开发的重中之重是地热发电，其能源利用系数远高于水力、风力、太阳能发电和地热直接利用(汪集旸, 2012).由于地热电站可不受季节、气候、昼夜等自然条件的影响而不间断运行，地热发电可以作为国家电网的基础载荷，也易于调峰和实施热电联供(庞忠和等, 2012; 王晓星等, 2012).我国地热能资源十分丰富，但目前的开发主要限于中低温地热能的直接利用，地热发电没有得到应有的发展(Huang, 2012; Feng et al., 2014; Lei et al., 2019; Pan et al., 2019).

地热资源按其成因和产出条件可分为水热型和干热岩型(汪集旸等, 2012; Jiang et al., 2014).水热型地热资源赋存于高渗透性的孔隙或裂隙介质中，与年轻火山活动或高热流背景相伴生形成高温水热系统，而处于正常或偏低热流背景下的地下水循环通常形成中低温水热系统.干热岩型地热资源赋存于地下较深处的高温但低渗透性岩体中，原则上不受地区分布限制，但其开采需要借助人工压裂进行储层改造(reservoir stimulation)，进而通过注水循环形成增强型地热系统(Enhanced Geothermal System, EGS)(Huang et al., 2017; Zhang et al., 2019; 丁军锋, 2018; 丁军锋和王世民, 2019).

EGS发电是一个有广泛应用前景、值得深入研究的方向.近年来EGS 的研究和开发广受关注，在美国、英国、法国、德国、瑞士、日本、澳大利亚等多个国家已经进行了EGS 工程试验(王晓星等, 2012; 许天福等, 2012)，但商业运行EGS 电站的成功案例还很少(汪集旸, 2012).我国EGS 的研发工作起步较晚，但针对EGS 资源评价(Liu et al., 2022; 汪集旸等, 2012; Wan et al., 2005)、工程选址(苏正等, 2014; 郭进京和周安朝, 2010; Feng et al., 2012, 2014)以及采热过程数值模拟(Jiang et al., 2013, 2014; Chen and Jiang, 2015; Luo et al., 2014; Zeng et al., 2014; Huang et al., 2014; Zhang et al., 2014)等方面的研究工作近年来愈益深入.

在中国大陆地区，高温水热系统主要分布在喜玛拉雅造山带附近，其他地区则以中低温水热系统为主.目前制约中低温地热发电的主要问题之一是钻井成本高(约占电站总成本的 60%以上)(李克文, 2010).利用油田伴生中低温地热发电，不但可以很好地降低成本，而且油田热储的回注或回灌系统大部分已经探明或者已经成功运行，是经济开发中低温地热发电的优选方向.结合深部钻井成本高昂与我国存在大量废弃油气井的实际情况，国内学者对改造废弃油气井进行地热发电开展了一系列模拟研究(李克文等, 2012; Gong et al., 2011, 2013; Xin et al., 2012; Sun et al., 2018b)，在华北潜山油田还实施了废弃油井地热发电的试验工程(李强和于红侠, 2012).

潜山油田试验电站在运行过程中发现存在三个主要问题(李强和于红侠, 2012).第一，管道和换热器中结垢严重.由于地热水矿化度较高且含有多种腐蚀离子，结垢后不但严重影响发电效率，长期运行还可能造成设备损坏甚至引发安全事故.第二，地热水流量不稳定且达不到额定流量，造成实际发电功率波动(160～240 kW)且远低于额定发电功率(360 kW).第三，电站自耗功率(80 kW)过高，高达实际净发电功率(80～160 kW)的50%～100%.

这些问题并不是潜山油田试验电站所独有，在其他地热电站实际运行中同样可能遇到.如西藏羊八井地热电站同样存在腐蚀结垢严重、地热水流量不稳定及压力损失较大等问题(周大吉, 2003)，而西藏郎久、那曲和台湾清水、土场的地热发电厂更由于结垢及流量不足等原因在运行不长时间后就停运(郑克棪等, 2010).未来投入商业运行的EGS 电站也会面临类似的问题，如法国苏尔士(Soultz-sous-Forets)EGS 试验工程在双工质循环发电设备的换热器系统和地热水回注系统均发现较严重的结垢和腐蚀现象(Mundhenk et al., 2013)，而能完全除垢或有效阻垢的工艺和系统还没能建立起来.因此，对这些在地热电站实际运行中具有普遍性的技术问题，需要进行深入的研究以找出切实可行的解决之道.

事实上，潜山试验电站运行中出现的几个主要问题都可归因于其采用的开放性采热模式，即通过开放系统的强迫对流方式，以压力驱动地热水流过地热储层的孔隙结构再循环到地表.在这种模式下，地热储层岩石孔隙中的水岩相互作用导致矿物成分溶解于地热水，而地表不同的温度、压力和化学条件可导致溶于水的矿物析出，造成结垢现象.由于地热储层的孔隙结构非常复杂，流体运移通道不稳定、渗透性较低且高度不均匀，驱动地热水在孔隙结构中流动不但要消耗可观的能量，而且流量难以稳定控制.另外，从废弃油气井中抽取的地热水常含有一定比例的油，除油过程需要消耗电能，增加了电站自耗功率.

然而，开放性采热模式不是地热电站的唯一选择(尽管这种采热模式被大量实际运行的地热电站普遍采用，并且也是文献中模拟未来EGS电站运行的主要采热模式).如果用封闭性采热模式取代开放性采热模式，前述潜山试验电站运行中出现的几个主要问题都可以从根本上得到解决.在封闭性采热模式下，传热工质在封闭的管道中循环.流经地热储层内的管道时，管道内的传热工质与管道外的地热储层不直接接触，管道内外只有能量交换但没有物质交换.这样，自然从源头上消除了结垢问题.同时，与驱动孔隙流动相比，封闭管道中的流体循环要求低得多的驱动压力，流量也非常易于控制.在封闭性采热模式下工作的油井地热电站当然也不需要消耗电能来除油.

与开放性采热模式相比，封闭性采热模式具有灵活选择传热工质的优势.开放性采热模式下运行的水热型地热电站采用水作为传热工质是必然的选择.现今运行的EGS 工程也普遍采用水作为传热工质，则很大程度上可归因于这些EGS均处于试验探索阶段，储层改造与采热开发往往需要多次交替进行.由于储层改造大都采用水力压裂的方式，其中水是压裂过程中的压力载体，在采热开发阶段选择水作为传热工质就是一个自然的选择.近年来，同样是采用开放性采热模式但以CO2作为传热工质的EGS模拟研究也不断发表于国内外文献(Luo et al., 2014)，引起一定程度的关注.采用封闭性采热模式为地热电站选择传热工质开辟了更广阔的空间.由于无需担心传热工质与储层岩石间的化学反应，在水与CO2之外，多种流体或流体混合物以及纳米流体也可供选择，只要其在地热电站工作条件下具备优良的热物理性质、与管道壁材料化学兼容并且价格可承受即可.

封闭性采热模式的另一个重要优势在于可以根据实际需要对传热流体工作状态进行设计和调整.这样的自由度在开放性采热模式下是不存在的，因为开放系统中传热流体压力是与地热储层的孔隙压力耦合的.由于储层孔隙结构的高度不均匀性和复杂性，人为控制地热水在地热储层中的热力学状态几乎是不可能的.然而在封闭性采热模式下，可以很容易控制传热流体的压力以满足实际需要.

封闭性采热模式还可以设计成两相流采热系统，为利用汽化热提高地热能开采效率提供了可能.众所周知，单相流只能通过流体升温来吸收热能，而气-液两相流还可以通过汽化热来吸收热能.例如，水在150 ℃的气化潜热为2113 kJ·kg-1(Cengel and Boles, 2008)，相当于液态水温度增加500 ℃所吸收的热量.应用于地热发电，利用气-液两相流开采同样多的地热能所需要的质量流量比单相流至少小一个数量级.而且，由于热流与温度差成正比，水在汽化过程中只吸热而不升温，有利于维持管道内外的温度差，从而提高换热效率.

封闭性采热模式实际应用于地热发电工程时，地热水在地热储层中的循环不再如开放模式那样由强迫对流驱动，而是由自然对流驱动.也就是说，封闭地热系统只能借助于自然对流强化采热，自然对流的流量将直接限制地热电站的采热功率与发电功率.因此, 自然对流如何影响采热效率对封闭地热系统来说具有根本的重要性.

虽然开放地热系统主要依靠强迫对流采热，前人研究表明地热储层内的自然对流对提高开放地热系统采热效率亦有显著作用.Bataillé 等(2006)对法国苏尔士EGS 地热储层中的自然对流和强迫对流模拟表明，自然对流的流量与强迫对流的流量属同一数量级，而且自然对流对维持地热田温度进而延长地热电站寿命有非常显著的作用，因为自然对流使得注入的冷流体更容易流向储层底部，增加了地热系统的有效采热面积.Cao等(2016)对比分析水工质和超临界CO2工质的开放地热系统，后者自然对流更强，对储层底部的采热能力也更强.Jiang等(2019)研究发现，自然对流发育的开放地热系统中，可以通过适当降低注水井位置并抬高生产井位置优化地热系统，提高采热效率.Liang等(2018)模拟具有断层发育储层的开放地热系统表明，采热效率并不总是随着储层渗透率增加而提高，渗透率过高反而导致地热系统的采热效率下降.这些基于开放地热系统的研究结果为自然对流能够提供足够的流量以满足地热电站在封闭性采热模式下工作提供了有力的支持.

封闭式地热系统具有多种分类方式.既可从规模角度分为大型系统和小型系统，从传热工质角度分成水工质系统、二氧化碳系统、氟利昂系统等，也可依据钻井数量分为单井型、多井型和井阵型，还可根据采热流体管道的形态和结构分成内外双层套管型、U型、多分支型、螺旋型等.

前人针对小型封闭地热系统已开展一系列研究工作.Gustafsson等(2010)对垂直单井U型管道封闭地热系统的研究表明，U型管道外井筒中的自然对流可以有效降低系统热阻，从而提高采热效率.Ghoreishi-Madiseh等(2013)对垂直U型封闭地热系统的研究发现，当渗透系数大于10-5m·s-1时，自然对流对采热有显著影响，并且自然对流对井阵式地热系统采热效率的提升效果优于对单井式地热系统.Bringedal等(2013)模拟垂直单井封闭地热系统发现，储层初始温度分布对自然对流强化采热效果有重要影响.Shi等(2018)研究垂直单井封闭地热系统表明，自然对流对小孔隙度储层采热有强化作用.Wang等(2021)对具有水平井段的单井内外双层套管封闭地热系统数值模拟表明，储层渗透率越高、水平井段长度越长，自然对流强化采热效果越好.此外，对于家庭供热制冷的小型单井型封闭地热系统，文献中已报告较多研究工作(Lyu et al., 2017, 2018; Song et al., 2017; Shi et al., 2018; Lee, 2019).

迄今为止，针对地热发电应用的大型封闭地热系统研究还很少.虽然德国联邦政府“投资未来项目”曾对“地下闭式循环热交换系统”地热发电从钻井与完井技术角度进行过探索和研究(郑秀华等, 2004; 周刚等, 2006)，但对地热电站在封闭性采热模式下工作所涉及的许多关键技术问题(例如，生产井和回注井的布局、地下换热系统结构设计与性能优化、采热系统运行过程模拟等)都还没有开展系统深入的研究工作.尽管针对封闭地热系统管道内流动如何影响采热效率已开展不少参数化(管道直径、长度、管内传热流体类型及注入温度、注入速率等)研究工作(Jiang et al., 2016; Cui et al., 2017; Wu et al., 2017; Song et al., 2018; Sun et al., 2018a,b, 2019; Wang et al., 2019)，但这些工作仅考虑储层中以传导方式传热.Oldenburg等(2016)研究了U型多井封闭系统自然对流强化采热问题，模拟了在水平井附近均匀提高储层渗透率对采热的影响，但没有系统分析储层厚度-宽度比、水平井井位、储层改造范围、非均匀渗透率强化改造等其它关键因素对采热的影响.

综上所述，封闭性采热模式应用于地热发电具有显著的技术优势和巨大的经济潜力，但对大型封闭地热系统采热性能的基础研究工作还很薄弱.本研究针对如图1所示适用于地热发电的大型封闭地热系统(多井U型封闭系统)(Schulz, 2008)开展有限元数值模拟工作.该地热系统的地下部分包括注水井、水平井和生产井，但采热主要由水平井段完成(Song et al., 2018; Sun et al., 2018a, 2019).地热采热过程是一个随时间演化的空间三维过程，采用三维模拟从理论上讲当然更理想.然而在数值模拟实践中，三维模拟常常由于计算量太大而达不到要求的精度(Sun et al., 2017; Ding and Wang, 2018).由于地热储层内的自然对流计算需要耦合流体流动和传热过程涉及的多物理场，而且计算区域具有水平井直径远远小于其长度的几何特征，要得到精确可靠的数值结果，计算所用网格尺寸必须足够小，时间步长必须足够短，导致三维模拟计算量巨大.事实上，对本文研究的问题来说，由于水平井直径远远小于其长度，温度场和流场在水平井走向上变化缓慢，而在垂直于水平井走向的截面内变化剧烈.因此，一个自然合理的选择是在垂直于水平井走向的截面内对耦合的渗流与传热过程进行二维模拟，并以水平井温度的不同取值代表从上游到下游不同的采热位置.这样处理，不但保留了自然对流影响水平井采热性能的基本物理要素，还可以节省计算资源以保证计算精度和全面系统地探查各种模型要素的影响.

图1 封闭地热系统示意图

本研究对封闭地热系统自然对流强化水平井采热性能的物理机制进行系统深入的二维数值模拟.研究中首先对一个参考模型进行详尽的分析，而后通过改变水平井温度、水平井位置、储层渗透率、水平井周围渗透率强化范围与储层厚度-宽度比，系统比较关键模型参数对水平井采热热流和采热年限的影响，从而为适用于地热发电的封闭地热系统设计和工程实践提供理论指导和基础数据.

1 有限元模型

本文参考模型计算采用的有限元模型如图2所示.设定水饱和地热储层厚度h=1000 m，宽度w=200 m，孔隙度φ=0.01，基础渗透率k=5×10-12m2.

储层中流体的密度随温度升高而降低，满足如下的状态方程(Nield and Bejan, 2017)：

ρf=ρ0[1-α(T-T0)],

(1)

其中，ρ0为参考温度T0对应的流体密度，α为流体热膨胀系数.

根据Oberbeck-Boussinesq近似(Nield and Bejan, 2017)，仅需在动量方程(达西定律)中的热浮力项考虑密度随温度变化，从而储层中的自然对流在图2所示的坐标系下满足如下形式的不可压缩连续性方程、达西定律和能量方程(Nield and Bejan, 2017)：

图2 本研究采用的有限元模型

(2)

(3)

(4)

其中，u和v分别是水平和垂向达西速度，p为孔隙压力，T为储层温度，g为重力加速度.流体和岩石骨架热物理性质的取值列于表1.储层的等效热容和等效热导率采用体积平均法(COMSOL, 2016)计算：

表1 模拟采用的流体和岩石骨架热物理性质取值

(ρCp)eff=φρ0Cpf+(1-φ)ρsCps,

(6)

λeff=φλf+(1-φ)λs,

(7)

下标f和s分别对应流体和岩石骨架.

参考模型模拟计算包括稳态背景自然对流计算和瞬态采热过程模拟两个步骤.首先计算在地热系统采热之前储层中的稳态背景自然对流.计算中设定储层四周均为不可渗透边界，储层底部和顶部温度分别固定为475 K和445 K而左右边界绝热，储层设为均匀渗透率5×10-12m2.这一步得到的结果为瞬态采热过程模拟提供初始条件.

瞬态采热通过位于储层中心对称线上且距离储层底部 900 m 处直径为0.3 m的水平井实现.瞬态模拟中，设定以水平井为中心、以30 m为半径的区域为储层渗透率强化改造区域.强化区内的渗透率指定为从强化区外的基础渗透率5×10-12m2线性递增到水平井壁处的10-9m2.这样的渗透率变化比Oldenburg 等(2016)采用的局部渗透率均匀提高的假设更符合实际情况.瞬态计算中，设定水平井井壁边界条件为恒定温度343 K且不可渗透，储层四周边界均设为不可渗透的绝热边界.由于水平井是通过自然对流方式冷却储层并从储层内采热，储层与围岩间的热传导对采热贡献甚微，因而储层四周设定为绝热条件是一个合理近似.

本研究采用多物理场耦合有限元软件COMSOL(COMSOL, 2016)进行数值模拟计算.为了保证自然对流的计算精度(Otero et al., 2004; Hewitt et al., 2012, 2013, 2014; Wen et al., 2015; Nield and Bejan, 2017)，在温度和达西速度梯度大的水平井附近(图2b)和储层边界区域(图2c)采用逐渐加密的有限元网格.参考模型有限元计算网格节点总数为328,544，有限单元总数为642,464，包括631,496个三角形单元和10,968个四边形单元.由于温度场和渗流场在采热开始阶段变化非常快，瞬态计算中在第1年内的采热模拟时间步长从万分之一年逐渐增加到最大时间步长0.01年.为了验证模拟结果的精度，我们在时间步长和网格尺寸减半后重新计算参考模型，得到的结果与本文报告的参考模型结果最大偏差在0.5%以内.

2 参考模型计算结果

参考模型稳态计算得到的背景自然对流结果如图2d所示.由于模型关于x=0左右对称，图中只绘出了温度和达西速度在中心对称线左侧的分布以及两者在中心对称线上的变化.图2d表明，在模拟的条件下，背景自然对流虽然发育并在储层顶部和底部边界附近呈现速度和温度边界层特征，但对流很弱，最大达西速度仅为10-7m·s-1量级.

参考模型瞬态计算模拟了100年的采热过程.将计算给出的水平井井壁热流和切向达西速度(法向不可渗透)沿井壁作积分平均，得到的井壁平均热流和平均达西速度结果由图3给出.图3右侧两个插入图分别给出热流和达西速度在采热0.2年、1年、5年、20年、40年和100年时从水平井底部(θ=0)沿左侧(或右侧)井壁到水平井顶部(θ=π)的变化.六个时刻的井壁热流和达西速度均呈现相似的总体趋势，即井壁热流自水平井底部向水平井顶部单调递增，而井壁切向达西速度从水平井底部的零值逐渐增加到圆心角稍小于π/2处的最大值、再逐渐下降到零值.这样的总体趋势可由水平井周围渗流场和温度场的总体特征解释，如图4所示采热1年时的情形.图4a的流线显示流体围绕水平井作自上而下的绕流，水平井顶部正对来流从而速度边界层最薄，随后速度边界层从上向下逐渐增厚，直到水平井底部达到最大厚度.这样的绕流特征控制了热边界层同样从上向下加厚，在图4b中表现为水平井顶部温度梯度最大，而水平井底部温度梯度最小.

图4 参考模型给出的采热1年后水平井周围流线(a)与温度分布(b)

在图3显示的模拟时间范围内，井壁平均热流先后经历了迅速增加、短暂下降、再次迅速增加而后缓慢上升、最终逐渐递减的过程.在采热27年时平均热流达到峰值(12.9 kW·m-2)，约为模拟结束时(采热100年)的两倍，而在采热的前47年内平均热流均超过10 kW·m-2.井壁平均热流与平均达西速度具有非常相似的演化趋势，说明传热过程受自然对流控制.如图3左侧插入图所示，井壁平均热流在采热0.24年到0.48年间有短暂下降，而平均达西速度下降先于平均热流下降发生，清楚地指示了自然对流对传热的控制作用.

图3 参考模型给出的水平井壁平均热流和达西速度随时间演化

水平井的冷却效应导致井下方的低温区具有较大的向下流动速度，形成一个冷柱.在采热0.2年、1年、5年、20年、40年和100年6个时刻，图5给出的储层温度分布及温度与达西速度在中心对称线上的变化清楚地揭示了冷柱形态的时间演化.在采热最初阶段，冷柱迅速向下扩展直到抵达渗透率强化区底部，形成井壁达西速度与热流最初的快速增加阶段.在采热0.24年到0.48年间，自然对流主要通过冷却渗透率强化区来增强停留在强化区内的冷柱，导致井壁达西速度与热流出现短暂下降.采热0.48年后，加强后的冷柱能够在强化区外持续扩展，使得井壁平均达西速度迅速增加到8.25×10-4m·s-1水平并维持20余年，直到采热30年后由于储层底部温度显著下降才导致冷柱强度降低，而储层内自然对流也逐渐减缓.相应地，井壁热流在采热0.48年后也经历了先迅速增加、再缓慢增加、最终逐渐衰减的过程.

图5中不同冷却时间对应的温度分布特征揭示出储层冷却的先后次序.首先冷却渗透率强化区内水平井下方区域，继而在水平井下方低渗透区域自上而下逐渐冷却，然后是储层两侧由下向上逐渐冷却，最后才是水平井上方从两侧向中心的冷却.这一顺序体现了自然对流区域逐渐拓展并从上游到下游区域渐次冷却储层的总体特征.

图5 参考模型给出的6个时刻储层左侧温度分布及中心对称线上温度与达西速度变化

参考模型结果表明，采热阶段的最大达西速度可达mm·s-1量级，比背景自然对流高4个量级.考虑到采热阶段132 K的系统总温差只有背景自然对流30 K总温差的4.4倍，驱动自然对流的温差与自然对流强度之间显然不是一个简单的线性关系，值得进一步深入研究.此外，水平井位置、储层厚度-宽度比、储层渗透率及其在水平井周围的强化改造范围，都会影响采热效率.为此，我们对各模型要素如何影响采热效率逐一进行讨论.

3 讨论

为方便讨论，后文将参考模型称作MR模型.表2提供了MR模型与其他14个模型的参数对比，表中将MR模型参数以及其他模型与MR不同的参数加黑表示.在MR模型基础上，通过改变水平井温度(MT1和MT2两个模型)、水平井高度(MH1和MH2两个模型)、储层基础渗透率(MK1和MK2两个模型)、渗透率强化区半径(MKR1、MKR2、MKR3和MKR4四个模型)和储层宽度(MW1、MW2、MW3和MW4四个模型)，形成其他14个模型.这15个模型计算结果的系统对比，可以揭示封闭采热系统中自然对流强化水平井采热的基本物理机制.

表2 计算模型总览

3.1 水平井壁温度的影响

封闭地热系统的传热流体在水平井段被逐渐加热，因而水平井壁温度从上游到下游逐渐升高.模型MR、MT1 和MT2分别设定井壁温度343 K、393 K和443 K，计算结果对比由图6给出.相对于采热开始时储层底部温度475 K，三个模型分别以132 K、82 K和32 K总温差驱动自然对流.图6a表明等差递减的驱动温差导致的水平井壁平均热流以非线性的方式递减：在采热的前40年内，MR热流大于12.0 kW·m-2，MT1热流略低于5.5 kW·m-2，而MT2热流仅有0.5 kW·m-2.图6b则揭示在采热40年后，MR和MT1均能通过有效的自然对流冷却水平井下方的整个储层，而MT2仅在渗透率强化区内发育.

图6 水平井壁温度对采热的影响

3.2 水平井高度的影响

在自然对流发育的封闭地热系统中，水平井优先冷却其下方的储层，与热传导向四周同步冷却的方式显著不同.这就提示在自然对流主导传热的储层中水平井高度会对系统采热产生影响.图7比较了水平井位高度分别为900 m(MR)、800 m(MH1)和700 m(MH2)三种情况的计算结果.图7a显示，水平井位越高，高热流持续时间越长，虽然热流峰值略有下降.高于10 kW·m-2的井壁平均热流在水平井高度900 m情况下可持续时间47年，比在水平井高度700 m情况下长10年.图7b显示采热40年时的温度和流场特征，显然水平井井位越高，系统有效采热面积越大，地热系统对储层高温区的利用程度越高，而主要集中在水平井上方的高温区域基本没有冷却.图7结果表明将水平井井位尽可能设置在储层高处，可以获得更持久的高热流，达到高效利用地热资源的效果.

图7 水平井位置对采热的影响

3.3 储层渗透率及其在水平井周围强化改造的影响

本文参考模型考虑了在储层基础渗透率5×10-12m2背景下，在水平井周围半径30 m的区域内渗透率从基础渗透率线性递增到水平井壁处的10-9m2.为了比较储层背景渗透率整体提高和水平井周围局部强化的不同影响，图8将基础渗透率分别为5×10-11m2(MK1)和5×10-12m2(MK2)但未进行局部强化改造的计算结果与参考模型对比.图8a显示在整个100年的模拟时间范围内，相对于基础渗透率低且未强化改造的MK2模型，基础渗透率整体提升一个量级的MK1模型虽然可将水平井井壁平均热流由1 kW·m-2水平提高到3 kW·m-2左右，但远没有只在水平井周围局部改造的参考模型采热效果好.这与前人的研究结论(Oldenburg et al., 2016)一致.图8b显示在采热40年后，参考模型对未改造的低渗透区域的冷却程度也远高于整体提升基础渗透率一个量级的MK1模型.这是因为水平井周围的局部高渗透区域可以有效帮助自然对流启动并加强，而一旦自然对流发展起来，就可以继续扩展到低渗透区，完成对低渗透区的采热.

图8 储层渗透率对采热的影响

图8结果对封闭地热系统工程选址和设计具有特别重要的指导意义.因为基础渗透率特别高的地热储层非常稀缺，整体提升基础渗透率的技术手段也非常有限，而通过水力压裂等技术大幅提高水平井周围局部渗透率的工程技术手段则很成熟.图8提示利用局部高渗透区启动并强化自然对流进而完成整个储层的采热是一个非常经济和高效的方案.

为了深入研究渗透率局部强化范围的影响，图9对比了强化半径分别取作1 m (MKR1)、10 m(MKR2)、30 m(MR)、50 m(MKR3)和100 m(MKR4)5个模型的计算结果.渗透率变化幅度与参考模型相同.图9a表明高渗透率区域半径为10 m (MKR2)和30 m(MR)两种情形采热效果各有优势，在采热27年到66年之间MKR2井壁热流高于参考模型MR，其他时间参考模型MR热流更高.这两个种模型均明显优于其他三个模型.相对而言，强化半径50 m的MKR3模型优于强化半径100 m的MKR4模型，而强化半径1 m的MKR1采热效果最差.这样的结果说明水平井周围渗透率局部强化改造存在最优范围，强化区域过小或过大都会使采热效果变差.图9b清楚地表明，强化半径1 m的MKR1模型因改造区过小，限制了自然对流在强化区内的启动和加强，导致在低渗透区域采热能力不足；而强化半径分别为50 m和100 m的MKR3和MKR4模型则因改造区过大，水平井正下方冷柱区域垂向流动过快、形成流体优势通道，导致冷却区域形成瓶颈形态，而在瓶颈外侧留下采热盲区.结合现今局部渗透率强化改造的施工能力，图9结果提示在水平井周围半径10 m范围内进行水力压裂提高渗透率应该是兼顾采热效果、技术难度和工程成本的最优设计.

图9 储层渗透率强化半径对采热的影响

3.4 储层宽度的影响

在给定储层厚度和温度梯度的前提下，储层宽度决定了储层内能否发生背景自然对流(Nield and Bejan, 2017).对本研究模拟的厚1 km温差30 K的地热储层，图10对比了厚度-宽度比分别为5(MR)、6(MW1)、7(MW2)、9(MW3)和10(MW4)5个模型的计算结果.其中，MR、MW1 和MW2三个模型以稳态背景自然对流作为计算初始条件；而MW3和MW4两个模型因储层厚度-宽度比太大背景自然对流不能发生，故采用稳态热传导解作为计算初始条件.厚度-宽度比等于8的储层内因背景自然对流不稳定，导致瞬态计算初始条件不能唯一给定，没有加入图10的模型对比.

图10 储层宽度对采热的影响

图10a显示随着储层宽度由200 m逐渐减小为100 m, 水平井壁热流相应递减，较高热流持续时间也逐渐递减.厚度-宽度比分别为5和6的MR和MW1模型在采热开始阶段呈现井壁热流快速增加的趋势，而背景自然对流不存在的MW3和MW4两个模型在采热开始阶段呈现井壁热流快速下降的趋势.两种趋势的截然差别反映了背景自然对流对于提高采热效率的贡献.

图10b几个模型结果的横向对比可以体现储层冷却的先后次序.由于不同宽度的储层地热储量不同，同样是冷却40年，储层最宽的模型还处于储层两侧由下而上的冷却阶段，而储层较窄的两个模型已经处于在水平井上方从两侧向中心冷却阶段.与图5结果相结合，图10b表明自然对流冷却储层的先后次序并不随储层宽度变化而改变.这样的冷却次序对封闭地热系统选址和设计具有指导意义.例如，对同样截面面积，高而窄的储层优于宽而薄的储层.增加储层厚度不但有利于增大储层内温差，强化背景自然对流，而且有效增加了自然对流优先采热的储层区域.

3.5 模型验证与参数取值讨论

验证数值模型的标准步骤是基于该模型求解经典基准(benchmark)问题，然后与已有的可靠理论解和数值解进行对比.为此，我们选择等温水平圆柱驱动的稳态孔隙自然对流这一基准问题考核本文有限元模型.在指定储层渗透率均匀且四周边界恒温的条件下，我们采用本文参考模型网格计算了水平井冷却储层的稳态自然对流问题.得到的有限元解可与前人给出的边界层理论解和有限差分数值解直接对比.对Rayleigh数Ra=10情形，三种不同方法给出的热流沿水平井井壁变化由图11绘出，其中井壁热流以无量纲Nusselt数Nu表达，而Ra与Nu的定义与Ingham和Pop(1987)一致.图11表明，在从水平井底部θ=0到顶部θ=π的整个井壁上，本文有限元解得到的Nu值处处介于前人边界层解(Merkin, 1979)与有限差分解(Ingham and Pop, 1987)之间，且与有限差分解的偏差非常小.对图11中沿水平井井壁变化的Nu作从θ=0到θ=π 的积分平均，本文有限元解给出井壁平均Nusselt数1.551，稍小于有限差分解给出的1.578，而两种数值解给出的平均Nusselt数均显著高于边界层解预测的1.2634.因为边界层解采用只有在Ra趋于无穷的极限情形才严格成立的近似，其应用于有限Ra自然对流会普遍倾向于低估热流(Nield and Bejan, 2017).在其他Ra数取值情况下，本文数值模型的可靠性同样得到了前人结果的验证.

图11 基于本文参考模型有限元网格计算Ra=10稳态自然对流得到的沿水平井井壁热流变化，及其与前人边界层解(Merkin, 1979)和有限差分解 (Ingham and Pop, 1987)的对比

本文模拟计算仅考虑垂直于水平井走向的二维平面内渗流与传热，忽略了沿水平井走向的三维效应.对这一简化近似的合理性有必要进行分析和讨论.由于水平段井筒内流体温度会从上游到下游逐渐升高，井筒周围储层温度也会相应地升高，导致自然对流包含一定程度的沿水平井走向渗流.以z和w分别表示从水平井上游到下游方向的坐标分量和达西速度分量，达西定律给出：

(8)

这说明水平井走向的达西速度大小正比于储层渗透率和z方向的孔隙压力梯度.水平井周围的渗透率强化将使得水平井走向的渗流主要集中于渗透率强化区，而将方程(3)、(4)和(8)表示的达西定律代入连续性方程(2)，可以得到孔隙压力满足泊松方程

(9)

方程(9)表明，储层内孔隙压力变化源于垂向温度梯度变化，正梯度区形成压力汇，而负梯度区形成压力源.在垂直于水平井走向的x-y平面内, 由于水平井周围渗透率强化区温度最低，导致在强化区上部(垂向温度梯度为正)形成压力汇，而在强化区下部形成压力源.随着水平井从上游到下游温度逐渐升高，压力源和压力汇的强度都将逐渐减弱，从而沿水平井走向的渗流在强化区上部表现为从下游弱汇向上游强汇的流动，而在强化区下部表现为从上游强源向下游弱源的流动.由此可以推断，平行水平井走向的渗流速度必然弱于垂直水平井走向的水平渗流，因为后者是压力源(汇)区与无源(汇)区之间的流动.上述定性分析的结论得到了定量计算结果的证实.图12a和图12b分别给出了采热1年时模型MR和MT1对应的x-y平面内达西速度大小分布.假设MR和MT1分别代表相距500 m的上、下游两个平面，且两平面之间相应于水平井温度变化50 K的孔隙压力变化可一级近似为沿水平井走向的线性变化，则可以通过模型MR和MT1对应的孔隙压力场之差估算水平井走向上的压力梯度，进而根据方程(8)计算达西速度分量w.图12c给出了这样得到的w大小随x和y坐标的变化，而图12d给出了在中心线x=0上w沿储层高度y的变化及其与MR、MT1垂向达西速度结果的对比.图12c和图12d清楚地表明，水平井走向的渗流主要集中于半径30 m的渗透率强化区内部，且呈现强化区下部向下游流动而强化区上部向上游流动的特征.对比图12给出的水平井走向达西速度w与x-y面内速度分布，可以看到平行水平井走向的渗流比x-y面内渗流弱得多，达西速度的峰值相差三个量级，其在渗透率强化区内的平均值相差两个量级，而在强化区外也有量级上的差别.因此，定性分析和定量对比一致表明，本文采用的二维模型是对实际三维问题的合理近似.

图12 不同方向达西速度分量(m·s-1)的对比

本文模拟中设定储层四周为绝热边界条件，忽略了围岩向储层的热传导效应.这样的边界条件设定被数值模拟地热能开采的前人研究普遍采用(如，Bataillé et al., 2006)，而且储层与围岩耦合传热的数值模拟研究(丁军锋和王世民, 2019)也直接证实围岩向储层的热传导对地热开采贡献甚微.本研究采用的储层初始温度、水平井温度以及储层渗透率值，是在综合前人同类研究采用的参数取值基础上设定的，以便于进行结果对比.需要说明的是，储层与采热井之间的温差决定了地热开采热流，而同步提高或降低储层和采热井的温度对模拟结果没有任何影响.本文模型采用的储层基础渗透率以及在水平井周围渗透率的强化程度，均与前人研究量级相当.例如，Oldenburg等(2016)设定高于基础渗透率100倍的均匀分布强化渗透率为10-10m2.这与本文设定渗透率从强化区边界到水平井井壁线性递增达200倍的平均效果接近.虽然本文考虑的渗透率强化在实际工程中需要通过压裂等储层改造手段实现，但并不需要对压裂裂缝内的流动进行模拟.原因在于由遍布储层且间隔不大的大量微裂隙形成的流体运移通道网才能起到强化对流采热的作用，而少量大尺度裂缝的存在将为流体提供流动“捷径”(优势通道)，造成“短路”效应，反而不利于流体遍历储层采热.本文考虑压裂微裂隙网络强化渗透率的情形，基于达西定律构建渗流模型，无需考虑微裂隙的具体形态和裂隙内流动.至于实际存在大尺度裂缝的特定地热储层，仅由基于达西定律的渗流模型不足以精确描述大尺度裂缝的作用，必须结合大裂缝的具体形态构建裂缝内流动与缝间渗流的耦合模型来研究.这样的耦合模型构建和数值研究在理论探索和实际应用上都有重要意义，是一个值得深入研究的课题.