swirl.h

Go to the documentation of this file.

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