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

Training structures in flow stream play an important role in shaping flow and bed properties. To investigate consequences of introducing training elements like groins or dikes into the river stream one can successfully use numerical modelling. This paper is the second part of the research concentrated on utilizing some hydrodynamic models for resolving the problem. In the first part some hydrodynamic background has been described and here the analysis of applying two-dimensional depth-averaged model for straight rectangular channel with a groyne is discussed. Three models of eddy viscosity were applied to investigate their influence on results of simulation and to attempt choosing the most suitable method of modelling the flow phenomena around the groin. Results of simulation in steady flow conditions show that from the hydrodynamic point of view the more profitable method of turbulence description is the mixed-length or even the parabolic method instead of the broadly suggested k-e model. They allow producing the most reliable vortexes in the shear layer of flow behind the training structure end well preserve the mass balance contrary to the results obtained with k-e model of eddy viscosity.


Introduction
This is the second part of the research effort concentrated on the numerical modelling of flow near the groyne. In the first part some theoretical background of hydro-dynamic description of the flow phenomenon has been discussed. Here the main goal is to discuss results of solution the problem using 2D versions of Reynolds equations with three methods of turbulence closure. Several papers documented that it is profitable and sufficiently adequate to use 2D depth-averaged models to find solutions of practical problems with acceptable correctness. Results of numerical simulations based on CCHE2D model, reported in [1], showed good approximation of cross-sectional depth-averaged velocity with measured 3-D profiles made in laboratory hydraulic model. Another paper [2] has 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. In the experimental field, one worth to mention the pioneering Rajnataram's work [3] for steady flow in rectangular laboratory channel with a groyne. Paper [3] contains detailed measurement and analysis results for flow parameters, shear stress fields and other quantities as functions of groyne length. Results from that paper were used by other researches for verification their own models. In the present study, the simulation results for flow in a groyne neighbourhood also are referenced to that obtained by Rajnataram's experiments.
More recent experimental results were found in the paper [4], where profiles and contours of velocity, depth and other flow quantities are presented for single groyne and two symmetrical groynes. Also, important results of applying numerical methods for solution flow field around spur dike are presented in [5]. The compared inviscid and viscous solutions for the groyne in channel flow were the numerical analogous of the Rajaratnam's laboratory flume. Another work [6] also referenced Rajaratnam's results with depth-averaged simulation of hydrodynamic flow field in a groyne proximity. In [6], the mixed-length turbulence model for eddy viscosity (EV) is applied. In [7], simulations of steady flow in natural meandering channels compared to available measurements showed that the 2D solution can be used for real flow situations but in some cases (where vertical velocity profiles are significant) the need of 3D solutions became evident. In [8], two different methods for the turbulence closure are introduced -the depth-averaged parabolic and modified mixed-length model but without any comparison between solutions based on that turbulence models. Another reviewed paper [9] contains more advanced applications of 2D depth-averaged hydrodynamic model for river flow with series of permeable pile groins using k-e eddy viscosity model. Simulation results are there successfully compared with measurements for one selected case of flow but no other EV models were used. 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. Thus, the work presented here is an attempt to achieve some valuable answers for that question.

Simulation of flow in channel with groyne
Lack of available measurements data was one of the reasons for choosing a model of real, small, laboratorysize channel in simulation experiments of this study, for which available results were published in the scientific literature. In this context a choice of numerical equivalent of the laboratory channel in University of Alberta (Canada), used by Rajaratnam and Nwachukwu for experiments in [3], seemed to be the most proper. For that channel with a spur dike, they carried out measurements given in [3] and their results are here selected to verify the simulation of hydrodynamic flow field in the numerical model of the channel. They measured flow parameters of a series of experiments in a straight rectangular channel with bed slope S 0 =0.001, which was 120 ft.

