耦合VOF的非线性k-ε模型三维溃坝数值模拟

丁伟业，金 生，艾丛芳

(大连理工大学 海岸和近海工程国家重点实验室, 大连 116024)

耦合VOF的非线性k-ε模型三维溃坝数值模拟

丁伟业，金 生，艾丛芳

(大连理工大学 海岸和近海工程国家重点实验室, 大连 116024)

溃坝是一种伴随有巨大破坏力的灾害性的水流现象。三维溃坝的数值模拟能够为溃坝水流的研究分析提供科学的依据。针对瞬间全溃模型水流的复杂性，考虑水库下游障碍物的滞水及阻水作用，利用开源程序OpenFOAM建立了耦合VOF对流方程的三维非线性(NL)k-ε紊流数学模型求解雷诺时均纳维斯托克斯方程。采用有限体积法进行离散，利用PISO算法进行数值求解。通过模拟下游为三角形障碍物、矩形障碍物以及90°弯道的溃坝实验，并将模拟结果与实验结果进行对比分析，初步证明了该模型的准确性和鲁棒性。

VOF法；NLk-ε模型；溃坝；OpenFOAM；三维数值模拟

2016年1月14日11时许，湖北省恩施龙凤镇三河村白庙突发堰塞湖溃坝事故。2016年2月3日，伊拉克摩苏尔水坝面临即将溃坝的风险，成千上万的人处在危险之中，工作人员对水坝进行加固。溃坝是指坝体溃决，是一种低频率高危害的社会致灾因素，会对下游地区人民生命财产安全造成灾难性的破坏。因此进行溃坝洪水数值计算与模拟分析具有重要的理论与现实意义[1]。研究溃坝问题的方法主要有三种：理论分析、模型试验和数值模拟。三种研究方法各有其优缺点：当处理十分简单的模型问题时理论分析方法能够得到较为精确的解析解，而对于复杂问题并不能得到其控制方程的解析解；模型试验方法是研究水流流动机理、分析水流流动现象、探讨并获得水流流动新概念的主要手段，然而模型试验受经费、时间及观测精度的限制，并且难以摆脱原型和模型之间不完全相似的困扰；数值模拟方法能够解决复杂的水流流动问题，同时不受经费和时间限制[2]。

近年来随着计算技术的不断发展，数值模拟越来越多的应用于溃坝以及溃坝水流对下游复杂地形冲刷的研究当中。虽然1-D和水深2-D溃坝水流模型能够模拟溃坝水流的演进过程，但是由于对控制方程进行了如静水压强分布假设，以及假设无明显的垂向加速度和自由表面曲率等，使得这些模型在溃坝的初始阶段、溃坝波前端和靠近结构物水流等方面均不适用了[3-4]。垂向2D和3D模型能够更精细地模拟溃坝水流的复杂流态及细节，尤其对于自由水面变化强烈的溃坝水流，3D模型所得到的丰富信息为研究局部水流特性、溃口发展机理以及水工建筑物结构设计等提供了依据[5]。对于自由水面的追踪捕捉，Lucy LB[6], Gingold和Monaghan[7]提出了利用流体粒子运动代替计算网格来求解控制方程的光滑粒子运动方法(SPH)。近些年三维非静压数值模型在自由表面的问题研究中具有较快的发展，从求解控制丰城数值方法分类主要有：显示投影法[8]；半隐分步法[9～10]；全隐式法[11]。由Hirt和Nichols[12]提出的VOF方法则更为研究者们所熟知。体积分数F介于0和1之间，当F等于0或1时，则单元内没有交界面。本文利用开源计算软件OpenFOAM研究3D溃坝流模型。利用VOF方法耦合NLk-ε紊流模型求解雷诺时均N-S方程。采用六面体网格单元建立网格模型，有限体积法进行离散，PISO算法进行求解。每一个PISO算法的内迭代步中的压力梯度通过压力泊松方程来求解，并用Rhie-Chow插值原理来消除压力波[13]。通过与实验案例进行对比验证来证明模型的可靠性。

1 数学模型

溃坝波的运动主要是受重力作用,水流和气体有明显的分界面,因此可以用水气两相流的分层模型VOF[14～16]来进行模拟自由液面。自由水面的运动通过流体体积函数的对流输运方程来捕捉，可以较为精细地描述溃坝洪水中的水面变化，克服了静压假定和刚盖假定对变化剧烈的自由水面的限制和导致的压力场失真。

VOF方法是将每个单元中一项流体体积与单元体积之比定义为流体体积函数，即

