1、非承压含水层定水头抽水两区井流数值模型研究史鹏钰1,2,3，宗一杰1,2,3，滕开庆1,2,3，刘健军1,2,3，肖良1,2,3(1.广西大学 土木建筑工程学院，广西 南宁 530004；2.广西防灾减灾与工程安全重点实验室，广西 南宁 530004；3.广西大学 工程防灾与结构安全教育部重点实验室，广西 南宁 530004)摘要:为了揭示在非承压含水层中定水头抽水试验引起的达西非达西两区流动机理，提出基于有限差分法的地下水定水头抽水井流数值模型。该模型根据抽水的流态特征将含水层分为 2 个区域：靠近抽水井的有限非达西渗流区域和远离抽水井的半无限达西渗流区域，其中非达西流区域流态的模拟基于 I
2、zbash 方程实现。通过与 COMSOL Multiphysics 的有限元数值解进行比较，验证了所提出数值解的可靠性。最后，研究有限非达西流效应对水头和抽水井抽水速率的影响以及井内水头对抽水井抽水速率的影响。研究表明：在抽水试验中非达西区域的影响不可忽略，湍流会分别导致两区流中水头较纯非达西流和纯达西流的水头偏大和偏小，且随抽水时间的增加逐渐变大；通过减小抽水井井内水头或增大非达西系数可提高抽水速率，但该影响会随抽水时间的增加而逐渐减弱；断面流量随径向距离的增大而不断减小，断面流量与径向距离曲线下降速率不断减小，且在转换界面处会出现转折点。该模型为定量研究在非达西流和达西流耦合作用下抽水井
3、附近的井流水头特征提供了一种简洁的方法，并为调查定水头抽水测试期间的抽水速率提供理论依据。关键词：非承压含水层；定水头抽水；两区流动；非达西效应；有限差分解；数值模型中图分类号：P345 文献标志码：A 文章编号：1001-1986(2023)10-0124-10Numericalmodeloftwo-regionflowbyconstant-headpumpinginanunconfinedaquiferSHI Pengyu1,2,3,ZONG Yijie1,2,3,TENG Kaiqing1,2,3,LIU Jianjun1,2,3,XIAO Liang1,2,3(1.College of
4、 Civil Engineering and Architecture,Guangxi University,Nanning 530004,China;2.Guangxi Key Laboratory ofDisaster Prevention and Engineering Safety,Nanning 530004,China;3.Key Laboratory of Disaster Prevention andStructural Safety of Ministry of Education,Guangxi University,Nanning 530004,China)Abstrac
5、t:In order to reveal the mechanism of two-region Darcian and non-Darcian flow induced by constant-headpumping tests in unconfined aquifers,a numerical model of well flow by constant-head pumping underground based onfinite difference solution was proposed.In the model,it is assumed that the aquifer i
6、s divided into two regions accordingto the pumping flow characteristics:the finite non-Darcian flow region near the pumping well and the semi-infinite Dar-cian flow region far away from the pumping well.Specifically,the flow regime in the non-Darcian flow region was sim-ulated by the Izbashs equatio
7、n,and the reliability of the proposed solution was verified by comparison with finite ele-ment numerical solutions by COMSOL Multiphysics.Finally,the influence of finite non-Darcian effect on hydraulichead and pumping rate,as well as the influence of hydraulic head in the pumping well on pumping rat
8、e,was especiallystudied.The results show that the influence of non-Darcian flow region in the pumping test cannot be ignored.The tur-bulent flow makes the hydraulic head of the two-region flow larger than that of the pure non-Darcian flow and smallerthan that of pure Darcian flow.Besides,such differ
9、ence of the hydraulic head is increased with the pumping time.Thepumping rate can be increased by reducing the hydraulic head in the pumping well or increasing the non-Darcian coeffi-cient,but the effect is gradually decreased as the pumping continues.The flow rate at the cross-section is decreased
10、withthe increase of radial distance,the gradient of the flow rate-radial distance curve is decreased with time,and a turning 收稿日期:2022-12-14；修回日期:2023-07-02基金项目:国家自然科学基金青年基金项目(41807197)；广西自然科学基金项目(2018GXNSFAA138042)；河北高层次人才资助项目(B2018003016)第一作者:史鹏钰，1997 年生，男，河南周口人，硕士，从事水文地质与地下水动力学研究.E-mail：通信作者:肖良，1
11、985 年生，男，广西崇左人，博士，副教授，从事水文地质与地下水动力学研究.E-mail： 第 51 卷 第 10 期煤田地质与勘探Vol.51 No.102023 年 10 月COAL GEOLOGY&EXPLORATIONOct.2023史鹏钰，宗一杰，滕开庆，等.非承压含水层定水头抽水两区井流数值模型研究J.煤田地质与勘探，2023，51(10)：124133.doi:10.12363/issn.1001-1986.22.12.0944SHI Pengyu，ZONG Yijie，TENG Kaiqing，et al.Numerical model of two-region flow b
12、y constant-head pumping in an uncon-fined aquiferJ.Coal Geology&Exploration，2023，51(10)：124133.doi:10.12363/issn.1001-1986.22.12.0944point will appear at the conversion interface between the two regions.The proposed model provides a simple method forquantitative studies on the characteristics of hyd
13、raulic head near the pumping well under the coupling of non-Darcianand Darcian effects,and provides a theoretical basis for determining the pumping rate during constant-head pumpingtest.Keywords:unconfined aquifer;constant-head pumping;two-region flow;non-Darcian effect;finite difference solution;nu
14、merical model 定水头抽水试验是一种特殊的抽水试验，其井中水头在抽水时需保持稳定1，常被应用于渗透系数较低的地层以及抽水量难以维持稳定的工况2。定水头抽水试验对于确定水力参数或解决垃圾填埋场引起的环境破坏问题具有重要意义2-3。自 1950 年以来，水文地质学家一直在研究由定水头抽水试验引起的抽水井流问题，C.E.Jacob 等4第一个提出了用于研究定水头抽水试验期间的达西流流动特性的解析解。此后，水文学家们在定水头抽水井流的理论和物理模型方面进行了大量有价值的工作2,5-8。一般情况下，抽水引起的地下水流动遵循达西定律，符合线性流动假设9。然而，许多学者发现，当水力梯度较大或较小
15、时，流量与水力梯度会呈非线性变化关系。这种非线性抽水井流被定义为非达西流10-11，通常由 Izbash 方程或 Forchheimer 方程描述，并根据现场条件在 2 个方程之间进行选择12-13。在过去10 年中，大量学者基于上述 2 种方法，通过原位试验研究了定水头抽水引起的非达西流流动特性问题5,8,12,14-15。例如，P.M.Quinn 等15在加拿大安大略省圭尔夫市附近的裂缝性白云岩和砂岩中进行了定水头原位试验，并指出使用线性特征的达西流来拟合实验数据可能会产生很大的误差。Dan Hancheng 等6基于定水头试验的降深数据，对无黏结级配骨料的参数进行评估，发现即使在低水力梯
16、度下非达西流现象也很普遍。Liu Zhongyu 等8使用改进的渗透率测试装置在饱和黏土含水层中进行定水头测试，发现抽水试验会导致水流具有明显的非线性特性。这些研究认为，忽略非达西效应会使抽水井流模拟和确定含水层参数产生较大的误差。然而，现有的定水头抽水模型大多是基于现场或室内试验数据拟合得出，适用地质范围较为狭窄，且不具备理论依据。到目前为止，关于定水头抽水引起的非达西流数值模型研究很少，提出可通用于常见含水层系统的相关数值模型是具备现实价值的。由文献评述笔者发现，现有定水头抽水试验研究通常基于抽水含水层是承压层的假设上提出的。然而，世界上许多抽水井都处于非承压含水层中16-17，并且关于非
17、承压含水层定水头抽水井流的研究一直较为缺乏。到目前为止，现有的非承压含水层抽水井流模型主要分为考虑和不考虑垂向流分量 2 类。在考虑垂直流分量模型中，重力释水不是瞬间释放的，这进一步导致水头上方的非饱和带会对下方饱和带进行补给18。基于含水层重力储水量远大于弹性储水量的事实，在抽水井内降深较小或抽水时间足够长的定水头试验中，垂直流的影响有限，可以合理使用不考虑垂向流分量的 Dupuit 假设，该假设认为重力储存为瞬时释放19-21。根据 Dupuit 的调查，当非承压含水层降深变化较缓时，可以忽略垂直流的影响，并假设所有补给均来自水平方向。因此，它可用于简化无约束流动问题，使其在解析上易于处理
18、。在非承压含水层中，当水位降深相对初始水位较小时，可以优先考虑使用 Dupuit 假设来模拟流动，以满足实际工程的需要2,13,20,22-24。从理论上说，由于抽水井附近的流速足够大，会对相应区域的流态造成影响。足够大的流速可以使雷诺数大于临界雷诺数，从而导致抽水井附近水流呈现非达西流态。而随着径向距离的增加，雷诺数随流速逐渐减小。因此，它可以进一步引起非达西流向达西流的转变，导致含水层中呈现出两区流特性25-29。迄今为止，在非承压含水层中通过定水头试验引发的这种两区流动机制仍不清楚，相关研究不足。针对上述问题，笔者提出一种非承压含水层定水头抽水试验导致的非达西达西两区渗流数值模型。当含水
19、层储水系数较小时，通过 Izbash 方程描述比流量和水力梯度之间的关系，以进行数学建模。基于Dupuit 假设，利用有限差分法推导出易于实际操作的简洁数值模型。所提出模型通过与 COMSOL Mul-tiphysics 的数值解进行比较，验证渗流模型的可靠性，以期为两区流特性分析奠定基础。1两区流数学模型图 1 为在非承压含水层中定水头抽水试验导致的两区渗流模型示意图，模型考虑了非达西和达西 2 个区域的水流流动，其数学模型的假设包括：(1)该含水层是非承压、均质、各向同性的非承压含水层，初始水头为一定值；(2)抽水井完全穿透非承压含水层，抽水井内的水头自抽水开始后保持不变；(3)非达西区域
20、(简称 N 区)为有限区域，主要发生在抽水井周围；达西区域(简称 D 区)为半无限区域，位于 N 区外围；(4)抽水井半径与有效半径相同，为一个固定正数23。第 10 期史鹏钰等：非承压含水层定水头抽水两区井流数值模型研究 125 基于质量守恒定律和 Dupuit 假设，根据 Boussinesq方程将 N 区和 D 区中的抽水流量分别用下式描述。q1r+q1r=S1h1h1t(1)q2r+q2r=S2h2h2t(2)q1q2h1h2trS1S2式中：和 分别为 N 区和 D 区中的比流量，m/s；和 分别为 N 区和 D 区中的水头，m；上述参数均为时间 和水平距离(指某一点与抽水井中心的水
21、平距离)的函数；和 分别为 N 区和 D 区的储水系数。r=R假设 N 区和 D 区之间的转换界面位于 处，根据流体连续性假设得到边界处的边界条件如下：limrRq1=limrRq2(3)h1|rR=h2|rR(4)IqReRe 10)的流动为湍流。对于抽水井附近的流动，水力梯度通常很高，流速随之增大31。1.1非达西区域IqnkcI=cqnc1 n 2Izbash 方程可用于描述非达西流的 和 之间的关系。使用该非线性函数的一个假设是，至少在某些流动阶段，所涉及的非达西系数 和准渗透系数 与空间和时间无关32。其公式表示为：,为常数，当参数 时，该公式用于表征因显著惯性效应n=10 n 10
22、当 时，可认为该区域处于 N 区，其抽水井流为非达西渗流25，N 区的比流量可用Izbash 方程表示。q1=(K1h1r)1n1(5)K1(m/s)n1n1n1 11 n1 2n1=1K1式中：为 N 区中的准渗透系数，(n1条件下)；为 N 区中的非达西湍流系数。当 时渗流被称为前达西流；当 时为后达西流；当，式(5)即为达西定律，此时 为渗透系数10。模型的初始条件是：h1(r,0)=b(6)在定水头抽水试验中，代表完整井的第一类边界条件如下：limrrwh1=hw(7)limrrw2rhw(K1h1r)1n1=Q(8)Qm3/s式中：为抽水流量，。将式(5)代入式(1)可得：2h1r2
23、+n1rh1r=S1h1n1K11n1(h1r)n11n1h1t(9)已知 N 区中的流量等于比流量乘以过流断面面积，如下式所述：Q Aq1=2rh1(K1h1r)1n1(10)A=2rh1m2h1式中：为井壁处过流断面面积，；在抽水井中 为已知参数。bS1h1n1K11n1S1bn1K11n1假设非承压含水层的降深远远小于，则可以基于 Dupuit 假设，在式(8)和式(9)中近似引入 来线性化一小部分推导过程。研究表明采用该近似方法可使建模过程大为简化，同时该近似方法所引起的计算误差相对较小，可以忽略不计34-35。因此，式(9)和式(10)可以近似改写为：2h1r2+n1rh1r=S1b
24、n1K11n1(h1r)n11n1h1t(11)h1r(Q2rb)n1K1(12)结合式(11)与式(12)，N 区中的控制方程可表示 定水头抽水井地表达西区非达西区不透水层初始水头bRrwhwrw抽水井有效半径，m；hw抽水井中水头，m；RN 区和 D区转换界面到抽水井中心的距离，m；b含水层厚度，m图 1 考虑两区转换的非承压含水层定水头抽水模型Fig.1 Schematic of the constant-head test in an unconfinedaquifer considering the two-region transform 126 煤田地质与勘探第 51 卷如下：2
25、h1r2+n1rh1r=1r1n1h1t(13)1=S1bn1K1(Q2b)n11其中：。1.2达西区域1 Re R 的径向节点上利用式(25)计算得出的结果，调用 Matlab 中名为 Equations and System Solver 的内置模块，即可得到在假设的定水头试验条件下任意点的流量及水头情况，N 区与 D 区求解矩阵分别如下：(hm+110hm+111hm+112 hm+11Jhm+11J+1)A1100B11A120C11B1200C120.00A1J00A1J00A1J=(hm10hm11 hm1J)(31)(hm+120hm+121hm+122 hm+12Jhm+12J
26、+1)A2100B21A220C21B2200C220.00A2J00A2J00A2J=(hm20hm21 hm2J)(32)hw定水头抽水试验可以使抽水井附近的水头在抽水早期迅速下降，并在短时间内接近。在此之后，水头变化的幅度会大大减小。因此，除了在抽水刚开始时的短时间内，其余阶段基于 Dupuit 假设提出的解是可以接受的。ReRecRe=Rec|rRRR=rwR=Re=qd/dm2/s1106m2/sRec在定水头抽水试验中，由于流动的湍流程度会随着与泵井径向距离的增加而减弱，因此可将 N 区定义为 大于临界雷诺数 的区域，并且通过判断 来确定该时刻的 值。此外，还应该注意的是，所提出的
27、解也可通过假设 来评估纯达西流或设定 调查纯非达西流的降深。计算雷诺数的公式25为，其中 为介质的特征直径，m；为水的运动黏度，其值通常为。可以得到达西区和非达西区转换界面处的 表达式如下：Rec=dK2h2j+1h2jr(33)Rec=10hRh2根据相关研究，通常认为。基于计算得到的水头 值，可以通过 Matlab 进而推求两区流界面位置 值。拟通过计算得出的达西区水头 进行迭代计算，其具体步骤如下：tm=0m=0R=Rm=0(1)选取初始值。假设在抽水前，井流为达西流，即 时，时间步长数，两区流交界面位置为。m+1tm+1h2m+1jRe(2)计算流速。计算第 个时间步长 时的水头，并计
28、算每个节点处的流速及雷诺数，然RecRmc后根据临界雷诺数 确定两区流交界面位置。|RmRmc|mR=Rmm=m+1m Rm=Rmc(3)收敛性判断。为一个给定的判别标准，若，则认为第 个时间步长的两区流交界面位置。然后计算下一时间步长两区流交界面位置：令，若，返回步骤(2)计算下一时刻的流速，反之则计算终止。若，令。RecRmc 5%(4)在新交界面处计算流速。根据计算得出的水头值，计算每个时间节点处的流速与雷诺数，由 确定两区流交界面位置，返回步骤(3)进行判断，进入下一步迭代计算，直至满足 的精度要求，迭代终止。3模型验证与对比鉴于目前缺少非承压含水层定水头抽水计算模型的相关研究，难以找
29、到相关工况完全一致的模型研究或实测案例，为验证本文提出有限差分数值模型的可靠性，本章将提出的非承压定水头抽水井流模型退化为承压定水头达西抽水井流模型，并与现有的经典承压定水头达西抽水井流解析模型(Jacob 和 Lohman 模型)进行对比分析4。同时，通过假设案例，对所提出有限差分模型与基于 COMSOL 有限元模型的计算结果进行对比分析，以验证所提出的有限差分模型的可靠性。为方便计算，本文基于 Matlab 软件实现对提出模型的自动化求解。128 煤田地质与勘探第 51 卷3.1与 Jacob 和 Hantush 解析解比较n1=1qt在抽水过程中，通过控制含水层压力水头始终大于含水层厚度
30、，将模型退化成承压含水层定水头抽水井流模型。然后，再令 将模型退化成达西流，采用 Jacob 和 Lohman 模型的参数，将比流量 及抽水时间 无量纲处理并代入本模型中，可得到对比结果如图 2 所示。从图中可以看出，本文在定水头抽水的承压含水层达西流情况下的有限差分解与 Jacob 和Lohman 模型解高度吻合，说明在模型推导过程中由数值差分带来的计算误差可忽略不计38。3.2与 COMSOL 有限元解比较非承压含水层非达西定水头抽水径流模型采用有限元模型基于 COMSOL Multiphysics 软件建立，该软件已被证明是最可靠的多物理场仿真软件之一39。COMSOL 模型是基于 N
31、区的式(6)式(9)及 D 区的式(15)式(17)推导而来。r=1 000 m0.010.0315 m/db=12 m hw=9 m需要注意的是，式(16)中无穷远边界的距离在COMSOL 解中用 表示，假设非承压含水层是砂质介质，基于现有研究可知非承压含水层的储水系数为，非承压含水层中细砂的渗透系数为 20。因此，案例 A 参数选取见表 1。考虑到 Dupuit 假设的要求，抽水井内部降深需要取一个相对较小的值，故取,。表 1 中的案例B、C 和 D 用于下文讨论部分的参数选取。RR=rw=0.2 mR=R=1 000 mK1S1K2S2图 3 展示了两区流(为变量)、纯达西流()和纯非达
32、西流()的水头时间曲线对比。纯非达西流 COMSOL 解近似为。应该注意到，和 用于纯非达西流动的模拟，而 和 用于纯达西流动的模拟。这一对比清楚地表明，在所有 3 种流动条件下，模拟结果与有限差分解之间的偏差都很小。在抽水开始后第一个小时内出现一定偏差，可认为是忽略了所提出的有限差分解的垂直流动而引起，这种偏差随着时间的推移逐渐消失，模拟结果进一步验证了所提出的有限差分模型用于两区流动的水头模拟的有效性。表1假设案例的参数值Table1Parametersofeachhypotheticalcases案例hw/mK1/(md1)n1K2/(md1)n1S1S2b/mrw/mr/md/mmA9
33、1.11.11.10.010.01120.230.3B7,8,91.11.11.20.010.01100.250.3C91.11.11.0,1.2,1.40.010.01120.220.3D91.11.11.20.010.01120.220.3 R此外，两区流的水头明显小于纯非达西流的水头，而大于纯达西流的水头(图 3)，说明非达西流对水头下降有负影响。这是因为 值越大表示非达西范围越大，抽水井附近的非承压含水层湍流范围更广。相比于纯达西流，纯非达西流的湍流程度更高。紊流会导致过水断面流速变大，相应比流量变大，含水层从外界得到水的能力增强，无穷远处对抽水井附近含水层补给变大，进而延缓因抽水导致
34、的含水层水位下降40。因此，在抽水后期，纯非达西流的水头能够明显大于纯达西流的水头。R-t图 4 展示了基于不同模型 曲线的对比。结果表明，不同时间下本文提出的有限差分模型解与COMSOL 有限元模型解拟合程度较高，验证了通过迭 100101102103104105106t/d101100q/(ms1)Jacob 和 Lohman 模型有限差分模型图 2 有限差分解与 Jacob 和 Lohman 解对比Fig.2 Comparison between the finite difference solution withJacob and Lohman solution 1021011001
35、01t/d8.08.28.48.68.89.09.29.49.69.810.0h/mR=R为时变函数R=0.2 mR=R为时变函数R=0.2 m有限差分解COMSOL 解 图 3 不同 R 值下有限差分解与数值解的水头时间曲线对比Fig.3 Comparison of head-time curves between finite differencesolution and numerical solution for different R第 10 期史鹏钰等：非承压含水层定水头抽水两区井流数值模型研究 129 RRhw代法计算出的 值的可靠性。而 值随时间推移不断减小，说明非达西流对两区
36、流动的影响随着抽水的继续而减弱。在抽水开始后，非承压含水层水头逐渐下降，接近于，它进一步导致水力梯度和比流量的下降。比流量越小，流速越慢，雷诺数越小，最终导致N 区消失。当抽水时间足够长时，N 区可以完全消失，两区流可以转化为纯达西流。103102101100101t/d0.60.70.80.91.01.11.21.31.4R/mCOMSOL 解有限差分解图 4 有限差分解界面距离时间曲线与数值解的比较Fig.4 Interface distance-time curve of finite difference solution vsnumerical solution 3.3与实测数据比较
37、0.0140.0420.2410.362 m/dS1=S2=0.031 K1=0.285 m/d K2=0.320 m/dn1=1.3L.Jones 等41在美国威斯康星州开展了定水头抽水试验，所调查的非承压含水层由厚度为 2.5 m 的含泥量低、含砂量高的浅风化冰碛层组成。抽水井有效半径为 0.051 m，且完全贯穿非承压含水层，其底部由砾石充填。抽水共进行24 h，抽水井内定水头为1.5 m42。两组降深数据分别在径向距离为 0.87、3.62 m 的观测井中采集。根据 C.S.Chen 等2的研究，水平方向的给水度和渗透系数分别为 和。因此，基于式(19)、式(25)经对比分析，选取参数
38、，时，可得到实测数据与本文模型解的对比曲线拟合度良好，如图 5 所示。但与图 3 的结果类似，有限差分解模拟曲线在抽水期间与实测数据的偏差很小，仅在抽水中期会出现一定程度的误差，且这种偏差随着时间的增加逐渐消失。图 5 进一步验证了所提出的有限差分模型在实际工程中的适用性。4讨论在实际抽水试验中，一般认为抽水附近区域流态为达西流，以方便概化模型及预测井内降深。本研究提出的定水头抽水两区流数值解，将抽水井附近的有限区域视为非达西区域，远离抽水井的半无限区域视为达西区域，从而将非承压含水层的地下水渗流根据流态区分为 2 个流域。相比于纯达西流模型，该两区Ren1hwQn1Rhw流模型更符合实际工况
39、，计算结果精度更高，但相应计算过程也较复杂。该模型适用于沙土或碎石比例较大、渗透性较强的含水层，其渗流流速相对较快，相对较高，在抽水井附近极易形成非达西流。但在黏土等渗透能力弱的含水层，渗流流速较小，非达西范围相对小，该模型将不再适用。在实际工程应用中，通常会着重观测井内水头及抽水量的动态变化。因此，该模型有必要讨论偏离达西程度对流态的影响，进而判断不同非达西湍流系数 对水头及抽水井抽水量的影响，并讨论井内水头 对抽水井抽水速率的影响。准确评估抽水流量 的动态变化对于设计定水头抽水试验至关重要。在下文中，基于假设案例 B、C 和 D(参数见表 1)揭示非达西流的、两区流转换界面位置 和井内水头
40、 对抽水速率动态变化的影响，以及径向距离对断面流量的影响。hwQ-tQhwQQhwQ-tQhwhw图 6 展示了案例 B 中不同时刻 的 曲线的模拟结果。结果表明，固定时刻下 随着 的增加不断降低，在抽水早期不同曲线对应的 之间偏差较大，随着抽水的继续，流量偏差逐渐减小。的变化与 呈负相关关系，在抽水早期 曲线下降速度较快，随着抽水的继续 下降的趋势越来越平缓，符合 的小幅降低会使流速大幅度升高的事实。流量变化实际上是由抽水井内部水头下降而引起的水力梯度变化造成的，因此，较大的 意味着较小的比流量，从而导致抽水流量较小。r=2 mQ-tt=1 dh-rQn1n1图 7 展示了案例 C 中非达西
41、湍流系数对泵送流量的影响，其中图 7a 和图 7b 分别为固定径向距离 时的 曲线和相同抽水时间 时的 曲线。从图 7a 中可以发现，在抽水早期，固定时间点的 随着 的增加而增加。而随着抽水的继续，不同 对应的流量偏差不断减小。从理论上讲，这是可以理解的。抽水井附近含水层补给水量以重力释水为主，如式(10)所示，假设其等于比流量乘以过流断面面积，00.20.40.60.81.0t/d2.102.152.202.252.302.352.402.452.50h/mr=0.87 mr=3.62 mr=0.87 mr=3.62 m有限差分解实测数据图 5 有限差分解与实测数据对比Fig.5 Finit
42、e difference solution vs measured data 130 煤田地质与勘探第 51 卷n1Q-t比流量越大，抽水流量越大。由此可见，在抽水早期，抽水井附近区域含水层的重力储水逐渐释放并补充到抽水井内，越大，湍流程度越大，导致非承压含水层储水的释放越快。随着抽水的继续，N 区重力储水的释放逐渐趋于稳定，抽水井中的水主要通过 D 区的水流进行补给，导致 曲线逐渐趋同。t=1 drhn1图 7b 体现了 时不同径向距离的非达西湍流的发展特征。结果表明，固定 下 随着 的增大rhrhrh-rrn1而增大，在 较小时不同曲线对应的 之间的偏差较大，随着 的增大，偏差逐渐减小。非
43、达西系数与水头之间存在正相关关系，在 较小时 曲线上升的速度较快，随着 的增大，曲线斜率逐渐减小。越大，表示水力梯度越大，对应流量越大。在抽水井附近，水头波动较大，但没有出现不合理的突变，证明了所提有限差分解的适用性。t=0.1 dt=100.0 dQ1-rt=0.1 dQ1rQ1t=100.0 dRR-tQ1rQ1RQ1Q1r图 8 展示了径向距离对断面流量的影响，其中图 8a 和图 8b 分别用于表示抽水前期()和抽水后期()的 曲线(Q1为断面流量)，其水力参数见表 1 的案例 D。对比图 8a 和图 8b 可以发现，时 随 的增加而快速下降，其 较 时更大。随着抽水时间的增加，逐渐变小
44、，该结论符合图 4 中 曲线的趋势。随着 的增大不断减小，在 N 区部分 曲线随距离下降的速度较快，并在经过转换界面 后出现突变，在 D 区 曲线随距离下降的速率明显变慢。从理论上来说，定水头抽水过程中，在井壁处存在变化的最大值，随着 的增大，距离抽水井越远，对抽水井补给量越小，断面流量越小。在抽水前中期，含水层中距离抽水井较远的部分无法及时对抽水井附近进行充足的补给，抽 102101100101102t/d510152025303540Q/(m3d1)hw=9 mhw=8 mhw=7 m案例 B图 6 泵井内不同水头下流量的动态变化Fig.6 Dynamic development of t
45、he flow rates with differenthydraulic head inside the pumping well 020406080100120140Q/(m3d1)n1=1.4n1=1.2n1=1.0案例 Cn1=1.4n1=1.2n1=1.0案例 C05101520r/m9.59.010.010.511.011.512.0h/m102103101100101102t/d(a)r=2 m 时 Q-t 曲线(b)t=1 d 时 h-r 曲线图 7 不同非达西湍流系数对泵送流量的影响Fig.7 Effects of different non-Darcy coefficien
46、ts on pumpingflow 101100101r/m(a)t=0.1 d4045505560Q1/(m3d1)Q1/(m3d1)N 区D 区R101100101r/m(b)t=100.0 dR1012141618202224N 区D 区图 8 不同时间下径向距离对断面流量的影响Fig.8 Influence of radial distance on cross-section flowat different time第 10 期史鹏钰等：非承压含水层定水头抽水两区井流数值模型研究 131 Q1-rr=R水井附近区域含水层的重力储水会逐渐释放并补充到抽水井内。随着抽水时间的增长，在抽
47、水后期，远处的水流逐渐补充到抽水井附近，抽水井周边区域的水位和断面流量的下降速度逐渐变小。相比于 N 区，D 区水力梯度更小，其比流量随径向距离的下降速率更慢，其 曲线的转折点出现在 处。5结论a.所提出的有限差分解可用于在由定水头抽水试验引起的两区流、纯达西流和纯非达西流的情况下进行降深模拟。在两区流动的情况下，所提出的解也可用于评价抽水流量和非达西区域的动态发展。b.忽略非达西区域的影响可能会给两区流动的水头模型带来误差。有限非达西区域的湍流会分别导致两区流中水头较纯非达西流和纯达西流的水头偏大和偏小，且随抽水时间的增加逐渐变大。hwn1hwn1c.利用 Izbash 方程，发现瞬态抽水速
48、率时间变化关系可以由抽水井内水头 和非达西滞流系数 显著控制。在抽水过程中，减小 或增大 可以提高抽水速率，这种影响随着抽水时间的增加会不断减弱并在抽水结束时消失。Q1rr=RQ1RQ1d.在抽水过程中，随着 的增大不断减小，并在 处出现转折点，具体表现为：在 N 区中 曲线随距离下降的速度较快，在经过转换界面 处后出现突变，在 D 区中 曲线随距离下降的速度明显变慢。参考文献(References)MARKLE J M，ROWE R K，NOVAKOWSKI K S.A model forthe constanthead pumping test conducted in verticall
49、y fracturedmediaJ.International Journal for Numerical and AnalyticalMethods in Geomechanics，1995，19(7)：457473.1 CHEN C S，CHANG C C.Well hydraulics theory and data ana-lysis of the constanthead test in an unconfined aquifer with theskin effectJ.Water Resources Research，2003，39(5)：11211135.2 柯瀚，胡杰，吴小雯
50、，等.竖井抽水下垃圾填埋场渗滤液运移规律研究J.岩土工程学报，2018，40(5)：786793.KE Han，HU Jie，WU Xiaowen，et al.Investigation into leachatetransport in MSW landfills under pumping of vertical wellsJ.Chinese Journal of Geotechnical Engineering，2018，40(5)：786793.3 JACOB C E，LOHMAN S W.Nonsteady flow to a well of con-stant drawdown i
