«上一篇
文章快速检索     高级检索
下一篇»
  应用科技  2018, Vol. 45 Issue (5): 6-10, 15  DOI: 10.11991/yykj.201712007
0

引用本文  

胡健, 刘立超. 二维溃坝的数值模拟研究[J]. 应用科技, 2018, 45(5), 6-10, 15. DOI: 10.11991/yykj.201712007.
HU Jian, LIU Lichao. Numerical simulation of 2D dam-break flow[J]. Applied Science and Technology, 2018, 45(5), 6-10, 15. DOI: 10.11991/yykj.201712007.

基金项目

国家自然科学基金项目(51678045,51579052,11102048)

通信作者

刘立超,E-mail:liulichao2015@126.com

作者简介

胡健(1979−),男,副教授,博士;
刘立超(1994−),男,硕士研究生

文章历史

收稿日期:2017-12-13
网络出版日期:2018-02-23
二维溃坝的数值模拟研究
胡健, 刘立超    
哈尔滨工程大学 船舶工程学院,黑龙江 哈尔滨 150001
摘要:为了研究溃坝问题中流体的运动和自由液面的演变,针对简化的块状水坍塌进行了数值模拟。其中,流体控制方程通过有限体积法进行离散,方程的对流项和扩散项分别采用了一阶迎风格式和中心差分格式;对压力与速度耦合的离散方程以IDEAL算法进行求解;以流体体积分数法和PLIC方法进行自由液面捕捉和重构;对两相流问题中的表面张力通过连续表面力(CSF)模型进行转化。基于上,分别对有限空间内的二维溃坝和带障碍物溃坝问题进行了数值模拟,对不同时刻的自由液面进行了捕捉并与实验结果进行了对比,证明了所采用方法的精确性。
关键词二维溃坝    障碍物    N-S方程    有限体积法    流体体积法    PLIC法    IDEAL算法    连续表面力模型    
Numerical simulation of 2D dam-break flow
HU Jian, LIU Lichao    
College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China
Abstract: In order to solve fluid movement and evolution of the free surface of liquid in the dam-break problem, this paper presents a numerical simulation on the simplified block water collapse. The simulation uses the Finite Volume Method (FVM) to discretize the fluid control equation, applies the first order upwind scheme and center difference scheme to the convective term and diffusive term of the equation respectively; solves the discrete equations of the coupling pressure and speed by IDEAL; captures and reconstructs the free surface of liquid by the Volume of Fluid (VOF) method and the procedural language for integrity constraints (PLIC) method; and transforms the surface tension in the two-phase flow problem by the continuous surface force (CSF) model. On this basis, the two-dimensional dam-breaking problem and the problem of dam-breaking with barrier in a finite space were simulated respectively, and the free surfaces of liquid at several moments were captured and compared with the experiment result, proving the accuracy of the methods mentioned above.
Keywords: 2-dimensional dam-break    obstacle    N-S equation    finite volume method    volume of fluid    procedural language for integrity constraints    IDEAL    continuous surface force    

溃坝的发生往往会给社会及人民造成巨大的灾难,对于溃坝问题的研究一直在水利工程中有着重要的地位。而另一方面,溃坝流动是带自由液面的基本问题,此类问题因初始条件及边界条件简单,却能充分反映出翻卷、入水、破碎等复杂的自由液面现象而引起众多流体力学人士的关注。

Deborah M. Greaves[1]基于VOF方法,以自适应网格对方柱水的坍塌及带障碍物坍塌过程进行了模拟,并进行了试验对比。Abdolmaleki等[2]分别基于VOF方法和FLUENT软件对二维溃坝问题进行了计算。Jaswant等[3]基于一种显式中心迎风格式进行了一维和二维的溃坝模拟计算以及带障碍物溃坝的计算。Brufau等[4]基于Delaunay非结构网格对二维浅水方程进行求解,对二维溃坝问题进行了模拟分析。

本文通过求解N-S方程得出溃坝问题中流体的流态,通过VOF方法追踪自由液面位置,分析了溃坝问题中的液体流动、自由液面的演变和溃坝造成的冲击力,得到了与实验结果高度吻合的结果,证明了本文方法的精确性。

1 理论基础 1.1 控制方程

由流体的连续性方程可知,封闭计算域内流体的质量只能由通过边界的流体流动改变,其数学表达式可以写为