(1)

并且单元内的物理特性为

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

(2)

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

(3)

其中l与g分别为气相和液相。则流体体积分数α满足对流输运方程[17]

(4)

模型的基本微分方程包括连续性方程、动量方程、以及NLk-ε方程，分别表示如下。

连续性方程

·u=0

(5)

雷诺时均动量方程

(6)

式中：u=(u,v,w)为笛卡尔坐标系下流速；p*=p-ρg·x为修正压力；ρ=ρ(x)、μ=μ(x)由式(2)、式(3)得到；σT和kα=分别为气液交界面处的表面张力和曲率。τNL为非线性雷诺应力张量

(7)

(8)

式中：Cτ1=-4.0,Cτ2=13.0,Cτ3=-2.0和A2=1 000均为常数。η=SKε由Yakhot et al. (1992)获得，且S为流体的平均应变率，有

(9)

k方程

(10)

方程

(11)

其中

(12)

(13)

式中：αξ=0.9和A1=1.25为常数，且ξ=Ωkε，且

(14)

NLk-ε中生成项Pk同样由非线性描述表达，有

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

(15)

式中：经验常数取值为C1=1.44,C2=1.92,σε=0.77，σk=1.

图1 模型三维尺寸图Fig.1 Three-dimensional size of model

2 模型验证

2.1下游为三角形障碍物溃坝模拟

模型为长38 m、宽1 m的矩形河道。河道上游为长15.5 m，初始水深 为0.75 m的水库。障碍物位于河道起点下游25.5 m位置，障碍物与水库间为干燥的河床，障碍物以下河床水深为0.15 m，具体尺寸如图1所示。图中测点G4、G10、G13、G20分别位于河道中部距离水库下游4 m、10 m、13 m和20 m[3,18,19]。对于水平网格：障碍物上游选择3 cm，障碍物附近选择1 cm，障碍物下游选择2 cm。在宽和高方向均采用1 cm网格。

图2 水位测点实验测量值与计算模拟值对比Fig.2 Comparison of water depth between experiment and computed results at gauging points G4, G10, G13 and G20

图2为各测点测得水位实验值与计算值的对比结果。由图可知，G4点整体吻合较好，但是计算模拟的反射波与实验值相比略微滞后；G10点计算模拟的最大水位略高于测量水深；G13点计算模拟上游溃坝来流时间与实验结果相比稍有提前，且在反射波过后的t=30～34 s，计算值与测量值差异较大；计算模拟溃坝水流抵达G20点时有冲刷，因而计算模拟的水位低于实验结果。虽然各测点的模拟值与实验值均存在差异，但模拟值与实验测量值整体吻合还是较好的。

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为水位测点，p1、p3、p5和p7分别为压力测点[3,20]。采用2 cm×2 cm×1 cm网格划分模型。

图3 模型三维尺寸图Fig.3 Three-dimensional size of model

采用k-ε和NLk-ε两种紊流方程对模型2进行数值模拟。图3为H4与H2水位测点数值模拟与实验测量和K.M.T. Kleefsman(2005)[20]结果的对比。图4为t=0.5 s,0.8 s,1.4 s时的流域特征量。图4-a为流速矢量压力云图，图4-b为相应时刻的河道中心剖面体积分数图，图4-c为3D流场示图。t=0.5 s时，溃坝波抵达障碍物，水流由中心向两侧及上方绕流经过障碍物，因此在障碍物前端形成高压区，H2点测得的水深为0.08 m。t=0.8 s时，溃坝波仍然作用于障碍物上，障碍物后方形成负压区，H2点测得的水深为0.14 m。t=1.4 s时，反射波经过障碍物与上游溃坝波相遇，此时H2点测得的水深为0.22 m。

4-a压力流速场 4-b中心剖面 4-c 3D流场图4模型2在t=0.5 s, 0.8 s和1.4 s时的内部流场Fig.4 Internal flow field of model 2 at t=0.5 s, 0.8 s and 1.4 s

图5 水位测点的数值模拟与实验测量结果Fig.5 Comparison of water depth between experiment and computed results at gauging points H4 and H2

图6 模型压力测点的数值模拟与实验测量结果的对比Fig.6 Comparison of pressure between experiment and computed results of NLk-ε model

