Abstract

A numerical research with different turbulence models for shallow water equations was carried out. This was done in order to investigate which model has the ability to reproduce more accurately the wakes produced by the shock of the water hitting a submerged island inside a canal. The study of this phenomenon is important for the numerical methods application advancement in the simulation of free surface flows since these models involve a number of simplifications and assumptions that can have a significant impact on the numerical solutions quality and thus can not reproduce correctly the physical phenomenon. The numerical experiments were carried out on an experimental case under controlled conditions, consisting of a channel with a submerged conical island. The numerical scheme is based on the Eulerian-Lagrangian finite volume method with four turbulence models, three mixing lengths (ml), and one joining on the horizontal axis with a mixing-length model (ml) on the vertical axis. The experimental results show that a with ml turbulence model makes it possible to approach the experimental results in a more qualitative manner. We found that when using only a model in the vertical and horizontal direction, the numerical results overestimate the experimental data. Additionally the computing time is reduced by simplifying the turbulence model.

1. Introduction

In coastal and estuary regions large-scale vortices in the water occur due to the effects caused by the water hitting structures such as breakwaters and islands [1]. Studying these vortices is of great interest to environmental engineers because they have great influence on silt and other solutes transportation. For instance, the study of vortices influenced by effluents, and other thermal and radioactive waste originating released from hydroelectric or nuclear-electric plants into coastal waters is important for marine animal life [2, 3].

Attemping to understand water movement using shallow water equations means that the magnitude of the vertical scale is much smaller than the horizontal scale. In other words, if the body of water (river or lagoon) is located in a cartesian plane, the and scales are bigger than the scale (depth of the body of water). Nonetheless, it has been documented [4] that a more accurate simulation of vortices formation requires a 3D or multilayered approach to the body of water and that this is tightly connected to the turbulence model used. In this study we experimented with four turbulence models. The first three are mixing-length (ml) based, and the last one consists of joining a turbulence model with a mixing-length (ml) one; the results show that this last combination produces a more qualitative solution that is not overestimated compared to the one using only a in the horizontal and vertical axes. The comparison and validation of our numerical experiments are done against the experimental and numerical model proposed by [5, 6], which consists of a canal with an island in front of it to produce the vorticity in which the depth of the water is smaller than the height of the island whilst the flux is subcritical. The velocities in this controlled experiment are achieved by means of tracking particles using the velocity system known as PIV, which creates velocity field maps at water surface level, as well as measuring the velocity of water at specific points of interest (numerical viewers). These velocity maps provide information on the evolution of transport and the scales relative to the formation of vortices within the flux, and they are useful to validate computational models.

In research conducted by Lloyd and Stansby [5, 6], a number of controlled experiments consisting of a canal with several types of conical islands to produce wakes were carried out; the numerical results show that the best reproduction of the results was achieved using a 3D or multilayered model. However, a mixing-length turbulence model was used, producing numerical values below experimental ones. In a later study, Stansby and Lloyd's [6], using a semi-implicit multilayered method with a turbulence model, spurious or higher results were found compared to the experimental ones. Further, improved numerical methods of second-order accuracy resulted in a closer solution compared to the formation of vortices [6] and better numerical behavior in mixing-length models [4]. More recent described the limitations of 2D models integrated into the vertical axis to reproduce turbulence wakes [7].

In studies by Alamatian and Jaefarzadeh [8], several numerical experiments were conducted using the finite triangular volume method with three turbulence models: mixing length, , and an algebraic stress model, which were more accurate than experimental values. Similar studies by Akilli et al. [9] and Fu and Rockwell [10, 11] found that the model overestimates the speeds. Similarly, Sadeque et al. report similar results to those of Fu [12, 13].

In experiments by Rodi and Scheuerer [14], the model and a model consisting of a single equation were used to predict the adverse pressure gradients on border layers. The single equation model showed good results but the model produced systematic discrepancies, excessive friction due to coefficients especially in relatively simple fluxes. Lam and Bremhorst showed that some simplifications of the model provide better results for semilaminar and laminar fluxes in pipes compared to previous forms proposed by the model [15].

Rodri­guez et al. [16, 17] studied mixing-length models focusing on the shallow waters vortices formation dynamics. Some of the applications show the quality and the capabilities of the developed numerical model, and the mixing-length model reproduced the behavior of dispersion of pollutants in water patterns satisfactorily albeit with some simplifications [1, 18–20]. From the literature we can conclude that a 3D volume control model is appropriate to reproduce wakes accurately, and more accurately describe the solution of shallow waters equations [21, 22].

