Application of a Primitive Variable Newton's Method for the Calculation of an Axisymmetric Laminar Diffusion Flame

Application of a Primitive Variable Newton's Method for the Calculation of an Axisymmetric Laminar Diffusion Flame

JOURNAL OF COMPUTATIONAL PHYSICS 104, 99-109(1993) Application of a Primitive Variable Newton's Method for the Calculation of an Axisymmetric Lamina...

853KB Sizes 0 Downloads 2 Views


104, 99-109(1993)

Application of a Primitive Variable Newton's Method for the Calculation of an Axisymmetric Laminar Diffusion Flame YUENONG





Department of Mechanical Engineering, Yale University, P.O. Box 2157, Yale Station, New Haven, Connecticut 06520-2157 Received April12, 1991; revised March 12, 1992

laminar diffusion flames requires the calculation of a large system of highly nonlinear elliptic equations. Newton's method has been proven to be a robust and efficient approach in the numerical calculation of such a system [17]. Recently, a few numerical solutions of axisymmetric laminar diffusion flames with different levels of chemical kinetic models have been reported [12, 17]. The numerical approaches used in these calculations were limited to the stream function-vorticity formulation. Such a formulation was adopted since the determination of the pressure in these systems is difficult due to its first-order nature in the momentum equations and its absence. in the continuity equation. One advantage of the stream function-vorticity formulation is that it eliminates the pressure as a dependent variable from the governing equations. However, this reduction brings with it some side effects, The correct boundary conditions for the vorticity are difficult to determine in the case of complex flows [14]. The vorticity boundary conditions specified from the stream function, e.g., at a solid wall, often cause trouble in obtaining a converged solution, The uncertainty of the vorticity boundary conditions at the inlet of the computational domain may result in a rough approximation of the true solution, thus altering the solution of the stream function whose derivatives are used in forming the velocity, In diffusion flames, combustion is controlled primarily by the rate at which the fuel and oxidizer are brought together in stoichiometric proportion. Thus, an accurate representation of the flow field is a precondition for the overall accuracy of the solution. In addition, the eliminated pressure may be desired as a solution during the iteration to evaluate thermodynamic properties, transport coefficients, and reaction rates (e.g., Lindemann fall-off effect) in some cases. Furthermore, in many problems, it is easier to estimate a reasonable velocity and pressure field rather than a stream function and vorticity field. Finally, a stream function-vorticity approach intrinsically limits the modeling of reacting flows to two-dimensional configurations including axisymmetric problems. While a vector-

In this paper we present a primitive variable Newton-based solution method with a block-line linear equation solver for the calculation of reacting flows. The present approach is compared with the stream function-vorticity Newton's method and the SIMPLER algorithm on the calculation of a system of fully elliptic equations governing an axisymmetric methane-air laminar diffusion flame. The chemical reaction is modeled by the flame sheet approximation. The numerical solution agrees well with experimental data in the major chemical species. The .comparison of three sets of numerical results indicates that the stream function-vorticity solution using the approximate boundary conditions reported in the previous calculations predicts a longer flame length and a broader flame shape. With a new set of modified vorticity boundary conditions, we obtain agreement between the primitive variable and stream function-vorticity solutions. The primitive variable Newton's method converges much faster than the other two methods. Because of much less computer memory required for the blockline tridiagonal solver compared to a direct solver, the present approach makes it possible to calculate multidimensional flames with detailed reaction mechanisms. The SIMPLER algorithm shows a slow convergence rate compared to the other two methods in the present calculation. © 1993 Academic Press. Inc.


Many gas turbines and commercial burners employ diffusion flames as their primary flame type. The ability to predict the coupled effects of transport phenomena with complex chemical process in these systems is critical in improving engine efficiency and in understanding the emission process. The solution of laminar flames can also be employed in the calculation of turbulent reacting flows using the laminar flamelet model, which interprets the turbulent flame as being comprised of an ensemble, of laminar flamelets, The equations governing laminar flames are strongly coupled together as a result of the interaction between the fluid transport, the heat transfer, and the chemical processes, The equations are also characterized by the presence of stiff source terms [13]. The solution of multidimensional 99

0021-9991/93 $5.00 Copyright © 1993 by Academic Press, Inc. All rights of reproduction in any form reserved.



