薄夫萍，吳海濤，王 星，張云云，顧聲龍

(青海大學水利電力學院，青海 西寧 810016)

水躍是底流消能工中常用的消能形式。經(jīng)過水躍這一局部水力現(xiàn)象，水流流態(tài)發(fā)生急劇變化，水躍區(qū)段內(nèi)底部主流區(qū)和上部旋滾區(qū)相互摻混、碰撞、摩擦，下泄水流消耗大部分動能從而使流速銳減、流態(tài)趨于平穩(wěn)，以達到消能效果。但當消能率較低時，下泄水流經(jīng)消力池后仍攜帶大量余能，致使下游河床遭受嚴重沖刷甚至造成壩體損壞[1]。因此，提高消力池的消能率是目前水工結(jié)構(gòu)設(shè)計中迫切需要解決的問題。

正弦形消力池底板最初由Ead等[2]提出，其本質(zhì)是將光滑的消力池底板改為正弦形以增加壁面粗糙度，使水流之間的紊動更為劇烈，水躍特性急劇改變，從而增加消能率。目前對此種消力池的研究方法主要有試驗研究、理論分析、數(shù)值模擬。Hughes等[3]通過對消力池底板加設(shè)不同礫石的試驗研究表明，隨著底板粗糙度的增加，水躍躍長明顯縮短，消能率有所提高；Tokyay[4]、Izadjoo等[5]、Abbaspour等[6]、Elsebaie等[7]、Samadi等[8]分別通過大量試驗研究了正弦形、梯形、三角形底板上的水躍，指明共軛水深是與弗勞德數(shù)(Fr)有關(guān)的函數(shù)，且正弦形底板在消能方面具有很大潛力。傅銘煥等[9]、張志昌等[10-11]結(jié)合斷面流速分布對數(shù)規(guī)律及動量方程給出了共軛水深、水躍長度等理論計算方法，并指出壁面阻力對波狀底板水躍特性有很大的影響，是解決實際問題中不可忽略的重要因素。由于試驗研究受到儀器、測量方法等因素的限制，很難準確測量復(fù)雜水力現(xiàn)象中急劇變化的參數(shù)；而理論分析需對研究問題做大量假設(shè)以簡化控制方程便于求解，但與實際現(xiàn)象偏離且求解過程繁瑣，且參數(shù)選取主要基于經(jīng)驗公式。數(shù)值模擬則可以彌補試驗研究及理論分析的不足，記錄現(xiàn)象演變的全過程，為所研究問題提供完整、具體的數(shù)據(jù)，便于細部結(jié)構(gòu)的觀察，有利于揭示現(xiàn)象變化的本質(zhì)。目前研究正弦形底板水躍現(xiàn)象采用的數(shù)值模擬方法主要為網(wǎng)格方法：魏文禮等[12-13]、程香菊等[14]、Abbaspour等[15]采用網(wǎng)格法建立不同形式的湍流模型，對波狀底板上的水躍進行數(shù)值模擬，得到水躍的演變過程及近壁面處存在的小渦分布。但網(wǎng)格法由于受到網(wǎng)格節(jié)點的限制，在模擬自由表面、運動交界面、大變形等引起網(wǎng)格嚴重變形的多尺度問題時精度和準確性差強人意。而水躍現(xiàn)象中水流紊動劇烈并伴隨空氣不斷地摻入[16]，自由表面存在旋滾，因而水面波動大且水體破碎嚴重，因此考慮采用新的數(shù)值模擬方法即光滑質(zhì)點水動力學法(SPH方法)進行研究。SPH方法是純拉格朗日形式的無網(wǎng)格粒子法之一，對需要計算的問題域采用任意分布的粒子框架進行描述，在處理大變形等問題上具有獨特優(yōu)勢，目前已廣泛地應(yīng)用于滑坡涌浪[17]、多相流[18]、潰壩[19]、溢洪道[20]等工程流體動力學問題的研究，并取得了重要的成果。

