邵家儒,楊 瑜

(重慶理工大學(xué)機(jī)械工程學(xué)院,重慶 400054)

三維洪水流動(dòng)的SPH數(shù)值模型

邵家儒,楊 瑜

(重慶理工大學(xué)機(jī)械工程學(xué)院,重慶 400054)

對(duì)光滑粒子動(dòng)力學(xué)(SPH)方法進(jìn)行了改進(jìn),提出了高精度的耦合動(dòng)力學(xué)邊界處理方法,保證粒子的無(wú)穿透條件,通過(guò)液面粒子的支持域求和準(zhǔn)確施加自由液面條件,準(zhǔn)確預(yù)測(cè)水流的翻卷、破碎效應(yīng)。建立洪水流動(dòng)的簡(jiǎn)化模型,分析防波堤對(duì)洪水的阻礙作用,準(zhǔn)確預(yù)測(cè)洪水流動(dòng)中各種三維特征。應(yīng)用SPH模型對(duì)三維圓柱形水體坍塌過(guò)程進(jìn)行模擬,結(jié)果表明,SPH數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合很好,從而驗(yàn)證了SPH模型模擬劇烈自由液面流動(dòng)的有效性。

光滑粒子動(dòng)力學(xué);自由液面;洪水;耦合動(dòng)力邊界

洪水是一種嚴(yán)重的自然災(zāi)害,根據(jù)其形成原因可以分為暴雨洪水、潰壩洪水、融雪洪水等。洪水發(fā)生時(shí),水流會(huì)向下游急速傳播,具有很強(qiáng)的突發(fā)性和危害性,嚴(yán)重影響人們的生命和財(cái)產(chǎn)安全[1-2]。洪水流動(dòng)涉及水體自由表面變化、翻卷、破碎,復(fù)雜的紊流與漩渦,以及強(qiáng)烈的液固耦合作用。數(shù)值模擬方法花費(fèi)少、可重復(fù)性好,沒(méi)有實(shí)驗(yàn)測(cè)試中可能出現(xiàn)的各種風(fēng)險(xiǎn),并且可以研究流動(dòng)中的強(qiáng)非線性特征,已經(jīng)逐漸成為研究自由液面流動(dòng)及其與結(jié)構(gòu)耦合作用問(wèn)題的重要方法[3-6]。

傳統(tǒng)基于網(wǎng)格的方法,如有限差分方法[7]、有限體積方法,已經(jīng)廣泛用于計(jì)算流體力學(xué)的研究中,并逐漸發(fā)展了一些自由表面追蹤技術(shù),如VOF、 Level set等[8-9]。然而在不規(guī)則的復(fù)雜區(qū)域上構(gòu)造規(guī)則的網(wǎng)格是比較困難的,常常需要復(fù)雜的數(shù)學(xué)變換,同時(shí)網(wǎng)格發(fā)生畸變時(shí)也會(huì)大大降低該方法的精度。近年來(lái),很多學(xué)者將研究方向轉(zhuǎn)移到光滑粒子動(dòng)力學(xué)(smoothed particle hydrodynamics,SPH)方法[10-11]上來(lái),并取得了很多成果。Shao等[12]對(duì)二維潰壩問(wèn)題進(jìn)行了數(shù)值模擬,應(yīng)用密度修正和核梯度修正消除了流場(chǎng)中的數(shù)值震蕩,得到了光滑的壓力場(chǎng)。Marrone等[13]應(yīng)用邊界虛粒子法對(duì)潰壩相關(guān)問(wèn)題進(jìn)行了研究;楊秀峰等[14]應(yīng)用修正的SPH方法在水波破碎方面進(jìn)行了相關(guān)研究;Liu等[15]提出了一種改進(jìn)的耦合動(dòng)力邊界條件,提高了潰壩流固體壁面上的插值精度,比較準(zhǔn)確地預(yù)測(cè)了壓力響應(yīng)。已有的工作大多針對(duì)二維模型展開研究,本文針對(duì)兩種典型的三維潰壩情況進(jìn)行數(shù)值模擬,以準(zhǔn)確預(yù)測(cè)潰壩流動(dòng)及有無(wú)防波堤對(duì)潰壩的影響。

1 SPH 方法

1.1 SPH 控制方程

潰壩流動(dòng)問(wèn)題主要由N-S方程控制,其具體形式如下:

(1)

(2)

