Hey,
I am kind of new to LBM. Currently I am working on a code for single disc particle sedimentation in a stationary fluid.
Gravity to the left. Right end is constant pressure boundary condition. Left, top bottom is wall.
g=0.001 in lattice units
The body force value need not be added to the fluid in the domain, right? Only the particle acceleration is update by this g.
I am using the momentum exchange algorithm by Ladd to simulate it. However, the velocity of the particle just keeps fluctuating between 0.0019 to 0.0023.(lattice units)
These values aren’t equal to the terminal velocity I determined by using Cd(Drag Coeffecient) given by Sucker and Brauer (1975).
I am stuck with this for a month now. I have searched but not been able to find a way. Please help.
Just for clarification the following is how I get the force:
if(isSolid(node[j][i-1],centroid))
{
xb=(node[j][i].x + c[3][0]/2.) - centroid.x;
yb=(node[j][i].y + c[3][1]/2.) - centroid.y;
ub=centroid.u - (wyb);
vb=centroid.v + (wxb);
fx=2*(node[j][i].f[3] - (6.*wt[3]*node[j][i].rho)*(c[3][0]*ub + c[3][1]*vb) - node[j][i-1].f[1])*c[3][0];
fy=2*(node[j][i].f[3] - (6.*wt[3]*node[j][i].rho)*(c[3][0]*ub + c[3][1]*vb) - node[j][i-1].f[1])*c[3][1];
Fx+=fx;
Fy+=fy;
T+=( (fy*xb) - (fx*yb) );
node[j][i].f[1]=node[j][i].f[3] - (6.*wt[3]*node[j][i].rho)*(c[3][0]*ub + c[3][1]*vb);
node[j][i-1].f[3]=node[j][i-1].f[1] + (6.*wt[3]*node[j][i].rho)*(c[3][0]*ub + c[3][1]*vb);
}
Similar code is for all the other neighbouring nodes of a particular fluid boundary node
The velocity update for the disc is done in the following snippet:
void updatingSolid(props ¢roid,int solid_count,long double Fx,long double Fy,long double T,long double &w,long double del)
{
long double k,M,I;
k=sqrt( 0.5*pow(a,2) );
M=(solid_rho)*PI*pow((a),2);
I=M*pow(k,2);
//Displacing
centroid.x+=( centroid.u + 0.5*(Fx/M - g) );
centroid.y+=( centroid.v + 0.5*(Fy/M) );
//Updating velocities
centroid.u+=(Fx/M - g);
centroid.v+=(Fy/M);
w+=(T/I);
}
Wherever it is written node[j] please replace with node(j)(i]).
Indexing is j increases downwards, i increases rightwards.
However, X is rightwards, Y is upwards.
fx, fy are forces at a particular node j,i along a particular direction of d2q9 lattice
F is total force, T is total torque, w= angular velocity, a is radius of disc, k is radius of gryration of disc, g is gravity.
xb,yb are coordinates of boundary node(halfway between solid boundary node and fluid boundary node) with respect to the disc centroid.
ub,vb are the corresponding velocities of boundary node
c[k][0] is lattice velocity of kth d2q9 direction along X.
c[k][1] is lattice velocity of kth d2q9 direction along Y.
I hope this clarifies the code. Please help with the code because quantitatively my results are wrong. It is very urgent.
The body force value need not be added to the
fluid in the domain, right? Only the particle
acceleration is update by this g.
Wrong! To obtain correct sedimentation behaviour you have to apply a net force to the particle which is equal to the difference of the particles gravity and buoyancy. The force that you measure after a terminal velocity is reached is the drag force.
In your code, I miss a loop going over all the lattice directions for the force calculation. Am I wrong or is there really only direction 3 that you include.
This is just a snippet of the code that checks for nodes where c[3] collides with the particle and bounces bacck…The check for solid node runs for all directions.
Also, I have included buoyancy while updating the solid particle’s velocity.
What I meant to ask was if any gravity force density(rho*g) term is to be added in the collision eqn like so:
f +=[ 1/tau(feq - f) + (gravity force term) ]
I have seen a case where body force due to free convection is added into collision term in a similar fashion as above.
Now I get your point! Including a body force (for gravitation) in the LB scheme would lead to a pressure-gradient in your domain. You could do that instead of adding a constant buoyancy force on the particle.
However, in doing so there are at least 2 additional points to consider:
- You have to preinitialize your lattice with the correct pressure (i.e. density) gradient. Otherwise you get shock waves.
- Because a pressure gradient is a density gradient in LBM you have to watch out for compressibility errors. (See also Buick & Greated paper on gravity in LB).
Hope that helps
Hi, you can refer to the paper : [Wen et. al. PRE 85, 016704(2012)].