Reference solutioncase of a straight rectangular channel without a groyne
The same artificial channel has been "constructed" as a finite element mesh for hydrodynamic simulation -in version without a groyne (reference polygon) and with the groyne of the size and shape similar to the laboratory case. In Fig. 1 some details of the numerical channel with groyne are shown. The calculation domain consists of 70 longitudinal profiles and near 1000 crosswise profiles that create near 70000 finite elements, particularly densely distributed in the proximity of the groyne location and near the channel walls. For testing purposes, first experiments were carried out for the test channel in absence of a groyne -for steady flow conditions. All three turbulence models discussed in the paper have been applied, giving similar solutions of velocity field distributions across the channel (Fig.2) only for parabolic and mixed-length models. Model k-e models produced different solutions (Fig. 3) for each water flow parameters. The obtained solution (Fig.3) were stable in time, and, using cross-section integration they gave values of total discharge and mean channel flow velocity close to the assumed boundary condition Q 0 =0.045 m 3 /s and U 0 =0.4 m/s (equal to flow parameters used by Rajaratnam -at least for solutions obtained for parabolic and mixedlength EV methods). In the case of k-e method there were significant differences from Q 0 and U 0 , especially near the walls of the channel. Solutions of flow parameters for mixed-length, parabolic and k-e methods are taken further as a reference background for analysis of simulation results for flow in channel with a groyne installed.

Calculation of steady flow in the channel with a groyne
In the middle of the channel a groyne has been built, sticking perpendicularly from the right bank into the flow stream -causing the FEM mesh to be reconstructed -the RL (Ryskin and Leal) orthogonal mesh ( Fig.4) with smoothness control has been applied. The groyne resembles the spur dike used by Rajaratnam (although not perfectly, because of the limitation of the FEM grid technology applied. Due to specific requirements of numerical methods applied to governing flow equations with dispersion terms used accordingly to the method described in the hydrodynamic background [10], the preferable way to achieve final state of steady flow was the unsteady flow regime calculation with steady boundary conditionssteady flow at upper entry to channel and steady or normal-depth condition at the downstream end of the channel. Fig. 5 shows progress with achieving the steady flow solution by model inside the channel that starts its runs from the "stream at the rest" condition with zero flow fields at initial moment of simulation. Analysing time processes of several model's runs, it appeared that in most cases around 54 seconds is consumed by the process of spinning up the model to obtain the solution that can be considered as the final state of the given steady flow. It is roughly twice the translation time (26.5 s) for shallow wave travelling along the test channel (c=(gh)=1.4 m/s). After that time the unsteady solution becomes the steady one and one can analyse its parameters, e.g. checking the flow balance (mass balance control), or, comparing the obtained solution (velocity field, shear stresses, etc). Many runs of the model have been conducted and they confirmed that the speed of convergence of model simulation to the given state of uniform flow, controlled by eqs. (12) in [10], depends on several factors, e.g. time step of integration schemas, geometric properties of channel (e.g. with/without training structures present), roughness of the channel bed. All next model runs utilized that way in order to effectively eliminate the inflow of non-steady disturbances on obtained results.

Verification of flow simulation results in the groyne area
Confrontation of the simulation results with Rajaratnam's experiment data is presented here in order to assess the model correctness against some experimental results. Rajaratnam carried out steady flow experiments in their test channel with a groyne installed in the middle of channel by the right bank (with the groyne of 0.15 m sticking perpendicularly into channel stream. In [3] are presented data about longitudinal profiles of shear stress (SS) in the vicinity of the groyne. It is noted in [3] the characteristic presence of SS peaks just behind the groyne positions. Similar shape of longitudinal SS profiles have been obtained in this simulation ( Fig.6 with x-axis given in ft. (x=0 at the groyne location) and for bed shear stress recalculated to psf units -for comparison the results with Rajaratnam's data). One can notice on the graphs in Fig.6 that, as the tip of the groyne is approached, the total bed shear stress (TSS) increases rather rapidly from the upstream level of  0 to the maximum of about to 0.012 -0.013 psf. Although shapes of the longitudinal profiles differ slightly from those measured by Rajaratnam, the peak values agreed well with the measurements presented in [3] -registered values of  0 are equal to 0.013 psf (0.62 Pa) on a profile y=0.5 ft. In Fig. 6, y is horizontal axis perpendicular to the x longitudinal channel axis and it corresponds to the z-axis in the Rajaratnam's work. From Fig. 6 it appears that immediately behind the groyne the bed shear stress diminishes rapidly in the area close to the right channel bank (y=0.15 m and y=0.23 m). This is caused by backward-flow region that is formed downstream the groyne by the right bank. This profile change also agrees with Rajaratnam's experiments. As is shown in Fig. 7, Rajaratnam's flume slope S 0 for uniform flow conditions has also been maintained in this model as well as the Froude number profiles (Rajaratnam's F r values were in the range of 0.2 -0.25 in channel far from the groyne, here the F r values varies between 0.2 and 0.27. In flow area affected by the groyne a characteristic change of water level is noticeable what also is found in the Rajarantam's paper. Profiles of F r show that it rises to its local maximum in the middle and left side of the channel and it decreases rapidly just behind the groyne -in the wake of the groyne near right bank (backward flow circulation). It agrees with the results of experiments, described in [3] and in [6]. Rajaratnam measured flow velocity along crosswise lines in different locations downstream the groyne. The measurements are used in this research for checking the correctness of simulated velocity fields (Fig.8). On Fig.8 they are shown (as black triangles) against the solutions from three different calculation cases: Case 1mixed-length method for EV, Case 2 -parabolic method and Case 3 -k-e equations. The velocity profiles are taken in 5 relative (x rg ) locations of cross-sections (x rg =0 = groyne, x rg =2 means 2ˑL g location downstream, L ggroyne length). Fig.8 shows that the calculated velocity profiles agreed reasonably with the experimental results, at least for mixed-length and parabolic methods. In case of k-e method there are more substantial discrepancies.

