Level Set method-based two-dimensional numerical model for simulation of nonuniform open-channel flow

Summary

The capture precision of the free surface of an open-channel with a water-air interface directly affects the calculation precision of flow field characteristics and general characteristics of the flow. Significant research effort has been devoted to Level Set since its creation, although the relevant research is mainly limited to bubble or droplet movement. In this paper, Level Set method is applied to a two-dimensional numerical simulation of open-channel turbulence, while a new numerical model is proposed and multispot synchronized experimental data are applied to the validation of numerical model. In addition, the model is used to study the flow field characteristics and general characteristics of open-channel flow, which have a water-level lowering curve. The study shows that (1) a semilogarithm zone of vertical distribution of longitudinal velocity is still present amid the transition of flow from nonuniform to uniform, and the depth-averaged velocity and wall shear stress increase along the flowing path. (2) both the energy loss coefficient and roughness coefficient of the flow at nonuniform flow region are greater than the respective values at uniform flow region, and the magnitude of the deviation is relevant to the magnitude of the flow deviation from uniform flow stage.

Keywords:

Simulation and modeling – Motion – Velocity – Shear stresses – Turbulence – Fluid flow – Flow field – Kinetics

Introduction

The flow in natural river course and artificial open-channel belongs to open-channel flow, which is characterized by the existence of free surface (water-air interface) [1,2]. Both the flood water-level in water resources and lowest navigable water-level in navigation are characteristic quantity representing free surface position. They are important parameters for measuring flood discharge capacity and navigation capacity of rivers in reality[3,4].

Research regarding open-channel flow through experimental research and numerical calculations has been completed by many researchers. In terms of experimental research, Nezu and Sanjou [5] adopted a laser Doppler anemometer (LDA) to study uniform open-channel flow turbulence, concluding that only near-wall flow velocity could be described using a logarithmic law, whereas a systematic deviation occurs between the far-wall mainstream flow velocity distribution that can be described using the logarithmic law and Coles wake flow law; T. Song and Y. M. Chiew [6] conducted a nonuniform flow experiment in a long water channel to study the impact of nonuniform flow on water movement and concluded that there is a relationship between the distribution of the parameters such as the flow velocity and turbulence intensity and the uniformity of the flow; A. H. Cardoso et al. [7] studied the distribution law of flow velocity and wall shear stress in accelerated flow. In terms of numerical calculations, K. Shiono et al. [8] studied flow in prism open-channel with complicated cross-section and deduced a depth-averaged velocity formula and wall shear stress transverse distribution equation; Arun Kamath et al. [9] established a numerical model based on the Reynolds Averaged Navier-Stokes (RANS) equations to simulate the complex surface flows over weirs around the bridge piers; Krishna Chandran et al. [10]studied the dam break flow, two-dimensional cavity filling, etc. by the volume of fluid (VOF) method to demonstrate the validity of the surface flow model; A. Khosronejad et al. [11] carried out the study of open-channel flow with a wall-mounted bridge abutment and compared the numerical simulation results of rigid-lid and level set method.

Three core issues exist in numerical simulation of open-channel turbulence, i.e., the choice of the turbulence model, discretization of the governing equation and treatment of the free surface [12]. A good turbulence model is particularly important for complicated boundary flow. But for a flow with simple boundary, such as the flow in a long straight open-channel, even the common k-ε model can basically reflect its average flow field characteristics. For the discretization of the governing equation, the finite difference method (FDM) or finite control volume is currently used, and many research results are available with relatively mature study methods.

In relative terms, the free surface treatment of open-channel flow remains to be improved in numerical simulations [13,14]. Rigid-lid and VOF are the two common methods of water surface handling at present [15,16,17]. The former is adopted to address a free surface problem of steady uniform-flow, but some researchers [11] pointed out that the reliable results may not be obtained by rigid-lid under certain conditions. The latter is used in mainstream business computing software for performing computations of the free surface of a nonuniform flow [18,19]. So we’ll choose the VOF method as the comparison object for Level Set method. In comparison, Level Set has the following advantages: (1) for the Level Set method, the free surface can be directly obtained by its definition, but for the VOF, additional geometric reconstruction must be carried out, and the results are different with different reconstruction strategies [20,21]; (2) most of the free surfaces obtained by Level Set are clear and sharp, while the results obtained by VOF may be blurred and serrated[22,23]; (3) it is easy to get a complete image by Level Set method, but the using of VOF may lead to numerical crushing, resulting in scattered images[24,25].

