光滑粒子流體動力學方法在拱壩中孔泄流沖擊水墊塘模擬中的應用

張建偉， 主攀， 陳海舟， 杜宇

(1.華北水利水電大學 水利學院，河南 鄭州 450046； 2.上海道盾科技股份有限公司，上海 200131)

對拱壩中孔泄流過程進行分析，對優(yōu)化泄流方式、減少壩體空蝕和磨蝕的產(chǎn)生、削弱河床和岸坡沖刷具有重要意義。拱壩中孔泄流時具有水頭高、流速大、下泄功率大的特點[1]，物理模型試驗獲取準確的觀測資料難度較大、投資大、周期長，每次物理模型試驗只能針對較少的組合工況進行模擬，具有局限性[2]。數(shù)值模擬具有建模方便、仿真時間短、成本低的特點，適用于大壩泄流模擬，可作為物理模型試驗研究的補充手段。2009年,孫建等[3]在剛蓋假定的前提下，采用紊流數(shù)學模型模擬了溪洛渡高拱壩水墊塘內三維流場變化過程；2019年,周振興等[4]采用UG-FLUENT軟件對二河新泄洪閘泄流時的下游流場進行了模擬。由于泄流過程自由液面會發(fā)生大變形，傳統(tǒng)的基于網(wǎng)格的數(shù)值方法如有限差分法(FDM)和有限體積法(FVM)只能對下游水墊塘進行模擬，很難對泄流過程進行仿真，尋找新的數(shù)值模擬方法已成為目前泄流仿真研究的熱點。

光滑粒子流體動力學(Smoothed Particle Hydrodynamics，SPH)方法是一種基于拉格朗日形式的純粒子法，最初用于解決天體物理學中的三維問題，后由MONAGHAN J[5]對其進行修改，用于解決流體計算中的自由表面流問題。SPH方法可以用來處理自由表面大變形問題，能夠模擬出流場的整體形態(tài)和局部的破碎、飛濺、卷曲等細節(jié)。目前,SPH方法已被廣泛應用于潰壩、晃蕩、液滴撞擊和滑坡涌浪等問題研究中。李大鳴等[6]應用SPH方法并使用邊壁加載和相對加速度加載兩種加載方式模擬二維矩形容器內液體的縱蕩，發(fā)現(xiàn)在不同水深下邊壁加載和相對加速度加載兩種方式下的計算結果較一致；陳海舟[7]應用SPH方法模擬了潰壩問題、容器內液體晃蕩問題和液滴沖擊液面問題，在SPH基本方程的基礎上推導出了適用于廣義流體力學的SPH控制方程，為SPH方法在流體力學中的廣泛應用奠定了基礎；張建偉等[8]應用SPH方法模擬了潰壩沖擊水墊塘的過程，對沖擊過程中的流速和壓強進行了分析，提出SPH方法未來可應用于消能研究。

文中應用SPH方法模擬了拱壩中孔泄流沖擊水墊塘過程，并對泄流過程中流速、壓強等進行分析，同時研究了下游水墊塘對泄流消能的影響。

1 SPH方法原理

1.1 基本原理

SPH方法是一種純拉格朗日粒子法，用該方法做相關分析時先進行核近似，然后再對粒子近似[9]。通過積分表示法對場函數(shù)進行近似，核近似方程為：

(1)

式中：f為場函數(shù)；Ω為計算域；x為坐標向量；h為光滑長度；W為核函數(shù)。在SPH方法中習慣用角括弧對核近似算子進行標記。

應用粒子近似法對核近似方程式(1)進一步近似，近似后的函數(shù)表達式為：

(2)

式中：ρ為密度；m為質量；i、j分別為粒子序號。

1.2 控制方程

采用無黏假定將二維Navier-Stokes方程離散為SPH形式的流體控制方程。

(3)

動量守恒方程為：

(4)

能量守恒方程為：

(5)

式中:P為壓強；v為速度；vij=vi-vj；e為能量；N為粒子支持域中的粒子數(shù)。

1.3 人工黏度

