Numerical Study of 2D and 3D Separation Phenomena in the Dam-Break Flow Interacting with a Triangular Obstacle

Dam-break turbulent flow interacting with obstacles is simulated with the VOF method implemented in an in-house unstructured-grid finite-volume Navier-Stokes code. A special attention is paid to prediction of separation phenomena using low-Re computational grids that provide full resolution of viscous sublayers on the bottom and side confining walls, if any. Some original developments aimed at improvement of the VOF method robustness for such kind of flows are presented. The test case considered is interaction of the dam-break induced water stream with a triangular obstacle. Computations under conditions of experiments by Soares-Frazao (2007) have been carried out on the base of 2D and 3D formulations. It is shown that action of the bottom wall friction leads to formation of one or two separation “bubbles”, depending on the flow development phase, and to occurrence of associated hills at the free surface, which are observed in experimental photos as well. Taking into account presence of side walls of the experimental channel results in solutions with a considerably 3D shape of the computed free surface, and its side view much better agrees with the experimental photos than that given by 2D solutions. Moreover, local-in-time separation of the flow from the side walls is predicted with the 3D formulation.


Introduction
Dam-break flows or other impact single-wave flows may cause damage of civil and industrial buildings and constructions.For the sake of the damage mitigation or minimization such kind of flows must be well studied including all the aspects of their interaction with different obstacles.
To date experimental data is available for several model dambreak flows that interact with obstacles having different shape (triangular - [1,2], trapezoidal - [3,4], vertical wall - [5,6] and others).As a rule, experimental works provide a set of photos showing free surface configuration for different time instants.
Currently experimental studies are being displaced more and more by numerical simulation applicable to a wide variety of flows.Different numerical models have different computational costs and areas of applicability.According to this trend, experiment is increasingly used as a data source for validation of mathematical models and computational tools.Sure, this statement is completely applicable to the problems attributed to the dam-break flows.
Traditionally numerical modeling of a single wave propagation and its interaction with different obstacles is performed using a method of shallow water equations (SWE) (see for example [1,[7][8][9][10]), which provides a proper prediction of unsteady areas of water raise before the obstacles (positions of these areas and their height).However, the SWE method is not capable to reproduce details of wave-obstacle interaction, such as overturning (breaking) of a negative wave, which appears when a stream is fully or partly reflected by an obstacle, and fails to predict accurately pressure distribution over walls of the obstacle (especially in case of bluff obstacles).As well, as pointed in some of the aforementioned contributions, the SWE method does not predict accurately the propagation speed of the main intensive wave and reflected waves, and does not reproduce details of their free surface shape.
Alternative to the SWE approach is solving the full threedimensional (3D) or two-dimensional (2D) Navier-Stokes (NS) equations, or the Reynolds Averaged Navier-Stokes equations (RANS), if turbulence is considered as a significant factor.Capabilities of the SWE and NS approaches are compared, in particular, in [3] and [11].Computational results obtained via solving the Navier-Stokes equations for dam-break flow interacting with different obstacles can be found in [2-5, 12, 13].
Separation zone may occur at front side of an obstacle.Such a zone can be seen, in particular, in the 2D velocity field computed with code Flow-3D for a model dam-break flow interacting with a trapezoidal obstacle [3].Apparently, appearance of this zone is caused by an adverse pressure gradient occurring when the flow goes up at a front side of the obstacle.
Application of numerical methods for analysis of such a complicated flow like the intensive single-wave interaction with an obstacle requires a justified assurance that no scheme factors are affecting considerably a solution.This requirement is related to both the schemes used for prediction of free surface evolution and to the approaches employed for near-wall layer treatment.
Typically, the flows under consideration are characterized by high Reynolds numbers, and are often (in most cases) modeled using one or another turbulence model.Here it should be noted that the majority of the previous computations of free surface flow developing over a wall were performed with application of the standard wall function technique (using socalled high-Re computational grids) to satisfy no-slip condition on the solid wall or even with prescribing slip condition.It is well known, however, that the standard wall functions do not perform well in case of boundary layer separation caused by adverse pressure gradient.
Nowadays the most popular approach for free surface tracking is the Volume-of-Fluid (VOF) method [14].It allows performing computations in cases of strong deformations of free surface including the case of wave breaking.In this method free surface position is determined by space distribution of so-called markerfunction presenting volume fraction of fluid.Marker-function evolution is governed by a convective transport equation, and the quality of numerical schemes used for approximation of this equation affects directly the free surface artificial smearing or/ and deformation.
It is well known that application of conventional schemes for convective flux evaluation is not suitable for solution of the marker-function equation as it leads to smearing of a transitional area, where the fluid volume fraction must vary rapidly from zero to unity, over plenty computational cells.There are several specialized (so-called "compressive") numerical schemes proposed in literature for approximation of this equation, for example HRIC [15], CICSAM [16] and M-CICSAM scheme [17].
The present work, extending studies reported in [18], is aimed to development and application of 3D numerical techniques for simulation of a dam-break turbulent flow interacting with an obstacle.Some original developments aimed at improvement of the VOF method robustness for such kind of flows are presented.A special attention is paid to prediction of separation phenomena using low-Re computational grids that provide full resolution of viscous sublayers on the bottom and confining side walls.Results of 2D and 3D computations performed for a test configuration with a triangular obstacle are presented and discussed.

Computational method
Present developments and computations are based on using an in-house unstructured-grid finite-volume Navier-Stokes code called Flag-S.In this code the VOF method is used for free surface tracking, and the SIMPLEC algorithm is used for pressure-velocity coupling in computations of incompressible fluid motion.Note also that all computational results presented below were obtained without taking into account effects of surface tension on the gas-liquid interface.

VOF method
As mentioned in Introduction, in the framework of the VOF method [14] free surface position is determined by spatial distribution of so-called marker-function C, which in fact is volume fraction of fluid in a computational grid cell: С = 1cell contains only liquid, С = 0 -cell contains only gas.It is assumed that the free surface coincides with an isosurface C = 0.5.Herewith, liquid and gas can be treated as a single fluid having variable material properties defined as follows (this approach is called one fluid formulation): Governing equations for this effective fluid motion are solved throughout the entire computational domain and no boundary conditions at the gas-liquid interface are needed.
To apply the finite-volume discretization method, the governing equations are written in a conservative form.Conservative form of the momentum equation is given by: The equation governing convective transport of the fluid volume fraction can be also transformed to a conservative form (4) using continuity equation (5) applicable for incompressible liquid and gas.
According to the finite-volume technique, discretized forms of these equations are as follows (summation is performed over all faces of a control volume being a computational cell): Different values of parameter β correspond to different time-discretization schemes.Note that form (6) ensures conservation of the marker-function total amount in the computational domain, i.e. ensures liquid and gas global conservation.
As seen from ( 6), values of the marker-function at computational cell faces, C f , are required.They should be evaluated via interpolation of C-values from neighboring cells.It is known that conventional upwind schemes most commonly used for discretization of convective terms (in momentum equation, for example) are inappropriate for Eq. ( 6) as they lead to smearing of the gas-liquid interface over many computational cells.
In the literature, several specialized (so-called "compressive") numerical schemes are proposed for approximation of C f .A detailed comparative study of performance of two popular schemes HRIC [15] and CICSAM [16], and also a promising modification of CICSAM scheme, M-CICSAM [17] was conducted in work [19].The latter one has demonstrated its superiority over the other schemes examined: less dependence on quality of computational grid and time step size.According to these findings, the M-CICSAM scheme was chosen for the computations presented below.It was found also that the Crank-Nicolson scheme (corresponding to β = 0.5) is preferable when solving Eq. ( 6) [19].
In case of taking into account turbulence effects, the set of above given governing equations is added by transport equations of turbulence parameters that define the eddy viscosity, and the latter is simply added to the molecular one.For the present computations, the Menter's two-equation SST turbulence model was used, since typically this model shows superiority over other eddy-viscosity turbulence models in case of separation flow computations [20].Transport equations of this model were solved throughout the computational domain, i.e. with no boundary conditions at the gas-liquid interface.

Approximation of convective part of the momentum equation
Convective flux of momentum at a cell face is evaluated via application of an interpolation procedure using density and velocity values from neighbouring cells.At that, a specific issue arises when treating the area corresponding to the gas-liquid interface since fluid density in this area changes rapidly by several orders.As a result of this rapid change, density approximation method has a dramatic effect on the cell-face momentum flux evaluated.As pointed in [21], writing momentum equation in the above given conservative form (3) implies that at the stage of their discretization the density interpolation procedure must provide those density values at cell faces that (together with the applied time approximation scheme) ensure fulfilment of discretized conservative form of non-stationary continuity equation for the effective fluid.This requirement may be satisfied by using the same time approximation for the momentum equation and for Eq. ( 4), and evaluating cell-face density values directly using C f values that are computed at solving Eq. ( 6) [21].
However, the above described approach prevents implementation of special numerical "tricks" that can improve effectiveness and robustness of the VOF method and quality of the computed marker-function field.One of such "tricks", aimed at extra sharpening of the gas-liquid interface, is reported in [22].Another motivation for rejection of the approach elaborated in [21] is due to difficulties with its implementation in the framework of fractional step strategy, which is rather efficient and implies that several time steps are done at solving Eq. ( 6) within one time step at solving the fluid dynamics governing equations.
To avoid the necessity of usage of fully consistent discrete approximations for Eqs.(3) and ( 4) one could, in particular, employ the momentum equation written in the conventional non-conservative form (8a), as it was done, for instance, in [23], or in the "partially" conservative formulation (8b) used in [24]: Discretized forms of these equations given by (9a) and (9b) illustrate that no density face values are needed: These two variants were tried at implementation of the VOF method in code Flag-S (note that for time discretization of momentum equation the Crank-Nicolson scheme was used for all the computations presented below).It has been established, however, that both schemes (9) result in dramatic false deformation of gas-liquid interface even in relatively simple cases.As an example, Fig. 1 depicts results of 2D test computations carried out with schemes (9a) and (9b) for the gravity-induced free falling of a circular liquid cylinder, diameter 0.04 m, surrounded by air (the latter remains at rest far from the cylinder).
Obviously, that at very short times (when the cylinder velocity is small and liquid-air interaction is insignificant) the shape of the cylinder has to remain practically the same as at the initial instant.Figure 1 shows, however, that both schemes (9) produce large deformations of the cylinder shape even for a time interval of 0.07 s that is sufficiently small (the cylinder passes less than one diameter).Figure 1 shows that even in this case large deformations of the cylinder shape are produced by both schemes (9).
In the present work, it is suggested to rewrite the momentum equation in a novel form given by (10) and to use its discrete analog (11).This form provides a proper reduction of errors originated from numerical violation of mass balance in computational cells.
For the test case with falling liquid cylinder, results of computations preformed with scheme (11) are also shown in Fig. 1.One can see a radical improvement of the prediction quality.Note that a form analogous to (11) was used for discretization of the convective terms in the turbulence model transport equations.

Approximation of pressure gradient
The finite-volume discretization for pressure gradient in momentum Eq. ( 3) is given by: Here again pressure values p f at cell face centers should be evaluated via interpolation of pressure values from neighboring cells.Unfortunately, conventional linear interpolation technique (13) works incorrectly in case of resting liquid and gas, which is characterized by discontinuity in pressure gradient.For the sake of clarity, consider a case where a mesh cell face coincides with a free surface (see Fig. 2).In this case, usage of correct pressure values at centers of neighboring cells 1 and 2 will produce an overestimated pressure value at the face.As a result, for cells 1 and 2 application of (12) gives a pressure gradient value that is not balanced with the gravity force.Consequently, momentum equation will not be satisfied in the case of zero fluid velocities (as it should be), and non-physical, increasing in time local oscillations of velocity and pressure fields will arise.It can be easily obtained that a correct face pressure value is generated when using a density-weighted interpolation given by: It has been established, however, that usage of scheme (14) in case of liquid and gas motion leads to occurrence of intensive even-odd pressure oscillations near the gas-liquid interface, in particular, in the above presented case of a free falling liquid cylinder.To overcome this issue, a combination of schemes ( 13) and ( 14) is proposed where the weight ξ is dependent on the angle, α, between the normal of the current cell face and the gravity vector:

, cos
Test computations performed in the present work with scheme (15) have shown that it works correctly in both the cases of resting and free falling liquid.

Problems of accurate resolution of viscous nearbottom layer in case of liquid spreading along a dry wall
As mentioned above, free surface flows, and in particular dam-break flows, are most commonly modeled using high-Re computational grids and the standard wall function technique for determination of wall friction.However, accurate prediction of flows with occurrence of near-bottom separation zones requires accurate resolution of near-wall viscous layer that implies using a low-Re grid strongly clustered near the wall (with normalized distance, Y + , of the first computational point to the wall less or about unity).First our attempts to use such a grid for simulation of dam-break flow developing along a dry bottom have highlighted a serious issue.Consider the simplest 2D dam-break test case, where after sudden removal of the retaining wall the fluid floods the dry horizontal wall due to gravity (see Fig. 3a).The computations with a low-Re grid have shown that the air initially located in computational cells in vicinity of the dry bottom can not be properly displaced by the spreading fluid because velocity values in these cells are very low.As a result, a thin elongated non-physical air layer occurs at the bottom (Fig. 3b).
To make sure that this artefact is not caused by any mistakes in code Flag-S, a computational run for this test case was also performed with the commercial CFD code ANSYS Fluent-14.0using the same grid, and results are shown in Fig. 3c.One can see that in the solution obtained with Fluent a non-physical near-bottom air layer is observed as well (the reason of spatial oscillations seen in the Fluent-solver solution remains unclear).Obviously that this non-physical air layer leads to a radical under-estimation of wall friction.As a result, boundary layer separation phenomenon can not be predicted accurately.So a special numerical technique should be introduced to overcome this issue.In the present work, an artificial diffusion term acting in a thin near-wall region was added to the marker-function transport equation as follows: Here χ max and d cell are user-defined parameters, and d cell defines distance from the wall within which the diffusion term is active.Typically, it can be set to one half of computational cell size in the flow core.With such a choice, d cell is sufficiently large to cover the area of the non-physical effect under consideration and at the same time it is sufficiently small to avoid noticeable influence on the flow field prediction.χ max is a diffusion coefficient.For the present computations it was set to 0.0001(gH 3 ) 1/2 , where H is the initial height of the water.This relatively small value was sufficient to completely avoid occurrence of the non-physical air layer.

Results of computations for a triangular obstacle 3.1 2D computations
Experimental study of dam-break water flow interaction with a triangular obstacle was conducted in [1] under geometrical and initial conditions defined in Fig. 4. In the experiments, the flow was treated as nominally two-dimensional.Confining side walls were transparent, and the free-surface shape at various time instants was fixed with a photo camera.Basic two-dimensional computations were performed with no-slip conditions on the initially dry bottom of the channel and the obstacle surface.The computational grid had 25 000 quadrangle cells and was clustered near the wall so that Y + -values were less than unity.Additional run was done using slip conditions, both on the bottom and the obstacle surface.For this run, the grid had 16 000 cells, with no clustering near the wall.For both the runs, the SST turbulence model [20] was employed.
Free surface shapes computed with two types of wall boundary conditions are compared in Fig. 5.The computational results are superposed onto experimental photos given in [1] for two time instants counted from the instant of sudden removal of retaining partition.One can see that the free surface shape predicted with the slip condition (Fig. 5a,c) disagrees with experimental observations dramatically.An accurate resolution of near-wall viscous effects results in a considerable improvement of the flow prediction quality.The most important difference between the solutions is attributed to manifestation of separation phenomena in the case of noslip conditions.Figure 5 shows that firstly one (at t = 3.0 s) and then two (t = 3.7 s) relatively large separation zones occur in front of and at the obstacle.Occurrence of these zones leads to formation of «hills» on the free surface, visible in the experimental photos as well.
(16) (17) In more details flow patterns predicted in case of no-slip conditions are illustrated in Fig. 6 where momentum vector fields for the same time instants are given (momentum vectors in every third computational cell along both the directions are displayed).It can be seen that (i) near-wall velocity gradients are well resolved, (ii) reverse-flow zones characterized also by intensive vortex motion are quite large in size, (iii) within the left separation zone, at t = 3.7 s, there is an additional smaller-size vortex adjacent to the bottom.
Two more runs with accurate resolution of viscous sublayer were performed reducing the molecular viscosity by 10 and 100 times, and keeping the sizes of the channel as in the experiments [1].Note that after rescaling of the problem according to the similarity theory, one can easily conclude that these runs are equivalent to cases of water flow developing after the dam breaks with the initial water levels of 55 cm and 240 cm correspondingly.Despite a considerable increase in the Reynolds number, in both the solutions the separation effects remain significant, and the separation zone size is reduced less than two times even in the case of the largest Reynolds number.

3D computations
Three-dimensional computations were performed to analyse the effect of the side walls confining the experimental channel in the spanwise direction (see Fig. 4).Due to the symmetry of the experimental configuration with respect to the middle plane, the computational domain covers only a half of the channel.The computational grid had one million hexahedral cells and was clustered both near the bottom and near the side wall, with the same evaluation of Y + -value as in the 2D runs.As previously, the SST turbulence model was used to introduce the eddy viscosity effects.
Figure 7 illustrates 3D free surface shape computed for the same time instants as in the experiments and 2D computations.As one can see in the plot, a considerable deviation of developing waves from a 2D form takes place in the area covering more that one third of the channel.3D effects are even more pronounced when considering the separation zones.For instance, vorticity lines originated from the middle-plane core of the separation zone(s) deviate from a straight line (that would be in a 2D flow) in the major part of the flow domain.
For the time instant of t = 3.0 s, Fig. 8 shows patterns of computed fluid velocity vectors projected on the 3D free surface and on the symmetry plane.One can see that a rather large recirculation occurs by this time near the side wall, namely at the area where the water starts to interact with the obstacle.This "local-in-time" recirculation zone evolves with time considerably.
For the same time instant, Fig. 9 partially illustrates wall friction distribution over the channel bottom, namely over the part that includes the obstacle front side (between lines A-A and B-B) and the adjoined horizontal segment of the bottom.The map given in Fig. 9a shows that in the region of nearly two-dimensional flow the wall friction streamwise component, τ x w , reaches extreme negative values twice: first near the front line of the obstacle and second at (approximately) one third of the obstacle front side.The second extreme is of higher absolute value and corresponds to position of the core of the separation vortex depicted in Fig. 6a.
The wall friction vector pattern given in Fig. 9b illustrates three-dimensional separation occurring in the region adjoined to the front edge point (being the point of intersection of three surfaces: sidewall, obstacle front side and horizontal bottom).Note also that in this region the wall friction reaches a maximum absolute value.
Figure 10 presents a comparison of side view of the computed 3D free surface with the experimental photos from [1].Due to 3D shape of the free surface, this view is not a line as in the 2D case, but looks like a complicated band.Notably that such a "band" is seen in the experimental photos as well: at each photo one can clearly see a dark area along the free surface.Comparing Figs. 5 and 10, one can conclude that a much better agreement with the experiments is achieved on the base of the 3D formulation taking into account viscosity effects near the side walls.

Conclusion
On the base of the VOF method, 3D computational techniques have been developed providing accurate treatment of both the free surface evolution and viscous near-bottom layer in the flow developing after a dam break and interacting with obstacles.Computations under conditions of experiments with a triangular obstacle have been carried out on the base of 2D and 3D formulations.Action of the bottom wall friction leads to formation of one or two separation "bubbles", depending on the flow development phase, and to occurrence of associated hill-like waves at the free surface.Taking into account presence of the experimental channel side walls gives a solution with a considerably 3D shape of the computed free surface, and its side view much better agrees with the experimental photos than that given by 2D solutions.As well, local-in-time separation of the flow from the side walls is predicted with the 3D formulation.

Fig. 2
Fig. 2 Scheme of pressure distribution in the vicinity of gas-liquid interface in case of resting liquid and gas, and interpolated pressure values at a mesh cell face coinciding with the interface.

Fig. 3
Fig. 3 Occurrence of non-physical near-bottom air layer in computations of water flow spreading along a dry wall: (a) general view, (b,c) view of nearbottom zone (stretched in the vertical direction) from results obtained with (b) code Flag-S and (c) ANSYS Fluent.

Fig. 5
Fig. 5 Effect of bottom friction on the free-surface shape computed: (a,c)slip condition, (b,d) -no-slip (separation zone is shown as well).Simulation results are superposed on the experimental photos given in [1] for (a,b) t = 3.0 s and (c,d) t = 3.7 s.

Fig. 6
Fig. 6 Near-obstacle momentum vector patterns predicted in case of no-slip condition for (a) t = 3.0 s and (b) t = 3.7 s.Free-surface shapes are shown by solid lines.

Fig. 7
Fig. 7 3D free surface shape computed at (a) t = 3.0 s and (b) t = 3.7 s for flow over the triangular obstacle under conditions of experiments [1].Vorticity lines originated from the middle-plane core of the separation zone are shown as well (dashed lines).

Fig. 8
Fig. 8 Patterns of computed fluid velocity vectors on the 3D free surface and in the symmetry plane, t = 3.0 s.

Fig. 9
Fig. 9 Wall friction distribution over the obstacle front side (bounded by lines A-A and B-B) and over adjoining part of the channel horizontal bottom: (a) streamwise component map computed at t = 3.0 s, (b) wall friction vector pattern in the vicinity of the front edge point.

Fig. 10
Fig. 10 Comparison of side views of the 3D free surface computed at (a) t = 3.0 s and (b) t = 3.7 s with experimental photos from [1].