zhyang-dev/CFDEMcoupling-PFM
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Some cleaning up of fines models.

Browse Source
This commit is contained in:
Thomas Lichtenegger
2016-07-05 10:27:41 +02:00
parent 1be85551c1
commit 2da9631f20
7 changed files with 569 additions and 142 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C
View File

@ -188,7 +188,7 @@ void ErgunStatFines::setForce() const
if(voidfraction > switchingVoidfraction_) //dilute
{
Rep=dSauterMix*voidfraction*magUr/nuf;
CdMagUrLag = (24.0*nuf/(ddSauterMix*voidfraction)) //1/magUr missing here, but compensated in expression for betaP!
CdMagUrLag = (24.0*nuf/(dSauterMix*voidfraction)) //1/magUr missing here, but compensated in expression for betaP!
*(scalar(1.0)+0.15*Foam::pow(Rep, 0.687));
betaP = 0.75* localPhiP * (

33
src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesDynFanning.C → src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C
View File

@ -17,7 +17,7 @@ License
#include "error.H"
#include "FinesDynFanning.H"
#include "FanningDynFines.H"
#include "addToRunTimeSelectionTable.H"
#include "averagingModel.H"
@ -28,12 +28,12 @@ namespace Foam
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(FinesDynFanning, 0);
defineTypeNameAndDebug(FanningDynFines, 0);
addToRunTimeSelectionTable
(
forceModel,
FinesDynFanning,
FanningDynFines,
dictionary
);
@ -41,7 +41,7 @@ addToRunTimeSelectionTable
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
FinesDynFanning::FinesDynFanning
FanningDynFines::FanningDynFines
(
const dictionary& dict,
cfdemCloud& sm
@ -56,14 +56,9 @@ FinesDynFanning::FinesDynFanning
UsFieldName_(propsDict_.lookup("granVelFieldName")),
UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
UDyn_(sm.mesh().lookupObject<volVectorField> ("uDyn")),
alphaDyn_(sm.mesh().lookupObject<volScalarField> ("alphaDyn")),
FanningCoeff_(sm.mesh().lookupObject<volScalarField> ("FanningCoeff")),
alphaP_(sm.mesh().lookupObject<volScalarField> ("alphaP")),
dHydMix_(sm.mesh().lookupObject<volScalarField> ("dHydMix")),
dSauter_(sm.mesh().lookupObject<volScalarField> ("dSauter")),
Froude_(sm.mesh().lookupObject<volScalarField> ("Froude")),
rhoDyn_(readScalar(propsDict_.lookup ("rhoDyn"))),
prefactor_(readScalar(propsDict_.lookup ("prefactor"))),
exponent_(readScalar(propsDict_.lookup ("exponent"))),
scaleDia_(1.),
scaleDrag_(1.)
{
@ -87,13 +82,13 @@ FinesDynFanning::FinesDynFanning
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
FinesDynFanning::~FinesDynFanning()
FanningDynFines::~FanningDynFines()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void FinesDynFanning::setForce() const
void FanningDynFines::setForce() const
{
vector UDyn(0,0,0);
vector drag(0,0,0);
@ -104,11 +99,7 @@ void FinesDynFanning::setForce() const
scalar ds(0);
scalar ds_scaled(0);
scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_;
scalar magUr(0);
scalar Fk(0);
vector dragExplicit(0,0,0);
scalar dragCoefficient(0);
@ -118,7 +109,6 @@ void FinesDynFanning::setForce() const
{
cellI = particleCloud_.cellIDs()[index][0];
drag = vector(0,0,0);
dragExplicit = vector(0,0,0);
UDyn = vector(0,0,0);
dragCoefficient = 0;
@ -127,16 +117,13 @@ void FinesDynFanning::setForce() const
UDyn = UDyn_[cellI];
Us = UsField_[cellI];
Ur = UDyn-Us;
magUr = mag(Ur);
ds = 2*particleCloud_.radius(index);
ds_scaled = ds/scaleDia_;
Fk = prefactor_*Foam::pow(Froude_[cellI], exponent_);
dragCoefficient = alphaDyn_[cellI] * rhoDyn_ * magUr * Fk / (2 * dHydMix_[cellI] )
dragCoefficient = FanningCoeff_[cellI];
// calc particle's drag
dragCoefficient *= M_PI/6 * ds_scaled * ds_scaled / alphaP_[cellI] * dSauter_[cellI] *scaleDia3*scaleDrag_;
dragCoefficient *= M_PI/6 * ds_scaled * ds_scaled / alphaP_[cellI] * dSauter_[cellI] * scaleDia3 * scaleDrag_;
if (modelType_=="B")
dragCoefficient /= voidfraction_[cellI];
@ -144,7 +131,7 @@ void FinesDynFanning::setForce() const
}
// write particle based data to global array
forceSubM(0).partToArray(index,drag,dragExplicit);
forceSubM(0).partToArray(index,drag,vector::zero);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

30
src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesDynFanning.H → src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.H
View File

@ -17,15 +17,15 @@ Description
fines influence on particles
SourceFiles
FinesDynFanning.C
FanningDynFines.C
\*---------------------------------------------------------------------------*/
#ifndef FinesDynFanning_H
#define FinesDynFanning_H
#ifndef FanningDynFines_H
#define FanningDynFines_H
#include "forceModel.H"
#include "interpolationCellPoint.H"
#include "FinesDynFanningFields.H"
#include "FanningDynFines.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -33,10 +33,10 @@ namespace Foam
{
/*---------------------------------------------------------------------------*\
Class FinesDynFanning Declaration
Class FanningDynFines Declaration
\*---------------------------------------------------------------------------*/
class FinesDynFanning
class FanningDynFines
:
public forceModel
{
@ -58,21 +58,11 @@ private:
const volVectorField& UDyn_;
const volScalarField& alphaDyn_;
const volScalarField& FanningCoeff_;
const volScalarField& alphaP_;
const volScalarField& dHydMix_;
const volScalarField& dSauter_;
const volScalarField& Froude_;
scalar rhoDyn_;
scalar prefactor_;
scalar exponent_;
mutable scalar scaleDia_;
@ -84,13 +74,13 @@ private:
public:
//- Runtime type information
TypeName("FinesDynFanning");
TypeName("FanningDynFines");
// Constructors
//- Construct from components
FinesDynFanning
FanningDynFines
(
const dictionary& dict,
cfdemCloud& sm
@ -98,7 +88,7 @@ public:
// Destructor
~FinesDynFanning();
~FanningDynFines();
// Member Functions

272
src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C
View File

@ -40,60 +40,29 @@ FinesFields::FinesFields
:
particleCloud_(sm),
propsDict_(dict.subDict(typeName + "Props")),
velFieldName_(propsDict_.lookup("velFieldName")),
velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
// this is probably really bad
voidfraction_(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_))),
UsFieldName_(propsDict_.lookup("granVelFieldName")),
UsFieldName_(propsDict_.lookupOrDefault<word>("granVelFieldName","Us")),
UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
dSauter_(sm.mesh().lookupObject<volScalarField> ("dSauter")),
dSauterMix_
pFieldName_(propsDict_.lookupOrDefault<word>("pFieldName","p")),
p_(sm.mesh().lookupObject<volScalarField> (pFieldName_)),
rhoGFieldName_(propsDict_.lookupOrDefault<word>("rhoGFieldName","rho")),
rhoG_(sm.mesh().lookupObject<volScalarField> (rhoGFieldName_)),
dSauter_(sm.mesh().lookupObject<volScalarField> ("dSauter")),
alphaG_
( IOobject
(
"dSauterMix",
"alphaG",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0)
),
dHydMix_
( IOobject
(
"dHydMix",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0)
),
alphaP_
( IOobject
(
"alphaP",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0)
),
alphaSt_
( IOobject
(
"alphaSt",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
sm.mesh()
dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0)
),
alphaDyn_
( IOobject
@ -106,10 +75,22 @@ FinesFields::FinesFields
),
sm.mesh()
),
uDyn_
alphaP_
( IOobject
(
"uDyn",
"alphaP",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0)
),
alphaSt_
( IOobject
(
"alphaSt",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::MUST_READ,
@ -117,17 +98,53 @@ FinesFields::FinesFields
),
sm.mesh()
),
Sds_
dHydMix_
( IOobject
(
"Sds",
"dHydMix",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0)
dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0)
),
DragCoeff_
( IOobject
(
"DragCoeff",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0), 0)
),
dSauterMix_
( IOobject
(
"dSauterMix",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,1,0,0,0), 0)
),
FanningCoeff_
( IOobject
(
"FanningCoeff",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0), 0)
),
Froude_
( IOobject
@ -139,24 +156,59 @@ FinesFields::FinesFields
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0)
dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0)
),
Sds_
( IOobject
(
"Sds",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0)
),
uDyn_
( IOobject
(
"uDyn",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
sm.mesh()
),
estimatedVelFrac_(readScalar(propsDict_.lookup ("estimatedVelFrac"))),
dFine_("dFine", dimensionSet(0,1,0,0,0), 0.0),
depRate_(readScalar(propsDict_.lookup ("depRate"))),
rhoDyn_(readScalar(propsDict_.lookup ("rhoDyn"))),
nCrit_(readScalar(propsDict_.lookup ("nCrit"))),
nuAve_("nuAve",dimensionSet(0,2,-1,0,0),1.6e-5),
rhoFine_("rhoFine",dimensionSet(1,-3,0,0,0),0.0),
g_("g",dimensionSet(0,1,-2,0,0),vector(0,0,-9.81)),
alphaMax_(readScalar(propsDict_.lookup ("alphaMax"))),
critVoidfraction_(readScalar(propsDict_.lookup ("critVoidfraction"))),
g("g",dimensionSet(0,1,-2,0,0),9.81)
depRate_(readScalar(propsDict_.lookup ("depRate"))),
exponent_(-1.33),
nCrit_(readScalar(propsDict_.lookup ("nCrit"))),
prefactor_(14.98)
{
dFine_.value()=readScalar(propsDict_.lookup ("dFine"));
Sds_.write();
// init force sub model
// setForceSubModels(propsDict_);
// define switches which can be read from dict
if (propsDict_.found("prefactor"))
prefactor_=readScalar(propsDict_.lookup ("prefactor"));
if (propsDict_.found("exponent"))
exponent_=readScalar(propsDict_.lookup ("exponent"));
if (propsDict_.found("dFine"))
dFine_.value()=readScalar(propsDict_.lookup ("dFine"));
else
FatalError <<"Please specify dFine.\n" << abort(FatalError);
if (propsDict_.found("rhoFine"))
rhoFine_.value()=readScalar(propsDict_.lookup ("rhoFine"));
else
FatalError <<"Please specify rhoFine.\n" << abort(FatalError);
if (propsDict_.found("nuAve"))
nuAve_.value()=readScalar(propsDict_.lookup ("nuAve"));
}
@ -172,30 +224,23 @@ FinesFields::~FinesFields()
void FinesFields::update()
{
// the sequence is important: for the source terms, the Sauter mean diameter of the large particles
// is needed, for the force calculation, the static hold-up's contribution needs to be included
alphaP_ = 1.0 - voidfraction_;
updateAlphaP();
updateAlphaG();
calcSource();
updateDSauter();
dHydMix_ = 2*(1 - alphaP_ - alphaSt_) / (3*(alphaP_ + alphaSt_) ) * dSauterMix_;
Froude_ = alphaDyn_*mag(uDyn_ - UsField_) / Foam::sqrt(dHydMix_*g_);
updateDHydMix();
updateFroude();
updateFanningCoeff();
updateUDyn();
integrateFields();
updateVoidfraction();
// update voidfraction, probably really bad...
voidfraction_ = alphaG_;
}
void FinesFields::updateDSauter()
{
dSauterMix_ = (alphaP_ + alphaSt_) / ( alphaP_/dSauter_ + alphaSt_/dFine_ );
// manipulate voidfraction field
}
void FinesFields::calcSource()
{
Sds_=0;
label cellI=0;
scalar f(0.0);
scalar critpore(0.0);
scalar deltaAlpha(0.0);
@ -230,17 +275,8 @@ void FinesFields::calcSource()
Sds_[cellI] = depRate_ * deltaAlpha;
}
}
}
void FinesFields::updateUDyn()
{
// more implementation to follow
}
void FinesFields::integrateFields()
{
@ -266,21 +302,87 @@ void FinesFields::integrateFields()
alphaDynEqn.solve();
}
void FinesFields::updateVoidfraction()
void FinesFields::updateAlphaG()
{
forAll(voidfraction_,cellI)
forAll(alphaG_,cellI)
{
if(voidfraction_[cellI] - alphaSt_[cellI] > critVoidfraction_)
{
voidfraction_[cellI] -= alphaSt_[cellI];
alphaG_[cellI] = voidfraction_[cellI] - alphaSt_[cellI];
}
else
{
voidfraction_[cellI] = critVoidfraction_;
alphaG_[cellI] = critVoidfraction_;
}
}
}
void FinesFields::updateAlphaP()
{
alphaP_ = 1.0 - voidfraction_;
}
void FinesFields::updateDHydMix()
{
dHydMix_ = 2*(1 - alphaP_ - alphaSt_) / (3*(alphaP_ + alphaSt_) ) * dSauterMix_;
}
void FinesFields::updateDragCoeff()
{
volScalarField beta = 0.75 * rhoG_ * alphaDyn_ / dFine_ * mag(U_ - uDyn_) * Foam::pow(alphaG_,-4.65);
volScalarField Ref = dFine_ * alphaG_ / nuAve_ * mag(U_ - uDyn_);
scalar Cd(0.0);
scalar Ref1(0.0);
forAll(DragCoeff_,cellI)
{
Ref1 = Ref[cellI];
if(Ref1 <= SMALL)
Cd = 24.0 / SMALL;
else if(Ref1 <= 1.0)
Cd = 24.0 / Ref1;
else if(Ref1 <= 1000)
Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
else
Cd = 0.44;
DragCoeff_ = Cd * beta[cellI];
}
}
void FinesFields::updateDSauter()
{
dSauterMix_ = (alphaP_ + alphaSt_) / ( alphaP_/dSauter_ + alphaSt_/dFine_ );
}
void FinesFields::updateFanningCoeff()
{
FanningCoeff_ = alphaDyn_ * rhoFine_ * mag(uDyn_ - UsField_) * prefactor_ * Foam::pow(Froude_, exponent_) / (2 * dHydMix_);
}
void FinesFields::updateFroude()
{
Froude_ = alphaDyn_ * mag(uDyn_ - UsField_) / Foam::sqrt(dHydMix_*mag(g_));
}
void FinesFields::updateUDyn()
{
volVectorField num = rhoFine_ * alphaDyn_ * g_ + FanningCoeff_ * UsField_ + DragCoeff_ * U_;
if (p_.dimensions()==dimensionSet(0,2,-2,0,0) )
num -= alphaDyn_ * rhoG_ * fvc::grad(p_) ;
else
num -= alphaDyn_ * fvc::grad(p_) ;
volScalarField denom = FanningCoeff_ + DragCoeff_;
uDyn_ = num / denom;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

58
src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H
View File

@ -54,46 +54,78 @@ private:
const volVectorField& UsField_;
word pFieldName_;
const volScalarField& p_;
word rhoGFieldName_;
const volScalarField& rhoG_;
const volScalarField& dSauter_;
volScalarField dSauterMix_;
volScalarField alphaG_;
volScalarField dHydMix_;
volScalarField alphaDyn_;
volScalarField alphaP_;
volScalarField alphaSt_;
volScalarField alphaDyn_;
volScalarField dHydMix_;
volVectorField uDyn_;
volScalarField DragCoeff_;
volScalarField Sds_;
volScalarField dSauterMix_;
volScalarField FanningCoeff_;
volScalarField Froude_;
volScalarField Sds_;
volVectorField uDyn_;
dimensionedScalar dFine_;
scalar rhoDyn_;
dimensionedScalar nuAve_;
scalar depRate_;
dimensionedScalar rhoFine_;
scalar nCrit_;
const dimensionedVector g_;
scalar alphaMax_;
scalar critVoidfraction_;
const dimensiondScalar g_;
scalar depRate_;
scalar exponent_;
scalar nCrit_;
scalar prefactor_;
void calcSource();
void integrateFields();
void updateAlphaG();
void updateAlphaP();
void updateDHydMix();
void updateDragCoeff();
void updateDSauter();
void updateFanningCoeff();
void updateFroude();
void updateUDyn();
void integrateFields();
public:
//- Runtime type information
@ -111,7 +143,7 @@ public:
// Destructor
~FinesFields();
virtual ~FinesFields();
// Member Functions

218
src/lagrangian/cfdemParticle/subModels/forceModel/Fines/gradPStatFines.C Normal file
View File

@ -0,0 +1,218 @@
/*---------------------------------------------------------------------------*\
License
This 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.
This code 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 this code. If not, see <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
Class
ErgunStatFines
SourceFiles
ErgunStatFines.C
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "gradPStatFines.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(gradPStatFines, 0);
addToRunTimeSelectionTable
(
forceModel,
gradPStatFines,
dictionary
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
gradPStatFines::gradPStatFines
(
const dictionary& dict,
cfdemCloud& sm
)
:
forceModel(dict,sm),
propsDict_(dict.subDict(typeName + "Props")),
pFieldName_(propsDict_.lookup("pFieldName")),
p_(sm.mesh().lookupObject<volScalarField> (pFieldName_)),
velocityFieldName_(propsDict_.lookup("velocityFieldName")),
U_(sm.mesh().lookupObject<volVectorField> (velocityFieldName_)),
alphaP_(sm.mesh().lookupObject<volScalarField> ("alphaP")),
alphaSt_(sm.mesh().lookupObject<volScalarField> ("alphaSt")),
useRho_(false),
useU_(false),
addedMassCoeff_(0.0)
{
// init force sub model
setForceSubModels(propsDict_);
// define switches which can be read from dict
forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
forceSubM(0).setSwitchesList(1,true); // activate treatForceDEM switch
forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
// read those switches defined above, if provided in dict
forceSubM(0).readSwitches();
if (modelType_ == "B")
{
FatalError <<"using model gradPStatFines with model type B is not valid\n" << abort(FatalError);
}else if (modelType_ == "Bfull")
{
if(forceSubM(0).switches()[1])
{
Info << "Using treatForceDEM false!" << endl;
forceSubM(0).setSwitches(1,false); // treatForceDEM = false
}
}else // modelType_=="A"
{
if(!forceSubM(0).switches()[1])
{
Info << "Using treatForceDEM true!" << endl;
forceSubM(0).setSwitches(1,true); // treatForceDEM = true
}
}
if (propsDict_.found("useU")) useU_=true;
if (propsDict_.found("useAddedMass"))
{
addedMassCoeff_ = readScalar(propsDict_.lookup("useAddedMass"));
Info << "gradP will also include added mass with coefficient: " << addedMassCoeff_ << endl;
Info << "WARNING: use fix nve/sphere/addedMass in LIGGGHTS input script to correctly account for added mass effects!" << endl;
}
if(p_.dimensions()==dimensionSet(0,2,-2,0,0))
useRho_ = true;
particleCloud_.checkCG(true);
particleCloud_.probeM().initialize(typeName, "gradP.logDat");
particleCloud_.probeM().vectorFields_.append("gradPStatFines"); //first entry must the be the force
particleCloud_.probeM().scalarFields_.append("Vs");
particleCloud_.probeM().scalarFields_.append("rho");
particleCloud_.probeM().writeHeader();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
gradPStatFines::~gradPStatFines()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void gradPStatFines::setForce() const
{
volVectorField gradPField = fvc::grad(p_);
/*if (useU_)
{
// const volScalarField& voidfraction_ = particleCloud_.mesh().lookupObject<volScalarField> ("voidfraction");
volScalarField U2 = U_&U_;// *voidfraction_*voidfraction_;
if (useRho_)
gradPField = fvc::grad(0.5*U2);
else
gradPField = fvc::grad(0.5*forceSubM(0).rhoField()*U2);
}*/
vector gradP;
scalar Vs;
scalar rho;
vector position;
vector force;
label cellI;
scalar volcorr(0.0);
interpolationCellPoint<vector> gradPInterpolator_(gradPField);
#include "setupProbeModel.H"
for(int index = 0;index < particleCloud_.numberOfParticles(); index++)
{
//if(mask[index][0])
//{
force=vector(0,0,0);
cellI = particleCloud_.cellIDs()[index][0];
if (cellI > -1) // particle Found
{
position = particleCloud_.position(index);
if(forceSubM(0).interpolation()) // use intepolated values for alpha (normally off!!!)
{
gradP = gradPInterpolator_.interpolate(position,cellI);
}else
{
gradP = gradPField[cellI];
}
Vs = particleCloud_.particleVolume(index);
rho = forceSubM(0).rhoField()[cellI];
// calc particle's pressure gradient force
if (useRho_)
force = -Vs*gradP*rho*(1.0+addedMassCoeff_);
else
force = -Vs*gradP*(1.0+addedMassCoeff_);
if( alphaP_[cellI] < SMALL )
{
volcorr = 1.0;
}
else
{
volcorr = (alphaP_[cellI] + alphaSt_[cellI]) / alphaP_[cellI];
}
force *= volcorr;
if(forceSubM(0).verbose() && index >=0 && index <2)
{
Info << "index = " << index << endl;
Info << "gradP = " << gradP << endl;
Info << "force = " << force << endl;
}
//Set value fields and write the probe
if(probeIt_)
{
#include "setupProbeModelfields.H"
vValues.append(force); //first entry must the be the force
sValues.append(Vs);
sValues.append(rho);
particleCloud_.probeM().writeProbe(index, sValues, vValues);
}
}
// write particle based data to global array
forceSubM(0).partToArray(index,force,vector::zero);
//}
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

98
src/lagrangian/cfdemParticle/subModels/forceModel/Fines/gradPStatFines.H Normal file
View File

@ -0,0 +1,98 @@
/*---------------------------------------------------------------------------*\
License
This 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.
This code 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 this code. If not, see <http://www.gnu.org/licenses/>.
Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
Class
gradPStatFines
SourceFiles
gradPStatFines.C
\*---------------------------------------------------------------------------*/
#ifndef gradPStatFines_H
#define gradPStatFines_H
#include "forceModel.H"
#include "interpolationCellPoint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class gradPStatFines Declaration
\*---------------------------------------------------------------------------*/
class gradPStatFines
:
public forceModel
{
private:
dictionary propsDict_;
word pFieldName_;
const volScalarField& p_;
word velocityFieldName_;
const volVectorField& U_;
const volScalarField& alphaP_;
const volScalarField& alphaSt_;
bool useRho_;
bool useU_; // if false: substitution p=0.5*rho*U^2
mutable double addedMassCoeff_; //added mass coefficient
public:
//- Runtime type information
TypeName("gradPStatFines");
// Constructors
//- Construct from components
gradPStatFines
(
const dictionary& dict,
cfdemCloud& sm
);
// Destructor
~gradPStatFines();
// Member Functions
void setForce() const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //