吳家陽，李安強，程永光
(1.長江勘測規劃設計研究有限責任公司，湖北 武漢 430010；
2.武漢大學水資源與水電工程科學國家重點實驗室，湖北 武漢 430072)
布置于明渠上的水電站在發生增、甩負荷等瞬變過程時往往會在上游水庫引發較大涌浪，產生機組出力和引航道水位劇烈波動等問題，對水電站的運行安全產生較大危害[1]。采用數值或物理模型手段研究涌浪在明渠內的傳播規律和控制措施，對減緩機組瞬變過程中的安全危害意義重大。目前，明渠瞬變流問題大都采用斷波法[2]或數值求解一維圣維南方程組和二維淺水方程組[3-4]解決。當瞬變過程發生在較為寬闊的河道型水庫(如三峽水庫)或湖泊(如鄱陽湖水利樞紐)內時，為了準確刻畫機組附近復雜的橫向和垂向流動過程，往往需要應用高維模型模擬，反之為了獲得較高的計算效率，則需要盡可能采用低維模型計算。不加區分地采用低維或高維模型模擬樞紐附近和庫尾河道等不同位置的瞬變過程，既不現實也無必要，亟需根據計算域幾何特征和流場特點，有的放矢地選取不同維度的模型模擬不同計算區域的瞬變過程，在不犧牲計算精度的同時保證計算效率。
多維耦合模型在平原河網、入海河灣和大型河庫等區域的水動力學求解中應用較為廣泛。Chen等[5]、王智勇等[6]采用水位預測矯正法，構建了一維-二維耦合的河湖系統整體水動力模型，并采用特征理論證明了該方法的收斂性；
諸裕良等[7]建立了一種一維、二維全隱河網海灣水動力聯網數學模型，采用接口斷面法實現一維、二維模型間的耦合；
Morales-Hernndez等[8]提出了一種適用于明渠淺水流動的一維-二維耦合方案，其中一維、二維模型分別采用隱式有限差分法和有限體積法數值求解；
陳文龍等[9]建立了基于側向聯解的一維-二維耦合水動力學模型，有效克服了基于堰流公式的傳統方法難以處理模型間動量交換的缺點；
顧巍巍等[10]建立了一維-二維耦合的洪水演進數學模型，并采用堰流公式實現一維、二維模型連接處的水流交互；
Liu等[11]提出了一種適用于復雜地形和非規則邊界的一維-二維水動力學耦合模型，并應用于河道-蓄滯洪區系統的水動力求解；
王秀杰等[12]考慮洪水和風暴潮的共同作用，構建了天然河道漫潰堤洪水在防洪保護區的一維-二維水動力耦合模型；
程濤等[13]利用MIKE軟件構建了適用于不規則計算域和地形的潰壩洪水一維-二維耦合模型。可見，目前多維耦合模型主要應用于一維和二維模型間的耦合求解，而三維模型因其模型實現難度大、計算效率較低和自由液面處理復雜等因素在模型耦合中尚不多見。
本文提出一種一維-三維耦合的浸沒邊界-格子玻爾茲曼明渠瞬變流模擬方法。其中，瞬變流場中的狹長河道采用一維圣維南方程組描述，水庫、寬闊湖面、機組進水口、泄水閘等高維流動特征明顯的水域采用三維Navier-Stokes(N-S)方程組描述；
一維圣維南方程組和三維N-S方程組分別采用隱式Pressimann差分方法和格子Boltzmann方法求解，而一維-三維模型間的耦合采用水位預測矯正法[5,7]實現。
1.1 一維明渠圣維南方程組及求解方法
1.1.1 控制方程
在假定明渠底坡較緩、渠道斷面壓力沿垂線按靜水壓力分布、渠道斷面流速均勻分布等條件下，一維明渠水流可由圣維南方程組描述：
(1)
(2)
式中：Z和Q分別為明渠的水位和流量；
A和B分別為明渠斷面的過水面積和寬度；
q為旁側入流或出流的單寬流量；
n為明渠河床糙率系數；
R為水力半徑；
g為重力加速度；
x為沿明渠水流方向的河道里程；
t為時間。式(1)和式(2)適用于一般非棱柱型明渠的水動力學模擬。
1.1.2 數值解法
采用廣泛應用的四點Pressimann隱格式求解一維圣維南方程組，時間偏導數取相鄰節點時間偏導數的平均值，空間偏導數取相鄰節點一階向前差分的加權平均值。控制方程離散后得到的非線性方程組采用Newton-Raphson迭代法求解。Pressimann隱格式的詳細推導過程詳見文獻[12]。
1.2 三維自由液面N-S方程組及求解方法
1.2.1 控制方程
三維黏性不可壓縮流體的控制方程為N-S方程組：
?·u=0
(3)
(4)
式中：u、p、f為分別流速、壓強、外力；
ν為流體的運動黏度；
ρ為流體密度。
1.2.2 數值解法
本文中，N-S方程組由格子Boltzmann方法數值求解。在格子Boltzmann方法中，最基本的物理量為分布函數fα(x,t)，其統計解釋為在t時刻x坐標處沿著α方向運動的粒子質量密度。流體密度和動量可由分布函數的各階矩計算：
(5)
(6)
式中：eα為離散速度方向矢量。
fα(x,t)滿足如下離散格子Boltzmann方程：
(7)
平衡態分布函數一般取為麥克斯韋統計分布，D3Q19離散速度模型的平衡態分布函數可表示為
(8)
式中：ωα為權系數，取值可參考文獻[14]。
1.2.3 三維自由液面處理方法
為了追蹤自由液面，參考文獻[15]中的方法，將三維網格點劃分為流體節點、氣體節點和液面節點3類，其中，流體節點完全被水體填充，氣體節點完全被氣體填充，而液面節點可視為被水體部分填充。液面節點在計算過程中的任一時刻均能形成封閉的曲面，將流體節點和氣體節點完全分隔開，如圖1所示。追蹤自由液面的運動過程分為追蹤液面節點的運動、處理液面節點邊界條件和更新流場節點類別3個步驟。
圖1 三維格子Boltzmann方法中3類網格節點示意Fig.1 Schematic diagram of three types of grid node in three-dimensional lattice Boltzmann method
本文通過計算每一網格點所含水體質量實現自由液面的追蹤。對于每一網格點，定義其質量(m)和體積分數(ε)，體積分數定義為網格點所含水體質量與水體密度之比(ε=m/ρ)。根據該定義，氣體節點ε=0，流體節點ε=1，自由液面節點0<ε<1。與Volume-of-Fluid(VOF)方法類似，液面追蹤通過計算任意網格點與周圍網格點的質量交換實現。考慮到格子Boltzmann方法中分布函數fα(x,t)的物理意義，可以利用分布函數在相鄰網格點間的遷移實現質量交換。基于上述思想，對于t+Δt時刻位于x位置的網格點，相鄰網格點既可能是流體節點，也可能是自由液面節點，當其與相鄰流體節點進行質量交換時，沿x方向的質量變化(Δmα(x，t+Δt))可表示為
(9)
(10)
m(x,t+Δt)=m(x,t)+∑Δmα(x,t+Δt)
(11)
式中：m(x,t+Δt)為自由液面節點下一時步的質量，可根據其計算結果實現節點類別的轉變。若m(x,t+Δt)>(1+κ)ρ(x,t+Δt)，則轉變為流體節點；
若m(x,t+Δt)<-κρ(x,t+Δt)，則轉變為氣體節點。其中，閾值κ=10-3，以防止新生成的自由液面節點在下一時步再次改變節點類型。
需要說明的是，對于新生成的流體節點，其質量可能大于水體密度，導致ε>1；
對于新生成的氣體節點，其質量可能為負值，導致ε<0。因此，需進行質量的二次分配，在保證系統質量守恒前提下，確保所有網格點的體積分數在[0,1]區間內。對于上述2類特殊節點，需進行二次分配的質量(mex)分別為mex=m-ρ和mex=m，為了在二次質量分配后不再出現體積分數超出[0,1]現象，mex將均勻分配給周圍新生成的自由液面節點。
1.3 一維-三維計算方法的耦合
本文將水位預測矯正法[5,7]推廣至一維和三維模型間的耦合。如圖2所示，若某一汊點有M個一維和三維河段與之相連，則在該汊點需滿足如下連接條件：
圖2 一維、三維模型耦合示意Fig.2 Schematic diagram of one-dimensional and three-dimensional model coupling
(12)
Z1=Z2=…=ZM-1=ZM
(13)
式中：Qi為流入和流出汊點的凈流量(流入為正，流出為負)；
Zi為與汊點相連的第i個河段的水位，i=1,2，…，M。
可見，汊點連接條件實質上為各連接河段提供邊界條件。在水位預測矯正方法中，所有模型在汊點處均采用水位邊界條件。若整個明渠系統在n時刻的水位、流量已知，則可假定汊點在n+1時刻的水位為Z′，從而計算整個明渠系統在n+1時刻的水位、流量，進而計算出流入、流出汊點的凈流量。對于非恒定流，汊點凈流量為Z′的函數，若能找到某一迭代方法逐步修正Z′以使得汊點凈流量趨于0，則可實現一維、三維模型間的耦合。根據文獻[10-11]，Z′可按式(14)修正：
(14)
若給定n+1時刻耦合節點的水位初始估計值(Zk，k為迭代步)，一維、三維模型耦合步驟總結如下(δ為某一小量)：
(1) 將明渠流場演化至k迭代步，并計算所有耦合節點的凈流量；
(2) 針對某一耦合汊點，若汊點凈流量小于δ，則停止該耦合節點的迭代，耦合汊點水位邊界條件取為Zk，否則執行步驟(3)；
(3) 根據式(14)計算所有耦合汊點的ΔZ′；
(4) 矯正所有耦合汊點的水位，并更新與耦合汊點相連的所有河段的水位邊界條件；
(5) 將迭代步更新為k=k+1，重復執行步驟(1) 。
采用潮汐波浪運動、局部潰壩波和順直明渠瞬變流3個算例分別驗證一維圣維南方程組求解方案、三維自由液面N-S方程組求解方法和一維-三維耦合的河庫瞬變流模擬方法的有效性和準確性。
2.1 潮汐波浪運動算例
潮汐所引起的波浪運動在海岸工程中十分常見，本小節采用順直河槽中的潮汐波浪運動算例驗證一維圣維南方程組數值解法的有效性和準確性。本例中，順直明渠總長L=14 000 m，渠底高程(H(x))由式(15)給定：
(15)
式中：x為從河槽左側端點起算的明渠里程。計算初始條件和邊界條件分別由式(16)和式(17)給定：
h(x,0)=H(x)u(x,0)=0
(16)
(17)
式中：h(x,t)和u(x,t)分別為明渠水深和斷面平均流速。在上述條件下，明渠內潮汐波浪的波長足夠短，本例存在如下漸進理論解。
(18)
(19)
本例中，計算空間步長Δx=70 m，時間步長Δt=10 s。t=9 120 s時刻明渠沿程水深和流速的計算結果以及漸進理論解如圖3所示。可見，模擬結果與漸進理論解吻合程度較好，沿程水深的最大統計誤差均在1%以內，沿程斷面平均流速的最大統計誤差也在5%以內，從而驗證了一維圣維南方程組四點Pressimann隱式差分解法的有效性和準確性。
圖3 潮汐波浪運動t=9 120 s計算結果與漸進理論解的比較Fig.3 Comparison between the simulated results and the asymptotic theoretical solutions of tidal wave motion at t=9 120 s
2.2 局部潰壩波算例
采用三維局部潰壩算例驗證三維N-S方程組格子Botlzmann數值解法的有效性和準確性，該算例廣泛應用于MacCormack格式、Pressimann隱格式等明渠非恒定流數值解法的驗證。計算條件設置如圖4所示，計算域為200 m×200 m的正方形區域，河床床底高程均為0 m，初始時刻被概化為一條直線的豎直壩體布置在計算域的中心線位置，計算域上游和下游水深分別為10 m和5 m。計算開始時，壩體在圖4所示的位置發生潰壩，潰壩段長度為75 m，潰口一直延伸至河床床底。計算參數如下：重力加速度g=9.8 m/s2，水的運動黏度ν=1.0×10-6m2/s，水體密度ρ=1.0×103kg/m3。
圖4 局部潰壩波算例計算條件設置Fig.4 Simulation setups of partial dam break wave case
圖5為不同時刻水深沿順河向潰壩段中心線的分布，圖6為不同時刻潰壩自由液面的形態。為了與其他數值格式計算結果進行對比，圖5也繪制了液面梯度法[16]和二維淺水格子Boltzmann方法(SWE-LBM)[17]的計算結果。可見，本文方法與其他數值格式的計算結果總體吻合較好，考慮到液面梯度法忽略了湍流效應，并忽略了控制方程中的擴散項，導致波前存在一定程度的不吻合現象。
圖5 潰壩段中心線的沿程水深分布Fig.5 Water depth distribution along the center line of dam break breach
圖6 不同時刻局部潰壩波自由液面形態Fig.6 Free surface configuration of partial dam break wave case at different time
2.3 順直明渠瞬變流算例
為驗證一維-三維耦合的河庫瞬變流模擬方法的有效性和準確性，本節計算簡單順直明渠中的瞬變流過程。計算域為一假想水電站的尾水矩形明渠，明渠長和寬分別為100 m和10 m，渠底底坡為0。尾水渠道下游水深恒定為15 m，上游取為流量邊界條件，且流量過程由水電站機組導葉控制。本例模擬水電站機組增、甩負荷過程中尾水明渠內的涌波過程，并與Delft 3D軟件的模擬結果進行對比。在水電站機組增、甩負荷過程中，假設機組導葉按照線性規律開啟或關閉，開啟或關閉歷時均為8 s，機組額定工況條件下應用流量為600 m3/s。增負荷工況下，初始時刻尾水明渠水體維持靜止，之后導葉開啟，瞬變過程開始。直至尾水明渠內水流達到恒定流狀態，該恒定流狀態也將作為甩負荷工況的計算初始條件。
為了驗證一維-三維耦合的河庫瞬變流模擬方法，本例將機組尾水明渠一分為二，上游50 m渠段采用一維圣維南方程組的Pressimann隱式差分法計算，下游50 m渠段采用三維N-S方程組的格子Boltzmann方法計算。一維、三維模型計算空間步長分別為Δx=0.25 m和Δx=0.20 m，計算時間步長均為t=0.01 s。
圖7為水電站機組增、甩負荷過程中上游水位和下游流量的歷時曲線。可見，本文方法計算結果與Delft 3D軟件模擬結果總體吻合較好，且本文方法捕捉涌波陡峭拐點的能力略強于商用軟件，說明其具有較小的數值擴散性。
圖7 機組尾水明渠上游水位和下游流量的歷時曲線Fig.7 Histories of upstream water level and downstream discharge of tailrace open channel
葛洲壩水利樞紐興建于20世紀七八十年代，是一座低水頭、大流量、徑流式水電站。作為上游三峽水利樞紐工程的反調節航運樞紐，與三峽樞紐工程聯合運行，可消除三峽電站日調節時的下游不穩定流，以改善河道航運條件，并利用該河段水位差發電。
由于泥沙淤積，在葛洲壩樞紐壩址形成了2座島嶼，將壩址江面分為大江、二江和三江3條支汊。樞紐工程共布置21臺機組，其中，大江廠房14臺，額定流量均為825 m3/s;余下7臺機組布置在二江廠房，包括2臺額定流量為1 130 m3/s和5臺額定流量為825 m3/s的機組;三江僅具備沖沙和航運功能。機組以額定工況恒定運行時，庫前水位維持正常蓄水位66 m。
3.1 瞬變流計算條件
本例模擬大江14臺機組同時甩負荷而二江7臺機組以額定流量運行時，上游江面的涌浪過程。大江機組甩負荷時間為33 s，機組應用流量假設按線性規律從滿發流量11 550 m3/s(14×825 m3/s)減小至0 m3/s。
計算范圍為葛洲壩樞紐壩址至三峽樞紐工程共40 km的庫區和河道，其中葛洲壩樞紐壩址上溯5 km范圍采用三維N-S方程組的格子Boltzmann方法模擬，而余下的35 km河道則采用一維圣維南方程組的Pressimann隱式差分法計算。上游邊界設為定流量邊界，下游邊界設為流量邊界條件，且出流過程由21臺機組的應用流量控制。為簡化起見，本例不考慮樞紐沖沙閘和船閘過流能力，且認為機組進水口之間沒有間隙，同時不考慮支流及區間入匯。河道糙率一律取為0.012，計算時間步長取為0.172 8 s，一維、三維模型的計算空間步長分別為11.67 m和8.64 m，網格總數達到千萬量級，恒定流計算耗時14.4 h。
為了對葛洲壩水庫各處的水位進行監測，本文將在圖8所示位置布置流場測點，其中，A測點位于大江沖沙閘前；
B點位于大江機組進水口前；
C點位于二江泄水閘前；
D點位于二江機組進水口前；
E點位于三江沖沙閘前；
F點位于大江防洪淤堤上游端點；
G點位于三江防洪淤堤上游端點；
H點位于支流河口。
圖8 流場監測點布置位置示意Fig.8 Layout of flow field monitoring points
3.2 大江機組甩負荷涌浪模擬結果
各時刻葛洲壩水庫水位和水平流速計算結果如圖9(分圖左為庫水位分布，分圖右為水平流速分布)所示。由于江面起伏相較于計算域水平尺度太小，為了更好地顯示江面起伏情況，圖中水位均沿縱向放大了100倍。可見，大江機組甩負荷后，涌浪隨即以機組為原點以弧形方式向四周傳播。當t=1.0 min時，涌浪首次來到三江防洪淤堤的右岸并發生反射，致使二江機組進水口水位開始上漲。隨后，涌浪向大江防洪淤堤傳播，并在隨后一段時間內在在大江和三江防洪淤堤之間往復來回行進，該現象在模擬初期尤為明顯，在二江左右岸間存在明顯的水位高低交替變化。
圖9 大江機組甩負荷后庫區流場變化示意Fig.9 Flow field in the load rejection process of Dajiang hydro-plant units
涌浪在大江、三江防洪淤堤間往復傳播的同時也在向上游行進。t=2.0 min時，大江沖沙閘首次受到涌波的影響，由于大江河道較為狹窄，其水位上漲非常明顯。t=8.0 min時，大江沖沙閘水位首次到達波谷；
t=3.5 min和t=4.5 min時，大江和二江機組分別達到波谷。
圖10為各流場監測點的水位歷時曲線，其中時間按小時計，且自大江機組甩負荷開始計時。不同頻率涌浪波動之間相互疊加的現象在圖中清晰可見。從圖10中可以得到樞紐各部位受涌浪影響的起始時間、最大和最小涌浪波高、涌浪完全衰減耗時等對樞紐安全運行至關重要的參數。例如，對于三江沖沙閘，最大和最小涌浪分別出現在機組甩負荷后1.2 h和2.55 h左右，最大和最小涌浪水位分別為68.0 m和65.7 m。
圖10 各流場監測點水位過程Fig.10 Water level histories of flow field monitoring points
圖11將大江和二江機組進水口水位進行疊加，高頻振動成分在圖中清晰可見，大江水位的波峰和波谷正好分別對應二江水位的波谷和波峰。如前所述，該現象源自于在大江、三江防洪淤堤間往復傳播的涌浪，因為二江左右岸距離較短，所以該振動以高頻形式附加在低頻基波上，本例中的低頻基波源自于大江機組進水口和上游流量邊界間往復傳播的涌浪。
圖11 大江、二江機組進水口水位—歷時曲線Fig.11 Water level histories of Dajiang and Erjiang hydro-plant unit′s inlets
利用快速傅里葉變換對大江機組、二江機組、大江沖沙閘和三江沖沙閘的水位—歷時曲線進行頻譜分析，區分水位波動的各種頻率成分并對其成因進行解釋。頻譜分析結果表明，大江機組、二江機組、大江沖沙閘均包含4種頻率成分，分別為1.92×10-4Hz(頻率1)、8.32×10-4Hz(頻率2)、2.50×10-3Hz(頻率3)、7.78×10-3Hz(頻率4)；
三江沖沙閘包含2種頻率成分，分別為1.92×10-4Hz(頻率1)、8.32×10-4Hz(頻率2)。本例中河道水深大致在28～31 m范圍內波動，根據波速公式可得涌浪傳播速度約為16 m/s。大江進水口至上游邊界和支流頂端邊界的距離分別為40 km和9 km左右。根據波速公式，往返于兩者間的涌浪周期分別為5 000 s和1 125 s，對應頻率分別為2.00×10-4Hz和8.89×10-4Hz，而這與頻率1和頻率2基本吻合;其次，大江進水口至三江沖沙閘的距離為3.55 km，往返于其間的涌浪周期為443.7 s，對應頻率為2.25×10-3Hz，這與頻率3基本一致；
最后，對于往返于大江、三江防洪淤堤間的涌浪(淤堤間距1.05 km)，其傳播周期為131.3 s，對應頻率為7.62×10-3Hz，這與頻率4基本一致。可見，計算結果基本符合物理規律，說明一維-三維耦合的河庫瞬變流模擬方法基本具備模擬實際明渠瞬變流的能力。
采用水位預測矯正法，提出了一種一維-三維耦合的明渠瞬變流計算模型，其中一維模型采用傳統有限差分法計算，三維模型采用基于自由液面的格子Boltzmann方法模擬，主要結論如下：
(1) 本模型可模擬具有任意地形、邊界條件的實際明渠瞬變流過程，可計算瞬變流過程中的水位、流量和流速等水力變量。
(2) 根據模擬需求，通過在計算域的不同位置運用不同維度的模型進行計算，本模型在不損失整體計算精度的前提下，大幅提高實際明渠瞬變流模擬效率。
(3) 采用潮汐波浪運動、局部潰壩波和順直明渠瞬變流等算例驗證了方法的有效性和正確性，模擬結果與漸進理論解、商用軟件模擬結果均吻合較好。
(4) 將本文方法應用在實際工程中，計算了葛洲壩水利樞紐瞬變過程中庫區的涌浪過程，分析了瞬變過程中的最大、最小涌浪波高和涌浪的衰減規律。通過對涌浪水位過程進行頻譜分析，得到樞紐工程各部位水位波動的頻率組成，并解釋了各頻率成分的形成原因，說明了本方法能準確模擬實際工況下的明渠瞬變流過程。
今后應重點研究本方法的并行加速策略，以大幅提高實際問題的模擬效率。
猜你喜歡變流明渠大江地鐵供電系統中雙向變流裝置設計方案對比研究都市快軌交通(2022年2期)2022-06-28雙向變流裝置在城市軌道交通中的多場景應用研究電力設備管理(2022年7期)2022-05-31百萬雄師過大江湘潮(上半月)(2021年3期)2021-07-20心中的大江音樂教育與創作(2020年2期)2020-06-19農田灌溉明渠水量計量方式分析中華建設(2019年7期)2019-08-27歡迎訂閱《管道系統瞬變流》水利水電工程設計(2019年4期)2019-02-14搞笑秀意林·全彩Color(2018年12期)2019-01-02平移式噴灌機供水明渠底坡設計節水灌溉(2017年2期)2017-03-21大江和堤岸小學閱讀指南·低年級版(2016年10期)2016-09-10北疆第四系覆蓋層上明渠結構優化研究中國水能及電氣化(2016年11期)2016-02-28