CFD Online Discussion Forums - Difficulties with premixed laminar methane flame in openFoam

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM (https://www.cfd-online.com/Forums/openfoam/)

- - Difficulties with premixed laminar methane flame in openFoam (https://www.cfd-online.com/Forums/openfoam/179526-difficulties-premixed-laminar-methane-flame-openfoam.html)

Difficulties with premixed laminar methane flame in openFoam

Hello Foamers,

I have spent the last month trying to get the correct laminar flame speeds for 1D premixed methane fuel and am stuck in two areas. Attached is my current problem setup and some images of the results.

I compared both 0D reactor and 1D flame simulations between cantera and OpeFoam using GriMech, BFER (2 equation), and MP1 (1 Equation) mechanisms. The results of this comparison are shown in Compare_Ignition.jpg and Compare_Flame.jpg (the results for MP1 are similar to the BFER one).
Attachment 51417
Attachment 51418

Since the ignition delay results are matching (exact match for BFER and MP1 where temperature, species, ect overlap) the kinetics seems to be working correctly. However, as you can see the flame speeds are totally not correct. This leaves transport properties and the fluid solution as the culprit.

I know I am having some problems with the boundary conditions. The temperature and velocity field throughout my simulation is shown in OF_Flame_Temperature.jpg and OF_Flame_Velocity.jpg. This leads to my first question. Can anyone explain why the velocity is increasing and actually flowing out of the tube? I think I must have the boundary over defined, but am trying to get an inflow close to the flame velocity (although the flame velocity at the moment is way too high). Can someone take a look at my boundaries and give me a hint what the problem may be?

Attachment 51419
Attachment 51420

The second question may be a little more difficult but I hope the discussion will be useful for others trying to do this kind of simulation.

Assuming that fixing the boundary conditions from above does not greatly improve my flame speeds (which I do not think it will as I have done simulations with constant pressure put into the solver both with and without velocity profile, and the flame speed for BFER and MP1 was still around 1 m/s), I originally thought it must be the assumptions that are in OF for diffusivity and conductivity that are the problem. However, from comparing the flame profile between cantera and OF, the transport parameters look to be too close to me to be the issue.

The temperature, density, and specific heat fields compare very well between the two. Thermal conductivity is also very similar (ambient value is within 5%). Diffusivity is definitely different (for openFoam I plotted mu/rho) however, to me it looks like it is in the wrong direction. The diffusivity is lower than the multidiffusivities of cantera which should make the flame slower not faster?

Attachment 51423
Attachment 51424
Attachment 51425
Attachment 51426
Attachment 51427

So besides the boundary conditions this is where I am stuck. I know other people have removed the assumptions inherent in openFoam (Schmidt number and lewis number) and/or have coupled directly with cantera. However, I would like to prove to myself that this is going to work before going through all that effort. Furthermore, without finding the answer with the assumptions it will be hard to gauge the impact of adding the more detailed approach.

Does anyone have an idea how close the flame speeds should be to experiential values with the constant lewis/schmidt number and thermal conductivity and diffusivity assumptions inherent in OF? I was expecting to get within 10-20% and then was going to look at adding the detailed approach if needed.

Any insight that anyone could provide would be much appreciated.

Also, for anyone running these types of simulations I would recommend learning cantera as it has been very useful for testing to me. If anyone wants any of the simulation files let me know and I will send them or the documentation that I used to get up and running.

Thanks!

Chris

Just an update on this. I realized that I made a mistake on my boundary conditions. Now my inflow has velocity, temperature, and species specified with pressure as zero gradient and outlet has pressure defined and species, temperature, and velocity as zero gradient. This is what I meant to setup in the first place but had the boundary names wrong in my pressure file.

This makes the result look better and the flame speed is more in line with cantera (about 47 cm/s compared 37 cm/s). There still seems to be some issue as the outflow velocity is too high. I believe u_out should be u_in*(rho_in/rho_out) but the outflow velocity is higher than that.

Attachment 51496

Does anyone have any thoughts on why this is the case? I am going to try commenting some parts of the solver and reviewing the boundary code to see if something is going over the whole mesh and scaling the velocity or something. i am not sure what the issue is yet but will post when I figure it out.

Hi, are you using PDRFoam for simulating this? I am actually working on a similar problem. Can you just provide me with any tutorials which I can refer to?

I am actually having difficulty running the PDRFoam tutorial in this link-

https://github.com/OpenFOAM/OpenFOAM...nWithObstacles

I do not understand the READ ME section, which says-

PDR test case Step to introduce the PDR fields: 1) Create zero-size patches for wall or/and coupled baffles in the boundary file. 2) Specify the boundary contitions for these patches in the fields. 3) Create the new PDR mesh using the PDRMesh utility.

