CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

What do you all do to stabilize reactingFoam (or in general)?

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   December 6, 2016, 11:48
Default What do you all do to stabilize reactingFoam (or in general)?
  #1
Member
 
Anonymouse
Join Date: Dec 2015
Posts: 73
Rep Power: 3
KarenRei is on a distinguished road
I've been working on this problem for the past several weeks. I have a simulation modeling the combustion of oxygen, paraffin and aluminum with several intermediary species. It's built atop a stable rhoSonicFoam run, dealing with it without combustion before switching over to reactingFoam.

It runs for a fair while, but almost inevitably, eventually "something happens" that leads to a blowup. In general, it appears to usually be along the lines of turbulence somewhere causing a localized patch of heating, which increases the reaction rate, increasing heating faster, etc such that temperatures skyrocket or go negative and things blow up. Yet it also happens (faster, actually!) if I set the turbulence model to laminar and/or the combustion turbulence to laminar.

I can generally prevent these blowups by dropping the timestep way down. But then it runs at a crawl. It's like 99.999999% of all calculations are fine, but that last 0.000001% blows things up, and to prevent that I have to tremendously slow things down for the 99.999999%. It would be lovely if chemistry could limit the timestep the same way as the courant number does... but I see nothing in either the docs nor the code to allow that.

* I make use of a temperature limit in fvOptions, but it doesn't seem to do anything useful.
* All of my fvSchemes are first-order that can be, most commonly Gauss linearUpwind grad(U).
* I use no relaxation parameter in fvSolution because they don't seem to help much, only delaying the inevitable
* My experimentation is mixed as to whether varying tolerances in fvSolution helps any or make things worse. And if they're high enough then strangeness crops up in the run, so for now I'm just keeping them low.

What do you do in such a situation to try to prevent these blowups? I've even considered a cheap hack, running reactingFoam via a looping external script that automatically lowers the timestep (then raises it after a period) if things blow up. But often when I restart reactingFoam the existing cell contents change for some bizarre reason - for example, all of the fuel might disappear and have to be repopulated, which of course screws up the burn, introduces more turbulence and just makes the problem worse.
KarenRei is offline   Reply With Quote

Old   December 6, 2016, 17:56
Default ODE solver choice
  #2
Member
 
Arvind Jay
Join Date: Sep 2012
Posts: 48
Rep Power: 6
arvindpj is on a distinguished road
What ODE solver are you using? Your reactions could be stiff.
arvindpj is offline   Reply With Quote

Old   December 6, 2016, 19:07
Default
  #3
Member
 
Anonymouse
Join Date: Dec 2015
Posts: 73
Rep Power: 3
KarenRei is on a distinguished road
Quote:
Originally Posted by arvindpj View Post
What ODE solver are you using? Your reactions could be stiff.
I'm using EulerImplicit, with parameters:

cTauChem 1.0;
equilibriumRateLimiter off;
KarenRei is offline   Reply With Quote

Old   December 7, 2016, 10:43
Default
  #4
Member
 
Arvind Jay
Join Date: Sep 2012
Posts: 48
Rep Power: 6
arvindpj is on a distinguished road
Quote:
Originally Posted by KarenRei View Post
I'm using EulerImplicit, with parameters:

cTauChem 1.0;
equilibriumRateLimiter off;
EulerImplict is 1st order method. Try using other higher-order ODE solvers. All depends on the stiffness of the reactions. Since you are using many reactions with different time scales, I am assuming your system of ODE is stiff. Try:

odeCoeffs
{
solver seulex;
absTol 1e-12;
relTol 1e-1;
}
arvindpj is offline   Reply With Quote

Old   December 7, 2016, 17:49
Default
  #5
Member
 
Anonymouse
Join Date: Dec 2015
Posts: 73
Rep Power: 3
KarenRei is on a distinguished road
I got an error that doesn't show up when I google it:


Quote:
--> FOAM FATAL ERROR:
Integration steps greater than maximum 10000xStart = 0, xEnd = 4.173913e-06, x = 3.1435182e-06, dxDid = 1.7084929e-13

From function virtual void Foam::ODESolver::solve(Foam::scalar, Foam::scalar, Foam::scalarField&, Foam::scalar&) const
in file ODESolvers/ODESolver/ODESolver.C at line 191.

FOAM exiting
Having not used this solver before, I have no clue what that means / how to fix it. I have no parameter called integration steps (or "maxSteps" as it appears to be in the code), nor do I know if it's even reasonable that it's that high.
KarenRei is offline   Reply With Quote

Old   December 8, 2016, 11:12
Default
  #6
Member
 
Arvind Jay
Join Date: Sep 2012
Posts: 48
Rep Power: 6
arvindpj is on a distinguished road
Quote:
Originally Posted by KarenRei View Post
I got an error that doesn't show up when I google it:




Having not used this solver before, I have no clue what that means / how to fix it. I have no parameter called integration steps (or "maxSteps" as it appears to be in the code), nor do I know if it's even reasonable that it's that high.
This means the ODE Solver is not converging.
Could you:
  • try reducing Co between 0.2-0.4
  • try different ODE solver
  • try reduced mechanisms.


As an exercise, you could try running chemFoam for various mole fractions to identify the bottleneck and pick the best performing ODE solver.
arvindpj is offline   Reply With Quote

Old   December 8, 2016, 16:49
Default
  #7
Member
 
Anonymouse
Join Date: Dec 2015
Posts: 73
Rep Power: 3
KarenRei is on a distinguished road
Quote:
Originally Posted by arvindpj View Post
This means the ODE Solver is not converging.
Could you:
  • try reducing Co between 0.2-0.4
  • try different ODE solver
  • try reduced mechanisms.


As an exercise, you could try running chemFoam for various mole fractions to identify the bottleneck and pick the best performing ODE solver.
I never run with a courant number greater than 0.2; my maxCo is 0.2.
I tried all of the ODEs in the list. The only ones that don't immediately crash are Euler and EulerSI. They eventually crash, like using EulerImplicit..

I am using reduced mechanisms, at least for part of the reactions.

Some of the reaction rates need to be determined empirically (my objective in this project involves determining them, but before that I need to set up the combustion model). Those reactions are the partial oxidation of the paraffin ("C30H62^0.25 + 15O2^1.5 = 30CO + 62H") - using A 3e15, beta 0, and Ta 15000.0 to start - and a variety of aluminum reactions. The exponents are from a paper, but A and Ta are initial guesses.

The aluminum is modeled as Al10000 burning to Al20000O30000, aka a high molecular weight so that neither the metal nor the oxide contribute relevantly to the gas pressure (emulating the combustion of particles to other particles, without having to actually model them). The aluminum reactions are all based on variants of "0.0002Al10000 + 1.5O2^0.2 = 0.0001Al20000O30000" (using A 0.3, beta 0.2, Ta 4000 for now), with the various alternatives for the O2 having different values for A. beta=0.2 and the relative burn rates between different aluminum reactions are from a paper, but the A=0.3 and Ta=4000 for the baseline case are initial guesses.

All other reactions are either from papers ("2AlOH + O^0.2 = 0.0001Al20000O30000 + 2H" A 5.0e11 beta 0 Ta 50000) or from OpenFOAM tutorial files, involving the various reactions between O2, O, N2, CO, H2O, O, CO2, OH, H, and H2. The total number of reactions is 18:

reaction "C30H62^0.25 + 15O2^1.5 = 30CO + 62H"
reaction "CO + O = CO2";
reaction "CO + OH = CO2 + H";
reaction "CO + O2 = CO2 + O";
reaction "O + H2 = H + OH";
reaction "OH + H2 = H + H2O";
reaction "H2 = H + H";
reaction "O + OH = O2 + H";
reaction "OH + OH = O + H2O";
reaction "H + OH = H2O";
reaction "H + O = OH";
reaction "O + O = O2";
reaction "0.0001Al10000 + H2O^0.2 = AlOH + 2H";
reaction "0.0001Al10000 + OH^0.2 = AlOH";
reaction "0.0002Al10000 + 1.5O2^0.2 = 0.0001Al20000O30000";
reaction "0.0002Al10000 + 3CO2^0.2 = 0.0001Al20000O30000 + 3CO";
reaction "0.0002Al10000 + 3O^0.2 = 0.0001Al20000O30000";
reaction "2AlOH + O = 0.0001Al20000O30000 + 2H";

Last edited by KarenRei; December 9, 2016 at 08:53.
KarenRei is offline   Reply With Quote

Old   December 11, 2016, 15:34
Default
  #8
Member
 
Anonymouse
Join Date: Dec 2015
Posts: 73
Rep Power: 3
KarenRei is on a distinguished road
No further thoughts on this?
KarenRei is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
error with reactingFoam BakedAlmonds OpenFOAM Running, Solving & CFD 4 June 22, 2016 02:21
How to switch off combustion and reaction in reactingFoam shenzhou1987 OpenFOAM Running, Solving & CFD 12 June 20, 2016 07:23
calculate flame speed using reactingFoam IColin OpenFOAM Running, Solving & CFD 0 February 4, 2014 16:14
reactingFoam wedge handling wrong U dhondupant OpenFOAM Bugs 1 December 9, 2010 08:34
reactingFoam - turbulent reacting flow hamburgFoam OpenFOAM 0 December 7, 2009 13:57


All times are GMT -4. The time now is 10:17.