基于分段函数的潜水井降深曲线在工程实践中的应用

2022-06-10朱悦璐熊俊飞

人民黄河 2022年6期

朱悦璐，熊俊飞，邓炜

（南昌工程学院水利与生态工程学院，江西南昌 330099）

1 引言

地下水向井运动在矿山排水、灌区农业灌溉、地下水资源评价中经常会遇到［1-3］。自Dupuit［4］推导出潜水向井稳定流公式以来，计算降深曲线、预测地下水位变动等，一直是该领域的热点［5-7］。中国地质大学陈崇希教授团队与南京大学薛禹群院士团队在2010—2014 年就稳定流Dupuit 公式渗流场影响半径问题进行多次商榷［8-10］。常云霞［11］利用Visual Modflow 软件，模拟分析了矿井周边区域地下水降深曲线变化趋势。兰盈盈［12］采用GMS（Groundwater Model System）模拟软件，对2015—2028 年大坳灌区地下水位进行动态预测，结果表明新建灌区对地下水位影响较小，不会造成土壤盐渍化问题。闫峭［13］利用FlowHeat 1.0，对西安市某小区热力井在各种布井方案和抽灌条件下的渗流场、温度场进行了模拟，并优化了监测井网的布局方案。 Pulat 等［14］、Michopoulos 等［15］模拟了抽水井热源泵系统在不同降深下的水热平衡问题。

地下水向井运动的研究成果已较为丰富，但仍存在局限性，如潜水井空间径向流场的复杂分布，很难求得解析解。为解决该问题，研究者大多假设向井流动的潜水流近似水平，忽略垂向分速度，这一近似处理手段导致地下水计算软件的底层框架大多演变为传统的Dupuit 假设，将潜水运动方程简化为一维或二维Boussinesq 方程。在该假设下经验公式存在的问题：①靠近井壁处，降深曲线计算误差较大，不能用于实际工程；②当r ＞1.5H0（r为与抽水井的距离，H0为潜水初始水头）时，理论公式才与实际相符。对于①，有学者提出井壁区域渗流场计算的相关问题，例如应用有限元三角形相似原理［16］，在稳定渗流场中确定一个能量损失率最小点，取井壁处的固定水头点为溢出点，求解r ＜1.5H0这一区域的降深曲线，该方案迭代速度快，但仍有计算精度不足等缺陷；对于②，当r ＞1.5H0时，实际潜水降深与理论计算值也有很大差距，Dupuit公式仍不可直接应用，但目前对此区间计算精度进行修正的研究较少。

基于上述问题，笔者利用16 口观测井实测数据，在稳定抽水条件下拟合方程，生成计算曲线，并采用试算法求得Dupuit 理论公式在实际工程中应用的精度分界点，以该点为区间临界点，将计算曲线与理论曲线组成分段函数，最终得出经过修正后的潜水渗流场。该方案可有效减小传统降深曲线在全区间的计算误差，为同类场地条件下确定计算分界点、优化观测井布置等提供参考。

2 传统Dupuit 公式及其不足

潜水井向井运动是潜水含水层中地下水运动的一个特例，因此其动力学特性满足潜水含水层中地下水运动的基本微分方程及假设条件，即Dupuit 假设：对于垂直的二维平面上任意一点（x，z），当潜水流垂向坡角θ很小（潜水位变化很小）时，可以忽略垂向分速度vz而仅考虑水平分速度vx（见图1）。

图1 Dupuit 假设

潜水井地下水向井运动的Dupuit 公式为

式中：H为潜水水头；h为潜水含水层厚度；x、z分别为水平向、垂向坐标。

式（1）为直角坐标系下的标准微分方程，为方便与抽水井实际工况匹配，常用的方法是将式（1）转化为柱坐标：

大多数研究者认为Dupuit 假设没有考虑水流在入井时的能量损失，因此在抽水井壁附近的理论降深和实测降深差异较大［17-19］。当边界条件为：r ＝rw、h ＝hw和r ＝Rw、h ＝H0（rw为抽水井半径，Rw为抽水井影响半径，hw为抽水井内水头）时，按文献［17］对式（2）进行积分运算，得到潜水井Dupuit 理论曲线表达式：

潜水含水层厚度h可表示为观测井到抽水井距离r的函数h ＝f（r），当潜水初始水头H0已知时，由Dupuit 公式得出的潜水井附近的水位降深计算公式为

式中：s为水位降深。

很多实际工程算例表明，直接使用式（4）会导致计算误差偏大。笔者在充分考虑前人研究成果的基础上，以抽水井实测数据为基础，拟合计算降深曲线，并通过相对误差分析确定计算精度分界点，在分界点以内采用拟合曲线，在分界点以外采用理论曲线，二者组成分段函数，建立全区间符合实际工况的潜水渗流场。

3 实例应用

3.1 工程概况

