灌溉渠系非恒定流数值模拟及闸门运行设计

古 玉，方天钰，陈 建，刘晓东，韩 宇

(1.中国农业大学水利与土木工程学院，北京 100083；2.中国农业大学工学院，北京 100083)

0 引 言

目前我国水资源短缺问题十分严重，已成为制约我国国民经济和社会发展的重要因素之一。我国现有灌溉面积0.664 亿hm2，农业用水更是占到用水总量的62.4%[1]，因此为应对水资源危机，采用农业用水节水措施是必要的。灌溉渠系进行输配水时，需要根据灌区用水量要求通过闸门进行调节，此时渠道中的水流受到干扰必定会出现非恒定流，使渠道中的水位、流量既沿程发生变化，又随时间发生变化。由于无法预测该过渡过程，常发生重复调度，导致渠道配水的可靠性及准确性很差，甚至跑水和弃水，造成水量浪费，降低灌溉质量。为使整个渠系引水、输水、配水过程科学合理，减少渠道输水损失，提高水的利用效率，达到节水的目的，数值模拟是最为有效可行的方法。

明渠非恒定流的理论研究起源于18世纪后期Laplace和Lagrange对浅水波的研究，之后Lagrange给出了其波速公式。到1871年，Saint-Venant在做了大量研究的基础上，提出了明渠非恒定流的理论和通用方程，即Saint-Venant方程，使明渠非恒定流现象得到了数学描述，推动了通过数学方法研究这一问题的发展。之后国内外学者不断深入研究，提出了有限差分法、特征线法和有限单元法等多种求解方法，其中有限差分法又包括直接解法、分级解法、单元划分法、混合解法以及松弛迭代法等。相关理论的不断发展完善的基础上，随着计算机的应用，明渠非恒定流的数值模拟也得到了充分发展。法国的CARIMA[2]模型以及荷兰的MODIS[3]、DUFLOW[4]模型都在田间试验中得到了较好验证。Gary P Merkley[5]教授开发出了一个可视化的模拟软件CanalMan，它允许用户自己建立所需渠道并定义相关参数，进一步修改后形成RootCanal水力模拟软件，但在后来的使用中发现此软件也存在着一些问题。到现在为止，应用较为广泛且效果较好的一个软件是丹麦水力研究所(DHI)开发的Mike11。Mike11虽然模拟效果较好，但其无法模拟闸门开度调控过程，且费用很高。在国内，20世纪80年代起，我国的非恒定流数值模拟开始有了进展，针对灌溉渠网，从1982年到2000年，王长德[6-9]建立了上、下游常水位自动控制渠道以闸门运动为动态边界条件的渠道非恒定流计算模型，并将P+PR算法(比例+比例微分算法)与比威尔算法相结合。2001年，赵竞成等[10]建立了较完整的渠系运行模型，编制了具有一定通用性和可扩充性的计算机模拟软件。到2010年，韩宇等[11]模拟了上游和下游常水深两种运行方式下梯形断面渠段由于不同闸门调控方式所引起的非恒定流过渡过程并进行了比较，相对明确地模拟了闸门调节过程。2013年，章少辉等[12]把渠系输配水控制方程变形处理，在具备高精度和分辨率的标量有限体积法的基础上建立了末级渠系非恒定流数值模拟模型并验证。

发展至今，非恒定流模拟的相关理论已经较为完善，Preissmann四点隐格式和特征线最为常用。在数值模拟方面，河网非恒定流数值模拟技术已经较好，Mike软件得到普遍使用且模拟效果较好，但灌溉渠网中非恒定流数值模拟仍缺乏较好的相关模型或程序。且国内外相较来说，我国针对非恒定流的研究开始的较晚，虽然理论方面发展较完善，但多基于南水北调工程，忽略某一段面的改变对其上下游稳态水流的干扰，并且由于灌溉渠道中闸门的调控方式及其形成的边界条件不容易合理概化，针对河网的数值模拟模型无法直接用于灌溉渠道，所以在灌溉渠网非恒定流数值模拟技术方面仍有所欠缺。另外，对于数值模拟可视化的研究较少，且多针对河网，2000年，徐小明等[13]用Fortran语言，对计算出的水位变化过程作了可视化显示研究；2002年，马洪明等[14]，运用计算机信息平台，建立了一套图形显示水力模型的支持系统；2008年，陈栋[15]基于松弛迭代法，运用混合编程技术，开发了河网非恒定流数值模拟的可视化系统；但这些成果都不够成熟故没有得到推广使用。随着计算机技术的不断发展成熟，将计算机技术与非恒定流模拟相结合，建立起一个较为简便且适用性强的灌溉渠系非恒定流模拟软件，对渠道各处的水力要素进行模拟，将能够使工作人员及时了解渠道水流的动态变化过程，给出适当的闸门调控方法，实现对水资源的科学合理调配，从而提高渠系输配水的利用率，大大减少农业用水的浪费。

