退水渠道水力特性模擬的SPH方法及工程應用

高仁祖，馬曉泉，李向東，顧聲龍

(1.青海省水利水電勘測設計研究院有限公司，青海 西寧 810001；2.青海大學水利電力學院，青海 西寧 810016；3.黃河上游生態保護與高質量發展實驗室，青海 西寧 810016)

0 引言

實現生態文明、綠色發展是我國西北社會經濟發展的重要目標。隨著蘭-西城市群發展規劃理念的提出，加快推進“一優兩高”戰略以及“大西寧”的規劃構建，“引黃濟寧”工程作為解決區域發展不平衡不充分問題、實施扶貧開發和改善民生具有重要意義[1]。“引黃濟寧”工程分為引水工程和供水工程兩部分，工程規劃從黃河干流上游龍羊峽水庫引水，經隧洞穿越拉脊山自流輸入到湟水南岸，向湟水干流及南岸、城鄉生活、工業園區、農業生產、生態環境建設保護供水，以提高水安全保障能力。

在供水工程中，為保證渠道及重要建筑物的運行安全和渠道分段檢修，滿足渠道安全運行等需要，結合實際地形溝道設置相應的退水建筑物。由于退水渠的特殊性，退水流量較大，坡度較陡，為防止水流對退水出口處地面造成沖刷，退水渠末端設置消力池進行消能。因此，退水渠的水力特性研究對工程設計和運行管理具有十分重要的意義。通常規模以上重要退水渠的水力特性研究一般以物理模型試驗[2]和數值模擬[3]為主，由于資金有限以及場地的限制，所以本次研究使用數值模擬方法。本文基于傳統的拉格朗日流體力學光滑質點流體動力學(SPH)方法，對退水渠的水力特性進行數值模擬。

SPH(Smoothed Particles Hydrodynamics)方法首次被Gingold和Monaghan提出[4]，近30多年的發展，在流體動力學領域應用廣泛，比如，潰壩流[5]、自由表面流動[6]、多相流[7]等，相應成熟的SPH程序軟件有DualSPHysics[8]、SPHinXsys[9]、GPUSPH[10]等。DualSPHysics開源程序包已得到國內外許多學者的驗證[11-12]以及應用[13-14]，2019年張云云[15]等數值模擬陡坡水躍，結果與文獻中的試驗結果基本吻合，驗證了SPH出入流方法對陡坡水躍的數值模擬可行性。然而，國內外對退水渠開展的數值模擬以及物理模型試驗甚少，相關成果不多，尤其通過SPH方法數值模擬退水渠并未給出相應合理的研究成果。

因此，本文使用DualSPHysics 5.0開源程序代碼對退水渠進行三維數值模擬，在4種流量的工況下，分析退水渠水流的水面線變化、流速變化以及渠道末端消力池的消能效果，并得出一些趨勢性的結論，更好地指導工程設計和運行管理。

1 模擬方法

1.1 SPH方法的基本方程

SPH方法認為，如果場變量在計算區域中是連續且光滑的，則粒子上的場變量可以由周圍粒子的場變量值核估計而來。SPH方程的建立有2個重要的步驟。第一步是核計算，第二步是粒子估計。核估計完成場變量或場變量梯度的插值，而粒子估計則完成對核估計積分表達式中的粒子離散。

1.2 核估計與粒子估計

對于一個連續光滑函數f(x)，在定義域Ω上一點的函數值可以表達為：

(1)

其中，

(2)

式中，x—空間位置矢量；δ(x-x′)—狄拉克δ函數。

若δ函數用函數W(x-x′，h)替代[16]，則f(x)可近似為：

(3)

式中，W(x-x′，h)—SPH中的核函數，它的值取決于兩點之間的距離|x-x′|和光滑長度h，它與光滑因子k一同決定了光滑函數影響域的大小。

對于(3)式，進行粒子離散化為：

(4)

式中，mj—配對粒子j的質量；ρj—配對粒子j的密度。

對粒子i處場函數的粒子估計最終為：

(5)

其中，

Wij=W(xi-xj，h)

(6)

粒子i處場函數空間導數的粒子估計最終[17]為：

(7)

