CN106055774A - 地面灌溉中地表水流运动的模拟方法 - Google Patents

Abstract

本发明公开了一种地面灌溉中地表水流运动的模拟方法，属于地面灌溉领域。该模拟方法包括：首先，选择向量形式的全水动力学方程作为模拟目标畦田地面灌溉中地表水流运动的地面灌溉模型。然后，对向量形式的全水动力学方程进行时空离散化处理，得到全水动力学方程的时空离散表达式；并将时空离散表达式中的和分别用q和h来表达。最后，将全水动力学方程的时空离散表达式整理成以q和h为未知数的地表水流运动模拟方程，通过求解该方程，获取目标畦田中每一测点处的q和h，从而对目标畦田地面灌溉中地表水流运动进行模拟。该模拟方法计算简便，有效缩短了模拟计算的时间，且计算精度高，适用于大规模地面灌溉中地表水流运动的模拟。

Description

地面灌溉中地表水流运动的模拟方法

技术领域

本发明涉及地面灌溉领域，特别涉及一种地面灌溉中地表水流运动的模拟方法。

背景技术

地面灌溉是一种常见的农田灌溉方法，其利用灌水沟、畦或格田等，将水引入农田，并使水在重力和毛细管作用下渗入土壤，为田间作物补充水分。进行地面灌溉时，为保证田间作物得到充分灌溉，并尽可能节约水资源，需要对水流在目标畦田中的灌溉性能进行预先评价，以对灌入水量及流速等进行调整，避免水资源的浪费。因此，提供一种地面灌溉的评价方法是十分必要的。

现有技术提供了一种地面灌溉中地表水流运动的模拟方法，该方法采用具有向量形式的全水动力学方程作为地面灌溉模型，根据此来描述地面灌溉中地表水流运动过程。在求解具有向量形式的全水动力学方程时，首先对其进行时空离散化处理，得到时空离散表达式：

$U_i^{n_t+1}=U_i^{n_t}+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(F_{i+1/2}^{n_t+1}-F_{i-1/2}^{n_t+1})+{\mathrm{\ΔtS}}_i^{n_t}$

式中，F为数值通量，U为因变量向量，S为源项向量；Δt和Δx分别为数值计算的时间步长和空间步长；n_t为时间离散迭代步；i为空间离散节点；i+1/2为空间离散节点i与i+1之间的中点；i-1/2为空间离散节点i与i-1之间的中点。

其中，在求解全水动力学方程的上述时空离散表达式时，关键是计算和现有技术通常采用经典Roe模拟方法进行计算，该经典Roe模拟方法的表达式如下所示：

$F_{i+1/2}^{n_t+1}=\frac{1}{2}\[F_{i+1}^{n_t}+F_i^{n_t}-\left(\left|{\mathrm{\λ}}_{1,i+1/2}^{n_t}\right|{\mathrm{\α}}_{1,i+1/2}^{n_t}{e}_{1,i+1/2}^{n_t}+\left|{\mathrm{\λ}}_{2,i+1/2}^{n_t}\right|{\mathrm{\α}}_{2,i+1/2}^{n_t}{e}_{2,i+1/2}^{n_t}\right)\]$

$F_{i-1/2}^{n_t+1}=\frac{1}{2}\[F_{i-1}^{n_t}+F_i^{n_t}-\left(\right|{\mathrm{\λ}}_{1,i-1/2}^{n_t}|{\mathrm{\α}}_{1,i-1/2}^{n_t}{e}_{1,i-1/2}^{n_t}+|{\mathrm{\λ}}_{2,i-1/2}^{n_t}\left|{\mathrm{\α}}_{2,i-1/2}^{n_t}{e}_{2,i-1/2}^{n_t}\right)\]$

其中，

${\mathrm{\λ}}_{2,i+1/2}^{n_t}=\frac{\sqrt{h_i^{n_t}}{u}_i^{n_t}+\sqrt{h_{i+1}^{n_t}}{u}_{i+1}^{n_t}}{\sqrt{h_i^{n_t}}+\sqrt{h_{i+1}^{n_t}}}-\sqrt{\frac{g}{2}(h_i^{n_t}+h_{i+1}^{n_t})}$

$e_{1,i+1/2}^{n_t}=\left(\begin{array}{c}1\\ \frac{\sqrt{h_i^{n_t}}{u}_i^{n_t}+\sqrt{h_{i+1}^{n_t}}{u}_{i+1}^{n_t}}{\sqrt{h_i^{n_t}}+\sqrt{h_{i+1}^{n_t}}}+\sqrt{\frac{g}{2}(h_i^{n_t}+h_{i+1}^{n_t})}\end{array}\right)$

$e_{2,i+1/2}^{n_t}=\left(\begin{array}{c}1\\ \frac{\sqrt{h_i^{n_t}}{u}_i^{n_t}+\sqrt{h_{i+1}^{n_t}}{u}_{i+1}^{n_t}}{\sqrt{h_i^{n_t}}+\sqrt{h_{i+1}^{n_t}}}+\sqrt{\frac{g}{2}(h_i^{n_t}+h_{i+1}^{n_t})}\end{array}\right)$

$\begin{array}{c}{\mathrm{\α}}_{2,i+1/2}^{n_t}=\frac{1}{2}(h_{i+1}^{n_t}-h_i^{n_t})\\ +\frac{1}{2\sqrt{g\left(h_i^{n_t}+h_{i+1}^{n_t}\right)/2}}\[\left(h_{i+1}^{n_t}{u}_{i+1}^{n_t}-h_i^{n_t}{u}_i^{n_t}\right)-\frac{\sqrt{h_i^{n_t}}{u}_i^{n_t}+\sqrt{h_{i+1}^{n_t}}{u}_{i+1}^{n_t}}{\sqrt{h_i^{n_t}}+\sqrt{h_{i+1}^{n_t}}}(h_{i+1}^{n_t}-h_i^{n_t})\]\\ {\mathrm{\α}}_{2,i+1/2}^{n_t}=\frac{1}{2}(h_{i+1}^{n_t}-h_i^{n_t})\\ -\frac{1}{2\sqrt{g\left(h_i^{n_t}+h_{i+1}^{n_t}\right)/2}}\[\left(h_{i+1}^{n_t}{u}_{i+1}^{n_t}-h_i^{n_t}{u}_i^{n_t}\right)-\frac{\sqrt{h_i^{n_t}}{u}_i^{n_t}+\sqrt{h_{i+1}^{n_t}}{u}_{i+1}^{n_t}}{\sqrt{h_i^{n_t}}+\sqrt{h_{i+1}^{n_t}}}(h_{i+1}^{n_t}-h_i^{n_t})\]\end{array}$

${\mathrm{\λ}}_{1,i-1/2}^{n_t}=\frac{\sqrt{h_{i-1}^{n_t}}{u}_{i-1}^{n_t}+\sqrt{h_i^{n_t}}{u}_i^{n_t}}{\sqrt{h_{i-1}^{n_t}}+\sqrt{h_i^{n_t}}}+\sqrt{\frac{g}{2}\left(h_{i-1}^{n_t}+h_i^{n_t}\right)}$

