馬 斌， 葉星宇， 緱文娟， 王孝群，2，3， 王宇航

(1.天津大學 水利工程仿真與安全國家重點實驗室， 天津 300350； 2.河北工程大學 水利水電學院，河北 邯鄲 056038； 3.河北省智慧水利重點實驗室， 河北 邯鄲 056038；)

1 研究背景

挑流消能具有施工方便，造價及維護成本低的優(yōu)點，在我國高壩工程建設中被廣泛使用，然而下泄水體能量集中而帶來的巨大沖擊力與強烈的水流脈動作用也是此類消能方式的特點[1]，多數工程通過挑坎的體型變化分散水舌、引導水舌、控制水舌的挑射距離等來優(yōu)化防沖消能效果。近年來雙曲面窄縫坎[2]、差動坎[3]、貼角扭面坎[4-5]、舌瓣坎[6]、翻轉坎[7]等已經在國內重大工程上得以運用，有的則作為試驗成果中的推薦體型用以解決工程上的局部動水荷載過大、岸坡沖刷和泄洪霧化等突出問題。對于高水頭大單寬流量的高壩工程，水舌的分散有助于解決下游荷載集中與消能水體利用不充分的問題。根據Deng等[8]的研究，燕尾坎對于水舌的縱向拉伸效果明顯，其團隊與國內學者對這一新型挑坎的水力特性進行了相關的試驗和研究[9-12]，在國內外一些重大工程的消能防沖設計中燕尾坎已經取得初步的實踐效果[13-14]。總體來講燕尾坎的研究和運用歷史尚短，結合大型重力壩壩身泄洪應用的研究則更少，因此對其水力特性的研究仍具有學術價值。經過長期的探索和發(fā)展數值模擬技術已經成為泄洪消能研究的重要手段，燕尾坎出挑水舌流態(tài)復雜，可以借助數值模擬的可視化優(yōu)勢來詳細了解其水力特性。從Xu等[15]的研究結果來看，雷諾時均法的k-ε模型對高壩泄洪中典型淹沒射流的模擬是適用的；郭文思等[16]將k-ε模型模擬的射流場結果與PIV(particle image velocimetry)測量結果進行了比對，進一步驗證了該湍流模型模擬此類淹沒射流的可行性。Yakhot等[17]將重整化群RNG (renormalization group)理論引入湍流領域，在Standardk-ε雙方程模型的基礎上建立了RNGk-ε湍流模型，此后該模型便廣泛應用于求解大彎曲流動問題。從文獻 [18]～[22]的研究成果來看，RNGk-ε湍流模型對于扭面坎、窄縫坎、斜切坎、單孔雙窄縫坎和多股水舌等復雜水舌的泄洪過程均能取得理想的模擬效果。

以上研究成果多以求解的水舌參數、流場、時均壓強等作為檢驗數值模擬精度的指標或作為求解目標，但運用該湍流模型求解水墊塘邊壁脈動壓強的研究很少。雷諾時均法在時間域上對湍流場中的物理量進行了平均，從而失去了部分脈動信息。結合梁在潮等[23]、崔廣濤等[24]、劉昉[25]的研究來看，大渦旋的紊動作用是水墊塘沖擊區(qū)邊壁脈動的主要成因，壓力脈動呈現出幅值大、頻率低的特點，低頻脈動的能量貢獻具有突出優(yōu)勢。雷諾時均法將每一步迭代計算的結果進行時均化處理，但在模擬時間序列中每一步計算結果時仍然具有瞬態(tài)性，若迭代計算的時間頻率遠大于模擬流場中水流的脈動頻率，那么來描述該頻率范圍內低頻的水流脈動理論上是可行的。

本文對連續(xù)坎和燕尾坎的泄洪過程進行了模型試驗和數值模擬，對比了兩種挑坎的水力特性和水墊塘底板動水荷載的分布特性；同時通過對兩種挑坎泄洪過程的數值模擬驗證RNGk-ε湍流模型求解水墊塘沖擊區(qū)低頻脈動的可行性。

2 研究方法

2.1 物理模型試驗

模型試驗依托黃河古賢水利工程，在天津大學高速水流實驗室進行，模型按重力相似準則設計，幾何比尺為1∶60，泄水建筑物與防沖建筑物制作材料為亞克力板，糙率為0.007，如圖1所示。

圖1 試驗模型及燕尾坎消能工體型裝置(單位：m)

