sandbox/Antoonvh/tube.c

    Antique cars often use a so-called thermosiphon instead of a water pump. Hence the large radiators. Image of a Ford Model T via Wikipedia

    Antique cars often use a so-called thermosiphon instead of a water pump. Hence the large radiators. Image of a Ford Model T via Wikipedia

    Thermal circulation

    We model a circulation in a torus with major radius R_t and minor radius r_t. The tube is diferentially heated so that a buoyancy-driven flow starts.

    #include "grid/octree.h"
#include "embed.h"
#include "navier-stokes/centered.h"
#include "tracer.h"
#include "diffusion.h"
#include "view.h"
#define BVIEW 1
#include "particles.h"

scalar b[], *tracers = {b};
face vector av[], muc[];

    The torus is #defined using these variables and macros:

    double minor_rt = 1., major_rt = 3.;

#define R      (sqrt(sq(x) + sq(y)))
#define SINPSI (y/R)
#define COSPSI (x/R)
#define TUBE   (R <= 0.01 ? major_rt				\
		: (sqrt(sq(x - major_rt*COSPSI) +		\
			sq(y - major_rt*SINPSI) + sq(z))))

    The donut is a no-slip wall.

    u.n[embed] = dirichlet (0);
u.t[embed] = dirichlet (0);
#if (dimension == 3)
u.r[embed] = dirichlet (0);
#endif
b[embed]   = dirichlet (x);

int main() {
  L0 = 10;
  X0 = Y0 = Z0 = -L0/2.;
  N = 64;
  TOLERANCE = 1.e-4;
  mu = muc;
  a = av;
  run();
}

    The geometry is initialized and tracer particles are seeded.

    event init (t = 0) {
  vertex scalar phi[];
  foreach_vertex() 
    phi[] = minor_rt - TUBE;
  fractions (phi, cs, fs);
  
  n_part = 400;
  int j = 0;
  loc = malloc (n_part*sizeof(coord));
  do {
    double a = noise();
    double b = noise();
    if (sq(a) + sq(b) < sq(0.8*minor_rt)) {
      loc[j].x = a + major_rt;
      loc[j].y = 0;
      loc[j].z = b;
    } 
  } while (j++ < n_part);
}

event properties (i++) {   //diffusivity 
  foreach_face()
    muc.x[] = fm.x[]/200.;
  boundary ((scalar*){muc});
}

event acceleration (i++) { //Gravity
  foreach_face(y)
    av.y[] = (b[] + b[0,-1])/2.;
}

event tracer_diffusion (i++) 
  diffusion (b, dt, mu);

    The maximum resolution corresponds to a 64^3 cell equidistant grid.

    event adapt (i++) 
  adapt_wavelet ((scalar*) {cs, b, u},
		 (double[]){0.03, 0.4, 0.3, 0.3, 0.3}, 6);

event stop (t = 50.);

    Results

    We output a movie, displaying the buoyancy field, the velocity field, particles and a poor reconstruction of a part of the embedded boundary.

    Thermal convection in a tube

    The particles reveal an initial, short-lived recurring (i.e. compensating) flow in the centre of the tube. This is when a shallow ring of warm fluid rises arround the neutral fluid in the center. Furthermore, an attractive secondary circulation is visible!

    Naturally, we wish to study the evolution of the azimutal velocity (u_{\theta}), averaged over the minor-radial coordinate (r). The results are plotted below.

    set xr [0:1.1]
set yr [-0.1:1.5]
set xlabel 'minor radius'
set ylabel 'u_{theta}'
set size ratio 1
set grid
plot 'prof1' u 1:2 w l lw 2 t 't = 1',		\
     'prof3' u 1:2 w l lw 2 t 't = 3' ,		\
     'prof5' u 1:2 w l lw 2 t 't = 5' ,		\
     'prof10' u 1:2 w l lw 2 t 't = 10' ,	\
     'prof20' u 1:2 w l lw 2 t 't = 20' ,	\
     'prof35' u 1:2 w l lw 2 t 't = 35' ,	\
     'prof50' u 1:2 w l lw 2 t 't = 50'
    Radial profiles of the azimuthal velocity. (script)

    Radial profiles of the azimuthal velocity. (script)

    event movie (t += 0.1) {
#if (dimension == 3)
  scalar U[], t[], bd[];
  face vector tf[];
  double slice_a = -1;
  foreach_cell() {
    bd[] = b[];
    t[] = cs[];
    if (x > slice_a)
      t[] = 0;
    if (cm[] < 1e-6) {
      U[] = nodata;
      bd[] = nodata;
    }
    else {
      U[] = 0;
      foreach_dimension()
	U[] += sq(u.x[]);
      U[] = sqrt(U[]);
    }
  }
  foreach_face() {
    tf.x[] = fm.x[];
    if (x > slice_a)
      tf.x[] = 0;
  }
  boundary_flux ({tf});
  view (phi = 0.1, psi = 0., theta = 0.8);
  cells();
  squares ("bd");
  squares ("U", n = {1,0,0}, alpha = slice_a + 0.1, min = 0);
  draw_vof ("t", "tf");
  scatter (loc, s = 40);
  box();
#else
  draw_vof ("cs", "fs", filled = -1, fc = {1,1,1});
  cells();
  squares ("b");
#endif
  save ("mov.mp4");
}

#include "profile6.h"
event profs (t += 1.) {
  char fname[99];
  sprintf (fname, "prof%g", t);
  scalar Uth[];
  foreach() 
    Uth[] = (u.y[] * COSPSI - u.x[]*SINPSI)*cm[];
  boundary ({Uth});
  vertex scalar phi[];
  foreach_vertex()
    phi[] = TUBE;
  profiles ({Uth}, phi, fname, 1.1*minor_rt, 0, 0.25);
}