京 都 大 学 防 災 研 究 所 年 報 第48号 B 平成17年4月
Annuals of Disas. Prev. Res. Inst., Kyoto Univ., No. 48 B, April 2005, pp. 691-707

A RANS Solver Using a 3D Unstructured FVM Procedure

Hao ZHANG*, Hajime NAKAGAWA, Taisuke ISHIGAKI** and Yasunori MUTO*

Graduate School of Engineering, Kyoto University
Dept. of Civil & Environmental Engineering, Kansai University**

SynopsisA three-dimensional unstructured mesh based Reynolds-Averaged Navier-Stokes equation
solver is presented with a Finite Volume Method. A strict but effective data structure is introduced
to store the control volume information resulted from different meshing strategies. The power law
scheme is adopted for the spatial discretization, and the implicit Crank-Nicolson scheme isemployed for the temporal integral. In order to avoid the checkerboard phenomenon, the surface
fluxes are calculated from the Rhie-Chow interpolation. The final algebraic equation system issolved using a preconditioned GMRES solver incorporated with an ILUTP preconditioner. The
proposed solver has been applied to three kinds of turbulent flows in laboratorial flumes. It is foundthat the solver is able to reproduce the selected flows in different domains with a reasonableaccuracy and manifests a potential to apply to the actual river conditions.
Keywords: Three-dimensional, FVM, unstructured mesh, k- ε model

1.

Introduction

The two most commonly adopted meshing strategies
in CFD (Computational fluid dynamics) community are
generally referred to as structured mesh and unstructured
mesh. The signature feature of structured mesh is theimplicit grid structure, which may alleviate the need tostore the mesh connectivity and permit the constructionof rapid solution algorithms (Mavriplis, 1996). In orderto take full advantages of structured mesh, many
researchers have been devoting themselves to improving
existing algorithms and developing new schemes. (e.g.
Jongen, 1997; Kimura and Hosoda, 2003; Li and
Fleming, 2003; Johnston and Liu, 2004; Brüger et al.,2005) However, despite its widespread use, thestructured mesh methods are generally limited to solution
domains of a relatively simple shape. It is very difficultand sometimes impossible to generate structured meshfor complex geometries, especially in 3D space. On theother hand, flows in complex geometries are usuallyencountered and of great interest in the engineering
practices nowadays. For example, turbulent flows pass
the natural rivers with extremely irregular cross sections
or artificial channels with small but important hydraulicstructures. An accurate resolution of the study domain
and local mesh refinement necessitates the employmentof unstructured mesh methods.
Although the past several decades have witnessed thesignificant development in unstructured mesh methods,
most of the published literatures dealt with flow
phenomena in other fields such as aerodynamics. (e.g.
Mavriplis, 1997; Nakahashi et al., 1999; Hassan et al.,
2000; Nakahashi et al., 2003; Basara, 2004) Hydraulicengineering poses many special problems requiring
solution methods probably not included in other fields.
But the related research is quite few to the authors’knowledge and is generally confined to 2D flows. For
instance, Kim et al. (1997) reported a RANS
(Reynolds-averaged Navier-Stokes equation) solver
based on a 2D unstructured mesh, and the solver wasverified to be capable of predicting some wall-bounded
turbulent flows with a good accuracy. Nallapati andPerot (2000) developed a staggered unstructured meshmethod for the 2D calculation of free surface flows and it
was applied successfully to some benchmark problems.
Recently, Olsen (2003) constructed a 3D model tocalculate the formation of a meandering channel. The
turbulent flow was predicted with an unstructured meshbased k- ε model. Unfortunately, very few information
has been shown on the quantitative comparison between
the computational result and the physical model
measurement in the paper. Haque et al. (2005) also
carried out a 3D numerical simulation with a RANS
solver provided by the commercial software FLUENT.
The hybrid mesh adopted in the study was argued to beable to resolve the detailed hydraulic structures with arelatively small total number of cells. They compared thetemperature distribution within a dam and its forebay
situated on the Colombia River and found a reasonablesimilarity between the calculation and measurement. But
the velocity field was not validated due to the shortage offield data.
In this study, an unstructured mesh based RANS
solver is proposed based on a 3D FVM (Finite volume
method) procedure. The standard k- ε model and some
non-linear k- ε model have been integrated with thissolver. The wall function approach is employed toresolve the near-wall area. Hybrid polyhedral mesh up tosix faces may be used. This solver is designed as amodule in a morphological model and may serve as a
powerful tool for the investigation of river engineering.
The paper is organized as follows. In the succeeding
section, a set of governing equations is firstly introduced.
After that comes the detailed information on FVM
formulation in Section 3. From Section 4 to Section 6,
the discretization methods, the boundary conditions andthe solution algorithms for sparse equation systems arepresented one after another. In Section 7, three numericalexamples are tested and computational results arecompared with existing laboratorial measurements,
which is followed by the final conclusion in Section 8.
2.

Governing equations

An engineer is usually concerned with thetime-averaged properties of turbulences (Launder and
Spalding, 1972). This suggests ways of accounting forthe turbulences on the mean flow behavior by solving the

RANS rather than solving the NS (Navier-Stokes
equation) directly. In hydraulic engineering, as is
well-known, the k- ε model and its variants are thewidely used turbulence models based on the RANS. In
this kind of models, the mean quantities are directly
solved from the corresponding transport equations andthe time-averaged turbulences are modeled byintroducing some assumptions and simplifications.
2.1 Mean flow
The unsteady 3D RANS and continuity equation
expressed in a Cartesian coordinate system with thetensor notation are as follows.

∂u
∂2ui
∂ui
1 ∂p
1 ∂τ ij+ u j i = fi −
+v
+
∂t
∂x j
ρ ∂xi ∂x j ∂x j ρ ∂x j

(1)

∂ui
=0
∂xi

(2)

