3.3.2　数值法

3.3.2　数值法

1.基于等效渗流通道模型

（1）等效渗流通道理论模型

等效渗流通道模型假设地热对井之间存在若干优势渗透通道，其主导了热储中的地热水渗流（图3-6）。优势渗流通道可能是近乎垂直的断裂带的一部分，或是水平岩层的一部分，并且优势流动通道内的渗流可近似等效为一维达西渗流。

图3-6 等效渗流通道模型示意图

（2）溶质运移理论

示踪剂在一维等效渗流通道中的主要运移机理为对流和水动力弥散，如忽略通道与其他热储岩石间的扩散过程及通道内的横向水动力弥散作用，等效渗流通道内的示踪剂运移过程可由经典的多孔介质中对流-扩散方程来描述：

式中，C表示等效渗流通道中的示踪剂质量浓度，kg/m³；t表示时间，s；x表示通道轴向坐标，m；D表示通道轴向弥散系数，m²/s；v表示通道中的平均流速，m/s。

在t=0时刻从回灌井瞬时投入质量为m的示踪剂，其他初始条件和边界条件可以表示为

式中，y表示通道纵向坐标，m。

结合初始条件和边界条件，偏微分方程（式3-8）的解为

式中，m_c为进入渗流通道中的示踪剂质量，kg；A为渗流通道的截面积，m²；ϕ为孔隙度；q_c为渗流通道中的回灌率，kg/s；ρ为水的密度，kg/m³；α_L为纵向弥散度，m。

考虑到开采井中示踪剂的质量守恒，在开采率为Q（kg/s）的情况下，有下列关系成立：

式中，c为开采井中示踪剂的质量浓度，kg/m³。

如果地热对井之间存在n条流动通道，开采井中示踪剂质量浓度表达式为

其中

式中，x_i为第i条渗流通道的轴向坐标，m；q_i为第i条渗流通道中的回灌率，kg/s；q_in为回灌率，kg/s；v_i为第i条渗流通道中的流速，m/s；为第i条渗流通道中的纵向弥散度，m；A_i为第i条渗流通道的截面积，m²；M_i为第i条渗流通道中示踪剂的质量，kg。

（3）反演模型

式3-13中的未知参数（如通道长度x_i、平均流速v_i、弥散度以及通道截面积A_i等）无法直接测定，须通过示踪试验反演求得。本文引入移动近似法（method of moving asymptotes，MMA），建立热储参数反分析理论框架。目标函数即原始优化问题P可定义为

式中，h_i（x）为由式3-12计算的示踪剂质量浓度，kg/m³；为示踪试验中示踪剂质量浓度的测量值，kg/m³；和分别表示给定的下界和上界；x={x₁，x₂，…，x_n}^T为式3-12中的未知参数。

拟合优度R²是衡量拟合程度的重要指标，其定义如下：

式中，h_mi为第i组质量浓度测量值；h_ei为第i组质量浓度估计值；为测量质量浓度平均值。

R²的取值范围为[0，1]，R²的值越接近1，说明拟合程度越好；反之，R²的值越接近0，说明拟合程度越差。

MMA算法由Svanberg提出，是特别适合应用于多变量的优化问题的一种算法。MMA算法需要先指定优化变量的上下界，之后按照如下步骤进行优化计算。

①选择初始计算点x^（1）。

②计算在该迭代点处的目标函数和约束函数各自的函数值f_i[x^（k）]和梯度值。

③生成一个子问题P^（k）来近似原始优化问题P，基于步骤②的计算，原始优化问题P中的函数f_i被子问题P^（k）中的函数近似替代。

④求解子问题P^（k），并令所得到的最优解成为下一个迭代点x^（k+1），然后重复步骤②～④。

⑤满足收敛判据时，迭代停止。计算所得子问题P^（k）的最优解就是原始优化问题P的唯一最优解。

Svanberg给出了如下收敛判据的形式：

式中，ε对于所有的j=1，…，n是一个极小值。

MMA算法的优势在于可同时优化多个变量，对优化变量的上下界不敏感，且极少出现迭代不收敛的情况。

（4）对流传热理论