本文采用特征线法求解Saint-Venant方程，并建立简单渠系模型模拟引汤灌区北干渠及九、十支渠系闸后水流过渡过程，并与Mike11模拟结果进行比对并验证该模型准确性。其后基于该模拟结果利用闸下出流公式建立闸门开度模拟模型，同时对多级闸门同步与逐级调控两种情况进行模拟，讨论其优缺点，对渠系闸门调控进行了简单设计。

1 灌溉渠系非恒定流数学模型

1.1 控制方程

非恒定流是指过水断面上的水力要素(流量、水位及流速等)不仅随空间位置的改变而变化，而且随时间不断变化的水流。灌溉渠道中的非恒定流为明渠非恒定流，圣维南方程就是针对这一问题的函数关系式，即明渠非恒定流的控制方程，它由连续性方程和运动方程组成。

从有旁侧入流的明渠非恒定流中取一长为ds的微小流段，由质量守恒原理可推出连续性方程如下：

(1)

假定水面波动是渐变的，忽略局部水头损失，由动量定理或牛顿第二定律可推导出其运动方程如下：

(2)

式中：B为水面宽度，m；Z为水位，m；Q为流量，m3/s；s为流程，m；ql为单位长度上的侧向汇流流量[流入为正，流出为负，m3/(s·m)]；A为过水断面面积，m2；g为重力加速度，m/s2；C为谢才系数；R为水力半径；h为水深，m。

1.2 求解方法

本文采用特征线法，用矩形网格划分求解域从而把微分方程的连续问题离散化，以有限个网格点代替原连续区域，用数值积分逼近微分方程，建立网格函数的代数方程组，在离散点上求解网格函数的值，总的来说就是把圣维南偏微分方程组化为常微分方程即特征方程组进行求解[16, 17]：

(3)

(4)

采用矩形网格特征差分法进行离散：首先将求解域依据时间t和距离s进行矩形网格划分，如图1，网格的交点称为结点。其中，i=0,1,2,…,n表示距离s的段号；j=0,1,2,…,n表示时间t的段号。Δsi=si-si-1表示距离步长，即渠段分段长度，可以等步长也可变步长；Δtj=tj-tj-1表示时间步长，需要根据稳定格式的限制和实际要求确定。图1中，顺、逆特征线(ω±)在时层j+1上的交点为点P，与时层j的交点分别为L和R，其间的结点为M，左右两边点各为点D和点E，求解方程即求解交点P。当给定初始条件后，则根据时间步长进行循环计算求解。

图1 矩形网格图示意图Fig.1 The Schematic of Rectangle grid diagram

要求点P的水位和流量，则需先求出点L和R，这两点不在网格结点上，因此采用柯朗格式进行线性内插，之后利用积分和积分中值定理可得到离散结果如下：

(5)

(6)

(7)

同时，由于柯朗格式为一阶显格式，为使计算稳定，需要选择合理的距离步长和时间步长使其满足稳定性条件式(8)：

(8)

1.3 边界条件

灌溉渠道边界条件包括内、外两种边界条件，其中内边界条件是指闸门等控制建筑物在渠段中的调控所产生的边界处的变化，外边界条件是指入流及出流边界条件。

上述的计算格式只适用于内部点的计算，如图2所示，对于边界点来说：仅有一条特征线，即只有两个方程，但其未知数有3个，故需补充边界条件。

上、下游边界条件一般情况下有如下3种情况：①Z=Z(t)，即水位过程已知；②Q=Q(t)，即流量过程已知；③Q=f(Z)，即水位流量关系已知。

图2 上下游边界点示意图Fig.2 The schematic diagram of the upstream and downstream boundary points

1.4 汊点处条件

灌溉渠系多呈树状，因此在其布置形式中必然存在众多分岔。灌溉渠道分为干、支、斗、农四级，它们之间必然由一系列分支、汇交点连接，这些连接点即为汊点。在汊点处，水流需满足相应的连接条件，即流量相容和能量相容条件[18, 19]。

图3 分岔渠道示意图Fig.3 The Schematic diagram of the distribution channel

流量相容条件为：

(9)

