Commit 19cbd48b authored by Andrew Heather's avatar Andrew Heather
Browse files

INT: Integration of rhoPimpleAdiabaticFoam

Solver for low Mach no. flows with adiabatic thermodynamics and updated
pressure-velocity coupling given by the RCM interpolation procedure
described in

    \verbatim
        Knacke, T. (2013).
        Potential effects of Rhie & Chow type interpolations in airframe
        noise simulations. In: Schram, C., Dénos, R., Lecomte E. (ed):
        Accurate and efficient aeroacoustic prediction approaches for
        airframe noise, VKI LS 2013-03.
    \endverbatim

Original code supplied by Thilo Knacke, CFD E+F GmbH
contact: info@cfd-berlin.com

Integrated into OpenFOAM by OpenCFD Ltd.
parent 57daf418
{
volScalarField& he = thermo.he();
const tmp<volScalarField>& tCp = thermo.Cp();
const tmp<volScalarField>& tCv = thermo.Cv();
const volScalarField& Cp = tCp();
const volScalarField& Cv = tCv();
const scalar gamma = max(Cp/Cv).value();
if (mag(gamma - min(Cp/Cv).value()) > VSMALL)
{
notImplemented("gamma not constant in space");
}
const dictionary& thermoDict = thermo.subDict("mixture");
const dictionary& eosDict = thermoDict.subDict("equationOfState");
bool local = eosDict.lookupOrDefault<bool>("local", false);
// Evolve T as:
//
// T_1 = T_0 \frac{p}{p_0}^{\frac{\gamma - 1}{\gamma}}
if (!local)
{
const scalar T0 = readScalar(eosDict.lookup("T0"));
const scalar p0 = readScalar(eosDict.lookup("p0"));
he = thermo.he(p, pow(p/p0, (gamma - scalar(1))/gamma)*T0);
}
else
{
const volScalarField& T0 = T.oldTime();
const volScalarField& p0 = p.oldTime();
he = thermo.he(p, pow(p/p0, (gamma - scalar(1))/gamma)*T0);
}
thermo.correct();
psi = 1.0/((Cp - Cv)*T);
rho = thermo.rho();
rho.relax();
rho.writeMinMax(Info);
}
rhoPimpleAdiabaticFoam.C
EXE = $(FOAM_APPBIN)/rhoPimpleAdiabaticFoam
EXE_INC = \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
EXE_LIBS = \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lspecie \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-lfiniteVolume \
-lmeshTools \
-lsampling \
-lfvOptions
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);
tmp<fvVectorMatrix> tUEqn
(
fvm::ddt(rho, U) + fvm::div(phi, U)
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(U)
==
fvOptions(rho, U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();
fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
fvOptions.correct(U);
}
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<fluidThermo> pThermo
(
fluidThermo::New(mesh)
);
fluidThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField& p = thermo.p();
volScalarField& T = thermo.T();
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
thermo.rho()
);
Info<< "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(rho)*linearInterpolate(U) & mesh.Sf()
);
Info<< "Calculating face flux field phiByRho\n" << endl;
surfaceScalarField phiByRho
(
IOobject
(
"phiByRho",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
phi/linearInterpolate(rho)
);
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
mesh.setFluxRequired(p.name());
Info<< "Creating field dpdt\n" << endl;
volScalarField dpdt
(
IOobject
(
"dpdt",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
);
#include "createMRF.H"
Info<< "Creating compressibility field psi\n" << endl;
volScalarField psi("psi", 1.0/((thermo.Cp() - thermo.Cv())*T));
psi.oldTime() = 1.0/((thermo.Cp() - thermo.Cv())*T.oldTime());
psi.oldTime().oldTime() = 1.0/((thermo.Cp()-thermo.Cv())*T.oldTime().oldTime());
{
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
// Define coefficients and pseudo-velocities for RCM interpolation
// M[U] = AU - H = -grad(p)
// U = H/A - 1/A grad(p)
// H/A = U + 1/A grad(p)
surfaceScalarField rhorAUf
(
"rhorAUf",
fvc::interpolate(rho)/fvc::interpolate(UEqn.A())
);
surfaceVectorField rhoHbyAf
(
"rhoHbyAf",
fvc::interpolate(rho)*fvc::interpolate(U)
+ rhorAUf*fvc::interpolate(fvc::grad(p))
);
#include "resetBoundaries.H"
if (pimple.nCorrPISO() <= 1)
{
tUEqn.clear();
}
if (pimple.transonic())
{
FatalError
<< "\nTransonic option not available for " << args.executable()
<< exit(FatalError);
}
else
{
// Rhie & Chow interpolation (part 1)
surfaceScalarField phiHbyA
(
"phiHbyA",
(
(rhoHbyAf & mesh.Sf())
+ rhorAUf*fvc::interpolate(rho)*fvc::ddtCorr(U, phiByRho)
+ fvc::interpolate(rho)
* fvc::alphaCorr(U, phiByRho, pimple.finalInnerIter())
)
);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
// Pressure corrector
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.finalInnerIter())));
// Rhie & Chow interpolation (part 2)
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
phiByRho = phi/fvc::interpolate(rho);
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
}
rho = thermo.rho();
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
{
// Keep standard formulation on domain boundaries to ensure compatibility
// with existing boundary conditions
const Foam::FieldField<Foam::fvsPatchField, scalar> rhorAUf_orig
(
fvc::interpolate(rho.boundaryField()*rAU.boundaryField())
);
const Foam::FieldField<Foam::fvsPatchField, vector> rhoHbyA_orig
(
fvc::interpolate(rho.boundaryField()*HbyA.boundaryField())
);
surfaceScalarField::Boundary& rhorAUfbf = rhorAUf.boundaryFieldRef();
surfaceVectorField::Boundary& rhoHbyAfbf = rhoHbyAf.boundaryFieldRef();
forAll(U.boundaryField(), patchi)
{
if (!U.boundaryField()[patchi].coupled())
{
rhorAUfbf[patchi] = rhorAUf_orig[patchi];
rhoHbyAfbf[patchi] = rhoHbyA_orig[patchi];
}
}
}
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Application
rhoPimpleAdiabaticFoam
Description
Transient solver for laminar or turbulent flow of weakly compressible
fluids for low Mach number aeroacoustic applications.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient simulations. The RCM interpolation is used as in
\verbatim
Knacke, T. (2013).
Potential effects of Rhie & Chow type interpolations in airframe
noise simulations. In: Schram, C., Dénos, R., Lecomte E. (ed):
Accurate and efficient aeroacoustic prediction approaches for
airframe noise, VKI LS 2013-03.
\endverbatim
Contact: info@cfd-berlin.com
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "fluidThermo.H"
#include "turbulentFluidThermoModel.H"
#include "bound.H"
#include "pimpleControl.H"
#include "fvOptions.H"
#include "ddtScheme.H"
#include "fvcCorrectAlpha.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
if (pimple.nCorrPIMPLE() <= 1)
{
#include "rhoEqn.H"
}
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
U.storePrevIter();
rho.storePrevIter();
phi.storePrevIter();
phiByRho.storePrevIter();
#include "UEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
#include "EEqn.H"
if (pimple.turbCorr())
{
turbulence->correct();
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
......@@ -401,6 +401,7 @@ finiteVolume/fvc/fvcFlux.C
finiteVolume/fvc/fvcMeshPhi.C
finiteVolume/fvc/fvcSmooth/fvcSmooth.C
finiteVolume/fvc/fvcReconstructMag.C
finiteVolume/fvc/fvcCorrectAlpha.C
general = cfdTools/general
$(general)/findRefCell/findRefCell.C
......
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -98,7 +98,19 @@ public:
backwardDdtScheme(const fvMesh& mesh, Istream& is)
:
ddtScheme<Type>(mesh, is)
{}
{
if (is.good() && !is.eof())
{
this->ddtPhiCoeff_ = readScalar(is);
}
// Ensure the old-old-time cell volumes are available
// for moving meshes
if (mesh.moving())
{
mesh.V00();
}
}
// Member Functions
......
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -156,13 +156,36 @@ tmp<surfaceScalarField> ddtScheme<Type>::fvcDdtPhiCoeff
const fluxFieldType& phiCorr
)
{
tmp<surfaceScalarField> tddtCouplingCoeff = scalar(1)
- min
tmp<surfaceScalarField> tddtCouplingCoeff
(
new surfaceScalarField
(
IOobject
(
"ddtCouplingCoeff",
U.mesh().time().timeName(),
U.mesh()
),
U.mesh(),
dimensionedScalar("one", dimless, 1.0)
)
);
if (ddtPhiCoeff_ < 0)
{
tddtCouplingCoeff.ref() =- min
(
mag(phiCorr)
/(mag(phi) + dimensionedScalar("small", phi.dimensions(), SMALL)),
scalar(1)
);
}
else
{
tddtCouplingCoeff.ref() =
dimensionedScalar("ddtPhiCoeff", dimless, ddtPhiCoeff_);
}
surfaceScalarField& ddtCouplingCoeff = tddtCouplingCoeff.ref();
......
......@@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -30,6 +30,7 @@ Group
Description
Abstract base class for ddt schemes.
SourceFiles
ddtScheme.C
......@@ -76,6 +77,9 @@ protected:
const fvMesh& mesh_;
//- Input for fvcDdtPhiCoeff (-1 default)
scalar ddtPhiCoeff_;
// Private Member Functions
......@@ -109,13 +113,15 @@ public:
//- Construct from mesh
ddtScheme(const fvMesh& mesh)
:
mesh_(mesh)
mesh_(mesh),
ddtPhiCoeff_(-1)
{}
//- Construct from mesh and Istream
ddtScheme(const fvMesh& mesh, Istream&)
ddtScheme(const fvMesh& mesh, Istream& is)
:
mesh_(mesh)
mesh_(mesh),
ddtPhiCoeff_(-1)
{}
......