Karst地區(qū)明渠伏流交替河道水動力學模型及應用

李松萍，陶 銘，金中武，胡德超

(1.華中科技大學 水電與數(shù)字化工程學院，武漢 430074; 2.長江科學院 河流研究所，武漢 430010)

1 研究背景

在我國西南Karst地區(qū)具有特殊的地理地質條件，常常出現(xiàn)明渠與伏流交替型河道。該類型河道的水流運動屬于自由水面流動、有壓流動的混合流運動，水流形態(tài)十分復雜，目前人們對于它們的認識還比較少。與此同時，研究混合流的運動規(guī)律也是研究明渠與伏流交替河道物質輸移、泥沙沖淤的基礎，具有重要科學意義。隨著我國西南Karst地區(qū)的水資源、水電開發(fā)工程逐漸增多，研究Karst地區(qū)明渠與伏流交替河道的水流規(guī)律可為工程建設提供設計依據(jù)和技術支撐，具有重要的工程應用價值。

混合流研究常采用基于圣維南方程或其變化形式的一維水動力模型模擬的方法。已開展的研究多集中在城市地下管網(wǎng)輸水領域，現(xiàn)有的混合流模型分為2類。第一類方法采用雙控制方程分別描述混合流中自由表面水流和有壓水流運動；第二類方法用單控制方程對混合流運動進行統(tǒng)一描述。Bousso等[1]等就上述2類模型優(yōu)缺點進行過比較，雙控制方程模型可以顯式追蹤水流與空氣的分界面，這類模型實現(xiàn)起來較為復雜且不具有適用于各種流態(tài)的通用性。單控制方程模型較為簡單，但常常需要借助輔助技術(Preissmann slot技術[2]、壓力分割技術[3])去描述混合流運動。此外，Casulli和Stelling[4]根據(jù)測壓管水頭在明渠、承壓水流中具有相同的表達形式，實現(xiàn)了混合流的統(tǒng)一描述和模擬。

以上混合流模型主要適用于城市管網(wǎng)中規(guī)則過流斷面(例如矩形、圓形)混合流的模擬。而真實河道斷面十分復雜，所以當這些模型用于模擬Karst河流明渠與伏流交替河道時常常遇到困難。目前，對于真實明渠與伏流交替型河道的水流模型的研究還很少。在僅有的文獻報道中，Zhang等[5]通過用雙控制方程描述的方法建立了明渠伏流一維水流模型。模型通過互提邊界條件法實現(xiàn)自由表面水流與有壓水流計算的耦合，在伏流段采用恒定流模型計算。他們使用所建立的模型模擬了真實河道不規(guī)則地形條件下的混合流。經(jīng)過文獻調研，真實河道不規(guī)則地形條件下混合流的統(tǒng)一型模型，至今尚未有文獻報道。

本文參考文獻[4]根據(jù)測壓管水頭的雙重含義實現(xiàn)混合流的統(tǒng)一描述；使用局部線性化水流控制方程代替文獻[4]的控制方程，大幅簡化了數(shù)值離散和求解；在前期規(guī)則封閉斷面混合流模型[6]的基礎上，提出三角形子斷面法描述伏流非規(guī)則封閉斷面，進而建立了真實明渠伏流交替河道的統(tǒng)一型混合流模型。以貴州六沖河夾巖水庫為背景，研究了明渠伏流交替河道縱向沿程的流速變化規(guī)律以及伏流河段前布置人工隧洞的分流規(guī)律。

2 模型簡介

2.1 控制方程

使用圣維南方程組作為描述明渠與伏流交替水流運動的控制方程組，其基本形式如式(1)、式(2)所示，其中式(1)為連續(xù)性方程，式(2)為動量方程：

(1)

(2)

式中:η表示測壓管水頭，對于封閉斷面過流通道中的承壓水流使用其本意，對于自由水面(明渠)水流而言測壓管水頭的含義即為自由水面高度[4](m)；u為斷面平均流速(m/s)；B為斷面水面寬度(m)；A為過水斷面面積(m2)；g為重力加速度(m/s2)；nm為河道綜合糙率值；R為斷面水力半徑(m)；t為時間(s)；x為距離(m)。

