非承压含水层定水头抽水两区井流数值模型研究

2023-11-06史鹏钰宗一杰滕开庆刘健军

煤田地质与勘探 2023年10期

史鹏钰，宗一杰，滕开庆，刘健军，肖良

(1.广西大学土木建筑工程学院，广西南宁 530004；2.广西防灾减灾与工程安全重点实验室，广西南宁 530004；3.广西大学工程防灾与结构安全教育部重点实验室，广西南宁 530004)

定水头抽水试验是一种特殊的抽水试验，其井中水头在抽水时需保持稳定[1]，常被应用于渗透系数较低的地层以及抽水量难以维持稳定的工况[2]。定水头抽水试验对于确定水力参数或解决垃圾填埋场引起的环境破坏问题具有重要意义[2-3]。自1950 年以来，水文地质学家一直在研究由定水头抽水试验引起的抽水井流问题，C.E.Jacob 等[4]第一个提出了用于研究定水头抽水试验期间的达西流流动特性的解析解。此后，水文学家们在定水头抽水井流的理论和物理模型方面进行了大量有价值的工作[2,5-8]。

一般情况下，抽水引起的地下水流动遵循达西定律，符合线性流动假设[9]。然而，许多学者发现，当水力梯度较大或较小时，流量与水力梯度会呈非线性变化关系。这种非线性抽水井流被定义为非达西流[10-11]，通常由Izbash 方程或Forchheimer 方程描述，并根据现场条件在2 个方程之间进行选择[12-13]。在过去10 年中，大量学者基于上述2 种方法，通过原位试验研究了定水头抽水引起的非达西流流动特性问题[5,8,12,14-15]。例如，P.M.Quinn 等[15]在加拿大安大略省圭尔夫市附近的裂缝性白云岩和砂岩中进行了定水头原位试验，并指出使用线性特征的达西流来拟合实验数据可能会产生很大的误差。Dan Hancheng 等[6]基于定水头试验的降深数据，对无黏结级配骨料的参数进行评估，发现即使在低水力梯度下非达西流现象也很普遍。Liu Zhongyu 等[8]使用改进的渗透率测试装置在饱和黏土含水层中进行定水头测试，发现抽水试验会导致水流具有明显的非线性特性。这些研究认为，忽略非达西效应会使抽水井流模拟和确定含水层参数产生较大的误差。然而，现有的定水头抽水模型大多是基于现场或室内试验数据拟合得出，适用地质范围较为狭窄，且不具备理论依据。到目前为止，关于定水头抽水引起的非达西流数值模型研究很少，提出可通用于常见含水层系统的相关数值模型是具备现实价值的。

由文献评述笔者发现，现有定水头抽水试验研究通常基于抽水含水层是承压层的假设上提出的。然而，世界上许多抽水井都处于非承压含水层中[16-17]，并且关于非承压含水层定水头抽水井流的研究一直较为缺乏。到目前为止，现有的非承压含水层抽水井流模型主要分为考虑和不考虑垂向流分量2 类。在考虑垂直流分量模型中，重力释水不是瞬间释放的，这进一步导致水头上方的非饱和带会对下方饱和带进行补给[18]。基于含水层重力储水量远大于弹性储水量的事实，在抽水井内降深较小或抽水时间足够长的定水头试验中，垂直流的影响有限，可以合理使用不考虑垂向流分量的Dupuit 假设，该假设认为重力储存为瞬时释放[19-21]。根据Dupuit 的调查，当非承压含水层降深变化较缓时，可以忽略垂直流的影响，并假设所有补给均来自水平方向。因此，它可用于简化无约束流动问题，使其在解析上易于处理。在非承压含水层中，当水位降深相对初始水位较小时，可以优先考虑使用Dupuit 假设来模拟流动，以满足实际工程的需要[2,13,20,22-24]。

