承压含水层大口井不稳定非达西流简化解析解

肖良，柳岳辉，陈立华，于江伟

(1.广西大学 土木建筑工程学院，广西 南宁 530004；2.广西防灾减灾与工程安全重点实验室，广西 南宁 530004；3.中天建设集团有限公司广西分公司，广西 南宁 530004)

0 引言

承压含水层抽水井流模型是地下水渗流理论的重要内容，抽水试验作为直接获取研究区水文地质参数的有效途径，一直被学者和工程师广泛采用。经典泰斯模型给出了水位降深(s)与其自变量(时间t和井距r)的解析解，从而能粗略估计研究地水文地质参数。该模型主要适用于完整井附近的达西渗流，即抽水井渗流场的雷诺系数介于1至10之间。根据室内试验和工程实践发现，抽水井附近渗流容易出现非达西现象，即抽水井渗流场的雷诺系数小于1或者大于10，因此基于达西定律的泰斯抽水井模型在非达西渗流应用上具有一定的局限性。此外，在发展中国家，大口井(指抽水条件下井中的储水效应不能忽略的抽水井)已被广泛应用于地下水开采中。大口井的直径一般在2～8 m，具有进水断面大、构造简单、取材容易、使用年限长等特点[1]。如何合理开发利用地下水、正确评价地下水资源、防止地下水污染等问题，已成为学者们的研究热点[2-4]。因此在考虑井储效应的基础上，开展完整井附近非达西渗流模型研究对于正确认识抽水试验具有一定理论支撑作用。

当抽水速率足够大时，在孔隙较大的多孔介质含水层或含有较多裂隙的含水层中，地下水渗流速度与水力梯度之间不再满足线性关系，对于此类不满足达西线性定律的地下渗流称之为非达西流[5]。由于达西定律不能满足所有地下水渗流情况，地下水流非线性运动机理的研究自20世纪以来便成为水文地质领域的研究热点，广大学者对不同介质下的非达西流进行试验及理论研究，并提出相应的非达西渗流公式。其中Forchheimer于1901年提出的Forchheimer公式和Izbash在1931年提出的Izbash公式在非达西流建模中应用最为广泛，具体选用哪个公式一般取决现场水文地质条件[6]。

目前关于大口井水力学问题的研究大多基于渗流满足达西定律而展开[7-14]。对于考虑井储效应下非达西井流建模研究，土耳其学者SEN[15]运用Forchheimer公式描述非达西流，通过Boltzmann变换得到完整井附近非达西流近似解析公式，并利用特征曲线对公式中的相关参数进行特征分析。WEN等[6]采用Laplace变换以及线性化手段求得单一抽水井附近非达西流降深半解析解，通过建立数值解对该解析解进行对照验证，发现线性化方法用于描述非达西流而解得的近似解在抽水前期的计算结果与数值解存在一定误差。陈建华等[16]用渗流理论对大口井地层气体进行模拟，将井筒内视为“达西-非达西”两区模型，并通过雷诺数来判定其分界点，得到井底压力分布图。上述研究为大口井附件非达西渗流模拟提供了一定的理论基础，但是其方法计算过于复杂，依赖于计算机程序以及特征曲线拟合等，不便于水文地质工作者直接使用。

因此，本文拟在无限承压含水层的等速抽水试验中，提出考虑井储效应的简化非达西完整井流模型。笔者采用Izbash公式来描述非线性非达西流，采用Boltzmann变换和线性化近似相结合的方法推导得其简化解析模型，并通过绘制降深-时间曲线对解析模型进行参数分析，为利用大口井开采的相应地下工程实施提供理论依据及参考。

1 数学模型与方法

承压含水层大口井附近水动力模型与文献[13]中物理模型一致。如图1所示，假设含水层中各向同性且均质，含水层厚度不变且水平无限长，且抽水速率始终保持恒定。在基于泰斯的假设条件之外，考虑大口井自身井储的影响，建立起数学模型如下文。

图1 无限承压含水层大口径井渗流模型

流体连续性方程：

(1)

式中，r为距抽水井井中心的距离,m；t表示为抽水时间,s；q(r,t)为单位水流通量,m/s；S为储水系数，用来表征含水层储水能力的无量纲系数;b为含水层厚度,m;s(r,t)表示为水位降深,m。