The third step is the creation of the mesh using pdrmesh, that,
I can understand. However, I am not able to understand the first
two steps. Can you help me as to how I can run this test case?

Thank you! :)

Hi, I am using coalChemistryFoam, which is derived from reactingFoam. There is not an exact tutorial but you can download the files that I attached to my previous post.

There is a batchRun file that will create the mesh and run the simulation. You will need to edit it to call coalChemistryFoam instead of myCoalChemistryFoam.

Let me know if you have any issues. I am just trying to debug the velocity issue and will post once I figure it out.

On the PDRFoam tutorial, I have not looked at this but maybe someone else can respond?

Hi DustExplosion!

Thanks for the reply. I am currently facing some issues with my system. I will try this in 2-3 days and get back to you.

Thanks! :)

From looking into the results further, I realized that the outflow velocity is correct. It is faster than cantera as my flame speed is higher. If I set the inflow velocity to 0 m/s (let the flame propagate into a stationary mixture) the velocity-continuity relation from the previous post is satisfied for both OpenFoam and Cantera. Attached is an updated run directory with the inflow velocity set to zero. Note that I am using cantera to generate the initial conditions for the openFoam simulation.

Now I am trying to figure out why the flamespeeds are off. I tryed the MP1, BFER, and GRI-MECH reaction mechanisms. The flame speed results look like this:

Attachment 51615

To figure out why the flamespeed is off I looked at all the transport and thermophysical variables. Specific heats are very close. Thermal conductivity seems to be within 5% even with the assumptions used in openFoam. I think the the issue is the constant Lewis number assumption in the diffusivity calculation.

Here is the species profile and diffusivities for the MP1 case:
Attachment 51616
Attachment 51617

It looks like for cantera the H2O is diffusing into the preheat zone and the zone where the maximum reaction is occurring (could not upload this plot). This reduces the maximum reaction rate, but this feature is not captured in OpenFoam. Similarly with the GRI-MECH reaction mechanism the highly diffusive H2 species is diffusing upstream:

Attachment 51620

In this case the H2 increases the reaction rate, which is not captured in OpenFoam (I Think).

So that is where I am at. Next, I would like to implement a constant lewis number on a species basis (difference constant lewis number for each specie) to see if it improves the flamespeed.

Does anyone have a reference for the different constant Lewis numbers that are typically used? Preferable for many species such as the GRI mechanism.

I know there are more detailed approaches but I would like to keep it lightweight as I will be moving to 2D and multiphase in the near future. Also a good reference for laminar premixed simulations with different diffusion approaches would be a useful contribution to this thread if anyone has it.

Thanks!

I realized that I made a mistake in my Sl vs phi plot from before. Attached is the updated version.

Attachment 51673

The gri-mech results compare well for lean concentrations with the standard coalChemFoam solver. I have not looked at it extensively yet, but my guess is the H2 and H diffusion are not as important for fuel lean conditions (likely because the concentrations and therefore concentration gradients are not as high?).

Ok, I am going to try and implement variable lewis number diffusion coefficients and see if that makes the flame speeds better. I think I have the specie equation derivation correct but I have some questions for the energy equation. I will post the derivation here and put any of my questions in red.

