# Development of a Nonlineark

Authors: Md. Shahjahan Ali; Takashi Hosoda; Ichiro Kimura
Journal: International Journal of Partial Differential Equations
 (2015)

---

## Abstract

The standardk
-
ε 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. 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.


---

## Body

## 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 k
-
ε model is the most frequently adopted one. However, the standard k
-
ε 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 nonlineark
-
ε 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 k
-
ε model, the value of the coefficient of eddy viscosity (c
μ) 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 c
μ 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 k
-
ε 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 ak
-
ε model for an unsteady incompressible flow are as follows.Continuity equation:(1)
∂
U
i
∂
x
i
=
0
,
momentum equation:(2)
∂
U
i
∂
t
+
∂
U
j
U
i
∂
x
j
=
g
i
-
1
ρ
∂
P
∂
x
i
+
∂
∂
x
j
-
u
i
u
j
¯
+
ν
∂
2
U
i
∂
x
j
2
,
k-equation:
(3)
∂
k
∂
t
+
∂
k
U
j
∂
x
j
=
-
u
i
u
j
¯
∂
U
i
∂
x
j
+
∂
∂
x
j
ν
t
σ
k
+
ν
∂
k
∂
x
j
-
ε
,
ɛ-equation:
(4)
∂
ε
∂
t
+
∂
ε
U
j
∂
x
j
=
-
c
ε
1
ε
k
u
i
u
j
¯
∂
U
i
∂
x
j
+
∂
∂
x
j
ν
t
σ
ε
+
ν
∂
ε
∂
x
j
-
c
ε
2
ε
2
k
,where x
i is the spatial coordinates, U
i and u
iare the average and turbulent velocities, respectively, in x
i direction, P is the pressure, ρ is the density of fluid, ν is the kinematic viscosity, k is the averaged turbulent energy, ɛ is the averaged turbulent energy dissipation rate, and ν
t is the eddy viscosity. σ
k, σ
ε,  c
ε
1, and c
ε
2 are the model constants; standard values (σ
k
=
1.0, σ
ε
=
1.3, c
ε
1
=
1.44, and c
ε
2
=
1.92) are used for these constants.


### 2.2. Turbulence Model

(a) Standard
k
-
ε
Model. In the standard k
-
ε model, the Reynolds stress tensor u
i
u
j
¯ is solved by linear constitutive equations derived from Boussinesq eddy viscosity concept, which does not take into account the anisotropy effect:
(5)
-
u
i
u
j
¯
=
ν
t
S
i
j
-
2
3
k
δ
i
j
,
S
i
j
=
∂
U
i
∂
x
j
+
∂
U
j
∂
x
i
.
ν
t is determined from the dimensional consideration of k and ɛ and is approximated by
(6)
ν
t
=
c
μ
k
2
ε
.
Here, c
μ bears a constant value of 0.09.(b) Nonlinear
k
-
ε
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:
(7)
-
u
i
u
j
¯
=
ν
t
S
i
j
-
2
3
k
δ
i
j
-
k
ε
v
t
∑
β
=
1
3
c
β
S
β
i
j
-
1
3
S
β
α
α
δ
i
j
,
(8)
S
1
i
j
=
∂
U
i
∂
x
γ
∂
U
j
∂
x
γ
,
S
2
i
j
=
1
2
∂
U
γ
∂
x
i
∂
U
j
∂
x
γ
+
∂
U
γ
∂
x
j
∂
U
i
∂
x
γ
,
S
3
i
j
=
∂
U
γ
∂
x
i
∂
U
γ
∂
x
j
.
Here, c
β (= c
1, c
2, c
3) is the coefficient of nonlinear quadratic term. In this equation, ν
t is also determined by (6), but c
μ is not a constant. Its functional form is derived as follows.The quadratic term in (7) is equivalent to the following nonlinear viscosity model:
(9)
α
1
S
i
l
Ω
l
j
+
Ω
l
i
S
j
l
+
α
2
S
i
l
S
l
j
-
1
3
S
k
m
S
m
k
δ
i
j
+
α
3
Ω
i
l
Ω
l
j
-
1
3
Ω
k
m
Ω
m
k
δ
i
j
,
where the strain and rotation tensors are defined as
(10)
S
i
j
=
∂
U
i
∂
x
j
+
∂
U
j
∂
x
i
,
Ω
i
j
=
∂
U
i
∂
x
j
-
∂
U
j
∂
x
i
.
From (7) and (9), the relations between the coefficients of two forms of nonlinear terms can be derived as
(11)
c
1
=
-
2
α
1
+
α
2
-
α
3
,
c
2
=
2
(
α
2
+
α
3
)
,
c
3
=
2
α
1
+
α
2
-
α
3
.
The relation inferred that the coefficient of eddy viscosity c
μ is a function of strain and rotation parameters. The strain parameter (S) and rotation parameter (Ω) are defined in the following equation:
(12)
S
=
k
ε
1
2
S
i
j
S
i
j
,
Ω
=
k
ε
1
2
Ω
i
j
Ω
i
j
.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 ofc
β gave better results instead of taking their constant values (shown later). Therefore, the coefficient of quadratic term c
β is also considered as a function of strain and rotation parameters.Therefore, the model used in this study differs from the standardk
-
ε model in two important ways: (i) the eddy viscosity coefficient (c
μ) 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 c
μ. 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 S or Ω). In this case, the value of eddy viscosity changes suddenly around the region of S
=
Ω and a sharp edge is observed in 2D profile of c
μ in S
-
Ω plane. This feature seems to be physically unsound. A typical flow which satisfies S
=
Ω is a simple shear flow. Such model is therefore not suitable for flows with vortex formation from a simple shear layer due to K
-
H instability. Authors have proposed more generalized functional form for these coefficients considering the effect of both parameters. The proposed functional forms are as follows:
(13)
c
μ
=
c
μ
0
1
+
c
n
s
S
2
+
c
n
Ω
Ω
2
·
1
+
c
d
s
S
2
+
c
d
Ω
Ω
2
+
c
d
s
Ω
S
Ω
+
c
d
s
1
S
4
+
c
d
Ω
1
Ω
4
+
c
d
s
Ω
1
S
2
Ω
2
h
h
h
+
c
d
Ω
1
Ω
4
+
c
d
s
Ω
1
S
2
Ω
2
-
1
,
(14)
c
β
=
c
β
0
1
+
m
d
s
S
2
+
m
d
Ω
Ω
2
,
β
=
1
,
2
,
3
.
Here, c
n
s, c
n
Ω, c
d
s, c
d
Ω, c
d
s
Ω, c
d
s
1, c
d
Ω
1, c
d
s
Ω
1, m
d
s, and m
d
Ω are the model constants. More simple functional forms of c
μ 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, S
=
Ω
=
0, c
μ becomes equal to the standard value of 0.09. Neglecting the quadratic term, the model reduces to the standard k
-
ε model.If the effect of strain and rotation parameters onc
β is neglected, (14) gives c
β
=
c
β
0. The values of c
β
0 used in this study are given as follows:
(15)
c
1
0
=
0.40
,
c
2
0
=
0
,
c
3
0
=
-
0.13
.


