ENH: rhePimpleFoam: Merged dynamic mesh functionality of rhoPimpleDyMFoam into rhoPimpleFoam

and replaced rhoPimpleDyMFoam with a script which reports this change.

The rhoPimpleDyMFoam tutorials have been moved into the rhoPimpleFoam directory.

This change is the first of a set of developments to merge dynamic mesh
functionality into the standard solvers to improve consistency, usability,
flexibility and maintainability of these solvers.

Henry G. Weller
CFD Direct Ltd.

rhoReactingFoam: Updated for changes to rhoPimpleFoam files
This commit is contained in:
Henry Weller 2017-11-23 12:13:37 +00:00 committed by Andrew Heather
parent 81cea09983
commit 2f888d1684
57 changed files with 218 additions and 359 deletions

View File

@ -53,10 +53,11 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUfIfPresent.H"
#include "createTimeControls.H"
turbulence->validate();

View File

@ -7,6 +7,8 @@ EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude
EXE_LIBS = \
-lcompressibleTransportModels \
@ -17,4 +19,7 @@ EXE_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lsampling \
-lfvOptions
-lfvOptions \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-ldynamicMesh

View File

@ -6,6 +6,6 @@ CorrectPhi
rho,
psi,
dimensionedScalar("rAUf", dimTime, 1),
divrhoU,
divrhoU(),
pimple
);

View File

@ -20,9 +20,10 @@ surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::interpolate(rho)*fvc::flux(HbyA)
+ rhorAUf*fvc::ddtCorr(rho, U, phi)
+ rhorAUf*fvc::ddtCorr(rho, U, phi, rhoUf)
);
fvc::makeRelative(phiHbyA, rho, U);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
@ -103,7 +104,15 @@ if (pressureControl.limit(p))
thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
rho = thermo.rho();
// Correct rhoUf if the mesh is moving
fvc::correctRhoUf(rhoUf, rho, U, phi);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
if (mesh.moving())
{
dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
}
}

View File

@ -1,3 +1,5 @@
InfoInFunction << "consistent" << endl;
if (!pimple.SIMPLErho())
{
rho = thermo.rho();
@ -21,10 +23,11 @@ surfaceScalarField phiHbyA
"phiHbyA",
(
fvc::interpolate(rho)*fvc::flux(HbyA)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+ fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi, rhoUf)
)
);
fvc::makeRelative(phiHbyA, rho, U);
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
volScalarField rhorAtU("rhorAtU", rho*rAtU);
@ -114,7 +117,15 @@ if (pressureControl.limit(p))
thermo.correctRho(psi*p - psip0, rhoMin, rhoMax) ;
rho = thermo.rho();
// Correct rhoUf if the mesh is moving
fvc::correctRhoUf(rhoUf, rho, U, phi);
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
if (mesh.moving())
{
dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
}
}

View File

@ -0,0 +1,10 @@
#include "readTimeControls.H"
bool correctPhi = pimple.dict().lookupOrDefault<Switch>
(
"correctPhi",
!isA<staticFvMesh>(mesh)
);
bool checkMeshCourantNo =
pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);

View File

@ -1,3 +0,0 @@
rhoPimpleDyMFoam.C
EXE = $(FOAM_APPBIN)/rhoPimpleDyMFoam

View File

@ -1,26 +0,0 @@
EXE_INC = \
-I.. \
-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 \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude
EXE_LIBS = \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lspecie \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-lfiniteVolume \
-lmeshTools \
-lsampling \
-lfvOptions \
-ldynamicFvMesh \
-ltopoChangerFvMesh \
-ldynamicMesh

View File

@ -1,11 +0,0 @@
#include "createTimeControls.H"
bool correctPhi
(
pimple.dict().lookupOrDefault("correctPhi", true)
);
bool checkMeshCourantNo
(
pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
);

View File

@ -1,121 +0,0 @@
if (!pimple.SIMPLErho())
{
rho = thermo.rho();
}
// Thermodynamic density needs to be updated by psi*d(p) after the
// pressure solution
const volScalarField psip0(psi*p);
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
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::interpolate(rho))*phiHbyA
);
phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(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 = phiHbyA + pEqn.flux();
}
}
}
else
{
while (pimple.correctNonOrthogonal())
{
// Pressure corrector
fvScalarMatrix pEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p))
+ fvc::div(phiHbyA)
- fvm::laplacian(rhorAUf, p)
==
fvOptions(psi, p, rho.name())
);
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();
U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
if (pressureControl.limit(p))
{
p.correctBoundaryConditions();
}
thermo.correctRho(psi*p - psip0, rhoMin, rhoMax);
rho = thermo.rho();
{
rhoUf = fvc::interpolate(rho*U);
surfaceVectorField n(mesh.Sf()/mesh.magSf());
rhoUf += n*(fvc::absolute(phi, rho, U)/mesh.magSf() - (n & rhoUf));
}
if (thermo.dpdt())
{
dpdt = fvc::ddt(p);
if (mesh.moving())
{
dpdt -= fvc::div(fvc::meshPhi(rho, U), p);
}
}

