耦合VOF的非線性k-ε模型三維潰壩數(shù)值模擬

丁偉業(yè)，金 生，艾叢芳

(大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室, 大連 116024)

耦合VOF的非線性k-ε模型三維潰壩數(shù)值模擬

丁偉業(yè)，金 生，艾叢芳

(大連理工大學(xué) 海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室, 大連 116024)

潰壩是一種伴隨有巨大破壞力的災(zāi)害性的水流現(xiàn)象。三維潰壩的數(shù)值模擬能夠?yàn)闈嗡鞯难芯糠治鎏峁┛茖W(xué)的依據(jù)。針對瞬間全潰模型水流的復(fù)雜性，考慮水庫下游障礙物的滯水及阻水作用，利用開源程序OpenFOAM建立了耦合VOF對流方程的三維非線性(NL)k-ε紊流數(shù)學(xué)模型求解雷諾時(shí)均納維斯托克斯方程。采用有限體積法進(jìn)行離散，利用PISO算法進(jìn)行數(shù)值求解。通過模擬下游為三角形障礙物、矩形障礙物以及90°彎道的潰壩實(shí)驗(yàn)，并將模擬結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行對比分析，初步證明了該模型的準(zhǔn)確性和魯棒性。

VOF法；NLk-ε模型；潰壩；OpenFOAM；三維數(shù)值模擬

2016年1月14日11時(shí)許，湖北省恩施龍鳳鎮(zhèn)三河村白廟突發(fā)堰塞湖潰壩事故。2016年2月3日，伊拉克摩蘇爾水壩面臨即將潰壩的風(fēng)險(xiǎn)，成千上萬的人處在危險(xiǎn)之中，工作人員對水壩進(jìn)行加固。潰壩是指壩體潰決，是一種低頻率高危害的社會致災(zāi)因素，會對下游地區(qū)人民生命財(cái)產(chǎn)安全造成災(zāi)難性的破壞。因此進(jìn)行潰壩洪水?dāng)?shù)值計(jì)算與模擬分析具有重要的理論與現(xiàn)實(shí)意義[1]。研究潰壩問題的方法主要有三種：理論分析、模型試驗(yàn)和數(shù)值模擬。三種研究方法各有其優(yōu)缺點(diǎn)：當(dāng)處理十分簡單的模型問題時(shí)理論分析方法能夠得到較為精確的解析解，而對于復(fù)雜問題并不能得到其控制方程的解析解；模型試驗(yàn)方法是研究水流流動機(jī)理、分析水流流動現(xiàn)象、探討并獲得水流流動新概念的主要手段，然而模型試驗(yàn)受經(jīng)費(fèi)、時(shí)間及觀測精度的限制，并且難以擺脫原型和模型之間不完全相似的困擾；數(shù)值模擬方法能夠解決復(fù)雜的水流流動問題，同時(shí)不受經(jīng)費(fèi)和時(shí)間限制[2]。