## 2.1. Governing Equations

The basic equations in ak
-
ε model for an unsteady incompressible flow are as follows.Continuity equation:(1)
∂
U
i
∂
x
i
=
0
,
momentum equation:(2)
∂
U
i
∂
t
+
∂
U
j
U
i
∂
x
j
=
g
i
-
1
ρ
∂
P
∂
x
i
+
∂
∂
x
j
-
u
i
u
j
¯
+
ν
∂
2
U
i
∂
x
j
2
,
k-equation:
(3)
∂
k
∂
t
+
∂
k
U
j
∂
x
j
=
-
u
i
u
j
¯
∂
U
i
∂
x
j
+
∂
∂
x
j
ν
t
σ
k
+
ν
∂
k
∂
x
j
-
ε
,
ɛ-equation:
(4)
∂
ε
∂
t
+
∂
ε
U
j
∂
x
j
=
-
c
ε
1
ε
k
u
i
u
j
¯
∂
U
i
∂
x
j
+
∂
∂
x
j
ν
t
σ
ε
+
ν
∂
ε
∂
x
j
-
c
ε
2
ε
2
k
,where x
i is the spatial coordinates, U
i and u
iare the average and turbulent velocities, respectively, in x
i direction, P is the pressure, ρ is the density of fluid, ν is the kinematic viscosity, k is the averaged turbulent energy, ɛ is the averaged turbulent energy dissipation rate, and ν
t is the eddy viscosity. σ
k, σ
ε,  c
ε
1, and c
ε
2 are the model constants; standard values (σ
k
=
1.0, σ
ε
=
1.3, c
ε
1
=
1.44, and c
ε
2
=
1.92) are used for these constants.


## 2.2. Turbulence Model

(a) Standard
k
-
ε
Model. In the standard k
-
ε model, the Reynolds stress tensor u
i
u
j
¯ is solved by linear constitutive equations derived from Boussinesq eddy viscosity concept, which does not take into account the anisotropy effect:
(5)
-
u
i
u
j
¯
=
ν
t
S
i
j
-
2
3
k
δ
i
j
,
S
i
j
=
∂
U
i
∂
x
j
+
∂
U
j
∂
x
i
.
ν
t is determined from the dimensional consideration of k and ɛ and is approximated by
(6)
ν
t
=
c
μ
k
2
ε
.
Here, c
μ bears a constant value of 0.09.(b) Nonlinear
k
-
ε
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:
(7)
-
u
i
u
j
¯
=
ν
t
S
i
j
-
2
3
k
δ
i
j
-
k
ε
v
t
∑
β
=
1
3
c
β
S
β
i
j
-
1
3
S
β
α
α
δ
i
j
,
(8)
S
1
i
j
=
∂
U
i
∂
x
γ
∂
U
j
∂
x
γ
,
S
2
i
j
=
1
2
∂
U
γ
∂
x
i
∂
U
j
∂
x
γ
+
∂
U
γ
∂
x
j
∂
U
i
∂
x
γ
,
S
3
i
j
=
∂
U
γ
∂
x
i
∂
U
γ
∂
x
j
.
Here, c
β (= c
1, c
2, c
3) is the coefficient of nonlinear quadratic term. In this equation, ν
t is also determined by (6), but c
μ is not a constant. Its functional form is derived as follows.The quadratic term in (7) is equivalent to the following nonlinear viscosity model:
(9)
α
1
S
i
l
Ω
l
j
+
Ω
l
i
S
j
l
+
α
2
S
i
l
S
l
j
-
1
3
S
k
m
S
m
k
δ
i
j
+
α
3
Ω
i
l
Ω
l
j
-
1
3
Ω
k
m
Ω
m
k
δ
i
j
,
where the strain and rotation tensors are defined as
(10)
S
i
j
=
∂
U
i
∂
x
j
+
∂
U
j
∂
x
i
,
Ω
i
j
=
∂
U
i
∂
x
j
-
∂
U
j
∂
x
i
.
From (7) and (9), the relations between the coefficients of two forms of nonlinear terms can be derived as
(11)
c
1
=
-
2
α
1
+
α
2
-
α
3
,
c
2
=
2
(
α
2
+
α
3
)
,
c
3
=
2
α
1
+
α
2
-
α
3
.
The relation inferred that the coefficient of eddy viscosity c
μ is a function of strain and rotation parameters. The strain parameter (S) and rotation parameter (Ω) are defined in the following equation:
(12)
S
=
k
ε
1
2
S
i
j
S
i
j
,
Ω
=
k
ε
1
2
Ω
i
j
Ω
i
j
.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 ofc
β gave better results instead of taking their constant values (shown later). Therefore, the coefficient of quadratic term c
β is also considered as a function of strain and rotation parameters.Therefore, the model used in this study differs from the standardk
-
ε model in two important ways: (i) the eddy viscosity coefficient (c
μ) 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 c
μ. 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 S or Ω). In this case, the value of eddy viscosity changes suddenly around the region of S
=
Ω and a sharp edge is observed in 2D profile of c
μ in S
-
Ω plane. This feature seems to be physically unsound. A typical flow which satisfies S
=
Ω is a simple shear flow. Such model is therefore not suitable for flows with vortex formation from a simple shear layer due to K
-
H instability. Authors have proposed more generalized functional form for these coefficients considering the effect of both parameters. The proposed functional forms are as follows:
(13)
c
μ
=
c
μ
0
1
+
c
n
s
S
2
+
c
n
Ω
Ω
2
·
1
+
c
d
s
S
2
+
c
d
Ω
Ω
2
+
c
d
s
Ω
S
Ω
+
c
d
s
1
S
4
+
c
d
Ω
1
Ω
4
+
c
d
s
Ω
1
S
2
Ω
2
h
h
h
+
c
d
Ω
1
Ω
4
+
c
d
s
Ω
1
S
2
Ω
2
-
1
,
(14)
c
β
=
c
β
0
1
+
m
d
s
S
2
+
m
d
Ω
Ω
2
,
β
=
1
,
2
,
3
.
Here, c
n
s, c
n
Ω, c
d
s, c
d
Ω, c
d
s
Ω, c
d
s
1, c
d
Ω
1, c
d
s
Ω
1, m
d
s, and m
d
Ω are the model constants. More simple functional forms of c
μ 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, S
=
Ω
=
0, c
μ becomes equal to the standard value of 0.09. Neglecting the quadratic term, the model reduces to the standard k
-
ε model.If the effect of strain and rotation parameters onc
β is neglected, (14) gives c
β
=
c
β
0. The values of c
β
0 used in this study are given as follows:
(15)
c
1
0
=
0.40
,
c
2
0
=
0
,
c
3
0
=
-
0.13
.


