From 6a636658aa58ad5d99a86e3fb6ba48ac606492d0 Mon Sep 17 00:00:00 2001 From: Paul Kieckhefen Date: Tue, 16 May 2017 13:55:12 +0200 Subject: [PATCH] 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;