高仁祖， 吳海濤， 顧聲龍*， 田麗蓉， 鄭溫剛

(1.青海大學(xué) 水利電力學(xué)院，西寧 810016；2.黃河上游生態(tài)保護(hù)與高質(zhì)量發(fā)展實驗室，西寧 810016)

1 引 言

溢洪道水流具有流速高和流量大等特點(diǎn)，為了保障大壩及下游安全，需要削減水流的動能。當(dāng)水流流經(jīng)臺階式溢洪道臺階表面時，水流會充分摻氣并發(fā)生旋滾，水流處于湍流狀態(tài)使得大小渦結(jié)構(gòu)之間存在剪切作用，從而消殺水流所具有的巨大能量，達(dá)到消能目的。

為了全面評估臺階表面的復(fù)雜水流狀況，并為設(shè)計實用提供可靠的信息，臺階式溢洪道的研究分為試驗研究和數(shù)值模擬。趙相航[1]通過不同體型臺階式溢洪道的研究，分析了臺階上的水面線、流速以及壓強(qiáng)等水力特性的分布規(guī)律。

相關(guān)研究表明，臺階式溢洪道消能顯著，臺階上的流動結(jié)構(gòu)比較復(fù)雜，自由表面的大變形以及水流破碎是研究臺階溢洪道的一大難題。光滑質(zhì)點(diǎn)水力學(xué)(SPH)方法處理大變形以及水流破碎方面，尤其在捕捉自由水面方面比較有優(yōu)勢。SPH(Smoothed Particles Hydrodynamics)方法首次由Gingold等[2]提出，并成功運(yùn)用到了流體動力學(xué)中。近20多年來，許多研究人員進(jìn)行了完善和改進(jìn)[3]，由于SPH方法在解決復(fù)雜的自由表面流上有獨(dú)特的優(yōu)勢，因此，應(yīng)用于許多流體動力學(xué)領(lǐng)域。如SPH方法應(yīng)用于研究自由表面流[4]、氣液固三相相互作用[5]及近海岸波浪破碎[6]等。

Husain等[7]采用開源的SPHysics程序及大水箱模型，研究了不同流量條件下臺階段上游水流區(qū)的壓力分布。Gu等[8]利用ParellelSPHysics開源程序?qū)Σ煌w型的臺階式溢洪道進(jìn)行了數(shù)值研究，采用推板模型穩(wěn)定入流流量，研究了臺階上的水深、速度和壓力變化。但是對于高雷諾數(shù)的湍流流動，需要考慮湍流帶來的諸多影響，需要在SPH算法中加入一定的湍流模型來優(yōu)化SPH算法。目前，常見的湍流模型有兩方程k-ε模型[6,9]和大渦模型[10,11]。

Gotoh等[10]首次成功地把大渦模型加入不可壓SPH方法中，并應(yīng)用于近岸孤立波的數(shù)值模擬。Ting等[11]對后臺階流動進(jìn)行了大渦模擬。Issa等[9]將RANS和LES模型分別加入SPH方法中，成功模擬了潰壩和近岸孤立波。隨后，Di Mascio等[12]在 δ -SPH 方法的基礎(chǔ)上提出了 δ -LES -SPH 模型。Meringolo等[13]運(yùn)用 δ -LES -SPH 方法數(shù)值模擬了重力波，分析了波的生成、傳播以及能量耗散隨時間的變化。Meringolo等[14]又進(jìn)一步對能量耗散進(jìn)行了研究分析，數(shù)值模擬了水片入水池和潰壩。Tripepi等[15]用 δ -LES -SPH 方法研究了孤立波與水下方形障礙物之間的關(guān)系。

相關(guān)研究表明，具有湍流封閉模型的 δ -LES -SPH 格式可以改善數(shù)值耗散，并能得到更準(zhǔn)確的計算結(jié)果。因此本文實現(xiàn) δ -LES -SPH 方法的C++程序，并利用經(jīng)典潰壩和寬頂堰進(jìn)行驗證，最后將 δ -LES -SPH 方法首次成功應(yīng)用于臺階溢洪道水流的水力特性研究，這對消能防沖工的設(shè)計具有工程指導(dǎo)性的意義，為具體的工程實踐提供有價值的參考。

