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

Minor changes.

Browse Source
This commit is contained in:
Thomas Lichtenegger
2021-05-19 11:54:47 +02:00
parent 04ff11b51c
commit 377d5dc665
3 changed files with 68 additions and 64 deletions
Download Patch File Download Diff File Expand all files Collapse all files

47
src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.C
View File

@ -64,10 +64,11 @@ terminalVelocity::terminalVelocity
// turbKineticEnergyFieldName_(propsDict_.lookupOrDefault<word>("turbKineticEnergyFieldName","")),
// turbDissipationRateFieldName_(propsDict_.lookupOrDefault<word>("turbDissipationRateFieldName","")),
// turbKineticEnergy_(NULL),
// turbDissipationRate_(NULL),
turbDissipationRate_(NULL),
// recurrenceBaseName_(propsDict_.lookupOrDefault<word>("recurrenceBaseName","")),
// recurrenceBase_(NULL),
// existturbDissipationRateInObjReg_(false),
existturbDissipationRateInObjReg_(false),
turbulenceCorrection_(propsDict_.lookupOrDefault<bool>("turbulenceCorrection",false)),
// delta( IOobject
// (
// "delta",
@ -90,8 +91,8 @@ terminalVelocity::terminalVelocity
sm.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
),
// liquidViscosity_(propsDict_.lookupOrDefault<scalar>("liquidViscosity", 1.0)),
// dragReductionFactor_(propsDict_.lookupOrDefault<scalar>("dragReductionFactor", 0.0)),
liquidViscosity_(propsDict_.lookupOrDefault<scalar>("liquidViscosity", 1.0)),
dragReductionFactor_(propsDict_.lookupOrDefault<scalar>("dragReductionFactor", 0.0)),
gravityFieldName_(propsDict_.lookupOrDefault<word>("gravityFieldName","g")),
g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
{
@ -122,7 +123,7 @@ terminalVelocity::terminalVelocity
recurrenceBase_ = &rBase_;
existturbDissipationRateInObjReg_ = true;
}
*/
turbDissipationRate_ = new volScalarField
(
IOobject
@ -136,7 +137,7 @@ terminalVelocity::terminalVelocity
mesh_,
dimensionedScalar("zero", dimensionSet(0,2,-3,0,0), 0)
);
/*
if (turbKineticEnergyFieldName_ != "")
{
turbKineticEnergy_ = new volScalarField
@ -176,7 +177,7 @@ terminalVelocity::terminalVelocity
terminalVelocity::~terminalVelocity()
{
// if (!existturbDissipationRateInObjReg_) delete turbDissipationRate_;
if (!existturbDissipationRateInObjReg_) delete turbDissipationRate_;
}
// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * //
@ -191,13 +192,15 @@ bool terminalVelocity::ignoreCell(label cell) const
void terminalVelocity::setForce() const
{
// updateEpsilon();
updateEpsilon();
vector position(0,0,0);
// scalar Urising=0.0;
label cellI=-1;
scalar radius=0.0;
// scalar epsilon=0.0;
scalar epsilon=0.0;
scalar dLambda = 0.0;
scalar velReductionFactor = 0.0;
vector Uparticle(0,0,0);
label patchID = -1;
@ -206,7 +209,7 @@ void terminalVelocity::setForce() const
vector faceINormal = vector::zero;
word patchName("");
// interpolationCellPoint<scalar> turbDissipationRateInterpolator_(*turbDissipationRate_);
interpolationCellPoint<scalar> turbDissipationRateInterpolator_(*turbDissipationRate_);
for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
{
cellI = particleCloud_.cellIDs()[index][0];
@ -217,20 +220,23 @@ void terminalVelocity::setForce() const
if (interpolate_)
{
position = particleCloud_.position(index);
// epsilon = turbDissipationRateInterpolator_.interpolate(position,cellI);
epsilon = turbDissipationRateInterpolator_.interpolate(position,cellI);
}
else
{
// epsilon = (*turbDissipationRate_)[cellI];
epsilon = (*turbDissipationRate_)[cellI];
}
position = particleCloud_.position(index);
radius = particleCloud_.radius(index);
if (turbulenceCorrection_)
{
// d * kolmogorov length scale
// scalar dLambda = 2*radius*pow(epsilon,0.25)/pow(liquidViscosity_,0.75);
// Urising = terminalVel_ / Foam::sqrt(1+(dragReductionFactor_*pow(dLambda,3)) );
dLambda = 2*radius*pow(epsilon,0.25)/pow(liquidViscosity_,0.75);
velReductionFactor = Foam::sqrt( 1 + ( dragReductionFactor_*pow(dLambda,3)));
terminalVel_ = terminalVel_ / velReductionFactor;
}
// particleCloud_.particleConvVels()[index][] += Urising_;
// read the new particle velocity
@ -267,21 +273,19 @@ void terminalVelocity::setForce() const
}
}
}
}
if (forceSubM(0).verbose() && index >0 && index <2)
{
Pout << "cellI = " << cellI << endl;
Pout << "index = " << index << endl;
// Pout << "epsilon = " << epsilon << endl;
Pout << "epsilon = " << epsilon << endl;
Pout << "rising velocity = " << terminalVel_ << endl;
}
}
}
/*
void terminalVelocity::updateEpsilon() const
{
if (!existturbDissipationRateInObjReg_)
@ -289,7 +293,7 @@ void terminalVelocity::updateEpsilon() const
Info << "epsilon is calculated from the turbulence model. " << endl;
*turbDissipationRate_ = particleCloud_.turbulence().epsilon()();
}
/*
else if (recurrenceBaseName_ != "")
{
Info << "updating epsilon from the database. " << endl;
@ -311,8 +315,9 @@ void terminalVelocity::updateEpsilon() const
{
FatalError <<"turbulent dissipation enery must be calculated either from the turbulence model or recurrence dataBase!\n" << abort(FatalError);
}
}
*/
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