## 3. Estimation of Model Constants

The values of model constants are given in Table1, 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 S
=
Ω, and therefore the model constants are reduced to 4.Table 1
Values for the coefficients ofc
μ and c
β.
Model const.
c
μ
0
c
n
s
c
n
Ω
c
d
s
c
d
s
Ω
c
d
Ω
c
d
s
1
c
d
Ω
1
c
d
s
Ω
1
m
d
s
m
d
Ω
Values
0.09
0.005
0.0068
0.008
−0.003
0.004
0.00005
0.00005
0.00025
0.01
0.003


### 3.1. Tuning of Coefficients inc
β

The coefficients of nonlinear quadratic term,c
β  (c
1, c
2, c
3) 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(16)
U
1
=
U
1
y
,
U
2
=
U
3
=
0
,
d
U
1
d
y
>
0
.
The strain and rotation parameters are defined as
(17)
S
=
Ω
=
M
=
k
ε
d
U
1
d
y
.
The diagonal components of the nondimensional Reynolds stress tensors are
(18)
u
1
u
1
¯
k
=
2
3
+
c
μ
2
c
1
-
c
3
3
M
2
,
u
2
u
2
¯
k
=
2
3
+
c
μ
2
c
3
-
c
1
3
M
2
,
u
3
u
3
¯
k
=
2
3
-
c
μ
c
3
+
c
1
3
M
2
.
The anisotropic tensor m
i
j are defined by
(19)
m
i
j
=
u
i
u
j
¯
k
-
2
3
δ
i
j
.
The ratio of the turbulent energy production term P
k to the dissipation rate ε is
(20)
α
=
P
k
ε
=
-
u
1
u
2
¯
ε
d
U
1
d
y
=
c
μ
M
2
.
The comparison between the experimental results and analytical solutions for the anisotropy of turbulent normal stresses (m
i
i) 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 c
β  (14), which gives better comparison with experiments than (15). Equation (15) shows the values of c
β bearing constant values, where the effects of strain and rotation parameters are neglected.Figure 1
Comparison of anisotropy of turbulence using (14) and (15) for coefficient c
β.


