Session 6


Chairman: G. Raithby


I.S. Choi and B. E. Milton

School of Mechanical and Manufacturing Engineering

The University of New South Wales, Sydney, Autralia 2052

Dual-fuel engines have existed for many years as a means of using alternative fuels as a distillate replacement. In general, the alternative fuels with the greatest potential are the gases such as propane, (predominant in LPG), methane (the largest constituent of both natural gas and bio-gas) and the alcohols, methanol and ethanol. Natural gas, being the most likely medium term fuel is that under study here. All these fuels have very low cetane numbers. They can therefore be used easily in spark ignition (SI) engines but not in compression ignition (CI) types. While medium (truck) size SI engines exist, there are problems in developing such engines for large cylinder volumes and many of the advantages of CI engines are therefore lost. In addition, during a fuel replacement phase, dedicated alternative fuel engines can only be used on limited routes where the fuelling infrastructure exists and an engine which can operate on both fuels is a preferred option. The solution is to dual-fuel with the alternative fuel being introduced in lean proportions much as in the SI types but with a pilot distillate spray to provide the ignition. One problem with such a system is that the turndown ratio of the distillate injector from full load, diesel only operation to minimum load, pilot operation is well outside the range of conventional injectors. However, new types of electro-hydraulic injectors are now being developed such as that described by Yudanov1, which may help to solve this problem.

The natural gas (or any of the vaporised low cetane number fuels) is normally mixed in the inlet manifold by what is essentially a gas carburettor although, to overcome defects of carburetted natural gas systems, some developments are considering direct injection of the gaseous fuel into the cylinder. While the former is the simplest and most convenient method, it has been noted by, for example, Karim2 that high pre-combustion gas levels can lead to end-gas autoignition type knock similar to that which exists in SI engines in addition to the diesel knock in CI engines which results from the ignition delay of the pilot fuel. The latter is generally thought to be exacerbated by the addition of the gas (Karim et al.3, Nielsen et al.4 and Milton5) although some differences exist on its extent. In addition, the slow burning of methane gas needs consideration in gas injection to vehicle engines and direct gas injection has been restricted to large marine engines.

With all recent development in internal combustion engines, a wide range of numerical models have been of immense importance. The difference between the conventional diesel and the dual-fuel engine lies in the mixed combustion process. Models of both the pre-mixed and injected gas dual-fuel combustion are therefore required. These developments are reported here. In order to better understand and calibrate these models, constant volume combustion bomb tests have also been carried out. This paper deals predominantly with the models developed for the constant volume process as is now described.

Models used with gasoline (SI) and diesel (CI) engines are either thermodynamic, quasi-dimensional types as discussed by Heywood6, with coupled phenomenological descriptions of the more important processes or, more rarely, fully multi-dimensional models. The latter are more demanding computationally but allow greater flexibility in that they can more readily handle such items as the heterogeneous gaseous fuel/air mixture surrounding the diesel spray in a dual-fuel engine. Hence, the dual-fuel combustion model developed here is based on a simplified multi-dimensional, finite volume system which provides a compromise between accuracy and computational time.

The core of the process is the diffusion model which describes the time dependent mixing of the gas phase substances (fuel, air and combustion products in differing proportions) throughout the spatial grid which is based here, for simplicity and relevance to the combustion bomb used in the tests, on cylindrical co-ordinates. The numerical solution was determined using the differential conservation equations of mass, energy and linear momentum, the turbulence of the flow being described by means of the k-e eddy diffusivity model of Launder and Spalding7.


  1. Yudanov, S.V., Development of the Hydraulically Activated Electronically Controlled Unit Injector for Diesel Engines, SAE Paper 952057, 1995.
  2. Karim, G.A., Knock in Dual-Fuel Engines, Proc. Inst. Mech. Engrs., Vol. 181, Pt. 1, No. 20, pp. 453-466, 1967.
  3. Karim, G.A., Jones, W. and Raine, R.R., An Examination of the Ignition Delay Period in Dual-fuel Engines, SAE Int. Fuel & Lubricants Meeting and Expo., Baltimore, Maryland, SAE Paper 892140, 1989.
  4. Nielsen, O.B., Qvale, B. and Sorenson, S.C., Ignition Delay in the Dual-fuel Engine, SAE 870589, 1987.
  5. Milton, B.E., Improving the Performance of Small Dual-Fuelled Engines, Paper 22, 1st NGV Int. Conf. and Exhib., Sydney, Australia, 1988.


J. I. Ramos

Departamento de Lenguajes y Ciencias de la Computación

