sandbox/Antoonvh/stokes-particles.h
Stokes’ drag law for spherical particles
The acceleration of a small particle \mathbf{a}_p in a fluid medium can be described by,
\displaystyle \mathbf{a}_p = \frac{\mathbf{F}_p}{m_p}
where F_p are the forces acting on the particle and m_p is the mass of the particle. For spherical particles, we consider a drag force due the the velocity difference of the particle and the fluid medium.
\displaystyle \mathbf{F}_d = C_d r_p \mu\left(\mathbf{u}_f - \mathbf{v}_p\right),
where C_d = 6\pi, the Stokes drag coefficient, r_p the particle radius , \mu is the fluid medium’s dynamic viscosity and \left(\mathbf{u}_f - \mathbf{v}_p\right) is the volocity difference between the undisturbed fluid flow and the particle. Furthermore, the acceleration of gravity (G
) is considered to act on the particle.
\displaystyle \mathbf{F}_p = \mathbf{F}_d + \mathbf{F}_G,
\displaystyle \mathbf{F}_G = G \times V_p \left(\rho_p - \rho_f\right),
with V_p the volume of the particle. This gives a formulation for one-way coupled partile-laden flow. It can be casted in an inertial particle formalism.
Implementation
It can be interesting to compute “inertial-particles
” paths, alongside flow tracer-particles. Therfore, it is chosen to reuse the names of tracer-particle members for the mass density (u2.x
) and the radius (u2.y
). Note if tracer-particles.h are also used, that headerfile should be included first. Finally, if u2.y
is not set u2.z
can be used to set the relaxation time scale \tau_p. By default it is computed from using u2.y
as r_p
\displaystyle \tau = \frac{\rho_p 4r_p^2}{18\mu}
#ifndef ADD_PART_MEM
#define ADD_PART_MEM coord u; coord u2; long unsigned int tag;
#endif
#define PA_INP particle * pp; double dt;
#define INP (PA_inp){p(), &p(), 2.*dt}
#include "inertial-particles.h"
For two-way couling, the forces on the particles are stored in a field.
extern vector u; //Fluid medium flow
extern face vector mu; //Fluid medium dynamic viscosity
extern scalar rho; //Fluid medium density
coord G; //Gravity acceleration vector (should it be external) ?
scalar * _automatics_ = NULL;
coord p_acc (PA_inp inp) {
particle pa = inp.p;
double xp = pa.x; double yp = pa.y; double zp = pa.z;
Implicity time discretization for v_p
Since the flow-adjustment timescale \tau can be smaller than the solver time step, it makes sense to use an implicit equation for v_p,
\displaystyle \frac{v_p^{n+1}-v_p^n}{\mathtt{dt}} = \frac{1}{\tau}\left(u - v_p^{n+1}\right) + g',
which yields a weighted averaged between v_p^n and the ternimal velocity for v_p^{n+1}.
\displaystyle v_p^{n+1} =\frac{\mathtt{dt}^{-1}v_p^n + \tau^{-1}u + g'}{\mathtt{dt}^{-1} + \tau^{-1}}.
Rather than computing the acceleration, the velocity is directly updated in this function.
double muc = is_constant (mu.x) ? constant(mu.x) : interpolate (mu.x, xp, yp, zp);
double rp = pa.u2.y;
double rhop = pa.u2.x;
double tau = rp ? (muc > 0 ? rhop * sq(2*rp)/(18*muc) : HUGE) : pa.u2.z ? pa.u2.z : HUGE;
double itau = tau > 0 ? 1/tau : HUGE;
double idt = inp.dt > 0 ? 1./inp.dt: HUGE;
foreach_dimension()
inp.pp->u.x = (idt*pa.u.x + itau*interpolate(u.x, xp, yp ,zp))/(idt + itau);
// Gravity
double rhof = is_constant(rho) ? constant(rho) : interpolate (rho, xp, yp, zp);
foreach_dimension()
inp.pp->u.x += rhop ? (G.x*(rhop - rhof)/rhop)/(idt + itau) : 0;
// No additional acceleration
return (coord){0., 0., 0.};
}
#if TWO_WAY
event defaults (i = 0) {
if (is_constant(a.x)) {
a = new face vector;
foreach_face()
a.x[] = 0.;
boundary ((scalar *){a});
}
foreach_dimension()
_automatics_ = list_add (_automatics_ , u.x);
if (!is_constant (rho))
_automatics_ = list_add (_automatics_, rho);
if (!is_constant (mu.x)) {
foreach_dimension()
_automatics_ = list_add (_automatics_, mu.x);
}
}
void add_force (Particles P, vector F) {
particle_boundary (P);
double rhof = constant(rho);
boundary ((scalar*){u});
foreach_particle_in(P) {
double muc = is_constant (mu.x) ? constant(mu.x) : interpolate (mu.x, x, y, z);
double rp = p().u2.y;
double pref = 6*pi*muc*rp;
double xp = x, yp = y , zp = z;
Point point = locate (x,y,z);
double rhop = p().u2.x;
foreach_dimension()
F.x[] += pref*(p().u.x - interpolate(u.x, xp, yp, zp));
pref = 4./3.*pi*cube(rp)*(rhop-rhof);
}
}
event acceleration (i++) {
face vector av = a;
vector Fp[];
foreach()
foreach_dimension()
Fp.x[] = 0;
foreach_P_in_list (inertial_particles)
add_force (P, Fp);
double rhof = constant(rho);
boundary ((scalar*){Fp});
foreach_face()
av.x[] += face_value(Fp.x, 0)/(dv()*rhof);
}
#endif
Tests
Usage
Todo
Examplify with usagesTwo-way coupling with the centered solver.Test ifinertial-particles
can indeed co-exist withtracer-particles
- Critical tests