• CFD, Fluid Flow, FEA, Heat/Mass Transfer
  • +91 987-11-19-383
  • amod.cfd@gmail.com

Combustion Simulation in OpenFOAM

Combustion Modeling - Theory and Numerical Simulation

Combustion is one of the most complex phenomena in engineering which requires to deal with turbulent fluid flow, radiative and convective heat transfer, chemical reaction and pollutant formation, mixing and quenching, particle break-up, evaporation ans mass transfer. Internal combustion engines, industrial furnaces used in metal and cement industries can greatly benefit with increased insight into the combustion system used there. Three primary mode of combustion can be nicely summarized in terms of slides shown below.


Disclaimer:

This content is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OpenFOAM and OpenCFD trademarks.


Some of the resource from the web are listed below. IPR: the ownership and copyright of these documents lie with the author mentioned in the document.


CFD Simulation and Combustion Modeling in OpenFOAM

In addition to standard files and dictionaries required to model any CFD simulation in OpenFOAM, following additional files are required exclusively to model combustion phenomena.
  • Flame Wrinking Parameter Xi refers to the ratio [St/Su] where St is Turbulent Flame Speed and the object Su is Laminar Flame Speed (3 different models: unstrained, equilibrium, transport)
  • All species fractions are given on mass-basis.
  • constant/chemistryProperties - This file under 'constant' folder of the case lists initial chemical time step, ODE solver and solver properties to be used.
  • constant/combustionProperties - Here turbulence chemistry interaction model and model constants are listed. Parameters such as fuel name [e.g. Propane, Methane, IsoOctane] laminar flame speed, equivalence ratio, XiModel, GuldersCoeffs, spark or ignition zone and its location / size / strength / duration, fuel fraction in ignition site ... are specified through this file.
  • constant/reactions - The detailed reaction mechanism in native OpenFOAM format (m3, kmole, Joules, K) or CHEMKIN format (mol, cm3, s, K; cal) is specified through this file. The utility chemkinToFoam can be used to convert CHEMKINIII thermodynamics and reaction data files into native OpenFOAM format.
  • constant/thermo.compressibleGas - Thermophysical data for the chemical species involved in reactant and product are listed here.
  • The initial and boundary conditions, in addition to standard U, p, T files, are specified by additional files such s CH4, CO2, H2O, N2, O2, Ydefault in 0/ directory to store chemical species mass fractions.
  • Ydefault specifies initial and boundary conditions for all species other than O2, N2 and fuel such as CH4, C7H16... that appear in reaction mechanism either on the reactants or products side.
  • ft is the fuel mass fraction i.e. fuel/(fuel + oxidant [usually air]), fu is unburnt fuel fraction.
  • Mean reaction regress variable 'b' in file 0/b specifies non-dimensional temperature boundary condition at all flow inlets: fresh gas [unburnt mixture] b = 1.0 and burnt gas [burnt mixture] b = 0.0. This is based on assumption that the rection is effectively instantaneous so the value jumps from 1 to 0 at the flame front.
  • The object b is defined as b = [TBURNT-GAS - T] / [TBURNT-GAS - TFRESH-GAS]
  • Temperature of unburnt gas at the boundaries and internal field is specified by object Tu in the file 0/Tu.
CHEMKIN to OpenFOAM Format: Note that the chemkinToFoam utility updates the data for Specific Heat Capacity only and the transport data for temperature coefficients of viscosity needs to be updated manually.
CHEMKIN to OpenFOAM Format

Combustion Solvers in OpenFOAM

The solvers related to combustion and chemical reaction are reactingFoam, XiFoam, rhoReactingFoam, PDRFoam (Porosity/Distributed Resistance solver such as for flame propagation over obstacles in HRSG - its development was funded by Shell), sprayFoam (earlier known as dieselFoam), reactingMultiphaseEulerFoam, reactingTwoPhaseEulerFoam. sprayFoma models flows involving compressible reacting flow with Lagrangian evaporating particles (such as diesel engines) in a three-dimensional domain.