There are also cases in which Level Set is adopted to study the free surface and the formation of bubbles in water and during the disengagement process [24,26,27], and it has been successfully used in numerical simulations of general free surface flow[28,29,30]; however, there are few numerical studies regarding open-channel turbulence that adopt Level Set. Moreover, it is particularly important to analyze parameters such as wall shear stress, roughness coefficient and mechanical energy loss coefficient while studying the open-channel flow evolution from nonuniform flow to uniform flow. Relevant research has not been extensively conducted due to the limit of many conditions.

Therefore, on the basis of the existing research work, this paper has carried out the following original work: (1) Level Set method was recommended to capture the free surface of open-channel flow; (2) a vertical two-dimensional numerical model based on Level Set method was established; (3) an experiment of open channel flow was carried out in a long straight tank, and the experimental results were obtained by, which were employed to test the numerical model; (4) the numerical results were used to study the variation tendency of hydraulic characteristics, during the transition of open-channel flow from non-uniform to uniform.

It’s worth noting that there are many different types of flows in the open-channel flow, such as backwater flow, dropdown flow, critical flow, etc. We chose one of the most typical flow, which have a water-level lowering curve, as the research object according to the experimental conditions.

Mathematical model

Free surface equation

Level Set was first proposed by Osher and others in 1988 [20]. The basic idea behind Level Set is to translate m-dimensional curve plane motion evolution to m+1-dimensional function projection and implicitly solve the high-dimensional function, thus precisely describing the low-dimensional curve plane topological structure.

As shown in Fig 1 above, the open-channel profile of the two-dimensional flow is zoned airflow V₁, water flow V₂ and interface S(x,t) between V₁ and V₂. The constructor φ(x,y,t) makes the free surface S(x,t) the exact zero isoline of function φ(x,y,t) at any moment, and the starter of φ(x,y,t) meets the conditions for normal direction monotony in the vicinity of free surface S(x, t). It is apparent that the requirements can be met when φ(x,y,t) is taken as the symbolic distance between point (x,y) and interface S(x,t), i.e.,

There is always φ(x, y, t) = 0 at any moment for a point on the free surface, leading to dφdt=0. The following equation is obtained thereafter:

Motion equation

The calculation zones Ω consist of airflow zone V₁ and water flow zone V₂. The RANS equation is used to describe the variation in the amount of the RANS flow characteristic in the main zones, and the standard k-ε model is used to compute the Reynolds stress. The corresponding equations are as follows:

Continuity equation:

Momentum equations:

Turbulent kinetic energy (k) equation:

Turbulent kinetic energy dissipation rate (ε) equation:

Numerical method

A numerical simulation calculation usually comprises grid division, discretization of the equation, boundary condition settings, and the equation solution. They are elaborated in this section.

The rectangular structure grid is used to divide the calculated field in the calculation.

Discretization of the motion equation

The motion equation may be commonly expressed by the following pattern:

This motion equation(RANS) is discretized by using a finite volume method(FVM)[32], and the SIMPLE algorithm[33] is used to solve the problem of pressure-velocity coupling. We use the second-order QUICK scheme for the convective terms in the RANS equation and the k and ε equations[34]. A central differencing scheme is used to discretize the viscous terms. The detailed discretization of RANS equation are given in S1 Appendix.

Discretization of level set equation

It is not a simple task to solve the Level Set equation directly, because the Level Set function is implicit. Therefore, some efficient measures must be taken to deal with the equation. In reference [35], the Eq (2) is rewritten as the spatial derivative operator below:

We use different schemes to discretize the left and right sides of the rewriting equation above as follows:

At the left end of the equation, third-order Runge-Kutta is used for time discretization:

At the right end of the equation, the fifth-order WENO combined structured template is used to directly compute the values of φ_x and φ_y, obtaining the L(φ) value.

For example, using the x direction, determine two derivatives from the left end and right end for any point i₀, i.e., the left derivative φ- x and right derivative φ+ x. To solve the left derivative, as shown in Fig 2, choose six node templates {φi0−3, φi0−2, φi0−1, φi0, φi0+1, φi0+2}, and obtain the five first-order Newton difference quotients, N₁, N₂, N₃, N₄, N₅.

To obtain the right derivative, choose another six node templates {φi0−2, φi0−1, φi0, φi0+1, φi0+2, φi0+3}, and similarly obtain the five first-order Newton difference quotients, N₁, N₂, N₃, N₄, N₅.

Following the WENO method, the convective term discretization pattern through weighting is obtained:

With the upwind scheme concept, φ_x is equal to φ- x when u_i0>0 or φ+ x when u_i0<0 for node i₀. The calculation in the y direction is consistent with the x direction, obtaining (φ± y) in the same manner(S2 Appendix).

Boundary conditions

As shown in Fig 1 above, the boundaries of the calculation fields include the inlet boundary, outlet boundary, pressure boundary and wall boundary.

(I) Inlet boundary

The inlet boundary gives the flow velocity at the intake and assumes a uniform distribution, i.e., U = U₀, on the basis of available conditions.

(II) Outlet boundary

The pressure boundary is used as the outlet boundary. The pressure is static and therefore used as a constant. The atmospheric pressure is used as the pressure boundary in the calculation since this boundary is open to the atmosphere, i.e., P = P_atm.

(III) Wall boundary conditions

There are no penetration and no slip conditions applied on the wall boundary, i.e., u = 0, v = 0. The wall function is applied in the calculation for the wall side and vicinity, and the flow velocity satisfies u+=1κln(Ey+) and y+=Δyp(Cμ1/4kp1/2)μ, where u⁺,y⁺ is the dimensionless parameters, Δy_p is the distance of the near-wall node P to the solid surface, Von Karman’s constant κ = 0.41, wall roughness parameter E = 9.8. See reference [37] for the selection of the relevant parameters and calculations.

Main numerical step

The principal calculation program consists of two modules, i.e., the motion equation solution module and the free surface solution module. They are mobilized for the circulation calculation after the initial assignment for the resulting output based on the condition of convergence. The calculation procedure is given in Fig 3.

• Initial assignment: set the boundary condition and assign starters for the entire field.

• Motion equation solution:

  • Solve the motion equation coefficient and motion equation.

  • Solve the k and ε equation.

  • Obtain the pressure intensity correction equation coefficient, solve the corrected value, and update the pressure intensity and flow velocity.

  • Repeat steps (b)-(d) until all parameters are converged, and then proceed to the next module.

• Free surface equation solution:

  • Obtain φ_t with reference to the obtained flow field (u, v) and last moment φ_t-1.

  • Reinitialize φ_t, keep the symbolic characteristics, and obtain output result.

  • Determine whether to proceed to next stage calculation or end it based on the time calculation.

Key issue handling

The following issues are addressed in order to guarantee precision in the calculations:

(a) the water flow-air density ratio ρ_w/ρ_a = 828.4 and viscosity coefficient ratio μ_w/μ_a = 55.5 may destabilize the calculation given the high ρ-μ difference on both sides of the free surface. Therefore, the physical property parameters in the vicinity of the interface must be smoothly treated. Consequently, the Heaviside function H is defined as follows:

(b) After a certain time increment, function φ(x, y) no longer meets Eq (1) defined as the symbolic distance characteristics in the numerical calculation. In this regard, φ(x, y) at the end of each stage is initialized to meet the requirement, a solution to the following initial value problem:

The computational details of the reinitialization are given in S3 Appendix.

Experimental study and validation of numerical model

Experimental equipment and methods

The experiment is conducted in a 28-m long, 0.56-m wide and 0.7-m high sloped channel. The sides and bottom of the channel are made of glass. A laser water-level gauge is used to measure the water-level, and an acoustic Doppler velocity meter (ADV) is used to measure the flow velocity.

