A Computational Study of Interaction of Main Channel and Floodplain: Open Channel Flows

The assumption of steady uniform flow permits the computation of the velocity isoline, secondary current and turbulent statistics in open channel flows. However, it becomes important to choose appropriate turbulence models to capture the length scale of turbulence near the interfacial zone of compound channels. This paper not only focusses on capturing the longitudinal vortex and primary mean velocity but also extrapolates the results of numerical analysis to understand the interaction between the main channel and floodplain in asymmetric compound channels. The results of computational fluid dynamics simulation showed that the velocity isoline bulging near the bed of the floodplain and sidewall at the junction, due to high-momentum transport by secondary current, can be captured with Reynolds stress model. Furthermore, by applying the three different cases of channels with varying geometrical aspects, the maximum velocity simulated showed similar results to the experiments where the structure of primary mean velocity is seen to be influenced by momentum transport due to the secondary current.


Introduction
The turbulent structures in compound open channels are characterized by large shear layers generated by the difference of velocity between main channel and floodplain flows.Differential flow velocities in these two different sections create vortices with vertical and longitudinal axis.Due to the in homogeneity of turbulence created in the shear layer region, these induced vortices are latterly defined as plan form vortices and secondary currents, respectively.The secondary cur-rents are found to influence the primary mean velocity field, which was experimentally explored by Tominaga and Nezu [1].Many other researchers lately have investigated on this topic and clarified the mean velocity and boundary shear stress in compound channel flows [2]- [7].In retrospect, the idea of mapping the channel cross section and defining these sub-regions by bounding orthogonal to the isovels is attributed to Leighly [8].
In experiments, it is quite difficult to obtain accurate secondary velocity components in open channel flows, because their order of magnitude is only 1% -3% of the primary mean velocity [1].However, the understanding of experimental techniques has evolved tremendously in these time using non-intrusive techniques of measurements.Nevertheless, it comes with costs and the technical superiority of intertwined interdisciplinary knowledge.In these cases computational fluid dynamics plays a pivotal role.Previously, researchers have detailed studies on turbulence closure models, which have eased the process of modelling [9] [10] [11] [12].Furthermore, recent application of these models on the open channel flow has intimidated researchers with producing better results for a wide range of channels including complex geometries like meandering, gravel bed and complex flows in vegetation [13] [14] [15] [16] [17].
In this study, the main focus was to simulate three conditions of asymmetric compound channels taken from Tominaga and Nezu [1].These channels are smooth homogeneous and have the bank full height of 20 mm (S1), 40 mm (S2) and 60 mm (S3).The first model which was tested to verify the primary velocity isolines and secondary current using the k-ϵ model, which was repeatedly simulated is S1.However, the results of the two aforementioned characteristic hydrodynamic variables were found to be sublimed.Basically, the modelling results showed that the model was unable to capture the bulging near the junction of two sections, and that the secondary currents were not visible in the contour plots.Nevertheless, the maximum mainstream velocity U max obtained in the simulation for S1 was 0.411 m/s, which was experimentally reported as 0.409 m/s.The percentage is overestimated around 1.2% which seemingly outcasts the above argument.Thus, a question raised here is whether or not the results obtained in these models are permissible without identifying the critical mechanism of flow?
To answer this question, another model Reynolds stress turbulence (RSS) was used and found to give a better outlook of quasi-2D structures in the cross section for better understanding.It has already been justified and proposed that these seven-equation models are more robust on the modelling of open channel flows through commercial generic CFD packages [18].Thus the key idea of this study is to capture secondary shear effects closely and then reveals these simulated results to understand the interaction between the main channel and floodplain.
The main aim of this paper is to study three dimensional structures of mean velocity and turbulence via CFD, in which the RSS model was used to produce the distribution of Reynolds shear stress, turbulent kinetic energy, and turbulence intensity.This analysis will overlay the impact of 3D structures in compound channel infused with insight of computational ability for simulation in ANSYS [18].The simulations showed good results for the mean velocity isoline and secondary currents for S1, S2 and S3.In addition, results of the turbulent statistics were discussed keeping the aforementioned corroboration as zero ground.