Barrios-Pina et al. [23] verified the applicability of a simpler turbulence model than the commonly used model to predict the mean flow through vegetation in shallow waters.

However, despite the great number of versions of turbulence models based mainly on mixing length, , and algebraic, none can be used universally for all kinds of problems, as the scale of the study domain and the computer equipment available must be considered.

2. Experimental Model Setup

Physical experiments against which the obtained numerical results using our model were compared against the numerical results obtained by Lloyd et al. [5, 24] using a canal in controlled conditions at the University of Manchester. The canal is shown in Figure 1 and its dimensions were 9.75 m × 1.52 m. An island built of aluminum was placed inside the channel; its size and dimensions are shown in Figure 2.

Figure 1 

Experimental model setup. Conical island inside the canal.

Figure 2 

Island dimensions (unscaled).

PIV was used to measure and the design was intended to measure velocities over relatively long sections, which is common in hydraulic experiments. The accuracy of PIV has been tested using a Doppler laser and the results show that the PIV has a margin of error below 10% [25].

The initial experimental conditions are an inflow to the canal at a velocity of 0.088 m/s and a depth of  m. As seen in Figure 2 the island is at a height of 0.066 m and thus remains at the same level as the initial water depth.

Figure 3 shows the location of numerical displays with the canal at a length of 0.73 m from the island and at 1.74 m. The physical display from which the data was obtained is therefore the only point where we compared numerical results against experimental results.

Figure 3 

Location of the viewers within the study domain.

3. Governing Equations

The equations upon which our computational implementation is based to model free surfaces flows are equations for shallow waters, where the velocities on the horizontal plane are considered to be greater than those on the vertical plane, and the scales of horizontal lengths are far greater than those of vertical lengths; it means . Thus, these are the equations that govern the flux are the following.

Equations for velocities on the horizontal plane are

Continuity equation is

Free surface equation is

The approach to shallow waters makes it acceptable not to solve the movement equation for the vertical speed, which is obtained from the continuity equation (2):

4. Turbulence Models

In this section the turbulence models used are described. Let us start with the mixing-length turbulence models used referred to as ml-03, ml-04, ml-05, and the hybrid + ml-06 model.

4.1. Mixing-Length Model ml-03

In the case of the ml-03 model taken from [4], we considered , where and (von Karman constant) and this is the vertical distance from the bottom to the surface of the free surface of the water. Thus, it is defined as where the distance between the wall to the point where the movement is being calculated. For the horizontal mixing, the following is the definition: where ; in general, this coefficient is calculated by means of experimental measurements. The turbulence diffusion coefficient is determined; thus,

Finally, the diffusion coefficient on the vertical plane is determined as and the diffusion coefficient on the horizontal plane is expressed as

4.2. Mixing-Length Model ml-04

In the case of the ml-04 mixing model taken from [26], let us take the vertical diffusion coefficient, just as we did in (7)   , and then change only the diffusion coefficient on the horizontal plane, calculated as where the chezy number is 55.33, is the vertical distance, and is the gravity.

4.3. Mixing-Length Model ml-05

In this mixing-length model the parameters used are the same as those of ml-01:

4.4. Hybrid with ml-06 Model

In this model the vertical diffusion coefficient is obtained using (12) from the ml-04 model and the horizontal diffusion coefficient is calculated as where , , , , , and is the production component defined as

5. Numerical Solution Scheme

The developed numerical model can work on three dimensions or on multilayers. The mesh used is a staggered cell. In Figure 4, the elements of the mesh are shown, where the black dot at the center of each cell represents any scale, such as temperature and salinity, and the circles on the faces indicate the velocity components position , and .

Figure 4 

Control volume discretization in 2D and 3D.

In terms of the vertical discretization, constant or variable layers can be considered. The variables distribution is shown in Figure 5.

Figure 5 

Discretization on the vertical plane.

