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

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

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

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

文中應(yīng)用SPH方法模擬了拱壩中孔泄流沖擊水墊塘過(guò)程，并對(duì)泄流過(guò)程中流速、壓強(qiáng)等進(jìn)行分析，同時(shí)研究了下游水墊塘對(duì)泄流消能的影響。

1 SPH方法原理

1.1 基本原理

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

(1)

式中：f為場(chǎng)函數(shù)；Ω為計(jì)算域；x為坐標(biāo)向量；h為光滑長(zhǎng)度；W為核函數(shù)。在SPH方法中習(xí)慣用角括弧對(duì)核近似算子進(jìn)行標(biāo)記。

應(yīng)用粒子近似法對(duì)核近似方程式(1)進(jìn)一步近似，近似后的函數(shù)表達(dá)式為：

(2)

式中：ρ為密度；m為質(zhì)量；i、j分別為粒子序號(hào)。

1.2 控制方程

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

(3)

動(dòng)量守恒方程為：

(4)

能量守恒方程為：

(5)

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

1.3 人工黏度

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

(6)

其中：

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

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

1.4 密度近似法

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

(7)

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

1.5 邊界處理

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

第一種虛粒子稱為固壁虛粒子，可以對(duì)邊壁附近的實(shí)粒子施加一個(gè)強(qiáng)排斥力，該排斥力的表達(dá)式為：

(8)

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

第二種虛粒子稱為鏡像虛粒子，其密度、光滑長(zhǎng)度和速度大小與對(duì)應(yīng)的實(shí)粒子相同，速度方向相反。邊壁粒子布置方案如圖1所示。

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

2 數(shù)值計(jì)算

2.1 模型驗(yàn)證

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

圖2 潰壩水位模型

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

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

圖3 不同時(shí)刻的水體流態(tài)與文獻(xiàn)[10]的結(jié)果對(duì)比

選取了0.2、0.4、0.8 s三個(gè)時(shí)刻，繪制出物理模型試驗(yàn)結(jié)果中的自由液面曲線和等壓線圖與本文的模擬結(jié)果進(jìn)行對(duì)比，如圖4所示。

圖4 自由液面比較曲線

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

表2 自由液面的相似度

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

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

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

圖5 中孔泄流水位模型

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

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

圖7 粒子的密度分布情況

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

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

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

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

3 結(jié)論

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

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