# 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 ${\partial }_{x}{u}_{x}+{\partial }_{y}{u}_{y}+\frac{{u}_{y}}{y}=0$ ${\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 \left(\mu y\nabla {u}_{x}\right)$ ${\partial }_{t}{u}_{y}+{u}_{x}{\partial }_{x}{u}_{y}+{u}_{y}{\partial }_{y}{u}_{y}-\frac{{w}^{2}}{y}=-\frac{1}{\rho }{\partial }_{y}p+\frac{1}{\rho y}\left(\nabla \cdot \left(\mu y\nabla {u}_{y}\right)-2\mu \frac{{u}_{y}}{y}\right)$ ${\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 \left(\mu y\nabla w\right)-w\left(\frac{\mu }{y}+{\partial }_{y}\mu \right)\right]$ 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 ${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.;
boundary ((scalar *){a});
}
}``````

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 $\theta {\partial }_{t}w=\nabla \cdot \left(D\nabla w\right)+\beta w$ Identifying with the diffusion part of the equation for $w$ above, we have $\begin{array}{ccc}\hfill \theta & \hfill =\hfill & \rho y\hfill \\ \hfill D& \hfill =\hfill & \mu y\hfill \\ \hfill \beta & \hfill =\hfill & -\left(\rho {u}_{y}+\frac{\mu }{y}+{\partial }_{y}\mu \right)\hfill \end{array}$ 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 β[], θ[];
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, μ, θ = θ, β = β);
}``````