# src/examples/atomisation.c

# Atomisation of a pulsed liquid jet

A dense cylindrical liquid jet is injected into a stagnant lighter phase (density ratio 1/27.84). The inflow velocity is modulated sinusoidally to promote the growth of primary shear instabilities. Surface tension is included and ultimately controls the characteristic scale of the smallest droplets.

We solve the two-phase Navier–Stokes equations with surface tension. We need the *tag()* function to count the number of droplets.

```
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
#include "tag.h"
```

We define the radius of the jet, the initial jet length, the Reynolds number and the surface tension coefficient.

```
#define radius 1./12.
#define length 0.025
#define Re 5800
#define SIGMA 3e-5
```

The default maximum level of refinement is 10 and the error threshold on velocity is 0.1.

```
int maxlevel = 10;
double uemax = 0.1;
```

To impose boundary conditions on a disk we use an auxilliary volume fraction field *f0* which is one inside the cylinder and zero outside. We then set an oscillating inflow velocity on the left-hand-side and free outflow on the right-hand-side.

```
scalar f0[];
u.n[left] = dirichlet(f0[]*(1. + 0.05*sin (10.*2.*π*t)));
u.t[left] = dirichlet(0);
p[left] = neumann(0);
f[left] = f0[];
u.n[right] = neumann(0);
p[right] = dirichlet(0);
```

The program can take two optional command-line arguments: the maximum level and the error threshold on velocity.

```
int main (int argc, char * argv[])
{
if (argc > 1)
maxlevel = atoi (argv[1]);
if (argc > 2)
uemax = atof (argv[2]);
```

The initial domain is discretised with ${64}^{3}$ grid points. We set the origin and domain size.

```
init_grid (64);
origin (0, -1.5, -1.5);
size (3.);
```

We set the density and viscosity of each phase as well as the surface tension coefficient and start the simulation.

```
rho1 = 1., rho2 = 1./27.84;
mu1 = 2.*radius/Re*rho1, mu2 = 2.*radius/Re*rho2;
f.σ = SIGMA;
run();
}
```

## Initial conditions

```
event init (t = 0) {
if (!restore (file = "dump")) {
```

We use a static refinement down to *maxlevel* in a cylinder 1.2 times longer than the initial jet and twice the radius.

` refine (x < 1.2*length && sq(y) + sq(z) < 2.*sq(radius) && level < maxlevel);`

We initialise the auxilliary volume fraction field for a cylinder of constant radius.

```
fraction (f0, sq(radius) - sq(y) - sq(z));
f0.refine = f0.prolongation = fraction_refine;
restriction ({f0}); // for boundary conditions on levels
```

We then use this to define the initial jet and its velocity.

```
foreach() {
f[] = f0[]*(x < length);
u.x[] = f[];
}
boundary ({f,u.x});
}
}
```

## Outputs

We log some statistics on the solver.

```
event logfile (i++) {
if (i == 0)
fprintf (ferr,
"t dt mgp.i mgpf.i mgu.i grid->tn perf.t perf.speed\n");
fprintf (ferr, "%g %g %d %d %d %ld %g %g\n",
t, dt, mgp.i, mgpf.i, mgu.i,
grid->tn, perf.t, perf.speed);
}
```

If gfsview is installed on the system, we can also visualise the simulation as it proceeds.

```
#if 0
event gfsview (i += 10) {
static FILE * fp = popen ("gfsview3D atomisation.gfv", "w");
output_gfs (fp, list = {u,p,f});
}
#endif
```

To generate an animation, uncomment this and use the `run.sh`

script below.

```
#if 0
event movie (t += 2e-2)
{
scalar pid[], ff[];
foreach() {
pid[] = fmod(pid()*(npe() + 37), npe());
ff[] = f[] < 1e-4 ? 0 : f[] > 1. - 1e-4 ? 1. : f[];
}
boundary ({pid,ff});
output_gfs (fout);
fprintf (fout, "Save stdout { format = PPM width = 1600 height = 1200 }\n");
}
#endif
```

We save snapshots of the simulation at regular intervals to post-process with gfsview.