从理论上说，由于抽水井附近的流速足够大，会对相应区域的流态造成影响。足够大的流速可以使雷诺数大于临界雷诺数，从而导致抽水井附近水流呈现非达西流态。而随着径向距离的增加，雷诺数随流速逐渐减小。因此，它可以进一步引起非达西流向达西流的转变，导致含水层中呈现出两区流特性[25-29]。迄今为止，在非承压含水层中通过定水头试验引发的这种两区流动机制仍不清楚，相关研究不足。

针对上述问题，笔者提出一种非承压含水层定水头抽水试验导致的非达西-达西两区渗流数值模型。当含水层储水系数较小时，通过Izbash 方程描述比流量和水力梯度之间的关系，以进行数学建模。基于Dupuit 假设，利用有限差分法推导出易于实际操作的简洁数值模型。所提出模型通过与COMSOL Multiphysics 的数值解进行比较，验证渗流模型的可靠性，以期为两区流特性分析奠定基础。

1 两区流数学模型

图1 为在非承压含水层中定水头抽水试验导致的两区渗流模型示意图，模型考虑了非达西和达西2 个区域的水流流动，其数学模型的假设包括：(1) 该含水层是非承压、均质、各向同性的非承压含水层，初始水头为一定值；(2) 抽水井完全穿透非承压含水层，抽水井内的水头自抽水开始后保持不变；(3) 非达西区域(简称N 区)为有限区域，主要发生在抽水井周围；达西区域(简称D 区)为半无限区域，位于N 区外围；(4) 抽水井半径与有效半径相同，为一个固定正数[23]。

图1 考虑两区转换的非承压含水层定水头抽水模型Fig.1 Schematic of the constant-head test in an unconfined aquifer considering the two-region transform

基于质量守恒定律和Dupuit 假设，根据Boussinesq方程将N 区和D 区中的抽水流量分别用下式描述。

式中：q1和q2分别为N 区和D 区中的比流量，m/s；h1和h2分别为N 区和D 区中的水头，m；上述参数均为时间t和水平距离r(指某一点与抽水井中心的水平距离)的函数；S1和S2分别为N 区和D 区的储水系数。

假设N 区和D 区之间的转换界面位于r=R处，根据流体连续性假设得到边界处的边界条件如下：

达西定律描述了水力坡度(I)和比流量(q)之间的线性关系。对于多孔介质流，研究表明其线性关系在以下2 种情况下是不成立的。一种是在较低水力梯度下通过低渗介质的流动，另一种则是在较高水力梯度下通过粗粒高渗介质的流动。J.Bear[30]提到，基于平均直径的雷诺数(Re)在Re<1∼10 时，达西定律是适用的。在此范围内，渗流流动可以被归纳为层流，称为达西流，否则即可被称为非达西流。此外，Re较高时(Re>10)的流动为湍流。对于抽水井附近的流动，水力梯度通常很高，流速随之增大[31]。

1.1 非达西区域

Izbash 方程可用于描述非达西流的I和q之间的关系。使用该非线性函数的一个假设是，至少在某些流动阶段，所涉及的非达西系数n和准渗透系数kc与空间和时间无关[32]。其公式表示为：I=cqn,c为常数，当参数 1

当Re>10 时，可认为该区域处于N 区，其抽水井流为非达西渗流[25]，N 区的比流量可用Izbash 方程表示。

式中：K1为N 区中的准渗透系数，(m/s)n1(n1条件下)；n1为N 区中的非达西湍流系数。当n1<1 时渗流被称为前达西流；当 1

模型的初始条件是：

在定水头抽水试验中，代表完整井的第一类边界条件如下：

式中：Q为抽水流量，m3/s。

将式(5)代入式(1)可得：

已知N 区中的流量等于比流量乘以过流断面面积，如下式所述：

式中：A=2πrh1为井壁处过流断面面积，m2；在抽水井中h1为已知参数。

假设非承压含水层的降深远远小于b，则可以基于Dupuit 假设，在式(8)和式(9)中近似引入来线性化一小部分推导过程。研究表明采用该近似方法可使建模过程大为简化，同时该近似方法所引起的计算误差相对较小，可以忽略不计[34-35]。因此，式(9)和式(10)可以近似改写为：

