2021, Vol. 53引用本文 |
近年来，化石资源临近枯竭，环境问题日益突出[1]。干热岩资源因储量巨大、清洁环保，决定了其必然会成为未来重要能源之一[2]。增强型地热系统（EGS）作为开发干热岩资源的有效手段，受到了广泛关注[3-5]。
研究EGS取热过程时，若利用建造干热岩电站的方法获得实际工程的地质资料、工程经验和数据，对研究干热岩储层热提取极为有利，但目前的干热岩储层改造技术尚不完善，投资巨大、周期长。因此，现阶段研究大多基于数值模拟方法建立数学模型，利用多场耦合过程对EGS取热过程进行评估，指导EGS的设计与优化[6-8]。在热储改造和EGS运行过程中，热储岩体的物理化学性质都会随时间和空间不断变化[9-10]，改造后热储尺寸多以km为量级，内部裂隙通常以mm为量级，在利用精细化网格直接模拟取热过程时计算量巨大[11]，需要对热储内的裂隙进行简化。目前，常用的简化模型有单孔隙模型[12]、双孔隙模型[13-14]和多孔隙模型。其中，双孔隙模型和多孔隙模型因精度要求较高，需要大量试验数据支撑，在实际中应用困难[15]。
干热岩热量资源预测方面，目前主要有两种预测方式：一种是利用目标热储体积及热储温度和某一确定温度（如地表温度）的差值计算，按百分比的可开采量预测[16-19]，主要有体积法、蒙特卡罗法及地表流量法等[20-23]，但这些方法在应用于EGS产能预测时，热储开采前后温差的设定值偏大，导致预测不准确；另一种是数值法[24]，根据详细的热储参数计算和评价勘查程度较高、具有一定开采历史的地热田[23]，此种方法较为繁琐，计算周期长，需要详细的热储参数，不适合未开采的热储产能的预测。在EGS运行系统可控因素的研究方面，目前的研究多集中于以不同的热媒介质注入温度、不同的热媒介质注入速度及不同井间距等人为可控因素控制EGS的运行[8,25]。以注水速度控制时，随着EGS的运行，注入井的压力不断产生变化，大大增加了EGS运行时的不可控性[26]。EGS取热效果方面，大多数研究将热储视为单孔隙模型，针对热储改造后的孔隙比、渗透率、注入井温度、布井方式等多种因素分析EGS取热效果[27-28]。如：王昌龙[29]考虑了11种因素对EGS取热效果的影响。在这些研究中，大多局限于单因素的定性研究，较少考虑多因素组合的定量研究，且鲜有将大地热流的热补偿作为初始条件加入到模型之中[30-33]。有研究发现，热储周围的大地热流对热储层长期积累的热补偿热量会对EGS运行寿命产生不可忽略的影响[34]。
以吸放热公式[35]和Dupuit公式为基础，建立EGS产能和寿命控制方程；依据产能和热储参数，通过调节热媒介质注入温度和注采压差两个可控因素取值预测EGS寿命，增加EGS运行时的可控性。同时，根据地下水流动方程和热传导方程分析出热储初始温度、热媒介质注入温度、EGS寿命、热储体积和热储比表面积5个因素对发生热突破（生产井产出温度低于423.16 K[36]）时对控制方程中热储减少热量、大地热补偿热量、生产井平均产出温度和热储形状系数4个未知参数的影响。在考虑大地热补偿的情况下，利用双井EGS的3维数值模型，对不同工况的EGS在发生热突破时的4个未知参数的变化规律进行研究，得到4个未知参数随5个影响因素变化的具体预测公式，实现EGS产能预测和多因素定量化分析。
1 流热耦合下EGS换热过程分析如图1所示，假设热储周围围岩渗透率为0，热媒介质在流动过程中无热量损失，人工储层原有热量Qm和大地补偿的热量Qt完全由热媒介质携带并提取，其中，热媒介质携带热量由Qw表示。
 |
