This article was downloaded by: [Moskow State Univ Bibliote]

On: 07 November 2013, At: 21:22

Publisher: Taylor & Francis
Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House,
37-41 Mortimer Street, London W1T 3JH, UK

Numerical Heat Transfer, Part A: Applications: An

International Journal of Computation and Methodology
Publication details, including instructions for authors and subscription information:
http://www.tandfonline.com/loi/unht20

CONTROL VOLUME FINITE-ELEMENT SOLUTION OF A

CONFINED TURBULENT DIFFUSION FLAME
a a a
D. Elkaim , M. Reggio & R. Camarero
a
Department of Mechanical Engineering , École Polytechnique de Montreal , C.P 6079,
Succ,A, Montréal, H3C 3A7, Canada
Published online: 05 Apr 2007.

To cite this article: D. Elkaim , M. Reggio & R. Camarero (1993) CONTROL VOLUME FINITE-ELEMENT SOLUTION OF A CONFINED
TURBULENT DIFFUSION FLAME, Numerical Heat Transfer, Part A: Applications: An International Journal of Computation and
Methodology, 23:3, 259-279, DOI: 10.1080/10407789308913672

To link to this article: http://dx.doi.org/10.1080/10407789308913672

PLEASE SCROLL DOWN FOR ARTICLE

Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained
in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no
representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the
Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and
are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and
should be independently verified with primary sources of information. Taylor and Francis shall not be liable for
any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever
or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of
the Content.
This article may be used for research, teaching, and private study purposes. Any substantial or systematic
reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any
form to anyone is expressly forbidden. Terms & Conditions of access and use can be found at http://
www.tandfonline.com/page/terms-and-conditions
 Numerical Heat Transfer, Part A, vol. 23, pp. 259-279, 1993

CONTROL VOLUME FINITE-ELEMENT SOLUTION

OF A CONFINED TURBULENT DIFFUSION FLAME

D. Elkaim, M. Reggio, and R. Camarero

Depanment of Mechanical Engineering, ~ c o l Polytechnique
e de Montrial, C.R
6079, Succ. A, Montrial, CaMda H3C 3A7
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

The control volume finite-element method in conjunction with the vom'ciry stream func-
tion fonnuhtion is used to produce a numerical solution of the confined turbulent diffu-
sion JIame. The turbulence is described by the k-c model with wall functions near solid
boundaries. The cornbuslion kinetics are defined by a one-step infhite mte chemical
readion. The validiry of the numerical method is assessed by a comparison between
present numerical results and experimental datafrom the litemturefor a hbomtory turbu-
lent diffusion W e .

INTRODUCTION

Combustion chambers may involve fine geometrical details, which numerical