式中：n为汊点处的渠道数目；Qi为第i个渠道在汊点处的流量；当渠道水流流入汊点时sli取“+1”，流出时取“-1”。对图3所示的分岔渠道，其流量相容条件可写为：

Q1=Q2+Q3

(10)

能量相容条件为：

首先将汊点概化为一个几何点，如图3，断面A-A、B-B、C-C围出了一个小流域，由于只须计算其边界断面上的流量或水位，故可不研究此内的非恒定流状态，假设流域的面积为零。在理想情况下，即不考虑局部水头损失，同时一般情况下不存在水位突变，故有：

Z1=Z2=Z3

(11)

式中：Zi(i=1,2,3)为第i个渠道在汊点处的水位。对于灌溉渠道，汊点处经常安有闸门等水工建筑物，故考虑局部水头损失，可列出能量方程如下：

(12)

(13)

式中：αi(i=1,2,3)为能量修正系数且α1=0。

结合上述汊点处条件并参考文献[17],本文对分岔渠系的模拟采用逐级追赶的方法。顾名思义，就是对渠系中的各级渠道逐次求解，以前一级渠道模拟得到结果作为下一级渠道的初始条件，分岔处则根据汊点处的流量相容条件，用干渠所得结果减去支渠的结果作为下一级渠道的初始条件，这样逐级追赶，得出所需结果。这种方法思路清晰，易于理解与掌握，在编程过程中只需注重汊点处的条件变化设置，可以大大减少编程工作量。在实际应用中，若汊点较少且流量变化率不大，则可直接按照此法逐级求解。

2 渠系数据及模拟条件

引汤灌区位于黑龙江省汤原县，灌区的水源来自汤旺河，设计灌溉区域2.68 万hm2，工程总的设计引流量为160 m3/s。本文选取部分引汤北干渠及其九、十支渠作为分岔渠系模拟模型，相关模型数据见表1、表2及表3。

表1 北干渠数据表Tab.1 North trunk channel data table

表2 九支渠数据表Tab.2 Data table for nine canals

表3 十支渠数据表Tab.3 Data table for ten canals

模拟条件：模拟总时长200 min，时间步长取值根据式(8)计算得到：干渠1，Δt≤62 s；干渠2，Δt≤61 s；干渠3，Δt≤62 s；九支渠，Δt≤108 s；十支渠，Δt≤112 s，故各渠道模拟时间步长均取为60 s。各渠道均采用下游常水深运行方式，其中干渠1上游流量在20 min内从25.28 m3/s线性增加到30.11 m3/s，九支渠上游流量在20 min内从1.02 m3/s线性增加到1.32 m3/s，干渠2的上游断面的流量变化则为干渠1模拟得到的下游流量减去九支渠上游断面的流量。十支渠上游流量在20 min内从0.79 m3/s线性增加到1.02 m3/s，干渠3上游断面的流量变化则为干渠2模拟得到的下游流量减去十支渠上游断面的流量。

3 模拟结果分析

在上述工况下，可模拟得到各渠道上游断面水位变化和下游断面流量变化，将所得过渡过程数值结果与用Mike11软件模拟得到的结果进行对比，来验证其可靠性，见图4～图8。

图4 干渠1上下游断面水位流量变化过程Fig.4 Process of water level flow change in the upper and lower sections of the main canal 1

由以上模拟结果可知，本研究中所建立模型与Mike11模型模拟结果大致吻合，说明该模型模拟结果较为准确，并由此得出具体结果如下。

(1)在同时调节干渠一、九支渠以及十支渠上游闸门使其按照模拟条件产生相应流量时，各渠道水位、流量在50 min之后都基本达到了稳定。

(2)干渠一的下游流量在一次波动后达到稳定，在流量线性增大的20 min内，上游水位也表现为线性增长，经过小幅波动后过渡到稳定状态。

(3)在经过之前的传播及九支渠的分流后，干渠二的下游流量波动更加明显，同时干渠二水位变化表现为非线性，在经历一定波动后达到稳定；干渠三的流量水位变化规律与干渠二类似，只是在经过进一步传播后，波动较干渠二更为明显。

(4)在给定模拟条件下，九支渠的流量变化幅度很小，故其下游断面流量波动很小，随后逐渐达到稳定。同时在流量线性增大的20 min内，其上游水位也表现为线性增长，之后趋于稳定；十支渠的流量水位变化与九支渠类似，故不再赘述。

图5 干渠2上下游断面水位流量变化过程Fig.5 Process of water level flow change in the upper and lower sections of the main canal 2