The numerical solution is based on a spatial discretization of second order (finite differences). It is also considered a temporal Eulerian-Lagrangian formulation where the highly nonlinear movement equations (1) are resolved by using the characteristics method (Langrangian part), using a second order of accuracy, which has a smaller numerical diffusion level compared to other schemes in the same order, whilst it avoids numerical changes in models of higher orders caused by strong gradients. On the other hand, the diffusion components are resolved by means of an implicit scheme (Eulerian part). To approximate the free surface equation (3), it is necessary to solve a linear pentadiagonal (or heptadiagonal in 3D) system, where the Conjugated Gradient method is used. More details on the numerical solution algorithm and the computational model can be found in [16, 17, 19–21].

6. Numerical Experiments

The mesh used to conduct the numerical simulations consists of discrete control volumes sized m, m, and m. The number of cells is 336, 60, and 10, in , , and , respectively, which results in a total of cells in the computational domain (see Figure 6).

Figure 6 

Computational domain discretization.

Simulations were conducted with all the turbulence models described previously for 300 s of real simulation, that is, 300,000 iterations in time; a small  seg was selected to avoid the possibility of introducing numerical oscillations. Figure 7 shows the field speed vectors obtained using the and ml-06 turbulence model when the flux is stable.

Figure 7 

Velocity field vectors when the flow has been stabilized using the + ml-06 turbulence model.

Figure 7 shows large-scale turbulence structures present in the surface. In the region close to the island, a very slow speed region is formed, which results in the formation of horizontal speed layers due to the interaction of the accelerated flux region around the island with a velocity field close to zero. These layers can be clearly observed in Figure 14, where the bottom friction causes a reduction in velocity.

For a better perspective on the solution quality, let us analyse the simulated results against the experimental data obtained using the viewers shown in Figure 3. This is an important analysis because it deals with the numerical model validation, which consists in verifying the implementation using measurements obtained using the viewers. The interpretation of a hydraulic physics phenomenon involves measurements taken with instruments across different time and space scales. The dynamics of these experiments entail an analysis of these scales using specific methods, both statistical and deterministic. Thus, calibration and the validation of the numerical model are quite complex, because, even though the experiments are controlled in the laboratory, the devices and mechanisms used to measure have a range of uncertainty and carry an intrinsic degree of error with them. In addition, human error might also play a part, as interpretation and perception add extra uncertainty. To verify the quality of the numerical solution in relation to the observed data the Nash-Sutcliffe efficiency [3, 27, 28] is used, defined by the following expression: where is an observed or measured experimental data and is the data obtained through numerical simulation, at the same point in space and time. The parameter is the measure of the observed data. Obviously, if the value of is close to the unit (1.0), it is possible to consider that the simulation data is quite close to the experimental data; if this shows that the predictions are as accurate as the average of the observed data; when the average of the data is the best estimator for the model. Figure 8 shows experimental and simulated data obtained by Stansby that have been compared to numerical results obtained in this research.

Figure 8 

Comparing measurement data (+) and simulated data () of the velocity component on the viewer obtained by Stansby and Lloyd [6]. The last 90 seg of simulation are shown when the flux is considered to be stable.

We proceed to show the results of the relevant simulations for each of the turbulence models (ml-03, ml-04, ml-05, and + ml-06) obtained on the viewer , as well as the corresponding efficiency factor to compare the results of each model with experimental data. Figure 9 shows data obtained using numerical simulation compared to experimental results using the turbulence model ml-03. Similarly, Figure 10 shows data obtained using the ml-04 model, and Figure 11 shows results of the ml-05 model. Finally, Figure 12 shows data from model + ml-06. Data shown has been obtained from 200 seg of real simulation up to 300 seg. Only the last 90 seg of simulation when the flux is stable are shown. Table 1 summarizes the efficiency obtained with each of the turbulence models used.

Figure 9 

Comparing measurement data (+) and simulated data () of the velocity component on the viewer obtained using the ml-03 model.

Figure 10 

Comparing measurement data (+) and simulated data () of the velocity component on the viewer obtained using the ml-04 model.

Figure 11 

Comparing measurement data (+) and simulated data () of the velocity component on the viewer obtained using the ml-05 model.

Figure 12 

Comparing measurement data (+) and simulated data () of the velocity component on the viewer obtained using the + lm06.

Even though following efficiency criteria the closest model is ml-04, this model completely destroys vorticity and is thus a bad model to be used to simulate turbulence wakes; amongst the models that reproduce vorticity best, (Stansby) and + ml-06, the latter seems to be the most efficient because it improves the solution quality without destroying vorticity. For a better perspective on the numerical solution quality let us evaluate the ranges associated with each frequency using a Fourier analysis to determine the evolution of a velocity field over time. To this effect we analyzed the frequency domain signals, subtracting their average value so that the values vary around zero. The average values found are shown in Table 2.

