src/navier-stokes/swirl.h
Azimuthal velocity for axisymmetric flows
The centered Navier–Stokes solver can be combined with the axisymmetric metric but assumes zero azimuthal velocity (“swirl”). This file adds this azimuthal velocity component: the w field.
Assuming that x is the axial direction and y the radial direction, as in axi.h, the incompressible, variable-density and viscosity, axisymmetric Navier–Stokes equations (with swirl) can be written \displaystyle \partial_x u_x + \partial_y u_y + \frac{u_y}{y} = 0 \displaystyle \partial_t u_x + u_x \partial_x u_x + u_y \partial_y u_x = - \frac{1}{\rho} \partial_x p + \frac{1}{\rho y} \nabla \cdot (2 \mu y \nabla \mathbf{D}_x) \displaystyle \partial_t u_y + u_x \partial_x u_y + u_y \partial_y u_y - {\color{blue}\frac{w^2}{y}} = - \frac{1}{\rho} \partial_y p + \frac{1}{\rho y} \left( \nabla \cdot (2 \mu y \nabla \mathbf{D}_y) - 2 \mu \frac{u_y}{y} \right) \displaystyle {\color{blue} \partial_t w + u_x \partial_x w + u_y \partial_y w + \frac{u_y w}{y} = \frac{1}{\rho y} \left[ \nabla \cdot (\mu y \nabla w) - w \left( \frac{\mu}{y} + \partial_y \mu \right) \right] } where the terms in blue are the missing “swirl” terms. We will thus need to solve an advection-diffusion equation for w.
#include "tracer.h"
#include "diffusion.h"
scalar w[], * tracers = {w};The azimuthal velocity is zero on the axis of symmetry (y=0).
w[bottom] = dirichlet(0);We will need to add the acceleration term w^2/y in the evolution equation for u_y. If the acceleration field is not allocated yet, we do so.
event defaults (i = 0)
{
if (is_constant(a.x)) {
a = new face vector;
foreach_face()
a.x[] = 0.;
}
}The equation for u_y is solved by the centered Navier–Stokes solver combined with the axisymmetric metric, but the acceleration term w^2/y is missing. We add it here, taking care of the division by zero on the axis, and averaging w from cell center to cell face.
event acceleration (i++)
{
face vector av = a;
foreach_face (y)
av.y[] += y > 0. ? sq(w[] + w[0,-1])/(4.*y) : 0.;
}The advection of w is done by the tracer solver, but we need to add diffusion. Using the diffusion solver, we solve \displaystyle \theta\partial_tw = \nabla\cdot(D\nabla w) + \beta w Identifying with the diffusion part of the equation for w above, we have \displaystyle \begin{aligned} \theta &= \rho y \\ D &= \mu y \\ \beta &= - \left(\rho u_y + \frac{\mu}{y} + \partial_y\mu \right) \end{aligned} Note that the rho and mu fields (defined by the Navier–Stokes solver) already include the metric (i.e. are \rho y and \mu y), which explains the divisions by y in the code below.
event tracer_diffusion (i++)
{
scalar beta[], theta[];
foreach() {
theta[] = rho[];
double muc = (mu.x[] + mu.x[1] + mu.y[] + mu.y[0,1])/4.;
double dymu = (mu.y[0,1]/fm.y[0,1] - mu.y[]/fm.y[])/Delta;
beta[] = - (rho[]*u.y[] + muc/y)/y - dymu;
}
diffusion (w, dt, mu, theta = theta, beta = beta);
}