式中:ρ——密度;t——時(shí)間;v——速度;P——壓強(qiáng);g——重力加速度;μ——?jiǎng)恿︷ば韵禂?shù)。為了將該方程應(yīng)用于數(shù)值模擬,需對(duì)其進(jìn)行數(shù)值離散。SPH方程的離散過(guò)程分為兩步,即核近似和粒子近似,離散后的SPH形式的控制方程描述如下:

(3)

(4)

其中vij=vi-vjxij=xi-xj

式中:m——粒子質(zhì)量;W——光滑函數(shù)[16];h——光滑長(zhǎng)度;x——粒子位置;r——粒子間的距離;i、j——粒子編號(hào)。光滑函數(shù)不僅決定函數(shù)近似形式、定義粒子支持域的尺寸,而且還決定核近似和粒子近似的緊致性和精度。本文模擬應(yīng)用了3次樣條函數(shù),其對(duì)應(yīng)的支持域半徑為2倍光滑長(zhǎng)度,具有較好的緊致性。其表達(dá)式如下:

(5)

式中:R——粒子間距與光滑長(zhǎng)度之比。

1.2 核梯度

(6)

(7)

在每個(gè)時(shí)間步粒子搜索結(jié)束后,只需進(jìn)行一次核梯度相關(guān)矩陣的計(jì)算就可以將結(jié)果覆蓋于所有的場(chǎng)變量,從而提高所有導(dǎo)數(shù)的求解精度。

每個(gè)詩(shī)節(jié)充滿了跨詩(shī)行連續(xù)的句子，這也是詩(shī)人精心的編排。這種跨詩(shī)行連續(xù)句的使用，能夠產(chǎn)生驅(qū)動(dòng)效果。不完整句子的詩(shī)行自然吸引讀者去補(bǔ)全缺失，尋找完整的句子和意義，讀者在驅(qū)動(dòng)下順著詩(shī)人的指引去找回消逝的童年和童年記憶中的美好生活。另外，跨行的連續(xù)句在詩(shī)歌中占有一定的比重，有助于營(yíng)造一幅連續(xù)動(dòng)態(tài)伊甸園般的田園景象，如同在讀者面前緩緩展開畫卷——那里有屋旁的蘋果樹，幽谷上的星空，成群的馬車，田野里的雛菊大麥和閃耀著蘋果光輝的小河……這種連續(xù)的句子把大自然的美一一展現(xiàn)，帶給讀者無(wú)盡的遐想。

1.3 改進(jìn)的耦合動(dòng)力學(xué)邊界

本文提出一種改進(jìn)的耦合動(dòng)力學(xué)邊界處理方法,如圖1所示。固壁邊界由兩類固定的虛粒子構(gòu)成,第一類虛粒子稱為排斥力粒子,均勻分布在固壁邊界位置,它對(duì)鄰近邊界的粒子作用排斥力,阻止粒子的非物理穿透;第二類虛粒子均勻分布在邊界外,其層數(shù)由所應(yīng)用核函數(shù)的支持域半徑來(lái)確定,從而保證真實(shí)粒子的支持域內(nèi)不會(huì)出現(xiàn)粒子缺失。

圖1 改進(jìn)的耦合動(dòng)力學(xué)邊界粒子分布示意圖Fig.1 Schematics of particle distributions in the modified coupled dynamic boundary

對(duì)于邊界粒子的壓力和速度,應(yīng)用如下形式進(jìn)行計(jì)算：

(8)

(9)

該計(jì)算形式可以恢復(fù)核函數(shù)的歸一化性質(zhì),提高邊界粒子相關(guān)物理量的計(jì)算精度,同時(shí)可以保證邊界上的無(wú)滑移條件。

1.4 自由液面判定

自由液面附近的粒子支持域上會(huì)存在一定的粒子缺失,根據(jù)該特征可以根據(jù)粒子支持域上核函數(shù)插值公式建立判定條件,即

(10)

判定為自由液面粒子后則可以設(shè)定此類壓力為零,其密度為初始參考密度。

2 數(shù) 值 算 例

2.1 圓柱形水體坍塌過(guò)程數(shù)值模擬

為了驗(yàn)證三維SPH模型的有效性,本文首先對(duì)圓柱形水體在圓柱形容器中的坍塌過(guò)程進(jìn)行了數(shù)值模擬。模型數(shù)據(jù)來(lái)源于Maschek等[17]的試驗(yàn),如圖2所示,Dc=0.11 m,Hc=0.2 m,Dv=0.44 m。模擬過(guò)程中初始粒子間距為0.002 5 m,流體密度為1 000 kg/m3,流體使用了約12萬(wàn)個(gè)粒子來(lái)代替,圓柱形剛體壁面由3層粒子構(gòu)造,使用了約20萬(wàn)個(gè)粒子。

