helical_rings_A015.c

Quasi-steady solutions for a vortex filament in shaped like a torus unknot, \mathcal{U}_{p,1} for p=2,4,8 and A=0.15R
References

Quasi-steady solutions for a vortex filament in shaped like a torus unknot, \mathcal{U}_{p,1} for p=2,4,8 and A=0.15R

Vortex rings and vortex helices are well-known steady solutions of the Biot-Savart law. Vortex filaments in the shape of a torus knots and unknots may also have steady solutions. In this example, we consider a vortex filament in the shape of a torus unknot, \mathcal{U}_{p,1}, which is can be seen as a helical coil wrapped around a torus with a single poloidal revolution and p toroidal revolutions. We can also think of this curve as a thin vortex ring superimposed with a Kelvin wave of azimuthal wavenumber m=p and finite-amplitude. Here, we consider 3 cases with p=2, p=4 and p=8, and equal amplitude A=0.15R of the Kelvin wave, where R is the radius of the undeformed vortex ring. The shape of the vortex filaments are shown below.

Steady vortex filaments in the shape of a torus unknot, \mathcal{U}_{2,1}, \mathcal{U}_{4,1}, and \mathcal{U}_{8,1}

The vortex filaments are initialized from a file containing the coordinates of the filament centerline, the core size and the volume of each segment.

The vortices move as material lines in the fluid according to: \begin{aligned} \frac{d\vec{x}_c}{dt} = \vec{U}_{ind} + \vec{U}_{\infty} \end{aligned} where \vec{x}_c is the position vector of vortex filament, \vec{U}_{ind} the velocity induced by the vortex filament and \vec{U}_{\infty} an external velocity field. The induced velocity \vec{U}_{ind} is given by the Biot-Savart law.

In this example, the three vortex filaments are independent but are shown together for visualization purposes. The vortex filaments rotate with a constant angular velocity around the vertical axis and translate with a constant velocity along the vertical axis. The motion of the vortex filaments is shown below.

#include "grid/octree.h"
#include "run.h"
#include "view.h"
#include "acastillo/input_fields/filaments.h"
#include "acastillo/input_fields/draw_filaments.h"
#include "acastillo/filaments/biot-savart.h"

In this example, we use a as the characteristic length scale and rescale all other lengths by 1/a. This ways, the velocity induced by a circular vortex ring of radius R and circulation \Gamma=2\pi remains of order unity.

#define SCALE (2e-2)

int nseg ;
double R = 1.0/SCALE;
double a = 2e-2;
double dtmax = 0.02;
double tend = 6000.0;
struct vortex_filament helical_ring1;
struct vortex_filament helical_ring2;
struct vortex_filament helical_ring3;

coord Uinfty = {0., 0., 0.};

The main time loop is defined in run.h.

int minlevel = 4;
int maxlevel = 4;
vector u[];
int main(){
  L0 = 9.0/SCALE;
  X0 = Y0 = Z0 = -L0 / 2;
  
  N = 1 << minlevel;  
  init_grid(N);  
  run();
}

Initial conditions

We read the vortex filaments from a file. The files must have the following format:

nseg
a
phi_0  x_0  y_0  z_0  dl_0
phi_1  x_1  y_1  z_1  dl_1
...
phi_(nseg-1)  x_(nseg-1)  y_(nseg-1)  z_(nseg-1)  dl_(nseg-1)

where nseg is the number of segments, a the core size (assumed constant), phi_i the azimuthal angle of the segment, (x_i,y_i,z_i) the coordinates of the center of the segment and dl_i the length of the segment. Viewed from the top, the vortex filaments look like planar curves, but are actually 3D curves.

import numpy as np
import matplotlib.pyplot as plt
	
ax = plt.figure()
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5), sharey=True, sharex=True)
data = np.loadtxt('./helical_rings_A015.filament1.txt', delimiter=' ', usecols=[0,1,2,3], skiprows=2)
ax1.plot(data[:,1], data[:,2], label=r'$p=2$')
ax2.plot(data[:,3], data[:,2], label=r'$p=2$')

data = np.loadtxt('./helical_rings_A015.filament2.txt', delimiter=' ', usecols=[0,1,2,3], skiprows=2)
ax1.plot(data[:,1], data[:,2], label=r'$p=4$')
ax2.plot(data[:,3], data[:,2], label=r'$p=4$')

data = np.loadtxt('./helical_rings_A015.filament3.txt', delimiter=' ', usecols=[0,1,2,3], skiprows=2);
ax1.plot(data[:,1], data[:,2], label=r'$p=8$');
ax2.plot(data[:,3], data[:,2], label=r'$p=8$');

phi = np.linspace(0,2*np.pi)
ax1.plot(np.cos(phi), np.sin(phi), '--', color='grey', lw=0.5)

ax1.legend(loc='best')
ax1.axis('image')
ax1.set_xlabel(r'Coordinate $x$')
ax1.set_ylabel(r'Coordinate $y$')

ax2.axis('image')
ax2.set_xlabel(r'Coordinate $z$')
plt.tight_layout()
plt.savefig('plot_init.svg')

