diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C index 73ec3fca..8f1a8534 100644 --- a/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C +++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C @@ -81,7 +81,7 @@ int main(int argc, char *argv[]) { particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces()); } - + Info << "update Ksl.internalField()" << endl; Ksl = particleCloud.momCoupleM(0).impMomSource(); Ksl.correctBoundaryConditions(); diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H index 606d9fe2..0f1e6aff 100644 --- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H +++ b/applications/solvers/cfdemSolverRhoPimple/EEqn.H @@ -39,7 +39,7 @@ == fvOptions(rho, he) ); - + EEqn.relax(); @@ -53,7 +53,8 @@ Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl; - particleCloud.clockM().start(31,"postFlow"); - particleCloud.postFlow(); - particleCloud.clockM().stop("postFlow"); + particleCloud.clockM().start(31,"energySolve"); + particleCloud.solve(); + particleCloud.clockM().stop("energySolve"); + } diff --git a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C index b9cf2fcb..39fa85cc 100644 --- a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C +++ b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C @@ -139,6 +139,10 @@ int main(int argc, char *argv[]) } } + particleCloud.clockM().start(31,"postFlow"); + particleCloud.postFlow(); + particleCloud.clockM().stop("postFlow"); + runTime.write(); diff --git a/applications/solvers/cfdemSolverRhoSimple/EEqn.H b/applications/solvers/cfdemSolverRhoSimple/EEqn.H new file mode 100644 index 00000000..3dbb5c3f --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/EEqn.H @@ -0,0 +1,57 @@ +// contributions to internal energy equation can be found in +// Crowe et al.: "Multiphase flows with droplets and particles", CRC Press 1998 +{ + // dim he = J / kg + volScalarField& he = thermo.he(); + particleCloud.energyContributions(Qsource); + particleCloud.energyCoefficients(QCoeff); + + //thDiff=particleCloud.thermCondM().thermDiff(); + thCond=particleCloud.thermCondM().thermCond(); + + addSource = + ( + he.name() == "e" + ? + fvc::div(phi, K) + + fvc::div + ( + fvc::absolute(phi/fvc::interpolate(rho), voidfraction*U), + p, + "div(phiv,p)" + ) + : fvc::div(phi, K) + ); + + Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp(); + + + fvScalarMatrix EEqn + ( + fvm::div(phi, he) + + addSource + - Qsource + - fvm::Sp(QCoeff/Cpv, he) + - fvm::laplacian(voidfraction*thCond/Cpv,he) + == + fvOptions(rho, he) + ); + + + EEqn.relax(); + + fvOptions.constrain(EEqn); + + EEqn.solve(); + + fvOptions.correct(he); + + thermo.correct(); + + Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl; + + + particleCloud.clockM().start(31,"energySolve"); + particleCloud.solve(); + particleCloud.clockM().stop("energySolve"); +} diff --git a/applications/solvers/cfdemSolverRhoSimple/Make/files b/applications/solvers/cfdemSolverRhoSimple/Make/files new file mode 100644 index 00000000..d9db0744 --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/Make/files @@ -0,0 +1,3 @@ +cfdemSolverRhoSimple.C + +EXE=$(CFDEM_APP_DIR)/cfdemSolverRhoSimple diff --git a/applications/solvers/cfdemSolverRhoSimple/Make/options b/applications/solvers/cfdemSolverRhoSimple/Make/options new file mode 100644 index 00000000..0377ece5 --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/Make/options @@ -0,0 +1,32 @@ +include $(CFDEM_ADD_LIBS_DIR)/additionalLibs + +PFLAGS+= -Dcompre + +EXE_INC = \ + $(PFLAGS) \ + -I$(CFDEM_OFVERSION_DIR) \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ + -I$(LIB_SRC)/finiteVolume/cfdTools \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \ + -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \ + +EXE_LIBS = \ + -L$(CFDEM_LIB_DIR)\ + -lcompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lspecie \ + -lturbulenceModels \ + -lcompressibleTurbulenceModels \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling \ + -lfvOptions \ + -l$(CFDEM_LIB_COMP_NAME) \ + $(CFDEM_ADD_LIB_PATHS) \ + $(CFDEM_ADD_LIBS) diff --git a/applications/solvers/cfdemSolverRhoSimple/UEqn.H b/applications/solvers/cfdemSolverRhoSimple/UEqn.H new file mode 100644 index 00000000..8fbd30f2 --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/UEqn.H @@ -0,0 +1,30 @@ +// Solve the Momentum equation +particleCloud.otherForces(fOther); + +tmp tUEqn +( + fvm::div(phi, U) + + particleCloud.divVoidfractionTau(U, voidfraction) + + fvm::Sp(Ksl,U) + - fOther + == + fvOptions(rho, U) +); +fvVectorMatrix& UEqn = tUEqn.ref(); + +UEqn.relax(); + +fvOptions.constrain(UEqn); + +if (modelType=="B" || modelType=="Bfull") +{ + solve(UEqn == -fvc::grad(p)+ Ksl*Us); +} +else +{ + solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us); +} + +fvOptions.correct(U); + +K = 0.5*magSqr(U); diff --git a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C new file mode 100644 index 00000000..4e1821fc --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C @@ -0,0 +1,142 @@ +/*---------------------------------------------------------------------------*\ +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 + +Application + cfdemSolverRhoPimple + +Description + Transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE) + algorithm. + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x, + where additional functionality for CFD-DEM coupling is added. +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "psiThermo.H" +#include "turbulentFluidThermoModel.H" +#include "bound.H" +#include "simpleControl.H" +#include "fvOptions.H" +#include "localEulerDdtScheme.H" +#include "fvcSmooth.H" + +#include "cfdemCloudEnergy.H" +#include "implicitCouple.H" +#include "clockModel.H" +#include "smoothingModel.H" +#include "forceModel.H" +#include "thermCondModel.H" +#include "energyModel.H" + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createTimeControls.H" + #include "createRDeltaT.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "createFieldRefs.H" + #include "createFvOptions.H" + + // create cfdemCloud + #include "readGravitationalAcceleration.H" + cfdemCloudEnergy particleCloud(mesh); + #include "checkModelType.H" + + turbulence->validate(); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (simple.loop()) + { + particleCloud.clockM().start(1,"Global"); + + Info<< "Time = " << runTime.timeName() << nl << endl; + + // do particle stuff + particleCloud.clockM().start(2,"Coupling"); + bool hasEvolved = particleCloud.evolve(voidfraction,Us,U); + + if(hasEvolved) + { + particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces()); + } + + Info << "update Ksl.internalField()" << endl; + Ksl = particleCloud.momCoupleM(0).impMomSource(); + Ksl.correctBoundaryConditions(); + + //Force Checks + vector fTotal(0,0,0); + vector fImpTotal = sum(mesh.V()*Ksl.primitiveFieldRef()*(Us.primitiveFieldRef()-U.primitiveFieldRef())); + reduce(fImpTotal, sumOp()); + Info << "TotalForceExp: " << fTotal << endl; + Info << "TotalForceImp: " << fImpTotal << endl; + + #include "solverDebugInfo.H" + particleCloud.clockM().stop("Coupling"); + + particleCloud.clockM().start(26,"Flow"); + + volScalarField rhoeps("rhoeps",rho*voidfraction); + // Pressure-velocity SIMPLE corrector + + #include "UEqn.H" + + + // besides this pEqn, OF offers a "simple consistent"-option + #include "pEqn.H" + rhoeps=rho*voidfraction; + + #include "EEqn.H" + + turbulence->correct(); + + particleCloud.clockM().start(32,"postFlow"); + if(hasEvolved) particleCloud.postFlow(); + particleCloud.clockM().stop("postFlow"); + + runTime.write(); + + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + + particleCloud.clockM().stop("Flow"); + particleCloud.clockM().stop("Global"); + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H b/applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H new file mode 100644 index 00000000..5842906a --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H @@ -0,0 +1,2 @@ +const volScalarField& T = thermo.T(); +const volScalarField& psi = thermo.psi(); diff --git a/applications/solvers/cfdemSolverRhoSimple/createFields.H b/applications/solvers/cfdemSolverRhoSimple/createFields.H new file mode 100644 index 00000000..a9225229 --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/createFields.H @@ -0,0 +1,242 @@ +Info<< "Reading thermophysical properties\n" << endl; + + autoPtr pThermo + ( + psiThermo::New(mesh) + ); + psiThermo& thermo = pThermo(); + thermo.validate(args.executable(), "h", "e"); + volScalarField& p = thermo.p(); + + Info<< "Reading field rho\n" << endl; + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + thermo.rho() + ); + + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl; + volScalarField voidfraction + ( + IOobject + ( + "voidfraction", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + volScalarField addSource + ( + IOobject + ( + "addSource", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + Info<< "\nCreating fluid-particle heat flux field\n" << endl; + volScalarField Qsource + ( + IOobject + ( + "Qsource", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0) + ); + + Info<< "\nCreating fluid-particle heat flux coefficient field\n" << endl; + volScalarField QCoeff + ( + IOobject + ( + "QCoeff", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0) + ); + + Info<< "\nCreating thermal conductivity field\n" << endl; + volScalarField thCond + ( + IOobject + ( + "thCond", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0) + ); + + Info<< "\nCreating heat capacity field\n" << endl; + volScalarField Cpv + ( + IOobject + ( + "Cpv", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0) + ); + + Info<< "\nCreating body force field\n" << endl; + volVectorField fOther + ( + IOobject + ( + "fOther", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero) + ); + + Info<< "Reading/calculating face flux field phi\n" << endl; + surfaceScalarField phi + ( + IOobject + ( + "phi", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + linearInterpolate(rho*U*voidfraction) & mesh.Sf() + ); + + dimensionedScalar rhoMax + ( + dimensionedScalar::lookupOrDefault + ( + "rhoMax", + simple.dict(), + dimDensity, + GREAT + ) + ); + + dimensionedScalar rhoMin + ( + dimensionedScalar::lookupOrDefault + ( + "rhoMin", + simple.dict(), + dimDensity, + 0 + ) + ); + + Info<< "Creating turbulence model\n" << endl; + autoPtr turbulence + ( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) + ); + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, simple.dict(), pRefCell, pRefValue); + + mesh.setFluxRequired(p.name()); + + Info<< "Creating field dpdt\n" << endl; + volScalarField dpdt + ( + IOobject + ( + "dpdt", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("dpdt", p.dimensions()/dimTime, 0) + ); + + Info<< "Creating field kinetic energy K\n" << endl; + volScalarField K("K", 0.5*magSqr(U)); + + Info<< "\nReading momentum exchange field Ksl\n" << endl; + volScalarField Ksl + ( + IOobject + ( + "Ksl", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0) + ); + + + Info<< "Reading particle velocity field Us\n" << endl; + volVectorField Us + ( + IOobject + ( + "Us", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + +//=============================== diff --git a/applications/solvers/cfdemSolverRhoSimple/pEqn.H b/applications/solvers/cfdemSolverRhoSimple/pEqn.H new file mode 100644 index 00000000..5b2dbcbd --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/pEqn.H @@ -0,0 +1,81 @@ +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); + +volScalarField rAU(1.0/UEqn.A()); +surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU)); +if (modelType=="A") +{ + rhorAUf *= fvc::interpolate(voidfraction); +} +volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); + +surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf()); + + +if (simple.transonic()) +{ +// transonic version not implemented yet +} +else +{ + surfaceScalarField phiHbyA + ( + "phiHbyA", + ( + fvc::flux(rhoeps*HbyA) + ) + ); + + // flux without pressure gradient contribution + phi = phiHbyA + phiUs; + + // Update the pressure BCs to ensure flux consistency + constrainPressure(p, rhoeps, U, phi, rhorAUf); + + while (simple.correctNonOrthogonal()) + { + // Pressure corrector + fvScalarMatrix pEqn + ( + fvc::div(phi) + - fvm::laplacian(rhorAUf, p) + == + fvOptions(psi, p, rho.name()) + ); + + pEqn.setReference(pRefCell, pRefValue); + + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) + { + phi += pEqn.flux(); + } + } +} + +// Explicitly relax pressure for momentum corrector +p.relax(); + +// Recalculate density from the relaxed pressure +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); +Info<< "rho max/min : " << max(rho).value() + << " " << min(rho).value() << endl; + +if (modelType=="A") +{ + U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us); +} +else +{ + U = HbyA - rAU*(fvc::grad(p)-Ksl*Us); +} + +U.correctBoundaryConditions(); +fvOptions.correct(U); +K = 0.5*magSqr(U); \ No newline at end of file diff --git a/doc/CFDEMcoupling_Manual.html b/doc/CFDEMcoupling_Manual.html index c100495e..eb31ac01 100644 --- a/doc/CFDEMcoupling_Manual.html +++ b/doc/CFDEMcoupling_Manual.html @@ -215,10 +215,11 @@ listing below of styles within certain commands. dataExchangeModel_twoWayMPIdataExchangeModel_twoWayMany2Many forceModelforceModel_Archimedes forceModel_ArchimedesIBforceModel_DiFeliceDrag -forceModel_GidaspowDragforceModel_KochHillDrag -forceModel_LaEuScalarTempforceModel_MeiLift -forceModel_SchillerNaumannDragforceModel_ShirgaonkarIB -forceModel_fieldStoreforceModel_gradPForce +forceModel_dSauterforceModel_GidaspowDrag +forceModel_KochHillDragforceModel_LaEuScalarTemp +forceModel_MeiLiftforceModel_SchillerNaumannDrag +forceModel_ShirgaonkarIBforceModel_fieldStore +forceModel_pdCorrelationforceModel_gradPForce forceModel_noDragforceModel_particleCellVolume forceModel_virtualMassForceforceModel_viscForce forceSubModelforceSubModel_ImEx diff --git a/doc/CFDEMcoupling_Manual.txt b/doc/CFDEMcoupling_Manual.txt index 3493934c..5a004c21 100644 --- a/doc/CFDEMcoupling_Manual.txt +++ b/doc/CFDEMcoupling_Manual.txt @@ -260,6 +260,7 @@ listing below of styles within certain commands. "forceModel_Archimedes"_forceModel_Archimedes.html, "forceModel_ArchimedesIB"_forceModel_ArchimedesIB.html, "forceModel_DiFeliceDrag"_forceModel_DiFeliceDrag.html, +"forceModel_dSauter"_forceModel_dSauter.html, "forceModel_GidaspowDrag"_forceModel_GidaspowDrag.html, "forceModel_KochHillDrag"_forceModel_KochHillDrag.html, "forceModel_LaEuScalarTemp"_forceModel_LaEuScalarTemp.html, @@ -267,6 +268,7 @@ listing below of styles within certain commands. "forceModel_SchillerNaumannDrag"_forceModel_SchillerNaumannDrag.html, "forceModel_ShirgaonkarIB"_forceModel_ShirgaonkarIB.html, "forceModel_fieldStore"_forceModel_fieldStore.html, +"forceModel_pdCorrelation"_forceModel_pdCorrelation.html, "forceModel_gradPForce"_forceModel_gradPForce.html, "forceModel_noDrag"_forceModel_noDrag.html, "forceModel_particleCellVolume"_forceModel_particleCellVolume.html, diff --git a/doc/Eqs/d32.png b/doc/Eqs/d32.png new file mode 100644 index 00000000..c8f91641 Binary files /dev/null and b/doc/Eqs/d32.png differ diff --git a/doc/Eqs/pdCorrelation.png b/doc/Eqs/pdCorrelation.png new file mode 100644 index 00000000..161dc5e8 Binary files /dev/null and b/doc/Eqs/pdCorrelation.png differ diff --git a/doc/forceModel.html b/doc/forceModel.html index 0e71fa58..fa2d2f78 100644 --- a/doc/forceModel.html +++ b/doc/forceModel.html @@ -41,7 +41,7 @@