自由水面高度(水位)的物理本質就是測壓管水頭，前者是后者在明渠水流中的一種具體表達形式。根據(jù)測壓管水頭在明渠、承壓水流中具有相同的表達形式，即可使用式(1)、式(2)對混合流的自由表面部分和有壓部分進行統(tǒng)一描述。在明渠河段水流為自由水面流動，η為斷面水位；在伏流段(充滿)時河道為有壓流動，η為測壓管水頭。

2.2 計算網(wǎng)格與明渠及伏流的模型描述

以斷面為中心構建控制體(也稱單元)，在控制體中心布設水位變量，在控制體界面(邊)處布設流量變量。ne和ns分別用于表示單元和邊的數(shù)量，每個單元/邊都有一個唯一的全局單元/邊索引。與計算網(wǎng)格(如圖1所示)相關聯(lián)的變量定義如下。對于計算單元i，Δxi表示計算單元的長度；Bi是布置在計算單元中心的橫斷面的濕表面寬度，它與單元的水位有關并隨時間變化；Ai代表位于單元i的橫斷面的輸運面積。對于連接單元i和i+1的邊i+1/2，Δxi+1/2是單元i和i+1中心之間的距離；Ai+1/2和Ri+1/2分別為邊i+1側的輸運面積和水力半徑，由單元i和i+1的值插值得到；i+1/2和i-1/2分別表示第i個單元的后一條邊和前一條邊。

圖1 明渠與伏流交替河道計算網(wǎng)格Fig.1 Computational grids of open channeland sinking stream

伏流和明渠區(qū)域中的計算單元用標記(CLO)區(qū)分。對于伏流范圍內的計算單元，CLO設置為1；對于明渠范圍的計算單元，CLO設置為0。在伏流區(qū)域，計算單元根據(jù)其是否充滿進一步區(qū)分(用“PRE”標記表示)。對于伏流區(qū)域內承壓流動的計算單元，PRE設置為1；對于明渠段自由表面流動狀態(tài)的計算單元，PRE設置為0。在該計算范圍中，所有自由表面和承壓區(qū)域的計算單元同步進行求解。

圖2 伏流洞非規(guī)則封閉斷面的描述方法Fig.2 Descriptions of the irregular closed cross-sectionof a sinking stream

伏流非規(guī)則封閉斷面的描述方法。使用傳統(tǒng)的子斷面法描述非規(guī)則封閉斷面(圖2(a))，需要進行分區(qū)描述，這種分區(qū)描述十分復雜，且分區(qū)可能需要隨著洞內水位的變動而調整；當封閉斷面左右兩側轉折點ZfL、ZfR處于不同高度時，會引起描述困難和不準確。初步測試表明，基于傳統(tǒng)子斷面分區(qū)描述法的非規(guī)則封閉斷面混合流模型，常常出現(xiàn)計算不穩(wěn)定或不準確的問題。本文提出三角形封閉斷面描述法，來解決這一問題。當過流斷面處于充滿狀態(tài)時，過流斷面面積等于斷面最大面積；當處于非充滿狀態(tài)時，斷面過流面積的計算方法(圖2(b))：在斷面最低點處豎直向上引一條線段與水面相交，得到一個交點，依次將水下所有斷面點與該交點連接，于是形成一系列三角形，通過計算各三角形面積并進行求和，即可得到過流面積。

2.3 數(shù)學模型框架

混合使用有限差分法和有限體積法對模型的控制方程進行離散。使用半隱方法離散壓力梯度項，使用歐拉-拉格朗日法[7](ELM)求解控制方程中的對流項；聯(lián)合動量方程和連續(xù)性方程構建三對角陣求解各控制體壓力；將壓力代入動量方程確定最終流速或流量。

2.3.1 控制方程的離散

