明渠非恒定流的一維淺水波與三維VOF耦合模擬*

張曉曦，張春澤，陳秋華

(1. 廈門理工學院 環境科學與工程學院，福建 廈門 361024； 2. 重慶交通大學 重慶西南水運工程科學研究所，重慶 400074；3. 廈門理工學院 土木工程與建筑學院，福建 廈門 361024)

明渠非恒定流的一維淺水波與三維VOF耦合模擬*

張曉曦1，張春澤2，陳秋華3

(1. 廈門理工學院 環境科學與工程學院，福建 廈門 361024； 2. 重慶交通大學 重慶西南水運工程科學研究所，重慶 400074；3. 廈門理工學院 土木工程與建筑學院，福建 廈門 361024)

提出了一種一維淺水波與三維VOF模型耦合的明渠非恒定流數值計算方法。此方法通過建立兩個計算模型間變量的數學關系，實現由三維區域為一維區域提供流量，一維區域為三維區域提供壓強和水體體積分數分布。用耦合方法模擬了矩形斷面渠道增、甩負荷兩種典型非恒定工況，將結果與一維淺水波模型計算結果對比，表明本方法在定性上正確，在定量上準確度也較高。最后利用這一方法模擬了某水電站帶引水渠壓力前池的非恒定流現象，獲得了最大溢流量、最高和最低水位等重要參數。耦合方法同時結合了一維模型計算較細長明渠漸變流的效率和三維方法模擬局部明渠急變流的準確性，可為分析工程中復雜的明渠流動問題提供更優選擇。

水利工程；明渠非恒定流；數值模擬；一維與三維耦合

0 引 言

明渠內任意一點的水位和流量隨時間變化的流動被稱為明渠非恒定流，如水電站工況調節在引水渠中激發的涌浪。這類涌浪對水電站的運行一般是有害的，如降水幅度過大可能導致與明渠銜接的壓力管道進氣，涌水幅度過大可能引起管道壓力超標，而水位和流量的交替變化則可能激發攔污柵和閘門等水力設施的振動等。因此，準確計算明渠非恒定流對水電站的設計和安全性評估具有重要意義。

到目前為止，一維和二維淺水波模型是計算明渠非恒定流的主要模型，也是研究水電站無壓引水系統涌浪的重要工具。例如：崔偉杰等[1]利用一維淺水波模型計算了某水電站長引水明渠的水力過渡過程并探討了壓力前池的數值計算方法；王明疆等[2]將這一模型線性化并計算了某水電站尾水明渠的小波動過渡過程；ZHANG Chuhze等[3]利用二維淺水波模型求解了梯級電站的銜接明渠段在電站工況轉換時的非恒定流現象，分析了部分機組甩負荷對其他機組的水力干擾。綜上，淺水波模型的應用形式多樣，計算效率高，基本可滿足大部分工程設計的需要。但這一類模型在推導過程中用到了一系列的簡化和假設，如只能用于渠道底坡接近于水平的情況、假設流速沿水深均勻分布以及壓力沿水深為靜壓分布等漸變流條件[4]，而不適用于諸如大底坡、急轉彎以及存在垂向流動等急變流情況。

為了研究明渠非恒定急變流，近年有學者開始嘗試用三維計算流體力學方法中能模擬自由水面的VOF(volume of fluid，流體體積分數)模型模擬此類問題。例如：CHENG Yongguang等[5]模擬了某變頂高尾水洞明滿流非恒定過程，探索了明滿流界面處滯留氣團生成、發展及潰滅機制；楊建東等[6]模擬了某水電站尾水洞與導流洞結合段的流動，分析了非恒定過程中尾水洞內水流中斷機理；朱方磊等[7]模擬了明渠銜接豎井的水電站的明渠非恒定流現象，探討了明渠斷流和明滿流發生的條件。綜上，三維方法適用于各種復雜流動條件且可獲得較準確的流場細節，應用潛力較傳統一維方法大。但此類方法的短板同樣顯著，即所耗費的計算機資源(內存容量和CPU計算時間)較大，目前尚難以全面應用于解決工程問題。