Evolution of the axial coordinate of the filaments (script)

void read_filament_from_file(const char *filename, struct vortex_filament *filament, double R=1.0, double xshift=0.0, double yshift=0.0, double zshift=0.0, double scale=1.0) {
  // Debug: Print current working directory and list files
  fprintf(stderr, "Attempting to open: %s\n", filename);
  fprintf(stderr, "Current working directory contents:\n");
  system("pwd > log2 2>&1");
  system("ls -la >> log2 2>&1");

  FILE *fp = fopen(filename, "r");
  if (!fp) {
    fprintf(stderr, "Error: could not open %s\n", filename);
    fprintf(stderr, "Errno: %s\n", strerror(errno));
    exit(1);
  }  

  int nseg;
  double a;
  fscanf(fp, "%d", &nseg);
  fscanf(fp, "%lf", &a);
  
  double phi[nseg];
  double a1[nseg];
  double vol1[nseg];
  coord C1[nseg];

  for (int i = 0; i < nseg; i++) {
    double x, y, z, dl;
    if (fscanf(fp, "%lf %lf %lf %lf %lf", &phi[i], &x, &y, &z, &dl) != 5) {
      fprintf(stderr, "Error: invalid format in %s at line %d\n", filename, i+2);
      exit(1);
    }
    C1[i] = (coord){x * scale + xshift, y * scale + yshift, z * scale + zshift};
    a1[i] = a * scale;
    vol1[i] = pi * sq(a) * dl * cube(scale);
  }
  fclose(fp);

  allocate_vortex_filament_members(filament, nseg);
  memcpy(filament->phi, phi, nseg * sizeof(double));
  memcpy(filament->C, C1, nseg * sizeof(coord));
  memcpy(filament->a, a1, nseg * sizeof(double));  
  memcpy(filament->vol, vol1, nseg * sizeof(double));
}

event init (t = 0) {

  read_filament_from_file("./helical_rings_A015.filament1.txt", &helical_ring1, R, xshift=-L0/3, zshift=-L0/2, scale=1/SCALE);
  read_filament_from_file("./helical_rings_A015.filament2.txt", &helical_ring2, R, xshift= 0.0,  zshift=-L0/2, scale=1/SCALE);
  read_filament_from_file("./helical_rings_A015.filament3.txt", &helical_ring3, R, xshift= L0/3, zshift=-L0/2, scale=1/SCALE);

  local_induced_velocity(&helical_ring1);  
  local_induced_velocity(&helical_ring2);
  local_induced_velocity(&helical_ring3);
  
  view (camera="front", fov=6, width=2560, height=800);
  draw_tube_along_curve(helical_ring1.nseg, helical_ring1.C, helical_ring1.a, scale=1);
  draw_tube_along_curve(helical_ring2.nseg, helical_ring2.C, helical_ring2.a, scale=1);
  draw_tube_along_curve(helical_ring3.nseg, helical_ring3.C, helical_ring3.a, scale=1);
  save ("prescribed_curve1.png");
  clear();
  
  view (camera="top", fov=25, width=1920, height=1920);
  draw_tube_along_curve(helical_ring1.nseg, helical_ring1.C, helical_ring1.a, scale=1);
  draw_tube_along_curve(helical_ring2.nseg, helical_ring2.C, helical_ring2.a, scale=1);
  draw_tube_along_curve(helical_ring3.nseg, helical_ring3.C, helical_ring3.a, scale=1);
  save ("prescribed_curve2.png");
  clear();
}

We release the vortex filaments at the end of the simulation.

event finalize(t = end){
  free_vortex_filament_members(&helical_ring1);
  free_vortex_filament_members(&helical_ring2);
  free_vortex_filament_members(&helical_ring3);
}

A Simple Time Advancing Scheme

Time integration is done in two steps. First, we evaluate the (self-)induced velocity at \vec{x}_c

event evaluate_velocity (i++) {

  memcpy(helical_ring1.Uprev, helical_ring1.U, helical_ring1.nseg * sizeof(coord));
  local_induced_velocity(&helical_ring1);
  for (int j = 0; j < helical_ring1.nseg; j++) {
    helical_ring1.Uauto[j] = nonlocal_induced_velocity(helical_ring1.C[j], &helical_ring1);
    foreach_dimension(){
      helical_ring1.U[j].x = helical_ring1.Uauto[j].x + helical_ring1.Ulocal[j].x - Uinfty.x;      
    }
  }

  memcpy(helical_ring2.Uprev, helical_ring2.U, helical_ring2.nseg * sizeof(coord));
  local_induced_velocity(&helical_ring2);
  for (int j = 0; j < helical_ring2.nseg; j++) {
    helical_ring2.Uauto[j] = nonlocal_induced_velocity(helical_ring2.C[j], &helical_ring2);
    foreach_dimension(){
      helical_ring2.U[j].x = helical_ring2.Uauto[j].x + helical_ring2.Ulocal[j].x - Uinfty.x;      
    }
  }

  memcpy(helical_ring3.Uprev, helical_ring3.U, helical_ring3.nseg * sizeof(coord));
  local_induced_velocity(&helical_ring3);
  for (int j = 0; j < helical_ring3.nseg; j++) {
    helical_ring3.Uauto[j] = nonlocal_induced_velocity(helical_ring3.C[j], &helical_ring3);
    foreach_dimension(){
      helical_ring3.U[j].x = helical_ring3.Uauto[j].x + helical_ring3.Ulocal[j].x - Uinfty.x;      
    }
  }
}