| 图1 EGS换热过程的流热耦合机制 Fig. 1 Flow-heat coupling mechanism of EGS heat transfer process |
根据热量平衡，由吸放热公式（式（1））结合Dupuit公式（式（2））可得热媒介质携带热量（EGS产出热量），如式（3）所示：
| ${{{Q}}_{\rm{w}}} = {c_{\rm{w}}}{{{m}}_{\rm{w}}}\Delta {T_{\rm{1}}}$ | (1) |
| ${{{v}}_{\rm{w}}} = k\frac{P}{{{\rho _{\rm{w}}}gL}}$ | (2) |
| ${{{Q}}_{\rm{w}}} = {c_{\rm{w}}}k\frac{P}{{gL}}{A_{\rm{1}}}t\Delta {T_{\rm{1}}}$ | (3) |
式（1）～（3）中：t为EGS寿命；k为热储渗透率；
令
若热储在tn时刻发生热突破，tm时刻的产出温度Tout(tm)可作为Tout(t)在t0至tn时刻的平均值，即：
| ${T_{{\rm{out}}}}({t_m}) = \frac{1}{{{t_n} - {t_0}}}\int_{{t_0}}^{{t_n}} {{T_{{\rm{out}}}}(t)} {\rm{d}}t$ | (4) |
因此，可以利用 EGS运行周期内的平均产出温度
整理式（3）和（4），可得EGS产能和寿命控制方程：
| $\frac{{{Q_{\rm{m}}} + {Q_{\rm{t}}}}}{t} = \alpha k({\overline T _{{\rm{out}}}} - {T_{{\rm{in}}}})P$ | (5) |
由式（5）可知，只需确定出热储在发生热突破时热储减少的热量Qm、大地热补偿热量Qt、热媒介质产出平均温度
利用达西定律和多孔介质传热方程描述流热耦合过程，如式（6）～（9）所示：
连续性方程：
| $\frac{{\partial {\rho _{\rm{w}}}}}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}({\rho _{\rm{w}}}{v_i}) = 0,\;\;\;i = 1,2,3$ | (6) |
流体动量方程：
| ${\rho _{\rm{w}}}\frac{{\partial {v_i}}}{{\partial t}} = {F_i} - \frac{{\partial {{P}}}}{{\partial {x_i}}} + \eta {\nabla ^{\rm{2}}}{v_i} = 0,\;\;i = 1,2,3$ | (7) |
流体能量方程：
| $ {\rho _{\rm{w}}}{c_{\rm{w}}}\left( {\frac{{\partial T}}{{\partial t}} + {{{v}}_{{i}}}\frac{{\partial T}}{{\partial {x_i}}}} \right) = {\lambda _{\rm{w}}}{\nabla ^2}T,\;\;\;{{i}} = 1,2,3 $ | (8) |
周围围岩热补偿能量方程：
| ${\rho _{\rm{m}}}{{{c}}_{\rm{m}}}\frac{{\partial T}}{{\partial t}} = {\lambda _{\rm{p}}}{\nabla ^2}T - h{A^{\rm{2}}}(T - {T_\infty })$ | (9) |
牛顿冷却公式：
| $\varPhi = {A_3}h({T_{\rm{s}}} - {T_{\rm{w}}})$ | (10) |
式（5）～（10）中，t为热媒介质流动时间，xi为热媒介质在i方向上的位移，ρw、ρm分别为流体、围岩的热储密度，vi为热储中流体的流速，Fi为体积力，
利用式（6）～（9），结合牛顿冷却公式（10），对式（5）中的4个参数进行分析，可确定出影响4个参数的主要因素，如表1所示。
| 表1 4个未知变量的影响因素汇总 Tab. 1 Summary of influencing factors of four unknown parameters |
 |