筆者旨在提出一種將一維淺水波與三維VOF模型耦合的明渠非恒定流計算方法，充分結合一維方法模擬細長渠道漸變流的效率和三維方法模擬局部急變流的準確性，為解決工程問題提供更優選擇。

1 數學模型及驗證

1.1 計算區域劃分

用耦合方法模擬明渠非恒定流的基本思路是將計算區域按流態分為兩部分，漸變流區域(如棱柱形渠道、無壓管道或隧洞等)采用一維方法求解，急變流區域(如壓力前池、進/出水口和渠道分岔口等區域)采用三維方法求解，兩個計算區域通過耦合邊界實時傳遞數據，如圖1。這里三維區域采用VOF模型，一維區域采用淺水波模型，下面具體介紹這兩個模型和耦合格式。

圖1 一維與三維耦合計算區域示意Fig.1 Schematic diagram of 1D and 3D coupling computational domain

1.2 三維VOF模型

VOF(volume of fluid)模型[8]是一種模擬多種互不穿透流體界面的方法。對于水氣兩相流而言，其基本思想是定義函數α表示計算網格內水流所占的比例，若α=1表示該網格完全被水充滿，α=0表示該網格完全被氣充滿，0<α<1則表示該網格內存在水氣界面。其連續性方程為

(1)

式中：u為流速；t為時間；x為坐標。在這一模型中，水相和氣相有相同的速度場和壓力場，動量方程為

(2)

以上即為VOF模型的基本控制方程，通過現有的算法(如SIMPLE、SIMPLEC或PISO算法等)，再結合模擬雷諾應力項的湍流模型(如k-ε、k-ω或LES模型等)，即可方便求解。

至于邊界條件，這里假設上游邊界流量(流速u)變化已知，下游(耦合邊界)則由耦合格式給出壓力p和水體體積分數α的分布規律。

1.3 一維淺水波模型

一維淺水波模型是描述一維明渠非恒定流的基本模型，包括連續性方程和動量方程[4]：

(3)

(4)

式中：Z為水位；Q為流量；B為渠道水面寬度；s為渠道長度坐標；A為過流斷面面積；g為重力加速度；C為謝才系數；R為水力半徑。

筆者采用顯格式有限差分法求解，差分網格如圖2。對于內節點，空間項采用二階中心差分格式，時間項采用一階差分格式，可將式(3)和式(4)轉化為差分方程式(5)和式(6)。

圖2 一維有限差分網格示意Fig.2 Schematic of 1D finite difference mesh

(5)

(6)

式中：α為權重系數(這里取為0.33)；Δt表示時間步長；Δs表示空間步長；is表示渠道坡度，變量的上標表示時間，下標表示差分網格節點編號。

這里假設下游邊界(節點N)水位條件已知，空間項采用一階向后差分得流量的差分格式如式(7)，上游(耦合邊界)邊界由耦合格式給出水位Z0和流量Q0的表達式。

(7)

1.4 耦合格式

通過以上對一維和三維區域的模型和算法分析可知，耦合邊界需實時給出三維區域下游壓強p和水體體積分數α的分布，以及一維區域上游水位Z0和流量Q0。筆者給出一種計算格式，流程圖如圖3。此格式基本思想是三維區域為一維區域提供流量，一維區域為三維區域提供壓強和水體體積分數分布。需要注意的是，耦合邊界必須同時滿足一維與三維的流動特點，即耦合邊界為漸變流態。

注：變量上標表示時間，下標表示邊界(三維區域)或節點編號(一維區域)圖3 耦合格式流程Fig.3 Flow chart of the coupling scheme

1.5 數值驗證

為了驗證上述耦合格式，這里采用耦合方法模擬圖1所示細長渠道的非恒定流，并將結果與一維淺水波模型的計算結果對比。計算區域總長100 m，斷面為寬4 m的矩形，底坡為平坡，其中上游側50 m用三維VOF模型模擬，下游側50 m用一維淺水波模型計算。渠道初始為恒定流，流量為10 m3/s，水位為5 m。這里模擬兩個工況，甩負荷和增負荷工況：甩負荷工況時，上游閘門逐漸關閉，流量在10 s內線性降到0；增負荷工況時，上游閥門逐漸開啟，流量在10 s內線性增加到20 m3/s，在這兩個過程中下游水位始終保持5 m不變。

