基于SPH法的正弦形消力池底板水躍現(xiàn)象數(shù)值模擬

薄夫萍，吳海濤，王 星，張?jiān)圃?，顧聲?/p>

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

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

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

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

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

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

(1)

(2)

(3)

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

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

(4)

(5)

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

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

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

(6)

(7)

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

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

(8)

(9)

(10)

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

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

(11)

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

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

(12)

(13)

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

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

圖1 試驗(yàn)布置示意圖[2]

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

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

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

表1 工況參數(shù)

3 結(jié)果分析

3.1 模型驗(yàn)證

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

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

圖3 閘口處流量對(duì)比

圖4 工況A水面線對(duì)比

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

圖5 各工況水面線對(duì)比

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

3.2 速度云圖、矢量圖

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

圖6 工況A速度云圖

圖7 工況A速度矢量

圖8 工況B速度矢量

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

3.3 躍長

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

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

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

3.4 躍后水深

圖10 共軛水深、躍后水深對(duì)比

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

3.5 消能率

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

圖11 流速分析

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

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

4 結(jié) 論

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

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

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

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

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