Analysis of flow simulation near groyne with different Eddy Viscosity Models
Three series of simulations has been conducted in the frame of this research and they were grouped in three different cases: Case 1 comprises all simulations for mixed-length model of eddy viscosity, Case 2 gathers all solution with parabolic EV model, and, Case 3 contains all results obtained with k-e model of EV. In every case the same steady flow conditions analogous to that from Rajaratnam's experiment were used: discharge Q=0.045 m 3 /s at the upper end of channel, undisturbed depth of flow = 0.19 m, normal flow condition at lower end of channel. Manning's roughness has been estimated (by Strickler formula (4c), [10]) as n=0.018 to resemble friction parameters of Rajaratnam's smooth bed (roughness height of 0.6 mm).

Case 1. Flow simulation using Mixed Length Method for Eddy Viscosity.
In Case 1 all runs of hydrodynamic model have been performed with the mixing-length model of EV (eq.(8), [10]). For each solution, the following parameters were analysed: velocity magnitude (U) together with velocity vectors, eddy viscosity (EV), water depth h and F r , in the channel near the groyne location. Only U and EV fields are shown in Fig.9. As shown in Fig.9, presence of the groyne affected the flow in the whole cross-section, especially downstream of the groyne, while in the upstream direction smaller flow area is influenced. Downstream of the groyne, a zone of re-circulating flow currents is reproduced. As is expected, the groyne forced the flow toward the central part of the channel, increasing its velocity due to flow field contraction. This change of velocity caused local depression of water surface near the groyne and influences distribution of EV i F r fields. In Fig. 10 dimensionless graphs of water surface and velocity magnitude are shown (with reference data of Z w0 and U 0 taken from hydrodynamic solution without the groyne (Fig.2). Fig.10(top) shows that water surface increased by 6 % at the groyne tip and decreased by about 3 % in the backward zone behind the groyne. In the short distance from the groyne nose flow velocity experienced high jump up to 150 % of its referenced value ( Fig.10(down)). More distant from groyne, the velocity increased by maximum 30 % in the contraction part of the groyne cross-section. In the zone of recirculation flow, in wake of the groyne, the velocity decreased by -45% due to negative velocities in this area. In the crosswise direction the distribution of velocity components u x is shown in Fig. 11(top). The velocity profiles were taken in six different relative locations (groyne position, x rg =0) and in five locations along the channel (x rg =(x c -x g )/L g =-25, -1, +0.7, +1.6 and + 22). The essential growth of flow complexity near the location of the groyne is noticeable for profiles x rg =-1, 0, 0.7 and 1.6. Crosswise profiles of u y velocity component (not shown) explicitly confirm the two-dimensional properties of the flow field in the proximity of the groyne. ). As expected, far away from the groyne (x rg =-25 and +22) the flow is uniform with u y =0. It is visible that near the channel centreline the flow velocities (top) generally increase due the groyne presence (stream contraction effect), and, in thin regions near left bank and right bank they decrease and become negative in the circulation zone. In Fig.11 (down) dimensionless crosswise profiles of eddy viscosity are shown to make the picture of flow more complementary. The zig-zag shape of the EV profiles for x rg =-1 means that in close distance upstream from the groyne there are intensive variability of the turbulent viscosity (by 100 -150 % of reference level).

