sandbox/piersonj/turbulence2D.c
(Naive) Forced turbulence in a 2D-periodic box with non-coalescing droplets
We investigate the dynamics of droplets in a two-dimensional periodic box. To induce turbulence, we employ a velocity linearly dependent forcing term, similar to the approach described in http://basilisk.fr/src/examples/isotropic.c. It should be noted that the intention is not to specifically generate two-dimensional isotropic turbulence with this forcing term, but rather to provide a crude way to make the droplets collide. To prevent numerical coalescence of droplets, we utilize the no-coalescence module.
#include "grid/multigrid.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
#include "no-coalescence.h"
We need to store the variable forcing term.
face vector av[];
int maxlevel = 8;
int main (int argc, char * argv[])
{
//maxruntime (&argc, argv);
if (argc > 1)
maxlevel = atoi(argv[1]);
The domain is (2\pi)^3 and doubly-periodic.
L0 = 2.*pi;
foreach_dimension()
periodic (right);
The acceleration is defined below. The level of refinement is maxlevel.
double rhoInt = 1.;
double muInt = 0.01;
rho1 = rhoInt, mu1 = muInt;
rho2 = rhoInt, mu2 = muInt;
f.sigma = 4.;
a = av;
N = 1 << maxlevel;
run();
}
Initial positions of the drop
double geometry (double x, double y) {
double inout=0.,i1=0.,i2=0.,i3=0.,i4=0.,Circle1=0.,Circle2=0.,Circle3=0.,Circle4=0.;
Circle1 = -sq(x-pi/2)-sq(y-pi/2)+sq(1);
Circle2 = -sq(x-pi/2)-sq(y-3.5*pi/2)+sq(1);
Circle3 = -sq(x-3.5*pi/2)-sq(y-3*pi/2)+sq(1);
Circle4 = -sq(x-3*pi/2)-sq(y-pi/2)+sq(1);
if(Circle1 > 0.) i1=1;
if(Circle2 > 0.) i2=1;
if(Circle3 > 0.) i3=1;
if(Circle4 > 0.) i4=1;
return(i1+i2+i3+i4);
}
Initial conditions
The initial condition is a periodic flow with white noise.
event init (i = 0) {
if (!restore (file = "restart"))
fraction(f, geometry(x,y));
foreach() {
u.x[] = (1-f[])*cos(pi*y) + noise()*0.1;
u.y[] = (1-f[])*sin(pi*x) + noise()*0.1;
}
}
Linear forcing term
We compute the average velocity and add the corresponding linear forcing term.
event acceleration (i++) {
coord ubar;
foreach_dimension() {
stats s = statsf(u.x);
ubar.x = s.sum/s.volume;
}
foreach_face()
av.x[] += 0.1*(1-f[])*((u.x[] + u.x[-1])/2. - ubar.x);
}
event logfile(i++) {
printf ("i = %d t = %g\n", i,t);
fflush(stdout);
}
We generate movies of the vortices, volume fraction and velocity (along x).
event output (t += 0.025; t <= 15.)
{
scalar omega[];
vorticity (u, omega);
output_ppm (omega, min = -1, max = 1, file = "omega.mp4");
output_ppm (f, file = "fraction.mp4");
output_ppm (u.x, file = "u.mp4");
}