potential vortIcIty model can be formulated for threedimensional flows, they introduce additional dependent variables [1]. The application of Newton's method to a primitive variable formulation of two-dimensional non-reacting viscous flows has been reported previously [22-24]. In these studies, the linear Newton equations were solved by a direct matrix solver. The reason for this is that some zero diagonal elements were present in the Jacobian matrix as a result of the absence of the pressure in the continuity equation. This prevented the use of a block-line iterative solver without pivoting in the diagonal blocks. Computations similar to those reported in [22-24] require a great amount of main memory and are possible only on a supercomputer. In the numerical modeling of multidimensional reacting flows, the computer memory size and time demanded by Newton's method, even with an iterative linear equation solver, is already near the computational limits of some supercomputers. For example, the number of nonzero storage locations in the Jacobian matrix formed in a Newton iteration in a nine-point finite difference scheme is given by 9xNxNxNODE, where Nis the number of unknowns and NODE is the number of mesh points in the computational domain. It is clear from this discussion that it is essentially infeasible to apply the Newton direct solver approach in modeling multidimensional reacting flows even with a simplified chemistry model. The difficulties described so far can be resolved by introducing Newton's method with a primitive variable formulation in such a way that the Jacobian matrix can still be solved by a block-line tridiagonal iterative method. The advantages of this numerical approach are: (1) Newton's method is particularly robust in handling the nonlinearity and the coupled effects among the equations of reacting flows; (2) the primitive variable formulation is suitable for complex flows and threedimensional problems; (3) the use of a block-line tridiagonal iterative method for the Jacobian matrix greatly reduces the requirement on the computer memory and thus makes it possible for large scale computations, e.g., threedimensional problems. In this paper we employ a primitive variable Newton's method to calculate a fully coupled elliptic model of an axisymmetric, methane-air, laminar diffusion flame in a confined co-flowing jet (see Fig. 1). The chemistry is approximated by a flame sheet model [2, 17]. The transport coefficients are evaluated with simple empirical formulae. The purpose of using a flame sheet model in the present work is twofold: the performance of the present numerical method is compared with Newton's method in the stream function-vorticity formulation and the SIMPLER algorithm [12] on the basis of this model; the solution of the flame sheet model can be used as an initial estimate for calculations using detailed chemistry models [11,17,25]. The Patankar-Spalding SIMPLER algorithm

FIG. 1.

Schematic of a diffusion flame in a confined co-flowing burner.

solves the primitive variable equations and has been quite successful in solving non-reacting flows, but it often fails or becomes prohibitively slow in calculating reacting flows associated with high nonlinearity in the source terms. In the next section, the flame sheet problem is formulated using the primitive variables and the stream function-vorticity approach. The numerical methods are presented in Section 3 and the numerical results are discussed in Section 4. Some final comments are contained in Section 5. 2. Problem Formulation

The laminar diffusion flame considered in the present work is produced by an axisymmetric co-flowing jet in which the fuel flows in the center jet with the air stream from the surrounding jet (see Fig. 1). The jet is confined by a shield wall. A series of extensive experimental studies have been made on this burner [12, 13], which makes it a good test case for the calculation of two-dimensional laminar diffusion flames. The previous numerical modeling of similar burners ranged from a flame sheet model to a complicated chemical reaction mechanism involving 16 chemical species equations [12, 17]. All previous elliptic solutions for this burner configuration were obtained using the stream function-vorticity formulation.



Our purpose in the present work is to demonstrate that Newton's me~hod in the primitive variable setting is computationally feasible and can perform well against Newton's method when applied to the stream function-vorticity formulation and the SIMPLER algorithm. In the chemical reaction model that follows, a flame sheet approximation is employed. 2.1. Primitive Variable Formulation of the Flame Sheet Model

The flame sheet concept was proposed by Burke and Schumann [2]. The flame sheet model of a laminar diffusion flame describes the chemical reaction with a one-step global irreversible reaction. The reaction is assumed infinitely fast and limited to a very thin exothermic reaction zone which separates the fuel and oxidizer. In the reaction zone, the fuel and oxidizer react in stoichiometric proportion. With such an approximation, no fuel is present on the oxidizer side and vice versa. In the presence of an inert gas (N), the reaction of fuel (F), and oxidizer (X) to form product (P) can be written in the form vFF+ vxX + N -+ vpP+ N, (2.1 ) where vF, vx, and vp are the stoichiometric coefficients of the reaction. To further simplify the governing equations [17], we assume that (1) thermal diffusion is negligible, (2) the specific heats, cp and CPk (k = F, X, P, N), are constant, (3) the ordinary mass diffusion velocities obey Fick's law, and (4) the Lewis numbers of all the species are equal to unity. With these approximations, the energy equation and the species equations become similar mathematically. By introducing Shvab-Zeldovich variables, the solutions of the temperature and the species can be formulated from the conserved scalar solution of a source-free convectivediffusive equation [17]. The primitive variable governing equations for a flame sheet model of an axisymmetric, laminar diffusion flame in cylindrical coordinates are Continuity,