泄水建筑物擬采用挑流設計結合重力壩身布置，表孔進口高程為612 m，挑流鼻坎末端高程為480 m。兩種挑坎挑角及反弧半徑相同。下游設消能防沖水墊塘，水墊塘底板高程為450 m，水墊塘長度為300 m，末端二道壩凈高為22 m；以水墊塘底板為基準，校核水位為178.75 m，本文模型試驗與數值模擬均在該水位下進行。

泄水建筑物與下游防沖水墊塘沿表孔中軸線和水墊塘中軸線對稱，物理模型與數值模型建立相同的坐標系，坐標原點設置在挑坎末端，即水墊塘起始點x=0，水墊塘末端x=300 m，數值模型和物理模型關于xoz面對稱。

時均壓強采用測壓排測量，脈動壓強采用YPS300微型數字壓力傳感器和SDA2000數據采集分析系統(tǒng)測量。根據奈奎斯特采樣定理[26]，采樣頻率必須大于脈動頻率的2倍。根據前人研究成果[27-29]可知，水墊塘中的水流脈動頻率主要集中于5 Hz以內，傳感器采樣頻率和數值模擬壓力采樣信號頻率均取100 Hz，采樣時長為180 s。

動水壓強以1 m高度水柱的壓強為計量單位，即9.8 kPa=1 m。時均壓強測點沿水墊塘底板中軸線布置，從67.0～297.5 m處依次布置32個測點，水舌沖擊區(qū)測點間距為6 m，沖擊區(qū)上、下游間距為12 m。脈動壓力傳感器根據不同挑坎體型其位置做出相應調整。數值模型中監(jiān)測點布置于水墊塘底板上表面0.01 m處，中軸線測點與試驗測點保持一致；在水舌沖擊區(qū)y向布置8列脈動壓強測線，每列21個測點，共計168個，如圖2所示。

圖2 數值模型脈動壓強測點布置

2.2 數值模型仿真

2.2.1 數學方程 RNGk-ε湍流模型紊流控制系列方程如下。

連續(xù)方程：

(1)

動量方程：

(2)

k方程：

(3)

ε方程：

(4)

(5)

(6)

式中：ρ為體積分數平均密度，kg/m3；μ為分子黏性系數；g為重力加速度，m/s2；p為修正壓強，Pa；k為湍動能，J；ε為湍動能耗散率；μt為紊流黏性系數；Gk為平均流速梯度變化引起的紊動能項；σk和σε分別為k方程和ε方程的紊流普朗特數；C1ε和C2ε為ε方程常數，ε方程常數值如表1所示。

表1 RNG k-ε控制方程各常數值

2.2.2 計算域及網格劃分處理 模擬范圍涵蓋上游庫區(qū)，其范圍為x向100 m、y向80 m，水庫底部距進口高程100 m，中表孔泄槽、下游水墊塘、下游20 m河床及自由液面以上為水氣流動區(qū)域。整體網格均采用六面體結構化劃分，網格正交性與連續(xù)性良好。網格全局尺寸為1.0 m；下游水墊塘水流行進方向網格尺寸設置為1.0 m，垂直于水流向網格尺寸為0.8 m；流道彎曲段和水舌沖擊區(qū)域以0.5～0.8 m進行局部加密；水舌運動區(qū)域水體與空氣發(fā)生強烈的卷氣和能量交換，其模擬精度關系到水舌的沖擊效果，該部分網格與周邊空氣運動區(qū)獨立劃分，尺寸為0.5 m；水舌運動區(qū)與周邊空氣運動區(qū)采用interface連接，整體網格數量為236×104個，流體域網格劃分如圖3所示。

圖3 數值仿真流體域網格劃分

2.2.3 求解器及邊界條件設置 數值模擬采用壓力基求解器，計算時間為瞬態(tài)，壓力速度耦合求解器選擇PISO (pressure-implicit with splitting of operators)，動量方程離散取二階迎風格式，湍動能及耗散率均為一階迎風格式，收斂條件保持默認設置。

上游庫區(qū)進口設置為壓力進口給出等水位的壓力水頭，流速水頭設置為0；下游出口設置為壓力出口；自由液面、水舌運動區(qū)域和燕尾坎鏤空段設置為壓力進口，設置水的體積分數為0；壁面均設置為無滑移固壁邊界，粗糙度設置為0.5 m，粗糙高度設置為0.001 m。整個模擬過程在標準大氣壓101.32 kPa環(huán)境下進行，空氣密度設置為1.225 kg/m3。

3 水舌形態(tài)及挑坎水流運動特性

3.1 水舌形態(tài)