Qm由热储终温Tml和热储初始温度Tm控制，所以在预测时仅需考虑Tml即可，整理后如式（11）所示。
| $\left\{\!\!\!\! \begin{array}{l} {T_{{\rm{ml}}}} = {f_{\rm{1}}}({T_{\rm{m}}},{T_{{\rm{in}}}},V,{S\!_{{\rm{Sa}}}}{{,t}}), \\ {Q_{\rm{t}}} = {f_{\rm{2}}}({T_{\rm{m}}},{T_{{\rm{in}}}},V,{S\!_{{\rm{Sa}}}},t), \\ {\overline T _{{\rm{out}}}} = {f_{\rm{3}}}({T_{\rm{m}}},{T_{{\rm{in}}}},V,{S\!_{{\rm{Sa}}}}), \\ \alpha = {f_{\rm{4}}}(V,{S\!_{{\rm{Sa}}}}) \end{array} \right.$ | (11) |
若利用式（5）达到讨论Tin和P与t之间关系的目的，可利用数值模拟预测4个参数的具体表达式，并代入式（5）中得到Tin、P、t三者关系。
2 EGS取热模型概述 2.1 地质模型构建如图2所示，EGS地下部分由人工热储、热储周围的围岩、注入井和生产井组成。为了更准确地模拟EGS地下热交换过程，在模型中作出主要假设：1）热储裂隙小而均匀，可看做单孔隙率的多孔介质；2）热媒介质流动为单相流；3）热储中热媒介质流动是饱和层流流动（0.003 m/s以下），初始时刻热储层中充满了与周围岩石温度相同的热媒介质；4）周围围岩渗透率为0，热媒介质在流动过程中无损失；5）热储及周围围岩不会产生形变和各种化学反应。
 |
| 图2 EGS地下部分概念图 Fig. 2 Concept picture of EGS underground part |
为了模拟EGS的取热过程，探究单孔隙双井EGS的Tml、Qt、
 |
| 图3 双井EGS模型3维图 Fig. 3 Three-dimensional pictures of the double well EGS model |
如图3所示，热储周围包覆200 m厚的花岗岩围岩，在外表面施加以3 K/100 m随深度线性增加的温度边界，模拟大地热流给予热储的热量补给；分别建立不同边长的正六面体，研究体积对4个参数的影响，如图3（a）所示；分别建立同体积下不同形状的热储，研究比表面积对4个参数的影响，如图3（b）所示；其余模型边界条件如表2所示。模型采用“下注上采”的开采方式，最大模拟时间50 a，采用瞬态求解方式，计算结果每年采集一次，且最终结果取符合EGS运行25 a，生产井出水速率达到80 kg/s以上，保证EGS商业化利用[35]。
| 表2 模型边界条件取值 Tab. 2 Values of model boundary conditions |
 |
2.2 模型热物性参数取值
| 表3 模型热物性参数取值 Tab. 3 Values of model thermal property parameters |
 |
2.3 网格独立性验证
本文主要对计算结果中的生产井平均产出温度和平均速度、热储平均温度和热储中流体平均速度进行独立性检验。
进行网格无关性检验时，按时间步长0.01 s，对不同自由四面体的网格数量进行无关性检验。采用热储边长为500 m正六面体，热储初始温度–热媒介质注入温度–注采压差为565.66 K–353.16 K–20 MPa的工况，分别选取第15年和第30年的计算结果进行比较（图4）。由图4可知：网格数量对计算结果有着较大的影响，当网格数量小于25万个时，计算结果存在较大差异；当网格数量大于30万个时，模拟结果变化不大。
 |
| 图4 网格无关性检验结果 Fig. 4 Results of the grid independence verification |
进行时间独立性检验时，以第12种（36.76万个）网格数量为网格模型，时间步长为0.007 5、0.010 0、0.012 5和0.015 0 s，其余条件与网格独立性检验时相同，结果如表4所示。由表4可知，当网格数量达到30万以后，时间步长在0.012 5 s以下时对模拟结果影响不大。
| 表4 时间独立性检验结果 Tab. 4 Results of the time independence verification |
 |
由图4和表4可知，在计算时，最大网格取40 m，最小网格取0.3 m，网格增长率取2.4，单元格数量41.21万个，时间步长选择0.01 s较为合理。
3 结果分析经过计算，表1中Tm、Tin、V、SSa和t这5个因素对Tml、Qt、
 |
|
图5 Tm、Tin、V、SSa、t对Tml、Qt、 |
由图5（a）、（c）、（d）、（e）可知：当Tm和 Tin增加时，Tml逐渐上升，且Tm对Tml的影响要低于Tin对Tml的影响；随着SSa的增加，Tml 也会逐渐升高；随着V的增大，Tml表现出逐渐减小的趋势；以Tm=565.66 K、Tin=293.16 K为例，EGS不同寿命所对应的Tml相差不大，近似服从高斯分布，与时间项无关。由于热突破条件并非与热媒介质注入温度相同，所以在注入井与生产井之间必定存在温度梯度，在均匀介质渗流场中，两井之间必存在流速最快的最优通路。在最优通路间，温度梯度介于注入温度和热突破温度之间，在最优通路外其温度梯度上限与热储初始温度有关。因此，随着Tm的升高，最优通路外的热储温度升高；随着Tin的升高，整个热储温度梯度下限升高；随着SSa的增加，热储与围岩接触面积增加，热补偿的热量增多，都会导致Tml的升高。此外，其他条件相同时，V的增加使得热突破时低温区域增大，Tml逐渐减小。
由图5（b）、（c）、（d）可知：随着Tm和V的增加，
由图5（f）、（g）、（h） 可知：受围岩热补偿的影响，25年后 ，Qt随着EGS运行时间逐渐增加；随着Tm、V和SSa的增大，Qt也逐渐增加；对于Tin对Qt的影响，由图5（f）分析可以发现，除热储初始温度处于485.66 K时，Qt随热媒介质注入温度逐渐增加，其余工况均无明显规律。分别取第27年和第29年Qt的值，利用极差分析法将Tm与Tin对Qt的影响进行比较，如表5和6所示，发现热储初始温度Tm对Qt的影响要远远大于Tin对Qt的影响，可近似将Tin对Qt的影响忽略。因此，可以确定，Qt主要受到Tm、Tin、V和SSa的影响。
| 表5 第27年Tm与Tin极差分析 Tab. 5 Range analysis of Tm and Tin in the 27th year |
 |
| 表6 第29年Tm与Tin极差分析 Tab. 6 Range analysis of Tm and Tin in the 29th year |
 |
由图5（g）、（h）可知：随着V的增大，α表现出先减小后增大的趋势，其最低点大致出现在边长为400 m的正六面体热储附近；随着SSa的增大，α表现出先增大后减小的趋势，其最高点出现在正六面体热储附近。
经上述分析，将表1中的影响因素予以修正，如表7所示，同时将式（11）修正为式（12）：
|
表7 修正后的Tml、Qt、 |
 |
| $\left\{\!\!\!\! \begin{array}{l} {T_{{\rm{ml}}}} = {f_{\rm{1}}}({T_{\rm{m}}},{T_{{\rm{in}}}},V,{S\!_{{\rm{Sa}}}}), \\ {Q_{\rm{t}}} = {f_{\rm{2}}}({T_{\rm{m}}},V,{S\!_{{\rm{Sa}}}},t), \\ {\overline T _{{\rm{out}}}} = {f_{\rm{3}}}({T_{\rm{m}}},{T_{{\rm{in}}}},V,{S\!_{{\rm{Sa}}}}), \\ \alpha = {f_{\rm{4}}}(V,{S\!_{{\rm{Sa}}}}) \end{array} \right.$ | (12) |
为得到式（12）的具体表现形式，需要对上述计算结果继续分析。
3.1 热储终温Tml 变化规律为探究表7中4个因素与Tml的定量关系，以Tin为353.16 K时为基准，确定出Tml353.16随Tm的变化；然后，基于此分别讨论Tin、SSa与V对Tml的增量的变化，如图6所示；最后，将上述4项进行线性叠加，如式（13）所示：
 |
| 图6 热储终温Tml计算图 Fig. 6 Calculation of final storage temperature (Tml) |
| ${T_{{\rm{ml}}}} = {T_{{\rm{ml353}}{\rm{.16}}}} + \Delta {T_{{\rm{ml}}T{\rm{in}}}} + \Delta {T_{{\rm{ml}}S}} + \Delta {T_{{\rm{ml}}V}}$ | (13) |
如图6（a）所示，极易确定出Tml353.16随Tm变化公式。如图6（b）所示，不同Tm的∆TmlTin随着Tin呈一致的线性负相关关系，因此，∆TmlTin不受Tm的影响。取6种不同Tm的∆TmlTin的平均值，得到∆TmlTin随Tin的变化公式，如图6（c）所示。由于比表面积在不同体积时难以一致表示，因此，预测比表面积引起的热储终温增量∆TmlS时，应以在同体积下与比表面积变化规律一致的量—各表面积与相应体积下正六面体表面积的比值S/SR表示，以消除热储体积对∆TmlS的影响。如图6（d）所示，得到了∆TmlS随S/SR的变化规律。考虑∆TmlV随不同体积热储与边长为500 m的正六面体热储体积比值V/V500的变化规律，用以量化体积对Tml的影响，如图6（e）所示。所以得到式（13）中各分项的具体表达式为：
| $\left\{\!\!\!\! \begin{array}{l} {T_{{\rm{ml}}{\rm{353}}{\rm{.16}}}} = {\rm{0}}{\rm{.052\;18}}{T_{\rm{m}}} + {\rm{373}}{\rm{.845}}, \\ \Delta {T_{{\rm{ml}}{{{T}}_{{\rm{in}}}}}} = {\rm{0}}{\rm{.381\;2}}{T_{{\rm{in}}}} - {\rm{134}}{\rm{.571\;9}}, \\ \Delta {T_{{\rm{ml}}S}} = {\rm{ - 10}}{\rm{.205\;77}}\exp ({\rm{ - 0}}{\rm{.5}} \times ((S/{S_{\rm{R}}}- \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ 0}}{\rm{.943\;91}})/{\rm{0}}{\rm{.031\;14}}{)^2} + {\rm{1}}{\rm{.530\;14}}), \\ \Delta {T_{{\rm{ml}}V}} = {\rm{14}}{\rm{.493\;2}} \times {\rm{0}}{\rm{.006\;4}}{{\rm{8}}^{V/{V_{500}}}} + {\rm{0}}{\rm{.411\;2}} \end{array} \right.$ | (14) |
首先，得到前25年的补偿热量随热储初始温度Tm的变化量Qt25，基于此分别讨论Qt随t和SSa的增量变化∆Qtt和∆QtS，将上述3个量进行线性叠加。如图5（g）所示，当热储体积为0时，Qt也一定为0。令不同体积得到的补偿热量与边长为500 m的正六面体热储体积得到的补偿热量的比值为ηV，将ηV作为体积系数与Qt25、∆Qtt和∆QtS叠加得到的补偿热量相乘，以考虑体积因素对Qt的影响，如式（15）所示：
| $Q_{\rm{t}} = {\eta _{{V}}}({Q_{{\rm{t25}}}} + {\rm{\Delta }}{Q_{{{{\rm{t}}t}}}} + {\rm{\Delta }}{Q_{{\rm{t}}V}})$ | (15) |
Qt25与Tm的关系、∆QtS与S/SR的关系及ηV与V/V500的关系如图7（a）、（b）和（d）所示。由图7（e）可知：∆Qtt随着时间的延长不断增加，近似呈线性变化；但不同Tm对应的∆Qtt的增量不同，且随着Tm的上升，∆Qtt随着时间变化的线性曲线斜率不断增加，截距不断减小。故在考虑∆Qtt随着时间的变化时，还要考虑Tm的影响。分别计算出斜率c和截距d随Tm的变化，如图7（c）所示。最终得到式（15）中各分项的具体表达式为：
 |
| 图7 大地热补偿热量Qt计算 Fig. 7 Calculation of geothermal compensation heat (Qt) |
| $\left\{\!\!\!\! \begin{array}{l} {Q_{{\rm{t25}}}} = {\rm{ - 9}}{\rm{.696\;8}} \times {\rm{1}}{{\rm{0}}^{\rm{2}}}T_{\rm{m}}^{\rm{2}}{\rm{ + 1}}{\rm{.190\;9}} \times {\rm{1}}{{\rm{0}}^{{\rm{15}}}}{T_{\rm{m}}}- \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ 3}}{\rm{.338\;3}} \times {\rm{1}}{{\rm{0}}^{{\rm{17}}}} ,\\ {\rm{\Delta }}{{\rm{Q}}_{{{{\rm{t}}t}}}} = ({\rm{3}}{\rm{.494\;6}} \times {\rm{1}}{{\rm{0}}^{{\rm{10}}}}T_{\rm{m}}^{\rm{2}}{\rm{ - 3}}{\rm{.408\;3}} \times {\rm{1}}{{\rm{0}}^{{\rm{13}}}}{T_{\rm{m}}} +\\ \;\;\;\;\;\;\;\;\;\;\;\;\; {\rm{ 8}}{\rm{.325\;9}} \times {\rm{1}}{{\rm{0}}^{{\rm{15}}}})t{\rm{ - 7}}{\rm{.187\;4}} \times {\rm{1}}{{\rm{0}}^{{\rm{11}}}}T_{\rm{m}}^{\rm{2}}+ \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ 6}}{\rm{.934\;8}} \times {\rm{1}}{{\rm{0}}^{{\rm{14}}}}{T_{\rm{m}}}{\rm{ - 1}}{\rm{.677\;7}} \times {\rm{1}}{{\rm{0}}^{{\rm{17}}}}, \\ {\rm{\Delta }}{{\rm{Q}}_{\rm{t}}}_S = {\rm{1}}{\rm{.702\;21}} \times {\rm{1}}{{\rm{0}}^{{\rm{17}}}}{(S/{S\!_{\rm{R}}})^{\rm{2}}}{\rm{ - 3}}{\rm{.186\;03}}\times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ 1}}{{\rm{0}}^{{\rm{17}}}}S/{S\!_{\rm{R}}}{\rm{ + 1}}{\rm{.487\;24}} \times {\rm{1}}{{\rm{0}}^{{\rm{17}}}}, \\ {\eta _{{V}}} = {\rm{0}}{\rm{.934}}\;{{\rm{5}}^{V/{V_{500}}}}{\rm{ + 0}}{\rm{.055\;5}} \end{array} \right.$ | (16) |
以Tin为353.16 K时为基准，确定出Tout353.16随Tm的变化；基于此，分别讨论Tin、V和SSa对
| ${\overline T _{{\rm{out}}}} = {T_{{\rm{out353}}{\rm{.16}}}} + \Delta {T_{{\rm{out}}{T{{\rm{in}}}}}} + \Delta {T_{{\rm{out}}S}} + \Delta {T_{{\rm{out}}V}}$ | (17) |
式中，Tout353.16、∆ToutS和∆ToutV极易通过数值计算结果通过拟合曲线得到，如图8（a）、（c）和（d）。如图8（e）所示，随着Tin的减小，∆ToutTin逐渐增加，可近似看做线性曲线，且不同Tm的∆ToutTin在同一注入温度的增量不同，随着注入温度的上升而不断增加。因此，在计算∆ToutTin时，需要考虑Tm的影响。与∆Qtt计算方法类似，如图8（b）所示，分别计算出∆ToutTin随Tin变化线性曲线的斜率a和截距b随Tm的变化，即可得到∆ToutTin的计算值，则式（17）中各分项的具体表达式为：
 |
