Merge branch 'MPPIC'

applications/solvers/lagrangian/DPMFoam/Allwclean

+#!/bin/sh
+
+cd ${0%/*} || exit 1
+set -x
+
+wclean DPMTurbulenceModels
+wclean
+wclean MPPICFoam

applications/solvers/lagrangian/DPMFoam/Allwmake

+#!/bin/sh
+
+cd ${0%/*} || exit 1
+set -x
+
+wmake DPMTurbulenceModels
+wmake
+wmake MPPICFoam

applications/solvers/lagrangian/DPMFoam/CourantNo.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+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;
+
+// ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/DPMFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+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;
+}
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/DPMTurbulenceModels.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#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);
+
+// ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/files

+DPMTurbulenceModels.C
+
+LIB = $(FOAM_LIBBIN)/libDPMTurbulenceModels

applications/solvers/lagrangian/DPMFoam/DPMTurbulenceModels/Make/options

+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

applications/solvers/lagrangian/DPMFoam/MPPICFoam/MPPICFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+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"
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/files

+MPPICFoam.C
+
+EXE = $(FOAM_APPBIN)/MPPICFoam

applications/solvers/lagrangian/DPMFoam/MPPICFoam/Make/options

+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

applications/solvers/lagrangian/DPMFoam/Make/files

+DPMFoam.C
+
+EXE = $(FOAM_APPBIN)/DPMFoam

applications/solvers/lagrangian/DPMFoam/Make/options

+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

applications/solvers/lagrangian/DPMFoam/UcEqn.H

+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()
+        )
+    );
+}

applications/solvers/lagrangian/DPMFoam/continuityErrs.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <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;
+}
+
+// ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/createFields.H

+    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<PhaseIncompressibleTurbulenceModel<singlePhaseTransportModel> >
+    continuousPhaseTurbulence
+    (
+        PhaseIncompressibleTurbulenceModel<singlePhaseTransportModel>::New
+        (
+            alphac,
+            Uc,
+            alphaPhic,
+            phic,
+            continuousPhaseTransport
+        )
+    );

applications/solvers/lagrangian/DPMFoam/pEqn.H

+{
+    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<fixedFluxPressureFvPatchScalarField>
+    (
+        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"

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

src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "MPPICCloud.H"
+#include "PackingModel.H"
+#include "ParticleStressModel.H"
+#include "DampingModel.H"
+#include "IsotropyModel.H"
+#include "TimeScaleModel.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::MPPICCloud<CloudType>::setModels()
+{
+    packingModel_.reset
+    (
+        PackingModel<MPPICCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+    dampingModel_.reset
+    (
+        DampingModel<MPPICCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+    isotropyModel_.reset
+    (
+        IsotropyModel<MPPICCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::MPPICCloud<CloudType>::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<CloudType>::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<class CloudType>
+Foam::MPPICCloud<CloudType>::MPPICCloud
+(
+    MPPICCloud<CloudType>& c,
+    const word& name
+)
+:
+    CloudType(c, name),
+    packingModel_(c.packingModel_->clone()),
+    dampingModel_(c.dampingModel_->clone()),
+    isotropyModel_(c.isotropyModel_->clone())
+{}
+
+
+template<class CloudType>
+Foam::MPPICCloud<CloudType>::MPPICCloud
+(
+    const fvMesh& mesh,
+    const word& name,
+    const MPPICCloud<CloudType>& c
+)
+:
+    CloudType(mesh, name, c),
+    packingModel_(NULL),
+    dampingModel_(NULL),
+    isotropyModel_(NULL)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::MPPICCloud<CloudType>::~MPPICCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::MPPICCloud<CloudType>::storeState()
+{
+    cloudCopyPtr_.reset
+    (
+        static_cast<MPPICCloud<CloudType>*>
+        (
+            clone(this->name() + "Copy").ptr()
+        )
+    );
+}
+
+
+template<class CloudType>
+void Foam::MPPICCloud<CloudType>::restoreState()
+{
+    this->cloudReset(cloudCopyPtr_());
+    cloudCopyPtr_.clear();
+}
+
+
+template<class CloudType>
+void Foam::MPPICCloud<CloudType>::evolve()
+{
+    if (this->solution().canEvolve())
+    {
+        typename parcelType::template
+            TrackingData<MPPICCloud<CloudType> > td(*this);
+
+        this->solve(td);
+    }
+}
+
+
+template<class CloudType>
+template<class TrackData>
+void Foam::MPPICCloud<CloudType>::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<class CloudType>
+void Foam::MPPICCloud<CloudType>::info()
+{
+    CloudType::info();
+
+    tmp<volScalarField> 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<label>());
+
+    if (nOutside > 0)
+    {
+        Info<< "    Number of parcels outside tets  = " << nOutside << endl;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloud.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::MPPICCloud
+
+Description
+    Adds collisions to kinematic clouds
+
+SourceFiles
+    MPPICCloudI.H
+    MPPICCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MPPICCloud_H
+#define MPPICCloud_H
+
+#include "particle.H"
+#include "Cloud.H"
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "fvMesh.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+template<class CloudType>
+class PackingModel;
+
+template<class CloudType>
+class DampingModel;
+
+template<class CloudType>
+class IsotropyModel;
+
+/*---------------------------------------------------------------------------*\
+                       Class MPPICCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class MPPICCloud
+:
+    public CloudType
+{
+public:
+
+    // Public typedefs
+
+        //- Type of cloud this cloud was instantiated for
+        typedef CloudType cloudType;
+
+        //- Type of parcel the cloud was instantiated for
+        typedef typename CloudType::parcelType parcelType;
+
+        //- Convenience typedef for this cloud type
+        typedef MPPICCloud<CloudType> MPPICCloudType;
+
+
+private:
+
+    // Private data
+
+        //- Cloud copy pointer
+        autoPtr<MPPICCloud<CloudType> > cloudCopyPtr_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        MPPICCloud(const MPPICCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const MPPICCloud&);
+
+
+protected:
+
+    // Protected data
+
+        // References to the cloud sub-models
+
+            //- Packing model
+            autoPtr<PackingModel<MPPICCloud<CloudType> > > packingModel_;
+
+            //- Damping model
+            autoPtr<DampingModel<MPPICCloud<CloudType> > >
+                dampingModel_;
+
+            //- Exchange model
+            autoPtr<IsotropyModel<MPPICCloud<CloudType> > >
+                isotropyModel_;
+
+
+        // Initialisation
+
+            //- Set cloud sub-models
+            void setModels();
+
+
+public:
+
+    // Constructors
+
+        //- Construct given carrier gas fields
+        MPPICCloud
+        (
+            const word& cloudName,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const volScalarField& mu,
+            const dimensionedVector& g,
+            bool readFields = true
+        );
+
+        //- Copy constructor with new name
+        MPPICCloud
+        (
+            MPPICCloud<CloudType>& c,
+            const word& name
+        );
+
+        //- Copy constructor with new name - creates bare cloud
+        MPPICCloud
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const MPPICCloud<CloudType>& c
+        );
+
+        //- Construct and return clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType> > clone(const word& name)
+        {
+            return autoPtr<Cloud<parcelType> >
+            (
+                new MPPICCloud(*this, name)
+            );
+        }
+
+        //- Construct and return bare clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType> > cloneBare(const word& name) const
+        {
+            return autoPtr<Cloud<parcelType> >
+            (
+                new MPPICCloud(this->mesh(), name, *this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~MPPICCloud();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return a reference to the cloud copy
+            inline const MPPICCloud& cloudCopy() const;
+
+            //- Return const access to the packing model
+            inline const PackingModel<MPPICCloud<CloudType> >&
+                packingModel() const;
+
+            //- Return a reference to the packing model
+            inline PackingModel<MPPICCloud<CloudType> >& packingModel();
+
+            //- Return condt access to the damping model
+            inline const DampingModel<MPPICCloud<CloudType> >&
+                dampingModel() const;
+
+            //- Return a reference to the damping model
+            inline DampingModel<MPPICCloud<CloudType> >& dampingModel();
+
+            //- Return condt access to the isotropy model
+            inline const IsotropyModel<MPPICCloud<CloudType> >&
+                isotropyModel() const;
+
+            //- Return a reference to the isotropy model
+            inline IsotropyModel<MPPICCloud<CloudType> >& isotropyModel();
+
+
+        // Cloud evolution functions
+
+            //- Store the current cloud state
+            void storeState();
+
+            //- Reset the current cloud to the previously stored state
+            void restoreState();
+
+            //- Evolve the cloud
+            void evolve();
+
+            //- Particle motion
+            template<class TrackData>
+            void motion(TrackData& td);
+
+
+        //- I-O
+
+            //- Print cloud information
+            void info();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "MPPICCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "MPPICCloud.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/MPPICCloud/MPPICCloudI.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+inline const Foam::MPPICCloud<CloudType>&
+Foam::MPPICCloud<CloudType>::cloudCopy() const
+{
+    return cloudCopyPtr_();
+}
+
+
+template<class CloudType>
+inline const Foam::PackingModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::packingModel() const
+{
+    return packingModel_();
+}
+
+
+template<class CloudType>
+inline Foam::PackingModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::packingModel()
+{
+    return packingModel_();
+}
+
+
+template<class CloudType>
+inline const Foam::DampingModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::dampingModel() const
+{
+    return dampingModel_();
+}
+
+
+template<class CloudType>
+inline Foam::DampingModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::dampingModel()
+{
+    return dampingModel_();
+}
+
+
+template<class CloudType>
+inline const Foam::IsotropyModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::isotropyModel() const
+{
+    return isotropyModel_();
+}
+
+
+template<class CloudType>
+inline Foam::IsotropyModel<Foam::MPPICCloud<CloudType> >&
+Foam::MPPICCloud<CloudType>::isotropyModel()
+{
+    return isotropyModel_();
+}
+
+
+// ************************************************************************* //