[Peoban2]

Hi

I have tried to modify reactingFoam in order to make it a steady state solver, so I took as base file rhoSimpleFoam and added the YEqn term on the main .C file of the solver and also the corresponding terms on theEEqn. The problem is that when I run it, the Temperature increase as well as the CO2 production are much lower than the results obtained using Fluent, which I was using in order to validate the solver, or even the same case using reactingFoam after a long time.

Actually I don´t know where the error is, whether in the solver itself, in the relaxation parameters or tolerances and numerical solvers defined in fvSolution. I´ve tryed to check as much as I can, but I have no clue.

Here I attach the solver, I couldn´t attach the case I´m modelling since it´s a little heavy, but I would have no problem also sending it.

I´m new in the whole CFD world and this is my first customized solver so would really appreciate any help I could get.

[wyldckat]

Greetings Pedro,

If the case is too big, then it's best that you create a much smaller test case that you can share, because it will also make it easier for anyone to help you.

Beyond this, I'm not so certain that it's possible to create a steady-state reaction solver, since it might not be able to release the heat fast enough, therefore the temperature would continue to rise continually. Keep in mind that steady-state simulations are almost an averaging of the flow over a long time frame.

Here are a few more ideas/hints:

• PIMPLE is a PISO-SIMPLE hybrid. And yes, SIMPLE can also be used in transient simulations. Here are a few wiki pages on these topics:
  • http://openfoamwiki.net/index.php/Ma...ransientSimple
  • http://openfoamwiki.net/index.php/Op...hm_in_OpenFOAM
    • http://openfoamwiki.net/index.php/PISO
  • http://openfoamwiki.net/index.php/Op...hm_in_OpenFOAM
• Have a look at the following tutorials, specifically at the "system/fvSchemes" file:
  

  tutorials/incompressible/simpleFoam/pitzDaily
  tutorials/incompressible/icoFoam/cavity
  tutorials/combustion/reactingFoam/ras/counterFlowFlame2D

  The main difference is that many of the transient tutorial cases use this:
  

  ddtSchemes
  {
      default         Euler;
  }

  While steady-state solvers use this:
  

  ddtSchemes
  {
      default         steadyState;
  }

  More information is given in the user guide: http://www.openfoam.org/docs/user/fvSchemes.php

Best regards,
Bruno

[Peoban2]

Hi Bruno

Thank you for replying, as a simpler case here I attach the modified counterFlowFlame2D case in order to be solved with the steady state solver that I previously uploaded. When comparing the results obtained with this solver when it converges (according to the tolerances that I stated that I´m also not sure if are ok) with the ones obtained with the transient solver after 10s (when I assumed that it was long enough as to have reached steady state) the differences in the obtained temperatures and CO2 concentrations are huge.

I have found this steady state solver with chemical reactions
https://github.com/Unofficial-Extend...rnateChemistry
It´s just that its in an older OF version which I haven´t installed so I couldn´t test it, however by checking the code, it seemed to me that it followed pretty much the same logic as the modifications I tried to implement, except that on the Yeqn they include a some controling loops in order to stabilize the solution.

thanks

Pedro

[wyldckat]

Greetings Pedro,

Sorry for the reaaaallly late reply, but only today did I finally managed to look into this.

Technically I'm not an expert on this topic, but as I stated in the previous post, the issue seems to be related to problems on how the heat is injected and extracted from the domain. Transient simulations are time dependant and steady-state simulations usually are not. And heat transfer is usually very time dependant, which is very likely why your simulations were either crashing or increasing too much in temperature.

And from what I could figure out, you used OpenFOAM 2.2 (not sure which minor version).

In file "EEqn.H", you have this line:

        reaction->Sh()  // THIS IS THE TERM THAT BLOWS IT OFF!!!!!

Indeed, this term has the following units:

dimEnergy/dimVolume/dimTime

which means that it's strongly dependant on time. My deduction is that you'll need to somehow neutralize the time dependence for this and other terms that depend directly on time.

If we look at the source code for "AlternateChemistry" repository you pointed out, you'll notice a few important details:

1. The question on the transient "h" equation is this:
   

       solve
       (
           fvm::ddt(rho, h)
         + mvConvection->fvmDiv(phi, h)
         - fvm::laplacian(turbulence->alphaEff(), h)
        ==
           DpDt
       );

2. While on the steady-state, the respective piece of code is this:
   

   fvScalarMatrix hEqn
       (
           fvm::div(phi, h)
         - fvm::Sp(fvc::div(phi), h)
         - fvm::laplacian(turbulence->alphaEff(), h)
        ==
           fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,p)"))
         - p*fvc::div(phi/fvc::interpolate(rho))
       );
   
       hEqn.relax();
   
       if(fixedCells.size()>0) {
           Info << "Fixing enthalpy in " << fixedCells.size() << " cells" << endl;
   
           scalarField hFix=dynamic_cast<basicThermo&>(thermo()).h(fixedTemp,fixedCells);
           hEqn.setValues(fixedCells,hFix);
       }
   
       eqnResidual = hEqn.solve().initialResidual();