結(jié)合水躍的特點及SPH方法的優(yōu)勢，筆者將使用SPHYSICS_2D軟件對正弦形消力池底板上水躍現(xiàn)象建模計算，以分析其水面線、流速分布、躍長、共軛水深、消能率等水躍特性相較于光滑底板的變化趨勢。

1 經(jīng)SPH方法離散后的N-S方程

流體動力學中通常使用邊界條件、初始狀態(tài)結(jié)合質(zhì)量、動量、能量三大守恒定律建立控制方程，以確定流體系統(tǒng)的運動情況。對無窮小流體單元采用拉格朗日描述法建立控制方程可得偏微分形式下的連續(xù)性方程、動量方程、能量方程，即N-S方程：

(1)

(2)

(3)

式中：α、β表示坐標方向；ρ為流體密度；t為時間；V為流體速度向量；X為流體位置向量；g為重力加速度；e為流體能量；σαβ為總應(yīng)力張量；P各向同性壓力；θ為黏性項。

SPH方法中將任意函數(shù)改寫為SPH形式需要通過核近似、粒子近似2個關(guān)鍵步驟。經(jīng)核近似后得場函數(shù)積分表示(式(4))，經(jīng)粒子近似后得函數(shù)的離散形式(式(5))：

(4)

(5)

其中Wij=W(ri-rj,h)

式中：A(r)為任意函數(shù)；W(r-r′,h)為光滑函數(shù)；h為光滑長度，用以控制計算域;r和r′代表粒子位置；N為粒子總數(shù)；mj和ρj分別為j粒子質(zhì)量和密度；Wij為i、j粒子間的光滑函數(shù)。

SPH方法中采用亞粒子模型(SPS)建立湍流模型，該模型由大渦模型(LES)和亞網(wǎng)格尺寸模型(SGS)組合而成。對黏性項θ的處理采用SPS方法，可得動量方程(式(6))；在SPH方法中流體被認為是弱可壓縮性的，因此為控制流體的壓縮性引入了狀態(tài)方程(7)，通過改變聲速調(diào)節(jié)數(shù)值模型時間步長以滿足壓縮性調(diào)節(jié)要求。但是聲速也不能過小，至少是流體最大流速的10倍，這樣才能確保粒子密度變化率小于1%，控制流體為弱可壓縮狀態(tài)。

(6)

(7)

式中：ρ0為預(yù)設(shè)密度，取值為1 000 kg/m3；c0為預(yù)設(shè)密度對應(yīng)的聲速;?0為層流動力黏度，取10-6m2/s；τ為黏性應(yīng)力。

綜上所述，經(jīng)SPH方法離散后的N-S方程分別為以下形式：

(8)

(9)

(10)

式中：下標a、b表示計算域內(nèi)任意2個粒子；mb為b粒子的質(zhì)量；Vab為a、b粒子之間的相對速度向量；Wab為a、b粒子間核函數(shù)梯度；ψab為經(jīng)SPH方法離散后的黏性項。

水躍消能率的計算通常使用躍前斷面能量E1與躍后斷面能量E2之差作為能量損失值ΔE，因此可以得到躍后能量損失為

(11)

式中：Z1、Z2分別為躍前斷面、躍后斷面水位高度；α1、α2分別為躍前斷面、躍后斷面的動能修正系數(shù);U1、U2分別為躍前斷面、躍后斷面平均流速；P1、P2分別為躍前斷面、躍后斷面壓強。

當渠底水平時，根據(jù)總流的動量方程可以推導(dǎo)出單寬流量q與躍前水深y1、躍后水深y2之間的關(guān)系,從而得到光滑底板上能量損失ΔE及消能率η：

(12)

(13)

2 水躍現(xiàn)象數(shù)值模型

