含水层水力联系温度-渗流耦合模型研究

刘殿德 陈志峰 糜行军

（滕州郭庄矿业有限责任公司，山东 滕州 277500）

探明含水层间水力联系是地下水资源和水文地质勘探研究的基础[1]，是预防矿井水害的重要工作之一[2]。含水层之间水力联系是一个复杂的系统工程[3]，需要综合多因素多方法确定，包括含水层水化学分析、含水层水位标高变化分析、构造赋存情况等因素。对于水化学分析方法，对比不同含水层的地下水矿物离子成分及其变化，有时会结合现场示踪试验方法，来判断含水层之间的水力联系[4]；对于含水层水位标高变化分析方法，可以对含水层进行抽水试验，研究不同含水层之间的水头高度变化，判断同一含水层不同区域以及不同含水层之间的水力联系。

在一些工程现场进行抽水试验时，通常还会利用水温变化判断含水层之间的水力联系。由于地温梯度的影响，含水层之间的温度也会不同，通常情况下深度越大，温度越高，可以利用温度来判断地下水所处的层位。但是，在含水层水力联系研究方面，从理论上对地温场的变化规律研究较少，缺少相关物理力学模型。针对上述存在的问题，基于多孔介质渗流、热力学相关理论，建立温度-渗流耦合条件下含水层水力联系数值模拟模型，研究不同含水层水力联系状态、不同回采距离条件下渗流场与温度场的变化规律。

1 控制方程

1856 年，法国水力学家达西（H.Darcy）通过大量的实验，得到线性渗透定律，即著名的达西定律：

式中：Q为通过渗流流量，即实验过程中通过砂柱各断面的流量，m3/s；K为渗透系数，m/s；A为过水断面，m2；I为水力梯度，等于过水断面之间的距离与水头损失的比值。可以看出，渗透系数可以理解为水力梯度为1 时（单位水力梯度）的渗流速度。同时应注意到达西定律实验得出的是渗透系数，而不是渗透率。

上述达西定律的形式还不能够进行数值模拟计算，需要将其转换成微分形式。达西定律也可以用压力P和流体的密度ρf表示，即可以写成如下形式[5]：

式中：K′即为渗透性常数，m3·s/kg；ρf为地下水的密度，kg/m3；P为水的压力，Pa；h为渗流断面之间的高度差，m。可以看出，其量纲包含了质量、长度和时间单位，表示了某种介质对某特定的流体的渗透能力，它的大小由介质和流体两者的性质而定[5]。但是，这种形式还是无法将多孔介质的影响与流体的影响分开，于是，便进一步引入一个新的物理量，其除以流体的黏度的值等于渗透性常数，即：

式中：κ为多孔介质的渗透率，m2；μ为地下水的动力黏度，Pa·s。渗透率的性质取决于多孔介质的内部结构，与流体的性质无关。

接下来将式（3）中的渗透性常数替换为渗透率，同时将压力差与h的比值表示为压力梯度，就可以将达西定律表示为常见的微分形式，即：

式中：v为地下水的渗流速度，m/s；p为孔隙水压力，Pa；g为重力加速度，m/s2。

地下水在地下岩体中的流动还需要满足连续性方程，即：

上式用来描述地下水在含水层中的流动问题，多孔介质内传热问题通常用基于局部热平衡理论和体积平均方法导出的多孔介质传热方程描述，即[6]：

式中：cf为地下水的恒压比热容，J/（kg·K）；cs为固体的恒压比热容，J/（kg·K）；qeff为多孔介质导热通量，J/（s·m2）；（ρcp）eff为等效单位体积热容，J/（m3·K）；ks为固体的导热系数，W/（m·K）；kf为流体的导热系数，W/（m·K）；keff为多孔介质的等效导热系数，W/（m·K）；T为温度，K。对于多孔介质，固体骨架和流体同时存在于一个体积空间内，两者热力学特征存在差异，需要使用等效的热力学参数将两者统一起来，得到多孔介质的能量方程。在累积项中，多孔介质内的等效单位体积热容可以使用式（10）进行计算。

2 数值模型建立