E. T. S. Ingenieros Industriales, Universidad de Málaga

Plaza El Ejido, s/n, 29013-Málaga, Spain

It is well known that the use of q-methods for the solution of multidimensional advection-diffusion- reaction equations such as those occurring in heat and mass transfer, combustion, fluid dynamics, etc., provides a large system of non-linearly coupled algebraic equations whose solution requires iterative techniques. Iterations may, however, be eliminated by employing linearized q-methods which are based on the discretization of the time derivative and the time-linearization of the nonlinear terms; however, a large system of linearly coupled algebraic equations still results upon the discretization of the spatial derivatives. In order to eliminate the coupling between different spatial directions in linearized q-methods, approximate factorization techniques which reduce a multi-dimensional problem to the solution of a sequence of one-dimensional ones have been proposed. The linear one- dimensional operators may be easily and inexpensively solved by means of a block tridiagonal solver. Unfortunately, these approximate factorization techniques introduce a factorization error due to the factorization of the three-dimensional operator into one-dimensional ones, and the magnitude of this error may be very large unless small time and space steps are employed in the calculations.

In this paper, an iterative, exact-factorization method based on linearized q-techniques is proposed. Since it is an exact factorization method, it does not have approximate factorization errors. Furthermore, since a linearized q-technique is employed, the factorized, one-dimensional operators in each spatial direction are linear and may be easily and inexpensively solved. However, the exact factorization introduces higher-order derivatives and a coupling between the different spatial directions. This coupling has been accounted for by means of a predictor-corrector strategy as follows. In the predictor step, an approximate factorization method which introduces factorization errors is employed, whereas the factorization errors are accounted for in an iterative manner in the corrector step. Moreover, the mixed fourth-order spatial derivatives have been eliminated by using the one-dimensional operators so that only second-order derivatives appear in the one-dimensional operators which are solved in an iterative manner. Furthermore, in order to accurately resolve the steep fronts or boundary layers that may occur, the spatial derivatives have been discretized by means of three-point, compact difference expressions which are fourth-order accurate in space, so that the resulting method is second-order accurate in time and fourth-order accurate in space, and the difference equations are formulated in delta form so that, for problems which have steady state solutions, the asymptotic numerical solution is independent of the time step employed in the calculations if approximate factorization techniques are employed to obtain the numerical solution. Since the approximate factorization errors depend on the curvature of the solution, they are expected to be largest where the curvature of the dependent variables is largest. In addition, since these errors depend on the time step, the steady state solution, if any, depends on the second power of the time step.

The use of three-point, compact difference expressions for the discretization of the spatial derivatives requires that, in addition to the unknown dependent variables, their first- and second-order spatial derivatives be considered as unknowns. This implies that, if U is the vector of dependent variables and n denotes the number of spatial directions, the blocks of the one-dimensional finite difference operators which result upon the factorization of compact, linearized q-methods are 3Nx3N for advection-diffusion-reaction equations compared with NxN for the case that the spatial derivatives are discretized with second-order accurate formula, where N denotes the number of dependent variables. The large increase in the dimensionality of the matrices is compensated by both the ease and inexpensiveness with which one-dimensional operators can be solved and the values of the first- and second-order spatial derivatives that they provide. These derivatives may be used to concentrate the grid points were steep gradients occur. Moreover, the use of the one-dimensional operators at three successive grid points together with the fourth-order accurate finite difference expressions which relate the first- and second-order derivatives to the discrete dependent variables provides seven equations for nine unknowns, i.e., the dependent variables and its first- and second-order derivatives; therefore, the spatial derivatives may be eliminated to obtain a tridiagonal matrix for the dependent variables.

The exact-factorization, compact method proposed here has been used to study the combustion of a mixture which is governed by a one-step irreversible chemical reaction in two-dimensions, i.e., n = 2, by solving the energy and mass conservation equations, i.e., N = 2, subject to homogeneous Dirichlet boundary conditions at the four boundaries as a function of time starting with a uniform mixture composition and a temperature spike at the center of the computational domain which serves as an ignition source. Calculations have been performed with second- and fourth-order accurate discretizations of the spatial derivatives, and with exact and approximate factorization methods in order to determine the influence of the spatial discretization and the factorization errors on the numerical solution. Some sample results are shown in Figures 1 and 2 which show the temperature (u) and fuel mass fraction (v) distributions at different times (t) and the errors that result from the approximate factorization of the two-dimensional operator into a sequence of two one-dimensional ones. Although the approximate factorization errors are second-order in time, their largest magnitude occurs at the steep fronts of both u and v.

