src/layered/dr.h
Boussinesq buoyancy
The momentum equation of the multilayer solver becomes \displaystyle \begin{aligned} \partial_t \left( h \mathbf{u} \right)_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \mathbf{u} \right)_k & = - gh_k \mathbf{{\nabla}} (\eta) {\color{blue} - \mathbf{{\nabla}} (h q)_k + \left[ q \mathbf{{\nabla}} z \right]_k} \end{aligned} where the terms in blue have been added and q is the (hydrostatic) pressure deviation due to (small) density variations.
The density variations are described by a scalar field T, which is advected by the flow, and an associated “equation of state” \Delta\rho(T) defined by the user. This gives the additional equations \displaystyle \begin{aligned} \partial_t \left( h T \right)_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} T \right)_k & = 0,\\ q(z) & = \int_0^z g \Delta \rho(T) dz \end{aligned}
scalar T;
event defaults (i = 0)
{
T = new scalar[nl];
tracers = list_append (tracers, T);
}
event cleanup (t = end, last)
{
delete ({T});
}The pressure gradient terms in blue are added to the acceleration of the multilayer solver.
event acceleration (i++)
{The hydrostatic pressure deviation q is stored on the interfaces between layers (consistently with the Keller box scheme discretisation, see Popinet, 2020). This gives the following vertical discrete integration scheme.
scalar q = new scalar[nl];
foreach() {
double ph = 0.;
for (point.l = nl - 1; point.l >= 0; point.l--) {
ph += G*drho(T[])*h[];
q[] = ph;
}
}Once the pressure deviation is known, the terms in blue above are added to the face acceleration field ha, using the pressure-gradient macro.
foreach_face()
hpg (pg, q, 0,
ha.x[] += pg);
delete ({q});
}