结合式(11)与式(12)，N 区中的控制方程可表示如下：

1.2 达西区域

当 1

式中：K2为D 区中的渗透系数，m/s。

初始条件为：

无穷远处的边界条件为：

将式(14)代入式(2)可得：

同理，为了简化建模过程，假设非承压含水层的降深远远小于b，基于Dupuit 假设，引入近似关系来线性化式(17)，可得D 区中的控制方程可表示如下：

2 基于有限差分的数值解

有限差分法简称FDM，常用于求解大型非线性地下水流运动问题，FDM 在不规则含水层边界和含水层内区域参数的紧密空间近似方面具有灵活性[36]。因此，本研究采用FDM 基本思想按时间步长和空间步长对定解区域进行网格划分[37]，把从无穷远处到抽水井处的空间区域及从抽水开始到抽水结束的时间区域用有限个网格节点代替。

采用FDM 推求上述方程的数值解，需先对抽水井流模型进行如下两步预处理。第一步，针对无穷远边界，本文拟采用恒定不变且取值较大的半径值re来代替，使抽水井流模型可被认为是一个因变量为半径r和时间t的一维非稳定流模型。第二步，将时间区间[0,t] 离散为N个步长，径向空间 [rw,re] 离散为J个步长。其中，任何一个时间坐标上节点m处的时间值tm必满足 0 ≤tm≤t,m=0,1,···,N，任何一个径向坐标上节点j处的径向距离rj需满足rw≤rj≤re,j=0,1,···,J。其步长大小的具体值表示为：Δt=t/N，Δr=(re-rw)/J。由于在抽水井附近水位降深变化较大，因此在离抽水井较近的地方空间步长应尽可能小，而离抽水井较远的地方空间步长可以适当加大。对N 区和D 区的控制方程及其所有边界条件进行离散化，采用向后差分的方法，依次得出N 区和D 区的水头。

2.1 非达西区域

N 区控制方程式(13)为一个二阶齐次微分方程，基于FDM 的基本理论，将控制方程离散化，得到离散化后的控制方程如下：

其中：

对于N 区的边界条件式(7)和式(8)，离散化后的边界条件如下：

2.2 达西区域

同理，D 区控制方程式(18)为一个二阶齐次微分方程，基于FDM 的基本理论离散化后可得：

其中：

对于D 区的边界条件式(15)和式(16)，离散化后的边界条件如下：

利用Matlab 软件，将式(19)-式(24)与式(25)-式(30)分别导入，输入合适的步长，即可由软件自动计算并分别得出N 区及D 区各个节点上的流量及水头值。通过在rj≤R的径向节点上利用式(19)以及在rj>R的径向节点上利用式(25)计算得出的结果，调用Matlab 中名为Equations and System Solver 的内置模块，即可得到在假设的定水头试验条件下任意点的流量及水头情况，N 区与D 区求解矩阵分别如下：

定水头抽水试验可以使抽水井附近的水头在抽水早期迅速下降，并在短时间内接近hw。在此之后，水头变化的幅度会大大减小。因此，除了在抽水刚开始时的短时间内，其余阶段基于Dupuit 假设提出的解是可以接受的。

在定水头抽水试验中，由于流动的湍流程度会随着与泵井径向距离的增加而减弱，因此可将N 区定义为Re大于临界雷诺数Rec的区域，并且通过判断Re=Rec|r→R来确定该时刻的R值。此外，还应该注意的是，所提出的解也可通过假设R=rw来评估纯达西流或设定R=∞ 调查纯非达西流的降深。计算雷诺数的公式[25]为Re=qd/υ，其中d为介质的特征直径，m；υ 为水的运动黏度，m2/s，其值通常为 1 ×106m2/s。可以得到达西区和非达西区转换界面处的Rec表达式如下：