水體坍塌過(guò)程如圖3所示,初始時(shí)刻水體保持靜止,在重力的作用下下端水體沿著底部壁面向周圍蔓延(如0.1 s時(shí)刻),并逐漸鋪滿整個(gè)底面(如0.2 s時(shí)刻)。之后底部水體與圓柱形壁面發(fā)生撞擊,產(chǎn)生較大的沖擊,沿著剛體壁面向上運(yùn)動(dòng)(如0.3 s時(shí)刻),當(dāng)水體達(dá)到一定高度后,動(dòng)能衰減為零,開始下落回流(如0.4 s時(shí)刻),并在容器底部聚集(如0.6 s時(shí)刻),此時(shí)容器中心處的水體受到擠壓,從底部的形心處向上流動(dòng),形成噴射狀的水柱。通過(guò)與Maschek等[17]的試驗(yàn)結(jié)果對(duì)比,可以看到,SPH模擬結(jié)果與試驗(yàn)結(jié)果吻合較好,證明改進(jìn)的SPH方法可以很好地模擬三維自由表面液體流動(dòng)問(wèn)題。

圖3 SPH結(jié)果與試驗(yàn)結(jié)果對(duì)比Fig.3 Comparison between SPH and experimental results

2.2 洪水流動(dòng)對(duì)建筑的沖擊

近年來(lái),由強(qiáng)暴雨導(dǎo)致的潰壩及潰堤事件頻發(fā),潰壩及潰堤洪水給人們的生命財(cái)產(chǎn)構(gòu)成了巨大威脅,尤其對(duì)位于洪流演進(jìn)路線上及周邊的城鎮(zhèn)危害最大。筆者針對(duì)水流對(duì)城區(qū)建筑物的沖擊問(wèn)題建立了簡(jiǎn)化模型,并進(jìn)行了數(shù)值模擬。

復(fù)雜街區(qū)模型如圖4所示,街區(qū)有兩排結(jié)構(gòu)相同、建筑高度相同的長(zhǎng)方體形構(gòu)成。街區(qū)位于一個(gè)長(zhǎng)18 m、寬5 m的長(zhǎng)方體空間內(nèi),建筑物的尺寸均為5.0 m×1.0 m×1.0 m,街區(qū)(第一排建筑物)距離水壩的距離為7.5 m,同一排建筑物間相距0.5 m,第一排與第二排建筑物間相距1 m,左側(cè)水體長(zhǎng)4.5 m、寬5 m、高2.5 m,初始時(shí)刻保持靜止。模型粒子間距為0.1 m,其中水體粒子約59 000個(gè),防波堤由3層粒子組成,粒子數(shù)約為1 600個(gè),長(zhǎng)方體空間的邊界為剛體壁面,粒子數(shù)約為58 000個(gè)。

圖4 洪水流動(dòng)模型示意圖Fig.4 Schematics of the flooding flow model

圖5 水流沖擊過(guò)程及速度場(chǎng)變化示意圖Fig.5 The impact process and velocity field of the flood flow

從圖5可以看到,受到復(fù)雜街區(qū)影響的洪水水流呈現(xiàn)出復(fù)雜的水流特征,如水躍、水波渦旋、水流間斷等,這些特性具有很強(qiáng)的三維特性(u為流體的水平速度)。洪水水流在每個(gè)建筑的后方都形成了湖區(qū),同時(shí)在建筑物的阻礙作用下水波發(fā)生破碎,速度方向發(fā)生改變,反向流動(dòng)速度超過(guò)2 m/s(如3.0 s時(shí)刻)。受復(fù)雜街區(qū)影響,洪水水流演進(jìn)到建筑物后方的空地上后形成了很多間斷水渦。上述水流特征都是典型的三維流態(tài),可見三維SPH數(shù)值模型能比較合理地模擬復(fù)雜的水流形態(tài)(圖6)。

圖6 洪水流動(dòng)俯視圖Fig.6 Top view of the flood flow

2.3 防波堤對(duì)洪水水流的阻礙作用