$\frac{{\partial \rho }}{{\partial t }} + \nabla \cdot (\rho {{V}}) = 0$

式中: $\rho $ 表示流体的密度, ${{V}} = (u,v,w)$ 代表速度。对于不可压缩流体,可以对连续方程做如下简化:

${\bf {\nabla}} \cdot {{V}} = 0$

流体的动量守恒方程也称作N-S方程,该定律的表述为,微元体中的动量对时间的变化率等于外界作用在该微元体上的各种力之和。这实际上是牛顿第二定律。

N-S方程有以下形式:

$\frac{{\partial \left( {\rho {{V}}} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho {{V}} \otimes {{V}}} \right) = \nabla \cdot \left( {\mu \nabla \otimes {{V}}} \right) - \nabla P + {{g}} + {{{f}}_\sigma } + {{S}}$ (1)

式中: $\mu $ 为液–气的混合动力黏性系数, $P$ 为压力, ${{g}} = (0,0, - g)$ 为重力矢量, ${{{f}}_{{\sigma }}} = \sigma \kappa $ 为表面张力, ${{S}} = ({S_x},{S_y},{S_z})$ 为方程源项[5]。式(1)左侧第2项称为对流项,右侧第1项称为扩散项。

1.2 控制方程离散及求解

采用有限体积法在交错网格之上离散控制方程,交错网格中,控制体的速度存储在网格边界上,其他标量存储在网格中心,这种变量存储方式有效避免了棋盘状压力分布。

对不可压连续方程在控制体上积分并应用高斯定理,得到连续方程的离散形式:

$\sum\limits_{nf = 1}^n {{{{u}}_f} \cdot {{{A}}_f} = 0} $

式中: ${{{u}}_f}$ 指界面处的速度向量, ${{{A}}_f}$ 指界面的面积矢量。

对于N-S方程,在控制体上积分之后,同样对对流项和扩散项使用高斯积分将体积分变换为面积分,而后分别对其采用一阶迎风格式和中心差分进行离散。

最终得到离散的N-S方程:

$\frac{{{{V}}_P^{t + \delta t} - {{V}}_P^t}}{{\delta t}}{V_P} + \sum\limits_{f = 1}^n {{\rho _f}{F_f}{{V}}_P^{t + \delta t}} = \sum\limits_{F = 1}^N {{\nu _f}{A_f} \cdot \left( {\nabla \otimes {{V}}} \right)_f^{t + \delta t} + {{{S}}_{{{{V}}_P}}}{V_P}} $

式中下标pf分别代表交错网格的中心和边界。

N-S方程的求解是通过IDEAL[6]算法实现的,其全称为内部双迭代高效算法。IDEAL算法基于semi-implicit method for pressure-linked equations (SIMPLE)开发而成。它通过对压力方程进行2次内部迭代计算,几乎完全舍弃了SIMPLE算法的2个基本假设,速度和压力得到了充分耦合,大大提高了算法的收敛性。

1.3 VOF和PLIC方法

VOF方法是基于流体体积分数 $F$ 这一概念提出的。流体体积分数为标量函数,定义为控制体积内特征流体所占体积与控制体体积的比值。以水为例,若控制体内没有水, $F$ 为0;若控制体被水充满, $F$ 为1;若有 $0 < F < 1$ ,则表示该控制体内有自由液面存在[7]

根据流体体积分数的定义,我们可以得出混合密度和混合动力黏性系数的表达式:

$\left\{\begin{aligned}& \rho = F{\rho _L} + (1 - F){\rho _G}\\& \mu = F{\mu _L} + (1 - F){\mu _G}\end{aligned}\right.$

流体体积分数的输运方程可以写为

$\frac{{\partial F}}{{\partial t}} + {{V}} \cdot \nabla F = 0$

由于VOF方法追踪的是网格中的流体体积,而不是追踪流体质点的运动,采用的施主–受主方法不仅要考虑本网格的流体体积分数,而且要用到相邻网格上的流体体积分数,以此来提高分辨率。因而具有容易实现、计算量小和精度高等优点,可以处理类似溃坝水流的自由面折叠的强间断问题[8]

在VOF方法中,在每个网格中定义流体体积分数 $F$ ,根据 $F$ 的值来重构自由液面的位置。PLIC方法将自由液面用网格内的一条斜线近似,这种方法较用垂直或水平线来近似自由液面的simple line interface calculation(SLIC)方法能更精确地表示出自由液面的形状,界面重构示例如图1所示[5]

Download:
图 1 PLIC界面重构的4种基本情况

由中心网格及其相邻的8个网格的流体体积分数,可确定自由液面的法向量[9-10],网格内用来近似自由液面位置的斜线就可以由此法向量和网格内的流体体积分数两者唯一确定。

1.4 表面张力的CSF模型

互不掺混的流体界面上一定存在表面张力,单位长度的表面张力称为表面张力系数,记为 $\sigma $ 。对于互不相溶流体来说 $\sigma > 0$ ;对于互溶的流体来说 $\sigma < 0$

两相流体间的这种张力可用Young-Laplace方程来表示:

$\Delta P = {P_I} - {P_O} = \sigma \kappa $

式中: $\sigma $ 为表面张力系数, $\kappa $ 为液面曲率,P为压强,下标IO分别代表曲面的凹形一侧和凸形侧。结构网格中曲率 $\kappa $ 可由自身及周围8个网格的流体体积分数确定[11-12]

2 数值模型

本文对二维溃坝的模拟采用文献[13]中的模型,示例见图2。图中,尺寸为0.6 m×1.2 m的水柱被一块挡板固定在长3.22 m,宽1.8 m的水箱左侧,在t=0时刻迅速抽去挡板,水柱开始向右坍塌。以此模型模拟溃坝水流的情形,观察过程中自由液面的演变。

Download:
图 2 二维溃坝模型

基于正交结构网格进行计算中网格数量为161×90。箱体四周均为不可滑移壁面边界,压力初场由静压条件给出,并指定壁面上的初始压力为0。计算中,取水的密度 ${\rho _W} = 998.2\; {\rm{kg/}}{{\rm{m}}^{\rm{3}}}$ ,水的动力黏性系数 ${\mu _W} = 1.003 \times {10^{ - 3}}\; {\rm{N}} \cdot {\rm{S/}}{{\rm{m}}^{\rm{2}}}$ ;空气密度 ${\rho _A} = 1.225\; {\rm{kg/}}{{\rm{m}}^{\rm{3}}}$ ,空气的动力黏性系数 $ {\mu _A} ={\rm{1}}{\rm{.789}} {\rm{\,\,4}}\times $ $ {\rm{1}}{{\rm{0}}^{ - 5}} \;{\rm{N}} \cdot {\rm{S/}}{{\rm{m}}^{\rm{2}}}$

带障碍物的二维溃坝模拟采用文献[14]中的模型,示意图如图3所示。计算域为0.584 m×0.584 m的方舱,尺寸为0.294 m×0.1461 m的水柱被挡在方舱的左侧,在横坐标0.31 m右侧有高0.048 m、宽0.024 m的障碍物,在t=0时刻迅速抽去挡板,观察过程中自由液面的变化。在此部分数值模拟中用到的常数与上一部分一致。箱体四周均为不可滑移壁面边界,压力初场由静压条件给出,并指定壁面上的初始压力为0。障碍物表面设置为无滑移壁面,同时为保证障碍物的速度为0,应采用以下方法:1)在每一层次迭代前令障碍物内速度为0;2)令障碍物速度离散方程主元系数为一极大值,以保证速度预估值为0;3)令障碍物速度修正值计算式系数为一极小值,以保证速度修正值为0[15]

