[universez]

Hi,
I try to study rhoSimpleFoam. In the pEqn.H of rhoSimpleFoam, there is a line "if (transonic)" see below for the file. I understand the function of the that is using incompressible SIMPLE method for low mach number cases and adding convertive term for compressible flow. (by the way, the compressible part in the rhoSimpleFoam of OpenFoam-dev1.5 is different).

I did output the the value transonic by adding
"Info<<"transonic = "<<transonic<<endl;"
in the rhoSimpleFoam.

the output of the is always "0" (see below)

My question are :

1, how to control transonic? it is specified at somewhere or it is defined by mach number?

2, if transonic = 0, it is always solving an incompressible problem but not a compressible problem.

Thank
yu

==============================================

Time = 172
smoothSolver: Solving for Ux, Initial residual = 0.130169, Final residual = 6.88999e-05, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.120254, Final residual = 9.18545e-05, No Iterations 2
DILUPBiCG: Solving for h, Initial residual = 0.0154733, Final residual = 4.72774e-05, No Iterations 1
transonic = 0
GAMG: Solving for p, Initial residual = 0.277657, Final residual = 0.00683687, No Iterations 6
time step continuity errors : sum local = 5.41896, global = 0.830172, cumulative = 3.09577
rho max/min : 1.79102 1.26744
ExecutionTime = 53.24 s ClockTime = 68 s

===============================================





===========================================
below is the modified pEqu file in the folder of rhoSimpleFoam
================================================
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
volScalarField rUA = 1.0/UEqn().A();
U = rUA*UEqn().H();
UEqn.clear();
bool closedVolume = false;
Info<<"transonic = "<<transonic<<endl; // output transonic
if (transonic)
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)*(fvc::interpolate(U) & mesh.Sf())
);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::div(phid, p)
- fvm::laplacian(rho*rUA, p)
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax(mesh.relaxationFactor("pEqn"));
pEqn.setReference(pRefCell, pRefValue);
// retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pEqn.solve();
}
if (nonOrth == nNonOrthCorr)
{
phi == pEqn.flux();
}
}
}
else
{
phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
closedVolume = adjustPhi(phi, U, p);
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rho*rUA, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
// Retain the residual from the first iteration
if (nonOrth == 0)
{
eqnResidual = pEqn.solve().initialResidual();
maxResidual = max(eqnResidual, maxResidual);
}
else
{
pEqn.solve();
}
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
}

#include "incompressible/continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/fvc::domainIntegrate(psi);
}
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;

Yu