結果表明文中一維淺水波與三維VOF耦合方法可以較為準確的模擬明渠涌浪。圖4為甩負荷和增負荷工況渠道中點(耦合邊界)流量和水位變化過程，其中實線為一維計算結果，虛線為耦合計算結果，兩者吻合較好。表1對比了兩種方法所得的涌浪周期、第一周期內水位極值和流量極值，可見除了甩負荷過程最大流量誤差稍大外(相對誤差為14.32%)，其他關鍵數據的相對誤差均小于5%。另外，以上結果也顯示出耦合方法所得波動的衰減要稍強于一維方法，這種現象一方面是由于三維方法在離散時引入的數值耗散要高于一維方法，另一方面則是因為一維水流在簡化過程中忽略了豎向和橫向能量耗散，阻尼稍小。綜上，文中耦合方法模擬細長渠道非恒定流的結果在定性上與一維方法一致，定量上與一維方法的差別也不大。考慮到一維方法已經過大量的工程實踐檢驗，這里認為本耦合方法也是正確且較準確的。

圖4 兩種計算方法所得明渠流量和水位隨時間變化過程對比Fig.4 Comparison of the flow and water level of open-channelobtained by two calculation methods changing with time

T/sZmax/mZmin/mQmax/(m3·s-1)Qmin(m3·s-1)一維淺水波29．1/41．06．30/5．864．13/4．797．47/22．26-7．20/19．37耦合方法27．8/42．86．13/5．834．14/4．796．40/21．60-7．15/18．83絕對誤差-1．3/1．8-0．17/-0．030．01/0-1．07/-0．660．05/-0．54相對誤差/%4．68/4．212．70/0．510．24/014．32/2．960．69/2．79

注：表中T表示涌浪周期，Zmax和Zmin分別表示第一周期內最高和最低水位，Qmax和Qmin分別表示第一周期內最大和最小流量，絕對誤差是耦合算法與一維方法結果的差值，相對誤差為絕對誤差的絕對值與相應一維結果的比值。

2 水電站帶引水渠的壓力前池非恒定流模擬

壓力前池是引水式水電站引水明渠與壓力鋼管結合處加寬、加深的水池，在池壁上還設置有溢流堰。這一結構一方面可以降低攔污柵處的水流流速，起到減小水頭損失和沉淀來流雜質的作用，更重要的是在電站發生過渡過程時可以宣泄多余水流，控制壓力鋼管的最大壓力。圖5為某電站的壓力前池及其附屬結構示意圖，這里由引水渠從上游水庫引水，來水經彎道擴散段進入壓力前池內，在池壁設有WES型溢流堰。其中引水渠長81 m，上游水庫額定水位為7.5 m，溢流堰頂高程為8.5 m。這里需要計算在最不利的雙機同時甩負荷的情況下壓力前池的水位波動過程，以判斷在最高水位時側堰溢流量是否滿足設計要求，壓力鋼管強度是否達標，以及最低水位是否會導致壓力鋼管進氣等。

這里采用一維與三維耦合的方法模擬壓力前池在全甩負荷工況下的非恒定流。一維計算區域包括上游水庫邊界和約51 m長引水渠。三維計算區域包括壓力前池、擴散段和部分引水渠(為保證一維與三維耦合邊界為漸變流態，三維計算區域保留約30 m長引水渠，以盡量減弱擴散段對耦合邊界流態的影響)，如圖5。

圖5 某水電站壓力前池體型示意Fig.5 Schematic diagram of the pressure fore bayof a hydropower station

全計算區域采用六面體和棱柱體混合非結構化網格離散，通過恒定流條件下的網格敏感性分析確定最終計算網格數約為110 000個。計算中用顯格式VOF模型模擬自由水面，湍流模擬采用RNGk-ε雙方程模型，計算格式采用模擬非恒定流較有優勢的PISO格式，時間步長按照庫朗穩定條件選為0.005 s。這里模擬兩機全甩負荷這一最不利工況，即1#和2#壓力鋼管流量在5 s內由54 m3/s線性降低到0，在此過程中上游水庫水位保持7.5 m不變。

