zhyang-dev/CFDEMcoupling-PFM
1
0
Fork 0
mirror of https://github.com/ParticulateFlow/CFDEMcoupling-PFM.git synced 2025-12-08 06:37:44 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge pull request #99 from ParticulateFlow/feature/fines

Browse Source 
Feature/fines
This commit is contained in:
tlichtenegger
2020-04-23 16:02:54 +02:00
committed by GitHub
parent 63f1f8fb68 52bcfc5dd3
commit 116fd3b810
18 changed files with 384 additions and 161 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
View File

@ -92,7 +92,7 @@ int main(int argc, char *argv[])
particleCloud.clockM().start(2,"Coupling");
bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
if(hasEvolved)
if(hasEvolved && smoothenForces)
{
particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
}

11
applications/solvers/cfdemSolverRhoPimple/createFields.H
View File

@ -177,6 +177,17 @@ Info<< "Reading thermophysical properties\n" << endl;
)
);
bool smoothenForces
(
pimple.dict().lookupOrDefault<bool>
(
"smoothenForces",
false
)
);
if (smoothenForces) Info << "Smoothening implicit particle forces.\n" << endl;
else Info << "Not smoothening implicit particle forces.\n" << endl;
Info<< "Creating turbulence model\n" << endl;
autoPtr<compressible::turbulenceModel> turbulence
(

38
src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C
View File

@ -68,6 +68,18 @@ twoWayMPI::twoWayMPI
:
dataExchangeModel(dict,sm),
propsDict_(dict.subDict(typeName + "Props")),
DEMVariableNames_(propsDict_.lookupOrDefault<wordList>("DEMVariables",wordList(0))),
DEMVariables_
(
IOobject
(
"DEMVariables",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
)
),
lmp(NULL)
{
Info << "Starting up LIGGGHTS for first time execution" << endl;
@ -100,7 +112,30 @@ twoWayMPI::~twoWayMPI()
delete lmp;
}
// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * //
void twoWayMPI::updateDEMVariables()
{
scalar variablevalue = 0.0;
word variablename("");
DEMVariables_.clear();
for (label i=0; i<DEMVariableNames_.size(); i++)
{
variablename = DEMVariableNames_[i];
variablevalue = DEMVariableValue(variablename);
DEMVariables_.append(variablevalue);
}
}
double twoWayMPI::DEMVariableValue(word variablename)
{
return liggghts_get_variable(lmp,variablename.c_str());
}
// * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * //
void twoWayMPI::getData
(
word name,
@ -355,6 +390,9 @@ bool twoWayMPI::couple(int i)
particleCloud_.clockM().start(4,"LIGGGHTS_reallocArrays");
particleCloud_.reAllocArrays();
particleCloud_.clockM().stop("LIGGGHTS_reallocArrays");
// retrieve DEM variables if present
updateDEMVariables();
}
return coupleNow;

9
src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H
View File

@ -44,6 +44,7 @@ SourceFiles
#include <mpi.h>
#include "dataExchangeModel.H"
#include "scalarIOList.H"
//=================================//
//LAMMPS/LIGGGHTS
@ -71,7 +72,14 @@ private:
dictionary propsDict_;
MPI_Comm comm_liggghts;
wordList DEMVariableNames_;
scalarIOList DEMVariables_;
// private member functions
double DEMVariableValue(word);
void updateDEMVariables();
protected:
LAMMPS_NS::LAMMPS *lmp;
@ -146,7 +154,6 @@ public:
int getNumberOfTypes() const;
double* getTypeVol() const;
scalar getCG() const { return lmp->force->cg(); }
};

2
src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C
View File

