黃逸軍，劉韓生

（西北農林科技大學水利與建筑工程學院，陜西楊凌 712100）

TVD格式在上游調壓室水擊數(shù)值模擬中的應用

黃逸軍，劉韓生

（西北農林科技大學水利與建筑工程學院，陜西楊凌 712100）

為了模擬含上游調壓室水電站的水力過渡過程，采用TVD格式求解水擊方程，并結合調壓室系統(tǒng)的基本方程組和給定的邊界條件來模擬水擊壓力變化。將此方法用于工程實例，其計算所得的水頭與實測水頭相吻合;并將計算所得的水頭與MacCormark格式和Lax-Friedrichs格式的計算結果比較，結果表明:TVD格式能有效模擬含上游調壓室的水擊現(xiàn)象，其結果具有耗散性低，不產生虛假數(shù)值震蕩，計算精度高等優(yōu)點。

水擊數(shù)值模擬;TVD格式;水力學;調壓室

壓力水管中因水體流速瞬時變化（如突然關閉或打開閥門）引起管內水體壓力波動，產生很大的撞擊力，這種瞬時壓力波動現(xiàn)象稱為水擊，調壓室普遍應用于水電站中就是用以減小水擊壓強的。水電站閥門的開閉往往引起下游壓力管道、上游引水隧洞以及調壓室等中水體的非恒定流動。壓力管道內的水擊方程屬于擬線性雙曲型偏微分方程，TVD格式是一種高分辨率差分格式，能夠很好地求解水擊方程?，F(xiàn)采用TVD格式求解水擊方程，再結合調壓室基本方程組計算，這樣就能較為完整地模擬上游調壓室系統(tǒng)的水力過渡過程。

1 二階TVD格式

在眾多的水擊計算數(shù)值解法中，TVD格式是較為成功的方法。定義函數(shù)u的總變差為:TV（un）=，如果某種差分格式滿足TV（un+1）≤TV （un），則稱之為TVD格式。TVD格式在間斷附近不會產生虛假數(shù)值震蕩，且對間斷有高分辨率。因為TVD格式屬于激光捕捉法，故而需要有統(tǒng)一的守恒格式。以往的研究對TVD格式做了詳細的介紹，并且將水擊控制方程推導出適用于TVD格式的嚴格守恒形式［1］。本文在此基礎上，進一步將TVD格式求解水擊方程應用于上游帶調壓室的電站水擊計算中。

考慮標量形式的雙曲型守恒律方程，TVD的具體形式為［2］

式中f為物理通量。式（1）的一般顯式格式可寫成

式中，時空步長比λ=Δt/Δx;數(shù)值通量h的形式為

式中:r=Δi+1/2-σu/Δi+1/2u，Δi+1/2u=ui+1-ui，σ=sgn （ai+1/2）;φ（r）為通量限制函數(shù)。

本文選用Roe的Superbee型限量函數(shù):φ（r）= max［0，min（2r，1），min（r，2）］進行計算。

2 水擊方程的二階TVD格式

水擊控制方程是一組擬線性雙曲型方程組:

式中:H為測壓管水頭;v為斷面平均流速;θ為管道傾斜度;J=λv2/2D，其中λ為沿程水頭損失系數(shù)，D為管道直徑;c為水擊波波速。

將水擊控制方程轉化為嚴格守恒形式［3］:

上述標量方程的TVD格式是在標量常系數(shù)方程上定義的，而要將此格式運用到方程組式（6）上，就要經過特征理論解耦［4］，即對方程組的系數(shù)加以局部凍結，然后應用單一方程的TVD格式進行推廣。

現(xiàn)假設Y=R-1U，則式（6）解耦為

用αl取代式（3）中的r的表達式中的u，即得式（13）中r的表達式?，F(xiàn)在將式（7）至式（10）以及式（6）中的U，S和F的表達式代入式（11），式（12）以及式（13）中就得到了水擊方程的二階TVD顯式格式，再結合一定的邊界條件和初始條件就可以進行求解。