Combustion of homogeneoug fuel-air mixture:
Combustion of homogeneous mixture

Combustion of fuel droplets:
Combustion of fuel droplets

Schematic of diffusion flames:
Combustion - Diffusion Flames


fireFoam -- simulation of transport of heat and smoke in fire incidents. Excerpts from <<Verification, validation and evaluation of FireFOAM as a tool for performance design>> - "There are mainly two versions of FireFOAM code. One is released as a solver for transient fire and diffusion flame simulation by Open CFD (an official release). The other is an extended version of the official release, and it is maintained by FM Global consisting of modified libraries, solvers and cases for fire research."

 Fires are essentially diffusion flames where fuel and air are mixing while burning as the fire propagates and increase in size and intesity. Thus, fireFoam is a turbulent diffucion combustion model. The OpenFOAM version includes fireFoam solver with "LES turbulence model" and has capability of simulating Lagrangian sprays such water sprinkler sprays for fire suppression. E.g. combustion>>fireFoam>>les>>flameSpreadWaterSuppressionPanel. As per the source code fireFoam.C, fireFoam is a "Transient solver for fires and turbulent diffusion flames with reacting particle clouds, surface film and pyrolysis modelling." Finite volume discrete ordinates model (fvDOM) is to solve radiation heat transfer equation.

  The desired objective to simulate fire phenomena is to get insight into fire detection, heat transfer and temperature rise of structures and smoke expansion / filling rates. Due to shear scale and hazard of fire, actual testing and measurements are not only expensive sometime impossible.

This simple tutorial case ia a 2D domain where methane is supplied through a samll opening at the bottom face. The combustion model used is infinitelyFastChemistry and no pyrolysisZone is active.
fireFoam smallPool tutorial

The chemical reaction is specified through the file 'reactions'.
fireFoam smallPool reaction CH4


fireFoam smallPool reaction C3H8

In addition to a p_rgh file, this solver needs another similar looking file ph_rgh. As per phrghEqn.H, p = ph_rgh + rho*gh + pRef. The 0/ directory also has a special file named FSDOmega. This is applicable to combustion model FSD which stands for "Flame Surface Density". A detailed description can be found in file FSD.H.
The variation of concentration of methane and CO2 with time and space are shown in the following video.

A tutorial case with 3D version of this phenomena is modeled using a cubical box as shown below.
fireFoam smallPool 3D

Another interesting example under fireFoam is modeling pyrolysis of wood. The computational domain used in tutorial combustion>>fireFoam>>>LES>>flameSpreadWaterSuppressionPanel is described below.
fireFoam pyrolysis
fireFoam pyrolysis


sprayFoam: Diesel Spray Combustion with Lagrangian method
The link of PDF tutorial is given at the beginninng of this page. Two sample contour plot for concentration of heptane at two time steps are shown below.
XiFoam Computational Domain

XiFoam Computational Domain Mesh


XiFoam - premixed and partially pre-mixed turbulent combustion

XiFoam: Homogeneous Combustion of Propane

The computational domain is shown by following images.
XiFoam Computational Domain

XiFoam Computational Domain Mesh

Temperature at various time [in seconds] are shown below.

Temperature: t = 0.005
XiFoam Temperature t = 0.005 [s]

Temperature: t = 0.010
XiFoam Temperature t = 0.010 [s]

Temperature: t = 0.015
XiFoam Temperature t = 0.015 [s]

Temperature: t = 0.020
XiFoam Temperature t = 0.020 [s]


