From 6a636658aa58ad5d99a86e3fb6ba48ac606492d0 Mon Sep 17 00:00:00 2001 From: Paul Kieckhefen Date: Tue, 16 May 2017 13:55:12 +0200 Subject: [PATCH 01/16] added grid coarsening/particle coarsening corrections to Beetstra drag model --- .../forceModel/BeetstraDrag/BeetstraDrag.C | 209 +++++++++++++++++- .../forceModel/BeetstraDrag/BeetstraDrag.H | 29 +++ .../CFD/constant/couplingProperties | 11 +- 3 files changed, 236 insertions(+), 13 deletions(-) diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C index 779a7a42..268a58cf 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C @@ -56,7 +56,16 @@ BeetstraDrag::BeetstraDrag UsFieldName_(propsDict_.lookup("granVelFieldName")), UsField_(sm.mesh().lookupObject (UsFieldName_)), scaleDia_(1.), - scaleDrag_(1.) + 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, "BeetstraDrag.logDat"); @@ -83,6 +92,30 @@ BeetstraDrag::BeetstraDrag if (propsDict_.found("scaleDrag")) scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag"))); + 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; + } + } + + + } @@ -122,6 +155,8 @@ void BeetstraDrag::setForce() const scalar magUr(0); scalar Rep(0); scalar localPhiP(0); + scalar GCcorr(1.); + scalar PCcorr(1.); vector dragExplicit(0,0,0); scalar dragCoefficient(0); @@ -163,7 +198,7 @@ void BeetstraDrag::setForce() const Ur = Ufluid-Us; magUr = mag(Ur); ds = 2*particleCloud_.radius(index); - ds_scaled = ds/scaleDia_; + ds_scaled = ds/scaleDia_; rho = rhoField[cellI]; nuf = nufField[cellI]; @@ -172,14 +207,29 @@ void BeetstraDrag::setForce() const // 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=ds_scaled*voidfraction*magUr/(nuf+SMALL); + dragCoefficient = F(Rep, voidfraction) + *3*M_PI*ds_scaled*nuf*rho*voidfraction + *scaleDia3*scaleDrag_; + + // calculate filtering corrections + if (useGC_) + { + GCcorr = 1.-h(localPhiP) + /( a(localPhiP) + *Lc2_ + /Foam::pow(U_.mesh().V()[cellI],.33333333) + + 1. + ); + if (usePC_) + { + PCcorr = Foam::exp(k_*(1.-scaleDia_)); + } + } + + // 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; @@ -200,6 +250,8 @@ void BeetstraDrag::setForce() const Pout << "nuf = " << nuf << endl; Pout << "voidfraction = " << voidfraction << endl; Pout << "Rep = " << Rep << endl; + Pout << "GCcorr = " << GCcorr << endl; + Pout << "PCcorr = " << PCcorr << endl; Pout << "drag = " << drag << endl; } @@ -220,6 +272,145 @@ 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; + } + return a0m + a1m*(phiP-phipam) + a2m*pow(phiP-phipam,2.) + a3m*pow(phiP-phipam,3.); +} + +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; + } + return h0m + h1m*(phiP-phiphm) + h2m*pow(phiP-phiphm,2) + h3m*pow(phiP-phiphm,3); +} + +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; +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H index a36ab894..ec59213f 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H @@ -59,6 +59,35 @@ private: 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_; + +protected: + 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/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties index 862de403..1216be50 100644 --- a/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties +++ b/tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties @@ -162,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; From 771dd7e6d05ac13da35e7751184b6ac88acf5f65 Mon Sep 17 00:00:00 2001 From: tlichtenegger Date: Thu, 12 Apr 2018 11:06:46 +0200 Subject: [PATCH 02/16] Initial commit for cfdemSolverRhoSimple. New heat transfer model to treat granular temperature as field to be solved on a grid. --- .../solvers/cfdemSolverRhoPimple/EEqn.H | 4 - .../cfdemSolverRhoPimple.C | 4 + .../solvers/cfdemSolverRhoSimple/EEqn.H | 57 +++++ .../solvers/cfdemSolverRhoSimple/Make/files | 3 + .../solvers/cfdemSolverRhoSimple/Make/options | 32 +++ .../solvers/cfdemSolverRhoSimple/UEqn.H | 32 +++ .../cfdemSolverRhoSimple.C | 141 ++++++++++ .../cfdemSolverRhoSimple/createFieldRefs.H | 2 + .../cfdemSolverRhoSimple/createFields.H | 242 ++++++++++++++++++ .../solvers/cfdemSolverRhoSimple/pEqn.H | 80 ++++++ etc/solver-list.txt | 1 + src/lagrangian/cfdemParticle/Make/files | 3 +- .../cfdemCloudEnergy/cfdemCloudEnergy.C | 6 + .../cfdemCloudEnergy/cfdemCloudEnergy.H | 2 + .../energyModel/energyModel/energyModel.H | 2 + .../heatTransferGunn/heatTransferGunn.C | 184 +++++++++++-- .../heatTransferGunn/heatTransferGunn.H | 74 ++++-- .../heatTransferGunnPartField.C | 149 +++++++++++ .../heatTransferGunnPartField.H | 93 +++++++ .../ZehnerSchluenderThermCond.C | 149 +++++++++++ .../ZehnerSchluenderThermCond.H | 107 ++++++++ src/lagrangian/cfdemParticleComp/Make/files | 3 +- 22 files changed, 1308 insertions(+), 62 deletions(-) create mode 100644 applications/solvers/cfdemSolverRhoSimple/EEqn.H create mode 100644 applications/solvers/cfdemSolverRhoSimple/Make/files create mode 100644 applications/solvers/cfdemSolverRhoSimple/Make/options create mode 100644 applications/solvers/cfdemSolverRhoSimple/UEqn.H create mode 100644 applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C create mode 100644 applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H create mode 100644 applications/solvers/cfdemSolverRhoSimple/createFields.H create mode 100644 applications/solvers/cfdemSolverRhoSimple/pEqn.H create mode 100644 src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C create mode 100644 src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H create mode 100644 src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C create mode 100644 src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H index 606d9fe2..fc770162 100644 --- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H +++ b/applications/solvers/cfdemSolverRhoPimple/EEqn.H @@ -52,8 +52,4 @@ thermo.correct(); Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl; - - particleCloud.clockM().start(31,"postFlow"); - particleCloud.postFlow(); - particleCloud.clockM().stop("postFlow"); } 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..c4a58ad2 --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/UEqn.H @@ -0,0 +1,32 @@ +// 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); + + fvOptions.correct(U); + K = 0.5*magSqr(U); +} +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..1b020f9f --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C @@ -0,0 +1,141 @@ +/*---------------------------------------------------------------------------*\ +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 PIMPLE corrector loop + { + #include "UEqn.H" + #include "EEqn.H" + + // besides this pEqn, OF offers a "simple consistent"-option + #include "pEqn.H" + rhoeps=rho*voidfraction; + + turbulence->correct(); + } + + particleCloud.clockM().start(32,"postFlow"); + 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..a85e7f44 --- /dev/null +++ b/applications/solvers/cfdemSolverRhoSimple/pEqn.H @@ -0,0 +1,80 @@ +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); diff --git a/etc/solver-list.txt b/etc/solver-list.txt index 9d9911fa..60c0eb39 100644 --- a/etc/solver-list.txt +++ b/etc/solver-list.txt @@ -1,5 +1,6 @@ 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 7cc1c491..b1b8e9e8 100644 --- a/src/lagrangian/cfdemParticle/Make/files +++ b/src/lagrangian/cfdemParticle/Make/files @@ -32,12 +32,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 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("expNusselt",false)), interpolation_(propsDict_.lookupOrDefault("interpolation",false)), verbose_(propsDict_.lookupOrDefault("verbose",false)), + implicit_(propsDict_.lookupOrDefault("implicit",true)), QPartFluidName_(propsDict_.lookupOrDefault("QPartFluidName","QPartFluid")), QPartFluid_ ( IOobject @@ -58,13 +60,26 @@ heatTransferGunn::heatTransferGunn sm.mesh(), dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0) ), + 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", sm.mesh().time().timeName(), sm.mesh(), - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, IOobject::NO_WRITE ), sm.mesh(), @@ -108,6 +123,8 @@ heatTransferGunn::heatTransferGunn ), partRefTemp_("partRefTemp", dimensionSet(0,0,0,1,0,0,0), 0.0), calcPartTempField_(propsDict_.lookupOrDefault("calcPartTempField",false)), + calcPartTempAve_(propsDict_.lookupOrDefault("calcPartTempAve",false)), + partTempAve_(0.0), tempFieldName_(propsDict_.lookupOrDefault("tempFieldName","T")), tempField_(sm.mesh().lookupObject (tempFieldName_)), voidfractionFieldName_(propsDict_.lookupOrDefault("voidfractionFieldName","voidfraction")), @@ -121,32 +138,47 @@ heatTransferGunn::heatTransferGunn partTemp_(NULL), partHeatFluxName_(propsDict_.lookup("partHeatFluxName")), partHeatFlux_(NULL), + partHeatFluxCoeff_(NULL), partRe_(NULL), partNu_(NULL) { - allocateMyArrays(); + allocateMyArrays(); if (propsDict_.found("maxSource")) { maxSource_=readScalar(propsDict_.lookup ("maxSource")); Info << "limiting eulerian source field to: " << maxSource_ << endl; } + if (calcPartTempField_) { + calcPartTempAve_ = true; if (propsDict_.found("partRefTemp")) + { partRefTemp_.value()=readScalar(propsDict_.lookup ("partRefTemp")); + } partTempField_.writeOpt() = IOobject::AUTO_WRITE; partRelTempField_.writeOpt() = IOobject::AUTO_WRITE; partTempField_.write(); partRelTempField_.write(); Info << "Particle temperature field activated." << endl; } + + if (!implicit_) + { + QPartFluidCoeff_.writeOpt() = IOobject::NO_WRITE; + } + if (verbose_) { ReField_.writeOpt() = IOobject::AUTO_WRITE; NuField_.writeOpt() = IOobject::AUTO_WRITE; ReField_.write(); NuField_.write(); + if (expNusselt_) + { + FatalError <<"Cannot read and create NuField at the same time!\n" << abort(FatalError); + } } } @@ -159,6 +191,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 * * * * * * * * * * * * * // @@ -168,6 +204,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_) { @@ -227,6 +267,7 @@ void heatTransferGunn::calcEnergyContribution() scalar Rep(0); scalar Pr(0); scalar Nup(0); + scalar Tsum(0.0); interpolationCellPoint voidfractionInterpolator_(voidfraction_); @@ -265,13 +306,12 @@ void heatTransferGunn::calcEnergyContribution() Nup = Nusselt(voidfraction, Rep, Pr); - + Tsum += partTemp_[index][0]; scalar h = kf0_ * Nup / ds; scalar As = ds * ds * M_PI; // surface area of sphere // calc convective heat flux [W] heatFlux(index, h, As, Tfluid); - heatFluxCoeff(index, h, As); if(verbose_) { @@ -295,6 +335,16 @@ void heatTransferGunn::calcEnergyContribution() } } + // 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_, @@ -305,6 +355,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; @@ -329,22 +394,22 @@ 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 @@ -352,6 +417,14 @@ void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const Qsource += QPartFluid_; } +void heatTransferGunn::addEnergyCoefficient(volScalarField& Qsource) const +{ + if(implicit_) + { + Qsource += QPartFluidCoeff_; + } +} + scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) const { return (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) * @@ -362,22 +435,77 @@ scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) con 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) + scalar hAs = h * As; + partHeatFlux_[index][0] = - hAs * partTemp_[index][0]; + if(!implicit_) { - Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl; - - particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_); + partHeatFlux_[index][0] += hAs * Tfluid; } + else + { + partHeatFluxCoeff_[index][0] = hAs; + } +} + +void heatTransferGunn::giveData() +{ + Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl; + + particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_); +} + +void heatTransferGunn::postFlow() +{ + if(implicit_) + { + label cellI; + scalar Tfluid(0.0); + scalar Tpart(0.0); + interpolationCellPoint TInterpolator_(tempField_); + + 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 + ); + + dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), partTempAve_); + partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_); } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H index 4862e5ed..aa6facf7 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,27 +44,39 @@ class heatTransferGunn protected: dictionary propsDict_; - + + bool expNusselt_; + bool interpolation_; - + bool verbose_; - + + bool implicit_; + word QPartFluidName_; - + volScalarField QPartFluid_; - - volScalarField partTempField_; - + + word QPartFluidCoeffName_; + + volScalarField QPartFluidCoeff_; + + mutable volScalarField partTempField_; + volScalarField partRelTempField_; - + volScalarField ReField_; - + volScalarField NuField_; - + dimensionedScalar partRefTemp_; - + bool calcPartTempField_; + bool calcPartTempAve_; + + scalar partTempAve_; + word tempFieldName_; const volScalarField& tempField_; // ref to temperature field @@ -78,33 +90,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 @@ -128,10 +142,14 @@ 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..9ac88758 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C @@ -0,0 +1,149 @@ +/*---------------------------------------------------------------------------*\ +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), + fvOptions(fv::options::New(sm.mesh())), + thermCondModel_ + ( + thermCondModel::New + ( + propsDict_, + sm + ) + ) +{ + if(!implicit_) + { + FatalError << "heatTransferGunnPartField requires implicit heat transfer treatment." << abort(FatalError); + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +heatTransferGunnPartField::~heatTransferGunnPartField() +{ +} + +// * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * * * Member Fct * * * * * * * * * * * * * * * // +void heatTransferGunnPartField::calcEnergyContribution() +{ + heatTransferGunn::calcEnergyContribution(); + // if heat sources in particles present, pull them here +} + +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_); + + 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() +{ +//use gasT and QCoeff to solve particle temp eqn + volScalarField Qsource = QPartFluidCoeff_*tempField_; + volScalarField alphaP = 1.0 - voidfraction_; + volScalarField thCond = thermCondModel_().thermCond(); + + fvScalarMatrix partTEqn + ( + Qsource + - fvm::Sp(QPartFluidCoeff_, partTempField_) + - fvm::laplacian(alphaP*thCond,partTempField_) + == + fvOptions(rho_, 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(); + + fvOptions.correct(partTempField_); + + 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..12ae0d82 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.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 + 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" +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +class thermCondModel; +/*---------------------------------------------------------------------------*\ + Class heatTransferGunnPartField Declaration +\*---------------------------------------------------------------------------*/ + +class heatTransferGunnPartField +: + public heatTransferGunn +{ +private: + + fv::options& fvOptions; + + autoPtr thermCondModel_; + + 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/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C new file mode 100644 index 00000000..b0241cf2 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C @@ -0,0 +1,149 @@ +/*---------------------------------------------------------------------------*\ +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")), + voidfractionFieldName_(propsDict_.lookupOrDefault("voidfractionFieldName","voidfraction")), + voidfraction_(sm.mesh().lookupObject (voidfractionFieldName_)), + ks0_(transportProperties_.lookup("ks")), + 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 (wallQFactor_.headerOk()) + { + Info << "Found field for scaling wall heat flux.\n" << endl; + hasWallQFactor_ = true; + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +ZehnerSchluenderThermCond::~ZehnerSchluenderThermCond() +{} + + +// * * * * * * * * * * * * * * * 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) + ) + ); + + volScalarField& svf = tvf.ref(); + + scalar A = ks0_.value()/kf0_.value(); + scalar B = 0.0; + scalar C = 0.0; + scalar k = 0.0; + scalar OnemBoA = 0.0; + scalar voidfraction = 0.0; + scalar w = 7.26e-3; + + forAll(svf, cellI) + { + voidfraction = voidfraction_[cellI]; + B = 1.25 * Foam::pow((1 - voidfraction) / voidfraction, 1.11); + OnemBoA = 1.0 - B/A; + C = (A - 1) / (OnemBoA * OnemBoA) * B/A * log(A/B) - (B - 1)/OnemBoA - 0.5 * (B + 1); + C *= 2.0 / OnemBoA; + k = Foam::sqrt(1 - voidfraction) * (w * A + (1 - w) * C) * kf0_.value(); + svf[cellI] = k / (1 - voidfraction); + } + + // 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); +} + + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // 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..12320d75 --- /dev/null +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H @@ -0,0 +1,107 @@ +/*---------------------------------------------------------------------------*\ +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" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class ZehnerSchluenderThermCond Declaration +\*---------------------------------------------------------------------------*/ + +class ZehnerSchluenderThermCond +: + public thermCondModel +{ + +private: + + dictionary propsDict_; + + word voidfractionFieldName_; + + const volScalarField& voidfraction_; + + dimensionedScalar ks0_; + + word wallQFactorName_; + + // ratio of half-cell-size and near-wall film + mutable volScalarField wallQFactor_; + + bool hasWallQFactor_; + +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/cfdemParticleComp/Make/files b/src/lagrangian/cfdemParticleComp/Make/files index b119a3b3..0613c1cf 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 From 76ea1af1abb02951aed8671352084a431a650089 Mon Sep 17 00:00:00 2001 From: tlichtenegger Date: Fri, 13 Apr 2018 11:19:08 +0200 Subject: [PATCH 03/16] CFDEMcloud can now hold particle types and densities. Heat transfer model needs to know heat capacities to solve T eqn. --- .../cfdemSolverRhoSimple.C | 17 +++--- .../cfdemParticle/cfdemCloud/cfdemCloud.C | 46 ++++++++++++++- .../cfdemParticle/cfdemCloud/cfdemCloud.H | 19 +++++- .../cfdemParticle/cfdemCloud/cfdemCloudI.H | 36 ++++++++++++ .../heatTransferGunnPartField.C | 58 ++++++++++++++++--- .../heatTransferGunnPartField.H | 10 +++- .../ZehnerSchluenderThermCond.C | 25 +++++--- 7 files changed, 181 insertions(+), 30 deletions(-) diff --git a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C index 1b020f9f..e960b40b 100644 --- a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C +++ b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C @@ -105,17 +105,16 @@ int main(int argc, char *argv[]) particleCloud.clockM().start(26,"Flow"); volScalarField rhoeps("rhoeps",rho*voidfraction); - // --- Pressure-velocity PIMPLE corrector loop - { - #include "UEqn.H" - #include "EEqn.H" + // Pressure-velocity SIMPLE corrector - // besides this pEqn, OF offers a "simple consistent"-option - #include "pEqn.H" - rhoeps=rho*voidfraction; + #include "UEqn.H" + #include "EEqn.H" - turbulence->correct(); - } + // besides this pEqn, OF offers a "simple consistent"-option + #include "pEqn.H" + rhoeps=rho*voidfraction; + + turbulence->correct(); particleCloud.clockM().start(32,"postFlow"); particleCloud.postFlow(); diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C index 89e2cead..5d132ddd 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C @@ -83,6 +83,8 @@ cfdemCloud::cfdemCloud ignore_(false), allowCFDsubTimestep_(true), limitDEMForces_(false), + getParticleDensities_(couplingProperties_.lookupOrDefault("getParticleDensities",false)), + getParticleTypes_(couplingProperties_.lookupOrDefault("getParticleTypes",false)), modelType_(couplingProperties_.lookup("modelType")), positions_(NULL), velocities_(NULL), @@ -127,6 +129,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_ ( @@ -358,6 +373,8 @@ cfdemCloud::~cfdemCloud() dataExchangeM().destroy(particleWeights_,1); dataExchangeM().destroy(particleVolumes_,1); dataExchangeM().destroy(particleV_,1); + if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1); + if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1); } // * * * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * * // @@ -369,6 +386,9 @@ 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(getParticleTypes_) dataExchangeM().getData("type","scalar-atom",particleTypes_); } void cfdemCloud::giveDEMdata() @@ -446,6 +466,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; @@ -598,7 +637,8 @@ bool cfdemCloud::evolve clockM().stop("setvoidFraction"); // set average particles velocity field - clockM().start(20,"setVectorAverage"); + clockM().start(20,"setAverages"); + setScalarAverages(); setVectorAverages(); @@ -612,7 +652,7 @@ bool cfdemCloud::evolve if(!treatVoidCellsAsExplicitForce()) smoothingM().smoothenReferenceField(averagingM().UsNext()); - clockM().stop("setVectorAverage"); + clockM().stop("setAverages"); } //============================================ @@ -703,6 +743,8 @@ 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(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 0f74cc1f..1b3f78cf 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H @@ -96,6 +96,10 @@ protected: bool limitDEMForces_; + bool getParticleDensities_; + + bool getParticleTypes_; + scalar maxDEMForce_; const word modelType_; @@ -122,6 +126,10 @@ protected: mutable int **cellIDs_; + mutable double **particleDensities_; + + mutable int **particleTypes_; + mutable double **particleWeights_; mutable double **particleVolumes_; @@ -162,6 +170,8 @@ protected: mutable volScalarField ddtVoidfraction_; + mutable volScalarField particleDensityField_; + mutable Switch checkPeriodicCells_; const turbulenceModel& turbulence_; @@ -205,6 +215,8 @@ protected: virtual void setParticleForceField(); + virtual void setScalarAverages(); + virtual void setVectorAverages(); public: @@ -323,9 +335,10 @@ 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;} + virtual inline double ** particleDensity() const; + virtual inline double particleDensity(label index) 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;} diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H index 2fbe1744..56ddfa48 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H @@ -216,6 +216,42 @@ inline double cfdemCloud::d32(bool recalc) return d32_; } +inline double ** cfdemCloud::particleDensity() const +{ + if(!getParticleDensities_) return NULL; + else + { + return particleDensities_; + } +} + +inline double cfdemCloud::particleDensity(label index) const +{ + if(!getParticleDensities_) return -1.0; + else + { + return particleDensities_[index][0]; + } +} + +inline int ** cfdemCloud::particleTypes() const +{ + if(!getParticleTypes_) return NULL; + else + { + 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/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C index 9ac88758..6e354810 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C @@ -44,7 +44,22 @@ heatTransferGunnPartField::heatTransferGunnPartField ) : heatTransferGunn(dict,sm), - fvOptions(fv::options::New(sm.mesh())), + 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")), + partCp_(NULL), thermCondModel_ ( thermCondModel::New @@ -52,12 +67,15 @@ heatTransferGunnPartField::heatTransferGunnPartField propsDict_, sm ) - ) + ), + fvOptions(fv::options::New(sm.mesh())) { if(!implicit_) { FatalError << "heatTransferGunnPartField requires implicit heat transfer treatment." << abort(FatalError); } + + allocateMyArrays(); } @@ -65,14 +83,36 @@ heatTransferGunnPartField::heatTransferGunnPartField 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 + partCpField_.primitiveFieldRef() = 0.0; + particleCloud_.averagingM().resetWeightFields(); + particleCloud_.averagingM().setScalarAverage + ( + partCpField_, + partCp_, + particleCloud_.particleWeights(), + particleCloud_.averagingM().UsWeightField(), + NULL + ); + } void heatTransferGunnPartField::addEnergyContribution(volScalarField& Qsource) const @@ -114,18 +154,20 @@ void heatTransferGunnPartField::postFlow() void heatTransferGunnPartField::solve() { -//use gasT and QCoeff to solve particle temp eqn - volScalarField Qsource = QPartFluidCoeff_*tempField_; + Info << "patch types of partTemp boundary: " << partTempField_.boundaryField().types() << endl; + volScalarField Qsource = QPartFluidCoeff_*tempField_/partCpField_; volScalarField alphaP = 1.0 - voidfraction_; - volScalarField thCond = thermCondModel_().thermCond(); + volScalarField partRhoEff = alphaP*partRhoField_; +// volScalarField thCond = thermCondModel_().thermCond(); + fvScalarMatrix partTEqn ( Qsource - - fvm::Sp(QPartFluidCoeff_, partTempField_) - - fvm::laplacian(alphaP*thCond,partTempField_) + - fvm::Sp(QPartFluidCoeff_/partCpField_, partTempField_) + // - fvm::laplacian(alphaP*thCond/partCpField_,partTempField_) == - fvOptions(rho_, partTempField_) + fvOptions(partRhoEff, partTempField_) ); // if transient add time derivative - need particle density and specific heat fields // if sources activated add sources diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H index 12ae0d82..a52cdec0 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H @@ -44,10 +44,18 @@ class heatTransferGunnPartField { private: - fv::options& fvOptions; + volScalarField partCpField_; + + const volScalarField& partRhoField_; + + mutable double **partCp_; autoPtr thermCondModel_; + fv::options& fvOptions; + + void allocateMyArrays() const; + void giveData(); public: diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C index b0241cf2..35f816dd 100644 --- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C @@ -99,7 +99,8 @@ tmp ZehnerSchluenderThermCond::thermCond() const IOobject::NO_WRITE ), voidfraction_.mesh(), - dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0) + dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0), + "zeroGradient" ) ); @@ -115,15 +116,25 @@ tmp ZehnerSchluenderThermCond::thermCond() const forAll(svf, cellI) { + // debugging + Pout << "calculating field in cell " << cellI << endl; voidfraction = voidfraction_[cellI]; - B = 1.25 * Foam::pow((1 - voidfraction) / voidfraction, 1.11); - OnemBoA = 1.0 - B/A; - C = (A - 1) / (OnemBoA * OnemBoA) * B/A * log(A/B) - (B - 1)/OnemBoA - 0.5 * (B + 1); - C *= 2.0 / OnemBoA; - k = Foam::sqrt(1 - voidfraction) * (w * A + (1 - w) * C) * kf0_.value(); - svf[cellI] = k / (1 - voidfraction); + if(voidfraction > 1.0 - SMALL) svf[cellI] = 0.0; + else + { + B = 1.25 * Foam::pow((1 - voidfraction) / voidfraction, 1.11); + OnemBoA = 1.0 - B/A; + C = (A - 1) / (OnemBoA * OnemBoA) * B/A * log(A/B) - (B - 1)/OnemBoA - 0.5 * (B + 1); + C *= 2.0 / OnemBoA; + k = Foam::sqrt(1 - voidfraction) * (w * A + (1 - w) * C) * kf0_.value(); + svf[cellI] = k / (1 - voidfraction); + } } + // debugging + Pout << "patch types of svf boundary: " << svf.boundaryField().types() << endl; + svf.correctBoundaryConditions(); + // if a wallQFactor field is present, use it to scale heat transport through a patch if (hasWallQFactor_) { From 24c5c25f7c83eab845f4abcb50a53286bb95fe52 Mon Sep 17 00:00:00 2001 From: tlichtenegger Date: Mon, 16 Apr 2018 11:44:48 +0200 Subject: [PATCH 04/16] Specific heat and thermal conductivities for various species. --- .../heatTransferGunn/heatTransferGunn.C | 5 +- .../heatTransferGunnPartField.C | 34 +++++++-- .../heatTransferGunnPartField.H | 3 + .../ZehnerSchluenderThermCond.C | 75 +++++++++++++++++-- .../ZehnerSchluenderThermCond.H | 11 ++- 5 files changed, 110 insertions(+), 18 deletions(-) diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C index cd945d30..47530705 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C @@ -76,14 +76,15 @@ heatTransferGunn::heatTransferGunn partTempField_ ( IOobject ( - "particleTemp", + "partTemp", sm.mesh().time().timeName(), sm.mesh(), IOobject::READ_IF_PRESENT, IOobject::NO_WRITE ), sm.mesh(), - dimensionedScalar("zero", dimensionSet(0,0,0,1,0,0,0), 0.0) + dimensionedScalar("zero", dimensionSet(0,0,0,1,0,0,0), 0.0), + "zeroGradient" ), partRelTempField_ ( IOobject diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C index 6e354810..17be0c17 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C @@ -59,6 +59,7 @@ heatTransferGunnPartField::heatTransferGunnPartField "zeroGradient" ), partRhoField_(sm.mesh().lookupObject("partRho")), + typeCp_(propsDict_.lookupOrDefault("specificHeatCapacities",scalarList(1,-1.0))), partCp_(NULL), thermCondModel_ ( @@ -70,11 +71,16 @@ heatTransferGunnPartField::heatTransferGunnPartField ), fvOptions(fv::options::New(sm.mesh())) { - if(!implicit_) + 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); + } + allocateMyArrays(); } @@ -102,6 +108,19 @@ void heatTransferGunnPartField::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 @@ -155,19 +174,20 @@ void heatTransferGunnPartField::postFlow() void heatTransferGunnPartField::solve() { Info << "patch types of partTemp boundary: " << partTempField_.boundaryField().types() << endl; - volScalarField Qsource = QPartFluidCoeff_*tempField_/partCpField_; + volScalarField Qsource = QPartFluidCoeff_*tempField_; volScalarField alphaP = 1.0 - voidfraction_; - volScalarField partRhoEff = alphaP*partRhoField_; -// volScalarField thCond = thermCondModel_().thermCond(); + volScalarField partCpEff = alphaP*partRhoField_*partCpField_; + volScalarField thCondEff = alphaP*thermCondModel_().thermCond(); +// thCondEff.correctBoundaryConditions(); fvScalarMatrix partTEqn ( Qsource - - fvm::Sp(QPartFluidCoeff_/partCpField_, partTempField_) - // - fvm::laplacian(alphaP*thCond/partCpField_,partTempField_) + - fvm::Sp(QPartFluidCoeff_, partTempField_) + - fvm::laplacian(thCondEff,partTempField_) == - fvOptions(partRhoEff, partTempField_) + fvOptions(partCpEff, partTempField_) ); // if transient add time derivative - need particle density and specific heat fields // if sources activated add sources diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H index a52cdec0..a321be7e 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H @@ -29,6 +29,7 @@ License #include "cfdemCloudEnergy.H" #include "heatTransferGunn.H" #include "fvOptions.H" +#include "scalarList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam @@ -48,6 +49,8 @@ private: const volScalarField& partRhoField_; + scalarList typeCp_; + mutable double **partCp_; autoPtr thermCondModel_; diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C index 35f816dd..7b8b9a95 100644 --- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C @@ -50,9 +50,24 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond : 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_)), - ks0_(transportProperties_.lookup("ks")), + typeKs_(propsDict_.lookupOrDefault("thermalConductivities",scalarList(1,-1.0))), + partKs_(NULL), wallQFactorName_(propsDict_.lookupOrDefault("wallQFactorName","wallQFactor")), wallQFactor_ ( IOobject @@ -68,6 +83,17 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond ), 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; @@ -79,9 +105,17 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond // * * * * * * * * * * * * * * * * 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 @@ -106,7 +140,8 @@ tmp ZehnerSchluenderThermCond::thermCond() const volScalarField& svf = tvf.ref(); - scalar A = ks0_.value()/kf0_.value(); + calcPartKsField(); + scalar A = 0.0; scalar B = 0.0; scalar C = 0.0; scalar k = 0.0; @@ -116,12 +151,11 @@ tmp ZehnerSchluenderThermCond::thermCond() const forAll(svf, cellI) { - // debugging - Pout << "calculating field in cell " << cellI << endl; voidfraction = voidfraction_[cellI]; - if(voidfraction > 1.0 - SMALL) svf[cellI] = 0.0; + 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); OnemBoA = 1.0 - B/A; C = (A - 1) / (OnemBoA * OnemBoA) * B/A * log(A/B) - (B - 1)/OnemBoA - 0.5 * (B + 1); @@ -131,8 +165,6 @@ tmp ZehnerSchluenderThermCond::thermCond() const } } - // debugging - Pout << "patch types of svf boundary: " << svf.boundaryField().types() << endl; svf.correctBoundaryConditions(); // if a wallQFactor field is present, use it to scale heat transport through a patch @@ -151,8 +183,35 @@ tmp ZehnerSchluenderThermCond::thermDiff() const FatalError << "ZehnerSchluenderThermCond does not provide thermal diffusivity." << abort(FatalError); } +void ZehnerSchluenderThermCond::calcPartKsField() const +{ + if (typeKs_.size() <= 1) return; + 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 index 12320d75..b3cafd7c 100644 --- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.H @@ -34,6 +34,7 @@ SourceFiles #define ZehnerSchluenderThermCond_H #include "thermCondModel.H" +#include "scalarList.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -52,12 +53,16 @@ class ZehnerSchluenderThermCond private: dictionary propsDict_; + + mutable volScalarField partKsField_; word voidfractionFieldName_; const volScalarField& voidfraction_; - dimensionedScalar ks0_; + scalarList typeKs_; + + mutable double **partKs_; word wallQFactorName_; @@ -66,6 +71,10 @@ private: bool hasWallQFactor_; + void allocateMyArrays() const; + + void calcPartKsField() const; + public: //- Runtime type information From d8682c6a79eb7823ab3775dc647c30363b5908be Mon Sep 17 00:00:00 2001 From: tlichtenegger Date: Thu, 19 Apr 2018 15:07:24 +0200 Subject: [PATCH 05/16] Some experiments to get more stability. --- .../solvers/cfdemSolverRhoPimple/EEqn.H | 7 ++++- .../solvers/cfdemSolverRhoSimple/UEqn.H | 19 +++++++++----- .../cfdemSolverRhoSimple.C | 4 ++- .../solvers/cfdemSolverRhoSimple/pEqn.H | 10 +++++++ .../heatTransferGunnPartField.C | 26 ++++++++++++++++--- 5 files changed, 55 insertions(+), 11 deletions(-) diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H index fc770162..0f1e6aff 100644 --- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H +++ b/applications/solvers/cfdemSolverRhoPimple/EEqn.H @@ -39,7 +39,7 @@ == fvOptions(rho, he) ); - + EEqn.relax(); @@ -52,4 +52,9 @@ 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/UEqn.H b/applications/solvers/cfdemSolverRhoSimple/UEqn.H index c4a58ad2..45b78abe 100644 --- a/applications/solvers/cfdemSolverRhoSimple/UEqn.H +++ b/applications/solvers/cfdemSolverRhoSimple/UEqn.H @@ -19,14 +19,21 @@ fvOptions.constrain(UEqn); if (modelType=="B" || modelType=="Bfull") { solve(UEqn == -fvc::grad(p)+ Ksl*Us); - - fvOptions.correct(U); - K = 0.5*magSqr(U); } else { solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us); - - fvOptions.correct(U); - K = 0.5*magSqr(U); } + +fvOptions.correct(U); + +// limit vel + +forAll(U,cellI) +{ + scalar mU = mag(U[cellI]); + if (mU > 200.0) U[cellI] *= 200.0/mU; +} + + +K = 0.5*magSqr(U); diff --git a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C index e960b40b..b505252e 100644 --- a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C +++ b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C @@ -108,12 +108,14 @@ int main(int argc, char *argv[]) // Pressure-velocity SIMPLE corrector #include "UEqn.H" - #include "EEqn.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"); diff --git a/applications/solvers/cfdemSolverRhoSimple/pEqn.H b/applications/solvers/cfdemSolverRhoSimple/pEqn.H index a85e7f44..c032cbeb 100644 --- a/applications/solvers/cfdemSolverRhoSimple/pEqn.H +++ b/applications/solvers/cfdemSolverRhoSimple/pEqn.H @@ -75,6 +75,16 @@ else { U = HbyA - rAU*(fvc::grad(p)-Ksl*Us); } + +// limit vel + +forAll(U,cellI) +{ + scalar mU = mag(U[cellI]); + if (mU > 200.0) U[cellI] *= 200.0/mU; +} + + U.correctBoundaryConditions(); fvOptions.correct(U); K = 0.5*magSqr(U); diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C index 17be0c17..2dac2222 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C @@ -81,6 +81,8 @@ heatTransferGunnPartField::heatTransferGunnPartField FatalError << "heatTransferGunnPartField: provide list of specific heat capacities." << abort(FatalError); } + partTempField_.writeOpt() = IOobject::AUTO_WRITE; + allocateMyArrays(); } @@ -173,18 +175,29 @@ void heatTransferGunnPartField::postFlow() void heatTransferGunnPartField::solve() { - Info << "patch types of partTemp boundary: " << partTempField_.boundaryField().types() << endl; - volScalarField Qsource = QPartFluidCoeff_*tempField_; + // 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(QPartFluidCoeff_, partTempField_) + - fvm::Sp(correctedQPartFluidCoeff, partTempField_) - fvm::laplacian(thCondEff,partTempField_) == fvOptions(partCpEff, partTempField_) @@ -199,8 +212,15 @@ void heatTransferGunnPartField::solve() partTEqn.solve(); + partTempField_.relax(); + fvOptions.correct(partTempField_); + dimensionedScalar pTMin("pTMin",dimensionSet(0,0,0,1,0,0,0),293.0); + dimensionedScalar pTMax("pTMin",dimensionSet(0,0,0,1,0,0,0),3000.0); + partTempField_ = max(partTempField_, pTMin); + partTempField_ = min(partTempField_, pTMax); + Info<< "partT max/min : " << max(partTempField_).value() << " " << min(partTempField_).value() << endl; } From bdfdb23cddc2fb2b9426a5bacf73eef5289179c4 Mon Sep 17 00:00:00 2001 From: tlichtenegger Date: Mon, 23 Apr 2018 14:15:34 +0200 Subject: [PATCH 06/16] Some debugging (initialization of arrays etc.). --- src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C | 6 +++++- src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H | 2 ++ src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H | 10 ++++++++++ .../ZehnerSchluenderThermCond.C | 5 +++++ 4 files changed, 22 insertions(+), 1 deletion(-) diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C index 5d132ddd..2eb91237 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C @@ -97,6 +97,8 @@ cfdemCloud::cfdemCloud radii_(NULL), voidfractions_(NULL), cellIDs_(NULL), + particleDensities_(NULL), + particleTypes_(NULL), particleWeights_(NULL), particleVolumes_(NULL), particleV_(NULL), @@ -743,8 +745,10 @@ bool cfdemCloud::reAllocArrays() dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle()); dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle()); dataExchangeM().allocateArray(particleV_,0.,1); + Info << "\nnow allocating particle densities array.\n" << endl; if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1); - if(getParticleTypes_) dataExchangeM().allocateArray(particleTypes_,0.,1); + Info << "\nnow allocating particle types array.\n" << endl; + 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 1b3f78cf..e736c974 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H @@ -335,8 +335,10 @@ public: virtual inline int maxType() {return -1;} virtual inline bool multipleTypesDMax() {return false;} virtual inline bool multipleTypesDMin() {return false;} + inline bool getParticleDensities() const; virtual inline double ** particleDensity() const; virtual inline double particleDensity(label index) const; + inline bool getParticleTypes() const; virtual inline int ** particleTypes() const; virtual inline label particleType(label index) const; diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H index 56ddfa48..0db389e0 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H @@ -216,6 +216,11 @@ inline double cfdemCloud::d32(bool recalc) return d32_; } +inline bool cfdemCloud::getParticleDensities() const +{ + return getParticleDensities_; +} + inline double ** cfdemCloud::particleDensity() const { if(!getParticleDensities_) return NULL; @@ -234,6 +239,11 @@ inline double cfdemCloud::particleDensity(label index) const } } +inline bool cfdemCloud::getParticleTypes() const +{ + return getParticleTypes_; +} + inline int ** cfdemCloud::particleTypes() const { if(!getParticleTypes_) return NULL; diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C index 7b8b9a95..1d1cec33 100644 --- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C @@ -187,6 +187,11 @@ 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; From 948a909da83919a13edd18123911c66fbd08809f Mon Sep 17 00:00:00 2001 From: tlichtenegger Date: Wed, 25 Apr 2018 15:28:03 +0200 Subject: [PATCH 07/16] Adapted Sauter mean diameter and void fraction models for eff. particle volume factor in case of polydisp. parcels. --- .../cfdemParticle/cfdemCloud/cfdemCloud.C | 6 ++- .../cfdemParticle/cfdemCloud/cfdemCloud.H | 9 +++++ .../cfdemParticle/cfdemCloud/cfdemCloudI.H | 14 +++++++ .../subModels/forceModel/dSauter/dSauter.C | 40 +++++++++---------- .../subModels/forceModel/dSauter/dSauter.H | 6 +-- .../GaussVoidFraction/GaussVoidFraction.C | 3 +- .../bigParticleVoidFraction.C | 3 +- .../centreVoidFraction/centreVoidFraction.C | 3 +- .../dividedVoidFraction/dividedVoidFraction.C | 1 + .../trilinearVoidFraction.C | 3 +- .../voidFractionModel/voidFractionModel.C | 2 + .../voidFractionModel/voidFractionModel.H | 7 ++++ 12 files changed, 68 insertions(+), 29 deletions(-) diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C index 2eb91237..680658c0 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C @@ -98,6 +98,7 @@ cfdemCloud::cfdemCloud voidfractions_(NULL), cellIDs_(NULL), particleDensities_(NULL), + particleEffVolFactors_(NULL), particleTypes_(NULL), particleWeights_(NULL), particleVolumes_(NULL), @@ -376,6 +377,7 @@ cfdemCloud::~cfdemCloud() dataExchangeM().destroy(particleVolumes_,1); dataExchangeM().destroy(particleV_,1); if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1); + if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1); if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1); } @@ -390,6 +392,7 @@ void cfdemCloud::getDEMdata() 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_); } @@ -745,9 +748,8 @@ bool cfdemCloud::reAllocArrays() dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle()); dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle()); dataExchangeM().allocateArray(particleV_,0.,1); - Info << "\nnow allocating particle densities array.\n" << endl; if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1); - Info << "\nnow allocating particle types array.\n" << endl; + 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 e736c974..aeee157f 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H @@ -98,6 +98,8 @@ protected: bool getParticleDensities_; + bool getParticleEffVolFactors_; + bool getParticleTypes_; scalar maxDEMForce_; @@ -128,6 +130,8 @@ protected: mutable double **particleDensities_; + mutable double **particleEffVolFactors_; + mutable int **particleTypes_; mutable double **particleWeights_; @@ -335,9 +339,14 @@ public: virtual inline int maxType() {return -1;} virtual inline bool multipleTypesDMax() {return false;} virtual inline bool multipleTypesDMin() {return false;} + inline bool getParticleDensities() const; virtual inline double ** particleDensity() 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; diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H index 0db389e0..685a8558 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H @@ -239,6 +239,20 @@ inline double cfdemCloud::particleDensity(label index) const } } +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_; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C b/src/lagrangian/cfdemParticle/subModels/forceModel/dSauter/dSauter.C index 081742c3..e7e6c3d8 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",labelList(1,1))), 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,34 @@ void dSauter::allocateMyArrays() const void dSauter::setForce() const { - if (scaleDia_ > 1) + if (typeCG_.size()>1 || typeCG_[0] > 1) { - 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 cg = 1; + label partType = 1; + 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 +181,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..1b969b9c 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 + labelList typeCG_; mutable volScalarField d2Field_; 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 3368fc7a..424b8eef 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 877b49bf..9a760b31 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_; @@ -132,6 +134,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 From 500db10717af5672686035ce0d61fe2e832335fd Mon Sep 17 00:00:00 2001 From: tlichtenegger Date: Mon, 30 Apr 2018 15:04:21 +0200 Subject: [PATCH 08/16] Forgot to set bool getParticleEffVolFactors. --- src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C | 1 + 1 file changed, 1 insertion(+) diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C index 680658c0..9f5d8955 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C @@ -84,6 +84,7 @@ cfdemCloud::cfdemCloud 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), From 7b59ed2ffa76eadfb1337a32eb81e3d6df0b76a9 Mon Sep 17 00:00:00 2001 From: Thomas Lichtenegger Date: Thu, 17 May 2018 14:13:13 +0200 Subject: [PATCH 09/16] Fix sequence of arguments in function calculating drag. --- .../subModels/forceModel/BeetstraDrag/BeetstraDrag.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C index 8bf3ed97..b771773b 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C @@ -236,7 +236,7 @@ void BeetstraDrag::setForce() const Rep=ds_scaled*voidfraction*magUr/nuf + SMALL; - dragCoefficient = F(Rep, voidfraction) + dragCoefficient = F(voidfraction, Rep) *3*M_PI*nuf*rho*voidfraction *effDiameter(ds_scaled, voidfraction, cellI, index) *scaleDia3*scaleDrag_; From e0530dcf5bf8c11ddc9f28f57a3a7398e1a0ff70 Mon Sep 17 00:00:00 2001 From: Thomas Lichtenegger Date: Fri, 18 May 2018 12:04:57 +0200 Subject: [PATCH 10/16] Cleaned up solver, removed testing code. --- applications/solvers/cfdemSolverRhoSimple/UEqn.H | 9 --------- .../cfdemSolverRhoSimple/cfdemSolverRhoSimple.C | 2 +- applications/solvers/cfdemSolverRhoSimple/pEqn.H | 11 +---------- 3 files changed, 2 insertions(+), 20 deletions(-) diff --git a/applications/solvers/cfdemSolverRhoSimple/UEqn.H b/applications/solvers/cfdemSolverRhoSimple/UEqn.H index 45b78abe..8fbd30f2 100644 --- a/applications/solvers/cfdemSolverRhoSimple/UEqn.H +++ b/applications/solvers/cfdemSolverRhoSimple/UEqn.H @@ -27,13 +27,4 @@ else fvOptions.correct(U); -// limit vel - -forAll(U,cellI) -{ - scalar mU = mag(U[cellI]); - if (mU > 200.0) U[cellI] *= 200.0/mU; -} - - K = 0.5*magSqr(U); diff --git a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C index b505252e..4e1821fc 100644 --- a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C +++ b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C @@ -119,7 +119,7 @@ int main(int argc, char *argv[]) turbulence->correct(); particleCloud.clockM().start(32,"postFlow"); - particleCloud.postFlow(); + if(hasEvolved) particleCloud.postFlow(); particleCloud.clockM().stop("postFlow"); runTime.write(); diff --git a/applications/solvers/cfdemSolverRhoSimple/pEqn.H b/applications/solvers/cfdemSolverRhoSimple/pEqn.H index c032cbeb..5b2dbcbd 100644 --- a/applications/solvers/cfdemSolverRhoSimple/pEqn.H +++ b/applications/solvers/cfdemSolverRhoSimple/pEqn.H @@ -76,15 +76,6 @@ else U = HbyA - rAU*(fvc::grad(p)-Ksl*Us); } -// limit vel - -forAll(U,cellI) -{ - scalar mU = mag(U[cellI]); - if (mU > 200.0) U[cellI] *= 200.0/mU; -} - - U.correctBoundaryConditions(); fvOptions.correct(U); -K = 0.5*magSqr(U); +K = 0.5*magSqr(U); \ No newline at end of file From 19dd9c96469d561c39eed7f8274327c8c7a751cb Mon Sep 17 00:00:00 2001 From: Thomas Lichtenegger Date: Fri, 18 May 2018 12:08:56 +0200 Subject: [PATCH 11/16] Beetstra drag model. For polydisp case, BeetstraDragPoly overrides functions for diameter. Effect of fines phase in terms of static hold-up field is included in BeetstraDragPoly. --- .../forceModel/BeetstraDrag/BeetstraDrag.C | 7 ++-- .../forceModel/BeetstraDrag/BeetstraDrag.H | 4 +- .../BeetstraDragPoly/BeetstraDragPoly.C | 40 ++++++++++++++----- .../BeetstraDragPoly/BeetstraDragPoly.H | 12 ++++-- 4 files changed, 44 insertions(+), 19 deletions(-) diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C index b771773b..98182357 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C @@ -163,6 +163,7 @@ void BeetstraDrag::setForce() const vector Ur(0,0,0); scalar ds(0); scalar ds_scaled(0); + scalar dSauter(0); scalar scaleDia3 = typeCG_[0]*typeCG_[0]*typeCG_[0]; scalar nuf(0); scalar rho(0); @@ -225,20 +226,20 @@ void BeetstraDrag::setForce() const magUr = mag(Ur); ds = 2*particleCloud_.radius(index); 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; + Rep=dSauter*voidfraction*magUr/nuf + SMALL; dragCoefficient = F(voidfraction, Rep) *3*M_PI*nuf*rho*voidfraction - *effDiameter(ds_scaled, voidfraction, cellI, index) + *effDiameter(ds_scaled, cellI, index) *scaleDia3*scaleDrag_; // calculate filtering corrections diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H index 8b8c0055..c2128855 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H @@ -85,7 +85,9 @@ protected: virtual void adaptVoidfraction(double&, label) const {} - virtual scalar effDiameter(double d, double voidfraction, label cellI, label index) const {return d;} + 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; diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C index fe4024cd..3fb66969 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C @@ -49,14 +49,26 @@ BeetstraDragPoly::BeetstraDragPoly : BeetstraDrag(dict,sm), fines_(propsDict_.lookupOrDefault("fines",false)), - dSauter_(sm.mesh().lookupObject ("dSauter")) + dFine_(propsDict_.lookupOrDefault("dFine",1.0)) { + // if fines are present, take mixture dSauter, otherwise normal dSauter if (fines_) { + if (!propsDict_.found("dFine")) + { + FatalError << "forceModel BeetstraDragPoly: Define fines diameter." << abort(FatalError); + } + volScalarField& alphaP(const_cast(sm.mesh().lookupObject ("alphaP"))); + alphaP_.set(&alphaP); volScalarField& alphaSt(const_cast(sm.mesh().lookupObject ("alphaSt"))); alphaSt_.set(&alphaSt); - volScalarField& dSauterMix(const_cast(sm.mesh().lookupObject ("dSauterMix"))); - dSauterMix_.set(&dSauterMix); + volScalarField& dSauter(const_cast(sm.mesh().lookupObject ("dSauterMix"))); + dSauter_.set(&dSauter); + } + else + { + volScalarField& dSauter(const_cast(sm.mesh().lookupObject ("dSauter"))); + dSauter_.set(&dSauter); } } @@ -75,24 +87,30 @@ void BeetstraDragPoly::adaptVoidfraction(double& voidfraction, label cellI) cons if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_; } -scalar BeetstraDragPoly::effDiameter(double d, double voidfraction, label cellI, label index) const +scalar BeetstraDragPoly::effDiameter(double d, label cellI, label index) const { - scalar dS = dSauter_[cellI]; - scalar effD = d*d/dS; + scalar dS = dSauter_()[cellI]; + scalar effD = d*d / dS + 0.064*d*d*d*d / (dS*dS*dS); + if (fines_) { - scalar pureVoidfraction = voidfraction_[cellI]; - scalar dSmix = dSauterMix_()[cellI]; - effD *= pureVoidfraction / voidfraction * (1 - voidfraction) / (1 - pureVoidfraction); - effD *= dS * dS / (dSmix * dSmix); + 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]; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H index 12bde2a3..88965010 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H @@ -43,16 +43,20 @@ class BeetstraDragPoly protected: const bool fines_; + + scalar dFine_; + autoPtr alphaP_; + autoPtr alphaSt_; - const volScalarField& dSauter_; - - autoPtr dSauterMix_; + autoPtr dSauter_; void adaptVoidfraction(double&, label) const; - scalar effDiameter(double, double, label, label) const; + scalar effDiameter(double, label, label) const; + + scalar meanSauterDiameter(double, label) const; public: From 05fb204dc397d7a716630237f8604e46257b5453 Mon Sep 17 00:00:00 2001 From: Thomas Lichtenegger Date: Fri, 18 May 2018 14:10:35 +0200 Subject: [PATCH 12/16] Some improvements based on helpful comments. --- .../cfdemParticle/cfdemCloud/cfdemCloudI.H | 12 ++------- .../forceModel/BeetstraDrag/BeetstraDrag.C | 25 +++++++++++-------- .../BeetstraDragPoly/BeetstraDragPoly.C | 7 ++---- 3 files changed, 19 insertions(+), 25 deletions(-) diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H index 685a8558..187d043f 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H @@ -223,11 +223,7 @@ inline bool cfdemCloud::getParticleDensities() const inline double ** cfdemCloud::particleDensity() const { - if(!getParticleDensities_) return NULL; - else - { - return particleDensities_; - } + return particleDensities_; } inline double cfdemCloud::particleDensity(label index) const @@ -260,11 +256,7 @@ inline bool cfdemCloud::getParticleTypes() const inline int ** cfdemCloud::particleTypes() const { - if(!getParticleTypes_) return NULL; - else - { - return particleTypes_; - } + return particleTypes_; } inline label cfdemCloud::particleType(label index) const diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C index 98182357..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 \*---------------------------------------------------------------------------*/ @@ -214,7 +215,7 @@ void BeetstraDrag::setForce() const // in case a fines phase is present, void fraction needs to be adapted adaptVoidfraction(voidfraction, cellI); - if (multiTypes_) + if (multiTypes_) { partType = particleCloud_.particleType(index); cg = typeCG_[partType - 1]; @@ -241,14 +242,14 @@ void BeetstraDrag::setForce() const *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_ - /Foam::pow(U_.mesh().V()[cellI],.33333333) + GCcorr = 1.-h(localPhiP) + /( a(localPhiP) + *Lc2_ + /std::cbrt(U_.mesh().V()[cellI]) + 1. ); if (usePC_) @@ -311,8 +312,8 @@ void BeetstraDrag::setForce() const 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)) + 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 @@ -362,7 +363,9 @@ double BeetstraDrag::a(double phiP) const { a0m = 1.79; } - return a0m + a1m*(phiP-phipam) + a2m*pow(phiP-phipam,2.) + a3m*pow(phiP-phipam,3.); + 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 @@ -405,7 +408,9 @@ double BeetstraDrag::h(double phiP) const { h0m = 0.680; h1m = -2.340; h2m = -225.200; phiphm = 0.55; } - return h0m + h1m*(phiP-phiphm) + h2m*pow(phiP-phiphm,2) + h3m*pow(phiP-phiphm,3); + 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 diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C index 3fb66969..0d886657 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C @@ -49,15 +49,12 @@ BeetstraDragPoly::BeetstraDragPoly : BeetstraDrag(dict,sm), fines_(propsDict_.lookupOrDefault("fines",false)), - dFine_(propsDict_.lookupOrDefault("dFine",1.0)) + dFine_(1.0) { // if fines are present, take mixture dSauter, otherwise normal dSauter if (fines_) { - if (!propsDict_.found("dFine")) - { - FatalError << "forceModel BeetstraDragPoly: Define fines diameter." << abort(FatalError); - } + dFine_ = readScalar(propsDict_.lookup("dFine")); volScalarField& alphaP(const_cast(sm.mesh().lookupObject ("alphaP"))); alphaP_.set(&alphaP); volScalarField& alphaSt(const_cast(sm.mesh().lookupObject ("alphaSt"))); From 554577d0bab3a7000ee8ba3abf0323e6fb221c33 Mon Sep 17 00:00:00 2001 From: Thomas Lichtenegger Date: Fri, 18 May 2018 14:30:42 +0200 Subject: [PATCH 13/16] Some more improvements. --- src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H | 2 +- src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H | 2 +- .../ZehnerSchluenderThermCond.C | 10 +++++----- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H index aeee157f..6f46acd1 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H @@ -341,7 +341,7 @@ public: virtual inline bool multipleTypesDMin() {return false;} inline bool getParticleDensities() const; - virtual inline double ** particleDensity() const; + virtual inline double ** particleDensities() const; virtual inline double particleDensity(label index) const; inline bool getParticleEffVolFactors() const; diff --git a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H index 187d043f..02b3a7dc 100644 --- a/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H +++ b/src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H @@ -221,7 +221,7 @@ inline bool cfdemCloud::getParticleDensities() const return getParticleDensities_; } -inline double ** cfdemCloud::particleDensity() const +inline double ** cfdemCloud::particleDensities() const { return particleDensities_; } diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C index 1d1cec33..e7647b7a 100644 --- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C @@ -51,7 +51,7 @@ ZehnerSchluenderThermCond::ZehnerSchluenderThermCond thermCondModel(dict,sm), propsDict_(dict.subDict(typeName + "Props")), partKsField_ - ( + ( IOobject ( "partKs", @@ -145,7 +145,7 @@ tmp ZehnerSchluenderThermCond::thermCond() const scalar B = 0.0; scalar C = 0.0; scalar k = 0.0; - scalar OnemBoA = 0.0; + scalar InvOnemBoA = 0.0; scalar voidfraction = 0.0; scalar w = 7.26e-3; @@ -157,9 +157,9 @@ tmp ZehnerSchluenderThermCond::thermCond() const { A = partKsField_[cellI]/kf0_.value(); B = 1.25 * Foam::pow((1 - voidfraction) / voidfraction, 1.11); - OnemBoA = 1.0 - B/A; - C = (A - 1) / (OnemBoA * OnemBoA) * B/A * log(A/B) - (B - 1)/OnemBoA - 0.5 * (B + 1); - C *= 2.0 / OnemBoA; + 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); } From c1cb77f08c6d3afab4f4f432c9aef029dc229709 Mon Sep 17 00:00:00 2001 From: Thomas Lichtenegger Date: Tue, 22 May 2018 09:42:24 +0200 Subject: [PATCH 14/16] Replaced hardcoded temperature limits with user-defined ones. --- .../heatTransferGunnPartField.C | 21 ++++++++++++++----- .../heatTransferGunnPartField.H | 4 ++++ 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C index 2dac2222..d4a116d0 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C @@ -56,11 +56,13 @@ heatTransferGunnPartField::heatTransferGunnPartField ), sm.mesh(), dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0), - "zeroGradient" + "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 @@ -80,6 +82,16 @@ heatTransferGunnPartField::heatTransferGunnPartField { 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; @@ -216,10 +228,9 @@ void heatTransferGunnPartField::solve() fvOptions.correct(partTempField_); - dimensionedScalar pTMin("pTMin",dimensionSet(0,0,0,1,0,0,0),293.0); - dimensionedScalar pTMax("pTMin",dimensionSet(0,0,0,1,0,0,0),3000.0); - partTempField_ = max(partTempField_, pTMin); - partTempField_ = min(partTempField_, pTMax); + 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; diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H index a321be7e..296c9a49 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.H @@ -53,6 +53,10 @@ private: mutable double **partCp_; + dimensionedScalar pTMax_; + + dimensionedScalar pTMin_; + autoPtr thermCondModel_; fv::options& fvOptions; From 0d7916746ab834af414d6ed6f71a81882eceb76a Mon Sep 17 00:00:00 2001 From: Thomas Lichtenegger Date: Tue, 22 May 2018 15:23:35 +0200 Subject: [PATCH 15/16] Added return NULL to silence compiler. --- .../ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C | 1 + 1 file changed, 1 insertion(+) diff --git a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C index e7647b7a..7e4ca613 100644 --- a/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C +++ b/src/lagrangian/cfdemParticle/subModels/thermCondModel/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C @@ -181,6 +181,7 @@ tmp ZehnerSchluenderThermCond::thermCond() const tmp ZehnerSchluenderThermCond::thermDiff() const { FatalError << "ZehnerSchluenderThermCond does not provide thermal diffusivity." << abort(FatalError); + return tmp(NULL); } void ZehnerSchluenderThermCond::calcPartKsField() const From 06633e9e4d5ba4c4e6d65e38468f2943bdc8480b Mon Sep 17 00:00:00 2001 From: Thomas Lichtenegger Date: Tue, 22 May 2018 15:24:36 +0200 Subject: [PATCH 16/16] Reset array where necessary. --- .../heatTransferGunnPartField/heatTransferGunnPartField.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C index d4a116d0..2d8af0e6 100644 --- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C +++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C @@ -164,6 +164,8 @@ void heatTransferGunnPartField::postFlow() 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];