swirl.h

📄 View in API Reference (Doxygen) · View on basilisk.fr

Test cases (1): swirl

Examples (1): yang

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 (2 \mu y \nabla \mathbf{D}_x)

\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)

{\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" [api]
#include "diffusion.h" [api]

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

\theta\partial_tw = \nabla\cdot(D\nabla w) + \beta w

Identifying with the diffusion part of the equation for $w$ above, we have

\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);
}

Basilisk CFD

Navier–Stokes

Saint-Venant / Shallow Water

Multilayer

Compressible

Advection & VOF

Interfacial Forces

Two-Phase

Elliptic Solvers

Electrohydrodynamics

Viscoelasticity

Grid System

Core

Coordinates & Metrics

Input / Output

Utilities

Other

swirl.h

Azimuthal velocity for axisymmetric flows