sandbox/Antoonvh/core.c

Liquid core convection in a planet.

A sphere with radius R contains a boussinesq fluid with viscosity \nu and reference density \rho_0. The fluid is subject to its own gravity potential \Phi and therefore has a quadratic variation with the radial coordinate (r), \Phi = A r^2, with A a constant (2\pi G\rho_0/dimension). For generaltiy we use a b (\langle b\rangle = t^{-2}) field to describe the effect of gravity:

\displaystyle \mathbf{F_g} = b \mathbf{r},

which eliminates \rho_0 as a system parameter.

The fluid is cooled at the sphere’s edge and thereby becomes slighly more dense. We use the (nagative) b-flux (B) to describe the cooling. The medium has a b-diffussivity \kappa.

With system parameters R, \nu, B,\kappa and the dimensionality we can construct the following dimensionless groups:

\displaystyle \Pr = \frac{\nu}{\kappa},

\displaystyle Re = \frac{\|B\|^{1/3}R^{8/3}}{\nu},

\displaystyle \Pi = \texttt{dimension}.

We choose Pr = 1 and Re = 30.

The 2D solution for the b field and some particle flow tracers

The 3D solution for the b field (on a hemisphere) and some particle flow tracers

Note that the earth has a solid core in a liquid mantle.

#include "embed.h"
#include
#include "tracer.h"
#include
//#include "view2.h"     //does not draw_vof facets for z > 0 in 3D
#include "view.h"
#define BVIEW 1
#include

#define RADIUS (sqrt(sq(x) + sq(y) + sq(z)))

scalar b[], * tracers = {b};
face vector muv[], av[];
double muf = 1./30., R = 1., B = -1.;

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

// Relaxation time (non physical)
double tau = 10.;

int main() {
L0 = R*2.1;
X0 = Y0 = Z0 = -L0/2.;
mu = muv;
a = av;
run();
}

event init (t = 0) {
scalar phi[];
foreach_vertex()
fractions (phi, cs, fs);
init_particles_2D_square_grid (6, 0, 0, R);
}

event properties (i++) {
foreach_face()
muv.x[] = fs.x[]*muf;
boundary ((scalar*){muv});
}

event acceleration (i++) {
foreach_face() {
coord f = {x, y, z};
av.x[] = face_value(b,0)*f.x;
}
}

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

We modify the buoyancy to omit large negative values. This “heating” does not influence the dynamics and should balance the cooling in a quasi steady state.

event no_run_away (i++) {
double tendency = dimension*B/R;
double mean = statsf(b).sum/statsf(b).volume;
foreach()
b[] -= (cs[] > 0)*dt*mean/tau;
boundary ({b});
printf ("%g %g cf.: %g\n",t, mean, tendency*tau);
}

#if TREE
adapt_wavelet ({b, u}, (double[]){0.4, 0.2, 0.2, 0.2}, 6, 4);
#endif

event movie (t += .1) {
view (fov = 19);
#if dimension == 2
scatter (loc);
draw_vof ("cs", "fs", filled = -1, fc = {1,1,1});
squares ("b", linear = true);
#else
view (theta = sin(t/10.), phi = cos(t/15));
scatter (loc, s = 30);
scalar phi[];
foreach_vertex()