### 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:
(21)
u
i
u
i
¯
≥
0
,
(22)
u
i
u
i
¯
·
u
j
u
j
¯
≥
u
i
u
j
¯
2
i
≠
j
,
(23)
det
⁡
u
1
u
1
¯
u
1
u
2
¯
u
1
u
3
¯
u
2
u
1
¯
u
2
u
2
¯
u
2
u
3
¯
u
3
u
1
¯
u
3
u
2
¯
u
3
u
3
¯
≥
0
.
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 c
μ 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):
(24)
u
1
u
1
¯
k
=
2
3
+
c
μ
2
c
1
-
c
3
3
M
2
≥
0
,
(25)
u
2
u
2
¯
k
=
2
3
+
c
μ
2
c
3
-
c
1
3
M
2
≥
0
.
For plane shear layer, M
=
S
=
Ω.Applying (22), the following inequality equation can be derived for Reynolds stress component, u
1
u
2
¯ (Schwarz inequality condition):
(26)
c
μ
2
9
M
2
+
2
c
1
2
+
2
c
3
2
-
5
c
1
c
3
M
4
-
2
c
μ
c
1
+
c
3
M
2
-
4
≤
0
.
Since the value of c
1 is positive and c
3 is negative, (24) is satisfied regardless of M. Thus, the restrictions on c
μ, derived from (25) and (26), are as follows:
(27)
c
μ
≤
2
c
1
-
2
c
3
M
2
,
c
μ
≤
(
c
1
+
c
3
)
M
c
1
+
c
3
2
M
2
+
4
9
+
(
2
c
1
2
+
2
c
3
2
-
5
c
1
c
3
)
M
2
h
+
c
1
+
c
3
2
M
2
+
4
9
+
(
2
c
1
2
+
2
c
3
2
-
5
c
1
c
3
)
M
2
·
9
M
+
2
c
1
2
+
2
c
3
2
-
5
c
1
c
3
M
3
-
1
.
The model constants in the functional form of c
μ 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 c
μ  (16) are plotted in Figure 2. The calculations are made with the following values of model constants:
(28)
c
μ
0
=
0.09
,
c
n
s
+
c
n
Ω
=
0.0118
,
c
d
s
+
c
d
Ω
+
c
d
s
Ω
=
0.009
,
c
d
s
1
+
c
d
Ω
1
+
c
d
s
Ω
1
=
0.00035
,
m
d
s
+
m
d
Ω
=
0.013
.
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.Figure 2
Relation betweenc
μ and M in a simple shear layer.In the log-law region, the assumed functional form ofc
μ shows almost a constant value of 0.09. It can be noted that, instead of functional form, if a constant value of c
μ (=
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 c
μ in S
-
Ω plane.Figure 3
Profile ofc
μ (S
,
Ω) on an S
-
Ω plane (13).


## 3.1. Tuning of Coefficients inc
β

The coefficients of nonlinear quadratic term,c
β  (c
1, c
2, c
3) 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(16)
U
1
=
U
1
y
,
U
2
=
U
3
=
0
,
d
U
1
d
y
>
0
.
The strain and rotation parameters are defined as
(17)
S
=
Ω
=
M
=
k
ε
d
U
1
d
y
.
The diagonal components of the nondimensional Reynolds stress tensors are
(18)
u
1
u
1
¯
k
=
2
3
+
c
μ
2
c
1
-
c
3
3
M
2
,
u
2
u
2
¯
k
=
2
3
+
c
μ
2
c
3
-
c
1
3
M
2
,
u
3
u
3
¯
k
=
2
3
-
c
μ
c
3
+
c
1
3
M
2
.
The anisotropic tensor m
i
j are defined by
(19)
m
i
j
=
u
i
u
j
¯
k
-
2
3
δ
i
j
.
The ratio of the turbulent energy production term P
k to the dissipation rate ε is
(20)
α
=
P
k
ε
=
-
u
1
u
2
¯
ε
d
U
1
d
y
=
c
μ
M
2
.
The comparison between the experimental results and analytical solutions for the anisotropy of turbulent normal stresses (m
i
i) 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 c
β  (14), which gives better comparison with experiments than (15). Equation (15) shows the values of c
β bearing constant values, where the effects of strain and rotation parameters are neglected.Figure 1
Comparison of anisotropy of turbulence using (14) and (15) for coefficient c
β.


## 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:
(21)
u
i
u
i
¯
≥
0
,
(22)
u
i
u
i
¯
·
u
j
u
j
¯
≥
u
i
u
j
¯
2
i
≠
j
,
(23)
det
⁡
u
1
u
1
¯
u
1
u
2
¯
u
1
u
3
¯
u
2
u
1
¯
u
2
u
2
¯
u
2
u
3
¯
u
3
u
1
¯
u
3
u
2
¯
u
3
u
3
¯
≥
0
.
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 c
μ 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):
(24)
u
1
u
1
¯
k
=
2
3
+
c
μ
2
c
1
-
c
3
3
M
2
≥
0
,
(25)
u
2
u
2
¯
k
=
2
3
+
c
μ
2
c
3
-
c
1
3
M
2
≥
0
.
For plane shear layer, M
=
S
=
Ω.Applying (22), the following inequality equation can be derived for Reynolds stress component, u
1
u
2
¯ (Schwarz inequality condition):
(26)
c
μ
2
9
M
2
+
2
c
1
2
+
2
c
3
2
-
5
c
1
c
3
M
4
-
2
c
μ
c
1
+
c
3
M
2
-
4
≤
0
.
Since the value of c
1 is positive and c
3 is negative, (24) is satisfied regardless of M. Thus, the restrictions on c
μ, derived from (25) and (26), are as follows:
(27)
c
μ
≤
2
c
1
-
2
c
3
M
2
,
c
μ
≤
(
c
1
+
c
3
)
M
c
1
+
c
3
2
M
2
+
4
9
+
(
2
c
1
2
+
2
c
3
2
-
5
c
1
c
3
)
M
2
h
+
c
1
+
c
3
2
M
2
+
4
9
+
(
2
c
1
2
+
2
c
3
2
-
5
c
1
c
3
)
M
2
·
9
M
+
2
c
1
2
+
2
c
3
2
-
5
c
1
c
3
M
3
-
1
.
The model constants in the functional form of c
μ 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 c
μ  (16) are plotted in Figure 2. The calculations are made with the following values of model constants:
(28)
c
μ
0
=
0.09
,
c
n
s
+
c
n
Ω
=
0.0118
,
c
d
s
+
c
d
Ω
+
c
d
s
Ω
=
0.009
,
c
d
s
1
+
c
d
Ω
1
+
c
d
s
Ω
1
=
0.00035
,
m
d
s
+
m
d
Ω
=
0.013
.
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.Figure 2
Relation betweenc
μ and M in a simple shear layer.In the log-law region, the assumed functional form ofc
μ shows almost a constant value of 0.09. It can be noted that, instead of functional form, if a constant value of c
μ (=
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 c
μ in S
-
Ω plane.Figure 3
Profile ofc
μ (S
,
Ω) on an S
-
Ω plane (13).


## 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,∂
P
/
∂
x
≈
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 (y- and z-coordinates) is much larger than the diffusion in the direction of flow (x-coordinate).Thus, the momentum equation in x-direction as well as the k and ε equations can be presented in a simplified form as follows.The momentum equation inx-direction:
(29)
U
x
∂
U
x
∂
x
+
U
y
∂
U
x
∂
y
+
U
z
∂
U
x
∂
z
=
∂
∂
x
-
u
x
u
x
¯
+
∂
∂
y
-
u
x
u
y
¯
+
∂
∂
z
-
u
x
u
z
¯
,
k-equation:
(30)
U
x
∂
k
∂
x
+
U
y
∂
k
∂
y
+
U
z
∂
k
∂
z
=
-
u
x
u
x
¯
∂
U
x
∂
x
-
u
x
u
y
¯
∂
U
x
∂
y
-
u
x
u
z
¯
∂
U
x
∂
z
-
u
y
u
x
¯
∂
U
y
∂
x
-
u
y
u
y
¯
∂
U
y
∂
y
-
u
y
u
z
¯
∂
U
y
∂
z
-
u
z
u
x
¯
∂
U
z
∂
x
-
u
z
u
y
¯
∂
U
z
∂
y
-
u
z
u
z
¯
∂
U
z
∂
z
+
∂
∂
y
ν
t
σ
k
∂
k
∂
y
+
∂
∂
z
ν
t
σ
k
∂
k
∂
z
-
ε
,
ε-equation:
(31)
U
x
∂
ε
∂
x
+
U
y
∂
ε
∂
y
+
U
z
∂
ε
∂
z
=
c
ε
1
ε
k
-
u
x
u
x
¯
∂
U
x
∂
x
-
u
x
u
y
¯
∂
U
x
∂
y
-
u
x
u
z
¯
∂
U
x
∂
z
h
h
h
h
h
h
h
-
u
y
u
x
¯
∂
U
y
∂
x
-
u
y
u
y
¯
∂
U
y
∂
y
-
u
y
u
z
¯
∂
U
y
∂
z
h
h
h
h
h
h
h
-
u
z
u
x
¯
∂
U
z
∂
x
-
u
z
u
y
¯
∂
U
z
∂
y
-
u
z
u
z
¯
∂
U
z
∂
z
+
∂
∂
y
ν
t
σ
ε
∂
k
∂
y
+
∂
∂
z
ν
t
σ
ε
∂
k
∂
z
-
c
ε
2
ε
2
k
.Figure4 shows the definition sketch of a swirl jet with Cartesian and Cylindrical coordinate systems. Consider that U
x, U
y, and U
z are the jet velocities in axial (x), lateral (y), and transverse (z) directions in Cartesian coordinate system and U
x, U
r, and U
θ are the velocities in axial (x), radial (r), and azimuthal (θ) directions in Cylindrical coordinate system, respectively.Figure 4
Definition sketch of swirl jet (Pratte and Keffer [8]).


### 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]:
(32)
B
∝
x
,
U
m
∝
1
x
,
k
0
∝
1
x
2
,
ε
0
∝
1
x
4
.
Here, B is the jet width; U
m is the centerline maximum velocity of the jet; k
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 andk
-
ε distributions, which are compatible with the decaying power law of velocity and k
-
ε along the centerline of jets:
(33)
U
x
=
a
1
x
exp
⁡
-
y
2
+
z
2
γ
2
x
2
+
a
3
x
y
2
+
z
2
γ
2
x
2
exp
⁡
-
y
2
+
z
2
γ
2
x
2
,
(34)
U
y
=
b
x
y
γ
x
exp
⁡
-
y
2
+
z
2
γ
2
x
2
+
α
z
γ
x
3
exp
⁡
-
y
2
+
z
2
γ
2
x
2
,
(35)
U
z
=
b
x
z
γ
x
exp
⁡
-
y
2
+
z
2
γ
2
x
2
-
α
y
γ
x
3
exp
⁡
-
y
2
+
z
2
γ
2
x
2
,
(36)
k
=
1
x
2
k
0
+
k
2
y
2
+
z
2
γ
2
x
2
exp
⁡
-
2
y
2
+
z
2
γ
2
x
2
,
(37)
ε
=
1
x
4
ε
0
+
ε
2
y
2
+
z
2
γ
2
x
2
exp
⁡
-
3
y
2
+
z
2
γ
2
x
2
,
where a
1, a
3, b, α, and γ are the unknown coefficients of the velocity profile and k
0, k
2, ε
0, and ε
2 are those of the k
-
ε 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 (W) and angular momentum (T) are preserved. They are defined by
(38)
W
=
∫
0
∞
U
x
2
-
U
θ
2
2
r
d
r
,
(39)
T
=
∫
0
∞
r
2
U
x
U
θ
d
r
,
respectively.ConsideringU
x
≫
U
θ,  (38) can be reduced to
(40)
W
=
∫
0
∞
U
x
2
r
d
r
.
Substituting the assumed velocity distributions into (39) and (40), the following relations are obtained (U
θ is calculated from U
y and U
z):
(41)
W
=
γ
2
8
2
a
1
2
+
2
a
1
a
3
+
a
3
2
,
T
=
α
γ
3
8
a
1
+
a
3
.
A nondimensional parameter combining W and T is used to represent the relative amount of swirl present in the flow; it is called swirl number (S
N):
(42)
S
N
=
T
W
r
0
,
where r
0 is the radius of the jet nozzle.Using (41) and (42), the swirl number can be defined as
(43)
S
N
=
6
α
γ
2
25
b
r
0
.
Hence,
(44)
α
=
25
S
N
b
r
0
6
γ
2
.


### 4.3. Derived Solutions

Substituting the assumed velocity distributions from (33) to (35) into the continuity equation, the following algebraic relations are derived:(45a)
a
1
=
2
b
γ
,
(45b)
a
3
=
2
3
b
γ
.The integral equation for conservation of momentum flux also results in the following equation:(46)
M
0
=
U
0
2
A
0
=
2
π
a
1
2
γ
2
4
+
2
a
1
a
3
γ
2
8
+
a
3
2
γ
2
8
,
where A
0, U
0, and M
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 “b” is determined as a function of inlet conditions:
(47)
b
=
3
5
M
0
π
.For the assumed velocity andk
-
ε distributions, using the Reynolds equation in the x-direction and k
-
ε equations, the following three algebraic expressions are derived as the relations of the lowest order with respect to the power of 1
/
x.Reynolds equation inx-direction:
(48)
1
x
3
:
-
a
1
2
ε
0
2
+
1
2
c
d
s
k
0
2
4
a
1
2
+
8
b
2
γ
2
=
2
c
μ
0
ε
0
k
ε
0
2
+
1
2
c
n
s
k
0
2
ε
0
2
4
a
1
2
+
8
b
2
γ
2
·
2
a
3
γ
2
-
2
a
1
γ
2
-
2
b
γ
0
2
,
k-equation:
(49)
1
x
4
:
-
2
a
1
k
0
+
ε
0
ε
0
2
+
1
2
c
d
s
k
0
2
4
a
1
2
+
8
b
2
γ
2
=
c
μ
0
σ
k
ε
0
k
0
2
4
k
2
γ
2
-
8
k
0
γ
2
·
ε
0
2
+
1
2
c
n
s
k
0
2
ε
0
2
4
a
1
2
+
8
b
2
γ
2
,
ε-equation:
(50)
1
x
6
:
-
4
a
1
ε
0
k
0
+
c
ε
2
ε
0
2
ε
0
2
+
1
2
c
d
s
k
0
2
4
a
1
2
+
8
b
2
γ
2
=
c
μ
0
σ
ε
ε
0
k
0
3
4
ε
2
γ
2
-
8
ε
0
γ
2
ε
0
2
+
1
2
c
n
s
k
0
2
ε
0
2
4
a
1
2
+
8
b
2
γ
2
.The values ofb, a
1, and a
3 are determined by (45a), (45b), and (47). Substituting the values of k
0 and ε
0, which are the k and ε values at centerline, the development rate of swirl jet (γ) can be estimated by (48). For any known value of swirl number S
N, the value of “α” can be determined using (44). The coefficients of k 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. Table1 shows the obtained values of model constants. Eliminating swirl parameter, that is, taking swirl number (S
N) 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 k
-
ε 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.Figure 5
Comparison of self-similar mean velocity profile.Figure 6
Centerline mean velocity along the axial distance.Figure 7
Distribution of turbulent intensity of the axial velocity in a round jet without swirl.Figure 8
Distribution of turbulent intensity of the radial velocity in a round jet without swirl.Figure 9
Distribution of turbulent intensity of the circumferential velocity in a round jet without swirl.Figure 10
Radial distribution of turbulent shear stress in a round jet without swirl.The distribution of turbulent kinetic energy (k) normalized by U
m
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 (ε
x
/
U
m
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 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 k
0 and ε
0, which are the centerline values of k and ε, are estimated by tuning their values for the best-fitted comparison with experimental results.Figure 11
Radial distribution of turbulent kinetic energy in a round jet without swirl.Figure 12
Radial distribution of turbulent kinetic energy dissipation rate in a round jet without swirl.Approximate solutions are also derived by using standardk
-
ε 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 S
N
=
0.3. 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 q
2
¯
1
/
2
/
U
m (where q
2
¯ is twice the turbulent kinetic energy) is favorably compared with the previous experiment. Figure 17 shows the radial distribution of turbulent Reynolds stress, u
x
u
r
¯; due to unavailability of experimental data, this result is not compared with the experiments.Figure 13
Distribution of turbulent intensity of the axial velocity in a swirl jet.Figure 14
Distribution of turbulent intensity of the radial velocity in a swirl jet.Figure 15
Distribution of turbulent intensity of the circumferential velocity in a swirl jet.Figure 16
Radial distribution of turbulent kinetic energy in a swirl jet.Figure 17
Radial distribution of turbulent shear stress in a swirl jet.The figures show that the nonlinear model is well compared with the experiments results. In Figure15, the measured data appeared to show a peak of turbulent intensity of the circumferential velocity at the jet edge (r
/
x~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 k
-
ε 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 ofS
N
=
0.3, 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.Figure 18
Distribution of turbulent intensity of the axial velocity for different swirl numbers (nonlineark
-
ε model).Figure 19
Distribution of turbulent intensity of the radial velocity for different swirl numbers (nonlineark
-
ε model).Figure 20
Distribution of turbulent intensity of the circumferential velocity for different swirl numbers (nonlineark
-
ε model).


## 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]:
(32)
B
∝
x
,
U
m
∝
1
x
,
k
0
∝
1
x
2
,
ε
0
∝
1
x
4
.
Here, B is the jet width; U
m is the centerline maximum velocity of the jet; k
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 andk
-
ε distributions, which are compatible with the decaying power law of velocity and k
-
ε along the centerline of jets:
(33)
U
x
=
a
1
x
exp
⁡
-
y
2
+
z
2
γ
2
x
2
+
a
3
x
y
2
+
z
2
γ
2
x
2
exp
⁡
-
y
2
+
z
2
γ
2
x
2
,
(34)
U
y
=
b
x
y
γ
x
exp
⁡
-
y
2
+
z
2
γ
2
x
2
+
α
z
γ
x
3
exp
⁡
-
y
2
+
z
2
γ
2
x
2
,
(35)
U
z
=
b
x
z
γ
x
exp
⁡
-
y
2
+
z
2
γ
2
x
2
-
α
y
γ
x
3
exp
⁡
-
y
2
+
z
2
γ
2
x
2
,
(36)
k
=
1
x
2
k
0
+
k
2
y
2
+
z
2
γ
2
x
2
exp
⁡
-
2
y
2
+
z
2
γ
2
x
2
,
(37)
ε
=
1
x
4
ε
0
+
ε
2
y
2
+
z
2
γ
2
x
2
exp
⁡
-
3
y
2
+
z
2
γ
2
x
2
,
where a
1, a
3, b, α, and γ are the unknown coefficients of the velocity profile and k
0, k
2, ε
0, and ε
2 are those of the k
-
ε 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 (W) and angular momentum (T) are preserved. They are defined by
(38)
W
=
∫
0
∞
U
x
2
-
U
θ
2
2
r
d
r
,
(39)
T
=
∫
0
∞
r
2
U
x
U
θ
d
r
,
respectively.ConsideringU
x
≫
U
θ,  (38) can be reduced to
(40)
W
=
∫
0
∞
U
x
2
r
d
r
.
Substituting the assumed velocity distributions into (39) and (40), the following relations are obtained (U
θ is calculated from U
y and U
z):
(41)
W
=
γ
2
8
2
a
1
2
+
2
a
1
a
3
+
a
3
2
,
T
=
α
γ
3
8
a
1
+
a
3
.
A nondimensional parameter combining W and T is used to represent the relative amount of swirl present in the flow; it is called swirl number (S
N):
(42)
S
N
=
T
W
r
0
,
where r
0 is the radius of the jet nozzle.Using (41) and (42), the swirl number can be defined as
(43)
S
N
=
6
α
γ
2
25
b
r
0
.
Hence,
(44)
α
=
25
S
N
b
r
0
6
γ
2
.


## 4.3. Derived Solutions

Substituting the assumed velocity distributions from (33) to (35) into the continuity equation, the following algebraic relations are derived:(45a)
a
1
=
2
b
γ
,
(45b)
a
3
=
2
3
b
γ
.The integral equation for conservation of momentum flux also results in the following equation:(46)
M
0
=
U
0
2
A
0
=
2
π
a
1
2
γ
2
4
+
2
a
1
a
3
γ
2
8
+
a
3
2
γ
2
8
,
where A
0, U
0, and M
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 “b” is determined as a function of inlet conditions:
(47)
b
=
3
5
M
0
π
.For the assumed velocity andk
-
ε distributions, using the Reynolds equation in the x-direction and k
-
ε equations, the following three algebraic expressions are derived as the relations of the lowest order with respect to the power of 1
/
x.Reynolds equation inx-direction:
(48)
1
x
3
:
-
a
1
2
ε
0
2
+
1
2
c
d
s
k
0
2
4
a
1
2
+
8
b
2
γ
2
=
2
c
μ
0
ε
0
k
ε
0
2
+
1
2
c
n
s
k
0
2
ε
0
2
4
a
1
2
+
8
b
2
γ
2
·
2
a
3
γ
2
-
2
a
1
γ
2
-
2
b
γ
0
2
,
k-equation:
(49)
1
x
4
:
-
2
a
1
k
0
+
ε
0
ε
0
2
+
1
2
c
d
s
k
0
2
4
a
1
2
+
8
b
2
γ
2
=
c
μ
0
σ
k
ε
0
k
0
2
4
k
2
γ
2
-
8
k
0
γ
2
·
ε
0
2
+
1
2
c
n
s
k
0
2
ε
0
2
4
a
1
2
+
8
b
2
γ
2
,
ε-equation:
(50)
1
x
6
:
-
4
a
1
ε
0
k
0
+
c
ε
2
ε
0
2
ε
0
2
+
1
2
c
d
s
k
0
2
4
a
1
2
+
8
b
2
γ
2
=
c
μ
0
σ
ε
ε
0
k
0
3
4
ε
2
γ
2
-
8
ε
0
γ
2
ε
0
2
+
1
2
c
n
s
k
0
2
ε
0
2
4
a
1
2
+
8
b
2
γ
2
.The values ofb, a
1, and a
3 are determined by (45a), (45b), and (47). Substituting the values of k
0 and ε
0, which are the k and ε values at centerline, the development rate of swirl jet (γ) can be estimated by (48). For any known value of swirl number S
N, the value of “α” can be determined using (44). The coefficients of k 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. Table1 shows the obtained values of model constants. Eliminating swirl parameter, that is, taking swirl number (S
N) 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 k
-
ε 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.Figure 5
Comparison of self-similar mean velocity profile.Figure 6
Centerline mean velocity along the axial distance.Figure 7
Distribution of turbulent intensity of the axial velocity in a round jet without swirl.Figure 8
Distribution of turbulent intensity of the radial velocity in a round jet without swirl.Figure 9
Distribution of turbulent intensity of the circumferential velocity in a round jet without swirl.Figure 10
Radial distribution of turbulent shear stress in a round jet without swirl.The distribution of turbulent kinetic energy (k) normalized by U
m
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 (ε
x
/
U
m
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 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 k
0 and ε
0, which are the centerline values of k and ε, are estimated by tuning their values for the best-fitted comparison with experimental results.Figure 11
Radial distribution of turbulent kinetic energy in a round jet without swirl.Figure 12
Radial distribution of turbulent kinetic energy dissipation rate in a round jet without swirl.Approximate solutions are also derived by using standardk
-
ε 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 S
N
=
0.3. 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 q
2
¯
1
/
2
/
U
m (where q
2
¯ is twice the turbulent kinetic energy) is favorably compared with the previous experiment. Figure 17 shows the radial distribution of turbulent Reynolds stress, u
x
u
r
¯; due to unavailability of experimental data, this result is not compared with the experiments.Figure 13
Distribution of turbulent intensity of the axial velocity in a swirl jet.Figure 14
Distribution of turbulent intensity of the radial velocity in a swirl jet.Figure 15
Distribution of turbulent intensity of the circumferential velocity in a swirl jet.Figure 16
Radial distribution of turbulent kinetic energy in a swirl jet.Figure 17
Radial distribution of turbulent shear stress in a swirl jet.The figures show that the nonlinear model is well compared with the experiments results. In Figure15, the measured data appeared to show a peak of turbulent intensity of the circumferential velocity at the jet edge (r
/
x~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 k
-
ε 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 ofS
N
=
0.3, 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.Figure 18
Distribution of turbulent intensity of the axial velocity for different swirl numbers (nonlineark
-
ε model).Figure 19
Distribution of turbulent intensity of the radial velocity for different swirl numbers (nonlineark
-
ε model).Figure 20
Distribution of turbulent intensity of the circumferential velocity for different swirl numbers (nonlineark
-
ε model).


## 5. Turbulent Structures in Stuart Vortex

Figure21 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 (S) characterized by negligible spanwise vorticity and the center of vortex (C) 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) Sketch of a coherent vortex. (b) Approximate contour of turbulent stresses at vortex center (Hussain [9]).
(a)
(b)The equation for the stream function of Stuart vortex can be expressed as follows:(51)
ψ
=
log
⁡
cosh
⁡
y
+
A
cos
⁡
x
,
U
x
=
∂
ψ
∂
y
=
sinh
⁡
y
cosh
⁡
y
+
A
cos
⁡
x
,
U
y
=
-
∂
ψ
∂
x
=
A
sin
x
cosh
⁡
y
+
A
cos
⁡
x
.
Here, 0
≥
A
≥
-
1 is constant and indicates the eccentricity of the elliptic streamline of the vortex. In this study, A
=
-
0.5 is used. For this velocity field, the k 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.Figure 22
Streamlines of periodic Stuart vortices.The periodic turbulent structures simulated by nonlinear model are shown in Figure23. It is observed that the turbulent kinetic energy (k), dissipation rate (ε), and turbulent normal stresses in x, y, and z directions (expressed as u
1
u
1
-, u
2
u
2
-, and u
3
u
3
-, resp.) show the contours in elliptical shape at the vortex center. On the other hand, the turbulent shear stress in x
y plane (u
1
u
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.Turbulent structures of Stuart vortex: (a) turbulent kinetic energy,k, (b) turbulent energy dissipation rate, ε, (c) turbulent normal stress in x direction, u
1
u
1
¯, (d) turbulent normal stress in y direction, u
2
u
2
¯, (e) turbulent normal stress in z direction, u
3
u
3
¯, and (f) turbulent shear stress in x
y plane, u
1
u
2
¯.
(a)
(b)
(c)
(d)
(e)
(f)However, in the standardk
-
ε model, coefficient c
μ 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 k
-
ε model causes an inconsistent prediction of turbulent structures at the center of vortex. Figure 24 shows the results simulated by standard k
-
ε model. The model failed to generate the elliptical profiles (focus pattern) for k, ε, 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 c
μ.Turbulent structures by numerical simulation using standardk
-
ε model: (a) k, (b) ε, (c) u
1
u
1
¯, (d) u
2
u
2
¯, (e) u
3
u
3
¯, and (f) u
1
u
2
¯.
(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 fork 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 k 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, x), 50 in transverse (widthwise, y), and 11 in vertical (depthwise, z) directions, for single flood plain case. For double flood plain case, 72 grids are used in transverse direction.Table 2
Hydraulic conditions for numerical simulation of compound channel flow.
Cases
B (cm)
b (cm)
H (cm)
h (cm)
Long. slope
Single flood plain
315
180
17.6
15.0
1/1000
Double flood plain
315
180
17.7
15.0
1/1000Figure 25
Cross section of compound open channel.


### 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.Figure 26
Comparison of averaged cross-sectional velocity profile in a compound open channel; Cal. (1) and Cal. (2) represent the prediction with and without considering horizontal vortices.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)).Comparison of spatial averaged secondary current in a compound open channel.
(a)
Experiment (Knight, 1991)
(b)
Simulation


## 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 fork 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 k 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, x), 50 in transverse (widthwise, y), and 11 in vertical (depthwise, z) directions, for single flood plain case. For double flood plain case, 72 grids are used in transverse direction.Table 2
Hydraulic conditions for numerical simulation of compound channel flow.
Cases
B (cm)
b (cm)
H (cm)
h (cm)
Long. slope
Single flood plain
315
180
17.6
15.0
1/1000
Double flood plain
315
180
17.7
15.0
1/1000Figure 25
Cross section of compound open channel.


## 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.Figure 26
Comparison of averaged cross-sectional velocity profile in a compound open channel; Cal. (1) and Cal. (2) represent the prediction with and without considering horizontal vortices.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)).Comparison of spatial averaged secondary current in a compound open channel.
(a)
Experiment (Knight, 1991)
(b)
Simulation


## 7. Conclusion

In this study, a second-order nonlineark
-
ε 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 (c
μ) and the coefficient of nonlinear quadratic term (c
1, c
2, and c
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 k
-
ε 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 k
-
ε 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 nonlineark
-
ε model can successfully predict the turbulent structures at vortex center. However, the standard k
-
ε 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.



---

*Source: 105809-2015-02-19.xml*