Turbulence Closure Models
The governing Reynolds averaged Navier-Stokes equation (RANS) is derived from fundamental fluid flow equation, which is theoretically and physically based equation of governing fluid flow.Through the control volume and force balance for an incline bed of open-channel in the stream wise direction, RANS equation can be combined with the continuity equation to give: Thus, Equation ( 1) is a force balance equation where the secondary flows are considered to equal the weight and gradient of (lateral + vertical) Reynolds stresses forces.Here {U, V, W} are the velocity components in the {x-longitudinal, y-lateral, z-vertical} directions respectively with respect to bed.Furthermore, ρ, g, o S are the fluid density, gravitational acceleration and channel bed slope respectively.Second order tensor { , ij yx zx τ τ τ = } are the Reynolds stresses in the direction x with its plane perpendicular to y and z, respectively.

K-ε Model
The K-ε is the most famous two-equation turbulence closure model, which is found in many CFD commercial codes, e.g. in ANSYS-Fluent.High-frequency fluctuation can be seen in the steady state analysis, commonly.To encapsulate this issue, time-averaging methodology is used, which gives additional terms that need empiricism.For the closure of solution, these additional terms are to be expressed as empirical approximations.The transport equations for the turbulent-kinetic energy k and its dissipation rate ε based semi-empirical model is expressed by the following equations: ( ) where k is the turbulent kinetic energy, defined as is the fluctuating quantity of velocity; p  is the instantaneous pressure; u′ is treated as the time and average velocity component; u' u' is the Reynolds stress and t ν is the eddy viscosity.The eddy viscosity employed in the calculation is quantified as Equation ( 4): The exact k-equation does not help for the closure of the model since other unknown correlations appear in the turbulent transport and dissipation terms, which lead to more unknown beside 10 unknowns of Navier-Stokes equation itself (three velocity components, one pressure term and six Reynolds stresses).To counter this problem empiricism and assumptions are used for the parameters.

Reynolds Shear Stress Turbulence Model (RSS)
To account for turbulence, the momentum and continuity equations are time averaged and the velocity as a whole is apportioned as fluctuation and sum of a mean component, which is established under Reynolds-averaging process.This leads to an extra term in the momentum equation as a consequence of the non-linearity of the advection term.This extra term is identified as the Reynolds stress which needs a mathematical closure together with additional transport equations.Thus, the RSS turbulence closure model consists of six Reynolds stresses (due to symmetry and one equation for the dissipation of turbulent energy making it a seven-equation model).
In the present analysis, Speziale, Sarkar and Gatski [19] (SSG) Reynolds shear stress turbulence model was used.The transport equation for the model is: where ij δ is the kronecker delta ( the right hand side of the equation has stress production, diffusion, pressure-strain, dissipation terms with µ and t µ being dynamic and eddy viscos- ity, respectively.In this specific model, the pressure strain term is modeled as: where ij a is the Reynolds-stress anisotropy term, ij S the strain rate, and ij Ω is vorticity, which are defined as: and (1 / 2) kk P P = , "k" is the turbulent kinetic energy and ε is the dissipation rate computed from the additional transport Equation ( 8).

(
) ( ) The constants in the SSG model are listed in Table 1.The value of C µ is taken as 0.09 for estimation of the turbulent eddy viscosity.The model used in the current prismatic compound channels is adequately sufficient to capture the secondary shear stress, since the vorticity equation dominates the anisotropy in the normal Reynolds stresses.

Model Setup and Solver Setting
Pre-processing of these simulations is started through creating the geometry of the three cases aforementioned as S1, S2 and S3 (Figure 1).ANSYS-Fluent [18] relies on the finite element based on the finite volume method with the solver Algebraic Multi-Grid (AMG) acceleration techniques.In between pre-processing and solving, there comes meshing part, which includes definition of cell sizes.In this particular instance three meshes with different element sizes were generated in the ICEM CFD [18].All these mesh configurations were treated with the same solving criteria such as second order discretization technique, a residual target of 10 −5 -10 −7 was used to capture secondary current terms, thus leading many thousands of iteration to obtain full convergence.All these parameters were chosen to facilitate the convergence indicators in the simulation undertaken here, while the intensity of secondary current terms are typically two order of magnitude smaller than downstream mean hydrodynamic quantities.

Boundary Conditions
Four boundary conditions were used here: 1) velocity inlet, 2) pressure outlet, 3) symmetry for surface flow, 4) no slip conditions for wall and bed.