2 δ -LES -SPH 模型

2.1 弱可壓縮的N-S方程

在拉格朗日型式下的弱可壓牛頓流體的N-S方程[14]表示為

(1)

(2)

Dr/Dt=u

(3)

(4,5)

本文是基于粒子法，因此方程(5)為SPS(Sub -Particle -Scale)湍流模型，也稱為Smagorinsky模式，在以往的大渦模擬[10]中，方程(5)用來封閉方程(4)，其中為應(yīng)變率張量，=(u+uT)/2，為Kronecker張量，λ和為粘性系數(shù)，tr()為對應(yīng)變率張量矩陣取跡，CS為Smagorinsky常數(shù)，Δ為粒子間距。弱可壓縮SPH模型中，狀態(tài)方程為[14]

(6)

(7)

2.2 δ -LES -SPH 模型

本文采用 δ -LES -SPH 模型，其中動量方程與連續(xù)性方程和文獻(xiàn)[16]的 δ -SPH 格式相近，就是在連續(xù)性方程中加入了擴(kuò)散項，動量方程中加入粘性項，目的是解決密度不連續(xù)帶來的數(shù)值計算不穩(wěn)定問題。δ -LES -SPH相比傳統(tǒng)的SPH以及 δ -SPH 格式有許多優(yōu)點(diǎn)，如δ -SPH中應(yīng)用人工粘度，而 δ -LES -SPH 中應(yīng)用了流體真實粘度；在 δ -LES -SPH 的連續(xù)性方程中，擴(kuò)散項系數(shù)和粘性項系數(shù)的求解與SPS湍流模型建立關(guān)系。δ -LES -SPH方法的N-S方程[14]為

(8)

(9)

(10)

式中下標(biāo)a和b為計算域內(nèi)任意2個粒子，a為相對于粒子a坐標(biāo)的梯度算子，Wa b為光滑核函數(shù)，本文使用Wendland C2核函數(shù)[17](光滑長度h=2dx,dx為粒子間距)，并將該核函數(shù)作為LES的濾波函數(shù)，a，ua，pa，ga和ra分別為第a個粒子的密度、速度矢量、壓力、重力加速度和位置矢量，Va為a粒子影響域的面積(本文計算域Ω為二維)。相比于 δ -SPH 格式，擴(kuò)散系數(shù)δ和粘性項系數(shù)α不再是常數(shù)，而是隨a粒子和b粒子變化的變量，δa b可以表達(dá)為

δa b=2 [δaδb/(δa+δb)]

(11)

(12)

(13)

式中Cδ為一個常數(shù)，取1.5，lLES為SPH濾波過程的一個參考量(lLES=4dx，d/dx為分辨率，d為初始自由水面的高度)，應(yīng)變率張量的計算公式為

(14)

(15)

(rb-ra)](rb-ra)/|rb-ra|2

(16)

(17)

動量方程(9)中，粘性項πa b表達(dá)為

πa b=[(ub-ua)·(rb-ra)]/(rb-ra)2

(18)

方程(9)的粘性系數(shù)αa b為

αa b=α+2 [αaαb/(αa+αb)]

(19)

(20)

(21)

式中CS為Smagorinsky常數(shù)，CS=0.12，ν為液體的動力粘滯系數(shù)。用式(12，20，21)來封閉 δ -LES -SPH 模型，當(dāng)研究劇烈自由表面流問題時，應(yīng)變率張量在h趨近于0時，流體互相碰撞時會出現(xiàn)無窮值，為了使數(shù)值模擬過程穩(wěn)定，根據(jù)文獻(xiàn)[14]取αa<0.2，δa<0.2，其物理意義在于消除弱可壓SPH數(shù)值模擬時帶來的虛假數(shù)值振蕩。

2.3 程序流程

本次開發(fā)的C++ δ -LES -SPH 模型程序的基本流程如圖1所示。

圖1 C++ δ -LES -SPH模型程序的基本流程

(1) 讀取幾何模型，并對幾何模型(計算域)進(jìn)行粒子化處理。

(2) 給每個粒子進(jìn)行粒子配對，其鄰近粒子搜索用到了鏈表搜索法。

