Development of a Nonlinear k-ε Model Incorporating Strain and Rotation Parameters for Prediction of Complex Turbulent Flows

The standard k-ε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 k-ε 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 (c μ ) 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 k-ε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 k-εmodel can successfully predict the turbulent structures at vortex center, while the standard k-εmodel fails. Themodel 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.


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 2 International Journal of Partial Differential Equations 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.

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.  ,   ,  1 , and  2 are the model constants; standard values (  = 1.0,   = 1.3,  1 = 1.44, and  2 = 1.92) are used for these constants.

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,   (=  1 ,  2 ,  3 ) 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,   ,  Ω ,   ,  Ω ,  Ω ,  1 ,  Ω1 ,  Ω1 ,   , 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,  = Ω = 0,   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   =   0 .The values of  0 used in this study are given as follows:

Estimation of Model Constants
The values of model constants are given in 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.

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 Eq.( 15) Eq. ( 14) 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,  1  2 (Schwarz inequality condition): Since the value of  1 is positive and  3 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   (= 0.09) 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.

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, / ≈ 0.
(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.

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]: International Journal of Partial Differential Equations Here,  is the jet width;   is the centerline maximum velocity of the jet;  0 and  0 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  1 ,  3 , , , and  are the unknown coefficients of the velocity profile and  0 ,  2 ,  0 , and  2 are those of the - distributions.

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 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  0 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  0 ,  0 , and  0 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 1/.
Reynolds equation in -direction: -equation:  -equation: The values of ,  1 , and  3 are determined by (45a), (45b), and (47).Substituting the values of  0 and  0 , 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.

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   2 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 (/  3 ), 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 International Journal of Partial Differential Equations turbulent kinetic energy ( 2 ) 1/2 /  (where  2 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  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.The equation for the stream function of Stuart vortex can be expressed as follows:

Turbulent Structures in Stuart Vortex
Here, 0 ≥  ≥ −1 is constant and indicates the eccentricity of the elliptic streamline of the vortex.In this study,  = −0.5 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  1  1 ,  2  2 , and  3  3 , resp.) show the contours in elliptical shape at the vortex center.On the other hand, the turbulent shear stress in  plane ( 1  2 ) 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.
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   .

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.the interface region.The comparison shows that the wellagreed 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

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 ( 1 ,  2 , and  3 ) 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 wellagreed 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.

Figure 2 :Figure 3 :
Figure 2: Relation between   and  in a simple shear layer.

Figure 5 :
Figure 5: Comparison of self-similar mean velocity profile.

Figure 6 :
Figure 6: Centerline mean velocity along the axial distance.

Figure 7 :
Figure 7: Distribution of turbulent intensity of the axial velocity in a round jet without swirl.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :
Figure 8: Distribution of turbulent intensity of the radial velocity in a round jet without swirl.

Figure 12 :
Figure 12: Radial distribution of turbulent kinetic energy dissipation rate in a round jet without swirl.

Figure 13 :
Figure 13: Distribution of turbulent intensity of the axial velocity in a swirl jet.

Figure 14 :Figure 15 :Figure 16 :Figure 17 :
Figure 14: Distribution of turbulent intensity of the radial velocity in a swirl jet.

Figure 21 Figure 18 :Figure 19 :
Figure21shows 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]

InternationalFigure 20 :Figure 21 :
Figure 20: Distribution of turbulent intensity of the circumferential velocity for different swirl numbers (nonlinear - model).

Figure 27 :
Figure 27: Comparison of spatial averaged secondary current in a compound open channel.

𝐴 0 :
Area of the jet at inlet  1 ,  3 , : Coefficients of the assumed velocity profiles : J e t w i d t h  1 ,  2 : - model constants   ,   ( 1 ,  2 ,  3 ): N o n l i n e a r- model constants  0 ,   ,  Ω ,   ,  Ω ,  Ω ,  1 ,  Ω1 ,  Ω1 : Model constants for    0 : Th e m o d e l c o n s t

Table 1 :
Values for the coefficients of   and   .

Table 2 :
Hydraulic conditions for numerical simulation of compound channel flow.
[20]n 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 Figure26, 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 International Journal of Partial Differential Equations a n t s f o r   : Turbulent energy  0 ,  2 : Coefficients of the assumed  distribution  0 : Initial momentum flux of the jet   ,  Ω ,   ,  Ω : Average velocity in   direction   ,   ,   : Velocities in Cartesian coordinate system   ,   ,   : Velocities in Cylindrical coordinate system   : T u r b u l e n t v e l o c i t y i n   direction     : Reynolds stress tensor : S p a t i a l c o o r d i n a t e s , , : Axial, lateral, and transverse directions in Cartesian coordinate system