為驗證SPH方法研究正弦形底板上水躍現(xiàn)象的可行性及適用性，選取文獻[2]中的試驗數(shù)據(jù)與SPH模擬值作對比。文獻[2]中試驗?zāi)Ｐ褪褂玫木匦嗡鄄捎糜袡C玻璃制作，寬0.446 m、高0.6 m、長7.6 m，并在水槽底部鋪設(shè)鋁制正弦形波狀底板。試驗布置如圖1所示，圖中Lj為躍長；Lrj為旋滾區(qū)長度；Um為斷面最大流速；S為波長，T為波幅。研究的波形有2種：S均為68 mm，T分別為13 mm、22 mm。由于正弦形底板凹凸不平，易引發(fā)空蝕問題，Ead等[21]研究發(fā)現(xiàn)將正弦形底板設(shè)置在光滑水平面以下，波峰與光滑底板位于同一高度時，能有效避免壁面凹凸不平引起的空蝕問題。

圖1 試驗布置示意圖[2]

根據(jù)試驗?zāi)Ｐ筒捎肧PHYSICS_2D軟件建立相應(yīng)的數(shù)值模型，如圖2所示。模型主要結(jié)構(gòu)包括蓄水池、水槽、尾水收集池。正弦邊界采用公式y(tǒng)=Asin(π/0.034x)建立，當參數(shù)A(A為波幅的1/2)分別為0.006 5和0.011時可得T分別為13 mm、22 mm，S為68 mm的正弦形底板。模型由2種粒子構(gòu)成，分別為可流動的水粒子和固定邊界粒子。各工況閘孔自由出流流量調(diào)節(jié)通過蓄水池的水位調(diào)節(jié)實現(xiàn)；流經(jīng)水槽的水流暫存在收集池以防水粒子運動到計算域以外，對數(shù)值計算產(chǎn)生干擾。建模時為減小運行時間預(yù)先在整個水槽段添加一定深度的水域。

圖2 數(shù)值模型示意圖

對文獻[2]中5組工況進行數(shù)值模擬，具體參數(shù)如表1所示，表中y2*為光滑壁面躍后水深。數(shù)值運算時選取的主要參數(shù)如下：核函數(shù)為5次樣條函數(shù)；黏性項處理采用亞粒子湍流模型；邊界采用動力邊界條件。模型粒子總數(shù)為96 476個，其中邊界粒子4 610個，初始粒子間距為5×10-3m，光滑長度為3.68×10-3m,時間步長為0.01 s，模擬時長為10 s。

表1 工況參數(shù)

3 結(jié)果分析

3.1 模型驗證

由于本數(shù)值模型只能在初始狀態(tài)添加水域進行計算，故為滿足試驗工況所需的恒定入流條件，數(shù)值建模時在水池中添加一個推板，即推板模型[22]。通過設(shè)定推板的速度對閘口附近水域進行補給，以保證蓄水池中水位恒定。分析得知各試驗工況閘口處流量的SPH模擬值與物理模型試驗值間的誤差分別為3.9%、6.6%。由流量對比(圖3)可得：推板模型能夠在短時間內(nèi)達到各工況所需的穩(wěn)定流量且流量值出現(xiàn)波動的情況較少，表明推板模型能夠有效解決正弦形底板水躍研究中的流量恒定問題。

將工況A的SPH水面線模擬值和已有文獻試驗值及模擬值作對比,如圖4所示。由圖4可得各工況水面變化趨勢基本保持一致，SPH數(shù)值模擬結(jié)果能夠充分展現(xiàn)水面線形態(tài)的變化。

圖3 閘口處流量對比

圖4 工況A水面線對比