根据相关研究，通常认为Rec=10。基于计算得到的水头h值，可以通过Matlab 进而推求两区流界面位置R值。拟通过计算得出的达西区水头h2进行迭代计算，其具体步骤如下：

(1) 选取初始值。假设在抽水前，井流为达西流，即tm=0 时，时间步长数m=0，两区流交界面位置为R=Rm=0。

(2) 计算流速。计算第m+1 个时间步长tm+1时的水头并计算每个节点处的流速及雷诺数Re，然后根据临界雷诺数Rec确定两区流交界面位置Rmc。

(3) 收敛性判断。δ为一个给定的判别标准，若 |Rm-Rmc|<δ，则认为第m个时间步长的两区流交界面位置R=Rm。然后计算下一时间步长两区流交界面位置：令m=m+1，若mδ，令Rm=Rmc。

(4) 在新交界面处计算流速。根据计算得出的水头值，计算每个时间节点处的流速与雷诺数，由Rec确定两区流交界面位置Rmc，返回步骤(3)进行判断，进入下一步迭代计算，直至满足δ<5% 的精度要求，迭代终止。

3 模型验证与对比

鉴于目前缺少非承压含水层定水头抽水计算模型的相关研究，难以找到相关工况完全一致的模型研究或实测案例，为验证本文提出有限差分数值模型的可靠性，本章将提出的非承压定水头抽水井流模型退化为承压定水头达西抽水井流模型，并与现有的经典承压定水头达西抽水井流解析模型(Jacob 和Lohman 模型)进行对比分析[4]。同时，通过假设案例，对所提出有限差分模型与基于COMSOL 有限元模型的计算结果进行对比分析，以验证所提出的有限差分模型的可靠性。为方便计算，本文基于Matlab 软件实现对提出模型的自动化求解。

3.1 与Jacob 和Hantush 解析解比较

在抽水过程中，通过控制含水层压力水头始终大于含水层厚度，将模型退化成承压含水层定水头抽水井流模型。然后，再令n1=1 将模型退化成达西流，采用Jacob 和Lohman 模型的参数，将比流量q及抽水时间t无量纲处理并代入本模型中，可得到对比结果如图2 所示。从图中可以看出，本文在定水头抽水的承压含水层达西流情况下的有限差分解与Jacob 和Lohman 模型解高度吻合，说明在模型推导过程中由数值差分带来的计算误差可忽略不计[38]。

图2 有限差分解与Jacob 和Lohman 解对比Fig.2 Comparison between the finite difference solution with Jacob and Lohman solution

3.2 与COMSOL 有限元解比较

非承压含水层非达西定水头抽水径流模型采用有限元模型基于COMSOL Multiphysics 软件建立，该软件已被证明是最可靠的多物理场仿真软件之一[39]。COMSOL 模型是基于N 区的式(6)-式(9)及D 区的式(15)-式(17)推导而来。

需要注意的是，式(16)中无穷远边界的距离在COMSOL 解中用r=1 000 m 表示，假设非承压含水层是砂质介质，基于现有研究可知非承压含水层的储水系数为 0 .01∼0.03，非承压含水层中细砂的渗透系数为 1∼5 m/d[20]。因此，案例A 参数选取见表1。考虑到Dupuit 假设的要求，抽水井内部降深需要取一个相对较小的值，故取b=12 m,hw=9 m。表1 中的案例B、C 和D 用于下文讨论部分的参数选取。

表1 假设案例的参数值Table 1 Parameters of each hypothetical cases

图3 展示了两区流(R为变量)、纯达西流(R=rw=0.2 m) 和纯非达西流(R=∞)的水头-时间曲线对比。纯非达西流COMSOL 解近似为R=1 000 m。应该注意到，K1和S1用于纯非达西流动的模拟，而K2和S2用于纯达西流动的模拟。这一对比清楚地表明，在所有3 种流动条件下，模拟结果与有限差分解之间的偏差都很小。在抽水开始后第一个小时内出现一定偏差，可认为是忽略了所提出的有限差分解的垂直流动而引起，这种偏差随着时间的推移逐渐消失，模拟结果进一步验证了所提出的有限差分模型用于两区流动的水头模拟的有效性。

