化石燃料需求增加导致供需压力逐年上升，环境污染也日益严重，可再生能源的开发利用是应对上述问题的有效途径^[1]。2019年，可再生能源消费量增长3.2×10¹⁸ J，占全球一次能源增长的40%以上，其在能源结构中占比由2018年的4.5%增至5.0%^[2]。然而，作为可再生能源的重要组成部分，地热能消费增长量不足1.0×10¹⁷ J，原因在于高昂的建井成本。钻完井成本占水热型地热系统取热总成本的30%～60%^[3]，高投入限制了地热井大范围推广利用。

2023年，国家能源局发布了《加快油气勘探开发与新能源融合发展行动方案（2023—2025年）》，建议将地热能等新能源设施作为油气开发项目的产能建设配套。在利好政策引导下，利用废弃井开发油区地热资源具备较好的发展前景。油田地区通常具有较高的地温梯度^[4]，且目前全球有2 000万至3 000万口废弃油井，早在2005年我国就存在16万口油气井，其中近8万口油气井因各种原因被废弃，约占总数的32%^[5]。靠近城镇油田所产热量可用于城市供暖、融雪和农作物养殖等，如胜利油田、华北油田和冀东油田等；远离城镇油田所产热量可用于油田伴热开发和油田建筑供暖等。

为了获取废弃井地热资源，学者们进行了大量研究。T. Kujawa等人^[6]使用双管换热器从废弃油井中提取热量，是较早利用废弃井进行地热能开发的研究。R. Caulk等人^[7]利用数模方法，分析了在废弃井中安装井下换热器的可行性。宋先知等人^[8–9]在河北雄县开展了单井取热现场试验，研制了高导热水泥和保温内管，形成了较为完备的单井开发利用理论体系。地热闭式系统改造和维护成本低，但取热效率相对较低；井网系统控制储层的体积大，加之油田地面配套设施齐全，开发油区地热具备较大优势^[10]。中国学者Gong Bin等人 ^[11]和美国学者D.Westphal等人^[12]已开展了固定井型的井网现场取热性能评估，证明了井网取热的可行性。

针对油田多相流研究，L.Nghiem等人^[13]以油藏水平井为对象开展了油水两相、气油两相数值建模，分析了油相产量和水油比等参数。朱广海等人^[14]建立了稠油热采过程中井筒内气液两相流温度场的计算方法，井口温度模拟结果与实测值的相对误差仅3.10%。蔡少斌等人^[15]考虑热流固耦合作用，开展了多孔介质孔隙尺度的两相流动模拟。

综合分析可知，尽管目前多相流研究已较为成熟，且有涉及温度场演化的文献报道，但针对大型油田井网多相流热流耦合的研究较少，且未见基于多相流热流耦合井网生产性能对比的研究。为此，笔者建立了考虑油水两相流的井网热流取热模型，并利用有限元软件进行求解，对比了不同井型下地热系统的热提取性能。该研究结果不仅适用于废弃井，同样适用于高/特高含水油藏取热，对油田地热开发利用具有一定的指导作用。

1.   模型描述

1.1   取热方法

以油田废弃井网为对象，进行地热井取热改造。考虑到油井长期生产过程中井底易堵塞，井筒内易结垢，故首先对井筒进行除垢和射孔改造。然后，将保温管、封隔器和电潜泵下到生产井中，并根据地下水位确定其下深。整个井段采用保温管进行绝热处理，以减少热流体在上返过程中的热量损失。最后，循环工作流体自注入井中泵入，与储层进行充分换热后从生产井中回收使用，后再次泵入储层完成循环。废弃井网热提取系统如图1所示。需要说明的是，本文只研究储层中的流动换热过程，井筒内流动和传热过程的研究后续进行。

图  1  废弃井网热提取系统示意

Figure  1.  Heat extraction system of abandoned well pattern

下载: 全尺寸图片 幻灯片

1.2   控制方程

水为热量提取的主要循环工质，其物性随温度发生变化^[16]：

式中：${\rho }_{{\mathrm{w}}}$为水相密度，kg/m³；$T$为温度，K；${c}_{p,{\mathrm{w}}}$为水相比热容，J/（kg·K）；${\lambda }_{{\mathrm{w}}}$为水相热导率，W/（m·K）；${\mu }_{{\mathrm{w}}}$为水相动力黏度，Pa·s。