For anyone who is interested there is also another post looking at more detailed implementation of variable diffusivity (http://www.cfd-online.com/Forums/openfoam/103227-multi-species-mass-transport-library-update.html) which includes source code. Unfortunately, the source code did not want to download from cfd-online over the last few days so I have not had a chance to look at it yet.

Specie Transport
The general transport equation is

$\frac{\partial \rho Y_{k}}{\partial t} + \frac{\partial}{\partial x}\left( \rho Y_{k} [u + V_{k}] \right) = S_{Y_{k}}$

With the diffusion velocity
$V_{k} =V_{k}^{*}+Vc$

and the correction velocity
$Vc = -\sum_{k=1}^{N}Y_{k}V_{k}^{*}$

Assuming Ficks Law
$V_{k}^{*} = -\frac{D_{k,mix}\nabla Y_{k}}{Y_{k}}$

and using a constant lewis number approach for diffusion coefficient

$D_{k,mix} = \frac{\lambda}{\rho C_{p}}\frac{1}{Le_{k}} = \frac{\alpha_{mix}}{Le_{k}}$

Current OpenFoam Treatment
Currently in coalChemistryFoam all lewis numbers are constant and equal to 1.

Therefore:

$D_{k,mix} =\frac{\lambda}{\rho C_{p}}\$

and

$V_{k}^{*} = -\frac{\lambda}{\rho C_{p}}\frac{\nabla Y_{k}}{Y_{k}}$

The correction velocity goes to zero

$Vc = -\sum_{k=1}^{N}Y_{k}\left(-\frac{\lambda}{\rho C_{p}}\frac{\nabla Y_{k}}{Y_{k}}\right) = \frac{\lambda}{\rho C_{p}}\sum_{k=1}^{N}\nabla Y_{k} = 0$

and the governing equation becomes:
$\frac{\partial \rho Y_{k}}{\partial t} + \frac{\partial}{\partial x}\left(\rho Y_{k} u\right) - \frac{\partial}{\partial x}\left(\frac{\lambda}{C_{p}}\nabla Y_{k}\right)= S_{Y_{k}}$

Variable Lewis Number Treatment
With the variable lewis number the diffusion velocity and correction velocity are:

$V_{k}^{*} = -\frac{\lambda}{\rho C_{p}}\frac{1}{Le_{k}} \frac{\nabla Y_{k}}{Y_{k}} = \frac{-\alpha}{Le_{k}} \frac{\nabla Y_{k}}{Y_{k}}$

$V_{c} = -\sum_{k=1}^{N}\frac{-\alpha}{Le_{k}}\nabla Y_{k}$

and the governing equation is
$\frac{\partial \rho Y_{k}}{\partial t} + \frac{\partial}{\partial x}\left(\rho Y_{k} \left(u + \sum_{k=1}^{N}\frac{\alpha}{Le_{k}}\nabla Y_{k}\right)\right) - \frac{\partial}{\partial x}\left(\frac{\lambda}{ C_{p}Le_{k}}\nabla Y_{k}\right)= S_{Y_{k}}$

Everything seems to make sense there. I am not that familiar with modifying the governing equations in open foam so my first question is can I just compute Vc before calling YiEqn.solve and then add it to phi in the definition like:

Code:

 fvScalarMatrix YiEqn

            (

                fvm::ddt(rho, Yi)

              + mvConvection->fvmDiv(phi + rho*Vc, Yi)

              - fvm::laplacian(turbulence->muEff(), Yi)

              ==

                coalParcels.SYi(i, Yi)

              + combustion->R(Yi)

              + fvOptions(rho, Yi)

            );



            YiEqn.relax();



            fvOptions.constrain(YiEqn);



            YiEqn.solve(mesh.solver("Yi"));

Enthalpy Transport

I have some more questions on the energy side. In my thermopysicalProperties files I am using sensibleEnthalpy

Code:

thermoType

{

    type            hePsiThermo;

    mixture         reactingMixture;

    transport       sutherland;

    thermo          janaf;

    energy          sensibleEnthalpy;

    equationOfState perfectGas;

    specie          specie;

}

and the energy equation defintion looks like

Code:

    fvScalarMatrix EEqn

    (

        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)

      + fvc::ddt(rho, K) + fvc::div(phi, K)

      + (

            he.name() == "e"

          ? fvc::div

            (

                fvc::absolute(phi/fvc::interpolate(rho), U),

                p,

                "div(phiv,p)"

            )

          : -dpdt

        )

      - fvm::laplacian(turbulence->alphaEff(), he)

     ==

        rho*(U&g)

      + combustion->Sh()

      + coalParcels.Sh(he)

      + limestoneParcels.Sh(he)

      + radiation->Sh(thermo)

      + fvOptions(rho, he)

    );

My second question is which "version" of the enthalpy equation is being used here? I am working fro Kuo's Fundamentals of Turbulent and Multiphase Combustion and I think the "version" is total nonchemical enthalpy (sensible + kinetic or $h_{tnc}$ on page 39 in Kuo). I think this is the case as the kinetic energy K is included and we explicitly have a reaction source term. Am I correct that if we were solving total enthalpy including chemical energy there would be no reaction source term?

If I am correct than the general governing equation for total nonchemical enthalpy $(h = h_{tnc} = h_{s} + \frac{u^{2}}{2})$ is:

$\frac{\partial \rho h }{\partial t} + \frac{\partial}{\partial x}\left(\rho h u\right) - \frac{\partial}{\partial x}\left(\lambda \frac{\partial T}{\partial x}\right) - \frac{\partial \tau_{ij}u_{i}}{\partial x} -\frac{\partial p}{\partial t} + \frac{\partial}{\partial x}\left(\rho \sum_{k=1}^{N}h_{s,k}Y_{k}V_{k}\right)= S_{h}$

So I have a couple of questions about this. I think the dp/dt term is included because he.name() != "e". Is this correct?

I feel like I should know the answer to this but where did the $\tau_{i,j}$ term go in OpenFoam?

This is another basic one but does the conversion from laplacian of temperature to convert to the laplacian of enthalpy by the following?:

$\frac{\partial}{\partial x}\left(\lambda \frac{\partial T}{\partial x}\right) =$ $\frac{\partial}{\partial x}\left(\alpha \rho C_{p}\frac{\partial T}{\partial x}\right) = \frac{\partial}{\partial x}\left(\alpha \frac{\partial T \rho C_{p}}{\partial x}\right) = \frac{\partial}{\partial x}\left(\alpha\frac{\partial h_{s}}{\partial x}\right)$

which would mean that $h_{s} = T \rho C_{p}$

And then for the final term on the right hand side, how does it dissapear with constant lewis number? The best I can do is

$\sum_{k=1}^{N}h_{s,k} Y_{k} V_{k} =$ $-\sum_{k=1}^{N}h_{s,k}\frac{\lambda}{\rho C_{p}}\nabla Y_{k} = -\frac{\lambda}{\rho C_{p}}\sum_{k=1}^{N}h_{s,k}\nabla Y_{k}$

but I cannot see how we know the summation of the product of enthalphy and the gradient goes to zero $\left(\sum_{k=1}^{N} h_{s,k} \nabla Y_{k} = 0\right)$

And my last question (for now at least!) is if I want to include variable diffusivities in the enthalpy equation what is the best way to include them in the equation defintion in OpenFoam? For example would I compute:

[math]A = \sum_{k=1}^{N}\left(\frac{-\alpha}{Le_{k}}\nabla Y_{k} \rho T C_{p,k}\right)[/imath]

and include it in the equation definition something like the following?

Code:

    fvScalarMatrix EEqn

    (

        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)

      + fvc::ddt(rho, K) + fvc::div(phi, K)

      + (

            he.name() == "e"

          ? fvc::div

            (

                fvc::absolute(phi/fvc::interpolate(rho), U),

                p,

                "div(phiv,p)"

            )

          : -dpdt

        )

      - fvm::laplacian(turbulence->alphaEff(), he)

   

     + fvc::ddx(rho, A) // New term?



     ==

        rho*(U&g)

      + combustion->Sh()

      + coalParcels.Sh(he)

      + limestoneParcels.Sh(he)

      + radiation->Sh(thermo)

      + fvOptions(rho, he)

    );

Thank you in advance for anyone who can provide some of the answers to my questions.

Hi DustExplosion!

Sorry, I am a beginner and not in a position to answer your questions.

By the way, in the PDRFoam problem I was talking about, the Allrun file gives these four commands-

runApplication blockMesh
runApplication changeDictionary
runApplication topoSet

runApplication PDRMesh -overwrite

I tried this. However, in the 'changeDictionary' step, I am facing the following error. Can you help me out?

balavignesh@bumblebee:~/OpenFOAM/balavignesh-4.1/run/flamePropagationWithObstacles$ changeDictionary
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : 4.1-1e03d68d4f4e
Exec : changeDictionary
Date : Nov 14 2016
Time : 14:52:26
Host : "bumblebee"
PID : 10607
Case : /home/balavignesh/OpenFOAM/balavignesh-4.1/run/flamePropagationWithObstacles
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Read dictionary changeDictionaryDict with replacements for dictionaries 1(dictionaryReplacement)
Reading polyMesh/boundary file to extract patch names
Loaded dictionary boundary with entries
6
(
outer
ground
blockedFaces
baffleWall
baffleCyclic_half0
baffleCyclic_half1
)

Extracted patch groups:
group cyclic with patches
2
(
baffleCyclic_half0
baffleCyclic_half1
)

group wall with patches
3
(
ground
blockedFaces
baffleWall
)

Replacing entries in dictionary dictionaryReplacement
Loading dictionary dictionaryReplacement

--> FOAM FATAL IO ERROR:
cannot find file

file: /home/balavignesh/OpenFOAM/balavignesh-4.1/run/flamePropagationWithObstacles/0/dictionaryReplacement at line 0.

