diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.C b/src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.C index cdfaa55f..98064a76 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.C @@ -64,10 +64,11 @@ terminalVelocity::terminalVelocity // turbKineticEnergyFieldName_(propsDict_.lookupOrDefault("turbKineticEnergyFieldName","")), // turbDissipationRateFieldName_(propsDict_.lookupOrDefault("turbDissipationRateFieldName","")), // turbKineticEnergy_(NULL), -// turbDissipationRate_(NULL), + turbDissipationRate_(NULL), // recurrenceBaseName_(propsDict_.lookupOrDefault("recurrenceBaseName","")), // recurrenceBase_(NULL), -// existturbDissipationRateInObjReg_(false), + existturbDissipationRateInObjReg_(false), + turbulenceCorrection_(propsDict_.lookupOrDefault("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("liquidViscosity", 1.0)), -// dragReductionFactor_(propsDict_.lookupOrDefault("dragReductionFactor", 0.0)), + liquidViscosity_(propsDict_.lookupOrDefault("liquidViscosity", 1.0)), + dragReductionFactor_(propsDict_.lookupOrDefault("dragReductionFactor", 0.0)), gravityFieldName_(propsDict_.lookupOrDefault("gravityFieldName","g")), g_(sm.mesh().lookupObject (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 turbDissipationRateInterpolator_(*turbDissipationRate_); + interpolationCellPoint 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 diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.H b/src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.H index b5fe43b9..8b5d1338 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/terminalVelocity/terminalVelocity.H @@ -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_; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C index a4fb742e..a9d05327 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C @@ -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]; } } }