其中，

(8)

rij=|xi-xj|

(9)

式中，rij—粒子間距；xi—i粒子的矢量位置；xj—j粒子的矢量位置。

1.3 控制方程

在流體動力學中，流體運動的控制方程(Navier-Stokes方程)表示為：

(10)

(11)

式中，ρ、v—流體的密度、速度；P、fs—流體的靜水壓力、表面張力；μ—流體的動力黏性系數；t—時間。

在不可壓縮流動問題中，利用弱可壓縮[18]狀態方程表示如下：

(12)

(13)

式中，B—參考壓強；γ—常數，在水中γ=7；ρ0—水的密度1000kg/m3；c0—聲速，應在流體最大速度的10倍以上，來保證模擬流場的不可壓縮性。

SPH方法通過核函數和粒子近似，將拉格朗日形式的N-S方程式(10)和式(11)離散為：

(14)

(15)

式中，Wij—核函數；Πij—人工黏性項[19]。

2 模型布置與方法驗證

2.1 幾何模型布置

根據“引黃濟寧”工程中的眾多退水渠原型，如圖1所示建立了1∶1的退水渠數值幾何模型，其中主干渠是“引黃濟寧”工程中的供水渠道，分水池的建立是為了調取部分主干渠的水量對周邊縣級區域起到灌溉作用，本文是研究退水渠的水力特性，因此不再考慮分水池調水灌溉這一作用。圖1中節制閘與渠道入口的距離為10.6m，退水渠的渠道中心與渠道入口的垂直距離為5.95m，圓形分水池的半徑為3.8m，坡度為1∶6，消力池底坎高度為1m，其它渠段尺寸見表1。在數值模擬中，粒子間距dx為0.1m，時間步長為1.5×10-4s，模擬總時長80s，使用出入流開放邊界條件[18]來穩定流量，入流面積為2m×2.1m(寬×高)，4種退水流量分別為5.7m3/s(最小流量)、5.99m3/s(設計流量)、7.29m3/s(加大流量)和8.4m3/s(限制流量)。節制閘的開閉狀態在數值模擬時為完全閉合。

表1 幾何模型各部分的尺寸大小

圖1 退水渠三維幾何模型布置圖及坐標原點分布

2.2 本文的SPH方法在陡坡水躍中的驗證

為了驗證DualSPHysics方法在數值模擬陡坡水力特性方面的準確性及可行性，選取文獻中陡坡水躍研究的方案，文獻中葛旭峰[20]研究不同坡度不同流量的陡坡水躍特性，給出了試驗值與流體軟件Fluent模擬值一致的結論。如圖2所示，參考了文獻[20]中的幾何模型尺寸建立了數值幾何模型，選取文獻[20]中坡度為30°、陡坡長度為0.6m、流量Q分別為42.17、25.18L/s的工況，且臨界水深hk分別為0.1043、0.074m。在模擬過程中選擇出入流開放邊界條件[15]，粒子間距dx=0.01m，時間步長為1.5×10-4s，模擬時長為50s。圖2中x表示從坡底的折坡點開始x方向上的距離，即離消力池起點的絕對距離。

圖2 陡坡三維幾何模型布置圖

如圖3所示，給出了不同流量的試驗值[20]與DualSPHysics方法的模擬值的水面線變化曲線圖，x/hk表示離消力池起點的相對位置。從圖3中分析得出，模擬值變化與試驗值變化趨勢一致，然而此處為水流從急變流向緩流變化的轉折處，水面線起伏較大，流態處于復雜且隨機的湍流狀態，因此，此處水面線的模擬值與試驗值之間的誤差較大，其誤差大于其余段的水面線誤差。當流量Q=25.18L/s時，模擬值與試驗值[20]的平均誤差為2.47%；流量Q=42.17L/s時，模擬值與試驗值[20]的平均誤差為-8.7%。在水躍旋滾劇烈處，水面波動明顯，數值模擬情況較差，但水面線總體趨勢一致，通過水面線的平均誤差分析，得出DualSPHysics方法在數值模擬陡坡水力特性方面，具有一定的準確性以及可行性。

圖3 水面線對比曲線圖

3 計算結果及分析

