Applying of depth-averaged hydrodynamic modeling with different methods of turbulence closure for flow near channel groyne. Part I – hydrodynamic background

Training structures in flow stream play an important role in shaping flow and bed properties. Planning to introduce such training elements like groins or dikes into the river stream one need to know consequences they may introduce into flow field and bed shear stresses. These consequences can be investigated by laboratory experiments on hydraulic models or by numerical modelling using hydrodynamic simulation models. In the paper the second possibility is exploited by applying two-dimensional depthaveraged model for straight rectangular channel with a groyne. This paper contains the first part of the research results and it describes hydrodynamic background of the flow phenomenon, concentrating on hydrodynamic equations for depth-averaged flow, types of eddy viscosity method used and kind of boundary conditions applied. Based on the hydrodynamic descriptions, different simulation experiments have been conducted for the flow problem and the whole analysis of simulation results for flow in channel near groyne is contained in the second part of the research activity (Part II = Analysis of simulation).


Introduction
Engineering structures for controlling river flow, such as spur dikes or groins, are effective methods for protecting river banks, stabilizing its navigation capacity and enhancing aquatic habitat.
Restorations of good flow conditions in natural channels often contain building of construction of spur dikes or groins, which extend from the bank projecting into the flow stream with different angles, usually in normal direction to the expected flow. The dikes should redirect flow into the main part of channel cross-section, protecting the riverbank from erosion and making stable pools for aquatic live with deposition of suspended sediment in recirculation areas behind groynes.
Deep view into the processes ongoing in space and time around such structures requires applying effective numerical models for simulation hydrodynamics and sediment transport, accompanied with some field measurements.
Flow in open channels around a training structure like groin or dike has a three-dimensional (3D) nature because of: • Irregular bed stresses -due to exchange of mass and transport phenomena in bed layer.