View File

@ -1,6 +0,0 @@
#include "readTimeControls.H"
correctPhi = pimple.dict().lookupOrDefault("correctPhi", true);
checkMeshCourantNo =
pimple.dict().lookupOrDefault("checkMeshCourantNo", false);

View File

@ -1,168 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2017 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
rhoPimpleDyMFoam
Group
grpCompressibleSolvers grpMovingMeshSolvers
Description
Transient solver for turbulent flow of compressible fluids for HVAC and
similar applications, with optional mesh motion and mesh topology changes.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient simulations.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.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"
#include "fvcSmooth.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUf.H"
#include "createControls.H"
turbulence->validate();
if (!LTS)
{
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readControls.H"
{
// Store divrhoU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
volScalarField divrhoU
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
);
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
}
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// Store momentum to set rhoUf for introduced faces.
volVectorField rhoU("rhoU", rho*U);
// Do any mesh changes
mesh.update();
if (mesh.changing())
{
MRF.update();
if (correctPhi)
{
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & rhoUf;
#include "correctPhi.H"
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
#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();
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -29,7 +29,7 @@ Group
Description
Transient solver for turbulent flow of compressible fluids for HVAC and
similar applications.
similar applications, with optional mesh motion and mesh topology changes.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient simulations.
@ -37,11 +37,13 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.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"
#include "fvcSmooth.H"
@ -54,12 +56,13 @@ int main(int argc, char *argv[])
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createTimeControls.H"
#include "initContinuityErrs.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUfIfPresent.H"
#include "createTimeControls.H"
turbulence->validate();
@ -75,21 +78,69 @@ int main(int argc, char *argv[])
while (runTime.run())
{
#include "readTimeControls.H"
#include "readControls.H"
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
}
// Store divrhoU from the previous mesh so that it can be mapped
// and used in correctPhi to ensure the corrected phi has the
// same divergence
autoPtr<volScalarField> divrhoU;
if (correctPhi)
{
divrhoU = new volScalarField
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
);
}
runTime++;
if (LTS)
{
#include "setRDeltaT.H"
}
else
{
#include "compressibleCourantNo.H"
#include "setDeltaT.H"
}
Info<< "Time = " << runTime.timeName() << nl << endl;
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
// Store momentum to set rhoUf for introduced faces.
autoPtr<volVectorField> rhoU;
if (rhoUf.valid())
{
rhoU = new volVectorField("rhoU", rho*U);
}
// Do any mesh changes
mesh.update();
#include "updateRhoUf.H"
if (mesh.changing())
{
MRF.update();
if (correctPhi)
{
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & rhoUf();
#include "correctPhi.H"
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
}
if (pimple.nCorrPIMPLE() <= 1)
{

View File

@ -1,6 +1,6 @@
EXE_INC = \
-I.. \
-I../../rhoPimpleFoam/rhoPimpleDyMFoam \
-I../../rhoPimpleFoam \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

View File

@ -51,7 +51,7 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createControls.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUf.H"

View File

@ -1,7 +1,7 @@
EXE_INC = \
-I.. \
-I../../reactingParcelFoam \
-I../../../compressible/rhoPimpleFoam/rhoPimpleDyMFoam \
-I../../../compressible/rhoPimpleFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I${LIB_SRC}/meshTools/lnInclude \
-I${LIB_SRC}/sampling/lnInclude \

View File

@ -54,7 +54,7 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createControl.H"
#include "createControls.H"
#include "createTimeControls.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "createRhoUf.H"

1
bin/rhoPimpleDyMFoam Symbolic link
View File

@ -0,0 +1 @@
mergedDyM

View File

@ -0,0 +1,51 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 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/>.
Global
createRhoUf
Description
Creates and initialises the velocity field rhoUf if present.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
autoPtr<surfaceVectorField> rhoUf;
IOobject rhoUfHeader
(
"rhoUf",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
);
if (rhoUfHeader.typeHeaderOk<surfaceVectorField>(true))
{
Info<< "Reading face momentum rhoUf\n" << endl;
rhoUf = new surfaceVectorField(rhoUfHeader, mesh);
}
// ************************************************************************* //

View File

@ -0,0 +1,55 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2017 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/>.
Global
updateRhoUf
Description
Constructs the face momentum field rhoUf if not already constructed.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
if (mesh.changing())
{
if (!rhoUf.valid())
{
Info<< "Constructing face momentum rhoUf" << endl;
rhoUf = new surfaceVectorField
(
IOobject
(
"rhoUf",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(rho*U)
);
}
}
// ************************************************************************* //

View File

@ -15,7 +15,7 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application rhoPimpleDyMFoam;
application rhoPimpleFoam;
startFrom latestTime;