河渠水位线性变化条件下河渠 - 潜水 非稳定流模型及其解 *

吴 丹, 陶月赞, 林 飞

(合肥工业大学 土木与水利工程学院, 合肥 230009)

摘要 在河渠水位迅速变化后再缓慢变化的条件下,建立了河渠半无限潜水含水层中非稳定渗流模型 利用Boussinesq第一线性化方法及Laplace变换,并注意应用Laplace变换中的“积分性质”,给出形式相对简单、由常用函数表达的解,阐述特定解及其相应的物理意义 由解所揭示的潜水位变化规律表明,含水层任一点处潜水位变动速度的时间变化曲线形态是固定的,与河渠边界水位变动速率 λ 无关;潜水最大变速发生的时间,随 λ 呈非线性位移 依据潜水位变化规律,建立利用潜水位变动速度求含水层参数的方法,并用实例演示了拐点法求参数的过程

潜水非稳定流; 河渠边界; 水位线性变化; Laplace变换; 积分性质

引 言

河渠附近潜水非稳定渗流问题,是地下水渗流力学中最为经典的问题之一 [1-3] 其中,被直线型河渠完全切割的半无限水平潜水含水层中的潜水一维渗流模型(河渠-半无限潜水渗流模型),是最基本同时也是最重要的河渠-潜水渗流模型 [1-4] 大量的研究文献以河渠-潜水渗流为基础,针对不同的水文地质条件,对模型中的初始条件、边界条件、源汇项等进行相应地修改 [4-10] ,以获得相应问题的解,从而为灌排水渠系统设计 [1-4] 、地表水体附近或河间地块中的潜水非稳定变化过程与渗流规律研究 [5-10] 、地表水与地下水之间相互作用的定量评价 [7-9] 、河道水文过程研究 [11-12] 、潜水含水层参数确定 [13-15] 等,提供了基本的理论工具

根据边界(视为一类边界的河渠完全切割潜水含水层)条件和有无源汇项,河渠-半无限潜水渗流模型,多假定河渠水位变化瞬时完成之后,在很长的计算期内,水位保持稳定 [1,3-9,13] 实际上,水位变化多是快速变化后再连接一相对缓慢的变化过程,如节制闸关闸蓄水,节制闸附近的渠水位快速上升后,因水流惯性而导致渠道有一较长时间的水位持续缓慢上升过程 为解决这类水位变动的一类边界问题,现有研究主要有2种方法:1) 对水位变化过程进行离散,使得 “水位分段水平”,再采取水流叠加原理,进行求解 这种方法相对简单,但实际应用中,当水位变化过程较复杂时,需要将水位过程分许多段,才能基本满足“水位分段水平”;分段过多,计算繁琐 2) 对水位变化过程进行相应的离散处理,不要求“水位分段水平”,由此求出的解,往往是不常用的特殊函数 [1,4] 此时,需要对特殊函数建立相应的计算方法,这为方法推广带来不便

对河流水位迅速变动后,水位呈缓慢的线性变化时的半无限潜水渗流模型,在利用Laplace变换时,应用Laplace变换的“积分性质”,可获得形式相对比较简单的解

1 渗 流 模 型

一顺直河渠,如图1所示,其所处地段的水文地质条件,可概括为:

图1 河渠附近潜水渗流场
Fig. 1 The phreatic seepage field by a stream

① 潜水含水层均质各向同性、具水平的隔水底板,在平面上无限延伸;

② 一条在剖面上基本完整切割含水层的河流,河流水位初期迅速升高Δ H 0 后,水位呈缓慢的线性变化;

③ 潜水初始水位 h ( x ,0)水平;

④ 潜水水流可视为一维流

对应经典的Ferris模型,上述水文地质概念模型中:条件②,由“水位迅速升高Δ H 0 后、水位保持不变”,修改为“水位迅速升高Δ H 0 后、水位呈缓慢的线性变化”

该问题的数学模型可写成

式中, μ 是含水层给水度; k 是含水层渗透系数,m/d; h 是地下水水位,m;Δ H 0 是河水迅速变化初期段的水位变化幅度,m; λ 是水位在计算期间随时间缓慢变化的速率,m/d; x 是计算点距边界的距离,m

