Commit 042354a8 authored by Henry's avatar Henry
Browse files

rhoPimpleDyMFoam: New version of rhoPimpleFoam for moving-mesh cases

parent c205857f
......@@ -5,5 +5,6 @@ set -x
wmake
wmake rhoPimplecFoam
wmake rhoLTSPimpleFoam
wmake rhoPimpleDyMFoam
# ----------------------------------------------------------------- end-of-file
rhoPimpleDyMFoam.C
EXE = $(FOAM_APPBIN)/rhoPimpleDyMFoam
EXE_INC = \
-I.. \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
-I$(LIB_SRC)/finiteVolume/cfdTools \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
EXE_LIBS = \
-lfluidThermophysicalModels \
-lspecie \
-lcompressibleTurbulenceModel \
-lcompressibleRASModels \
-lcompressibleLESModels \
-lfiniteVolume \
-lmeshTools \
-lsampling \
-lfvOptions \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-ldynamicMesh \
-lmeshTools
{
if (mesh.changing())
{
forAll(U.boundaryField(), patchi)
{
if (U.boundaryField()[patchi].fixesValue())
{
U.boundaryField()[patchi].initEvaluate();
}
}
forAll(U.boundaryField(), patchi)
{
if (U.boundaryField()[patchi].fixesValue())
{
U.boundaryField()[patchi].evaluate();
phi.boundaryField()[patchi] =
rho.boundaryField()[patchi]
*(
U.boundaryField()[patchi]
& mesh.Sf().boundaryField()[patchi]
);
}
}
}
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("pcorr", p.dimensions(), 0.0),
pcorrTypes
);
dimensionedScalar Dp("Dp", dimTime, 1.0);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(Dp, pcorr) == fvc::div(phi) - divrhoU
);
pcorrEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phi -= pcorrEqn.flux();
}
}
}
wordList pcorrTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
for (label i=0; i<p.boundaryField().size(); i++)
{
if (p.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
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();
if (pimple.nCorrPISO() <= 1)
{
UEqn.clear();
}
if (pimple.transonic())
{
surfaceScalarField phid
(
"phid",
fvc::interpolate(psi)
*(
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
)
);
fvOptions.relativeFlux(fvc::interpolate(psi), phid);
volScalarField Dp("Dp", rho*rAU);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvm::div(phid, p)
- fvm::laplacian(Dp, p)
==
fvOptions(psi, p, rho.name())
);
fvOptions.constrain(pEqn);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi == pEqn.flux();
}
}
}
else
{
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)
*(
(fvc::interpolate(HbyA) & mesh.Sf())
- fvc::meshPhi(rho, U)
+ fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
)
);
fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);
volScalarField Dp("Dp", rho*rAU);
while (pimple.correctNonOrthogonal())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvm::ddt(psi, p)
+ fvc::div(phiHbyA)
- fvm::laplacian(Dp, p)
==
fvOptions(psi, p, rho.name())
);
fvOptions.constrain(pEqn);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phi = phiHbyA + pEqn.flux();
}
}
}
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
// 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);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
}
#include "readTimeControls.H"
bool correctPhi =
pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
bool checkMeshCourantNo =
pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\/ 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
rhoPimpleFoam
Description
Transient solver for laminar or turbulent flow of compressible fluids
for HVAC and similar applications.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient simulations.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "psiThermo.H"
#include "turbulenceModel.H"
#include "bound.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
#include "readControls.H"
#include "createFields.H"
#include "createFvOptions.H"
#include "createPcorrTypes.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
// Create old-time absolute flux for ddtPhiCorr
surfaceScalarField phiAbs("phiAbs", phi);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readControls.H"
#include "compressibleCourantNo.H"
// Make the fluxes absolute before mesh-motion
fvc::makeAbsolute(phi, rho, U);
// Update absolute flux for ddtPhiCorr
phiAbs = phi;
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
{
// Store divrhoU from the previous time-step/mesh for the correctPhi
volScalarField divrhoU(fvc::div(phi));
// Do any mesh changes
mesh.update();
if (mesh.changing() && correctPhi)
{
#include "correctPhi.H"
}
}
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);
if (mesh.changing() && checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
if (pimple.nCorrPIMPLE() <= 1)
{
#include "rhoEqn.H"
Info<< "rhoEqn max/min : " << max(rho).value()
<< " " << min(rho).value() << endl;
}
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
#include "UEqn.H"
#include "EEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.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;
}
// ************************************************************************* //
EXE_INC = \
-I.. \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/fvOptions/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
EXE_LIBS = \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-ldynamicMesh \
-lmeshTools \
-lincompressibleTransportModels \
-lincompressibleTurbulenceModel \
-lincompressibleRASModels \
-lincompressibleLESModels \
-lfiniteVolume \
-lfvOptions \
-lsampling
-lsampling \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-ldynamicMesh \
-lmeshTools
// Solve the Momentum equation
tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
==
fvOptions(U)
);
UEqn().relax();
fvOptions.constrain(UEqn());
rAU = 1.0/UEqn().A();
if (pimple.momentumPredictor())
{
solve(UEqn() == -fvc::grad(p));
fvOptions.correct(U);
}
......@@ -22,20 +22,6 @@
}
}
wordList pcorrTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
forAll(p.boundaryField(), patchI)
{
if (p.boundaryField()[patchI].fixesValue())
{
pcorrTypes[patchI] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField pcorr
(
IOobject
......@@ -51,11 +37,13 @@
pcorrTypes
);
dimensionedScalar Dp("Dp", dimTime, 1.0);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::laplacian(rAU, pcorr) == fvc::div(phi)
fvm::laplacian(Dp, pcorr) == fvc::div(phi)
);
pcorrEqn.setReference(pRefCell, pRefValue);
......@@ -68,6 +56,4 @@
}
}
phi.oldTime() = phi;
#include "continuityErrs.H"
......@@ -40,19 +40,3 @@
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
Info<< "Reading field rAU if present\n" << endl;
volScalarField rAU
(
IOobject
(
"rAU",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
runTime.deltaT(),
zeroGradientFvPatchScalarField::typeName
);
wordList pcorrTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
for (label i=0; i<p.boundaryField().size(); i++)
{
if (p.boundaryField()[i].fixesValue())
{
pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
}
}
......@@ -10,13 +10,9 @@ surfaceScalarField phiHbyA
(
"phiHbyA",
(fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, U, phiAbs)
);
if (ddtPhiCorr)
{
phiHbyA += fvc::ddtPhiCorr(rAU, U, phi);
}
if (p.needReference())
{
fvc::makeRelative(phiHbyA, U);
......
......@@ -33,9 +33,9 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "dynamicFvMesh.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
......@@ -44,15 +44,21 @@ Description
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "initContinuityErrs.H"
pimpleControl pimple(mesh);
#include "createFields.H"
#include "createFvOptions.H"
#include "readTimeControls.H"
#include "createPcorrTypes.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"
pimpleControl pimple(mesh);
// Create old-time absolute flux for ddtPhiCorr
surfaceScalarField phiAbs("phiAbs", phi);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -66,6 +72,9 @@ int main(int argc, char *argv[])
// Make the fluxes absolute
fvc::makeAbsolute(phi, U);
// Update absolute flux for ddtPhiCorr
phiAbs = phi;
#include "setDeltaT.H"
runTime++;
......
#include "readTimeControls.H"
const dictionary& pimpleDict = pimple.dict();
const bool correctPhi =
pimpleDict.lookupOrDefault("correctPhi", false);
pimple.dict().lookupOrDefault("correctPhi", false);