zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

rhoPimpleFoam: Added support for transonic flow of liquids and real gases

Browse Source 
Both stardard SIMPLE and the SIMPLEC (using the 'consistent' option in
fvSolution) are now supported for both subsonic and transonic flow of all
fluid types.

rhoPimpleFoam now instantiates the lower-level fluidThermo which instantiates
either a psiThermo or rhoThermo according to the 'type' specification in
thermophysicalProperties, see also commit a1c8cde310
This commit is contained in:
Henry Weller
2017-02-28 11:14:59 +00:00
parent 7515ea7550
commit 50516486a4
21 changed files with 295 additions and 343 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
applications/solvers/compressible/rhoPimpleFoam/createFields.H
View File

@ -2,11 +2,11 @@
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<psiThermo> pThermo
autoPtr<fluidThermo> pThermo
(
psiThermo::New(mesh)
fluidThermo::New(mesh)
);
psiThermo& thermo = pThermo();
fluidThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField& p = thermo.p();
@ -40,27 +40,7 @@ volVectorField U
#include "compressibleCreatePhi.H"
dimensionedScalar rhoMax
(
dimensionedScalar::lookupOrDefault
(
"rhoMax",
pimple.dict(),
dimDensity,
GREAT
)
);
dimensionedScalar rhoMin
(
dimensionedScalar::lookupOrDefault
(
"rhoMin",
pimple.dict(),
dimDensity,
0
)
);
pressureControl pressureControl(p, rho, pimple.dict(), false);
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence

65
applications/solvers/compressible/rhoPimpleFoam/pEqn.H
View File

@ -1,8 +1,3 @@
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@ -12,55 +7,54 @@ if (pimple.nCorrPISO() <= 1)
tUEqn.clear();
}
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ 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);
if (pimple.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)
*(
fvc::flux(HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
)
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
MRF.makeRelative(fvc::interpolate(psi), phid);
phiHbyA -= fvc::interpolate(p)*phid;
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi == pEqn.flux();
phi = phiHbyA + pEqn.flux();
}
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ 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
@ -87,19 +81,20 @@ else
// Explicitly relax pressure for momentum corrector
p.relax();
// Recalculate density from the relaxed pressure
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
Info<< "rho max/min : " << max(rho).value()
<< " " << min(rho).value() << endl;
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
pressureControl.limit(p);
p.correctBoundaryConditions();
rho = thermo.rho();
if (!pimple.transonic())
{
rho.relax();
}
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);

76
applications/solvers/compressible/rhoPimpleFoam/pcEqn.H
View File

@ -1,8 +1,3 @@
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn.A());
volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1()));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@ -12,72 +7,64 @@ if (pimple.nCorrPISO() <= 1)
tUEqn.clear();
}
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
volScalarField rhorAtU("rhorAtU", rho*rAtU);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
if (pimple.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)
*(
fvc::flux(HbyA)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
/fvc::interpolate(rho)
)
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
MRF.makeRelative(fvc::interpolate(psi), phid);
surfaceScalarField phic
(
"phic",
phiHbyA +=
fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf()
);
- fvc::interpolate(p)*phid;
HbyA -= (rAU - rAtU)*fvc::grad(p);
volScalarField rhorAtU("rhorAtU", rho*rAtU);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
+ fvm::div(phid, p)
+ fvc::div(phic)
- fvm::laplacian(rhorAtU, p)
==
fvOptions(psi, p, rho.name())
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi == phic + pEqn.flux();
phi = phiHbyA + pEqn.flux();
}
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
(
fvc::flux(rho*HbyA)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU)*fvc::grad(p);
volScalarField rhorAtU("rhorAtU", rho*rAtU);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
@ -109,19 +96,16 @@ U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
// Recalculate density from the relaxed pressure
pressureControl.limit(p);
p.correctBoundaryConditions();
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
if (!pimple.transonic())
{
rho.relax();
}
Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl;
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}

64
applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
View File