From function regIOobject::readStream()
in file db/regIOobject/regIOobjectRead.C at line 72.

FOAM exiting

Thanks in advance!

Hi Balav,

I just copied the tutorial from OpenFoam-3.0.1/tutorial/combustion/PDRFoam/flamePropagationWithObstacles and ran ./Allcean and ./Allrun and it seems to be going fine.

On your latter question I believe that the 0 directory cannot be found. My guess is you are not running the "cp -r 0.org 0" command from the Allrun script.

On your original question, I am not sure exactly what the readMe is trying to say but it sounds like you are supposed to "hack" the mesh and break some connectivity to get the "obstacles" in place. Unfortunately, I am not sure where to do that. Maybe someone else who has used this tutorial can help? From the simulation results it looks like two box shaped obstacles are there, so I am not sure what else is needed.

Hi DustExplosion!

I was able to successfully run the tutorial. Sorry for the late reply. Thank you! :)

Hi DustExplosion:
What you have done is very impressive! But I don't think you can get a constant lewis number on a species basis by fix the conservation equations. OF do not use the transport data. http://cfd.direct/openfoam/energy-equation/ may help you understand the energy-equation, so according you works,the laminar flame speed calculated by OF is wrong? If you have some progress on this,please let me know. Many thanks!
Zhong

Hi!

I successfully managed to setup my case using XiFoam. However, after a few timesteps, I am getting the janaf (temperature out of bounds) warning. Eventually, after some timesteps, the simulation stops.

I searched through the posts and found some threads. But, I am not able to understand as to how I can resolve this error.

Can someone help in this regard?

Thanks!

Hi Zong,

Sorry for the late reply on this, I missed the email notification.

Thank you for the reference. It actually looks like it may provide some very good information that I was missing.

On "OF does not use transport data" - This is not really true, although there are approximations to transport data used. In the link you provided the "transport data" is hidden in q. A good description of some of the assumptions are given on page 34 of https://publikationen.bibliothek.kit...037446/2873782

Note that there are also other assumptions on how viscosity and thermal conductivity are calculated. One option to fix this is to couple to cantera (as is done in the thesis at the link). A second option is to create an openFoam library (see this thread for more info https://www.cfd-online.com/Forums/op...-update-5.html)

I would say that the laminar flame velocity calculated by openFoam is actually quite good considering the simplifications used (and how difficult doing it full from first principals is). I have decided to keep the code as it is for now, but am only studying up to an equivalence ratio of 0.8 (above this point the mass diffusion simplifications start to come in, but I think you could make an argument that the result is still qualitatively correct).

As an aside I also found several mechanisms (LuLaw30, DRM22, and DRM19) that can reproduce the GRI-Mech result that looks good up to 0.8, but runs faster.

Balav, Normally the jannaf errors happen when the boundarys are not set correctly or the simulation is unstable. Try reducing the timestep and checking your boundaries first. After that you may need to start removing reactions and stuff untill you can figure out what the issue it.

Another thought/question on this.

Does anyone know if changing the unity Schmidt number assumption to a unity lewis number assumption in Yenq, will have an negative effects?

I have modified the transport equation to use alpha (equal to conducitivty/Cp) instead of viscosity for my flame simulations, as it seems to produce better results for simple 1 and 2 step reaction models.

Code:

 fvScalarMatrix YiEqn

            (

                fvm::ddt(rho, Yi)

              + mvConvection->fvmDiv(phi, Yi)

              //- fvm::laplacian(turbulence->muEff(), Yi)

        - fvm::laplacian(turbulence->alphaEff(), Yi)

              ==

                coalParcels.SYi(i, Yi)

              + combustion->R(Yi)

              + fvOptions(rho, Yi)

            );

The part I am not sure about is on page 34/35 of this thesis: https://publikationen.bibliothek.kit...037446/2873782. The authors say that a unity Lewis number is used to derive the energy equation. However, the thermal diffusion term appears to not have any assumptions about Lewis number (it uses cond/Cp, which is the correct "fundamentally" derived coefficient.

In that case I am not sure why the authors bring up the Lewis number at all unless it is to say that the multicomponent part of the diffusion term can be taken away (but this happens whenever D is constant, irregardless of whether the Lewis number is unity).

It seems to me that the OF equations on use a constant Schmidt number assumption, but not constant Lewis number assumption. Does this makes sense?

Quote:

Originally Posted by DustExplosion (Post 623738)

Hello, Chris. I have looked your case setup files, and want to know how the file 'flameOutput.txt' is generated? Thanks in advance for your help.