rhoPorousMRFSimpleFoam: Changed to rhoThermo

Also renamed addEnthalpySource -> addEnergySource
This commit is contained in:
Henry
2012-07-16 09:49:04 +01:00
parent 734f853ca9
commit 09aa9fbeb2
18 changed files with 430 additions and 259 deletions

View File

@ -1,98 +1,102 @@
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn().A());
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
UEqn.clear();
bool closedVolume = false;
if (simple.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::div(phid, p)
- fvm::laplacian(rho*rAU, p)
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi == pEqn.flux();
}
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
);
closedVolume = adjustPhi(phiHbyA, U, p);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rho*rAU, p)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
#include "incompressible/continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
U = HbyA - rAU*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);
if (!simple.transonic())
{
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
}
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
volScalarField rAU(1.0/UEqn().A());
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn().H();
UEqn.clear();
bool closedVolume = false;
if (simple.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::div(phid, p)
- fvm::laplacian(rho*rAU, p)
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi == pEqn.flux();
}
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf())
);
closedVolume = adjustPhi(phiHbyA, U, p);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::div(phiHbyA)
- fvm::laplacian(rho*rAU, p)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (simple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
#include "incompressible/continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
U = HbyA - rAU*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);
if (!simple.transonic())
{
rho.relax();
}
Info<< "rho max/min : "
<< max(rho).value() << " "
<< min(rho).value() << endl;
}

View File

@ -0,0 +1,61 @@
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<rhoThermo> pThermo
(
rhoThermo::New(mesh)
);
rhoThermo& thermo = pThermo();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
thermo.rho()
);
volScalarField& p = thermo.p();
volScalarField& e = thermo.he();
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, simple.dict(), pRefCell, pRefValue);
dimensionedScalar rhoMax(simple.dict().lookup("rhoMax"));
dimensionedScalar rhoMin(simple.dict().lookup("rhoMin"));
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::RASModel> turbulence
(
compressible::RASModel::New
(
rho,
U,
phi,
thermo
)
);
dimensionedScalar initialMass = fvc::domainIntegrate(rho);

View File

@ -1,4 +1,5 @@
{
// Kinetic + pressure energy
volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho);
fvScalarMatrix eEqn
@ -10,7 +11,7 @@
fvc::div(phi)*Ekp - fvc::div(phi, Ekp)
);
//pZones.addEnergySource(thermo, rho, eEqn);
pZones.addEnergySource(thermo, rho, eEqn);
eEqn.relax();
eEqn.solve();

View File

@ -1,52 +1,24 @@
volVectorField HbyA("HbyA", U);
if (pressureImplicitPorosity)
{
HbyA = trTU() & UEqn().H();
}
else
{
HbyA = trAU()*UEqn().H();
}
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
UEqn.clear();
volVectorField HbyA("HbyA", U);
bool closedVolume = false;
if (simple.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf())
);
mrfZones.relativeFlux(fvc::interpolate(psi), phid);
while (simple.correctNonOrthogonal())
if (pressureImplicitPorosity)
{
tmp<fvScalarMatrix> tpEqn;
if (pressureImplicitPorosity)
{
tpEqn = (fvc::div(phid, p) - fvm::laplacian(rho*trTU(), p));
}
else
{
tpEqn = (fvc::div(phid, p) - fvm::laplacian(rho*trAU(), p));
}
tpEqn().setReference(pRefCell, pRefValue);
tpEqn().solve();
if (simple.finalNonOrthogonalIter())
{
phi == tpEqn().flux();
}
HbyA = trTU() & UEqn().H();
}
}
else
{
else
{
HbyA = trAU()*UEqn().H();
}
UEqn.clear();
bool closedVolume = false;
surfaceScalarField phiHbyA
(
"phiHbyA",
@ -79,34 +51,37 @@ else
phi = phiHbyA - tpEqn().flux();
}
}
#include "incompressible/continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
if (pressureImplicitPorosity)
{
U = HbyA - (trTU() & fvc::grad(p));
}
else
{
U = HbyA - trAU()*fvc::grad(p);
}
U.correctBoundaryConditions();
// For closed-volume cases adjust the pressure and density levels
// to obey overall mass continuity
if (closedVolume)
{
const volScalarField& psi = thermo.psi();
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;
}
#include "incompressible/continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
if (pressureImplicitPorosity)
{
U = HbyA - (trTU() & fvc::grad(p));
}
else
{
U = HbyA - trAU()*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;

View File

@ -32,7 +32,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "psiThermo.H"
#include "rhoThermo.H"
#include "RASModel.H"
#include "MRFZones.H"
#include "thermalPorousZones.H"

View File

@ -16,7 +16,6 @@ porousPhi =
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(porousRho*rAUPorous, porousP) == fvc::div(porousPhi)