where ui= time-averaged velocity; xi= Cartesian
coordinate component; ρ = density of the fluid; fi=body force; p= time-averaged pressure; v= molecular
kinematic viscosity of the fluid; τ ij = − ρ ui' u 'j , are theReynolds stress tensors, and ui' is the fluctuating velocity
component. The equation system is not closed, and the
Reynolds stress tensors have to be evaluated in someother ways.
2.2 Turbulence closure
In the standard k- ε model, the Reynolds tensors are
acquired through the linear constitutive equation.

2−ui' u 'j = 2vt Sij − kδ ij
3(3)

where k= turbulence kinetic energy; δ ij = theKronecker delta; vt = eddy viscosity and Sij= thestrain-rate tensor, the latter three are expressed by

i= j
i≠ j

1 if
0 if

δ ij = ν t = Cµ

k2

ε

(4)

1  ∂u ∂u j 
Sij =  i +2  ∂x j ∂xi in which Cµ is a coefficient, and is usually set to be aconstant and equal to 0.09, ε is the dissipation rate of
the turbulence kinetic energy k. Two transport equations
as described below are employed to estimate k and ε ,respectively.
vt  ∂k 
∂k
∂k
∂ 
+uj
=
 v + 
 +G −ε∂t∂x j ∂x j  σ k  ∂x j 

(5)vt  ∂ε 
∂ε
∂ε
∂ 
ε
+uj
=
 + (C1ε G − C2ε ε ) v + ∂t
∂x j ∂x j  σ ε  ∂x j 
k

(6)

∂ui
∂x j

(7)

and the model constants suggested by Rodi (1980)
generally take the universal values as below.

σk =1.0 σε =1.3 C1ε =1.44 C2ε =1.92(8)The standard k- ε model suffers from some inherent
defects in particular the omission of any anisotropic eddy
viscosity effects. This may be partially cured byintroducing a non-linear constitutive relation between the
turbulence stresses and the mean strain rate satisfying
certain tensorial properties. Efforts and achievements
have been make by many research groups such as
Rubinstein and Barton, 1990; Gatski and Speziale, 1993;
Shih et al, 1995 and Kimura and Hosoda, 2003. A
general form for a quadratic constitutive equation can besummarized as2
ui'u'j = kδij − 2vt Sij
33 
k ∂u ∂u ∂u j ∂ul 2 ∂uk ∂um δ +a1 2  i l +−ε  ∂xl ∂x j ∂xl ∂xi 3 ∂xm ∂xk ij 
k 3  ∂u ∂u j 1 ∂um ∂um δ +a2 2  i−ε  ∂xl ∂xl 3 ∂xk ∂xk ij 
+a3

(9)

k 3  ∂ul ∂ul 1 ∂um ∂um δ −
ε 2  ∂xi ∂x j 3 ∂xk ∂xk ij where a1, a2 and a3 are coefficients, which may havedifferent values due to different models.
3.

∂
φ dV+∫ φu⋅ ndS=∫ Γ∇φ⋅ ndS+∫bdV
S
S
V
∂t ∫V

(10)

where V= the volume of the CV; S= the CV surface witha unit normal vector n directing outwards; φ = general

where G= the rate-of-production of the turbulencekinetic energy k, is defined as

G = −ui' u 'j

into a number of continuous polyhedral CVs (Control
volumes). If one integrates the governing PDEs (Partial
differential equations) over a CV, the general form reads

Finite volume formulation
In an FVM procedure, the study domain is divided

conserved quantity representing either scalars or vector
and tensor field components; u= the fluid velocity vector
whose Cartesian components are ui or (u, v, w); Γ =diffusion coefficient and b= the volumetric source of thequantity φ . It is noted that Eq. 10 has a left side of atransient term and a convective term balanced by adiffusive term and a source term on the right. The
equation is mesh-independent and is valid for arbitrary
polyhedral CVs. As a result, an FVM procedure is able
to take full advantages of an unstructured mesh system.
Moreover, summation of the equations for all CVs leads
to the global conservation equation because the inner CV
faces will cancel out. This is also one of the most
attractive points to use FVM in engineering practices.
3.1 Data structure
In an unstructured mesh method, the meshconnectivity is not implicitly known. Since the CVs may
be numbered in any order and have any number ofneighbors, an extra space in the memory should bemaintained for that kind of information. This may
sometimes result in a considerable consumption ofcomputer resources.
As is known, any kind of polyhedral mesh may be
used in the simulation theoretically, but the generally
adopted unstructured mesh in engineering practices isconfined to tetrahedra, pyramids, prisms andhexahedra. A partial explanation lies in the difficulty ofgenerating polyhedral mesh with more than six faces.
Amongst all the above mesh types, the hexahedron
may be accounted as the general case, and thequadrilateral is the general case for the CV facecorrespondingly. It indicates that a hybrid mesh system
in practice is possible to be considered as a hexahedralmesh provided that the polyhedra with less than six
faces are assigned some nominal ones to satisfy theconvention. This provides a way to simplify the data
structure in storage.
In order to define the connectivity of the CVs, there
are a lot of alternatives. The most general data structure

5

5

8

6

7

6

5,6 ,7,8

5,6 ,7,8

7,8
1

1

1

4

2

4

1
2

3

Hexahedron

3 ,4

2

Prism

3

2

3 ,4

Pyramid

Tetrahedron

Fig.1 CV and its vertices
contains all the information of the CV explicitly. That
is: a CV is defined by its six faces, faces by the lists ofedges and edges by their corresponding vertices.
However this method is not preferred to reduce thenumber of arrays needed for the definition of meshconnectivity. A simplified data structure is introduced
here. This method is also employed in some CFD
codes and recommended by Ferziger and Perić (2002).
As indicated in Fig. 1, the CV is defined by a list of
its eight vertices in a counter-clockwise order. With
such definition, the faces enclosing the CV and theedges forming the faces are also uniquely identified.
They may be implicitly organized in an ordered way
without occupying any computer memory. Moreover,since the arbitrary polyhedral mesh is stored as ahexahedral mesh, the number of the neighboring CVs
is fixed to be six. Hence the neighboring CV whichshares a common face with the current CV may also beindexed in the same ordered way as that for the faces.
For nominal hexahedra, due to the existence of nominal
faces, some of the vertices will be repeated duringstoring as also shown in Fig. 1. Using this kind oftreatment, the input of the mesh system includes only alist of nodal coordinates, a list of vertices of the CVs
and a list of the neighboring CVs.

