一維-三維耦合的明渠瞬變流模擬方法與應(yīng)用

吳家陽，李安強(qiáng)，程永光

(1.長江勘測規(guī)劃設(shè)計(jì)研究有限責(zé)任公司，湖北 武漢 430010；2.武漢大學(xué)水資源與水電工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072)

布置于明渠上的水電站在發(fā)生增、甩負(fù)荷等瞬變過程時往往會在上游水庫引發(fā)較大涌浪，產(chǎn)生機(jī)組出力和引航道水位劇烈波動等問題，對水電站的運(yùn)行安全產(chǎn)生較大危害[1]。采用數(shù)值或物理模型手段研究涌浪在明渠內(nèi)的傳播規(guī)律和控制措施，對減緩機(jī)組瞬變過程中的安全危害意義重大。目前，明渠瞬變流問題大都采用斷波法[2]或數(shù)值求解一維圣維南方程組和二維淺水方程組[3-4]解決。當(dāng)瞬變過程發(fā)生在較為寬闊的河道型水庫(如三峽水庫)或湖泊(如鄱陽湖水利樞紐)內(nèi)時，為了準(zhǔn)確刻畫機(jī)組附近復(fù)雜的橫向和垂向流動過程，往往需要應(yīng)用高維模型模擬，反之為了獲得較高的計(jì)算效率，則需要盡可能采用低維模型計(jì)算。不加區(qū)分地采用低維或高維模型模擬樞紐附近和庫尾河道等不同位置的瞬變過程，既不現(xiàn)實(shí)也無必要，亟需根據(jù)計(jì)算域幾何特征和流場特點(diǎn)，有的放矢地選取不同維度的模型模擬不同計(jì)算區(qū)域的瞬變過程，在不犧牲計(jì)算精度的同時保證計(jì)算效率。

多維耦合模型在平原河網(wǎng)、入海河灣和大型河庫等區(qū)域的水動力學(xué)求解中應(yīng)用較為廣泛。Chen等[5]、王智勇等[6]采用水位預(yù)測矯正法，構(gòu)建了一維-二維耦合的河湖系統(tǒng)整體水動力模型，并采用特征理論證明了該方法的收斂性；諸裕良等[7]建立了一種一維、二維全隱河網(wǎng)海灣水動力聯(lián)網(wǎng)數(shù)學(xué)模型，采用接口斷面法實(shí)現(xiàn)一維、二維模型間的耦合；Morales-Hernndez等[8]提出了一種適用于明渠淺水流動的一維-二維耦合方案，其中一維、二維模型分別采用隱式有限差分法和有限體積法數(shù)值求解；陳文龍等[9]建立了基于側(cè)向聯(lián)解的一維-二維耦合水動力學(xué)模型，有效克服了基于堰流公式的傳統(tǒng)方法難以處理模型間動量交換的缺點(diǎn)；顧巍巍等[10]建立了一維-二維耦合的洪水演進(jìn)數(shù)學(xué)模型，并采用堰流公式實(shí)現(xiàn)一維、二維模型連接處的水流交互；Liu等[11]提出了一種適用于復(fù)雜地形和非規(guī)則邊界的一維-二維水動力學(xué)耦合模型，并應(yīng)用于河道-蓄滯洪區(qū)系統(tǒng)的水動力求解；王秀杰等[12]考慮洪水和風(fēng)暴潮的共同作用，構(gòu)建了天然河道漫潰堤洪水在防洪保護(hù)區(qū)的一維-二維水動力耦合模型；程濤等[13]利用MIKE軟件構(gòu)建了適用于不規(guī)則計(jì)算域和地形的潰壩洪水一維-二維耦合模型?？梢?，目前多維耦合模型主要應(yīng)用于一維和二維模型間的耦合求解，而三維模型因其模型實(shí)現(xiàn)難度大、計(jì)算效率較低和自由液面處理復(fù)雜等因素在模型耦合中尚不多見。

本文提出一種一維-三維耦合的浸沒邊界-格子玻爾茲曼明渠瞬變流模擬方法。其中，瞬變流場中的狹長河道采用一維圣維南方程組描述，水庫、寬闊湖面、機(jī)組進(jìn)水口、泄水閘等高維流動特征明顯的水域采用三維Navier-Stokes(N-S)方程組描述；一維圣維南方程組和三維N-S方程組分別采用隱式Pressimann差分方法和格子Boltzmann方法求解，而一維-三維模型間的耦合采用水位預(yù)測矯正法[5,7]實(shí)現(xiàn)。