14
src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.H
View File

@ -66,21 +66,23 @@ private:
// volScalarField* turbKineticEnergy_;
// volScalarField* turbDissipationRate_;
volScalarField* turbDissipationRate_;
// mutable volScalarField delta;
mutable volScalarField wallIndicatorField_;
// bool existturbDissipationRateInObjReg_;
bool existturbDissipationRateInObjReg_;
// scalar liquidViscosity_;
bool turbulenceCorrection_;
scalar liquidViscosity_;
// drag reduction factor in turbulent flow
// scalar dragReductionFactor_;
scalar dragReductionFactor_;
// terminal rising velocity of the bubble in the resting liquid
vector terminalVel_;
mutable vector terminalVel_;
word gravityFieldName_;
@ -88,7 +90,7 @@ private:
bool ignoreCell(label) const;
// void updateEpsilon() const;
void updateEpsilon() const;
// recBase* recurrenceBase_;

71
src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C
View File

@ -198,47 +198,44 @@ void turbulentDispersion::setForce() const
cellI = particleCloud_.cellIDs()[index][0];
if (cellI > -1 && !ignoreCell(cellI))
{
// particles in dilute regions follow fluid without fluctuations
if (voidfraction_[cellI] < critVoidfraction_)
if (interpolate_)
{
if (interpolate_)
position = particleCloud_.position(index);
D = nutInterpolator_.interpolate(position,cellI) / turbulentSchmidtNumber_;
}
else
{
D = nut_[cellI] / turbulentSchmidtNumber_;
}
// include concentration dependence on the diffusivity at this point if necessary
flucU=unitFlucDir()*Foam::sqrt(6.0*D/dt_);
// prevent particles being pushed through walls by regulating velocity fluctuations
// check if cell is adjacent to wall and remove corresponding components
if (wallIndicatorField_[cellI] > 0.5)
{
const cell& faces = mesh_.cells()[cellI];
forAll (faces, faceI)
{
position = particleCloud_.position(index);
D = nutInterpolator_.interpolate(position,cellI) / turbulentSchmidtNumber_;
}
else
{
D = nut_[cellI] / turbulentSchmidtNumber_;
faceIGlobal = faces[faceI];
patchID = mesh_.boundaryMesh().whichPatch(faceIGlobal);
if (patchID < 0) continue;
patchName = mesh_.boundary()[patchID].name();
if (patchName.rfind("procB",0) == 0) continue;
faceINormal = mesh_.Sf()[faceIGlobal];
faceINormal /= mag(faceINormal);
flucProjection = faceINormal&flucU;
if (flucProjection > 0.0) flucU -= flucProjection*faceINormal;
}
}
flucU=unitFlucDir()*Foam::sqrt(6.0*D/dt_);
// prevent particles being pushed through walls by regulating velocity fluctuations
// check if cell is adjacent to wall and remove corresponding components
if (wallIndicatorField_[cellI] > 0.5)
{
const cell& faces = mesh_.cells()[cellI];
forAll (faces, faceI)
{
faceIGlobal = faces[faceI];
patchID = mesh_.boundaryMesh().whichPatch(faceIGlobal);
if (patchID < 0) continue;
patchName = mesh_.boundary()[patchID].name();
if (patchName.rfind("procB",0) == 0) continue;
faceINormal = mesh_.Sf()[faceIGlobal];
faceINormal /= mag(faceINormal);
flucProjection = faceINormal&flucU;
if (flucProjection > 0.0) flucU -= flucProjection*faceINormal;
}
}
for(int j=0;j<3;j++)
{
particleCloud_.particleFlucVels()[index][j] += flucU[j];
}
for(int j=0;j<3;j++)
{
particleCloud_.particleFlucVels()[index][j] += flucU[j];
}
}
}