显式Runge-Kutta法求解明槽恒定渐变流的稳定性研究

周 斌

(汕尾市水利水电规划设计院，广东 汕尾 516600)

明槽恒定渐变流是水力学的经典问题，明槽通过一定流量时，由于底坡、上下游进出流边界条件差异及明槽内建筑物所形成的控制水深不同，明槽中的水流可以形成多种形式的水面线，并可采用微分方程表述[1]。明槽恒定渐变流微分方程工程上常采用隐式尤拉公式差分方程求解，即为应用广泛的能量方程[2]。尽管众多学者对能量方程进行了广泛的研究以提高其计算精度，但由于尤拉公式固有的低精度，使得算法精度提升潜力有限。显式Runge-Kutta法是一种在水利工程计算中应用较广泛的常微分方程数值解法，国内有不少学者尝试将其用于求解明槽恒定渐变流微分方程[3-5]。显式Runge-Kutta法是一种常微分方程的数值解法，具有高精度、计算简便的特点[6]，有一定的推广价值，可作为常用的能量方程算法的一个有益补充。显式Runge-Kutta法是条件稳定的[6]，若未严格控制计算步长，成果容易失真。有学者在显式Runge-Kutta算法引入了自动选择步长算法[6]以防计算成果失真[7]。究其原因，成果失真大多是由于计算失稳产生。因而查明算法的稳定性对明槽恒定渐变流显式Runge-Kutta法的推广运用是十分必要的。采用多元函数泰勒公式可推导显式Runge-Kutta法误差传播方程[8]，得到明槽恒定渐变流显式Runge-Kutta法的误差传播方程，以验证和控制算法的稳定性。

1 明槽恒定渐变流水面线微分方程[2]和4阶显式Runge-Kutt公式

明槽恒定渐变流满足如下微分方程：

(1)

式中s——流向方向的坐标；h——水深；i——明槽底坡；J——水力坡降；Fr——弗汝得数。

为方便讨论，将i、J、Fr写为函数的形式，有：

(2)

其中：

(3)

(4)

式中n(s)——s处的糙率；Q——流量；A(s,h)——s处水深h的过水面积；R(s,h)——s处水深h的水力半径；V(s,h)——s处水深h的流速；B(s,h)——s处水深h的水面宽度；g——重力加速度。

引入4阶显式Runge-Kutta公式，对明槽恒定渐变流分段求解，对任意计算段Δs有[3-5]：

(5)

k1=k(s,h)|s=sn,h=hn

(6)

(7)

(8)

k4=k(s,n)|s=sn+1,h=hn+Δs·k3

(9)

2 4阶显式Runge-Kutt公式的误差传播方程

式(5)为显式公式，已知hn可直接解得hn+1。若已知的hn存在误差εn，则式(5)的参数k1、k2、k3、k4必将产生相应的误差(分别记为εnk1、εnk2、εnk3、εnk4)，并在式(5)的左边产生相应的误差εn+1，则式(5)可写为：

(10)

式(10)中消去式(5)有：

(11)

对式(6)—(9)引入多元函数泰勒公式并略去高阶微量，则有：

(12)

(13)

(14)

(15)

(16)

则求解式(16)所需的其余各相关公式可罗列如下：

J(s,h)

3 4阶显式Runge-Kutta算法的稳定性分析

要使得水面线计算过程稳定，需使得|η|≤1[6]，即：

(17)

综上所述，基于计算的稳定性，显式Runge-Kutta法计算明槽恒定渐变流时，缓流宜从下游向上游推算，急流时宜从上游向下游推算。显式Runge-Kutta法的计算方向与稳定性的关系与应用广泛的能量方程法[9]是完全一致的。

从上述推导过程也易知，“计算明槽恒定渐变流时，缓流宜从下游向上游推算，急流时宜从上游向下游推算”的结论对一般显式Runge-Kutta法均能成立。常见各阶显式Runge-Kutta公式和误差传播方程见表1。

表1 常见显式Runge-Kutta公式和误差传播方程

4 算例

某一矩形泄槽，底宽为5 m，前段底坡为i=0.001，后段底坡i=0.25，糙率n=0.015，流量Q=30 m3/s。则其临界水深为hk=1.542 m，前段(i=0.001)的均匀流水深h0=2.465 m、后段(i=0.250)的均匀流水深为h0=0.379 m。上段取段长|Δs|=5 m(个别误差过大者酌减)试向上、下游方向各推算几组不同水深的单段水面线及其误差传播系数η，成果见表2；下段取段长|Δs|=0.05 m试向上、下游方向推算几组不同水深的单段水面线及其误差传播系数η，成果见表3。

表2 上段不同水深推算单段水面线及误差传播系数η

表3 下段不同水深推算单段水面线及误差传播系数η

从表2、3可见，缓流从下游向上游计算水面线和急流从上游向下游计算水面线时，均有|η|≤1，计算是稳定的；反之，均有|η|>1，计算是不稳定的。

5 结语

显式Runge-Kutta法也可以作为水面线计算的一种数值方法，近年来得到了一些学者的研究和推广。作为一种数值算法，其计算应保持稳定，以防止误差噪音“淹没”真解。采用显式Runge-Kutta法计算水面线时，按照流态选择合适的计算方向，结合误差传播系数η合理选择计算步长|Δs|，可以严格控制计算的稳定性；同时，当沿着显式Runge-Kutta法适宜的计算方向计算时，随着|Δs|的增加，误差传播系数η会出现1→0→-∞的变化规律。

显式Runge-Kutta法是一种高精度算法，不需要迭代试算就能得到高质量的解答，特别适合渠道等断面尺寸变化不大的明槽水面线计算，可供同行参考选用。