1 一維-三維耦合的明渠瞬變流模擬方法

1.1 一維明渠圣維南方程組及求解方法

1.1.1 控制方程

在假定明渠底坡較緩、渠道斷面壓力沿垂線按靜水壓力分布、渠道斷面流速均勻分布等條件下，一維明渠水流可由圣維南方程組描述：

(1)

(2)

式中：Z和Q分別為明渠的水位和流量；A和B分別為明渠斷面的過水面積和寬度；q為旁側(cè)入流或出流的單寬流量；n為明渠河床糙率系數(shù)；R為水力半徑；g為重力加速度；x為沿明渠水流方向的河道里程；t為時間。式(1)和式(2)適用于一般非棱柱型明渠的水動力學(xué)模擬。

1.1.2 數(shù)值解法

采用廣泛應(yīng)用的四點(diǎn)Pressimann隱格式求解一維圣維南方程組，時間偏導(dǎo)數(shù)取相鄰節(jié)點(diǎn)時間偏導(dǎo)數(shù)的平均值，空間偏導(dǎo)數(shù)取相鄰節(jié)點(diǎn)一階向前差分的加權(quán)平均值。控制方程離散后得到的非線性方程組采用Newton-Raphson迭代法求解。Pressimann隱格式的詳細(xì)推導(dǎo)過程詳見文獻(xiàn)[12]。

1.2 三維自由液面N-S方程組及求解方法

1.2.1 控制方程

三維黏性不可壓縮流體的控制方程為N-S方程組：

?·u=0

(3)

(4)

式中：u、p、f為分別流速、壓強(qiáng)、外力；ν為流體的運(yùn)動黏度；ρ為流體密度。

1.2.2 數(shù)值解法

本文中，N-S方程組由格子Boltzmann方法數(shù)值求解。在格子Boltzmann方法中，最基本的物理量為分布函數(shù)fα(x,t)，其統(tǒng)計(jì)解釋為在t時刻x坐標(biāo)處沿著α方向運(yùn)動的粒子質(zhì)量密度。流體密度和動量可由分布函數(shù)的各階矩計(jì)算：

(5)

(6)

式中：eα為離散速度方向矢量。

fα(x,t)滿足如下離散格子Boltzmann方程：

(7)

平衡態(tài)分布函數(shù)一般取為麥克斯韋統(tǒng)計(jì)分布，D3Q19離散速度模型的平衡態(tài)分布函數(shù)可表示為

(8)

式中：ωα為權(quán)系數(shù)，取值可參考文獻(xiàn)[14]。

1.2.3 三維自由液面處理方法

為了追蹤自由液面，參考文獻(xiàn)[15]中的方法，將三維網(wǎng)格點(diǎn)劃分為流體節(jié)點(diǎn)、氣體節(jié)點(diǎn)和液面節(jié)點(diǎn)3類，其中，流體節(jié)點(diǎn)完全被水體填充，氣體節(jié)點(diǎn)完全被氣體填充，而液面節(jié)點(diǎn)可視為被水體部分填充。液面節(jié)點(diǎn)在計(jì)算過程中的任一時刻均能形成封閉的曲面，將流體節(jié)點(diǎn)和氣體節(jié)點(diǎn)完全分隔開，如圖1所示。追蹤自由液面的運(yùn)動過程分為追蹤液面節(jié)點(diǎn)的運(yùn)動、處理液面節(jié)點(diǎn)邊界條件和更新流場節(jié)點(diǎn)類別3個步驟。

圖1 三維格子Boltzmann方法中3類網(wǎng)格節(jié)點(diǎn)示意Fig.1 Schematic diagram of three types of grid node in three-dimensional lattice Boltzmann method