${\mathrm{\λ}}_{2,i-1/2}^{n_t}=\frac{\sqrt{h_{i-1}^{n_t}}{u}_{i-1}^{n_t}+\sqrt{h_i^{n_t}}{u}_i^{n_t}}{\sqrt{h_{i-1}^{n_t}}+\sqrt{h_i^{n_t}}}+\sqrt{\frac{g}{2}\left(h_{i-1}^{n_t}+h_i^{n_t}\right)}$

$e_{1,i-1/2}^{n_t}=\left(\begin{array}{c}1\\ \frac{\sqrt{h_{i-1}^{n_t}}{u}_{i-1}^{n_t}+\sqrt{h_i^{n_t}}{u}_i^{n_t}}{\sqrt{h_{i-1}^{n_t}}+\sqrt{h_i^{n_t}}}+\sqrt{\frac{g}{2}(h_{i-1}^{n_t}+h_i^{n_t})}\end{array}\right)$

$e_{2,i-1/2}^{n_t}=\left(\begin{array}{c}1\\ \frac{\sqrt{h_{i-1}^{n_t}}{u}_{i-1}^{n_t}+\sqrt{h_i^{n_t}}{u}_i^{n_t}}{\sqrt{h_{i-1}^{n_t}}+\sqrt{h_i^{n_t}}}+\sqrt{\frac{g}{2}(h_{i-1}^{n_t}+h_i^{n_t})}\end{array}\right)$

$\begin{array}{c}{\mathrm{\α}}_{1,i-1/2}^{n_t}=\frac{1}{2}(h_i^{n_t}-h_{i-1}^{n_t})\\ +\frac{1}{2\sqrt{g\left(h_{i-1}^{n_t}+h_i^{n_t}\right)/2}}\[\left(h_i^{n_t}{u}_i^{n_t}-h_{i-1}^{n_t}{u}_{i-1}^{n_t}\right)-\frac{\sqrt{h_{i-1}^{n_t}}{u}_{i-1}^{n_t}+\sqrt{h_i^{n_t}}{u}_{i1}^{n_t}}{\sqrt{h_{i-1}^{n_t}}+\sqrt{h_i^{n_t}}}(h_i^{n_t}-h_{i-1}^{n_t})\]\end{array}$

$\begin{array}{c}{\mathrm{\α}}_{2,i-1/2}^{n_t}=\frac{1}{2}(h_i^{n_t}-h_{i-1}^{n_t})\\ -\frac{1}{2\sqrt{g\left(h_{i-1}^{n_t}+h_i^{n_t}\right)/2}}\[\left(h_i^{n_t}{u}_i^{n_t}-h_{i-1}^{n_t}{u}_{i-1}^{n_t}\right)-\frac{\sqrt{h_{i-1}^{n_t}}{u}_{i-1}^{n_t}+\sqrt{h_i^{n_t}}{u}_{i1}^{n_t}}{\sqrt{h_{i-1}^{n_t}}+\sqrt{h_i^{n_t}}}(h_i^{n_t}-h_{i-1}^{n_t})\]\end{array}$

q＝hu

其中，u为水流沿所述目标畦田畦长方向的垂向均布流速，单位为m/s；h为地表水深，单位为m；q为水流沿所述目标畦田畦长方向的单宽流量，单位为m³/(s·m)。

根据上述各计算公式能够求解得到关于q和h的方程，通过该方程绘制h相对于q的变化曲线，即获得地表水深相对于水流沿目标畦田畦长方向的单宽流量变化情况，以完成地面灌溉中地表水流运动的模拟，从而实现对地面灌溉的评价。

发明人发现现有技术至少存在以下技术问题：

采用现有技术对地面灌溉中地表水流运动进行模拟计算时，虽然经典Roe模拟方法的模拟精度较高，但其表达式极为复杂，导致求解效率较低，无法有效完成对大规模地面灌溉水流运动的模拟。

发明内容

本发明实施例所要解决的技术问题在于，提供了一种计算简单快捷、且精度高的地面灌溉中地表水流运动的模拟方法，具体技术方案如下：

本发明实施例提供了一种地面灌溉中地表水流运动的模拟方法，该模拟方法包括以下步骤：

步骤a、选择向量形式的全水动力学方程作为模拟目标畦田地面灌溉中地表水流运动的地面灌溉模型。

所述向量形式的全水动力学方程为：

$\frac{\∂U}{\∂t}+\frac{\∂F}{\∂x}=S$

$U=\left(\begin{array}{c}h\\ q\end{array}\right)$

$F=\left(\begin{array}{c}q\\ qu+\frac{1}{2}{\mathrm{gh}}^2\end{array}\right)$

$S=\left(\begin{array}{c}{i}_c\\ -gh\frac{\∂b}{\∂x}-gh\frac{n^{2}u\left|u\right|}{h^{4/3}}\end{array}\right)$

其中，F为物理通量，U为因变量向量，S为源项向量，t为所述目标畦田以灌溉起始时刻作为原点的时间坐标，单位为s，x为所述目标畦田中任意测点与所述目标畦田的畦首之间的距离，单位为m；u为水流沿所述目标畦田畦长方向的垂向均布流速，单位为m/s；h为地表水深，单位为m；q为水流沿所述目标畦田畦长方向的单宽流量，单位为m³/(s·m)，且q＝hu；g为重力加速度，单位为m/s²；i_c为实测得到的所述目标畦田地表水入渗率；b为实测得到的所述目标畦田地表相对高程，单位为m；n为已知的所述目标畦田的畦面糙率系数，单位为s/m^1/3。

步骤b、对所述向量形式的全水动力学方程进行时空离散化处理，得到所述全水动力学方程的时空离散表达式。

所述全水动力学方程的时空离散表达式为：

$U_i^{n_t+1}=U_i^{n_t}+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(F_{i+1/2}^{n_t+1}-F_{i-1/2}^{n_t+1})+{\mathrm{\ΔtS}}_i^{n_t};$

其中，Δt为针对所述t得到的时间步长；Δx为针对所述x得到的空间步长；n_t为针对所述t得到的时间离散节点；i为针对所述x得到的空间离散节点；i+1/2为空间离散节点i与i+1之间的中点；i-1/2为空间离散节点i与i-1之间的中点。

步骤c、提供利用所述q和所述h来表达所述的第一表达式。

所述第一表达式的为：

$S_i^{n_t}=\left(\begin{array}{c}{\left(i_c\right)}_i^{n_t}\\ -{\mathrm{gh}}_i^{n_t}\frac{b_{i+1}-b_{i-1}}{2\mathrm{\Δ}x}-g\frac{n^2{u}_i^{n_t}\left|u_i^{n_t}\right|}{{\left(h^{1/3}\right)}_i^{n_t}}\end{array}\right).$

步骤d、分别提供利用所述q和所述h来表达所述和的第二表达式和第三表达式。