油相动力黏度随温度发生变化：温度为283.15 K时，动力黏度为0.220 Pa·s；温度为293.15 K时，动力黏度为0.150 Pa·s；温度为303.15 K时，动力黏度为0.091 Pa·s；温度为313.15 K时，动力黏度为0.061 Pa·s；温度为323.15 K时，动力黏度为0.041 Pa·s；温度为333.15 K时，动力黏度为0.029 Pa·s；温度为343.15 K时，动力黏度为0.019 Pa·s；温度为353.15 K时，动力黏度为0.011 Pa·s；温度为363.15 K时，动力黏度为0.006 Pa·s；温度为373.15 K时，动力黏度为0.004 Pa·s。其他物性参数均设置为定值：密度为870 kg/m³，比热容为2 000 J/（kg·K），热导率为0.2 W/（m·K）。以上相关取值来源于多个油田数据和已发表文献^[17]。为便于开展研究，进行以下合理性假设：1）岩石基质水平渗透率为垂直渗透率的10倍^[18]，岩石物理性质恒定；2）模型采用局部热平衡假设^[19]；3）地层压力始终高于20 MPa，故储层中仅存在油水两相，不存在气相^[20]；4）油相和水相均为不可压缩流体。

储层中油水两相采用质量守恒和扩展达西定律进行求解：

式中：$t$为时间，s；$\phi$为孔隙度；${\rho }_{i}$为i相密度，kg/m³；${S}_{i}$为i相饱和度；${{{\boldsymbol{u}}}}_{i}$为i相速度，m/s；${Q}_{i}$为i相质量流量，kg/（m³·s）；${K}_{{\mathrm{r}}i}$为i相相对渗透率；${\mu }_{i}$为i相动力黏度，Pa·s；$K$为多孔介质渗透率，m²；${p}_{i}$为i相压力，Pa；g为重力加速度，m/s²。

水相和油相体积分数之和为1。将水相压力设置为润湿相，表示为p_w。油相压力p_o由毛细管压力进行计算，表达式如下：

式中：p_o为油相压力，Pa；p_c为毛细压力，Pa；p_w为水相压力，Pa；${S}_{{\mathrm{w}}}$和${S}_{{\mathrm{o}}}$分别为水相和油相饱和度。

混合流体物性由水相和油相物性加权得到，在此不再赘述^[21]。

毛细管力的计算采用了布鲁克斯−科里模型^[22]，其结果取决于饱和度：

式中：${p}_{{\mathrm{ec}}}$为入口毛细管压力，Pa；${f}_{{\mathrm{p}}}$为孔隙分布指数；${S}_{{\mathrm{rw}}}$为水相残余饱和度；${S}_{{\mathrm{ro}}}$为油相残余饱和度；${K}_{{\mathrm{rw}}}$和${K}_{{\mathrm{ro}}}$分别为水相和油相的相对渗透率。

传热方程表达式如下^[23]：

式中：${\left({\rho c}_{p}\right)}_{{\mathrm{eff}}}$为恒压下的有效体积热容，J/（m³·K）；${\rho }_{{\mathrm{f}}}$为流体密度，kg/m³；${c}_{p,{\mathrm{f}}}$为恒压下流体比热容，J/（kg·K）；${{\boldsymbol{u}}}$为速度场，m/s；${\lambda }_{{\mathrm{eff}}}$为有效导热系数，W/（m·K）；$Q$为热源，W/m³；${\rho }_{{\mathrm{s}}}$为恒压下岩石基质密度，kg/m³；${c}_{p,{\mathrm{s}}}$为恒压下岩石基质比热容，J/（kg·K）；${\lambda }_{{\mathrm{f}}}$为流体热导率，W/（m·K）；${\lambda }_{{\mathrm{s}}}$为岩石基质热导率，W/（m·K）。

1.3   几何模型

计算域长、宽、高分别为3 000，3 000和1 000 m，储层位于计算域中心位置，为边长500 m的立方体，如图2（a）所示。上述尺寸设计，基于边界效应分析和计算效率而确定。储层分为未改造段和改造段，后者指利用射孔、压裂等工艺进行处理。井筒以一维直线表示，生产井和注入井的长度均为300 m，对应储层的改造段，井筒直径均为0.152 4 m。以基本算例为例，蓝线段和红线段分别代表注入井和生产井，注采井水平间距分别为176.8，176.8和250.0 m，如图2（b）所示。