2 模 型 的 解

对于模型 ,当 h ( x , t )- h ( x ,0)≤0.1 h m ( h m 为潜水流的平均厚度,这在实际中大都能满足)时,可利用Boussinesq方程第一线性化方法,令 u ( x , t )= h ( x , t )- h ( x ,0);再对模型 求关于 t 的Laplace变换,可得模型

(1)

(2)

(3)

式中, u 关于 t 的Laplace变换的象函数, s 为Laplace算符, a = kh m / μ a 为潜水含水层的导压系数,m 2 /d

模型 中,式(1)的通解为

由上述通解,结合边界条件(2)、(3),可得模型 的定解:

(4)

对式(4)进行Laplace逆变换,对1/ s 2 项求逆变换时,注意Laplace变换中的“积分性质”,即 是关于 f ( t )的Laplace算符 注意 u ( x , t )= h ( x , t )- h ( x ,0),可获得上述模型的解:

(5)

式中, f erc ( z )是余误差函数,

式(5)是在河渠边界控制下的半无限潜水含水层中,在水位迅速升高Δ H 0 后、水位呈缓慢线性变化的河渠水位影响下,潜水非稳定渗流过程的解析解 该式是由比较常用的函数表达的;与现有类似模型的解相比,形式也比较简单

3 解的数理特征

3 . 1 特定解及其水文地质意义

3.1.1 λ =0

λ =0,对应模型中的条件②,修改为“水位迅速升高Δ H 0 后、水位保持不变”;此时,式(5)转化为式(6):

h *

(6)

式(6)是经典的Ferris模型解;为便于以下讨论,式(6)计算的 h ( x , t ),记为 h * ( x , t )

3.1.2 x →0

此时,对应的是边界点上的地下水水位,与河渠水位变化同步、同幅度,这与地下水渗流力学的基本原理相符

3.1.3 x →∞

此条件对应的是距离边界无穷远处的地下水水位,与河渠水位变化无关,这与地下水渗流力学的基本原理相符;实际情况中,一般距离边界有一定远距离时,河渠水位变动对地下水水位的影响就基本可以忽略了

3 . 2 水位复杂变化条件下的解

当河渠水位 H ( t )变化过程复杂,此时,可根据 H ( t )实测过程,采用线性插值方法,先将评价期划分为若干个计算时段,在每个时段内, H ( t )被描述为线性变化 则有

(7)

γ( t - t i -1 )为Heaviside函数, 具有如下性质:当 t < t i -1 ,γ( t - t i -1 )=0;当 t t i -1 ,γ( t - t i -1 )=1 则有

(8)

依据式(5),将边界条件(2)变换为式(8),根据叠加原理,可得问题的解:

(9)

式(9)是在河渠水位采用线性插值离散条件下,潜水非稳定渗流过程的解析解

3 . 3 潜水位变动规律与应用

潜水位变动速度 φ ( x , t )=∂ h ( x , t )/∂ t ,由式(5)可得

(10)

(11)

3.3.1 λ 形成的潜水位变速的时间变化规律

由式(10),潜水位变速是Δ H 0 λ 两种作用的合成

当Δ H 0 =0时,潜水位变速就完全由 λ 作用形成;由式(10), λ 作用形成的潜水位变速为式(10)右端第2项,是河渠水位后期变化速率 λ f erc ( z )的乘积

对于一特定含水层( a 值确定)中的某一点( x 值确定),其潜水位变动速度 φ ( x , t )随时间 t 的曲线形态是固定的,与河渠边界水位变动速率 λ 无关

图2为 x =60 m时,在不同 a 值条件下的 φ ( x , t )- t 曲线

φ ( x , t )- t 的影响规律, x 等效

图2 不同a值时λ在x处形成的潜水位变速
Fig. 2 The variable speed of the phreatic level formed by λ at point x under different values of a

3.3.2 λ 对潜水位最大变速的影响规律

φ ( x , t )- t 上的拐点,为潜水获得最大变速的时刻,是研究潜水渗流问题的一个重要变量

当Δ H 0 ≠0∩ λ =0时,解转化为Ferris模型解 此时,将 φ ( x , t )记为 φ 0 ( x , t ); φ 0 ( x , t )- t 存在一拐点、拐点处时间为 t k 由式(11),有