@ -1,8 +1,3 @@
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
@ -12,55 +7,53 @@ if (pimple.nCorrPISO() <= 1)
tUEqn.clear();
}
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(rho*HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
);
fvc::makeRelative(phiHbyA, rho, U);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
if (pimple.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)
*(
fvc::flux(HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, rhoUf)/fvc::interpolate(rho)
)
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
fvc::makeRelative(phid, psi, U);
MRF.makeRelative(fvc::interpolate(psi), phid);
phiHbyA -= fvc::interpolate(p)*phid;
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
// Relax the pressure equation to ensure diagonal-dominance
pEqn.relax();
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi == pEqn.flux();
phi = phiHbyA + pEqn.flux();
}
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(rho*HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, rhoUf)
);
fvc::makeRelative(phiHbyA, rho, U);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
while (pimple.correctNonOrthogonal())
{
// Pressure corrector
@ -88,19 +81,20 @@ else
// Explicitly relax pressure for momentum corrector
p.relax();
// Recalculate density from the relaxed pressure
rho = thermo.rho();
rho = max(rho, rhoMin);
rho = min(rho, rhoMax);
rho.relax();
Info<< "rho max/min : " << max(rho).value()
<< " " << min(rho).value() << endl;
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
pressureControl.limit(p);
p.correctBoundaryConditions();
rho = thermo.rho();
if (!pimple.transonic())
{
rho.relax();
}
{
rhoUf = fvc::interpolate(rho*U);
surfaceVectorField n(mesh.Sf()/mesh.magSf());

10
applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -22,10 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
rhoPimpleFoam
Group
grpCompressibleSolvers grpMovingMeshSolvers
rhoPimpleDyMFoam
Description
Transient solver for turbulent flow of compressible fluids for HVAC and
@ -38,10 +35,11 @@ Description
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "psiThermo.H"
#include "fluidThermo.H"
#include "turbulentFluidThermoModel.H"
#include "bound.H"
#include "pimpleControl.H"
#include "pressureControl.H"
#include "CorrectPhi.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"

5
applications/solvers/compressible/rhoPimpleFoam/rhoPimpleFoam.C
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -34,10 +34,11 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "psiThermo.H"
#include "fluidThermo.H"
#include "turbulentFluidThermoModel.H"
#include "bound.H"
#include "pimpleControl.H"
#include "pressureControl.H"
#include "fvOptions.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"

4
applications/solvers/compressible/rhoSimpleFoam/createFields.H
View File

@ -7,6 +7,8 @@ autoPtr<fluidThermo> pThermo
fluidThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField& p = thermo.p();
volScalarField rho
(
IOobject
@ -20,8 +22,6 @@ volScalarField rho
thermo.rho()
);
volScalarField& p = thermo.p();
Info<< "Reading field U\n" << endl;
volVectorField U
(

18
applications/solvers/compressible/rhoSimpleFoam/pEqn.H
View File

@ -6,16 +6,18 @@
bool closedVolume = false;
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
if (simple.transonic())
{
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
surfaceScalarField rhof(fvc::interpolate(rho));
MRF.makeRelative(rhof, phiHbyA);
surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/rhof)*phiHbyA
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
phiHbyA -= fvc::interpolate(p)*phid;
@ -49,14 +51,8 @@
}
else
{
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
closedVolume = adjustPhi(phiHbyA, U, p);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn

28
applications/solvers/compressible/rhoSimpleFoam/pcEqn.H
View File

@ -5,16 +5,20 @@ tUEqn.clear();
bool closedVolume = false;
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
volScalarField rhorAtU("rhorAtU", rho*rAtU);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
if (simple.transonic())
{
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
surfaceScalarField rhof(fvc::interpolate(rho));
MRF.makeRelative(rhof, phiHbyA);
surfaceScalarField phid
(
"phid",
(fvc::interpolate(psi)/rhof)*phiHbyA
(fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
);
phiHbyA +=
@ -23,14 +27,12 @@ if (simple.transonic())
HbyA -= (rAU - rAtU)*fvc::grad(p);
volScalarField rhorAtU("rhorAtU", rho*rAtU);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::div(phid, p)
+ fvc::div(phiHbyA)
fvc::div(phiHbyA)
+ fvm::div(phid, p)
- fvm::laplacian(rhorAtU, p)
==
fvOptions(psi, p, rho.name())
@ -55,19 +57,11 @@ if (simple.transonic())
}
else
{
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(rho*HbyA));
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
closedVolume = adjustPhi(phiHbyA, U, p);
phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
HbyA -= (rAU - rAtU)*fvc::grad(p);
volScalarField rhorAtU("rhorAtU", rho*rAtU);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF);
while (simple.correctNonOrthogonal())
{
fvScalarMatrix pEqn

59
applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createFields.H
View File

@ -1,59 +0,0 @@
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<fluidThermo> pThermo
(
fluidThermo::New(mesh)
);
fluidThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
thermo.rho()
);
volScalarField& p = thermo.p();
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "compressibleCreatePhi.H"
pressureControl pressureControl(p, rho, simple.dict());
mesh.setFluxRequired(p.name());
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
dimensionedScalar initialMass = fvc::domainIntegrate(rho);
#include "createMRF.H"