diff --git a/applications/solvers/lagrangian/DPMFoam/Allwclean b/applications/solvers/lagrangian/DPMFoam/Allwclean new file mode 100755 index 0000000000..25d0b2955f --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/Allwclean @@ -0,0 +1,8 @@ +#!/bin/sh + +cd ${0%/*} || exit 1 +set -x + +wclean DPMTurbulenceModels +wclean +wclean MPPICFoam diff --git a/applications/solvers/lagrangian/DPMFoam/Allwmake b/applications/solvers/lagrangian/DPMFoam/Allwmake new file mode 100755 index 0000000000..6308a7052b --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/Allwmake @@ -0,0 +1,8 @@ +#!/bin/sh + +cd ${0%/*} || exit 1 +set -x + +wmake DPMTurbulenceModels +wmake +wmake MPPICFoam diff --git a/applications/solvers/lagrangian/DPMFoam/CourantNo.H b/applications/solvers/lagrangian/DPMFoam/CourantNo.H new file mode 100644 index 0000000000..9ff53ed401 --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/CourantNo.H @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Global + CourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers. + +\*---------------------------------------------------------------------------*/ + +scalar CoNum = 0.0; +scalar meanCoNum = 0.0; + +if (mesh.nInternalFaces()) +{ + scalarField sumPhi + ( + fvc::surfaceSum(mag(phic))().internalField() + ); + + CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + + meanCoNum = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); +} + +Info<< "Courant Number mean: " << meanCoNum + << " max: " << CoNum << endl; + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/DPMFoam/DPMFoam.C b/applications/solvers/lagrangian/DPMFoam/DPMFoam.C new file mode 100644 index 0000000000..9c97968c1c --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/DPMFoam.C @@ -0,0 +1,143 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Application + DPMFoam + +Description + Transient solver for the coupled transport of a single kinematic particle + could including the effect of the volume fraction of particles on the + continuous phase. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "singlePhaseTransportModel.H" +#include "PhaseIncompressibleTurbulenceModel.H" +#include "pimpleControl.H" +#include "fixedFluxPressureFvPatchScalarField.H" + +#ifdef MPPIC + #include "basicKinematicMPPICCloud.H" + #define basicKinematicTypeCloud basicKinematicMPPICCloud +#else + #include "basicKinematicCollidingCloud.H" + #define basicKinematicTypeCloud basicKinematicCollidingCloud +#endif + +int main(int argc, char *argv[]) +{ + argList::addOption + ( + "cloudName", + "name", + "specify alternative cloud name. default is 'kinematicCloud'" + ); + + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readGravitationalAcceleration.H" + #include "createFields.H" + #include "initContinuityErrs.H" + + pimpleControl pimple(mesh); + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + #include "CourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + continuousPhaseTransport.correct(); + muc = rhoc*continuousPhaseTransport.nu(); + + Info<< "Evolving " << kinematicCloud.name() << endl; + kinematicCloud.evolve(); + + // Update continuous phase volume fraction field + alphac = max(1.0 - kinematicCloud.theta(), alphacMin); + alphac.correctBoundaryConditions(); + alphacf = fvc::interpolate(alphac); + alphaPhic = alphacf*phic; + + fvVectorMatrix cloudSU(kinematicCloud.SU(Uc)); + volVectorField cloudVolSUSu + ( + IOobject + ( + "cloudVolSUSu", + runTime.timeName(), + mesh + ), + mesh, + dimensionedVector + ( + "0", + cloudSU.dimensions()/dimVolume, + vector::zero + ), + zeroGradientFvPatchVectorField::typeName + ); + + cloudVolSUSu.internalField() = - cloudSU.source()/mesh.V(); + cloudVolSUSu.correctBoundaryConditions(); + cloudSU.source() = vector::zero; + + // --- Pressure-velocity PIMPLE corrector loop + while (pimple.loop()) + { + #include "UcEqn.H" + + // --- PISO loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + continuousPhaseTurbulence->correct(); + } + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMTurbulenceModels.C b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMTurbulenceModels.C new file mode 100644 index 0000000000..1aeca77957 --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMTurbulenceModels.C @@ -0,0 +1,61 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "PhaseIncompressibleTurbulenceModel.H" +#include "singlePhaseTransportModel.H" +#include "addToRunTimeSelectionTable.H" +#include "makeTurbulenceModel.H" + +#include "laminar.H" +#include "RASModel.H" +#include "LESModel.H" + +makeBaseTurbulenceModel +( + volScalarField, + geometricOneField, + incompressibleTurbulenceModel, + PhaseIncompressibleTurbulenceModel, + singlePhaseTransportModel +); + +#define makeRASModel(Type) \ + makeTemplatedTurbulenceModel \ + (singlePhaseTransportModelPhaseIncompressibleTurbulenceModel, RAS, Type) + +#define makeLESModel(Type) \ + makeTemplatedTurbulenceModel \ + (singlePhaseTransportModelPhaseIncompressibleTurbulenceModel, LES, Type) + +#include "kEpsilon.H" +makeRASModel(kEpsilon); + +#include "Smagorinsky.H" +makeLESModel(Smagorinsky); + +#include "kEqn.H" +makeLESModel(kEqn); + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/files b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/files new file mode 100644 index 0000000000..3b6b48c2c4 --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/files @@ -0,0 +1,3 @@ +DPMTurbulenceModels.C + +LIB = $(FOAM_LIBBIN)/libDPMTurbulenceModels diff --git a/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options new file mode 100644 index 0000000000..716929b562 --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -I$(LIB_SRC)/foam/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude diff --git a/applications/solvers/lagrangian/DPMFoam/MPPICFoam/MPPICFoam.C b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/MPPICFoam.C new file mode 100644 index 0000000000..5a036f8f50 --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/MPPICFoam.C @@ -0,0 +1,40 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +Application + MPPICFoam + +Description + Transient solver for the coupled transport of a single kinematic particle + could including the effect of the volume fraction of particles on the + continuous phase. Multi-Phase Particle In Cell (MPPIC) modeling is used to + represent collisions without resolving particle-particle interactions. + +\*---------------------------------------------------------------------------*/ + +#define MPPIC + +#include "DPMFoam.C" + + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/files b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/files new file mode 100644 index 0000000000..d311f5a33c --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/files @@ -0,0 +1,3 @@ +MPPICFoam.C + +EXE = $(FOAM_APPBIN)/MPPICFoam diff --git a/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options new file mode 100644 index 0000000000..948df09838 --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options @@ -0,0 +1,34 @@ +EXE_INC = \ + -I.. \ + -I../DPMTurbulenceModels/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \ + -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ + -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -llagrangian \ + -llagrangianIntermediate \ + -lthermophysicalFunctions \ + -lspecie \ + -lradiationModels \ + -lincompressibleTransportModels \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lDPMTurbulenceModels \ + -lregionModels \ + -lsurfaceFilmModels \ + -lsampling diff --git a/applications/solvers/lagrangian/DPMFoam/Make/files b/applications/solvers/lagrangian/DPMFoam/Make/files new file mode 100644 index 0000000000..c4d729205d --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/Make/files @@ -0,0 +1,3 @@ +DPMFoam.C + +EXE = $(FOAM_APPBIN)/DPMFoam diff --git a/applications/solvers/lagrangian/DPMFoam/Make/options b/applications/solvers/lagrangian/DPMFoam/Make/options new file mode 100644 index 0000000000..8a4d6cc23b --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/Make/options @@ -0,0 +1,33 @@ +EXE_INC = \ + -I./DPMTurbulenceModels/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \ + -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ + -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -llagrangian \ + -llagrangianIntermediate \ + -lthermophysicalFunctions \ + -lspecie \ + -lradiationModels \ + -lincompressibleTransportModels \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lDPMTurbulenceModels \ + -lregionModels \ + -lsurfaceFilmModels \ + -lsampling diff --git a/applications/solvers/lagrangian/DPMFoam/UcEqn.H b/applications/solvers/lagrangian/DPMFoam/UcEqn.H new file mode 100644 index 0000000000..0bd4a7146d --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/UcEqn.H @@ -0,0 +1,32 @@ +fvVectorMatrix UcEqn +( + fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc) + - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc) + + continuousPhaseTurbulence->divDevRhoReff(Uc) + == + (1.0/rhoc)*cloudSU +); + +UcEqn.relax(); + +volScalarField rAUc(1.0/UcEqn.A()); +surfaceScalarField rAUcf("Dp", fvc::interpolate(rAUc)); + +surfaceScalarField phicForces +( + (fvc::interpolate(rAUc*cloudVolSUSu/rhoc) & mesh.Sf()) + + rAUcf*(g & mesh.Sf()) +); + +if (pimple.momentumPredictor()) +{ + solve + ( + UcEqn + == + fvc::reconstruct + ( + phicForces/rAUcf - fvc::snGrad(p)*mesh.magSf() + ) + ); +} diff --git a/applications/solvers/lagrangian/DPMFoam/continuityErrs.H b/applications/solvers/lagrangian/DPMFoam/continuityErrs.H new file mode 100644 index 0000000000..f07d42b39a --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/continuityErrs.H @@ -0,0 +1,48 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +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; +} + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/DPMFoam/createFields.H b/applications/solvers/lagrangian/DPMFoam/createFields.H new file mode 100644 index 0000000000..16ed9fa919 --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/createFields.H @@ -0,0 +1,166 @@ + Info<< "\nReading transportProperties\n" << endl; + + IOdictionary transportProperties + ( + IOobject + ( + "transportProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + IOobject::NO_WRITE + ) + ); + + word contiuousPhaseName(transportProperties.lookup("contiuousPhaseName")); + + dimensionedScalar rhocValue + ( + IOobject::groupName("rho", contiuousPhaseName), + dimDensity, + transportProperties.lookup + ( + IOobject::groupName("rho", contiuousPhaseName) + ) + ); + + volScalarField rhoc + ( + IOobject + ( + rhocValue.name(), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + rhocValue + ); + + Info<< "Reading field U\n" << endl; + volVectorField Uc + ( + IOobject + ( + IOobject::groupName("U", contiuousPhaseName), + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "Reading field p\n" << endl; + volScalarField p + ( + IOobject + ( + "p", + runTime.timeName(), + 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", contiuousPhaseName), + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + linearInterpolate(Uc) & mesh.Sf() + ); + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue); + + Info<< "Creating turbulence model\n" << endl; + + singlePhaseTransportModel continuousPhaseTransport(Uc, phic); + + volScalarField muc + ( + IOobject + ( + IOobject::groupName("mu", contiuousPhaseName), + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + rhoc*continuousPhaseTransport.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", contiuousPhaseName), + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("0", dimless, 0) + ); + + word kinematicCloudName("kinematicCloud"); + args.optionReadIfPresent("cloudName", kinematicCloudName); + + Info<< "Constructing kinematicCloud " << kinematicCloudName << endl; + basicKinematicTypeCloud kinematicCloud + ( + kinematicCloudName, + rhoc, + Uc, + muc, + g + ); + + // Particle fraction upper limit + scalar alphacMin + ( + 1.0 + - readScalar + ( + kinematicCloud.particleProperties().subDict("constantProperties") + .lookup("alphaMax") + ) + ); + + // Update alphac from the particle locations + alphac = max(1.0 - kinematicCloud.theta(), alphacMin); + alphac.correctBoundaryConditions(); + + surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac)); + surfaceScalarField alphaPhic("alphaPhic", alphacf*phic); + + autoPtr > + continuousPhaseTurbulence + ( + PhaseIncompressibleTurbulenceModel::New + ( + alphac, + Uc, + alphaPhic, + phic, + continuousPhaseTransport + ) + ); diff --git a/applications/solvers/lagrangian/DPMFoam/pEqn.H b/applications/solvers/lagrangian/DPMFoam/pEqn.H new file mode 100644 index 0000000000..34cce4f166 --- /dev/null +++ b/applications/solvers/lagrangian/DPMFoam/pEqn.H @@ -0,0 +1,52 @@ +{ + volVectorField HbyA("HbyA", Uc); + HbyA = rAUc*UcEqn.H(); + + surfaceScalarField phiHbyA + ( + "phiHbyA", + ( + (fvc::interpolate(HbyA) & mesh.Sf()) + + alphacf*rAUcf*fvc::ddtCorr(Uc, phic) + + phicForces + ) + ); + + // Update the fixedFluxPressure BCs to ensure flux consistency + setSnGrad + ( + p.boundaryField(), + ( + phiHbyA.boundaryField() + - (mesh.Sf().boundaryField() & Uc.boundaryField()) + )/(mesh.magSf().boundaryField()*rAUcf.boundaryField()) + ); + + // Non-orthogonal pressure corrector loop + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::laplacian(alphacf*rAUcf, p) + == + fvc::ddt(alphac) + fvc::div(alphacf*phiHbyA) + ); + + pEqn.setReference(pRefCell, pRefValue); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phic = phiHbyA - pEqn.flux()/alphacf; + + p.relax(); + + Uc = HbyA + + rAUc*fvc::reconstruct((phicForces - pEqn.flux()/alphacf)/rAUcf); + Uc.correctBoundaryConditions(); + } + } +} + +#include "continuityErrs.H" diff --git a/src/lagrangian/intermediate/Make/files b/src/lagrangian/intermediate/Make/files index 6d4e26a441..9cadf55117 100644 --- a/src/lagrangian/intermediate/Make/files +++ b/src/lagrangian/intermediate/Make/files @@ -19,6 +19,8 @@ KINEMATICPARCEL=$(DERIVEDPARCELS)/basicKinematicParcel $(KINEMATICPARCEL)/defineBasicKinematicParcel.C $(KINEMATICPARCEL)/makeBasicKinematicParcelSubmodels.C + +/* kinematic colliding parcel sub-models */ KINEMATICCOLLIDINGPARCEL=$(DERIVEDPARCELS)/basicKinematicCollidingParcel $(KINEMATICCOLLIDINGPARCEL)/defineBasicKinematicCollidingParcel.C $(KINEMATICCOLLIDINGPARCEL)/makeBasicKinematicCollidingParcelSubmodels.C @@ -42,6 +44,12 @@ $(REACTINGMPPARCEL)/defineBasicReactingMultiphaseParcel.C $(REACTINGMPPARCEL)/makeBasicReactingMultiphaseParcelSubmodels.C +/* kinematic MPPIC parcel sub-models */ +KINEMATICMPPICPARCEL=$(DERIVEDPARCELS)/basicKinematicMPPICParcel +$(KINEMATICMPPICPARCEL)/defineBasicKinematicMPPICParcel.C +$(KINEMATICMPPICPARCEL)/makeBasicKinematicMPPICParcelSubmodels.C + + /* bolt-on models */ RADIATION=submodels/addOns/radiation $(RADIATION)/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C @@ -71,6 +79,24 @@ $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphase $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphaseParcelInjectionDataIO.C $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphaseParcelInjectionDataIOList.C +MPPICPARTICLESTRESS=submodels/MPPIC/ParticleStressModels +$(MPPICPARTICLESTRESS)/ParticleStressModel/ParticleStressModel.C +$(MPPICPARTICLESTRESS)/HarrisCrighton/HarrisCrighton.C +$(MPPICPARTICLESTRESS)/Lun/Lun.C +$(MPPICPARTICLESTRESS)/exponential/exponential.C + +MPPICCORRECTIONLIMITING=submodels/MPPIC/CorrectionLimitingMethods +$(MPPICCORRECTIONLIMITING)/CorrectionLimitingMethod/CorrectionLimitingMethod.C +$(MPPICCORRECTIONLIMITING)/noCorrectionLimiting/noCorrectionLimiting.C +$(MPPICCORRECTIONLIMITING)/absolute/absolute.C +$(MPPICCORRECTIONLIMITING)/relative/relative.C + +MPPICTIMESCALE=submodels/MPPIC/TimeScaleModels +$(MPPICTIMESCALE)/TimeScaleModel/TimeScaleModel.C +$(MPPICTIMESCALE)/equilibrium/equilibrium.C +$(MPPICTIMESCALE)/nonEquilibrium/nonEquilibrium.C +$(MPPICTIMESCALE)/isotropic/isotropic.C + /* integration schemes */ IntegrationScheme/makeIntegrationSchemes.C @@ -81,8 +107,13 @@ phaseProperties/phaseProperties/phaseProperties.C phaseProperties/phaseProperties/phasePropertiesIO.C phaseProperties/phasePropertiesList/phasePropertiesList.C -/* Additional helper classes */ + +/* additional helper classes */ clouds/Templates/KinematicCloud/cloudSolution/cloudSolution.C +/* averaging methods */ +submodels/MPPIC/AveragingMethods/makeAveragingMethods.C + + LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate diff --git a/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C new file mode 100644 index 0000000000..f22cac7b0f --- /dev/null +++ b/src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C @@ -0,0 +1,299 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 . + +\*---------------------------------------------------------------------------*/ + +#include "MPPICCloud.H" +#include "PackingModel.H" +#include "ParticleStressModel.H" +#include "DampingModel.H" +#include "IsotropyModel.H" +#include "TimeScaleModel.H" + +// * * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * // + +template +void Foam::MPPICCloud::setModels() +{ + packingModel_.reset + ( + PackingModel >::New + ( + this->subModelProperties(), + *this + ).ptr() + ); + dampingModel_.reset + ( + DampingModel >::New + ( + this->subModelProperties(), + *this + ).ptr() + ); + isotropyModel_.reset + ( + IsotropyModel >::New + ( + this->subModelProperties(), + *this + ).ptr() + ); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::MPPICCloud::MPPICCloud +( + const word& cloudName, + const volScalarField& rho, + const volVectorField& U, + const volScalarField& mu, + const dimensionedVector& g, + bool readFields +) +: + CloudType(cloudName, rho, U, mu, g, false), + packingModel_(NULL), + dampingModel_(NULL), + isotropyModel_(NULL) +{ + if (this->solution().steadyState()) + { + FatalErrorIn + ( + "Foam::MPPICCloud::MPPICCloud" + "(" + "const word&, " + "const volScalarField&, " + "const volVectorField&, " + "const volScalarField&, " + "const dimensionedVector&, " + "bool" + ")" + ) << "MPPIC modelling not available for steady state calculations" + << exit(FatalError); + } + + if (this->solution().active()) + { + setModels(); + + if (readFields) + { + parcelType::readFields(*this); + } + } +} + + +template +Foam::MPPICCloud::MPPICCloud +( + MPPICCloud& c, + const word& name +) +: + CloudType(c, name), + packingModel_(c.packingModel_->clone()), + dampingModel_(c.dampingModel_->clone()), + isotropyModel_(c.isotropyModel_->clone()) +{} + + +template +Foam::MPPICCloud::MPPICCloud +( + const fvMesh& mesh, + const word& name, + const MPPICCloud& c +) +: + CloudType(mesh, name, c), + packingModel_(NULL), + dampingModel_(NULL), + isotropyModel_(NULL) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +template +Foam::MPPICCloud::~MPPICCloud() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +void Foam::MPPICCloud::storeState() +{ + cloudCopyPtr_.reset + ( + static_cast*> + ( + clone(this->name() + "Copy").ptr() + ) + ); +} + + +template +void Foam::MPPICCloud::restoreState() +{ + this->cloudReset(cloudCopyPtr_()); + cloudCopyPtr_.clear(); +} + + +template +void Foam::MPPICCloud::evolve() +{ + if (this->solution().canEvolve()) + { + typename parcelType::template + TrackingData > td(*this); + + this->solve(td); + } +} + + +template +template +void Foam::MPPICCloud::motion(TrackData& td) +{ + // Kinematic + // ~~~~~~~~~ + + // force calculation and tracking + td.part() = TrackData::tpLinearTrack; + CloudType::move(td, this->db().time().deltaTValue()); + + + // Preliminary + // ~~~~~~~~~~~ + + // switch forces off so they are not applied in corrector steps + this->forces().setCalcNonCoupled(false); + this->forces().setCalcCoupled(false); + + + // Damping + // ~~~~~~~ + + if (dampingModel_->active()) + { + // update averages + td.updateAverages(*this); + + // memory allocation and eulerian calculations + dampingModel_->cacheFields(true); + + // calculate the damping velocity corrections without moving the parcels + td.part() = TrackData::tpDampingNoTrack; + CloudType::move(td, this->db().time().deltaTValue()); + + // correct the parcel positions and velocities + td.part() = TrackData::tpCorrectTrack; + CloudType::move(td, this->db().time().deltaTValue()); + + // finalise and free memory + dampingModel_->cacheFields(false); + } + + + // Packing + // ~~~~~~~ + + if (packingModel_->active()) + { + // same procedure as for damping + td.updateAverages(*this); + packingModel_->cacheFields(true); + td.part() = TrackData::tpPackingNoTrack; + CloudType::move(td, this->db().time().deltaTValue()); + td.part() = TrackData::tpCorrectTrack; + CloudType::move(td, this->db().time().deltaTValue()); + packingModel_->cacheFields(false); + } + + + // Isotropy + // ~~~~~~~~ + + if (isotropyModel_->active()) + { + // update averages + td.updateAverages(*this); + + // apply isotropy model + isotropyModel_->calculate(); + } + + + // Final + // ~~~~~ + + // update cell occupancy + this->updateCellOccupancy(); + + // switch forces back on + this->forces().setCalcNonCoupled(true); + this->forces().setCalcCoupled(this->solution().coupled()); +} + + +template +void Foam::MPPICCloud::info() +{ + CloudType::info(); + + tmp alpha = this->theta(); + + Info<< " Min cell volume fraction = " + << gMin(alpha().internalField()) << endl; + Info<< " Max cell volume fraction = " + << gMax(alpha().internalField()) << endl; + + label nOutside = 0; + + forAllIter(typename CloudType, *this, iter) + { + typename CloudType::parcelType& p = iter(); + const tetIndices tetIs(p.cell(), p.tetFace(), p.tetPt(), this->mesh()); + nOutside += !tetIs.tet(this->mesh()).inside(p.position()); + } + + reduce(nOutside, plusOp