0 前言
赋存于地下3~10 km范围内的干热岩(hot dry rock,HDR)地热能具有巨大的储量和广泛的分布[1]，增强型地热系统(enhanced geothermal systems,EGS)旨在将HDR地热能提取至地上并加以利用。EGS的原理是通过水力激发等手段，在低渗透性的HDR内产生具有一定连通性的裂隙网络，形成人工热储；然后由注入井灌注冷流体工质，流体工质流过地下裂隙时获取HDR热量，热流体经由生产井开采出来后用于地面发电，发电后的流体工质经进一步梯级利用降温后再回注至地下热储，从而形成循环生产。地下采热过程是EGS的关键，直接影响EGS的产能和寿命。采热过程机理复杂，包含了热(T)、水力(H)、固体力学(M)和化学(C)等过程的综合作用，采热工质在岩体裂隙中的流动与换热则是EGS采热的核心作用过程，通常采用数值模拟方法研究和揭示相关的多场耦合机制，预测评估采热性能，为EGS的合理设计和优化提供基础条件和技术支持。随着模型的逐渐发展完善以及地质和其他相关基础数据的日渐丰富，数值模拟结果的准确性也将越来越高。Yang等[2]以及Gelet等[3]采用了双孔隙率的方法来区分岩石和流道的物性。为了更精确地模拟热储层的各向异性效果，Kalinina等[4]采用了多孔隙率的多孔介质模型。Jiang等[5]采用局部非热平衡热模型，实现了对热储层中岩石和流体之间局部换热过程的模拟。Shaik等[6]采用边界元和有限元耦合的方法计算裂隙内流动传热。Pruess等[7, 8]首次采用数值方法研究了以CO2为循环工质的EGS采热过程，比较了以CO2和水为工质的EGS采热过程，认为CO2 EGS具有较大的采热优势。
由于循环工质通常是以低温注入热储内，通过裂隙网络渗流并与高温岩体进行热交换获得热能，同时其动量产生损耗，因此在地下热开采过程中循环工质的温度和压力经历较大范围变化。温度和压力的变化将改变循环工质的热物性，从而影响采热过程中工质流体的输运和岩石-流体换热，也会影响EGS的寿命、出力等生产指标。本文在前期EGS采热过程数值模型工作基础上[5, 9, 10, 11, 12]考虑工质流体(水及二氧化碳)的变物性，并对EGS采热过程进行模拟，旨在探究工质的物性变化对EGS采热性能的影响。
1 EGS地下热开采过程建模
1.1 主控方程
EGS的地下部分包括井通道、热储和热储周围岩体3个区域。模型将EGS的地下部分处理为3个性质不同的子区域：①开放流道性质的注入井和生产井; ②多孔介质性质的热储; ③渗透性可忽略不计的热储周围岩体。模型假设单相流体流动，不考虑循环流体与岩石的化学、力学作用，主控方程[5]如下：
循环工质的连续性方程为
循环工质的动量守恒方程为
岩石区能量守恒方程为
流体工质内能量守恒方程为
其中： u 、p、T为待求解量，分别表示速度矢量、压力、温度；ρ、cp、λ、μ 为热物性参数，分别表示密度、比定压热容、导热系数和黏度系数，流体工质的热物性参数受工质的温度和压力决定，将通过后续变物性模型给出其与温度和压力的关系；下标f和s分别表示参数属于流体工质或岩石，导热系数λ的上标eff表示计算中采用有效导热系数，这里采用修正因子为1.5的Bruggeman关系式来确定，即λeffs=λs(1-ε)1.5和λefff=λfε1.5； ε和K分别为热储的孔隙率和渗透率；t为时间； g 为重力加速度；ha表示岩石-流体对流换热基础条件参数，本文取ha =1.0 W/(m3·K)。模型基于局部非热平衡思想，采用两个能量方程分别描述热储内流体和岩石骨架的温度场，可方便地处理采热过程中实际存在的岩石-流体换热过程。
1.2 水的变物性模型
笔者采用国际水及蒸汽物性组织(IAPWS)[13]的公式建立水工质的变物性模型。IAPWS模型针对水不同的相建立了相应的计算公式，考虑到EGS系统中水的温度和压力范围(压力低于100 MPa，温度一般低于623.15 K)，仅选用液相区域进行建模。该区域采用Gibbs自由能描述水的状态：
其中：G(p,T)为Gibbs自由能，单位为J/kg；Rw为水的比气体常数，Rw =461.526 J/(kg·K)；γ(π,τ)、π与τ分别为无量纲化的Gibbs自由能、压力和温度，π=p/p*，τ=T*/T，参考压力 p*=16.53 MPa，参考温度T*= 1 386 K；ni，Ii，Ji均为常数，可参考文献[13]。利用Gibbs自由能的偏导数，就可得到随温度和压力变化的密度及比定压热容 其中：γπ表示γ(π,τ)关于π的一阶偏导数；γττ表示γ(π,τ)关于τ的二阶偏导数。 式中：参考值λ*=1.0×10-3 W/(m·K)，μ*=1.0×10-6 Pa·s；
及δ分别为无量纲化的温度和密度，=T/T*，δ＝ρ/ρ*，T*=647.096 K，参考密度ρ*=322.0 kg/m3；
分别对应稀释气体项(仅与温度相关)、有限密度项(与温度和密度相关)和临界点修正项，详细模型可参看文献[14, 15]。由于热储内的温度没有达到水的临界点，因此本文计算中
。
1.3 超临界CO2变物性模型
在EGS开采过程中，二氧化碳在热储高温高压条件下处于超临界状态[16]。笔者采用改进的Redlich-Kwong状态方程[17]进行超临界二氧化碳(supercritical carbon dioxide，SCCO2)的密度计算，该模型可由压缩因子的多项式表达：
其中：Z为压缩因子；参数A和B的计算式可参看文献[17]。SCCO2的比定压热容采用Helmholtz能量进行建立，详细模型可参考文献[18]。Helmholtz能量公式为 ρ*=467.6 kg/m3，T*=304.13 K。等式右端由理想气体项φo(δ,τ)和残余项φr(δ,τ)组成，比定压热容则可由Helmholtz能量的偏导数表述： 式中，φo及φr的下角标τ和δ分别表示该函数关于τ或δ的一阶及二阶偏导数。 其中，C0——C8以及D0——D8均为常数参数，取值列于表 1中。
| 参数 | 取值 | 参数 | 取值 |
| C0 | 1.492 882×101 | D0 | -1.146 067×10-1 |
| C1 | 2.625 411×10-3 | D1 | 6.978 380×10-7 |
| C2 | 8.778 046×10-6 | D2 | 3.976 765×10-10 |
| C3 | -5.114 246 | D3 | 6.336 120×10-2 |
| C4 | 4.377 109×10-1 | D4 | -1.166 119×10-2 |
| C5 | 2.114 051×10-5 | D5 | 7.142 596×10-4 |
| C6 | -4.730 357×10-1 | D6 | 6.519 333×10-6 |
| C7 | 7.366 358×10-2 | D7 | -3.567 559×10-1 |
| C8 | -3.763 399×10-3 | D8 | 3.180 473×10-2 |
2 算例及分析
2.1 计算模型及参数
图 1为假设的5口井(1注入井、4生产井)布局的EGS，由于对称性取1/4进行计算。人工热储为500 m×500 m×500 m的立方体，注入井和生产井均为0.2 m×0.2 m的方形通道。热储周围包覆足够体积的不可渗透岩石，在所考察的EGS运行期内，热储内岩石温度降低不会波及计算域边界，避免了人为设定边界条件带来的误差。计算域尺寸如图 1a所示。岩石的初始温度按照4 K/hm的地温梯度随深度线性增加，地表温度设为300 K，初始温度分布如图 1b所示。计算域最外侧均采用绝热边界；人工热储与周围岩石间为计算域内部界面，热流和温度在界面处连续。本文的热边界条件对计算结果的影响在文献[10]中已讨论。热储中裂隙流体与当地岩石温度相同，所有与流体接触的壁面均为非滑移边界，均采用定压力边界，注入井与生产井压差取10 MPa，热储内环境压力取40 MPa。注入温度为343 K，热储孔隙率为0.01，渗透率为1.0×10-14 m2。
 |
