sandbox/Antoonvh/tube.c
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 #define
d 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.
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'
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);
}