图3 不同R 值下有限差分解与数值解的水头-时间曲线对比Fig.3 Comparison of head-time curves between finite difference solution and numerical solution for different R

此外，两区流的水头明显小于纯非达西流的水头，而大于纯达西流的水头(图3)，说明非达西流对水头下降有负影响。这是因为R值越大表示非达西范围越大，抽水井附近的非承压含水层湍流范围更广。相比于纯达西流，纯非达西流的湍流程度更高。紊流会导致过水断面流速变大，相应比流量变大，含水层从外界得到水的能力增强，无穷远处对抽水井附近含水层补给变大，进而延缓因抽水导致的含水层水位下降[40]。因此，在抽水后期，纯非达西流的水头能够明显大于纯达西流的水头。

图4 展示了基于不同模型R-t曲线的对比。结果表明，不同时间下本文提出的有限差分模型解与COMSOL 有限元模型解拟合程度较高，验证了通过迭代法计算出的R值的可靠性。而R值随时间推移不断减小，说明非达西流对两区流动的影响随着抽水的继续而减弱。在抽水开始后，非承压含水层水头逐渐下降，接近于hw，它进一步导致水力梯度和比流量的下降。比流量越小，流速越慢，雷诺数越小，最终导致N 区消失。当抽水时间足够长时，N 区可以完全消失，两区流可以转化为纯达西流。

图4 有限差分解界面距离-时间曲线与数值解的比较Fig.4 Interface distance-time curve of finite difference solution vs numerical solution

3.3 与实测数据比较

L.Jones 等[41]在美国威斯康星州开展了定水头抽水试验，所调查的非承压含水层由厚度为2.5 m 的含泥量低、含砂量高的浅风化冰碛层组成。抽水井有效半径为0.051 m，且完全贯穿非承压含水层，其底部由砾石充填。抽水共进行24 h，抽水井内定水头为1.5 m[42]。两组降深数据分别在径向距离为0.87、3.62 m 的观测井中采集。根据C.S.Chen 等[2]的研究，水平方向的给水度和渗透系数分别为 0.014∼0.042 和0.241∼0.362 m/d。因此，基于式(19)、式(25)经对比分析，选取参数S1=S2=0.031，K1=0.285 m/d，K2=0.320 m/d，n1=1.3时，可得到实测数据与本文模型解的对比曲线拟合度良好，如图5 所示。但与图3 的结果类似，有限差分解模拟曲线在抽水期间与实测数据的偏差很小，仅在抽水中期会出现一定程度的误差，且这种偏差随着时间的增加逐渐消失。图5 进一步验证了所提出的有限差分模型在实际工程中的适用性。

图5 有限差分解与实测数据对比Fig.5 Finite difference solution vs measured data

4 讨论

在实际抽水试验中，一般认为抽水附近区域流态为达西流，以方便概化模型及预测井内降深。本研究提出的定水头抽水两区流数值解，将抽水井附近的有限区域视为非达西区域，远离抽水井的半无限区域视为达西区域，从而将非承压含水层的地下水渗流根据流态区分为2 个流域。相比于纯达西流模型，该两区流模型更符合实际工况，计算结果精度更高，但相应计算过程也较复杂。该模型适用于沙土或碎石比例较大、渗透性较强的含水层，其渗流流速相对较快，Re相对较高，在抽水井附近极易形成非达西流。但在黏土等渗透能力弱的含水层，渗流流速较小，非达西范围相对小，该模型将不再适用。在实际工程应用中，通常会着重观测井内水头及抽水量的动态变化。因此，该模型有必要讨论偏离达西程度对流态的影响，进而判断不同非达西湍流系数n1对水头及抽水井抽水量的影响，并讨论井内水头hw对抽水井抽水速率的影响。准确评估抽水流量Q的动态变化对于设计定水头抽水试验至关重要。在下文中，基于假设案例B、C 和D(参数见表1)揭示非达西流的n1、两区流转换界面位置R和井内水头hw对抽水速率动态变化的影响，以及径向距离对断面流量的影响。

