# Difficulties with premixed laminar methane flame in openFoam

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

November 1, 2016, 19:38
Difficulties with premixed laminar methane flame in openFoam
#1
Member

Chris Cloney
Join Date: Jun 2016
Posts: 62
Rep Power: 5
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).
Compare_Ignition.jpg
Compare_Flame.jpg

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?

OF_Flame_Temperature.jpg
OF_Flame_Velocity.jpg
Attached Files
 OF_Flame.tar.gz (119.3 KB, 110 views)

 November 1, 2016, 19:59 #2 Member   Chris Cloney Join Date: Jun 2016 Location: Halifax, Canada Posts: 62 Rep Power: 5 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? Flame_Structure_Temperature.jpg Flame_Structure_Species.jpg Flame_Structure_Conductivity.jpg Flame_Structure_SpecificHeat.jpg Flame_Structure_Diffusivity.jpg 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 Thamali likes this.

November 4, 2016, 16:46
#3
Member

Chris Cloney
Join Date: Jun 2016
Posts: 62
Rep Power: 5
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.

MethaneFlame.jpg

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.
Attached Files
 methaneCase.tar.gz (115.8 KB, 96 views)

 November 5, 2016, 03:26 #4 New Member   Bala Join Date: Jul 2015 Posts: 18 Rep Power: 6 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!

 November 6, 2016, 13:33 #5 Member   Chris Cloney Join Date: Jun 2016 Location: Halifax, Canada Posts: 62 Rep Power: 5 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?

 November 9, 2016, 11:33 #6 New Member   Bala Join Date: Jul 2015 Posts: 18 Rep Power: 6 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!

November 9, 2016, 14:23
#7
Member

Chris Cloney
Join Date: Jun 2016
Posts: 62
Rep Power: 5
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:

OFFigure1.jpg

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:
heat_release.jpg
diffusivities.jpg

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:

species_GRIMECH.jpg

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!
Attached Files
 BFER_Mechanism_Phi1.0.tar.gz (87.1 KB, 60 views)

 November 13, 2016, 14:13 #8 Member   Chris Cloney Join Date: Jun 2016 Location: Halifax, Canada Posts: 62 Rep Power: 5 I realized that I made a mistake in my Sl vs phi plot from before. Attached is the updated version. OFFigure1.jpg 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?).

 November 13, 2016, 16:20 #9 Member   Chris Cloney Join Date: Jun 2016 Location: Halifax, Canada Posts: 62 Rep Power: 5 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 With the diffusion velocity and the correction velocity Assuming Ficks Law and using a constant lewis number approach for diffusion coefficient Current OpenFoam Treatment Currently in coalChemistryFoam all lewis numbers are constant and equal to 1. Therefore: and The correction velocity goes to zero and the governing equation becomes: Variable Lewis Number Treatment With the variable lewis number the diffusion velocity and correction velocity are: and the governing equation is 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 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 is: 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 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?: which would mean that 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 but I cannot see how we know the summation of the product of enthalphy and the gradient goes to zero 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. Tobi, arvindpj and Eldor like this.

 November 15, 2016, 09:53 #11 Member   Chris Cloney Join Date: Jun 2016 Location: Halifax, Canada Posts: 62 Rep Power: 5 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.

 November 21, 2016, 03:29 #12 New Member   Bala Join Date: Jul 2015 Posts: 18 Rep Power: 6 Hi DustExplosion! I was able to successfully run the tutorial. Sorry for the late reply. Thank you!

 November 22, 2016, 04:58 #13 New Member   Join Date: Jan 2016 Posts: 15 Rep Power: 5 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

 November 23, 2016, 02:38 #14 New Member   Bala Join Date: Jul 2015 Posts: 18 Rep Power: 6 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!

 December 6, 2016, 13:03 #15 Member   Chris Cloney Join Date: Jun 2016 Location: Halifax, Canada Posts: 62 Rep Power: 5 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. balav94 likes this.

 December 6, 2016, 13:23 #16 Member   Chris Cloney Join Date: Jun 2016 Location: Halifax, Canada Posts: 62 Rep Power: 5 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? arvindpj likes this.

August 9, 2019, 05:31
#17
New Member

Yujuan Luo
Join Date: Jul 2019
Posts: 3
Rep Power: 2
Quote:
 Originally Posted by DustExplosion 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

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.

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post SBusch OpenFOAM 22 December 26, 2016 14:24 AleDR OpenFOAM Running, Solving & CFD 22 October 12, 2016 03:12 cfd.direct OpenFOAM Announcements from Other Sources 0 September 14, 2016 03:19 fidy FLUENT 0 June 28, 2014 23:45 Cristian_Alex CFX 7 November 27, 2012 10:22

All times are GMT -4. The time now is 19:49.

 Contact Us - CFD Online - Privacy Statement - Top