所述第二表达式为：

$F_{i+1/2}^{n_t+1}=0.5\[F_{i+1}^{n_t}+F_i^{n_t}-0.5{\mathrm{\Λ}}_{i+1/2}^{n_t}(U_{i+1}^{n_t+1}-U_i^{n_t})\]$

其中，

所述第三表达式为：

$F_{i-1/2}^{n_t+1}=0.5\[F_i^{n_t}+F_{i-1}^{n_t}-0.5{\mathrm{\Λ}}_{i-1/2}^{n_t}(U_i^{n_t+1}-U_{i-1}^{n_t})\]$

${\mathrm{\Λ}}_{i-1/2}^{n_t}=\left(\begin{array}{cc}0.5\left(u_{i-1}^{n_t}+u_i^{n_t}-\sqrt{{\mathrm{gh}}_{i-1}^{n_t}}-\sqrt{{\mathrm{gh}}_i^{n_t}}\right)& 0\\ 0& 0.5\left(u_{i-1}^{n_t}+u_i^{n_t}+\sqrt{{\mathrm{gh}}_{i-1}^{n_t}}+\sqrt{{\mathrm{gh}}_i^{n_t}}\right)\end{array}\right).$

步骤e、将所述第一表达式、所述第二表达式和所述第三表达式代入所述全水动力学方程的时空离散表达式，得到以所述q和所述h作为未知数的地表水流运动的模拟方程，通过对所述地表水流运动的模拟方程进行求解，获取所述目标畦田中每一测点处的所述q和所述h，从而对地面灌溉中地表水流运动进行模拟。

具体地，作为优选，所述目标畦田地表水入渗率i_c通过Kostiakov入渗公式计算得到，所述Kostiakov入渗公式为：

$i_c=k_{in}{\mathrm{\α\τ}}_{in}^{\mathrm{\α}-1}$

其中，k_in为土壤入渗参数，单位为m/s^α；α为无量纲入渗指数参数；τ_in为入渗受水时间，单位为s。

具体地，作为优选，所述畦面糙率系数n通过在所述目标畦田的畦首实测畦首地表水深值，并根据所述畦首地表水深值估算得到。

具体地，作为优选，所述时间步长为1-10s。

具体地，作为优选，所述空间步长为1-10m。

具体地，作为优选，在所述步骤e中，在求解所述地表水流运动的模拟方程时，将水流首次进入畦田时的单宽流量作为求解所述q的初始条件。

本发明实施例提供的技术方案带来的有益效果是：

本发明实施例提供的地面灌溉中地表水流运动的模拟方法，通过对向量形式的全水动力学方程进行时空离散，获取该全水动力学方程的时空离散表达式，并将该表达式转换成以单宽流量q和地表水深h作为未知数的地表水流运动的模拟方程。通过求解该模拟方程，即可获取目标畦田中每一测点处的q和h，进而可获得目标畦田中所有测点处的q和h，进而通过所求解得到的数据来对地面灌溉中地表水流运动进行模拟。本发明实施例所使用的上述时空离散表达式简便可靠，能在保证计算精确度的前提下，有效缩短q和h的计算时间，提高了地面灌溉中地表水流运动模拟数据的计算效率，而且，该方法能够适用于农田等流速较缓的地面灌溉中地表水流运动的模拟，并可用于大规模地面灌溉中地表水流运动的模拟。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例提供的地面灌溉中某一时刻下地表水流运动的状态示意图；

图2是本发明又一实施例提供的河北冶河灌区、灌溉时间为40min时的地面灌溉中地表水流运动的模拟结果示意图；

图3是本发明又一实施例提供的河北冶河灌区、灌溉时间为80min时的地面灌溉中地表水流运动的模拟结果示意图；

图4是本发明又一实施例提供的新疆建设兵团222团、灌溉时间为40min时的地面灌溉中地表水流运动的模拟结果示意图；

图5是本发明又一实施例提供的新疆建设兵团222团、灌溉时间为80min时的地面灌溉中地表水流运动的模拟结果示意图；

图6是本发明又一实施例提供的山东簸箕李灌区、灌溉时间为40min时的地面灌溉中地表水流运动的模拟结果示意图；

图7是本发明又一实施例提供的山东簸箕李灌区、灌溉时间为80min时的地面灌溉中地表水流运动的模拟结果示意图。

附图标记分别表示：

h 地表水深，

d 入渗水深，

b 地表高程。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将对本发明实施方式作进一步地详细描述。

本发明实施例提供了一种地面灌溉中地表水流运动的模拟方法，该模拟方法包括以下步骤：

步骤101、选择向量形式的全水动力学方程作为模拟目标畦田地面灌溉中地表水流运动的地面灌溉模型。

向量形式的全水动力学方程为：

$\frac{\∂U}{\∂t}+\frac{\∂F}{\∂x}=S$

$U=\left(\begin{array}{c}h\\ q\end{array}\right)$

$F=\left(\begin{array}{c}q\\ qu+\frac{1}{2}{\mathrm{gh}}^2\end{array}\right)$

$S=\left(\begin{array}{c}{i}_c\\ -gh\frac{\∂b}{\∂x}-gh\frac{n^{2}u\left|u\right|}{h^{4/3}}\end{array}\right)$

其中，F为物理通量，U为因变量向量，S为源项向量，t为目标畦田以灌溉起始时刻作为原点的时间坐标，单位为s，x为目标畦田中任意测点与目标畦田的畦首之间的距离，单位为m；u为水流沿目标畦田畦长方向的垂向均布流速，单位为m/s；h为地表水深，单位为m；q为水流沿目标畦田畦长方向的单宽流量，单位为m³/(s·m)，且q＝hu；g为重力加速度，单位为m/s²；i_c为实测得到的目标畦田地表水入渗率；b为实测得到的目标畦田地表相对高程，单位为m；n为已知的目标畦田的畦面糙率系数，单位为s/m^1/3。

步骤102、对向量形式的全水动力学方程进行时空离散化处理，得到全水动力学方程的时空离散表达式。

全水动力学方程的时空离散表达式为：

$U_i^{n_t+1}=U_i^{n_t}+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(F_{i+1/2}^{n_t+1}-F_{i-1/2}^{n_t+1})+{\mathrm{\ΔtS}}_i^{n_t};$

其中，Δt为针对t得到的时间步长；Δx为针对x得到的空间步长；n_t为针对t得到的时间离散节点；i为针对x得到的空间离散节点；i+1/2为空间离散节点i与i+1之间的中点；i-1/2为空间离散节点i与i-1之间的中点；

步骤103、提供利用q和h来表达的第一表达式。

第一表达式的为：

$S_i^{n_t}=\left(\begin{array}{c}{\left(i_c\right)}_i^{n_t}\\ -{\mathrm{gh}}_i^{n_t}\frac{b_{i+1}-b_{i-1}}{2\mathrm{\Δ}x}-g\frac{n^2{u}_i^{n_t}\left|u_i^{n_t}\right|}{{\left(h^{1/3}\right)}_i^{n_t}}\end{array}\right).$