图6 展示了案例B 中不同时刻hw的Q-t曲线的模拟结果。结果表明，固定时刻下Q随着hw的增加不断降低，在抽水早期不同曲线对应的Q之间偏差较大，随着抽水的继续，流量偏差逐渐减小。Q的变化与hw呈负相关关系，在抽水早期Q-t曲线下降速度较快，随着抽水的继续Q下降的趋势越来越平缓，符合hw的小幅降低会使流速大幅度升高的事实。流量变化实际上是由抽水井内部水头下降而引起的水力梯度变化造成的，因此，较大的hw意味着较小的比流量，从而导致抽水流量较小。

图6 泵井内不同水头下流量的动态变化Fig.6 Dynamic development of the flow rates with different hydraulic head inside the pumping well

图7 展示了案例C 中非达西湍流系数对泵送流量的影响，其中图7a 和图7b 分别为固定径向距离r=2 m 时的Q-t曲线和相同抽水时间t=1 d 时的h-r曲线。从图7a 中可以发现，在抽水早期，固定时间点的Q随着n1的增加而增加。而随着抽水的继续，不同n1对应的流量偏差不断减小。从理论上讲，这是可以理解的。抽水井附近含水层补给水量以重力释水为主，如式(10)所示，假设其等于比流量乘以过流断面面积，比流量越大，抽水流量越大。由此可见，在抽水早期，抽水井附近区域含水层的重力储水逐渐释放并补充到抽水井内，n1越大，湍流程度越大，导致非承压含水层储水的释放越快。随着抽水的继续，N 区重力储水的释放逐渐趋于稳定，抽水井中的水主要通过D 区的水流进行补给，导致Q-t曲线逐渐趋同。

图7 不同非达西湍流系数对泵送流量的影响Fig.7 Effects of different non-Darcy coefficients on pumping flow

图7b 体现了t=1 d 时不同径向距离的非达西湍流的发展特征。结果表明，固定r下h随着n1的增大而增大，在r较小时不同曲线对应的h之间的偏差较大，随着r的增大，h偏差逐渐减小。非达西系数与水头之间存在正相关关系，在r较小时h-r曲线上升的速度较快，随着r的增大，曲线斜率逐渐减小。n1越大，表示水力梯度越大，对应流量越大。在抽水井附近，水头波动较大，但没有出现不合理的突变，证明了所提有限差分解的适用性。

图8 展示了径向距离对断面流量的影响，其中图8a 和图8b 分别用于表示抽水前期 (t=0.1 d) 和抽水后期 (t=100.0 d) 的Q1-r曲线(Q1为断面流量)，其水力参数见表1 的案例D。对比图8a 和图8b 可以发现，t=0.1 d 时Q1随r的增加而快速下降，其Q1较t=100.0 d 时更大。随着抽水时间的增加，R逐渐变小，该结论符合图4 中R-t曲线的趋势。Q1随着r的增大不断减小，在N 区部分Q1曲线随距离下降的速度较快，并在经过转换界面R后出现突变，在D 区Q1曲线随距离下降的速率明显变慢。从理论上来说，定水头抽水过程中，在井壁处存在Q1变化的最大值，随着r的增大，距离抽水井越远，对抽水井补给量越小，断面流量越小。在抽水前中期，含水层中距离抽水井较远的部分无法及时对抽水井附近进行充足的补给，抽水井附近区域含水层的重力储水会逐渐释放并补充到抽水井内。随着抽水时间的增长，在抽水后期，远处的水流逐渐补充到抽水井附近，抽水井周边区域的水位和断面流量的下降速度逐渐变小。相比于N 区，D 区水力梯度更小，其比流量随径向距离的下降速率更慢，其Q1-r曲线的转折点出现在r=R处。