modifying fix_viscous

Submitted by Somnath on Wed, 10/05/2011 - 02:58

Hi,

I would like to model a viscous drag force (stokes drag) for a polydispersed simulation where the particle size distribution is specified using fix_particletemplate/sphere and fix_particledistribution/discrete. Since, all the particles have similar properties, I would also like to have one atom_type in the simulations.

However, as it is right now, fix_viscous is defined in force/velocity units and to scale the viscous force for different particles (with different radii), different atom_types must be used. This becomes cumbersome while defining material properties for all the different atom_type combinations, especially when the material properties are the same. Is there a way around this?

Or, is it possible to modify fix_viscous such that the viscosity of the fluid is specified instead of gamma (=force/velocity) and the drag force is calculated as, F_drag=viscosity*radius[i]*v[i]

Has anyone tried this? I would be grateful, if someone can help me with the code.

Regards
Som

Forums:

ckloss's picture

ckloss | Wed, 10/05/2011 - 07:59

Sure, that is possible. You can implement it like

double *radius = atom->radius;
double *v = atom->v;
F_drag=viscosity*radius[i]*v[i]

Cheers, Christoph

Somnath | Wed, 10/05/2011 - 20:54

Hey Christoph,

Thanks for your reply. Just wanted to make sure that you meant modifying fix_viscous like this:


void FixViscous::post_force(int vflag)
{
// apply drag force to atoms in group
// direction is opposed to velocity vector
// magnitude depends on atom type

double **v = atom->v;
double **f = atom->f;
double *radius = atom->radius
int *mask = atom->mask;
int *type = atom->type;
int nlocal = atom->nlocal;

double drag;

for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
drag = gamma[type[i]];
f[i][0] -= drag*radius[i]*v[i][0];
f[i][1] -= drag*radius[i]*v[i][1];
f[i][2] -= drag*radius[i]*v[i][2];
}
}

Regards
Som

Matteo | Mon, 06/06/2016 - 11:42

Hello,

Probably this thread is a bit old, but I think it helped me in clarifying my doubt.

I struggled a bit to understand what LIGGGHTS wants as gamma value when using "fix viscous".

Can you please confirm that I should set this value just computing the product 3*pi*eta*d? This would mean that LIGGGHTS does not take in consideration, during the computation, the effect of the particle diameter for the calculation of the viscous force. And, in case I use a Gaussian distribution to describe the particle sizes, this would introduce some approximations for the calculation of the viscous force.
In case this is true, I guess I have to modify the fix_viscous as Somnath described above?

Thanks

Best regards
Matteo

Matteo | Wed, 06/15/2016 - 16:36

Hi Christoph,

Thanks for your answer.
I studied C at high school but I have no idea what the symbol -> means. IS that part of C++?
Looking at the code, anyway, it seems that the radius of the particles is not considered, as probably contained into the gamma variable, as guessed from the previous posts/manual.

It would be useful if in the next release of LIGGGHTS this would be changed.
At this point I guess I have to modify the code and compile it again?
Should I copy the code suggested by Somnath?

Thanks.

Best regards
Matteo

ckloss's picture

ckloss | Tue, 06/28/2016 - 17:41

Hi Matteo,

it's correct the implementation of this fix is quite basic, not including particle radii. If you don't like it, you can change it as mentioned above (or ask someone with C++ experience to do it for you)

Best wishes
Christoph