近年來隨著計(jì)算技術(shù)的不斷發(fā)展，數(shù)值模擬越來越多的應(yīng)用于潰壩以及潰壩水流對下游復(fù)雜地形沖刷的研究當(dāng)中。雖然1-D和水深2-D潰壩水流模型能夠模擬潰壩水流的演進(jìn)過程，但是由于對控制方程進(jìn)行了如靜水壓強(qiáng)分布假設(shè)，以及假設(shè)無明顯的垂向加速度和自由表面曲率等，使得這些模型在潰壩的初始階段、潰壩波前端和靠近結(jié)構(gòu)物水流等方面均不適用了[3-4]。垂向2D和3D模型能夠更精細(xì)地模擬潰壩水流的復(fù)雜流態(tài)及細(xì)節(jié)，尤其對于自由水面變化強(qiáng)烈的潰壩水流，3D模型所得到的豐富信息為研究局部水流特性、潰口發(fā)展機(jī)理以及水工建筑物結(jié)構(gòu)設(shè)計(jì)等提供了依據(jù)[5]。對于自由水面的追蹤捕捉，Lucy LB[6], Gingold和Monaghan[7]提出了利用流體粒子運(yùn)動代替計(jì)算網(wǎng)格來求解控制方程的光滑粒子運(yùn)動方法(SPH)。近些年三維非靜壓數(shù)值模型在自由表面的問題研究中具有較快的發(fā)展，從求解控制豐城數(shù)值方法分類主要有：顯示投影法[8]；半隱分步法[9～10]；全隱式法[11]。由Hirt和Nichols[12]提出的VOF方法則更為研究者們所熟知。體積分?jǐn)?shù)F介于0和1之間，當(dāng)F等于0或1時(shí)，則單元內(nèi)沒有交界面。本文利用開源計(jì)算軟件OpenFOAM研究3D潰壩流模型。利用VOF方法耦合NLk-ε紊流模型求解雷諾時(shí)均N-S方程。采用六面體網(wǎng)格單元建立網(wǎng)格模型，有限體積法進(jìn)行離散，PISO算法進(jìn)行求解。每一個(gè)PISO算法的內(nèi)迭代步中的壓力梯度通過壓力泊松方程來求解，并用Rhie-Chow插值原理來消除壓力波[13]。通過與實(shí)驗(yàn)案例進(jìn)行對比驗(yàn)證來證明模型的可靠性。

1 數(shù)學(xué)模型

潰壩波的運(yùn)動主要是受重力作用,水流和氣體有明顯的分界面,因此可以用水氣兩相流的分層模型VOF[14～16]來進(jìn)行模擬自由液面。自由水面的運(yùn)動通過流體體積函數(shù)的對流輸運(yùn)方程來捕捉，可以較為精細(xì)地描述潰壩洪水中的水面變化，克服了靜壓假定和剛蓋假定對變化劇烈的自由水面的限制和導(dǎo)致的壓力場失真。

VOF方法是將每個(gè)單元中一項(xiàng)流體體積與單元體積之比定義為流體體積函數(shù)，即

(1)

并且單元內(nèi)的物理特性為

ρ=ρ1α+(1-α)ρg

(2)

μ=μ1α+(1-α)μg

(3)

其中l(wèi)與g分別為氣相和液相。則流體體積分?jǐn)?shù)α滿足對流輸運(yùn)方程[17]

(4)

模型的基本微分方程包括連續(xù)性方程、動量方程、以及NLk-ε方程，分別表示如下。

連續(xù)性方程

·u=0

(5)

雷諾時(shí)均動量方程

(6)

式中：u=(u,v,w)為笛卡爾坐標(biāo)系下流速；p*=p-ρg·x為修正壓力；ρ=ρ(x)、μ=μ(x)由式(2)、式(3)得到；σT和kα=分別為氣液交界面處的表面張力和曲率。τNL為非線性雷諾應(yīng)力張量

(7)

(8)

式中：Cτ1=-4.0,Cτ2=13.0,Cτ3=-2.0和A2=1 000均為常數(shù)。η=SKε由Yakhot et al. (1992)獲得，且S為流體的平均應(yīng)變率，有

(9)

k方程

(10)

方程

(11)

其中

(12)

(13)

式中：αξ=0.9和A1=1.25為常數(shù)，且ξ=Ωkε，且

(14)

NLk-ε中生成項(xiàng)Pk同樣由非線性描述表達(dá)，有

Pk=ρ(μtS:u-τNL:u)

(15)

式中：經(jīng)驗(yàn)常數(shù)取值為C1=1.44,C2=1.92,σε=0.77，σk=1.

圖1 模型三維尺寸圖Fig.1 Three-dimensional size of model

2 模型驗(yàn)證

2.1下游為三角形障礙物潰壩模擬

