利用变分迭代方法的新排水间距公式

《Irrigation and Drainage》：New Drain Spacing Formulas Using the Variational Iteration Method

　　**摘要**
本研究的目的是利用变分迭代方法（VIM）对线性化的Boussinesq方程进行求解，以计算排水管间距。通过在三种不同的初始条件情景下最多应用两次VIM方法，得出了三个用于计算排水管间距的方程。与Glover–Dumm方程相比，这些方程预测的排水管间距值通常更接近van Schilfgaarde方程的结果。具体来说，当参数λ的值在0.42–0.9范围内时，其中一个方程与van Schilfgaarde方程的相对误差（RE）小于10%；而当λ在0.51到0.75之间时，相对误差小于5%。与Glover–Dumm方程的比较表明，在λ值的范围更窄的情况下也能得到相同的相对误差值。在m值较大的情况下，地下水位下降的速度低于van Schilfgaarde方程的预测值；而在后期阶段，这一趋势则相反。

**引言**
农业用地的排水涉及人工去除多余的地表水和地下水，以实现适当的土壤湿度、通气和温度，特别是在根区，从而促进作物生长并保持土壤生产力（ICID 1978；Ziccarelli和Valore 2018）。不正确的灌溉水管理可能导致土壤盐碱化问题，以及灌溉农业中地下水位控制不足（Sallam 2017）。地下排水设计的一个关键方面是计算最佳的排水管间距，以确保有效的地下水位控制和盐分去除。用于确定排水管间距的方程通常基于Dupuit–Forchheimer假设，这些假设将流动域简化为具有平行水平流线的一维系统。这些方程主要分为两类：一类来自稳态理论，另一类基于非稳态（瞬态）理论。根据Ritzema（2007）的观点，稳态方法仅描述了地下水位深度与流量之间的简化、恒定关系。Hooghoudt（1940）被认为是稳态排水理论的先驱（Tzimopoulos等人2023）。他假设流动域由两个区域组成：一个是水流主要为水平的区域，另一个是靠近排水管的区域，那里会发生径向流动。通过引入等效深度的概念，Hooghoudt将水平流动和径向流动的结合转化为等效的水平流动。后来，研究人员开发了不同的稳态排水方程，这些方程要么包含了Dupuit–Forchheimer假设，要么对其进行了修正，排水管位于两层界面的不同深度（Kirkham 1958；Ernst 1962；Toksoz和Kirkham 1971；Terzides和Karamouzis 1986；Ritzema 1994；Smedema等人2004；Mishra和Singh 2007；Ostad-Ali-Askari和Shayannejad 2021）。然而，在实际中，地下水位补给随时间变化，导致流向排水管的水流速度不恒定。为了描述随时间变化的地下水位波动，采用了瞬态（非稳态）方法。稳态和瞬态方法都基于相同的假设，唯一的区别是在瞬态条件下补给量随时间变化。这种现象在灌溉区或大雨事件后尤为明显，此时地下水位会突然上升，然后在某个时间t逐渐恢复到灌溉前的水平。通常使用Glover–Dumm方程（Dumm 1954）来描述在排水管安装在不透水层上方后地下水位下降的情况，假设土壤是均匀的（图1）。