In the experiment, the channel bottom slope i₀ is taken to be 0.002, and the intake flow Q is 0.06 m³/s. The purpose of this experiment is to simultaneously measure the water-level and the vertical distribution of the longitudinal velocity when the open-channel flow steadily develops to a uniform flow (the normal depth is h₀) from a nonuniform flow. In the experiment, seven gauging points (#1-#7) are set up longitudinally (point #1 is 2 m away from the intake, and the other points are at 4 m intervals) and four velocity measuring cross-sections are at points #2, #3, #5, and #6. See Fig 4 below for the plan.

We can get a water-level lowering flow under the conditions mentioned above, and the experiment date will be applied to test the numerical model.

Validation of numerical model

The rectangular grids are applied to the mesh generation of calculation zones, as shown in Fig 5. Forty grids and 400 grids are distributed vertically and longitudinally, respectively. The aforementioned two-dimensional numerical model of Level Set-based open-channel turbulence is applied to simulate an open-channel flow under certain experimental conditions.

In this numerical computation, the total calculation time depends on the time required for the flow reaches steady state. We employ the following criteria to determine whether the flow has reached steady state: (1) There is an obvious uniform flow section existing downstream of the open-channel flow; (2) The variation value of the length of uniform flow section in adjacent interval times is less than 5×10⁻³. In each time step, the calculation converges when the average residual error (the arithmetic mean value of the sum of all-unit residual errors) is less than 10⁻³. The value of time step is Δt = 0.02s, which is associated with Courant-Friedrichs-Lewy (CFL) number less than 1.

The following sections will provide a description of the computational validation.

Water surface profile

The water surface is the water-air interface line at different locations in the longitudinal direction. Improving water surface capture precision is the purpose of the Level Set calculation in this paper. Fig 6 provides the water surface experiment and the verifying calculation results. They are well matched, as shown in the figure below.

Flow velocity

Fig 7 below provides the longitudinal flow experiment and verifying calculation results. The y-coordinate is u+=uu* and the x-coordinate is y+=u*yν in the figure. The following conclusions can be determined from the figure: (1) the longitudinal flow velocity experiment and verifying calculation results are well matched, and (2) the longitudinal velocity is still subject to semilogarithmic distribution within vertical range y⁺>30 from the wall during the evolution of open-channel flow from nonuniform flow to uniform flow.

Efficiency of computation

In this numerical work, the factors that affect the computational efficiency are as follows:

• The reinitialization of Level Set equation. The solution of Level Set equation is the key to this numerical computation, and it is also a challenging work because the Level Set function φ(x, y,) is implicit. And the reinitialization of Level Set equation is an important step in the process of solving. Because of the numerical effect, after the accumulation of several time steps, Level Set function no longer maintain the symbol distance characteristics in calculation. Therefore, the Level Set function must be reinitialized after each time step, which is the main reason for the heavy calculation workload. It can be argued that the reinitialization is directly related to computational efficiency.

• The time required for the flow to reach steady state. Although we are concerned about the characteristics of steady flow, in order to get the flow most close to the actual flow, we establish an unsteady numerical model to simulate the flow developing from the initial state to the final steady state. It is obviously that the time required for flow in different conditions to reach steady state is different, that is the total time of computation, so we can easily find that this time is another factor to affect the computational efficiency.

• The mesh size. For the same case, a more accurate free surface can be obtained by grid refinement, but the computational effort will increase greatly, thus the computational efficiency will be reduced.

From the above analysis, we can see that there are many factors that affect the computational efficiency, but the reinitialization computation is the most important factor.

Application

The same numerical model and mesh grids as the verification calculation are applied to simulate an open-channel flow in the aforesaid experimental water tank, in order to have a better understanding of general characteristics (including the water depth h, depth-averaged velocity U, and wall shear stress τ_b) and variation in the flow evolution from nonuniform flow to uniform flow and to acquire the variation features of general flow characteristic parameters (energy loss coefficient λ and roughness coefficient n) under nonuniform flow. The generalization of our research content is shown in Fig 8.

Variation of water depth and depth-averaged velocity

The water depth and depth-averaged velocity when the open-channel flow is finally in the uniform state are expressed as h₀ and U₀, respectively, and h₀ is used to nondimensionalize the x-coordinate. Fig 9A shows the h/h0—x/h₀ and U/U0—x/h₀ numerical simulation results. The following conclusions can be obtained from the figure: (1) under simulated working conditions, the water depth and depth-averaged velocity of nonuniform open-channel flow approximately change to the corresponding value of uniform flow at x/h₀ = 6~8 of longitudinal location and (2) under the same bottom slope conditions, the higher the Froude number is, the shorter the longitudinal length is required, which is the length required for nonuniform flow reach to uniform state.

Variation of wall shear stress

Fig 9B provides the wall shear stress τ_b process variation numerical simulation results. τ_b in the figure has been nondimensionalized by the final uniform flow wall shear stress τ_b0. The following can be concluded from the figure: (1) for an open-channel flow with water-level lowering, as shown in Fig 8, when the flow tends to be uniform, wall shear stress rises along the flowing path, and wall shear stress change amplitude τ_b/τ_b0 is the magnitude of 0.5~1 while water depth change amplitude h/h₀ is 1~1.6 and (2) under the same bottom slope conditions, the higher the Froude number is, the shorter the longitudinal flow length is required.

Variation of energy loss coefficient

The energy equation of open-channel steady flow is

In Eq (18), h_w is the energy loss between two sections and may be calculated using the Darcy-Weisbach formula hwL=λU28gh. Fig 9C shows the numerical calculation results of the energy loss coefficient λ process variation. The following can be concluded from the figure: (1) λ decreases along the flowing path amid amid the transition of flow from uniform to nonuniform; i.e., the nonuniform flow energy loss coefficient is higher than the uniform flow coefficient under corresponding conditions. From the quantitative point of view, the energy loss coefficient change amplitude λ/λ₀ varies by the magnitude of 1~1.4 while the water depth change amplitude h/h₀ is 1~1.6 and (2) under the same bottom slope conditions, the higher the Froude number is, the shorter the longitudinal flow length is required.

Variation of roughness coefficient

Roughness coefficient n is an important characteristic parameter for open-channel flow, and its relationship with the energy loss coefficient λ is expressed by the equation below:

Fig 9D provides the numerical simulation results of roughness coefficient n process variation. The following can be determined from the figure: (1) n decreases along the flowing path, i.e., the roughness coefficient of open-channel nonuniform flow is higher than that of uniform flow under corresponding conditions. Quantitatively, the roughness coefficient change amplitude n/n₀ varies by the magnitude of 1~1.9 while water depth change amplitude h/h₀ is 1~1.6. (2) under the same bottom slope conditions, the higher Froude number is, the shorter the longitudinal flow length required.

Conclusion

Level Set method is an excellent free surface capturing method, but the research on it has not been widely carried out. So we established the vertical two-dimensional numerical based on Level Set method. In order to test the validation of the numerical model, the open-channel flow experiment was conducted in a long straight tank by measuring the water-level and flow velocity synchronously, and the experimental results were employed in the verification computation. Finally, the numerical computation was carried out to study the flow characteristics of open-channel flow under different flow conditions.

Based on experimental research and numerical calculation, the following conclusions are made in this paper:

• The results of verification computation show that it is possible to precisely capture water surface profile of nonuniform flow by applying Level Set method, and the vertical distribution of the longitudinal velocity can agree with the experimental results well.

• A semilogarithmic zone of vertical distribution of longitudinal velocity exists amid the transition of flow from nonuniform state to uniform state;

• For an open-channel flow with a water-level lowering curve, the depth-averaged velocity and wall shear stress increase along the flowing path, but the energy loss coefficient and roughness coefficient decrease by a big margin.

Supporting information

S1 Appendix [docx]
Discretization.

S2 Appendix [docx]
Source code.

S3 Appendix [docx]
Reinitialization.

S1 Table [docx]
Notation.