模型為長38 m、寬1 m的矩形河道。河道上游為長15.5 m，初始水深 為0.75 m的水庫。障礙物位于河道起點(diǎn)下游25.5 m位置，障礙物與水庫間為干燥的河床，障礙物以下河床水深為0.15 m，具體尺寸如圖1所示。圖中測點(diǎn)G4、G10、G13、G20分別位于河道中部距離水庫下游4 m、10 m、13 m和20 m[3,18,19]。對于水平網(wǎng)格：障礙物上游選擇3 cm，障礙物附近選擇1 cm，障礙物下游選擇2 cm。在寬和高方向均采用1 cm網(wǎng)格。

圖2 水位測點(diǎn)實(shí)驗(yàn)測量值與計(jì)算模擬值對比Fig.2 Comparison of water depth between experiment and computed results at gauging points G4, G10, G13 and G20

圖2為各測點(diǎn)測得水位實(shí)驗(yàn)值與計(jì)算值的對比結(jié)果。由圖可知，G4點(diǎn)整體吻合較好，但是計(jì)算模擬的反射波與實(shí)驗(yàn)值相比略微滯后；G10點(diǎn)計(jì)算模擬的最大水位略高于測量水深；G13點(diǎn)計(jì)算模擬上游潰壩來流時(shí)間與實(shí)驗(yàn)結(jié)果相比稍有提前，且在反射波過后的t=30～34 s，計(jì)算值與測量值差異較大；計(jì)算模擬潰壩水流抵達(dá)G20點(diǎn)時(shí)有沖刷，因而計(jì)算模擬的水位低于實(shí)驗(yàn)結(jié)果。雖然各測點(diǎn)的模擬值與實(shí)驗(yàn)值均存在差異，但模擬值與實(shí)驗(yàn)測量值整體吻合還是較好的。

2.2下游為矩形障礙物潰壩模擬

模型為長3.22 m、寬1 m的矩形河道，河道上游水庫長1.23 m，水庫初始水深h0為0.55 m，下游河床干燥無水。矩形障礙物尺寸為0.403×0.161×0.161，位于水庫下游1.17 m。圖中H4和H2為水位測點(diǎn)，p1、p3、p5和p7分別為壓力測點(diǎn)[3,20]。采用2 cm×2 cm×1 cm網(wǎng)格劃分模型。

圖3 模型三維尺寸圖Fig.3 Three-dimensional size of model

采用k-ε和NLk-ε兩種紊流方程對模型2進(jìn)行數(shù)值模擬。圖3為H4與H2水位測點(diǎn)數(shù)值模擬與實(shí)驗(yàn)測量和K.M.T. Kleefsman(2005)[20]結(jié)果的對比。圖4為t=0.5 s,0.8 s,1.4 s時(shí)的流域特征量。圖4-a為流速矢量壓力云圖，圖4-b為相應(yīng)時(shí)刻的河道中心剖面體積分?jǐn)?shù)圖，圖4-c為3D流場示圖。t=0.5 s時(shí)，潰壩波抵達(dá)障礙物，水流由中心向兩側(cè)及上方繞流經(jīng)過障礙物，因此在障礙物前端形成高壓區(qū)，H2點(diǎn)測得的水深為0.08 m。t=0.8 s時(shí)，潰壩波仍然作用于障礙物上，障礙物后方形成負(fù)壓區(qū)，H2點(diǎn)測得的水深為0.14 m。t=1.4 s時(shí)，反射波經(jīng)過障礙物與上游潰壩波相遇，此時(shí)H2點(diǎn)測得的水深為0.22 m。

4-a壓力流速場 4-b中心剖面 4-c 3D流場圖4模型2在t=0.5 s, 0.8 s和1.4 s時(shí)的內(nèi)部流場Fig.4 Internal flow field of model 2 at t=0.5 s, 0.8 s and 1.4 s

圖5 水位測點(diǎn)的數(shù)值模擬與實(shí)驗(yàn)測量結(jié)果Fig.5 Comparison of water depth between experiment and computed results at gauging points H4 and H2