Figure 1. u (top left), v (top right), and errors in u (bottom left) and v (bottom right) at t =10. (The top figures were obtained with the implicit, linearized compact method presented here Dt=0.04, Dx=Dy=0.5 and q=0.5. The bottom figures correspond to the difference between the results of a second-order accurate approximate factorization method in both space and time and the compact, linearized technique presented here.)

Numerical experiments indicate that the exact factorization technique proposed here is able to accurately resolve the steep gradients which occur in both the temperature and fuel mass fraction as the combustion front propagates through the combustible mixture and approaches the boundaries of the domain. The differences between the results obtained with second- and fourth-order accurate discretizations of the spatial derivatives are largest at the steep fronts; in fact, fourth-order accurate discretizations provide steeper fronts than second-order ones. Moreover, the iterative, exact- factorization method presented here converges in only two or three iterations and is very efficient because the one-dimensional finite difference operators provide block tridiagonal matrices and their solution was obtained by means of LU decomposition.

Numerical experiments have also been performed in order to determine the effects of q, i.e., the implicitness parameter, and the allocation of the reaction terms to the one-dimensional operators which result from the factorization. These experiments indicate that the factorization errors of the fully-implicit method, i.e., q=1, are second-order in time which is the same order as that of the Crank-Nicolson method, i.e., q=1/2; however, the discretization errors of the former, i.e., first- order accuracy in time, are larger than those of the latter, i.e., second-order, and the fully-implicit method predicts a slightly faster flame speed than the Crank-Nicoloson technique when the combustion front is sufficiently far from the boundaries.

The allocation of the reaction terms to the one-dimensional operators results in errors which are also second-order in time and the differences between different allocations and the one which splits the reaction terms in equal amounts amongst the different spatial operators are largest at the steep combustion front. However, equal allocation of the reaction terms to the one-dimensional operators is computationally more demanding that those which assign the source terms to only one-dimensional operator because the Jacobian matrix of the reaction terms has to be evaluated for the two one- dimensional operators. Furthermore, for symmetric problems such as the one considered here, the unequal allocation of the reaction terms to the one-dimensional operators may result in solutions which are not symmetric because the non-symmetric distribution of the reaction terms increases the splitting errors associated with the factorization of the two-dimensional reaction-diffusion operator into one-dimensional ones, although, if the source terms are allocated to only one-dimensional operator, the other one is of the diffusion type.

Figure 2. u (top left), v (top right), and errors in u (bottom left) and v (bottom right) at t = 20. (The top figures were obtained with the implicit, linearized compact method presented here Dt=0.04, Dx=Dy=0.5 and q=0.5. The bottom figures correspond to the difference between the results of a second-order accurate approximate factorization method in both space and time and the compact, linearized technique presented here.)


This work has been supported by Project No. PB94-1494 from the DGICYT of Spain.


D. Morvan, B. Porterie, M. Larini and J.C. Loraud

IRPHE UMR CNRS 138 60 rue J. Curie Technopôle de Château Gombert

13453 Marseille cedex 13 France email:

IUSTI UMR CNRS 139 5 rue E. Fermi Technopôle de Château Gombert

13453 Marseille cedex 13 France email:


This paper reports numerical simulations of unconfined pool fires represented as a Methane/Air turbulent diffusion flame which develops from a porous burner. The modelling of turbulent combustion is solved using an eddy-dissipation combustion model coupled with a RNG k-e statistical model for the turbulent flow. Radiation heat transfer and soot formation have been taken into account using P1-approximation and transport submodels which permit to simulate main phenomena (nucleation, coagulation, surface growth ...) present in such system. The set of coupled transport equations is solved numericaly using a finite volume method, the velocity-pressure coupling is treated with a projection technique.


The progress performed for the modelling of turbulent combustion permit now to simulate with a more realistic manner pool fires which develop in industrial or natural environment. Field models represent a promissing way to improve knowledge of fundamental physical mechanisms which are associated with fire spread. They must permit for example to understand the interaction between gas flow dynamics (turbulence, hydrodynamics instabilities ...), heat flow (by convection, diffusion and radiation) and chemical kinetics (reaction rate, chemical species mixing ...) [1]. Because of high fluctuation levels observed in turbulent flame, the classical moments method generally used in turbulent flows modelling cannot be used for the calculation of chemical species in turbulent combustion. Various models based generally on physical considerations permit to represent multiple interactions between turbulent structures and reaction rate. The reaction rate for example could be evaluated from characteristic turbulent time scale (calculated with turbulent kinetic energy and its dissipation rate) and the fluctuation of fuel mass fraction mixture (Eddy Break Up model for premixed flames) or the average mass fraction of the more deficient chemical species (Eddy Dissipation Combustion model for diffusion flames). An other approach based on probability density function (Pdf) permits also to take into account of finite rate chemical kinetics, the average mass fractions of chemical species are calculated from simple quadrature between instantaneous species mass fractions exprimed as a function of instantaneous mixing fraction and the Pdf generally evaluated from average mixing fraction and its fluctuation [2,3].