假设水-岩界面温度相等，等效渗流通道中的对流换热控制微分方程为

式中，T为温度，K；（ρc）_f为渗流通道中材料的体积热容，J/（m³·K）；k_r为热储岩石的导热系数，W/（m·K）；b为渗流通道的宽度，m；ρ_w为水的密度，kg/m³；ρ_r为热储岩石密度，kg/m³；c_w为水的质量热容，J/（kg·K）；c_r为热储岩石的质量热容，J/（kg·K）；ϕ为岩石孔隙度。

式3-17可进一步化简为如下形式：

令

式中，h为渗流通道的高度，m。

则式3-18可化简为

当h≫b时，初始条件和边界条件为

所以，渗流通道中水的温度表达式为

式中，erfc（ ）表示误差补函数，它与误差函数erf（ ）的关系为

所以，式3-22可以写为

式中，（ρc）_f表示流动通道中材料的体积热容，这其中的材料不仅有水，还有多孔的热储岩石，因此，k_r/（ρc）_f这一项其实是考虑了流动通道中材料的热扩散系数，而k_r/（ρ_rc_r）只是考虑了热储岩石的热扩散系数，整个热储层是由热储岩石和地下水共同组成的，因此，在理论上，考虑整个热储层的热扩散系数要更严谨一些。考虑到开采率为Q，则在不同时刻开采井中的水温T_p为

如果两井之间有i条流动通道，则开采井中的水温T_p可以通过以下表达式计算：

（5）算例

为了验证MMA算法是否适合优化多参数问题，假设表达式3-13中的所有参数都是已知的，则参数取值如表3-10所示。通过计算，可以得到不同时间点的质量浓度随时间分布的散点图，将其作为外部数据输入到模型中进行多参数的优化计算，最后比较原始给定的参数取值与优化参数值之间的误差。试验数据点分布如图3-7所示。

表3-10 表达式参数取值

图3-7 试验数据点分布图

MMA算法支持约束，因此，可以指定优化参数的取值范围，优化参数的取值范围如表3-11所示。

表3-11 优化参数的取值范围

优化拟合结果如图3-8所示，参数优化结果与原始参数值对比及误差如表3-12所示。

表3-12 优化结果与初始取值对比

图3-8 优化拟合结果示意图

由表3-12可以看出，优化结果与初始取值非常接近，误差控制在2%以内。因此，可以证实，MMA算法不仅能优化多参数问题，其误差也在合理的范围之内。

（6）数值模拟结果分析

①参数反演结果分析

根据鲁北院办公地示踪试验所取得的数据，钼元素质量浓度随时间变化的散点分布图如图3-9所示。

图3-9 钼元素质量浓度随时间变化的散点分布图

通过分析所测得的示踪剂质量浓度随时间分布的散点分布图（图3-10），发现：在中后期，由于取样频率过高，加大了反分析的计算量。因此，为了能更好地分析示踪试验数据，对测得的数据进行人为处理；为了提高优化程序的计算效率，将取样频率处理为每天取样一次，对于每天多次取样的数据点，取其平均值作为该天的质量浓度。考虑到储层的底值质量浓度较高，约为1.4×10^-5kg/m³，不利于反分析，将前60d的质量浓度均取为0kg/m³，处理后的示踪剂质量浓度随时间变化的散点分布图如图3-10所示。

图3-10 处理后的示踪剂质量浓度随时间变化的散点分布图

由于示踪试验数据点分布出现了一个明显的质量浓度峰值，因而采用了之前所提出的基于等效渗流通道模型的热储参数反演理论框架对本次示踪试验数据进行反分析，拟合所得曲线如图3-11所示，热储参数计算结果如表3-13所示。计算所得的拟合优度R²=0.603＞0.600，拟合程度较好，具有一定的参考价值，为未来进一步研究德州市的砂岩热储资源奠定了理论基础。

图3-11 示踪试验数据点与拟合曲线示意图

表3-13 热储参数计算结果