(3) 遍歷每一個粒子，通過式(17)進(jìn)行密度梯度修正，以確保在自由液面處核函數(shù)粒子缺失條件下，密度梯度具有二階精度，計算的結(jié)果參與步驟(6)的計算。

(4) 通過虛粒子法進(jìn)行固壁邊界處理，其中遍歷每一個虛粒子，通過狀態(tài)方程反推出其密度。

(5) 遍歷每一個流體粒子，通過式(11，12，19～21)，求解出每個流體粒子的δa b和αa b，并參與當(dāng)前時間步內(nèi)的每個流體粒子密度變化量以及加速度的求解計算。

(6) 求解式(8，9)的左項(注：此程序中，沒有提前遍歷求解出每個粒子的核函數(shù)以及偏導(dǎo)數(shù)，而是每次遍歷粒子需要時，才去求解該粒子的核函數(shù)以及偏導(dǎo)數(shù))。

(7) 時間進(jìn)程采用蛙跳法推進(jìn)每一步求解，并且每500步輸出一次文件。最終需要進(jìn)行計算時間的判斷，直到計算達(dá)到要求。

2.4 邊界條件的處理

SPH方法中，靠近邊界的流體粒子支持域內(nèi)只有部分粒子，在固壁邊界外無流體粒子，出現(xiàn)了積分邊界截斷的問題，造成邊界附近的粒子穿過邊界，因此需對固壁邊界進(jìn)行處理，以獲得正確的計算結(jié)果。本文使用Adami等[18]提出的一種在流體外布置多層虛粒子實施固壁邊界的方法，此方法直接將內(nèi)部流體粒子的壓力插值到虛粒子上，如式(22)，再利用狀態(tài)方程反推出虛粒子的密度，虛粒子通過壓力梯度對流體施加邊界力，該邊界處理方法可以很好地適應(yīng)復(fù)雜幾何邊界問題。

(22)

式中ab為虛粒子的加速度。

3 算例驗證

3.1 經(jīng)典潰壩模型

本文采用經(jīng)典潰壩算例進(jìn)行驗證，與文獻(xiàn)[14,19]的模擬結(jié)果和Lobovsky等[20]的試驗結(jié)果作了對比分析。如圖2所示，初始的水箱長度L=5.367 m，布置L1為2 m和H1為1 m的矩形流體區(qū)域。流體粒子總數(shù)8萬個，粒子間距為0.005 m，時間步長為10-4s，模擬總時長為10 s。

圖2 潰壩數(shù)值模擬模型

圖3分別給出了無量綱時間為t(g/d)1/2=1.41，5.24和6.41時的流態(tài)。圖3(a)為文獻(xiàn)[20]的試驗流態(tài)，圖3(b)為文獻(xiàn)[18]的模擬流態(tài)，圖3(c)為本文的模擬流態(tài)。通過 δ -LES -SPH 數(shù)值模擬的流態(tài)與試驗的流態(tài)以及傳統(tǒng)SPH模擬得到的流態(tài)比較分析，發(fā)現(xiàn)數(shù)值計算過程穩(wěn)定，δ -LES -SPH 相比傳統(tǒng)的SPH方法，其數(shù)值模擬的流場更加明顯地反映了潰壩真實流動情況，說明 δ -LES -SPH 可精確捕捉自由水面。

圖3 t (g/d)1/2=1.41，5.24和6.41(從上到下)時的流態(tài)

圖4給出了t(g/d)1/2=9.45的αa值和δa值的模擬結(jié)果，其中Ma=0.05，在圖3的流體區(qū)域發(fā)生了翻滾水流的猛烈撞擊，耗散項主要作用于自由表面重聯(lián)和破裂而產(chǎn)生的渦旋周圍，此時αa值和δa值表現(xiàn)得很高，在大部分流體區(qū)域，αa值和δa值都接近于零。這些現(xiàn)象與文獻(xiàn)[14]的結(jié)果一致。在3.1節(jié)中，知道αa項和δa項都與應(yīng)變率張量成正比關(guān)系，因此，這兩項的變化值都出在了同一個圖中。整個計算域中，αa和δa的平均值為0.005和0.1。圖5給出了t(g/d)1/2=20時的渦量圖，并以無量綱渦量進(jìn)行了分析，圖5出現(xiàn)了諸多湍流漩渦斑點(diǎn)。