將5種工況水面線的SPH模擬值與文獻[2]試驗值作對比,如圖5所示，可知兩者基本吻合，僅在某些工況，如工況B、C旋滾區(qū)段水面有一定的偏差。考慮到水躍是一種較為復(fù)雜的流態(tài)，在水躍區(qū)段內(nèi)主流區(qū)和旋滾區(qū)之間的水質(zhì)點不斷摻混，伴隨大量的氣泡產(chǎn)生或摻入，導(dǎo)致水面劇烈波動，部分粒子飛濺，造成試驗值存在較大誤差。工況D和E流量相同、波形不同時，水面線變化趨勢基本相同，水面線基本重合，說明波形的改變對水面線高度的影響較小，分析其原因可能是由于波狀底板位于水平面下方，沒有伸入流體內(nèi)部；且水面的高度遠大于波狀底板波幅值，故導(dǎo)致正弦形底板上波形的變化對水面線影響相對較小。

圖5 各工況水面線對比

從5種工況水面線的整體情況來看，SPH數(shù)值模擬結(jié)果和試驗結(jié)果吻合較好，說明SPH方法能夠較為準確地反映水躍過程中水面的變化趨勢及形態(tài)，在處理波狀消力池底板上的水躍特性方面具有適用性及可靠性。

3.2 速度云圖、矢量圖

由工況A速度分布(圖6)得：水躍段包括速度為正的順流及速度為負的逆流，其中逆流主要聚集在水流表面和底部的主流處，表面由于部分水流回流形成逆流，底部因正弦形底板的波狀結(jié)構(gòu)形成波內(nèi)旋渦導(dǎo)致流速為負；水躍段存在表面旋滾，且有大量氣體摻入，速度分布波動較大，當?shù)装鍨檎倚螘r波內(nèi)形成旋渦，相當于增加了邊界層的厚度和粗糙度，近壁面及表面流速較小，中間層流速較大；當水流流經(jīng)明渠段時受邊界層影響，表面流速最大，近壁面處最小。

圖6 工況A速度云圖

圖7 工況A速度矢量

圖8 工況B速度矢量

結(jié)合工況A、B的速度矢量(圖7、圖8)可得：水躍段內(nèi)有明顯水流破碎和摻氣現(xiàn)象；在閘口開度相同時，隨著Fr的增大，水流破碎、摻氣現(xiàn)象加劇并且存在明顯渦結(jié)構(gòu)，分析其原因主要是水躍段流速分布不均勻，流層間的水質(zhì)點有相對運動，導(dǎo)致流層間產(chǎn)生的內(nèi)摩擦切應(yīng)力在水流表面為順流方向，而在水流底部為逆流方向，故形成了順時針流向的旋滾。對比兩工況旋滾得：同等條件下旋滾傾斜程度隨Fr的增大逐漸增加，即形成旋滾的長度逐漸縮短，旋滾內(nèi)的渦結(jié)構(gòu)愈發(fā)明顯，渦團尺寸逐漸增加，這將意味著旋滾的影響范圍逐漸變大。

3.3 躍長

水躍段水流劇烈翻滾摻混使水面形態(tài)急劇變化，旋滾的尾部水流不斷匯入主流并被逐漸淹沒，直到水流表面不會因為水流的翻滾而產(chǎn)生較大波動時流態(tài)演變?yōu)槊髑鲃?。根?jù)圖7可得，水躍前部1—1斷面至形成旋滾的2—2斷面之間的距離為旋滾區(qū)長度Lrj；1—1斷面至水流經(jīng)過旋滾又回到較為穩(wěn)定的自由水面3—3斷面之間的距離為Lj。

將5種工況Lrj、Lj模擬值與試驗值及光滑底板躍長試驗值Lj*對比得出,Lrj、Lj模擬值相比對應(yīng)試驗值的誤差分別在7%、8%以內(nèi)；Lj試驗值較光滑底板減小30%左右，表明正弦形底板上的躍長較光滑底板明顯縮短，若應(yīng)用在實際工程中能夠有效縮短消力池長度。結(jié)合表1對比工況A、B、C的Lrj、Lj如圖9所示，表明試驗值和模擬值變化趨勢是一致的；且同等條件下，Lrj、Lj均隨著Fr的遞增呈遞增趨勢。