圖6 模型壓力測點(diǎn)的數(shù)值模擬與實(shí)驗(yàn)測量結(jié)果的對比Fig.6 Comparison of pressure between experiment and computed results of NLk-ε model

由圖5可知，從整體來看k-ε和NLk-ε兩種紊流模型都能反映出潰壩發(fā)生后水庫及河道的水位變化。當(dāng)潰壩水流遇到下游障礙物形成的反射波傳遞到H4點(diǎn)時(shí)NLk-ε模型的水位低于實(shí)驗(yàn)測量結(jié)果且存在瞬間的水位下降，而k-ε模型計(jì)算模擬的反射波要滯后于測量結(jié)果近0.5 s。當(dāng)水流從下游邊界反射回上游到達(dá)H4點(diǎn)時(shí)，NLk-ε模型的水位低于實(shí)驗(yàn)測量結(jié)果，而k-ε模型計(jì)算模擬的反射波仍然存在著滯后。H2點(diǎn)k-ε和NLk-ε兩紊流模型數(shù)值模擬與實(shí)驗(yàn)測值要好于H4點(diǎn)。下游反射波到達(dá)H2點(diǎn)時(shí)，NLk-ε模型的水位高于測量值，模擬結(jié)果沒有k-ε模型模擬結(jié)果好。當(dāng)上游水流再次流過H2點(diǎn)形成涌水時(shí)，兩模型均滯后于測量值。本文所提出的模型在與K.M.T. Kleefsman的結(jié)果對比當(dāng)中表現(xiàn)出了良好的準(zhǔn)確性。

圖7 90°彎道模型尺寸及各測點(diǎn)位置分布Fig.7 Dimensions of physical model with 90° bend and locations of P1,P2,P3,P4,P5and P6

圖6為NLk-ε模型壓力測點(diǎn)的數(shù)值模擬與實(shí)驗(yàn)測量結(jié)果的對比。從圖中可知，NLk-ε模型能夠較好的模擬流場內(nèi)的壓力。但數(shù)值模擬的測點(diǎn)p1和p3的壓力峰值均小于測量值。且當(dāng)反射波流過測點(diǎn)時(shí)，兩測點(diǎn)模擬值與測量值相比均存在滯后。p5和p7測點(diǎn)的模擬結(jié)果與測量值對比均較好。

由以上數(shù)值模擬與實(shí)驗(yàn)測量結(jié)果可知，NLk-ε模型能夠較好的反應(yīng)潰壩水流的演進(jìn)過程，并且流場內(nèi)的各物理特性均能反應(yīng)出來。在潰壩模擬中，NLk-ε模型比k-ε模型表現(xiàn)的要好。

2.3下游為90°直角彎道潰壩模擬

模型上游為長×寬×高=244 cm×244 cm×63 cm的矩形水庫。水庫初始時(shí)刻水深h0為53 cm，在水庫中安置P1號水位測點(diǎn)。下游為寬49.5 cm的90°直角彎道，沿河道分別布置P2～P6號水位測點(diǎn)。水庫底部與河道底部存在33 cm的高差。模型具體尺寸以及各水位測點(diǎn)分布見圖7。

圖8 各測點(diǎn)水位計(jì)算值與實(shí)驗(yàn)值隨時(shí)間變化的對比Fig.8 Comparison of water surface profiles between measured and calculated results at different gauging points