| 图 1 模型几何尺寸及初始温度分布 Fig. 1 Geometry dimensions and initial temperature of the quintuplet EGS |
| |
2.2 水工质变物性对EGS采热的影响
为了研究热物性变化对EGS采热性能的影响，笔者分别计算了以水为工质的常物性条件和变物性条件下的热开采过程。由于密度包含在所有控制方程中，是实现流动传热双向耦合的关键，同时也是求解变黏度和变导热系数的基础，因此本文在变物性模拟算例中均包含变密度条件，计算方案如表 2所示。常数物性算例1的取值为水在温度450 K、压力40 MPa状态下的值。
| ρ/(kg/m3) | cp/(J/(kg·K)) | λ/(W/(m·K)) | μ/(Pa·s) | |
| 算例1 | 913.8 | 4 254.2 | 0.70 | 1.63×10-4 |
| 算例2 | ρ(p,T) | 4 254.2 | 0.70 | 1.63×10-4 |
| 算例3 | ρ(p,T) | cp(p,T) | 0.70 | 1.63×10-4 |
| 算例4 | ρ(p,T) | 4 254.2 | λ(T,ρ) | 1.63×10-4 |
| 算例5 | ρ(p,T) | 4 254.2 | 0.70 | μ(T,ρ) |
| 算例6 | ρ(p,T) | cp(p,T) | λ(T,ρ) | μ(T,ρ) |
图 2显示了算例6中开采至第10年时热储内工质温度、压力以及各物性的三维等值面分布。密度、黏度系数与温度成反相关关系，比定压热容则与温度成正相关关系，而导热系数在考察温度范围内存在一个峰值，随着温度由高到低导热系数先增大后减小。在热储四周边缘区，储热量几乎未被开采，温度还大致遵从初始分布：沿深度方向存在约4 K/hm的温度梯度。在温度场分布的影响下，流体热物性参数在热储四周边缘区域同样表现出随深度方向的梯度分布，如图 2所示。计算物性相对于其在考察区间内最小值的相对变化率，则密度的相对变化率为10.00%，比定压热容为4.39%，导热系数为3.81%，黏度系数为150.0%；可见，随着开采的进行，黏度系数的变化极为显著，这将直接影响工质在热储内的流动阻力，进而影响生产井的采出质量和出力。
 |