本文通過計(jì)算每一網(wǎng)格點(diǎn)所含水體質(zhì)量實(shí)現(xiàn)自由液面的追蹤。對于每一網(wǎng)格點(diǎn)，定義其質(zhì)量(m)和體積分?jǐn)?shù)(ε)，體積分?jǐn)?shù)定義為網(wǎng)格點(diǎn)所含水體質(zhì)量與水體密度之比(ε=m/ρ)。根據(jù)該定義，氣體節(jié)點(diǎn)ε=0，流體節(jié)點(diǎn)ε=1，自由液面節(jié)點(diǎn)0<ε<1。與Volume-of-Fluid(VOF)方法類似，液面追蹤通過計(jì)算任意網(wǎng)格點(diǎn)與周圍網(wǎng)格點(diǎn)的質(zhì)量交換實(shí)現(xiàn)?？紤]到格子Boltzmann方法中分布函數(shù)fα(x,t)的物理意義，可以利用分布函數(shù)在相鄰網(wǎng)格點(diǎn)間的遷移實(shí)現(xiàn)質(zhì)量交換?；谏鲜鏊枷耄瑢τ趖+Δt時刻位于x位置的網(wǎng)格點(diǎn)，相鄰網(wǎng)格點(diǎn)既可能是流體節(jié)點(diǎn)，也可能是自由液面節(jié)點(diǎn)，當(dāng)其與相鄰流體節(jié)點(diǎn)進(jìn)行質(zhì)量交換時，沿x方向的質(zhì)量變化(Δmα(x，t+Δt))可表示為

(9)

(10)

m(x,t+Δt)=m(x,t)+∑Δmα(x,t+Δt)

(11)

式中：m(x,t+Δt)為自由液面節(jié)點(diǎn)下一時步的質(zhì)量，可根據(jù)其計(jì)算結(jié)果實(shí)現(xiàn)節(jié)點(diǎn)類別的轉(zhuǎn)變。若m(x,t+Δt)>(1+κ)ρ(x,t+Δt)，則轉(zhuǎn)變?yōu)榱黧w節(jié)點(diǎn)；若m(x,t+Δt)<-κρ(x,t+Δt)，則轉(zhuǎn)變?yōu)闅怏w節(jié)點(diǎn)。其中，閾值κ=10-3，以防止新生成的自由液面節(jié)點(diǎn)在下一時步再次改變節(jié)點(diǎn)類型。

需要說明的是，對于新生成的流體節(jié)點(diǎn)，其質(zhì)量可能大于水體密度，導(dǎo)致ε>1；對于新生成的氣體節(jié)點(diǎn)，其質(zhì)量可能為負(fù)值，導(dǎo)致ε<0。因此，需進(jìn)行質(zhì)量的二次分配，在保證系統(tǒng)質(zhì)量守恒前提下，確保所有網(wǎng)格點(diǎn)的體積分?jǐn)?shù)在[0,1]區(qū)間內(nèi)。對于上述2類特殊節(jié)點(diǎn)，需進(jìn)行二次分配的質(zhì)量(mex)分別為mex=m-ρ和mex=m，為了在二次質(zhì)量分配后不再出現(xiàn)體積分?jǐn)?shù)超出[0,1]現(xiàn)象，mex將均勻分配給周圍新生成的自由液面節(jié)點(diǎn)。

1.3 一維-三維計(jì)算方法的耦合

本文將水位預(yù)測矯正法[5,7]推廣至一維和三維模型間的耦合。如圖2所示，若某一汊點(diǎn)有M個一維和三維河段與之相連，則在該汊點(diǎn)需滿足如下連接條件：

圖2 一維、三維模型耦合示意Fig.2 Schematic diagram of one-dimensional and three-dimensional model coupling

(12)

Z1=Z2=…=ZM-1=ZM

(13)

式中：Qi為流入和流出汊點(diǎn)的凈流量(流入為正，流出為負(fù))；Zi為與汊點(diǎn)相連的第i個河段的水位，i=1,2，…，M。

可見，汊點(diǎn)連接條件實(shí)質(zhì)上為各連接河段提供邊界條件。在水位預(yù)測矯正方法中，所有模型在汊點(diǎn)處均采用水位邊界條件。若整個明渠系統(tǒng)在n時刻的水位、流量已知，則可假定汊點(diǎn)在n+1時刻的水位為Z′，從而計(jì)算整個明渠系統(tǒng)在n+1時刻的水位、流量，進(jìn)而計(jì)算出流入、流出汊點(diǎn)的凈流量。對于非恒定流，汊點(diǎn)凈流量為Z′的函數(shù)，若能找到某一迭代方法逐步修正Z′以使得汊點(diǎn)凈流量趨于0，則可實(shí)現(xiàn)一維、三維模型間的耦合。根據(jù)文獻(xiàn)[10-11]，Z′可按式(14)修正：