**排水系统符号说明**
用于描述排水到沟渠或管状排水系统的符号。两个排水管分别位于距离不透水层D的垂直距离x=0和x=L处。排水管轴线在x=L/2处的初始地下水位高度用m0表示，Z(x,t)表示地下水位距离不透水层的高度。Glover–Dumm方程（Dumm 1954）是基于Boussinesq方程（Boussinesq 1903；方程（1））的线性化，并考虑了Dupuit–Forchheimer假设推导出来的。这是一个二阶非线性偏微分方程，属于抛物线类型。其解利用了最初为热传导问题开发的方法（Carslaw和Jaeger 1959）。排水管间距的计算适用于以下情况：在某一时刻，地下水位表面与排水管水平平行，或者初始剖面呈四次抛物线形状。美国垦务局（US Bureau of Reclamation，Dumm 1960；USBR 1993；Ayars和Evans 2015）使用Glover–Dumm方程。Van Schilfgaarde（1963, 1964）在不进行线性化的情况下得到了Boussinesq方程的解（方程（1））。针对特定排水问题，不同研究人员提供了在不同初始和边界条件下的方程（1）的解析解（Terzidis 1968；Fath-Allah 2025）。除了解析解之外，还开发了数值近似方法，可以在任意初始和边界条件下以及对于非均匀流动域得到非线性微分方程的近似解（Tzimopoulos和Terzidis 1975；Skaggs和Tang 1979；Fipps等人1986）。联合国粮农组织（FAO，Van der Molen等人2007）提供了广泛采用的排水管间距计算软件，如DRAINMOD（Skaggs 1999）和Modflow（McDonald和Harbaugh 1988）。然而，在某些情况下（例如Modflow），需要在整个域和边界条件下详细定义初始地下水位，这使得其在实际排水设计中的应用具有挑战性（Tzimopoulos等人2023）。因此，实践者和排水工程师继续广泛使用经典模型，特别是Glover–Dumm方程。在本研究中，我们提出了使用变分迭代方法（VIM），这是一种基于拉格朗日乘数原理的相对较新的方法，最初由He（1997, 1998）引入。这种方法是一种用于求解线性和非线性微分方程的替代迭代技术（He和Wu 2007）。VIM因其简单性和能够在不需要离散化的情况下为广泛问题提供准确的近似解而受到广泛接受。该方法已成功应用于多个领域，如流体动力学（Liu和Kurra 2011）、热传递（Slota 2009）、波传播（Hemeda 2008）和多孔介质流动（Kargas和Mindrinos 2024）。然而，据我们所知，它尚未在排水研究中应用，特别是在求解用于确定排水管间距的Boussinesq方程方面。本研究的目的是将VIM应用于线性化的Boussinesq方程，以在三种不同的初始条件下计算均匀土壤的排水管间距。我们假设初始条件可以由以下形式描述：（i）具有椭圆形状的三角函数；（ii）具有抛物线形状的四次多项式函数；（iii）具有椭圆形状的幂律函数。我们还将结果与经典的Glover–Dumm和Van Schilfgaarde方程进行了比较。

**问题表述与迭代求解**
排水问题通过一维Boussinesq方程（Boussinesq 1903）进行建模，这是一个基于Dupuit–Forchheimer假设的非线性偏微分方程（PDE）：
(1)
其中，h是地下水位距离不透水层的距离，k是饱和水力传导率，γ是均匀土壤的特定产出率。方程（1）结合以下初始和边界条件进行考虑：
(2a)
(2b)
(2c)
其中d表示排水管距离不透水层的深度，分别位于x1和x2处。Boussinesq（Boussinesq 1903）在没有线性化的情况下求解了方程（1），假设其变量是可分离的。函数f可以是常数（本工作中未探讨）、三角函数或多项式。随后考虑了三种不同的情况，并检验了所提出方法的适用性。方程（1）可以重写为：
(3)

如果f相对于其幅度的变化很小，即Δf与h相比很小，则第一项分母中的f可以视为在平均值h处的常数。此外，如果f相对较大且斜率Δf/h很小，则可以忽略第二项。在这些假设下，我们得到以下形式的线性扩散型PDE：
(4)
其中h是排水管轴线处的地下水位高度。在一维区间[0, L]上，方程（4）可以使用多种经典理论和现代数值技术进行求解。解析方法包括分离变量法和傅里叶级数法（适用于均匀边界条件，Carslaw和Jaeger 1959），特征函数展开法（适用于更一般的不均匀问题，Strauss 2007），以及用于时变强迫的格林函数方法（Debnath 2016）。数值离散化方案包括显式有限差分法、隐式方案和Crank–Nicolson方法，后者在时间和空间上都能达到二阶精度（Smith 1985；Morton和Mayers 2005）。我们建议使用VIM来求解方程（4），该方法以级数形式提供解，并在几次迭代后收敛。有关该方法的基本原理，请参见He（1997, 1998）；关于其在渗透问题中的最新应用，请参见Kargas和Mindrinos（2024）。使用VIM对方程（4）的修正函数表示为：
(5)
其中λ是一阶时间导数的拉格朗日乘数，由Wazwaz（2014）给出。我们设置λ的值，具体将在后面指定。