Case 2. Simulation with Parabolic Method for EV.
In Case 2 all runs of hydrodynamic model have been conducted with parabolic eddy viscosity model, given by eq. (7) in [10]. As in the first case, results of the same basic parameters (u x , u y , U, TSS, F r , EV) were analysed. Results of that simulation are shown in Fig. 12, both for velocity and eddy viscosity. Comparing them with results for mixed-length method (Fig.9), it is visible that parabolic method gives more complicated and enlarged backward-flow zone and it is very similar to the Shear Layer given in [3]. The computed field of eddy viscosity, shown in the lower part of Fig. 12 follows in great extent the velocity field and it is agreeable with the parabolic theoretical background [10]. Dimensionless quantities, shown in Fig.13 discover more detailed pattern of relative changes of water level and flow velocity magnitude near the groyne influence, more complicated in the wake region behind the groyne than in the Case 1. Fig. 13. Ddimensionless x-profiles of Z w and u. P model.
Crosswise profiles of dimensionless velocity, computed by using parabolic method, are shown in Fig. 14 (left), and, as compared to those for mixed-length method, they show bigger changes in central part of channel (up to 20%). But, in the recirculating area smaller changes in flow velocity are visible as compared to M-L method. Contrary to the EV profiles for Case 1, there are no zigzag shapes in the EV profiles close to the groyne, and, the calculated EV change in central part of channel is higher (20 -30 %) than for the case M-L model.