(https://www.daowen.com)

②热突破预测结果与分析

回灌的“冷水”通过渗流通道被加热，这取决于渗流通道的表面积而不是体积，因此，必须对流动通道的几何特征做出一定的假设。考虑到初始条件和边界条件存在的前提是h≫b，结合示踪试验反分析得出的参数取值，对开采井中水的温度变化进行预测，相关参数取值如表3-14所示。

表3-14 相关参数取值

通过分析示踪试验数据可以得到示踪剂的回收率M_i/M，计算所得示踪剂回收率为0.3982%，计算所得开采水温在100a内的变化曲线如图3-12所示。本次回灌示踪试验的开采量及回灌量较小，且示踪剂回收率较低，说明采灌井之间的水力联系较差，大部分的回灌水通过其他流动通道流向热储层中，仅有一小部分回灌水通过流动通道到达了开采井，开采井温度下降幅度很小，仅为0.2℃。因此，基于目前的研究预测，该地热对井系统在长期回灌条件下，在100a的使用期限内开采水温变化不明显，不会发生热突破。但是，该地热对井运行时间只有2a，相关数据很少，这方面的研究有待以后进行。

图3-12 开采水温在100a内的变化曲线

③参数敏感性分析

为了讨论渗流通道宽度b、回灌量q和对井间距x对热突破的影响程度，试验将用单通道模型的相关数据对以上参数进行敏感性分析，为进一步优化地热回灌提供理论依据。

a.当b取0.1m、0.2m、0.3m，对应的h取6.4m、3.2m、2.13m时，开采井中温度随时间的变化规律如图3-13所示。随着b的不断增大，开采井温度的下降幅度也不断增大，这是由于该模型只考虑上下平板的热传导作用；b越大，上下平板的间距就越大，在相同时间内通过渗流通道的水被加热得就越不充分，开采井温度下降的幅度也就越大。

图3-13 开采井温度随时间的变化规律

b.开采量保持65m³/h（18.06kg/s）不变，回灌量分别取20m³/h（5.56kg/s）、40m³/h（11.12kg/s）、60m³/h（16.68kg/s）时，开采井温度在10a内的变化曲线如图3-14所示。随着回灌量的不断增大，开采井温度下降幅度也在不断增大（表3-15）。这是由于采灌井之间流动通道中的回灌水流速大小主要由回灌量决定，由式3-17、式3-19可知，在开采量保持不变的情况下，回灌量q越大，则β越大，因此，t就越小，回灌水到达开采井的时间就越短，换热时间就越短，故温度下降幅度也就越大。

图3-14 不同回灌量下开采井温度在10 a内的变化曲线

表3-15 不同回灌量下温度下降幅度

c.开采量保持不变65m³/h（18.06kg/s），回灌量保持20m³/h（5.56kg/s）不变，对井间距x分别取300m、600m、900m时，开采井温度在10a内的变化曲线如图3-15所示。随着对井间距的不断增大，开采井温度变化的幅度在不断减小。根据式3-24得出，在β保持不变的前提下，x越大，t就越大，因此，回灌水到达开采井的时间就越长，在流动通道中被加热的时间就越长，最后开采井温度变化幅度也就越小（表3-16）。

图3-15 不同对井间距下开采井温度在10a内的变化曲线

表3-16 不同对井间距下开采井温度下降幅度

（7）小结

①将理论模型应用到工程实例中，结合示踪试验数据，进行了参数反演计算，计算得到拟合优度R²=0.603＞0.600，说明此优化结果与测量值有较强的相关性，拟合程度较好。

②结合示踪试验反分析的热储参数取值，对该地热对井系统进行了热突破预测。结果表明示踪剂回收率偏低，说明采灌井之间的水力联系较差，在100a的使用期限里，开采温度变化幅度很小，不会发生热突破现象。

③对主要影响开采井热突破的参数进行了敏感性分析，发现：渗流通道宽度b越大，开采井地热水温度下降的幅度就越大；回灌量q越大，回灌地热尾水到达开采井的时间就越短；井间距x越大，回灌水到达开采井的时间就越长。因此，在实际工程中，砂岩热储回灌量不宜过大，回灌地热尾水温度不宜过低，对热量不宜“吃干榨净”；同时，应尽量设计较大的采灌井间距。

2.地热回灌中井壁与流体换热现象的研究模型

等效渗流通道模型对示踪试验数据进行了合理的解释，但是，此模型在其理论的假设基础上只能计算出开采井井底的温度变化情况，而开采的热水从井底抽到井口，必然会有井中热水与井壁周围岩石的热交换现象，这可能会导致开采热水到达井口的温度不同于井底温度。因此，在这一部分，建立了一个从井口到井底的完整模型来研究井中热水在抽采过程中与井壁周围岩石的热交换现象以及该现象对开采水温的影响程度。

（1）理论推导

回灌井中的回灌流体从井口注入井底的过程、开采井中的开采“热水”从井底抽采到井口的过程，井中流体与周围岩石会发生热交换。描述这一物理过程，需要对井中流体的能量、动量和质量守恒方程，以及井壁附近岩石的能量守恒方程进行求解。

传热的基本定律是热力学第一定律，通常也被称为能量守恒原理。然而，内能U是一个难以测量的量，因此，传热方程通常改写成用温度T来表示的形式，对于井中的流体而言，忽略黏滞耗散，则传热方程可以表示为

对于井壁附近的岩石而言，其传热方程为

式中，ρ_f、ρ_s为流体密度和固体密度，kg/m³；c_p，f、c_p，s分别为流体、固体的常压比热容，J/（kg·K）；T为温度，K；v为速度矢量，m/s；k_f、k_s分别为流体、固体的导热系数，W/（m·K）；Q为热源，W/m³。

流体传热方程中的速度矢量v依赖于流体的流动状态，由于井径很小，因此假设井中流体的流动状态为层流，则该过程可以由纳维-斯托克斯（Navier-Stokes）方程来描述，其连续性方程和动量方程分别为

式中，p为压力，Pa；τ为黏性应力张量，Pa，对于可压缩流体而言，τ可以表示为

式中，μ为动力黏度，Pa·s；F表示体积力，N/m³。

在深层热储中，含水层是完全饱和的，那么饱和压力流可以用达西定律来描述：

式中，v′为达西流速，m/s；κ为多孔介质渗透率，m²；p_f为流体压力，Pa。

流体连续性方程可以表示为

式中，Q_m为质量源项，kg/（m³·s）。

回灌流体进入热储层后的传热过程可以用以下多孔介质传热方程来描述：

式中，（ρc_p）_eff为常压有效体积热容，对于多孔介质而言，可以表示为

式中，θ_s为固体材料的体积分数，与液体体积分数θ_f（或孔隙度）的关系为

多孔介质的有效导热系数k_eff与固体材料的导热系数k_s和流体材料的导热系数k_f有关，可以表示为k_s和k_f的加权算术平均值：

（2）三维数值计算模型

①模型定义

根据钻孔资料显示，地层自上而下：0～270.00m为第四纪平原组，厚270m；270.00～1290.40m为新近纪明化镇组，厚1020.40m；1290.40～1494.40m为新近纪馆陶组，厚204m；1494.40～1603.10m为古近纪东营组，厚108.70m（未揭穿）。取水层主要为新近纪馆陶组和古近纪东营组的砂砾岩和细砂岩。因此，数值计算模型根据地质钻孔资料将地层简化为5层，三维数值计算模型尺寸为800m×800m×1500m，回灌井与开采井在三维数值计算模型中的相对位置如图3-16所示。将模型自由剖分成四面体网格，并细化井中及井附近区域的网格，共剖分为144917个网格单元，如图3-17所示。

图3-16 回灌井与开采井在三维数值计算模型中的相对位置

图3-17 三维数值计算模型网格剖分

②初始条件及边界条件

根据资料显示：地层的温度分布随深度的变化存在一个温度梯度，经过计算可得到该温度梯度为0.0311℃/m。因此，初始条件设置为在整个地层施加一个初始温度分布T₀（y）=25-0.02y。

水力边界条件为在回灌井井口设置流量流入边界，回灌流量为10.36m³/h，在井底取水段设置流量流出边界；在开采井井底取水段设置流量流入边界，在井口设置流量流出边界，流入等于流出，开采流量为10.53m³/h；对于含水层，在回灌井井底取水段设置流量流入边界，在开采井井底取水段设置流量流出边界，含水层左边和右边均表示无限远处，因此直接指定其水位为0m。

温度边界条件为在回灌井井口设置温度流入边界，在井底取水段设置温度流出边界；在开采井井底取水段设置温度流入边界，该温度值等于热储层中开采井取水段附近的平均温度，在井口设置温度流出边界；对于热储层，在回灌井井底取水段设置温度流入边界，假设该温度值等于回灌井流出温度的平均值；在开采井井底取水段设置温度流出边界，模型左边和右边均表示无限远处，并将其边界设置为开边界。

③数值模型计算结果

在三维数值计算模型中，图3-18表示含水层渗透系数服从统计规律时，不同时刻地层的温度场和渗流场分布。从流线分布可以看出，在抽水和注水的过程中，井中流体的温度会影响井壁周围岩石的温度变化，井壁周围岩石的温度分布是以井为圆心的同心圆区域。考虑到本次回灌量与开采量偏小以及大地热流的影响，结合如图3-19所示的开采井温度变化曲线可以看到，开采井温度的变化幅度很小，不大于0.55℃。因此，可以预测该地热对井系统在50a的使用年限里并不会影响地热井的开发利用。

图3-18 不同时刻地层的温度场和渗流场分布云图

图3-19 开采井温度变化曲线

（3）小结

通过建立三维模型，揭示了流体在热储层中的渗流规律、井中的流动规律以及热储层温度和开采水温的变化规律。在开采过程中，开采井井壁附近的岩石逐渐被开采热水加热，因此，随着开采时间的不断增加，井中流体与井壁岩石之间的热交换逐渐减小，故通过分析开采井井底温度的变化情况来预测开采水温的变化规律是可行的。

3.采灌井间距的研究

地热回灌过程中应避免开采井发生热突破，因此需要确定采灌井的合理井间距。本文根据回灌数值模型开发App，通过对采灌井间距、开采及回灌水温、开采及回灌水量等参数的控制，模拟一定时期内的热突破曲线，进而分析、确定采灌井间距。

地热回灌过程主要涉及渗流场和温度场的耦合，并且按回灌井和开采井所处的热储层性质将回灌分为同层回灌和异层回灌，该工程实例属于同层回灌，因此，为了更好地求解多物理场耦合问题，数值计算模型在实际工程的问题上进行了相应的简化处理，将研究重点区域放在取水段热储含水层，模型主要运用式3-32～式3-37，其概念模型示意图如图3-20所示。

图3-20 概念模型示意图

考虑到地热供暖系统一年里只有冬季的4个月运行，其余8个月停运的情况，可以采用一个矩形周期函数来描述这种系统的运行和停运，其作用类似于开关，该函数图像如图3-21所示，横坐标表示时间，周期为1a；纵坐标的1表示系统处于运行状态，0表示系统处于停运状态。

图3-21 矩形周期函数图像

数值模型中考虑了大地热流密度的影响，因此，在模型底部施加了一个边界热源，模型的初始及边界条件如图3-22所示。根据鲁北砂岩热储的岩性特征，模型参数取值如表3-17所示。

图3-22 模型的初始及边界条件

表3-17 模型参数取值列表

模型的渗流场与温度场分布云图如图3-23所示。从图中可以看出：回灌冷水从回灌井进入热储中，在渗流场的作用下，回灌冷水逐渐向开采井方向运移，并导致开采温度的降低，但是由于一年里只有4个月的运移时间，因此，有效的开采时间其实只有16.67a左右。给定参数为开采量与回灌量均为60m³/h，结合如图3-24所示的开采井温度变化曲线可以看到，50a开采井温度的降低幅度约为1℃。以地热井的一般使用寿命30a计，其温度降低幅度约为0.65℃。因此，预测该地热对井系统在采灌间距为300m，采灌量同为60m³/h的情况下，50a的使用年限里不会影响其开发利用。

图3-23 模型的渗流场与温度场分布云图

图3-24 开采井温度变化曲线

---

[1] 指山东省鲁北地质工程勘察院，下同。