Then, we advect the filament segments using an explicit Adams-Bashforth scheme

event advance_filaments (i++) {      
  for (int j = 0; j < helical_ring1.nseg; j++) {
    foreach_dimension(){
      helical_ring1.C[j].x += dt*(3.*helical_ring1.U[j].x - helical_ring1.Uprev[j].x)/2.;
    }
  } 

  for (int j = 0; j < helical_ring2.nseg; j++) {
    foreach_dimension(){
      helical_ring2.C[j].x += dt*(3.*helical_ring2.U[j].x - helical_ring2.Uprev[j].x)/2.;
    }
  }
  
  for (int j = 0; j < helical_ring3.nseg; j++) {
    foreach_dimension(){
      helical_ring3.C[j].x += dt*(3.*helical_ring3.U[j].x - helical_ring3.Uprev[j].x)/2.;
    }
  }
}

Finally, we advance by one time-step

event advance_filaments (i++, last) {      
  dt = dtnext (dtmax);
}

Outputs

Viewed from the top, we can see the elliptical ring seems to be rotating with a constant angular velocity around the vertical axis and translating with a constant velocity along the vertical axis. However, these velocities depend on the amplitude of the Kelvin wave.

import numpy as np
import matplotlib.pyplot as plt

a = 2e-2
ax = plt.figure()
n = 64*2
data = np.loadtxt('curve1.txt', delimiter=' ', usecols=[0,2,3,4])
data[:,1] = data[:,1] + 9/a/3
r = np.sqrt(data[:,1]**2 + data[:,2]**2)
plt.plot(data[::n,3]*a, r[::n]*a, label='$p=2$')

n = 64*4
data = np.loadtxt('curve2.txt', delimiter=' ', usecols=[0,2,3,4])
r = np.sqrt(data[:,1]**2 + data[:,2]**2)
plt.plot(data[::n,3]*a, r[::n]*a, label='$p=4$')

n = 64*8
data = np.loadtxt('curve3.txt', delimiter=' ', usecols=[0,2,3,4])
data[:,1] = data[:,1] - 9/a/3
r = np.sqrt(data[:,1]**2 + data[:,2]**2)
plt.plot(data[::n,3]*a, r[::n]*a, label='$p=8$')

plt.xlabel(r'Coordinate $z/R$')
plt.ylabel(r'Coordinate $r/R$')
plt.legend()
plt.tight_layout()
plt.savefig('plot_xy_vs_t.svg')

Evolution of the x- and y-coordinates of the first filament segment (script)

event movie (t += 10.00){
  view (camera="front", fov=6, width=2560, height=800);
  draw_tube_along_curve(helical_ring1.nseg, helical_ring1.C, helical_ring1.a, scale=1);
  draw_tube_along_curve(helical_ring2.nseg, helical_ring2.C, helical_ring2.a, scale=1);
  draw_tube_along_curve(helical_ring3.nseg, helical_ring3.C, helical_ring3.a, scale=1);
  save("top_view.mp4");

  view (camera="top", fov=25, width=1920, height=1920);
  draw_tube_along_curve(helical_ring1.nseg, helical_ring1.C, helical_ring1.a, scale=1);
  draw_tube_along_curve(helical_ring2.nseg, helical_ring2.C, helical_ring2.a, scale=1);
  draw_tube_along_curve(helical_ring3.nseg, helical_ring3.C, helical_ring3.a, scale=1);
  save("side_view.mp4");
}

event final (t += 10.00, t <= tend){
  if (pid() == 0){
    FILE *fp;
    fp = fopen("curve1.txt", "a"); 
    write_filament_state(fp, &helical_ring1); 
    fclose(fp);

    fp = fopen("curve2.txt", "a"); 
    write_filament_state(fp, &helical_ring2); 
    fclose(fp);

    fp = fopen("curve3.txt", "a"); 
    write_filament_state(fp, &helical_ring3); 
    fclose(fp);
  }
}

References

[castillo2022]	A. Castillo-Castellanos and S. Le Dizès. Closely spaced co-rotating helical vortices: long-wave instability. Journal of Fluid Mechanics, 946:A10, 2022. [ DOI ]
[castillo2021]	A. Castillo-Castellanos, S. Le Dizès, and E. Durán Venegas. Closely spaced corotating helical vortices: General solutions. Phys. Rev. Fluids, 6:114701, Nov 2021. [ DOI \| http ]
[duran2019]	E. Durán Venegas and S. Le Dizès. Generalized helical vortex pairs. Journal of Fluid Mechanics, 865:523–545, 2019. [ DOI ]