|
图8 平均产出温度
|
| $\left\{\!\!\!\! \begin{array}{l} {{{T}}_{{\rm{out353}}{\rm{.16}}}} = {\rm{ - 8}}{\rm{.784}} \times {\rm{1}}{{\rm{0}}^{{\rm{ - 4}}}}T_{\rm{m}}^{\rm{2}}{\rm{ + 1}}{\rm{.411\;4}}{T_{\rm{m}}}- \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{ 15}}{\rm{.128\;9}}, \\ {\rm{\Delta }}{{{T}}_{{\rm{out}}T{\rm{in}}}} = ({\rm{ - 0}}{\rm{.002\;16}}{T_{\rm{m}}} + {\rm{0}}{\rm{.986\;3}}){T_{{\rm{in}}}}+ \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\,\;\;\;\,{\rm{ 0}}{\rm{.769\;22}}{{{T}}_{\rm{m}}}{\rm{ - 353}}{\rm{.693\;5}}, \\ {\rm{\Delta }}{T_{{\rm{out}}S}} = {\rm{222}}{\rm{.877\;19}}{(S/{S\!_{\rm{R}}})^{\rm{2}}}- \\ \;\;\;\;\;\;\;\;\;\;\; \;\;\;\;\,{\rm{ 503}}{\rm{.588\;06}}S/{S\!_{\rm{R}}} + {\rm{279}}{\rm{.462\;73}}, \\ {\rm{\Delta }}{T_{{\rm{out}}V}} = {\rm{ - 34}}{\rm{.275\;6}} \times {\rm{0}}{\rm{.004\;0}}{{\rm{4}}^{V/{V_{{\rm{500}}}}}}{\rm{ - 1}}{\rm{.520\;1}} \end{array} \right.$ | (18) |
由表7可知，系数α只与V和SSa有关。为得到三者间的定量关系，将α分为αV和∆αS。首先，得到αV随V的变化，如图9（a）所示；然后，基于αV计算由SSa引起的增量∆αS，如图9（b）所示；最后，将αV与∆αS进行线性叠加得到系数α的计算式，如式（19）所示：
 |