步骤104、分别提供利用q和h来表达和的第二表达式和第三表达式。

第二表达式为：

$F_{i+1/2}^{n_t+1}=0.5\[F_{i+1}^{n_t+1}+F_i^{n_t+1}-0.5{\mathrm{\Λ}}_{i+1/2}^{n_t+1}(U_{i+1}^{n_t+1}-U_i^{n_t+1})\]$

其中，

第三表达式为：

$F_{i-1/2}^{n_t+1}=0.5\[F_i^{n_t+1}+F_{i-1}^{n_t+1}-0.5{\mathrm{\Λ}}_{i-1/2}^{n_t+1}(U_i^{n_t+1}-U_{i-1}^{n_t+1})\]$

${\mathrm{\Λ}}_{i-1/2}^{n_t+1}=\left(\begin{array}{cc}0.5\left(u_{i-1}^{n_t+1}+u_i^{n_t+1}-\sqrt{{\mathrm{gh}}_{i-1}^{n_t+1}}-\sqrt{{\mathrm{gh}}_i^{n_t+1}}\right)& 0\\ 0& 0.5\left(u_{i-1}^{n_t+1}+u_i^{n_t+1}+\sqrt{{\mathrm{gh}}_{i-1}^{n_t+1}}+\sqrt{{\mathrm{gh}}_i^{n_t+1}}\right)\end{array}\right).$

步骤105、将第一表达式、第二表达式和第三表达式代入全水动力学方程的时空离散表达式，得到以q和h作为未知数的地表水流运动的模拟方程，通过对地表水流运动的模拟方程进行求解，获取目标畦田中每一测点处的q和h，从而对地面灌溉中地表水流运动进行模拟。

本发明实施例提供的地面灌溉中地表水流运动的模拟方法，通过对向量形式的全水动力学方程进行时空离散，获取该全水动力学方程的时空离散表达式，并将该表达式转换成以单宽流量q和地表水深h作为未知数的地表水流运动的模拟方程。通过求解该模拟方程，即可获取目标畦田中每一测点处的q和h，进而可获得目标畦田中所有测点处的q和h，进而通过所求解得到的数据来对地面灌溉中地表水流运动进行模拟。本发明实施例所使用的上述时空离散表达式简便可靠，能在保证计算精确度的前提下，有效缩短q和h的计算时间，提高了地面灌溉中地表水流运动模拟数据的计算效率，而且，该方法能够适用于农田等流速较缓的地面灌溉中地表水流运动的模拟，并可用于大规模地面灌溉中地表水流运动的模拟。

步骤101中，向量形式的全水动力学方程中的t为目标畦田以灌溉起始时刻作为原点的时间坐标，目标畦田空间内所有测点采用统一的时间坐标，计算时需要对t进行时间离散；x为目标畦田任意测点与目标畦田的畦首之间的距离，即将目标畦田的畦首位置定义为0m，则目标畦田中任意一个测点所在位置与畦首之间的距离则定义为xm，该距离可以通过实测得到。根据向量形式的全水动力学方程，对于因变量向量U求其在t轴方向上的偏导数，对于物理通量F求其在x轴方向上的偏导数，得到源项向量S，以对目标畦田中每个测点处的q和h相对于t和x的变化情况进行模拟计算。

可以理解的是，u为水流沿目标畦田畦长方向的垂向均布流速，其指的是与目标畦田畦长方向垂直的水流速度；如附图1所示，h为地表水深，即目标畦田中某个测点处土壤表面以上的灌溉水深度；q为水流沿目标畦田畦长方向的单宽流量，即目标畦田中单位面积和单位时间内流过的水的体积，因此，可以得出q＝hu。在向量形式的全水动力学方程中，上述q、h、u三个量均为未知量，其中u可用q和h来表达，进而可得出仅仅关于q和h的方程，以对地面灌溉中地表水流的运动进行模拟。

具体地，g为重力加速度，可以根据目标畦田所在的地理位置，确定重力加速度的值，一般选取9.8m/s²。

i_c为实测得到的目标畦田地表水入渗率，即目标畦田中地表水渗入土壤的速度，通过地表水入渗率，即可求出某一水流推进时间内，目标畦田内某测点处水渗入土壤的深度，即地表水的入渗水深d，其等于i_c乘以t，从而对地面灌溉的情况进行模拟和预测。如附图1所示，b为实测得到的目标畦田在某一测点处的地表高程b，即目标畦田相对于地表的高度。地表高程b可以采用水准仪在间隔预设距离上进行测量，其中预设距离可以为5m或者其他数值，例如4m、6m、8m、10m等，本领域技术人员可以根据实际情况进行设定。地表高程b可以作为地表水深h以及入渗水深d的基准线，便于更直观地得到以地表高程b、地表水深h、入渗水深d随畦田长度的变化关系来表示地面灌溉中地表水流运动的模拟曲线。

n为目标畦田的畦面糙率系数，在本发明实施例中，其作为已知值，通常情况下，畦面糙率系数多在0.04-0.3之间选择，其具体值可通过如下方法来确定：对目标畦田畦首处的地表水深h进行观测，获取目标畦田畦首处的地表水深h的观测数据，根据该观测数据和Chézy公式估算畦面糙率系数。本领域技术人员可以理解的是，目标畦田不同位置处的糙率系数可能不相同，可以以不同位置处的糙率系数的平均值作为该目标畦田的糙率系数。通过上述实测或估算得到的参数来辅助向量形式的全水动力学方程的求解，以实现地面灌溉中对地表水流的运动的模拟。

虽然在本发明实施例中，向量形式的全水动力学方程也是关于q和h的表达式，但是该表达式为一个连续方程，无法得到目标畦田中每一测点处的q和h。所以在求解每一测点处的q和h时，需要对向量形式的全水动力学方程进行离散，然后通过它的离散表达式进行求解。

向量形式的全水动力学方程的离散表达式为：

$U_i^{n_t+1}=U_i^{n_t}+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(F_{i+1/2}^{n_t+1}-F_{i-1/2}^{n_t+1})+{\mathrm{\ΔtS}}_i^{n_t}$

具体地，本发明实施例提供的地面灌溉中地表水流运动的模拟方法，由于其向量形式的全水动力学方程是对时间和空间(即长度)进行的偏导，所以需要对向量形式的全水动力学方程进行时空离散化处理。进行时空离散化处理时，t为目标畦田沿水流方向上的水流推进时间，t被离散成0、1、2、3、.....、n_t、n_t+1等有限个时间节点，时间节点之间的时间差就是针对t得到的时间步长Δt，Δt可以选择1-10s中的任意数值，例如选择1s、2s、5s、10s等；x为目标畦田水流方向上的长度，x也被离散成0、1、2、3、.....、i、i+1等有限个空间节点，i+1/2即表示空间离散节点i与i+1之间的中点；i-1/2为空间离散节点i与i-1之间的中点，空间节点之间的距离就是针对x得到的空间步长Δx，Δx可以选择1-10m中的任意数值，例如选择1m、2m、5m、10m等。Δt和Δx值是对时间和空间进行不同程度的离散，Δt和Δx的取值越小，说明离散程度越高，计算结果的精确度也越高。

