周 斌

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

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

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

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

(1)

式中s——流向方向的坐標(biāo)；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公式，對(duì)明槽恒定漸變流分段求解，對(duì)任意計(jì)算段Δ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)生相應(yīng)的誤差(分別記為εnk1、εnk2、εnk3、εnk4)，并在式(5)的左邊產(chǎn)生相應(yīng)的誤差εn+1，則式(5)可寫為：

(10)

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

(11)

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

(12)

(13)

(14)

(15)

(16)

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

J(s,h)

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

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

(17)

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

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

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

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

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

5 結(jié)語

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

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