t k = x 2 /(6 a )

(12)

当Δ H 0 ≠0∩ λ ≠0时, φ ( x , t )- t 也存在一个拐点,拐点处的时间为 t g ;由式(11),有

(13)

式(13)表明,在 λ 作用的影响下, φ ( x , t )- t 的拐点时间发生时间位移

由式(10), φ ( x , t )= φ 0 ( x , t )+ λf erc ( z ), f erc ( z )为单调递增函数, λf erc ( z )将导致拐点时间位移有如下规律(如图3):相对于 λ =0时的拐点时间 t k λ >0时拐点发生时间 t g 向后延时,即 t g > t k λ 越大延时越长; λ <0时延时规律相反;不同 λ 值对应的拐点发生时间不在一条直线上,也即,潜水最大变速发生的时间随 λ 呈非线性位移

图3 潜水变速拐点随λ位移
Fig. 3 The inflection point of the phreatic variable speed vs. displacement λ

3.3.3 λ 对潜水位计算的影响

对于水位持续变化的河渠,如未考虑河渠水位变动的影响,而直接用Δ H 0 来计算潜水水位,形成的计算误差为Δ h ,是式(10)右端第2项;Δ h 随时间的变化,记为 V Δ h

(14)

V Δ h 为相对式(6)计算的地下水水位变幅,有

V Δ h /[ h * ( x , t )- h ( x ,0)]=

λ H 0

(15)

式(14)表明:不考虑河渠水位变化导致的潜水变幅计算误差,其误差的变动速度 V Δ h ,相比Δ H 0 独立形成的潜水水位变幅,为 λ H 0 是一常数 这表明,Δ H 0 引起的地下水变幅、 λt H 0 之后的水位变化) 引起的潜水变幅,两者对潜水水位的作用相对独立,潜水水位由这两者合成而成;这与式(5)所表达的数学规律相一致,也符合地下水流叠加原理

3.3.4 应用

由上述规律,可直接应用于实际工作中

适线法 当Δ H 0 =0∩ λ ≠0时,根据 λ 形成的 φ ( x , t )- t 规律,可建立求解含水层参数 a 的适线法 x 监测点,构建与不同 a 值相对应的 f erc ( z )- t 理论曲线族;由 x 点实测的潜水水位动态,获得 φ ( x , t )- t 曲线;当形成 φ ( x , t )- t 的含水层参数 a 值,与 f erc ( z )- t 理论曲线图族中某曲线 a 值相等时,则 φ ( x , t )- t f erc ( z )- t 的曲线形态完全相同,仅相差一个常数 λ 倍,也即两条曲线完全重合

因此,根据潜水水位动态测验数据,形成 φ ( x , t )- t 曲线,将之与 f erc ( z )- t 理论曲线族进行适线,可确定出含水层的 a 值(如图2)

拐点法 当Δ H 0 ≠0∩ λ ≠0时,根据 λ 形成的 φ ( x , t )- t 曲线拐点,可建立求解含水层参数 a 的拐点法 x 监测点,由 x 点实测的潜水水位动态,绘出 φ ( x , t )- t 曲线;利用曲线的拐点 t g ,以及河渠水位实际变动数据而获得 H 0 λ ;由式(13),可计算出模型参数 a (此时, H 0 λ x 都是已知值)

4 实 例 研 究

安徽淮北平原中部,以粉细砂为主的潜水含水层发育广泛,厚度8 m左右,底部一般发育有不完全连续黏性土层;由于潜水位埋深浅,为2.5~3.0 m;区内农田灌溉渠系统比较完善,干渠基本深切至隔水底板、干渠的渠间距为2 km左右,干渠渠首多有节制闸控制;一口国家级地下水位自记观测井,距离干渠直线距离为60 m处,观测井附近地面标高31.02 m

图4 拐点法求a
Fig. 4 Calculation of a by means of the inflection point method

