sandbox/cselcuk/balistic-sphere.c
Balistic test case, a sphere with initial transversal velocity and subject to gravity
# include "grid/octree.h"
# define DLM_Moving_particle 1
# define NPARTICLES 1
# define DLM_alpha_coupling 1
Parameters related to the spatial resolution
# define LEVEL 4
# define adaptive 1
# if adaptive
# define MAXLEVEL (LEVEL + 8)
# endif
Physical parameters
The characteristic velocity scale is U_c = \sqrt(gD) with D the sphere’s diameters (characteristic length scale).The Reynolds number based on U_c is Re = 25.
# define rhoval 1000. // fluid density
# define rhosolid 1140. // particle's density
# define gravity_y -9.81 // gravity acceleration
# define Lc (0.01) // characteristic length - particle diameter
# define tval 0.125284 // fluid dynamic viscosity, here Re = 25
# define Uc 0.31321 // sqrt(grav*Lc)
# define Tc (Lc/Uc) // characteristic time
/* output and numerical parameters */
# define mydt (Tc/50) //time-step
# define maxtime (100*mydt) // simulation time
Fictitious domain implementation with toy model granular solver
# include "dlmfd-toygs.h"
# include "navier-stokes/perfs.h"
double deltau;
scalar un[];
int main () {
origin (0., 0., 0.);
L0 = 120*Lc;
/* set time step */
DT = mydt;
/* initialise a uniform grid */
init_grid (1 << LEVEL);
/* boundary conditions */
u.t[left] = dirichlet(0.); //v
u.r[left] = dirichlet(0.); //w
u.n[left] = dirichlet(0.); //u
uf.t[left] = dirichlet(0.); //v
uf.r[left] = dirichlet(0.); //w
uf.n[left] = dirichlet(0.); //u
u.t[right] = dirichlet(0.); //v
u.r[right] = dirichlet(0.); //w
u.n[right] = dirichlet(0.); //u
uf.t[right] = dirichlet(0.); //v
uf.r[right] = dirichlet(0.); //w
uf.n[right] = dirichlet(0.); //u
u.n[bottom] = dirichlet(0.); //v
u.t[bottom] = dirichlet(0.); //w
u.r[bottom] = dirichlet(0.); //u
uf.n[bottom] = dirichlet(0.); //v
uf.t[bottom] = dirichlet(0.); //w
uf.r[bottom] = dirichlet(0.); //u
u.n[top] = dirichlet(0.); //v
u.t[top] = dirichlet(0.); //w
u.r[top] = dirichlet(0.); //u
uf.n[top] = dirichlet(0.); //v
uf.t[top] = dirichlet(0.); //w
uf.r[top] = dirichlet(0.); //u
u.n[front] = dirichlet(0.); //w
u.t[front] = dirichlet(0.); //u
u.r[front] = dirichlet(0.); //v
uf.n[front] = dirichlet(0.); //w
uf.t[front] = dirichlet(0.); //u
uf.r[front] = dirichlet(0.); //v
u.n[back] = dirichlet(0.); //w
u.t[back] = dirichlet(0.); //u
u.r[back] = dirichlet(0.); //v
uf.n[back] = dirichlet(0.); //w
uf.t[back] = dirichlet(0.); //u
uf.r[back] = dirichlet(0.); //v
// Convergence criteria for the N-S solver
TOLERANCE = 1e-4;
NITERMIN = 2;
run();
}
event init (i = 0) {
particle * p = particles;
/* Initial condition for the fluid */
foreach() {
foreach_dimension() {
u.x[] = 0.;
}
un[] = u.x[];
}
double diam = Lc;
for (int k = 0; k < NPARTICLES; k++) {
/* initial condition: particle's position */
GeomParameter gp = {0};
gp.center.x = 20*diam;
gp.center.y = L0-20*diam;
gp.center.z = L0/2;
gp.radius = diam/2.;
p[k].g = gp;
/* initial condition: particle's velocities */
coord c = {0., 0., 0.};
c.x = 1.; p[k].U = c;
c.x = 0.; p[k].w = c;
}
}
event output_data (t += Tc; t <= maxtime)
{
dump(file = "dump");
particle * p = particles;
dump_particle(p);
}
event last_output (t=end) {
p.nodump = false;
dump("dump");
}
Results
set grid
show grid
Lc = 0.01
Uc = 0.31321
tc = Lc/Uc
set xlabel "t/t_c"
set ylabel "t_c*\omega"
plot "particle-data-0" u ($1/tc):($8*tc) w l title "t_c\omega^p_x", "particle-data-0" u ($1/tc):($9*tc) w l title "\omega^p_y t_c", "particle-data-0" u ($1/tc):($10*tc) w l title "\omega^p_z t_c"
reset
set grid
show grid
Lc = 0.01
Uc = 0.31321
tc = Lc/Uc
set xlabel "t/t_c"
set ylabel "u/Uc"
plot "particle-data-0" u ($1/tc):($5/Uc) w l title "U_x/Uc","particle-data-0" u ($1/tc):($6/Uc) w l title "U_y/Uc", "particle-data-0" u ($1/tc):($7/Uc) w l title "U_z/Uc"
reset
diam = 0.01
show grid
set xlabel "x/diam"
set ylabel "y/diam"
plot "particle-data-0" u ($2/diam):($3/diam) w l title "trajectory"