圖4 t (g/d)1/2=9.45的αa值和δa值的模擬結(jié)果

圖5 t(g/d)1/2=20時的渦量

圖6是t(g/d)1/2=1.079和1.75時的不同粒子間距下的自由水面曲線，選擇這兩時刻的原因是此時的潰壩流態(tài)比較穩(wěn)定。其中dx表示粒子間距，縱橫坐標(biāo)表示自由水面上粒子的水平和垂直距離。可以看出，當(dāng)粒子間距變小，粒子數(shù)量變多時，同一時間下，自由水面明顯呈收斂趨勢，圖6中dx=0.005和0.0025的自由水面基本一致，然而在同一工況下，dx= 0.005的模擬時長需要3天，而 dx=0.0025的模擬時長需要10天，最后兩者的模擬結(jié)果基本一致，為了節(jié)省計算時間，dx=0.005 作為最優(yōu)粒子間距進(jìn)行后面的案例分析。

圖6 t (g/d)1/2=1.079和1.75時的不同粒子間距下的自由水面曲線

綜上所述，本節(jié)模擬的結(jié)果以及相應(yīng)的對比分析，說明了自主實現(xiàn)的 δ -LES -SPH 算法具有穩(wěn)定性以及可靠性。

3.2 寬頂堰的驗證

圖7 寬頂堰數(shù)值模擬模型

圖8 流量隨時間變化的曲線

圖9給出了t=6 s時沿水平方向的流速分布云圖，可以看出，水流受到寬頂堰方形堰角的影響，在堰頂上游局部區(qū)域的水平速度為負(fù)，說明該區(qū)域發(fā)生了回流現(xiàn)象；水流流過堰頂流入消力池中，水舌下面區(qū)域的流速分布極不均勻，說明水流發(fā)生了旋滾以及回流現(xiàn)象且自由水面波動較大，在水舌下面區(qū)域的水流受到漩渦以及水質(zhì)點(diǎn)相對運(yùn)動的影響，加劇了水流能量耗散，使得此區(qū)域的水流流速急劇下降。

圖9 t =6 s時沿水平方向的流速分布云圖

圖10 沿堰頂表面水平方向的壓強(qiáng)變化曲線

通過寬頂堰的數(shù)值模擬以及結(jié)果分析可知,本次對 δ -LES -SPH 模型的實現(xiàn)，考慮了湍流這部分因素對計算帶來的影響，有效改善了數(shù)值計算結(jié)果。

3.3 臺階式溢洪道的數(shù)值模擬

表1 模型相關(guān)參數(shù)

圖11 臺階式溢洪道的數(shù)值模型

圖12給出了31階溢洪道在t=8.4 s時的水流流態(tài)，可以看出，臺階上游段的水面基本平穩(wěn)，臺階下游段水面破碎嚴(yán)重，水深略有緩慢上升趨勢，這與試驗觀測的水面現(xiàn)象相近。水流過臺階端角時出現(xiàn)了局部空化現(xiàn)象，而空腔上部水流為挑射流，如圖12(b)[1]試驗中也能觀察到局部空化現(xiàn)象，與試驗中觀察到的吻合，模擬中的31階溢洪道上的水流流態(tài)為跌落水流，因此會對臺階造成較大的沖擊和沖刷。

圖12 31階溢洪道在t =8.8 s時的數(shù)值模擬水流流態(tài)和31階溢洪道試驗時的水流流態(tài)