@ -218,8 +218,6 @@ void BeetstraDrag::setForce() const
voidfraction = voidfraction_[cellI];
Ufluid = U_[cellI];
}
// in case a fines phase is present, void fraction needs to be adapted
adaptVoidfraction(voidfraction, cellI);
if (multiTypes_)
{

2
src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H
View File

@ -90,8 +90,6 @@ protected:
bool usePC_;
virtual void adaptVoidfraction(double&, label) const {}
virtual scalar effDiameter(double d, label cellI, label index) const {return d;}
virtual scalar meanSauterDiameter(double d, label cellI) const {return d;}

11
src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C
View File

@ -59,19 +59,15 @@ BeetstraDragPoly::BeetstraDragPoly
{
dFine_ = readScalar(propsDict_.lookup("dFine"));
volScalarField& alphaP(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaP")));
// alphaP_.set(&alphaP);
alphaP_ = &alphaP;
volScalarField& alphaSt(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaSt")));
// alphaSt_.set(&alphaSt);
alphaSt_ = &alphaSt;
volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauterMix")));
// dSauter_.set(&dSauter);
dSauter_ = &dSauter;
}
else
{
volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauter")));
// dSauter_.set(&dSauter);
dSauter_ = &dSauter;
}
}
@ -84,13 +80,6 @@ BeetstraDragPoly::~BeetstraDragPoly()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void BeetstraDragPoly::adaptVoidfraction(double& voidfraction, label cellI) const
{
if (fines_) voidfraction -= (*alphaSt_)[cellI];
if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
}
scalar BeetstraDragPoly::effDiameter(double d, label cellI, label index) const
{
scalar dS = (*dSauter_)[cellI];

2
src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H
View File

@ -52,8 +52,6 @@ protected:
volScalarField* dSauter_;
void adaptVoidfraction(double&, label) const;
scalar effDiameter(double, label, label) const;
scalar meanSauterDiameter(double, label) const;

5
src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.C
View File

@ -50,6 +50,11 @@ Fines::Fines
forceModel(dict,sm),
finesFields_(dict,sm)
{
// safety check: fines should be first force model in list because they correct the voidfraction field
if (particleCloud_.forceModels()[0] != "dSauter" || particleCloud_.forceModels()[1] != "Fines")
{
FatalError <<"Force models dSauter and Fines need to be first in list.\n" << abort(FatalError);
}
}

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

@ -45,11 +45,15 @@ FinesFields::FinesFields
propsDict_(dict.subDict("FinesFieldsProps")),
smoothing_(propsDict_.lookupOrDefault<bool>("smoothing",false)),
verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
clogKin_(propsDict_.lookupOrDefault<bool>("kineticClogging",false)),
clogStick_(propsDict_.lookupOrDefault<bool>("stickyClogging",false)),
movingBed_(propsDict_.lookupOrDefault<bool>("movingBed",true)),
fluxFieldName_(propsDict_.lookupOrDefault<word>("fluxFieldName","phi")),
phi_(sm.mesh().lookupObject<surfaceScalarField> (fluxFieldName_)),
velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
// this is probably really bad
voidfraction_(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_))),
voidfraction_(sm.mesh().lookupObjectRef<volScalarField> (voidfractionFieldName_)),
UsFieldName_(propsDict_.lookupOrDefault<word>("granVelFieldName","Us")),
UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
pFieldName_(propsDict_.lookupOrDefault<word>("pFieldName","p")),
@ -200,7 +204,6 @@ FinesFields::FinesFields
),
sm.mesh(),
dimensionedScalar("zero", dimensionSet(1,0,-1,0,0), 0)
//dimensionedVector("zero", dimensionSet(1,-2,-1,0,0), vector::zero)
),
uDyn_
( IOobject
@ -219,54 +222,114 @@ FinesFields::FinesFields
rhoFine_("rhoFine",dimensionSet(1,-3,0,0,0),0.0),
g_("g",dimensionSet(0,1,-2,0,0),vector(0,0,-9.81)),
alphaDynMax_(0.1),
alphaMax_(readScalar(propsDict_.lookup ("alphaMax"))),
critVoidfraction_(readScalar(propsDict_.lookup ("critVoidfraction"))),
depRate_(readScalar(propsDict_.lookup ("depRate"))),
alphaMax_(propsDict_.lookupOrDefault<scalar>("alphaMax",0.95)),
alphaMinClog_(propsDict_.lookupOrDefault<scalar>("alphaMinClog",0.3)),
critVoidfraction_(propsDict_.lookupOrDefault<scalar>("critVoidfraction", 0.05)),
deltaT_(voidfraction_.mesh().time().deltaTValue()),
depositionLength_(readScalar(propsDict_.lookup ("depositionLength"))),
exponent_(-1.33),
nCrit_(readScalar(propsDict_.lookup ("nCrit"))),
poresizeWidth_(readScalar(propsDict_.lookup ("poresizeWidth"))),
nCrit_(0.0),
poresizeWidth_(0.0),
prefactor_(10.5),
ratioHydraulicPore_(1.5)
ratioHydraulicPore_(1.5),
tauRelease_(readScalar(propsDict_.lookup ("tauRelease"))),
uBindHigh_(propsDict_.lookupOrDefault<scalar>("uBindHigh",0.0)),
uBindLow_(propsDict_.lookupOrDefault<scalar>("uBindLow",0.0)),
uMin_(0.001)
{
Sds_.write();
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("diffCoeff"))
diffCoeff_.value()=readScalar(propsDict_.lookup ("diffCoeff"));
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"));
if (propsDict_.found("alphaDynMax"))
alphaDynMax_=readScalar(propsDict_.lookup ("alphaDynMax"));
if(verbose_)
{
alphaG_.writeOpt() = IOobject::AUTO_WRITE;
alphaG_.write();
alphaP_.writeOpt() = IOobject::AUTO_WRITE;
alphaP_.write();
deltaAlpha_.writeOpt() = IOobject::AUTO_WRITE;
deltaAlpha_.write();
dHydMix_.writeOpt() = IOobject::AUTO_WRITE;
dHydMix_.write();
DragCoeff_.writeOpt() = IOobject::AUTO_WRITE;
DragCoeff_.write();
dSauterMix_.writeOpt() = IOobject::AUTO_WRITE;
dSauterMix_.write();
FanningCoeff_.writeOpt() = IOobject::AUTO_WRITE;
FanningCoeff_.write();
Froude_.writeOpt() = IOobject::AUTO_WRITE;
Froude_.write();
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("diffCoeff"))
{
diffCoeff_.value()=readScalar(propsDict_.lookup ("diffCoeff"));
}
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"));
}
if (propsDict_.found("alphaDynMax"))
{
alphaDynMax_=readScalar(propsDict_.lookup ("alphaDynMax"));
}
if (clogKin_)
{
nCrit_ = readScalar(propsDict_.lookup ("nCrit"));
poresizeWidth_ = readScalar(propsDict_.lookup ("poresizeWidth"));
}
if (clogStick_)
{
if (uBindHigh_ - uBindLow_ < SMALL)
{
FatalError <<"No reasonable values for critical binding velocities.\n" << abort(FatalError);
}
uReconstructed_.set
(
new volVectorField
(
IOobject
(
"uReconstructed",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
U_//sm.mesh()
)
);
}
if (verbose_)
{
alphaG_.writeOpt() = IOobject::AUTO_WRITE;
alphaG_.write();
alphaP_.writeOpt() = IOobject::AUTO_WRITE;
alphaP_.write();
deltaAlpha_.writeOpt() = IOobject::AUTO_WRITE;
deltaAlpha_.write();
dHydMix_.writeOpt() = IOobject::AUTO_WRITE;
dHydMix_.write();
DragCoeff_.writeOpt() = IOobject::AUTO_WRITE;
DragCoeff_.write();
dSauterMix_.writeOpt() = IOobject::AUTO_WRITE;
dSauterMix_.write();
FanningCoeff_.writeOpt() = IOobject::AUTO_WRITE;
FanningCoeff_.write();
Froude_.writeOpt() = IOobject::AUTO_WRITE;
Froude_.write();
if (clogStick_)
{
uReconstructed_().writeOpt() = IOobject::AUTO_WRITE;
uReconstructed_().write();
}
}
}
@ -281,76 +344,93 @@ FinesFields::~FinesFields()
void FinesFields::update()
{
if(verbose_) Info << "FinesFields: Updating alphaP.\n" << endl;
if (verbose_) Info << "FinesFields: Updating alphaP.\n" << endl;
updateAlphaP();
if(verbose_) Info << "FinesFields: Updating alphaG.\n" << endl;
if (verbose_) Info << "FinesFields: Updating alphaG.\n" << endl;
updateAlphaG();
if(verbose_) Info << "FinesFields: Calculating source terms.\n" << endl;
if (clogStick_)
{
if (verbose_) Info << "FinesFields: Updating uReconstructed.\n" << endl;
updateUReconstructed();
}
if (verbose_) Info << "FinesFields: Calculating source terms.\n" << endl;
calcSource();
if(verbose_) Info << "FinesFields: Updating dSauter.\n" << endl;
if (verbose_) Info << "FinesFields: Updating dSauter.\n" << endl;
updateDSauter();
if(verbose_) Info << "FinesFields: Updating dHydMix.\n" << endl;
if (verbose_) Info << "FinesFields: Updating dHydMix.\n" << endl;
updateDHydMix();
if(verbose_) Info << "FinesFields: Updating Froude.\n" << endl;
if (verbose_) Info << "FinesFields: Updating Froude.\n" << endl;
updateFroude();
if(verbose_) Info << "FinesFields: Updating FanningCoeff.\n" << endl;
if (verbose_) Info << "FinesFields: Updating FanningCoeff.\n" << endl;
updateFanningCoeff();
if(verbose_) Info << "FinesFields: Updating DragCoeff.\n" << endl;
if (verbose_) Info << "FinesFields: Updating DragCoeff.\n" << endl;
updateDragCoeff();
if(verbose_) Info << "FinesFields: Updating uDyn.\n" << endl;
if (verbose_) Info << "FinesFields: Updating uDyn.\n" << endl;
updateUDyn();
if(verbose_) Info << "FinesFields: Integrating alphas.\n" << endl;
if (verbose_) Info << "FinesFields: Integrating alphas.\n" << endl;
integrateFields();
if(verbose_) Info << "FinesFields: Update finished.\n" << endl;
if (verbose_) Info << "FinesFields: Update finished.\n" << endl;
}
void FinesFields::calcSource()
{
if (!clogKin_ & !clogStick_) return;
Sds_.primitiveFieldRef() = 0;
deltaAlpha_.primitiveFieldRef() = 0.0;
scalar f(0.0);
scalar critpore(0.0);
scalar dmean(0.0);
scalar d1(0.0);
scalar d2(0.0);
scalar fKin = 0.0;
scalar fStick = 0.0;
scalar critpore = 0.0;
scalar dmean = 0.0;
scalar d1 = 0.0;
scalar d2 = 0.0;
scalar magU = 0.0;
scalar tauDeposition = 0.0;
forAll(Sds_,cellI)
{
// calculate everything in units auf dSauter
critpore = nCrit_*dFine_.value()/dSauter_[cellI];
// pore size from hydraulic radius
dmean = 2 * (1 - alphaP_[cellI]) / ( (1 + poresizeWidth_*poresizeWidth_/3) * 3 * alphaP_[cellI] );
// Sweeney and Martin, Acta Materialia 51 (2003): ratio of hydraulic to pore throat radius
dmean /= ratioHydraulicPore_;
d1 = dmean * (1 - poresizeWidth_);
d2 = dmean * (1 + poresizeWidth_);
fKin = 0.0;
fStick = 0.0;
if (clogKin_ && alphaP_[cellI] > alphaMinClog_) // no kinetic cloggig in dilute regions
{
// calculate everything in units auf dSauter
critpore = nCrit_*dFine_.value()/dSauter_[cellI];
// pore size from hydraulic radius
dmean = 2 * (1 - alphaP_[cellI]) / ( (1 + poresizeWidth_*poresizeWidth_/3) * 3 * alphaP_[cellI] );
// Sweeney and Martin, Acta Materialia 51 (2003): ratio of hydraulic to pore throat radius
dmean /= ratioHydraulicPore_;
d1 = dmean * (1 - poresizeWidth_);
d2 = dmean * (1 + poresizeWidth_);
f = (critpore*critpore*critpore - d1 * d1 * d1) / (d2 * d2 * d2 - d1 * d1 * d1);
if (f < 0)
{
f = 0.0;
fKin = (critpore*critpore*critpore - d1 * d1 * d1) / (d2 * d2 * d2 - d1 * d1 * d1);
if (fKin < 0) fKin = 0.0;
else if (fKin > 1.0) fKin = 1.0;
}
else if (f > 1.0)
if (clogStick_)
{
f = 1.0;
magU = mag(uReconstructed_()[cellI]); // use U reconstructed from phi to suppress oscillations at interfaces
// fStick = 1.0 / ( 1.0 + magU/uBind_) * alphaP_[cellI] / 0.65;
if (magU < uBindLow_) fStick = 1.0;
else if (magU > uBindHigh_) fStick = 0.0;
else fStick = 1.0 - (magU - uBindLow_) / (uBindHigh_ - uBindLow_);
// if (fStick > 1.0) fStick = 1.0;
fStick *= alphaP_[cellI] / 0.65;
}
// at this point, voidfraction is still calculated from the true particle sizes
deltaAlpha_[cellI] = f * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI];
// too much volume occupied: release it (50% per time step)
deltaAlpha_[cellI] = max(fKin,fStick) * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI];
// too much volume occupied: release it
if (deltaAlpha_[cellI] < 0.0)
{
Sds_[cellI] = 0.5*deltaAlpha_[cellI];
}
// volume too occupy available: deposit at most 80% of dyn hold up
else if (depRate_ * deltaAlpha_[cellI] > 0.8 * alphaDyn_[cellI])
{
Sds_[cellI] = 0.8 * alphaDyn_[cellI];
Sds_[cellI] = deltaAlpha_[cellI] / tauRelease_;
}
// volume to occupy available
else
{
Sds_[cellI] = depRate_ * deltaAlpha_[cellI];
tauDeposition = depositionLength_ / max(mag(U_[cellI]),uMin_);
if (tauDeposition < 4.0*deltaT_) tauDeposition = 4.0*deltaT_; // do not deposit all dyn hold-up within less than 3 time steps
Sds_[cellI] = alphaDyn_[cellI] / tauDeposition;
}
}
}
@ -364,11 +444,12 @@ void FinesFields::integrateFields()
fvScalarMatrix alphaStEqn
(
fvm::ddt(alphaSt_)
+ fvm::div(phiSt,alphaSt_)
fvm::ddt(alphaSt_)
==
Sds_
);
if (movingBed_) alphaStEqn += fvm::div(phiSt,alphaSt_);
fvScalarMatrix alphaDynEqn
(
fvm::ddt(alphaDyn_)
@ -377,10 +458,11 @@ void FinesFields::integrateFields()
==
-Sds_
);
alphaStEqn.solve();
alphaDynEqn.solve();
if(smoothing_)
if (smoothing_)
particleCloud_.smoothingM().smoothen(alphaDyn_);
// limit hold-ups, should be done more elegantly
@ -422,13 +504,14 @@ void FinesFields::integrateFields()
}
void FinesFields::updateAlphaG()
void FinesFields::updateAlphaG() // called after updateAlphaP() - correct voidfraction by removing space occupied by fines
{
alphaG_ = max(voidfraction_ - alphaSt_ - alphaDyn_, critVoidfraction_);
voidfraction_ = alphaG_;
}
void FinesFields::updateAlphaP()
void FinesFields::updateAlphaP() // called first in the update cycle - voidfraction_ is current with DEM data
{
alphaP_ = 1.0 - voidfraction_ + SMALL;
}
@ -439,7 +522,7 @@ void FinesFields::updateDHydMix()
forAll(dHydMix_,cellI)
{
scalar aPSt = alphaP_[cellI] + alphaSt_[cellI];
if(aPSt < SMALL || aPSt > 1 - SMALL)
if (aPSt < SMALL || aPSt > 1 - SMALL)
dHydMix_[cellI] = SMALL;
else
dHydMix_[cellI] = 2*(1 - aPSt) / (3*aPSt ) * dSauterMix_[cellI];
@ -459,11 +542,11 @@ void FinesFields::updateDragCoeff()
forAll(DragCoeff_,cellI)
{
Ref1 = Ref[cellI];
if(Ref1 <= SMALL)
if (Ref1 <= SMALL)
Cd = 24.0 / SMALL;
else if(Ref1 <= 1.0)
else if (Ref1 <= 1.0)
Cd = 24.0 / Ref1;
else if(Ref1 <= 1000)
else if (Ref1 <= 1000)
Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
else
Cd = 0.44;
@ -475,11 +558,11 @@ void FinesFields::updateDragCoeff()
forAll(DragCoeff_.boundaryField()[patchI], faceI)
{
Ref1 = Ref.boundaryField()[patchI][faceI];
if(Ref1 <= SMALL)
if (Ref1 <= SMALL)
Cd = 24.0 / SMALL;
else if(Ref1 <= 1.0)
else if (Ref1 <= 1.0)
Cd = 24.0 / Ref1;
else if(Ref1 <= 1000)
else if (Ref1 <= 1000)
Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
else
Cd = 0.44;
@ -496,9 +579,9 @@ void FinesFields::updateDSauter()
{
scalar aP = alphaP_[cellI];
scalar aSt = alphaSt_[cellI];
if(aSt < SMALL)
if (aSt < SMALL)
dSauterMix_[cellI] = dSauter_[cellI];
else if(aP < SMALL)
else if (aP < SMALL)
dSauterMix_[cellI] = dFine_.value();
else
dSauterMix_[cellI] = (aP + aSt) / (aP / dSauter_[cellI] + aSt / dFine_.value() );
@ -541,7 +624,7 @@ void FinesFields::updateUDyn()
{
scalar mU(mag(U_[cellI]));
scalar muDyn(mag(uDyn_[cellI]));
if(muDyn > mU && muDyn > SMALL)
if (muDyn > mU && muDyn > SMALL)
{
uDyn_[cellI] *= mU / muDyn;
}
@ -552,13 +635,24 @@ void FinesFields::updateUDyn()
{
scalar mU(mag(U_.boundaryField()[patchI][faceI]));
scalar muDyn(mag(uDyn_.boundaryField()[patchI][faceI]));
if(muDyn > mU && muDyn > SMALL)
if (muDyn > mU && muDyn > SMALL)
{
uDyn_.boundaryFieldRef()[patchI][faceI] *= mU / muDyn;
}
}
}
void FinesFields::updateUReconstructed()
{
if (phi_.dimensions() == dimensionSet(1, 0, -1, 0, 0)) // compressible
{
uReconstructed_() = fvc::reconstruct(phi_) / (rhoG_ * voidfraction_);
}
else
{
uReconstructed_() = fvc::reconstruct(phi_) / voidfraction_;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

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

@ -46,6 +46,16 @@ private:
bool verbose_;
bool clogKin_;
bool clogStick_;
bool movingBed_;
word fluxFieldName_;
const surfaceScalarField& phi_;
word velFieldName_;
const volVectorField& U_;
@ -90,11 +100,12 @@ private:
volScalarField Sds_;
//volVectorField massFluxDyn_;
surfaceScalarField massFluxDyn_;
volVectorField uDyn_;
autoPtr<volVectorField> uReconstructed_; // velocity field can show oscillations at porosity steps
dimensionedScalar dFine_;
dimensionedScalar diffCoeff_;
@ -109,9 +120,13 @@ private:
scalar alphaMax_;
scalar alphaMinClog_;
scalar critVoidfraction_;
scalar depRate_;
scalar deltaT_;
scalar depositionLength_;
scalar exponent_;
@ -123,6 +138,14 @@ private:
scalar ratioHydraulicPore_;
scalar tauRelease_;
scalar uBindHigh_;
scalar uBindLow_;
scalar uMin_;
void calcSource();
void integrateFields();
@ -143,6 +166,8 @@ private:
void updateUDyn();
void updateUReconstructed();
public:
//- Runtime type information

6
src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C
View File

@ -69,7 +69,8 @@ forceModel::forceModel
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero) // N
dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero),
"zeroGradient"
),
expParticleForces_
( IOobject
@ -81,7 +82,8 @@ forceModel::forceModel
IOobject::AUTO_WRITE
),
sm.mesh(),
dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero) // N
dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero),
"zeroGradient"
),
coupleForce_(true),
modelType_(sm.modelType()),