取數值模擬結果與試驗結果對比兩種挑坎的水舌形態(tài)及射流流態(tài)，進而研究其不同的水舌運動特性和沖擊特性。兩種挑坎的水舌形態(tài)及射流流態(tài)如圖4和5所示，水舌形態(tài)參數試驗值與模擬值匯總見表2。

圖4 兩種挑坎的水舌形態(tài)對比

圖5 數值模擬兩種挑坎射流流態(tài)縱剖面圖(t=60s)

表2 兩種挑坎的水舌形態(tài)參數匯總 m

由圖4和5可看出，連續(xù)坎水舌出挑后有部分橫向擴散，經過最高點后舌心區(qū)周邊水體開始分散，至末端部分水體已經破碎。水墊塘上游0～170 m段水面無明顯波動，強烈的紊動發(fā)生在沖擊區(qū)及下游。燕尾坎泄洪過程中下泄水體行進至燕尾開口處時，水體下邊界失去約束水體先從中間缺口處起挑，另一部分水體從兩側底板起挑。相比于連續(xù)坎，其水舌后緣大幅度向上游移動，前緣挑距也有明顯的減小，水舌縱向獲得大幅度拉伸，側面呈扇形。水舌上緣與下緣之間水體摻氣效果明顯，從水墊塘上游x=70 m處水體開始受水舌后緣水體沖擊，沖擊區(qū)范圍大幅度增加，水墊塘水體的紊動程度減弱且紊動區(qū)范圍擴大。

兩種挑坎水舌橫斷面的沿程變化如圖6所示。由圖6可見，連續(xù)坎水舌出挑后向兩側擴散，中間水體沿z向上凸，沿x向運動趨勢增強，水體主要分布于斷面中上部；燕尾坎水舌出挑后斷面則呈下凹狀，沿x向運動，兩側水體向射流中軸線集中，至x=80 m處斷面開始向z向上凸，相比于連續(xù)坎水舌其橫向擴散幅度減小，水體集中于射流上部，射流下部摻氣明顯。

圖6 兩種挑坎水舌橫斷面沿程變化

3.2 挑坎水流運動特性

為了進一步探究兩種挑坎的水舌擴散特性，對兩種挑坎的水流運動特性進行分析。如圖7所示，截取t=60 s時刻挑坎x=-17 m、x=-11 m、x=-5 m處3個橫斷面，提取各橫斷面上的水氣二相分布、壓強和不同方向流速分量u、v、w的時均分布云圖，u、v、w分別為x、y、z方向的流速分量。

圖7 兩種挑坎水流運動特性研究斷面選取

3.2.1 挑坎斷面水氣二相及壓強分布 圖8、9分別為兩種挑坎各選取斷面的水氣二相分布和壓強分布情況。

圖8 兩種挑坎各斷面水氣二相分布

圖9 兩種挑坎各斷面壓強分布

如圖8、9所示，連續(xù)坎斷面上水體的自由液面靠近擋墻的兩側，略高于流道中間部分，整體呈水平分布，沿x方向沒有明顯的起伏變化；燕尾坎上水體行進至開口處時液面已經有下凹跡象，沿x方向過水斷面呈“漏斗狀”向底板缺口處拉伸。連續(xù)坎上壓強在z方向近似線性分布；燕尾坎壓強兩側大中間小，在y方向形成明顯的梯度變化；兩種挑坎上的斷面壓強分布均沿x向減小。

3.2.2 挑坎斷面流速分布 兩種挑坎各選取斷面流速分布如圖10所示。分析圖10可知，連續(xù)坎上水流流速分布較為均勻，除固壁邊界的粘滯作用影響整體分布均勻，沒有明顯的梯度變化，流速u沿x向運動過程中受挑坎阻力影響流速減?。涣魉賤沿x向運動過程中沒有明顯的梯度變化；由于挑坎上挑方向與流速w保持一致，因此流速w沿x向則有明顯的增大。燕尾坎流速u的核心區(qū)集中于中間缺口處，沿x向運動過程中向下移動；流速v在沿x向運動過程中兩側流速增勢明顯，在y向與中間部分形成較大流速差。

圖10 兩種挑坎各斷面x、y、z向的流速分量分布

如圖11所示，進一步對兩種挑坎末端x=-5 m斷面上的流速u、v、w分量按監(jiān)測面Ⅰ、Ⅱ分塊進行局部監(jiān)測，連續(xù)坎和燕尾坎監(jiān)測面Ⅰ、Ⅱ流速u、v、w均值匯總分別如表3、4所示，其中，u+w為流速u和流速w在平行于平面xoz的縱斷面上的合速度；v+w為流速v和流速w在平行于平面yoz的橫斷面上的合速度；θ11和θ21分別為監(jiān)測面Ⅰ、Ⅱ上平行于xoz的縱斷面上合速度與x軸正方向的夾角；θ12和θ22分別為監(jiān)測面Ⅰ、Ⅱ上平行于yoz的橫斷面上合速度與z軸正方向的夾角。

