Abstract
The standard model has the deficiency of predicting swirling and vortical flows due to its isotropic assumption of eddy viscosity. In this study, a second-order nonlinear model is developed incorporating some new functions for the model coefficients to explore the models applicability to complex turbulent flows. Considering the realizability principle, the coefficient of eddy viscosity () is derived as a function of strain and rotation parameters. The coefficients of nonlinear quadratic term are estimated considering the anisotropy of turbulence in a simple shear layer. Analytical solutions for the fundamental properties of swirl jet are derived based on the nonlinear model, and the values of model constants are determined by tuning their values for the best-fitted comparison with the experiments. The model performance is examined for two test cases: (i) for an ideal vortex (Stuart vortex), the basic equations are solved numerically to predict the turbulent structures at the vortex center and the (ii) unsteady 3D simulation is carried out to calculate the flow field of a compound channel. It is observed that the proposed nonlinear model can successfully predict the turbulent structures at vortex center, while the standard model fails. The model is found to be capable of accounting the effect of transverse momentum transfer in the compound channel through generating the horizontal vortices at the interface.
1. Introduction
It is well known that the RANS (Reynolds Averaged Navier Stokes) type turbulence models, such as two-equation model or Reynolds stress model, are the most popular tool used for practical engineering applications [1], because it requires less CPU time and computer memory compared to LES and DNS. Therefore, the clarification of the possibility, the limitation, and areas of improvement of RANS models should be still paid attention to. To resolve the Reynolds stress term that appeared in the averaged Navier-Stokes equations, the model is the most frequently adopted one. However, the standard model cannot produce satisfactory results for the flow field having high rate of strain and rotation because of its isotropic assumption of eddy viscosity [2]. On the other hand, a nonlinear model can predict superior result by capturing the anisotropy of turbulence.
The nonlinear model is a generalized eddy viscosity model, where additional nonlinear terms of mean strain rate are added in the Reynolds stress equation. Generally, in a model, the value of the coefficient of eddy viscosity () is assumed to be constant throughout the turbulent flow field that overpredicts the value of eddy viscosity, especially in the case of large rate of strain and rotation. If the strain is sufficiently large, the model may produce negative normal stresses. Moreover, comparing the nonlinear quadratic term in the Reynolds stress equation proposed by Yoshizawa [3] with its other forms, it is inferred that is a function of strain and rotation parameters (shown later). The nonlinear model used in previous studies (Craft et al. [4], Cotton and Ismail [5], Kato and Launder [6], and Kimura and Hosoda [7]) mainly considers only strain rate, and the effect of rotation rate is neglected. In this study a modified and generalized nonlinear model is proposed, where both the strain and rotation parameters are taken into account.
The content of this paper can be classified into twofold. In the first part, the model development and evaluation of model constants are explained considering the analytical solution of some basic turbulent flows. In the second part, the model performance is examined considering the predictability of (i) turbulent structures at the center of an ideal vortex (2D simulation) and (ii) flow field in a compound channel (unsteady 3D simulation).
The model constants are initially estimated considering the realizability conditions and the anisotropy of turbulence in simple shear flows. Then, the approximate solutions are derived for the fundamental properties of a swirl jet. Neglecting the swirl parameters from the derived solutions, the turbulent properties are also calculated for a round jet without swirl. The model constants are finally determined by tuning their values for best fitting the approximate solutions with the previous experimental results.
The numerical simulation is carried out to solve the turbulent field of an ideal vortex street using standard and nonlinear models. The Stuart vortex is considered for this analysis. The model performance is evaluated considering the predictability of turbulent structures at vortex center.
In a compound channel, large-scale horizontal vortices are generated at the interface of main channel and flood plain due to velocity gradient. The vortex formation causes momentum transfer in the lateral direction and forms a complex flow field. Using the proposed model, unsteady 3D simulation was carried out, and the models capability is examined considering both the averaged velocity and secondary flow predictions.
2. Mathematical Formulation
2.1. Governing Equations
The basic equations in a model for an unsteady incompressible flow are as follows. Continuity equation: momentum equation: -equation: -equation: where is the spatial coordinates, and are the average and turbulent velocities, respectively, in direction, is the pressure, is the density of fluid, is the kinematic viscosity, is the averaged turbulent energy, is the averaged turbulent energy dissipation rate, and is the eddy viscosity. , , , and are the model constants; standard values (, , , and ) are used for these constants.
2.2. Turbulence Model
(a) Standard Model. In the standard model, the Reynolds stress tensor is solved by linear constitutive equations derived from Boussinesq eddy viscosity concept, which does not take into account the anisotropy effect: is determined from the dimensional consideration of and and is approximated by Here, bears a constant value of 0.09.
(b) Nonlinear Model. Including the nonlinear anisotropy term in the Reynolds stress equation introduced by Yoshizawa [3], the constitutive equation can be expressed in the following form: Here, (= , , ) is the coefficient of nonlinear quadratic term. In this equation, is also determined by (6), but is not a constant. Its functional form is derived as follows.
The quadratic term in (7) is equivalent to the following nonlinear viscosity model: where the strain and rotation tensors are defined as From (7) and (9), the relations between the coefficients of two forms of nonlinear terms can be derived as The relation inferred that the coefficient of eddy viscosity is a function of strain and rotation parameters. The strain parameter () and rotation parameter (Ω) are defined in the following equation:
Comparing the analytical results for diagonal components of the anisotropic tensors with those of experiments for simple shear flows, it is observed that the functional form of gave better results instead of taking their constant values (shown later). Therefore, the coefficient of quadratic term is also considered as a function of strain and rotation parameters.
Therefore, the model used in this study differs from the standard model in two important ways: (i) the eddy viscosity coefficient () is not a constant but a function of strain and rotation parameters and (ii) nonlinear terms are added in the Reynolds stress equation to account the anisotropy of turbulence. Many kinds of model functions have been proposed for coefficient . Most of them considered only strain parameter and rotation parameter is neglected (such as Cotton and Ismail [5] and Kato and Launder [6]). Craft et al. [4] considered one dominant parameter of two (either or ). In this case, the value of eddy viscosity changes suddenly around the region of and a sharp edge is observed in 2D profile of in plane. This feature seems to be physically unsound. A typical flow which satisfies is a simple shear flow. Such model is therefore not suitable for flows with vortex formation from a simple shear layer due to instability. Authors have proposed more generalized functional form for these coefficients considering the effect of both parameters. The proposed functional forms are as follows: Here, , , , , , , , , , and are the model constants. More simple functional forms of used in the algebraic stress model can be obtained from the above equation just neglecting some of the terms. Moreover, when the strain and rotation effects are neglected, that is, , becomes equal to the standard value of 0.09. Neglecting the quadratic term, the model reduces to the standard model.
If the effect of strain and rotation parameters on is neglected, (14) gives . The values of used in this study are given as follows:
3. Estimation of Model Constants
The values of model constants are given in Table 1, which are obtained through some theoretical considerations as explained below. It is a tedious job to determine 10 model constants. That is why, firstly, the model is applied to plain shear layer, where , and therefore the model constants are reduced to 4.
3.1. Tuning of Coefficients in
The coefficients of nonlinear quadratic term, (, , ) in (6), should be determined carefully because they are expected to influence the physical accuracy and numerical performance of the model. In this study, these model constants are tuned considering the anisotropy of turbulent normal stresses in simple shear flows.
For simple shear flow, the flow can be described as The strain and rotation parameters are defined as The diagonal components of the nondimensional Reynolds stress tensors are The anisotropic tensor are defined by The ratio of the turbulent energy production term to the dissipation rate is The comparison between the experimental results and analytical solutions for the anisotropy of turbulent normal stresses () in plane shear layer is shown in Figure 1 ( = the ratio of turbulent production and dissipation rate). In the figure, CHC and HGC denote the experimental results by Champagne et al. [10] and Harris et al. [11], respectively. The bold line indicates the functional form for (14), which gives better comparison with experiments than (15). Equation (15) shows the values of bearing constant values, where the effects of strain and rotation parameters are neglected.
3.2. Consideration of Realizability for Plane Shear Layer
Realizability can be defined as the requirement of the nonnegativity of turbulent normal stresses and Schwarz inequality between any turbulent velocity correlations. It is a basic physical and mathematical principle that the solution of any turbulence model equation should obey [12]. The realizability inequalities for 3D turbulent flows are as follows: Einstein’s summation rule is not applied in the above equations. In a two-dimensional averaged flow, (22) coincides with (23). In this study, the restrictions on are derived from the mentioned realizability equations for simple shear flow.
Applying (21) to plane shear layer, the following two equations are derived for the diagonal components of the Reynolds stress tensor (nonnegativity condition): For plane shear layer, .
Applying (22), the following inequality equation can be derived for Reynolds stress component, (Schwarz inequality condition): Since the value of is positive and is negative, (24) is satisfied regardless of . Thus, the restrictions on , derived from (25) and (26), are as follows: The model constants in the functional form of in (13) are tuned to satisfy the realizability conditions derived in (27).
For plane shear layer, the realizability conditions (27) as well as the proposed functional form of (16) are plotted in Figure 2. The calculations are made with the following values of model constants: Using these estimated values, the approximate solutions for swirl jet are compared with the previous experimental results and the values of model constants are finally determined by tuning their values for the best-fitted comparison. Table 1 shows the final values of constants obtained by such a trial and error method, and Figure 2 confirms that the model obeys the realizability conditions with these values of constants.
In the log-law region, the assumed functional form of shows almost a constant value of 0.09. It can be noted that, instead of functional form, if a constant value of () is used throughout the turbulent flow field, it fails to satisfy the realizability conditions. Figure 3 shows the distribution of assumed functional form of in plane.
4. Analytical Solution of Swirling Jet
To obtain the basic equations of turbulent shear flows, the following assumptions are made.(a)The pressure gradient is small; that is, .(b)Viscous shear stress is much smaller than the turbulent shear stress and can be neglected.(c)Diffusion in the direction normal to the flow (- and -coordinates) is much larger than the diffusion in the direction of flow (-coordinate).Thus, the momentum equation in -direction as well as the and equations can be presented in a simplified form as follows. The momentum equation in -direction: -equation: -equation:
Figure 4 shows the definition sketch of a swirl jet with Cartesian and Cylindrical coordinate systems. Consider that , , and are the jet velocities in axial (), lateral (), and transverse () directions in Cartesian coordinate system and , , and are the velocities in axial (), radial (), and azimuthal () directions in Cylindrical coordinate system, respectively.
4.1. Assumed Profiles
In the self-similar region, the following relations can be obtained for the attenuation rate of hydraulic variables of an axisymmetric swirl jet [13, 14]: Here, is the jet width; is the centerline maximum velocity of the jet; and are the centerline values of turbulent kinetic energy and its dissipation rate, respectively.
The following assumptions are made for the functional forms of velocity and distributions, which are compatible with the decaying power law of velocity and along the centerline of jets: where , , , , and are the unknown coefficients of the velocity profile and , , , and are those of the distributions.
4.2. Swirl Parameter
To define the swirl parameter, consider the Cylindrical coordinate system. Rajaratnam [13] showed that, for the circular jet with swirl, integral of the pressure plus momentum () and angular momentum () are preserved. They are defined by respectively.
Considering , (38) can be reduced to Substituting the assumed velocity distributions into (39) and (40), the following relations are obtained ( is calculated from and ): A nondimensional parameter combining and is used to represent the relative amount of swirl present in the flow; it is called swirl number (): where is the radius of the jet nozzle.
Using (41) and (42), the swirl number can be defined as Hence,
4.3. Derived Solutions
Substituting the assumed velocity distributions from (33) to (35) into the continuity equation, the following algebraic relations are derived:
The integral equation for conservation of momentum flux also results in the following equation: where , , and are the area, velocity, and initial momentum flux of the jet at inlet, respectively.
Substituting (45a) and (45b) into (46), the coefficient of attenuation of radial velocity “” is determined as a function of inlet conditions:
For the assumed velocity and distributions, using the Reynolds equation in the -direction and equations, the following three algebraic expressions are derived as the relations of the lowest order with respect to the power of . Reynolds equation in -direction: -equation: -equation:
The values of , , and are determined by (45a), (45b), and (47). Substituting the values of and , which are the and values at centerline, the development rate of swirl jet () can be estimated by (48). For any known value of swirl number , the value of “” can be determined using (44). The coefficients of and distributions are determined by solving (49) and (50).
The radial distributions of turbulent intensities as well as turbulent shear stresses are derived by constitutive equations.
4.4. Results
Tuning the model constants (estimated from realizability considerations) the approximate solutions are compared with the previous experimental results and the values of model constants are determined for the best-fitted comparison. Table 1 shows the obtained values of model constants. Eliminating swirl parameter, that is, taking swirl number () and hence the swirl parameter () as zero, the derived solutions for swirl jet are applied to a round jet without swirl. The comparisons of the results for swirl and nonswirl jets with experiments are presented below.
(i) Round Jet without Swirl. The self-similar velocity profile predicted by the nonlinear model is compared with the experimental data of Wygnanski and Fielder [15] in Figure 5. In Figure 6, the normalized centerline velocity decay along the centerline is compared with the experimental result reported in Pope [16]. The predictions agree well with the previous experimental studies. Figures 7, 8, and 9 show the radial distributions of turbulent intensities for axial, radial, and circumferential velocities. The results are compared with the experimental data by Wygnanski and Fielder [15] and Wang and Law [17]. In Figure 10, the calculated shear stress profile is compared with the range of experimental results by Fukushima et al. [18], where two experimental profiles represent the lower and the upper boundary of a number of measured profiles for different downstream distances.
The distribution of turbulent kinetic energy () normalized by is favorably compared with the experiment of Wygnanski and Fielder and is shown in Figure 11. The approximate solution for normalized turbulent energy dissipate rate (), compared with the experimental data of Wygnanski and Fielder as well as with that of Sami [19], is shown in Figure 12. The available experimental data on turbulent kinetic energy are quite scattered. The present prediction coincides with the Wygnanski and Fielder data near the central region of jet, and it is close to Sami’s data for higher radial distances. In comparison to experimental results, a little deviation is observed in the spreading of profile; it can be noted that this deficiency is due to the simple assumptions made in the exponential part of distributions in (37). The values of and , which are the centerline values of and , are estimated by tuning their values for the best-fitted comparison with experimental results.
Approximate solutions are also derived by using standard model and shown in the figures by dotted line.
(ii) Swirl Jet. The turbulent intensities for axial, lateral, and circumferential velocities are calculated for a moderate swirl number of . The results are compared with the experimental data by Pratte and Keffer [8] in Figures 13, 14 and 15. In Figure 16, the distribution of normalized turbulent kinetic energy (where is twice the turbulent kinetic energy) is favorably compared with the previous experiment. Figure 17 shows the radial distribution of turbulent Reynolds stress, ; due to unavailability of experimental data, this result is not compared with the experiments.
The figures show that the nonlinear model is well compared with the experiments results. In Figure 15, the measured data appeared to show a peak of turbulent intensity of the circumferential velocity at the jet edge (~0.1) which is not observed in the calculated results. This deficiency may be due to neglecting the higher order terms in the analytical solution for the sake of simplicity. In comparison to experimental results as well as nonlinear model prediction, the standard model underpredicts the turbulent intensities for axial direction and overpredicts for radial and circumferential directions. The discrepancy is higher for swirl jet case. Reasonably, the nonlinear model showed better comparison than standard one.
The distributions of turbulent intensities are found to be changed significantly depending on swirl number. For different swirl numbers of , 1.5, and 2.5, the turbulent intensities are calculated. The radial distributions of turbulent intensities for axial, radial, and circumferential velocities are shown in Figures 18, 19, and 20, respectively. It is observed that, with the increasing swirl number, the turbulent intensity is decreased for axial velocity component and increased for lateral and circumferential velocity components.
5. Turbulent Structures in Stuart Vortex
Figure 21 shows an example of the spanwise cut of coherent structure, in particular in the plane mixing layer (Hussain [9]). The outer contour of coherent vorticity denotes the structure boundary. It shows that there are two critical points in the structure: the saddle () characterized by negligible spanwise vorticity and the center of vortex () characterized by peak spanwise vorticity. Based on the experimental results on turbulent structures of a plane shear layer, Hussain [9] reported that the structures of turbulent normal stresses at vortex point are elliptical in shape; on the other hand, the turbulent shear stresses show hyperbolic profile. The same patterns are also confirmed in further studies with turbulent axisymmetric jets and wakes. In this section, 2D numerical simulations were carried out for the turbulent characteristics of an ideal vortex street (Stuart vortices) using the proposed nonlinear model; the predicted turbulent structures for vortex center are compared with those of Hussain [9] and the model’s performance is assessed. The results of the standard model will also be presented for comparison.
(a)
(b)
The equation for the stream function of Stuart vortex can be expressed as follows: Here, is constant and indicates the eccentricity of the elliptic streamline of the vortex. In this study, is used. For this velocity field, the and equations (3) and (4) are solved with the finite volume method based on a staggered grid system. Turbulent stresses are calculated using both standard and nonlinear model equations expressed in (5) and (7), respectively. As shown in Figure 22, two vortices are considered in the flow domain with a periodic boundary condition at upstream and downstream end.
The periodic turbulent structures simulated by nonlinear model are shown in Figure 23. It is observed that the turbulent kinetic energy (), dissipation rate (), and turbulent normal stresses in , , and directions (expressed as , , and , resp.) show the contours in elliptical shape at the vortex center. On the other hand, the turbulent shear stress in plane () shows hyperbolic profile (saddle pattern). The results of turbulent structural topology at vortex center are found compatible with those of previous experimental investigations [9] of coherent structures in free shear flows.
(a)
(b)
(c)
(d)
(e)
(f)
However, in the standard model, coefficient has no dependency on the rate of strain and rotation and bears a constant value (= 0.09) throughout the turbulent flow field. This deficiency of standard model causes an inconsistent prediction of turbulent structures at the center of vortex. Figure 24 shows the results simulated by standard model. The model failed to generate the elliptical profiles (focus pattern) for , , and turbulent normal stresses. Changing the values of model constants in (13), it is also observed that the turbulent structures at vortex center are sensitive to the functional form of .
(a)
(b)
(c)
(d)
(e)
(f)
6. 3D Unsteady Simulation of Compound Channel
In a compound channel, a shear layer is developed at the interface of the main channel and flood plain due to velocity gradient and large-scale vortices are generated with the vertical axis. The vortex formation causes momentum transfer in the lateral direction from the main channel to the flood plains and reduces the channel conveyance, and the flow resistance in the channel is increased. The vortices also significantly distort the cross-sectional velocity profile and secondary flow pattern. Therefore, for modeling the flow in a compound channel, the model should be capable of accounting the effect of horizontal vortices to predict the flow field accurately. Since the secondary current is generated, due to anisotropy of turbulence, the model constants have high influence on it. In this section, 3D unsteady simulation is carried out for compound channel, and the model performance is examined comparing the simulated results with the previous experiments.
6.1. Flow Configuration and Computational Domain
The governing equations for mean velocities and turbulent flows are discretized with the finite volume method based on a staggered grid system. For the momentum equation, convective and diffusive fluxes are approximated with Quick and central difference schemes, respectively. The hybrid central upwind scheme is used for and equations. Time advancement is achieved by Adam-Bashforth scheme of second-order accuracy in each equation. The basic equations are discretized as fully explicit forms and are solved successively with the time increment step by step. The pressure field is solved using iterative procedure at each time step. The free surface elevation is solved by continuity equation integrated over the control volume of the surface layer. The wall functions are employed as the wall boundary conditions for and . The frictions near the bed and side walls are estimated by log-law. Periodic boundary conditions are used in upstream and downstream ends of the flow domain.
The experimental results reported by Knight [20] are used for the comparison of the simulated result. The cross section of the flow domain for single flood plain is shown in Figure 25, and the hydraulic conditions are shown in Table 2. The computational domain consists of 200 grids in longitudinal (streamwise, ), 50 in transverse (widthwise, ), and 11 in vertical (depthwise, ) directions, for single flood plain case. For double flood plain case, 72 grids are used in transverse direction.
6.2. Results
The calculated depth averaged velocity profile is compared with the laboratory experiments of Knight [20]. In Figure 26, Cal. (1) represents the predicted profile with the proposed nonlinear model at time = 300 sec. At this stage, the horizontal vortices at the interface of main channel and flood plain are fully developed. On the other hand, Cal. (2) is the profile where the effect of horizontal vortices is neglected. It is observed that the velocity profile is highly influenced by the horizontal vortices, especially near the interface region. The comparison shows that the well-agreed profiles are obtained when the effects of vortices are considered. Therefore, the model is capable of accounting the effect of transverse momentum transfer in the compound channel through generating the horizontal vortices at the interface.
Since the pattern of secondary current is influenced by the model constants, the simulated secondary flow is compared with experiments of Knight [20] in Figure 27. It shows that the patterns as well as the magnitude of simulated secondary flow (Figure 27(b)) are well comparable with the experiments (Figure 27(a)).
(a) Experiment (Knight, 1991)
(b) Simulation
7. Conclusion
In this study, a second-order nonlinear model is proposed incorporating some new functions for the model coefficients to explore the model applicability to complex turbulent flows. The model constants are determined considering some basic turbulent flows, and the model performance is examined considering the predictability of simulated results for (i) turbulent structures at the center of an ideal vortex (Stuart vortex) and (ii) the flow field of a compound open channel.
The coefficient of eddy viscosity () and the coefficient of nonlinear quadratic term (, , and ) are derived as a function of strain and rotation parameters. The model constants are estimated considering the realizability conditions and the anisotropy of turbulence in a simple shear layer. Approximate solutions for the fundamental properties of swirl and nonswirl jets are derived based on the nonlinear model. Tuning the model constants (estimated from realizability considerations), the approximate solutions are compared with the previous experimental results, and the values of model constants are determined for the best-fitted comparison. The solutions of standard model are also presented. Reasonably, the nonlinear model showed better comparison than the standard one.
From the simulated results of an ideal vortex street, it is observed that the present nonlinear model can successfully predict the turbulent structures at vortex center. However, the standard model failed to generate the elliptical structures of turbulent normal stresses at vortex center observed in previous experimental studies.
For the compound channel, it is observed that the model is capable of accounting the transverse momentum transfer in the compound channel through generating the horizontal vortices at the interface. The velocity profile is found to be highly influenced by the horizontal vortices, and the well-agreed profiles are obtained when the effects of vortices are considered. The patterns and the magnitude of simulated secondary flow are also found to be well comparable with experiments.
Notation of Symbols
| : | Area of the jet at inlet |
| , , : | Coefficients of the assumed velocity profiles |
| : | Jet width |
| , : | model constants |
| , : | Nonlinear model constants |
| , , , , , , , , : | Model constants for |
| : | The model constants for |
| : | Turbulent energy |
| , : | Coefficients of the assumed distribution |
| : | Initial momentum flux of the jet |
| , , , : | Model constants for |
| : | Pressure |
| : | Radius of the jet nozzle |
| : | Strain parameter defined by (15) |
| : | Swirl number |
| : | Strain rate tensors |
| : | Velocity of the jet at inlet |
| : | Centerline maximum velocity of the jet |
| : | Average velocity in direction |
| , , : | Velocities in Cartesian coordinate system |
| , , : | Velocities in Cylindrical coordinate system |
| : | Turbulent velocity in direction |
| : | Reynolds stress tensor |
| : | Spatial coordinates |
| , , : | Axial, lateral, and transverse directions in Cartesian coordinate system |
| , , : | Axial, radial, and azimuthal directions in Cylindrical coordinate system |
| , : | Coefficients of the assumed velocity profiles |
| : | Turbulent energy dissipation rate |
| , : | Coefficients of the assumed distribution |
| : | Density of fluid |
| : | Molecular dynamic viscosity |
| : | Eddy viscosity |
| , : | model constants |
| : | Rotation parameter defined by (12) |
| : | Rotation rate tensor. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.