Case 3. Simulation with k-e method for EV.
In Case 3 all runs of hydrodynamic model have been conducted with k-e model, defined by eqs. (10 (a-c) in [10]. In the case of k-e model the solver required very small time steps (of 0.01 s). As in Case 1 and Case 2, results of the same basic flow parameters (u x , u y , U, TSS, F r , EV) were taken into account. Case 3 results are shown in Fig. 15 (velocity, EV). Comparing these fields with solutions obtained for mixed-length and parabolic solutions it becomes evident that k-e solution differs in great extent from the previous two cases. The maximum velocity area is far downstream from the groyne comparing the left-upper picture with its equivalents of Case1 (Fig. 9) and Case 2 (Fig.12). Even more sharp differences exist between the k-e EV field and its counterparts from Fig.9 and Fig.12. For flow depth the k-e solution demonstrates local rises and depressions (Fig.15,down). Dimensionless values of velocity change in this case (Fig.16, top) are similar in their ranges (2-9%) for the change obtained in Case 1 and Case 2, but they are characterized by sharp irregularities upstream to the groyne, and, in the region of groyne they resemble profiles obtained for mixed-length method (Fig.10). In this case variability of relative velocity is smaller near the groyne than the same change for parabolic case (Fig.  14) and for mixed-length (Fig. 13). In Fig.17 (top) are shown crosswise profiles of velocity components in both horizontal directions and it is clearly visible that the distribution of u x across the channel in the groyne area and downstream is completely different from those obtained for Case 1 and Case 2. The u x profiles across the channel resemble triangles with peak value near the central line of the channel. The specific distribution of u x across the channel is tightly connected with eddy viscosity field, as is shown in Fig.17 (down) where it is noticeable that the lowest levels of EV are present in the central part of channel and they rise substantially towards both banks.

Comparison of solutions for eddy viscosity models for flow near groyne
Comparison of solution for three EV Methods can be conducted by using dimensionless quantities. In Fig.18 crosswise profiles of dimensionless velocities (u x ) are presented in six different locations along the channel: x rg =x/L g =-25, -1, 0.7, 0.7, 1.4, 22. Vertical axes of the graphs are also expressed relatively to the groyne length (from 0 at right bank up to 6L g to the right bank). It can notice from the Fig.18 that starting from the profile similar to mixed-length solution in upstream part of the channel, the velocity profile calculated from k-e method differs greatly from the profiles obtained for mixed-length and parabolic methods and u x from k-e method are generally smaller (-100% till to -50%) and has completely different shape of triangle distribution with maximum in the centerline of the channel.
To check the performance of each model one can make tests about the fulfilment of the basic physics laws that should be meet accordingly to assumed mathematical model given in [10]. One of the criterion can be the mass balance, taken here as the balance of calculated total discharge in several cross-sections along the channel to the presumed inflow into the channel (upper boundary condition).The discharge at each test cross-section has been obtained from numerical integration of unit discharge field accordingly to the formula: (1) where: A -cross-section area, q -unit discharge [m 2 /s], q x,j -unit discharge component in x direction in j-th cell , y j -y-width of the j-th cell (j runs in y-direction).
Applying the integration procedure given by (1) for the same cross-sections as above have produced discharges that are gathered altogether in the Table 1. Data gathered in Table 1 show clearly that both mixedlength model and parabolic model conserve balance of flow very well. Deviations from 100% do not exceed 2.5% for both models and they probably are unavoidable errors due to numerical integration properties.
In the case of k-e model a rising loss of mass balance in downstream direction occurs (from 2.3% near upstream channel end up to 40% loss of the input discharge at downstream location far away from the groyne). Hydrodynamic model that uses the k-e model for closing its Reynolds equations simply does not preserve the basic physical law, mass conservation. Taking into account that all numerical experiments conducted in this research have utilized the same FEM network, used the same boundary and initial conditions and the same numerical parameters, it should be stated here that for such kind of flow modeling the k-e model should be rather avoided. There is no direct explanation of this peciular behaviour of k-e method. Some suggestions in literature state that for the flow-with-groyne the parabolic model over-predicts the strength of backwater flow, the mixed-length model provides better prediction for that flow, the k-e model under-predicts the length of the recirculation zone and produces viscosity field larger from any other turbulence models.

Conclusions
Three turbulence models (parabolic, mixed-length, k-e) were used for simulation flow in channel with a groyne. The hydrodynamic model has been verified with experimental data of the Rajaratnam's channel. Only k-e results substantially differed from measured values. Comparison between results of the three turbulence model showed that the flow velocity other flow parameters, predicted by parabolic model and mixedlength model are similar. Results of k-e model were significantly different from those produced by parabolic and mixing-length models. The balance of mass revealed that both parabolic model and mixed-length model fulfilled the criterion very well, contrary to the k-e solutions that suffered loss of it. The discrepancy in mass balance suggest using the most advanced turbulence model with great caution with deeper searching for reasons of such peculiar behaviour. 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.