Download:
图 3 带障碍物溃坝模型
3 数值结果

本文模拟结果与文献[13]的模拟结果对比如图4(左侧一列为文献[13]的结果,右侧一列为模拟结果,下同)。

Download:
图 4 自由液面演变与文献[13]的对比

文献[13]的结果是基于光滑粒子流体动力学方法(SPH) 方法得到的。在t=1.19 s左右,液面由于重力开始翻卷下落;t=1.415 s时,本文的模拟结果更能突显出由于重力作用驱使的水头向下的趋势,而这一差别在t=1.824 s时得到了更明显的体现。

对于带障碍物的二维溃坝模拟,自由液面的演变及与文献[14]中的实验结果对比见图5

Download:
图 5 带障碍物二维溃坝自由液面演变

对比试验结果和文中的计算结果,由于抽去挡板时对水柱造成了扰动,在0.1 s时实验结果水柱侧面有明显的不光滑,而在数值模拟中不会有这种现象。在0.2 s和0.3 s时,水舌的出现和甩出的水滴都得到了很好的模拟,然而数值模拟中并没有出现严重的液体飞溅。0.4 s时,由于数值模拟中水箱高度要大于试验,故模拟中水舌接触到侧壁后沿侧壁向上爬升而不是如实验中那样接触到顶壁然后翻卷,但基本的自由液面形状是一致的。