某灌区试验场抽水井与周边16 口观测井W1、W2、…、W16的距离分别为43、60、90、120、135、140、170、220、270、480、510、590、650、700、730、780 m。布井场地范围内地质条件差异不大，可按各向同性考虑。抽水井是一个浅层开放性地下水系统，可接收大气降水、地表水、灌溉回归水等垂直入渗补给，通过浅层潜水蒸发、人工开采侧向径流等输出，地下水水力特性属于潜水-微承压水。

3.2 实测降深曲线数学模型

由于各井实际工况不同，因此判断抽水井在长时间定流量抽水后是否在井附近形成稳定降落漏斗，目前尚无统一标准。一般方法是选取一个抽水时段Δt，如果在该时段观测井水位基本无变化，则认为达到稳定状态。本研究选取Δt ＝200 min，监测抽水10、210、1 600、1 800 min 时各观测井的水位降深（见表1）。

表1 不同抽水时间观测井水位降深实测值

将上述4 个抽水时间的水位降深曲线绘制在同一张图上进行对比（见图2），可以看出：在抽水开始t＝10 min 和t＝210 min 时同一观测井的水位降深差异较大，在距抽水井较近的范围（ 0 ≤r≤300 m）差异较为显著；当t＝1 600 min和t＝1 800 min 时，同一观测井水位降深差异很小，表明此时抽水范围内已形成相对稳定的渗流场，降落漏斗形状和水位线保持稳定。

由图2 可知，t＝1 800 min 时，各观测井实测降深为稳定流数据。对该组数据应用Matlab 软件进行拟合，根据曲线形态和工程经验，选取幂函数、指数函数、有理数逼近以及傅里叶函数等4 种函数进行拟合，结果见图3。

图2 不同抽水时间各观测井水位降深曲线

图3 不同函数拟合效果对比

以幂函数拟合的降深曲线表达式为

以指数函数拟合的降深曲线表达式为

以有理数逼近拟合的降深曲线表达式为

以傅里叶函数拟合的降深曲线表达式为

其中：a0＝4.39×107，a1＝-4.39×107，b1＝1.593×104，w ＝-5.602×10-7。

采用误差平方和（SSE）、决定系数（R2）、校正决定系数（Adjusted R2）、均方根误差（RMSE）对上述4 种函数拟合程度的优劣进行判断，结果见表2。

表2 拟合程度指标

SSE越接近0，说明模型选择越合理、拟合越好、数据预测越成功；R2越接近1，说明拟合效果越好；Adjusted R2判断标准与决定系数相同；RMSE越接近0拟合效果越好。根据上述判断标准，由表2 可知，通过有理数逼近拟合的水位降深曲线在4 种函数中效果最好，基本可以反映稳定抽水时真实的水位降深状态，其函数表达式为

3.3 潜水渗流场及其母线方程

以抽水井底部中心为坐标原点，井轴线为z轴，建立直角坐标系。以式（9）为母线绕井轴线旋转，此时在潜水含水层中形成一个旋转曲面，该旋转曲面即为潜水含水层中的实际降落漏斗，其表达式为

式（10）即为研究工况下地下水潜水位分布的空间方程，该式所表达的渗流场基于工程实测资料生成，因此在靠近井壁区域，即当r≤1.5H0时，同样符合实际情况。

传统的此类研究，由于井壁附近边界条件复杂，其概念化的微分方程甚至不能写出显性表达式，更难于求出解析解，因此工程意义有限；而本研究采用实测资料拟合方案所推导出的方程，过程清晰、物理概念明确，降落漏斗空间形态是由真实降深曲线为母线旋转而来，且具有数值解，相比微分方程（例如二维潜水运动的Boussinesq 方程），更方便实际应用。

考虑到抽水井实际的坐标系模型，将式（10）通过坐标变换：x ＝rcosφ，y ＝rsinφ，z ＝z，以柱坐标形式表示，其方程可转化为式（9）。

4 结果与讨论

4.1 计算曲线与理论曲线对比

将本研究方案计算生成的降深曲线与按式（4）理论降深曲线在同一坐标系下对比，见图4（图4（b）中理论降深曲线仅为一般经验形态，无具体坐标数值）。可以看出，计算曲线和理论曲线在趋势上一致，整体上呈现离抽水井越远降深越小的趋势，当r大于一定值时，降深几乎不再变化；但在距抽水井较近位置时，二者形态有较大差异。本研究方案降深曲线与实测数据吻合较好，而理论降深曲线在井附近变动趋势更为剧烈，存在明显的“陡降”。一般认为，造成该现象的主要原因有二：一是由流速造成地下水入井水头损失，二是地下水触井后由近似水平运动转为垂直运动造成能量损失。

图4 计算生成的降深曲线与理论降深曲线对比

4.2 理论曲线分界点计算结果

利用式（4）计算得到的降深曲线见图5。可以看出，理论计算结果与实测结果偏差较大。传统理论认为“当r≥1.5H0时，理论降深曲线与实际降深曲线一致”，但该结论与本研究实际情况不符。研究区初始潜水含水层厚度为23.25 m，按照传统理论，当与观测井距离r≥1.5H0（为方便计算，不妨取40 m）时，理论曲线与实际曲线应当一致，而事实上直至r＝300 m 时理论值与实测值才较为接近（抽水井与最近的观测井距离设计为43 m）。

