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

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.


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 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 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 wallmounted 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 openchannel 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.

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 lowdimensional curve plane topological structure.
As shown in Fig 1 above, the open-channel profile of the two-dimensional flow is zoned airflow V 1 , water flow V 2 and interface S(x,t) between V 1 and V 2 . 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: where u and v are the RANS flow velocity in the x and y directions, respectively. Eq (2) is called the Level Set equation, in which the symbolic distant function φ(x, y, t) is the Level Set function.

Motion equation
The calculation zones O consist of airflow zone V 1 and water flow zone V 2 . 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: where u and v are flow velocity in the x and y directions, respectively; u 0 ,v 0 are fluctuating velocity components; À ru 0 v 0 is Reynolds stress; p is the pressure intensity; G k is the produced turbulent kinetic energy; θ is the dip angle of open-channel; ρ and μ are the density and dynamic viscosity coefficient of air and water, respectively; ν t is the turbulent viscosity coefficient; C μ , σ k , σ C ε1 and C ε2 ,are all empirical coefficients, which are valued for the constant in the standard k-ε model [31]: C μ = 0.09, σ k = 1.00, σ ε = 1.30, C ε1 = 1.44, C ε2 = 1.92(S1 Table).

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: where F is the generalized variable, Γ is the generalized diffusive coefficient, and S 0 (F) is the source term. These variables have different meanings in different equations. 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: where L(φ) is the spatial operator, φ is the Level Set function. 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 0 , i.e., the left derivative φ-x and right derivative φ+ x. To solve the left derivative, as shown in Fig 2, To obtain the right derivative, choose another six node templates and similarly obtain the five first-order Newton difference quotients, N 1 , N 2 , N 3 , N 4 , N 5 .
Following the WENO method, the convective term discretization pattern through weighting is obtained: where parameter ω is the weighting coefficient. It can be obtained it by the reference [36]: where a is the transition coefficient, they are given as follows: where ε 0 is an extremely small number which is used to avoid the denominator of a becoming zero. In our later computation, the value of ε 0 is 10 −7

.
With the upwind scheme concept, φ x is equal to φ-x when u i0 >0 or φ+ x when u i0 <0 for node i 0 . 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 0 , 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 k lnðEy þ Þ and y þ ¼

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. i. Initial assignment: set the boundary condition and assign starters for the entire field. iii. Free surface equation solution: a. Obtain φ t with reference to the obtained flow field (u, v) and last moment φ t-1 .
b. Reinitialize φ t , keep the symbolic characteristics, and obtain output result.
c. 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: where ε 1 is the low amount rectifying parameter, usually ε 1 = 1.5O(Δx,Δy), we take ε 1 = 1.5Δy in this calculation. ρ and μ in the whole field can be expressed as rðx; yÞ ¼ r a þ ðr w À r a ÞHðφÞ mðx; yÞ ¼ m a þ ðm w À m a ÞHðφÞ ð16Þ ( where the subscripts w and a represent water and air, respectively, H(φ) is smoothed Heaviside function.
(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 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 0 is taken to be 0.002, and the intake flow Q is 0.06 m 3 /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 0 ) 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 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 −3 . 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 −3 . 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. Efficiency of computation. In this numerical work, the factors that affect the computational efficiency are as follows: 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.
2. 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.
3. 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

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 0 and U 0 , respectively, and h 0 is used to nondimensionalize the xcoordinate. Fig 9A shows (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 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 . 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 λ/λ 0 varies by the magnitude of 1~1.4 while the water depth change amplitude h/h 0 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 0 varies by the magnitude of 1~1.9 while water depth change amplitude h/h 0 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: 1. 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.

2.
A semilogarithmic zone of vertical distribution of longitudinal velocity exists amid the transition of flow from nonuniform state to uniform state; 3. 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.