79
src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C
View File

@ -65,14 +65,70 @@ constDiffSmoothing::constDiffSmoothing
propsDict_(dict.subDict(typeName + "Props")),
lowerLimit_(readScalar(propsDict_.lookup("lowerLimit"))),
upperLimit_(readScalar(propsDict_.lookup("upperLimit"))),
smoothingLength_(dimensionedScalar("smoothingLength", dimLength, readScalar(propsDict_.lookup("smoothingLength")))),
smoothingLengthReferenceField_(dimensionedScalar("smoothingLengthReferenceField", dimLength, readScalar(propsDict_.lookup("smoothingLength")))),
DT_("DT", dimensionSet(0,2,-1,0,0), 0.),
smoothingLength_(propsDict_.lookupOrDefault<scalar>("smoothingLength", -1.0)),
smoothingLengthReference_(propsDict_.lookupOrDefault<scalar>("smoothingLengthReference",smoothingLength_)),
smoothingLengthFieldName_(propsDict_.lookupOrDefault<word>("smoothingLengthFieldName","smoothingLengthField")),
smoothingLengthField_
( IOobject
(
smoothingLengthFieldName_,
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("smoothingLength", dimensionSet(0,1,0,0,0,0,0), smoothingLength_),
"zeroGradient"
),
smoothingLengthReferenceFieldName_(propsDict_.lookupOrDefault<word>("smoothingLengthReferenceFieldName","smoothingLengthReferenceField")),
smoothingLengthReferenceField_
( IOobject
(
smoothingLengthReferenceFieldName_,
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("smoothingLengthReference", dimensionSet(0,1,0,0,0,0,0), smoothingLengthReference_),
"zeroGradient"
),
DT_
( IOobject
(
"DT",
sm.mesh().time().timeName(),
sm.mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
sm.mesh(),
dimensionedScalar("DT", dimensionSet(0,2,-1,0,0,0,0), 0.0),
"zeroGradient"
),
verbose_(propsDict_.found("verbose"))
{
if(propsDict_.found("smoothingLengthReferenceField"))
// either use scalar or field parameters for smoothing
if (smoothingLength_ > 0.0 || smoothingLengthReference_ > 0.0)
{
smoothingLengthReferenceField_.value() = double(readScalar(propsDict_.lookup("smoothingLengthReferenceField")));
if (smoothingLengthField_.headerOk() || smoothingLengthReferenceField_.headerOk())
{
FatalError <<"constDiffSmoothing: Either use scalar or field parameter for smoothing.\n" << abort(FatalError);
}
}
if (smoothingLength_ < 0.0 && !smoothingLengthField_.headerOk())
{
FatalError <<"constDiffSmoothing: Provide scalar or field parameter for smoothing.\n" << abort(FatalError);
}
// if no scalar length for smoothing wrt reference field is provided and no
// such field, use smoothingLengthField
if (smoothingLengthReference_ < 0.0 && !smoothingLengthReferenceField_.headerOk())
{
smoothingLengthReferenceField_ = smoothingLengthField_;
}
checkFields(sSmoothField_);
@ -104,8 +160,8 @@ void constDiffSmoothing::smoothen(volScalarField& fieldSrc) const
sSmoothField.oldTime()=fieldSrc;
sSmoothField.oldTime().correctBoundaryConditions();
double deltaT = sSmoothField.mesh().time().deltaTValue();
DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT;
dimensionedScalar deltaT = sSmoothField.mesh().time().deltaT();
DT_ = smoothingLengthField_ * smoothingLengthField_ / deltaT;
// do smoothing
solve
@ -145,8 +201,8 @@ void constDiffSmoothing::smoothen(volVectorField& fieldSrc) const
vSmoothField.oldTime()=fieldSrc;
vSmoothField.oldTime().correctBoundaryConditions();
double deltaT = vSmoothField_.mesh().time().deltaTValue();
DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT;
dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
DT_ = smoothingLengthField_ * smoothingLengthField_ / deltaT;
// do smoothing
solve
@ -183,9 +239,8 @@ void constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc) const
double sourceStrength = 1e5; //large number to keep reference values constant
dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
DT_.value() = smoothingLengthReferenceField_.value()
* smoothingLengthReferenceField_.value() / deltaT.value();
DT_ = smoothingLengthReferenceField_ * smoothingLengthReferenceField_ / deltaT;
tmp<volScalarField> NLarge
(
new volScalarField

20
src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.H
View File

@ -60,11 +60,25 @@ class constDiffSmoothing
private:
dictionary propsDict_;
const scalar lowerLimit_;
const scalar upperLimit_;
const dimensionedScalar smoothingLength_;
dimensionedScalar smoothingLengthReferenceField_;
mutable dimensionedScalar DT_;
const scalar smoothingLength_;
const scalar smoothingLengthReference_;
word smoothingLengthFieldName_;
mutable volScalarField smoothingLengthField_;
word smoothingLengthReferenceFieldName_;
mutable volScalarField smoothingLengthReferenceField_;
mutable volScalarField DT_;
const bool verbose_;
public:

11
tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/constant/couplingProperties
View File

@ -38,7 +38,7 @@ meshMotionModel noMeshMotion;
regionModel allRegion;
IOModel "basicIO";
IOModel "off";
dataExchangeModel twoWayMPI;//twoWayM2M;//twoWayFiles;//oneWayVTK;//
@ -55,7 +55,7 @@ forceModels
dSauter
Fines
FanningDynFines
ErgunStatFines
BeetstraDragPoly
gradPForce
viscForce
);
@ -140,13 +140,14 @@ viscForceProps
interpolation;
}
ErgunStatFinesProps
BeetstraDragProps
{
fines true;
dFine 0.000388;
velFieldName "U";
granVelFieldName "Us";
densityFieldName "rho";
voidfractionFieldName "voidfraction";
phi 1;
}
FanningDynFinesProps
@ -207,7 +208,7 @@ dividedProps
{
alphaMin 0.25;
scaleUpVol 1.0;
weight 1.0; //1.33;
weight 1.0;
verbose;
}

2
tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/blockMeshDict
View File

@ -13,7 +13,7 @@ FoamFile
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "../../geometry"
#include "../geometry"
convertToMeters 1;
vertices #codeStream

4
tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/controlDict
View File

@ -59,7 +59,7 @@ functions
volInt
{
type volRegion;
type volFieldValue;
functionObjectLibs ("libfieldFunctionObjects.so");
writeControl timeStep;
log true;
@ -76,7 +76,7 @@ functions
inflow
{
type surfaceRegion;
type surfaceFieldValue;
libs ("libfieldFunctionObjects.so");
writeControl timeStep;
log true;

20
tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/fvOptions
View File

@ -15,20 +15,6 @@ FoamFile
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
/*
source1
{
type temperatureLimitsConstraint;
selectionMode all;
active true;
temperatureLimitsConstraintCoeffs
{
Tmin 288;
Tmax 298;
}
}
*/
source1
{
type limitTemperature;
@ -37,7 +23,9 @@ source1
{
active yes;
selectionMode all;
Tmin 288;
Tmax 298;
Tmin 288; // OF4
Tmax 298; // OF4
min $Tmin; // OF5,OF6
max $Tmax; // OF5,OF6
}
}