reactingFoam
  • This is a transient, incompressible solver based on PIMPLE algorithm for non-premixed turbulent combustion of gaseous hydrocarbons.
  • The turbulence-chemistry interaction models available for this solver are [a] infinitelyFastChemistry - a single step combustion reaction and [b] diffusion based combustion model - also based on single step combustion reaction. [c] PaSR: Partially Stirred Reactor - a finite rate chemistry model based on both turbulence and chemistry time scales.
  • The compressible version is rhoReactingFoam
  • 0/Ydefault file is used in case many type of species share a common set of boundary conditions i.e. a default set which is useful in case there are many species involved.

Mixing and Combustion of Methane and Oxygen

The computational domain is shown by following images.
reactingFoam Computational Domain

Velocity Magnitude: t = 0.5 [s]
reactingFoam velocity Magnitude

Mass Fraction of Methane: t = 0.5 [s]
Fuel Mass Fraction - raectingFoam

Mass Fraction of Oxygen: t = 0.5 [s]
reactingFoam - Oxygen Mass Fraction

Temperature: t = 0.5 [s]
reactingFoam CO2


combustionProperties file with comments
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1606+                                |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile  {
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      combustionProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
/* Comments based on A Brief tutorial for XiFoam in OpenFOAM 1.7.x by Ehsan 
Yasari and Simulating the combustion of gaseous fuels 6th OpenFoam Workshop 
Training Session by Dominik Christ and files in tutorial folders             */

// Gulders, GuldersEGR, constant, SCOPE (PDRFoam);
laminarFlameSpeedCorrelation Gulders;

//----------------------------------------------------------------------------
/* Partially Stirred Reactor [PaSR] model - 
   coalChemistryFoam, sprayFoam, simpleReactingParcelFoam                   */
/*
combustionModel  PaSR;
active  true;    //False
PaSRCoeffs  {
    Cmix                1.0;
    turbulentReaction   on;     // sprayFoam: set to 'yes';
	
	useReactionRate     true;   // simpleReactingParcelFoam
}
*/
//----------------------------------------------------------------------------
// must be specified, if Gulders/GuldersEGR is selected above
fuel                    Propane;
//fuelFile              "fuels/propane.dat";  //PDRFoam

//Laminar Flame Speed - if laminarFlameSpeedCorrelation = constant
Su                      Su [0 1 -1 0 0 0 0] 0.434;

//model  to calculate laminar flame speed: unstrained, equilibrium, transport
SuModel                 unstrained;

// The ratio [fuel-to-oxidizer ratio] / [stoichiometric fuel-to-oxidizer ratio]
equivalenceRatio equivalenceRatio [0 0 0 0 0 0 0] 1;

// The strain rate at extinction which obtained from the Markstein length by 
// extrapolating to Su --> 0
sigmaExt                sigmaExt [0 0 -1 0 0 0 0] 100000;

/*Methods  to calculate the Xi parameters: [1] fixed (Xi = constant), 
  [2]algebraic: Xi == scalar(1)+
     (scalar(1)+(2*XiShapeCoef)*(scalar(0.5)-b))*XiCoef*sqrt(up/(Su+SuMin))*Reta;
  [3]transport: solve a transport equation for Xi                            */
XiModel                 transport;

//Coefficients used in algebraic model for Xi
//Xi = 1 + [1 + 2 * XiShapeCoef *(0.5 - b)] * XiCoef * sqrt[u' / (Su + Su_Min)] * R_eta
//R_eta = u' / sqrt(epsilon * tau_eta)
//tau_eta = sqrt(mu_unburnt / epsilon / rho_unburnt)
XiCoef                  XiCoef [0 0 0 0 0 0 0] 0.62;
XiShapeCoef             XiShapeCoef [0 0 0 0 0 0 0] 1;

//calculation of velocity fluctuation: u(t) = u_mean + u'
//u' = uPrimeCoef * sqrt(2*k/3), k = TKE (Turbulent Kinetic Energy)
uPrimeCoef              uPrimeCoef [0 0 0 0 0 0 0] 1;

//Coefficients for calculation of laminar flame speed according to the Gulders 
//formulation for specific fuel.
GuldersCoeffs {
    Methane     {
        W               0.422;
        eta             0.15;
        xi              5.18;
        alpha           2;
        beta            -0.5;
        f               2.3;
    }

    Propane     {
        W               0.446;
        eta             0.12;
        xi              4.95;
        alpha           1.77;
        beta            -0.2;
        f               2.3;
    }

    IsoOctane    {
        W               0.4658;
        eta             -0.326;
        xi              4.48;
        alpha           1.56;
        beta            -0.22;
        f               2.3;
    }
}

RaviPetersenCoeffs  {
    HydrogenInAir      {
        TRef            320;
        pPoints         ( 1.0e05 5.0e05 1.0e06 2.0e06 3.0e06 );
        EqRPoints       (0.5 2.0 5.0);
        alpha           ( ( (-0.03   -2.347  9.984  -6.734  1.361)
                            ( 1.61   -9.708 19.026 -11.117  2.098)
                            ( 2.329 -12.287 21.317 -11.973  2.207)
                            ( 2.593 -12.813 20.815 -11.471  2.095)
                            ( 2.728 -13.164 20.794 -11.418  2.086) )
                          ( ( 3.558   0.162 -0.247   0.0253 0    )
                            ( 4.818  -0.872 -0.053   0.0138 0    )
                            ( 3.789  -0.312 -0.208   0.028  0    )
                            ( 4.925  -1.841  0.211  -0.0059 0    )
                            ( 4.505  -1.906  0.259  -0.0105 0    ) ) );
        beta            ( ( ( 5.07  -6.42  3.87  -0.767)
                            ( 5.52  -6.73  3.88  -0.728)
                            ( 5.76  -6.92  3.92  -0.715)
                            ( 6.02  -7.44  4.37  -0.825)
                            ( 7.84 -11.55  7.14  -1.399) )
                          ( ( 1.405  0.053 0.022  0    )
                            ( 1.091  0.317 0      0    )
                            ( 1.64  -0.03  0.07   0    )
                            ( 0.84   0.56  0      0    )
                            ( 0.81   0.64  0      0    ) ) );
    }
}

// yes or no
ignite          yes;    
ignitionSites
(
    {
        location        (0.0 0.0 0.0005);
        diameter        0.003;
        start           0.0;
        duration        0.003;
        strength        2.0;
    }
);

ignitionSphereFraction                1;
ignitionThickness ignitionThickness   [0 1 0 0 0 0 0] 0.001;
ignitionCircleFraction                0.5;
ignitionKernelArea ignitionKernelArea [0 2 0 0 0 0 0] 0.001;

// ************************************************************************* //

thermophysicalProperties file with comments
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1606+                                |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile  {
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      thermophysicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

thermoType  {
    type            heheuPsiThermo; 
    mixture         homogeneousMixture;   
	//singleStepReactingMixture; reactingMixture [sprayFoam]; inhomogeneousMixture;
    
	transport       sutherland;
    thermo          janaf;
    equationOfState perfectGas;
    specie          specie;
    energy          absoluteEnthalpy;     //sensibleEnthalpy;
}

stoichiometricAirFuelMassRatio stoichiometricAirFuelMassRatio [0 0 0 0 0 0 0] 15.675;

//homogeneousMixture: reactants, inhomogeneousMixture: fuel, oxidant
reactants {
    specie   {                            //keyword
        nMoles          24.8095;
        molWeight       29.4649;
    }
    thermodynamics     {
        Tlow            200;              //Lower, Upper and Common temperatures
        Thigh           5000;
        Tcommon         1000;
		// Cp / R = a1 + a2 * T + a3 * T^2 + a4 * T^4 + a5 * T^5, a6:enthalpy offset,a7:entropy offset
        highCpCoeffs    (3.28069 0.00195035 -6.53483e-07  1.0024e-10 -5.64653e-15 -1609.55 4.41496);
        lowCpCoeffs     (3.47696 0.00036750  1.84866e-06 -9.8993e-10 -3.10214e-14 -1570.81 3.76075);
    }
    transport  {                  //Sutherland coefficient, mu = As /[1+Ts/T] * [Ts]^0.5
        As              1.67212e-06;
        Ts              170.672;
    }
}

// inhomogeneousMixture only
/*
oxidant   {
    specie     {
        nMoles          1;
        molWeight       28.8504;
    }
    thermodynamics    {
        Tlow            200;
        Thigh           6000;
        Tcommon         1000;
        highCpCoeffs    ( 3.10131 0.00124137 -4.18816e-07 6.64158e-11 -3.91274e-15 -985.266 5.35597 );
        lowCpCoeffs     ( 3.58378 -0.000727005 1.67057e-06 -1.09203e-10 -4.31765e-13 -1050.53 3.11239 );
    }
    transport    {
        As              1.67212e-06;
        Ts              170.672;
    }
}
*/ inhomogeneousMixture only

//homogeneousMixture: products, inhomogeneousMixture: burntProducts, 
products  {
    specie     {
        nMoles          1;
        molWeight       28.3233;
    }
    thermodynamics     {
        Tlow            200;
        Thigh           5000;
        Tcommon         1000;
        highCpCoeffs    (3.1060 0.00179682 -5.94382e-07 9.04998e-11 -5.08033e-15 -11003.7 5.11872);
        lowCpCoeffs     (3.4961 0.00065036 -2.08029e-07 1.22910e-09 -7.73697e-13 -11080.3 3.18978);
    }
    transport     {
        As              1.67212e-06;
        Ts              170.672;
    }
}

// singleStepReactingMixture; reactingMixture;
/*
inertSpecie              N2;
fuel                     CH4;
chemistryReader          foamChemistryReader;
foamChemistryFile        "$FOAM_CASE/constant/reactions";
foamChemistryThermoFile  "$FOAM_CASE/constant/thermo.compressibleGas";       */

//sprayFoam
//CHEMKINTransportFile     "$FOAM_CASE/chemkin/transportProperties";

//sprayFoam
/*
liquids  {
    C7H16     {
        defaultCoeffs   yes;
    }
}

solids  {
    // none
}
*/
// ************************************************************************* //

Combustion Simulation in IC Engines

coldEnginengineFoam, engineFoam, moveEngineMesh, sprayEngineFoam, engineCompRatio, engineSwirl are some of the utilities available in OpenFOAM. moveEngineMesh only move the mesh and intended to ensure the mesh motions are working properly and to determine if the mesh quality will be acceptable during the compression/expansion strokes. coldEnginengineFoam is engineFoam without combustion. In the solver engineFoam, "no special pre-processing for movement of pistons is required because mesh motion is integrated in the solver where at every time step mesh nodes are moved and mesh topology eventually changed when needed. Thus, a mesh is used for a certain crank angle interval (CA) only. However, no new layer of mesh is generated with setting "engineMesh layered;" in engineGeometry dictionary. The mesh is compressed and expanded as the piston moves towards TDC and BDC respectively.

The main consideration of an engine mesh is to have separate patches for cylinderHead, liner and piston and the appropriate values in the dictionary engineGeometry. engineMesh will still work with additional patches other than the three mentioned above but only the liner and piston patches will move, everything else will stay stationary. Other patches can be moved using utility fvMotionSolverEngineMesh.

engineCompRatio Calculate the geometric compression ratio. As described in the source code, "Note that if you have valves and/or extra volumes it will not work, since it calculates the volume at BDC and TCD".

Geometry associated with tutorial kivaTest is shown below.
Geometry kivaTest engineFoam

The animation of temperature profile is embedded below.
engineGeometry Dictionary
/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1606+                                |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile  {
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      engineGeometry;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Parameters specified in this file are for adjusting the mesh parameters for 
// the mesh motion when piston moves from Top Dead Centre [TDC] to Bottom Dead 
// Centre [BDC] for every crank ange step size specified. The parameters are
// read by engineTime.H

// How engine mesh are to be treated when piston moves?
engineMesh      layered;   //layeredAR;

// L
conRodLength    conRodLength [0 1 0 0 0 0 0] 0.147;

// D
bore            bore [0 1 0 0 0 0 0] 0.092;

// S
stroke          stroke [0 1 0 0 0 0 0] 0.08423;

// c
clearance       clearance [0 1 0 0 0 0 0] 0.00115;

// N
rpm             rpm [0 0 -1 0 0 0 0] 1500;

// Intake swirl specification - helpful when valve geometry is not modelled.
swirlAxis       (0.0 0.0 1.0);

swirlCenter     (0.0 0.0 0.056);

swirlRPMRatio   swirlRPMRatio [ 0 0 0 0 0 0 0 ] 2.0; 

swirlProfile    swirlProfile [ 0 0 0 0 0 0 0 ] 3.11;

/* ----------------------------------------------------------------------------
            |   |    |  |
            |   |    |  |    
    *    |```   ``````  ```| '''''' Clearance = c
    |    |                 |-------
    |    | _______________ |   .
    |    |!               !|  /|\  Piston Motion
   \|/   |!       O       !|   |
    x    |!........\......!|   |
         |          \      |  \|/
         |           \     |   `
         |            \    |
                   .   \
                  /|\   \
                   |    O
                   |_  /
                   | \/ Crank Angle, q: -180 when piston is at BDC.
                   | /   
                   O--------------------->   Crank Centre
   r = S/2, x = c + r *(1 + L/r - cos(q) - L/r * [1 - (sin(q) * r/L)^2]^0.5  */

// ************************************************************************* //
  

Radiation Modeling

Radiative heat transfer is a predominate mode of heat transfer in combustion. The dictionaries Qr, IDefault and G are used for this purpose. Qr is for viewFactor radiation model only [this method cannot handle participating media and is used when just surface to surface radiation is to be calculated]. IDefault and G are for other radiation models. Some versions of OF doesn't support radiative transport within a solid media such as radiation through semi-transparent solid. In such cases, radiation can only be defined for fluid region(s).


Particle injection

The files specifying injection of particles in some olders versions of OpenFOAM are [a] /constant/injectorProperties, [b] /constant/sprayProperties and in version 2.0 or v1606+ it is /constant/sprayCloudProperties. For (solid) particle injection as in multi-phase flows, the dictionary file kinematicCloudProperties is used.
Spray phenomena involves [a] particle breakup - breakupModel / atomizationModel, [b] evaporation and heat transfer -evaporationModel / heatTransferModel, [c] droplet drag - dragModel, [d] spray induced turbulence - dispersionModel, [e] combustion.

Definitions
  • parcel: collection of droplets and is a representation of a gathering of real particles. This consruction is plainly made because that it is almost always too computational demanding to simulate all the real particles.

combustion >> coalChemistryFoam >> simplifiedSiwekAs per the source code, coalChemistryFoam is a transient solver for: [a] compressible, [b] turbulent flow with [I] coal and limestone parcel injections, [II] energy source and [III] combustion. The coal particles are injected as per specified positions mentioned in the file constant/coalCloud1Positions and constant/limestonePositions.
Coal Cloud Positions

Limestone Cloud Positions


Types of Combustion

Combustion: Types


Non-Premixed Combustion

Non-Premixed Combustion


Premixed Combustion

Premixed Combustion


Partially Premixed Combustion

Partially Premixed Combustion


Contact us
Disclaimers and Policies

The content on CFDyna.com is being constantly refined and improvised with on-the-job experience, testing, and training. Examples might be simplified to improve insight into the physics and basic understanding. Linked pages, articles, references, and examples are constantly reviewed to reduce errors, but we cannot warrant full correctness of all content.