yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12669
[Branch ~yade-pkg/yade/git-trunk] Rev 3863: HydroForceEngine: modify turbulent fluctuation model.
------------------------------------------------------------
revno: 3863
committer: Raphael Maurin <raph_maurin@xxxxxxxxxxx>
timestamp: Thu 2016-05-19 16:19:12 +0200
message:
HydroForceEngine: modify turbulent fluctuation model.
Add a fluctuation along the spanwise (y) direction.
Reset the fluctuation to zero when the particle is out of the flow (fix a bug)
Comment a bit more the turbulentFluctuation function
modified:
pkg/common/ForceEngine.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== modified file 'pkg/common/ForceEngine.cpp'
--- pkg/common/ForceEngine.cpp 2016-02-09 17:55:09 +0000
+++ pkg/common/ForceEngine.cpp 2016-05-19 14:19:12 +0000
@@ -97,7 +97,7 @@
if ((p<nCell)&&(p>0)) {
Vector3r liftForce = Vector3r::Zero();
Vector3r dragForce = Vector3r::Zero();
- Vector3r vRel = Vector3r(vxFluid[p]+vFluctX[id],0.0,vFluctZ[id]) - b->state->vel;//fluid-particle relative velocity
+ Vector3r vRel = Vector3r(vxFluid[p]+vFluctX[id],vFluctY[id],vFluctZ[id]) - b->state->vel;//fluid-particle relative velocity
//Drag force calculation
if (vRel.norm()!=0.0) {
dragForce = 0.5*densFluid*Mathr::PI*pow(sphere->radius,2.0)*(0.44*vRel.norm()+24.4*viscoDyn/(densFluid*sphere->radius*2))*pow(1-phiPart[p],-expoRZ)*vRel;
@@ -155,7 +155,7 @@
// Relative fluid/particle velocity using also the associated fluid vel. fluct.
if ((Np>=0)&&(Np<nCell)){
- uRel = Vector3r(vxFluid[Np]+vFluctX[b->id], 0.0,vFluctZ[b->id]) - b->state->vel;
+ uRel = Vector3r(vxFluid[Np]+vFluctX[b->id], vFluctY[b->id],vFluctZ[b->id]) - b->state->vel;
// Drag force with a Dallavalle formulation (drag coef.) and Richardson-Zaki Correction (hindrance effect)
fDrag = 0.5*Mathr::PI*pow(s->radius,2.0)*densFluid*(0.44*uRel.norm()+24.4*viscoDyn/(densFluid*2.0*s->radius))*pow((1-phiPart[Np]),-expoRZ)*uRel;
}
@@ -287,10 +287,12 @@
if(size<scene->bodies->size()){
size=scene->bodies->size();
vFluctX.resize(size);
+ vFluctY.resize(size);
vFluctZ.resize(size);
}
/* reset stored values to zero */
memset(& vFluctX[0],0,size);
+ memset(& vFluctY[0],0,size);
memset(& vFluctZ[0],0,size);
/* Create a random number generator rnd() with a gaussian distribution of mean 0 and stdev 1.0 */
@@ -301,6 +303,7 @@
double rand1 = 0.0;
double rand2 = 0.0;
+ double rand3 = 0.0;
/* Attribute a fluid velocity fluctuation to each body above the bed elevation */
FOREACH(Body::id_t id, ids){
Body* b=Body::byId(id,scene).get();
@@ -310,37 +313,53 @@
if (sphere){
Vector3r posSphere = b->state->pos;//position vector of the sphere
int p = floor((posSphere[2]-zRef)/deltaZ); //cell number in which the particle is
- // if the particle is inside the water and above the bed elevation, so inside the turbulent flow, evaluate a turbulent fluid velocity fluctuation which will be used to apply the drag.
- if ((p<nCell)&&(posSphere[2]-zRef>bedElevation)) {
+ // If the particle is inside the water and above the bed elevation, so inside the turbulent flow, evaluate a turbulent fluid velocity fluctuation which will be used to apply the drag.
+ // The fluctuation magnitude is linked to the value of the Reynolds stress tensor at the given position, a kind of local friction velocity ustar
+ // The fluctuations along wall-normal and streamwise directions are correlated in order to be consistent with the formulation of the Reynolds stress tensor and to recover the result
+ // that the magnitude of the fluctuation along streamwise = 2*along wall normal
+ if ((p<nCell)&&(posSphere[2]-zRef>bedElevation)) { // Remove the particles outside of the flow and inside the granular bed, they are not submitted to turbulent fluctuations.
Real uStar2 = simplifiedReynoldStresses[p];
if (uStar2>0.0){
Real uStar = sqrt(uStar2);
rand1 = rnd();
- rand2 = -rand1 + rnd();
+ rand2 = rnd();
+ rand3 = -rand1 + rnd();// x and z fluctuation are correlated as measured by Nezu 1977 and as expected from the formulation of the Reynolds stress tensor.
vFluctZ[id] = rand1*uStar;
- vFluctX[id] = rand2*uStar;
+ vFluctY[id] = rand2*uStar;
+ vFluctX[id] = rand3*uStar;
}
}
+ else{
+ vFluctZ[id] = 0.0;
+ vFluctY[id] = 0.0;
+ vFluctX[id] = 0.0;
+
+ }
}
}
}
-/* Alternative Velocity fluctuation determination, fluctuation time step depend on z. Should be executed at a period dtFluct corresponding to the smallest fluctTime */
+/* Alternative Velocity fluctuation model, same as turbulentFluctuation model but with a time step associated with the fluctuation generation depending on z */
+/* Should be executed in the python script at a period dtFluct corresponding to the smallest value of the fluctTime vector */
/* Should be initialized before running HydroForceEngine */
void HydroForceEngine::turbulentFluctuationBIS(){
int idPartMax = vFluctX.size();
double rand1 = 0.0;
double rand2 = 0.0;
+ double rand3 = 0.0;
/* Create a random number generator rnd() with a gaussian distribution of mean 0 and stdev 1.0 */
/* see http://www.boost.org/doc/libs/1_55_0/doc/html/boost_random/reference.html and the chapter 7 of Numerical Recipes in C, second edition (1992) for more details */
static boost::minstd_rand0 randGen((int)TimingInfo::getNow(true));
static boost::normal_distribution<Real> dist(0.0, 1.0);
static boost::variate_generator<boost::minstd_rand0&,boost::normal_distribution<Real> > rnd(randGen,dist);
+ //Loop on the particles
for(int idPart=0;idPart<idPartMax;idPart++){
- fluctTime[idPart]-=dtFluct; //define dtFluct in global
- if (fluctTime[idPart]<=0){
- fluctTime[idPart] = 10*dtFluct;
+ //Remove the time ran since last application of the function (dtFluct define in global)
+ fluctTime[idPart]-=dtFluct;
+ //If negative, means that the time of application of the fluctuation is over, generate a new one with a new associated time
+ if (fluctTime[idPart]<=0){
+ fluctTime[idPart] = 10*dtFluct; //Initialisation of the application time
Body* b=Body::byId(idPart,scene).get();
if (!b) continue;
if (!(scene->bodies->exists(idPart))) continue;
@@ -350,15 +369,26 @@
Vector3r posSphere = b->state->pos;//position vector of the sphere
int p = floor((posSphere[2]-zRef)/deltaZ); //cell number in which the particle is
if (simplifiedReynoldStresses[p]>0.0) uStar = sqrt(simplifiedReynoldStresses[p]);
- // if the particle is inside the water and above the bed elevation, so inside the turbulent flow, evaluate a turbulent fluid velocity fluctuation which will be used to apply the drag.
+ // Remove the particles outside of the flow and inside the granular bed, they are not submitted to turbulent fluctuations.
if ((p<nCell)&&(posSphere[2]-zRef>bedElevation)) {
rand1 = rnd();
- rand2 = -rand1 + rnd();
+ rand2 = rnd();
+ rand3 = -rand1 + rnd(); // x and z fluctuation are correlated as measured by Nezu 1977 and as expected from the formulation of the Reynolds stress tensor.
vFluctZ[idPart] = rand1*uStar;
- vFluctX[idPart] = rand2*uStar;
+ vFluctY[idPart] = rand2*uStar;
+ vFluctX[idPart] = rand3*uStar;
+ // Limit the value to avoid the application of fluctuations in the viscous sublayer
const Real zPos = max(b->state->pos[2]-zRef-bedElevation,11.6*viscoDyn/densFluid/uStar);
+ // Time of application of the fluctuation as a function of depth from kepsilon model
if (uStar>0.0) fluctTime[idPart]=min(0.33*0.41*zPos/uStar,10.);
}
+ else{
+ vFluctZ[idPart] = 0.0;
+ vFluctY[idPart] = 0.0;
+ vFluctX[idPart] = 0.0;
+ fluctTime[idPart] = 0.0;
+
+ }
}
}
}