图5 理论曲线与实测值对比

为定量分析不同区间理论降深和实测值的差异，按照40＜r ＜120 m、120 m＜r ＜400 m、400 m＜r ＜800 m 将观测井分为3 个区间进行试算，通过考察不同区间理论降深和实测值的相对误差，确定符合实际工况的计算分界点。相对误差计算公式为

Er＝（ss- sm）／sm（11）

式中：Er为相对误差；ss为理论降深；sm为实测降深。

40 m＜r＜120 m 区间理论降深和实测值的相对误差见表3。可以看出，该区间4 组数据平均相对误差为63.75%，最大相对误差为65%；应用同样方法可得，120 m＜r ＜400 m 区间4 组数据平均相对误差为17%，最大相对误差为21%。上述说明理论曲线在本研究40 m＜r ＜120 m 和120 m＜r ＜400 m 这两个区间都不适用，尽管该距离已完全满足r≥1.5H0。当观测井处于400 m＜r ＜800 m 这一区间时，平均相对误差为7%、最大相对误差为9%，这表明直到观测半径r＞400 m 时，理论曲线精度才可满足工程需求，因此Dupuit 公式在实际应用中需十分慎重。

表3 40 m＜r＜120 m 区间理论降深和实测值的相对误差

4.3 计算曲线-理论曲线组合方程

在实际工程中一般无法大量设置观测井，因此判断精度分界点的位置显得尤为重要。由上述理论曲线计算结果可知，r＝400 m 可作为本研究计算误差分界点，在该范围之外潜水降深理论值与实测值一致，可不设置观测井；在该范围之内理论值与实测值误差较大，需设置观测井，不可用理论公式直接计算。显然对于具体工程，采用本文的计算分界点方案，比直接以r ＝1.5H0作为分界点更为精确且更具有实际意义。

将式（4）与式（9）联立，在实际分界点之内采用本研究构建的有理数逼近方程，在实际分界点之外采用Dupuit 方程推导出的理论降深公式，这样即可建立适用于本研究对象的自井附近至无限远处全区间范围降深分段函数：

4.4 分段函数工程验证

将观测井与抽水井的距离及其他相关数据代入式（12）进行计算，得出分段降深曲线（见图6）。以该曲线为母线绕井轴旋转所形成的降落漏斗即在本研究条件下的渗流场。为验证式（12）在实际工程中的适用性，停止抽水一段时间，待试验场地潜水分布回归稳定状态后重新进行独立抽水试验。重新抽水后，取t＝2 000 min，对比各观测井实测降深与计算降深。

图6 分段降深曲线

以理论分界点r＝1.5H0（取40 m）和计算分界点r＝400 m 将试验场地分为r＜40 m、40 m＜r＜400 m、r＞400 m 3 个区域，将各区间降深实测值、分段函数计算值、Dupuit 理论曲线值三者对比，以判断计算曲线和理论曲线在实际应用中的差异。

在r＜40 m 区间，选择井壁附近5、20、30 m 处3 个观测孔，分别对比实测值、计算值、理论值，见表4。可以看出，本文所构建的分段函数计算值皆优于Dupuit理论曲线计算值，尤其在靠近井壁区域。当观测井距抽水井为5 m 时，Dupuit 理论曲线值与实测值相差2.67 m，相对误差为56%，理论值与实测值偏差较大，Dupuit 公式显然无法应用于实际工程；相同位置处，本文构建的函数曲线计算结果与实测值绝对误差仅为0.16 m，相对误差仅为3%，其计算精度远远高于理论曲线的。

表4 井壁附近观测降深与计算降深对比

在40 m＜r＜400 m 区间，选取6 个观测点，进行理论值、计算值与实测值对比，见图7。可以看出，实测降深与理论值在该区间仍有较大偏差，绝对误差最大值出现在r＝135 m 处，其值为0.31 m，最大相对误差为15%；实测值与计算值趋势则较为吻合，绝对误差最大值出现在r＝265 m 处，其值为0.15 m，最大相对误差为8%。

图7 较远处理论值、计算值与实测值对比

进一步，选用建模序列对比中常用的总量偏差BLAS、均方根纳什效率系数SNSE两个校验指标进行判断，其中BLAS越小则效果越好，SNSE越接近1 效果越好，其数学意义见参考文献［20］，二者计算公式分别为

式中：ssn为降深计算值；son为降深观测值；N为观测数据的数量。

理论值与实测值对比检验结果为BLAS＝0.17、SNSE＝0.54，表明理论值不能满足工程精度［20］要求；实测值与计算值对比检验结果为BLAS＝0.08、SNSE＝0.86，表明计算值可以满足工程精度要求。本组对比结果再次表明，Dupuit 理论曲线即便在较远处与实测值仍有一定误差，直接使用Dupuit 公式需谨慎，这一结论与大量前人研究结果相吻合［21-24］。