Related commands:

-

Archimedes, DiFeliceDrag, gradPForce, viscForce +

Archimedes, DiFeliceDrag, gradPForce, viscForce">dSauter

Note: This examples list may be incomplete - please look for other models (forceModel_XY) in this documentation.

diff --git a/doc/forceModel.txt b/doc/forceModel.txt index 4ae3a1a1..2ada849a 100644 --- a/doc/forceModel.txt +++ b/doc/forceModel.txt @@ -39,7 +39,7 @@ None. [Related commands:] -"Archimedes"_forceModel_Archimedes.html, "DiFeliceDrag"_forceModel_DiFeliceDrag.html, "gradPForce"_forceModel_gradPForce.html, "viscForce"_forceModel_viscForce.html +"Archimedes"_forceModel_Archimedes.html, "DiFeliceDrag"_forceModel_DiFeliceDrag.html, "gradPForce"_forceModel_gradPForce.html, "viscForce"_forceModel_viscForce.html,"dSauter"_forceModel_dSauter.html Note: This examples list may be incomplete - please look for other models (forceModel_XY) in this documentation. diff --git a/doc/forceModel_dSauter.html b/doc/forceModel_dSauter.html new file mode 100644 index 00000000..4f8050e3 --- /dev/null +++ b/doc/forceModel_dSauter.html @@ -0,0 +1,47 @@ + +
CFDEMproject WWW Site - CFDEM Commands +
+ + + + +
+ +

forceModel_dSauter command +

+

Syntax: +

+

Defined in couplingProperties dictionary. +

+
forceModels
+(
+    dSauter
+);
+dSauterProps
+{
+    coarseGrainingFactors
+    (
+        X Y Z
+    );
+}; 
+
+
  • coarseGrainingFactors = list of coarse graining factors by type, separated by whitespace, optional + + +
+

Description: +

+

This "forceModel" does not influence the particles or the flow - it calculates the Sauter diameter +

+
+
+

. +

+

Restrictions: +

+

none. +

+

Related commands: +

+

forceModel

+ diff --git a/doc/forceModel_dSauter.txt b/doc/forceModel_dSauter.txt new file mode 100644 index 00000000..478c8dbd --- /dev/null +++ b/doc/forceModel_dSauter.txt @@ -0,0 +1,44 @@ +"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c + +:link(lws,http://www.cfdem.com) +:link(lc,CFDEMcoupling_Manual.html#comm) + +:line + +forceModel_dSauter command :h3 + +[Syntax:] + +Defined in couplingProperties dictionary. + +forceModels +( + dSauter +); +dSauterProps +\{ + coarseGrainingFactors + ( + X Y Z + ); +\}; :pre + +{coarseGrainingFactors} = list of coarse graining factors by type, separated by whitespace, optional :ulb,l + +:ule + +[Description:] + +This "forceModel" does not influence the particles or the flow - it calculates the Sauter diameter + +:c,image(Eqs/d32.png) +. + + +[Restrictions:] + +none. + +[Related commands:] + +"forceModel"_forceModel.html \ No newline at end of file diff --git a/doc/forceModel_pdCorrelation.html b/doc/forceModel_pdCorrelation.html new file mode 100644 index 00000000..c838e5e2 --- /dev/null +++ b/doc/forceModel_pdCorrelation.html @@ -0,0 +1,58 @@ + +
CFDEMproject WWW Site - CFDEM Commands +
+ + + + +
+ +

forceModel_pdCorrelation +

+

Syntax: +

+

Defined in couplingProperties dictionary. +

+
forceModels
+(
+    pdCorrelation
+);
+pdCorrelationProps
+{
+    coarseGrainingFactors
+    (
+        X Y Z
+    );
+    particleDensities
+    (
+        A B C
+    );
+    runOnWriteOnly true;
+}; 
+
+
  • coarseGrainingFactors = list of coarse graining factors by type, separated by whitespace, optional + +
  • particleDensities = list of particle densities by type, separated by whitespace, optional + +
  • runOnWriteOnly = switch if this should be executed on write, optional (default: false - execute every coupling step). + + +
+

Description: +

+

This "forceModel" does not influence the particles or the flow - it calculates the particle momentum-diameter correlation +

+
+
+

where delta is the type-specific coarsegraining factor. +

+

This model is sensitive to additionally pulled particle type info, and can either use the type-specific densities from the dict or those pulled from LIGGGHTS. +

+

Restrictions: +

+

none. +

+

Related commands: +

+

forceModel