**2.1 三角函数**
我们考虑以下形式的三角函数：
(6)
该函数满足方程（2）中的边界条件，并且满足其他条件，从而得到所需的解。鉴于在边界处具有均匀Dirichlet条件的热方程解的级数表示（Fokas和Kaxiras 2023），这是一个自然的选择。这个初始条件描述了一个最大高度为H的椭圆形地下水位。然后，从方程（5）出发，我们得到：
执行第二步并使用上述公式，我们得到：

**2.2 四次函数**
设f(x)满足方程（2），则该四次抛物线由Dumm（1968）提出。应用VIM（见方程（5）后，经过两步得到：

**2.3 幂函数**
受Tzimopoulos等人（2025）及其参考文献的启发，我们考虑以下平方根加权的多项式（幂律函数）：
(10)
其中系数选择满足条件（2）。初始函数方程（10）描述了一个最大高度为H的椭圆形地下水位。为了得到解的紧凑表示，我们仅应用VIM的一次迭代，得到：

**2.4 结果比较**
我们使用方程（6）得到第二步后的近似解，其最终形式为：

我们观察到该解满足方程（2）中提出的初始和边界条件：

**2.5 四次多项式的进一步处理**
设f(x)满足方程（2），则该四次抛物线由Dumm（1968）提出。应用VIM（见方程（5）后，经过两步得到：

我们用φ表示经过一步迭代后的近似解，其形式为：

我们观察到该迭代解不满足方程（2）中的边界条件：

**结论**
方程（7）表示的解是线性化方程（4）对于形式为方程（6）的初始条件的近似解。这个选择是基于以下部分的目标，其中需要一个简单的显式表达式。我们观察到，对于某些情况，近似解是没有定义的，但这不是一个实际问题，因为解关于某个点是对称的，这意味着它满足对称性质；因此，我们可以考虑该函数在该点及其反射点的值。在边界条件和齐次诺伊曼条件下的情况，见方程（2c），并没有得到满足，如下所示：
（12）这意味着随着时间的增加，我们在该点的值会发生变化。

3 确定排水间距的新公式
与排水问题相关的有两个有趣的逆问题：

估计在给定时间内，由于地下水位（高度）在中点下降，排水管之间的间距。根据排水管之间的距离，估计中间高度下降所需的时间。第一个问题在数学上表述如下：找到满足以下条件的解
（13）对于某个给定的时间和高度。

3.1 三角函数
对于形式如方程（6）的初始函数，由方程（7）和（13）给出的解可以表示为
（14）
我们定义，并且方程（14）可以表示为
这是一个关于的二次方程，其解为
或者等价地
（15）这只有在中间地下水位的高度大于其初始值的一半时才有意义。鉴于此，两个解都是正的。接下来，我们展示其中只有一个解在物理上是有意义的。如果，我们得到一个唯一的正解，形式为

第二个逆问题也是通过方程（13）来表述的，现在已知，而（用替换）是未知函数，从而得到
（16）
这是一个关于的二次方程，其解为
（17）
这里同样要求。

3.2 四次函数
作为另一种方法，可以使用方程（9）来求解逆问题（13），对于初始函数（8），我们得到唯一的正解
（21a）
（21b）

3.3 幂函数
同样地，我们可以使用方程（11）中的表达式来求解逆问题（13），基于初始函数（10），其形式为：
（22a）
（22b）

4 近似解析解
线性化的方程（4）也可以通过分离变量来求解，得到级数表示（Carslaw和Jaeger 1959）。Dumm（1964）推导出了解，其形式为
（23），假设并且初始函数的形式为（8）。如果，级数表示中的第二项及后续项会变得非常小，可以忽略不计。因此，只考虑方程（23）中的第一项并求解，得到
（24），这被称为Glover–Dumm方程。注意，当差异（）很大时，由给出。除了之前提到的需要用替换之外，还应强调，当差异变大时，假设常数是值得怀疑的。此外，Van Schilfgaarde提出了非线性方程（1）的解（Van Schilfgaarde 1963, 1964）。初始条件是一个最大高度为的椭圆形地下水位。提出的计算排水间距的公式为
（25）

方程（24）和（25）分别考虑了线性化的方程（4）和非线性Boussinesq方程（1），将在下一节中用于与我们的公式进行比较。我们观察到，在所有公式中都出现了项。由于，我们设置，其中，我们定义比率
其中 trig, qua, pol 分别对应于方程（20a）、（21a）和（22a），以及（24）和（25），以便它可以表示为的函数。我们得到
（26a）
（26b）
（26c）
（26d）
（26e）
（26f）

