A flame sheet calculation of a confined buoyancy laminar diffusion flame

CEJ-a

MATHEMATICS

AND COMPUTERS IN SIMULATION

Mathematics

ELSEVIER

and Computers

in Simulation 36 (1994) 203-208

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 fuel-oxygen interface and, hence, the points of heat release. The numerically simulated steady-state 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 methane-air diffusion flame.

1. Introduction The governing equations that describe a reacting flow are highly coupled and non-linear. 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 fluid-dynamic 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)E0048-A

Science

B.V. All rights reserved

204

R. Villaseiior /Mathematics

and Computers in Simulation 36 (1994) 203-208

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 Navier-Stokes 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 atmospheric-pressure 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 non-linear 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 ideal-gas equation of state serves to close the system of equations. The solution to the above set of non-linear, 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, “g-G;=;

+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) 203-208

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 finite-difference procedure, Eq. (1) can be expressed in terms of algebraic equations at the i-th nodal point as follows AI-J”;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 + Aqi-l/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. (3-6) 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 non-linear 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 (VAX-3100). 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) 203-208

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 energy-absorbing 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 fuel-oxygen 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) 203-208

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 non-premixed

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 general-purpose problem-independent, Report, 80-8003 (1980). and combustion fuel droplet, in:

fortran

208

R. Villaserior /Mathematics

and Computers

in Simulation

36 (1994) 203-208

[5] D. Keyes, D. Philbin and M. Smooke, Modification and improvement of software for modeling multidimensional reacting fuel flows, Technical Report, WRDC-TR-89-2056 (1989) l-92. [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) 227-244. [7] S.V. Patankar and D.B. Spalding, Heat and Mass Transfer in Boundary Layers, 2nd. Ed. (Intertext Books, London, 1970).