| 图 2 算例6热储内温度压力及各物性参数分布(10 a后) Fig. 2 Temperature,pressure and property distributions in the heat reservoir of Case 6 after 10 a operation |
| |
图 3显示了算例1至算例6 EGS采出温度(production temperature)和工质的质量流率(mass flow rate)随时间变化的曲线。可以看出，热物性的变化对开采寿命、采出质量具有显著的影响。以采出温度降低10 K作为系统废止的条件，即将本文中生产井采出温度降至450 K的时刻作为该系统的开采寿命。由图 3a可见，与给定的常物性条件相比:引入变密度条件后开采寿命有所降低(约为9.0 a);在变密度基础上增加变比定压热容条件发现开采寿命进一步降低至7.5 a。这是因为相比常物性条件而言，由变密度模型计算的密度随着开采的进行逐渐增大，由变比定压热容模型获取的比定压热容值则在温度为450~460 K时高于恒定比定压热容条件;由能量守恒方程可知,较大的比定压热容和密度能够在单位时间带走更多的热量。在变密度基础上引入变导热系数，发现采出温度曲线无明显变化；表明导热系数对热开采过程影响较小，对流换热是热开采的主导因素。引入变黏度模型发现获得的开采寿命显著增大(约为18.0 a)。结合图 3b可知，由于后续注入热储的工质温度较低，黏度系数较大，在注入泵功恒定的条件下质量流率将逐渐降低，生产井出力逐渐减少，因此开采寿命有所延长。
 |
| 图 3 不同物性条件的采出温度曲线(a)和质量流率曲线(b) Fig. 3 EGS production temperature curves (a) and mass flow rate curves (b) for cases with different differing property schemes assumed |
| |
2.3 变物性条件下水及SCCO2采热性能对比
为了比较水及SCCO2的采热性能，针对SCCO2 EGS计算了表 2描述的算例6变物性条件算例，物性参数由公式(11)——(14)给出。计算模型、初始条件及边界条件均按照2.1节设定。
开采至第10年计算获得的各变量分布如图 4所示。可以看出，SCCO2的密度变化具有与水比拟的液相流体特征，而黏度系数则具有气相流体特征，远低于水的黏度。此外与水工质不同的是，在热储内的温度和压力变化范围内，SCCO2的比定压热容容与温度呈负相关关系。
 |