(14)

若給定n+1時刻耦合節(jié)點(diǎn)的水位初始估計(jì)值(Zk，k為迭代步)，一維、三維模型耦合步驟總結(jié)如下(δ為某一小量)：

(1) 將明渠流場演化至k迭代步，并計(jì)算所有耦合節(jié)點(diǎn)的凈流量；

(2) 針對某一耦合汊點(diǎn)，若汊點(diǎn)凈流量小于δ，則停止該耦合節(jié)點(diǎn)的迭代，耦合汊點(diǎn)水位邊界條件取為Zk，否則執(zhí)行步驟(3)；

(3) 根據(jù)式(14)計(jì)算所有耦合汊點(diǎn)的ΔZ′；

(4) 矯正所有耦合汊點(diǎn)的水位，并更新與耦合汊點(diǎn)相連的所有河段的水位邊界條件；

(5) 將迭代步更新為k=k+1，重復(fù)執(zhí)行步驟(1) 。

2 方法驗(yàn)證

采用潮汐波浪運(yùn)動、局部潰壩波和順直明渠瞬變流3個算例分別驗(yàn)證一維圣維南方程組求解方案、三維自由液面N-S方程組求解方法和一維-三維耦合的河庫瞬變流模擬方法的有效性和準(zhǔn)確性。

2.1 潮汐波浪運(yùn)動算例

潮汐所引起的波浪運(yùn)動在海岸工程中十分常見，本小節(jié)采用順直河槽中的潮汐波浪運(yùn)動算例驗(yàn)證一維圣維南方程組數(shù)值解法的有效性和準(zhǔn)確性。本例中，順直明渠總長L=14 000 m，渠底高程(H(x))由式(15)給定：

(15)

式中：x為從河槽左側(cè)端點(diǎn)起算的明渠里程。計(jì)算初始條件和邊界條件分別由式(16)和式(17)給定：

h(x,0)=H(x)u(x,0)=0

(16)

(17)

式中：h(x,t)和u(x,t)分別為明渠水深和斷面平均流速。在上述條件下，明渠內(nèi)潮汐波浪的波長足夠短，本例存在如下漸進(jìn)理論解。

(18)

(19)

本例中，計(jì)算空間步長Δx=70 m，時間步長Δt=10 s。t=9 120 s時刻明渠沿程水深和流速的計(jì)算結(jié)果以及漸進(jìn)理論解如圖3所示?？梢?，模擬結(jié)果與漸進(jìn)理論解吻合程度較好，沿程水深的最大統(tǒng)計(jì)誤差均在1%以內(nèi)，沿程斷面平均流速的最大統(tǒng)計(jì)誤差也在5%以內(nèi)，從而驗(yàn)證了一維圣維南方程組四點(diǎn)Pressimann隱式差分解法的有效性和準(zhǔn)確性。

圖3 潮汐波浪運(yùn)動t=9 120 s計(jì)算結(jié)果與漸進(jìn)理論解的比較Fig.3 Comparison between the simulated results and the asymptotic theoretical solutions of tidal wave motion at t=9 120 s

2.2 局部潰壩波算例

采用三維局部潰壩算例驗(yàn)證三維N-S方程組格子Botlzmann數(shù)值解法的有效性和準(zhǔn)確性，該算例廣泛應(yīng)用于MacCormack格式、Pressimann隱格式等明渠非恒定流數(shù)值解法的驗(yàn)證。計(jì)算條件設(shè)置如圖4所示，計(jì)算域?yàn)?00 m×200 m的正方形區(qū)域，河床床底高程均為0 m，初始時刻被概化為一條直線的豎直壩體布置在計(jì)算域的中心線位置，計(jì)算域上游和下游水深分別為10 m和5 m。計(jì)算開始時，壩體在圖4所示的位置發(fā)生潰壩，潰壩段長度為75 m，潰口一直延伸至河床床底。計(jì)算參數(shù)如下：重力加速度g=9.8 m/s2，水的運(yùn)動黏度ν=1.0×10-6m2/s，水體密度ρ=1.0×103kg/m3。

圖4 局部潰壩波算例計(jì)算條件設(shè)置Fig.4 Simulation setups of partial dam break wave case