图6 干渠3上下游断面水位流量变化过程Fig.6 Process of water level flow change in the upper and lower sections of the main canal 3

图7 九支渠上下游断面水位流量变化过程Fig.7 Process of water level flow change in the upper and lower sections of the nine branch canal

图8 十支渠上下游断面水位流量变化过程Fig.8 Process of water level flow change in the upper and lower sections of the ten branch canal

4 闸门运行设计

4.1 闸门开度调节设计

实现灌溉渠道自动化是提高渠系输配水效率，建设节水农业的重要内容。而其自动化设计中的一个重要环节就是合理的闸门调控设计。灌溉渠道运行时，往往需要根据灌溉流量需求来进行水量输配，所需水流从水源到田间的过程中受各种因素干扰，引起取水口水深变化，就需要通过相应闸门来进行控制，要得到合理的闸门调控措施，就需要知道相应波动下闸门的开度随时间的变化。根据本文建立的渠道模型可以求出相应于闸门断面的水位或流量变化，从而可以根据闸孔出流的公式添加闸门开度计算模型，得出闸门的调控过程，为相关工作人员提供决策依据。

闸孔出流可分为自由出流和淹没出流，如图9和图10所示，它们的流量计算都有对应的经验公式，可统一为以下格式：

(14)

式中：Qz表示通过闸门的流量；bz表示闸孔宽度；Go表示闸门开度；hu表示闸门上游水深；cd表示流量系数，对于平板闸门，其取值可按下式计算：

图9 淹没出流Fig.9 Flooded outflow

图10 自由出流Fig.10 Free flow

自由出流：

(15)

淹没出流：

(16)

式中：hd表示闸门下游水深。

在上述引汤灌区分岔渠系模型中，渠道采用下游常水深运行方式，其给定的初始条件中的上游流量的线性变化应该是在下游水位保持不变的情况下由上游闸门调控所产生的，因此根据模拟得到的上游断面的水位变化过程(即闸后水深变化)，反推出上游闸门在流量线性变化情况下对应的开度随时间的调节过程。

对于九支渠，首先通过试算得出其初始闸门开度：已知初始闸门上、下游水深分别为1.39 m和0.66 m，闸孔宽度为0.8 m，先假定一个闸门开度，判断其属于自由出流还是淹没出流，根据相应公式计算出其流量系数，之后代入式(14)求出流量并与已知初始流量对比，若不一致则继续试算至一致为止，通过计算得到九支渠初始闸门开度为0.47 m。得到初始闸门开度后，以此为起点，利用上述试算过程针对上游断面每一时刻都进行求解，则可得到九支渠上游闸门开度调节过程，如图11(a)所示。由图可看出：在前5 min闸门开度不变(水深较小，闸门对水流没有控制)之后闸门开度随时间线性增加，到20min后保持不变，其线性调控过程可拟合为以下函数：y=0.007 6x+0.419 8；随后在本研究中又模拟该渠道在流量从1.32减少到1.02时的闸门开度变化过程，如图11(b)所示，从图中可以看到，在前7 min开度保持不变(水深较大，不需调节闸门)，随后线性减小，到20 min时达到稳定，其线性调控过程可拟合为y=-0.007 9x+0.677 1。

图11 流量增加或减少时九支渠上游闸门调节过程Fig.11 Regulating process of the upstream gate of the nine branch canal when the flow increases or decreases

同理，对于十支渠也可求出其闸门调控过程。十支渠初始闸门上、下游水深分别为2.49 m和0.63 m，闸孔宽度为0.47 m，试算得其初始闸门开度为0.44 m。十支渠闸门调控过程如图12(a)所示，在前7 min闸门开度不变，之后闸门开度随时间线性增加，到20 min后保持不变，其线性调控过程可拟合为以下函数：y=0.006 7x+0.381 9；随后又模拟了十支渠流量从1.02减少到0.79时的闸门调控过程，如图12(b)所示，从图12中可以看到，在前7 min闸门开度保持不变，随后线性减小，到20 min时达到稳定，其线性调控过程可拟合为y=-0.007 9x+0.677 1，十支渠总体规律与九支渠类似。

图12 流量增加或减少时十支渠上游闸门调节过程Fig.12 Adjustment process of the upstream gate of the ten branch channel when the flow rate increases or decreases

4.2 多级闸门运行设计