Figure 13 shows the velocity Fourier velocity transformations norms, where similarities and differences between the signals present in the frequency function can be observed. As can be seen, the approach each signal shows to the experimental value in the frequency domain, for the lowest value, is within the acceptable range and that as the frequency increases the simulated signals move away from experimental results. Thus, any model within the Fourier domain is quite close to experimental results. In fact, if we calculate FFT and we establish the correlation to experimental signals all the numerical models are quite close to experimental results; in other words, when the numerical data increases so does the experimental model in a constant proportion (see Table 3).

Figure 13 

Velocity Fourier transformations norms.

(a)

(b)

(c)

(d)

(e)

(f)

(a)
(b)
(c)
(d)
(e)
(f)

Figure 14 

Velocity field on the , plane obtained at different depths; cross-section (a) shows the velocity field on the surface, cross-section (b) at 0.011 m, cross-section (c) at 0.023 m, cross-section (d) at 0.046 m, cross-section (e) at 0.58 cm, and (f) at the bottom.

Even though the best efficiency is shown by the ml-04 model (close to zero), and this indicates that the obtained results using this model are close to the average of the data observed in Figure 10, where the result was almost a straight line on average going through the velocity variation on , this shows that the model is efficient but it adequately reproduces the vorticity. Both an accurate efficiency and the adequate reproduction of the vorticity are desired. The model that meets both requirements is hybrid model + ml-06. The results of the reproduction of turbulence trails can be seen in Figure 14 and the experimental in Figure 19. Figure 15 shows the right combination of layers, with the speed field on the , plane half way through the canal, which indicates that the mixing-length model is working accurately on the vertical plane. The vorticity fileds obtained numerically and experimentally using the hybrid model are depicted in the Figures 17 and 18, respectively.

Figure 15 

Velocity field on the , plane and half way through the plane.

It is important to clarify that the problem of simulating turbulence wakes in a channel laboratory of relatively small dimensions means that the absolute values of the fluctuations of the free surface are too small in order to be analyzed. One of the main reasons is that if it was possible to measure experimentally these fluctuations, it would be necessary to consider the effect of the capillary that fluids yield when they are in contact with a solid surface, for this reason it is not meaningful to compare experimental measurements against the numerical values. In Figure 16 is depicted the free surface elevation at different simulation times.

(a)

(b)

(c)

(a)
(b)
(c)

Figure 16 

Distribution of free surface elevation with respect to a reference level (m); incise (a) at 350 seg, incise (b) at 400 seg, and incise (c) at 450 seg of real time simulation.

Figure 17 

Surface vorticity field numerically obtained.

Figure 18 

Experimental flow around island, surface vorticity field.

Figure 19 

Experimental surface velocity vector fields of wake of island.

7. Conclusions

A number of numerical experiments have been conducted to test various turbulence models for shallow fluxes with recirculation of wakes produced by an island with a slight slope. The obtained results were compared to experimental data obtained using a PTV and against a numerical model. The hybrid turbulence model proposed was found to reproduce the wakes formed by the clash of the flux with the island better than others. The mixing-length models showed better efficiency but could not reproduce vorticity like the model; that is, they cause more dissipation. We think that a 3D scheme using a mixing-length model comes quite close to the flux speeds produced by the wakes and that a model overestimates the speeds. We found that using a on the horizontal plane and mixing-length model on the horizontal place has a positive impact on the quality of the solution because it reproduces patterns of wakes and it produces similar characteristics in the speed change periods; in addition, it requires less computational time showing better efficiency than the model. Finally, the computational time using the hybrid turbulence model is lower by approximately 40% compared to the model and an algebraic model. Thus it can be used to reproduce turbulence patterns in industrial applications at a lower computational cost with an adequate reproduction of turbulence wakes.

Nomenclature

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors wish to thank Secretariat of Public Education (Mexico) specifically for the program (PROMED) of the professors improving the support in order to publish the results. Also, they are thankful to the Autonomous University of San Luis Potosí for giving them the facilities to develop the research in its laboratories and to CONACYT-ESIA-IPN for the grant given to the second author.