COMSOL 博客

传热与地下多孔介质流的耦合仿真

作者 Phillip Oberdorfer

2014年 4月 24日

本文是地热能系列博客的第二篇，将重点探讨传热与地下水流耦合的过程，并据此确定地下储藏的热量是否足够多，是否值得开采以获得地热能。文中，我们将通过一个地下水换热回灌系统示例模型演示这一耦合过程。

深层地热能：潜能巨大，风险也巨大

地热能开采的最大挑战之一就是如何将勘查风险降到最低。例如，如何确定所选开采点的热能足够开采 30 年？通常情况下，几乎没有关于当地地下结构及水质等特性的信息，是否值得开采，不确定性相当大。

最近几十年中，数值模型成为评估这一风险的重要工具，通过这种方法可以将不确定性限制在一个合理的范围，然后再进行参数化研究。今天，我将简单地介绍地下水流与传热耦合问题的数学描述，这个问题在许多地热应用中都会遇到。我还将向您展示如何利用 COMSOL 软件研究和预测（水利）地热系统的性能。

地下水换热系统中的方程

可以用下方程描述地下传热：

(1)

(\rho C_p)_{eq} \frac{\partial T}{\partial t} + \rho C_p {\bf u } \cdot \nabla T = \nabla \cdot (k_{eq} \nabla T ) + Q + Q_{geo}

热量经传导和对流过程实现平衡，并且可以通过在源项 Q 中定义热量来表示热量的生成或流失。多孔介质传热 接口中有一个特殊的地热采暖 功能，表示为域条件：Q_{geo}。

COMSOL 中还有一个功能可以使地热能仿真工程师的工作稍微更轻松些。目前可以做到的是将热力参数的平均表示实现为一个权重因子，这些热力参数包括岩石基质和地下水，用矩阵体积分数 \theta 表示。对于几种固定的固体和流体，可以选择使用体积平均方法，也可以采用幂律平均方法。

要使用体积平均方法，热传导方程中的体积比热变为：

(2)

(\rho C_p )_{eq} = \sum_{i} ( \theta_{pi}\rho_{pi}C_{p,pi})+(1-\sum_{i}\theta_{pi})\rho C_p

热导率变为：

(3)

k_{eq}=\sum_{i} \theta_{pi} k_{pi} + ( 1-\sum_{i} \theta_{pi} ) \rho C_p

要正确求解传热，还要考虑与流场相结合。一般而言，地下情况复杂多变，要求我们使用不同的方法对地下水流进行数学描述。如果侧重于微尺度，解决孔隙空间内水的流动，则需要求解蠕动流方程或 Stokes 流方程。对于非饱和区，需要求解 Richard 方程，它常用于有关环境污染的研究中（例如，可以参阅我们之前发表的博客文章：模拟农药径流，涕灭威的影响）。

而深层地热层中的饱和流及主要由压力驱动的流动则完全可以通过 Darcy 定律求解：

(4)

{\mathbf u} = -\frac{\kappa}{\mu} \nabla p

其中，速度场 \mathbf{u} 取决于渗透率 \kappa，流体的动态粘度 \mu 由压力梯度 p 驱动。然后，将 Darcy 定律并入以下连续性方程：

(5)

\frac{\partial}{\partial t} (\rho \epsilon_p) + \nabla \cdot ( \rho {\bf u} ) = Q_m

如果地质年代较为久远，则由于水流具有储能效应，有关时间的项可以忽略不计。由此，上述方程左边第一项为 0，因为密度 \rho 和孔隙率 \epsilon_p 都可以看作常数。通常，水力属性随温度的变化可以忽略不计。因此，（稳态）流动方程与（时间相关的）传热方程互不相关。在有些情况下，特别是当自由度数量相当大时，利用这种无关性将此问题拆分为一个稳态研究步骤和一个时间相关的研究步骤是很有意义的。

裂隙流和多孔弹性

裂隙流是局部地热系统中主要的水流动态，例如在卡斯特含水层系统中。地下水流模块中包含的裂隙流接口可用于表示二维裂隙和裂缝中的 Darcy 流场。

地下水换热提取系统通常包含一个或多个回灌井和生产井。很多情况下这些井都是单独钻取的，但目前更多的做法是钻取一个（或多个）分支井。还有人提出仅钻取单个井，再分出单独的回灌区和生产区。

需要注意的是，由于水的回灌和抽取而产生的人为压力变化会对多孔介质的结构产生影响，并产生水力压裂。如果考虑这些因素，则还需要进行多孔弹性分析，但在这里我们暂不考虑。

地下水换热应用的 COMSOL 模型：地热回灌

在 COMSOL Multiphysics 中可以非常简单地建立一个对地下水-地热应用作长期预测的模型。

模型区包含三层具有不同热性能和水力属性的地质层，位于一个体积约为 500 立方米的箱体中。该箱体表示处于一个大断裂带上的地热开采点的一部分，层高为源自外部数据集的插值函数。相关含水层已完全饱和，顶部和底部均为弱透水层（隔水层）。温度分布通常是一个不确定性因素，但我们可以假设地温梯度为 0.03 [°C/m]，故初始温度分布：T₀(z)=10 [°C] – z·0.03 [°C/m]。

使用 COMSOL Multiphysics 创建的分层地下区域中地下水换热对井回灌系统的几何
分层地下区域中的地下水换热回灌系统，处于一个断裂带上。边缘长度约为 500m。左侧为回灌井，右侧为生产井。两井之间的横向距离约 120 m。