甩負荷發生時，壓力管道流量減小并在前池出現向上游傳播的雍水波。在此過程中若池內水位高于溢流堰頂高程，則會發生溢流，以減緩池內水位上升。同時，由于過流斷面在橫向和縱向的擴散，池內流動的三維特性顯著。圖6為最高水位時前池內流線分布，可以看出此時前池內出現多個大范圍回流區，這種流動用三維方法模擬顯然比用傳統一維方法更為準確。圖7為甩負荷過程壓力前池及其附屬結構3個特征斷面處水位和WES堰溢流量隨時間變化過程，從中可以看出壓力前池在恒定及非恒定工況下的流動特征。恒定流時，由于流動斷面逐漸擴大，擴散段入口、前池入口和壓力鋼管入口的水位也逐漸升高，此時壓力前池內的水位約為7.6 m。甩負荷后，壓力鋼管入口的水位首先上升并在池內形成雍水波。隨著雍水波向上游傳播，壓力前池入口和擴算段入口的水位依次上升，在此過程中當池內水位高于堰頂高程時(約第16 s時刻)開始溢流(這里溢流量由溢流堰下部統計得出，因此變化略晚于池內水位)。之后在上游水庫反射波到達前池后，在池內形成降水波，并向上游傳播，出現周期性涌浪。通過本次模擬得出，在最不利甩負荷工況下壓力鋼管入口最大動水壓力水頭約為5.6 m(以鋼管軸線高程為基準)，最小淹沒水深約為1.6 m(以進水口頂板高程為基準)，溢流堰最大溢流量約為24 m3/s，均符合設計要求。

圖7 壓力前池特征斷面水位和溢流量隨時間變化過程Fig.7 Time history of water level and dischargein characteristic section of pressure fore bay

綜上，由于壓力前池體型特殊，內部水流流態，尤其是非恒定流態呈典型三維特征，用傳統一維或二維方法無法準確模擬其水動力特性。筆者提出的一維淺水波與三維VOF模型耦合的方法既可準確考慮非恒定過程中壓力前池內局部擴散、回流和溢流等復雜三維流態，又可計入引水渠中的水流慣性對局部流動的影響，所得特征參數的可信度更高。同時這一方法所消耗計算機資源和時間遠小于全三維模擬，特別適用于解決實際工程問題，具有廣闊的應用前景。

3 結 語

提出了一種一維淺水波與三維VOF模型耦合的明渠非恒定流數值計算方法，并在驗證了其準確性的基礎上模擬了某電站壓力前池甩負荷涌浪過程。

這一方法可以充分結合傳統一維模型計算較細長型明渠中漸變流的效率和三維方法模擬局部帶自由水面急變流的準確性，為分析工程中復雜的明渠流動問題提供了新的途徑。但筆者使用經過簡化的數值解驗證所提出的方法，不夠精確，后續應通過更符合實際情況的高精度水力學模型試驗對其進行更全面的驗證。

[1] 崔偉杰，張健，陳勝. 某水電站壓力前池不同模型過渡過程的計算對比[J]. 三峽大學學報(自然科學版)，2016，38(4)：36-39.

CUI Weijie，ZHANG Jian，CHEN Sheng. Calculation and comparison of transients of forebay with different models for a hydropower station[J].JournalofChinaThreeGorgesUniversity(NaturalSciences)，2016，38(4)：36-39.

[2] 王明疆，楊建東，王煌. 含明渠尾水系統小波動調節穩定性分析[J]. 水力發電學報，2015，34(1)：161-168.

WANG Mingjiang，YANG Jiandong，WANG Huang. Analysis of stability of tailrace open channel in small fluctuation governing system[J].JournalofHydroelectricEngineering，2015，34(1)：161-168.