圖9 躍長、旋滾區(qū)長度對比

3.4 躍后水深

圖10 共軛水深、躍后水深對比

將圖7中水面變化較穩(wěn)定的3—3斷面處水深值作為躍后水深y2。為比較共軛水深與躍后水深的關(guān)系，將A、B、C這3種工況的共軛水深值統(tǒng)一縮小至原來的1/10與其躍后水深y2繪制在圖10中，對比y2模擬值與試驗值及光滑底板躍后水深y2*得：共軛水深、躍后水深的變化趨勢相同；y2模擬值與試驗值之間的誤差在15%以內(nèi)，y2試驗值較y2*減少30%左右；比較工況A、B、C得同等條件下y2、y2*隨Fr的增加而增加；結(jié)合表1中工況D、C可發(fā)現(xiàn)波形對躍后水深影響較小。采用無量綱分析法探究躍后水深、共軛水深與躍長間的相互關(guān)系。結(jié)果表明：y2和Lj變化趨勢呈正比關(guān)系；當共軛水深呈遞增趨勢時，躍長呈減小趨勢，兩者呈反比關(guān)系，表明波狀底板雖然能夠減小躍長及躍后水深，但是不能同時減小共軛水深和躍長。

3.5 消能率

水躍消能是泄水建筑物中常用的消能形式之一，高速水流經(jīng)過消力池后形成水躍，躍后流速的急劇衰減消耗了上游水流的大部分能量，從而起到消能的目的。因此，在工程造價預(yù)算范圍內(nèi)對消能工進行合理設(shè)計，以盡可能地提高消能率。

圖11 流速分析

將5組工況正弦形底板消能率模擬值、試驗值及光滑底板消能率進行對比，見表2，可知當Fr=4～8時，光滑消力池底板的消能率范圍為39.36%～66.12%，正弦形消力池底板的消能率范圍為47.08%～73.12%，因此正弦形底板較光滑底板消能率明顯提高了10%左右，且5種工況的消能率均在47%以上。由表2還可得知，消能率試驗值和模擬值吻合度較高，變化規(guī)律相似。在所研究的Fr范圍里，光滑底板和正弦形底板的消能率均隨Fr的增加而增加；對比工況A、B、C可得，當波形相同時消能率隨Fr的增加而增加；對比D、E兩種工況可得當Fr相同時，正弦波形不同時消能率較為接近，說明其他條件相同時，正弦波形對消能率的影響不明顯，正弦形底板消能效果明顯。

表2 各工況不同F(xiàn)r時的消能率

4 結(jié) 論

a. 推板模型能夠獲取試驗所需的恒定流量條件，使模擬流量誤差在7%以內(nèi)；流量能在較短的時間內(nèi)達到穩(wěn)定且很少出現(xiàn)數(shù)值波動。

b. 數(shù)值模擬水面線與試驗值基本吻合，僅存在較小波動，驗證了SPH方法研究正弦形底板水躍現(xiàn)象的可行性及適用性。

c. 水躍段流速分布不均勻且波動較大；相同條件下，隨著Fr的增加摻氣量及水面破碎現(xiàn)象越來越明顯，旋滾的影響范圍逐漸變大。

d. 正弦形底板躍長、躍后水深較光滑底板均較小30%左右，表明波狀底板能夠有效降低躍后水深，減小消力池護坦高度；對比模擬工況發(fā)現(xiàn)水面線、躍長受波形影響較小；躍長、躍后水深的變化趨勢一致，但共軛水深與躍長變化趨勢相反，故不能同時減小躍長和共軛水深。

e. 所模擬工況斷面平均流速、最大流速分布規(guī)律具有相似性，即在水躍段存在較大波動、衰減，明渠段基本保持穩(wěn)定；正弦形底板消能率較光滑底板提高10%左右，且5種工況消能率均在47%以上。