對于明渠和伏流這2種流動形態(tài)，采用θ半隱式方法[8-9]對動量方程的水位梯度項進行離散在邊i+1/2處，動量方程可統(tǒng)一離散為

式中:θ為隱式因子；Δt為時間步長；上標n表示n時刻；ubt為用ELM求解動量方程對流項后得到的結果。河床阻力項采用ubt和un+1進行離散，以提高計算穩(wěn)定性。

當離散的動量方程中所有顯式項合并后，僅剩下未知項(n+1時刻的壓力ηn+1)，進一步將式(3)轉換為

(4)

(5)

式中Δxi表示第i個單元的長度。

2.3.2 速度-壓力耦合的求解

采用ELM求解得到動量方程對流項和局部線性化連續(xù)性方程對流項的方法，連續(xù)性方程和動量方程的耦合解只需要求解一個線性問題，將式(4)代入式(5)進行速度-壓力耦合，得到

(6)

其中，

式中a、b、c分別表示對應單元i-1、i、i+1的系數(shù)項。

2.3.3 枝狀、環(huán)狀河網(wǎng)的預測-校正法

在實際中，河道常常并不是單一的而是多條縱橫交錯構成河網(wǎng)，呈枝狀、環(huán)狀形態(tài)。此時，采用預測校正法[10-11]實現(xiàn)河網(wǎng)各支汊的耦合求解，在保證計算精度的前提下能有效提高計算速度。河網(wǎng)模擬的預測-校正法以汊點為界將河網(wǎng)分成若干段，分為預測步和校正步進行求解。

2.3.4 測試案例

該試驗中水平管道連接2個水庫，有壓流動發(fā)生在水平管道中，如圖3所示。管道長L=400 m，正方形截面面積a=1 m2。管道中間設有一個關閉的閘門，2個水庫的水位分別為η1和η2且保持不變，水位差η1-η2=1 m。

圖3 水平有壓管段連接2個水庫Fig.3 Two reservoirs connected by a horizontalpressurized pipe

當t=0時，閘門突然打開。假定管道無摩擦，流速和管道內壓力用解析解描述為：

u(x,t)=u0tanh(t/t0) ,

(7)

在該算例中，忽略水的可壓縮性和管壁彈性，即上游到下游的流量相等。在逐漸縮小的均勻網(wǎng)格上對模型進行測試，網(wǎng)格尺度dx依次為40、20、16、10、5 m，建立解的網(wǎng)格獨立性。該模型穩(wěn)定地完成了所有的模擬。隨著網(wǎng)格尺度變小，計算結果逐漸收斂于解析解，如圖4。當網(wǎng)格尺度≤16 m時，計算結果穩(wěn)定趨近于解析解，體現(xiàn)了網(wǎng)格無關性。

圖4 管道內第一個網(wǎng)格點不同網(wǎng)格尺度下速度和壓力的模擬結果Fig.4 Simulated results of velocity and pressure at thefirst grid point inside the pipe (using different grid scales)

3 模型應用

夾巖水利樞紐工程位于貴州省畢節(jié)市與遵義市境內，主要涉及的河流有六沖河以及3條支流：大河、挖嘎河和引底河。大壩位于烏江支流六沖河中游河段上，該段河谷深切，巖溶相當發(fā)育，屬于典型的Karst地貌。擬建夾巖水庫庫區(qū)河流明暗相間(如圖5)，有多處伏流段(大、中、小天橋伏流段)，為典型的明渠伏流交替河道。應用本文建立的明渠與伏流交替河道水動力學模型，研究不同工況條件下六沖河干流沿程的流速變化規(guī)律以及人工隧洞的分流能力。

圖5 夾巖水庫干支流、伏流及人工隧洞布置Fig.5 Arrangement of main stream,tributary,sinkingstream and artificial tunnels in Jiayan Reservoir

3.1 夾巖水庫明渠伏流水動力模型建模與測試

