DPMDyMFoam, MPPICDyMFoam: New dynamic mesh versions of DPMFoam and MPPICFoam

supporting both mesh morphing and topology change.

applications/solvers/compressible/rhoPimpleFoam/pEqn.H

@@ -10,10 +10,8 @@ if (pimple.nCorrPISO() <= 1)
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
-    (
-        fvc::flux(rho*HbyA)
-      + rhorAUf*fvc::ddtCorr(rho, U, phi)
-    )
+    fvc::flux(rho*HbyA)
+  + rhorAUf*fvc::ddtCorr(rho, U, phi)
 );
 
 MRF.makeRelative(fvc::interpolate(rho), phiHbyA);

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -67,7 +67,6 @@ int main(int argc, char *argv[])
     {
         #include "readControls.H"
         #include "CourantNo.H"
-
         #include "setDeltaT.H"
 
         runTime++;

applications/solvers/lagrangian/DPMFoam/Allwclean

 #!/bin/sh
-
-cd ${0%/*} || exit 1
+cd ${0%/*} || exit 1    # Run from this directory
 
 wclean libso DPMTurbulenceModels
+
 wclean
 wclean MPPICFoam
+wclean DPMDyMFoam
+wclean DPMDyMFoam/MPPICDyMFoam
+
+#------------------------------------------------------------------------------

applications/solvers/lagrangian/DPMFoam/Allwmake

 #!/bin/sh
-cd ${0%/*} || exit 1
+cd ${0%/*} || exit 1    # Run from this directory
 
 # Parse arguments for library compilation
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
@@ -8,3 +8,7 @@ wmake $targetType DPMTurbulenceModels
 
 wmake $targetType
 wmake $targetType MPPICFoam
+wmake $targetType DPMDyMFoam
+wmake $targetType DPMDyMFoam/MPPICDyMFoam
+
+#------------------------------------------------------------------------------

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/DPMDyMFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    DPMDyMFoam
+
+Description
+    Transient solver for the coupled transport of a single kinematic particle
+    cloud including the effect of the volume fraction of particles on the
+    continuous phase, with optional mesh motion and mesh topology changes.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "singlePhaseTransportModel.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
+#include "pimpleControl.H"
+#include "CorrectPhi.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 "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createControls.H"
+    #include "createFields.H"
+    #include "createUcf.H"
+    #include "initContinuityErrs.H"
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        mesh.update();
+
+        // Calculate absolute flux from the mapped surface velocity
+        phic = mesh.Sf() & Ucf;
+
+        if (mesh.changing() && correctPhi)
+        {
+            #include "correctPhic.H"
+        }
+
+        // Make the flux relative to the mesh motion
+        fvc::makeRelative(phic, Uc);
+
+        if (mesh.changing() && checkMeshCourantNo)
+        {
+            #include "meshCourantNo.H"
+        }
+
+        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,
+                Zero
+            ),
+            zeroGradientFvPatchVectorField::typeName
+        );
+
+        cloudVolSUSu.primitiveFieldRef() = -cloudSU.source()/mesh.V();
+        cloudVolSUSu.correctBoundaryConditions();
+        cloudSU.source() = 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/DPMDyMFoam/MPPICDyMFoam/MPPICDyMFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    MPPICDyMFoam
+
+Description
+    Transient solver for the coupled transport of a single kinematic particle
+    cloud 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,
+    with optional mesh motion and mesh topology changes.
+
+\*---------------------------------------------------------------------------*/
+
+#define MPPIC
+
+#include "DPMDyMFoam.C"
+
+
+// ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/files

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

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/MPPICDyMFoam/Make/options

+EXE_INC = \
+    -I.. \
+    -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)/transportModels/compressible/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 \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lspecie \
+    -lradiationModels \
+    -lincompressibleTransportModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lDPMTurbulenceModels \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lsampling \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/files

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

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/Make/options

+EXE_INC = \
+    -I.. \
+    -I../DPMTurbulenceModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/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 \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude
+
+EXE_LIBS = \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -llagrangianTurbulence \
+    -lspecie \
+    -lradiationModels \
+    -lincompressibleTransportModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lDPMTurbulenceModels \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lsampling \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/correctPhic.H

+CorrectPhi
+(
+    Uc,
+    phic,
+    p,
+    dimensionedScalar("rAUf", dimTime, 1),
+    geometricZeroField(),
+    pimple
+);
+
+#include "continuityErrs.H"

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/createControls.H

+#include "createControl.H"
+#include "createTimeControls.H"
+
+bool correctPhi
+(
+    pimple.dict().lookupOrDefault("correctPhi", false)
+);
+
+bool checkMeshCourantNo
+(
+    pimple.dict().lookupOrDefault("checkMeshCourantNo", false)
+);

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/createUcf.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    createUcf
+
+Description
+    Creates and initialises the velocity velocity field Ucf.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Info<< "Reading/calculating continuous phase face velocity Ucf\n" << endl;
+
+surfaceVectorField Ucf
+(
+    IOobject
+    (
+        "Ucf",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    fvc::interpolate(Uc)
+);
+
+// ************************************************************************* //

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/pEqn.H

+{
+    volVectorField HbyA("HbyA", Uc);
+    HbyA = rAUc*UcEqn.H();
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+           fvc::flux(HbyA)
+         + alphacf*rAUcf*fvc::ddtCorr(Uc, Ucf)
+         + phicForces
+        )
+    );
+
+    if (p.needReference())
+    {
+        fvc::makeRelative(phiHbyA, Uc);
+        adjustPhi(phiHbyA, Uc, p);
+        fvc::makeAbsolute(phiHbyA, Uc);
+    }
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, Uc, phiHbyA, rAUcf);
+
+    // 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();
+
+            {
+                Ucf = fvc::interpolate(Uc);
+                surfaceVectorField n(mesh.Sf()/mesh.magSf());
+                Ucf += n*(phic/mesh.magSf() - (n & Ucf));
+            }
+
+            // Make the fluxes relative to the mesh motion
+            fvc::makeRelative(phic, Uc);
+        }
+    }
+}
+
+#include "continuityErrs.H"

applications/solvers/lagrangian/DPMFoam/DPMDyMFoam/readControls.H

+#include "readTimeControls.H"
+
+correctPhi = pimple.dict().lookupOrDefault("correctPhi", false);
+
+checkMeshCourantNo = pimple.dict().lookupOrDefault("checkMeshCourantNo", false);