[3] ZHANG Chunze，CHENG Yongguang，WU Jiayang，et al. Lattice boltzmann simulation of the open channel flow connecting two cascaded hydropower stations[J].JournalofHydrodynamicsSer.B，2016，28(3)：400-410.

[4] 程永光，楊建東，賴旭，等. 實用水力瞬變過程[M]. 3版.北京：中國水利水電出版社，2015.

CHENG Yongguang, YANG Jiandong, LAI Xu, et al.AppliedHydraulicTransients，ThirdEdition[M]. 3th ed. Beijing：China Water Power Press，2015.

[5] CHENG Yongguang，LI Jinping，YANG Jiandong. Free surface-pressurized flow in ceiling-sloping tailrace tunnel of hydropower plant：Simulation by VOF model[J].JournalofHydraulicResearch，2007，45(1)：88-99.

[6] 楊建東，李玲，周俊杰，等. 尾導結合的尾水系統水流中斷的機理分析[J]. 水動力學研究與進展A輯，2012，27(4)：394-400.

YANG Jiandong，LI Ling，ZHOU Junjie，et al. Study on the water flow disruption in the system of tailrace tunnel combined with diversion tunnel[J].ChineseJournalofHydrodynamics，2012，27(4)：394-400.

[7] 朱方磊，程永光. 明渠接豎井水電站引水系統的瞬變流規律[J]. 武漢大學學報(工學版)，2015，48(6)：744-750.

ZHU Fanglei，CHENG Yongguang. Transient flows in headrace system combining open channel and shaft for hydropower stations[J].EngineeringJournalofWuhanUniversity，2015，48(6)：744-750.

[8] HIRT C W，NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J].JournalofComputationalPhysics，1981，39(1)：201-225.

1D Shallow Water Wave and 3D VOF Coupling Simulationof Transient Open-Channel Flow

ZHANG Xiaoxi1，ZHANG Chunze2，CHEN Qiuhua3

(1. School of Environmental Science and Engineering，Xiamen University of Technology，Xiamen 361024，Fujian，P. R. China; 2. Chongqing Southwest Scientific Research Institute for Water Transportation Engineering，Chongqing Jiaotong University，Chongqing 400074，P. R. China; 3. School of Civil Engineering and Architecture，Xiamen University of Technology，Xiamen 361024，Fujian，P. R. China)

A numerical method for simulating the transient open-channel flow by coupling 1D shallow water wave and 3D VOF models was proposed. Firstly，by establishing the mathematical relationship between the variables of the two calculation models，the proposed method realized that 3D region provided flow for 1D region and 1D region provided the pressure and volume fraction distribution of water for 3D region. Then，two typical transient conditions of rectangular cross-section channel during load rejection and increase were simulated by coupling method. Comparing the above results with the calculation results obtained by 1D shallow water wave model，it shows that the proposed method is qualitatively correct，and the accuracy is high in quantitative. Finally，the transient open-channel flow in a pressure fore bay of a hydropower station with inlet pipes was simulated by the proposed method，and the critical parameters，such as the maximal overflow discharge，the maximal and minimal water levels were obtained. The coupling method combines the efficiency of 1D model in the calculation of a fine long canal and the accuracy of 3D model in the simulation of local open-channel blast flow，which can provide a better choice for analyzing complex open-channel flow problems in engineering.

hydraulic engineering; transient open-channel flow; numerical simulation; 1D and 3D coupling

10.3969/j.issn.1674-0696.2017.12.10

2016-10-12；

2016-12-02

廈門市科技計劃項目(3502Z20150051)；廈門理工學院高層次人才項目(YKJ15042R，YKJ14025R)；福建省中青年教師教育科研項目(JA15384)；重慶市科委自然科學基金計劃資助項目(cstc2016jcyjA1935)；重慶市教委自然科學基金計劃資助項目(KJ1600514)

張曉曦(1986—)，男，湖北襄陽人，講師，博士，主要從事水力學方面的研究。

張春澤(1986—)，男，寧夏中衛人，助理研究員，博士，主要從事計算流體力學方面的研究。E-mail: zhangchunze@whu.edu.cn。

TV131.2

A

1674-0696(2017)12-053-05

朱漢容)