`event snapshot (t = 0.1; t += 0.1; t <= 2.5) {`

We also dump the entire simulation to be able to restart.

```
dump (file = "dump");
char name[80];
sprintf (name, "snapshot-%g.gfs", t);
scalar pid[], ff[];
foreach() {
pid[] = fmod(pid()*(npe() + 37), npe());
ff[] = f[] < 1e-4 ? 0 : f[] > 1. - 1e-4 ? 1. : f[];
}
boundary ({pid,ff});
output_gfs (file = name);
}
```

## Counting droplets

The number and sizes of droplets generated by the atomising jet is a useful statistics for atomisation problems. This is not a quantity which is trivial to compute. The *tag()* function is designed to solve this problem. Any connected region for which *f[] > 1e-3* (i.e. a droplet) will be identified by a unique “tag” value between 0 and *n-1*.

```
event droplets (t += 0.1)
{
scalar m[];
foreach()
m[] = f[] > 1e-3;
int n = tag (m);
```

Once each cell is tagged with a unique droplet index, we can easily compute the volume *v* and position *b* of each droplet. Note that we use *foreach_leaf()* rather than *foreach()* to avoid doing a parallel traversal when using OpenMP. This is because we don’t have reduction operations for the *v* and *b* arrays (yet).

```
double v[n];
coord b[n];
for (int j = 0; j < n; j++)
v[j] = b[j].x = b[j].y = b[j].z = 0.;
foreach_leaf()
if (m[] > 0) {
int j = m[] - 1;
v[j] += dv()*f[];
coord p = {x,y,z};
foreach_dimension()
b[j].x += dv()*f[]*p.x;
}
```

When using MPI we need to perform a global reduction to get the volumes and positions of droplets which span multiple processes.

```
#if _MPI
MPI_Allreduce (MPI_IN_PLACE, v, n, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce (MPI_IN_PLACE, b, 3*n, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#endif
```

Finally we output the volume and position of each droplet to standard output.

```
for (int j = 0; j < n; j++)
fprintf (fout, "%d %g %d %g %g %g\n", i, t,
j, v[j], b[j].x/v[j], b[j].y/v[j]);
fflush (fout);
}
```

## Mesh adaptation

We adapt the mesh according to the error on the volume fraction field and the velocity.

```
event adapt (i++) {
adapt_wavelet ({f,u}, (double[]){0.01,uemax,uemax,uemax}, maxlevel);
}
```

## Running in parallel

To run in 3D on occigen, we can do on the local machine

```
local% qcc -source -grid=octree -D_MPI=1 atomisation.c
local% scp _atomisation.c occigen.cines.fr:
```

and on occigen (to run on 64*24 = 1536 cores, with 12 levels of refinement)

`occigen% sbatch --nodes=64 --time=5:00:00 run.sh`

with the following `run.sh`

script

```
#!/bin/bash
#SBATCH -J basilisk
#SBATCH --nodes=1
#SBATCH --constraint=HSW24
#SBATCH --ntasks-per-node=24
#SBATCH --threads-per-core=1
#SBATCH --time=10:00
#SBATCH --output basilisk.output
#SBATCH --exclusive
LEVEL=12
module purge
module load openmpi
module load intel
NAME=atomisation
mpicc -Wall -std=c99 -O2 _$NAME.c -o $NAME -lm
srun --mpi=pmi2 -K1 --resv-ports -n $SLURM_NTASKS ./$NAME $LEVEL \
2> log-$LEVEL-$SLURM_NTASKS | \
sed 's/VariableTracerVOF f/VariableTracerVOF ff/g' | \
$HOME/local/bin/gfsview-batch3D $NAME.gfv \
> out-$LEVEL-$SLURM_NTASKS.ppm
```

Note that this assumes that `gfsview-batch`

is installed on the machine to generate the animation. The animated GIF below can then be created using

`ppm2gif --scale 0.5 --delay 10 < out-12-1536.ppm > atomisation.gif`

## Results

The gain in number of grid points of the adaptive simulation (relative to a constant 4096^{3} Cartesian grid discretisation) is illustrated below as well as the computational speed (in points.timesteps/sec).