yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 4144
committer: Raphael Maurin <raph_maurin@xxxxxxxxxxx>
timestamp: Tue 2014-08-19 15:59:12 +0200
message:
  New force engine to couple Yade with a 1D RANS code.
  Add a new force engine applying the main hydrodynamical forces in function of a 1D average fluid velocity vector which depends only on the depth.
  The engine is calculating at each time step the drag, lift and buoyant forces for each particle.
  Complete the references for the documentation associated to the engine.
modified:
  doc/references.bib
  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 'doc/references.bib'
--- doc/references.bib	2014-07-28 06:24:36 +0000
+++ doc/references.bib	2014-08-19 13:59:12 +0000
@@ -735,3 +735,50 @@
 	pages={447-454},
 	language={English}
 }
+
+@article{Richarson1954,
+  author = {Richardson, J. F., and W. N. Zaki},
+  title = {Sedimentation and fluidization: Part i},
+  journal = {Trans. Instn. Chem. Engrs},
+  year = {1954},
+  volume = {32}
+}
+
+@article{RevilBaudard2013,
+  author = {Revil-Baudard, T. and Chauchat, J.},
+  title = {A two-phase model for sheet flow regime based on dense granular flow
+	rheology},
+  journal = {Journal of Geophysical Research: Oceans},
+  year = {2013},
+  volume = {118},
+  pages = {619--634},
+  number = {2}
+}
+
+@article{Schmeeckle2007,
+  author = {Schmeeckle, Mark W. and Nelson, Jonathan M. and Shreve, Ronald L.},
+  title = {Forces on stationary particles in near-bed turbulent flows},
+  journal = {Journal of Geophysical Research: Earth Surface},
+  year = {2007},
+  volume = {112},
+  number = {F2},
+  doi = {10.1029/2006JF000536},
+  url = {http://dx.doi.org/10.1029/2006JF000536}
+}
+
+@article{Wiberg1985,
+  author = {Wiberg, Patricia L. and Smith, J. Dungan},
+  title = {A theoretical model for saltating grains in water},
+  journal = {Journal of Geophysical Research: Oceans},
+  year = {1985},
+  volume = {90},
+  pages = {7341--7354},
+  number = {C4}
+}
+@book{Dallavalle1948,
+  title = {Micrometrics : The technology of fine particles},
+  publisher = {Pitman Pub. Corp},
+  year = {1948},
+  author = {J. M. DallaValle},
+  volume = {2nd edition}
+}

=== modified file 'pkg/common/ForceEngine.cpp'
--- pkg/common/ForceEngine.cpp	2011-04-20 15:52:11 +0000
+++ pkg/common/ForceEngine.cpp	2014-08-19 13:59:12 +0000
@@ -1,6 +1,6 @@
 // 2004 © Janek Kozicki <cosurgi@xxxxxxxxxx> 
 // 2009 © Václav Šmilauer <eudoxos@xxxxxxxx> 
-
+// 2014 © Raphael Maurin <raphael.maurin@xxxxxxxxx> 
 
 #include"ForceEngine.hpp"
 #include<yade/core/Scene.hpp>
@@ -11,7 +11,11 @@
 #include<yade/core/IGeom.hpp>
 #include<yade/core/IPhys.hpp>
 
-YADE_PLUGIN((ForceEngine)(InterpolatingDirectedForceEngine)(RadialForceEngine)(DragEngine)(LinearDragEngine));
+#include <boost/random/linear_congruential.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <boost/random/variate_generator.hpp>
+
+YADE_PLUGIN((ForceEngine)(InterpolatingDirectedForceEngine)(RadialForceEngine)(DragEngine)(LinearDragEngine)(HydroForceEngine));
 
 void ForceEngine::action(){
 	FOREACH(Body::id_t id, ids){
@@ -77,3 +81,90 @@
 		}
 	}
 }
