TUT: Adding new MPPICDyMFoam and uncoupledKinematicParcelDyMFoam tutorials

STY: Style and header-content changes

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

@@ -111,7 +111,6 @@ int main(int argc, char *argv[])
         continuousPhaseTransport.correct();
         muc = rhoc*continuousPhaseTransport.nu();
 
-        Info<< "Evolving " << kinematicCloud.name() << endl;
         kinematicCloud.evolve();
 
         // Update continuous phase volume fraction field

src/lagrangian/basic/particle/particleTemplates.C

@@ -470,13 +470,6 @@ void Foam::particle::hitCyclicAMIPatch
     const label sendFacei = cpp.whichFace(facei_);
     const label receiveFacei = cpp.pointFace(sendFacei, displacement, pos);
 
-    if (false)
-    {
-        Info<< "My new pos : " <<  pos << endl;
-        Info<< "Particle " << origId() << " crossing AMI from " << cpp.name()
-            << " to " << receiveCpp.name() << endl << endl;
-    }
-
     if (receiveFacei < 0)
     {
         // If the patch face of the particle is not known assume that the
@@ -495,9 +488,15 @@ void Foam::particle::hitCyclicAMIPatch
     vector displacementT = displacement;
     cpp.reverseTransformDirection(displacementT, sendFacei);
 
+    // NOTE: The ray used to find the hit location accross the AMI might not
+    // be consistent in the displacement direction. Therefore a particle can
+    // be looping accross AMI patches indefinitely. Advancing the particle
+    // trajectory inside the cell is a possible solution.
+    const vector dispDir = cpp.fraction()*displacementT;
+    stepFraction_ += cpp.fraction();
     locate
     (
-        pos,
+        pos + dispDir,
         &displacementT,
         mesh_.faceOwner()[facei_],
         false,
@@ -534,7 +533,6 @@ void Foam::particle::hitCyclicAMIPatch
 
     //if (onBoundaryFace())
     {
-        //DebugVar("On boudanry")
 //         vector receiveNormal, receiveDisplacement;
 //         patchData(receiveNormal, receiveDisplacement);
 //

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

@@ -758,12 +758,10 @@ void Foam::KinematicCloud<CloudType>::preEvolve
 
     if (this->dampingModel().active())
     {
-        DebugVar("dampingModel")
         this->dampingModel().cacheFields(true);
     }
     if (this->packingModel().active())
     {
-        DebugVar("packingModel")
         this->packingModel().cacheFields(true);
     }
 
@@ -895,6 +893,48 @@ void Foam::KinematicCloud<CloudType>::info()
     injectors_.info(Info);
     this->surfaceFilm().info(Info);
     this->patchInteraction().info(Info);
+
+    if (this->packingModel().active())
+    {
+        tmp<volScalarField> alpha = this->theta();
+
+        if (this->db().time().writeTime())
+        {
+            alpha().write();
+        }
+
+        const scalar alphaMin = gMin(alpha().primitiveField());
+        const scalar alphaMax = gMax(alpha().primitiveField());
+
+        Info<< "    Min cell volume fraction        = " << alphaMin << endl;
+        Info<< "    Max cell volume fraction        = " << alphaMax << endl;
+
+        if (alphaMax < SMALL)
+        {
+            return;
+        }
+
+        scalar nMin = GREAT;
+
+        forAll(this->mesh().cells(), celli)
+        {
+            const label n = this->cellOccupancy()[celli].size();
+
+            if (n > 0)
+            {
+                const scalar nPack = n*alphaMax/alpha()[celli];
+
+                if (nPack < nMin)
+                {
+                    nMin = nPack;
+                }
+            }
+        }
+
+        reduce(nMin, minOp<scalar>());
+
+        Info<< "    Min dense number of parcels     = " << nMin << endl;
+    }
 }
 
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -361,7 +361,6 @@ bool Foam::KinematicParcel<ParcelType>::move
 
         const scalar dt = (p.stepFraction() - sfrac)*trackTime;
 
-
         // Avoid problems with extremely small timesteps
         if (dt > ROOTVSMALL)
         {
@@ -386,12 +385,11 @@ bool Foam::KinematicParcel<ParcelType>::move
         {
             cloud.functions().postFace(p, ttd.keepParticle);
         }
-
         cloud.functions().postMove(p, dt, start, ttd.keepParticle);
 
         if (p.active() && p.onFace() && ttd.keepParticle)
         {
-            p.hitFace(f*s - d, cloud, ttd);
+            p.hitFace(s, cloud, ttd);
         }
     }
 

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/LocalInteraction/LocalInteraction.C

@@ -294,11 +294,11 @@ bool Foam::LocalInteraction<CloudType>::correct
                 // Calculate motion relative to patch velocity
                 U -= Up;
 
-                if (mag(U) < this->Urmax())
+                if (mag(Up) > 0 && mag(U) < this->Urmax())
                 {
                     WarningInFunction
-                    << "Particle U the same as patch "
-                    << "    The particle has been removed" << nl << endl;
+                        << "Particle U the same as patch "
+                        << "    The particle has been removed" << nl << endl;
 
                     keepParticle = false;
                     p.active(false);

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/PatchInteractionModel/PatchInteractionModel.C

@@ -162,7 +162,7 @@ Foam::PatchInteractionModel<CloudType>::PatchInteractionModel
     UName_(this->coeffDict().template getOrDefault<word>("U", "U")),
     escapedParcels_(0),
     escapedMass_(0.0),
-    Urmax_(this->coeffDict().template getOrDefault<scalar>("UrMax", 1e-4))
+    Urmax_(this->coeffDict().template getOrDefault<scalar>("UrMax", 0))
 {}
 
 

src/lagrangian/intermediate/submodels/Kinematic/PatchInteractionModel/StandardWallInteraction/StandardWallInteraction.C

@@ -211,11 +211,11 @@ bool Foam::StandardWallInteraction<CloudType>::correct
                 // Calculate motion relative to patch velocity
                 U -= Up;
 
-                if (mag(U) < this->Urmax())
+                if (mag(Up) > 0 && mag(U) < this->Urmax())
                 {
                     WarningInFunction
-                    << "Particle U the same as patch "
-                    << "    The particle has been removed" << nl << endl;
+                        << "Particle U the same as patch "
+                        << "    The particle has been removed" << nl << endl;
 
                     keepParticle = false;
                     p.active(false);

src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C

@@ -92,8 +92,6 @@ Foam::AveragingMethod<Type>::New
         dict.template getOrDefault<word>(typeName, "basic")
     );
 
-    //Info<< "Selecting averaging method " << modelType << endl;
-
     auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
     if (!cstrIter.found())

src/lagrangian/spray/parcels/derived/basicSprayParcel/makeBasicSprayParcelSubmodels.C

@@ -52,6 +52,11 @@ License
 #include "makeSprayParcelAtomizationModels.H"
 #include "makeSprayParcelBreakupModels.H"
 
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 makeParcelCloudFunctionObjects(basicSprayCloud);
@@ -78,5 +83,9 @@ makeParticleForceModelType(DistortedSphereDragForce, basicSprayCloud);
 makeSprayParcelAtomizationModels(basicSprayCloud);
 makeSprayParcelBreakupModels(basicSprayCloud);
 
+// MPPIC sub-models
+makeMPPICParcelDampingModels(basicSprayCloud);
+makeMPPICParcelIsotropyModels(basicSprayCloud);
+makeMPPICParcelPackingModels(basicSprayCloud);
 
 // ************************************************************************* //

src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C

@@ -1069,7 +1069,8 @@ const
     {
         const face& f = srcPatch[srcFacei];
 
-        pointHit ray = f.ray(tgtPoint, n, srcPoints);
+        pointHit ray =
+            f.ray(tgtPoint, n, srcPoints, intersection::algorithm::VISIBLE);
 
         if (ray.hit())
         {
@@ -1078,6 +1079,7 @@ const
         }
         else if (ray.distance() < nearest.distance())
         {
+
             nearest = ray;
             nearestFacei = srcFacei;
         }
@@ -1116,7 +1118,8 @@ const
     {
         const face& f = tgtPatch[tgtFacei];
 
-        pointHit ray = f.ray(srcPoint, n, tgtPoints);
+        pointHit ray =
+            f.ray(srcPoint, n, tgtPoints, intersection::algorithm::VISIBLE);
 
         if (ray.hit())
         {
@@ -1131,9 +1134,9 @@ const
             nearestFacei = tgtFacei;
         }
     }
-
     if (nearest.hit() || nearest.eligibleMiss())
     {
+
         srcPoint = nearest.rawPoint();
         return nearestFacei;
     }

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C

@@ -568,6 +568,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     coupledPolyPatch(name, size, start, index, bm, patchType, transform),
     nbrPatchName_(word::null),
     nbrPatchID_(-1),
+    fraction_(Zero),
     rotationAxis_(Zero),
     rotationCentre_(Zero),
     rotationAngleDefined_(false),
@@ -603,6 +604,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     nbrPatchName_(dict.getOrDefault<word>("neighbourPatch", word::null)),
     coupleGroup_(dict),
     nbrPatchID_(-1),
+    fraction_(dict.getOrDefault<scalar>("fraction", Zero)),
     rotationAxis_(Zero),
     rotationCentre_(Zero),
     rotationAngleDefined_(false),
@@ -707,6 +709,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     nbrPatchName_(pp.nbrPatchName_),
     coupleGroup_(pp.coupleGroup_),
     nbrPatchID_(-1),
+    fraction_(pp.fraction_),
     rotationAxis_(pp.rotationAxis_),
     rotationCentre_(pp.rotationCentre_),
     rotationAngleDefined_(pp.rotationAngleDefined_),
@@ -742,6 +745,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     nbrPatchName_(nbrPatchName),
     coupleGroup_(pp.coupleGroup_),
     nbrPatchID_(-1),
+    fraction_(pp.fraction_),
     rotationAxis_(pp.rotationAxis_),
     rotationCentre_(pp.rotationCentre_),
     rotationAngleDefined_(pp.rotationAngleDefined_),
@@ -784,6 +788,7 @@ Foam::cyclicAMIPolyPatch::cyclicAMIPolyPatch
     nbrPatchName_(pp.nbrPatchName_),
     coupleGroup_(pp.coupleGroup_),
     nbrPatchID_(-1),
+    fraction_(pp.fraction_),
     rotationAxis_(pp.rotationAxis_),
     rotationCentre_(pp.rotationCentre_),
     rotationAngleDefined_(pp.rotationAngleDefined_),
@@ -1153,6 +1158,8 @@ void Foam::cyclicAMIPolyPatch::write(Ostream& os) const
         os.writeEntry("tgtSize", tgtFaceIDs_.size());
         os.writeEntry("moveFaceCentres", moveFaceCentres_);
     }
+
+    os.writeEntryIfDifferent<scalar>("fraction", Zero, fraction_);
 }
 
 

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.H

@@ -97,6 +97,9 @@ protected:
         //- Index of other half
         mutable label nbrPatchID_;
 
+        //- Particle displacement fraction accross AMI
+        const scalar fraction_;
+
 
         // Transformations
 
@@ -357,6 +360,9 @@ public:
             //- Neighbour patch ID
             virtual label neighbPatchID() const;
 
+            //- Particle fraction increase between AMI pathces
+            inline scalar fraction() const;
+
             //- Does this side own the patch?
             virtual bool owner() const;
 

src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatchI.H

@@ -58,6 +58,13 @@ inline const Foam::word& Foam::cyclicAMIPolyPatch::neighbPatchName() const
     return nbrPatchName_;
 }
 
+
+inline Foam::scalar Foam::cyclicAMIPolyPatch::fraction() const
+{
+    return fraction_;
+}
+
+
 inline const Foam::scalarListList& Foam::cyclicAMIPolyPatch::weights() const
 {
     if (owner())

tutorials/lagrangian/MPPICDyMFoam/column/0/U.air

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      binary;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    sides
+    {
+        type            symmetry;
+    }
+
+    top
+    {
+        type            noSlip;
+    }
+
+    bottom
+    {
+        type            noSlip;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/0/p

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    sides
+    {
+        type            symmetry;
+    }
+
+    top
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+
+    bottom
+    {
+        type            fixedFluxPressure;
+        value           $internalField;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/constant/dynamicMeshDict

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicMotionSolverFvMesh;
+
+motionSolver    solidBody;
+
+solidBodyMotionFunction  oscillatingLinearMotion;
+
+amplitude       (3 0 0);
+omega           0.1; // rad/s  (.5 rps)
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/constant/g

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       uniformDimensionedVectorField;
+    object      g;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -2 0 0 0 0];
+value           (0 0 -9.81);
+
+
+// ************************************************************************* //

tutorials/lagrangian/MPPICDyMFoam/column/constant/kinematicCloudProperties

+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2006                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      kinematicCloudProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solution
+{
+    active          true;
+    coupled         true;
+    transient       yes;
+    cellValueSourceCorrection off;
+
+    interpolationSchemes
+    {
+        rho.air         cell;
+        U.air           cellPoint;
+        mu.air          cell;
+    }
+
+    averagingMethod dual;
+
+    integrationSchemes
+    {
+        U               Euler;
+    }
+
+    sourceTerms
+    {
+        schemes
+        {
+            U     semiImplicit 1;
+        }
+    }
+}
+
+constantProperties
+{
+    rho0            2526;
+    alphaMax        0.9;
+}
+
+subModels
+{
+    particleForces
+    {
+        ErgunWenYuDrag
+        {
+            alphac alpha.air;
+        }
+    }
+
+    injectionModels
+    {
+        model1
+        {
+            type            manualInjection;
+            massTotal       0;
+            parcelBasisType fixed;
+            nParticle       750;
+            SOI             0;
+            positionsFile   "kinematicCloudPositions";
+            U0              (0 0 0);
+            sizeDistribution
+            {
+                type        fixedValue;
+                fixedValueDistribution
+                {
+                    value   0.0003;
+                }
+            }
+        }
+    }
+
+    dispersionModel none;
+
+    patchInteractionModel localInteraction;
+
+    localInteractionCoeffs
+    {
+        patches
+        (
+            sides
+            {
+                type rebound;
+                e    1;
+                mu   0.09;
+            }
+            top
+            {
+                type rebound;
+                e    1;
+                mu   0.09;
+            }
+            bottom
+            {
+                type rebound;
+                e    1;
+                mu   0.09;
+            }
+        );
+    }
+
+    heatTransferModel none;
+
+    surfaceFilmModel none;
+
+    packingModel implicit;
+
+    explicitCoeffs
+    {
+        particleStressModel
+        {
+            type HarrisCrighton;
+            alphaPacked 0.6;
+            pSolid 10.0;
+            beta 2.0;
+            eps 1.0e-7;
+        }
+        correctionLimitingMethod
+        {
+            type absolute;
+            e 0.9;