# 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

# 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++)
(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 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)

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;
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);
}