為了研究防波堤對(duì)水流運(yùn)動(dòng)產(chǎn)生的影響,本文在距離左側(cè)8 m(與水體左側(cè)壁面間的距離為3.5 m)處設(shè)置一個(gè)高1.0 m的防波堤(圖4中的黑色區(qū)域)。在防波堤作用下,洪水呈現(xiàn)出了更為復(fù)雜的水流特征。如圖7所示,水流與防波堤撞擊之后會(huì)沿著防波堤向上運(yùn)動(dòng),此過(guò)程中水平流動(dòng)速度大大降低(如1.2 s時(shí)刻)。水體越過(guò)防波堤后回落,這一過(guò)程中會(huì)在防波堤的左側(cè)出現(xiàn)巨大的空腔,相互作用過(guò)程中,水波從破碎狀態(tài)逐漸恢復(fù)為較穩(wěn)定的流動(dòng)狀態(tài),直到與建筑物相遇。

圖7 有防波堤時(shí)的水流沖擊過(guò)程及速度場(chǎng)變化示意圖Fig.7 The impact process and velocity field of the flood flow with a bulwark

由圖5～8可知,防波堤對(duì)洪水的運(yùn)動(dòng)有巨大的阻礙作用。如T=1.8 s時(shí)刻,在無(wú)防波堤作用時(shí),洪水已經(jīng)蔓延到了第一街區(qū),并向第二街區(qū)繼續(xù)運(yùn)動(dòng),而有防波堤阻礙時(shí)洪水剛剛越過(guò)防波堤。從圖8中T=3.0 s及T=4.2 s時(shí)的速度分布可以看到,有防波堤作用時(shí)水流的流量和流速明顯降低,這表明防波堤對(duì)洪水起到很大的截流作用,街區(qū)中的水位變低,這將大大降低城市排水系統(tǒng)的排水壓力。

圖8 有防波堤時(shí)的洪水流動(dòng)俯視圖Fig.8 Top view of the flood flow with a bulwark

3 結(jié) 語(yǔ)

本文建立了洪水流動(dòng)問(wèn)題的三維SPH數(shù)值模型?；谔├占?jí)數(shù)展開的核梯度計(jì)算可以提高函數(shù)導(dǎo)數(shù)的計(jì)算精度,改進(jìn)的耦合動(dòng)力邊界處理能夠保證邊界粒子的規(guī)則分別及無(wú)穿透條件,結(jié)合自由液面條件可以準(zhǔn)確地預(yù)測(cè)水流的翻卷、破碎等效應(yīng)。圓柱體水流坍塌問(wèn)題中模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合良好,驗(yàn)證了改進(jìn)的SPH方法的有效性。有無(wú)防波堤作用下的洪水流動(dòng)模擬表明,SPH模型可以準(zhǔn)確預(yù)測(cè)洪水流動(dòng)中的各種復(fù)雜三維特征,對(duì)于工程問(wèn)題具有一定的指導(dǎo)作用。

[ 1 ] 王立輝,胡四一.潰壩問(wèn)題研究綜述.水利水電科技進(jìn)展[J].2007,27(1):80-85.(WANG Lihui,HU Siyi.Study on dam failure-related problems[J].Advances in Science and Technology of Water Resources,2007,27(1):80-85.(in Chinese))

[ 2 ] 周克發(fā),李雷.基于社會(huì)經(jīng)濟(jì)發(fā)展的潰壩洪水損失動(dòng)態(tài)預(yù)測(cè)評(píng)價(jià)模型[J].長(zhǎng)江流域資源與環(huán)境,2015(增刊1):147-150.(ZHOU Kefa,LI Lei.Dynamic forecasting evaluation model of flood loss due to dam breach based on socioeconomic development[J].Resources and Environment in the Yangtze Basin,2015(Sup1): 147-150.(in Chinese))