圖5為不同時刻水深沿順河向潰壩段中心線的分布，圖6為不同時刻潰壩自由液面的形態(tài)。為了與其他數(shù)值格式計(jì)算結(jié)果進(jìn)行對比，圖5也繪制了液面梯度法[16]和二維淺水格子Boltzmann方法(SWE-LBM)[17]的計(jì)算結(jié)果?？梢?，本文方法與其他數(shù)值格式的計(jì)算結(jié)果總體吻合較好，考慮到液面梯度法忽略了湍流效應(yīng)，并忽略了控制方程中的擴(kuò)散項(xiàng)，導(dǎo)致波前存在一定程度的不吻合現(xiàn)象。

圖5 潰壩段中心線的沿程水深分布Fig.5 Water depth distribution along the center line of dam break breach

圖6 不同時刻局部潰壩波自由液面形態(tài)Fig.6 Free surface configuration of partial dam break wave case at different time

2.3 順直明渠瞬變流算例

為驗(yàn)證一維-三維耦合的河庫瞬變流模擬方法的有效性和準(zhǔn)確性，本節(jié)計(jì)算簡單順直明渠中的瞬變流過程。計(jì)算域?yàn)橐患傧胨娬镜奈菜匦蚊髑?，明渠長和寬分別為100 m和10 m，渠底底坡為0。尾水渠道下游水深恒定為15 m，上游取為流量邊界條件，且流量過程由水電站機(jī)組導(dǎo)葉控制。本例模擬水電站機(jī)組增、甩負(fù)荷過程中尾水明渠內(nèi)的涌波過程，并與Delft 3D軟件的模擬結(jié)果進(jìn)行對比。在水電站機(jī)組增、甩負(fù)荷過程中，假設(shè)機(jī)組導(dǎo)葉按照線性規(guī)律開啟或關(guān)閉，開啟或關(guān)閉歷時均為8 s，機(jī)組額定工況條件下應(yīng)用流量為600 m3/s。增負(fù)荷工況下，初始時刻尾水明渠水體維持靜止，之后導(dǎo)葉開啟，瞬變過程開始。直至尾水明渠內(nèi)水流達(dá)到恒定流狀態(tài)，該恒定流狀態(tài)也將作為甩負(fù)荷工況的計(jì)算初始條件。

為了驗(yàn)證一維-三維耦合的河庫瞬變流模擬方法，本例將機(jī)組尾水明渠一分為二，上游50 m渠段采用一維圣維南方程組的Pressimann隱式差分法計(jì)算，下游50 m渠段采用三維N-S方程組的格子Boltzmann方法計(jì)算。一維、三維模型計(jì)算空間步長分別為Δx=0.25 m和Δx=0.20 m，計(jì)算時間步長均為t=0.01 s。

圖7為水電站機(jī)組增、甩負(fù)荷過程中上游水位和下游流量的歷時曲線?？梢?，本文方法計(jì)算結(jié)果與Delft 3D軟件模擬結(jié)果總體吻合較好，且本文方法捕捉涌波陡峭拐點(diǎn)的能力略強(qiáng)于商用軟件，說明其具有較小的數(shù)值擴(kuò)散性。

圖7 機(jī)組尾水明渠上游水位和下游流量的歷時曲線Fig.7 Histories of upstream water level and downstream discharge of tailrace open channel

3 葛洲壩水庫涌浪特性模擬

葛洲壩水利樞紐興建于20世紀(jì)七八十年代，是一座低水頭、大流量、徑流式水電站。作為上游三峽水利樞紐工程的反調(diào)節(jié)航運(yùn)樞紐，與三峽樞紐工程聯(lián)合運(yùn)行，可消除三峽電站日調(diào)節(jié)時的下游不穩(wěn)定流，以改善河道航運(yùn)條件，并利用該河段水位差發(fā)電。

由于泥沙淤積，在葛洲壩樞紐壩址形成了2座島嶼，將壩址江面分為大江、二江和三江3條支汊。樞紐工程共布置21臺機(jī)組，其中，大江廠房14臺，額定流量均為825 m3/s;余下7臺機(jī)組布置在二江廠房，包括2臺額定流量為1 130 m3/s和5臺額定流量為825 m3/s的機(jī)組;三江僅具備沖沙和航運(yùn)功能。機(jī)組以額定工況恒定運(yùn)行時，庫前水位維持正常蓄水位66 m。