SPH算法在模擬沖擊波問題時會產(chǎn)生非物理振蕩，采用MONAGHAN型人工黏度Πij能有效防止粒子相互接近時的非物理穿透。MONAGHAN人工黏度表達式為：

(6)

其中：

全國已建基層水利站29040個，12個省份實行了村級水管員制度。全國80%的縣成立了縣級農(nóng)村飲水安全工程管理機構。在全國27個省55個縣深入開展農(nóng)業(yè)水價綜合改革示范。取消4項行政審批事項，到2015年年底擬將水利部行政審批事項由48項減少到29項。55個全國示范縣市、150個省級試點縣市開展小型工程管理體制改革試點工作。

設置αΠ分別取1.0、1.2、1.4、1.6、1.8和2.0，通過試算發(fā)現(xiàn)，當αΠ取1.8時，既不會造成水流的間斷，同時也能很好地模擬水體大變形狀態(tài)，故αΠ設置為1.8。

1.4 密度近似法

在SPH方法中,光滑長度的變化和粒子的分配與粒子密度變化具有直接的關系。目前,關于粒子密度近似方面應用最廣泛的方法是密度求和法。任意粒子i的密度ρi的計算式為：

(7)

式中Wij為粒子j對粒子i影響的光滑函數(shù)。

1.5 邊界處理

SPH方法屬于無網(wǎng)格方法，只需要考慮固壁邊界條件的處理。為防止實粒子穿透邊壁，采用第一種虛粒子與第二種虛粒子相結合的方法來布置邊壁粒子。

第一種虛粒子稱為固壁虛粒子，可以對邊壁附近的實粒子施加一個強排斥力，該排斥力的表達式為：

(8)

式中：Fij為固壁粒子j對邊壁附近實粒子的排斥力；n1、n2分別取12和4；r0為截止半徑；D是由具體問題而定的參數(shù)，一般取與速度最大值的平方相等的量級值；rij為粒子i與粒子j之間的距離。

第二種虛粒子稱為鏡像虛粒子，其密度、光滑長度和速度大小與對應的實粒子相同，速度方向相反。邊壁粒子布置方案如圖1所示。

圖1 邊壁粒子布置方案示意圖

2 數(shù)值計算

2.1 模型驗證

拱壩中孔泄流仿真難點在于對自由液面大變形進行模擬，為驗證SPH方法是否能準確模擬流體大變形問題，選取文獻[10]中潰壩的算例進行驗證。采用FORTRAN語言編寫了SPH方法的計算程序，應用所開發(fā)的程序模擬了潰壩過程。模型為長0.584 m、高0.350 m的單寬矩形槽，上游水位設置為0.293 m，通過快速提閘模擬潰壩過程[11]，模型尺寸如圖2所示，計算參數(shù)見表1。

圖2 潰壩水位模型

表1 模擬潰壩過程的參數(shù)條件

圖3為模擬所得的不同時刻水體流態(tài)與文獻[10]的對比結果。

圖3 不同時刻的水體流態(tài)與文獻[10]的結果對比

選取了0.2、0.4、0.8 s三個時刻，繪制出物理模型試驗結果中的自由液面曲線和等壓線圖與本文的模擬結果進行對比，如圖4所示。

圖4 自由液面比較曲線

在圖4中物理模型試驗結果中的自由表面曲線和數(shù)值模擬結果中的自由表面曲線上分別取若干個點，計算點的位置距離差后與物理模型試驗值進行對比，得到各點的位置相似度，取平均值后得到自由液面曲線的相似度，結果見表2。

表2 自由液面的相似度

文中的模擬結果和文獻[10]的物理模型試驗結果十分接近，自由液面曲線綜合相似度為91.7%。相對于文獻[10]中MPS半隱式運動粒子法模擬結果，SPH方法模擬得到的液面流態(tài)更接近物理模型試驗的結果。同時，在圖4中的等壓線圖中可以看到，相對壓強為零的等壓線與自由液面曲線基本重合，說明壓強分布規(guī)律合理，驗證了文中所編寫的計算程序的準確性。在0.8 s時，潰壩洪水形成了漩渦，自由液體表面產(chǎn)生了大變形，說明SPH方法可以準確模擬自由液面大變形問題。