圖11 兩種挑坎x=-5 m斷面流速監(jiān)測面分塊

表3 兩種挑坎監(jiān)測面Ⅰ流速均值匯總表

表4 兩種挑坎監(jiān)測面Ⅱ流速均值匯總表

由表3、4可知，連續(xù)坎監(jiān)測面Ⅰ、Ⅱ合速度u+w分別為47.64、50.39 m/s，挑角分別為24.71°、26.00°，根據類平拋運動，監(jiān)測面Ⅱ對應的中間部分水體起挑速度和起挑角度更大，成為水舌外緣；監(jiān)測面Ⅰ對應的起挑水體成為水舌內緣，但兩者相差不大，這與連續(xù)坎水舌入水寬度集中的特點相符合。燕尾坎監(jiān)測面Ⅰ、Ⅱ合速度u+w分別為41.52、49.46 m/s，挑角分別為21.81°和9.90°，可見燕尾坎中部起挑水舌與兩側底板起挑水舌存在較大的挑角差，兩側水體獲得更大的挑角成為水舌外緣，中間水體成為水舌內緣；兩者合速度的巨大差異形成燕尾坎水舌縱向的大尺度拉伸。

對比監(jiān)測斷面Ⅰ連續(xù)坎和燕尾坎合速度u+w，兩種挑坎分別為47.64、41.52 m/s，挑角分別為24.71°、21.81°，由此可以解釋等反弧半徑挑角的燕尾坎前緣挑距小于連續(xù)坎。

對比兩挑坎監(jiān)測面Ⅰ、Ⅱ速度v均值的變化，兩者斷面中部的速度v近乎為零，兩側水體的流速燕尾坎和連續(xù)坎分別為2.07和0.18 m/s，可知燕尾坎y向流速梯度遠大于連續(xù)坎，在監(jiān)測面Ⅰ中v+w與z軸正方向的夾角θ12為7.65°，這也解釋了燕尾坎水舌在橫向上的收縮現象。

4 水墊塘底板壓強分布特性

4.1 時均壓強分布特性

兩種挑坎相應的水墊塘底板時均壓強試驗值與模擬值分布如圖12所示。由圖12可看出，連續(xù)坎泄洪時x=50～150 m區(qū)間時均壓強分布較為平緩；x=210～250 m為其沖擊區(qū)，該區(qū)間主流受水墊與底板的折沖作用流速急劇減小，動能轉化為勢能形成沖擊峰，峰值點在x=234.5 m；壁射流區(qū)水體附著于底板向上下游擴散，該部分水體動能增加、勢能減小、時均壓強在峰值點上下游呈下凹狀。燕尾坎時均壓強分布整體呈上升趨勢，x=90～200 m為水舌沖擊區(qū)，時均壓強在該區(qū)間形成兩個連續(xù)的沖擊峰，其局部峰值點在底板x=145.0 m與x=186.5 m，相比于連續(xù)坎時均荷載核心區(qū)向上游移動了70 m，沖擊區(qū)范圍大幅增加，形成的沖擊壓強減小。

連續(xù)坎峰值點試驗值與模擬值分別在x=234.5 m和x=240.5 m，其沖擊壓強值分別為3.09和3.67 m；燕尾坎時均壓強試驗值與模擬值的分布規(guī)律相似，兩者的數值誤差均在5%以內。兩種挑坎時均壓強的沿程分布規(guī)律與試驗結果吻合，其數值誤差均在10%以內。

4.2 脈動壓強均方根分布特性

泄洪過程中水墊塘的水流運動屬于紊流，因此流場中壓強的變化隨時間呈現瞬態(tài)性特征。圖13為連續(xù)坎泄洪過程中水墊塘底板沖擊區(qū)時均壓強峰值點去均值后的試驗實測和模擬壓強時程曲線，由圖13可看出，模擬壓強信號與實測壓強信號分布規(guī)律相似。

圖12 兩種挑坎相應的水墊塘底板時均壓強試驗值與模擬值分布

圖13 連續(xù)坎水墊塘底板沖擊區(qū)峰值點去均值后壓強時程曲線