在分岔渠系中，往往存在多级闸门，当要通过闸门调控来满足下游流量需求时，可以同时调节各级闸门，也可以一级一级分别调节，这样同步与逐级的两种调节方式下，渠道中的水流动态也会有所差别。在上述引汤灌区渠道模型中，模拟情况为同时调节干渠一、九支渠和十支渠的上游闸门，在此基础上，再模拟逐级调节情况，即干渠一达到稳定后调节九支渠，干渠二稳定后调节十支渠，由于干渠一和九、十支渠的模拟工况未变化，因此将模拟得到的干渠二与干渠三的流量变化过程进行对比，如图13、图14所示。

图13 干渠2下游流量变化过程对比Fig.13 Comparison of downstream flow changes in the main canal 2

图14 干渠3下游流量变化过程对比Fig.14 Comparison of downstream flow changes in the main canal 3

由以上两图可以看出：当各级闸门同步调节时，渠道中的水流波动不断传播，各级渠道中的波动幅度越来越明显，而各渠段达到稳定的时间虽有所延长，但却变化不大，均在50～60 min之后稳定，最后一级渠道即北干渠三在60 min之后达到稳定；当各级闸门逐级调节时，下一级闸门是在上一级渠道中水流达到稳定后开始调节，故渠道中的水流波动无法逐级传播，对于各级渠道来说，其波动幅度都比较小且稳定时间也比较短，在35～40 min之后达到稳定。但在这种情况下，最后一级渠道稳定到目标量的时间不能仅仅看这一渠道，而是应该从上游第一级闸门的调控开始考虑，于是北干渠三中水流波动达到稳定的时间为自身的40 min加上各上级渠道的稳定时间，则其最终的总稳定时间比同步调控下更长。

目前在我国的大部分灌区还无法实现同步调控，在传统的运行方式中采用逐级调控，即首先调节渠首的闸门，达到稳定后进行下级闸门调节，按此规律依次推进。在这种情况下，渠道中的水流量波动小，对渠道的冲刷也更小，运行也更安全，同时可降低渠道的超高标准，减少工程造价，但是各级闸门的启闭互相影响，不利于简化操作，并且渠系达到最终稳定所需的总时长也会比较长。若采用同步调控，则渠系到达最终稳定的总时长较短，且各级闸门可同时调控，不用以上游渠道的过渡过程为依据，相对独立，启闭操作简单方便，但这种情况下渠道中水流波动会比较大，对渠道产生冲刷，也不利于渠道的安全运行。在实际灌溉中，若时间允许，则推荐逐级调控的方式；但如果需要在尽可能短的时间内达到需求量，则推荐采用同步调控的方式。

5 结 语

本文总结了采用特征线法对圣维南方程进行离散的具体过程、边界条件、汊点处条件以及稳定性条件，建立了渠系非恒定流数值模拟模型以及闸控模型，且所建模型可以通过修改参数来实现不同渠道的模拟。

(1)本文基于特征线法及汊点处条件，建立了渠系模拟模型，并对引汤灌区北干渠及九、十支渠在给定模拟条件下渠道中所产生非恒定流过渡过程进行模拟。将该模拟结果与Mike11所获结果比较后发现较为吻合，证明本研究所建模型具有一定可靠性。

(2)依据该模型模拟结果，建立了相应闸控模型，结果表明在流量线性变化的条件下闸门开度变化也是线性的。在本文中通过计算得到九支、十支渠的闸门开度调节过程，得到拟合的闸控方程，作为相关工作人员决策依据，同时实现其智能控制提供依据。

(3)基于本文所建立的渠系模型，模拟各级闸门在同步调控与逐级调控两种不同方式下的非恒定流过渡过程，结果表明在同步调控时，总调控时间短且此时闸门的启闭操作方便，但渠道中的水流波动较大，不利于渠道安全；逐级调控时，其水流波动则会比较小且针对各部分渠道来说稳定时间较短，但此种情况下闸门启闭操作则需满足上级渠段稳定时长，同时渠系调控的总时长会变长。在实际中，若时间允许则可采用逐级调控，若需在尽量短的时间内得到目标流量则可采用同步调控，但目前我国还有很多灌区无法实现同步调控。

在实际当中，一个灌区的灌溉渠道系统多为由干渠、支渠、斗渠和农渠等多级渠道组成的树状网络，其连接非常复杂，因此，在对简单的分岔渠系模拟的基础上，进一步扩展到对更加复杂的大型灌溉渠网的模拟仍是未来工作的重要方向。同时，由于实测资料的欠缺，本文未能将模拟结果与实测资料进行对比验证。另外，本文所建模型未考虑渠道渗漏损失，对于渠段中有较大底坡跌落的模拟也会失真，这些问题都需要在今后的工作中进行不断修正与改进。