•
Curvilinear flow path and the secondary flowcirculation of water in transverse direction due to activity of the centrifugal force (near the water surface it moves toward the channel axis, and that near the bed moves toward the bank line). Consequently, the shear force, which has the same direction as the local flow close to the bed, deviates slightly from the direction of the mean flow ( [1], [2]).
For simulation such complicated flow phenomenon a full, good quality 3D hydrodynamic model would be the best choice, and, recently, several 3D models have been worked out in leading research centres. One of relatively advanced method of 3D simulation was presented by Kuhnle at al. in [3] for 3D flows near spur dikes. They confirmed results of simulations with measurements for physical test channel with the exception of significant discrepancy between simulation and measurement in the recirculation region downstream the groyne location. An interesting example of applying Large Eddy Simulation (LES) for turbulent flow around obstacles in the river flow path was work [4] of Teruzzi and Ballio, where they presented qualitative results for flow around the bridge abutment with complex highly 3D of vortex structures and shear stress fields tightly connected with scour shape around the abutment being formed.
Also, satisfactory results for the similar technique, large eddy turbulence model, were presented by Shahrokhi and Serveram in [5] for 3D simulation of flow around a groyne. Results of their simulation allowed for precise estimation of separation region behind the groyne in dependency on the groyne length and its inclination angle.
However, enormous efforts and costs are to be paid off to implement 3D models for real situations. Often, for practical engineering problems, such as maintaining navigational depths in rivers and channels, investigation of geomorphic processes in alluvial beds, estimating circulation zones behind spur dikes and similar, it is not effective and efficient applying the 3D models with very time consuming preparation of their numerical bathymetry, quality-degree computational meshes and substantial amount of input data. Much more effective is applying 2D versions of full Reynolds equations or other approximated methods of turbulence phenomenon treatment.
Several papers documented that it is profitable and sufficiently adequate to use 2D depthaveraged models to find solutions of practical problems with acceptable correctness. Results of numerical simulations based on CCHE2D model, reported in [6], showed good approximation of cross-sectional depthaveraged velocity with measured 3-D profiles made in laboratory hydraulic model. Another paper, [7], showed that this kind of CFD modelling can be successfully utilized for analysis of complicated flow fields in settling basins that in turn influenced sedimentation patterns at the bed of that basins.
But, the depth-averaged flow models should be applied with a deep cautiousness -they cannot properly transfer the momentum of flow because they are not able to include the secondary flow and the transfer should be corrected by some empirical functions. In several 2D hydrodynamic models some different solutions for the problem have been introduced, e.g. the prediction of the transverse component of velocity near the bed and surface, the use of dispersion terms in the momentum equations, the application of dispersion terms obtained from the streamwise and transverse velocity profiles in the curvilinear coordinate system.
One of pioneering researches in the applying numerical methods for solution flow field around spur dike was paper [7], where authors prepared depth-averaged hydrodynamic model with constant closure of turbulence and applied it to the flow near spur-dike. Also Duan and Nanda in [9] dealt with depth-averaged simulation of hydrodynamic flow field in a groyne proximity. For calculations they applied the mixed-length turbulence model for eddy viscosity (EV) both in the case of flow near spur dikes and in the case of suspended sediment concentration simulation in the area of groyne influence. Earlier, Duan applied the depth-averaged HD model with dispersion terms for momentum equations and parabolic EV closure for simulation of steady flow in natural meandering channels [10]. Their simulations compared to available measurements showed that the 2D solution should be used for real flow situations but in some cases (where vertical velocity profiles are significant) the need of 3D solutions became evident. An interesting example of depth-averaged flow-with-obstacle model is described in [11], where the concept of 'quadtree' grid for computational domain near obstacles (inlets, groins, etc.) was applied. In their description of hydrodynamic equations they introduced two different methods for the turbulence closure -the depth-averaged parabolic EVM and modified mixed-length model but they did not make any comparison between solutions based on that turbulence models. Instead of that, they showed (in the model testing part) results for steady flow around a spurdike (flow patterns, crosswise velocity profiles, velocity and water depth contours) that are similar to the results presented in this paper and can be used for comparison. One of the newer and more advanced applications of 2D depth-averaged hydrodynamic model is work done by Lee at al. [12] for the river flow with series of permeable pile groins where they applied k-e eddy viscosity model with Rodi formulation [13]. They successfully compared obtained numerical results with measurements in one selected case of flow but they did not utilize other descriptions of eddy viscosity field.
Reviewing sources connected with flow simulation fields near flow obstacles one can notice that small if none efforts were made in the scope of the influence of different closure method of eddy viscosity in depthaveraged Reynolds equations on the simulation results obtained. Only one paper [14], available to the author, linked the quality of numerical solutions with the turbulence model applied. Their authors considered five different EV closure models ( depth-averaged parabolic, modified mixing-length, standard k-e, non-equilibrium k-e and RNG k-e models) and investigated their inflows on solutions for three cases of flow: around a spur-dyke, in sudden-expanded flume and in natural Fall River meandering reach. Their conclusions differed substantially on which model is preferable for each of the three cases, and, especially in the case of spur-dyke, it seems that further research and comparisons would be profitable.
The work presented here is an attempt for preparation some hydraulic background for questions arisen while trying to achieve the following goals: • Investigating the influence of different turbulence closure models on simulation results obtained from 2D hydrodynamic model for mean steady flow conditions. • Choosing of such a turbulence model that would be possible to get flow velocities fulfilled in some determined limits, prescribed for given bathymetry and bed slope. That EV model should retrieve real circulation patterns in the regions determined geometrically and physically, such as areas near spur dikes, along stream bends, steady pools and sediment traps near the shorelines. • Estimation of 2D HD model performance about the capability of smoothly treating dry and wet nodes in the mesh of active flow area in order to produce velocity fields well fitted to actual morphometric and geometric conditions near flow training structures.

