yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #12530
[Branch ~yade-pkg/yade/git-trunk] Rev 3790: HydroForceEngine: modify the turbulent fluctuations formulation and add a new DRW model
------------------------------------------------------------
revno: 3790
committer: Raphael Maurin <raph_maurin@xxxxxxxxxxx>
timestamp: Tue 2016-02-09 18:45:15 +0100
message:
HydroForceEngine: modify the turbulent fluctuations formulation and add a new DRW model
The turbulent fluctuations model takes now the form of a function of HydroForceEngine and is not anymore called from a flag.
In addition, an alternative DRW model has been added.
modified:
pkg/common/ForceEngine.cpp
pkg/common/ForceEngine.hpp
--
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 2015-07-21 16:19:06 +0000
+++ pkg/common/ForceEngine.cpp 2016-02-09 17:45:15 +0000
@@ -82,54 +82,7 @@
}
}
-void HydroForceEngine::action(){
- /* Velocity fluctuation determination (not usually done at each dt, that is why it not placed in the other loop) */
- if (velFluct == true){
- /* check size */
- size_t size=vFluctX.size();
- if(size<scene->bodies->size()){
- size=scene->bodies->size();
- vFluctX.resize(size);
- vFluctZ.resize(size);
- }
- /* reset stored values to zero */
- memset(& vFluctX[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 */
- /* 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);
-
- double rand1 = 0.0;
- double rand2 = 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();
- if (!b) continue;
- if (!(scene->bodies->exists(id))) continue;
- const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
- 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)) {
- Real uStar2 = simplifiedReynoldStresses[p];
- if (uStar2>0.0){
- Real uStar = sqrt(uStar2);
- rand1 = rnd();
- rand2 = -rand1 + rnd();
- vFluctZ[id] = rand1*uStar;
- vFluctX[id] = rand2*uStar;
- }
- }
- }
-
- }
- velFluct = false;
- }
-
+void HydroForceEngine::action(){
/* Application of hydrodynamical forces */
if (activateAverage==true) averageProfile(); //Calculate the average solid profiles
@@ -177,7 +130,7 @@
Vector3r uRel = Vector3r::Zero();
Vector3r fDrag = Vector3r::Zero();
- int nMax = 2*nCell;
+ int nMax = nCell;
vector<Real> velAverageX(nMax,0.0);
vector<Real> velAverageY(nMax,0.0);
vector<Real> velAverageZ(nMax,0.0);
@@ -282,3 +235,91 @@
}
+/* Velocity fluctuation determination. To execute at a given (changing) period corresponding to the eddy turn over time*/
+/* Should be initialized before running HydroForceEngine */
+void HydroForceEngine::turbulentFluctuation(){
+ /* check size */
+ size_t size=vFluctX.size();
+ if(size<scene->bodies->size()){
+ size=scene->bodies->size();
+ vFluctX.resize(size);
+ vFluctZ.resize(size);
+ }
+ /* reset stored values to zero */
+ memset(& vFluctX[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 */
+ /* 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);
+
+ double rand1 = 0.0;
+ double rand2 = 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();
+ if (!b) continue;
+ if (!(scene->bodies->exists(id))) continue;
+ const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
+ 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)) {
+ Real uStar2 = simplifiedReynoldStresses[p];
+ if (uStar2>0.0){
+ Real uStar = sqrt(uStar2);
+ rand1 = rnd();
+ rand2 = -rand1 + rnd();
+ vFluctZ[id] = rand1*uStar;
+ vFluctX[id] = rand2*uStar;
+ }
+ }
+ }
+ }
+}
+
+/* Alternative Velocity fluctuation determination, fluctuation time step depend on z. Should be executed at a period dtFluct corresponding to the smallest fluctTime */
+/* Should be initialized before running HydroForceEngine */
+void HydroForceEngine::turbulentFluctuationBIS(){
+ int idPartMax = vFluctX.size();
+ double rand1 = 0.0;
+ double rand2 = 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);
+
+ for(int idPart=0;idPart<idPartMax;idPart++){
+ fluctTime[idPart]-=dtFluct; //define dtFluct in global
+ if (fluctTime[idPart]<=0){
+ fluctTime[idPart] = 10*dtFluct;
+ Body* b=Body::byId(idPart,scene).get();
+ if (!b) continue;
+ if (!(scene->bodies->exists(idPart))) continue;
+ const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
+ Real uStar = 0.0;
+ 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 (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.
+ if ((p<nCell)&&(posSphere[2]-zRef>bedElevation)) {
+ rand1 = rnd();
+ rand2 = -rand1 + rnd();
+ vFluctZ[idPart] = rand1*uStar;
+ vFluctX[idPart] = rand2*uStar;
+ const Real zPos = max(b->state->pos[2]-zRef-bedElevation,11.6*viscoDyn/densFluid/uStar);
+ if (uStar>0.0) fluctTime[idPart]=min(0.33*0.41*zPos/uStar,10.);
+ }
+ }
+ }
+ }
+}
+
+
+
+
=== modified file 'pkg/common/ForceEngine.hpp'
--- pkg/common/ForceEngine.hpp 2015-07-21 16:19:06 +0000
+++ pkg/common/ForceEngine.hpp 2016-02-09 17:45:15 +0000
@@ -73,6 +73,8 @@
class HydroForceEngine: public PartialEngine{
public:
void averageProfile();
+ void turbulentFluctuation();
+ void turbulentFluctuationBIS();
public:
virtual void action();
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(HydroForceEngine,PartialEngine,"Apply drag and lift due to a fluid flow vector (1D) to each sphere + the buoyant weight.\n The applied drag force reads\n\n $F_{d}=\\frac{1}{2} C_d A\\rho^f|\\vec{v_f - v}| \\vec{v_f - v}$ \n\n where $\\rho$ is the medium density (:yref:`densFluid<HydroForceEngine.densFluid>`), $v$ is particle's velocity, $v_f$ is the velocity of the fluid at the particle center(:yref:`vxFluid<HydroForceEngine.vxFluid>`), $A$ is particle projected area (disc), $C_d$ is the drag coefficient. The formulation of the drag coefficient depends on the local particle reynolds number and the solid volume fraction. The formulation of the drag is [Dallavalle1948]_ [RevilBaudard2013]_ with a correction of Richardson-Zaki [Richardson1954]_ to take into account the hindrance effect. This law is classical in sediment transport. It is possible to activate a fluctuation of the drag force for each particle which account for the turbulent fluctuation of the fluid velocity (:yref:`velFluct<HydroForceEngine.velFluct>`). The model implemented for the turbulent velocity fluctuation is a simple discrete random walk which takes as input the Reynolds stress tensor $R^f_{xz}$ as a function of the depth, and allows to recover the main property of the fluctuations by imposing $<u_x'u_z'> (z) = <R^f_{xz}>(z)/\\rho^f$. It requires as input $<R^f_{xz}>(z)/\\rho^f$ called :yref:`simplifiedReynoldStresses<HydroForceEngine.simplifiedReynoldStresses>` in the code. \n The formulation of the lift is taken from [Wiberg1985]_ and is such that : \n\n $F_{L}=\\frac{1}{2} C_L A\\rho^f((v_f - v)^2_{top} - (v_f - v)^2_{bottom})$ \n\n Where the subscript top and bottom means evaluated at the top (respectively the bottom) of the sphere considered. This formulation of the lift account for the difference of pressure at the top and the bottom of the particle inside a turbulent shear flow. As this formulation is controversial when approaching the threshold of motion [Schmeeckle2007]_ it is possible to desactivate it with the variable :yref:`lift<HydroForceEngine.lift>`.\n The buoyancy is taken into account through the buoyant weight : \n\n $F_{buoyancy}= - \\rho^f V^p g$ \n\n, where g is the gravity vector along the vertical, and $V^p$ is the volume of the particle. This engine also evaluate the average particle velocity, solid volume fraction and drag force depth profiles, through the function averageProfile. This is done as the solid volume fraction depth profile is required for the drag calculation, and as the three are required for the independent fluid resolution.",
@@ -105,9 +107,14 @@
((vector<Real>,vFluctZ,,,"Vector associating a normal fluid velocity fluctuation to each particle. Fluctuation calculated in the C++ code from the discrete random walk model"))
((vector<Real>,simplifiedReynoldStresses,,,"Vector of size equal to :yref:`nCell<HydroForceEngine.nCell>` containing the Reynolds stresses divided by the fluid density in function of the depth. simplifiedReynoldStresses(z) $= <u_x'u_z'>(z)^2$"))
((Real,bedElevation,,,"Elevation of the bed above which the fluid flow is turbulent and the particles undergo turbulent velocity fluctuation."))
+ ((vector<Real>,fluctTime,,,"Vector containing the time of life of the fluctuations associated to each particles."))
+ ((Real,dtFluct,,,"Execution time step of the turbulent fluctuation model."))
+
,/*ctor*/
,/*py*/
.def("averageProfile",&HydroForceEngine::averageProfile,"Compute and store the particle velocity (:yref:`vxPart<HydroForceEngine.vxPart>`, :yref:`vyPart<HydroForceEngine.vyPart>`, :yref:`vzPart<HydroForceEngine.vzPart>`) and solid volume fraction (:yref:`phiPart<HydroForceEngine.phiPart>`) depth profile. For each defined cell z, the k component of the average particle velocity reads: \n\n $<v_k>^z= \\sum_p V^p v_k^p/\\sum_p V^p$,\n\n where the sum is made over the particles contained in the cell, $v_k^p$ is the k component of the velocity associated to particle p, and $V^p$ is the part of the volume of the particle p contained inside the cell. This definition allows to smooth the averaging, and is equivalent to taking into account the center of the particles only when there is a lot of particles in each cell. As for the solid volume fraction, it is evaluated in the same way: for each defined cell z, it reads: \n\n $<\\phi>^z= \\frac{1}{V_{cell}}\\sum_p V^p$, where $V_{cell}$ is the volume of the cell considered, and $V^p$ is the volume of particle p contained in cell z.\n This function gives depth profiles of average velocity and solid volume fraction, returning the average quantities in each cell of height dz, from the reference horizontal plane at elevation :yref:`zRef<HydroForceEngine.zRef>` (input parameter) until the plane of elevation :yref:`zRef<HydroForceEngine.zRef>` plus :yref:`nCell<HydroForceEngine.nCell>` times :yref:`deltaZ<HydroForceEngine.deltaZ>` (input parameters). When the option :yref:`twoSize<HydroForceEngine.twoSize>` is set to True, evaluate in addition the average drag (:yref:`averageDrag1<HydroForceEngine.averageDrag1>` and :yref:`averageDrag2<HydroForceEngine.averageDrag2>`) and solid volume fraction (:yref:`phiPart1<HydroForceEngine.phiPart1>` and :yref:`phiPart2<HydroForceEngine.phiPart2>`) depth profiles considering only the particles of radius respectively :yref:`radiusPart1<HydroForceEngine.radiusPart1>` and :yref:`radiusPart2<HydroForceEngine.radiusPart2>` in the averaging.")
+ .def("turbulentFluctuation",&HydroForceEngine::turbulentFluctuation,"Apply turbulent fluctuation to the problem.")
+ .def("turbulentFluctuationBIS",&HydroForceEngine::turbulentFluctuationBIS,"Apply turbulent fluctuation to the problem with an alternative formulation.")
);
};
REGISTER_SERIALIZABLE(HydroForceEngine);