3.2 Surface and volume calculation
With the aforementioned data structure, it is no needto evaluate the geometrical elements for arbitrary
polyhedra. The calculation is limited to quadrilaterals in
2D and hexahedra in 3D. As the area and the center of atriangle is very easy to acquire. One can subdivide aquadrilateral face into two triangles and the vector area
is approximated by the summation of those of the two
sub-triangles. For instance, the vector area of thequadrilateral (1,2,3,4) in Fig. 2 is computed fromri = xi i + yi j+ zi k
S1234 = S123 + S134 =1
( r2 − r1 ) × ( r3 − r1 ) + ( r3 − r1 ) × ( r4 − r1 ) 
(12)2
where ri = the ith vertex of the CV; i, j, k = unit vectors
in the x, y, z coordinate directions; S123, S134= vector
areas of the sub-triangles ∆123 and ∆134, respectively andS1234= vector area of the quadrilateral surface (1,2,3,4).
The center of the surface is the mean value of thoseof the sub-triangles weighted by the correspondingarea, i.e.
c1234 =

S123c123 + S134 c134
S1234

8

5

(11)

(13)

1

7

6
z

C123
1

y

4
x

2

2

y

C134

x
3

3

Fig.2 Surface and volume calculation4

where c= the center of a surface and S= the area of asurface.
Extending the above logic from 2D to 3D, ahexahedron is divided into 6 tetrahedra. The evaluationof the volume is thus acquired by summation of thevolumes of all the sub-tetrahedra corresponding to their
partitioning. And the volume of a tetrahedron is a dot
product of two vectors, for example the volume oftetrahedron (1,5,6,7) is1
V1567 = (r5 − r1 ) ⋅ S567
3(14)

where V1567= the volume of tetrahedron (1,5,6,7). The
centroid of the CV is the mean value of those of the
sub-tetrahedra weighted by the corresponding volume.

cV =

∑V c
i

i

V

Discretization methods∂φ∑φ (u S ) = ∑Γ ∂n S + b − s φ
f

f⊥

f

f

φ f = α f φP + (1 − α f )φ A
where α f =

f

f

P

P P

(16)

f

where uf⊥= the fluid velocity normal to the surface; bP=part of the source term containing all the contributions
excluding unknown variables and − s pφ p = part of thesource term including the unknown variables which canbe treated implicitly. The subscript P= the present CV
and the subscript f= the face of the CV.
It is readily seen that Eq. 16 is not an explicit
expression of the variables defined at the centers of theCVs. Values at other locations have to be obtained bysome kind of interpolation methods. The diffusive term
contains the gradient of a quantity, which necessitatessome numerical differentiation techniques.

(17)

d Af
d Af + d fP

(18)

in which the subscript A= the adjacent CV, dfP and dAf
are the distances from the surface to the present CV andto the adjacent CV, respectively. If the line connecting
the two CV centers does not pass through the center ofthe common face center, a correction term may beintroduced. In this case

φ f = φ f ′ + (∇φ ) f ′ ⋅ (rf − rf ′ )(19)where f’= the intersection of the surface and the line
connecting the two neighboring CVs (see Fig. 3), thegradient at f’ is obtained by interpolating the cell-center
gradients at either side of the face. And the Gauss’
theorem may be employed to determine the cell-center
gradients.

A

In a collocated FVM procedure, the variables aredefined at the center of the CV. A second order midpoint
rule is generally used for the integral approximation. If
the transient term is absent for the time being, the control
volume equation (i.e. Eq. 10) may be discretized term byterm and written asf

the surface,(15)

where cV= the centroid of the CV; Vi= the volume of
sub-tetrahedron i (i ranges from 1 to 6) and ci= thecentroid of sub-tetrahedron i .
According the above analysis, the sub-routine for the
calculation of the geometry elements can be easily
coded.
4.

An arithmetic interpolation method is commonlyused to evaluate the surface values. For a quantity φ on

f’
f

P

Fig.3 Interpolation of the surface value
The discretization methods adopted here have beensuggested by some other publications such as Ferziger
and Perić (2002) and applied to various problems. It
serves as a basic scheme whenever a surface quantity
needs to be interpolated.
4.1 Spatial discretization
From the viewpoint of stability and ease of
programming, higher order scheme is avoided in thisstudy. The power law scheme is employed during the
spatial discretization. This scheme is relatively easy andhas been confirmed to be applicable in 3D calculations(Wilson et al., 2003; Olsen, 2004).
With the power law scheme, Eq. 16 can be finallyassembled as{} (20)( )∑f Df A Pf + max( −Ff ,0)(φP −φA ) + Ff φP= bP − sPφP

where the strength of the convection Ff, diffusion
conductance Df and the ratio of them are given as
Ff = u f ⊥ S f , D f =

and

( )

Γf S fd AP, Pf =[(A Pf = max 0, 1 − 0.1 Pf

Ff
Df

)]
5

(21)

ΓA ΓP
α f ΓP + (1 − α f )ΓA

(22)

(23)

A simple arithmetic mean for the surface flux maylead to checkerboard variable distribution, which has
caused the slow acceptance of the use of collocatedmesh. This problem can be cured by the interpolation
method proposed by Rhie and Chow (1983). The
method introduces an additional term related to thepressure gradient when calculating the fluxes on thesurface.
As is known, the unknown quantities of the present
CV can be finally expressed by all of its neighboring
CVs after discretization. For instance, the momentum
equations for u at present CV and one of its adjacentCVs are written as∂pdV P + bPV ∂x

a P u P = ∑ a nb u nb P − ∫
nb

V

nb

