顯式Runge-Kutta法求解明槽恒定漸變流的穩(wěn)定性研究

周 斌

(汕尾市水利水電規(guī)劃設計院，廣東 汕尾 516600)

明槽恒定漸變流是水力學的經(jīng)典問題，明槽通過一定流量時，由于底坡、上下游進出流邊界條件差異及明槽內(nèi)建筑物所形成的控制水深不同，明槽中的水流可以形成多種形式的水面線，并可采用微分方程表述[1]。明槽恒定漸變流微分方程工程上常采用隱式尤拉公式差分方程求解，即為應用廣泛的能量方程[2]。盡管眾多學者對能量方程進行了廣泛的研究以提高其計算精度，但由于尤拉公式固有的低精度，使得算法精度提升潛力有限。顯式Runge-Kutta法是一種在水利工程計算中應用較廣泛的常微分方程數(shù)值解法，國內(nèi)有不少學者嘗試將其用于求解明槽恒定漸變流微分方程[3-5]。顯式Runge-Kutta法是一種常微分方程的數(shù)值解法，具有高精度、計算簡便的特點[6]，有一定的推廣價值，可作為常用的能量方程算法的一個有益補充。顯式Runge-Kutta法是條件穩(wěn)定的[6]，若未嚴格控制計算步長，成果容易失真。有學者在顯式Runge-Kutta算法引入了自動選擇步長算法[6]以防計算成果失真[7]。究其原因，成果失真大多是由于計算失穩(wěn)產(chǎn)生。因而查明算法的穩(wěn)定性對明槽恒定漸變流顯式Runge-Kutta法的推廣運用是十分必要的。采用多元函數(shù)泰勒公式可推導顯式Runge-Kutta法誤差傳播方程[8]，得到明槽恒定漸變流顯式Runge-Kutta法的誤差傳播方程，以驗證和控制算法的穩(wěn)定性。

1 明槽恒定漸變流水面線微分方程[2]和4階顯式Runge-Kutt公式

明槽恒定漸變流滿足如下微分方程：

(1)

式中s——流向方向的坐標；h——水深；i——明槽底坡；J——水力坡降；Fr——弗汝得數(shù)。

為方便討論，將i、J、Fr寫為函數(shù)的形式，有：

(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)的參數(shù)k1、k2、k3、k4必將產(chǎn)生相應的誤差(分別記為εnk1、εnk2、εnk3、εnk4)，并在式(5)的左邊產(chǎn)生相應的誤差εn+1，則式(5)可寫為：

(10)

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

(11)

對式(6)—(9)引入多元函數(shù)泰勒公式并略去高階微量，則有：

(12)

(13)

(14)

(15)

(16)

則求解式(16)所需的其余各相關公式可羅列如下：

J(s,h)

3 4階顯式Runge-Kutta算法的穩(wěn)定性分析

要使得水面線計算過程穩(wěn)定，需使得|η|≤1[6]，即：

(17)

綜上所述，基于計算的穩(wěn)定性，顯式Runge-Kutta法計算明槽恒定漸變流時，緩流宜從下游向上游推算，急流時宜從上游向下游推算。顯式Runge-Kutta法的計算方向與穩(wěn)定性的關系與應用廣泛的能量方程法[9]是完全一致的。

從上述推導過程也易知，“計算明槽恒定漸變流時，緩流宜從下游向上游推算，急流時宜從上游向下游推算”的結(jié)論對一般顯式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(個別誤差過大者酌減)試向上、下游方向各推算幾組不同水深的單段水面線及其誤差傳播系數(shù)η，成果見表2；下段取段長|Δs|=0.05 m試向上、下游方向推算幾組不同水深的單段水面線及其誤差傳播系數(shù)η，成果見表3。

表2 上段不同水深推算單段水面線及誤差傳播系數(shù)η

表3 下段不同水深推算單段水面線及誤差傳播系數(shù)η

從表2、3可見，緩流從下游向上游計算水面線和急流從上游向下游計算水面線時，均有|η|≤1，計算是穩(wěn)定的；反之，均有|η|>1，計算是不穩(wěn)定的。

5 結(jié)語

顯式Runge-Kutta法也可以作為水面線計算的一種數(shù)值方法，近年來得到了一些學者的研究和推廣。作為一種數(shù)值算法，其計算應保持穩(wěn)定，以防止誤差噪音“淹沒”真解。采用顯式Runge-Kutta法計算水面線時，按照流態(tài)選擇合適的計算方向，結(jié)合誤差傳播系數(shù)η合理選擇計算步長|Δs|，可以嚴格控制計算的穩(wěn)定性；同時，當沿著顯式Runge-Kutta法適宜的計算方向計算時，隨著|Δs|的增加，誤差傳播系數(shù)η會出現(xiàn)1→0→-∞的變化規(guī)律。

顯式Runge-Kutta法是一種高精度算法，不需要迭代試算就能得到高質(zhì)量的解答，特別適合渠道等斷面尺寸變化不大的明槽水面線計算，可供同行參考選用。