基于SPH法的正弦形消力池底板水躍現象數值模擬

2019-03-26 08:56:16薄夫萍吳海濤張云云顧聲龍

水利水電科技進展 2019年2期

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

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

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

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

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

1 經SPH方法離散后的N-S方程

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

(1)

(2)

(3)

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

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

(4)

(5)

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

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

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

(6)

(7)

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

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

(8)

(9)

(10)

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

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

(11)

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

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

(12)

(13)

2 水躍現象數值模型

為驗證SPH方法研究正弦形底板上水躍現象的可行性及適用性，選取文獻[2]中的試驗數據與SPH模擬值作對比。文獻[2]中試驗模型使用的矩形水槽采用有機玻璃制作，寬0.446 m、高0.6 m、長7.6 m，并在水槽底部鋪設鋁制正弦形波狀底板。試驗布置如圖1所示，圖中Lj為躍長；Lrj為旋滾區長度；Um為斷面最大流速；S為波長，T為波幅。研究的波形有2種：S均為68 mm，T分別為13 mm、22 mm。由于正弦形底板凹凸不平，易引發空蝕問題，Ead等[21]研究發現將正弦形底板設置在光滑水平面以下，波峰與光滑底板位于同一高度時，能有效避免壁面凹凸不平引起的空蝕問題。

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

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

圖2 數值模型示意圖

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

表1 工況參數

3 結果分析

3.1 模型驗證

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

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

圖3 閘口處流量對比

圖4 工況A水面線對比

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

圖5 各工況水面線對比

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

3.2 速度云圖、矢量圖

由工況A速度分布(圖6)得：水躍段包括速度為正的順流及速度為負的逆流，其中逆流主要聚集在水流表面和底部的主流處，表面由于部分水流回流形成逆流，底部因正弦形底板的波狀結構形成波內旋渦導致流速為負；水躍段存在表面旋滾，且有大量氣體摻入，速度分布波動較大，當底板為正弦形時波內形成旋渦，相當于增加了邊界層的厚度和粗糙度，近壁面及表面流速較小，中間層流速較大；當水流流經明渠段時受邊界層影響，表面流速最大，近壁面處最小。

圖6 工況A速度云圖

圖7 工況A速度矢量

圖8 工況B速度矢量

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

3.3 躍長

水躍段水流劇烈翻滾摻混使水面形態急劇變化，旋滾的尾部水流不斷匯入主流并被逐漸淹沒，直到水流表面不會因為水流的翻滾而產生較大波動時流態演變為明渠流動。根據圖7可得，水躍前部1—1斷面至形成旋滾的2—2斷面之間的距離為旋滾區長度Lrj；1—1斷面至水流經過旋滾又回到較為穩定的自由水面3—3斷面之間的距離為Lj。

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

圖9 躍長、旋滾區長度對比

3.4 躍后水深

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

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

3.5 消能率

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

圖11 流速分析

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