or (rpv)


+ oz (rpu)=O;


Radial momentum,

o 0 oz (rpuv) + or (rpvv) _

3 or


(/1 OU) oz + r ~ oz (/1 OU)]; Or

(2.4 )

Conserved scalar, oS OS] - 0 ( rpDOS) [ rpv-+rpuor OZ or or

_~ (rPD OS) =0;




Equation of state,

pW p= RT'


k=F, X, P, N.


Of critical importance to the flame sheet model is the ability to recover the temperature and the major species profiles from the conserved scalar solution. If we denote variables at the flame front with the subscript f, then it can be shown that the location of the flame front at the axial coordinate z can be obtained using

(2.7) where the SUbscripts I and 0 refer to the inner jet and the outer jet, respectively. Using the result of Eq. (2.7) we can generate expressions for the temperature and species on the fuel and oxidizer sides of the flame as follows. On the fuel side we have (2.8 )


~ ~ (r/1 OU) _ ~ (r/1 OU) OZ


= _ r op + rpg + [_ ~ ~



3 or

o 0 oz (rpuu) + or (rpvu) 30Z


+ [_~ r ~

Axial momentum,


~ (r/1 OV) _ ~ ~ (r/1 OV)

3 oz



(/1 o(rv)) + ~ (r/1 OV)]; or or oz (2.3 )


(2.11 )

(2.12 )



On the oxidizer side we have (2.13 )

where Pr rer is a reference Prandtl number and A is the thermal conductivity of the mixture. Hence, determination of pD is reduced to the specification of a transport relation for the viscosity, which is approximated by the power law

(2.14 ) (2.23) (2.15 )

Pr rer = 0.75, r = 0.7, To = 298 K, and /10 = 10- 4 gm/cm are reference values for air [9]. The heat release parameter Q/cp in Eq. (2.8) can be determined from where 1.85


an estimate of the peak temperature (e.g., from an experiment) or from the heat of combustion and a representative heat capacity. In the case of two products, i.e.,

2.2. Stream Function-Vorticity Formulation of the Flame Sheet Model we can solve for Y PI and Y P, by forming

The formulation of the stream function and the vorticity equations eliminates the pressure as one of the dependent variables. The vorticity is defined such that

ov au




Y P,


W P2 v P , ) Y p. W PI V PI + W P2 V P,


In the above equations, v and u denote the velocity of the fluid mixture in the radial and axial directions, respectively: S the conserved scalar; Y k the mass fraction of the kth species; p the mass mixture density; W k the molecular weight of the kth species; W the mean molecular weight of the mixture; R the universal gas constant; /1 the mixture viscosity; D the mixture diffusivity; cp the specific heat capacity of the mixture; g the gravitational factor; and Q the heat release per unit mass of the fuel. Equations (2.2 )-(2.6) are solved for the velocities, the pressure, the conserved scalar, and the density. For a given profile of the conserved scalar, we then utilize the relations in (2.7)-(2.20) to obtain solutions for T, Y F , Y x , Y p , and YiI/. The viscosity and the diffusion coefficient depend on the temperature in the chosen transport model. If we introduce the Prandtl number

(2.24 )


The velocity is related to the stream function

IjJ by (2.25)

oljJ pru=a;.


We note that in this formulation the stream function satisfies the continuity equation identically. The governing equations of the flame sheet model in this formulation [7, 17J are:

Stream function, (2.27)

Vorticity, Pr = /1C p A '

(2.21 )

r and recall that all of the Lewis numbers are equal to one, we can write (2.22 )


[aoz (W-; a; OIjJ) a(W OIjJ)] - or -; oz - ~oZ [r 3~oz (!!:.r W)] - ~or [r 3~or (!!:.r W)] (2.28 )



Conserved scalar conservation, 0 [ oz


ot/l) or

0 or


Outer zone (r

v= u=O,

ot/l)] OZ





(2.29 )

where the components of the iso operator are given by (oploz, -oplor). The equations for recovering the temperature and the species are the same as in Eqs. (2.7)-(2.20).

2.3. Boundary Conditions for the Confined Diffusion Flame The computational domain in the co-flowing burner (Fig. 1) is enclosed by the inlet and exit in the axial direction, the symmetric centerline and the solid wall in the radial direction. The inlet of the burner is formed by an inner fuel jet and an outer air jet. The boundary conditions for the primitive variable and stream function-vorticity formulations under the assumption of plug flow for the inlet velocities are as follows: Inlet (z = 0):

r ~ RI ,

S= 1,

v=O, ou w=--


S= 1;





(2.31b) ou

Exit (z -+



00 ),

ou oS -=v=-=O oz oz'

P= Palm'

ot/l = ow = oS = O' oz oz oz ' Axis of symmetry (r


+ !Pouo(R~ - R;),



where the equations marked with "a" are for the primitive variable formulation, and those marked with "b" are for the stream function-vorticity formulation. The subscript np refers to grid points located next to the shield. The specified quantities u I , Uo , YkI' Y ko ' T I , To, PI' Po, and Palm represent the operating conditions of the burner. The primitive variable governing equations require only one pressure boundary condition in each direction. In addition, the vorticity boundary conditions are similar to those in the previous calculations [13, 17]. The effect of the vorticity boundary conditions on the solution is discussed in Section 4.


1 2 1 2 2 t/I="2PIUIRI+"2POuO(r -R I ), w= -


t/I = t/lmax = !PluIR;

~ (rPD OS) _ ~ (rPD OS) = 0,


= R o),

(2.32a) (2.32b)

= 0),