圖13給出了t=8.4 s的31階和8.8 s時的26階溢洪道的水深試驗值[1]和 δ -LES -SPH 大渦模擬值的沿程水深曲線，圖中橫坐標(biāo)無量綱化為x/L,坐標(biāo)L為溢洪道長度是70 m,x為監(jiān)測點(diǎn)到堰頂?shù)乃骄嚯x，h為沿程水深，結(jié)果表明,圖13(a)與無湍流模型的開源DualSPHysics的31階模擬值[22]相比，開源DualSPHysics的模擬值與試驗值的平均誤差為4.40%。溢流堰上的水深變化先逐漸降低，進(jìn)入臺階段水深又明顯增加，隨后水深變得平穩(wěn)。31階溢洪道臺階段上模擬的水流水深最大值為0.785 m，發(fā)生在臺階30附近，水流趨于平穩(wěn)后，水深達(dá)0.76 m，沿程的模擬水深值與試驗水深值的平均誤差約為4.20%，提高了4.55%的計算精度(精度計算是兩誤差之差除以原平均誤差)。其26臺階溢洪道上最大水深模擬值為0.869 m，水流趨于穩(wěn)定時，水深維持約在0.8 m左右，模擬水深與試驗水深的沿程平均誤差約為4.23%,這里提到的誤差是通過各個監(jiān)測點(diǎn)處的誤差求和再求平均得到。

兩種工況的對比結(jié)果表明，δ -LES -SPH 數(shù)值模擬和試驗值的結(jié)果吻合度較高，說明 δ -LES -SPH 能夠較為準(zhǔn)確地反映臺階式溢洪道上水面的變化趨勢以及形態(tài)。

圖13 31階和26階溢洪道的水深試驗值和模擬值的沿程曲線

圖14給出了t=8.4 s的31階和8.8 s時的26階臺階式溢洪道流速的試驗值[1]和模擬值的沿程曲線變化，圖14與無湍流模型的開源DualSPHysics的31階模擬值[22]都以相應(yīng)比例轉(zhuǎn)化為原型數(shù)據(jù)作了比較，分析得出，堰上的水流流速是先沿程增大，在臺階段增大到最大值又會沿程減小，隨后流速略有增大或減小的趨勢呈波浪狀分布，其試驗值和 δ -LES -SPH 模擬值以及DualSPHysics模擬的流速沿程分布趨勢一致。

圖14 31階和26階溢洪道的流速試驗值和模擬值的沿程曲線

試驗的最大流速為13.8 m/s,δ -LES -SPH 模擬的31階溢洪道的最大流速為13.9 m/s，試驗值和DualSPHysics的模擬值沿程平均誤差約為6.10%，其試驗值和本文的模擬值沿程平均誤差約為3.00%，提高了50.82%的精度；26階溢洪道模擬的最大流速為15.64 m/s，試驗中最大流速為14.993 m/s，其試驗值和模擬值沿程平均誤差約為4.60%。

比較 δ -LES -SPH 和DualSPHysics[22]的結(jié)果，發(fā)現(xiàn) δ -LES -SPH 的模擬結(jié)果與試驗結(jié)果吻合度較高，尤其在臺階下游段 δ -LES -SPH 試驗的流速都呈波浪形分布，而DualSPHysics的流速曲線表現(xiàn)較為平穩(wěn)，反映了 δ -LES -SPH 方法在曲線變化的拐點(diǎn)處具有很好的修正作用，在Tripepi等[15]用 δ -LES -SPH 方法研究孤立波與水下方形障礙物之間的關(guān)系時，也同樣提到了這一作用。

圖16對t=8.4 s的31階和8.8 s的26階臺階式溢洪道的消能率進(jìn)行研究分析，分析得出，消能率的變化從堰頂向下游沿程增加，其試驗值[1]和 δ -LES -SPH 模擬結(jié)果沿程分布趨勢一致,而且從整體趨勢來看，相比較開源DualSPHysics方法模擬得到的結(jié)果，31階中 δ -LES -SPH 模擬值與試驗值吻合較好，計算消能率的公式為

(23，24)

(26)

式中E1和E2分別為1-1斷面和2-2斷面的總能量,H0為2-2斷面到溢流堰堰頂?shù)拇怪本嚯x,H1為水庫水位到溢流堰堰頂?shù)拇怪本嚯x，h2和v2分別為2-2斷面的水深和流速，g為重力加速度,ΔE為損耗能量，α為溢洪道坡度，具體如圖15所示。

圖15 消能斷面選擇

