solvers::incompressibleDenseParticleFluid: New solver module for particle laden incompressible flow

executed with foamRun for single region simulations of foamMultiRun for
multi-region simulations.  Replaces denseParticleFoam and all the corresponding
tutorials have been updated and moved to
tutorials/modules/incompressibleDenseParticleFluid.

Class
    Foam::solvers::incompressibleDenseParticleFluid

Description

    Solver module for transient flow of incompressible isothermal fluids coupled
    with particle clouds including the effect of the volume fraction of
    particles on the continuous phase, with optional mesh motion and change.

    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
    pseudo-transient and steady simulations.

    Optional fvModels and fvConstraints are provided to enhance the simulation
    in many ways including adding various sources, constraining or limiting
    the solution.
This commit is contained in:
Henry Weller
2023-04-11 16:56:13 +01:00
parent 288fc74ff0
commit e40198353b
82 changed files with 962 additions and 738 deletions

View File

@ -1,46 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2013-2022 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
CourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar CoNum = 0;
{
const scalarField sumPhi(fvc::surfaceSum(mag(phic))().primitiveField());
CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
const scalar meanCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
Info<< "Courant Number mean: " << meanCoNum
<< " max: " << CoNum << endl;
}
// ************************************************************************* //

View File

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

View File

@ -1,28 +0,0 @@
fvVectorMatrix UcEqn
(
fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
- fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
+ continuousPhaseTurbulence->divDevTau(Uc)
);
UcEqn.relax();
fvConstraints.constrain(UcEqn);
if (pimple.momentumPredictor())
{
solve
(
UcEqn
==
fvc::reconstruct
(
Fgf - fvc::snGrad(p)*mesh.magSf()
- Dcf*(phic - phid)
)
+ Dc*fvc::reconstruct(phic - phid)
+ Fd - fvm::Sp(Dc, Uc)
);
fvConstraints.constrain(Uc);
}

View File

@ -1,48 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2013-2020 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
continuityErrs
Description
Calculates and prints the continuity errors.
\*---------------------------------------------------------------------------*/
{
volScalarField contErr(fvc::ddt(alphac) + fvc::div(alphacf*phic));
scalar sumLocalContErr = runTime.deltaTValue()*
mag(contErr)().weightedAverage(mesh.V()).value();
scalar globalContErr = runTime.deltaTValue()*
contErr.weightedAverage(mesh.V()).value();
cumulativeContErr += globalContErr;
Info<< "time step continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;
}
// ************************************************************************* //

View File

@ -1,20 +0,0 @@
// Calculate absolute flux from the mapped surface velocity
phic = mesh.Sf() & Ucf();
correctUphiBCs(Uc, phic, true);
fv::correctPhi
(
phic,
Uc,
p,
autoPtr<volScalarField>(),
autoPtr<volScalarField>(),
pressureReference,
pimple
);
#include "continuityErrs.H"
// Make the flux relative to the mesh motion
fvc::makeRelative(phic, Uc);

View File