ou oS op v=-=-=-=O or or or '


oS t/I=w= or =0;



The problem formulated in the above section is solved numerically on the computational domain using (1) Newton's method with the primitive variable formulation, (2) Newton's method with the stream function-vorticity formulation, and (3) the SIMPLER algorithm.

3.1. Discretization Method To solve the problem, we first discretize the governing equations on pre-assigned grid cells. The discretization method adopted in the present study is finite differences with a monotonicity preserving upwind scheme [21] for the convective terms on a nine-point stencil, which is formed on a nonuniform tensor product grid. The control volume method is adopted in the SIMPLER algorithm for the convenient arrangement of the weight coefficients in the pressure and pressure correction equations [15]. The control volume method is consistent with the finite difference approach on a nonuniform mesh if the first-order convective terms are discretized using a monotonicity preserving upwind difference method [21]. For example,

-max( -U i +l/Z, O)(Si+l-SJ],




where the max operator is equivalent to the intrinsic function AMAX 1 in FORTRAN.

the new grid to generate the initial solution estimate for the calculation on the new mesh. The grid refinement may be applied several times depending on the particular problem.

3.2. Newton's Method with AdaptiiJe Time Steps The system of discretized equations (see [25] for detail) including an unsteady term is given by (3.2) where D is a scaling matrix and Fis the steady-state residual function evaluated at the solution ¢J. If the unsteady term is implicitly differenced with a backward Euler method, we have

where rand r - I are the solutions at step nand n -1, respectively, and L1 t n = t n - t n - I. The above system of equations can be solved efficiently using a damped and modified Newton iteration [17,25], k=O, 1,2, ... ,


where ),k is the damping parameter at the kth Newton iteration. J(¢Jk) is the Jacobian matrix, (3.5) With the spatial discretizations described in Section 3.1, the Jacobian matrix in (3.5) can be written in block ninediagonal form. For problems involving complex transport and chemistry models, it is often more efficient to evaluate the Jacobian matrix numerically as opposed to analytically [3, 14]. Once the Jacobian matrix is formed, the Newton equations (3.4) can be solved by a matrix solver. Given an initial solution estimate, we follow the transient solution (Eq. (3.3)) until the Euclidean norm of the residual, IIF(¢Jn)112' is small enough so that Newton's method can be employed with the scaling matrix D equal to zero. The size of the time steps in Eq. (3.3) is adaptively chosen by monitoring the local truncation error of the time discretization process [17, 19]. Once the solution is converged on the initial grid, the grid is then either refined or res elected by equidistributing the positive weight functions of the gradient and the curvature of the converged solution. The variation of step sizes can be smoothed by bounding the local (adjacent) and global (overall) step size ratios [6, 10]. The grid refinement is performed independently in the axial and the radial directions [25]. The solution on the previous grid is interpolated onto