具体地，在求解上述全水动力学方程的时空离散表达式时，首先求解其中的的计算公式中涉及到目标畦田地表水入渗率i_c，i_c是通过Kostiakov入渗公式计算得到，Kostiakov入渗公式为：

$i_c=k_{in}{\mathrm{\α\τ}}_{in}^{\mathrm{\α}-1}$

其中，k_in为土壤入渗参数，单位为m/s^α；α为无量纲入渗指数参数；τ_in为入渗受水时间，单位为s。具体地，通过双环渗水试验，可以测得目标畦田的土壤入渗参数以及无量纲入渗指数参数。双环入渗水试验是一种常用的土壤入渗参数测量方法，其所用设备主要包括内外两个同心设置的圆环以及马利奥特瓶。将内环和外环插入土层一定深度，利用马利奥特瓶向内环注水，连续观测并记录马利奥特瓶内水位随时间的变化，即可获得目标畦田的土壤入渗参数以及无量纲入渗指数参数。入渗受水时间即试验土壤从有水渗入直至土壤表面没有水的时间，且目标畦田中每一个测点处的入渗受水时间是不同的。通过上述三个参数，可以计算得到目标畦田的地表水入渗率。

本领域技术人员可以理解的是，在一块目标畦田内，通常会设置多个土壤入渗参数测点，每一个测点最终得到的土壤入渗参数k_in以及无量纲入渗指数参数α可能会有所不同。当多个测点得到的土壤入渗参数k_in以及无量纲入渗指数参数α不相同时，在计算该目标畦田的地表水入渗率时所用的土壤入渗参数k_in以及无量纲入渗指数参数α的取值应当为各测点得到的土壤入渗参数k_in以及无量纲入渗指数参数α的平均值。

具体地，求解全水动力学方程的时空离散表达式的关键步骤是求解其中的迎风格式耗散项，即和将和用q和h表示出来，即用第二表达式和第三表达式来表示和现有技术中，一般通过经典Roe模拟方法对和进行求解，经典Roe模拟方法一般适用于洪水等湍急的水流模式，求解过程十分复杂，计算量庞大。而本发明实施例提供的地面灌溉中地表水流运动的模拟方法，对经典Roe模拟方法中求得的一系列F值进行了算数平均，使和的计算过程得到了简化。这种简化使是建立在地面灌溉模式下的，且模拟出的数值曲线与经典Roe模拟方法十分接近，说明该模拟方法准确度较高，其完全适用于农田地面灌溉中地表水流运动的模拟。采用本发明实施例提供的地面灌溉中地表水流运动的模拟方法，缩短了模拟计算的时间，提高了模拟计算的效率，使该模拟方法在进行大规模地面灌溉水流运动的模拟时，具有更突出的优势。

具体地，通过将第一表达式、第二表达式和第三表达式代入全水动力学方程的时空离散表达式，得到以q和h为未知数的地表水流运动的模拟方程，该模拟方程为：

$\begin{array}{c}{h}_i^{n_t+1}=h_i^{n_t}+0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\[q_{i+1}^{n_t+1}+q_i^{n_t+1}-0.25\left(u_{i+1}^{n_t+1}+u_i^{n_t+1}-\sqrt{{\mathrm{gh}}_{i+1}^{n_t+1}}-\sqrt{{\mathrm{gh}}_i^{n_t+1}}\right)\left(q_{i+1}^{n_t+1}-q_i^{n_t+1}\right)\]\\ -0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\[q_i^{n_t+1}-q_{i-1}^{n_t+1}-0.25\left(u_i^{n_t+1}+u_{i-1}^{n_t+1}-\sqrt{{\mathrm{gh}}_i^{n_t+1}}-\sqrt{{\mathrm{gh}}_{i-1}^{n_t+1}}\right)\left(q_i^{n_t+1}-q_{i-1}^{n_t+1}\right)\]\\ +\mathrm{\Δ}t\×k_{in}\mathrm{\α}{\left({\mathrm{\τ}}_{in}^{1-\mathrm{\α}}\right)}_i^{n_t}\end{array}$

$\begin{array}{c}{q}_i^{n_t+1}=q_i^{n_t}+0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\{{\left(qu\right)}_{i+1}^{n_t+1}+{\left(qu\right)}_i^{n_t+1}-0.25\left(u_{i+1}^{n_t+1}+u_i^{n_t+1}+\sqrt{{\mathrm{gh}}_{i+1}^{n_t+1}}+\sqrt{{\mathrm{gh}}_i^{n_t+1}}\right)\[{\left(qu\right)}_{i+1}^{n_t+1}-{\left(qu\right)}_i^{n_t+1}\]\}\\ -0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\[q_i^{n_t+1}-q_{i-1}^{n_t+1}-0.25\left(u_i^{n_t+1}+u_{i-1}^{n_t+1}+\sqrt{{\mathrm{gh}}_i^{n_t+1}}+\sqrt{{\mathrm{gh}}_{i-1}^{n_t+1}}\right)\[{\left(qu\right)}_i^{n_t+1}-{\left(qu\right)}_{i-1}^{n_t+1}\]\}\\ -0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}{h}_i^{n_t}\left(b_{i+1}-b_{i-1}\right)-\mathrm{\Δ}tg\frac{n^2{q}_i^{n_t}\left|q_i^{n_t}\right|}{{\left(h^{7/3}\right)}_i^{n_t}}\end{array}$

该模拟方程中b、n、k_in、α、τ_in、Δx、Δt为已知量，n_t和i根据迭代的步数具体进行取值，将q的初始条件设置为畦田内所有空间点均为无单宽流量值(即零值)，边界条件为畦田首部为实测获得的入畦单宽流量，该值采用超声波流量计进行测量。

由于u＝q/h，故通过联立求解上述两式，即可计算出每一测点处的q和h，并得到h关于q变化的曲线，以此判断地面灌溉中地表水流沿畦长方向各点的入渗情况，从而选择合适的灌溉水量，节约水资源。

具体地，由于本发明实施例提供的模拟方法适用于农田的地面灌溉，所以沿目标畦田畦长方向的单宽流量q的适用范围为0.002-0.01m³/(s·m)，对具有适宜q值的地面灌溉过程进行模拟，可以获得更加精准的模拟结果。

以下将通过具体实施例进行详细阐述：

实施例1

本实施例提供了一种地面灌溉中地表水流运动的模拟方法，通过该模拟方法对河北冶河灌区畦田中地面灌溉时地表水流的运动情况进行模拟，选取的河北冶河灌区畦田的基本情况为：畦长100m，通过超声波流量计测得水流首次进入畦田时的单宽流量为3L/s·m，以作为求解q的初始条件，通过秒表测得畦田的受水时间为100min，通过双环渗水试验测得畦田的土壤入渗参数k_in为0.00063m/s^α，无量纲入渗指数参数α为0.23，通过畦首实测的地表水深h估算得到畦面糙率系数n为0.08s/m^1/3。