4 结论

本文以方形水柱的自由坍塌模拟了二维溃坝和带障碍物溃坝问题,可得到以下结论:1)文中对自由液面的数值模拟与实验结果有很高的一致性,证明了所采用方法的精确度;2)本文对溃坝水流的预估模拟可以为现实中防灾减灾措施提供参考;3)文中以二维近似三维,将程序拓展到三维进行更接近现实的模拟预测,或者尝试对可移动边界或障碍物进行处理都是以后可参考的研究方向。

参考文献
[1] GREAVES D M. Simulation of viscous water column collapse using adapting hierarchical grids[J]. International journal for numerical methods in fluids, 2006, 50(6): 693-711. DOI:10.1002/(ISSN)1097-0363 (0)
[2] ABDOLMALEKI K, THIAGARAJAN K P, MORRIS-THMOAS M T. Simulation of the dam break peoblem and impact flows using a navier-stokes solver[C]//15th Australasian Fluid Mechanics Conference. Sydney, Australia, 2004. (0)
[3] SINGH J, ALTINAKAR M S, DING Yan. Two-dimensional numerical modeling of dam-break flows over natural terrain using a central explicit scheme[J]. Advances in water resources, 2011, 34(10): 1366-1375. DOI:10.1016/j.advwatres.2011.07.007 (0)
[4] BRUFAU P, GARCIA-NAVARRO P. Two-dimensional dam break flow simulation[J]. International journal for numerical methods in fluids, 2000, 33(1): 35-57. DOI:10.1002/(ISSN)1097-0363 (0)
[5] 陈朝晖. 固定方舱破损进水数值模拟[D]. 哈尔滨: 哈尔滨工程大学. 2015. (0)
[6] 孙东亮, 屈治国, 何雅玲, 陶文铨. 求解流动与传热问题的一种高效稳定的分离式算法—IDEAL[J]. 工程热物理学报, 2009, 30(08): 1369-1372. DOI:10.3321/j.issn:0253-231X.2009.08.031 (0)
[7] SEIFOLLAHI M, SHIRANI E, ASHGRIZ N. An improved method for calculation of interface pressure force in PLIC-VOF methods[J]. European journal of mechanics-B/fluids, 2008, 27(1): 1-23. DOI:10.1016/j.euromechflu.2007.01.002 (0)
[8] 史宏达, 刘臻. 溃坝水流数值模拟研究进展[J]. 水科学进展, 2006(01): 129-135. DOI:10.3321/j.issn:1001-6791.2006.01.021 (0)
[9] 赵大勇, 李维仲. VOF方法中几种界面重构技术的比较[J]. 热科学与技术, 2003(04): 318-323. DOI:10.3969/j.issn.1671-8097.2003.04.008 (0)
[10] 王宗鹏. 基于有限体积法的液滴撞击固体平壁数值模拟[D]. 大连: 大连理工大学, 2010. (0)
[11] MEIER M, YADIGAROGLU G, SMITH B L. A novel technique for including surface tension in PLIC-VOF methods[J]. European journal of mechanics-B/Fluids, 2002, 21(1): 61-73. DOI:10.1016/S0997-7546(01)01161-X (0)
[12] MARKUS M. Numerical and experimental study of large steam-air bubbles injected in a water pool[D]. Zurich: ETH Zürich, 1999. (0)
[13] GAO Zhiliang, VASSALOS D, GAO Qiuxin. Numerical simulation of water flooding into a damaged vessel’s compartment by the volume of fluid method[J]. Ocean engineering, 2010, 37(16): 1428-1442. DOI:10.1016/j.oceaneng.2010.07.010 (0)
[14] HÄNSCH S, LUCAS D, HÖHNE T, et al. Application of a new concept for multi-scale interfacial structures to the dam-break case with an obstacle[J]. Nuclear engineering and design, 2014, 279: 171-181. DOI:10.1016/j.nucengdes.2014.02.006 (0)
[15] 陶文铨. 数值传热学[M]. 2版. 西安: 西安交通大学出版社, 2001: 244-245. (0)