3.3. Primitive Variable Formulation with Newton's Method The difficulties with the primitive variable formulation lie in the treatment of the pressure terms in the momentum equations and the velocity components in the continuity equation. If the first-order derivatives of the pressure in the momentum equations are discretized by a central difference operator on a regular grid, a zigzag pressure field will be regarded as realistic by the momentum equations [15]. A similar kind of difficulty arises in discretizing the continuity equation. However, these difficulties can be resolved by employing the staggered grid technique. In the staggered grid, the axial velocity u i ,} is shifted half way from the regular grid point to the point (i + ~,j); the radial velocity Vi,} is similarly moved to the point (i,j + ~), and the rest of the .dependent variables are still kept at the regular point (i,j) (see Fig. 2); The structure of the present staggered grid requires specification of the pressure boundary condition at the outlet of the reactor, which reflects the physical characteristics of internal flows. It is worth mentioning that the shifting of grids for the velocity only moderately increases the programming complexity on a rectangular mesh. A byproduct of the approach is that it increases the accuracy of the discretization for some of the first-order derivatives. For example, the continuity equation can be discretized using central differences at adjacent points. In addition, it enhances the diagonal dominance of the Jacobian matrix due to the compact points at each grid cell. As we discussed in Section 1, a direct matrix solver is still infeasible for the calculation of reacting flows. Therefore, a key factor in the primitive variable Newton's method is the

L x

FIG. 2. Layout of the staggered grids, ..... : u, j: v, 0: other dependent unknowns.



successful application of a block-line iterative method for the solution of the linear Newton equations (3.4). The structure of the nine-strip blocks in the Jacobian matrix should be kept so that the block-line tridiagonal iteration method can be applied efficiently. While a direct matrix solver can be applied to any nonsingular matrix, a block-line tridiagonal iterative method without pivoting can only be applied to a nonsingular matrix without zero diagonal elements. Since the pressure does not explicitly appear in the continuity equation, some zero diagonal elements are present in the Jacobian matrix. As a result, we employ a block-line tridiagonal iterative method with partial pivoting [8]. In this way we only require that the overall matrix and each diagonal block be nonsingular.

Ro = 2.54 cm, and the length of the tubular pyrex shield Z = 30.0 cm. The flow conditions at the inlet (z = 0) are given by


= 0.0 cm/s,

y CH4 = 1.0,


= 0.0 cm/s,


= 4.50 cm/s,

T= 298 K,


Y02= Y H20 = Y C02 = Y N2 =0;

Vz =

9.88 cm/s,

T=298 K,

Y 02 = 0.232,

(4.3 )

Y CH4 = Y H20 = Y C02 =0.

3.4. Stream Function-Vorticity Formulation with Newton's Method

Newton's method applied to the stream functionvorticity formulation is presented in detail in [17]. In the present work, we solve for the vorticity w as one of the dependent variables instead of w/r as in [12, 17]. This gives the advantage of utilizing Dirichlet boundary conditions for the vorticity at the symmetric line and thus avoids the undefined value of wjr at the same boundary. 3.5. The SIMPLER Algorithm

The SIMPLER algorithm described in [15] is modified to calculate compressible reacting flows. Due to the large variation of the density occurring in reacting flows, the pressure and pressure correction equations must be solved exactly to ensure the conservation of mass before calculating the other equations, e.g., the conserved scalar equation (Eq. (2.5)). We initially encountered difficulties in calculating the pressure and pressure correction equations when using a simple tridiagonal method. The time step had to be kept very small in order for the iteration to converge. Due to the relatively small size, the very sparse nature, and the same structure of the pressure and pressure correction equations, a direct sparse matrix solver, YSMP [5], was employed to solve these two equations. As a result, mass is conserved over each grid cell after the velocity correction is made.