3.1 瞬變流計(jì)算條件

本例模擬大江14臺機(jī)組同時甩負(fù)荷而二江7臺機(jī)組以額定流量運(yùn)行時，上游江面的涌浪過程。大江機(jī)組甩負(fù)荷時間為33 s，機(jī)組應(yīng)用流量假設(shè)按線性規(guī)律從滿發(fā)流量11 550 m3/s(14×825 m3/s)減小至0 m3/s。

計(jì)算范圍為葛洲壩樞紐壩址至三峽樞紐工程共40 km的庫區(qū)和河道，其中葛洲壩樞紐壩址上溯5 km范圍采用三維N-S方程組的格子Boltzmann方法模擬，而余下的35 km河道則采用一維圣維南方程組的Pressimann隱式差分法計(jì)算。上游邊界設(shè)為定流量邊界，下游邊界設(shè)為流量邊界條件，且出流過程由21臺機(jī)組的應(yīng)用流量控制。為簡化起見，本例不考慮樞紐沖沙閘和船閘過流能力，且認(rèn)為機(jī)組進(jìn)水口之間沒有間隙，同時不考慮支流及區(qū)間入?yún)R。河道糙率一律取為0.012，計(jì)算時間步長取為0.172 8 s，一維、三維模型的計(jì)算空間步長分別為11.67 m和8.64 m，網(wǎng)格總數(shù)達(dá)到千萬量級，恒定流計(jì)算耗時14.4 h。

為了對葛洲壩水庫各處的水位進(jìn)行監(jiān)測，本文將在圖8所示位置布置流場測點(diǎn)，其中，A測點(diǎn)位于大江沖沙閘前；B點(diǎn)位于大江機(jī)組進(jìn)水口前；C點(diǎn)位于二江泄水閘前；D點(diǎn)位于二江機(jī)組進(jìn)水口前；E點(diǎn)位于三江沖沙閘前；F點(diǎn)位于大江防洪淤堤上游端點(diǎn)；G點(diǎn)位于三江防洪淤堤上游端點(diǎn)；H點(diǎn)位于支流河口。

圖8 流場監(jiān)測點(diǎn)布置位置示意Fig.8 Layout of flow field monitoring points

3.2 大江機(jī)組甩負(fù)荷涌浪模擬結(jié)果

各時刻葛洲壩水庫水位和水平流速計(jì)算結(jié)果如圖9(分圖左為庫水位分布，分圖右為水平流速分布)所示。由于江面起伏相較于計(jì)算域水平尺度太小，為了更好地顯示江面起伏情況，圖中水位均沿縱向放大了100倍。可見，大江機(jī)組甩負(fù)荷后，涌浪隨即以機(jī)組為原點(diǎn)以弧形方式向四周傳播。當(dāng)t=1.0 min時，涌浪首次來到三江防洪淤堤的右岸并發(fā)生反射，致使二江機(jī)組進(jìn)水口水位開始上漲。隨后，涌浪向大江防洪淤堤傳播，并在隨后一段時間內(nèi)在在大江和三江防洪淤堤之間往復(fù)來回行進(jìn)，該現(xiàn)象在模擬初期尤為明顯，在二江左右岸間存在明顯的水位高低交替變化。

圖9 大江機(jī)組甩負(fù)荷后庫區(qū)流場變化示意Fig.9 Flow field in the load rejection process of Dajiang hydro-plant units

涌浪在大江、三江防洪淤堤間往復(fù)傳播的同時也在向上游行進(jìn)。t=2.0 min時，大江沖沙閘首次受到涌波的影響，由于大江河道較為狹窄，其水位上漲非常明顯。t=8.0 min時，大江沖沙閘水位首次到達(dá)波谷；t=3.5 min和t=4.5 min時，大江和二江機(jī)組分別達(dá)到波谷。

圖10為各流場監(jiān)測點(diǎn)的水位歷時曲線，其中時間按小時計(jì)，且自大江機(jī)組甩負(fù)荷開始計(jì)時。不同頻率涌浪波動之間相互疊加的現(xiàn)象在圖中清晰可見。從圖10中可以得到樞紐各部位受涌浪影響的起始時間、最大和最小涌浪波高、涌浪完全衰減耗時等對樞紐安全運(yùn)行至關(guān)重要的參數(shù)。例如，對于三江沖沙閘，最大和最小涌浪分別出現(xiàn)在機(jī)組甩負(fù)荷后1.2 h和2.55 h左右，最大和最小涌浪水位分別為68.0 m和65.7 m。

