yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #13058
[Branch ~yade-pkg/yade/git-trunk] Rev 4030: extended version of ForceEngine with its own fluid solver (by J. Chauchat)
------------------------------------------------------------
revno: 4030
committer: bchareyre <bruno.chareyre@xxxxxxxxxxxxxxx>
timestamp: Wed 2017-04-05 19:18:06 +0200
message:
extended version of ForceEngine with its own fluid solver (by J. Chauchat)
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 2017-01-23 16:38:29 +0000
+++ pkg/common/ForceEngine.cpp 2017-04-05 17:18:06 +0000
@@ -82,7 +82,8 @@
}
}
-void HydroForceEngine::action(){
+
+void HydroForceEngine::action(){
/* Application of hydrodynamical forces */
if (activateAverage==true) averageProfile(); //Calculate the average solid profiles
@@ -98,6 +99,8 @@
Vector3r liftForce = Vector3r::Zero();
Vector3r dragForce = Vector3r::Zero();
Vector3r convAccForce = Vector3r::Zero();
+ //deterministic version
+// Vector3r vRel = Vector3r(vxFluid[p],0,0) - 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) {
@@ -455,4 +458,360 @@
}
}
-
+/* function declaration */
+void doubleq(double ddam1[],double ddam2[],double ddam3[],double ddbm[],double ddxm[],int n);
+void calbeta(int irheolf, double alphas[], double beta[], double alphasmax, unsigned long n);
+void calviscotlm(int iturbu, int ilm, double dz, double h,double ufn[], double alphas[], double alphasmax, double kappa, double lmExp, double viscoft[], unsigned long n);
+void fluidModel(double h,double sig[],double dsig[],double dp,double ufn[],double alphaf[],double rhof,double viscof,double usnp[],double alphas[],double rhos,double alphasmax,double dpdx,double slope, double gra, double tfin,double dt,double cfdpYade[], double ufnp[], double viscoft[]);
+
+void HydroForceEngine::updateVelocity() {
+ fluidModel(fluidHeight,&sig[0],&dsig[0],diameterPart,&vxFluid[0],
+ &phiPart[0],//FIXME: in the py script it's in fact (1-phiPart) which is passed to nsmp, wtf? and why the redundancy wrt [#] (below)
+ densFluid,
+ viscoDyn /*is it really the dynamic one here? else pass viscoDyn/densFluid (awkward anyway) */,
+ &vxPart[0],
+ &phiPartFluid[0] /*[#] WHY PASSING IT AGAIN?!! It seems to be equal to phiPart in practice */,
+ densPart, alphasmax, dpdx, slope,
+ gravity[2],/*FIXME: I'm assuming that gravity is along 2-axis (??), does it mean that users have to define gravity multiple times in one single script? Newton::gravity, HydroForceEngine::gravity? ugly... :-\ */
+ fluidResolPeriod,dtFluid,
+ &taufsi[0],
+ &vxFluid[0],&turbulentViscosity[0]); //<-------- output of the function
+}
+
+
+void fluidModel(double h,double sig[],double dsig[],double dp,double ufn[],double alphaf[],double rhof,double viscof,double usnp[],double alphas[],double rhos,double alphasmax,double dpdx,double slope, double gra, double tfin,double dt,double cfdpYade[], double ufnp[], double viscoft[])
+{
+ const unsigned& ndimz=HydroForceEngine::ndimz;
+ unsigned j;
+ int irheolf,idrag,iturbu,ilm,iusl,idrift;
+ double dummy,dn1,ds1,dn2,ds2,dz,dzn,dzs,dzm;
+ double alphafp,alphafwn,alphafws,rhop;
+ double kappa,lmExp,expoRZ;
+ double as,an,ap1,a[ndimz],b[ndimz],c[ndimz],s[ndimz],udrift[ndimz],beta[ndimz],cfdp[ndimz];
+ double Rep,cd;
+ double time=0;
+
+ // impose constant grid size
+ dz=dsig[0]*h;
+
+ //
+ // Option of the code
+ //
+ // irheolf = 0 : Viscosite du fluide pur
+ // 1 : Viscosite d'Einstein
+ // 2 : Viscosite de Graham
+ // 3 : Viscosite de Krieger-Dougherty / Ishii-Zuber
+ // 4 : Viscosite de Boyer et al.
+ irheolf = 0;
+
+ // idrag = 0 : Trainee de Dallavale
+ // 1 : Trainee de Schiller & Naumann
+ // 2 : Trainee de Clift & Gauvin
+ // 3 : Loi de Stokes
+ // 4 : Loi de Darcy
+ // 5 : Loi de Ergun
+ // 6 : Imposed from averaged drag force (Yade)
+ idrag = 0;
+
+ // exposant de Richardson-Zaki (fonction d'entravement)
+ expoRZ = -3.1;
+
+ // iturbu = 0 : Pas de turbulence
+ // 1 : Longueur de mélange
+ // ilm = 0 : longueur de melange de Prandtl
+ // 1 : longueur de melange de Prandtl avec effet de Surface libre
+ // 2 : Li and Sawamoto (1995)
+ iturbu = 1;
+ ilm = 2;
+
+ kappa = 0.41;
+ lmExp = 1;
+
+ // iusl = 0 : Condition de Dirichlet (u=0 en z=h)
+ // 1 : Condition de Neumann (dudz=0 en z=h)
+ iusl = 1;
+
+ // idrift = 0 : Sans vitesse de dispersion
+ // 1 : Avec vitesse de dispersion
+ idrift = 0;
+
+
+ // Time loop
+ while (time < tfin)
+ {
+ // Advance time
+ time = time + dt;
+ printf("t=%7.4f s\n",time);
+
+ // Viscosity amplification factor
+ calbeta(irheolf,alphas,beta,alphasmax,ndimz);
+
+ // Eddy viscosity
+ calviscotlm(iturbu,ilm,dz,h,ufn,alphas,alphasmax,kappa,lmExp,viscoft,ndimz);
+
+ // Bottom boundary condition: (always no-slip)
+ j=0;
+ a[j]=0.;
+ b[j]=1.;
+ c[j]=0.;
+
+ s[j]=0.;
+
+ // Top boundary condition: (0: no-slip / 1: zero gradient)
+ j=ndimz-1;
+ if (iusl==0)
+ {
+ a[j]=0.;
+ b[j]=1.;
+ c[j]=0.;
+ }
+ else if (iusl==1)
+ {
+ a[j]=-1.;
+ b[j]=1.;
+ c[j]=0.;
+ }
+ s[j]=0.;
+
+
+ //Main loop
+ for(j=1;j<=ndimz-2;j++)
+ {
+ // volume fraction interpolation (staggered grid)
+ alphafp = 0.5*(alphaf[j]+alphaf[j+1]);
+ if (j==1)
+ {
+ dzn = dz;
+ dzs = 1.5*dz;
+
+ alphafwn=0.5*(alphaf[j+2]+alphaf[j+1]);
+ alphafws=alphaf[j-1];
+ }
+ else if (j==ndimz-2)
+ {
+ dzn = 0.5*dz;
+ dzs = dz;
+
+ alphafwn=alphaf[j+1];
+ alphafws=0.5*(alphaf[j ]+alphaf[j-1]);
+ }
+ else
+ {
+ dzn = dz;
+ dzs = dz;
+
+ alphafwn=0.5*(alphaf[j+2]+alphaf[j+1]);
+ alphafws=0.5*(alphaf[j ]+alphaf[j-1]);
+ }
+ dzm = 0.5*(dzn+dzs);
+
+ // Drag model
+ if (idrag==6)
+ {
+ // read from file:
+ //cfdp[j] = cfdpYade[j];
+ cfdp[j] = max(0.,cfdpYade[j]/max(fabs(usnp[j]-ufn[j]),1e-5))/rhof;
+ //printf("%u\t%12.8f\t%12.8f\t%12.8f\n", j,ufn[j],usnp[j],cfdp[j]);
+ }
+ else if (idrag==0)
+ {
+ // Dallavalle + Richardson & Zaki
+ Rep=max(fabs(ufn[j]-usnp[j])*dp/viscof,1e-8);
+ cd=(24.4/Rep+0.4)*pow(alphafp,expoRZ);
+ cfdp[j]=0.75*(1-alphafp)/dp*cd*fabs(ufn[j]-usnp[j]);
+ //printf("%u\t%12.8f\t%12.8f\t%12.8f\n", j,ufn[j],usnp[j],cfdp[j]);
+ }
+ else
+ {
+ printf("idrag value undefined");
+ break;
+ }
+ // Diffusion coefficients
+ // Eddy viscosity terms
+ dummy=dt/dzm;
+ ds1=dummy*viscoft[j ]/dzn;
+ dn1=dummy*viscoft[j+1]/dzs;
+ // Viscous terms
+ ds2=dummy*viscof*beta[j ]/dz*alphafp;
+ dn2=dummy*viscof*beta[j+1]/dz*alphafp;
+
+ // Numerical scheme coefficient (diffussion only in 1DV)
+ an=dn1+(dn2)*alphafwn;
+ as=ds1+(ds2)*alphafws;
+
+ ap1=dn1+(dn2)*alphafp+ds1+(ds2)*alphafp;
+
+ // LHS: algebraic system coefficients
+ a[j] = - as;
+ b[j] = alphafp + ap1 + dt*cfdp[j];
+ c[j] = - an;
+
+ // RHS: unsteady, gravity, drag, pressure gradient
+ s[j]= alphafp*ufn[j] + alphafp*gra*sin(slope)*dt + dt*cfdp[j]*(usnp[j]+udrift[j]) - alphafp*dpdx/rhof*dt;
+ //printf("%u\t%12.8f\t%12.8f\t%12.8f\n", j,a[j],b[j],c[j]);
+ }
+ // Implicit solution using tridiag (useful because of potential very high viscosities)
+ doubleq( a, b , c, s , ufnp,ndimz);
+ // Update solution for next time step
+ for(j=0;j<=ndimz-1;j++)
+ {
+ ufn[j]=ufnp[j];
+ //printf("%u\t%12.8f\t%12.8f\t%12.8f\n", j,ufn[j],usnp[j],cfdp[j]);
+ }
+ }
+}
+
+void calbeta(int irheolf, double alphas[], double beta[], double alphasmax, unsigned long n)
+{
+ const unsigned& ndimz=HydroForceEngine::ndimz;
+ int j;
+ double ratio1,hsuram1;
+
+ // viscosity amplification factor
+ if (irheolf==0)
+ // 0 : Viscosite du fluide pur
+ {
+ for(j=0;j<=ndimz;j++)
+ beta[j]=1.;
+ }
+ else if (irheolf==1)
+ // 1 : Viscosite d'Einstein
+ {
+ for(j=0;j<=ndimz;j++)
+ beta[j]=1.+2.5*alphas[j];
+ }
+ else if (irheolf==2)
+ // 2 : Viscosite de Graham // this one is buged
+ {
+ for(j=0;j<=ndimz;j++)
+ {
+ ratio1=pow(min(alphas[j]/alphasmax,0.99),0.3333333);
+ hsuram1=0.5*max(ratio1/(1.-ratio1),1e-3);
+ // printf("%u\t%7.4f\t%7.4f\n",j,ratio1,hsura);
+ beta[j]=1.+2.5*alphas[j]+2.25*(1+0.5/hsuram1)*(hsuram1-pow(1.+1./hsuram1,-1)-pow(1+1./hsuram1,2));
+ }
+ }
+ else if (irheolf==3)
+ // 3 : Viscosite de Krieger-Dougherty / Ishii-Zuber
+ {
+ for(j=0;j<=ndimz;j++)
+ beta[j]=pow(1.-min(alphas[j]/alphasmax,0.99),-(2.5*alphasmax));
+ }
+ else if (irheolf==4)
+ // 4 : Viscosite de Boyer et al.
+ {
+ for(j=0;j<=ndimz;j++)
+ beta[j]=1. + 2.5*alphas[j]*pow(1.-min(alphas[j]/alphasmax,0.99),-1);
+ }
+}
+
+// ------------------------------------------------------------------------------//
+
+void calviscotlm(int iturbu, int ilm, double dz, double h,double ufn[], double alphas[], double alphasmax, double kappa, double lmExp, double viscoft[], unsigned long n)
+{
+ const unsigned& ndimz=HydroForceEngine::ndimz;
+ int j;
+ double lm[n],dist,sum,alphasmid,dzm,dudz;
+
+ // Eddy viscosity
+ if (iturbu==0)
+ // 0 : No turbulence
+ {
+ for(j=0;j<=ndimz;j++)
+ viscoft[j]=0.;
+ }
+ else if (iturbu==1)
+ // 1 : Turbulence activated
+ {
+ if (ilm==0)
+ // 0 : Prandtl mixing length
+ {
+ dist = 0;
+ lm[0]=0.;
+ for(j=1;j<=ndimz-1;j++)
+ {
+ lm[j]=kappa*dist;
+ dist = dist + dz;
+ }
+ }
+ else if (ilm==1)
+ // 1 : Parabolic profile (free surface flows)
+ {
+ dist = 0;
+ lm[0]=0.;
+ for(j=1;j<=ndimz-1;j++)
+ {
+ lm[j]=kappa*dist*sqrt(1.-dist/h);
+ dist = dist + dz;
+ }
+ lm[ndimz-1]=0.;
+ }
+ else if (ilm==2)
+ // 2 : Li and Sawamoto (1995) integral of concentration profile
+ {
+ dist = 0;
+ sum = 0.;
+ lm[0]=0.;
+ for(j=1;j<=ndimz-1;j++)
+ {
+ alphasmid=0.5*(alphas[j-1]+alphas[j]);
+ sum = sum + min(alphasmid/alphasmax,0.9999999)*dz;
+ lm[j]=kappa*(dist-pow(sum,lmExp));
+ // printf("%u\t%7.4f\t%7.4f\n",j,lm[j],dist*kappa);
+ dist=dist+dz;
+ }
+ lm[ndimz-1]=lm[ndimz-2];
+ }
+ // Compute the velocity gradient and the mixing length
+ for(j=1;j<ndimz-1;j++)
+ {
+ if (j==1)
+ dzm=1.5*dz;
+ else if (j==ndimz-1)
+ dzm=0.5*dz;
+ else
+ dzm=dz;
+ dudz=(ufn[j]-ufn[j-1])/dz;
+ viscoft[j]=(1.-alphas[j])*pow(lm[j],2)*fabs(dudz);
+ //printf("%u\t%7.4f\t%7.4f\t%7.4f\n",j,fabs(dudz),lm[j],viscoft[j]);
+ }
+ viscoft[ndimz-1]=viscoft[ndimz-2];
+ }
+}
+
+// ------------------------------------------------------------------------------//
+void doubleq(double ddam1[],double ddam2[],double ddam3[],double ddbm[],double ddxm[],int n)
+{
+ /* reosolution of tridiagonal system
+ am1,2,3[n] - Tridiagonal matrix coefficients
+ bm[n] - RHS vector
+ xm[n] - Solution vector
+ em[n], fm[n) - working arrays
+ n - algebraic system size
+ */
+
+ int i,ii;
+ double ddem[n],ddfm[n],dddiv;
+
+ // downward sweep
+ dddiv=ddam2[0];
+ ddem[0]=-ddam3[0]/dddiv;
+ ddfm[0]=ddbm[0]/dddiv;
+
+ for (i=1;i<=n-2;i++)
+ {
+ dddiv=ddam2[i]+ddam1[i]*ddem[i-1];
+ ddem[i]=-ddam3[i]/dddiv;
+ ddfm[i]=(ddbm[i]-ddam1[i]*ddfm[i-1])/dddiv;
+ }
+ // upward sweep
+ dddiv=ddam2[n-1]+ddam1[n-1]*ddem[n-2];
+ ddfm[n-1]=(ddbm[n-1]-ddam1[n-1]*ddfm[n-2])/dddiv;
+ ddxm[n-1]=ddfm[n-1];
+
+ for (ii=1;ii<=n-1;ii++)
+ {
+ i=n-1-ii;
+ ddxm[i]=ddem[i]*ddxm[i+1]+ddfm[i];
+ }
+}
\ No newline at end of file
=== modified file 'pkg/common/ForceEngine.hpp'
--- pkg/common/ForceEngine.hpp 2017-01-23 16:38:29 +0000
+++ pkg/common/ForceEngine.hpp 2017-04-05 17:18:06 +0000
@@ -1,6 +1,6 @@
// 2004 © Janek Kozicki <cosurgi@xxxxxxxxxx>
// 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
-// 2014 © Raphael Maurin <raphael.maurin@xxxxxxxxx>
+
#pragma once
@@ -72,13 +72,30 @@
class HydroForceEngine: public PartialEngine{
public:
+ static const unsigned ndimz = 900;
void averageProfile();
void turbulentFluctuation();
void turbulentFluctuationBIS();
void turbulentFluctuationFluidizedBed();
+ void updateVelocity();
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.",
+ /// BEGIN Transitory variables for migration fortran->c++
+ ((Real,fluidHeight,1.,,"Height of the flow from the bottom of the sample"))
+ ((vector<Real>,sig,vector<Real>(ndimz),,"???????"))
+ ((vector<Real>,dsig,vector<Real>(ndimz),,"???????"))
+ ((Real,diameterPart,0.,,"Reference particle diameter"))
+ ((Real,densPart,1.,,"mass density of the particles"))
+ ((Real,dpdx,0.,,"pressure gradient along streamwise direction"))
+ ((Real,slope,0.,,"Angle of the slope in radian"))
+ ((Real,fluidResolPeriod,0.,,"1/fluidResolPeriod = frequency of the calculation of the fluid profile"))
+ ((vector<Real>,taufsi,vector<Real>(ndimz),,"Create Taufsi/rhof = dragTerm/(rhof(vf-vxp)) to transmit to the fluid code"))
+ ((Real,dtFluid,0.,,"Time step for the fluid resolution"))
+ ((vector<Real>,turbulentViscosity,vector<Real>(ndimz),,"Turbulent viscocity"))
+ ((vector<Real>,phiPartFluid,vector<Real>(ndimz),,"???"))
+ ((Real,alphasmax, 0.61,,"????"))
+ //// END
((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, i.e. to the ground."))
@@ -119,10 +136,10 @@
((vector<Real>,convAcc,0,,"Convective acceleration, depth dependent"))
((Real,convAccOption,false,,"To activate the convective acceleration"))
((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("updateVelocity",&HydroForceEngine::updateVelocity,"Calculate the fluid velocity profile.")
.def("turbulentFluctuation",&HydroForceEngine::turbulentFluctuation,"Apply turbulent fluctuation to the problem.")
.def("turbulentFluctuationBIS",&HydroForceEngine::turbulentFluctuationBIS,"Apply turbulent fluctuation to the problem with an alternative formulation.")
.def("turbulentFluctuationFluidizedBed",&HydroForceEngine::turbulentFluctuationFluidizedBed,"Apply turbulent fluctuation to the problem with another alternative formulation.")
