/** ![Antique cars often use a so-called [thermosiphon](https://en.wikipedia.org/wiki/Thermosiphon) instead of a water pump. Hence the large radiators. Image of a Ford Model T via [Wikipedia](https://en.wikipedia.org/wiki/Antique_car)](https://upload.wikimedia.org/wikipedia/commons/thumb/0/00/1916_Ford_Model_T_touring_car.JPG/250px-1916_Ford_Model_T_touring_car.JPG) # 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](tube/mov.mp4) 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. ~~~gnuplot Radial profiles of the azimuthal velocity. 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); }