夾巖水庫數(shù)學模型的計算區(qū)域如圖5，包括六沖河(到壩址段)、支流(大河、挖嘎河和引底河)，其中六沖河干流的大、中、小天橋段為伏流段。平面上呈現(xiàn)為枝狀河網(wǎng)形態(tài)，在數(shù)學模型建模時，將計算區(qū)域分為7個河段。采用200 m尺度的計算網(wǎng)格對計算區(qū)域進行剖分，得到241個計算單元。利用散點地形對斷面地形進行插值，在伏流段的33個斷面采用設計單位提供的實測斷面地形資料。

采用實測水文資料和設計單位提供的河道糙率資料對模型進行測試。計算區(qū)域內有七星關水文站，位置見圖5。數(shù)學模型的入流邊界為4個：六沖河、大河、挖嘎河與引底河，出流邊界為一個，即壩址。由于模型的各支流入口均沒有實測入流資料，故采用如下方法設計測試的邊界條件。根據(jù)實測水文資料，將七星關站日流量實測數(shù)據(jù)采用同比放大法得到對應壩址日流量。壩址以上流域集水面積為4 312 km2，大河流域面積為1 195 km2，挖嘎河流域面積為1 195 km2，引底河流域面積為535 km2。采用流域面積比例法得到六沖河與大河流量分配系數(shù)分別為0.60與0.40；根據(jù)流量守恒定理與流域面積比例法得到挖嘎河與引底河流量分配系數(shù)分別為0.38與0.62。

采用2000年水文過程、設計單位提供的糙率進行模型測試。根據(jù)水文資料，河流多處于平水年，2000年為典型平水年，選擇該年具有普適性。河床糙率取值在0.055～0.06，在該糙率條件下，將七星關站的實測值與計算值的流量過程進行套繪，如圖6，實測和計算的流量過程未出現(xiàn)相位偏移，該站的流量相對誤差在1.6%以內。表明本文數(shù)學模型計算方法可行，參系數(shù)取值合理，可用于計算和研究該河段的水流特性。

圖6 七星關站模擬值與計算值的流量過程Fig.6 Flow rate processes at Qixingguan station frommodel test and computation

3.2 干流沿程的流速變化規(guī)律

由于夾巖水庫泥沙問題較為突出，設計部分對庫區(qū)伏流洞是否會被淤堵較為關心。水流流速是泥沙是否淤積的重要判斷指標，因此本節(jié)開展了六沖河干流沿程的流速變化規(guī)律模擬研究。采用恒定水流邊界來研究六沖河干流沿程流速變化規(guī)律，入流邊界選擇枯季流量、多年平均流量以及洪峰流量3種工況，壩前水位均為正常蓄水位1 323 m。第1種工況各開邊界流量分別為4.0、2.7、1.8、3.0 m3/s；第2種工況各開邊界流量分別為22.9、15.1、7.4、12.2 m3/s；第3種工況各開邊界流量分別為425.3、281.7、106.4、174.0 m3/s。3種工況七星關的流量分別為6.7、38.0、707.0 m3/s。

3種工況下干流縱向沿程的流速變化曲線如圖7。流速隨著河道斷面面積的變化而變化，在進入大、中、小天橋伏流段時由于斷面面積驟減而流速增加。由于該河流位于山區(qū)峽谷，過流斷面多變，流速呈現(xiàn)起伏變化。在第1種工況下，干流明渠段與伏流段流速均在0.02 m/s以下；在第2種工況下，干流明渠段流速在0.02 m/s以下，伏流段流速在0.1～0.2 m/s；在第3種工況下，明渠與伏流段流速差異顯著，明渠段流速一般在0.25 m/s以下，伏流段流速一般在1.0～3.0 m/s。

圖7 干流沿程的流速變化Fig.7 Variation of flow velocity along the main stream

3.3人工隧洞的分流能力

夾巖水庫庫尾共有3段伏流河段，分別為大、中、小天橋伏流段。由于伏流段可能存在斷面小過洪能力不足的情況，加之水庫蓄水后可能存在泥沙淤堵風險，產(chǎn)生不利影響。因此，設計單位布置了伏流段分洪隧洞以解決泄洪問題。

