diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C index 53fadacc..f62e1be7 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C @@ -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 * ( diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesDynFanning.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C similarity index 80% rename from src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesDynFanning.C rename to src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C index 9175e30c..8d4fb4dd 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesDynFanning.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C @@ -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 (UsFieldName_)), UDyn_(sm.mesh().lookupObject ("uDyn")), - alphaDyn_(sm.mesh().lookupObject ("alphaDyn")), + FanningCoeff_(sm.mesh().lookupObject ("FanningCoeff")), alphaP_(sm.mesh().lookupObject ("alphaP")), - dHydMix_(sm.mesh().lookupObject ("dHydMix")), dSauter_(sm.mesh().lookupObject ("dSauter")), - Froude_(sm.mesh().lookupObject ("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); } } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesDynFanning.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.H similarity index 82% rename from src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesDynFanning.H rename to src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.H index fef9faa8..ff1803bf 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesDynFanning.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.H @@ -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 diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C index 1be78714..ec833ec8 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C @@ -40,60 +40,29 @@ FinesFields::FinesFields : particleCloud_(sm), propsDict_(dict.subDict(typeName + "Props")), - velFieldName_(propsDict_.lookup("velFieldName")), + velFieldName_(propsDict_.lookupOrDefault("velFieldName","U")), U_(sm.mesh().lookupObject (velFieldName_)), - voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")), + voidfractionFieldName_(propsDict_.lookupOrDefault("voidfractionFieldName","voidfraction")), // this is probably really bad voidfraction_(const_cast(sm.mesh().lookupObject (voidfractionFieldName_))), - UsFieldName_(propsDict_.lookup("granVelFieldName")), + UsFieldName_(propsDict_.lookupOrDefault("granVelFieldName","Us")), UsField_(sm.mesh().lookupObject (UsFieldName_)), - dSauter_(sm.mesh().lookupObject ("dSauter")), - dSauterMix_ + pFieldName_(propsDict_.lookupOrDefault("pFieldName","p")), + p_(sm.mesh().lookupObject (pFieldName_)), + rhoGFieldName_(propsDict_.lookupOrDefault("rhoGFieldName","rho")), + rhoG_(sm.mesh().lookupObject (rhoGFieldName_)), + dSauter_(sm.mesh().lookupObject ("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 diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H index 59b9f161..b2c76968 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H @@ -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 diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/gradPStatFines.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/gradPStatFines.C new file mode 100644 index 00000000..05d8800a --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/gradPStatFines.C @@ -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 . + + 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 (pFieldName_)), + velocityFieldName_(propsDict_.lookup("velocityFieldName")), + U_(sm.mesh().lookupObject (velocityFieldName_)), + alphaP_(sm.mesh().lookupObject ("alphaP")), + alphaSt_(sm.mesh().lookupObject ("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 ("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 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 + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/gradPStatFines.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/gradPStatFines.H new file mode 100644 index 00000000..25948776 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/gradPStatFines.H @@ -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 . + + 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 + +// ************************************************************************* //