圖8為各測點(diǎn)水深數(shù)值模擬與實(shí)驗(yàn)測量的結(jié)果對比，從圖中可以看出，數(shù)值模擬的結(jié)果與實(shí)驗(yàn)測量結(jié)果整體擬合較好。在t=0 s時(shí)刻，撤去閘門，水庫內(nèi)水流入河道。t=0～2.5 s時(shí)間段，潰壩波進(jìn)入與水庫相連河道，除去P2測點(diǎn)水位峰值數(shù)值模擬值略低于實(shí)驗(yàn)測量值外，其它各水位測點(diǎn)均擬合的較好。t=2.5 s后潰壩水流分為兩部分在河道中流動，其中一部分經(jīng)河道壁面反射向水庫方向傳播，另一部分經(jīng)90°轉(zhuǎn)角沿河道向下游傳播。由圖8可知，數(shù)值模擬的反射波經(jīng)過P4測點(diǎn)的時(shí)間滯后于實(shí)驗(yàn)測量結(jié)果，且水位峰值高于實(shí)驗(yàn)測量結(jié)果，從而導(dǎo)致了向下游傳播的水量減少，這一點(diǎn)從P5和P6測點(diǎn)數(shù)值模擬的水位峰值低于實(shí)驗(yàn)測量結(jié)果就可以觀察到。數(shù)值模擬的反射波在經(jīng)過P3測點(diǎn)時(shí)與實(shí)驗(yàn)測量值有較好擬合，在經(jīng)過P2測點(diǎn)時(shí)有明顯的提前，且峰值高于實(shí)驗(yàn)測量值，進(jìn)而導(dǎo)致從水庫進(jìn)入河道水量減少，水庫內(nèi)水位數(shù)值模擬的水位高于實(shí)驗(yàn)測量結(jié)果，這一點(diǎn)從P1測點(diǎn)可以觀察到。

3 結(jié)論

本文利用開源計(jì)算軟件OpenFOAM中的VOF法與NLk-ε紊流模型進(jìn)行耦合來研究3-D潰壩問題。采用有限體積法進(jìn)行離散，利用PISO算法進(jìn)行求解。利用水庫下游存在障礙物以及存在90°轉(zhuǎn)角河道的三維潰壩模型來驗(yàn)證耦合模型的可靠性。從模擬結(jié)果中可知，模型能夠較好的模擬潰壩流場的演進(jìn)過程，為潰壩水流的研究工作提供了科學(xué)依據(jù)。但針對數(shù)值模擬與實(shí)驗(yàn)測量之間存在的水庫下游水位峰值差異以及反射波形成的滯后問題均需在以后的工作做進(jìn)一步的修正。本文只針對瞬間全潰模型進(jìn)行了驗(yàn)證，然而實(shí)際的潰壩過程要復(fù)雜的多，許多其他因素對于潰壩水流的運(yùn)動都會產(chǎn)生影響，例如：河床的地形變化、河床侵蝕、局部潰口及逐漸潰決等。耦合利用VOF法與3-D NLk-ε紊流模型對于其他條件下的潰壩水流運(yùn)動的模擬能力強(qiáng)弱需要再進(jìn)行專門的驗(yàn)證。

[1]WANG X L, ZHANG A L, CHEN H L. Three-dimensional numerical simulation of dam-break flood routing in the complex inundation areas[J]. Journal of Hydraulic Engineering, 2012(9): 1 025-1 033,1041.

[2]LIN J B, JIN S, DING W Y. Dam break water flow simulation based on HydroInfo software[J]. Journal of Water Resources and Architectural Engineering,2015(6):113-117.

[3]REZA M , WU W M. 3-D finite-volume model of dam-break flow over uneven beds based on VOF method[J]. Advances in Water Resources,2014, 70: 104-117.

[4]YANG C, LIN B, JIANG C. Predicting near-field dam-break flow and impact force using a 3D model[J]. Journal of Hydraulic Research,2010, 48(6):784-792.

[5]FU C J, LIAN J J. Three-dimensional numerical simulation of dam-break flow in complicated river sections[J]. Journal of Hy-draulic Engineering,2007,38(10): 151-157.

[6]LUCY L B. A numerical approach to the testing of the fission hypothesis[J]. The Astronomical Journal,1977,82:1 013-1 024.