Hydrodynamic background
This study employed the solutions of the depth-averaged Reynolds equations in Cartesian coordinate, proposed in the CCHE2D model [15]. To include in some degree the effects of real-world secondary flow, the dispersion terms have been added for each momentum equations. They are obtained from assumed streamwise and transverse velocity profiles and they are transformed to Cartesian coordinates. Additionally, the hydrodynamic model can adopt easily (in the momentum equations) the variable density of flow allowing for coupling the hydrodynamic model with model of sediment transport model in the next step of research program (mobile bed with mass entrainment and deposition).
The governing equations are the two momentum equations (in [kg/(ms 2 )] = [Pa]): The bed shear stress term in (1) and (2) can be prescribed by two methods.
Knowing the stress components one can calculate the dynamic (shear) velocity: One can note that in practical problems the coefficient n contains in its value the effects of bed irregularities, bed forms, channel geometry, sediment particles size, vegetation, etc.
Second method comes from the depth-integrated logarithmic law for velocity distribution [16]: Here, z 0 is the smoothness parameter of the depth, depending on different flow conditions, k s is the height of bed roughness and v -the kinematic viscosity.
The relationship for shear velocity u  is implicit and its value must be calculated iteratively. Then, one can directly determine the Darcy-Weisbach's f c coefficient, needed to the stresses calculation. Formulas given by Van Rijn in [16] are used for these calculations:  Using values of f c , one can evaluate finally the bed stress components : In case of detailed simulation and model verification based on measured data, the second method seems to be more "physically profound" (provided the k s data are available). On the other hand, the k s parameter will also contain the effects of sediment movement and bed properties.
When the necessity exists one can attempt to convert values of n and k s using the Strickler's formula: There is no way to calculate values of dispersion terms from (5) -actual temporal velocities are unknown. In order to make some practical approximations for these terms, it is necessary to utilize empirical relationships between depth-averaged and local variables. For the streamwise velocity u l one can assume that it satisfies the logarithmic law along vertical coordinate: where: z -coordinate in vertical direction from bed to surface, u l (z) -component of velocity vector in stream direction, u  -shear velocity, z 0 -near-bed level calculated for bed stresses.
The profile of transverse component u r of the secondary flow velocity is taken as linear one accordingly to the following form: Here u r (z) is transverse velocity component along the water depth column, u r means its depth-averaged value and u s means its value at the water surface.
Based on the empirical relationships (5a) and (5b), an estimation of the dispersion terms in (5) is performed here using the procedure formulated by Duan and Nanda in [9] Finally, the dispersion terms (5) in Cartesian system of coordinates are expressed in the following way:   Reynold's stress terms in eqs. (1) and (2) where:  t -Eddy Viscosity (EV) coefficient [m 2 /s].
EV is not a constant -it is very complicated function of flow field variables, subjected to different theories of approximation. Several of them have been chosen in this study as a one of major factor in the estimation of the influence of the turbulence model on properties of the flow field obtained from simulation.

