yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3569
committer: Raphael Maurin <raph_maurin@xxxxxxxxxxx>
timestamp: Tue 2015-01-27 16:52:46 +0100
message:
  HydroForceEngine: minor fix + add an averaging function
  Minor fix : modify some default value of the parameters, remove unecessary parameters.
  Add a 1D averaging function which evaluate the average depth profiles of the solid volume fraction, the solid velocity, and the drag force.
  The averaging is weighted by the volume occupied by the particles in the layer considered.
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	2014-12-12 17:01:09 +0000
+++ pkg/common/ForceEngine.cpp	2015-01-27 15:52:46 +0000
@@ -131,6 +131,8 @@
 	}
 	
 	/* Application of hydrodynamical forces */
+	if (activateAverage==true) averageProfile(); //Calculate the average fluid and velocity profile
+	
 	FOREACH(Body::id_t id, ids){
 		Body* b=Body::byId(id,scene).get();
 		if (!b) continue;
@@ -142,27 +144,20 @@
 			if ((p<nCell)&&(p>0)) {
 				Vector3r liftForce = Vector3r::Zero();
 				Vector3r dragForce = Vector3r::Zero();
-				Vector3r vFluid(vxFluid[p]+vFluctX[id],0.0,vFluctZ[id]); //fluid velocity at the particle's position
-				Vector3r vPart = b->state->vel;//particle velocity
-				Vector3r vRel = vFluid - vPart;//fluid-particle relative velocity
-
+				Vector3r vRel = Vector3r(vxFluid[p]+vFluctX[id],0.0,vFluctZ[id]) -  b->state->vel;//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();
+					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;
 				}
 				//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);
+					Real vRelTop = vxFluid[p+intRadius] -  b->state->vel[0]; // relative velocity of the fluid wrt the particle at the top of the particle
+					Real vRelBottom = vxFluid[p-intRadius] -  b->state->vel[0]; // same at the bottom
+					liftForce[2] = 0.5*densFluid*Mathr::PI*pow(sphere->radius,2.0)*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;
+				Vector3r buoyantForce = -4.0/3.0*Mathr::PI*pow(sphere->radius,3.0)*densFluid*gravity;
 				//add the hydro forces to the particle
 				scene->forces.addForce(id,dragForce+liftForce+buoyantForce);		
 			}
@@ -170,3 +165,81 @@
 	}
 }
 
+void HydroForceEngine::averageProfile(){
+	//Initialization
+	int Np;
+	int minZ;
+	int maxZ;
+	int numLayer;
+	Real deltaCenter;
+	Real zInf;
+	Real zSup;
+	Real volPart;
+	Vector3r uRel = Vector3r::Zero();
+	Vector3r fDrag  = Vector3r::Zero();
+
+	int nMax = 2*nCell;
+	vector<Real> velAverage(nMax,0.0);
+	vector<Real> phiAverage(nMax,0.0);
+	vector<Real> dragAverage(nMax,0.0);
+
+	//Loop over the particles
+	FOREACH(const shared_ptr<Body>& b, *Omega::instance().getScene()->bodies){
+		shared_ptr<Sphere> s=YADE_PTR_DYN_CAST<Sphere>(b->shape); if(!s) continue;
+		const Real zPos = b->state->pos[2]-zRef;
+		int Np = floor(zPos/deltaZ);	//Define the layer number with 0 corresponding to zRef. Let the z position wrt to zero, that way all z altitude are positive. (otherwise problem with volPart evaluation)
+
+		// 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;
+			// 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;
+		}
+		else fDrag = Vector3r::Zero();
+
+		minZ= floor((zPos-s->radius)/deltaZ);
+		maxZ= floor((zPos+s->radius)/deltaZ);
+		deltaCenter = zPos - Np*deltaZ;
+	
+		// Loop over the cell in which the particle is contained
+		numLayer = minZ;
+		while (numLayer<=maxZ){
+			if ((numLayer>=0)&&(numLayer<nMax)){ //average under zRef does not interest us, avoid also negative values not compatible with the evaluation of volPart
+				zInf=(numLayer-Np-1)*deltaZ + deltaCenter;
+				zSup=(numLayer-Np)*deltaZ + deltaCenter;
+				if (zInf<-s->radius) zInf = -s->radius;
+				if (zSup>s->radius) zSup = s->radius;
+
+				//Analytical formulation of the volume of a slice of sphere
+				volPart = Mathr::PI*pow(s->radius,2)*(zSup - zInf +(pow(zInf,3)-pow(zSup,3))/(3*pow(s->radius,2)));
+
+				phiAverage[numLayer]+=volPart;
+				velAverage[numLayer]+=volPart*b->state->vel[0];
+				dragAverage[numLayer]+=volPart*fDrag[0];
+			}
+			numLayer+=1;
+		}
+	}
+	//Normalized the weighted velocity by the volume of particles contained inside the cell
+	for(int n=0;n<nMax;n++){
+		if (phiAverage[n]!=0){
+			velAverage[n]/=phiAverage[n];
+			dragAverage[n]/=phiAverage[n];
+			//Normalize the concentration after
+			phiAverage[n]/=vCell;
+		}
+		else {
+			velAverage[n] = 0.0;
+			dragAverage[n] = 0.0;
+		}
+	}
+	//Assign the results to the global/public variables of HydroForceEngine
+	phiPart = phiAverage;
+	vxPart = velAverage;
+	averageDrag = dragAverage;
+
+	//desactivate the average to avoid calculating at each step, only when asked by the user
+	activateAverage=false; 
+}
+
+

=== modified file 'pkg/common/ForceEngine.hpp'
--- pkg/common/ForceEngine.hpp	2014-12-12 17:01:09 +0000
+++ pkg/common/ForceEngine.hpp	2015-01-27 15:52:46 +0000
@@ -71,20 +71,26 @@
 
 
 class HydroForceEngine: public PartialEngine{
+	private:
+		void averageProfile();
 	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"))
+	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.densFluid>`), $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. This engine also evaluate the average particle velocity, solid volume fraction and drag force depth profiles. 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, and C++ code is faster than python.",
+		((Real,densFluid,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"))
+                ((bool,lift,false,,"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"))
+		((Real,vCell,,,"Volume of averaging cell"))
+		((int,nCell,,,"number of cell in the height"))
+                ((Vector3r,gravity,Vector3r(0,0,-9.81),,"Gravity vector (may depend on the slope)."))
+		((vector<Real>,vxFluid,,,"Discretized streamwise fluid velocity depth profile"))
+		((vector<Real>,phiPart,,,"Discretized solid volume fraction depth profile"))
+		((vector<Real>,vxPart,,,"Discretized solid velocity depth profile"))
+		((vector<Real>,averageDrag,,,"Discretized drag depth profile"))
+		((bool,activateAverage,false,,"If true, activate the calculation of the average depth profiles of drag, solid volume fraction, and solid velocity for the application of the force (phiPart in hindrance function) and to use in python for the coupling with the fluid."))
 		((bool,velFluct,false,,"If true, activate the determination of turbulent fluid velocity fluctuation for the next time step only at the position of each particle, using a simple discrete random walk model based on the Reynolds stresses :yref:`turbStress<HydroForceEngine.squaredAverageTurbfluct>`"))
 		((vector<Real>,vFluctX,,,"Vector associating a x fluid velocity fluctuation to each particle. Fluctuation calculated in the C++ code"))
 		((vector<Real>,vFluctZ,,,"Vector associating a Z fluid velocity fluctuation to each particle. Fluctuation calculated in the C++ code"))