The pressure at the exit is p = 1 atm. The shield temperature is approximated at 298 K. The peak burning temperature is estimated to be 2050 K from the experimental data [12]. For simplicity, we denote p-u-v as the primitive variable solution and t{!-w as the stream function-vorticity solution with Newton's method in the following discussion. The numerical solutions are first computed on a 55 x 50 nonuniform grid with all three methods discussed in the previous section. This initial grid is constructed by using algebraic expressions [1] to cluster grid points in the regions of high spatial activity (e.g., near the solid wall, the fuel and air jet boundaries). The converged solutions are used to refine the grid. The grid refinement is applied two times to produce grids of 65 x 65 and finally 75 x 75 points. The solution on the new mesh is easily obtained with a small

1.0 , - - - , - - - , - - - - , , - - - - , - - - - , - - - - , 400




Grid: 75(z}-75(r}

200 >-


'0 o

Grid: lOO(z}-90(r}


:> 100


The overall reaction for a methane-air system with an inert gas (N 2) is 0.0

(4.1 ) The reactor configuration is such that the radius ofthe inner fuel jet R[ = 0.635 cm, the radius of the outer oxidizer jet

0~--1---~1~0===:;:::==::;;2r;;:0==:::::L===~30.0 z (em)

FIG. 3. Computed axial profiles of the conserved scalar and the axial velocity at the centerline (r = O) from two sets of grids, 75(z} x 75(r) grid (solid line), 100(z} x 90(r} grid (dashed line). (The lines are overlapped.)



number of Newton iterations with an initial solution estimate formed by interpolating the old solution onto the new mesh. The computations were all performed on a Multiflow Trace 14/300 computer (128 MB main memory, 107 mips/60 Mflops scalar system). The computation on the 75 x 75 grid requires a total of 8.2 MB, and the whole process of the computation with the grid refinement takes about 2 h of CPU time. Finally, we also obtain the primitive variable solution on a much finer grid (100 x 90), using Newton's method for the purpose of testing grid independence of the results. Figure 3 shows the comparisons of the profiles of the conserved scalar and the axial velocity at the centerline computed on two meshes. One can see that the solution on the grid (75 x 75) overlaps the solution on the finer grid. To compare the performance of the three different numerical approaches, the profiles of the residual norm versus computer CPU time for the calculation on the initial grid (55 x 50) are plotted in Fig. 4. Since different iteration schemes are involved, we used the computing time instead of the number of iterations. From the results in Fig. 4, we can clearly see that the p-u-v Newton's method converges monotonically at a much faster rate than the other two methods. Part of the reason is that the steady Newton iteration (i.e., D=O in Eq. (3.2» can be started after 40 time steps. At this point, the steady Newton's method converges quadratically. This is one of the advantages of the present method compared to the ljJ-w Newton's method in which the unsteady iteration must be kept for a larger number of time steps to maintain convergence stability. This behavior

10- 1


tll ~




¢ -


o -


x -


CH 4


CO 2

+ - N2

~ 10-2


- - p-u-v _ .. - SIMPLER




'IjI-r.J (new b.c.) - - - - 'IjI-r.J (old b.c.)

o FIG. 5. Radial profiles of the species mole fractions at z = 1.2 cm, experiment: CH 4 ( 0), 02( 0), H2 O( D), CO 2 ( x), N 2 ( +); computed: p-u-v Newton's method (solid line), SIMPLER (dot-dashed line), t/I-w Newton's method with the new b.c. (dotted line), t/I-w Newton's method with the old b.c. (dashed line). (The first three lines may be overlapped.)

with the ljJ-w approach was also reported in [12, 17]. The SIMPLER algorithm gives a smooth but a slow convergence rate. This phenomenon is caused by the large variation of the density present in the reacting flow and the equation by equation iteration nature of the SIMPLER algorithm.


10 1

Primitive variables


Stream function-Vorticity SIMPLER algorithm

tll ~





~ 10-2 r;..

.. ~:,. .~.

\"" ---.-.-~.;.;---":;: .....



+ - N2

............ 'IjI-r.J (new b.c.)

-...... .

,, ,


- .. - SIMPLER






---- p-u-v


-.. -....



x - C02

:;::: "

¢ -


- - - - 'IjI-r.J (old b,c.)


\ I I \

, \


Computer CPU time (sec) FIG. 4. Overall convergence history. p-u-v Newton's method (solid line), t/I-w Newton's method with the old b.c. (dashed line), SIMPLER algorithm (dotted line).

FIG. 6. Radial profiles of the species mole fractions at z = 2.4 em, experiment: CH 4 ( 0), 02( 0), H 20( D), CO 2 ( x), N 2 ( +); computed: p-u-v Newton's method (solid line), SIMPLER (dot-dashed line), t/I-w Newton's method with the new b,c, (dotted line), t/I-w Newton's method with the old b,e, (dashed line), (The first three lines may be overlapped.)




10- 1

experiment: 0 - CH4




o - H2O


.: 0

x - CO2 + - N2

:;:; ()




- - p-u-v








.... 'I/I-w (new h.c.) 'I/I-G) (old h.c.)


FIG. 7. Radial profiles of the species mole fractions at z = 5.0 cm, experiment: CH4(O), O 2 (0), H 2 0(O), CO 2 (x), N 2 (+); computed: p-u-v Newton's method (solid line), SIMPLER (dot-dashed line), JjJ-w Newton's method with the new b.c. (dotted line), Newton's method with the old h.c. (dashed line). (The first three lines may he overlapped:)


The radial profiles of the species mole fractions are ploted against the experimental data at three axial locations in Figs. 5-7. The flame sheet model compares favorably with the experimental data for most of the major species. The numerical solution tends to give somewhat higher CH 4

mole fractions in the fuel rich region. Since the penetration of oxidizer into the fuel side is not allowed in the flame sheet model, the discrepancy of O 2 on the fuel side is due to this approximation. The higher H 2 0 value in the experimental data is caused by the moisture carried by the re-entrant flow from the exit due to the recirculation in the burner. The tjJ-w solution with the boundary conditions described in Section 2 (which we denote by "old b.c."), predicts an even higher CH 4 and lower N2 and O 2 near or in the fuel side compared to the p-u-v solution. The radial temperature profiles at the two axial locations are shown in Figs. 8-9. Also plotted in Fig. 10 is the profile of the temperature along the symmetric line. The numerical solutions agree reasonably with the experimental data on the oxidizer side, but they fail to match the experimental temperature on the fuel side. As explained in [13], this may be caused by the infinite reaction approximation in the flame sheet model. The tjJ-w solution (old b.c.), however, gives lower temperatures on the fuel side. This lower temperature causes the higher CH 4 mole fraction mentioned above. The preliminary primitive variable solution of a complex chemistry model shows similar trends as in the flame sheet model-a higher temperature on the fuel side and a consistent set of major species mole fractions. Finally, all of the three approaches over-predict the flame length. The experimental data [12] reported a 5.8-cm flame length, while the p-u-v approach predicts a 7.6~cm flame length. As can be seen in Figs. 8-10, the solution of the tjJ-w approach (old b.c.) gives a much longer flame length (9.7 cm) and a broader flame shape. The flame length is defined as the axial




o experiment

experiment p-u-v






";::I., ..... "p.,

'I/I-w (new b.c.)


'I/I-w (old b.c.)

E 1000



'f/;-w (new






,.,III .....;::I ,.,<0



(old h.c.)






Radial temperature profiles at z = 1.2 em, experimental data Newton's method (solid line), SIMPLER (dot-dashed line), Newton's method with the new h.c. (dotted line), Newton's method with the old h.c. (dashed line). (The first three lines may he overlapped.) FIG. 8.

(0), computed: p-u-v


- -0 00


FIG. 9. Radial temperature profiles at z = 5.0 cm, experimental data computed: p-u-v Newton's method (solid line), SIMPLER (dot-dashed line), ljJ-w Newton's method with the new h.c. (dotted line), Method's method with the old h.c. (dashed line). (The first three lines may he overlapped.) (0),





/",,~o ~

," I I I



"8 :0


:0 :0



o experiment

p-u-v SIMPLER 'I{I-w (new b.e.)



'I{I-w (old b.c.)


z (em) FIG. 10. Axial temperature profiles at the centerline (r = 0), experimental data (0), computed: p-u-v Newton's method (solid line), SIMPLER (dot-dashed line), t/J-w Newton's method with the new b.c. (dotted line), t/J-w Newton's method with the old b.c. (dashed line). (The first three lines may be overlapped.)

location at which the maximum temperature occurs on the symmetric line and the flame shape is formed by the location of the maximum temperature isotherm. The disagreement is largely due to the simple chemistry and transport models employed in the present calculation. The exact causes can only be identified by a simulation using the detailed chemistry and complex transport models and a careful examination of the experimental procedures used to measure the temperature. The disagreement between the p-u-v solution and ljI-w solution with the old boundary conditions is due to the form of the inlet boundary conditions. In the p-u-v approach, the uniform axial velocity profile and ?ero radial velocity are assigned at the fuel and air jets as in the experiment. The variety of boundary conditions on the vorticity and the stream function at the inlet is carefully discussed in [16]. A convenient approximation of the vorticity boundary condition at the burner inlet is wlz=o=O as in [12,17] or those in Eqs. (2.30b) and (2.31 b). In the calculations of the present confined flame, the inlet boundary conditions played an important role in determining the flame shape. We found that the correct velocity boundary conditions can be accurately represented in the ljI-w formulation in the following way. We specify two boundary conditions for the stream function at the inlet, 81j1 8z = rpv==o = 0,

IjI =




pruz=o dr.

(4.4 )

This leaves no boundary condition for the vorticity. The vorticity at the inlet is then calculated by Eq. (2.27). The first boundary condition in Eq. (4.4) can be incorporated in the discretized form of Eq. (2.27) at the inlet. The same arrangement can also be done at the exit of the burner. With this correct boundary condition in the ljI-w approach, we obtain a new ljI-w solution with the iteration starting from the converged solution with the old boundary conditions. This new ljI-w solution is compared with the other solutions in Figs. 5-10. It is clear from these plots that the correct ljI-w boundary conditions result in a solution almost identical to the primitive variable solution. This observation is confirmed by the solution with the ljI-w boundary conditions at the inlet specified from the p-u-v solution. Finally, we note that the p-u-v Newton solution and the SIMPLER algorithm solution overlap each other in Figs. 5-10. This results from the fact that the finite difference scheme with a mono tonicity preserving upwind difference is identical to the control volume discretization scheme discussed in Section 3.


Application of Newton's method to the calculation of the primitive variable governing equations of reacting or nonreacting flows is presented. The difficulty associated with solving the large Jacobian based Newton equations is resolved by using the staggered grid technique and a blockline tridiagonal iterative method with partial pivoting. The present numerical approach is tested against Newton's method with the stream function-vorticity formulation and the SIMPLER algorithm on the calculation of an axisymmetric, methane-air, laminar diffusion flame. The new approach performs much better than the other numerical approaches. The primitive variable formulation yields an accurate solution since it eliminates the uncertainty of the vorticity boundary conditions. Agreement on the solutions between the primitive variable and the stream function-vorticity approach is obtained when a new set ofljl-w boundary conditions is applied. The numerical results of the flame sheet model agree well with the mole fraction profiles ofthe major species from the experimental data. Preliminary calculations indicate that, if the result of the flame sheet model is used to generate an initial solution estimate for a detailed chemical reaction model, significant gains in CPU time can be obtained.

ACKNOWLEDGMENTS This work was supported in part by the Air Force Office of Scientific Research and the U.S. Department of Energy. Discussions with Professor D. E. Keyes are gratefully acknowledged.




12. R. E. Mitchell, A. F. Sarofim, and L. A. Clomburg, Combust. Flame 37, 227 (1980).

D. A. Anderson, J. C. Tannehill, and R. H. Pletcher, Computational Fluid Mechanics and Heat Transfer (McGraw-Hill, New York, 1984). 2. S. P. Burke and T. E. W. Schumann, Indus. Eng. Chem. 20 (10), 998 (1928).

13. R. E. Mitchell, SC.D. thesis, MIT, 1975 (unpublished). 14. G. N. Newsam and J. D. Ramsdell, Harvard University Research Report TR-17-81, 1981 (unpublished).

3. A. R. Curtis and 1. K. Reid, 1. Inst. Math. Appl. 13,121 (1974).

15. S. V. Patankar, Numerical Heat Transfer and Fluid Flow (McGrawHill, New York, 1980), p. 131.

4. P. Deuflhard, Numer. Math. 22, 289 (1974). 5. S. C. Eisenstat, M. C. Gursky, M. H. Schulti, and A. H. Sherman, Yale University Research Report DOC 114, 1975 (unpublished). 6. V. Giovangigli and M. D. Smooke, Appl. Numer. Math. 5, 305 (1989). 7. A. D. Gosman, W. M. Pun, A. K. Runchal, D. B. Spalding, and M. Wolfshtein, Heat and Mass Transfer in Recirculating Flows (Academic Press, London, 1969). 8. A. C. Hindmarsh, Lawrence Livermore Laboratory Report UCID30150, 1977 (unpublished). 9. A. M. Kanury, Combustion Phenomena (Gordon & Breach, New York, 1982). 10. 1. Kautsky and N. K. Nichols, SIAM 1. Sci. Statist. Comput. 1, 499 (1980). 11. D. E. Keyes and M. D. Smooke, 1. Comput. Phys. 73, 267 (1987).

16. P. 1. Roache, Computational Fluid Dynamics (Hermosa, Albuquerque, NM, 1972), p. 139. 17. M. D. Smooke, R. E. Mitchell, and D. E. Keyes, Combust. Sci. Tech. 67, 85 (1989). 18. M. D. Smooke, AIChE 1.32, 1233 (1986). 19. M. D. Smooke, 1. Opt. Theory Appl. 39, 489 (1983). 20. M. D. Smooke, 1. Comput. Phys. 48,72 (1982). 21. G. A. Sod, Numerical Methods In Fluid Dynamics (Cambridge Univ. Press, London, 1985), p. 304. 22. S. C. P. Van Dam, M. Hafez, and J. Ahmad, AIAA 1. 28, 937 (1990). 23. S. P. Vanka, AIAA 1. 24, 462 (1986). 24. V. Venkatakrishnan, AIAA 1. 27,885 (1989). 25. Y. Xu, Ph.D. thesis, Yale University, 1991 (unpublished).