∂p
dV f + b f
∂x

(25)

Approximating the solution uf of Eq.25, theinformation from Eq.24 can be used. By using somelinear interpolation and simplification, the following
equation is obtained.

uf =uf +∂p1  ∂p∫dV f − ∫dV f 
VV∂xa f  ∂x


(26)

where

The interpolation of the diffusive coefficient on thesurface and the surface flux deserve special attention.
For the diffusive coefficient, the harmonic mean may
reflect more physics and reasonableness in particularnear the boundary. This results inΓf =a f u f = ∑ a nb u nb f − ∫u f = α f u P + (1 − α f )u A

a f = α f aP + (1 − α f )a A∂p

∫ ∂x dV∂p∫ ∂x dV

where a= coefficient for the unknown at the center of
the approximated CV and nb= the neighboring CV.
From the conservation principle of the FVM
formulation, the velocity at the common face of the twoneighboring CVs must also have a discretized
momentum equation of the similar form as that of
Eq.24, i.e.

fV

=α f ∫V∂p∂p
dV P + (1 − α f )∫dV A (27)V∂x∂x= S fx ( p A − pP )in which Sfx= projected area of the surface to the yz
plane (perpendicular to the x axis). The extension toother velocity components is straightforward.
4.2 Temporal integral
At the end of the spatial discretization, one can get∂∫ φ dV = F
∂t V

(28)

where F = ∑ anbφnb + bP − aPφP
nb

Taking the second order implicit Crank-Nicolson
scheme as follows

φ m+1 − φ m

F m +1 + F m
(29)
∆t
2
where the superscript m and m+1 stand for the previous
and the current time step. One can arrive at the final
algebraic equation set.
V=

a P φ Pm +1 = ∑ a nb φ m +1 + b P

(24)

∂p
dV A + bA
a Au A = ∑ a nb u nb A − ∫
V ∂x
nb

f

V

(30)

nb

in whichm +1
aPm +1 = ∑ anb+ ∑ F f + sPm +1
nb

aP = aPm +1 +bP = bPm +1 +

f

2V
∆t
2V m 

φP +  ∑ anbφnb + bP − aPφP 
∆t
 nb


(31)
m

It is seen from Eq. 31 that after the time integral, atime-related term is introduced to the coefficient of aP
comparing with the steady case. And besides atemporal term, there is another contribution from theprevious time step to the source term.

then introduced to ensure the continuity equation. An
equation set for the pressure correction will be obtained.

4.3 Pressure-velocity coupling
The solution procedure follows the SIMPLE(Semi-implicit method for pressure-linked equations)algorithm. The main concept of this method is to guess
the pressure and get a pressure correction with thecontinuity equation. The procedure is summarized here.
In the derivation, the guessed pressure and thevelocity field not satisfying the continuity equation are
denoted with an index *, the correction of the variable isdenoted with an index ‘, and the variable without any
superscript stands for the corrected value. Hence,

This equation set is closed and solution of theequations yields the pressure correction. The velocity
field is renewed and a new iteration starts untilconvergence.

ui = ui* + ui′
p = p* + p ′(32)

Extracting the pressure term from the source term in
the discretized momentum equation, the velocity
component in the x direction satisfies*
aP u*P = ∑ anb unb− ∑ A fx p*f + b p
nb

f

(33)

On the other hand, the same momentum equation
based on the corrected velocity component has thefollowing formaP u P = ∑ anb unb − ∑ A fx p f + b p
nb

f

(34)

If Eq. 33 is subtracted from Eq. 34 and Eq. 32 is taken
into account, the velocity correction is obtained.

u ′P =∑ anb u ′nb ∑ A fx Pf′
−
aP
aP

(35)

In a SIMPLE method, the first term on the right handside is omitted. Then the velocity correction in the CV
center has a simple relation with the pressure correction.
For the velocity correction on the surface, the same
relationship is assumed to be valid instead ofinterpolating the value from the CV centers. It meansu ′f = −1
A fx ( PA′ − PP′ )
aP

(36)

The correction of other velocity components may bederived in the same way. The corrected velocity field is
A2fi  
 ( PP′ − PA′ )  = 0∑  ∑ u fi Afi  +  ∑
f  i = x, y, z  i = x, y, z aPi 



5.

(37)

Boundary conditions and near-wall treatment

Appropriate boundary conditions should be specified
depending on the nature of the flow. The inlet boundaryis generally considered as a Dirichlet boundary and all
the quantities have to be prescribed. At the outlet
boundary, the flow information is usually little known. A
Neumann boundary with zero gradients may be assumed
if the outlet is set as far downstream of the study domain
as possible. The water surface is considered as asymmetrical plane for the time being. It may be
acceptable for many hydraulic problems. Special
attention has been paid for the impermeable wallboundaries.
The no-slip condition is the appropriate condition for
velocity components at both the riverbed and the sidewalls. However the wall function approach is preferredhere to avoid the possible integration through the viscous
sub-layer and implement the wall roughness more
flexibly. In the wall function approach, the near wall CV
velocity is assumed to be parallel to the wall and denotedby u//. Although it is not always the case, the treatment
can be simplified without significant influence on the
result.
With the definition of the dimensionless distance y+and dimensionless velocity u+ as follows

y+ =u* y ⊥

ν, u+ =

u //
u*

(38)

where u*= the friction velocity near the bed and y⊥=
the normal distance from the center of the near wall CV
to the wall surface, the universal wall function can be
expressed by

u+ =

1

κ

(

ln Ey +)

(39)

where κ = the van Karman constant (= 0.41) and E=
roughness parameter of the wall. Assuming that theflow is in local equilibrium, i.e. the production anddissipation rate of the turbulence are nearly equal, one
can obtainu* = Cµ1/ 4 k P1/ 2(40)

Then the wall shear stress is written asτ w = ρu*2 =ρCµ1 / 4 k P1/ 2u //
u+

(41)

