[ Pobierz całość w formacie PDF ]
solving the two partial differential equations for k and µ
(4.7.3) and (4.7.5). The terms on the left side of the equa-
tions and the first term on the right side form a convection-
diffusion equation, which can be solved using standard al-
40
CFD Algorithms for Hydraulic Engineering 4. The Navier-Stokes equations
gorithms (Olsen, 1999). Afterwards, Eq. 4.7.1 is solved to
obtain the eddy-viscosity.
The main work from a numerical point of view is the calcu-
lation of Pk and taking care of the boundary conditions.
Estimating Pk
Pk is calculated by integrating Eq. 3.4.5 over a control vol-
ume. The discretized version for a non-orthogonal 3D grid
is then (Melaaen, 1990):
"U "U "Uiöø 2 "Uj"Ui
-íø - + - = ½T ëø"U-jöø + -------- -
Pk= ½T--------jëø--------j --------øø ---------
íø--------øø
"xi "xi "xj "xi "xi "xj
(4.7.7)
The equation can be rewritten using Eq. 4.2.2:
½T
2
-
Pk = -----[(Éi ) + Éi Éj ] (4.7.8)
j j i
V2
or
½T
2 2 2
-
Pk = -----[2(É1) + 2(É2) + 2(É3) + (4.7.9)
1 2 3
V2
2 2 2
(É1 + É2) + (É1 + É3) + (É2 + É3) ]
2 1 3 1 3 2
Each É has three terms, so 27 terms are obtained altogeth-
er. If the lines in the · direction are identical to the z direc-
tion, A13 and A23 will be zero, so some terms will disappear.
The velocity gradients can be obtained by central differenc-
ing. This means for example:
"U3 1(Un
-
--------- = -- Us ) (4.7.10)
3 3
"¾2 2
where n denotes the cell north of the current cell, and s de-
notes the south cell.
Integration of the term over the cell means that the source
term in the k equation is given by Pk multiplied with the vol-
ume of the cell.
41
CFD Algorithms for Hydraulic Engineering 4. The Navier-Stokes equations
Boundary conditions
Zero gradients are used as boundary conditions close to
the water surface and at outlets. Alternatively, Rodi (1980)
gives a method to account for damping of k from the water
surface.
Boundary conditions are also needed for k and µ at the wall.
Formulas are given by Rodi (1980), based on the assump-
tion that turbulent production is equal to the dissipation of k
near the wall. A formula for µ is then
C3/4k3/2
µ
µ = ------------------ : (4.7.11)
º´n
The formula gives the dissipation in the bed cell directly. It
is therefore specified in the cell, and it would not be neces-
sary to solve the convection-diffusion equation for e for
these cells. However, the solvers usually need a standard
method will all grid cells. To make the solver give the result
of Eq. 3.4.12, the result is multiplied with 1030 and added to
the source, simultaneously as 1030 is added to ap. The
solver will thereby come up with the result of Equation
3.4.12 in the bed cell.
The source terms for k is obtained by integrating the source
of the k-equation (Eq. 4.7.3) over the bed cell:
ÄwUw ÁC3/4kw3/2Uw+
µ
-
Áµ) = ------------- ------------------------------------
(4.7.12)
+"(Pk
´n ´n
´V
The wall law is then used to find Äw and U+w:
ºUw
-
Uw+ = -----------------------
(4.7.13)
lnëø30´nöø
íø------------øø
ks
2
Äw = Á(u+) (4.7.14)
The velocity in the bed cell, Uw, is taken from the result of
the last iteration.
The terms in Eq. 4.7.13 are added to the source. The addi-
tion to ap coefficient becomes:
ÁC3/4kw1/2Uw+V
µ
ap = ap + -----------------------------------------
(4.7.15)
´n
After having used Eqs. 4.7.14 and 4.7.15.
The addition to the source becomes:
42
CFD Algorithms for Hydraulic Engineering 4. The Navier-Stokes equations
ÄwUwV
source = source -----------------
(4.7.16)
´n
Eqs. 4.7.14 and 4.7.15 has also been used.
Inflow boundary
On inflow boundaries, the values of k and µ have to be
specified. A procedure outlined by Olsen (1991) can be
used for free surface flows. It is based on the investigations
of the magnitude (Keefer, 1971) and vertical profile (Naas,
1977) of the eddy-viscosity in rivers. The average eddy-vis-
cosity is given by:
½T = 0.11u*h (4.7.17)
where u* is the shear velocity and h is the water depth. Note
this equation is based on natural rivers. For straight wide
channels with uniform flow, a different coefficient from 0.11
is obtained.
½T
The profile of the eddy-viscosity is given in Fig. 4.7.1. The
value of k at the bed can be estimated from the shear
0.3h
stress:
Ä
Figure 4.7.1 Vertical
k = -------------
(4.7.18)
profile of the eddy-vis- Á cµ
cosity according to
Naas (1977)
The value of µ at the boundary is given by Eq. 4.7.11.
Thereby the bed value of the eddy-viscosity in Fig. 4.7.1
can be computed, using Eq. 4.7.1. Assuming for example
a linear vertical profile of k, with a surface value equal to
half the bed value, Eq. 4.7.1 can be used to compute the
vertical profile of µ.
4.8 The stress terms
The Boussinesq approximation for the turbulent stress
term is given as:
1------ "U- + "Ujöø
"-
- -
uiuj = -- ½Tëø--------i --------øø
(4.8.1)
Á"xj íø "xj "xi
The first term on the right side is the ordinary diffusive term,
used for solving the convection-diffusion equation. The
second term is here denoted the stress term. It is often ne-
glected in CFD programs tailor-made for hydraulic engi-
43
CFD Algorithms for Hydraulic Engineering 4. The Navier-Stokes equations
neering, such as SSIIM and Telemac. The reasons are
given in the following.
In hydraulic engineering flow fields, there is usually one
dominant flow direction. Let us assume such a situation,
with the x-direction being in the flow direction in an orthog-
onal grid. The z direction is vertical and y is in the cross-
streamwise direction. The stress term for the velocity U in
the x direction can then be written:
1---- "Uj 1 " 1----- 1----- -
" " "
-- ½T--------- = -- ½T"U + -- ½T"V + -- ½T"W
- ------ ------- - ------ - -------
Á"j "x Á"x "x Á"x "x Á"x "x
(4.8.2)
[ Pobierz całość w formacie PDF ]