圖4 四種流量下，沿程水面線變化曲線圖

如圖5所示，給出了設計流量5.99m3/s和加大流量7.29m3/s的流速分布云圖，對退水渠的流場進一步分析，從圖5中分析得出：整個渠道中的水流流態穩定，無嚴重的水流沖擊、飛濺等對工程不利的現象。整個退水渠流速最大處分布在陡坡段末，然而，流速越大的位置，受到水流的沖刷越大，因此，在具體工程中應該對沖刷嚴重的部位加強鞏固。通過整個渠道流速色值來分析，陡坡色值表較大，而水流進入消力池后，色值變低，從而更加反映了消力池削減水流流速的效果。如圖6所示，給出了4種流量下沿退水渠的流速變化曲線，從圖6中分析得出：流速整體沿陡坡增大，圖6中色值明顯增加；在陡坡段末流速受到消力池旋滾區域的影響，水流流態處于湍流，水流之間存在著許多個隨機的大大小小的漩渦結構，漩渦結構之間的剪切作用使得流速趨勢變得急劇下降，最后沿消力池段變得平緩。最大流速6.53m/s發生在陡坡段末x=35m的位置。

圖5 兩種流量下，流速分布三維云圖

圖6 四種流量下，流速分布曲線圖

最后，如圖7所示，對沿退水渠變化的消能率進行分析，從圖7中分析得出：消能率沿退水渠增大，在消力池段，整體消能率的增長變得緩慢。在消力池段，消能率隨著流量的增大，消能效果反而下降。在x=35m至x=45m之間，由于水流流態復雜，處于湍流旋滾區域，藍點線和紅點線消能率曲線呈先減小后增大趨勢。然而，總體消能率的趨勢沿程呈增長趨勢，總體消能率在55%左右，消能效果顯著。其中消能率的計算公式[22]如下：

圖7 四種流量下，消能率分布曲線圖

(16)

式中，Ei—斷面單位重量液體所具有的總能量；i=1，2，分別表示躍前、后斷面；zib—渠道底面面標高；vi—斷面平均流速；hi—斷面平均水深。

4 討論與結論

退水渠是渠系建筑物不可缺少的部分，在“引黃濟寧”的供水工程各段分布眾多，退水渠承擔著排泄灌溉渠道內剩余水量或入渠洪水的責任，不同流量的退水渠會產生不同的水力特性，這些變量最終將會影響整個退水渠排水的功能，影響著下游原始河道防沖消能，因此，對退水渠的數值模擬其意義重大。DualSPHysics 5.0開源程序代碼是基于SPH粒子法比較成熟的流體軟件，在國內外應用廣泛[8，11-12]，在自由表面流方面的數值模擬技術趨于成熟，其次三維的數值模擬相比于以往的二維流體仿真[11-12]，在模擬整個流場方面會更加真實以及數值結果精度會提高很多。本文基于DualSPHysics方法對退水渠進行三維數值模擬，通過四種不同流量的工況，分析研究退水渠水流沿程水深、水流流速以及消能率的變化規律，得出以下結論：

(1)當退水渠道體型一定，入流流量越大，渠道沿程水深在消力池中也隨之增大。根據GB 50288—2018《灌溉與排水工程設計標準》[21]對分水池段的水面線超高驗算，得出加大流量為7.29m3/s時退水渠的邊墻高度設計合理。

(2)對退水渠中水流流速分析，在陡坡段末x=35m處的渠道水流流速最大，此位置受水流沖刷影響較大，需在具體工程中對沖刷嚴重的部位加強鞏固。通過陡坡段以及消力池段的流速云圖發現，色值變化明顯，極大的反應了消力池的效能效果，同時消力池中的流態為湍流，表現出流態穩定且無明顯不利工程運行的問題。

(3)退水渠道消能率隨著入流流量的減小而增大，即消力池在小流量的工況下消能效果更好，總體消能率達到了55%左右，消能效果明顯。

(4)本文只是對流量變化下的退水渠的水力特性進行了數值分析，得出了一些趨勢性的結論，但是對于不同尺寸退水渠的數值分析沒有系統的論證，希望在今后的時間里做出合理的成果。