| 图 4 热储内SCCO2工质温度压力及各物性参数分布(10 a后) Fig. 4 Temperature,pressure and property distributions in the heat reservoir of SCCO2 after 10 a operation |
| |
为了比较热储由不同工质采出热量的情况，定义采热比(heat extraction ratio)为
式中：Tg为地表温度；Vh为热储体积；Vb为基岩体积；Vh+Vb构成了本文中图 1所示的计算域体积。等式右端分母为以地表温度为参考的热储内总热能，分子包含了由热储内与基岩内已开采出的能量，该参数的实质是t时刻由EGS系统采出的无量纲化的总能量。计算获得的水及SCCO2采出温度和热开采率对比如图 5所示。 可以看出，SCCO2的开采寿命约为12.5 a，低于水的开采寿命(16.3 a)；在各自废止的时刻，采用SCCO2的EGS热开采率约为0.44，采用水的约为0.45，略高于前者，这是因为以水为工质的开采寿命较长，周围岩石依靠热传导对已开采的低温区域进行热补偿，从而提高了热储的整体热开采率；而在相同开采时刻，SCCO2的热开采率则高于水。虽然SCCO2的密度和比定压热容低于水，但结合图 6质量流率对比可知注入泵功恒定条件下单位时间采出的SCCO2质量约为水的4倍，在密度、比定压热容和黏度物性的综合作用下采用SCCO2工质的采热速率高于水，相同开采时刻的热开采率就相应地更高。
 |
| 图 5 水及SCCO2 EGS采出温度和热开采率对比 Fig. 5 Comparison of production temperature and heat extraction ratio between water and SCCO2 EGSs |
| |
 |
