基于SPH方法的二維寬頂堰溢流數值模擬

顧聲龍　吳玉帥　解宏偉　袁曉偉

摘要：光滑粒子水動力學方法（SPH）是一種基于純拉格朗日思想的無網格粒子方法，在眾多科學和工程領域得到了初步應用。采用SPH方法對二維寬頂堰溢流進行了數值模擬，并利用機群進行并行計算來提高計算效率，將模擬數據與試驗數據進行了對比。結果表明，SPH方法可以模擬出二維寬頂堰自由出流和淹沒出流時水流流態、水跌、回流區、水躍現象，計算出溢流流量、流速等水力參數，并且采用數據擬合得出與理論相符的斷面流速分布圖。結果表明了SPH方法可以較好模擬二維寬頂堰溢流中的流速分布以及水面曲線。

關鍵詞：SPH；二維寬頂堰溢流；數值模擬；自由面計算；并行計算

中圖分類號：TV652.1 文獻標志碼：A 文章編號：1672-1683（2017）00-0127-05

近年來SPH方法在自由表面流動模擬中取得了較大進展，能夠模擬傳統有網格算法較難處理的自由面翻卷破碎、液艙晃蕩、流固耦合等現象，得益于其拉格朗日粒子特性；SPH方法中粒子不僅用于估算場變與近似控制方程，也代表介質系統，具有諸如密度、壓力、速度、內能等宏觀物理量，相當于物質點。SPH方法既具有了拉格朗日方法描述物質界面準確的優勢，又兼備歐拉方法容易計算的長處。因此，SPH法已用來研究漬壩、自由表面流動、多相流等問題。

本文采用SPH方法來模擬二維寬頂堰溢流現象，利用寬頂堰溢流實驗實測數據，采用集群并行計算方式，來研究SPH方法在寬頂堰自由出流和淹沒出流水力特性方面的適用性及模擬的精度。

1SPH理論

1.1SPH方法的基本方程

SPH方法，本質上是把真正的空間近似劃分，最小的劃分空間由一系列表面絕對光滑的粒子填充。這些空間粒子常常用坐標、速度、質量、內能等參量表征。推導SPH的基本方程，分為兩步：第一步，給出SPH方法的積分表示法，即在整個變量場中，對任意位置應用核近似法進行積分近似；第二步給出其粒子近似值，對給出的函數積分表達式，通過對在其光滑長度范圍內的粒子疊加求近似，得到其粒子的近似公式。

1.2SPH方法核近似與粒子近似

在SPH方法中，函數f（x）的積分表示式如下：

（1）

（2）

（3）式中：h為光滑長度，用來定義光滑函數影響區域的光滑長度。由于W不是狄拉克函數，故式（3）的積分表示式只能是近似式，上式即為核近似方程的初始形式。

（4）

（5）

式（5）即為SPH方法的粒子近似法。在粒子i處的函數的粒子近似式最終可寫為

（6）

式（6）說明了粒子i處的任一函數值可通過應用光滑函數對其緊支域內所有粒子相對應的函數值進行加權平均近似求得。

1.3控制方程

將描述流體運動的質量守恒方程和動量守恒方程（Navier-Stokes方程）表示為

（7）

（8）

（9）

（10）

2數值模擬

2.1模型建立

本文選取寬頂堰溢流為試驗模型，見圖1，參數：堰高P=0.08 m，堰寬B=0.15 m，堰厚δ=0.4m，堰上水頭為H=0.065 m，δ/H=6。二維計算模型的幾何和物理參數與該試驗模型相同，由于SPH方法中入流問題仍是一個未解決的難題；故本文入流方式采用“大水箱”，該方法可以在一定范圍內保持模型的“恒水位，對計算結果影響不大。計算模型共計20 256個粒子，邊界粒子2 125個，本文利用集群進行并行計算；模型參數：核心函數為五次樣條核函數；時間步長算法為預估校正法；密度過濾為謝巴德濾波器方法；黏性處理方法為層流和亞粒子尺度紊流；狀態方程為微可壓；邊界類型為邊界排斥力型；時間步長為1×10^-5s，并且可變；柯朗弗里德里希一列維數為0.2；光滑長度系數為0.92；模擬溢流時長為10 s，集群共運行1.5 h。由于試驗模型的二維性較強以及構建三維計算模型的不易之處，所以本文采用二維模型；分析模擬結果時，將試驗模型數據轉換為單寬數據與模擬結果進行對比即可。

2.2結果及分析

2.2.1水流形態

本文試驗模型中h_s為下游水位高出堰項距離，H₀=H+α₀v²₀/2g；由圖2可知，寬頂堰溢流水面有兩次跌落；堰坎首端水面跌落是由于水流經過堰坎時，在縱向受到邊界的約束，過流斷面面積減小，流速增大，勢能減小。水面最大跌落處形成收縮斷面C-C；而后，由于堰頂阻力，使水面形成壅水曲線，逐漸接近堰頂斷面的臨界水深。下游水位較低，在堰坎末端再次出現跌落。當下游水位大于堰高后，隨著下游水位超過堰項水深幅度的增大，堰項在收縮斷面后發生波狀水躍，但這時下游水深并不影響收縮斷面水深；當h_s大于堰項收縮斷面水深的躍后共軛水深時，成為淹沒出流，下游水深將影響堰頂水位變化，堰上水位被壅高，堰頂呈緩流，見圖3。實驗研究得出h_s/H₀=0.8是寬頂堰是否淹沒的判別條件，在以下內容中會通過試驗數據得出這一比值。

