Turbulent dispersion model.

src/lagrangian/cfdemParticle/Make/files

@@ -91,6 +91,7 @@ $(forceModels)/pdCorrelation/pdCorrelation.C
 $(forceModels)/surfaceTensionForce/surfaceTensionForce.C
 $(forceModels)/gradPForceSmooth/gradPForceSmooth.C
 $(forceModels)/particleDeformation/particleDeformation.C
+$(forceModels)/turbulentDispersion/turbulentDispersion.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C

src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C

@@ -0,0 +1,207 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "turbulentDispersion.H"
+#include "addToRunTimeSelectionTable.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(turbulentDispersion, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    turbulentDispersion,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+turbulentDispersion::turbulentDispersion
+(
+    const dictionary& dict,
+    cfdemCloud& sm,
+    word type
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(type + "Props")),
+    interpolate_(propsDict_.lookupOrDefault<bool>("interpolation", false)),
+    ignoreCellsName_(propsDict_.lookupOrDefault<word>("ignoreCellsName","none")),
+    ignoreCells_(),
+    existIgnoreCells_(true),
+    wallIndicatorField_
+    (   IOobject
+        (
+            "wallIndicator",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+    ),
+    turbKinetcEnergyFieldName_(propsDict_.lookupOrDefault<word>("turbKinetcEnergyFieldName","k")),
+    turbKinetcEnergy_(sm.mesh().lookupObject<volScalarField> (turbKinetcEnergyFieldName_)),
+    minTurbKinetcEnergy_(propsDict_.lookupOrDefault<scalar>("minTurbKinetcEnergy", 0.0)),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    critVoidfraction_(propsDict_.lookupOrDefault<scalar>("critVoidfraction", 0.9)),
+    ranGen_(clock::getTime()+pid())
+{
+    if(ignoreCellsName_ != "none")
+    {
+       ignoreCells_.set(new cellSet(particleCloud_.mesh(),ignoreCellsName_));
+       Info << type << ": ignoring fluctuations in cellSet " << ignoreCells_().name() <<
+        " with " << ignoreCells_().size() << " cells." << endl;
+    }
+    else existIgnoreCells_ = false;
+
+    // define a field to indicate if a cell is next to boundary
+    // label cellI = -1;
+    // // set wall indicator field
+    // const fvMesh& mesh = sm.mesh();  
+    // forAll(mesh.boundary(),patchI)
+    // {
+    //     forAll(mesh.boundary()[patchI], faceI)
+    //     {
+    //         cellI = patch().faceCells()[faceI];
+    //         wallIndicatorField_[cellI] = 1.0;
+    //     }
+    // }
+
+    // make sure this is the last force model in list so that fluid velocity does not get overwritten
+    label numLastForceModel = sm.nrForceModels();
+    word lastForceModel = sm.forceModels()[numLastForceModel-1];
+    if (lastForceModel != "turbulentDispersion")
+    {
+        FatalError <<"Force model 'turbulentDispersion' needs to be last in list!\n" << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+turbulentDispersion::~turbulentDispersion()
+{
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+bool turbulentDispersion::ignoreCell(label cell) const
+{
+    if (!existIgnoreCells_) return false;
+    else return ignoreCells_()[cell];
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void turbulentDispersion::setForce() const
+{
+    vector position(0,0,0);
+    scalar k = 0.0;
+
+    vector flucU(0,0,0);
+    label cellI = 0;
+
+    interpolationCellPoint<scalar> turbKinetcEnergyInterpolator_(turbKinetcEnergy_);
+
+    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    {
+        cellI = particleCloud_.cellIDs()[index][0];
+        flucU = vector(0,0,0);
+        k=0.0;
+
+        if (cellI > -1 && !ignoreCell(cellI))
+        {
+            // particles in dilute regions follow fluid without fluctuations
+            if (voidfraction_[cellI] < critVoidfraction_)
+            {
+                if (interpolate_)
+                {
+                    position = particleCloud_.position(index);
+                    k = turbKinetcEnergyInterpolator_.interpolate(position,cellI);
+                }
+                else
+                {
+                    k = turbKinetcEnergy_[cellI];
+                }
+
+                if (k < minTurbKinetcEnergy_) k = minTurbKinetcEnergy_;
+
+                flucU=unitFlucDir()*Foam::sqrt(2.0*k);
+
+                // if particles are pushed through walls, the velocity fluctuations may be regulated here
+                // check if cell is adjacent to wall and invert corresponding component
+                // if (wallIndicatorField_[cellI] > 0.5)
+                // {
+                // }
+
+                for(int j=0;j<3;j++)
+                {
+                    particleCloud_.fluidVels()[index][j] += flucU[j];
+                }
+            }
+        }
+    }   
+}
+
+vector turbulentDispersion::unitFlucDir() const
+{
+    // unit random vector
+    // algorithm according to:
+    // Marsaglia. "Choosing a point from the surface of a sphere." The Annals of Mathematical Statistics 43.2 (1972): 645-646.
+    scalar v1(0.0);
+    scalar v2(0.0);
+    scalar s(10.0);
+    scalar s2(0.0);
+    vector rvec(0,0,0);
+    while(s>1.0)
+    {
+        v1=2*(ranGen_.scalar01()-0.5);
+        v2=2*(ranGen_.scalar01()-0.5);
+        s=v1*v1+v2*v2;
+    }
+    s2=Foam::sqrt(1-s);
+    rvec[0]=2*v1*s2;
+    rvec[1]=2*v2*s2;
+    rvec[2]=1-2*s;
+    return rvec;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.H

@@ -0,0 +1,113 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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.
+
+    CFDEMcoupling academic 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 CFDEMcoupling academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#ifndef turbulentDispersion_H
+#define turbulentDispersion_H
+
+#include "forceModel.H"
+#include "interpolationCellPoint.H"
+#include "Random.H"
+#include "autoPtr.H"
+#include "cellSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class turbulentDispersion Declaration
+\*---------------------------------------------------------------------------*/
+
+class turbulentDispersion
+:
+    public forceModel
+{
+protected:
+    dictionary propsDict_;
+
+    bool interpolate_;
+
+    // ignore fluctuations in region
+    word ignoreCellsName_;
+
+    autoPtr<cellSet> ignoreCells_;
+
+    bool existIgnoreCells_;
+
+    mutable volScalarField wallIndicatorField_;
+
+    word turbKinetcEnergyFieldName_;
+
+    const volScalarField& turbKinetcEnergy_;
+
+    scalar minTurbKinetcEnergy_;
+
+    word voidfractionFieldName_;
+
+    const volScalarField& voidfraction_;
+
+    scalar critVoidfraction_;
+
+    virtual vector unitFlucDir() const;
+
+    mutable Random ranGen_;
+
+    bool ignoreCell(label) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("turbulentDispersion");
+
+
+    // Constructors
+
+        //- Construct from components
+        turbulentDispersion
+        (
+            const dictionary& dict,
+            cfdemCloud& sm,
+            word type = "turbulentDispersion"
+        );
+
+    // Destructor
+
+        ~turbulentDispersion();
+
+
+    // Member Functions
+        void setForce() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticleComp/Make/files

@@ -82,6 +82,7 @@ $(forceModels)/Fines/FanningDynFines.C
 $(forceModels)/Fines/ErgunStatFines.C
 $(forceModels)/granKineticEnergy/granKineticEnergy.C
 $(forceModels)/particleDeformation/particleDeformation.C
+$(forceModels)/turbulentDispersion/turbulentDispersion.C
 
 $(forceSubModels)/forceSubModel/newForceSubModel.C
 $(forceSubModels)/forceSubModel/forceSubModel.C