图  2  模型计算域

Figure  2.  Model computation zone

下载: 全尺寸图片 幻灯片

储层渗透率远高于围岩。储层和盖层物性参数见表1，数据参考了相关文献^[24–25]。

表  1  储层和盖层物性参数

Table  1.  Physical properties of reservoir and cap rocks

类型	密度/ （kg·m−3）	比热容/ （J·（kg·K）−1）	热导率/ （W·（m·K）−1）	水平渗透率/ mD	孔隙度， %
储层	2 700	800	2.0	500.0	30.0
围岩	2 800	850	2.2	0.1	5.0

下载: 导出CSV 

| 显示表格

1.4   初始条件及边界条件

模型采用质量流量入口和压力出口，注入流体温度为303.15 K，注入流量为60 kg/s（为所有注入井的总流量），平均生产压力设定为23.0 MPa（等于储层中心点处原始压力）。储层初始含油饱和度、束缚水饱和度和残余油饱和度分别为0.30，0.22和0.20^[13]。模型顶部和底部边界设置为无流动和绝热边界，侧边界设置为恒压和恒温边界，用以模拟储层存在的流体和热量补充。模型顶部初始压力为20.00 MPa，压力梯度约为6 000 Pa/m，模型顶部初始温度为333.15 K，温度梯度为0.035 K/m^[26]。

1.5   网格剖分及计算

计算域的网格划分方案如图3所示。首先，在储层顶面生成三角形单元，并细化注入井和生产井周围网格。然后，沿垂直方向（z方向）扫掠网格至储层底面，以生成三角形棱柱体单元，剩余几何区域采用自由四面体网格进行划分。通过网格无关性分析，决定采用12.94万个网格和43.67万个自由度的划分方案进行计算，以确保具有较高的计算效率和求解精度。

图  3  网格划分方案示意

Figure  3.  Mesh division scheme

下载: 全尺寸图片 幻灯片

对于单相热流耦合已经进行充分证明，本文主要对油水两相渗流进行验证。现有文献使用不同数值软件（不包括COMSOL）来计算相同的原油生产问题（油水两相流动）^[13]，利用本模型对文献中2个算例进行了计算，结果如图4所示。产油量最大误差分别为19.8 和22.4 t/d，原因是不同软件的计算和求解方法不同，生产后期吻合度较好。此外，利用本文模型计算得到的生产水油比位于参考文献中生产水油比的上下限之间，再次证明了本文多相流计算的可靠性。

图  4  本文模型结果与现有文献结果对比^{[13, 19]}

Figure  4.  Comparison of results of proposed model with those of existing literature

下载: 全尺寸图片 幻灯片

模型研究周期为3 660 d，时间步长设为1 d。所有变量的相对容差，包括压力、温度、各相饱和度、因变量和状态变量等，均设为10⁻⁵，并视为收敛标准，绝对容差设置为0.1。

2.   数模结果及讨论

2.1   压力场

基于第2节所述基础算例，对压力场进行了分析，结果如图5所示。储层附近，特别是注入井附近存在明显高压，而生产井附近压力保持在23.0 MPa（恒压生产）。流体自注入井沿径向流动至储层其他位置，波及区域明显。此外，边界处压力几乎无变化，表明所选几何模型尺寸合理，可减少边界效应对结果的影响。

图  5  生产3 660 d后不同截面压力分布

Figure  5.  Pressure distribution at different sections after production of 3 660 days

下载: 全尺寸图片 幻灯片

为定量分析，在图5中选取AC和BD两条直线作为观测对象，所得压力分布如图6所示。图6中，E点和F点、G点和H点分别为BD线和AC线与储层边界的交点。由图6（a）可知，注入井和生产井间的压力梯度远大于注入井与其他区域的对应值。由此可推断出，注入井注入流体首先充满储层，导致储层压力升高，因盖层渗透率远低于储层，流体主要流向生产井，进而形成较大压力梯度。当流量达到稳定后，流体将逐渐流向渗透低的盖层。渗透率的巨大差异将导致储层和盖层交界处的压力急剧变化，此现象在靠近注入井一侧更为明显。流入盖层流体较少，对其压力分布影响不大，故压力梯度较小，盖层中压力变化相对平缓。在储层中，注入井a与生产井b的间距为176.8 m，而注入井b与生产井d的间距为250.0 m，故a井注入流体相比于b井更易流向生产井，所需生产压差较小，如图6（b）所示。

