cryabroad |
December 29, 2019 21:06 |
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)
Quote:
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
Quote:
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:D.
|