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

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

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

關(guān)鍵詞：SPH；二維寬頂堰溢流；數(shù)值模擬；自由面計(jì)算；并行計(jì)算

中圖分類號(hào)：TV652.1 文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：1672-1683（2017）00-0127-05

近年來SPH方法在自由表面流動(dòng)模擬中取得了較大進(jìn)展，能夠模擬傳統(tǒng)有網(wǎng)格算法較難處理的自由面翻卷破碎、液艙晃蕩、流固耦合等現(xiàn)象，得益于其拉格朗日粒子特性；SPH方法中粒子不僅用于估算場(chǎng)變與近似控制方程，也代表介質(zhì)系統(tǒng)，具有諸如密度、壓力、速度、內(nèi)能等宏觀物理量，相當(dāng)于物質(zhì)點(diǎn)。SPH方法既具有了拉格朗日方法描述物質(zhì)界面準(zhǔn)確的優(yōu)勢(shì)，又兼?zhèn)錃W拉方法容易計(jì)算的長處。因此，SPH法已用來研究漬壩、自由表面流動(dòng)、多相流等問題。

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

1SPH理論

1.1SPH方法的基本方程

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

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

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

（1）

（2）

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

（4）

（5）

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

（6）

式（6）說明了粒子i處的任一函數(shù)值可通過應(yīng)用光滑函數(shù)對(duì)其緊支域內(nèi)所有粒子相對(duì)應(yīng)的函數(shù)值進(jìn)行加權(quán)平均近似求得。

1.3控制方程

將描述流體運(yùn)動(dòng)的質(zhì)量守恒方程和動(dòng)量守恒方程（Navier-Stokes方程）表示為

（7）

（8）

（9）

（10）

2數(shù)值模擬

2.1模型建立

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

2.2結(jié)果及分析

2.2.1水流形態(tài)

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

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

2.2.2流量及流速

試驗(yàn)?zāi)Ｐ土髁客ㄟ^直角形三角薄壁量水堰來測(cè)量，公式如下：

Q=1.365H^2.48 （11）

其中：H=H₁-H₂。在與本文二維計(jì)算模型相同水頭下，試驗(yàn)?zāi)Ｐ蛿?shù)據(jù)為：Q=6.5×10^-4m³/s：?jiǎn)螌捔髁縌_單=Q/B=4.3×10^-3m³/s；下游水位超過堰項(xiàng)水深h_s=0.056 m，堰前斷面平均流速v₀=Q/A=0.19 m/s，故h_s/H₀=0.84，據(jù)此可說明發(fā)生了淹沒出流。

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

（12）

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

二維計(jì)算模型中，當(dāng)t<1.2 s，模型剛開始運(yùn)行，還未達(dá)到穩(wěn)定狀態(tài)，其流量變化幅度比較大，即圖6中N區(qū)；當(dāng)1.2 s6.9 s后，隨著模型運(yùn)行，計(jì)算模型并沒有維掙恒定水位，故計(jì)算流量有一定程度的下降，即圖6中0區(qū)。

圖7展示了堰前斷面6.75 m在t=2.66 s流速變化情況，由圖可知，在初始狀態(tài)中入口處水流速度很穩(wěn)定，所以從渠底到堰高這段距離內(nèi)，流速分布比較均勻；從堰高至水流表面，水流速度變化幅度很大。對(duì)數(shù)據(jù)進(jìn)行擬合，可得變化曲線為三次方程。由于計(jì)算粒子數(shù)目有限并且一些參數(shù)設(shè)定帶來的計(jì)算誤差等原因，并沒有得到流速分布的對(duì)數(shù)方程；不過斷面上流速的變化趨勢(shì)與對(duì)數(shù)方程很相似。圖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屬于堰前，流速變化不大，分布比較均勻，有稍微的增長趨勢(shì)；斷面7.0m至斷面7.4 m屬于堰上斷面，流速快速增長，變化趨勢(shì)為先增大后減小，最大值在斷面7.2 m上；原因是隨著堰上斷面面積變小，水流動(dòng)能增加，勢(shì)能減小，并且斷面7.2 m為收縮斷面，收縮斷面之后水深有一定的增長，水面線與堰項(xiàng)基本保持平行。斷面7.4 m至斷面7.5 m屬于堰后，流速變大，是由于水流勢(shì)能減小，動(dòng)能增加。從圖7、圖8、圖9可看出，流速存在一定的震蕩現(xiàn)象，主要原因是：（1）計(jì)算粒子數(shù)目不足，導(dǎo)致結(jié)果數(shù)據(jù)過少；（2）并沒有實(shí)現(xiàn)恒定流的入流方式，使水頭在逐漸漸小。（3）模型參數(shù)設(shè)置上還未達(dá)到最佳，還需進(jìn)一步調(diào)試。

3結(jié)語

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

在下一步的研究中，建立三維模型和解決入流問題，使模擬結(jié)果與實(shí)際更相近以及水流形態(tài)與實(shí)驗(yàn)觀察更吻合，使SPH方法在水動(dòng)力學(xué)方面得到了更廣泛的應(yīng)用。