抽水前天然状态下水力梯度为零，水平面保持一致且降深为零。模型初始条件为

s(r,0)=0。

(2)

在含水层无限远处的边界条件为

s(r→∞,t)=0。

(3)

考虑到抽水井自身井筒的储水，根据质量守恒定律，井内部的定解关系可描述为

(4)

式中，Q为抽水速率，m3/s;rw为抽水井滤网有效半径,m；rc为井套半径,m3；sw(t)为抽水井井内降深,m。

式(4)说明在大口井中，由于井径较大，其储水能力较强，抽水过程中井内流量来自于两方面：一是井筒内储水，即井储效应；二是含水层中对于井内的补给流量。

运用Izbash公式来描述非达西流，其水流通量与水力梯度关系为

(5)

式中，n、k为经验系数，k为准水力传导系数,(ms-1)n，是用来描述介质传输水流快慢的一个特征参数;n是用来描述非达西效应的无量纲系数，一般取值为1～2。当n等于1时，式(5)即为达西定律。

由式(5)求得q对于r的偏导数：

(6)

将式(6)代入式(1)可得

(7)

为除去式(7)中的非线性项，做线性化近似假设如下：

(8)

式(8)表示在抽水稳定时，单位时间内含水层中离井中心距离为r的任意过水断面的流量均为Q，即2πbrq=Q。文献[6]和文献[17]研究已表明该近似替换引起的模型计算误差非常小，可用于非达西流建模过程。

将式(8)代入式(7)可得非达西流线性近似为

(9)

将式(5)代入式(4)可得

(10)

Boltzmann变化在求解关于偏微分方程问题上是一种常用且便利的技术手段，因此采用Boltzmann变化将上述偏微分方程转换为常微分方程[18]，引入参量v(η)和η如下所示：

(11)

(12)

式中，v(η)为s和t的相似替代项；η为r和t的相似替代项;β和α为常数。

通过应用式(11)、式(12)，将式(7)中三个不同偏微分项可统一替换为关于η的方程，并如下所示：

(13)

(14)

(15)

假设β=2α以及γ=0，式(9)可表述为

(16)

由式(16)解得

(17)

将上述变化代入其初始条件及边界条件中，式(2)、式(3)、式(10)可改写为

v(η→∞)=0，

(18)

v(0)=0，

(19)

(20)

由式(12)、式(20)解得

(21)

从而求得

(22)

联立式(17)、式(22)可解得

(23)

将式(23)进行化简及积分求解可得

(24)

由式(18)、式(19)可知

(25)

因此联立式(24)、式(25)解得

(26)

(27)

将式(27)简化如下：

(28)

2 计算结果分析

2.1 模型准确性分析

本文将上述模型解与文献[1]研究成果进行比较验证提出模型的准确性。假设无井储效应，即取rc、rw为0，可得M=1。因此，不考虑井储效应前提下，由式(27)可将承压无限含水层等速非达西完整大口井流模型可简化为

(29)

由此可得，式(29)与文献[17]中式(11)表达形式完全一致，由此可证明本文推导模型的准确性。

2.2 井储效应影响分析

为研究井储效应对抽水降深的影响，本文假设在一承压无限含水层中，令抽水速率Q=0.108 m3/s，观测井r取值为10 m，井筒半径rc为1 m，滤网半径取rw为0.5 m，含水层厚度b设为50 m，储水系数S取值为0.003，将非达西系数n设为1.4，k取5×10-6(m/s)1.4，分别采用式(27)、式(29)计算降深-时间曲线，结果如图2所示。

图2 井储效应对于降深影响

由图2可知，井储效应的存在导致其降深-时间曲线负偏离无井储效应的降深-时间曲线，即因井储作用导致其在某一个抽水时间点的降深小于无井储作用下的降深。抽水过程中井储效应对于降深的影响一般发生在抽水前期，而随着抽水过程的进行，抽水井自身的储水将逐渐被抽干，降深变化主要受含水层补给控制，因此井储效应的影响随着抽水的进行逐渐消失。

2.3 准水力传导系数k对于降深的影响