In the momentum equations, the link with the wallis suppressed by setting it to zero and adding the wall
force in Eq.41 as a source term. The normal derivativeof k at the wall boundary CV is set to be zero in thek- ε equation, and the production in the wall region iscomputed fromGP =

τ w ∂u // τ w u //=ρ ∂n ρ y ⊥(42)

ε in the near wall CV is directly set toεP =C µ3 / 4 k P3 / 2
u*3
=
κy ⊥
κy ⊥

(43)

The roughness parameter in Eq. 39 provides a simple
way to take the bed roughness into account. Suggested
by some researchers (e.g. Wu et al., 2000; Salaheldin,2004), this parameter may be evaluated as below.

E = exp κ ( B − ∆B ) (44)

where B=5.2, is a constant; ∆B= roughness function
defining the shift of the intercept due to roughness effect,
it is evaluated from0ks+ < 2.25
∆B = Bm sin 0.4285 ln ks+ − 0.811  2.25 ≤ ks+ < 90 (45)




Bm
ks+ ≥ 90()

in which
Bm = B − 8.5 +k s+ =

Solution methods

7.

Applications

The proposed methodology is applied to three kindsof turbulent flows in the following sub-sections.
Different turbulence models and different meshingstrategies have been tested.
7.1 Flow in a straight rectangular channel
Imamoto et al. (1987) carried out a series ofexperiments concerning the flow in open channels. An
LDV (Laser Doppler velocimeter) system was used tomeasure the fully developed turbulent flow velocities in atransverse cross-section. In this section, one case wasselected to verify the proposed methods.

ln ks+

u* k sv

κTable 1 Experiment conditions(46)

and ks is the equivalent roughness height. For a hydraulicsmooth bed, ks=0, and for a rough bed, many empirical
relations are available in the literatures according to the
bed conditions (e.g. Chang, 1988; van Rijn, 1993).
6.

The equations are solved in an implicit decoupledway. The final algebraic equation systems areunder-relaxed before submitted to the equation solver.
The widely used method proposed by Patankar (1980)
is adopted here. It has been found to be very efficient
owing to its increasing of the diagonal dominance ofthe coefficient matrices.
Considering the sparseness and non-symmetry
characteristic of the coefficient matrices, Krylov
subspace iterative methods are preferred. The mostpopular and representative methods may include
GMRES (Generalized minimal residual method),
Bi-CGSTAB (Bi-conjugate gradient stabilized method)
and TFQMR (Transpose-free variant of quasi-minimal
residual method). In order to get practically useful, theiterative solver usually works with a suitablepreconditioner. In this study, a preconditioned GMRES
is employed together with an ILUTP (Incomplete LU
factorizationwiththresholdandpivoting)
preconditioner. A complete presentation on iterative
methods and preconditioning techniques can be foundin the book published by Saad (2003).

Width Depth Discharge Slope Reynolds Froude
B(cm) H (cm) Q (l/s)
Ienumber number
20.0

4.00

2.055

1/1,400

7,700

0.48

The experiment conditions are given in Table 1. The
laboratory flume is geometrically symmetrical, and only
half of the domain is calculated by employing asymmetrical plane boundary condition. A longitudinal

z/H
1.0

z/H

0.5

0.5

0

1.0

1.0
(a) Experimental result

2.0

0

y/H

z/H

z/H

1.0

1.0

0.5

0.5

0

1.0
2.0
(b) Linear model result

y/H

0

1.0
2.0
(c) Non-linear model result

y/H

1.0
2.0
(d) Structured mesh resulty/HFig. 4 Comparison of the longitudinal velocity u (cm/s)
distance of 140cm (=35H) is assumed to be long enoughto diminish the influence of inlet and outlet boundaries
and has been chosen for the simulation after some trial
computations. The standard k- ε model and a non-linear
k- ε model proposed by Kimura and Hosoda (2003) are
implemented on a hexahedral mesh with a total numberof 10,626, respectively. A logarithmic velocity profile isprescribed at the inlet boundary. The turbulent quantities
k and ε are specified corresponding to a viscosity ratioof 10.0 and taking the turbulence intensity 8%.
The computed longitudinal velocity and turbulencekinetic energy in the transverse cross-section comparedwith the experimental result are shown in Fig.4 and Fig.5,
respectively. Ishigaki (1993) simulated the same flow

z/H
1.0

z/H
1.0

0.5

0.5

0

1.0
2.0
(a) Experimental result

y/H 0

z/H
1.0

z/H

0.5

0.5

0

1.0
2.0
(c) Non-linear model result

y/H

1.0
2.0
(d) Structured mesh result

y/H

1.0

1.0
2.0
(b) Linear model resulty/H0Fig. 5 Comparison of turbulence kinetic energy normalized by squared friction velocity k/u*2

phenomenon with two ARS (Algebraic Reynolds stresssimulation) models based on a structured mesh. The
result of one of the models proposed by Naot and Rodi
(1982) is also shown in the figures.
A standard model does not include any informationof the vortices and the secondary currents, which inreturn affects the prediction of the mean flow pattern.
This may be evidently observed in the comparison of thelongitudinal velocity in Fig. 4. Both the non-linear modeland the ARS model give a closer result to the LDV
measurements than the linear one. It is also interesting tonote that the non-linear model in this selected case seems
to perform even better than the ARS model. The reason
is not much known, but it demonstrates that the−u ′v ′ / u*2

unstructured mesh method in this study is at least
comparable to a structured one.
In the comparison of the turbulence kinetic energy(the absolute value has been normalized by the squared
friction velocity), the non-linear model again manifests
itself to be the best solution. However, in allcomputational cases, the area near the wall boundaries
has been slightly over-estimated. A further comparison iscarried out in Fig. 6 between the experimental result andthe non-linear model for the lateral Reynolds stress
distribution at different water depths. The values areagain normalized by the squared friction velocity. The
over-estimation of the Reynolds stresses near theboundary area is observed. And away from the boundary,−v ′w′ / u*2− w′u ′ / u*2