2013年8月26日,干渠关闸蓄水灌溉 26日15时关闸至27日15时,渠道水位变化过程可划分为2个明显的时段:关闸15 min内,水位迅速上升,升幅为2.00 m(由26.80 m上升至28.80 m);之后,水位并不是保持不变,而是一直在缓慢上升,至27日15时升幅为0.21 m(上升至29.01 m) 例中,地下水观测孔水位按3 h摘录,如表1

该时段, H 0 =2.0 m、 λ =0.21 m/d,可采用拐点法求参数 如图4, t g =16.5 h=0.68 d;拐点法结果 a 值为860 m 2 /d

文献[13]中,以同一地段地下水动态数据获得的 a 值为855 m 2 /d,两结果基本一致

表1 潜水水位动态数据(2013-08-26~2013-08-27)与“拐点法”计算过程

Table 1 The dynamic data of water table (2013-08-26~2013-08-27) and the process of “inflection point method”

t/h36910111215182124h/m26.8426.8926.9527.0127.0827.1527.2327.3127.3927.46φ(x,t)/(m/h)1.33×10-21.67×10-22.00×10-22.00×10-22.33×10-22.33×10-22.67×10-22.67×10-22.57×10-22.43×10-2

初始水位为26.80 m,附近地面标高30.66 m

Note The initial water table is 26.80 m and the ground level is 30.66 m

5 讨论与结论

在上述研究过程中,形成以下结论:

1) 实际工作中,水闸关闸引起的河渠水位变化,通常是一迅速上升再连接一缓慢上升过程;本文给出的河渠-半无限潜水渗流模型,更为符合上述河渠水位变化过程

2) 借助Laplace变换,可获得形式相对比较简单、由比较常用的函数表达的解

3) 河渠初期快速上升的水位变幅Δ H 0 λ 形成的水位变化,两者对地下水水位的作用相对独立,地下水水位变化由这两者作用合成

4) 对于一特定含水层( a 值确定)中的某一点( x 值确定), 其潜水位变动速度 φ ( x , t )随时间 t 的曲线形态是固定的, 与河渠边界水位变动速率 λ 无关 φ ( x , t )- t 的影响规律, x 等效

5) 相对于Ferris模型解,潜水最大变速发生的时间随 λ 呈非线性位移, λ >0时拐点发生时间向后延时, λ 越大延时越长; λ <0时延时规律相反

参考文献 ( References ):

[1] 张蔚榛. 地下水非稳定流计算和地下水资源评价[M]. 北京: 科学出版社,1983.(ZHANG Weizhen. Calculation of Unsteady Flow of Groundwater and Evaluation of Groundwater Resources [M]. Beijing: Science Press, 1983.(in Chinese))

[2] BEAR J. 多孔介质流体力学[M]. 李竞生, 陈崇希, 译. 北京: 中国建筑工业出版社, 1983.(BEAR J. Dynamics of Fluids in Porous Media [M]. LI Jingsheng, CHEN Chongxi, transl. Beijing: China Architecture & Building Press, 1983.(Chinese version))

[3] 束龙仓, 陶月赞. 地下水水文学[M]. 北京: 中国水利水电出版社, 2009.(SHU Longcang, TAO Yuezan. Groundwater Hydrology [M]. Beijing: China Water & Power Press, 2009.(in Chinese))

[4] 陶月赞, 姚梅, 张炳峰. 时变垂向入渗影响下河渠-潜水非稳定渗流模型的解及应用[J]. 应用数学和力学, 2007, 28 (9): 1047-1053.(TAO Yuezan, YAO Mei, ZHANG Bingfeng. Solution and its application of transient stream/groundwater model subjected to time-dependent vertical seepage[J]. Applied Mathematics and Mechanics , 2007, 28 (9): 1047-1053.(in Chinese))

[5] MAHDAVI A. Transient-state analytical solution for groundwater recharge in anisotropic sloping aquifer[J]. Water Resources Management , 2015, 29 (10): 3735-3748.

[6] Bansal R K. Approximation of surface-groundwater interaction mediated by vertical stream bank in sloping terrains[J]. Journal of Ocean Engineering and Science , 2017, 2 (1): 18-27.

[7] HUANG F K, CHUANG M H, WANG G S, et al. Tide-induced groundwater level fluctuation in a U-shaped coastal aquifer[J]. Journal of Hydrology , 2015, 530 (1): 291-305.