| 图 6 水及SCCO2 EGS质量流率对比 Fig. 6 Comparison of mass flow rate between water and SCCO2 EGS |
| |
3 结语
笔者在基于局部非热平衡假设的EGS热开采过程三维计算模型基础上，引入了水和超临界二氧化碳的变物性模型，实现了求解过程的热流双向耦合。以水工质EGS为例分析了各物性参数变化对EGS采热性能的影响，并对变物性条件的水和超临界二氧化碳EGS采热进行了对比研究。计算结果显示，工质的密度和比定压热容越大则EGS开采寿命越短，黏度系数越大则EGS开采寿命越长，导热系数则对EGS采热性能无明显影响。注入压力一定的条件下以水为工质的EGS具有较长开采时间，相同时刻的质量流率和热开采率低于以二氧化碳EGS,但在运行终止时两者的热开采总量却大致相等。
本文的研究结果表明工质的热物性变化对EGS开采寿命、热开采率等具有重要影响，而热物性变化由工质的温度、压力状态所决定，因此本文研究结果可为工质的注入温度、压力参数设置、工质的选择方面提供理论依据。由于本文的计算模型采用了理想化的热储，仅从机理方面解释了工质变物性对EGS热开采的影响，在后续研究中有必要结合实际地质模型，研究实际EGS运行中工质物性的变化规律及其效应。
| [1] | Tester J W, Anderson B J, Batchelor A S, et al. The Future of Geothermal Energy: Impact of Enhanced Geothermal Systems (EGS) on the United States in the 21st Century[R]. Cambridge: Massachussets of Institute of Technology, 2006. |
| [2] | Yang Y, Yeh H. Modeling Heat Extraction From Hot Dry Rock in a Multi-Well System[J]. Applied Thermal Engineering, 2009, 29: 1676-1681. |
| [3] | Gelet R, Loret B, Khalili N. A Thermo-Hydro-Mechanical Coupled Model in Local Thermal Non-Equilibrium for Fractured HDR Reservoir with Double Porosity[J]. Journal of Geophysical Research, 2012, 117: 7205-7228. |
| [4] | Kalinina E, McKenna S, Hadgu T, et al. Analysis of the Effects of Heterogeneity on Heat Extraction in an EGS Represented with the Continuum Fracture Model[C]//Thirty-Seventh Workshop on Geothermal Reservoir Engineering, Stanford University. Stanford:CA, 2012:436-445. |
| [5] | Jiang F M, Luo L, Chen J L. A Novel Three-Dimensional Transient Model for Subsurface Heat Exchange in Enhanced Geothermal Systems[J]. International Communications in Heat and Mass Transfer, 2013, 41: 57-62. |
| [6] | Shaik A R, Rahman S S, Tran N H, et al. Numerical Simulation of Fluid-Rock Coupling Heat Transfer in Naturally Fractured Geothermal System[J]. Applied Thermal Engineering, 2011, 31: 1600-1606. |
| [7] | Pruess K. Enhanced Geothermal Systems (EGS) Using CO2 as Working Fluid:A Novel Approach for Generating Renewable Energy with Simultaneous Sequestration of Carbon[J]. Geothermics, 2006, 3: 51-67. |
| [8] | Pruess K. On Production Behavior of Enhanced Geothermal Systems with CO2 as Working Fluid[J]. Applied Thermal Engineering, 2008, 49: 1446-1454. |
| [9] | 陈继良, 蒋方明. 增强型地热系统热开采性能的数值模拟分析[J]. 可再生能源, 2013, 31(12): 111-117. Chen Jiliang, Jiang Fangming. A Numerical Study on Heat Extraction Performance of Enhanced Geothermal Systems[J]. Renewable Energy Resources, 2013, 31(12): 111-117. |
| [10] | 陈继良, 罗良, 蒋方明.热储周围岩石热补偿对增强型地热系统采热过程的影响[J]. 计算物理, 2013, 30(6): 862-870. Chen Jiliang, Luo Liang, Jiang Fangming. Thermal Compensation of Rocks Encircling Heat Reservoir in Heat Extraction of Enhanced Geothermal System[J]. Chinese Journal of Computational Physics, 2013, 30(6): 862-870. |
| [11] | 陈继良, 蒋方明, 罗良. 增强型地热系统地下渗流场的模拟和分析[J]. 计算物理, 2013, 30(6): 871-878. Chen Jiliang, Jiang Fangming, Luo Liang. Numerical Simulation of Down-Hole Seepage Flow in Enhanced Geothermal System[J]. Chinese Journal of Computational Physics, 2013, 30(6): 871-878. |
| [12] | Jiang F M, Chen J L, Huang W B, et al. A Three-Dimensional Transient Model for EGS Subsurface Thermo-Hydraulic Process[J]. Energy, 2014, 72: 300-310. |
| [13] | The International Association for the Properties of Water and Steam. Revised Release on the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam[R]. Lucerne:International Association for the Properties of Water and Steam, 2007. |
| [14] | The International Association for the Properties of Water and Steam. Lamda Release on the IAPWS Formulation 2011 for the Thermal Conductivity of Ordinary Water Substance[R]. Plzeň:International Association for the Properties of Water and Steam, 2011. |
| [15] | The International Association for the Properties of Water and Steam. Viscosity Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance[R]. Berlin:International Association for the Properties of Water and Steam, 2008. |
| [16] | 雷宏武, 李佳琦, 许天福, 等. 鄂尔多斯盆地深部咸水层二氧化碳地质储存热-水动力-力学(THM)耦合过程数值模拟[J]. 吉林大学学报:地球科学版, 2015, 45(2): 552-563. Lei Hongwu,Li Jiaqi,Xu Tianfu, et al. Numerical Simulation of Coupled Thermal-Hydrodynamic-Mechanical (THM) Processes for CO2 Geological Sequestration in Deep Saline Aquifers at Ordos Basin, China[J]. Journal of Jilin University: Earth Science Edition, 2015, 45(2): 552-563. |
| [17] | Heidaryan E, Jarrahian A. Modified Redlich-Kwong Equation of State for Supercritical Carbon Dioxide[J]. The Journal of Supercritical Fluids, 2013, 81: 92-98. |
| [18] | Span R, Wagner W. A New Equation of State for Carbon Dioxide Covering the Fluid Region from the Triple Point Temperature to 1 100 K at Pressures up to 800 MPa[J]. Journal of Physical and Chemical Reference Data, 1996, 25: 1509-1596. |
| [19] | Jarrahian A, Heidaryan E. A Novel Correlation Approach to Estimate Thermal Conductivity of Pure Carbon Dioxide in the Supercritical Region[J]. The Journal of Supercritical Fluids, 2012, 64:39-45. |
| [20] | Heidaryan E, Hatami T H, Rahimi M, et al. Viscosity of Pure Carbon Dioxide at Supercritical Region: Measurement and Correlation Approach[J]. The Journal of Supercritical Fluids, 2011, 56: 144-151. |