Pool fires and jet flames constitute two regimes which limit the behaviour of diffusion flames. They can be distinguished from dimensionless heat flux Q' and the dependence of flame height:

where rĄ, Cp, TĄ design respectively the density, the specific heat and the temperature of the ambient air, and , g, d the energy release by the chemical reaction, the gravity magnitude and the burner diameter.

If the dimensionless flame height L/d is ploted as a function of dimensionless heat flux Q', we can distinguish three zones:

-for Q’ <1 the flame height L varies inversely proportionally to d 2/3(pool fire regime)

-for 1 < Q’ < 100 the flame height L is independant of burner diameter d

-for 100 < Q’ the flame height varies linearly to the burner diameter d (jet flame regime)

Therefore pool fires regime can be characterized by a dependence of flame height with burner diameter and small dimensionless heat flux Q'. Because of large base dimensions and small fuel mass rate, wildeland fires and industrial fires are allways classified as pool fires.


Field modelling of pool fires is based on the resolution of conservation equations (mass, momentum, energy ...) inside the region including the fire and the solid oxidizer situated near this region. The buoyancy flow generated by hot gas expension due to heat release, is sufficiently important that we can suppose that the regime is locally turbulent. Average procedure, mass, momentum, species (fuel, oxidant, combustion products and inert gas) and energy conservation could be writen using classical Favre procedure [1]. The radiation heat transfer is treated using P1-approximation, which consists to consider the gas mixture as a gray medium [4] with a mean absorption coefficient as a function of mole fraction of combustion products, soot volume fraction and mean temperature [5]. The turbulent flow is evaluated using RNG k-e statistical model which permits to improve significantly the results obtained from classical high or low Reynolds number k- e model for the description of turbulent flow including weak turbulent and recirculating flow regions [6]. Modelling of turbulent combustion consists to define a relation between reaction rate, turbulent variables (k, e...) and average or mean square fluctuations of mass fractions (chemical species, mixing fraction ...). The present calculation is based on Eddy Dissipation Combustion Model (EDCM) which represents an adaptation of the Eddy Break Up model for diffusion flames [7]. In this case the mean reaction rate is evaluated from a characteristic turbulent time scale evaluated from the turbulent kinetics energy and its dissipation rate, and the average mass fraction of the deficient chemical species (fuel, oxidizer or combustion products).

The present study includes the numerical simulation of an unconfined turbulent pool fire in a configuration represented in Figure 1, the dimension of the computational domain is 4 x 4 m. From numerical results the effects of input fuel mass rate upon the gaz flow, the dimension of the flame are studied.


The European Economic Commission is gratefully acknowledged for providing partial funding for this research in the frame of the EFAISTOS project.


  1. Cox, G., Combustion fundamentals of fire. Academic Press, 1995.
  2. Borghi, R., Turbulent combustion modelling, Prog. Energy Combust. Sci., Vol. 14, pp 245-292, 1988.
  3. Libby, P.A., Williams, F.A., Turbulent Reacting Flows, Academic Press,1993.
  4. Fusegi, T., Farouk, B., Laminar and turbulent natural convection-radiation interaction in a square enclosure filled with a nongray gas, Numerical Heat Transfer, Part A, Vol. 15, pp 303-322, 1989.
  5. Kaplan, C.R., Shaddix, C.R., Smyth, K.C., Computations of enhanced soot production in time-varying CH4/Air diffusion flames, Combust. Flame, Vol. 106, pp 392-405, 1996.
  6. Zijlema, M., Segal, A., Wesseling, P., Finite volume computation of incompressible turbulent flows in general co-ordinates on staggered grids, Int. J. Numerical Methods Fluids, Vol. 20, pp 621-640.
  7. Magnussen, B.F., Hjertager, H., On mathematical modelling of turbulent combustion with special emphasis on soot formation and combustion, Proceedings 16th Symposium (International) on Combustion, The Combustion Institute, pp 719-729, 1976.


Future Meetings | Past Meetings | Proceedings on Sale | Related Links

ICHMT World Wide Web Administrator