伏流段分洪隧洞分布在大中天橋和小天橋河段，每段由2條隧洞并列組成，隧洞斷面均為12 m×16.5 m的城門洞，頂拱中心角120°。分別將2條隧洞合并布置，寬度加倍(24 m)，隧洞的平面布置見圖5。大中天橋隧洞長2 843 m，進口底板高程1 315.980 m，底坡α=0.008，出口高程1 313.706 m；小天橋隧洞長603 m，進口底板高程1 310.732 m，底坡α=0.008，出口高程1 310.250 m。隧洞按無壓流設計，過流能力按明渠均勻流計算。加上人工隧洞后，計算區(qū)域在平面上呈現(xiàn)為環(huán)狀、枝狀河網(wǎng)混合形態(tài)。

圖8 1971—1990年大中天橋、小天橋隧洞和伏流段流量過程Fig.8 Flow processes of artificial tunnels and sinkingstreams in front of big-medium and small overbridgesfrom 1971 to 1990

為闡明人工隧洞的分流能力，開展了河網(wǎng)水動力模擬研究。計算時，河道糙率使用設計單位提供的數(shù)據(jù)0.055～0.06，隧洞糙率取設計值0.014。1971—1990年期間水文資料完整，包含平水年、豐水年、枯水年等各種典型水文年份，能較全面地代表該河段水文條件，故采用1971—1990年作為人工隧洞分流過程計算的系列年。大中天橋、小天橋人工隧洞及大中天橋、小天橋伏流段的流量過程如圖8所示。大中天橋伏流前的人工隧洞的分流比為40.3%～42.8%，小天橋伏流前的人工隧洞的分流比為76.8%～78.6%。從20 a水文過程模擬的統(tǒng)計數(shù)據(jù)來看，大中天橋隧洞分流量為90 622.7 m3/s，占河道總流量的32.7%；小天橋隧洞分流量為197 150.3 m3/s，占河道總流量的70.7%。

4 結 論

(1)基于圣維南方程組，根據(jù)測壓管水頭在明渠、承壓水流中具有相同的表達形式，實現(xiàn)了自由水面水流、承壓水流的統(tǒng)一描述。提出三角形子斷面法描述伏流非規(guī)則封閉斷面,建立了適用于真實Karst地區(qū)明渠伏流交替河道的統(tǒng)一型混合流模型。采用半隱方法離散壓力梯度項，使用歐拉-拉格朗日法求解對流項，并使用預測-校正法實現(xiàn)了枝狀、環(huán)狀河網(wǎng)各支汊的耦合求解。

(2)開展了貴州六沖河夾巖水庫明渠伏流交替河段(包含有大中天橋、小天橋伏流)的數(shù)學模型建模。將計算區(qū)域劃分為7個河段，剖分成241個計算單元。采用2000年水文過程對數(shù)學模型進行了測試，水文站的流量計算誤差在1.6%以內，表明本文數(shù)學模型計算方法可行，參系數(shù)取值合理，可用于研究該河段的水流特性。

(3)采用典型水流條件(枯季、多年平均、洪峰流量)，研究了六沖河干流縱向沿程的流速變化規(guī)律。在洪峰流量條件下，明渠與伏流段流速差異顯著，明渠段流速一般在0.25 m/s以下，伏流段流速一般為1.0～3.0 m/s，有利于伏流過流通道的保持。

(4)采用1971—1990年水文過程，研究了伏流河段前布置人工隧洞的分流規(guī)律。大中天橋伏流前的人工隧洞的分流比為40.3%～42.8%，小天橋伏流前的人工隧洞的分流比為76.8%～78.6%。從20 a水文過程的模擬結果來看，大中天橋隧洞分流量占河道總流量的32.7%；小天橋前人工隧洞分流量占河道總流量的70.7%。由此可推測，小天橋隧洞分流將大幅減小小天橋伏流洞內的流速，加速泥沙淤堵。