Fig. 6 Comparison of the lateral Reynolds stress distribution (Left: Experiment; Right: Computation)the model result is very close to the measurements. The
employment of the wall function approach may be themain error source for this simulation.
7.2 Flow in a square embayment
The flow in an embayment is one of the mostimportant hydraulic phenomena in the river engineering
practice. Due to the complex flow structure, the balanceof the sediment transport is broken. It results in a great
number of problems, for instance, the local scouring, thechange of environmental parameters, etc. By using two
LDAs (Laser Doppler anemometers) and an EMC(Electromagnetic current meter), Muto et al. (2000)
experimentally investigated the flow exchanges between
the main channel and an embayment area. In this section,
the methodology is applied to one of this kind of flows
and the result is compared with the measurements.

profile in the vertical direction. At the beginning of thecalculation, this velocity profile is applied to the wholedomain.
Two kinds of meshing strategies have been tested
with a non-linear k- ε model under the same initial andboundary conditions. The mesh system near theembayment in either case is shown in Fig. 8 and Fig. 9,
respectively.
Around the embayment, the mesh is clustered for abetter resolution, and away from the embayment, themesh is relatively coarse to save computational time. The
hybrid mesh consists of both hexahedra and prisms
having a total mesh number of 13,086. And thehexahedral mesh is a bit finer with a total number of
19,683.

Flow

B
L
x
y

b

Fig.7 Top view of the experimental setup
Fig.8 Hybrid mesh
The laboratory flume consists of a main channel
which has a width B=16.0cm and a flood plain with awidth b=16.0cm. Part of the flood plain is removed
where forms an embayment with a length L=16.0cm.
The experiment set up is shown in Fig. 7. And theexperiment conditions are summarized in Table 2.
Table 2 Experiment conditions
Discharge
Q (l/s)

Water
depth
H (cm)

Friction
velocityu*(cm/s)

2.271

3.8

2.055

Reynolds Froudenumber number
9,650

0.74

In the computation, the inlet boundary is set at the
distance of about 10b from the upstream of theembayment and the outlet is 20b from the downstream ofthe embayment. This distance seems long enough andmay diminish the effect of the inlet and outlet boundaries.
The inlet flow is assumed to have a logarithmic velocity

Fig.9 Hexahedral mesh
The longitudinal velocity and the horizontal vortex at
half water depth are shown in Fig. 10 and Fig. 11,
respectively.
Both the hybrid mesh and the hexahedral mesh work
well in this test case. The flow pattern in the embayment
has been reproduced with a reasonable accuracy. A large

(a) Experiment result(b) Hexahedral mesh result
Fig. 10 Comparison of the longitudinal velocity profile(c) Hybrid mesh result

10cm /s(a) Experiment result(b) Hexahedral mesh result(c) Hybrid mesh result
Fig. 11 Comparison of the horizontal cross-sectional velocity profile (u, v)
circulating flow is induced in the embayment area asshown in Fig. 11. This circulation occupies almost thewhole region of the embayment with a maximum
velocity about a quarter of that in the main channel. At
the center of the circulation, the flow is stagnant.
Velocity gradients in the junction zone (i.e. the interface
between the main channel and the embayment) are quite
steep. Although there are some differences in the result
due to the two different mesh systems, it is very difficult
to distinguish the winner.

(a) Plan viewA
14.6cm
Flow

60cm

1.5cm

99cm

A

(b) Side view: cross-section A-A

5.6cm

7.3 Flow around spur dykes with local scour holes
The final application is the prediction of the flow
field around a series of impermeable spur dykes with

15cm

Fig. 12 Experiment setup

Table 3 Experiment conditions
Discharge Mean velocity
Q (l/s)u (cm)
10.52

165
150
135
120
y (cm)
105
90
75
60
45
30
15
0
100

18.98

Slope
Ie

Approach flow depth
H (cm)

Shear velocityu*(cm/s)

1/3,000

5.60

1.35

-13 -11 -9 -7 -5 -3 -1 1

150

200

250

300

350

400

450

3

Shear velocity Reynolds Froude
ratio (u*/ u*c) number number
0.915

10,682

Unit: cm

500

550

600

650

Fig. 13 Bed topography under equilibrium condition
local scour holes. Detailed measurements were reported
by Khaleduzzaman (2004) and Zhang et al. (2005). The
experiment is conducted in a straight tilting flume, whichis 20m long, 0.99m wide and 0.3m deep. The initialriverbed was covered with 15cm-thick fine sediment.
The sediment has a mean diameter of 196 µ m and ageometric standard deviation of 1.496.
The experiment setup and the experiment conditions
are given in Fig. 12 and Table 3, respectively. It satisfies
the condition of clear-water scouring. The final riverbed
topography after a continuous running of 893 hours isshown in Fig. 13. Severe local scour holes are observedat the toes of all the spur dykes. Around the first pair, all
the sediment has been eroded. The main channel area hasy (cm)

0.26

700 750x (cm)

been significantly degraded, while along the channelbank, obvious deposition occurs. The 3D flow velocity inthe equilibrium state is measured with an I-shape and an
L-shape electromagnetic velocity meters under dynamic
flow conditions. The average of 600 samples processedby a computerized data acquisition system at a frequencyof 10Hz is taken as the time-averaged value of eachmeasured quantity at each point.
The flow calculation is based on the final riverbed
with a k- ε turbulence model. A hexahedral mesh systemwith a total number of 31,808 is employed for the
calculation. Considering the geometrical symmetry, only
half of the flume is selected as the computational
domain.

10 cm/sy (cm)x (cm)x (cm)
Fig. 14 Stream-wise velocity profile (u, v) around the first five spur dykes at z=1.0cm(Experiment: top; Computation: bottom)

y (cm)
10 cm/sx (cm)y (cm)