圖10 各流場監(jiān)測點(diǎn)水位過程Fig.10 Water level histories of flow field monitoring points

圖11將大江和二江機(jī)組進(jìn)水口水位進(jìn)行疊加，高頻振動成分在圖中清晰可見，大江水位的波峰和波谷正好分別對應(yīng)二江水位的波谷和波峰。如前所述，該現(xiàn)象源自于在大江、三江防洪淤堤間往復(fù)傳播的涌浪，因?yàn)槎笥野毒嚯x較短，所以該振動以高頻形式附加在低頻基波上，本例中的低頻基波源自于大江機(jī)組進(jìn)水口和上游流量邊界間往復(fù)傳播的涌浪。

圖11 大江、二江機(jī)組進(jìn)水口水位—?dú)v時曲線Fig.11 Water level histories of Dajiang and Erjiang hydro-plant unit′s inlets

利用快速傅里葉變換對大江機(jī)組、二江機(jī)組、大江沖沙閘和三江沖沙閘的水位—?dú)v時曲線進(jìn)行頻譜分析，區(qū)分水位波動的各種頻率成分并對其成因進(jìn)行解釋。頻譜分析結(jié)果表明，大江機(jī)組、二江機(jī)組、大江沖沙閘均包含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范圍內(nèi)波動，根據(jù)波速公式可得涌浪傳播速度約為16 m/s。大江進(jìn)水口至上游邊界和支流頂端邊界的距離分別為40 km和9 km左右。根據(jù)波速公式，往返于兩者間的涌浪周期分別為5 000 s和1 125 s，對應(yīng)頻率分別為2.00×10-4Hz和8.89×10-4Hz，而這與頻率1和頻率2基本吻合;其次，大江進(jìn)水口至三江沖沙閘的距離為3.55 km，往返于其間的涌浪周期為443.7 s，對應(yīng)頻率為2.25×10-3Hz，這與頻率3基本一致；最后，對于往返于大江、三江防洪淤堤間的涌浪(淤堤間距1.05 km)，其傳播周期為131.3 s，對應(yīng)頻率為7.62×10-3Hz，這與頻率4基本一致。可見，計(jì)算結(jié)果基本符合物理規(guī)律，說明一維-三維耦合的河庫瞬變流模擬方法基本具備模擬實(shí)際明渠瞬變流的能力。

4 結(jié) 論

采用水位預(yù)測矯正法，提出了一種一維-三維耦合的明渠瞬變流計(jì)算模型，其中一維模型采用傳統(tǒng)有限差分法計(jì)算，三維模型采用基于自由液面的格子Boltzmann方法模擬，主要結(jié)論如下：

(1) 本模型可模擬具有任意地形、邊界條件的實(shí)際明渠瞬變流過程，可計(jì)算瞬變流過程中的水位、流量和流速等水力變量。

(2) 根據(jù)模擬需求，通過在計(jì)算域的不同位置運(yùn)用不同維度的模型進(jìn)行計(jì)算，本模型在不損失整體計(jì)算精度的前提下，大幅提高實(shí)際明渠瞬變流模擬效率。

(3) 采用潮汐波浪運(yùn)動、局部潰壩波和順直明渠瞬變流等算例驗(yàn)證了方法的有效性和正確性，模擬結(jié)果與漸進(jìn)理論解、商用軟件模擬結(jié)果均吻合較好。

(4) 將本文方法應(yīng)用在實(shí)際工程中，計(jì)算了葛洲壩水利樞紐瞬變過程中庫區(qū)的涌浪過程，分析了瞬變過程中的最大、最小涌浪波高和涌浪的衰減規(guī)律。通過對涌浪水位過程進(jìn)行頻譜分析，得到樞紐工程各部位水位波動的頻率組成，并解釋了各頻率成分的形成原因，說明了本方法能準(zhǔn)確模擬實(shí)際工況下的明渠瞬變流過程。

今后應(yīng)重點(diǎn)研究本方法的并行加速策略，以大幅提高實(shí)際問題的模擬效率。