取連續(xù)坎沖擊區(qū)時均壓強峰值點壓強時程信號進行傅里葉變換，得到的水墊塘沖擊區(qū)典型功率譜如圖14所示。分析圖14可知，試驗值與模擬值的主頻分別為0.016 7和0.013 8 Hz，優(yōu)勢頻率均集中于0～1 Hz；從頻譜的能量分布來看試驗值與模擬值的功率譜能量主要集中于0～2 Hz，由于真實的流場包含更小的渦旋結構， 因此試驗結果在更高頻率段也有一定的能量占比，模擬功率譜大于3 Hz的能量分布幾乎為0，這符合RNGk-ε湍流模型本身的低通濾波屬性；從整體分布趨勢看在0～2Hz能量均隨頻率的增加而遞減，這與前人對水墊塘沖擊區(qū)頻譜特性的研究結論是吻合的。

兩種挑坎相應的水墊塘底板脈動壓強均方根試驗值和模擬值的沿程分布如圖15所示。由圖15可看出，連續(xù)坎在x=50～180 m區(qū)間脈動壓強均方根均在1 m以內；x=180～280 m區(qū)間為波峰段，峰值點在x=234.5 m，峰值為5.92 m；峰值段靠近水墊塘下游，以峰值點為對稱點40 m范圍內下游均方根值大于上游。相比于連續(xù)坎，燕尾坎脈動壓強均方根從起始測點至末端測點整體分布呈波峰狀，在測線上的分布更為均勻，整體降幅明顯，其峰值點在x=207.5 m，相比連續(xù)坎向上游移動了27 m，峰值大小為2.18 m，同比降幅高達57.14%。

圖14 水墊塘底板脈動壓力功率譜密度

圖15 兩種挑坎相應的水墊塘底板脈動壓強均方根值沿程分布

連續(xù)坎試驗值與模擬值的脈動均方根峰值點分別在x=234.5 m和x=240.5 m，數值大小分別為5.92 m和6.84m；燕尾坎試驗值與模擬值峰值點分別在x=207.5 m和x=201.5 m，數值大小分別為2.18和2.41 m。對比模擬結果與實測結果，沖擊區(qū)脈動壓強均方根偏差均在20%以內；相比模型試驗中的脈動壓強均方根峰值降幅為63.18%，數值模擬中相應的降幅為64.77%，可見數值模擬能較好地反映脈動壓強均方根峰值隨挑坎體型的變化。

4.3 水舌沖擊特性

提取連續(xù)坎和燕尾坎兩種挑坎水墊塘沖擊區(qū)監(jiān)測點的脈動壓強均方根模擬值繪制其分布云圖，如圖16所示。由圖16可看出，連續(xù)坎脈動壓強均方根核心區(qū)峰值優(yōu)勢明顯(圖16(a))，從核心區(qū)沿x軸正方向脈動壓強均方根呈衰減趨勢，而在y軸方向擴散半徑沿程增大；沿x軸負方向脈動壓強均方根的傳遞減弱，在y軸方向擴散半徑也沿程減小。相比于連續(xù)坎，燕尾坎脈動壓強均方根整體大幅度減小(圖16(b))，沿核心區(qū)向x軸正方向及x軸負方向衰減速度也更均勻。兩者在y軸方向的衰減速度均大于x軸方向，這也意味著能量的消殺主要集中在射流軸線附近。

圖16 兩種挑坎相應的水墊塘沖擊區(qū)脈動壓強均方根值平面分布

5 結 論

(1)通過數值模擬結果與試驗結果的對比和分析可知，RNGk-ε湍流模型對連續(xù)坎和燕尾坎的泄洪過程模擬效果良好，水舌形態(tài)、水墊塘底板時均壓強分布的模擬結果與實測結果基本吻合，脈動均方根值的對比結果表明該湍流模型求解水墊塘大尺度水流脈動是可行的。

(2)燕尾坎上水流邊界的突變使斷面壓力和流速分布產生了明顯的橫向梯度，梯度指向斷面中軸線；挑坎中間缺口出挑水流流速大、挑角小，形成了水舌內緣，兩側出挑水體流速小、挑角大，沿程向射流軸線集中形成了水舌外緣，水舌橫向收縮明顯而縱向得到大幅度拉伸。

(3)與連續(xù)坎相比，燕尾坎泄洪時水舌沖擊區(qū)范圍及水墊塘水體紊動區(qū)范圍增大，水舌沖擊區(qū)紊動程度減弱；水墊塘底板的時均壓強和脈動壓強核心區(qū)分別向上游移動了70和27 m，沖擊壓強和脈動壓強均方根峰值均大幅減小，荷載能量集中程度減弱。