具体模拟步骤如下：

步骤1、选择向量形式的全水动力学方程作为模拟目标畦田地面灌溉中地表水流运动的地面灌溉模型，其中，向量形式的全水动力学方程为：

$\frac{\∂U}{\∂t}+\frac{\∂F}{\∂x}=S$

$U=\left(\begin{array}{c}h\\ q\end{array}\right)$

$F=\left(\begin{array}{c}q\\ qu+\frac{1}{2}{\mathrm{gh}}^2\end{array}\right)$

$S=\left(\begin{array}{c}{i}_c\\ -gh\frac{\∂b}{\∂x}-gh\frac{n^{2}u\left|u\right|}{h^{4/3}}\end{array}\right)$

其中，F为物理通量，U为因变量向量，S为源项向量，t为目标畦田沿水流方向上的水流推进时间，x为目标畦田水流方向上的长度；u为水流沿目标畦田畦长方向的垂向均布流速，单位为m/s；h为地表水深h，单位为m；q为水流沿目标畦田畦长方向的单宽流量，单位为m³/(s·m)，且q＝hu；g为重力加速度，选择9.8m/s²；i_c为实测得到的目标畦田地表水入渗率；b为实测得到的目标畦田地表高程b，单位为m；n为已知的目标畦田的畦面糙率系数，即0.08s/m^1/3。

其中，i_c通过Kostiakov入渗公式计算得到，在本实施例中i_c为：

$i_c=k_{in}{\mathrm{\α\τ}}_{in}^{\mathrm{\α}-1}=0.00063\×0.23\×100{{\mathrm{\τ}}_{in}}^{-0.77}$

步骤2、对向量形式的全水动力学方程进行时空离散化处理，得到全水动力学方程的时空离散表达式：

$U_i^{n_t+1}=U_i^{n_t}+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(F_{i+1/2}^{n_t+1}-F_{i-1/2}^{n_t+1})+{\mathrm{\ΔtS}}_i^{n_t}$

Δt为目标畦田受水时间变量所对应的时间步长，本实施例中Δt为1s；Δx为目标畦田的畦长变量所对应的空间步长，本实施例中Δx为1m；n_t为时间离散迭代步；i为空间离散节点；i+1/2为空间离散节点i与i+1之间的中点；i-1/2为空间离散节点i与i-1之间的中点。

步骤3、提供利用q和h来表达的第一表达式，具体地，第一表达式的为：

$S_i^{n_t}=\left(\begin{array}{c}{\left(i_c\right)}_i^{n_t}\\ -{\mathrm{gh}}_i^{n_t}\frac{b_{i+1}-b_{i-1}}{2\mathrm{\Δ}x}-g\frac{n^2{u}_i^{n_t}\left|u_i^{n_t}\right|}{{\left(h^{1/3}\right)}_i^{n_t}}\end{array}\right)=\left(\begin{array}{c}{\left(i_c\right)}_i^{n_t}\\ -{\mathrm{gh}}_i^{n_t}\frac{b_{i+1}-b_{i-1}}{2\mathrm{\Δ}x}-g\frac{n^2{q}_i^{n_t}\left|q_i^{n_t}\right|}{{\left(h^{7/3}\right)}_i^{n_t}}\end{array}\right).$

步骤4、分别提供利用q和h来表达和的第二表达式和第三表达式，其中，第二表达式为：

$F_{i+1/2}^{n_t+1}=0.5\[F_{i+1}^{n_t}+F_i^{n_t}-0.5{\mathrm{\Λ}}_{i+1/2}^{n_t}(U_{i+1}^{n_t+1}-U_i^{n_t})\]$

其中，

第三表达式为：

$F_{i-1/2}^{n_t+1}=0.5\[F_i^{n_t}+F_{i-1}^{n_t}-0.5{\mathrm{\Λ}}_{i-1/2}^{n_t}(U_i^{n_t+1}-U_{i-1}^{n_t})\]$

${\mathrm{\Λ}}_{i-1/2}^{n_t}=\left(\begin{array}{cc}0.5\left(u_{i-1}^{n_t}+u_i^{n_t}-\sqrt{{\mathrm{gh}}_{i-1}^{n_t}}-\sqrt{{\mathrm{gh}}_i^{n_t}}\right)& 0\\ 0& 0.5\left(u_{i-1}^{n_t}+u_i^{n_t}+\sqrt{{\mathrm{gh}}_{i-1}^{n_t}}+\sqrt{{\mathrm{gh}}_i^{n_t}}\right)\end{array}\right).$

步骤5、将第一表达式、第二表达式和第三表达式代入全水动力学方程的时空离散表达式，得到以q和h为未知数的地表水流运动的模拟方程，

$\begin{array}{c}{h}_i^{n_t+1}=h_i^{n_t}+0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\[q_{i+1}^{n_t+1}+q_i^{n_t+1}-0.25\left(u_{i+1}^{n_t+1}+u_i^{n_t+1}-\sqrt{{\mathrm{gh}}_{i+1}^{n_t+1}}-\sqrt{{\mathrm{gh}}_i^{n_t+1}}\right)\left(q_{i+1}^{n_t+1}-q_i^{n_t+1}\right)\]\\ -0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\[q_i^{n_t+1}-q_{i-1}^{n_t+1}-0.25\left(u_i^{n_t+1}+u_{i-1}^{n_t+1}-\sqrt{{\mathrm{gh}}_i^{n_t+1}}-\sqrt{{\mathrm{gh}}_{i-1}^{n_t+1}}\right)\[q_i^{n_t+1}-q_{i-1}^{n_t+1}\]\\ +\mathrm{\Δ}t\×k_{in}\mathrm{\α}{\left({\mathrm{\τ}}_{in}^{1-\mathrm{\α}}\right)}_i^{n_t}\end{array}$

$\begin{array}{c}{q}_i^{n_t+1}=q_i^{n_t}+0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\left\{{\left(qu\right)}_{i+1}^{n_t+1}+{\left(qu\right)}_i^{n_t+1}-0.25\left(u_{i+1}^{n_t+1}+u_i^{n_t+1}+\sqrt{{\mathrm{gh}}_{i+1}^{n_t+1}}+\sqrt{{\mathrm{gh}}_i^{n_t+1}}\right)\[{\left(qu\right)}_{i+1}^{n_t+1}-{\left(qu\right)}_i^{n_t+1}\]\right\}\\ -0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}\left\{q_i^{n_t+1}-q_{i-1}^{n_t+1}-0.25\left(u_i^{n_t+1}+u_{i-1}^{n_t+1}+\sqrt{{\mathrm{gh}}_i^{n_t+1}}+\sqrt{{\mathrm{gh}}_{i-1}^{n_t+1}}\right)\[{\left(qu\right)}_i^{n_t+1}-{\left(qu\right)}_{i-1}^{n_t+1}\]\right\}\\ -0.5\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}{h}_i^{n_t}\left(b_{i+1}-b_{i-1}\right)-\mathrm{\Δ}tg\frac{n^2{q}_i^{n_t}\left|q_i^{n_t}\right|}{{\left(h^{7/3}\right)}_i^{n_t}}\end{array}$

