
[Sponsors] 
Difficulties with premixed laminar methane flame in openFoam 

LinkBack  Thread Tools  Search this Thread  Display Modes 
November 1, 2016, 19:38 
Difficulties with premixed laminar methane flame in openFoam

#1 
Member
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
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 

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 1020% 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 

November 4, 2016, 16:46 

#3 
Member
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
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. 

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 zerosize 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 23 days and get back to you. Thanks! 

November 9, 2016, 14:23 

#7 
Member
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
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 velocitycontinuity 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 GRIMECH 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 GRIMECH 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! 

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 grimech 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.cfdonline.com/Forums/openfoam/103227multispeciesmasstransportlibraryupdate.html) which includes source code. Unfortunately, the source code did not want to download from cfdonline 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; } 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) ); 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) ); 

November 14, 2016, 04:34 

#10 
New Member
Bala
Join Date: Jul 2015
Posts: 18
Rep Power: 6 
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/balavignesh4.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.11e03d68d4f4e Exec : changeDictionary Date : Nov 14 2016 Time : 14:52:26 Host : "bumblebee" PID : 10607 Case : /home/balavignesh/OpenFOAM/balavignesh4.1/run/flamePropagationWithObstacles nProcs : 1 sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE). fileModificationChecking : Monitoring runtime modified files using timeStampMaster allowSystemOperations : Allowing usersupplied 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/balavignesh4.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! 

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 OpenFoam3.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/energyequation/ may help you understand the energyequation, 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.cfdonline.com/Forums/op...update5.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 GRIMech 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. 

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) ); 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? 

August 9, 2019, 05:31 

#17  
New Member
Yujuan Luo
Join Date: Jul 2019
Posts: 3
Rep Power: 2 
Quote:
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 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
OpenFOAM v3.0+ ??  SBusch  OpenFOAM  22  December 26, 2016 14:24 
Combustion modelling in OpenFOAM  Difficulties  AleDR  OpenFOAM Running, Solving & CFD  22  October 12, 2016 03:12 
OpenFOAM Training, London, Chicago, Munich, Houston 20162017  cfd.direct  OpenFOAM Announcements from Other Sources  0  September 14, 2016 03:19 
Modelling premixed laminar flame using Fluent  fidy  FLUENT  0  June 28, 2014 23:45 
Flame Lenght Overpredicted  Methane / Air Combustion  Cristian_Alex  CFX  7  November 27, 2012 10:22 