+
+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=vFluct.size();
+		if(size<scene->bodies->size()){
+			size=scene->bodies->size();
+			vFluct.resize(size);
+		}
+		/* reset stored values to zero */
+		memset(& vFluct[0],0,sizeof(Vector2r)*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();
+						vFluct[id] = Vector2r(rand1*uStar,rand2*uStar);
+					}
+				}
+			}
+			
+		}
+		velFluct = false;
+	}
+	
+	/* Application of hydrodynamical forces */
+	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 ((p<nCell)&&(p>0)) {
+				Vector3r liftForce = Vector3r::Zero();
+				Vector3r dragForce = Vector3r::Zero();
+				Vector2r fluctVelBody = vFluct[id];//fluid velocity fluctuation associated to the particle's position considered.
+				Vector3r vFluid(vxFluid[p]+fluctVelBody.x(),0.0,fluctVelBody.y()); //fluid velocity at the particle's position
+				Vector3r vPart = b->state->vel;//particle velocity
+				Vector3r vRel = vFluid - vPart;//fluid-particle relative velocity
+
+				//Drag force calculation
+				Real Rep = vRel.norm()*sphere->radius*2*rhoFluid/viscoDyn; //particles reynolds number
+				Real A = sphere->radius*sphere->radius*Mathr::PI;	//Crossection of the sphere
+				if (vRel.norm()!=0.0) {
+					Real hindranceF = pow(1-phiPart[p],-expoRZ); //hindrance function
+					Real Cd = (0.44 + 24.4/Rep)*hindranceF; //drag coefficient
+					dragForce = 0.5*rhoFluid*A*Cd*vRel.squaredNorm()*vRel.normalized();
+				}
+				//lift force calculation due to difference of fluid pressure between top and bottom of the particle
+				int intRadius = floor(sphere->radius/deltaZ);
+				if ((p+intRadius<nCell)&&(p-intRadius>0)&&(lift==true)) {
+					Real vRelTop = vxFluid[p+intRadius] - vPart[0]; // relative velocity of the fluid wrt the particle at the top of the particle
+					Real vRelBottom = vxFluid[p-intRadius] - vPart[0]; // same at the bottom
+					liftForce[2] = 0.5*rhoFluid*A*Cl*(vRelTop*vRelTop-vRelBottom*vRelBottom);
+				}
+				//buoyant weight force calculation
+				Vector3r buoyantForce = -4.0/3.0*Mathr::PI*sphere->radius*sphere->radius*sphere->radius*rhoFluid*gravity;
+				//add the hydro forces to the particle
+				scene->forces.addForce(id,dragForce+liftForce+buoyantForce);		
+			}
+		}
+	}
+}
+

=== modified file 'pkg/common/ForceEngine.hpp'
--- pkg/common/ForceEngine.hpp	2011-04-20 15:52:11 +0000
+++ pkg/common/ForceEngine.hpp	2014-08-19 13:59:12 +0000
@@ -1,5 +1,6 @@
 // 2004 © Janek Kozicki <cosurgi@xxxxxxxxxx> 
 // 2009 © Václav Šmilauer <eudoxos@xxxxxxxx> 
+// 2014 © Raphael Maurin <raphael.maurin@xxxxxxxxx> 
 
 #pragma once
 
@@ -67,3 +68,29 @@
 	);
 };
 REGISTER_SERIALIZABLE(LinearDragEngine);
+
+
+class HydroForceEngine: public PartialEngine{
+	private:
+		vector<Vector2r> vFluct;
+	public:
+		virtual void action();
+	YADE_CLASS_BASE_DOC_ATTRS(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.. math:: 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:`density<HydroForceEngine.rhoFluid>`), $v$ is particle's velocity,  $v_f$ is the velocity of the fluid at the particle center,  $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`). The model implemented for the turbulent velocity fluctuation is a simple discrete random walk which takes as input the reynolds stress tensor Re_{xz} in function of the depth and allows to recover the main property of the fluctuations by imposing <u_x'u_z'> (z) = <Re>(z)/rho^f. It requires as input <Re>(z)/rho^f called :yref:`simplifiedReynoldStresses` in the code. \n The formulation of the lift is taken from [Wiberg1985]_ and is such that : \n\n.. math:: 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`.\n The buoyancy is taken into account through the buoyant weight : \n\n.. math:: 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.",
+		((Real,rhoFluid,1000,,"Density of the fluid, by default - density of water"))
+		((Real,viscoDyn,1e-3,,"Dynamic viscosity of the fluid, by default - viscosity of water"))
+		((Real,zRef,,,"Position of the reference point which correspond to the first value of the fluid velocity"))
+		((Real,nCell,,,"Size of the vector of the fluid velocity"))
+		((Real,deltaZ,,,"width of the discretization cell "))
+		((Real,expoRZ,3.1,,"Value of the Richardson-Zaki exponent, for the correction due to hindrance"))
+                ((bool,lift,true,,"Option to activate or not the evaluation of the lift"))
+		((Real,Cl,0.2,,"Value of the lift coefficient taken from [Wiberg1985]_"))
+                ((Vector3r,gravity,,,"Gravity vector (may depend on the slope)."))
+		((vector<Real>,vxFluid,,,"Discretized streamwise fluid velocity profile in function of the depth"))
+		((vector<Real>,phiPart,,,"Discretized solid volume fraction profile in function of the depth"))
+		((bool,velFluct,false,,"If true, activate the determination of turbulent fluid velocity fluctuation at the position of each particle, using a simple discrete random walk model based on the Reynolds stresses :yref:`turbStress<HydroForceEngine.squaredAverageTurbfluct>`"))
+		((vector<Real>,simplifiedReynoldStresses,,,"Vector of size equal to :yref:`turbStress<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."))
+	);
+};
+REGISTER_SERIALIZABLE(HydroForceEngine);
+