其中，i＝0，1，2，3，4，.....，I(I为空间节点总数)，n_t＝0，1，2，3，4，.....，N(N为时间节点总数)。

求解时，初始条件q＝0，边界条件为q₀＝0.003。

通过求解目标畦田中每一测点处的q和h，从而获取目标畦田中所有测点处的q和h，以对地面灌溉中地表水流运动进行模拟。

同时，本实施例还采用经典Roe模拟方法对该同一个目标畦田进行模拟，并将本实施例提供的模拟方法与经典Roe模拟方法的模拟数值、模拟精度以及模拟计算所需时间进行对比，以对本实施例提供的模拟方法的模拟效果进行检验。

附图2和附图3分别为河北冶河灌区的畦田灌溉40min和80min时，采用本实施例提供的地面灌溉中地表水流运动的模拟方法与采用经典Roe模拟方法得到的地表水深h的数值模拟对比图，并通过畦田地表水入渗率与受水时间计算入渗水深d，绘制入渗水深d曲线。而单宽流量q与地表水深h相对应，当地表水深h的模拟值较为准确时，则单宽流量q的值也是较为准确的。因此，本领域技术人员在评价地面灌溉中地表水流运动的模拟方法的可行性时，只需参考地表水深h与入渗水深d之间的管线图即可，即如附图2、3中所示。从附图2和附图3中可以看出，通过上述两种方法模拟出的数值，无论是地表水深h还是入渗水深d都十分接近，表示这两者的曲线基本重合在一起。由于经典Roe模拟方法的已经得到广泛认可和应用，说明其模拟真实度较高，模拟效果较好，而通过本实施例提供的模拟方法得到的模拟数值与经典Roe模拟方法得到的数值十分相近，说明本实施例提供的模拟方法也具有较好的模拟效果，可以用来模拟地面灌溉中地表水流的运动情况。

具体地，采用本实施例提供的地面灌溉中地表水流运动的模拟方法测得的地表水深h和入渗水深d与采用经典Roe模拟方法得到的数值相比，在灌水时间为40min时，两种方法所得到的地表水深h之间的平均相对误差为0.12％，入渗水深d之间的平均相对误差为0.11％。在灌水时间为80min时，两种方法所得到的地表水深h之间的平均相对误差为0.09％，入渗水深d之间的平均相对误差为0.11％。可见，利用本实施例提供的方法与经典Roe模拟方法所模拟得到的数据之间的误差很小，模拟效果良好。

然而，实验证明，利用本实施例提供的模拟方法进行上述各数值的模拟计算的时间为2.92min，用经典Roe模拟方法进行模拟计算的时间则为9.42min，可见，采用本实施例提供的模拟方法进行地面灌溉中地表水流运动的模拟计算，可使计算时间缩短69％，从而提高了模拟计算的效率，更适用于大规模畦田的模拟计算。

实施例2

本实施例提供了一种地面灌溉中地表水流运动的模拟方法，该方法的操作过程与实施例1相同，所不同的是：

进行试验的畦田为新疆建设兵团222团畦田，该畦田的基本情况为：畦长150m，通过超声波流量计测得水流首次进入畦田时的单宽流量为6L/s·m(即q的初始条件q₀＝0.006)，通过秒表测得畦田的受水时间为160min，通过双环渗水试验测得畦田的土壤入渗参数k_in为0.00016m/s^α，无量纲入渗指数参数α为0.42，通过畦首实测的地表水深h估算得到畦面糙率系数n为0.11s/m^1/3。

同时，采用经典Roe模拟方法对该同一畦田进行模拟，并将本实施例提供的模拟方法与经典Roe模拟方法的模拟数值、模拟精度以及模拟计算所需时间进行对比，以对本实施例提供的模拟方法的模拟效果进行检验。

附图4和附图5分别为新疆建设兵团222团的畦田灌溉40min和80min时，采用本实施例提供的地面灌溉中地表水流运动的模拟方法与采用经典Roe模拟方法得到的地表水深h的数值模拟对比图，并通过畦田地表水入渗率与受水时间计算入渗水深d，绘制入渗水深d曲线。从附图4和附图5中可以看出，通过上述两种方法模拟出的数值，无论是地表水深h还是入渗水深d都十分接近。由于经典Roe模拟方法的已经得到广泛认可和应用，说明其模拟真实度较高，模拟效果较好，而通过本实施例提供的模拟方法得到的模拟数值与经典Roe模拟方法得到的数值十分相近，说明本实施例提供的模拟方法也具有较好的模拟效果，可以用来模拟地面灌溉中地表水流的运动情况。

进一步地，采用本实施例提供的地面灌溉中地表水流运动的模拟方法测得的地表水深h和入渗水深d与采用经典Roe模拟方法得到的数值相比，在灌水时间为40min时，地表水深h之间的平均相对误差为0.13％，入渗水深d之间的平均相对误差为0.09％，在灌水时间为80min时，地表水深h之间的平均相对误差为0.08％，入渗水深d之间的平均相对误差为0.08％。可见，利用本实施例提供的方法与经典Roe模拟方法所模拟得到的数据之间的误差很小，模拟效果良好。

然而，试验证明，采用本实施例提供的模拟方法进行模拟计算的时间为4.25min，用经典Roe模拟方法进行模拟计算的时间则为13.67min，可见，采用本实施例提供的模拟方法进行地面灌溉中地表水流运动的模拟计算，可使计算时间缩短69％，从而提高了模拟计算的效率，更适用于大规模畦田的模拟计算。

实施例3

本实施例提供了一种地面灌溉中地表水流运动的模拟方法，该方法的操作过程与实施例1相同，所不同的是：

进行试验的畦田为山东簸箕李灌区畦田，该畦田的基本情况为：畦长130m，通过超声波流量计测得水流首次进入畦田时的单宽流量为4L/s·m(即q的初始条件q₀＝0.004)，通过秒表测得畦田的受水时间为120min，通过双环渗水试验测得畦田的土壤入渗参数k_in为0.00026m/s^α，无量纲入渗指数参数α为0.38，通过畦首实测的地表水深h估算得到畦面糙率系数n为0.09s/m^1/3。

同时，采用经典Roe模拟方法对河北冶河灌区的同一畦田进行模拟，并将本实施例提供的模拟方法与经典Roe模拟方法的模拟数值、模拟精度以及模拟计算所需时间进行对比，以对本实施例提供的模拟方法的模拟效果进行检验。