由图5可知，从整体来看k-ε和NLk-ε两种紊流模型都能反映出溃坝发生后水库及河道的水位变化。当溃坝水流遇到下游障碍物形成的反射波传递到H4点时NLk-ε模型的水位低于实验测量结果且存在瞬间的水位下降，而k-ε模型计算模拟的反射波要滞后于测量结果近0.5 s。当水流从下游边界反射回上游到达H4点时，NLk-ε模型的水位低于实验测量结果，而k-ε模型计算模拟的反射波仍然存在着滞后。H2点k-ε和NLk-ε两紊流模型数值模拟与实验测值要好于H4点。下游反射波到达H2点时，NLk-ε模型的水位高于测量值，模拟结果没有k-ε模型模拟结果好。当上游水流再次流过H2点形成涌水时，两模型均滞后于测量值。本文所提出的模型在与K.M.T. Kleefsman的结果对比当中表现出了良好的准确性。

图7 90°弯道模型尺寸及各测点位置分布Fig.7 Dimensions of physical model with 90° bend and locations of P1,P2,P3,P4,P5and P6

图6为NLk-ε模型压力测点的数值模拟与实验测量结果的对比。从图中可知，NLk-ε模型能够较好的模拟流场内的压力。但数值模拟的测点p1和p3的压力峰值均小于测量值。且当反射波流过测点时，两测点模拟值与测量值相比均存在滞后。p5和p7测点的模拟结果与测量值对比均较好。

由以上数值模拟与实验测量结果可知，NLk-ε模型能够较好的反应溃坝水流的演进过程，并且流场内的各物理特性均能反应出来。在溃坝模拟中，NLk-ε模型比k-ε模型表现的要好。

2.3下游为90°直角弯道溃坝模拟

模型上游为长×宽×高=244 cm×244 cm×63 cm的矩形水库。水库初始时刻水深h0为53 cm，在水库中安置P1号水位测点。下游为宽49.5 cm的90°直角弯道，沿河道分别布置P2～P6号水位测点。水库底部与河道底部存在33 cm的高差。模型具体尺寸以及各水位测点分布见图7。

图8 各测点水位计算值与实验值随时间变化的对比Fig.8 Comparison of water surface profiles between measured and calculated results at different gauging points

图8为各测点水深数值模拟与实验测量的结果对比，从图中可以看出，数值模拟的结果与实验测量结果整体拟合较好。在t=0 s时刻，撤去闸门，水库内水流入河道。t=0～2.5 s时间段，溃坝波进入与水库相连河道，除去P2测点水位峰值数值模拟值略低于实验测量值外，其它各水位测点均拟合的较好。t=2.5 s后溃坝水流分为两部分在河道中流动，其中一部分经河道壁面反射向水库方向传播，另一部分经90°转角沿河道向下游传播。由图8可知，数值模拟的反射波经过P4测点的时间滞后于实验测量结果，且水位峰值高于实验测量结果，从而导致了向下游传播的水量减少，这一点从P5和P6测点数值模拟的水位峰值低于实验测量结果就可以观察到。数值模拟的反射波在经过P3测点时与实验测量值有较好拟合，在经过P2测点时有明显的提前，且峰值高于实验测量值，进而导致从水库进入河道水量减少，水库内水位数值模拟的水位高于实验测量结果，这一点从P1测点可以观察到。

3 结论

本文利用开源计算软件OpenFOAM中的VOF法与NLk-ε紊流模型进行耦合来研究3-D溃坝问题。采用有限体积法进行离散，利用PISO算法进行求解。利用水库下游存在障碍物以及存在90°转角河道的三维溃坝模型来验证耦合模型的可靠性。从模拟结果中可知，模型能够较好的模拟溃坝流场的演进过程，为溃坝水流的研究工作提供了科学依据。但针对数值模拟与实验测量之间存在的水库下游水位峰值差异以及反射波形成的滞后问题均需在以后的工作做进一步的修正。本文只针对瞬间全溃模型进行了验证，然而实际的溃坝过程要复杂的多，许多其他因素对于溃坝水流的运动都会产生影响，例如：河床的地形变化、河床侵蚀、局部溃口及逐渐溃决等。耦合利用VOF法与3-D NLk-ε紊流模型对于其他条件下的溃坝水流运动的模拟能力强弱需要再进行专门的验证。

[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中的实现[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

国家自然科学基金项目(51309052; 51479022)；中央高校基本科研业务费专项基金(DUT15LK01)；辽宁省教育厅重点实验室基础研究项目(LZ2015012) .

丁伟业(1988-)，男，黑龙江省人，博士研究生，主要从事流体自由表面研究。

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

TV 122

A

1005-8443(2017)05-0489-06