[8] MUNUSAMY S B, DHAR A. Homotopy perturbation method-based analytical solution for tide-induced groundwater fluctuations[J]. Groundwater , 2016, 54 (3): 440-447.

[9] SU N H. The fractional Boussinesq equation of groundwater flow and its applications[J]. Journal of Hydrology , 2017, 547 (2): 403-412.

[10] DAVID P V, SAHUQUILLO A, ANDREU J, et al. A general methodology to simulate groundwater flow of unconfined aquifers with a reduced computational cost[J]. Journal of Hydrology , 2007, 38 (1/2): 42-56.

[11] KIM K Y, KIM T, KIM Y, et al. A semi-analytical solution for groundwater responses to stream-stage variations and tidal fluctuations in a coastal aquifer[J]. Hydrological Processes , 2007, 21 (5): 665-674.

[12] SERRANO S E, WORKMAN S R, SRIVASTAVA K, et al. Models of nonlinear stream aquifer transients[J]. Journal of Hydrology , 2007, 336 (1): 199-205.

[13] 陶月赞, 曹彭强, 席道瑛. 垂向入渗与河渠边界影响下潜水非稳定流参数的求解[J]. 水利学报, 2006, 37 (8): 913-917.(TAO Yuezan, CAO Pengqiang, XI Daoying. Parameter estimation for semi-infinite phreatic aquifer subjected to vertical seepage and bounded by channel[J]. Journal of Hydraulic Engineering , 2006, 37 (8): 913-917.(in Chinese))

[14] MIZUMURA K. Approximate solution of nonlinear Boussinesq equation[J]. Journal of Hydrologic Engineering , 2009, 14 (10): 1156-1164.

[15] YOUNGS E G, RUSHTON K R. Dupuit-forchheimer analyses of steady-state water-table heights due to accretion in drained lands overlying undulating sloping impermeable beds[J]. Journal of Irrigation and Drainage Engineering , 2009, 135 (4): 467-473.

Solution of the Transient Stream-Groundwater Model With Linearly Varying Stream Water Levels

WU Dan, TAO Yuezan, LIN Fei

( School of Civil and Hydraulic Engineering , Hefei University of Technology , Hefei 230009, P . R . China )

Abstract : Based on the first linearized Boussinesq equation, the analytical solution of the transient groundwater model for description of phreatic flow in a semi-infinite aquifer bordered by a linear stream with linearly varying stream water levels, was derived through the Laplace transform and in view of the integral property of the Laplace transform. The solution is composed of some common functions and its expression form is relatively simple. According to the mathematical characteristics of the solution, its corresponding physical meaning was discussed. The variation rule of the phreatic level revealed by the solution shows that the temporal variation curve of the aquifer at any point is fixed and has nothing to do with the change rate of the water level of the river channel. The time of the maximum speed change of the phreatic aquifer nonlinearly varies with λ . Based on the variation rule of the phreatic level, the method determining the aquifer parameters with the changing velocity of the phreatic level was established, and the process of obtaining the parameter with the inflection point method was demonstrated through an example.

Key words: phreatic flow; channel boundary; linearly varying water level; Laplace transform; integral property

Foundation item: The National Natural Science Foundation of China(51309071)

ⓒ 应用数学和力学编委会,ISSN 1000-0887

文章编号 1000-0887(2018)09-1043-08

吴丹(1986—),女,博士生(E-mail: wudanyt@163.com);

陶月赞(1964—),男,教授,博士生导师(通讯作者. E-mail: taoyuezan@163.com).

基金项目 国家自然科学基金(51309071)

作者简介

* 收稿日期 2017-09-06;

修订日期: 2017-11-21

文献标志码: A

DOI: 10.21656/1000-0887.380250

中图分类号 P641.132

①引用本文 / Cite this paper: 吴丹, 陶月赞, 林飞. 河渠水位线性变化条件下河渠-潜水非稳定流模型及其解[J]. 应用数学和力学, 2018, 39 (7): 1043-1050.WU Dan, TAO Yuezan, LIN Fei. Solution of the transient stream-groundwater model with linearly varying stream water levels[J]. Applied Mathematics and Mechanics , 2018, 39 (7): 1043-1050.