+ diff --git a/doc/forceModel_pdCorrelation.txt b/doc/forceModel_pdCorrelation.txt new file mode 100644 index 00000000..34049c96 --- /dev/null +++ b/doc/forceModel_pdCorrelation.txt @@ -0,0 +1,53 @@ +"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c + +:link(lws,http://www.cfdem.com) +:link(lc,CFDEMcoupling_Manual.html#comm) + +:line + +forceModel_pdCorrelation :h3 + +[Syntax:] + +Defined in couplingProperties dictionary. + +forceModels +( + pdCorrelation +); +pdCorrelationProps +\{ + coarseGrainingFactors + ( + X Y Z + ); + particleDensities + ( + A B C + ); + runOnWriteOnly true; +\}; :pre + +{coarseGrainingFactors} = list of coarse graining factors by type, separated by whitespace, optional :ulb,l +{particleDensities} = list of particle densities by type, separated by whitespace, optional :l +{runOnWriteOnly} = switch if this should be executed on write, optional (default: false - execute every coupling step). :l + + +:ule + +[Description:] + +This "forceModel" does not influence the particles or the flow - it calculates the particle momentum-diameter correlation + +:c,image(Eqs/pdCorrelation.png) +where delta is the type-specific coarsegraining factor. + +This model is sensitive to additionally pulled particle type info, and can either use the type-specific densities from the dict or those pulled from LIGGGHTS. + +[Restrictions:] + +none. + +[Related commands:] + +"forceModel"_forceModel.html \ No newline at end of file diff --git a/etc/bashrc b/etc/bashrc index 00f84d39..159fdf83 100755 --- a/etc/bashrc +++ b/etc/bashrc @@ -17,7 +17,7 @@ #------------------------------------------------------------------------------ export CFDEM_PROJECT=CFDEM -export CFDEM_VERSION=17.08 +export CFDEM_VERSION=18.03 ################################################################################ # USER EDITABLE PART: Changes made here may be lost with the next upgrade diff --git a/etc/cshrc b/etc/cshrc index 3ae372da..e2ac96a0 100755 --- a/etc/cshrc +++ b/etc/cshrc @@ -15,7 +15,7 @@ #------------------------------------------------------------------------------ setenv CFDEM_PROJECT CFDEM -setenv CFDEM_VERSION 17.08 +setenv CFDEM_VERSION 18.03 ################################################################################ # USER EDITABLE PART: Changes made here may be lost with the next upgrade diff --git a/etc/solver-list.txt b/etc/solver-list.txt index 86e94167..f849a120 100644 --- a/etc/solver-list.txt +++ b/etc/solver-list.txt @@ -6,5 +6,7 @@ rtfmSolverSpecies/dir cfdemSolverRhoPimple/dir cfdemSolverPisoMS/dir cfdemSolverPiso/dir +cfdemSolverRhoPimple/dir +cfdemSolverRhoSimple/dir cfdemSolverIB/dir cfdemSolverPisoScalar/dir diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files index 631291a1..58721d89 100644 --- a/src/lagrangian/cfdemParticle/Make/files +++ b/src/lagrangian/cfdemParticle/Make/files @@ -33,12 +33,13 @@ $(cfdTools)/newGlobal.C $(energyModels)/energyModel/energyModel.C $(energyModels)/energyModel/newEnergyModel.C $(energyModels)/heatTransferGunn/heatTransferGunn.C -$(energyModels)/heatTransferGunnImplicit/heatTransferGunnImplicit.C +$(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C $(energyModels)/reactionHeat/reactionHeat.C $(thermCondModels)/thermCondModel/thermCondModel.C $(thermCondModels)/thermCondModel/newThermCondModel.C $(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C +$(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C $(thermCondModels)/noTherm/noThermCond.C $(forceModels)/forceModel/forceModel.C @@ -70,12 +71,14 @@ $(forceModels)/evaluateFluctuations/evaluateFluctuations.C $(forceModels)/directedDiffusiveRelaxation/directedDiffusiveRelaxation.C $(forceModels)/potentialRelaxation/potentialRelaxation.C $(forceModels)/BeetstraDrag/BeetstraDrag.C +$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C $(forceModels)/dSauter/dSauter.C $(forceModels)/Fines/Fines.C $(forceModels)/Fines/FinesFields.C $(forceModels)/Fines/FanningDynFines.C $(forceModels)/Fines/ErgunStatFines.C $(forceModels)/granKineticEnergy/granKineticEnergy.C +$(forceModels)/pdCorrelation/pdCorrelation.C $(forceModelsMS)/forceModelMS/forceModelMS.C $(forceModelsMS)/forceModelMS/newForceModelMS.C diff --git a/src/lagrangian/cfdemParticle/cfdTools/debugInfo.H b/src/lagrangian/cfdemParticle/cfdTools/debugInfo.H index 285e291f..f5968190 100755 --- a/src/lagrangian/cfdemParticle/cfdTools/debugInfo.H +++ b/src/lagrangian/cfdemParticle/cfdTools/debugInfo.H @@ -69,8 +69,8 @@ else { meanAlpha_field = 0; - meanU_field = vector(0,0,0); - meanUs_field = vector(0,0,0); + meanU_field = vector::zero; + meanUs_field = vector::zero; } meanUs_array /= numberOfParticles()+SMALL; meanR_array /= numberOfParticles()+SMALL; diff --git a/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H b/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H index bebaac8c..43650545 100755 --- a/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H +++ b/src/lagrangian/cfdemParticle/cfdTools/versionInfo.H @@ -34,8 +34,8 @@ Description #ifndef versionInfo_H #define versionInfo_H -word CFDEMversion="PFM 17.08"; -word compatibleLIGGGHTSversion="PFM 17.08"; +word CFDEMversion="PFM 18.03"; +word compatibleLIGGGHTSversion="PFM 18.03"; word OFversion="4.x"; Info << "\nCFDEMcoupling version: " << CFDEMversion << endl; diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C index 6828fd51..9147a02f 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C @@ -83,6 +83,9 @@ cfdemCloud::cfdemCloud ignore_(false), allowCFDsubTimestep_(true), limitDEMForces_(false), + getParticleDensities_(couplingProperties_.lookupOrDefault("getParticleDensities",false)), + getParticleEffVolFactors_(couplingProperties_.lookupOrDefault("getParticleEffVolFactors",false)), + getParticleTypes_(couplingProperties_.lookupOrDefault("getParticleTypes",false)), modelType_(couplingProperties_.lookup("modelType")), positions_(NULL), velocities_(NULL), @@ -95,6 +98,9 @@ cfdemCloud::cfdemCloud radii_(NULL), voidfractions_(NULL), cellIDs_(NULL), + particleDensities_(NULL), + particleEffVolFactors_(NULL), + particleTypes_(NULL), particleWeights_(NULL), particleVolumes_(NULL), particleV_(NULL), @@ -129,6 +135,19 @@ cfdemCloud::cfdemCloud mesh, dimensionedScalar("zero", dimensionSet(0,0,-1,0,0), 0) // 1/s ), + particleDensityField_ + ( + IOobject + ( + "partRho", + mesh.time().timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", dimensionSet(1,-3,0,0,0), 0.0) + ), checkPeriodicCells_(false), turbulence_ ( @@ -362,6 +381,9 @@ cfdemCloud::~cfdemCloud() dataExchangeM().destroy(particleV_,1); dataExchangeM().destroy(particleConvVel_,3); dataExchangeM().destroy(particleFlucVel_,3); + if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1); + if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1); + if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1); } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // @@ -373,6 +395,10 @@ void cfdemCloud::getDEMdata() if(impDEMdragAcc_) dataExchangeM().getData("dragAcc","vector-atom",fAcc_); // array is used twice - might be necessary to clean it first + + if(getParticleDensities_) dataExchangeM().getData("density","scalar-atom",particleDensities_); + if(getParticleEffVolFactors_) dataExchangeM().getData("effvolfactor","scalar-atom",particleEffVolFactors_); + if(getParticleTypes_) dataExchangeM().getData("type","scalar-atom",particleTypes_); } void cfdemCloud::giveDEMdata() @@ -452,6 +478,25 @@ void cfdemCloud::setParticleForceField() ); } +void cfdemCloud::setScalarAverages() +{ + if(!getParticleDensities_) return; + if(verbose_) Info << "- setScalarAverage" << endl; + + particleDensityField_.primitiveFieldRef() = 0.0; + averagingM().resetWeightFields(); + averagingM().setScalarAverage + ( + particleDensityField_, + particleDensities_, + particleWeights_, + averagingM().UsWeightField(), + NULL + ); + + if(verbose_) Info << "setScalarAverage done." << endl; +} + void cfdemCloud::setVectorAverages() { if(verbose_) Info << "- setVectorAverage(Us,velocities_,weights_)" << endl; @@ -604,7 +649,8 @@ bool cfdemCloud::evolve clockM().stop("setvoidFraction"); // set average particles velocity field - clockM().start(20,"setVectorAverage"); + clockM().start(20,"setAverages"); + setScalarAverages(); setVectorAverages(); @@ -618,7 +664,7 @@ bool cfdemCloud::evolve if(!treatVoidCellsAsExplicitForce()) smoothingM().smoothenReferenceField(averagingM().UsNext()); - clockM().stop("setVectorAverage"); + clockM().stop("setAverages"); } //============================================ @@ -711,6 +757,9 @@ bool cfdemCloud::reAllocArrays() dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle()); dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle()); dataExchangeM().allocateArray(particleV_,0.,1); + if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1); + if(getParticleEffVolFactors_) dataExchangeM().allocateArray(particleEffVolFactors_,0.,1); + if(getParticleTypes_) dataExchangeM().allocateArray(particleTypes_,0,1); arraysReallocated_ = true; return true; } diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H index c7ea502a..6ecae3e1 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H @@ -96,6 +96,12 @@ protected: bool limitDEMForces_; + bool getParticleDensities_; + + bool getParticleEffVolFactors_; + + bool getParticleTypes_; + scalar maxDEMForce_; const word modelType_; @@ -122,6 +128,12 @@ protected: mutable int **cellIDs_; + mutable double **particleDensities_; + + mutable double **particleEffVolFactors_; + + mutable int **particleTypes_; + mutable double **particleWeights_; mutable double **particleVolumes_; @@ -166,6 +178,8 @@ protected: mutable volScalarField ddtVoidfraction_; + mutable volScalarField particleDensityField_; + mutable Switch checkPeriodicCells_; const turbulenceModel& turbulence_; @@ -209,6 +223,8 @@ protected: virtual void setParticleForceField(); + virtual void setScalarAverages(); + virtual void setVectorAverages(); public: @@ -331,18 +347,26 @@ public: virtual inline int maxType() {return -1;} virtual inline bool multipleTypesDMax() {return false;} virtual inline bool multipleTypesDMin() {return false;} - virtual inline double ** particleDensity() const {return NULL;} - virtual inline int ** particleTypes() const {return NULL;} - virtual label particleType(label index) const {return -1;} + + inline bool getParticleDensities() const; + virtual inline double ** particleDensities() const; + virtual inline double particleDensity(label index) const; + + inline bool getParticleEffVolFactors() const; + virtual inline double particleEffVolFactor(label index) const; + + inline bool getParticleTypes() const; + virtual inline int ** particleTypes() const; + virtual inline label particleType(label index) const; //access to the particle's rotation and torque data virtual inline double ** DEMTorques() const {return NULL;} virtual inline double ** omegaArray() const {return NULL;} - virtual vector omega(int) const {return vector(0,0,0);} + virtual vector omega(int) const {return vector::zero;} //access to the particles' orientation information virtual inline double ** exArray() const {return NULL;} - virtual vector ex(int) const {return vector(0,0,0);} + virtual vector ex(int) const {return vector::zero;} //Detector if SRF module is enable or not virtual inline bool SRFOn(){return false;} diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H index 2edc3c1e..f1863bd4 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H @@ -226,6 +226,58 @@ inline double cfdemCloud::d32(bool recalc) return d32_; } +inline bool cfdemCloud::getParticleDensities() const +{ + return getParticleDensities_; +} + +inline double ** cfdemCloud::particleDensities() const +{ + return particleDensities_; +} + +inline double cfdemCloud::particleDensity(label index) const +{ + if(!getParticleDensities_) return -1.0; + else + { + return particleDensities_[index][0]; + } +} + +inline bool cfdemCloud::getParticleEffVolFactors() const +{ + return getParticleEffVolFactors_; +} + +inline double cfdemCloud::particleEffVolFactor(label index) const +{ + if(!getParticleEffVolFactors_) return -1.0; + else + { + return particleEffVolFactors_[index][0]; + } +} + +inline bool cfdemCloud::getParticleTypes() const +{ + return getParticleTypes_; +} + +inline int ** cfdemCloud::particleTypes() const +{ + return particleTypes_; +} + +inline label cfdemCloud::particleType(label index) const +{ + if(!getParticleDensities_) return -1; + else + { + return particleTypes_[index][0]; + } +} + inline int cfdemCloud::numberOfParticles() const { return numberOfParticles_; diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C b/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C index 1ed2727b..4a2f44d3 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C @@ -169,6 +169,12 @@ void cfdemCloudEnergy::postFlow() energyModel_[i]().postFlow(); } +void cfdemCloudEnergy::solve() +{ + for (int i=0;i> just_read; // skip text for dataType // need to distinguish between file formats from dump custom/vtk and from LPP if(just_read.compare("3") == 0) - { - input >> just_read; - input >> just_read; - } + { + input >> just_read; + input >> just_read; + } for (int index = 0; index> field[index][0] >> field[index][1] >> field[index][2]; diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic index 4035f2a6..b2c8df9c 100644 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic @@ -27,7 +27,7 @@ lib: $(OBJ) $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) clean: - rm $(LIB) *.o *.d + rm -f $(LIB) *.o *.d # Compilation rules diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic_debug b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic_debug index 28614b91..4749e990 100644 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic_debug +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic_debug @@ -27,7 +27,7 @@ lib: $(OBJ) $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) clean: - rm $(LIB) *.o *.d + rm -f $(LIB) *.o *.d # Compilation rules diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.g++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.g++ index 07e58076..cf4211d9 100644 --- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.g++ +++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.g++ @@ -26,7 +26,7 @@ lib: $(OBJ) $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) clean: - rm $(LIB) *.o *.d + rm -f $(LIB) *.o *.d # Compilation rules diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/energyModel/energyModel.H b/src/lagrangian/cfdemParticle/subModels/energyModel/energyModel/energyModel.H index eef829ed..32f44070 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/energyModel/energyModel.H +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/energyModel/energyModel.H @@ -104,6 +104,8 @@ public: virtual void calcEnergyContribution() = 0; virtual void postFlow() {} + + virtual void solve() {} scalar Cp() const; diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C index 3f7ad3d3..679552d1 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C @@ -46,7 +46,7 @@ heatTransferGunn::heatTransferGunn expNusselt_(propsDict_.lookupOrDefault("expNusselt",false)), interpolation_(propsDict_.lookupOrDefault("interpolation",false)), verbose_(propsDict_.lookupOrDefault("verbose",false)), - NusseltFieldName_(propsDict_.lookupOrDefault("NusseltFieldName","NuField")), + implicit_(propsDict_.lookupOrDefault("implicit",true)), QPartFluidName_(propsDict_.lookupOrDefault("QPartFluidName","QPartFluid")), QPartFluid_ ( IOobject @@ -61,13 +61,26 @@ heatTransferGunn::heatTransferGunn dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0), "zeroGradient" ), + QPartFluidCoeffName_(propsDict_.lookupOrDefault("QPartFluidCoeffName","QPartFluidCoeff")), + QPartFluidCoeff_ + ( IOobject + ( + QPartFluidCoeffName_, + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + sm.mesh(), + dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0) + ), partTempField_ ( IOobject ( - "particleTemp", + "partTemp", sm.mesh().time().timeName(), sm.mesh(), - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, IOobject::NO_WRITE ), sm.mesh(), @@ -100,6 +113,7 @@ heatTransferGunn::heatTransferGunn dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0), "zeroGradient" ), + NusseltFieldName_(propsDict_.lookupOrDefault("NusseltFieldName","NuField")), NuField_ ( IOobject ( @@ -130,10 +144,11 @@ heatTransferGunn::heatTransferGunn partTemp_(NULL), partHeatFluxName_(propsDict_.lookup("partHeatFluxName")), partHeatFlux_(NULL), + partHeatFluxCoeff_(NULL), partRe_(NULL), partNu_(NULL) { - allocateMyArrays(); + allocateMyArrays(); if (propsDict_.found("maxSource")) { @@ -154,7 +169,12 @@ heatTransferGunn::heatTransferGunn partRelTempField_.write(); Info << "Particle temperature field activated." << endl; } - + + if (!implicit_) + { + QPartFluidCoeff_.writeOpt() = IOobject::NO_WRITE; + } + if (verbose_) { ReField_.writeOpt() = IOobject::AUTO_WRITE; @@ -184,6 +204,10 @@ heatTransferGunn::~heatTransferGunn() particleCloud_.dataExchangeM().destroy(partHeatFlux_,1); particleCloud_.dataExchangeM().destroy(partRe_,1); particleCloud_.dataExchangeM().destroy(partNu_,1); + if (implicit_) + { + particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1); + } } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // @@ -193,6 +217,10 @@ void heatTransferGunn::allocateMyArrays() const double initVal=0.0; particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1); // field/initVal/with/lenghtFromLigghts particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1); + if(implicit_) + { + particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1); + } if(verbose_) { @@ -201,7 +229,68 @@ void heatTransferGunn::allocateMyArrays() const } } + +void heatTransferGunn::partTempField() +{ + partTempField_.primitiveFieldRef() = 0.0; + particleCloud_.averagingM().resetWeightFields(); + particleCloud_.averagingM().setScalarAverage + ( + partTempField_, + partTemp_, + particleCloud_.particleWeights(), + particleCloud_.averagingM().UsWeightField(), + NULL + ); + + dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), partTempAve_); + partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_); +} + + +scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) const +{ + return (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) * + (1 + 0.7 * Foam::pow(Rep,0.2) * Foam::pow(Pr,0.33)) + + (1.33 - 2.4 * voidfraction + 1.2 * voidfraction * voidfraction) * + Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33); +} + +void heatTransferGunn::giveData() +{ + Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl; + particleCloud_.clockM().start(30,"giveDEM_Tdata"); + particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_); + particleCloud_.clockM().stop("giveDEM_Tdata"); +} + +void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid) +{ + scalar hAs = h * As; + partHeatFlux_[index][0] = - hAs * partTemp_[index][0]; + if(!implicit_) + { + partHeatFlux_[index][0] += hAs * Tfluid; + } + else + { + partHeatFluxCoeff_[index][0] = hAs; + } +} + // * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * // +void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const +{ + Qsource += QPartFluid_; +} + +void heatTransferGunn::addEnergyCoefficient(volScalarField& Qsource) const +{ + if(implicit_) + { + Qsource += QPartFluidCoeff_; + } +} void heatTransferGunn::calcEnergyContribution() { @@ -292,7 +381,6 @@ void heatTransferGunn::calcEnergyContribution() // calc convective heat flux [W] heatFlux(index, h, As, Tfluid); - heatFluxCoeff(index, h, As); if(verbose_) { @@ -326,6 +414,16 @@ void heatTransferGunn::calcEnergyContribution() if(calcPartTempField_) partTempField(); + // gather particle temperature sums and obtain average + if(calcPartTempAve_) + { + reduce(Tsum, sumOp()); + partTempAve_ = Tsum / particleCloud_.numberOfParticles(); + Info << "mean particle temperature = " << partTempAve_ << endl; + } + + if(calcPartTempField_) partTempField(); + particleCloud_.averagingM().setScalarSum ( QPartFluid_, @@ -336,6 +434,21 @@ void heatTransferGunn::calcEnergyContribution() QPartFluid_.primitiveFieldRef() /= -QPartFluid_.mesh().V(); + if(implicit_) + { + QPartFluidCoeff_.primitiveFieldRef() = 0.0; + + particleCloud_.averagingM().setScalarSum + ( + QPartFluidCoeff_, + partHeatFluxCoeff_, + particleCloud_.particleWeights(), + NULL + ); + + QPartFluidCoeff_.primitiveFieldRef() /= -QPartFluidCoeff_.mesh().V(); + } + if(verbose_) { ReField_.primitiveFieldRef() = 0.0; @@ -360,82 +473,60 @@ void heatTransferGunn::calcEnergyContribution() ); } - // limit source term - forAll(QPartFluid_,cellI) + // limit source term in explicit treatment + if(!implicit_) { - scalar EuFieldInCell = QPartFluid_[cellI]; - - if(mag(EuFieldInCell) > maxSource_ ) + forAll(QPartFluid_,cellI) { - Pout << "limiting source term\n" << endl ; - QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_; + scalar EuFieldInCell = QPartFluid_[cellI]; + + if(mag(EuFieldInCell) > maxSource_ ) + { + Pout << "limiting source term\n" << endl ; + QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_; + } } } QPartFluid_.correctBoundaryConditions(); - - giveData(0); - } -void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const +void heatTransferGunn::postFlow() { - Qsource += QPartFluid_; -} - -scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) const -{ - return (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) * - (1 + 0.7 * Foam::pow(Rep,0.2) * Foam::pow(Pr,0.33)) + - (1.33 - 2.4 * voidfraction + 1.2 * voidfraction * voidfraction) * - Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33); -} - -void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid) -{ - partHeatFlux_[index][0] = h * As * (Tfluid - partTemp_[index][0]); -} - -void heatTransferGunn::heatFluxCoeff(label index, scalar h, scalar As) -{ - //no heat transfer coefficient in explicit model -} - -void heatTransferGunn::giveData(int call) -{ - if(call == 0) + if(implicit_) { - Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl; + label cellI; + scalar Tfluid(0.0); + scalar Tpart(0.0); + interpolationCellPoint TInterpolator_(tempField_); - particleCloud_.clockM().start(30,"giveDEM_Tdata"); - particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_); - particleCloud_.clockM().stop("giveDEM_Tdata"); + for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + { + cellI = particleCloud_.cellIDs()[index][0]; + if(cellI >= 0) + { + if(interpolation_) + { + vector position = particleCloud_.position(index); + Tfluid = TInterpolator_.interpolate(position,cellI); + } + else + { + Tfluid = tempField_[cellI]; + } + + Tpart = partTemp_[index][0]; + partHeatFlux_[index][0] = (Tfluid - Tpart) * partHeatFluxCoeff_[index][0]; + } + } } + giveData(); } scalar heatTransferGunn::aveTpart() const { return partTempAve_; } - -void heatTransferGunn::partTempField() -{ - partTempField_.primitiveFieldRef() = 0.0; - particleCloud_.averagingM().resetWeightFields(); - particleCloud_.averagingM().setScalarAverage - ( - partTempField_, - partTemp_, - particleCloud_.particleWeights(), - particleCloud_.averagingM().UsWeightField(), - NULL - ); - - //volScalarField sumTp (particleCloud_.averagingM().UsWeightField() * partTempField_); - // dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), gSum(sumTp) / particleCloud_.numberOfParticles()); - dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), partTempAve_); - partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_); -} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H index 63a6505f..fb077759 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H @@ -15,7 +15,7 @@ License along with this code. If not, see . Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria - + Description Correlation for Nusselt number according to Gunn, D. J. International Journal of Heat and Mass Transfer 21.4 (1978) @@ -44,29 +44,35 @@ class heatTransferGunn protected: dictionary propsDict_; - + bool expNusselt_; - + bool interpolation_; - + bool verbose_; - - word NusseltFieldName_; - + + bool implicit_; + word QPartFluidName_; - + volScalarField QPartFluid_; - - volScalarField partTempField_; - + + word QPartFluidCoeffName_; + + volScalarField QPartFluidCoeff_; + + mutable volScalarField partTempField_; + volScalarField partRelTempField_; - + volScalarField ReField_; - + + word NusseltFieldName_; + volScalarField NuField_; - + dimensionedScalar partRefTemp_; - + bool calcPartTempField_; bool calcPartTempAve_; @@ -86,35 +92,35 @@ protected: word velFieldName_; const volVectorField& U_; - + word densityFieldName_; - + const volScalarField& rho_; word partTempName_; - mutable double **partTemp_; // Lagrangian array + mutable double **partTemp_; word partHeatFluxName_; - mutable double **partHeatFlux_; // Lagrangian array - + mutable double **partHeatFlux_; + + mutable double **partHeatFluxCoeff_; + mutable double **partRe_; - + mutable double **partNu_; void allocateMyArrays() const; - + void partTempField(); - + scalar Nusselt(scalar, scalar, scalar) const; - - virtual void giveData(int); - + + virtual void giveData(); + virtual void heatFlux(label, scalar, scalar, scalar); - - virtual void heatFluxCoeff(label, scalar, scalar); - + public: //- Runtime type information @@ -138,11 +144,13 @@ public: // Member Functions void addEnergyContribution(volScalarField&) const; - - void addEnergyCoefficient(volScalarField&) const {} - + + void addEnergyCoefficient(volScalarField&) const; + void calcEnergyContribution(); + void postFlow(); + scalar aveTpart() const; }; diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C new file mode 100644 index 00000000..2d8af0e6 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C @@ -0,0 +1,244 @@ +/*---------------------------------------------------------------------------*\ +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 + +\*---------------------------------------------------------------------------*/ + +#include "error.H" +#include "heatTransferGunnPartField.H" +#include "addToRunTimeSelectionTable.H" +#include "thermCondModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(heatTransferGunnPartField, 0); + +addToRunTimeSelectionTable(energyModel, heatTransferGunnPartField, dictionary); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +heatTransferGunnPartField::heatTransferGunnPartField +( + const dictionary& dict, + cfdemCloudEnergy& sm +) +: + heatTransferGunn(dict,sm), + partCpField_ + ( + IOobject + ( + "partCp", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + sm.mesh(), + dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0), + "zeroGradient" + ), + partRhoField_(sm.mesh().lookupObject("partRho")), + typeCp_(propsDict_.lookupOrDefault("specificHeatCapacities",scalarList(1,-1.0))), + partCp_(NULL), + pTMax_(dimensionedScalar("pTMax",dimensionSet(0,0,0,1,0,0,0), -1.0)), + pTMin_(dimensionedScalar("pTMin",dimensionSet(0,0,0,1,0,0,0), -1.0)), + thermCondModel_ + ( + thermCondModel::New + ( + propsDict_, + sm + ) + ), + fvOptions(fv::options::New(sm.mesh())) +{ + if (!implicit_) + { + FatalError << "heatTransferGunnPartField requires implicit heat transfer treatment." << abort(FatalError); + } + + if (typeCp_[0] < 0.0) + { + FatalError << "heatTransferGunnPartField: provide list of specific heat capacities." << abort(FatalError); + } + + if (propsDict_.found("pTMax")) + { + pTMax_.value()=scalar(readScalar(propsDict_.lookup("pTMax"))); + } + + if (propsDict_.found("pTMin")) + { + pTMin_.value()=scalar(readScalar(propsDict_.lookup("pTMin"))); + } + + partTempField_.writeOpt() = IOobject::AUTO_WRITE; + + allocateMyArrays(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +heatTransferGunnPartField::~heatTransferGunnPartField() +{ + particleCloud_.dataExchangeM().destroy(partCp_,1); +} + +// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // + +void heatTransferGunnPartField::allocateMyArrays() const +{ + double initVal=0.0; + particleCloud_.dataExchangeM().allocateArray(partCp_,initVal,1); +} +// * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * // +void heatTransferGunnPartField::calcEnergyContribution() +{ + allocateMyArrays(); + heatTransferGunn::calcEnergyContribution(); + + // if heat sources in particles present, pull them here + + // loop over all particles to fill partCp_ based on type + label cellI=0; + label partType = 0; + for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + { + cellI = particleCloud_.cellIDs()[index][0]; + if(cellI >= 0) + { + partType = particleCloud_.particleType(index); + // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0 + partCp_[index][0] = typeCp_[partType - 1]; + } + } + + partCpField_.primitiveFieldRef() = 0.0; + particleCloud_.averagingM().resetWeightFields(); + particleCloud_.averagingM().setScalarAverage + ( + partCpField_, + partCp_, + particleCloud_.particleWeights(), + particleCloud_.averagingM().UsWeightField(), + NULL + ); + +} + +void heatTransferGunnPartField::addEnergyContribution(volScalarField& Qsource) const +{ + Qsource -= QPartFluidCoeff_*partTempField_; +} + +void heatTransferGunnPartField::giveData() +{ + particleCloud_.dataExchangeM().giveData(partTempName_,"scalar-atom", partTemp_); +} + +void heatTransferGunnPartField::postFlow() +{ + label cellI; + scalar Tpart(0.0); + interpolationCellPoint partTInterpolator_(partTempField_); + + particleCloud_.dataExchangeM().allocateArray(partTemp_,0.0,1); + + for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + { + cellI = particleCloud_.cellIDs()[index][0]; + if(cellI >= 0) + { + if(interpolation_) + { + vector position = particleCloud_.position(index); + Tpart = partTInterpolator_.interpolate(position,cellI); + } + else + { + Tpart = partTempField_[cellI]; + } + partTemp_[index][0] = Tpart; + } + } + + giveData(); +} + +void heatTransferGunnPartField::solve() +{ + // for some weird reason, the particle-fluid heat transfer fields were defined with a negative sign + + volScalarField alphaP = 1.0 - voidfraction_; + volScalarField correctedQPartFluidCoeff(QPartFluidCoeff_); + // no heattransfer in empty cells -- for numerical stability, add small amount + forAll(correctedQPartFluidCoeff,cellI) + { + if (-correctedQPartFluidCoeff[cellI] < SMALL) correctedQPartFluidCoeff[cellI] = -SMALL; + } + + volScalarField Qsource = correctedQPartFluidCoeff*tempField_; + volScalarField partCpEff = alphaP*partRhoField_*partCpField_; + volScalarField thCondEff = alphaP*thermCondModel_().thermCond(); +// thCondEff.correctBoundaryConditions(); + + + + fvScalarMatrix partTEqn + ( + // fvm::ddt(partCpEff, partTempField_) + // + Qsource + Qsource + - fvm::Sp(correctedQPartFluidCoeff, partTempField_) + - fvm::laplacian(thCondEff,partTempField_) + == + fvOptions(partCpEff, partTempField_) + ); + // if transient add time derivative - need particle density and specific heat fields + // if sources activated add sources + // if convection activated add convection + + partTEqn.relax(); + + fvOptions.constrain(partTEqn); + + partTEqn.solve(); + + partTempField_.relax(); + + fvOptions.correct(partTempField_); + + if (pTMax_.value() > 0.0) partTempField_ = min(partTempField_, pTMax_); + if (pTMin_.value() > 0.0) partTempField_ = max(partTempField_, pTMin_); + + + Info<< "partT max/min : " << max(partTempField_).value() << " " << min(partTempField_).value() << endl; + +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H new file mode 100644 index 00000000..296c9a49 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H @@ -0,0 +1,108 @@ +/*---------------------------------------------------------------------------*\ +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 + + Description + Correlation for Nusselt number according to + Gunn, D. J. International Journal of Heat and Mass Transfer 21.4 (1978) + +\*---------------------------------------------------------------------------*/ + +#ifndef heatTransferGunnPartField_H +#define heatTransferGunnPartField_H + +#include "fvCFD.H" +#include "cfdemCloudEnergy.H" +#include "heatTransferGunn.H" +#include "fvOptions.H" +#include "scalarList.H" +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +class thermCondModel; +/*---------------------------------------------------------------------------*\ + Class heatTransferGunnPartField Declaration +\*---------------------------------------------------------------------------*/ + +class heatTransferGunnPartField +: + public heatTransferGunn +{ +private: + + volScalarField partCpField_; + + const volScalarField& partRhoField_; + + scalarList typeCp_; + + mutable double **partCp_; + + dimensionedScalar pTMax_; + + dimensionedScalar pTMin_; + + autoPtr thermCondModel_; + + fv::options& fvOptions; + + void allocateMyArrays() const; + + void giveData(); + +public: + + //- Runtime type information + TypeName("heatTransferGunnPartField"); + + // Constructors + + //- Construct from components + heatTransferGunnPartField + ( + const dictionary& dict, + cfdemCloudEnergy& sm + ); + + + // Destructor + + virtual ~heatTransferGunnPartField(); + + + // Member Functions + void addEnergyContribution(volScalarField&) const; + + void calcEnergyContribution(); + + void postFlow(); + + void solve(); + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C index 58677ad3..59c39254 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C @@ -12,6 +12,7 @@ License along with this code. If not, see . Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria + Contributing author: Paul Kieckhefen, TU Hamburg, Germany \*---------------------------------------------------------------------------*/ @@ -49,14 +50,26 @@ BeetstraDrag::BeetstraDrag : forceModel(dict,sm), propsDict_(dict.subDict(typeName + "Props")), + multiTypes_(false), velFieldName_(propsDict_.lookup("velFieldName")), U_(sm.mesh().lookupObject (velFieldName_)), voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")), voidfraction_(sm.mesh().lookupObject (voidfractionFieldName_)), + minVoidfraction_(propsDict_.lookupOrDefault("minVoidfraction",0.1)), UsFieldName_(propsDict_.lookup("granVelFieldName")), UsField_(sm.mesh().lookupObject (UsFieldName_)), scaleDia_(1.), - scaleDrag_(1.) + typeCG_(propsDict_.lookupOrDefault("coarseGrainingFactors",scalarList(1,1.0))), + scaleDrag_(1.), + rhoP_(0.), + rho_(0.), + Lc2_(0.), + dPrim_(0.), + nuf_(0.), + g_(9.81), + k_(0.05), + useGC_(false), + usePC_(false) { //Append the field names to be probed particleCloud_.probeM().initialize(typeName, typeName+".logDat"); @@ -78,10 +91,42 @@ BeetstraDrag::BeetstraDrag forceSubM(0).readSwitches(); particleCloud_.checkCG(true); - if (propsDict_.found("scale")) + if (propsDict_.found("scale") && typeCG_.size()==1) + { + // if "scale" is specified and there's only one single type, use "scale" scaleDia_=scalar(readScalar(propsDict_.lookup("scale"))); + typeCG_[0] = scaleDia_; + } if (propsDict_.found("scaleDrag")) + { scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag"))); + } + + if (typeCG_.size()>1) multiTypes_ = true; + + if (propsDict_.found("useFilteredDragModel")) + { + useGC_ = true; + g_=propsDict_.lookupOrDefault("g", 9.81); + dPrim_=scalar(readScalar(propsDict_.lookup("dPrim"))); + rhoP_=scalar(readScalar(propsDict_.lookup("rhoP"))); + rho_=scalar(readScalar(propsDict_.lookup("rho"))); + nuf_=scalar(readScalar(propsDict_.lookup("nuf"))); + scalar ut = terminalVelocity(1., dPrim_, nuf_, rho_, rhoP_, g_); + scalar Frp = ut*ut/g_/dPrim_; + Lc2_ = ut*ut/g_*pow(Frp, -.6666667); // n is hardcoded as -2/3 + Info << "using grid coarsening correction with Lc2 = " << Lc2_ << " and ut = " << ut << " and Frp = " << Frp<< endl; + + if (propsDict_.found("useParcelSizeDependentFilteredDrag")) + { + usePC_ = true; + if (propsDict_.found("k")) + k_=scalar(readScalar(propsDict_.lookup("k"))); + Info << "using particle coarsening correction with k = " << k_ << endl; + } + } + + } @@ -96,9 +141,9 @@ BeetstraDrag::~BeetstraDrag() void BeetstraDrag::setForce() const { - if (scaleDia_ > 1) + if (typeCG_.size()>1 || typeCG_[0] > 1) { - Info << "Beetstra using scale = " << scaleDia_ << endl; + Info << "Beetstra using scale = " << typeCG_ << endl; } else if (particleCloud_.cg() > 1) { @@ -119,12 +164,18 @@ void BeetstraDrag::setForce() const vector Ur(0,0,0); scalar ds(0); scalar ds_scaled(0); - scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_; + scalar dSauter(0); + scalar scaleDia3 = typeCG_[0]*typeCG_[0]*typeCG_[0]; scalar nuf(0); scalar rho(0); scalar magUr(0); scalar Rep(0); scalar localPhiP(0); + scalar GCcorr(1.); + scalar PCcorr(1.); + + scalar cg = typeCG_[0]; + label partType = 1; vector dragExplicit(0,0,0); scalar dragCoefficient(0); @@ -137,52 +188,79 @@ void BeetstraDrag::setForce() const for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; - drag = vector(0,0,0); - dragExplicit = vector(0,0,0); - Ufluid =vector(0,0,0); + drag = vector::zero; + dragExplicit = vector::zero; + Ufluid = vector::zero; voidfraction=0; dragCoefficient = 0; if (cellI > -1) // particle found { - if( forceSubM(0).interpolation() ) + if ( forceSubM(0).interpolation() ) { position = particleCloud_.position(index); voidfraction = voidfractionInterpolator_.interpolate(position,cellI); Ufluid = UInterpolator_.interpolate(position,cellI); //Ensure interpolated void fraction to be meaningful // Info << " --> voidfraction: " << voidfraction << endl; - if(voidfraction>1.00) voidfraction = 1.0; - if(voidfraction<0.10) voidfraction = 0.10; + if (voidfraction > 1.00) voidfraction = 1.0; + if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_; } else { voidfraction = voidfraction_[cellI]; Ufluid = U_[cellI]; } + // in case a fines phase is present, void fraction needs to be adapted + adaptVoidfraction(voidfraction, cellI); + + if (multiTypes_) + { + partType = particleCloud_.particleType(index); + cg = typeCG_[partType - 1]; + scaleDia3 = cg*cg*cg; + } Us = particleCloud_.velocity(index); Ur = Ufluid-Us; magUr = mag(Ur); ds = 2*particleCloud_.radius(index); - ds_scaled = ds/scaleDia_; + ds_scaled = ds/cg; + dSauter = meanSauterDiameter(ds_scaled, cellI); + rho = rhoField[cellI]; nuf = nufField[cellI]; - Rep=0.0; localPhiP = 1.0f-voidfraction+SMALL; // calc particle's drag coefficient (i.e., Force per unit slip velocity and Stokes drag) - Rep=ds_scaled*voidfraction*magUr/nuf+SMALL; - dragCoefficient = 10.0*localPhiP/(voidfraction*voidfraction) + - voidfraction*voidfraction*(1.0+1.5*Foam::sqrt(localPhiP)) + - 0.413*Rep/(24*voidfraction*voidfraction)*(1.0/voidfraction+3*voidfraction*localPhiP+8.4*Foam::pow(Rep,-0.343))/ - (1+Foam::pow(10,3*localPhiP)*Foam::pow(Rep,-0.5*(1+4*localPhiP))); + Rep=dSauter*voidfraction*magUr/nuf + SMALL; + + dragCoefficient = F(voidfraction, Rep) + *3*M_PI*nuf*rho*voidfraction + *effDiameter(ds_scaled, cellI, index) + *scaleDia3*scaleDrag_; + + // calculate filtering corrections + if (useGC_) + { + GCcorr = 1.-h(localPhiP) + /( a(localPhiP) + *Lc2_ + /std::cbrt(U_.mesh().V()[cellI]) + + 1. + ); + if (usePC_) + { + PCcorr = Foam::exp(k_*(1.-cg)); + } + } + + // apply filtering corrections + dragCoefficient *= GCcorr*PCcorr; - // calc particle's drag - dragCoefficient *= 3*M_PI*ds_scaled*nuf*rho*voidfraction*scaleDia3*scaleDrag_; if (modelType_=="B") dragCoefficient /= voidfraction; @@ -198,11 +276,13 @@ void BeetstraDrag::setForce() const Pout << "Us = " << Us << endl; Pout << "Ur = " << Ur << endl; Pout << "ds = " << ds << endl; - Pout << "ds/scale = " << ds/scaleDia_ << endl; + Pout << "ds/scale = " << ds/cg << endl; Pout << "rho = " << rho << endl; Pout << "nuf = " << nuf << endl; Pout << "voidfraction = " << voidfraction << endl; Pout << "Rep = " << Rep << endl; + Pout << "GCcorr = " << GCcorr << endl; + Pout << "PCcorr = " << PCcorr << endl; Pout << "drag = " << drag << endl; } @@ -223,7 +303,149 @@ void BeetstraDrag::setForce() const } } +/********************************************************* + * "Drag Force of Intermediate Reynolds Number Flow Past * + * Mono- and Bidisperse Arrays of Spheres", eq. 16 * + * R Beetstra, M. A. van der Hoef, JAM Kuipers * + * AIChE Journal 53(2) (2007) * + *********************************************************/ +double BeetstraDrag::F(double voidfraction, double Rep) const +{ + double localPhiP = max(SMALL,min(1.-SMALL,1.-voidfraction)); + return 10.0*localPhiP/(voidfraction*voidfraction) + + voidfraction*voidfraction*(1.0+1.5*Foam::sqrt(localPhiP)) + + 0.413*Rep/(24*voidfraction*voidfraction) + *(1.0/voidfraction + +3*voidfraction*localPhiP + +8.4*Foam::pow(Rep,-0.343) + ) + /(1+Foam::pow(10,3*localPhiP) + *Foam::pow(Rep,-0.5*(1+4*localPhiP)) + ); +} + + +/********************************************************* + * "A drag model for filtered Euler-Lagange simulations * + * of clustered gas-particle suspension", * + * S. Radl, S. Sundaresan, * + * Chemical Engineering Science 117 (2014). * + *********************************************************/ +double BeetstraDrag::a(double phiP) const +{ + double a0m = 0.; + double a1m = 0.; + double a2m = 0.; + double a3m = 0.; + double phipam = 0.; + if (phiP < 0.016) + { + a0m = 21.51; + } + else if (phiP < 0.100) + { + a0m = 1.96; a1m = 29.40; a2m = 164.91; a3m = -1923.; + } + else if (phiP < 0.180) + { + a0m = 4.63; a1m = 4.68; a2m = -412.04; a3m = 2254.; phipam = 0.10; + } + else if (phiP < 0.250) + { + a0m = 3.52; a1m = -17.99; a2m = 128.80; a3m = -603.; phipam = 0.18; + } + else if (phiP < 0.400) + { + a0m = 2.68; a1m = -8.20; a2m = 2.18; a3m = 112.33; phipam = 0.25; + } + else + { + a0m = 1.79; + } + const double phiP_minus_phipam = phiP-phipam; + return a0m + a1m*phiP_minus_phipam + a2m*phiP_minus_phipam*phiP_minus_phipam + + a3m*phiP_minus_phipam*phiP_minus_phipam*phiP_minus_phipam; +} + +double BeetstraDrag::h(double phiP) const +{ + double h0m = 0.; + double h1m = 0.; + double h2m = 0.; + double h3m = 0.; + double phiphm = 0.; + + if (phiP < 0.03) + { + h1m = 7.97; + } + else if (phiP < 0.08) + { + h0m = 0.239; h1m = 4.640; h2m = -4.410; h3m = 253.630; phiphm = 0.03; + } + else if (phiP < 0.12) + { + h0m = 0.492; h1m = 6.100; h2m = 33.630; h3m = -789.600; phiphm = 0.08; + } + else if (phiP < 0.18) + { + h0m = 0.739; h1m = 5.010; h2m = -61.100; h3m = 310.800; phiphm = 0.12; + } + else if (phiP < 0.34) + { + h0m = 0.887; h1m = 1.030; h2m = -5.170; h3m = 5.990; phiphm = 0.18; + } + else if (phiP < 0.48) + { + h0m = 0.943; h1m = -0.170; h2m = -2.290; h3m = -9.120; phiphm = 0.34; + } + else if (phiP < 0.55) + { + h0m = 0.850; h1m = -1.350; h2m = -6.130; h3m = -132.600; phiphm = 0.48; + } + else + { + h0m = 0.680; h1m = -2.340; h2m = -225.200; phiphm = 0.55; + } + const double phiP_minus_phiphm = phiP-phiphm; + return h0m + h1m*phiP_minus_phiphm + h2m*phiP_minus_phiphm*phiP_minus_phiphm + + h3m*phiP_minus_phiphm*phiP_minus_phiphm*phiP_minus_phiphm; +} + +double BeetstraDrag::terminalVelocity(double voidfraction, double dp, double nuf, double rhof, double rhop, double g) const +{ + scalar u0(dp*dp*fabs(rhof-rhop)*g/18./rhof/nuf); + scalar Re(u0*dp/nuf); + scalar res(1.); + scalar u(u0); + scalar Fi(0); + scalar CdSt(0); + Info << "uo: " << u0< 1.e-6) && (i<100)) + { + Info << "Iteration " << i; + u0 = u; + Info << ", u0 = " << u0; + CdSt = 24/Re; + Info << ", CdSt = " << CdSt; + Fi = F(voidfraction, Re); + Info << ", F = "; + u = sqrt(1.333333333*fabs(rhof-rhop)*g*dp + /(CdSt*voidfraction*Fi*rhof) + ); + Info << ", u = " << u; + Re = fabs(u)*dp/nuf*voidfraction; + res = fabs((u-u0)/u); + Info << "Res: " << res << endl; + i++; + } + if (res >1.e-6) + FatalError << "Terminal velocity calculation diverged!" << endl; + + return u; +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H index f6e079ab..c2128855 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H @@ -40,9 +40,11 @@ class BeetstraDrag : public forceModel { -private: +protected: dictionary propsDict_; + bool multiTypes_; + word velFieldName_; const volVectorField& U_; @@ -51,14 +53,52 @@ private: const volScalarField& voidfraction_; + const scalar minVoidfraction_; + word UsFieldName_; const volVectorField& UsField_; mutable scalar scaleDia_; + scalarList typeCG_; + mutable scalar scaleDrag_; + mutable scalar rhoP_; + + mutable scalar rho_; + + mutable scalar Lc2_; + + mutable scalar dPrim_; + + mutable scalar nuf_; + + mutable scalar g_; + + mutable scalar k_; + + bool useGC_; + + 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;} + + double F(double, double) const; + + double terminalVelocity(double, double, double, double, double, double) const; + + double a(double) const; + + double h(double) const; + + + public: //- Runtime type information diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C new file mode 100644 index 00000000..04e12ed9 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C @@ -0,0 +1,116 @@ +/*---------------------------------------------------------------------------*\ +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 + +\*---------------------------------------------------------------------------*/ + +#include "error.H" + +#include "BeetstraDragPoly.H" +#include "addToRunTimeSelectionTable.H" +#include "averagingModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(BeetstraDragPoly, 0); + +addToRunTimeSelectionTable +( + forceModel, + BeetstraDragPoly, + dictionary +); + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +BeetstraDragPoly::BeetstraDragPoly +( + const dictionary& dict, + cfdemCloud& sm +) +: + BeetstraDrag(dict,sm), + fines_(propsDict_.lookupOrDefault("fines",false)), + dFine_(1.0) +{ + // if fines are present, take mixture dSauter, otherwise normal dSauter + if (fines_) + { + dFine_ = readScalar(propsDict_.lookup("dFine")); + volScalarField& alphaP(const_cast(sm.mesh().lookupObject ("alphaP"))); + alphaP_.set(&alphaP); + volScalarField& alphaSt(const_cast(sm.mesh().lookupObject ("alphaSt"))); + alphaSt_.set(&alphaSt); + volScalarField& dSauter(const_cast(sm.mesh().lookupObject ("dSauterMix"))); + dSauter_.set(&dSauter); + } + else + { + volScalarField& dSauter(const_cast(sm.mesh().lookupObject ("dSauter"))); + dSauter_.set(&dSauter); + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +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]; + scalar effD = d*d / dS + 0.064*d*d*d*d / (dS*dS*dS); + + if (fines_) + { + scalar fineCorr = dFine_*dFine_ / dS + 0.064*dFine_*dFine_*dFine_*dFine_ / (dS*dS*dS); + fineCorr *= d*d*d / (dFine_*dFine_*dFine_) * alphaSt_()[cellI] / alphaP_()[cellI]; + effD += fineCorr; + } + + if (particleCloud_.getParticleEffVolFactors()) + { + scalar effVolFac = particleCloud_.particleEffVolFactor(index); + effD *= effVolFac; + } + + return effD; +} + +scalar BeetstraDragPoly::meanSauterDiameter(double d, label cellI) const +{ + return dSauter_()[cellI]; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H new file mode 100644 index 00000000..88965010 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H @@ -0,0 +1,93 @@ +/*---------------------------------------------------------------------------*\ +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 + +Description + drag law for monodisperse systems according to + Beetstra et al. AIChE J 53.2 (2007) + +SourceFiles + BeetstraDragPoly.C +\*---------------------------------------------------------------------------*/ + +#ifndef BeetstraDragPoly_H +#define BeetstraDragPoly_H + +#include "BeetstraDrag.H" +#include "interpolationCellPoint.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class BeetstraDragPoly Declaration +\*---------------------------------------------------------------------------*/ + +class BeetstraDragPoly +: + public BeetstraDrag +{ +protected: + + const bool fines_; + + scalar dFine_; + + autoPtr alphaP_; + + autoPtr alphaSt_; + + autoPtr dSauter_; + + void adaptVoidfraction(double&, label) const; + + scalar effDiameter(double, label, label) const; + + scalar meanSauterDiameter(double, label) const; + +public: + + //- Runtime type information + TypeName("BeetstraDragPoly"); + + + // Constructors + + //- Construct from components + BeetstraDragPoly + ( + const dictionary& dict, + cfdemCloud& sm + ); + + // Destructor + + ~BeetstraDragPoly(); + + + // Member Functions +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C index c15da1e1..8a7deab9 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/DiFeliceDrag/DiFeliceDrag.C @@ -148,10 +148,10 @@ void DiFeliceDrag::setForce() const for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; - drag = vector(0,0,0); - dragExplicit = vector(0,0,0); + drag = vector::zero; + dragExplicit = vector::zero; dragCoefficient=0; - Ufluid =vector(0,0,0); + Ufluid = vector::zero; if (cellI > -1) // particle Found { diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C index e229c15e..08fbf8b3 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C @@ -164,10 +164,10 @@ void ErgunStatFines::setForce() const for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; - drag = vector(0,0,0); - dragExplicit = vector(0,0,0); + drag = vector::zero; + dragExplicit = vector::zero; betaP = 0; - Ufluid = vector(0,0,0); + Ufluid = vector::zero; voidfraction = 0; dragCoefficient = 0; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C index 29cbdd36..0dffe38b 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C @@ -120,8 +120,8 @@ void FanningDynFines::setForce() const for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; - drag = vector(0,0,0); - UDyn = vector(0,0,0); + drag = vector::zero; + UDyn = vector::zero; dragCoefficient = 0; if (cellI > -1) // particle found diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C index 57bb5d3c..826d2a76 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/GidaspowDrag/GidaspowDrag.C @@ -159,11 +159,11 @@ void GidaspowDrag::setForce() const //if(mask[index][0]) //{ cellI = particleCloud_.cellIDs()[index][0]; - drag = vector(0,0,0); - dragExplicit = vector(0,0,0); + drag = vector::zero; + dragExplicit = vector::zero; betaP = 0; Vs = 0; - Ufluid = vector(0,0,0); + Ufluid = vector::zero; voidfraction = 0; dragCoefficient = 0; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C index 7368d942..23f866ab 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillDrag/KochHillDrag.C @@ -157,12 +157,12 @@ void KochHillDrag::setForce() const for (int index = 0; index -1) // particle Found diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C index 616069c8..5303c7fd 100755 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/KochHillRWDrag/KochHillRWDrag.C @@ -200,12 +200,12 @@ void KochHillRWDrag::setForce() const //if (mask[index][0]) //{ cellI = particleCloud_.cellIDs()[index][0]; - drag = vector(0,0,0); - dragExplicit = vector(0,0,0); + drag = vector::zero; + dragExplicit = vector::zero; dragCoefficient = 0; betaP = 0; Vs = 0; - Ufluid = vector(0,0,0); + Ufluid = vector::zero; // Pout << "RW-TEST: cellI = " << cellI << endl; // TEST-Output if (cellI > -1) // particle Found diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C index 081742c3..9f1a555b 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C @@ -53,10 +53,10 @@ dSauter::dSauter : forceModel(dict,sm), propsDict_(dict.subDict(typeName + "Props")), + multiTypes_(false), d2_(NULL), d3_(NULL), - scaleDia_(1.0), - scaleDiaDist_(1.0), + typeCG_(propsDict_.lookupOrDefault("coarseGrainingFactors",scalarList(1,1.0))), d2Field_ ( IOobject ( @@ -95,17 +95,13 @@ dSauter::dSauter "zeroGradient" ) { + if (typeCG_.size()>1) multiTypes_ = true; allocateMyArrays(); dSauter_.write(); // init force sub model setForceSubModels(propsDict_); - - if (propsDict_.found("scaleCG")) - scaleDia_ = scalar(readScalar(propsDict_.lookup("scaleCG"))); - if (propsDict_.found("scaleDist")) - scaleDiaDist_ = scalar(readScalar(propsDict_.lookup("scaleDist"))); } @@ -131,30 +127,37 @@ void dSauter::allocateMyArrays() const void dSauter::setForce() const { - if (scaleDia_ > 1) + if (typeCG_.size()>1 || typeCG_[0] > 1.0) { - Info << "dSauter using scaleCG = " << scaleDia_ << endl; - } - else if (particleCloud_.cg() > 1) - { - scaleDia_ = particleCloud_.cg(); - Info << "dSauter using scaleCG from liggghts cg = " << scaleDia_ << endl; + Info << "dSauter using CG factor(s) = " << typeCG_ << endl; } allocateMyArrays(); - label cellI=0; - scalar ds(0); - scalar scale = scaleDiaDist_/scaleDia_; + label cellI = 0; + label partType = 1; + scalar cg = typeCG_[0]; + scalar ds = 0.0; + scalar effVolFac = 1.0; + for(int index = 0; index < particleCloud_.numberOfParticles(); ++index) { cellI = particleCloud_.cellIDs()[index][0]; if (cellI >= 0) { + if (particleCloud_.getParticleEffVolFactors()) + { + effVolFac = particleCloud_.particleEffVolFactor(index); + } + if (multiTypes_) + { + partType = particleCloud_.particleType(index); + cg = typeCG_[partType - 1]; + } ds = particleCloud_.d(index); - d2_[index][0] = ds*ds; - d3_[index][0] = ds*ds*ds; + d2_[index][0] = ds*ds*effVolFac*cg; + d3_[index][0] = ds*ds*ds*effVolFac; } } @@ -181,7 +184,7 @@ void dSauter::setForce() const { if (d2Field_[cellI] > ROOTVSMALL) { - dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI] * scale; + dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI]; } else { diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H index 7a81fab5..f1e1ca6e 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.H @@ -43,13 +43,13 @@ private: dictionary propsDict_; + bool multiTypes_; + mutable double **d2_; mutable double **d3_; - mutable scalar scaleDia_; // scaling due to coarse graining - - mutable scalar scaleDiaDist_; // scaling due to modified particle size distribution + scalarList typeCG_; mutable volScalarField d2Field_; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C index 0b593924..4abb55af 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C @@ -69,7 +69,7 @@ forceModel::forceModel IOobject::AUTO_WRITE ), sm.mesh(), - dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector(0,0,0)) // N + dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero) // N ), expParticleForces_ ( IOobject @@ -81,7 +81,7 @@ forceModel::forceModel IOobject::AUTO_WRITE ), sm.mesh(), - dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector(0,0,0)) // N + dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero) // N ), coupleForce_(true), modelType_(sm.modelType()), diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C index 586579b7..ed8481ad 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/gradPForce/gradPForce.C @@ -160,7 +160,7 @@ void gradPForce::setForce() const { //if(mask[index][0]) //{ - force = vector(0,0,0); + force = vector::zero; cellI = particleCloud_.cellIDs()[index][0]; if (cellI > -1) // particle Found diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.C b/src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.C index ad77200e..33b44627 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/interface/interface.C @@ -142,7 +142,7 @@ Info << "interface::setForce" << endl; } // Initialize an interfaceForce vector - vector interfaceForce = Foam::vector(0,0,0); + vector interfaceForce = Foam::vector::zero; // Calculate the interfaceForce (range of alphap needed for stability) diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C new file mode 100644 index 00000000..ff0de5d0 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.C @@ -0,0 +1,317 @@ +/*---------------------------------------------------------------------------*\ +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 . + +Description + calculates the correlation between particle momentum and diameter + +SourceFiles + pdCorrelation.C + +Contributing Author + 2018 Paul Kieckhefen, TUHH +\*---------------------------------------------------------------------------*/ + +#include "error.H" + +#include "pdCorrelation.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(pdCorrelation, 0); + +addToRunTimeSelectionTable +( + forceModel, + pdCorrelation, + dictionary +); + +const dimensionSet dimMomentum(dimForce*dimTime); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +pdCorrelation::pdCorrelation +( + const dictionary& dict, + cfdemCloud& sm +) +: + forceModel(dict,sm), + propsDict_(dict.subDict(typeName + "Props")), + d_(nullptr), + p_(nullptr), + d2_(nullptr), + pd_(nullptr), + cg3_(nullptr), + dField_ + ( IOobject + ( + "dField", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + sm.mesh(), + dimensionedScalar("zero", dimLength, 0), + "zeroGradient" + ), + pField_ + ( IOobject + ( + "pField", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + sm.mesh(), + dimensionedVector("zero", dimMomentum, vector::zero) + ), + d2Field_ + ( IOobject + ( + "d2Field", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + sm.mesh(), + dimensionedScalar("zero", dimLength*dimLength, 0) + ), + pdField_ + ( IOobject + ( + "pdCorrelation", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + sm.mesh(), + dimensionedVector("zero", dimMomentum*dimLength, vector::zero), + "zeroGradient" + ), + cg3Field_ + ( IOobject + ( + "cg3Field", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + sm.mesh(), + dimensionedScalar("zero", dimless, 0) + ), + typeCG_ + ( + propsDict_.lookupOrDefault + ( + "coarseGrainingFactors", + scalarList(1, sm.cg()) + ) + ), + particleDensities_ + ( + propsDict_.lookupOrDefault + ( + "particleDensities", + scalarList(1, -1.) + ) + ), + constantCG_(typeCG_.size() < 2), + CG_(!constantCG_ || typeCG_[0] > 1. + SMALL), + runOnWriteOnly_(propsDict_.lookupOrDefault("runOnWriteOnly", false)) +{ + if ((particleDensities_[0] < 0) && !particleCloud_.getParticleDensities()) + { + FatalError<< "Please set the particle density either in LIGGGHTS" + << "or the pdCorrelationPropsDict." + << "In the first case, set getParticleDensities in the couplingProperties." + << abort(FatalError); + } + + allocateMyArrays(); + + dField_.write(); + pdField_.write(); + + + // init force sub model + setForceSubModels(propsDict_); +} + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +pdCorrelation::~pdCorrelation() +{ + particleCloud_.dataExchangeM().destroy(cg3_, 1); + particleCloud_.dataExchangeM().destroy(d_, 1); + particleCloud_.dataExchangeM().destroy(p_, 3); + particleCloud_.dataExchangeM().destroy(d2_, 1); + particleCloud_.dataExchangeM().destroy(pd_, 3); +} + +// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // +void pdCorrelation::allocateMyArrays() const +{ + // get memory for 2d arrays + double initVal = 0.0; + particleCloud_.dataExchangeM().allocateArray(d_, initVal, 1); + particleCloud_.dataExchangeM().allocateArray(p_, initVal, 3); + particleCloud_.dataExchangeM().allocateArray(d2_, initVal, 1); + particleCloud_.dataExchangeM().allocateArray(pd_, initVal, 3); + particleCloud_.dataExchangeM().allocateArray(cg3_, initVal, 1); +} +// * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * // + +void pdCorrelation::setForce() const +{ + const fvMesh& mesh = particleCloud_.mesh(); + + if (runOnWriteOnly_ && !mesh.write()) return; // skip if it's not write time + + allocateMyArrays(); + + const Switch densityFromList + ( + particleCloud_.getParticleTypes() + && !particleCloud_.getParticleDensities() + ); + const scalar piOverSix = constant::mathematical::pi / 6.; + + dField_.primitiveFieldRef() = 0.0; + pField_.primitiveFieldRef() = vector::zero; + d2Field_.primitiveFieldRef() = 0.0; + pdField_.primitiveFieldRef() = vector::zero; + cg3Field_.primitiveFieldRef() = 0.; + + label celli(0); + scalar dp(0.); + scalar cg(typeCG_[0]); + scalar rhop(particleDensities_[0]); + label typei(0); + scalar particleV(0.); + for(int pi = 0; pi < particleCloud_.numberOfParticles(); ++pi) + { + celli = particleCloud_.cellIDs()[pi][0]; + if (celli < 0) continue; + + dp = particleCloud_.d(pi); + if (particleCloud_.getParticleTypes()) + { + typei = particleCloud_.particleTypes()[pi][0] - 1; + rhop = densityFromList ? particleDensities_[typei] + : particleCloud_.particleDensity(pi); + cg = constantCG_ ? typeCG_[0] : typeCG_[typei]; + } + + particleV = piOverSix * rhop * dp * dp * dp; + d_ [pi][0] = cg * cg * dp; + d2_[pi][0] = cg * dp * dp; + p_ [pi][0] = particleV * particleCloud_.velocities()[pi][0]; + p_ [pi][1] = particleV * particleCloud_.velocities()[pi][1]; + p_ [pi][2] = particleV * particleCloud_.velocities()[pi][2]; + pd_[pi][0] = p_[pi][0] * dp / cg; + pd_[pi][1] = p_[pi][1] * dp / cg; + pd_[pi][2] = p_[pi][2] * dp / cg; + + cg3_[pi][0] = CG_ ? cg * cg * cg : 1.; + } + + + particleCloud_.averagingM().setScalarSum + ( + dField_, + d_, + particleCloud_.particleWeights(), + nullptr + ); + particleCloud_.averagingM().setVectorSum + ( + pField_, + p_, + particleCloud_.particleWeights(), + nullptr + ); + particleCloud_.averagingM().setScalarSum + ( + d2Field_, + d2_, + particleCloud_.particleWeights(), + nullptr + ); + particleCloud_.averagingM().setVectorSum + ( + pdField_, + pd_, + particleCloud_.particleWeights(), + nullptr + ); + + particleCloud_.averagingM().setScalarSum + ( + cg3Field_, + cg3_, + particleCloud_.particleWeights(), + nullptr + ); + + scalar oneOverCG3(1.); + const scalar oneMinusSmall(1.-SMALL); + forAll(cg3Field_, celli) + { + if (cg3Field_[celli] < oneMinusSmall) continue; + oneOverCG3 = 1./cg3Field_[celli]; + + dField_ [celli] *= oneOverCG3; + pField_ [celli] *= oneOverCG3; + d2Field_[celli] *= oneOverCG3; + pdField_[celli] *= oneOverCG3; + } + + scalar denom(0.); + forAll(dField_, celli) + { + if (dField_[celli] < SMALL) + { + pdField_[celli] = vector::zero; + continue; + } + denom = d2Field_[celli] - dField_[celli] * dField_[celli]; + if (denom < SMALL) + { + pdField_[celli] = vector::zero; + continue; + } + + pdField_[celli] -= dField_[celli] * pField_[celli]; + pdField_[celli] /= denom; + } +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H new file mode 100644 index 00000000..7ea99050 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/pdCorrelation/pdCorrelation.H @@ -0,0 +1,102 @@ +/*---------------------------------------------------------------------------*\ +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 . + + +Description + calculates the correlation between particle diameter and momentum + +SourceFiles + pdCorrelation.C + +Contributing Author + 2018 Paul Kieckhefen, TUHH +\*---------------------------------------------------------------------------*/ + +#ifndef pdCorrelation_H +#define pdCorrelation_H + +#include "forceModel.H" +#include "averagingModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class pdCorrelation Declaration +\*---------------------------------------------------------------------------*/ + +class pdCorrelation +: + public forceModel +{ +private: + + dictionary propsDict_; + + mutable double **d_; + mutable double **p_; + mutable double **d2_; + mutable double **pd_; + mutable double **cg3_; + + mutable volScalarField dField_; + mutable volVectorField pField_; + mutable volScalarField d2Field_; + mutable volVectorField pdField_; + + mutable volScalarField cg3Field_; + + const scalarList typeCG_; + const scalarList particleDensities_; + const Switch constantCG_; + const Switch CG_; + const Switch runOnWriteOnly_; + + void allocateMyArrays() const; + +public: + + //- Runtime type information + TypeName("pdCorrelation"); + + // Constructors + + //- Construct from components + pdCorrelation + ( + const dictionary& dict, + cfdemCloud& sm + ); + + // Destructor + + ~pdCorrelation(); + + + // Member Functions + void setForce() const; + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/viscForce/viscForce.C b/src/lagrangian/cfdemParticle/subModels/forceModel/viscForce/viscForce.C index 30762ce8..23b3298a 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/viscForce/viscForce.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/viscForce/viscForce.C @@ -142,7 +142,7 @@ void viscForce::setForce() const { //if(mask[index][0]) //{ - force=vector(0,0,0); + force = vector::zero; cellI = particleCloud_.cellIDs()[index][0]; if (cellI > -1) // particle Found diff --git a/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C b/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C index 7cb38d76..1bebcc2a 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C @@ -156,7 +156,7 @@ void DiFeliceDragMS::setForce() const //{ cellI = cloudRefMS().cellIDCM(index); - drag = vector(0,0,0); + drag = vector::zero; if (cellI > -1) // particle Found { diff --git a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C index 46e99223..da5238ac 100644 --- a/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C +++ b/src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C @@ -73,7 +73,7 @@ explicitCouple::explicitCouple IOobject::NO_WRITE ), sm.mesh(), - dimensionedVector("zero", dimensionSet(1,-2,-2,0,0), vector(0,0,0)), // N/m3 + dimensionedVector("zero", dimensionSet(1,-2,-2,0,0), vector::zero), // N/m3 "zeroGradient" ), fNext_ @@ -86,7 +86,7 @@ explicitCouple::explicitCouple IOobject::NO_WRITE ), sm.mesh(), - dimensionedVector("zero", dimensionSet(1,-2,-2,0,0), vector(0,0,0)), // N/m3 + dimensionedVector("zero", dimensionSet(1,-2,-2,0,0), vector::zero), // N/m3 "zeroGradient" ), fLimit_(1e10,1e10,1e10) diff --git a/src/lagrangian/cfdemParticle/subModels/otherForceModel/expParticleForces/expParticleForces.C b/src/lagrangian/cfdemParticle/subModels/otherForceModel/expParticleForces/expParticleForces.C index e7818931..65bbe71d 100644 --- a/src/lagrangian/cfdemParticle/subModels/otherForceModel/expParticleForces/expParticleForces.C +++ b/src/lagrangian/cfdemParticle/subModels/otherForceModel/expParticleForces/expParticleForces.C @@ -69,18 +69,18 @@ tmp expParticleForces::exportForceField() particleCloud_.mesh(), dimensionedVector ( - "zero",dimensionSet(1, -2, -2, 0, 0),vector(0,0,0) + "zero",dimensionSet(1, -2, -2, 0, 0),vector::zero ) ) ); - + volVectorField& source = tsource.ref(); - + // negative sign in sum because force on particles = - force on fluid for(int i=0; i gravity::exportForceField() particleCloud_.mesh(), dimensionedVector ( - "zero",dimensionSet(1, -2, -2, 0, 0),vector(0,0,0) + "zero",dimensionSet(1, -2, -2, 0, 0),vector::zero ) ) ); - + volVectorField& source = tsource.ref(); - + source = rhoG_ * voidfraction_ * g_; - + return tsource; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/otherForceModel/gravity/gravity.H b/src/lagrangian/cfdemParticle/subModels/otherForceModel/gravity/gravity.H index b9e650b5..24efa9a7 100644 --- a/src/lagrangian/cfdemParticle/subModels/otherForceModel/gravity/gravity.H +++ b/src/lagrangian/cfdemParticle/subModels/otherForceModel/gravity/gravity.H @@ -48,14 +48,13 @@ protected: word voidfractionFieldName_; const volScalarField& voidfraction_; - + word rhoGFieldName_; const volScalarField& rhoG_; const dimensionedVector g_; - public: //- Runtime type information diff --git a/src/lagrangian/cfdemParticle/subModels/otherForceModel/otherForceModel/otherForceModel.H b/src/lagrangian/cfdemParticle/subModels/otherForceModel/otherForceModel/otherForceModel.H index 5b7c2ac5..0813ba38 100644 --- a/src/lagrangian/cfdemParticle/subModels/otherForceModel/otherForceModel/otherForceModel.H +++ b/src/lagrangian/cfdemParticle/subModels/otherForceModel/otherForceModel/otherForceModel.H @@ -47,7 +47,6 @@ protected: cfdemCloud& particleCloud_; - public: //- Runtime type information @@ -95,7 +94,6 @@ public: // Member Functions - virtual tmp exportForceField() = 0; }; diff --git a/src/lagrangian/cfdemParticle/subModels/otherForceModel/weightSecondaryPhase/weightSecondaryPhase.C b/src/lagrangian/cfdemParticle/subModels/otherForceModel/weightSecondaryPhase/weightSecondaryPhase.C index 6e6fba8d..69d5ef63 100644 --- a/src/lagrangian/cfdemParticle/subModels/otherForceModel/weightSecondaryPhase/weightSecondaryPhase.C +++ b/src/lagrangian/cfdemParticle/subModels/otherForceModel/weightSecondaryPhase/weightSecondaryPhase.C @@ -78,7 +78,7 @@ tmp weightSecondaryPhase::exportForceField() particleCloud_.mesh(), dimensionedVector ( - "zero",dimensionSet(1, -2, -2, 0, 0),vector(0,0,0) + "zero",dimensionSet(1, -2, -2, 0, 0),vector::zero ) ) ); diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C new file mode 100644 index 00000000..7e4ca613 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C @@ -0,0 +1,225 @@ +/*---------------------------------------------------------------------------*\ +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 + +\*---------------------------------------------------------------------------*/ + +#include "error.H" +#include "ZehnerSchluenderThermCond.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(ZehnerSchluenderThermCond, 0); + +addToRunTimeSelectionTable +( + thermCondModel, + ZehnerSchluenderThermCond, + dictionary +); + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +ZehnerSchluenderThermCond::ZehnerSchluenderThermCond +( + const dictionary& dict, + cfdemCloud& sm +) +: + thermCondModel(dict,sm), + propsDict_(dict.subDict(typeName + "Props")), + partKsField_ + ( + IOobject + ( + "partKs", + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + sm.mesh(), + dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0), + "zeroGradient" + ), + voidfractionFieldName_(propsDict_.lookupOrDefault("voidfractionFieldName","voidfraction")), + voidfraction_(sm.mesh().lookupObject (voidfractionFieldName_)), + typeKs_(propsDict_.lookupOrDefault("thermalConductivities",scalarList(1,-1.0))), + partKs_(NULL), + wallQFactorName_(propsDict_.lookupOrDefault("wallQFactorName","wallQFactor")), + wallQFactor_ + ( IOobject + ( + wallQFactorName_, + sm.mesh().time().timeName(), + sm.mesh(), + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + sm.mesh(), + dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 1.0) + ), + hasWallQFactor_(false) +{ + if (typeKs_[0] < 0.0) + { + FatalError << "ZehnerSchluenderThermCond: provide list of thermal conductivities." << abort(FatalError); + } + + if (typeKs_.size() > 1) allocateMyArrays(); + else + { + partKsField_ *= typeKs_[0]; + } + + if (wallQFactor_.headerOk()) + { + Info << "Found field for scaling wall heat flux.\n" << endl; + hasWallQFactor_ = true; + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +ZehnerSchluenderThermCond::~ZehnerSchluenderThermCond() +{ + if (typeKs_.size() > 1) particleCloud_.dataExchangeM().destroy(partKs_,1); +} + + +void ZehnerSchluenderThermCond::allocateMyArrays() const +{ + double initVal=0.0; + particleCloud_.dataExchangeM().allocateArray(partKs_,initVal,1); +} + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +tmp ZehnerSchluenderThermCond::thermCond() const +{ + tmp tvf + ( + new volScalarField + ( + IOobject + ( + "tmpThCond", + voidfraction_.instance(), + voidfraction_.mesh(), + IOobject::NO_READ, + IOobject::NO_WRITE + ), + voidfraction_.mesh(), + dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0), + "zeroGradient" + ) + ); + + volScalarField& svf = tvf.ref(); + + calcPartKsField(); + scalar A = 0.0; + scalar B = 0.0; + scalar C = 0.0; + scalar k = 0.0; + scalar InvOnemBoA = 0.0; + scalar voidfraction = 0.0; + scalar w = 7.26e-3; + + forAll(svf, cellI) + { + voidfraction = voidfraction_[cellI]; + if(voidfraction > 1.0 - SMALL || partKsField_[cellI] < SMALL) svf[cellI] = 0.0; + else + { + A = partKsField_[cellI]/kf0_.value(); + B = 1.25 * Foam::pow((1 - voidfraction) / voidfraction, 1.11); + InvOnemBoA = 1.0/(1.0 - B/A); + C = (A - 1) * InvOnemBoA * InvOnemBoA * B/A * log(A/B) - (B - 1) * InvOnemBoA - 0.5 * (B + 1); + C *= 2.0 * InvOnemBoA; + k = Foam::sqrt(1 - voidfraction) * (w * A + (1 - w) * C) * kf0_.value(); + svf[cellI] = k / (1 - voidfraction); + } + } + + svf.correctBoundaryConditions(); + + // if a wallQFactor field is present, use it to scale heat transport through a patch + if (hasWallQFactor_) + { + wallQFactor_.correctBoundaryConditions(); + forAll(wallQFactor_.boundaryField(), patchi) + svf.boundaryFieldRef()[patchi] *= wallQFactor_.boundaryField()[patchi]; + } + + return tvf; +} + +tmp ZehnerSchluenderThermCond::thermDiff() const +{ + FatalError << "ZehnerSchluenderThermCond does not provide thermal diffusivity." << abort(FatalError); + return tmp(NULL); +} + +void ZehnerSchluenderThermCond::calcPartKsField() const +{ + if (typeKs_.size() <= 1) return; + + if (!particleCloud_.getParticleTypes()) + { + FatalError << "ZehnerSchluenderThermCond needs data for more than one type, but types are not communicated." << abort(FatalError); + } + + allocateMyArrays(); + label cellI=0; + label partType = 0; + for(int index = 0;index < particleCloud_.numberOfParticles(); ++index) + { + cellI = particleCloud_.cellIDs()[index][0]; + if(cellI >= 0) + { + partType = particleCloud_.particleType(index); + // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0 + partKs_[index][0] = typeKs_[partType - 1]; + } + } + + partKsField_.primitiveFieldRef() = 0.0; + particleCloud_.averagingM().resetWeightFields(); + particleCloud_.averagingM().setScalarAverage + ( + partKsField_, + partKs_, + particleCloud_.particleWeights(), + particleCloud_.averagingM().UsWeightField(), + NULL + ); +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H new file mode 100644 index 00000000..b3cafd7c --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H @@ -0,0 +1,116 @@ +/*---------------------------------------------------------------------------*\ +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 + + +Description + thermal conductivity of PARTICLE phase in presence of a fluid according to + Zehner and Schluender (1970) + +Class + ZehnerSchluenderThermCond + +SourceFiles + ZehnerSchluenderThermCond.C + +\*---------------------------------------------------------------------------*/ + + +#ifndef ZehnerSchluenderThermCond_H +#define ZehnerSchluenderThermCond_H + +#include "thermCondModel.H" +#include "scalarList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class ZehnerSchluenderThermCond Declaration +\*---------------------------------------------------------------------------*/ + +class ZehnerSchluenderThermCond +: + public thermCondModel +{ + +private: + + dictionary propsDict_; + + mutable volScalarField partKsField_; + + word voidfractionFieldName_; + + const volScalarField& voidfraction_; + + scalarList typeKs_; + + mutable double **partKs_; + + word wallQFactorName_; + + // ratio of half-cell-size and near-wall film + mutable volScalarField wallQFactor_; + + bool hasWallQFactor_; + + void allocateMyArrays() const; + + void calcPartKsField() const; + +public: + + //- Runtime type information + TypeName("ZehnerSchluenderThermCond"); + + + // Constructors + + //- Construct from components + ZehnerSchluenderThermCond + ( + const dictionary& dict, + cfdemCloud& sm + ); + + // Destructor + + ~ZehnerSchluenderThermCond(); + + + // Member Functions + + tmp thermCond() const; + + tmp thermDiff() const; + + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C index 3964b074..3db3e053 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/GaussVoidFraction/GaussVoidFraction.C @@ -95,7 +95,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract scalar radius(-1); scalar volume(0); - scalar scaleVol= weight(); + scalar scaleVol = weight(); scalar scaleRadius = pow(porosity(),1./3.); for(int index=0; index< particleCloud_.numberOfParticles(); index++) @@ -115,6 +115,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract label particleCenterCellID=particleCloud_.cellIDs()[index][0]; radius = particleCloud_.radius(index); + if (multiWeights_) scaleVol = weight(index); volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol; radius *= scaleRadius; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C index d176668f..e7dc28bc 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/bigParticleVoidFraction/bigParticleVoidFraction.C @@ -94,7 +94,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi scalar radius(-1); scalar volume(0); - scalar scaleVol= weight(); + scalar scaleVol = weight(); scalar scaleRadius = pow(porosity(),1./3.); for(int index=0; index< particleCloud_.numberOfParticles(); index++) @@ -113,6 +113,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi //collecting data label particleCenterCellID=particleCloud_.cellIDs()[index][0]; radius = particleCloud_.radius(index); + if (multiWeights_) scaleVol = weight(index); volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol; radius *= scaleRadius; vector positionCenter=particleCloud_.position(index); diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C index 0f1fb2bf..6aa96bce 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/centreVoidFraction/centreVoidFraction.C @@ -85,7 +85,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac scalar radius(-1); scalar volume(0); scalar cellVol(0); - scalar scaleVol= weight(); + scalar scaleVol = weight(); for(int index=0; index< particleCloud_.numberOfParticles(); index++) { @@ -99,6 +99,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac if (cellI >= 0) // particel centre is in domain { + if (multiWeights_) scaleVol = weight(index); cellVol = voidfractionNext_.mesh().V()[cellI]; radius = particleCloud_.radius(index); volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C index 3066be22..6648d3c7 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/dividedVoidFraction/dividedVoidFraction.C @@ -200,6 +200,7 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra position = particleCloud_.position(index); cellID = particleCloud_.cellIDs()[index][0]; radius = particleCloud_.radius(index); + if (multiWeights_) scaleVol = weight(index); volume = Vp(index,radius,scaleVol); radius *= scaleRadius; cellVol = 0; diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C index 5e395e31..8413532a 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C @@ -94,7 +94,7 @@ void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidf scalar radius(-1.); scalar volume(0.); - const scalar scaleVol = weight(); + scalar scaleVol = weight(); vector partPos(0.,0.,0.); vector pt(0.,0.,0.); @@ -140,6 +140,7 @@ void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidf if (cellI >= 0) // particel centre is in domain { radius = particleCloud_.radius(index); + if (multiWeights_) scaleVol = weight(index); volume = constant::mathematical::fourPiByThree * radius * radius * radius * scaleVol; // store volume for each particle diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C index d711590c..a2587c36 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C @@ -57,6 +57,7 @@ voidFractionModel::voidFractionModel : dict_(dict), particleCloud_(sm), + multiWeights_(false), voidfractionPrev_ ( IOobject ( @@ -89,6 +90,7 @@ voidFractionModel::voidFractionModel porosity_(1.) { particleCloud_.dataExchangeM().allocateArray(cellsPerParticle_,1,1); + if (particleCloud_.getParticleEffVolFactors()) multiWeights_ = true; } diff --git a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H index 4bdbe4c2..beedb6e2 100644 --- a/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H +++ b/src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H @@ -63,6 +63,8 @@ protected: cfdemCloud& particleCloud_; + bool multiWeights_; + mutable volScalarField voidfractionPrev_; mutable volScalarField voidfractionNext_; @@ -140,6 +142,11 @@ public: inline scalar weight()const { return weight_; } + inline scalar weight(label index) const + { + return particleCloud_.particleEffVolFactor(index); + } + inline scalar porosity()const { return porosity_; } inline void checkWeightNporosity(dictionary& propsDict) const diff --git a/src/lagrangian/cfdemParticleComp/Make/files b/src/lagrangian/cfdemParticleComp/Make/files index 58f28aa1..eb2a0e9d 100644 --- a/src/lagrangian/cfdemParticleComp/Make/files +++ b/src/lagrangian/cfdemParticleComp/Make/files @@ -31,12 +31,13 @@ $(cfdTools)/newGlobal.C $(energyModels)/energyModel/energyModel.C $(energyModels)/energyModel/newEnergyModel.C $(energyModels)/heatTransferGunn/heatTransferGunn.C -$(energyModels)/heatTransferGunnImplicit/heatTransferGunnImplicit.C +$(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C $(energyModels)/reactionHeat/reactionHeat.C $(thermCondModels)/thermCondModel/thermCondModel.C $(thermCondModels)/thermCondModel/newThermCondModel.C $(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C +$(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C $(thermCondModels)/noTherm/noThermCond.C $(forceModels)/forceModel/forceModel.C @@ -64,6 +65,7 @@ $(forceModels)/isotropicFluctuations/isotropicFluctuations.C $(forceModels)/evaluateFluctuations/evaluateFluctuations.C $(forceModels)/directedDiffusiveRelaxation/directedDiffusiveRelaxation.C $(forceModels)/BeetstraDrag/BeetstraDrag.C +$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C $(forceModels)/dSauter/dSauter.C $(forceModels)/Fines/Fines.C $(forceModels)/Fines/FinesFields.C diff --git a/tutorials/.gitignore b/tutorials/.gitignore index 6154fdee..a41c02d3 100644 --- a/tutorials/.gitignore +++ b/tutorials/.gitignore @@ -6,3 +6,16 @@ log_* log.* *~ +**/DEM/post/* +**/CFD/processor*/ +**/CFD/0*/ +**/CFD/1*/ +**/CFD/2*/ +**/CFD/3*/ +**/CFD/4*/ +**/CFD/5*/ +**/CFD/6*/ +**/CFD/7*/ +**/CFD/8*/ +**/CFD/9*/ +!**/CFD/0/ diff --git a/tutorials/cfdemSolverIB/twoSpheresGlowinskiMPI/parCFDDEMrun.sh b/tutorials/cfdemSolverIB/twoSpheresGlowinskiMPI/parCFDDEMrun.sh old mode 100644 new mode 100755 index 967b517e..f496b3f1 --- a/tutorials/cfdemSolverIB/twoSpheresGlowinskiMPI/parCFDDEMrun.sh +++ b/tutorials/cfdemSolverIB/twoSpheresGlowinskiMPI/parCFDDEMrun.sh @@ -22,6 +22,7 @@ solverName="cfdemSolverIB" nrProcs="4" machineFileName="none" # yourMachinefileName | none debugMode="off" # on | off| strict +reconstructCase="false" # true | false testHarnessPath="$CFDEM_TEST_HARNESS_PATH" runOctave="true" postproc="false" @@ -30,6 +31,14 @@ postproc="false" #- call function to run a parallel CFD-DEM case parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode +#- case needs special reconstruction +if [ $reconstructCase == "true" ] + then + cd $casePath/CFD + reconstructParMesh -mergeTol 1e-06 + reconstructPar -noLagrangian +fi + if [ $runOctave == "true" ] then @@ -60,9 +69,6 @@ if [ $postproc == "true" ] read fi -#- copy log file to test harness -cp ../../$logfileName $testHarnessPath - #- clean up case echo "deleting data at: $casePath" source $WM_PROJECT_DIR/bin/tools/CleanFunctions diff --git a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties index fd901353..1216be50 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties @@ -33,7 +33,7 @@ couplingInterval 100; voidFractionModel divided;//centre;// -locateModel engine;//turboEngine;// +locateModel engine;//turboEngineM2M;// meshMotionModel noMeshMotion; @@ -49,7 +49,7 @@ averagingModel dense;//dilute;// clockModel off;//standardClock;// -smoothingModel off;// constDiffSmoothing; // +smoothingModel off;// localPSizeDiffSmoothing;// constDiffSmoothing; // forceModels ( @@ -61,6 +61,7 @@ forceModels GidaspowDrag //Archimedes //volWeightedAverage + //totalMomentumExchange particleCellVolume ); @@ -74,6 +75,14 @@ turbulenceModelType "turbulenceProperties"; //===========================================================================// // sub-model properties +localPSizeDiffSmoothingProps +{ + lowerLimit 0.1; + upperLimit 1e10; + dSmoothingLength 1.5e-3; + Csmoothing 1.0; +} + constDiffSmoothingProps { lowerLimit 0.1; @@ -123,7 +132,13 @@ volWeightedAverageProps lowerThreshold 0; verbose true; } - +totalMomentumExchangeProps +{ + implicitMomExFieldName "Ksl"; + explicitMomExFieldName "none"; + fluidVelFieldName "U"; + granVelFieldName "Us"; +} GidaspowDragProps { verbose true; @@ -147,12 +162,15 @@ KochHillDragProps BeetstraDragProps { velFieldName "U"; - gravityFieldName "g"; - rhoParticle 2000.; voidfractionFieldName "voidfraction"; + granVelFieldName "Us"; interpolation true; - useFilteredDragModel ; - useParcelSizeDependentFilteredDrag ; +// useFilteredDragModel; +// useParcelSizeDependentFilteredDrag; + g 9.81; + rhoP 7000.; + rho 10.; + nuf 1.5e-4; k 0.05; aLimit 0.0; // verbose true; diff --git a/tutorials/cfdemSolverPiso/ErgunTestMPI_restart/parCFDDEMrun.sh b/tutorials/cfdemSolverPiso/ErgunTestMPI_restart/parCFDDEMrun.sh index 2dd14810..691b9397 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestMPI_restart/parCFDDEMrun.sh +++ b/tutorials/cfdemSolverPiso/ErgunTestMPI_restart/parCFDDEMrun.sh @@ -22,10 +22,10 @@ solverName="cfdemSolverPiso" nrProcs="4" machineFileName="none" # yourMachinefileName | none debugMode="off" # on | off| strict -reconstuctCase="true" # true | false +reconstructCase="true" # true | false testHarnessPath="$CFDEM_TEST_HARNESS_PATH" #--------------------------------------------------------------------------------# #- call function to run a parallel CFD-DEM case -parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode $reconstuctCase +parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode $reconstructCase diff --git a/tutorials/cfdemSolverPiso/settlingTestMPI/parCFDDEMrun.sh b/tutorials/cfdemSolverPiso/settlingTestMPI/parCFDDEMrun.sh index 3e8fcd20..f0d2ca8c 100644 --- a/tutorials/cfdemSolverPiso/settlingTestMPI/parCFDDEMrun.sh +++ b/tutorials/cfdemSolverPiso/settlingTestMPI/parCFDDEMrun.sh @@ -22,6 +22,7 @@ solverName="cfdemSolverPiso" nrProcs="2" machineFileName="none" # yourMachinefileName | none debugMode="off" # on | off| strict +reconstructCase="true" # true | false testHarnessPath="$CFDEM_TEST_HARNESS_PATH" runOctave="true" cleanUp="true" @@ -29,7 +30,7 @@ postproc="false" #--------------------------------------------------------------------------------# #- call function to run a parallel CFD-DEM case -parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode "true" +parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode $reconstructCase if [ $runOctave == "true" ] then