[7]Gingold R A, Monaghan J J. Smoothed particle hydrodynamics: theory and application to non-spherical stars[J]. Monthly No-tices of the Royal Astronomical Society, 1977, 181(3):375-389.

[8]WANG K, SHENG J, GANG L. Numerical modeling of free-surface flows with bottom and surface-layer pressure treatment [J].Journal of Hydrodynamics, Ser. B, 2009, 21(3):352-359.

[9]AI C F, JIN S. Three-dimensional free-surface flow model for simulating water wave motions [J]. Journal of Hydrodynamics, Ser. A, 2008, 23(3):338-347.

[10]AI C F, JIN S, LV B. A new fully non-hydrostatic 3D free surface flow model for water wave motions[J].International Journal for Numerical methods in Fluids, 2011,66(11):1 354-1 370.

[11]YOUNG C C, WU C H, KUO J T, et al. A higher-order sigma-coordinate model for simulating non-linear refraction-diffraction of water waves[J].Costal Engineering, 2009,56(9): 919-930.

[12]Hirt C W, Nichols B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39:201-225.

[13]東岳流體.Rhie-Chow插值在OpenFOAM中的實(shí)現(xiàn)[EBOL]. http:www.dyfluid.comRhieChow.html.

[14]Torey M D, Cloutman L D, Mjilsness R C, et al. NASA-VOF2D: A computer program for incompressible flows with free sur-faces[R]. USA:Los Alamos Scientific Laboratory,1985.

[15]Rider W J, Kothe D B. Reconstructing volume tracking[J]. Journal of Computational Physics, 1998，141:112-152.

[16]WANG Z D, WANG D G. Free surface reconstruction with VOF method[J]. Journal of Hydrodynamic, 2003(1): 52-56.

[17]Rusche H. Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions[D]. London:Imperial College Lon-don (University of London), 2003.

[18]Marsooli R, Zhang M, Wu W. Vertical and Horizontal Two-Dimensional Numerical Modeling of Dam-Break Flow over Fixed Beds[J]. World Environmental and Water Resources Congress,2011，79:2 225-2 233.

[19]Zhou J G, Causon D M, Mingham C G, et al. Numerical prediction of dam-break flows in general geometries with complex bed topography[J]. Journal of Hydraulic Engineering,2004, 130(4):332-340.

[20]Kleefsman K M T, Fekken G, Veldman A E P, et al. A Volume-of-Fluid based simulation method for wave impact problems[J]. Journal of Computational Physics, 2005, 206(1):363-393.

3-D numerical simulation of dam-break flow based on VOF method with nonlineark-εturbulent model

DINGWei-ye,JINSheng,AICong-fang

(StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian116024,China)

Dam break is a kind of disastrous water phenomenon with great destructive power. 3-D numerical model can provide scientific basis for dam-break flow research and analysis. In view of the complexity of the instantaneous total dam-break flow model, considering the stagnant water and water blocking of the obstacle in the downstream of the reservoir, a 3-D NLk-εturbulence model coupled with the VOF method in OpenFOAM was established to solve the Reynolds-Averaged Navier-Stokes equations. The PISO algorithm was utilized for the pressure-velocity coupling and the equations were discretized using the finite volume method. The accuracy and robustness of the model has been tested comparing with the laboratory experiments of dam-break flow over a triangle obstacle, rectangle obstacle and 90 degree channel in this paper.

VOF；NLk-εModel; dam-break; OpenFOAM; 3-D numerical simulation

2017-03-27；

2017-05-05

國家自然科學(xué)基金項(xiàng)目(51309052; 51479022)；中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)基金(DUT15LK01)；遼寧省教育廳重點(diǎn)實(shí)驗(yàn)室基礎(chǔ)研究項(xiàng)目(LZ2015012) .

丁偉業(yè)(1988-)，男，黑龍江省人，博士研究生，主要從事流體自由表面研究。

Biography:DING Wei-ye(1988-),male,doctor student.

TV 122

A

1005-8443(2017)05-0489-06