3. There are 2 major differences between the two:
   1. This:
      

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

      was replaced by this:
      

              fvm::div(phi, h)
            - fvm::Sp(fvc::div(phi), h)

   2. This:
      

      DpDt

      was replaced by this:
      

              fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p, "div(U,p)"))
            - p*fvc::div(phi/fvc::interpolate(rho))

4. This essentially points out that the convection term has to be accounted a bit differently, as well as the pressure based changes (the second difference).

But then I took a look at the source code differences between buoyantSimpleFoam and rhoSimpleFoam, and the only difference between the two in the "EEqn.H" file is the source term "Sh" for the radiation, which has the same exact units as the other "Sh" term:

    fvScalarMatrix EEqn
    (
        fvm::div(phi, he)
      + (
            he.name() == "e"
          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
     ==
        radiation->Sh(thermo)
      + fvOptions(rho, he)
    );

... uhm, sorry, now I'm lost. I tried using OpenFOAM 2.2.x, but the solver is crashing regardless of what is done.

Weirder still is that I've taken a look at the source code for reactingParcelFoam and simpleReactingParcelFoam, which are essentially the same changes you've tried to do.
So I took another glance at the tutorials for these two solvers and found that there the chemistry solver cannot use the ODE option, in the file "constant/chemistryProperties".
By using the following in your test case, the simulation worked fine:

chemistryType
{
    //chemistrySolver   ode;
    chemistrySolver   noChemistrySolver;
    chemistryThermo   psi;
}

As for the results I got, they seem consistent albeit a bit lower in the maximum temperature. This is somewhat natural to occur, because not solving the chemistry part can lead to inaccuracies...

I then tried:

chemistryType
{
    //chemistrySolver   ode;
    chemistrySolver   sequential; //noChemistrySolver;
    chemistryThermo   psi;
}

chemistry           on;

initialChemicalTimeStep 1e-07;

sequentialCoeffs
{
    cTauChem        0.001;
    equilibriumRateLimiter on;
}

And it gave much better results!

In addition, I removed the residual stopping criteria and at around iteration 4410 it was as good as stable. I then compared with the run in the original tutorial "combustion/reactingFoam/ras/counterFlowFlame2D" and it gave closer results.

Attached is this last test case I used, in the package "steady_counterFlowFlame2D_v2.tar.gz".
The values I got:

• transient:
  

  min/max(T) = 788.691, 828.766

• steady-state:
  

  min/max(T) = 792.176, 829.846

These differences scale are normal to occur, since the resolution between the two methods are not identical. In addition, a mesh study should be done in order to access what is the exact origin of the differences.

By the way, I used your solver without additional modifications.

Best regards,
Bruno

[flowAlways]

Thank you very much Pedro and Bruno! I

Now have some good hopes with my masters thesis!!. The current simulation time for our furnace application using reactingFoam is in days. Hope we will bring it to hours using the steady state version(hope it works!!).

Earlier I was planning to accelerate the simulation using the idea of more outer correctors and under-relaxation with PIMPLE. But it is good to know that SIMPLE can be used too!

I will soon let you guys know how my results with the new solver compare to the existing reactingFoam results(after a long time).

[flowAlways]

Hello Bruno,

If you could find time please can you have a look. I tested the solver simpleReactingFoam in OF2.2.2 with the following problem.
1. Non-premixed turbulent gaseous combustion in a furnace.
2. Swirling Flow at Inlet.
3. Global 1 step irreversible reaction mechanism for CH4.

I used all the settings with the counterFlowFlame2D tutorial(except 0 and polyMesh) provided by you. The simulation in the beginning throws janafThermo error but later(after many iterations) the simulation proceeds without error. However the result compared to reactingFoam and Fluent is completely wrong. The max temperatures falls to 1200K where it should be ~2000K .

I dont know what am I doing wrong or the solver is not completely correct. ?? The mesh is very refined and therefore the simulation takes many hours(5-6) for convergence. I increased the max steps to 10000.

SIMPLE: no convergence criteria found. Calculations will run for 4410 steps.

Creating reaction model

Selecting combustion model PaSR<psiChemistryCombustion>
Selecting chemistry type
{
    chemistrySolver sequential;
    chemistryThermo psi;
}

Selecting thermodynamics package
{
    type            hePsiThermo;
    mixture         reactingMixture;
    transport       sutherland;
    thermo          janaf;
    energy          sensibleEnthalpy;
    equationOfState perfectGas;
    specie          specie;
}

Selecting chemistryReader foamChemistryReader
chemistryModel: Number of species = 5 and reactions = 1
Reading field U

Reading/calculating face flux field phi

Creating turbulence model.

Selecting turbulence model type RASModel
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    C3              -0.33;
    sigmak          1;
    sigmaEps        1.3;
    Prt             1;
}

Creating field kinetic energy K

No finite volume options present


Reading g
Courant Number mean: 109.439 max: 21454.1

