复杂输水系统水力过渡的数值解法比较及适用性分析
穆祥鹏
（1979—），男（回族），天津人，博士，从事长距离输水的水力瞬变及控制研究。
天津大学 建筑工程学院，天津 300072
练继建
天津大学 建筑工程学院，天津 300072
刘翰河
天津大学 建筑工程学院，天津 300072
不同数值解法在求解不同复杂输水系统水力过渡问题时，会体现出各自的特点及不同程度的适用性。因此，针对不同的输水方式和水力过渡类型来选择合适的数值解法，既可以保证计算精度，又可以提高计算速度。针对输水系统的无压流、有压流、无压接有压、明满流交替等不同过渡过程，对特征线法和隐式差分法进行了对比，分析了这两种方法的特点，研究和确定了这两种方法在不同输水方式和水流状态下的适用性。
基金项目:国家杰出青年科学基金（50725929）；国家自然科学基金项目（50538060；50679052）。
输水系统在发生水力过渡过程时，水力参数随时间发生急剧变化，往往超出恒定流的设计范围，影响到输水系统的设计和管理运行，越来越受到人们的重视。非恒定流基本微分方程包括运动方程和连续方程［1］，这是实现输水系统过渡过程模拟的理论基础。此二方程联立，加上相应的边界条件，辅之一定的数值求解方法就可求得水力过渡的变化过程。根据Preissmann窄缝原理［2］，无压流和有压流的过渡过程可用统一的控制方程来求解。
运动方程：
连续方程：
式中：y为对于无压流的水深，而对于有压流相当于输水管道中的相对于管道底高程的压力；c为对于无压流的水面波速，即
，而对于有压流则为水击波速。以上两个控制方程是拟线性双曲型偏微分方程组，无法得出解析解，一般采用差分方法得到近似的数值解，即用差商代替导数，将偏微分方程的求解问题转化为求解代数方程组的问题。同一个问题，采用不同的差分格式，将得到不同的代数方程组，对代数方程组求解就可得到定解在差分网格点上的近似解，从而实现对非恒定流问题的求解。
对拟线性双曲型偏微分方程组的数值求解一般有两种最常用的方法：一种是特征线法；另一种是隐式差分法。特征线法是通过特征线理论［4-5］是应用最为广泛的隐式差分格式。
笔者将针对代表复杂输水系统水力过渡的无压流、有压流、无压接有压、明满流交替等不同非恒定流过程，对特征线法和隐式差分法进行对比，并对其在不同输水方式和水流状态下的适用性进行分析和研究。
1 非恒定流基本方程的数值解法
1.1 特征线法
将非恒定流的控制方程用一个未知因子λ进行线性组合得到：
令λ=±1/c，可得到以下两组常微分方程：
以上的特征线方法在x-t平面上可以表示成如图1所示的样子，式
称为正负特征线，沿各自的特征线其相容性方程成立。对C+和C-方程组中的常微分方程沿各自的特征线积分，可得：
计算时假定初始时刻所有截面上的v和y值都是已知的，以上方程组连同端点的已知边界条件，就可以依次计算出Δt以后各个截面上的vP和yP值。为保持计算稳定性，方程必须满足库朗（Courant）稳定判据Δt≤Δx/（｜v｜+c）。
1.2 Preissmann四点隐格式
图1 特征线差分网格
Preissmann四点隐格式可表示成如图2的网格，整个管道分成m段、m+1个断面，n为已知时层，n+1为待求时层。差分格式中f代表y和v，f（x，t）是控制方程中偏导数系数和源项，θ=Δt′/Δt（0.5＜θ≤1.0）为时间加权因子。理论上隐格式是无条件稳定的，实际经验指出0＜θ≤0.5，难稳定；0.5＜θ≤0.6，弱稳定；0.6＜θ≤1.0，强稳定［6］。
将式（8）、式（9）、式（10）代入非恒定流控制方程式（1）和式（2），经整理可得：
图2 Preissmann四点隐格式
其中D=A/B，A为过水面积，B为水面宽。式中共4个未知量
、
、
、
，假设t时刻，各计算断面处的水力要素均已求出，循环进入t+Δt时刻，1，2，…，i，…，n+1各断面的流量与水深未知，故n+1个断面共有2n+2个未知数；第i个网格单元满足式（11）、式（12）两式，故n个单元共有2n个方程，而首、末端边界条件又可列出两个方程，共有2n+2个方程，方程封闭，可解出t+Δt时刻每一个断面处的流量与水深。
2 特征线法与Preissmann四点隐格式的比较与分析
特征线法和Preissmann四点隐格式作为计算非恒定流最常用的两种数值方法，在精度、稳定性、计算速度、程序实现难易程度上存在差异。特征线法在差分处理上是一种显式的方法，一般，特征线法和Preissmann四点隐格式有以下几点区别，这些区别同时也是显式差分法和隐式差分法的不同之处。
（1）稳定性：特征线法和显式差分法的稳定性均受库朗稳定判据限制，时间步长的选取受到了限制；而隐式法是无条件稳定的差分格式，时间步长选取不受稳定性的限制。
（2）计算速度：特征线法和显式差分方法每个时步都是求解简单的代数方程，故计算速度很快；而隐式差分法每个时步均需求解大型联立代数方程组，所花机时较长。
（3）计算机内存占用：因隐式法需要存储大型矩阵，故特征线法和显式差分方法占用计算机内存比隐式法少。
（4）程序设计：特征线法和显式差分法的程序设计比隐式差分法简单，这种程序设计上的便利更多体现在边界条件的处理上。
（5）计算精度：特征线法及一般的显式差分具有一阶精度，而Preissmann四点隐格式在加权因子取0.5时，具有二阶精度，但属于弱稳定状态，否则具有一阶精度。对于工程应用而言，一阶精度已可得到满意的计算结果。
以上这些特点是就一般而言的，不同数值方法在处理不同非恒定流问题时会体现出各自的特点及不同程度的适用性。下面将针对有压流、无压流、无压接有压和明满流交替等4种输水方式和输水形态，通过对特征线法和Preissmann四点隐格式在处理非恒定流问题时的比较，对这两种方法的特点及在求解不同输水形式时的适用性问题进行分析研究。
2.1 基于有压非恒定流的比较与分析
输水系统有压非恒定流的特点是在过渡过程中，输水管路始终处于有压状态，管路水击波速大，水力要素变化迅速，波动频率高，水力要素在波动过程中往往带有尖峰。本例设置的有压输水模型如下，一孔有压箱涵尺寸为4.4m×4.4m，长16114m，糙率0.0135，水击波速取800m/s。首端边界：首端流量12.33m3/s线性增为20.0m3/s，调节时间300s；末端边界：末端为堰，流量系数0.52，堰宽4.4m，堰顶高程4.1m，堰前池面积42.24m2。
图3 有压管道进口的压力波动
分别采用特征线法和Preissmann四点隐格式编制程序，对有压管路的过渡过程进行模拟。图3所示的是时步分别取1.0s和3.0s时，特征线法及Preissmann四点隐格式得到的管道进口压力波动过程。表1是不同时步下两种方法得到的压力最大值。从以上结果可看出随着时间步长的增加，压力波动衰减的速度增加，最大峰值逐渐消减。特征线法虽然不同时步的计算结果的波动衰减速度有差异，但在压力波动初期，不同时步计算结果基本一致，得到的压力波动的最大峰值相差很小，特别是时步0.5s和1.0s得到的压力最大值只相差3cm。在小时步时，Preissmann四点隐格式的计算结果与特征线法计算结果基本一致，但随时间步长的增加，压力波动衰减的速度较特征线法快。且可以看出随时间步长的增加，隐式法的计算结果在压力波动初期即出现了明显的削峰现象，影响了压力波动的最大值和最小值的计算精度。
表1 有压管道进口的压力最大值 单位：m
输水系统有压非恒定流的水力特点适合用较小时步来模拟，这种小时间步长很容易满足库朗准则，不会对特征线法或显式法的计算造成影响，且特征线法计算速度快的特点也能在此时发挥出来。而隐式法如采用较小时步，因其本身需求解大型联立方程组，故其计算速度将会变得很慢。另外，对于有压非恒定流而言是带有水力要素尖峰的过渡流，工程上所关心的往往是在过渡过程中所出现的水力要素的最大值和最小值。而从特征线法计算的结果可以看出，采用不同的时间步长，虽然水力波动衰减速度不同，但计算所得到的水力要素的最大值和最小值相差不大，在时步小于1s时一般不会出现削峰现象。但是，从Preissmann四点隐格式的计算结果来看，如果隐式法的时步稍稍取大的话，就将抹去尖峰，削峰现象非常明显，预测的水力波动的最值偏小，不利于工程安全，但如减小时步又会遇到计算速度的问题。
综上所述，对于输水系统有压非恒定流的模拟宜采用特征线法或显式差分法，在进行时间步长的选择时，可以根据具体的情况，通过试算选择满足精度要求的较大的时间步长。
2.2 基于无压非恒定流的比较与分析
输水系统的无压非恒定流特点是在过渡过程中，水流始终处于无压状态，无压流的水面波速相对于有压流的水击波速来说要小得多，加之渠道自身的调蓄作用，无压非恒定流的水力要素变化较有压流要平缓得多，波动频率较低。本例设置的无压输水模型条件如下，输水明渠长16114m，底宽8m，边坡系数3，糙率0.0135，底坡0.00009。
首端边界：首端流量由12.33m3/s减为8.33m3/s，调节时间为1200s。
末端边界：末端为堰，流量系数0.52，堰宽4.4m，堰顶高程0.0m，堰前池面积126.0m2。
图4 明渠进口的水位波动
分别采用特征线法和Preissmann四点隐格式法，对其过渡过程进行模拟。图4所示的是时步分别取4.0s和50.0s时，特征线法及Preissmann四点隐格式的计算结果，可以看出各种计算结果的水位波动趋势是相同的，但特征线法随时间步长的增加，沿线的水位、流量会出现数值上的衰减。这主要是因为特征线法中采用了插值网格造成的。如图1差分网格所示，明渠特征线无法相交到A、B两点，须通过插值求得交点的水力参数。插值法的好处在于可不用调整波速或管长而满足多管系统里共同的时间增量的要求。主要缺点是给结果带来了人为的数值衰减。从3.1节的计算结果来看，当用插值特征线法计算有压非恒定流时，未出现数值衰减的情况，主要是因为有压特征线方程，流速v在数值上与水击波速c相比，相差很大，特征线斜率主要由水击波速决定，特征线与差分网格的交点均接近于A、B两点，从而减小了插值误差。当采用插值特征线法计算无压非恒定流时，可采用增加分段数的方法，将插值误差降到最小，从图4结果来看，特征线法采用4.0s的时步所得到的水位波动与隐式法计算结果吻合很好。从Preissmann四点隐格式的计算结果看，时步4.0s和50.0s的计算结果基本一致，未出现时步增大，精度降低的情况。
综上所述，输水系统无压非恒定流的水力特点适合用较大的时步来模拟，而隐式法因其计算速度慢适合采用较大时间步长，且该法的空间步长不受时间步长的影响，时间步长的适当增大不会对计算精度产生影响，是模拟输水系统无压过渡过程的一种理想方法。而特征线法由于库朗准则的限制，在增大时间步长的同时，必须同时增大空间步长，且为避免插值所造成的数值衰减，其时间步长不能取得太大，但是由于特征线法或显式差分法计算过程非常简单，计算速度快，所以选择好合适的时间步长也可以得到令人满意的结果。
2.3 基于无压接有压的非恒定流的比较与分析
在实际工程中，由于地形和水力条件的限制，输水系统有时采取无压流和有压流相结合的输水形式，一般通过调节池将无压与有压部分衔接起来。这种输水形式的过渡过程包括无压流和有压流两部分，其特点在于有压流的水击波速通常是无压流水面波速的几十倍，乃至百倍。对于这种输水形式通常采用Preissmann窄缝假设对无压和有压两种流态统一进行求解。本例设置的无压接有压的非恒定流模型条件如下，无压段长4500.983m，箱涵尺寸4.3m×4.3m，3孔并联，底坡0.00007758，糙率0.0135；有压段长11987.0m，箱涵边长4.3m×4.3m，3孔并联，底坡0.00007758，糙率0.0135，水击波速取800m/s。
首端边界：首端流量由30.0m3/s减为20.0m3/s，调节时间为600s。
末端边界：末端有压暗渠出口是一个堰，流量系数0.45，堰宽16.667m，堰顶高程1.2m，堰前池面积166.667m2。
分别采用特征线法和Preissmann四点隐格式，对无压接有压的过渡过程进行模拟。图5、图6所示的是分别用特征线法和Preissmann四点隐格式得到的无压段进口水位以及有压段进口压力过程线，可以看到两种方法得到的结果是一致的。说明通过选择各自适宜的时间步长，两种方法都可求解无压接有压的非恒定流过程。但无论是特征线法或显式法还是隐式法，有时都会在处理无压接有压的非恒定流时，遇到时间步长选择的困难。因为无压流的水面波速小，水力要素波动缓慢，宜采用较大的时间步长，而有压流的水击波速大，水力要素变化剧烈，宜采用较小的时间步长。根据前面的分析结果已知，对于隐式法来说，如果采用较小的时间步长将会使得计算速度变得很慢，如果采用较大的时间步长又有可能使得有压瞬变的计算结果出现削峰的情况，影响精度。而对于特征线法或显式法而言，因为有压流的水击波速通常是无压流水面波速的几十倍，乃至百倍，如本例，初始时刻无压流的水面波速为4.54m/s，而有压流的水击波速为800m/s，二者相差170多倍，根据库朗准则的限制，如无压流与有压流采用统一的时间步长求解，那么有压段的空间步长就应该是无压段空间步长的100多倍，因此协调无压与有压二者的空间步长及统一的时间步长将变得比较困难。如果输水系统中存在长度较短的有压管道，就会使得有压管道的空间步长受限，从而也限制了其时间步长。
图5 用两种方法得到的无压进口水位
图6 用两种方法得到的有压进口压力
针对无压接有压的输水系统中，无压段与有压段的波速相差较大的特性，可通过变时间步长法［7］来处理。因为无压流宜采用较大的时步，而有压流宜采用较小的时步，故对于无压流可采用同一个时间步长Δt1和相应的空间步长Δx1，而有压流采用同一个时间步长Δt2和相应的空间步长Δx2，对于有压系统中长度较短的有压管道也可以单独采用较短的时间步长。这样对同一段箱涵内部的结点仍采用等时间步长和等空间步长的计算方法，而在无压与有压相结合、长管与短管相结合的变时间步长的地方需要进行特别处理。
如图7所示，在变时间步长处采用线性插值，在求下一时刻的待求结点的值时，如右边有压段S点的特征线正好处于结点上，而过待求点的左边特征线不在结点上，就做无压段的线性插值，得出相应R点水力参数，然后求解待求结点的有压段进口的水力参数。如果左边无压段R点的特征线正好处于结点上，而过待求点的右边特征线不在结点上，就做有压段的线性插值，得出相应S点水力参数，然后求解待求结点的无压段出口的水力参数。随着时间的推进及时间步长的交替，将完成整个过渡过程的求解。该方法很好地适应无压接有压的水力特性，提高了计算精度，并大大节省了计算时间。
图8是分别用特征线定时步方法和特征线变时步方法得到的有压段进口的水位过程线。定时步方法采用0.5s的时间步长，变时步方法的有压段仍采用0.5s的时间步长，而无压段采用3.0s的时间步长，通过比较可以看出两种方法计算的结果吻合很好，从而也证明了该方法处理无压接有压非恒定流问题的可行性。
综上所述，特征线法或显式法以及隐式法都能处理无压接有压的非恒定流问题，但这些方法有时还是会遇到时间步长选择的困难。这时采用变时步法，是一个很好的选择。对于瞬变过程缓慢的无压流采用较大的时间步长，而对于变化剧烈的有压流则采用较小的时间步长，通过对无压与有压衔接断面的处理，实现不同时间步长对无压接有压非恒定流的模拟，既保证了求解的精度，又提高了计算的速度。
图7 变时步插值网格图
图8 用两种方法得到的有压进口水位
2.4 基于明满流交替的比较与分析
明满流交替是指输水系统中同一断面在不同时刻可能出现无压流也可能出现有压流的水力现象。如某些无压暗渠，在发生过渡过程时，因流量增大，洞内水位超过了洞顶，水流由无压状态转变为有压状态，而后流量减小，洞内水头又低于洞顶，又出现有压流向无压流的转变。这种明满流交替的现象在水工隧洞、输水管道、水电站尾水系统均有可能发生，且其水力特性较一般的无压流或有压流都更为复杂。
早在20世纪30年代，人们就开始了对明满流交替的关注和研究。Priessman于1961年提出了虚拟窄缝法使得无压流和有压流可以统一描述，无需再追踪明满流的交界面，使得明满流交替模拟变得方便。
特征线法或显式法在模拟明满流交替时有其自身不可避免的局限性，这是因为在Priessman窄缝假设中，当管路的水位低于管顶高程时，就认为水流处于无压状态，在非恒定流的模拟中采用无压流的水面波速，而当水位高于管顶高程时，认为水流进入有压状态，则采用有压流的水击波速。水击波速通常是无压流水面波速的几十倍，乃至百倍，因库朗准则的限制，当某断面由无压变为有压流时，原本适宜的时间步长，将因水击波速被引入非恒定流控制方程，而减小为原来的几十分之一，甚至更小。虽然整个管路有些部分处于有压状态，有些处于无压状态，但都将在这样小的时间步长下运算，这将令计算速度变得很慢，有时甚至还会出现数值振荡的现象。对于插值特征线法还易出现插值所带来的数值衰减。
而隐式法因其自身无条件稳定的特点，使得其在求解明满流交替时有着显式法不可比拟的优势。管路的任意断面不论处于有压状态，还是无压状态，隐式法都将以其确定的时间步长和空间步长进行稳定运算。
本节采用2.1节所设计的管道模型，首端初始流量为12.33m3/s，在一定时间内线性减为零，如果流量减小得过快，那么，有压箱涵首端就会出现剧烈的水力波动，甚至会造成局部的短暂脱空，出现明满流交替的水力现象。用Preissmann四点隐格式计算的首端关闸的过渡过程，时间步长取1.0s，可以看出当流量减小时间为1800s时，在水力过渡过程中，水流始终处于有压状态；而当流量减小时间减为1200s时，水力波动变得剧烈，箱涵进口出现了明满流交替的水力现象。该计算实例表明，隐式差分法能很好地模拟明满流交替现象。
3 结论
（1）不同数值解法在求解不同复杂输水系统水力过渡问题时会体现出各自的特点及不同程度的适用性。本文针对代表复杂输水系统水力过渡过程的无压流、有压流、无压接有压、明满流交替等不同非恒定流过程，对特征线法和隐式差分法进行对比，并对其在不同输水方式和水流状态下的适用性进行了分析和研究。
（2）输水系统有压瞬变流水力要素变化快，宜用较小时步。特征线法或显式法可以发挥其计算速度快，编程方便的优势，而隐式法因其计算速度及存在的削峰现象表现出一定的局限性。
（3）输水系统无压瞬变流水力要素变化慢，宜采用较大时步。隐式法可以发挥其无条件稳定的优点，采用较大时步可取得理想的效果。特征线法采用较大时步，会出现插值带来的数值衰减，但是因其计算速度快，所以选择好合适的时间步长也可以得到令人满意的结果。
（4）特征线法或显式法，以及隐式法都可处理无压接有压的输水系统水力过渡问题，但均会遇到时间步长选择的困难。针对无压与有压的水力特点，采用变时步法，是一个很好的选择。
（5）明满流交替的无压与有压界面不固定，因库朗准则的限制，特征线法或显式法具有很大的局限性，而隐式法无条件稳定，在明满流的模拟上有着特征线法或显式法不可比拟的优势。
所得结论可以为研究复杂输水系统水力过渡问题模拟时数值方法的选择及步长的处理提供一定的理论依据和借鉴。
参考文献
［1］E.Benjamin Wylie，Victor L.Streeter.Fluid transients［M］.New York：McGraw-Hill International Book Co.，1978.
［2］Mei-Ling Chen，Georges.D.Nonlinear optimal control of an open-channel hydraulic system based on an infinite-dimensional model［C］//Decision and control：Proceedings of the 38th IEEE Conference.Phoenix，Arizona，1999，-5（12）：4313-4318.
［3］Streeter，V.L.，Wylie，E.B.，Hydraulic Transients［M］.New York：McGraw-Hill Book Co.，1967.
［4］Abbot M B，Havno K，Lindberg S.The fourth generation of numerical modeling in hydraulics［J］.Journal of Hydraulic Research，1991，29（5）：581-600.
［5］杨敏，李强，李琳，等.有压管道充水过程数值模拟［J］.水力学报，2007，38（2）：171-175.
［6］Abbot M B.Computational hydraulics［M］.London，UK：Pitman Publishing Ltd.，1979.
［7］练继建，王俊，万五一，等.变时步的特征线法计算复杂输水系统的水力过渡过程［J］.水利水电技术，2003，34（9）：12-14.