sandbox/Antoonvh/vraxi.c
A lab experiment by T.T. Lim and T.B. Nickels reveals the spectecular dynamics of a collision between two vortex rings. Animated .gif via imgur.
The collision of two vortex rings
Inspired by the experimental results of Lim and Nickels (1992), and the more recent hype generated by this youtube movie (why cannot I get 2.5 Million downloads on my scintific produce?), a study on the collision of two vortex rings is presented. The lab set-up aims to create an axisymmetric flow configuration where two toroidal vortices collide head-on. Different from a planar vortex collision, during the toroidal collision, the rings expand at a remarkably fast rate, after which they become unstable and seem to break-up in to smaller, perpendicularly-alligned vortex rings. On this page the initial (axisymmetric) phase of the experiment is studied. This is useful before one redirects attention towards a fully 3D flow simulation. More specifically, here the effect of the fluid’s viscosity (\nu) on the collision dynamics is studied.
Set-up
A vortex ring is characterized by a circular patch of vorticiy with (minor) radius r that revolves arround an axis of symmetry with major ring radius R (R>r). In the set-up a Gaussian distribution is used for the initial vorticity,
\displaystyle \omega = \omega_0 e^{-\frac{d^2}{\sigma^2}},
where d is the distance from a circle with distance R away from the symmtery axis. Introducing a time scale, \mathcal{T}=\omega_0^{-1}, and a length scale, \mathcal{L} = \sigma. In our set-up we use normalized values for \omega_0 and \sigma. This also warrants the introduction of a dimensionless parameter regarding the shape of the torus;
\displaystyle \Pi = \frac{R}{\sigma}=10.
One can express a velocity scale, U =\sigma \omega _0. This enables the expression of the relative importance of the fluid’s viscosity using a well-know non-dimensional number, the so-called Reynolds number,
\displaystyle Re = \frac{UR}{\nu}= \Pi \frac{\sigma^2 \omega_0}{\nu},
which we vary to assess the influence of the fluid’s viscosity on the dynamics. More specifically, we are interested in what role Re plays with respect to the flow-acceleration phase of the experiment by solving the relevant equations of motion under the axisymmetric approximation.
#include "axi.h"
#include "navier-stokes/centered.h"We set R = 10\sigma, and the initial separation distance between the two vortex-ring centres is dx = 20\sigma. The domain is a periodic channel of square size, 75\sigma \times 75\sigma.
double rp = 10, xp = 10., omg0 = 10, vis = 0.;
double xc[61], yc[61];
int j, l;
double tend = 25;
double deltat = 0.5;
int main(){
periodic(left);
L0 = 75.;
X0 = -L0/2.;
Y0 = 0.;
foreach_dimension()
u.x.refine = refine_linear;
init_grid(256);By iteratively increasing the viscosity, simulations are run with four different Reynolds numbers; Re = \{ 10000,\ 1000,\ 100,\ 10 \}.
vis = 0.0001;
for (j = 0; j < 4; j++){
vis *= 10.;
run();
}
}
event init(t = 0){
DT = 0.01;
const face vector muc[] = {vis, vis};
mu = muc;
l = 0;Initialization
We use an \omega \rightarrow \psi \rightarrow \mathbf{u} formulation to find the velocity components corresponding to the described vorticity distribution.
scalar omg[], psi[];
foreach_cell() // A zero guess prooved to be important.
psi[] = 0.;
psi[bottom] = dirichlet(0.); //anti-symmetry for psi.
TOLERANCE = 10E-6;
for (int g = 0; g < 10; g++){
foreach()
omg[] = omg0*(exp(-(sq(y - rp) + sq(x - xp))) - exp(-(sq(y - rp) + sq(x + xp))));
boundary({omg, psi});
poisson(psi, omg);
boundary({psi});
foreach(){
u.x[] = ((psi[0,1] - psi[0,-1])/(2*Delta));
u.y[] = -(psi[1,0] - psi[-1,0])/(2*Delta);
}
boundary((scalar*){u});
adapt_wavelet((scalar*){u}, (double[]){0.005, 0.005}, 11);
}
TOLERANCE = 10E-4;
}Adaptation
We use grid adaptivity. Since we do not know too much about the to-be-simulated dynamics, we only set a refinement criterion for the velocity components \zeta_u = U/20, and do not limit the maximum level of refinement.
event adapt(i++)
adapt_wavelet((scalar*){u}, (double[]){0.05 ,0.05}, 19);Output
The evolution of the location of the positive vortex patch patch is monitored along with the kinetic energy. Furthermore, two movies are produced for each run, to display both the grid structre and to the reveal vorticity dynamics.
event output (t += deltat){
scalar omg[], lev[];
vorticity(u, omg);
double xx = 0, yy = 0, ens = 0 , e = 0.;
int n = 0;
foreach(){
e += M_PI*sq(Delta)*y*(sq(u.x[]) + sq(u.y[]));
n++;
lev[] = level;
if (omg[] > 0){
xx += sq(omg[])*x*sq(Delta);
yy += sq(omg[])*y*sq(Delta);
ens += sq(omg[])*sq(Delta);
}
}
xc[l] = xx/ens;
yc[l] = yy/ens;
l++;
char fname[99];
sprintf(fname,"dataj=%d",j);
static FILE * fpw = fopen(fname, "w");
fprintf(fpw, "%d\t%g\t%g\t%d\t%g\n", i, t, e, depth(), (double)(1<<(depth()*dimension))/(double)n);
printf("%g\t%g\t%d\n", t, e, depth());
sprintf(fname,"ppm2mp4 vort%d.mp4", j);
static FILE * fp = popen(fname, "w");
output_ppm(omg, fp, n = 256, map = cool_warm);
sprintf(fname,"ppm2mp4 level%d.mp4", j);
static FILE * fpl = popen(fname, "w");
output_ppm(lev, fpl, n = 512, min = 2, max = depth());
}
event stop(t = tend){
char fname[99];
sprintf(fname, "velo%d", j);
FILE * fpv = fopen(fname, "w");
for (int m = 1; m < 50; m++){
double v = sqrt ((sq(xc[m+1] - xc[m-1]) + sq(yc[m+1] - yc[m-1])))/(2*deltat);
fprintf(fpv, "%g\t%g\t%g\t%g\n", 0.05*(double)m, v, xc[m], yc[m]);
}
return 1;
}Results
First, the vorticity dynamics are revealed for the highest and lowest Re values:
Re = 10000
Re = 10
The used grid structure for the highest-Reynolds-number case:
Grid for Re = 10000, More red colors correspond to a higher level of refinement
We have traced the vortex paths of the red-coloured vortex, so we plot them:
set terminal pngcairo enhanced font ',12'
set xlabel 'z / {/Symbol s}'
set ylabel 'r / {/Symbol s}'
set size square
set key top left
plot 'velo0' u 3:4 w l lw 3 t 'Re = 10000' , \
'velo1' u 3:4 w l lw 3 t 'Re = 1000' , \
'velo2' u 3:4 w l lw 3 t 'Re = 100' , \
'velo3' u 3:4 w l lw 3 t 'Re = 10' 
(script)
Next, The evolution of the vortex-centre velocity is shown.
set terminal pngcairo enhanced font ',12'
set xlabel 't{/Symbol s}'
set ylabel 'Velocity / U'
set size square
set key top left
plot 'velo0' u 1:2 w l lw 3 t 'Re = 10000' , \
'velo1' u 1:2 w l lw 3 t 'Re = 1000' , \
'velo2' u 1:2 w l lw 3 t 'Re = 100' , \
'velo3' u 1:2 w l lw 3 t 'Re = 10' 
(script)
It appears that the severity and duration of the acceleration phase is a function of the Reynoldsnumber. Also note that the initial translation speed is of the order U. Meaning that the experiment with Re=1000, apart from the 3D effects, approximately corresponds to the experiment of Lim and Nickels (1992).
We do a sanity check of our numerics and validate that the kinetic energy is monotomically decreasing:
set ylabel 'Kin. En.'
set key top right
plot 'dataj=0' u 2:3 w l lw 3 t 'Re = 10000' , \
'dataj=1' u 2:3 w l lw 3 t 'Re = 1000' , \
'dataj=2' u 2:3 w l lw 3 t 'Re = 100' , \
'dataj=3' u 2:3 w l lw 3 t 'Re = 10'
(script)
It appears that Re = 10000 is not resolved well enough and would require a lower error threshold. The other results do pass this particular check.
Finally, we check if the adaptive gridding approach made some sense by plotting the reduction factor of the number of grid cells compared to an equidistant grid apprach at the same depth of refinment. Noting that in our approach, the depth is diagnosed by the grid-adaptation algorithm. A non-trivial feature that the equidistant pre-tuned grid typically does not have.
set xlabel 'Iteration'
set ylabel 'Reduction factor'
set key bottom right
plot 'dataj=0' u 1:5 w l lw 3 t 'Re = 10000' , \
'dataj=1' u 1:5 w l lw 3 t 'Re = 1000' , \
'dataj=2' u 1:5 w l lw 3 t 'Re = 100' , \
'dataj=3' u 1:5 w l lw 3 t 'Re = 10'
(script)
We conclude that it is wise to use grid adaptivity if we extend the analysis to a full 3D experiment.
Reference
Lim, T. T.; Nickels, T. B. Instability and reconnection in the head-on collision of two vortex rings. Nature, 1992, 357.6375: 225.