5 数值实验
在本节中，将展示所提出的方程（20a）、（21a）和（22a）的结果，并将其与广泛使用的Glover–Dumm方程（24）和van Schilfgaarde方程（25）的结果进行比较，使用的数据来自Van Beers（1976）的非稳态流动示例。目的是为了相对全面地评估它们在这个示例中的性能，并说明它们的行为。然而，将结果与不同条件下的实验现场数据进行比较以得出更安全的结论的问题仍然悬而未决。排水管放置在土壤表面下2米的深度，不透水的土壤层位于7米的深度，这给出了。均匀土壤的水力传导率为，特定产率为。假设地下水位从土壤表面的初始距离为0.80米，结果为

在表1中，根据方程（20a）、（21a）、（22a）、（24）和（25）计算了在时间和10天，以及0.5米时的值。表1. 对应时间和的值的，以及方程（20a）、（21a）、（22a）、（24）和（25）的相对误差（RE）值，以及与方程（24）和（25）的比较。

从上述结果来看，方程（20a）、（21a）和（22a）通常预测的值比方程（24）的值要大，除了在方程（20a）的情况下。相比之下，方程（20a）、（21a）和（22a）与方程（25）相比，并不总是预测更大的值，如表1所示。然而，很明显，这些方程与方程（25）的偏差比方程（24）的偏差要小得多。为了更好地评估方程（20a）、（21a）和（22a）与方程（24）和（25）之间的差异，我们计算了相对误差值（见表1）。我们定义相对误差：

这些方程的行为对于值和0.5米以及固定值进行了检验（表1）。方程（22a）产生的值的相对误差小于方程（24），对于所有选定的值。方程（20a）只有在时才产生较小的误差（8%），而它不允许计算的值。方程（21a）仅在时产生相对较小的RE值。观察到方程（22a）的RE值对于所有三个值都相对恒定，而对于方程（21a），随着的减小而增加，而对于方程（20a），随着的减小而减小，直到值。应该注意的是，当差异（）变得足够大时，即，RE值会显著增加。这可能是因为在时，地下水位与不透水层的距离不能被认为是恒定的，这违反了与相比较小的假设。还应该注意的是，如果改变每个的（排水管与不透水层的距离）的值，同时保持其他参数不变，那么与方程（24）相比，三个方程的RE值保持不变。当改变或的值同时保持其他参数不变时，也会出现同样的情况。这是显而易见的，因为项在所有方程中都是共有的。如方程（26a）–（26c）所示，三个方程与方程（24）之间的差异仅取决于参数的值。如果应用更多的VIM迭代，RE值可能会减小。然而，在这种情况下，方程会变得更加复杂，应用也会更加困难。值得注意的是，方程（22a）仅进行一次VIM迭代，总体上给出了最佳结果。为了更全面地了解方程的行为，图2展示了三个方程（20a）、（21a）和（22a）以及方程（24）和（25）的预测值作为时间（天）的函数，对于和0.50的值。在纵轴上，使用米来消除土壤特性的影响（Van Schilfgaarde 1974）。从结果来看，在三种情况中的两种情况下，方程（24）的曲线位于其他方程曲线的左侧。这意味着对于相同的值，方程（24）预测地下水位下降到和0.50所需的时间更长。实际上，这表明与其它方程预测的相比，地下水位下降的速度较慢。特别是，方程（22a）在所有三种情况下都与方程（24）的行为一致，而方程（21a）的偏差随着的增加而增加。相比之下，方程（20a）预测的偏差随着的增加而减小。图2在图形查看器中打开

图2显示了使用方程（20a）、（21a）、（22a）、（24）和（25）的三个值和的与排水间距之间的关系。此外，表1展示了方程（20a）、（21a）和（22a）与方程（25）的RE值对于相同的值。然而，应该注意的是，在这种情况下，值不仅受影响，还受值的影响。对于特定的值，与方程（24）相比，RE值要小得多，除了在时的方程（20a）。在三个方程中，方程（22a）产生的RE值最低，且始终小于。最大的RE值再次出现在时的方程（21a）中。因此，方程（22a）的预测更接近方程（25）的预测。关于对RE的影响，还检查了和的额外情景。表2中的结果表明，对于相同的值，方程（22a）的RE随着的增加而减小。具体来说，对于，对于，以及对于。总体而言，的变化对的影响似乎并不显著。表2显示了使用方程（20a）、（21a）、（22a）和（25）的对应时间和的值的值，以及方程（20a）、（21a）和（22a）的相对误差（RE）值与方程（25）的比较。