|
图9 系数
|
| $\alpha = {\alpha _V} + \Delta {\alpha _S}$ | (19) |
如图9所示，最终得到式（19）中各分项的具体表达式为：
| $\left\{\!\!\!\! \begin{array}{l} {\alpha _V} = {\rm{0}}{\rm{.007\;5}}{V^{\rm{3}}}{\rm{ - 169\;948}}{\rm{.8}}{V^{\rm{2}}}- \\ \;\;\;\;\;\;\;\;\;\;{\rm{ 1}}{\rm{.082\;3}} \times {\rm{1}}{{\rm{0}}^{{\rm{14}}}}V{\rm{ + 1}}{\rm{.341\;1}} \times {\rm{1}}{{\rm{0}}^{{\rm{23}}}}, \\ {\rm{\Delta }}{\alpha _S} = {\rm{4}}{\rm{.833\;4}} \times {\rm{1}}{{\rm{0}}^{{\rm{21}}}}\exp ({\rm{ - 0}}{\rm{.5}}((S/{S_{\rm{R}}}- \\ \;\;\;\;\;\;\;\;\;\;\;\;{\rm{ 1}}{\rm{.033\;5}})/{\rm{0}}{\rm{.047\;6}}{{\rm{)}}^{\rm{2}}}){\rm{ - 4}}{\rm{.088\;2}} \times {\rm{1}}{0^{{\rm{21}}}} \end{array} \right.$ | (20) |
经上述分析，将式（13）、（15）、（17）、（19）代入式（12），得到Tml、Qt、
| $\left\{\!\!\!\! \begin{array}{l} {T_{{\rm{ml}}}} = {T_{{\rm{ml353}}{\rm{.16}}}} + {\rm{\Delta }}{T_{{\rm{ml}}T{\rm{in}}}} + {\rm{\Delta }}{T_{{\rm{ml}}S}} + {\rm{\Delta }}{T_{{\rm{ml}}V}}, \\ {Q_{\rm{t}}} = {\eta _V}({Q_{{\rm{t25}}}} + {\rm{\Delta }}{Q_{{{{\rm{t}}t}}}} + {\rm{\Delta }}{Q_{{\rm{t}}V}}), \\ {\overline T _{{\rm{out}}}} = {T_{{\rm{out353}}{\rm{.16}}}} + {\rm{\Delta }}{T_{{\rm{out}}T{\rm{in}}}} + {\rm{\Delta }}{T_{{\rm{out}}S}} + {\rm{\Delta }}{T_{{\rm{out}}V}}, \\ \alpha = {\alpha _V} + {\rm{\Delta }}{\alpha _S} \end{array} \right.$ | (21) |
式（21）中各分式可由式（14）、（16）、（18）、（20）表示。将式（21）代入式（5）中，结合勘测到的热储渗透率即可得到Tin、P、t三者关系，指导EGS建设。
3.5 计算方程的适用性验证目前，大多数研究主要利用工质循环流量表示热媒介质注入的量及热媒介质产出量，利用注采压差P与已有研究比较来验证计算公式适用性的方法难以直接实现。根据Dupuit公式可知，通过热储尺寸和渗透率可将注采压差P换算为热媒介质在生产井的平均流速，由生产井的尺寸即可得到热媒介质产出量，即工质循环流量。为验证计算方程的可靠性，将其应用至文献[15]的模型中。将文献[15]中的热储尺寸、热储初始温度、热媒介质注入温度和不同EGS运行时间代入式（5）中，得到不同EGS运行时间对应的注采压差P；然后，仅使用 COMSOL Multiphysics 中达西定律模块计算得到工质循环流量，所得结果如图10所示。研究发现，利用式（5）计算结果与文献[15]所得曲线结果一致，说明式（5）具有良好的适用性。
 |
| 图10 计算公式适用性验证 Fig. 10 Applicability verification of calculation formula |
4 结 论
利用理论分析与数值模拟相结合，描述EGS地下取热过程，主要结论如下：
1）将吸放热公式与Dupuit公式相结合，得到了双井EGS产能与寿命控制方程，揭示了EGS产能随热储渗透率、EGS寿命、热媒介质注入温度、注采压差和热储形状系数的变化关系。
2）通过理论分析分别得到热储减少热量Qm、大地热补偿热量Qt、EGS平均产出温度
3）将控制方程和预测公式与现有研究进行对比，控制方程和预测公式具有良好的适用性，可有效指导EGS建设，加强EGS运行期间的实时调控。
| [1] |
Crowley B. Statistical review of world energy[J]. Hydrocarbon Processing, 2000, 79(11): 23. |
| [2] |
Lu Chuan,Wang Guiling. Current status and prospect of hot dry rock research[J]. Science and Technology Review, 2015, 33(19): 13-21. [陆川,王贵玲. 干热岩研究现状与展望[J]. 科技导报, 2015, 33(19): 13-21. DOI:10.3981/j.issn.1000-7857.2015.19.00] |
| [3] |
Bertani R. Geothermal power generation in the world 2005—2010 update report[J]. Geothermics, 2012, 41: 1-29. DOI:10.1016/j.geothermics.2011.10.001 |
| [4] |
Bertani R. Geothermal power generation in the world 2010—2014 update report[J]. Geothermics, 2016, 60: 31-43. DOI:10.1016/j.geothermics.2015.11.003 |
| [5] |
Xu Tianfu,Hu Zixu,Li Shengtao,et al. Enhanced geothermal system:International progresses and research status of China[J]. Acta Geologica Sinica, 2018, 92(9): 1936-1947. [许天福,胡子旭,李胜涛,等. 增强型地热系统:国际研究进展与我国研究现状[J]. 地质学报, 2018, 92(9): 1936-1947. DOI:10.3969/j.issn.0001-5717.2018.09.01] |
| [6] |
Chen Jiliang,Jiang Fangming. Designing multi-well layout for enhanced geothermal system to better exploit hot dry rock geothermal energy[J]. Renewable Energy, 2015, 74: 37-48. DOI:10.1016/j.renene.2014.07.056 |
| [7] |
Asai P,Panja P,McLennan J,et al. Performance evaluation of enhanced geothermal system (EGS):Surrogate models,sensitivity study and ranking key parameters[J]. Renewable Energy, 2018, 122: 184-195. DOI:10.1016/j.renene.2018.01.098 |
| [8] |
Beardsmore G R. Thermal modelling of the hot dry rock geothermal resource beneath GEL99 in the Cooper Basin,South Australia[J]. Aseg Extended Abstracts, 2004(1): 1-4. DOI:10.1071/aseg2004ab009 |
| [9] |
Taron J,Elsworth D. Coupled mechanical and chemical processes in engineered geothermal reservoirs with dynamic permeability[J]. International Journal of Rock Mechanics and Mining Sciences, 2010, 47(8): 1339-1348. DOI:10.1016/j.ijrmms.2010.08.021 |
| [10] |
Wang Xiaoxing,Wu Nengyou,Zhang Keni,et al. Multi-field coupling for enhanced geothermal system development[J]. Hydrogeology and Engineering Geology, 2012, 39(2): 131-135. [王晓星,吴能友,张可霓,等. 增强型地热系统开发过程中的多场耦合问题[J]. 水文地质工程地质, 2012, 39(2): 131-135. DOI:CNKI:SUN:SWDG.0.2012-02-02] |
| [11] |
Sausse J,Dezayes C,Dorbath L,et al. 3D model of fracture zones at Soultz-sous-Forêts based on geological data,image logs,induced microseismicity and vertical seismic profiles[J]. Comptes Rendus Geosciences, 2010, 342(7/8): 531-545. DOI:10.1016/j.crte.2010.01.011 |
| [12] |
Arnaud B,Pierre G,Michel R,et al. Modeling the coupling between free and forced convection in a vertical permeable slot:Implications for the heat production of an enhanced geothermal system[J]. Geothermics, 2006, 35(5/6): 654-682. DOI:10.1016/j.geothermics.2006.11.008 |
| [13] |
He Luwu,Jin Zhihe. A local thermal none quilibrium poroelastic theory for fluid saturated porous media[J]. Journal of Thermal Stresses, 2010, 33(8): 799-813. DOI:10.1080/01495739.2010.482358 |
| [14] |
Karsten P. The TOUGH codes—A family of simulation tools for multiphase flow and transport processes in permeable media[J]. Vadose Zone Journal, 2004, 3: 738-746. DOI:10.2113/3.3.738 |
| [15] |
Chen Jiliang,Jiang Fangming. A numerical study on heat extraction performance of enhanced geothermal systems[J]. Renewable Energy Resources, 2013, 31(12): 111-117. [陈继良,蒋方明. 增强型地热系统热开采性能的数值模拟分析[J]. 可再生能源, 2013, 31(12): 111-117. DOI:10.3969/j.issn.1671-5292.2013.12.02] |
| [16] |
MIT.The future of geothermal energy:Impact of enhanced geothermal systems(EGS) on the United States in the 21st Century[R].Boston:Massachusetts Institute of Technology,2006.
|
| [17] |
Zhang Shengsheng,Zhang Lei,Cai Jingshou,et al. Preliminary estimation and evaluation of hot dry rock resources in Qiabuqia area of Gonghe basin[J]. Journal of Qinghai University, 2018, 36(4): 75-78. [张盛生,张磊,蔡敬寿,等. 共和盆地恰卜恰地区干热岩资源量初步估算及评价[J]. 青海大学学报, 2018, 36(4): 75-78. DOI:CNKI:SUN:QHXZ.0.2018-04-01] |
| [18] |
Zhai Haizhen,Su Zheng,Ling Lulu,et al. Study of EGS heat recovery by multi-parallel fracture model at desert peak field,USA[J]. Journal of Engineering Thermophysics, 2017, 38(1): 134-140. [翟海珍,苏正,凌璐璐,等. 基于平行多裂隙模型的美国沙漠峰地热田EGS热恢复研究[J]. 工程热物理学报, 2017, 38(1): 134-140. DOI:CNKI:SUN:GCRB.0.2017-01-02] |
| [19] |
Ren Guocheng. Talking about the calculation method of geothermal resources and reserves[J]. West-China Exploration Engineering, 2015, 27(11): 112-114. [任国澄. 浅谈地热资源及储量计算方法[J]. 西部探矿工程, 2015, 27(11): 112-114. DOI:10.3969/j.issn.1004-5716.2015.11.03] |
| [20] |
Muffler L J P,Christiansen R L. Geothermal resource assessment of the United States[J]. Pure and Applied Geophysics, 1978, 117(1/2): 160-171. |
| [21] |
Yan Yingjun.Application on Monte Carlo method in the typical geothermal assessment in China[D].Beijing:China University of Geosciences,Beijing,2009. 颜英军.蒙特卡罗法在中国典型地热资源评价中的应用[D].北京:中国地质大学(北京),2009. |
| [22] |
Wang Xinyi,Liao Zisheng,Han Xingxia,et al. Research on the method of Monte Carlo for assessing geothermal resources[J]. Hydrogeology and Engineering Geology, 2002, 29(4): 10-13. [王心义,廖资生,韩星霞,等. 地热资源评价的蒙特卡罗法[J]. 水文地质工程地质, 2002, 29(4): 10-13. DOI:10.3969/j.issn.1000-3665.2002.04.00] |
| [23] |
Li Tongbiao. A review of geothermal resource assessment methods[J]. Energy and Environment, 2015(5): 91-92. [李同彪. 地热资源评估方法综述[J]. 能源与环境, 2015(5): 91-92. DOI:10.3969/j.issn.1672-9064.2015.05.04] |
| [24] |
Noorollahi Y,Itoi R. Production capacity estimation by reservoir numerical simulation of northwest (NW) Sabalan geothermal field,Iran[J]. Energy, 2011, 36(7): 4552-4569. DOI:10.1016/j.energy.2011.03.046 |
| [25] |
Zhao Yangsheng,Feng Zijun,Feng Zengchao,et al. THM (Thermo-hydro-mechanical) coupled mathematical model of fractured media and numerical simulation of a 3D enhanced geothermal system at 573 K and buried depth 6 000–7 000 m[J]. Energy, 2015, 82: 193-205. DOI:10.1016/j.energy.2015.01.030 |
| [26] |
Gao Ping.Analysis of rock thermal physical parameters and research on multi-field thermal effect coupled model[D].Changchun:Jilin University,2015. 高平.岩石热物性参数分析及多场热效应耦合模型研究[D].长春:吉林大学,2015. |
| [27] |
Pandey S N,Vishal V. Sensitivity analysis of coupled processes and parameters on the performance of enhanced geothermal systems[J]. Scientific Reports, 2017, 7(1): 17057. DOI:10.1038/s41598-017-14273-4 |
| [28] |
Zeng Yuchao,Tang Liansheng,Wu Nengyou,et al. Analysis of influencing factors of production performance of enhanced geothermal system:A case study at Yangbajing geothermal field[J]. Energy, 2017, 127: 218-235. DOI:10.1016/j.energy.2017.03.100 |
| [29] |
Wang Changlong.Numerical simulation study of enhanced geothermal systems(EGS) considering fluid losses[D].Hefei:University of Science and Technology of China,2017. 王昌龙.考虑流体损失的增强型地热系统(EGS)数值模拟研究[D].合肥:中国科学技术大学,2017. |
| [30] |
Pruess K. On production behavior of enhanced geothermal systems with CO2 as working fluid
[J]. Energy Conversion and Management, 2008, 49(6): 1446-1454. DOI:10.1016/j.enconman.2007.12.029 |
| [31] |
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, 35(4): 351-367. DOI:10.1016/j.geothermics.2006.08.002 |
| [32] |
Spycher N,Pruess K. A phase-partitioning model for CO2—Brine mixtures at elevated temperatures and pressures:Application to CO2-enhanced geothermal systems
[J]. Transport in Porous Media, 2010, 82(1): 173-196. DOI:10.1007/s11242-009-9425-y |
| [33] |
White M D,Phillips B R.Code comparison study fosters confidence in the numerical simulation of enhanced geothermal systems[R].Richland:Pacific Northwest National Laboratory (PNNL),2015.
|
| [34] |
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. [陈继良,罗良,蒋方明. 热储周围岩石热补偿对增强型地热系统采热过程的影响[J]. 计算物理, 2013, 30(6): 862-870. DOI:10.3969/j.issn.1001-246X.2013.06.00] |
| [35] |
张三慧.大学物理学–第2册–热学[M].北京:清华大学出版社,1999.
|
| [36] |
Jelacic A,Fortuna R,LaSala R,et al.An evaluation of enhanced geothermal systems technology[R].Washington:Office of Energy Efficiency and Renewable Energy (EERE),2008.
|