多物理场耦合数值模拟软件COMSOL Multiphysics 能够将不同物理场的控制方程耦合，模拟研究多物理场耦合问题。使用COMSOL Multiphysics 多孔介质流动模块和多孔介质传热模块，将达西定律、多孔介质流动控制方程以及多孔介质传热方程耦合，建立含水层水力联系温度场-渗流场耦合模型，实现地下水在地下岩体多孔介质中的流动过程和传热过程的耦合。具体包括建立几何模型、设置物理场与控制方程、设置初始与边界条件以及制定研究方案求解温度场-渗流场耦合控制方程。

软件中的达西定律接口用于描述流体在完全饱和的多孔介质中通过间隙的流动，这种运动由压力梯度驱动，流体的剪切应力引起的动量传递可以忽略不计。可以使用达西定律接口计算压力，然后根据压力梯度、流体黏度和渗透率来确定速度场。其控制方程就是式（6）。

软件中的多孔介质传热接口是通过传导、对流和弥散进行的。弥散是由液体在多孔介质中的曲折路径引起的，并通过包含平均对流以外的其他效应进行描述。在许多情况下，固相可以由电导率不同的多种材料组成，也可能会涉及许多不同的流体。多孔介质传热接口可以自动分析这些因素，并提供用于计算有效传热属性的混合规则。孔隙空间中的流体也可以经历一个或多个相变，这在模拟冻土过程时很有意义。可以使用专门的相变材料特征通过指定两种材料和多个相变属性（相变温度、转变间隔和潜热）来模拟此过程和类似的过程。其控制方程就是式（8），其参数可以在相应的节点处定义。

2.1 几何模型

COMSOL Multiphysics 软件包提供了丰富的几何建模工具，支持使用实体对象、表面、曲线和布尔操作来创建零件。可以通过操作序列来定义几何，序列中的每个操作都可以接收输入参数，方便在多物理场模型中轻松进行编辑和参数化研究。几何与定义的物理场设置之间是完全相互关联的，只要几何发生变化，系统便会自动将相关变化传递到整个关联的模型设置中。可以将材料域和表面等几何实体进行分组，创建不同的选择，并在定义物理场、网格划分以及绘图等后续操作中使用这些选择。不仅如此，还可以使用一系列操作来创建参数化几何零件（包括相关选择），然后将它们存储到“零件库”中，以便在多个模型中重复使用。

为了研究不同回采距离下含水层之间的水力联系以及渗流场、温度场之间的关系，建立三维几何模型。模型长500 m、高400 m，沿煤层走向厚度600 m，煤层厚度5 m，工作面宽150 m（见图1）。

图1 三维几何模型

2.2 物理场设置

COMSOL 软件提供了一系列预定义的物理场接口，用于模拟各种物理现象，包括许多常见的由多个物理场共同作用引起的现象。每个物理场接口都提供专用于相关科学或工程领域的特定设置。当选择所需的物理场接口之后，软件会给出适用的研究类型建议，例如瞬态或稳态求解器。选定求解器后，将对数学模型、求解器序列以及可视化和结果设置进行适当的数值离散化。同时，可以对所有设置进行编辑操作。

渗流场选择达西定律，包含了达西定律方程和连续性方程，在其内部各节点可以设置流体的密度、动力黏度等参数，同时可以设置不同域对应的控制方程，本模型所有域均为达西定律方程；温度场选择多孔介质传热方程，在其内部各节点可以设置流体和固体的热力学参数，比如导热系数、比热容，并选择有效导热系数计算多孔介质的热力学参数，同时本模型所有域均设置为多孔介质传热控制方程。

2.3 初始条件与边界条件

分别对渗流场和温度场设置初始条件和边界条件。渗流场的边界条件设置分为两种情况：一种情况是设置垂直方向的压力梯度，模型上部边界为定压力边界，水压值根据研究方案确定，工作面边界设置为0.1 MPa，其余边界设置为无流动；另一种情况是设置水平方向的压力梯度，给定模型左侧边界为定压力边界，水压值根据研究方案确定，模型右侧边界设置为0.1 MPa，其余边界设置为无流动。对于温度场，整个区域的地温梯度设置为0.025 ℃/m，模型的上下边界设置为定温度边界，分别为25 ℃和35 ℃。

2.4 研究方案

COMSOL Multiphysics 软件可以提供瞬态、稳态或特征频率研究，同时可以进行扫描、优化和估计研究，通过参数化扫描来运行，其中可以包含模型中的一个或多个参数，涉及几何参数、物理场定义中的设置等。可以使用不同的材料及其定义的属性来执行扫描，也可以对一组定义的函数执行扫描。