如图2所示，方程（22a）在所有情况下都更接近方程（25）。然而，观察到对于相同的值，在较高的时，地下水位下降的速度低于方程（25），而在后期阶段，趋势则相反。总体而言，如表3所示，对于参数值，方程（22a）和方程（25）之间的RE小于10%，而对于方程（20a）和方程（21a），相应的参数范围分别是和。这表明方程（20a）和（21a）在达到可接受范围内的值范围较窄。表3。在方程（20a）、（21a）和（22a）中，达到最小RE的参数值与方程（24）和（25）进行了比较，同时也确定了RE偏差发生的参数值范围。

最小值 5%–10%
最小值 5%–10%
0%–5%

方程（20a） 0.50 7.5% [0.50, 0.59]
0.56 0.06%
0.51 [0.52, 0.90]

方程（21a） 0.87 8.0% [0.80, 0.90]
0.71 2.0% [0.67, 0.74]

方程（22a） 0.60 9.0% [0.50, 0.68]
0.62 0.03% [0.42, 0.50] [0.51, 0.75]

与方程（24）相比，这三个方程对应的参数值范围要窄得多，表明新方程，特别是方程（22a），更接近方程（25）。然而，公式（22a）源自方程（11），而方程（11）未能满足边界条件；参见方程（12）。Van Schilfgaarde（1974）报告称，方程（24）通常低估了该参数的值，而方程（25）则倾向于高估它。这实际上可能意味着所提出的方程，特别是方程（22a），能够预测出更接近他方程的值，这一点通过具体例子的分析得到了证实。因此，基于线性化Boussinesq方程的VIM解得到的方程（22a），其预测的值与van Schilfgaarde得到的非线性解相似。需要注意的是，在计算每个方程中的实际参数值时，必须考虑由于流线在排水口附近汇聚而造成的额外能量损失，这需要通过在相应方程中引入Hooghout的等效深度来替代原来的参数。此外，基于与方程（24）相同初始条件的方程（21a）出现较大的RE值，需要进一步研究。如上所述，通过应用超过两次VIM迭代可以改善这一问题。

6 结论

排水系统在农业生产中至关重要，尤其是在排水条件较差的地区。本研究尝试使用相对较新的VIM方法，从线性化Boussinesq方程计算排水间距。该方法应用于三种不同的初始条件情景，得到了三个不同的排水间距方程，每个方程最多进行了两次VIM迭代。然而，其中两个方程（21a）和（22a）的解不满足边界条件之一，而满足边界条件的第三个方程（20a）仅在某个参数值范围内有效。尽管存在这些限制，这些结果仍为地下排水系统设计、地下水控制以及洪水风险的有效管理提供了宝贵的见解。将这三个方程的结果与著名的Glover–Dumm方程和van Schilfgaarde方程进行比较后发现，总体而言，这三个方程预测的排水间距值都更接近van Schilfgaarde方程的结果。特别是具有幂律函数（平方根加权多项式）初始条件的方程，在较宽的参数值范围内，其预测的值与van Schilfgaarde方程的结果相似。此外，这三个方程预测的地下水位下降速度都比Glover–Dumm方程快；而与van Schilfgaarde方程相比，情况更为复杂：在较大的参数值下，地下水位下降速度低于方程（25）的预测值，但在后期阶段，趋势则相反。通过增加VIM迭代次数，这些方程的预测性能可能会得到提高；然而，由于方程的复杂性，应用起来会更加困难。为了对这三个方程得出可靠的结论，应将其结果与实际现场数据以及其他广泛使用的排水间距计算方程进行比较。未来的研究应侧重于长期监测和评估，以验证这些结果并进一步改进排水设计。此外，还需要进行现场研究，以确定哪些参数对设计最为关键，以及设计过程中每个参数所需的置信区间。

本文以开放获取（OA）模式发表，得到了HEAL-Link的财政支持。

作者声明没有利益冲突。

数据可应要求向作者索取。