附图6和附图7分别为河北冶河灌区的畦田灌溉40min和80min时，采用本实施例提供的地面灌溉中地表水流运动的模拟方法与采用经典Roe模拟方法得到的地表水深h的数值模拟对比图，并通过畦田地表水入渗率与受水时间计算入渗水深d，绘制入渗水深d曲线。从附图6和附图7中可以看出，通过上述两种方法模拟出的数值，无论是地表水深h还是入渗水深d都十分接近。由于经典Roe模拟方法的已经得到广泛认可和应用，说明其模拟真实度较高，模拟效果较好，而通过本实施例提供的模拟方法得到的模拟数值与经典Roe模拟方法得到的数值十分相近，说明本实施例提供的模拟方法也具有较好的模拟效果，可以用来模拟地面灌溉中地表水流的运动情况，且地面灌溉时间越长，模拟精度越高。

进一步地，采用本实施例提供的地面灌溉中地表水流运动的模拟方法测得的地表水深h和入渗水深d与采用经典Roe模拟方法得到的数值相比，在灌水时间为40min时，地表水深h之间的平均相对误差为0.12％，入渗水深d之间的平均相对误差为0.11％，在灌水时间为80min时，地表水深h之间的平均相对误差为0.10％，入渗水深d之间的平均相对误差为0.11％。可见，利用本实施例提供的方法与经典Roe模拟方法所模拟得到的数据之间的误差很小，模拟效果良好。

然而，试验证明，采用本实施例提供的模拟方法进行模拟计算的时间为3.61min，用经典Roe模拟方法进行模拟计算的时间则为11.58min，可见，采用本实施例提供的模拟方法进行地面灌溉中地表水流运动的模拟计算，可使计算时间缩短69％，从而提高了模拟计算的效率，更适用于大规模畦田的模拟计算。

以上所述仅为本发明的较佳实施例，并不用以限制本发明的保护范围，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种地面灌溉中地表水流运动的模拟方法，其特征在于，所述模拟方法包括以下步骤：

步骤a、选择向量形式的全水动力学方程作为模拟目标畦田地面灌溉中地表水流运动的地面灌溉模型；

所述向量形式的全水动力学方程为：

$\frac{\∂U}{\∂t}+\frac{\∂F}{\∂x}=S$

$U=\left(\begin{array}{c}h\\ q\end{array}\right)$

$F=\left(\begin{array}{c}q\\ qu+\frac{1}{2}{\mathrm{gh}}^2\end{array}\right)$

$S=\left(\begin{array}{c}{i}_c\\ -gh\frac{\∂b}{\∂x}-gh\frac{n^{2}u\left|u\right|}{h^{4/3}}\end{array}\right)$

其中，F为物理通量，U为因变量向量，S为源项向量，t为所述目标畦田以灌溉起始时刻作为原点的时间坐标，单位为s，x为所述目标畦田中任意测点与所述目标畦田的畦首之间的距离，单位为m；u为水流沿所述目标畦田畦长方向的垂向均布流速，单位为m/s；h为地表水深，单位为m；q为水流沿所述目标畦田畦长方向的单宽流量，单位为m³/(s·m)，且q＝hu；g为重力加速度，单位为m/s²；i_c为实测得到的所述目标畦田地表水入渗率；b为实测得到的所述目标畦田地表相对高程，单位为m；n为已知的所述目标畦田的畦面糙率系数，单位为s/m^1/3；

步骤b、对所述向量形式的全水动力学方程进行时空离散化处理，得到所述全水动力学方程的时空离散表达式；

所述全水动力学方程的时空离散表达式为：

$U_i^{n_t+1}=U_i^{n_t}+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(F_{i+1/2}^{n_t+1}-F_{i-1/2}^{n_t+1})+{\mathrm{\ΔtS}}_i^{n_t};$

其中，Δt为针对所述t得到的时间步长；Δx为针对所述x得到的空间步长；n_t为针对所述t得到的时间离散节点；i为针对所述x得到的空间离散节点；i+1/2为空间离散节点i与i+1之间的中点；i-1/2为空间离散节点i与i-1之间的中点；

步骤c、提供利用所述q和所述h来表达所述的第一表达式；

所述第一表达式的为：

$S_i^{n_t}=\left(\begin{array}{c}{\left(i_c\right)}_i^{n_t}\\ -{\mathrm{gh}}_i^{n_t}\frac{b_{i+1}-b_{i-1}}{2\mathrm{\Δ}x}-g\frac{n^2{u}_i^{n_t}\left|u_i^{n_t}\right|}{{\left(h^{1/3}\right)}_i^{n_t}}\end{array}\right);$

步骤d、分别提供利用所述q和所述h来表达所述和的第二表达式和第三表达式；

所述第二表达式为：

$F_{i+1/2}^{n_t+1}=0.5\[F_{i+1}^{n_t}+F_i^{n_t}-0.5{\mathrm{\Λ}}_{i+1/2}^{n_t}(U_{i+1}^{n_t+1}-U_i^{n_t})\]$

其中，

所述第三表达式为：

$F_{i-1/2}^{n_t+1}=0.5\[F_i^{n_t}+F_{i-1}^{n_t}-0.5{\mathrm{\Λ}}_{i-1/2}^{n_t}(U_i^{n_t+1}-U_{i-1}^{n_t})\]$

${\mathrm{\Λ}}_{i-1/2}^{n_t}=\left(\begin{array}{cc}0.5(u_{i-1}^{n_t}+u_i^{n_t}-\sqrt{{\mathrm{gh}}_{i-1}^{n_t}}-\sqrt{{\mathrm{gh}}_i^{n_t}}& 0\\ 0& 0.5(u_{i-1}^{n_t}+u_i^{n_t}+\sqrt{{\mathrm{gh}}_{i-1}^{n_t}}+\sqrt{{\mathrm{gh}}_i^{n_t}}\end{array}\right);$

步骤e、将所述第一表达式、所述第二表达式和所述第三表达式代入所述全水动力学方程的时空离散表达式，得到以所述q和所述h作为未知数的地表水流运动的模拟方程，通过对所述地表水流运动的模拟方程进行求解，获取所述目标畦田中每一测点处的所述q和所述h，从而对地面灌溉中地表水流运动进行模拟。

2.根据权利要求1所述的模拟方法，其特征在于，所述目标畦田地表水入渗率i_c通过Kostiakov入渗公式计算得到，所述Kostiakov入渗公式为：

$i_c=k_{in}{\mathrm{\α\τ}}_{in}^{\mathrm{\α}-1}$

其中，k_in为土壤入渗参数，单位为m/s^α；α为无量纲入渗指数参数；τ_in为入渗受水时间，单位为s。

3.根据权利要求2所述的模拟方法，其特征在于，所述畦面糙率系数n通过在所述目标畦田的畦首实测畦首地表水深值，并根据所述畦首地表水深值估算得到。

4.根据权利要求1所述的模拟方法，其特征在于，所述时间步长为1-10s。

5.根据权利要求1所述的模拟方法，其特征在于，所述空间步长为1-10m。

6.根据权利要求1所述的模拟方法，其特征在于，在所述步骤e中，在求解所述地表水流运动的模拟方程时，将水流首次进入畦田时的单宽流量作为求解所述q的初始条件。