3 邊界條件

3.1 上游調壓室邊界條件

壓力管道上游為調壓室系統(tǒng)。調壓室系統(tǒng)的基本方程為［5］

式中:z為調壓室水面相對基準面的高度，規(guī)定向上為正;AT，AS分別為上游引水隧洞的橫斷面積和調壓室的橫斷面積;Q為水輪機的引用流量;L為上游引水隧洞的長度;v為上游引水隧洞的水流流速;ζc為上游引水隧洞的水頭損失系數(shù);hf為阻抗孔處的水頭損失。

對于式（14）和式（15）可以采用四階龍格庫塔法進行求解。

3.2 下游邊界條件

給定下游的邊界條件為［6］

式中:vm為閥門全開時的管中流速;H0為恒定流時管道末端的作用水頭;τj為j時刻管道末端的閥門相對開度。根據(jù)特征線理論，沿著C+特征線有

4 初始條件

初始條件是指水擊現(xiàn)象發(fā)生前，恒定流動時管道中的水頭H0和流速v0，可以通過恒定流的水力計算確定。

5 工程實例

某水電站位于陜西省嵐皋縣境內的嵐河干流上，水庫正常水位為512.5 m。調壓室為阻抗式，壓力引水隧洞長2 715.213 m，過水洞徑6.0 m，調壓室內徑為14.0 m。阻抗孔口直徑為2.94 m，底板高程為460.524 m。調壓室后接壓力埋管，內徑為5.5 m，壓力鋼管主管長134.281 m，在進入主廠房約102 m處分為3條支管，支管內徑為3.0 m［7］。

現(xiàn)以3臺機組全部同時甩荷作為計算工況，機組引用流量由27.6 m3/s逐步減少到0。導葉關閉時間為7 s。將導葉關閉規(guī)律作簡化處理，近似地認為開度和時間的關系式為

分別用一階Lax-Friedrichs格式，二階Mac-Cormack格式以及上述TVD格式進行計算，計算所得的末端閥門水頭線與實測末端閥門水頭線的比較如圖1和圖2所示。

6 結果分析

從圖1與圖2中可以看到，TVD格式計算所得的末端閥門水頭線與工程的實測水頭線比較符合。從表1中可以看出，一階Lax-Friedrichs格式的計算所得的末端閥門水頭線存在較大的誤差，部分時段的誤差已超過10%。雖然二階MacCormack格式的計算結果與TVD格式的計算結果二者相對誤差相差不大，但從圖2中可以看到，二階MacCormack格式計算所得的末端閥門水頭存在明顯的虛假數(shù)值震蕩問題。表2中顯示，采用此3種格式計算所得的水頭在波峰處和波谷處的值的相對誤差相差不大。

圖1 L-F格式、TVD格式和實測三者末端閥門斷面水擊壓力水頭線Fig.1Curves of water hammer pressure head at the tail valve’s cross-section obtained from the measured data and the calculation by L-F scheme and TVD scheme

圖2 MC格式、TVD格式和實測三者末端閥門斷面水擊壓力水頭線Fig.2Curves of water hammer pressure head at the tail valve’s cross-section obtained from the measured data and the calculation by MC scheme and TVD scheme

表1 各格式計算所得水頭的誤差比較Table 1Comparison of errors in the calculations of different schemes

表2 t=7 s后各格式計算所得的水頭在波峰、波谷處的誤差比較Table 2Comparison of errors in the wave peaks and troughs in the calculations of different schemes（after t=7 s）

7 結語

（1）用TVD格式求解水擊方程，再結合調壓室基本方程和下游邊界條件，計算得到的水頭與實測水頭能較好地吻合，且誤差較小。這說明該方法能夠較好地模擬有上游調壓室的水電站的水擊壓力變化。