图  6  AC和BD线的压力分布

Figure  6.  Pressure distribution at Line AC and Line BD

下载: 全尺寸图片 幻灯片

2.2   温度场

与压力波及区域相比，温度场影响范围较小。开采1 860和3 660 d后储层温度场的分布如图7所示。从注入井到生产井，呈现明显的指进现象，生产后期发生热突破，即低温区开始波及到生产井处，生产温度从初始地层温度开始下降。就不同截面而言，a井（c井）的热突破时间早于b井，主要由前者井间距较小导致。此外，开采1 860 d后，低温区尚未扩散到生产井，3 660 d后生产井附近出现了明显的低温，主要是由a井和c井的来流造成的。从图7（d）和图7（e）可以得出，a井（c井）和d井之间的裂缝已完全被低温区所覆盖，而b井和d井之间尚未实现完全的热突破，如图7（f）所示。

图  7  开采1 860和3 660 d后不同截面的温度分布

Figure  7.  Temperature distribution at different sections after production of 1 860 and 3 660 days

下载: 全尺寸图片 幻灯片

开采360，2 960和3 660 d后b井和d井间的温度分布如图8所示。通过分析b井和d井间的温度分布，可对上述现象进行定量解释。

图  8  开采360，2 960和3 660 d后b井和d井间的温度分布

Figure  8.  Temperature distribution between Well b and Well d after production of 360, 2 960, and 3 660 days

下载: 全尺寸图片 幻灯片

从图8可以得出，注入井附近的温度已近乎等于注入水的温度，表明该区域已完全冷却。早期，生产井附近的初始储层温度保持不变，经过长时间开采，生产井附近温度同样开始下降。值得注意的是，以3 660 d的温度曲线为例，曲线在生产井处显著下降，而生产井附近的温度高于生产井的温度。说明此时b井和d井之间尚未实现完全热突破，生产井温度降是由注入井a和c的来水造成的，即注入井c和生产井d及注入井c和生产井间已发生完全热突破，导致注入水直接进入生产井。

2.3   油水两相分布

为分析不同初始含油饱和度（S₂）下井网取热系统的流动和传热特性差异，研究了生产温度、注采压差和生产质量流量随时间的变化规律，结果见图9。从图9可以看出，当初始含油饱和度为0.25时，算例的生产温度和注采压差均小于初始含油饱和度为0.30时。其中，生产温度差异不明显，因油相黏度远高于水相黏度，注采压差存在显著差异。上述规律对于纯水储层和初始含油饱和度为0.30的储层更为明显：开采3 660 d后注采压差和生产温度对应差分别为0.291 MPa和2.19 K。由于注入水部分保留在储层当中，导致失水现象发生，故各情况下生产质量流量均小于注入质量流量（60 kg/s）。黏度低的流体更容易流动，因此纯水储层中水损失较大，开采10年后初始含油饱和度为0.30的储层和纯水储层最大差值达0.19 kg/s。随着开采进行，由于不断注入水以补充地层压力，水损失程度逐渐降低。

图  9  不同初始含油饱和度下的生产特征

Figure  9.  Comparison of production characteristics under different initial oil saturation

下载: 全尺寸图片 幻灯片

为分析储层基质中油相饱和度变化，选取图6（b）中b井左侧5 m处的点为观测点，可知油相饱和度逐渐降低，水相饱和度逐渐增加，但变化幅度减小（见图10）。开采3660 d后油相饱和度由0.300变为0.251，水相饱和度则由0.700变为0.749，油水饱和度相加始终为1。

图  10  相饱和度变化曲线

Figure  10.  Phase saturation curve

下载: 全尺寸图片 幻灯片

2.4   井网对比

油田井网类型复杂，为对比不同井型组合、注采比下的取热性能差异，本文设置8个算例，对应8种不同井网类型，如图11所示。算例1是对照案例，不同井网的总注入流量和注入温度均相同（见前文）。不同井网注入井与生产井间距依次为：算例1注采井间距为176.8和250.0 m，算例2注采井间距为176.8 m，算例3注采井间距为176.8和250.0 m，算例4和算例5注采井间距为125.0 m，算例6注采井间距为125.0和176.8 m，算例7和算例8注采井间距为125.0 m。井网类型的选取，均基于现场实际井网或已有文献：通过调整油田区域井网注采比，可获得与算例1—算例6对应的井网，算例7和算例8基于已有文献进行设计^[27]。