Inlet (Velocity) and Outlet (Pressure)
At the inlet, turbulence properties i.e. k (turbulence kinetic energy) and (ε turbulence dissipation rate) must be specified.These were calculated as [20]: where I is the turbulence intensity, U is the mean value of stream wise velocity and 1000 t I µ µ = .At the outlet, the pressure condition was given as the boundary condition, and pressure was fixed at zero.

Wall, Bed and Surface
A no-slip boundary condition was applied at the walls.This means that the velocity components should be zero on the walls.Furthermore, the standard wall-function which uses log-law of the wall to compute the wall shear stress was used.
The water surface was defined as a plane of symmetry, which means that the normal velocity and normal gradients of all variables are zero at this plane.

Results
Table 2 shows the case definition and comparison of the simulated results with the experimental results.The simulated results showed a very close and consistent agreement with the experiment.Note that the results given in [1] are all normalized with U max and thus simulated results have been represented in the same order and scales.

Secondary Currents (Longitudinal Vortex)
Figure 2 & Figure 3 reveal the vector description of the secondary currents (U, V) simulated through two models discussed above.Figure 2(a) & Figure 2(b) interpret the secondary vector results of K-ε and RSS models.The incapability of K-ε model requires a higher order turbulence model to capture re-circulation cells, which identify the direction of secondary current in vector plots.
In Figure 2(b), inclined up flow at the junction was visible and more prominent in the floodplain region, similar to the experimental data.Moreover, a weak main channel vortex in the S1 case where H/h = 0.75 was visible as same that of   the experimental plot (Figure 2(c)).The floodplain vortex was so strong that it reached the free surface as also depicted in the experimental results.To understand the similarity index of both simulated and experimental results, case S2 should be a good example (Figure 3).As (H/h) decreased, pronounced effect of secondary cell at main channel and floodplain increased significantly.These longitudinal vortexes not only covered the junction region due to differential velocity but also occupied the far end of the floodplain channels, which showed anisotropy of turbulence due to wall effect, also evidently visible in straight rec-tangular channels [21].These secondary cells were easily captured in the simulated results, although they were missing in the experimental results due to the inability of the data from the weak spots.The two vortices that reached the free surface and covered the junction edge in Figure 3(a) are called main-channel and floodplain vortex, respectively.Experimentally, the region of vortex given in [1] was in between 1.6 / 3.0 y H ≤ ≤ which is the same region speculated in our simulation.These results showed the potential of the RSS over the K-ε model and also depicted that the results promise for further investigation on turbulence characteristics in depth for understanding the interaction of main channel and floodplain flow on simulation ground.

Distribution of Primary Mean Velocity
For clarity and affirmation on the simulated results for the smooth asymmetric channels, the contour mappings of primary mean velocity are illustrated in show that the isovel line bulges upward in the vicinity of junction edge as a result of the deceleration of velocity, thus pushing towards the bed of the floodplain and sidewall at the junction.This phenomenon was appropriately mimicked in the simulated results of the RSS model (Figure 4(b)).This is not surprising that the K-ε model in Figure 4(a) was not able to produce the bulging because the high momentum transport at the junction is secondary current induced, as not shown before.
In Figure 4(b), in spite of large flow depth at the floodplain, the primary mean velocity isolines are shown with the bulge upward at the junction, while the deceleration due to the low momentum transport from junction edge seems   to extend to the free surface like that in experimental cases (Figure 4(c)).However, these results show more conspicuous in the cases of S2 and S3 in Figure 4(a) & Figure 5(a), which is predictable due to the increasing momentum transfer rate for the low depth ratio.Furthermore, the deceleration due to the bottom vortex in the cases of S2 seems more clearly visible and is accurately captured in the simulated results.
For case S3 in Figure 5(a), interestingly as pointed in experimental results, velocity isolines no longer significantly bulge from the junction towards the free surface.Rather it bulges towards the sidewall of the main channel because of the very low depth ratio in rectangular open channel cases.Furthermore, the velocity dip phenomenon shows a remarkable distinction in case S3 as the depth ratio decreases.The isoline bulging in the plot of k is inclined outward the same direction as that of the experimental result given in Figure 7(c).This implies that the turbulence rises in the proximity of the junction of the two sections of flows.However, these isolines never overpower the primary velocity isolines as depicted by the experimental results.Furthermore, Figure 7(b) reveals the results of the turbulence intensity (I), which is normalized by the root-mean-square of the velocity fluctuations.Due to secondary currents, this plays a key role in understanding the connection between the primary mean energy and turbulent kinetic energy due to secondary currents.The isolines in Figure 7(a) have the same trend as those of the k plot in Figure 7(c).Such redistribution of the turbulent energy and turbulence intensity shows that the higher anisotropy occurs near the wall, bed and junction of the two sections, which was primarily revealed in the experimental data as well.The overall result in this category has 4.2% of overestimation near the boundary.However, in the experimental results the isolines close to the boundary are missing, which is considered to be weak spot for experimental data collections.These two terms significantly indicate that there exits the apparent shear stress, which is recently explored in depth in [22] [23] [24] [25] [26] studies.However, our simulated results make the picture more clear and profound.

