[cryabroad]

Hello Foamers,

As far as I understand, the two solvers reactingFoam and rhoReactingFoam differ ONLY in the way how the density is calculated. One is calculated from psi*p and one is directly accessed from the thermodynamics library (By the way is this true?). They are both pressure-based solvers, because they both solve a pressure equation to satisfy the continuity equation.

However, their pressure equations are different. I also notice that starting from OpenFOAM-5.x, the pressure equation in rhoReactingFoam is the same as the one in rhoPimpleFoam. But before that, the pressure equation in reactingFoam is the same as the one in rhoPimpleFoam. I basically have two questions: 1. Why the pressure equation is different in reactingFoam and rhoReactingFoam? Aren't they solving the same pressure equation? 2. How is rhoReactingFoam (or reactingFoam) connected to rhoPimpleFoam? I would think that with reaction turned off, these three should be the same thing.

Thanks in advance,
Ruiyan

[HPE]

Would you mind to give snippets of the differences, otherwise quite difficult to follow?

[cryabroad]

Hi HPE, thank you for your interest! I should've done that to make myself more clear. Take OpenFOAM-6 for example, the relevant code are as follows.

From the pEqn.H in rhoReactingFoam (and rhoPimpleFoam)

if (!pimple.simpleRho())
{
rho = thermo.rho();
}

// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
const volScalarField psip0(psi*p);

code that are not relevant

else
{
fvScalarMatrix pDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
==
fvOptions(psi, p, rho.name())
);

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));

pEqn.solve(mesh.solver(p.select(pimple.finalInnerI ter())));

if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}

// Thermodynamic density update
thermo.correctRho(psi*p - psip0);

and from the pEqn.H in reactingFoam

rho = thermo.rho();

code that are not relevant

else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi))
)
);

MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);

pEqn.solve(mesh.solver(p.select(pimple.finalInnerI ter())));

if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}

If I understand it correctly, the difference is the additional psi*correction(fvm::ddt(p)) term used in rhoReactingFoam. However, this seems to be coming from the fact that rho = psi * p (which is the basis of all psi-based thermodynamics in FOAM), and if we take psi as a constant within the iteration, then we have rho_new-rho_old = psi*(p_new-p_old).

The pressure equation in reactingFoam seems natural and easy to understand. It is just a Poisson's equation where rho does not appear explicitly because of the fact that rho = psi * p. However, the pressure euqation in rhoReactingFoam involves a fvc::ddt(rho) term, which is explicit, and the correction term I mentioned above. Why can't we use the same Poisson's equation as that in reactingFoam? I feel like I'm missing some obvious things here and it could be very obvious to other people.

And one more thing, if I run a same case with these two solvers, it should give me the exact same results right? In other words, in terms of the ability of modelling the physics, these two should be identical, or maybe one is better than the other under some extreme conditions? Sorry for my lengthy statements/questions above, as I'm getting picky about every detail since I started using FOAM, dare I say it's kind of addictive.

[cryabroad]

Anyone? Still waiting for some answers/discussions here.

[holmej]

Four years have passed, can anyone have a clearer understanding of this issue and share it? I would appreciate it.

[dlahaye]

One element in the discussion might be the fact that reactingFoam is able to solve flows that are not incompressible.

Understanding the pressure equation in reactingFoam might thus involve two steps:

1/ step 1/2: understand how the pressure equation changes when passing from incompressible constant density flow (say simpleFoam) to still incompressible, but now variable density flow (say buoyantSimpleFoam). The energy equation and the thermodynamics (equation of state) needs to be taken into account.

2/ step 2/2: understand how the pressure equation changes when passing from incompressible variable density flow to no longer incompressible flows. The book of e.g. Darwish-Mangani-Moukalled gives a good theoretical basis. Having settled that basis, one can discuss what happens in the implementation.

[holmej]

Originally Posted by dlahaye

One element in the discussion might be the fact that reactingFoam is able to solve flows that are not incompressible.

Understanding the pressure equation in reactingFoam might thus involve two steps:

1/ step 1/2: understand how the pressure equation changes when passing from incompressible constant density flow (say simpleFoam) to still incompressible, but now variable density flow (say buoyantSimpleFoam). The energy equation and the thermodynamics (equation of state) needs to be taken into account.

2/ step 2/2: understand how the pressure equation changes when passing from incompressible variable density flow to no longer incompressible flows. The book of e.g. Darwish-Mangani-Moukalled gives a good theoretical basis. Having settled that basis, one can discuss what happens in the implementation.

Thanks for your reply. But different version of the OF have distinguished equations about the reactingFoam and rhoReactingFoam, so I need more time to figure out which solver is suitable for compressible reaction flow.Thank you again.

[dlahaye]

In most recent versions of OF, reactingFoam and rhoReactingFoam have merged into reactingFoam.

If flow is compressible, make sure to set transonic on in system/controlDict.

[holmej]

Originally Posted by dlahaye

In most recent versions of OF, reactingFoam and rhoReactingFoam have merged into reactingFoam.

If flow is compressible, make sure to set transonic on in system/controlDict.

Yes, you are right. In OF V9, these two solvers had been merged. The "transonic" you mentioned seems to be effective in simple flow without combustion. I'm not familiar with OF. If I am wrong, please point it out.

[dlahaye]

What do you mean by "effective" in this context?

[holmej]

Originally Posted by dlahaye

What do you mean by "effective" in this context?

I'm sorry for unclear expression. I mean "transonic" switch only appear in the tutorials without combustion, like rhoCentralFoam. Tutorials of reactingFoam have no option about "transonic".

[dlahaye]

What follows could merely be my points of view. The transonic options belongs to a solver like reactingFoam, more than to a tutorial case.

I do know at Mach number e.g. the chokedNozzle or membrane tutorial case are set. One could merely change the inlet velocity, vary the Mach number and verify how reactingFoam behaves with transonic set either on or off.

Not sure whether this addresses your concern.

[holmej]

Originally Posted by dlahaye

What follows could merely be my points of view. The transonic options belongs to a solver like reactingFoam, more than to a tutorial case.

I do know at Mach number e.g. the chokedNozzle or membrane tutorial case are set. One could merely change the inlet velocity, vary the Mach number and verify how reactingFoam behaves with transonic set either on or off.

Not sure whether this addresses your concern.

Thank you for providing the idea. It is a way to verify whether the reactingFoam acts like compressible solvers with transonic or supersonic conditions.

[dlahaye]

Pls beware that changing the combustion/reactingFoam/RAS tutorial cases from incompressible to transonic (by increasing e.g. the mass flow rate of the fuel) might require changing the boundary conditions for the pressure. We have found the wave transmissive boundary conditions at the outlet to work well. Good luck.

[holmej]

Originally Posted by dlahaye

Pls beware that changing the combustion/reactingFoam/RAS tutorial cases from incompressible to transonic (by increasing e.g. the mass flow rate of the fuel) might require changing the boundary conditions for the pressure. We have found the wave transmissive boundary conditions at the outlet to work well. Good luck.

You are a kind man. Thank you for your advice!