Commit 5c58b07f authored by sergio's avatar sergio Committed by Andrew Heather
Browse files

ENH: Adding overBuoyantPimpleFoam and tutorial

parent aaf5d7be
overBuoyantPimpleDyMFoam.C
EXE = $(FOAM_APPBIN)/overBuoyantPimpleDyMFoam
EXE_INC = \
-I$(FOAM_SOLVERS)/heatTransfer/buoyantPimpleFoam \
-I$(FOAM_SOLVERS)/compressible/rhoPimpleFoam/overRhoPimpleDyMFoam \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/overset/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/transportModels/compressible/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude
EXE_LIBS = \
-lfiniteVolume \
-lsampling \
-lmeshTools \
-lfvOptions \
-lcompressibleTransportModels \
-lfluidThermophysicalModels \
-lradiationModels \
-lspecie \
-lturbulenceModels \
-lcompressibleTurbulenceModels \
-latmosphericModels \
-loverset \
-ldynamicFvMesh \
-ltopoChangerFvMesh
if (mesh.changing())
{
volVectorField::Boundary& bfld = U.boundaryFieldRef();
forAll(bfld, patchi)
{
if (bfld[patchi].fixesValue())
{
bfld[patchi].initEvaluate();
}
}
surfaceScalarField::Boundary& phiBfld = phi.boundaryFieldRef();
forAll(bfld, patchi)
{
if (bfld[patchi].fixesValue())
{
bfld[patchi].evaluate();
phiBfld[patchi] =
rho.boundaryField()[patchi]
* (
bfld[patchi]
& mesh.Sf().boundaryField()[patchi]
);
}
}
}
// Initialize BCs list for pcorr to zero-gradient
wordList pcorrTypes
(
p.boundaryField().size(),
zeroGradientFvPatchScalarField::typeName
);
// Set BCs of pcorr to fixed-value for patches at which p is fixed
forAll(p.boundaryField(), patchi)
{
if (p.boundaryField()[patchi].fixesValue())
{
pcorrTypes[patchi] = fixedValueFvPatchScalarField::typeName;
}
}
volScalarField pcorr
(
IOobject
(
"pcorr",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar(p.dimensions(), Zero),
pcorrTypes
);
mesh.setFluxRequired(pcorr.name());
{
dimensionedScalar rAUf("rAUf", dimTime, 1.0);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pcorrEqn
(
fvm::ddt(psi, pcorr)
+ fvc::div(phi)
- fvm::laplacian(rAUf, pcorr)
==
divrhoU()
);
//pcorrEqn.solve(mesh.solver(pcorr.select(pimple.finalInnerIter())));
//Bypass virtual layer
const dictionary& d = mesh.solver
(
pcorr.select
(
pimple.finalInnerIter()
)
);
mesh.fvMesh::solve(pcorrEqn, d);
if (pimple.finalNonOrthogonalIter())
{
phi += pcorrEqn.flux();
}
}
}
Info<< "Reading thermophysical properties\n" << endl;
autoPtr<rhoThermo> pThermo(rhoThermo::New(mesh));
rhoThermo& thermo = pThermo();
thermo.validate(args.executable(), "h", "e");
volScalarField rho
(
IOobject
(
"rho",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_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"
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
#include "readGravitationalAcceleration.H"
#include "readhRef.H"
#include "gh.H"
Info<< "Reading field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
// Force p_rgh to be consistent with p
p_rgh = p - rho*gh;
mesh.setFluxRequired(p_rgh.name());
label pRefCell = 0;
scalar pRefValue = 0.0;
if (p_rgh.needReference())
{
setRefCell
(
p,
p_rgh,
pimple.dict(),
pRefCell,
pRefValue
);
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
}
dimensionedScalar initialMass("initialMass", fvc::domainIntegrate(rho));
#include "createDpdt.H"
#include "createK.H"
#include "createMRF.H"
#include "createRadiationModel.H"
#include "createFvOptions.H"
//- Overset specific
// Add solver-specific interpolations
{
wordHashSet& nonInt =
const_cast<wordHashSet&>(Stencil::New(mesh).nonInterpolatedFields());
nonInt.insert("HbyA");
nonInt.insert("grad(p_rgh)");
nonInt.insert("surfaceIntegrate(phi)");
nonInt.insert("surfaceIntegrate(phiHbyA)");
nonInt.insert("cellMask");
nonInt.insert("cellDisplacement");
nonInt.insert("interpolatedCells");
nonInt.insert("cellInterpolationWeight");
}
// Mask field for zeroing out contributions on hole cells
#include "createCellMask.H"
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2019 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
overBuoyantPimpleDymFoam
Group
grpHeatTransferSolvers
Description
Transient solver for buoyant, turbulent flow of compressible fluids for
ventilation and heat-transfer with overset feature
Turbulence is modelled using a run-time selectable compressible RAS or
LES model.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "dynamicFvMesh.H"
#include "rhoThermo.H"
#include "turbulentFluidThermoModel.H"
#include "radiationModel.H"
#include "fvOptions.H"
#include "pimpleControl.H"
#include "pressureControl.H"
#include "CorrectPhi.H"
#include "cellCellStencilObject.H"
#include "localMin.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
argList::addNote
(
"Transient solver for buoyant, turbulent fluid flow"
" of compressible fluids, including radiation."
);
#include "postProcess.H"
#include "addCheckCaseOptions.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createDynamicFvMesh.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createFieldRefs.H"
#include "initContinuityErrs.H"
#include "createRhoUfIfPresent.H"
#include "createControls.H"
#include "compressibleCourantNo.H"
#include "setInitialDeltaT.H"
turbulence->validate();
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
while (runTime.run())
{
#include "readTimeControls.H"
#include "readControls.H"
#include "readDyMControls.H"
#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.reset
(
new volScalarField
(
"divrhoU",
fvc::div(fvc::absolute(phi, rho, U))
)
);
}
++runTime;
Info<< "Time = " << runTime.timeName() << nl << endl;
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
if (pimple.firstIter() || moveMeshOuterCorrectors)
{
// Do any mesh changes
mesh.update();
if (mesh.changing())
{
MRF.update();
#include "setCellMask.H"
const surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
// Zero Uf on old faceMask (H-I)
rhoUf() *= faceMaskOld;
//fvc::correctRhoUf(rhoUfint, rho, U, phi);
surfaceVectorField rhoUfint(fvc::interpolate(rho*U));
// Update Uf and phi on new C-I faces
rhoUf() += (1-faceMaskOld)*rhoUfint;
// Update Uf boundary
forAll(rhoUf().boundaryField(), patchI)
{
rhoUf().boundaryFieldRef()[patchI] =
rhoUfint.boundaryField()[patchI];
}
// Calculate absolute flux from the mapped surface velocity
phi = mesh.Sf() & rhoUf();
if (correctPhi)
{
#include "correctPhi.H"
}
// Zero phi on current H-I
const surfaceScalarField faceMask
(
localMin<scalar>(mesh).interpolate(cellMask)
);
phi *= faceMask;
U *= cellMask;
// Make the fluxes relative to the mesh-motion
fvc::makeRelative(phi, rho, U);
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
if (pimple.firstIter())
{
#include "rhoEqn.H"
}
#include "UEqn.H"
#include "EEqn.H"
// --- Pressure corrector loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.turbCorr())
{
turbulence->correct();
}
}
rho = thermo.rho();
runTime.write();
runTime.printExecutionTime(Info);
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //
dimensionedScalar compressibility = fvc::domainIntegrate(psi);
bool compressible = (compressibility.value() > SMALL);
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());
mesh.interpolate(rAU);
surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p_rgh));
surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(rho*HbyA) + phig
);
if (ddtCorr)
{
surfaceScalarField faceMaskOld
(
localMin<scalar>(mesh).interpolate(cellMask.oldTime())
);
phiHbyA +=
faceMaskOld*MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi));
}
MRF.makeRelative(fvc::interpolate(rho), phiHbyA);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF);
fvScalarMatrix p_rghDDtEqn
(
fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
+ fvc::div(phiHbyA)
==
fvOptions(psi, p_rgh, rho.name())
);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
p_rghDDtEqn
- fvm::laplacian(rhorAUf, p_rgh)
);
p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
// Calculate the conservative fluxes
phi = phiHbyA + p_rghEqn.flux();
// Explicitly relax pressure for momentum corrector
p_rgh.relax();
// Correct the momentum source with the pressure gradient flux
// calculated from the relaxed pressure
U =
cellMask*
(
HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf)
);
U.correctBoundaryConditions();
fvOptions.correct(U);
K = 0.5*magSqr(U);
}
}
p = p_rgh + rho*gh;
#include "rhoEqn.H"
#include "compressibleContinuityErrs.H"
if (p_rgh.needReference())
{
if (!compressible)
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
}
else
{
p += (initialMass - fvc::domainIntegrate(psi*p))
/compressibility;
thermo.correctRho(psi*p - psip0);
rho = thermo.rho();
p_rgh = p - rho*gh;
}
}
else
{
thermo.correctRho(psi*p - psip0);
}
rho = thermo.rho();