Turbulence Statistics
The characteristic behavior of the spanwise Reynold shear stress has a negative value which increases at the junction (y/H = 2.6) and the free surface, showing an inclined upflow, also discussed in the experimental results of [1].Furthermore, the vertical distribution of the Reynolds shear stress in Figure 8(b) has a concavity at the surface between the main channel and floodplain side walls.At the center of the junction (z/H= 0.2), the Reynolds stress has a minimum and then increases towards the free surface.In the experimental results, this behavior of the Reynolds stress implied the momentum transport from the

Conclusion
The numerical modelling in this study for the smooth asymmetric open channel flows showed consistent results with the experiments.We have discussed the hydrodynamic parameters like secondary current and primary velocity isolines, which showed near exceptional results in capturing major features of flow using the Reynolds shear stress model in CFD.Further analysis of turbulence statistics was given about the interaction of main channel and floodplain flow.The results obtained in this study have shown the capability of the RSS model in capturing the momentum transfer phenomenon at the junction of the edge of the two sections.However, overestimation in the turbulence characteristics has been identified over the bed and wall region.The overall results showed that the internal shear and secondary currents are not negligible and hence should be included in any kind of hydrodynamic studies in compound channels.The distribution of Reynolds shear stress notably indicated the increase in anisotropy at the junction edge, as pointed in the experimental analysis and replicated as well in our simulation.In addition, the values near the boundaries were noteworthy in the simulated results, which were not visible in the experimental analysis, since these points close to borders are identified and generally overlooked as weak spots.

Figure 1 .
Figure 1.Schematic cross sectional geometry, H is the total depth and h is the flow depth of flood plan from S1 to S3 but in the figure it only denotes the depth for S1 configuration.

Figure 3 .
Figure 3.Comparison of secondary current vectors for S2 for depicting both main channel and floodplain longitudinal vortex: (a) RSS results and (b) experimental results from Tominga and Nezu [1].

Figure 5 .
Figure 5.Comparison of velocity isoline results for S2 depicting pronounced bulging at the junction of main channel and floodplain: (a) RSS results and (b) experimental results from Tominga and Nezu [1].

Figure 7 (
Figure 7(a) & Figure 7(b) show the normalized turbulent kinetic energy (k) and turbulence intensity (I) by the experimental averaged friction velocity U * of 0.0164 m/s for S2 case.For comparison, the turbulence characteristic plot of k is depicted for case S2 only ecause it was given in [1].

Figure 8 (
Figure 8(a) & Figure 8(b) represent the shear stress due to spanwise and vertical Reynolds stress.The peak values of these terms are high at the junction of the main channel and floodplain, specifically for the spanwise Reynolds stress.

Figure 8 .
Figure 8.(a) Spanwise and (b) vertical Reynolds shear stress terms in case S1, illustrated to be strong at the junction between the sections and near walls.

P
. K. Singh et al.DOI: 10.4236/jamp.2020.8111882537 Journal of Applied Mathematics and Physics main channel to floodplain.Thus Figure 8 demonstrates that the momentum transport over the two sections was well captured in this study.

Table 2 .
Comparison between experimental and simulated results.