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 xux+yuy+uyy=0 tux+uxxux+uyyux=1ρxp+1ρy(μyux) tuy+uxxuy+uyyuyw2y=1ρyp+1ρy((μyuy)2μuyy) tw+uxxw+uyyw+uywy=1ρy[(μyw)w(μy+yμ)] We will thus need to solve and 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 w2/y in the evolution equation for uy. 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.;
    boundary ((scalar *){a});
  }
}

The equation for uy is solved by the centered Navier–Stokes solver combined with the axisymmetric metric, but the acceleration term w2/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 θtw=(Dw)+βw Identifying with the diffusion part of the equation for w above, we have θ=ρyD=μyβ=(ρuy+μy+yμ) Note that the rho and mu fields (defined by the Navier–Stokes solver) already include the metric (i.e. are ρy and μy), which explains the divisions by y in the code below.

event tracer_diffusion (i++)
{
  scalar β[], θ[];
  foreach() {
    θ[] = ρ[];
    double muc = (μ.x[] + μ.x[1] + μ.y[] + μ.y[0,1])/4.;
    double dymu = (μ.y[0,1]/fm.y[0,1] - μ.y[]/fm.y[])/Δ;
    β[] = - (ρ[]*u.y[] + muc/y)/y - dymu;
  }
  diffusion (w, dt, μ, θ = θ, β = β);
}

Usage

Examples

Tests