Starting time loop

Time = 1

smoothSolver:  Solving for Ux, Initial residual = 0.00530601, Final residual = 0.000446973, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 0.0424394, Final residual = 0.00302406, No Iterations 2
smoothSolver:  Solving for Uz, Initial residual = 0.00906213, Final residual = 0.000165917, No Iterations 4
DILUPBiCG:  Solving for O2, Initial residual = 1, Final residual = 6.42829e+09, No Iterations 1001
DILUPBiCG:  Solving for H2O, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for CH4, Initial residual = 1, Final residual = 4.10261e+10, No Iterations 1001
DILUPBiCG:  Solving for CO2, Initial residual = 0, Final residual = 0, No Iterations 0
DILUPBiCG:  Solving for h, Initial residual = 1, Final residual = 0.0107451, No Iterations 1
min/max(T) = 300, 2102.62
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.0378342, No Iterations 16
time step continuity errors : sum local = 0.124907, global = -0.123693, cumulative = -0.123693
rho max/min : 1.14168 0.174843
smoothSolver:  Solving for epsilon, Initial residual = 0.190803, Final residual = 0.0150082, No Iterations 2
smoothSolver:  Solving for k, Initial residual = 1, Final residual = 0.040656, No Iterations 4
ExecutionTime = 1.97 s  ClockTime = 2 s

Time = 2

smoothSolver:  Solving for Ux, Initial residual = 0.0950581, Final residual = 0.00133338, No Iterations 4
smoothSolver:  Solving for Uy, Initial residual = 0.156577, Final residual = 0.0023826, No Iterations 4
smoothSolver:  Solving for Uz, Initial residual = 0.0181092, Final residual = 0.000291685, No Iterations 4
--> FOAM Warning :
    From function janafThermo<EquationOfState>::limit(const scalar T) const
    in file /opt/apps/openfoam-2.2.2/OpenFOAM-2.2.2/src/thermophysicalModels/specie/lnInclude/janafThermoI.H at line 108
    attempt to use janafThermo<EquationOfState> out of temperature range 200 -> 5000;  T = 5748.93
--> FOAM Warning :
    From function janafThermo<EquationOfState>::limit(const scalar T) const
    in file /opt/apps/openfoam-2.2.2/OpenFOAM-2.2.2/src/thermophysicalModels/specie/lnInclude/janafThermoI.H at line 108
    attempt to use janafThermo<EquationOfState> out of temperature range 200 -> 5000;  T = 6743.5

[TBO]

Interesting discussion how to make reactingFoam steady state, looking forward to the results. Might be worth implementing it in OF 2.3 instead of 2.2. Some new solvers have been added in OF 2.3 which can speed up the chemistry calculations significantly. For some reason, OF 2.2 is relatively slow with chemistry (compared to OF 2.1 and OF 2.3). Running similations with chemFoam (which basically uses the same chemistry models and solvers as reactingFoam, but only for a single cell) are much faster in versions other than 2.2.

[wyldckat]

Greetings to all!

@Ali: Sorry for the very late reply, but only today did I finally manage to give at least a quick look into this. I guess that since so much time as passed, you've already moved on and/or solver your problem.

Either way, as pointed out by TBO in another thread: http://www.cfd-online.com/Forums/ope...tml#post542001 post #108 - the LTSReactingFoam solver would be a preferable approach. Because steady state heat transfer and chemical reactions, when done incorrectly, can result in having the whole fluid domain populated by the maximum temperature, due to averaging effects (not the actual technical term, but I can't remember right now which one is the correct term).

I'll answer your other question on the thread I mentioned.

Best regards,
Bruno

[TinaB]

Dear all,

I have been reading this (and the related) threads with great interest, since I also want to have a steady state solver for combustion. Following the suggestions posted above, I made some experiments with LTSReactingParcelFoam, based on the counterFlowFlame2D tutorial. At first, it looked really promising, but I found a problem. Maybe you can shed some light.

This is what I did:

Test 1:
No modifications, I just ran the tutorial. After 2000 iterations, the maximal temperature was 2432 K.

Test 2:
I changed the maximal Courant number from 1 to 5 in fvSolution:

    maxCo           5;

Resulting maximal temperature did not change significantly. The same happened with maxCO=10.

Test 3:
I used a Courant number smaller than one.

    maxCo           0.5;

The maximal temperature changed by 10 K to 2442 K.

Test 4:

    maxCo           0.1;

Maximal temperature: 2491 K.

Can anyone explain this behaviour? The Courant number gives - as far as I understand it - a maximum limit for the "time-step". Why does the temperature change if I make the time step smaller?

I suspect that there are problems in the solver (or the chemistryModel, or both) concerning the modeling of the species and energy source terms for combustion. The source terms seem to depend on the time step, even though the solver is steady-state.

A related question: Did anyone ever successfully use simpleReactingParcelFoam for a combustion case?

[ssss]

Because if you reduce the time step you reduce the time discretization errors and thus your solution is less diffusive