Fig. 15 Stream-wise velocity profile (u, v) around the first five spur dykes at z=3.25cm(Experiment: top; Computation: bottom)
The stream-wise velocity profiles (u, v) around thefirst several spur dykes at two different layers (z=1.0cm
and z=3.25cm) are shown in Fig. 14 and Fig. 15,
respectively.
Although the magnitude of the velocity vector in the
two layers is different, the flow pattern shows somesimilarity. The velocity diverts at the heads of the spurdykes. In the main channel area, the velocity isintensified, and there is a significant velocity reductionin the embayment areas formed by the consecutive spurdykes. Take a look at the experimental data, one mayfind that there is a wake vortex in each embayment area.
The center of the vortex is closer to the upstream spurdyke forming the embayment area. This phenomenon
has been reasonably reproduced by the computation.
A further inspection may lead to some moreinteresting discoveries. In the computational result, one
may find a small eddy at the downstream corner of each
embayment. These eddies are not found in the plot for
the experimental data. This may be due to the shortageof enough measured data. If the data has been collected
at more locations, the small eddies may be expected toappear. In fact, similar small eddies have been observedon the free surface during the surface flow visualization(Zhang et al., 2005). In the first embayment, thecomputed vortex system is a little different from theothers, which coincides with the experimentmeasurements. The obvious 3D characteristics of theflow structure due to the local scour hole in this area
may be responsible for this difference. As has
mentioned before, the first spur dyke is surrounded byx (cm)

the deepest scour hole. It has a great effect on the flow
field.
In the junction zone, the measured velocity seems a
bit larger than the computed one. This may be attributedz (cm)
10cm/sz (cm)y (cm)y (cm)
Fig. 16 Transverse velocity profile (v, w) at the head of thefirst spur dyke (Experiment: top; Computation: bottom)

to the coarse mesh adopted in the test case. As thevelocity gradient in this area is quite deep, a finer meshmay be suggested for a better resolution.
A comparison of the transverse cross-sectional
velocity profile (v, w) at the head of the first spur dyke
(i.e. x=150.0cm) is depicted in Fig. 16 finally. This figure
gives an intuitive image of the vertical flow pattern in thelocal scour hole. The basic aspects of the flow structure
have been reasonably captured with the numericalsimulation. The vertical velocity component w at thehead of the spur dyke is quite large. A big anti-clockwise
vortex occupies almost the whole scour area if one takes
a look from the downstream. No wonder the sediment
just in front of the spur dyke is completely eroded.
8.

Conclusions

A 3D unstructured mesh based RANS solver using
an FVM procedure has been presented in this paper.
With a strict but simple data structure, the solver can takefull advantages of arbitrary polyhedral mesh up to six
faces. Special techniques distinguished from those ofstructured mesh methods have been accentuatedincluding the discretization of the governing equations,
the wall boundary treatments, etc.
The solver has been incorporated into different kindsof turbulence models and applied to flows in different
domains with different meshing strategies. The
comparison of the computational result and theexperimental data demonstrates that the solver is able toreproduce the selected flows with a reasonable accuracy.
Furthermore, the stableness, effectiveness and simplicityof the solver are also observed during the model
verification. It may be accounted as a promising solution
for the turbulence modeling in hydraulic engineering
practices.
Nevertheless, some assumptions and simplifications
have been introduced in the solver, which may become
the bottleneck and limit the application of the proposed
method. These include the omission of the water surface
variation and the assumption made for the wall function
approach. Moreover, higher order discretization schemes
have been avoided in this study in order that the solver is
as stable as possible. This may sometimes become themain error source for the simulation. All these problems
should be taken into account in the future research.

Acknowledgements

This research is sponsored by the Monbukagakusho(Ministry of Education, Culture, Sports, Science andTechnology, Japan) under Grant No. 14350265. The
scholarship provided by the Monbukagakusho to the firstauthor is also gratefully acknowledged.

References

Basara, B. (2004): Employment of the second-momentturbulence closure on arbitrary unstructured grids,
International Journal for Numerical Methods in Fluids,
Vol.44, pp.377-407.
Brüger, A., Gustafsson, B, Lötstedt, P. and Nilsson, J.
(2005): High order accurate solution of theincompressible Navier-Stokes equations, Journal ofComputational Fluids, Vol. 203, pp. 47-91.
Chang, H.H. (1988): Fluvial Processes in River
Engineering, John Wiley & Sons, Inc.
Ferziger, J.H. and Peric, M. (2002): Computational
Methods for Fluid Dynamics, 3rd rev. ed.,Springer-Verlag.
Gatski, T.B. and Speziale, C.G. (1993): On explicit
algebraic stress models for complex turbulent flows,
Journal of Fluid Mechanics, Vol.254, pp.59-78.
Haque, M.D., Koken, M., Constantinescu, S.G. andWeber, L. (2005): Use of a 3D RANS model to predict
stratification effects related to fish passage athydropower dams, http://www.iihr.uiowa.edu/people/
resume/constantinescu/NEW_PAPERS_WEB/SEOUL_Haque.pdf
Hassan, O., Probert, E.J., Morgan, K. and Weatherill,
N.P. (2000): Unsteady flow simulation using
unstructured meshes, Computer Methods in Applied
Mechanics and Engineering, Vol. 189, pp. 1247-1275
Imamoto, H., Ishigaki, T. and Kajima, A. (1987):
Velocity vector measurement by Laser Doppler
Velocimeter in an open channel flow (3), Annuals ofDisas. Prev. Res. Inst., Kyoto Univ. No. 30 B (2), pp.
643-654. (in Japanese)
Ishigaki, T. (1993): Study on 3D Flow Structures in
Compound Open Channels, Ph. D Thesis, Kyoto Univ.
(in Japanese)