本模型采用稳态研究，为了研究不同回采距离范围内、不同水压、不同压力梯度方向条件下渗流场和温度场的空间分布，对开采距离、水压值进行参数化扫描，两种情况下渗流场上游水压设置为0.5～3 MPa。

3 结果

不同工作面回采距离会影响工作面围岩温度场的分布，不同水压值和压力梯度方向决定了渗流速度的变化，地下水渗透压和渗流速度会改变地下岩体多孔介质的热力学性质，进而决定了地温的变化。

3.1 垂直方向温度场分布

图2 展示了回采350 m 时，沿回采方向不同位置处、不同水压下温度场的分布，沿回采方向（Z轴方向）选取3 个截面，分别为Z=0 m、Z=300 m和Z=600 m。Z=600 m 截面表示未回采区域，可以看出，水平方向温度相同；Z=0 m 和Z=300 m 表示已回采区域的不同位置，Z=0 m 表示切眼位置，可以看出，工作面上方一定范围内温度低于两侧温度，相同区域内Z=0 m 截面的温度低于Z=300 m 截面的温度。

图2 垂直方向压力梯度条件下温度场分布

为了更清晰地展示沿回采方向不同位置处工作面上方的温度变化规律，选取工作面上方附近与Z轴平行的一条测线，其端点坐标为（300, 160, 0）和（300, 160, 600），测线上的温度变化见图3。可以看出，水压由0.5 MPa 增加至3 MPa 过程中，在切眼处（Z=0 m）温度从29.2 ℃降至25.8 ℃，在Z=600 m 处温度从29.4 ℃降至26.5 ℃，在Z 轴方向上距离切眼越远，工作面上方相同区域内温度越高，同时压力越大，工作面附近温度越低。产生以上现象的原因是，压力梯度的方向是垂直方向，地下水在含水层之间的联系也是主要发生在垂直方向，工作面上部区域的地下水补给主要来自上部含水层，上部含水层由于埋深较浅，地温较低，地下水由含水层流入工作面的过程中会与周围岩体发生热量交换，降低工作面附近的温度。由上述现象可以得出，若工作面附近温度降低，则反映了地下水来自赋存较浅的含水层，上覆浅部含水层和深部含水层之间发生水力联系。同时，对工作面上部区域进行温度监测，若工作面上部同一高度位置温度发生明显的温度变化，则说明该区域的渗流场发生变化，需要加强对水害的监测预警。

图3 不同水压下测线上的温度分布

3.2 水平方向渗流速度及温度场分布

图4 展示了回采350 m 时，沿回采方向不同位置处、不同水压下温度场的分布。沿回采方向同样选取上述3 个截面，分别为Z=0 m、Z=300 m 和Z=600 m。可以看出，3 个截面表示的区域温度分布相同，没有变化。同时，在不同水压和回采方向上，测线上的温度也没有发生变化。产生这种现象的原因是，压力梯度的方向是水平方向，地下水在含水层之间没有水力联系，水力联系主要发生在同一岩层内，地下水补给主要来自同一含水层，同一含水层在水平方向上由于所处的深度相同，温度也相同，地下水与其在水平方向上流动所经过的区域温度相同，不会发生热量交换，所以不会使流经区域的温度升高或者降低，因而，温度场不会发生变化。由上述现象可以得出，含水层之间没有水力联系，水力联系只发生在水平方向，不会引起温度场的变化。

图4 水平方向压力梯度条件下温度场分布

4 结论

1）温度场-渗流场耦合数值模拟模型能够模拟回采方向上不同位置、不同含水层水力联系情况下的温度场分布，温度场的变化规律能够用来判断含水层水力联系规律。

2）当垂直方向上含水层之间发生水力联系时，工作面附近温度会降低，且距离开切眼越远，温度越高；当地下水在含水层之间没有水力联系时，在同一含水层不同区域之间有水平方向的水力联系时，不会引起工作面附近温度场的变化。

3）工作面附近温度场的变化能够反映含水层之间的水力联系情况，可以用来判断顶板含水层不同区域之间以及不同含水层之间是否发生水力联系，进而可以用来预防矿井水害。