CEJa
MATHEMATICS
AND COMPUTERS IN SIMULATION
Mathematics
ELSEVIER
and Computers
in Simulation 36 (1994) 203208
A flame sheet calculation of a confined buoyancy laminar diffusion flame R. Villasefior Instituto de Inuestigaciones Electricas, Sistemas de Combustidn., Interior Internado Palmira, 62000 Cuernacaca, Mor., Mexico
Abstract A numerical model of an axisymmetric confined diffusion flame formed between a CH, jet and coflowing air is presented. The model utilizes the flame sheet concept to locate the stoichiometric fueloxygen interface and, hence, the points of heat release. The numerically simulated steadystate flame assumes a nonunity Lewis number. Despite the low Mach number of the flow, the viscous dissipation term is allowed in the simulation. Allowances are made for natural convection effects and variable thermodynamic and molecular transport properties. The model developed solves the boundary layer equations for describing the thermal and aerodynamic fields established in the confined, laminar diffusion flame under the restriction that the outer oxidizer tube diameter be very large in comparison with the fuel jet diameter. Numerical predictions of temperature, velocity, and stable species concentrations are compared with measured data by Mitchell et al. (1980) of a laminar methaneair diffusion flame.
1. Introduction The governing equations that describe a reacting flow are highly coupled and nonlinear. When a sufficiently good starting solution for the discretized equations can be known, then, a converging solution can be achieved in less iterations than if no initial solution is available. Finding a starting estimate solution with a simple but robust computer code is essential to make the modeling computationally less expensive. This practice is vital specially for multidimensional models combining both fluiddynamic effects with finite rate chemistry. The aim of the present work is to develop a computer program for reacting flows that can determine a sufficiently “good” initial solution estimate in two dimensions for velocity, temperature and major species. In the present computer code the set of parabolic equations is written first in matrix form and discretized with a standard finite difference technique. The main difference between the present method and the existing numerical schemes consists of the form in which the finite difference equations are solved when the downstream integration is carried out. For instance, in the methods of Patankar and Spalding [7] and Chen et al. [2], the solution to a given flow at 037%4754/94/$07.00 0 1994 Elsevier SSDI 037%4754(93)E0048A
Science
B.V. All rights reserved
204
R. Villaseiior /Mathematics
and Computers in Simulation 36 (1994) 203208
some axial location is obtained in two steps. In the first step the velocity field is determined separately from the scalar field. The second step requires the computed values of the velocity field for calculating the scalar field. In cases in which the number of dependent variables is very large numerical instability can arise using sequential solution techniques due to the uncoupling of the equations. The salient feature of the newly developed scheme is that the coupled effects are fully taken into account.
2. Model formulation Invoking the boundary layer approximations, the NavierStokes equations can be reduced to a system of parabolic partial differential equations (Kuo [4]) describing the conservation of mass, momentum, energy, and species equations. The transport equations are transformed by the Von Mises transformation (Patankar et al. [7]) with a normalized stream function !Jf, so that the computation can be carried out in a domain of fixed coordinates. For atmosphericpressure flames, the flame thickness is small relative to the flame radius and a reasonable representation of flame shape can be obtained by assuming that the reaction between fuel and air is infinitely fast (Mitchell et al. [6]). The infinitely fast chemistry assumption implies that the instantaneous species concentrations are functions only of a conserved scalar (Bilger [l]), the mixture fraction f. This assumption reduces the species conservation equations to one single nonlinear transport equation with no source term. The recovery of the major species profiles follows from the conserved scalar f. The mass fraction concentrations of stable species as a function of f for the CH,Air system have been derived by Keyes et al. [S] and will not be repeated here. The energy equation is expressed in terms of total enthalpy. The radiant energy flux and the energy flux due to mechanical driving forces (Dufour effect) are not considered in this work. To maintain generality in the newly developed numerical scheme the viscous dissipation and the interdiffusion terms, which are normally ignored for low Mach number flows, are included in the energy equation. The temperature profiles are recovered from total enthalpy and the mass fraction concentrations. The idealgas equation of state serves to close the system of equations. The solution to the above set of nonlinear, parabolic, partial differential equations leads to the determination of velocity, mixture fraction and total enthalpy. In matrix representation the axisymmetric flow equations are given by, “gG;=;
+F
Q;
i
(1)
i
Where U = [U fhrlT is the vector of dependent variables and F = [(p pJ/pu source vector. The elements of the diffusion matrix Q in Eq. (1) are written as r
a/J
O
O
1
0 01'is the
R. Villasetior / Mathematics
and Computers
in Simulation
36 (1994) 203208
205
Also in Eq. (l), H is the identity matrix and G is a diagonal matrix whose diagonal elements are all the same and equal to A = (a + bp). The parameters a, b are related to the entrainment rates (Patankar [7]) at the inner (streamline I,!J~)and outer (streamline +,> boundaries. The coefficient (Y in the diffusion matrix Q is defined as p~r’/($~  I,!J~)~. The finite difference representation of each term in Eq. (1) is obtained by considering an integration over a control area AhA = Ax * Aq. Between the grid points, a linear interpolation formula, similar to that proposed by Chen et al. [2] is used when the integration is carried out in the transverse direction. An implicit method is used for all the dependent variables with lagging of the elements in matrices G, Q and F. With the above finitedifference procedure, Eq. (1) can be expressed in terms of algebraic equations at the ith nodal point as follows AIJ”;1 + &!qn+l+ c;u,::’ = Dn. The coefficient matrices A:, I?:, C,? and the independent vector On are given by: Al
=
Aqi
4Ax
l/2
3 Api &! = 4 &J’” cy
=
ATi+
H:,
+ +(3G: + G;,)
q;_1,2,
 b(Gi”l  Gi”,,) + (L/2
I/*
4AX

Hi?
1
+(3G:
(3)
+ q1:&
+ Gi”+l) qi”+1,2,
(5)
1 0” = A~~+ 1/2( HU)~+ 1 + 3A1r,(HU)P + Aqil/z( HU)r I] 4AX [ + a[ A~i+,,,F,“+,
+ 3A~iFFi” + A1v,~,/,F:_I].
(6)
In the above equations Aqi+1,2 = vi+1  pi, and Aqi_,,, = pi  pi_,. The sum of these two delta increments gives Api. The diffusion matrices in Eqs. (36) are defined as, qr+1,2 = (Qi”,l + Qi”)/Api+1/2 and qr_1,2 = /Aqi_,/,*
3. Results and final remarks To demonstrate the importance of nonlinear coupling among the PDE’s, calculations of a laminar, axisymmetric, buoyancy flame for methane has been performed. The thermodynamic and transport properties are calculated with the CHEMKIN and transport subroutines (Kee et al. [3]). Measurements of a confined laminar diffusion flame by Mitchell [6] are used for verification of the model. Initial profiles for mixture fraction, total enthalpy, and velocity are assumed. Eighty grid points across the jet are selected to solve this flow. The downstream integration is carried out up to 3.94 jet diameters downstream. It takes 11 minutes of CPU time to solve completely this flow on a microcomputer (VAX3100). Fig. 1 shows a comparison between experimental data and predicted radial temperature distribution at x/d = 0.94 fuel jet diameters downstream. Analysis of the temperature profile in Fig. 1 shows that the model overpredicts temperature in the fuel rich zone and that the profile is shifted toward the jet centerline. This shift is attributed to the actual finite rates of
206
R. Wllasetior /Mathematics

and Computers in Simulation 36 (1994) 203208
I
0.0
015
110
115
r/d Fig. 1. A comparison x/d = 0.94.
between
experimental
data and numerical
predictions
of radial
temperature
and velocity
at
reaction in the reaction zone and hence, a spreading of the heat release due to chemical reaction. The overprediction of temperature inside the fuel rich zone is caused mainly to the neglect of energyabsorbing pyrolysis reactions. It is important to point out that the numerical results predict a much higher temperature value at the flame front than the measured temperature. The reason for that is due to the elimination of the radiation heat flux from the energy equation, hence, the predicted flame temperature is about 2300 K, a value that corresponds to the adiabatic flame temperature of methane. Also shown in Fig. 1 is the predicted velocity profile and that measured by Mitchell et al. [6]. The expansion during combustion explains the peak in velocity with a gradual decaying behaviour toward the outer flow velocity. Farther up the velocity profile continues to accelerate at the jet centerline due to the buoyancy effects. Fig. 2 compares radial experimental and computational profiles for the major species in the flame at 1.2 cm above the burner plate. Several features in Fig. 2 are worth noting. Clearly, from Fig. 2 the methane diffuses to the reaction zone where it is completely consumed. In the model, the only region for fuel consumption of methane is at the fueloxygen interface. Since allowance is not made for any fuel pyrolysis product, nor for the production and subsequent consumption of any intermediate oxidation product, the methane concentration is overpredieted. The oxygen concentration outside the flame region is near its inlet value and then drops nearly to zero in the reaction zone. The radial mole fractions of water and carbon dioxide maximize in the region of the peak temperature. Due to space limitations profiles for temperature, velocity and mass fractions of major species at other axial locations are not presented. However, a similar trend is observed at those stations.
R. Villaseiior /Mathematics
and Computers in Simulation 36 (1994) 203208
207
0.4
CH4
0.0
Fig. 2. Radial = 0.94.
concentration
distributions
1 .o
0.5
of the
predicted
and
measured
mole
fractions
of stable
x/d
model
developed
the thermal
and
and molecular transport properties been considered. The has been relaxed so that the flux due to concentration gradients is allowed the stimulation. Although the has a Mach number, the the viscous dissipation The advantages the present numerical scheme older methods the full coupling the equations and the flexibility this approach add more transport equations the structure the code.
References [ll R.W. Bilger, 121 J.Y.
Turbulent
flows with nonpremixed
in: P.A. and F.A. Williams, New York, and R.W. Dibble, Numerical for a staggered grid system, the Eighteenth Annual Pittsburgh
[3] R.J. Kee, J.A. Miller and code package, Sandia [4] Kuo, Gaseous diffusion 2nd. (Wiley, New York,
Jefferson, CHEMKIN: A generalpurpose problemindependent, Report, 808003 (1980). and combustion fuel droplet, in:
fortran
208
R. Villaserior /Mathematics
and Computers
in Simulation
36 (1994) 203208
[5] D. Keyes, D. Philbin and M. Smooke, Modification and improvement of software for modeling multidimensional reacting fuel flows, Technical Report, WRDCTR892056 (1989) l92. [6] R.E. Mitchell, A.F. Sarofim and L.A. Clomburg, Experimental and numerical investigation of confined laminar diffusion flames, Combustion and Flame 37 (1980) 227244. [7] S.V. Patankar and D.B. Spalding, Heat and Mass Transfer in Boundary Layers, 2nd. Ed. (Intertext Books, London, 1970).