Conceptions of Eddy Viscosity modeling for 2D simulation
Three methods for depth-averaged eddy viscosity calculation were applied in this study: • Parabolic Eddy Viscosity (EV) model • Mixing-Length EV model • k- (k-e) EV model.
The parabolic model is the simplest one with the direct implementation into the calculation modules of HD model. It allows rather a general pattern of turbulence to include into the hydrodynamic solution of flow field. Several papers ( [4], [14], [15] induced that using the more advanced mixed-length model gives best solutions for more complex flow fields.
Smagorinsky model is an option of closure the Reynolds 2-D equations promising more adequate simulation of real flows, although this model of EV is more often used for 3D simulation techniques based on LES concept [5], and it is not used in this paper.
For applications where more detailed and precise treatment of the turbulence is required, the two-equation k-ε model for depth-integrated flows can be applied. But, it is the most expansive and prone to any nonlinear behaviour of basic flow parameters.
The goals of this study comprised incorporation all these models into hydrodynamic simulation, investigation of the influence of the model applied on results of simulation and deduction which of the models will be the most appropriate for flow field calculations in channels with bed and bank obstacles.
The simplest way of modeling eddy viscosity is given by parabolic model of EV: with   Near the actual river bank parabolic model can produce very high eddy viscosity values and it needs some improvement in that areas. The improvement consists of applying a distance to the wall instead of depth of flow as the length scale. Here, Fig. 1 shows how the normal distance from a node to the wall (d w ) is taken into calculation in the CCHE 2D model -the ratio of d w /h < 0.21 (= relative distance where the parabolic profile is equal to its depth-averaged value) should be fulfilled to use the d w value as the length scale. Other HD models utilize similar procedures with slightly different ratio value.
The mixing-length model calculates the eddy viscosity as a product of mixing length and gradients of 2D velocities in X-Y space ( [9], [11]): The mean mixing length L in eq. 8 is taken as an integral of its 3D distribution along depth of water column: According to Rodi's concept [13] for depth-averaged flows, the standard k-e model comprises the equation for k in the following form: and the equation for turbulence dissipation  : In above equations P km means the production of kinetic energy k from mean flow and it is calculated from the following formula: In eqs. 10a, 10b additional terms of turbulent energy production and dissipation of the energy are presentfrom bed friction, valid in case of uniform flow. That production of k (P kV ) and  (P eV ) is given by the following expressions:

Boundary and initial conditions for the hydrodynamic model
To start hydrodynamic simulation, it is required to prepare and input into the model all necessary boundary conditions -along inlet boundary velocity vectors are to be prescribed, along outlet boundary water level hydrograms should be described and along solid wall boundaries condition of slip or no slip velocity and initial water surface elevation in the computational domain are necessary.
At the inlet boundary one can prescribed the following conditions: Total discharge, discharge per unit width, flow velocity, hydrograph Q(t), hydrograph formula, water surface elevation.
At the outlet boundary the following conditions are possible: Constant water surface elevation, Free surface elevation, digitized rating curve, hydrograph Z w (t), hydrograph formula. For steady flow simulation the conservation of fluid mass should be achieved and controlled, allowing to determine the velocity field at the outlet boundary.
On solid boundaries: For water surface elevation the following is assumed: no pressure gradients normal to the wall are present.
For velocity components along that boundaries it is possible to apply one of the three possibilities: non-slip condition, partial slip condition or condition of full slippery. In case where more detailed solution is needed near the wall, one can use the logarithmic law of velocity distribution, provided that the quality and density of mesh near the wall is sufficient.
Additionally, in case of using k-e model for calculating eddy viscosity, it is necessary to provide another set of boundary conditions -for inlet, outlet and along solid boundaries, respectively. They are indispensable, because the k-e model solves its own transport equations.
In order to estimate energy of incoming turbulence, according to [18[, one can use the following depthaveraged expressions along the inlet boundary: Here flow is uniform without wall effect and  is the dimensionless depth of flow along the inlet line.
Afterwards, the eddy viscosity along the inlet is calculated by using the following expression: h u t * 6 1    Then, the value of dissipation ε along the inlet boundary is taken from transformation of eq. (10c): Along the outlet boundary it is assumed that the turbulence energy is streaming out freely, accordingly to the following conditions

 
where: u  -shear velocity, y -the distance between the current wall line and the closest mesh node in the actual calculation domain.
Besides the boundary conditions defined, another set of input data about all flow variables (u x , u y , Z w ) is required at the beginning of simulation period (initial time). This initial condition can either be prescribed externally (Cold Start) or adopted from previous solutions of the same case (as Hot Start). Using Hot Start method substantially advances the process of reaching final stable, correct and balanced solution of the given problem.
Computationally, each hydrodynamic model must iteratively solve a huge non-linear set of algebraic equations, that arises from time and space discretization of the governing equations, boundary conditions and other relationships at ever time step of simulation period. Thus, checking of the convergence process at each step of simulation is required and it can be done successfully by observing and controlling the variations of increments of the computed flow variables u x , u y , Z w and h.
The following relative maximum increments can be investigated across the whole modelling domain for every iteration and every time step: Convergence of solution at each time step will be satisfied if the maximum relative increments of velocity and water level are less than the accuracy limits assumed for the given run.
For real situations, the physical domain of modelling can shrink or enlarge, depending on boundary conditions applied and the convergence process being worked out and lots of dry nodes may appear or disappear. In such cases values Max u tot and Max Z tot in (12) could not reach the presumed accuracy limits but, usually, the convergence of the main flow will be satisfied.
Another way for checking the quality of solution obtained from such kind of numerical calculations is the condition of preservation the mass balance along some selected cross-sections of computational domain. This is especially helpful for steady flow cases where the discharge balance preservation is theoretically predominant condition. This rule is utilized in the next part of the research for simulation steady flow in a channel with groyne using the presented above theoretical descriptions of the numerical model, based on the CCHE 2D FEM solution. This model will be adopted to the following channel finite element mesh (Fig.2): Steady flow in the channel given in Fig.2 will be calculated using the numerical counterparts of the equations presented here and an analysis of results obtained will be conducted in the Part II of the research.

Conclusions
Depth-averaged two-dimensional horizontally set of hydrodynamic equations have been characterized as a calculation model for channel flow with a groyne in its bathymetry. Three turbulence models for depth-averaged 2-D Reynolds equations -parabolic, mixed-length and k-e models -have been considered in this study for simulation of steady flow in the test channel with groyne. Short discussion of kinds of advisable and necessary boundary and initial conditions has been given in the paper. They will be practically utilized during preparation of input data for hydrodynamic model of channel with groyne in the second part of the research. The basic checking criteria, relative maximum increments (12) and the balance of mass (discharge) have been suggested for checking the performance of each turbulence model. The findings from this simple analysis may have practical value regarding the quality of simulation results for any hydrodynamic model -checking the model quality should be often performed with suitable and simple criterion. Finally, one can underline that the quality of simulation of flow parameters is essential for other research activities e.g. in the field of fluid-sediment dynamics, due to importance of bed shear stresses and any effort for better understanding the obtained HD solutions would probably pay off with better sediment movement simulation.