methods must represent with sufficient flexibility. Traditional finite-difference based
methods work with structured orthogonal meshes, which do not provide efficient means
for handling complex geometries or sufficient mesh flexibility. Finite-element methods
provide enough mesh flexibility, and some applications of the finite-element method in
combustion problems have been reported [ l , 21. Another way to handle complex geome-
tries is with a finite-volume or control volume based finite-element method [3], which
combines the easy physical interpretation of control volume based methods with the
inherent capabilities of the finite-element method.
Recently, the authors have shown [4] that it is worthwhile to use the control
volume based finite-element method in conjunction with the vorticity \ stream function
formulation to simulate reacting laminar flows. This interesting avenue deserves to be
continued, and this paper investigates the use of a control volume based finite-element
method for a numerical simulation of turbulent reacting flows. The turbulence will be
accounted for via the standard k-Emodel [5] together with wall functions for a more
adequate representation of near-wall phenomena. Particular attention will be given to the
vorticity boundary conditions as well.

The authors acknowledge the financial support provided for this research by the CRlM (Centre de
Recherche Informatique de Montrhl), the FP (Fondation de Polytechnique), Gaz Metropolitain of Quebec and
the NSERCC (National Science and Engineering Research Council of Canada).

Copyright @ 1993 lhylor & Francis

1040-7782193 $10.00 + .00
 260 D. ELKAIM ET AL.

NOMENCLATURE
constants temperature, 'K
constants velocity components in the x and y
mole fraction directions, respectively, d s
specific heat of species i at constant average velocity components in the x
pressure a n d y directions, respectively, d s
mixture fraction local velocity components in the X and
enthalpy, Ilkg Y directions. respectively, d s
heat of reaction, Ilkg Cartesian coordinates, m
turbulent kinetic energy, m2/s2 local coordinates, m
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

mass fraction of species i exchange (diffusion) coefficient

molecular weight of species i dissipation rate of k, m2/s3
pressure, ~ l m ' laminar viscosity, kgl(m s)
radial cwrdinaa, m mixture density, kglm3
universal gas constant (- 8.314). WI any iutlmown quantity
(mol OK) stream function
source term of the 6 transport equation vorticity

MATHEMATICAL FORMULATION
The equations governing steady two-dimensional axisymmetric turbulent reacting
flow may be divided into two groups for convenience: the fluid motion equations and the
species transport equations with the related combustion model. Although radiation heat
transfer may be important in such flows (especially for large industrial furnaces), it is
not taken into account here.

Fluid Flow Equations

The steady axisyrnmetric equations representing the conservation of mass and mo-
mentum are written via the vorticity stream function formulation (w, $). The vorticity is
expressed by

The stream function $ is expressed in such a way that the continuity equation is
identically satisfied. We then have

rpu - a$
-
ar

rpv - --a$
ax

Using Eqs. (1) and (2). the governing equations for the fluid flow problem become
[6] as follows.
 C O m D TURBULENT DIFFUSION FLAME 261

For the stream function

For the vorticity

a (rpuw) + -
- a (rpuw) -- -
ax ar
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

where a is the effective viscosity, equal to

where p, is the turbulent viscosity evaluated according to the k-e model of turbulence [5]
as follows.

nrbulent kinetic energy

a
- (rpuk)
a (rpvk) -
+- +r(Gk - pe) (6)
ax ar

Dissipation rate of k
a (rpue) + -
-
ax
a (rpve)
ar
-

Gk or c,Gk are generation terms, as opposed to the destruction terms -pe and - cfle.
Gk is given by

The turbulent viscosity is related to k and e via

This turbulence model has five constants, uk,a,, c,, c,, c,, for which the values
0.9, 1.22, 1.44, 1.92, and 0.09 are commonly used [5], respectively.
 262 D. ELKAIM ET AL.

In general, the k-E model is only valid in regions where the flow is entirely
turbulent. Close to the solid walls, viscous effects become dominant, and such a model
does not lead to acceptable predictions. Later in the present work, the wall function
method is used to adequately describe the flow in the vicinity of solid walls.

Combustion Model

We consider only diffusion flames and a fast chemistry with a one-step irreversible
chemical reaction. It is assumed that the chemistry is sufficiently fast and that intermedi-
ate species do not play a significant role. In these physically controlled diffusion flames,
the mixture composition can be related to one conserved scalar quantity, namely, the
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

mixture fraction. In a two-feed system, the mixture fraction is conserved under chemical
reactions and is defined by

where

and m, and m,, denote the mass fraction of fuel and oxidizer, respectively, i is the
stoichiometric oxidant required to bum 1 kg of fuel, and the subscripts F and 0 denote
the fuel and air stream conditions at the inlet.
From the assumption of chemical equilibrium, there will be no oxidant present for
a mixture fraction richer than stoichiometric and no fuel present when the mixture
fraction is below stoichiometric. Neglecting PV work and radiation, and for adiabatic
operation under the assumption of a unit Lewis number, the enthalpy equation is similar
to the mixture fraction equation, and the enthalpy can be calculated directly from the
mixture fraction and inlet enthalpy values. Thus for all the combustion-related variables,
i.e., the enthalpy, the mass fraction of the fuel, the oxygen, and the products, it is
sufficient to solve the transport equation for the mixture fraction f. All the subsequent
unknowns can be deduced from it [7].
Because we have chosen to put more emphasis on the numerical method than on
the combustion model, at the present stage, we will not take into account the chemistry-
turbulence interaction. However, for a more rigorous modeling of a combustion system,
the mixture fraction variance should also be calculated [7].
The variable f is obtained by solving the following transport equation:

where r - -
pJaf and of 0.9.
In order to close the system given by Eqs. (1)-(12). other equations for the ther-
modynamical properties are needed. These are as follows.
 CONFINED TURBULENT DIFFUSION FLAME

For the pressure

For the enthalpy

Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

where 112 v2 represents the mixture kinetic energy.

BOUNDARY CONDITIONS
To solve the governing equations, some boundary conditions are required at the
boundaries of the computational domain. For flow problems, the most common types of
boundaries are inlet, outlet, solid walls, and symmetry axis. For each type and for each
unknown, let us see what boundary conditions are to be specified.

Boundary Conditions for k-e

Solid walls. To adequately describe the steep gradients in the wall region would
require many mesh points. Moreover, the k-c model does not allow good predictions in
the vicinity of the wall because it is primarily suited for fully turbulent flow regions. A
well-known method for compensating this weakness is the wall function method, which
has been adopted here. The region close to solid walls can be divided into two sublayers
(Fig. 1): a laminar sublayer or viscous sublayer where purely viscous effects are domi-
nant, and a turbulent sublayer. Then it is supposed that the first computational point p
adjacent to the wall is in the turbulent sublayer (Fig. 1) and that the velocity vector at
this point is parallel to the wall. In this sublayer, the velocity has a known logarithmic

r main flow .

Fig. 1 Solid wall region.

 264 D. ELKAIM ET AL.

variation [8]. In this case, conditions on the wall are related to conditions at point p
through this law.
The computational procedure thus skips over the viscous sublayer, where many
grid points are necessary to adequately describe the steep gradients in this region. It is
only necessary to make sure that point p is indeed in the turbulent sublayer. Other
types of sublayers can be considered in building the law of the wall [9] with the point p
within any of these sublayers. Naturally, the velocity variation has to be changed
accordingly.
Here we consider only two sublayers (viscous and turbulent), with the first compu-
tational point p being in the turbulent sublayer. At this point the velocity Up is parallel to
the boundary and has a logarithmic variation ([8]):
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

where u* is the friction velocity and y;represents a dimensionless distance from point p
to the wall. They are defined by

parameter ( K - -
where 7, is the shear stress at the wall, K is the von Kannan constant, E is a roughness
0.4, E 9.7 [8]), and yp is the actual distance from point p to the wall.
It is the value of the dimensionless distance yithat sets the limits between the different
sublayers. For the turbulent layer, y; is approximately between 10 and 400 [9].
To obtain relations involving the values of k and E at point p, it is common to
suppose that the turbulent sublayer is in local equilibrium, so that the rate of production
of k is exactly equal to its destruction rate. Therefore, from Eqs. (6) and (8) and taking
into account only the planar two-dimensional terms gives

where y is the distance from the wall. Moreover, considering the zero velocity at the
wall and supposing that the velocity at point p is parallel to this (so that there is no
variation of the v velocity component in the normal direction to the wall), Eq. (18) at
point p becomes
 CONFINED TURBULENT DIFFUSION FLAME 265

Now let us further suppose that the shear stress is constant in the sublayer, so that

Then from Eq. (16),

replacing (aulay)lp from Eq. (21) into Eq. (19) and using the definition of p,, Eq. (9),
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

one finally derives the turbulent kinetic energy at point p:

On the other hand, the combination of Eq. (18) with the logarithmic variation of the
velocity, Eq. (IS), and with kp described by the above relation, yields
u*~
4'- (23)
KYP

Equations (22) and (23) give the values of k and E at point p without solving the
transport equations, Eqs. (6) and (7). In order to obtain u*, Eqs. (15) and (17) are
combined, which yields the following nonlinear equation:

The velocity Up is known from the solution of the vorticity and stream function
transport equations, Eqs. (3) and (4). The values of kp and ep are used to calculate the
turbulent viscosity and serve as boundary conditions (Dirichlet) for the rest of the do-
main (Fig. 2). It is not necessary to calculate k and e at the walls. The viscosity there is
set equal to the laminar viscosity.
Inlet. Values of k or E are not known at the inlet and, if they are not given by
experimental data, some reasonable assumptions can be made. The kinetic energy of
turbulence is estimated according to a certain percentage of the square of the average
inlet velocity:

where ii is the average inlet velocity and A is a percentage.

The dissipation is calculated according to the equation
 D. ELKAIM ET AL.

boundary poiyts for k and r

I computational domain for k and r

I
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

Fig. 2 Calculation domain for the law of the wall

where D is the inlet diameter. The values A

and may vary slightly in the literature.
- -
0.03 and a 0.005 are commonly used

Outlet. We suppose that the flow extends over a sufficiently long domain so that
it is fully developed at the exit section. Thus for any variable @ the condition is

Symmetry axis. Here the radial derivative of the variables is set equal to zero:

Boundary Conditions for the Vorticity and Stream Function

As a preamble to the calculation of the boundary condition for turbulent flows, a
brief account of its application for laminar flows is given. Although similar develop-
ments can be found in the literature [6, 101, this description is considered to be helpful
prior to the introduction of the law of the wall for the vorticity.
Solid walls.
Laminar flow. For the vorticity \ stream function formulation, there is no ex-
plicit Dirichlet-type boundary condition like a zero velocity at the walls imposed to
satisfy the no-slip condition. The value of w to be imposed is implicit because it depends
on the flow itself. Moreover, it is the vorticity that is produced on solid walls that affects
the flow field and that therefore has to be computed with care. For example, as has been
done recently [ll], the Poisson-type equation vZ@ -
- w can be used to compute the
vorticity on solid boundaries. In a more classical way [lo], the vorticity values at the
solid walls are deduced from a 'lkylor series expansion of the stream function $ around
the solid point. Then, depending upon the number of terms of the series that are taken
into account, one gets first-order or second-order or even higher order boundary condi-
tions. For a first-order boundary condition [lo],
 CONFINED TURBULENT DIFFUSION FLAME

while for the second-order boundary condition, the vorticity is given by

where points W and P and the distance Y, are indicated on Fig. 3. Note that Eq. (29) is
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

valid for orthogonal and nonorthogonal grids, so that irregular boundaries can be cor-
rectly treated. It is also valid for variable density flows.
Even if a second-order boundary condition is more accurate, numerical experi-
ments show that a higher order formulation can lead the overall computational procedure
to divergence as the Reynolds number increases. It is therefore more common to use a
first-order boundary condition, as is done here.
For the stream function, we impose a constant value on the walls because they are
also streamlines. In order to find this value, the velocity profile can be integrated at a
section where it is known (inlet), and the value of $ can be fixed on another boundary
(symmetry axis or solid). This allows the integration constant to be calculated.
Turbulent flow. For turbulent flows, the calculation of the stream function re-
mains the same as before. However, the condition for the vorticity has to be changed to
reflect the law of the wall. In this case, we must find a way to impose the shear stress 7,
at the wall according to the logarithmic distribution of the velocity close to the wall. In
order to do this we must first recall that the mesh point adjacent to the wall is in the
turbulent sublayer (Fig. l), and seek one way of forcing the logarithmic variation close
to the wall. When the primitive (u, u, p) formulation is used, a well-known alternative is
to compute a fictitious slip velocity on the wall. It is deduced from the logarithmic
variation as follows: First, Eq. (15) is derived with respect t o y , which here represents a
local coordinate normal to the wall:

then, at point P and using backward differences, the derivative aU/ay is evaluated:

and, finally, the slip velocity Uw is equal to

 D. ELKAIM ET AL.
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

Fig. 3 Calculation o f the vorticity at solid walls

and has the same sign as U p .

We now have to modify the calculation of the vorticity on solid walls by reworking
Eq. (29), taking into account the slip velocity Uw at the wall. This yields an implicit
Dirichlet-type boundary condition of the form (neglecting curvature effects)

Uw is the slip velocity relative to the local frame, illustrated in Fig. 3.

Inlet. At the inlet, thevelocity is known, so that the stream function and the
vorticity distributions are calculated from their own definitions, Eqs. (1) and (2).
Outlet. The type of conditions at the outlet are the same as for k or e , Eq. (24),
whether for $ or w .

ity o is equal to zero because of the conditions aulay - -

Symmetry axis. The stream function is assigned a fixed value, while the vortic-
0 and u 0 on the axis.
Finally, for the combustion-related variables, at the inlet, we specify (mi, T),from

condition aflan -
which the values of the mixture fraction f can be deduced. On solid boundaries, the
0 is imposed, n being the direction normal to the wall.

NUMERICAL SOLUTION

To solve the set of partial differential equations previously described, the control
volume based finite-element technique presented by Baliga and Patankar [3] has been
followed. Each equation in the system can be cast in the following generic form:
 CONFINED TURBULENT DIFFUSION FLAME 269

where 6 represents the scalar, which undergoes convection and is diffused through the
field, r the exchange coefficient, and S, a source or sink term. A description of the
discretization method by reference to this general transport equation follows. More
details may be found in Refs. [3, 12, 131.

Domain Discretization and Interpolation Function

The domain of interest is first divided into three-node triangular elements. Around
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

the computational point P, a control volume is created by joining the centroids of all
neighboring elements through the midpoints of the corresponding sides (Fig. 4).
Following the finite-volume framework, Eq. (35) is integrated over the control
+.
volume. This procedure requires an interpolation function for Baliga and Patankar [3]
developed an approach based on the idea of using the exact solution of the one-
dimensional convection-diffusion equation as the interpolation function, namely, the ex-
ponential scheme proposed earlier by Patankar [14]. As an extension, an interpolation
function is introduced that is as close as possible to the exact solution of the two-
dimensional convectiondiffusion equation. With the origin located at the centroid of the
element, a locally flow-oriented coordinate system (X, Y) is defined (Fig. 5). For each
triangular element, the interpolation function for 6 is given by

Fig. 4 The polygonal control volume

 D. ELKAIM ET AL.
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

Fig. S Definition of flow-oriented coordinate system

with

where X and Y are the local coordinates and Xm,

average local velocity component, equal to
- max (XI, X,, XI). U,,is the

where

r and p are average values of the exchange coefficient and the density that prevail over
 CONFINED TURBULENT DIFFUSION FLAME 271

the element. The values of A , B, and C are uniquely determined by the values of 6
pertaining to the three nodes 1 , 2, 3.
The beneficial characteristics of this shape function have been further demon-
strated by Prakash and Patankar [15] and Hookey and Baliga [13, 161 for fluid flow or
heat transfer problems using the primitive variable formulation or the vorticity \ stream
function formulation [ll]. It must be pointed out that because there are no convection
terms in the partial differential equation for the stream function +, the interpolation
function for this variable is bilinear and given by $ -
ax + by + c. The resulting
numerical scheme is then applied to solve reacting turbulent incompressible flows.

Discretized Equations
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

The discretization of the equations is carried out by integrating Eq. (35) over the
defined control volume. Using Simpson's rule through points a-r-o or o-r-c (Fig. 5)
and applying the divergence theorem, an equation of the form

is obtained for each computational point with a,, a,,, and dp called the discretization
coefficients. Details on how these coefficients can be computed may be found in Refs.
[12] and [13]. The subscript nb refers to neighboring points, and p stands for the compu-
tational point P.

Boundary Nodes and Source Terms

Boundary nodes. For nodes at Dirichlet boundaries, a, - -
1 , a,,
becomes the known value of 6.For Neumann boundaries, an equation of type Eq. (39)
0, and dp

has to be written for each point of the boundary. Consider the boundary point W together
with the triangle containing the normal to the boundary at that point, and suppose a
linear variation of the property 4 in that triangle. Then in the local X, Y frame of Fig. 6,
b has the form

Fig. 6 A boundary node and its related nomenclature

 D. ELKAIM ET AL.

From Eq. (40), the normal derivative is easily evaluated. In fact,

where N, are the base functions associated with the linear interpolation function of 6 in
the triangle containing the normal. Thus, if the normal gradient is known from the
boundary condition, then Eq. (41) is precisely of type Eq. (39). This is very useful when
computations are done in a totally implicit way, as is the case here.
Source terms. Source terms S(6) are first linearized in the following form:
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

where subscripts L and NL stand for the linear and nonlinear parts of the source term.
Then Eq. (42) is integrated over the control volume by assuming average values of S,
and S,, that prevail over the element. Thus

where Vi is the volume (an area in two-dimensions) of a neighbor element (nbe) i, and
the sum is taken over all the neighboring elements of the computational point. Assuming
that an average value mi within an element is given by

where the subscripts 1, 2, and 3 indicate local node numbers, it is possible to calculate
the contribution of the linearized source term to the coefficient a, of the computational
point P, and also to the coefficient a, of its neighboring nodes. This contribution of the i
element, denoted by as,,is given by

Node P and its neighbors will receive this contribution from all neighboring elements.
Note also that the minus sign in Eq. (45) appears because the discretization coefficients
of Eq. (39) are on the left-hand side, while the source term of Eq. (35) is on the right-
hand side.
Additional considerations have to be taken into account for the k and E source
terms in order to avoid overshoots in the solution (negative values of k or E). The
procedure consists in linearizing the negative part of the source term first and then
including it in the discretization coefficient a, on the left-hand side of Eq. (39). For a
reason that will become apparent later, even if the remaining positive part of the source
term can also be linearized and moved over to the left-hand side, it is left on the right-
hand side. The linearization of the negative part is done in the following way.
 CONFINED TURBULENT DIFFUSION FLAME 273

Let S, denote the negative part of the source term, 4 being k or E. S; may be
written as

where f(4) is a function of 4. The linearization is then taken to be

where n stands for the value at the previous iteration.

This practice of taking only the negative part of the source term and including it on
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

the left-hand side of Eq. (39) makes the resulting matrix more diagonally dominant, and
the iterative procedure converges. Obviously, the positive part must remain on the right-
hand side. Note that during the first iteration steps, the values of k and 6 remain nega-
tive, and after that they become admissible.

Calculation of the Velocities

Because the stream function is allowed one linear variation within the element, the
velocities that we can deduce from it are constant over the element:

rpu - a~ -
- b

rpu - --all.ax - -,
However, for comparison purposes, it is often desirable to compute the nodal
velocities. In this case, a good approximation is obtained by computing a weighted-area
nodal velocity at point p, according to

Again, the sum is over all the elements neighboring the computational point P. The u
velocity component is computed in a similar way.

Solution of the Discretized Equations

For each computational point, one equation of the Eq. (39) type is written. These
equations are then assembled to solve the entire field implicitly. The discretized equa-
tions for the vorticity and stream function are solved in a coupled manner (including
boundary conditions), while a segregated approach is used for the rest of the variables.
The resulting matrices are sparse and without any particular structure, and there are a
 274 D. ELKAIM ET AL.

number of suitable methods available to solve the discretized equations. We used a

sparse matrix solver from IBM's ESSL library [17]. The overall solution procedure can
be outlined as follows.

1. Guess the necessary variables.

2. Compute the slip velocity and boundary conditions for k-6.
3. Solve the stream function and vorticity transport equations.
4. Solve the mixture fraction transport equation.
5. Solve the transport equations for k-e.
6. Calculate auxiliary variables such as temperature and density from the asso-
ciated combustion model.
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

7. Update the turbulent viscosity.

8. Treat the updated values of all variables as improved guesses, and return to step
2 and repeat the process until convergence.

NUMERICAL RESULTS

The procedure and the method described above were applied to solve test cases for
(1) turbulent axisymmetric flow in a channel and (2) turbulent axisymmetric confined
diffusion flame. All test cases are accompanied by experimental or analytical data from
the literature, which is used to assess the validity and accuracy of the method.

Turbulent Flow in a Circular Pipe

'hrbulent flow in a circular pipe is a classical and yet essential test for validating
the setting of the turbulence model. We consider here a fluid that enters a circular
channel with a uniform velocity profile and a Reynolds number (based on the diameter
and the average inlet velocity) equal to 1.1 x 10' (Fig. 7). At the exit, the flow is
entirely developed, and we want to predict the radial velocity distribution at this station.
Figure 8 shows a comparison between our numerical results and Nikuradse's
quasi-analytical results [8]. The mesh consisted of 576 elements and 325 nodes uni-
fomdy distributed (13 x 25). Correspondence of the results is excellent.

Confined Diffusion Flame

This test case has been experimentally [I81 and numerically [19] studied. The
general characteristics of this flow are given in Fig. 9. Air and fuel, initially separated,
are injected into an axisymmetric burner, where they mix and react. Because of the
symmetry of the problem, the computational domain is only half the physical domain.
Comparisons are made between our numerical predictions and the experimental data of
the mixture fraction at several stations in the burner and of the combustion products near
the exit.
At all three first-comparison stations (Fig. lo), the value o f f is overpredicted,
especially in the mixing zone of the two jets. Part of this disagreement is attributed to the
difficulty in the adequate physical modeling of turbulent reacting flows. As stated ear-
lier, the turbulence-combustion interaction is not taken into account. Of course, at the
last comparison station for f, concordance is excellent because the mixing and reacting
 CONFINED TURBULENT DIFFUSION FLAME
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

Fig. 7 Global feaNres of the turbulent channel flow

zone is passed. Closer to the solid walls, predictions are better as well. All these results
are also corroborated by the numerical predictions of Ref. [19].
Figure 11 shows the comparison of the molar concentrations of the combustion
products. There are slight differences between the predictions and the experimental data
for H,O and CO, because of the oversimplified one-step five-species chemical reaction.
Indeed, it can be seen that small traces of Ar, CO, and H, are measured experimentally.
Those species are not accounted for by the present combustion model. For the nonreact-
ing species (N, in Fig. l l ) , concordance is excellent, as expected.

CONCLUSIONS
A numerical procedure has been developed for solving turbulent reacting flow.
The velocity-pressure coupling was addressed via the vorticity stream function formula-
tion. Turbulence has been modeled using the k-e model together with wall functions.
The control volume based finite-element method has been used as the numerical method
in conjunction with the exponential interpolation function. Boundary conditions for the
vorticity have been given careful consideration when used with wall functions. To avoid
undershoots, special treatments for the k and E source terms had to be done. Physical
modeling has to be improved to enhance the overall accuracy of the resulting numerical

Fig. 8 Radial velocity distribution in a circular channel. Solid line, present predictions;
squares, Ref. [a].
 D. ELKAIM ET AL.
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

Fig. 9 Global features of the turbulent confined diffusion flame and location of the comparison
stations (drawing scaled in the x direction).

procedure. Nevertheless, the study showed what w e believe is an interesting application

of the control volume finite-element method that also allows treatment of geometric
complexity with enough flexibility.

REFERENCES

I. F. Benkhaldoun and B. Larrouturou, Adaptive Calculations of Wrinkled Flames, Proc. Symp.

6th Int. Symp. on Finite Element Methods in Flow Problems, pp. 145-170, 1986.
2. A. C. Benim, Finite Element Solution of an Enclosed Turbulent Diffusion Flame, Int. I.
Numer. Methods Fluids, vol. 9, pp. 289-303, 1989.
3. B. R. Baliga and S. V. Patankar, A New Finite-Element Formulation for Convection-
Diffusion Problems, Numer. Heat Transfer, vol. 3, pp. 393-409, 1980.
4. D. Elkaim, M. Reggio, and R. Camarero, Numerical Solution of Reactive Laminar Flow by
a Control Volume Based Finite-Element Method and the Vorticity-Stream Function Fonnula-
tion, Numer. Heat 7ban.sfer, Pan B, vol. 20, pp. 223-240, 1991.
5. W. P. Jones and B. E. Launder, The Prediction of Laminarization with a Two-Equation Model
of 'hrbulence, Int. J. Heat Mass Transfer, vol. 15, pp. 301-314, 1972.
6. A. D. Gosman, W. M. Pun, A. K. Runchal, D. B. Spalding, and M. Wolfshtein, Heat and
Mass nansfer in Recirculating Flow, Academic, San Diego, Calif., 1969.
7. E. E. Khalil, Modeling of Furnaces and Combustors. Abacus, Tunbridge Wells, U.K., 1982.
8. H. Schlichting, Boundnry Layer nteory, 7th ed., Mechanical Engineering Ser., McGraw-
Hill, New York, 1979.
9. R. S. Ammo, Development of a lbrbulence Near-Wall Model and Its Application to Sepa-
rated and Reattached Flows, Numer. Heat nansfer, vol. 7, pp. 59-75, 1984.
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013
 278 D. ELKAIM ET AL.
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013

-
Fig. 11 Radial distribution of the species. Solid line, present predictions; squares and
circles. experimental data [18]; axial distance x 1.375 m.

10. 1. P. Roache, Compurario~lFluid Dynamics, Hermosa Publishers, 1972.

I I. C. F. Kettleborough, S. R. Husain, and C. Prakash, Solution of Fluid Flow Problems with the
Vorticity-Stream Function Formulation and Control Volume Based Finite-Element Method,
Numer. Heat Transfer, Part B, vol. 16, pp. 31-58, 1989.
12. B. R. Baliga and S. V. Patankar, A Control Volume Based Finite-Element Method for Two-
~imensionalFluid Flow and Heat Transfer, Numer. Heat Transfer, vol. 3, pp. 245-261,
1983.
13. N. A. Hookey and B. R. Baliga, Evaluation and Enhancements of Some Control Volume
Finite-Element Methods: Pan 1, Convection-Diffusion Problems, Numer. Heat Transfer, vol.
..
14., OD. 255-272. 1988.
14. S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, D.C.,
1980.
..
15. C. Prakash and S. V. Patankar, A Control Volume Based Finite-Element Method for Solving
the Navier-Stokes Equations Using Equal-Order Velocity-Pressure Interpolation, Numer. Heat
Transfer, vol. 8, pp. 259-280, 1985.
16. N. A. Hookey and B. R. Baliga, Evaluation and Enhancements of Some Control Volume
 CONFINED TURBULENT DIFFUSION FLAME 279

Finite-Element Methods: Pan 2, Incompressible Fluid Flow Problems, Nwner. Heat Transfer,
vol. 14, pp. 273-293, 1988.
17. IBM. Engineering and Scienrific Subroutine Library, 4th ed., 1983.
18. D. L. Smoot and H. M. Lewis, 'hrbulent Gaswus Combustion: Pan I, Local Species Con-
centration Measurements, Combust. Flame, vol. 42, pp. 183-196, 1981.
19. D. L. Smoot and P. J. Smith, 'hrbulent Gaseous Combustion: Part LI, Theory and Evaluation
of Local Properties, Combust. Flame, vol. 42, pp. 277-285, 1981.

Received 21 February 1991

Accepred 28 February 1992
Requests for reprints should be sent to M. Reggio.
Downloaded by [Moskow State Univ Bibliote] at 21:22 07 November 2013