图  11  不同井网类型示意

Figure  11.  Different well pattern types

下载: 全尺寸图片 幻灯片

以储层中心水平截面为研究平面，不同井网类型温度场和压力场分布如图12所示。由图12可知，在相同条件下，不同井网类型温压场分布存在明显差异，主要体现在低压−高压边界差异，低温区指进程度不同。上述不同是由井网控制体积不同造成的，为此，以生产温度、注采压差、生产流量和日产油量为指标，对不同井网的生产性能进行了对比，结果见图13—图16。

图  12  不同井网类型下温压场的对比

Figure  12.  Comparison of temperature and pressure fields in different well pattern types

下载: 全尺寸图片 幻灯片

图  13  不同井网类型下生产温度的对比

Figure  13.  Comparison of production temperatures in different well pattern types

下载: 全尺寸图片 幻灯片

图13为不同井网类型下生产温度随时间变化情况的对比。开采3 660 d后，不同井网类型下生产温度差异明显，最大差值达到22.56 K，其中算例2的生产温度达到344.72 K，为各算例最高值。算例4和算例5的生产温度相对较低，均在322 K左右。从生产温度变化趋势可以得出，整个生产过程中，算例5的生产温度低于算例4。

图14为不同井网类型下注采压差随时间的变化情况的对比。由于本文所设置的基质渗透率较大，算例1—6的注采压差一般在0.30～0.60 MPa，算例7和算例8的生产压差值可达1.0 MPa以上。此外，通过分析发现，含油储层注采压差的演变与纯水储层存在差异，并非持续呈上升趋势，而是存在短期下降阶段（算例2—4，算例6），甚至在研究期持续下降（算例7）。随着开采进行，纯水储层温度降低将导致水相黏度增加，进而增大注采压差。在含油储层中，影响注采压差的因素不仅包括温度降低导致流体黏度升高，还包括油相饱和度变化。在相同温度下，油相黏度远高于水相黏度。随着生产的进行，参与流动油相的饱和度持续降低，导致油水混合物黏度降低。注采压差曲线上升抑或是下降，关键在于影响混合物黏度2个因素的主次关系。在上述几个算例中，注采压差降低的原因，可认为是油相饱和度降低所造成的油水混合物黏度降低效应强于受温度影响油水混合物的黏度升高效应。

图  14  不同井网类型下注采压差的对比

Figure  14.  Comparison of pressure differences between injection and production in different well pattern types

下载: 全尺寸图片 幻灯片

图15和图16分别为不同井网类型下生产质量流量和日产油量随时间的变化情况的对比。经过10年的开采，不同井网类型下生产质量流量间最大差值为0.45 kg/s。除算例4和算例8外，其他算例的生产质量流量约为59.60 kg/s。此外，由于失水，所有情况下的生产质量流量均小于注入质量流量。一般来说，注入水损失越小，提取效果越好。日产油量的变化与生产温度相似，10年后日产油量的区间为1.29 ～2.15 t，最大差值达到0.87 t。

图  15  不同井网类型下生产质量流量的对比

Figure  15.  Comparison of production mass flows in different well pattern types

下载: 全尺寸图片 幻灯片

图  16  不同井网类型下日产油量的对比

Figure  16.  Comparison of daily oil production in different well pattern types

下载: 全尺寸图片 幻灯片

综合分析生产温度、注采压差、生产质量流量和日产油量，算例2的生产温度和日产油量最高，其生产压差相对较低。因此，在本文研究条件下，算例2可被视为8种井网类型中的最优井型。

3.   结　论

1）温度场影响范围小于压力场，生产压差随初始含油饱和度增加而增加，生产温度随初始含油饱和度变化较小，油相在含油储层取热研究中不可忽视。

2）与纯水储层不同，含油储层中的流体流动和传热受温度、油相饱和度的双重影响；在生产早期，温度降低占主导地位，而在生产后期，油相饱和度降低起主导作用，前者使流体黏度升高，后者使流体黏度降低，将对注采压差产生显著的影响。

3）与算例2对应的井网类型（生产井和注水井交叉对称分布）具有最高的生产温度、日产油量和较低的生产压差，在所研究算例，热采效果最好。不同井网类型具有不同的注采比和井距，因此生产特征差异明显，证明了开展废弃油井转地热井取热性能评价及井型对比的必要性。