Johnston, H. and Liu, J.G. (2004): Accurate, stable andefficient Navier-Stokes solvers based on explicit
treatment of the pressure term, Journal ofComputational Physics, Vol. 199, pp. 221-259
Jongen, T. (1997): Simulation and Modeling of
Turbulent Incompressible Fluid Flows, Ph. D Thesis,
Swiss Federal Institute of Technology.
Khaleduzzaman, A.T.M. (2004): Experimental Study onRiver Course Stabilization and Restoration by Using
Groin-like Structures, Master thesis, Kyoto Univ.
Kim, S. E., Mathur, S. R., Murthy, J. Y. and Choudhury,
D. (1997): A Reynolds-averaged Navier-Stokes solverusing an unstructured mesh based finite-volumescheme, Fluent Technical Notes, (TN117) Fluent, Inc.
Kimura, I. and Hosoda, T. (2003): A Non-linear k- ε
Model with Realizability for Prediction of Flowsaround Bluff Bodies, International Journal for
Numerical Methods in Fluids, Vol. 43, pp. 818-837.
Launder, B.E. and Spalding, D.B. (1972): Lectures inMathematical Models of Turbulence, Academic Press.
Li, B. and Fleming, C.A. (2003): Three-dimensionalhydrodynamic model for free surface flow, Journal of
Hydraulic Research, IAHR, Vol.41, No.4, pp. 367-377.
Mavriplis, D.J. (1996): Mesh generation and adaptivity
for complex geometries and flows, Handbook ofComputational Fluid Mechanics, Peyret, R. (ed.),Academic Press, pp.417-459.
Mavriplis, D.J. (1997): Unstructured grid techniques,
Annual Review of Fluid Mechanics, Vol. 29, pp.
473-514.
Muto, Y., Imamoto, H. and Ishigaki, T. (2002): Velocity
measurements in a straight open channel with arectangular embayment, Proc. of the 12th APD-IAHR,
Bangkok, Thailand.
Naot, D. and Rodi, W. (1982): Calculation of secondary
currents in channel flow, Journal of Hydraulics
Division, Proc. ASCE, Vol. 108, pp. 948-968.
Nakahashi, K., Sharov, D., Kano, S. and Kodera, M.
(1999): Applications of unstructured hybrid grid
method to high-Reynolds number viscous flows,
International Journal for Numerical Methods in Fluids,
Vol. 31, pp. 97-111.
Nakahashi, K., Ito, Y. and Togashi, F. (2003): Somechallenges of realistic flow simulations by unstructured
grid CFD, International Journal for Numerical Methods
in Fluids, Vol. 43, pp. 769-783.
Nallapati, R. and Perot, B.J. (2000) Numericalsimulation of free-surface flows using a moving

unstructured mesh, Proc. of ASME 2000 Fluids
Engineering Division Summer Meeting, June 11-25,
Boston, Massachusetts
Olsen, N.R.B. (2003): Three-dimensional CFD modelingof self-forming meandering channel, Journal of
Hydraulic Engineering, ASCE, Vol. 129 (5), pp.
366-372.
Patankar, S.V. (1980): Numerical Heat Transfer andFluid Flow, McGraw-Hill, New York.
Rhie, C.M. and Chow, W.L. (1983): Numerical study ofthe turbulent flow past an airfoil with trailing edge
separation, Journal of AIAA, Vol. 21, pp. 1525-1532.
Rodi, W. (1980): Turbulence Models and Their
Application in Hydraulics - a State of the Art Review,
University of Karlsruhe, Karlsruhe, Germany.
Rubinstein, R. and Barton, J. M. (1990): Nonlinear
Reynolds Stress Models and the Renormalization
Group, Physics of Fluids A 2 (8), pp. 1472-1476.
Shih, T.H., Zhu, J. and Lumley, J.L. (1995): A new
Reynolds stress algebraic equation model,
Computational Methods in Applied Mechanics andEngineering, Vol. 125, pp. 287-302.
Saad, Y. (2003): Iterative Methods for Sparse Linear
System (Second Edition), Society for Industrial &Applied Mathematics.
Salaheldin, T.M., Imran, J. and Chaudhry, M.H. (2004):
Numerical modeling of three-dimensional flow fieldaround circular piers, Journal of Hydraulic Engineering,
ASCE, Vol. 130, No. 2, pp. 91-100.
van Rijn, L.C. (1993): Principles of Sediment Transportin Rivers, Estuaries and Coastal Seas, Aqua
Publications.
Wilson, C.A.M.E., Boxall, J.B., Guymer, I. and Olsen,
N.R.B. (2003): Validation of a three-dimensional
numerical code in the simulation of pseudo-natural
meandering flows, Journal of Hydraulic Engineering,
ASCE, Vol. 129, No. 10, pp. 758-768.
Wu, W.M., Rodi, W. and Wenka, T. (2000): 3D
numerical modeling of flow and sediment transport inopen channels, Journal of Hydraulic Engineering,
ASCE, Vol. 126, No. 1, pp. 4-15.
Zhang, H., Nakagawa, H., Ishigaki, T., Muto, Y. andKhaleduzzaman, A.T.M. (2005): Flow and bed
deformation around a series of impermeable andpermeable spur dykes, Proc. of the international
conference on monitoring, prediction and mitigation ofwater-related disasters (MPMD2005), Jan. 12-15,
Kyoto, pp.197-202.

三次元非構造格子を用いたレイノルズ方程の有限体積解法について
張 浩*・中川一・石垣泰輔**・武藤裕則*

京都大学大学院工学研究科**

関西大学工学部都市環境工学科
要旨

本研究では，三次元非構造格子を用いたレイノルズ平均ナビエ・ストークス方程式の有限体積解法を示している．異なる格子構造を結合することにより生じるコントロール・ボリュームにおける保存則の問題を，厳密かつ有効なデータ構造とすることで解決する方法を示している．空間差分には指数法が，また時間差分には隠的Crank-Nicolson法が用いられている．表面フラックスの計算にあたっては，チェッカーボード現象を避けるた
めにRhie-Chowによる補間法を用いている．系を構成する方程式群は，最終的にILUTP 法による前処理を伴ったGMRES法を用いて解かれている．本研究による方法を３種の室内実験流れに対して適用したところ，それらの条件の異なるいずれに対してもほぼ満足な流況を再現することが示された．
キーワード:

三次元，有限体積法，非構造格子， k − ε モテル