[ 3 ] 王軍,梁忠民,施曄.基于GIS的水庫(kù)洪水風(fēng)險(xiǎn)圖編制[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,38(1):20-25.(WANG JUN,LIANG Zhongmin,SHI Ye.Mapping of flood risk of reservoirs using GIS technology[J].Journal of Hohai University(Natural Sciences),2010,38(1):20-25(in Chinese))

[ 4 ] 鄧鵬,李致家.3種水文模型在淮河息縣流域洪水模擬中的比較[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,41(5):377-382.(DENG Peng,LI Zhijia.Comparison of three hydrological models in flood simulation for Xixian Basin of Huaihe River[J].Journal of Hohai University(Natural Sciences),2013,41(5): 377-382.(in Chinese))

[5] 包紅軍,王莉莉,李致家,等.基于Holtan產(chǎn)流的分布式水文模型[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2016,44(4):340-346.( BAO Hongjun,WANG Lili,LI Zhijia,et al.A distributed hydrological model based on Holtan runoff generation theory[J].Journal of Hohai University(Natural Sciences), 2016,44(4):340-346.(in Chinese))

[6] 劉佳嘉,周祖昊,賈仰文,等.分布式水文模型子流域編碼方法對(duì)比分析[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2017,45(1):22-29.(LIU Jiajia,ZHOU,Zuhao,JIA Yangwen,et al.Comparison analysis of subwatershed codificationmethods for distributed hydrological model[J].Journal of Hohai University(Natural Sciences), 2017,45(1):22-29.(in Chinese))

[ 7 ] THOMAS J W.Numerical partial differential equations: finite difference methods[M].New York:Springer Science amp; Business Media,2013.

[ 8 ] RAEINIA Q,BLUNT M J,BIJELJIC B.Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method[J].Journal of Computational Physics,2012,231(17): 5653-5668.

[ 9 ] DOYEUX V,GUYOT Y,CHABANNES.Simulation of two-fluid flows using a finite element/level set method[J].Journal of Computational and Applied Mathematics,2013,246: 251-259.

[10 ] MONAGHAN J J.Smoothed particle hydrodynamics and its diverse applications[J].Annual Review of Fluid Mechanics,2012,44: 323-346.

[11] LIU Moubin,LIU Guirong.Smoothed particle hydrodynamics(SPH):an overview and recent developments[J].Archives of Computational Methods in Engineering,2010,17(1): 25-76.

[12] SHAO Jiaru,LIU Moubin,YANG Xiufeng,et al.Improved smoothed particle hydrodynamics with RANS for free surface flow problem[J].International Journal of Computational Methods,2012,9(1): 1240001.

[13] MARRONE S,ANTUONO M,COLAGROSSI.δ-SPH model for simulating violent impact flows[J].Computer Methods in Applied Mechanics and Engineering,2011,200(13): 1526-1542.

[14] 楊秀峰,彭世镠.潰壩流的光滑粒子法模擬[J].計(jì)算物理,2010,272(2):173-180.(YANG Xiufeng,PENG Shiliu.Simulation of dam-break flow with SPH method[J].Chinese Journal of Computational Physics,2010,272(2):173-180.(in Chinese))

[15] LIU Moubin,SHAO Jiaru,CHANG Jianzhong.On the treatment of solid boundary in smoothed particle hydrodynamics[J].Science China Technological Sciences,2012,55(1): 244-254.

[16] BONET J,KULASEGARAM S.Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations[J].International Journal for Numerical Methods in Engineering,2000,47(6): 1189-1214.

[17] MASCHEK W,MUNZ C D,MEYER L.Investigations of sloshing fluid motions in pools related to recriticalities in liquid-metal fast breeder reactor core meltdown accidents[J].Nuclear Technology,1992,98(1): 27-43.

Numericalstudyof3DfloodingflowusingSPHmethod

SHAOJiaru,YANGYu

(CollegeofMechanicalEngineering,ChongqingUniversityofTechnology,Chongqing400054,China)

In this paper, the SPH method is modified by incorporating a new coupled dynamic solid boundary with high calculation accuracy, in which non-penetrating condition is maintained within the particles and the free surface conditions are applied by using the sum of the support regions of corresponding particles. This method is expected to provide an accurate estimation of the rolling and broken effects of the water flow. A simplified flooding flow model is thus established to examine the effect of bulwark against the flood, in which various characteristics of the 3D flooding flow are well predicted. The collapse behavior of the cylindrical body of water is simulated using a 3D SPH model. The finding shows that the numerical results agree well with the experimental results, which demonstrates the effectiveness of the modified SPH method in modelling the free surface with violent flow behavior.

smoothed particle hydrodynamics (SPH), free surface, flood flow, coupled dynamic solid boundary

10.3876/j.issn.1000-1980.2017.06.007

2016-11-29

國(guó)家自然科學(xué)基金(11602045);重慶市基礎(chǔ)科學(xué)與前沿技術(shù)研究項(xiàng)目(cstc2016jcyjA0373);重慶市教育委員會(huì)科學(xué)技術(shù)研究項(xiàng)目(KJ1600918)

邵家儒(1986—),男,山東淄博人,副教授,博士,主要從事光滑粒子動(dòng)力學(xué)和流體力學(xué)研究。E-mail: shaojiaru@cqut.edu.cn

TV131.2

A

1000-1980(2017)06-0515-07