@ -1,162 +0,0 @@
#include "readGravitationalAcceleration.H"
word continuousPhaseName
(
IOdictionary
(
IOobject
(
"physicalProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ
)
).lookup("continuousPhaseName")
);
Info<< "Reading field U\n" << endl;
volVectorField Uc
(
IOobject
(
IOobject::groupName("U", continuousPhaseName),
runTime.name(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.name(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info<< "Reading/calculating continuous-phase face flux field phic\n"
<< endl;
surfaceScalarField phic
(
IOobject
(
IOobject::groupName("phi", continuousPhaseName),
runTime.name(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(Uc) & mesh.Sf()
);
pressureReference pressureReference(p, pimple.dict());
mesh.schemes().setFluxRequired(p.name());
Info<< "Creating turbulence model\n" << endl;
autoPtr<viscosityModel> continuousPhaseViscosity(viscosityModel::New(mesh));
dimensionedScalar rhocValue
(
IOobject::groupName("rho", continuousPhaseName),
dimDensity,
continuousPhaseViscosity->lookup
(
IOobject::groupName("rho", continuousPhaseName)
)
);
volScalarField rhoc
(
IOobject
(
rhocValue.name(),
runTime.name(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
rhocValue
);
volScalarField muc
(
IOobject
(
IOobject::groupName("mu", continuousPhaseName),
runTime.name(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
rhoc*continuousPhaseViscosity->nu()
);
Info << "Creating field alphac\n" << endl;
// alphac must be constructed before the cloud
// so that the drag-models can find it
volScalarField alphac
(
IOobject
(
IOobject::groupName("alpha", continuousPhaseName),
runTime.name(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimless, 0)
);
Info<< "Constructing clouds" << endl;
parcelClouds& clouds = parcelClouds::New(mesh, rhoc, Uc, muc, g);
// Particle fraction upper limit
scalar alphacMin
(
1 - mesh.solution().solverDict(alphac.name()).lookup<scalar>("max")
);
// Update alphac from the particle locations
alphac = max(1.0 - clouds.theta(), alphacMin);
alphac.correctBoundaryConditions();
surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
surfaceScalarField alphaPhic
(
IOobject::groupName
(
"alphaPhi",
continuousPhaseName
),
alphacf*phic
);
autoPtr<phaseIncompressible::momentumTransportModel> continuousPhaseTurbulence
(
phaseIncompressible::momentumTransportModel::New
(
alphac,
Uc,
alphaPhic,
phic,
continuousPhaseViscosity
)
);
#include "createFvModels.H"
#include "createFvConstraints.H"

View File

@ -1,208 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2013-2023 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
denseParticleFoam
Description
Transient solver for the coupled transport of particle clouds including the
effect of the volume fraction of particles on the continuous phase, with
optional mesh motion and mesh topology changes.
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "timeSelector.H"
#include "viscosityModel.H"
#include "phaseIncompressibleMomentumTransportModel.H"
#include "pimpleControl.H"
#include "pressureReference.H"
#include "findRefCell.H"
#include "constrainPressure.H"
#include "constrainHbyA.H"
#include "adjustPhi.H"
#include "uniformDimensionedFields.H"
#include "zeroGradientFvPatchFields.H"
#include "fvModels.H"
#include "fvConstraints.H"
#include "fvcDdt.H"
#include "fvcGrad.H"
#include "fvcSnGrad.H"
#include "fvcFlux.H"
#include "fvcMeshPhi.H"
#include "fvCorrectPhi.H"
#include "fvcReconstruct.H"
#include "fvmDdt.H"
#include "fvmDiv.H"
#include "fvmLaplacian.H"
#include "parcelClouds.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "postProcess.H"
#include "setRootCaseLists.H"
#include "createTime.H"
#include "createMesh.H"
#include "createDyMControls.H"
#include "createFields.H"
#include "createUcfIfPresent.H"
#include "initContinuityErrs.H"
Info<< "\nStarting time loop\n" << endl;
while (pimple.run(runTime))
{
#include "readDyMControls.H"
#include "CourantNo.H"
#include "setDeltaT.H"
// Update the mesh for topology change, mesh to mesh mapping
mesh.update();
runTime++;
Info<< "Time = " << runTime.userTimeName() << nl << endl;
// Store the particle positions
clouds.storeGlobalPositions();
// Move the mesh
mesh.move();
if (mesh.changing())
{
if (correctPhi)
{
#include "correctPhic.H"
}
if (checkMeshCourantNo)
{
#include "meshCourantNo.H"
}
}
continuousPhaseViscosity->correct();
muc = rhoc*continuousPhaseViscosity->nu();
clouds.evolve();
// Update continuous phase volume fraction field
alphac = max(1.0 - clouds.theta(), alphacMin);
alphac.correctBoundaryConditions();
alphacf = fvc::interpolate(alphac);
alphaPhic = alphacf*phic;
// Dispersed phase drag force
volVectorField Fd
(
IOobject
(
"Fd",
runTime.name(),
mesh
),
mesh,
dimensionedVector(dimAcceleration, Zero),
zeroGradientFvPatchVectorField::typeName
);
// continuous-dispersed phase drag coefficient
volScalarField Dc
(
IOobject
(
"Dc",
runTime.name(),
mesh
),
mesh,
dimensionedScalar(dimless/dimTime, Zero),
zeroGradientFvPatchVectorField::typeName
);
{
const fvVectorMatrix cloudSU(clouds.SU(Uc));
// Dispersed phase drag force
Fd.primitiveFieldRef() = -cloudSU.source()/mesh.V()/rhoc;
Fd.correctBoundaryConditions();
// Continuous phase drag coefficient
Dc.primitiveFieldRef() = -cloudSU.diag()/mesh.V()/rhoc;
Dc.correctBoundaryConditions();
}
// Face continuous-dispersed phase drag coefficient
const surfaceScalarField Dcf(fvc::interpolate(Dc));
// Effective flux of the dispersed phase
const surfaceScalarField phid
(
fvc::flux(Fd)/(Dcf + dimensionedScalar(Dc.dimensions(), small))
);
// Face buoyancy force
const surfaceScalarField Fgf(g & mesh.Sf());
// --- Pressure-velocity PIMPLE corrector loop
while (pimple.loop())
{
fvModels.correct();
#include "UcEqn.H"
// --- PISO loop
while (pimple.correct())
{
#include "pEqn.H"
}
if (pimple.correctTransport())
{
continuousPhaseTurbulence->correct();
}
}
runTime.write();
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,50 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 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
meshCourantNo
Description
Calculates and outputs the mean and maximum Courant Numbers.
\*---------------------------------------------------------------------------*/
scalar meshCoNum = 0.0;
scalar meanMeshCoNum = 0.0;
{
scalarField sumPhi
(
fvc::surfaceSum(mag(mesh.phi()))().primitiveField()
);
meshCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
meanMeshCoNum =
0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
}
Info<< "Mesh Courant Number mean: " << meanMeshCoNum
<< " max: " << meshCoNum << endl;
// ************************************************************************* //

View File

@ -1,89 +0,0 @@
{
const volScalarField rAUc(1.0/UcEqn.A());
const volScalarField r1ADUc(1/(1 + rAUc*Dc));
const surfaceScalarField rAUcf(fvc::interpolate(rAUc));
const surfaceScalarField r1ADUcf(1/(1 + rAUcf*Dcf));
const surfaceScalarField rADUcf("Dp", r1ADUcf*rAUcf);
volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
surfaceScalarField phiHbyAD
(
"phiHbyAD",
(
r1ADUcf
*(
fvc::flux(HbyA)
+ alphacf*rAUcf*fvc::ddtCorr(Uc, phic, Ucf)
)
)
);
if (p.needReference())
{
fvc::makeRelative(phiHbyAD, Uc);
adjustPhi(phiHbyAD, Uc, p);
fvc::makeAbsolute(phiHbyAD, Uc);
}
phiHbyAD += rADUcf*(Fgf + Dcf*phid);
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, Uc, phiHbyAD, rADUcf);
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(alphacf*rADUcf, p)
==
fvc::ddt(alphac)
+ fvc::div(alphacf*phiHbyAD)
);
pEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
pEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phic = phiHbyAD - pEqn.flux()/alphacf;
// Explicitly relax pressure for momentum corrector
p.relax();
Uc =
r1ADUc
*(
HbyA
+ rAUc
*(
fvc::reconstruct
(
Fgf - pEqn.flux()/alphacf/rADUcf
- Dcf*(phic - phid)
)
+ Dc*fvc::reconstruct(phic - phid)
+ Fd
)
);
Uc.correctBoundaryConditions();
fvConstraints.constrain(Uc);
// Correct Ucf if the mesh is moving
fvc::correctUf(Ucf, Uc, phic);
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phic, Uc);
}
}
}
#include "continuityErrs.H"

View File

@ -21,6 +21,7 @@ wmake $targetType multiphaseVoFSolver
wmake $targetType incompressibleMultiphaseVoF
wmake $targetType compressibleMultiphaseVoF
multiphaseEuler/Allwmake $targetType $*
wmake $targetType incompressibleDenseParticleFluid
isothermalFilm/Allwmake $targetType $*
film/Allwmake $targetType $*
wmake $targetType solid

View File

@ -0,0 +1,6 @@
moveMesh.C
momentumPredictor.C
correctPressure.C
incompressibleDenseParticleFluid.C
LIB = $(FOAM_LIBBIN)/libincompressibleDenseParticleFluid

View File

@ -1,26 +1,25 @@
EXE_INC = \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/parcel/lnInclude \
-I$(LIB_SRC)/physicalProperties/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-I$(LIB_SRC)/thermophysicalModels/multicomponentThermo/lnInclude \
-I$(FOAM_SOLVERS)/modules/fluidSolver/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
-I$(LIB_SRC)/MomentumTransportModels/phaseIncompressible/lnInclude \
-I$(LIB_SRC)/surfaceFilmModels/lnInclude \
-I$(LIB_SRC)/physicalProperties/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/lagrangian/parcel/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-llagrangian \
-llagrangianParcel \
-llagrangianParcelTurbulence \
-lphysicalProperties \
LIB_LIBS = \
-lfluidSolver \
-lmomentumTransportModels \
-lincompressibleMomentumTransportModels \
-lphaseIncompressibleMomentumTransportModels \
-lphysicalProperties \
-llagrangian \
-llagrangianParcel \
-llagrangianParcelTurbulence \
-lfiniteVolume \
-lfvModels \
-lfvConstraints \
-lsampling \
-lmeshTools

View File

@ -0,0 +1,138 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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/>.
\*---------------------------------------------------------------------------*/
#include "incompressibleDenseParticleFluid.H"
#include "constrainHbyA.H"
#include "constrainPressure.H"
#include "adjustPhi.H"
#include "fvcMeshPhi.H"
#include "fvcFlux.H"
#include "fvcDdt.H"
#include "fvcSnGrad.H"
#include "fvcReconstruct.H"
#include "fvmLaplacian.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::incompressibleDenseParticleFluid::correctPressure()
{
volScalarField& p(p_);
volVectorField& Uc(Uc_);
surfaceScalarField& phic(phic_);
fvVectorMatrix& UcEqn = tUcEqn.ref();
const volScalarField rAUc(1.0/UcEqn.A());
const volScalarField r1ADUc(1/(1 + rAUc*Dc()));
const surfaceScalarField rAUcf(fvc::interpolate(rAUc));
const surfaceScalarField r1ADUcf(1/(1 + rAUcf*Dcf()));
const surfaceScalarField rADUcf("Dp", r1ADUcf*rAUcf);
volVectorField HbyA(constrainHbyA(rAUc*UcEqn.H(), Uc, p));
surfaceScalarField phiHbyAD
(
"phiHbyAD",
(
r1ADUcf
*(
fvc::flux(HbyA)
+ alphacf*rAUcf*fvc::ddtCorr(Uc, phic, Ucf)
)
)
);
if (p.needReference())
{
fvc::makeRelative(phiHbyAD, Uc);
adjustPhi(phiHbyAD, Uc, p);
fvc::makeAbsolute(phiHbyAD, Uc);
}
// Face buoyancy force
const surfaceScalarField Fgf(g & mesh.Sf());
phiHbyAD += rADUcf*(Fgf + Dcf()*phid());
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, Uc, phiHbyAD, rADUcf);
// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvm::laplacian(alphacf*rADUcf, p)
==
fvc::ddt(alphac)
+ fvc::div(alphacf*phiHbyAD)
);
pEqn.setReference
(
pressureReference.refCell(),
pressureReference.refValue()
);
pEqn.solve();
if (pimple.finalNonOrthogonalIter())
{
phic = phiHbyAD - pEqn.flux()/alphacf;
// Explicitly relax pressure for momentum corrector
p.relax();
Uc =
r1ADUc
*(
HbyA
+ rAUc
*(
fvc::reconstruct
(
Fgf - pEqn.flux()/alphacf/rADUcf
- Dcf()*(phic - phid())
)
+ Dc()*fvc::reconstruct(phic - phid())
+ Fd()
)
);
Uc.correctBoundaryConditions();
fvConstraints().constrain(Uc);
// Correct Ucf if the mesh is moving
fvc::correctUf(Ucf, Uc, phic);
// Make the fluxes relative to the mesh motion
fvc::makeRelative(phic, Uc);
}
}
}
// ************************************************************************* //

View File

@ -0,0 +1,390 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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/>.
\*---------------------------------------------------------------------------*/
#include "incompressibleDenseParticleFluid.H"
#include "surfaceInterpolate.H"
#include "fvMatrix.H"
#include "fvcFlux.H"
#include "zeroGradientFvPatchFields.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace solvers
{
defineTypeNameAndDebug(incompressibleDenseParticleFluid, 0);
addToRunTimeSelectionTable
(
solver,
incompressibleDenseParticleFluid,
fvMesh
);
}
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::solvers::incompressibleDenseParticleFluid::correctCoNum()
{
fluidSolver::correctCoNum(phic);
}
void Foam::solvers::incompressibleDenseParticleFluid::continuityErrors()
{
fluidSolver::continuityErrors(phic);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::solvers::incompressibleDenseParticleFluid::
incompressibleDenseParticleFluid
(
fvMesh& mesh
)
:
fluidSolver(mesh),
continuousPhaseName
(
IOdictionary
(
IOobject
(
"physicalProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ
)
).lookup("continuousPhaseName")
),
p_
(
IOobject
(
"p",
runTime.name(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
pressureReference(p_, pimple.dict()),
g
(
IOobject
(
"g",
runTime.constant(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
Uc_
(
IOobject
(
IOobject::groupName("U", continuousPhaseName),
runTime.name(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
phic_
(
IOobject
(
IOobject::groupName("phi", continuousPhaseName),
runTime.name(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(Uc_) & mesh.Sf()
),
viscosity(viscosityModel::New(mesh)),
rhoc
(
IOobject
(
IOobject::groupName("rho", continuousPhaseName),
runTime.name(),
mesh
),
mesh,
dimensionedScalar
(
IOobject::groupName("rho", continuousPhaseName),
dimDensity,
viscosity->lookup
(
IOobject::groupName("rho", continuousPhaseName)
)
)
),
muc
(
IOobject
(
IOobject::groupName("mu", continuousPhaseName),
runTime.name(),
mesh
),
rhoc*viscosity->nu()
),
alphac_
(
IOobject
(
IOobject::groupName("alpha", continuousPhaseName),
runTime.name(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar(dimless, 0)
),
alphacMin
(
1 - mesh.solution().solverDict(alphac_.name()).lookup<scalar>("max")
),
alphacf("alphacf", fvc::interpolate(alphac_)),
alphaPhic
(
IOobject::groupName("alphaPhi", continuousPhaseName),
alphacf*phic_
),
momentumTransport
(
phaseIncompressible::momentumTransportModel::New
(
alphac_,
Uc_,
alphaPhic,
phic_,
viscosity
)
),
clouds
(
parcelClouds::New(mesh, rhoc, Uc_, muc, g)
),
p(p_),
Uc(Uc_),
phic(phic_),
alphac(alphac_)
{
mesh.schemes().setFluxRequired(p.name());
momentumTransport->validate();
// Update alphac from the particle locations
alphac_ = max(1 - clouds.theta(), alphacMin);
alphac_.correctBoundaryConditions();
alphacf = fvc::interpolate(alphac);
alphaPhic = alphacf*phic;
if (mesh.dynamic())
{
Info<< "Constructing face momentum Ucf" << endl;
// Ensure the U BCs are up-to-date before constructing Ucf
Uc_.correctBoundaryConditions();
Ucf = new surfaceVectorField
(
IOobject
(
"Ucf",
runTime.name(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(Uc)
);
}
correctCoNum();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::solvers::incompressibleDenseParticleFluid::
~incompressibleDenseParticleFluid()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::incompressibleDenseParticleFluid::preSolve()
{
// Read the controls
readControls();
// Store the particle positions
if (mesh.topoChanging())
{
clouds.storeGlobalPositions();
}
fvModels().preUpdateMesh();
correctCoNum();
// Update the mesh for topology change, mesh to mesh mapping
mesh_.update();
}
void Foam::solvers::incompressibleDenseParticleFluid::prePredictor()
{
if (pimple.firstIter())
{
clouds.evolve();
// Update continuous phase volume fraction field
alphac_ = max(1 - clouds.theta(), alphacMin);
alphac_.correctBoundaryConditions();
alphacf = fvc::interpolate(alphac);
// ... and continuous phase volumetric flux
alphaPhic = alphacf*phic;
// Update the continuous phase dynamic viscosity
muc = rhoc*viscosity->nu();
Fd = new volVectorField
(
IOobject
(
"Fd",
runTime.name(),
mesh
),
mesh,
dimensionedVector(dimAcceleration, Zero),
zeroGradientFvPatchVectorField::typeName
);
Dc = new volScalarField
(
IOobject
(
"Dc",
runTime.name(),
mesh
),
mesh,
dimensionedScalar(dimless/dimTime, Zero),
zeroGradientFvPatchVectorField::typeName
);
const fvVectorMatrix cloudSU(clouds.SU(Uc));
Fd().primitiveFieldRef() = -cloudSU.source()/mesh.V()/rhoc;
Fd().correctBoundaryConditions();
Dc().primitiveFieldRef() = -cloudSU.diag()/mesh.V()/rhoc;
Dc().correctBoundaryConditions();
Dcf = fvc::interpolate(Dc()).ptr();
phid =
(
fvc::flux(Fd())
/(Dcf() + dimensionedScalar(Dc().dimensions(), small))
).ptr();
}
fvModels().correct();
if (pimple.predictTransport())
{
momentumTransport->predict();
}
}
void Foam::solvers::incompressibleDenseParticleFluid::thermophysicalPredictor()
{}
void Foam::solvers::incompressibleDenseParticleFluid::pressureCorrector()
{
while (pimple.correct())
{
correctPressure();
}
tUcEqn.clear();
}
void Foam::solvers::incompressibleDenseParticleFluid::postCorrector()
{
if (pimple.correctTransport())
{
viscosity->correct();
momentumTransport->correct();
}
}
void Foam::solvers::incompressibleDenseParticleFluid::postSolve()
{
Fd.clear();
Dc.clear();
Dcf.clear();
phid.clear();
}
// ************************************************************************* //

View File

@ -0,0 +1,265 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2023 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/>.
Class
Foam::solvers::incompressibleDenseParticleFluid
Description
Solver module for transient flow of incompressible isothermal fluids coupled
with particle clouds including the effect of the volume fraction of
particles on the continuous phase, with optional mesh motion and change.
Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
pseudo-transient and steady simulations.
Optional fvModels and fvConstraints are provided to enhance the simulation
in many ways including adding various sources, constraining or limiting
the solution.
Reference:
\verbatim
Greenshields, C. J., & Weller, H. G. (2022).
Notes on Computational Fluid Dynamics: General Principles.
CFD Direct Ltd.: Reading, UK.
\endverbatim
SourceFiles
incompressibleDenseParticleFluid.C
See also
Foam::solvers::fluidSolver
Foam::solvers::incompressibleFluid
\*---------------------------------------------------------------------------*/
#ifndef incompressibleDenseParticleFluid_H
#define incompressibleDenseParticleFluid_H
#include "fluidSolver.H"
#include "viscosityModel.H"
#include "phaseIncompressibleMomentumTransportModel.H"
#include "pressureReference.H"
#include "uniformDimensionedFields.H"
#include "parcelClouds.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace solvers
{
/*---------------------------------------------------------------------------*\
Class incompressibleDenseParticleFluid Declaration
\*---------------------------------------------------------------------------*/
class incompressibleDenseParticleFluid
:
public fluidSolver
{
protected:
//- Name of the continuous phase
const word continuousPhaseName;
// Pressure
//- Pressure field
volScalarField p_;
//- Pressure reference
Foam::pressureReference pressureReference;
// Buoyancy
//- Acceleration due to gravity
uniformDimensionedVectorField g;
// Kinematic properties
//- Continuous phase velocity field
volVectorField Uc_;
//- Continuous phase flux field
surfaceScalarField phic_;
//- Continuous phase kinematic viscosity model
autoPtr<viscosityModel> viscosity;
// Continuous phase properties for Lagrangian cloud
//- Continuous phase density
volScalarField rhoc;
//- Continuous phase viscosity
volScalarField muc;
//- Continuous phase-fraction
volScalarField alphac_;
//- Minimum continuous phase-fraction
scalar alphacMin;
//- Interpolated continuous phase-fraction
surfaceScalarField alphacf;
//- Continuous phase volumetric-flux field
surfaceScalarField alphaPhic;
// Momentum transport
//- Pointer to the momentum transport model
autoPtr<phaseIncompressible::momentumTransportModel> momentumTransport;
// Dispersed phase Lagrangian clouds
parcelClouds& clouds;
// Cached temporary fields
//- Dispersed phase drag force
autoPtr<volVectorField> Fd;
// Continuous-dispersed phase drag coefficient
autoPtr<volScalarField> Dc;
//- Face continuous-dispersed phase drag coefficient
autoPtr<surfaceScalarField> Dcf;
//- Effective volumetric flux of the dispersed phase
autoPtr<surfaceScalarField> phid;
//- Pointer to the surface momentum field
// used to recreate the flux after mesh-change
autoPtr<surfaceVectorField> Ucf;
//- Cached momentum matrix
// shared between the momentum predictor and pressure corrector
tmp<fvVectorMatrix> tUcEqn;
// Protected Member Functions
//- Correct the cached Courant numbers
void correctCoNum();
//- Calculate and print the continuity errors
void continuityErrors();
//- Construct the pressure equation
// and correct the pressure and velocity
virtual void correctPressure();
public:
// Public Data
//- Reference to the pressure field
const volScalarField& p;
//- Reference to the continuous phase velocity field
const volVectorField& Uc;
//- Reference to the continuous phase volumetric-flux field
const surfaceScalarField& phic;
//- Reference continuous phase-fraction
const volScalarField& alphac;
//- Runtime type information
TypeName("incompressibleDenseParticleFluid");
// Constructors
//- Construct from region mesh
incompressibleDenseParticleFluid(fvMesh& mesh);
//- Disallow default bitwise copy construction
incompressibleDenseParticleFluid
(
const incompressibleDenseParticleFluid&
) = delete;
//- Destructor
virtual ~incompressibleDenseParticleFluid();
// Member Functions
//- Called at the start of the time-step, before the PIMPLE loop
virtual void preSolve();
//- Called at the start of the PIMPLE loop to move the mesh
virtual void moveMesh();
//- Called at the start of the PIMPLE loop
virtual void prePredictor();
//- Construct and optionally solve the momentum equation
virtual void momentumPredictor();
//- Construct and solve the energy equation,
// convert to temperature
// and update thermophysical and transport properties
virtual void thermophysicalPredictor();
//- Construct and solve the pressure equation in the PISO loop
virtual void pressureCorrector();
//- Correct the momentum and thermophysical transport modelling
virtual void postCorrector();
//- Called after the PIMPLE loop at the end of the time-step
virtual void postSolve();
// Member Operators
//- Disallow default bitwise assignment
void operator=(const incompressibleDenseParticleFluid&) = delete;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace solvers
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2017-2022 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,35 +21,56 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Global
createUcfIfPresent
Description
Creates and initialises the continuous phase face velocity field Ufc
if required.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "incompressibleDenseParticleFluid.H"
#include "fvcDdt.H"
#include "fvcSnGrad.H"
#include "fvcReconstruct.H"
#include "fvmDiv.H"
#include "fvmDdt.H"
autoPtr<surfaceVectorField> Ucf;
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
if (mesh.dynamic())
void Foam::solvers::incompressibleDenseParticleFluid::momentumPredictor()
{
Info<< "Constructing continuous phase face velocity Ucf\n" << endl;
volVectorField& Uc(Uc_);
Ucf = new surfaceVectorField
tUcEqn =
(
IOobject
(
"Ucf",
runTime.name(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(Uc)
fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
- fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
+ momentumTransport->divDevSigma(Uc)
==
fvModels().source(Uc)
);
fvVectorMatrix& UcEqn = tUcEqn.ref();
UcEqn.relax();
fvConstraints().constrain(UcEqn);
if (pimple.momentumPredictor())
{
// Face buoyancy force
const surfaceScalarField Fgf(g & mesh.Sf());
solve
(
UcEqn
==
fvc::reconstruct
(
Fgf - fvc::snGrad(p)*mesh.magSf()
- Dcf()*(phic - phid())
)
+ Dc()*fvc::reconstruct(phic - phid())
+ Fd() - fvm::Sp(Dc(), Uc)
);
fvConstraints().constrain(Uc);
}
}
// ************************************************************************* //

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2017-2022 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -21,29 +21,51 @@ License
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Global
createUcf
Description
Creates and initialises the velocity velocity field Ucf.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "incompressibleDenseParticleFluid.H"
#include "fvCorrectPhi.H"
#include "fvcMeshPhi.H"
#include "geometricZeroField.H"
Info<< "Reading/calculating continuous phase face velocity Ucf\n" << endl;
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::solvers::incompressibleDenseParticleFluid::moveMesh()
{
if (pimple.firstIter() || pimple.moveMeshOuterCorrectors())
{
// Move the mesh
mesh_.move();
if (mesh.changing())
{
if (correctPhi || mesh.topoChanged())
{
// Calculate absolute flux
// from the mapped surface velocity
phic_ = mesh.Sf() & Ucf();
correctUphiBCs(Uc_, phic_, true);
fv::correctPhi
(
phic_,
Uc,
p,
autoPtr<volScalarField>(),
autoPtr<volScalarField>(),
pressureReference,
pimple
);
// Make the flux relative to the mesh motion
fvc::makeRelative(phic_, Uc);
}
meshCourantNo();
}
}
}
surfaceVectorField Ucf
(
IOobject
(
"Ucf",
runTime.name(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
fvc::interpolate(Uc)
);
// ************************************************************************* //

View File

@ -92,7 +92,7 @@ protected:
//- Velocity field
volVectorField U_;
//- Mass-flux field
//- Volumetric-flux field
surfaceScalarField phi_;
@ -151,7 +151,7 @@ public:
//- Reference to the velocity field
const volVectorField& U;
//- Reference to the mass-flux field
//- Reference to the volumetric-flux field
const surfaceScalarField& phi;

46
bin/denseParticleFoam Executable file
View File

@ -0,0 +1,46 @@
#!/bin/sh
#------------------------------------------------------------------------------
# ========= |
# \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
# \\ / O peration | Website: https://openfoam.org
# \\ / A nd | Copyright (C) 2023 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/>.
#
# Script
# denseParticleFoam
#
# Description
# Script to inform the user that denseParticleFoam has been superseded and
# replaced by the more general incompressibleDenseParticleFluid solver
# module executed by the foamRun application.
#
#------------------------------------------------------------------------------
cat <<EOF
denseParticleFoam has been superseded and replaced by the more general
incompressibleDenseParticleFluid solver module executed by the foamRun application:
foamRun -solver incompressibleDenseParticleFluid
EOF
exec env foamRun -solver incompressibleDenseParticleFluid "$@"
#------------------------------------------------------------------------------

View File

@ -14,7 +14,9 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application denseParticleFoam;
application foamRun;
solver incompressibleDenseParticleFluid;
startFrom startTime;

View File

@ -64,9 +64,5 @@ PIMPLE
pRefValue 0;
}
relaxationFactors
{
}
// ************************************************************************* //

View File

@ -14,7 +14,9 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application denseParticleFoam;
application foamRun;
solver incompressibleDenseParticleFluid;
startFrom startTime;

View File

@ -72,9 +72,5 @@ PIMPLE
pRefValue 0;
}
relaxationFactors
{
}
// ************************************************************************* //

View File

@ -14,7 +14,9 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application denseParticleFoam;
application foamRun;
solver incompressibleDenseParticleFluid;
startFrom startTime;

View File

@ -51,7 +51,7 @@ solvers
relTol 0;
}
cloud:alpha
"cloud:alpha.*"
{
solver smoothSolver;
smoother symGaussSeidel;
@ -70,9 +70,5 @@ PIMPLE
pRefValue 0;
}
relaxationFactors
{
}
// ************************************************************************* //

View File

@ -14,7 +14,9 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application denseParticleFoam;
application foamRun;
solver incompressibleDenseParticleFluid;
startFrom latestTime;

View File

@ -53,7 +53,7 @@ solvers
relTol 0.1;
}
cloud:alpha
"cloud:alpha.*"
{
solver GAMG;
tolerance 1e-06;
@ -72,9 +72,5 @@ PIMPLE
pRefValue 0;
}
relaxationFactors
{
}
// ************************************************************************* //

View File

@ -14,7 +14,9 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application denseParticleFoam;
application foamRun;
solver incompressibleDenseParticleFluid;
startFrom startTime;

View File

@ -63,9 +63,5 @@ PIMPLE
pRefValue 0;
}
relaxationFactors
{
}
// ************************************************************************* //