對于31階溢洪道，試驗中最終的消能率為76.05%，通過開源DualSPHysics方法[22]模擬的結(jié)果分析得到的消能率為78.34%，與試驗值作比較，消能率的誤差約為3.01%；通過 δ -LES -SPH 數(shù)值模擬的消能率為75.89%，與試驗值作比較，消能率的誤差約為0.29%，在SPH方法中實現(xiàn)大渦模擬，使得模擬結(jié)果提高了90.37%的精度。對于26階的溢洪道，試驗中最終的消能率為 75.35%，通過 δ -LES -SPH 模擬的消能率為74.72%，與試驗值作比較，消能率的誤差約為0.26%。圖16中 δ -LES -SPH 數(shù)值模擬得到的消能率相比DualSPHysics方法得到的結(jié)果略偏小，這是因為DualSPHysics用的是人為調(diào)節(jié)擴(kuò)散項和耗散項的 δ -SPH 模型，而 δ -LES -SPH 用大渦湍流模型封閉了人為調(diào)節(jié)的擴(kuò)散項和耗散項，因此擴(kuò)散項和耗散項的變化與實際流場有關(guān)。由以上模擬結(jié)果來看，在計算中比 δ -SPH 模型要精細(xì)很多。DualSPHysics方法與 δ -LES -SPH 方法的邊界條件不同，然而邊界條件的給出是為了阻止流體粒子穿透固壁邊界，而對計算結(jié)果影響甚小。

綜上所述，通過曲線變化趨勢的比較以及相應(yīng)的平均誤差分析可知，在具有復(fù)雜湍流性質(zhì)的SPH數(shù)值模擬中加入 δ -LES -SPH 模型，優(yōu)化傳統(tǒng)SPH算法，使得數(shù)值模擬過程穩(wěn)定，且精確了模擬結(jié)果。以上結(jié)果也表明縮小比例尺后的數(shù)值試驗沒有對原有試驗產(chǎn)生影響。此次臺階式溢洪道數(shù)值模擬進(jìn)行28核并行處理計算，一組工況花費(fèi)了大約15天的時間,DualSPHysics方法一組工況在32核處理計算需要1天時間，相比之下添加復(fù)雜湍流模型，反而增加了所需模擬時間。

圖16 31階和26階溢洪道的試驗值和模擬值的消能率沿程曲線

4 結(jié) 論

本文就基于Mascio提出的大渦模型思想，自主實現(xiàn)了 δ -LES -SPH 程序，并通過潰壩案例與寬頂堰案例驗證了 δ -LES -SPH 程序，驗證結(jié)果表明，該算法在湍流流動數(shù)值模擬過程中具有較好的準(zhǔn)確性與穩(wěn)定性。最終，用 δ -LES -SPH 程序首次數(shù)值模擬兩種階數(shù)的臺階式溢洪道，并與現(xiàn)有文獻(xiàn)對比分析得出以下結(jié)論。

(1) δ -LES -SPH 方法可以很好地捕捉其自由水面，能夠準(zhǔn)確地反映其臺階上復(fù)雜變化的水流流態(tài)，臺階上的各斷面水深、各斷面的水流流速以及消能率的變化趨勢總體與試驗趨勢一致。31階的最大流速為13.9 m/s，出現(xiàn)在第20級臺階上，26階的最大流速為15.64 m/s，出現(xiàn)在第4級臺階上。

(2) 對31階臺階式溢洪道數(shù)值模擬，與無大渦模型的DualSPHysics方法比較，提高了沿程水深4.55%的精度、流速50.82%的精度、效能率90.37%的精度，而且在曲線拐點(diǎn)附近 δ -LES -SPH 的曲線變化更接近于試驗值，對曲線拐點(diǎn)處的數(shù)值結(jié)果有很好的修正作用，充分體現(xiàn)了 δ -LES -SPH 的優(yōu)勢。

(3) 由消能率的數(shù)值結(jié)果分析可知，臺階式溢洪道上的消能率沿程逐級增大，最后趨于穩(wěn)定，31階和26階的消能率分別達(dá)到了75.89%和74.72%，明顯31階在消殺下泄水流的能量方面優(yōu)于26階，充分體現(xiàn)了臺階式溢洪道的消能效果。這些數(shù)值結(jié)果對消能防沖工的設(shè)計具有指導(dǎo)性的意義，本文首次在臺階式溢洪道方面使用 δ -LES -SPH 方法，這是在具有復(fù)雜湍流流動研究方面的探索，也是創(chuàng)新，為未來開發(fā)高效的算法提供有價值的思路。