圖4（a）-圖4（d）四幅圖為計算模型1.04 s、2.66 s、3.98 s、9.98 s四個時刻的水流形態。二維寬頂堰溢流模型初始狀態中水面線高度為0.145m，在水流達到穩定狀態前，堰前水面線的變化可明顯體現出寬頂堰對水流的影響，見圖4（a）中斷面6.92 m處。當計算模型達到穩定狀態后，水面線在寬項堰前基本保持平行而且有壅高現象，寬項堰阻礙作用得到了體現，見圖4（b）和圖4（c）中斷面6.8～7.0 m處；堰項水流出現了第一次跌落之后，隨后的水面線與堰頂基本保持平行，計算模型在出堰后成功地展示了出堰流的第二次跌落和寬項堰末端有回流區現象的存在，見圖4（b）和圖4（c）中斷面7.2～7.6 m處。這些特性與已知的試驗模型寬頂堰堰流特性非常相符。并且從圖4（a）-圖4（c），可明顯看出自由出流的水躍現象；圖4（d）呈現出淹沒出流時的水流流態，由圖4（d）可得，淹沒出流的水面曲線變化平緩近乎和堰項平等，在堰出口處，下游水深稍有回升，所以堰下游水位稍高于堰頂水位；在斷面7.5 m處還有波狀水躍，由于計算模型入流方式為“大水箱”方式，在9.98 s時，水位有一定的下降，故波狀水躍現象呈現的并不十分明顯。

2.2.2流量及流速

試驗模型流量通過直角形三角薄壁量水堰來測量，公式如下：

Q=1.365H^2.48 （11）

其中：H=H₁-H₂。在與本文二維計算模型相同水頭下，試驗模型數據為：Q=6.5×10^-4m³/s：單寬流量Q_單=Q/B=4.3×10^-3m³/s；下游水位超過堰項水深h_s=0.056 m，堰前斷面平均流速v₀=Q/A=0.19 m/s，故h_s/H₀=0.84，據此可說明發生了淹沒出流。

計算模型流量采用積分求和方式，步驟如下。（1）首先從結果文件中選取計算斷面的水平流速，選取方法為：先從結果文件中讀取計算斷面上水平速度以及豎直坐標，讀取依據為水平坐標是否等于設定的計算斷面位置；（2）選取數據之后，將水平速度和豎直坐標按豎直坐標進行排序，然后根據公式（12）進行計算，見圖（5）。

（12）

流量計算斷面選為x=6.75 m，距離堰前端△x=0.25 m，△x/H=3.8；計算結果見圖（6）。

二維計算模型中，當t<1.2 s，模型剛開始運行，還未達到穩定狀態，其流量變化幅度比較大，即圖6中N區；當1.2 s6.9 s后，隨著模型運行，計算模型并沒有維掙恒定水位，故計算流量有一定程度的下降，即圖6中0區。

圖7展示了堰前斷面6.75 m在t=2.66 s流速變化情況，由圖可知，在初始狀態中入口處水流速度很穩定，所以從渠底到堰高這段距離內，流速分布比較均勻；從堰高至水流表面，水流速度變化幅度很大。對數據進行擬合，可得變化曲線為三次方程。由于計算粒子數目有限并且一些參數設定帶來的計算誤差等原因，并沒有得到流速分布的對數方程；不過斷面上流速的變化趨勢與對數方程很相似。圖8展示了堰上斷面7.2 m在t=2.66 s的流速變化情況，該斷面上流速明顯比6.75 m斷面上流速大，主要是由于過流斷面面積變小。

圖9展示了在t=2 66 s，從斷面6.75 m到斷面7.5 m上平均流速變化情況。由圖可得，斷面6.75 m至斷面7.0 m屬于堰前，流速變化不大，分布比較均勻，有稍微的增長趨勢；斷面7.0m至斷面7.4 m屬于堰上斷面，流速快速增長，變化趨勢為先增大后減小，最大值在斷面7.2 m上；原因是隨著堰上斷面面積變小，水流動能增加，勢能減小，并且斷面7.2 m為收縮斷面，收縮斷面之后水深有一定的增長，水面線與堰項基本保持平行。斷面7.4 m至斷面7.5 m屬于堰后，流速變大，是由于水流勢能減小，動能增加。從圖7、圖8、圖9可看出，流速存在一定的震蕩現象，主要原因是：（1）計算粒子數目不足，導致結果數據過少；（2）并沒有實現恒定流的入流方式，使水頭在逐漸漸小。（3）模型參數設置上還未達到最佳，還需進一步調試。

3結語

本文采用SPH方法對寬項堰溢流進行了數值模擬，得出了寬項堰溢流流態、流量、流速并且觀察到了水面跌落、回流區、水躍等現象，SPH方法模擬數值結果與試驗數據相近以及所得到的撞擊、翻卷及碎波的形態、大小和位置都與實驗觀察吻合，證明了SPH方法在模擬自由水面出流方面的獨特優勢。由于“入流”問題的存在，對模擬結果帶來了一定的誤差，在斷面流速變化圖上可以看到；若模型保持水位恒定，模擬所得的斷面流速變化曲線將滿足對數方程。本文利用集群進行并行計算，計算效率高，節約了計算時間。

在下一步的研究中，建立三維模型和解決入流問題，使模擬結果與實際更相近以及水流形態與實驗觀察更吻合，使SPH方法在水動力學方面得到了更廣泛的應用。