2.2 拱壩泄流沖擊水墊塘過程數(shù)值模擬

為模擬拱壩中孔泄流沖擊水墊塘過程，布置了一個長12.2 m、高3.5 m的單寬矩形水槽模型。距水槽左端4 m處設置一擋水建筑物，該擋水建筑物高2.5 m、寬0.2 m。距其底部1.1 m位置處設置一直徑為0.1 m的泄水孔用來模擬中孔泄流。下游布置有水墊塘。上游水位設置為2.5 m，下游水位設置為0.4 m。模型尺寸如圖5所示，計算參數(shù)見表3，初始粒子布置如圖6所示。初始粒子密度分布如圖7所示。

圖5 中孔泄流水位模型

表3 模擬拱壩泄流沖擊水墊塘過程的參數(shù)條件

圖6 粒子的初始布置位置示意圖

圖7 粒子的密度分布情況

中孔泄流沖擊水墊塘過程中的幾個關鍵時刻的壓強圖(左)和流速圖(右)如圖8所示。由圖8可以看出：上游水體從泄水孔流出后，由于中孔位置處受到的壓強較大，使得下泄水流的流速較快；在0.16 s時，上游孔口處粒子受到的最大壓強達到14 kPa，孔口的最大流速為3.2 m/s，出射水股具有極大的動能；從孔口流出后，下泄水體在空中形成射流，這一過程中粒子的勢能逐步轉化為動能；在0.46 s時，可以觀察到粒子的壓強進一步減小，流速增大，粒子的最大流速達到5.1 m/s；在0.49 s時，上游水體開始沖擊下游水墊塘，發(fā)生碰撞的一瞬間，由于受到下游水墊塘粒子的作用力，上游水體粒子的壓強突然增大到24 kPa，由于動量守恒，下泄水體粒子的速度減小為3.8 m/s；隨后在0.76～1.10 s時，由于上游水體的持續(xù)沖擊，下游水墊塘形成水舌，并逐漸擴大為涌浪，上下游水體粒子混摻在一起向下游邊界處運動；在1.32 s時，由于水體的持續(xù)沖擊，下游水墊塘再次出現(xiàn)水舌并持續(xù)翻滾形成漩渦；在1.49 s時可以清楚地觀察到漩渦中心處出現(xiàn)空穴，隨后上下游水體粒子再次混摻在一起形成涌浪，推高下游水位，說明出射水股動能逐漸轉化為水墊塘內水體的勢能。

圖8 拱壩中孔泄流沖擊水墊塘過程中水流的壓強圖(左)和流速圖(右)

為了更清晰地展示漩渦部分的流場分布，繪制出漩渦形成過程中幾個關鍵時刻的流速矢量圖(左)和局部放大圖(右)，如圖9所示。圖9中：在1.32 s時，由于上游水體的流速較大，受到?jīng)_擊的下游水墊塘形成水舌；隨后下游水墊塘的水舌進一步發(fā)展，上下游水體混摻后翻滾形成漩渦，在1.49 s時，在漩渦的中心處出現(xiàn)明顯的空穴，由于空穴處水體粒子相互擠壓，空穴處流速較大；隨后空穴部分逐漸被下泄水體填充，空穴處流速減小，實現(xiàn)消能，在1.73 s時，漩渦演變?yōu)橛坷?，在這一過程中水股動能逐漸轉化為內能，具有一定的消能效果。

圖9 拱壩中孔泄流沖擊水墊塘過程的流速矢量圖(左)及其局部放大圖(右)

3 結論

1)SPH方法將粒子近似法和拉格朗日公式有機結合在一起分析水流狀態(tài)，使得SPH方法中的粒子既可以看作插值近似點同時還攜帶著材料屬性，會在內部和外部荷載的相互作用下運動，可以很好地模擬出拱壩泄流過程中的自由表面大變形和沖擊過程。

2)拱壩中孔泄流過程中，出射水股沖擊下游水墊塘后上下游水體粒子混摻在一起，出現(xiàn)涌浪、漩渦和空穴等現(xiàn)象，說明下游設置水墊塘具有很好的消能效果。