（2）TVD格式相對于一階Lax-Friedrichs格式和二階MacCormack格式的優(yōu)點在于:其計算得到的水頭不發(fā)生較大的耗散，也不存在虛假數(shù)值震蕩。這說明TVD格式具有高分辨率、高精度、低耗散、無虛假震蕩的優(yōu)點，具有較強的激波捕捉能力，能夠有效地模擬水擊現(xiàn)象。

［1］劉韓生，樊書剛.TVD格式在水擊數(shù)值模擬中的應用［J］.水力發(fā)電學報，2010，29（4）:108-112.（LIU Han-sheng，F(xiàn)AN Shu-gang.Numerical Simulation of Water hammer with TVD Scheme［J］.Journal of Hydroelectric Engineering，2010，29（4）:108-112.（in Chinese））

［2］HARTEN A.High Resolution Schemes for Hyperbolic Conservation Laws［J］.Journal of Computational Physics，1983，49（3）:357-393.

［3］張丹，劉韓生，李順兵.水擊的數(shù)值模擬方法比較［J］.人民長江，2008，39（18）:75-76.（ZHANG Dan，LIU Han-sheng，LI Shun-bing.Comparison of Numerical Simulation Methods for Water-Hammer［J］.Yangtze River，2008，39（18）:75-76.（in Chinese））

［4］余德浩，湯華中.微分方程數(shù)值解法［M］.北京:科學出版社，2003:181-279.（YU De-hao，TANG Huazhong.Numerical Solution of Differential Equations［M］. Beijing:Science Press，2003:181-279.（in Chinese））

［5］呂宏興.水力學［M］.北京:中國農業(yè)出版社，2002: 296-335.（LV Hong-xing.Hydraulics［M］.Beijing: China Agriculture Press，2002:296-335.（in Chinese））

［6］彭守拙.上游氣墊調壓室系統(tǒng)過渡過程的計算［J］.水力發(fā)電學報，2002，（4）:62-70.（PENG Shou-zhuo. Calculation of Transient Process for Upstream Air Cushion Surge Chamber System［J］.Journal of Hydroelectric Engineering，2002，（4）:62-70.（in Chinese））

［7］張曉宏.氣墊式調壓室過渡過程的計算研究［D］.西安:西安理工大學，2007.（ZHANG Xiao-hong.Calculable Research on Fluid Transients of Air Cushion Surge Chamber［D］.Xi’an:Xi’an University of Technology，2007.（in Chinese））

（編輯:劉運飛）

Numerical Simulation of Water Hammer with TVD Scheme in the Upstream Surge Chamber

HUANG Yi-jun，LIU Han-sheng
（College of Water Resources and Architectural Engineering，Northwest A＆F University，Yangling712100，China）

In order to simulate the hydraulic transient of the hydropower station which has upstream surge chamber，a TVD scheme is employed to solve the water hammer equations.In association with the basic equations of the surge chamber system and the given boundary conditions，the variation of water hammer pressure is simulated.The approach is applied to engineering practice and the calculated waterhead is consistent with measured data.Furthermore，the calculation result is compared with those calculated by using MacCormark and Lax-Friedrichs schemes. It’s revealed that the TVD scheme can effectively simulate the water hammer effect in the upstream surge chamber，and the calculation is less dissipative，non-oscillatory，and highly precise.

numerical simulation of water hammer;TVD scheme;hydraulics;surge chamber

TV131.4

A

1001-5485（2013）05-0047-04

10.3969/j.issn.1001-5485.2013.05.011

2013，30（05）:47-50

2012-09-28;

2012-11-07

黃逸軍（1988-），男，江西撫州人，碩士研究生，從事水力學與水工建筑物方面的研究，電話（15891398953）電子信箱（huangyijun1988@163.com）。

劉韓生（1962-），男，陜西韓城人，教授，博士，從事水力學與水工建筑物方面的研究，（電話）13319231569（電子信箱）hanshengliu@ 126.com。