基于上述解析模型，本文将对准水力传导系数k和非达西系数n对特征曲线的影响进行分析。假设抽水泵抽水流量Q取值为0.008 m3/s，观测井r取值为10 m，井筒半径rc为1 m，滤网半径取rw为0.3 m，含水层厚度b设为50 m，储水系数S取值为0.003，将非达西系数n设为1.4，编写相关MATLAB程序，就k在取值为5×10-5、5×10-6、5×10-7(m/s)1.4三种情况下分别进行分析讨论。

如图3所示为不同k值情况下的观测井水位降深-时间曲线。结果表明，随着抽水过程的进行，承压含水层中水位降深均随时间增加而变大。在抽水前期，水位降深变化与准水力传导系数k之间呈正相关，而在中后期，呈负相关。一般而言，k值越小含水层中介质输送水流速度越慢。在保持抽水泵抽水速率恒定时，由于观测井周围含水层对于井内的补给速率不同会导致井内降深不同。因而，在抽水中后期，k值较大时，观测井内水位降深相对于k值较小时降深相比较小。

图3 不同k值下水位降深-时间曲线

2.4 非达西系数n对于降深影响

同理，取假设抽水泵抽水流量Q取值为0.008 m3/s，观测井r取值为10 m，井筒半径rc为1 m，滤网半径取rw为0.3 m，含水层厚度b设为50 m，储水系数S取值为0.003，准水力传导系数k取值为5×10-6(m/s)n，n取值分别为1.2、1.4、1.6三种情况进行讨论。

不同n取值下降深-时间曲线模拟结果如图4所示。由此可知，非达西系数n与准水力传导系数k在对水位降深影响的效应相同；即在前期表现为非达西系数n与准水力传导系数k之间呈现正相关，而在抽水中后期呈负相关。在水力学中，水头指的是含水层中单位液体所具有的机械能，包括位置水头、压强水头、流速水头。随着抽水的进行，井内发生能量转换，对于观测井而言，主要表现在井内降深(即水头)的增大。由于流体的连续性，抽水井与含水层之间进行着能量的转换。因此，井内水位降深变化，侧面反映着流体能量的交换。一般来说，当n值越大(n>1)，流体的紊动程度越大，流体由于惯性力占据主导，流体内包含的势能较大。紊流相对于的层流(n=1的达西流)，短期供水能力更强。因此，在同一抽水时刻，非达西系数n越大的流体，将自身势能转换补给井内水头能力越强，从而水位降深变化反而较小。此外，也可解释另一规律：当非达西系数n的取值越大，非达西流趋于稳定的时间越早。

图4 不同n值下降深―时间曲线

2.5 参数敏感性分析

敏感性分析是一种量化输入量与输出量之间的不确定性关系的手段，参考文献[19]提出的标准化公式，本文将作为标准化的敏感度自变量定义为

(30)

(31)

式中，ΔPj一般为10-2×Pj。运用式(31)对Izbash公式中两个常数(n和k)进行参数敏感性分析(图5、图6)，选取3.2小节中同一组参数分别做出k=5×10-5(m/s)1.4、n=1.1的NSC-时间曲线图。

图5 n值敏感性分析

图6 k值敏感性分析

由图5、图6可知，Izbash公式中的两个经验系数中非达西系数n相对于准水力传导系数k而言的敏感性要强得多。从Izbash公式的结构来看，n作为反映水流非达西效应强弱的主要系数，相对于k而言，敏感性更强。

4 结论

本文提出了一种简化数学解析模型用于求解非达西效应下承压无限含水层大口井附近的水位降深，通过对所得到的的简化解析解进行分析，得到如下结论：

① 井储效应仅出现在抽水前期，随着井筒内部储水的消耗，水位降深的变化仅来自于含水层的补给。

② 本文通过对Izbash公式的两个参数进行分析，抽水前期水位降深与准水力传导系数k和非达西系数n为正相关，而在抽水后期呈负相关。且当n的取值越大，非达西流趋于稳定的时间越早。

③ 通过对Izbash公式中的两个参数敏感性分析，相较于准水力传导系数k，非达西系数n敏感性更强。研究结果将对水文地质参数的获取、地下水资源评估等工作提供了一定的理论支撑。