src/axistream.h

Axisymmetric stream function

The axisymmetric or Stokes stream function ψ verifies ur=1rzψ uz=+1rrψ with ur and uz the radial and longitudinal velocity components of an incompressible axisymmetric flow.

Given the definition of the vorticity ω ω=zurruz ψ verifies the variable-coefficient Poisson equation z(1rzψ)+r(1rrψ)=ω This is not well-conditioned due to the divergence of the coefficients at r=0 and can be rewritten instead as 2ψz2+2ψr21rrψ=ωr which is better behaved.

This equation can be inverted using the multigrid solver combined with the relaxation and residual functions defined below.

#include "poisson.h"

static void relax_psi (scalar * al, scalar * bl, int l, void * data)
{
  scalar a = al[0], b = bl[0];
  foreach_level_or_leaf (l)
    a[] = (- sq(Δ)*b[] - Δ*(a[0,1] - a[0,-1])/(2.*y) +
	   a[1] + a[-1] + a[0,1] + a[0,-1])/4.;
}

static double residual_psi (scalar * al, scalar * bl, scalar * resl, void * data)
{
  scalar a = al[0], b = bl[0], res = resl[0];
  double maxres = 0.;
#if TREE
  /* conservative coarse/fine discretisation (2nd order) */
  face vector g[];
  foreach_face()
    g.x[] = face_gradient_x (a, 0);
  boundary_flux ({g});
  foreach (reduction(max:maxres)) {
    res[] = b[] + (a[0,1] - a[0,-1])/(2.*y*Δ);
    foreach_dimension()
      res[] -= (g.x[1] - g.x[])/Δ;
#else // !TREE
  /* "naive" discretisation (only 1st order on trees) */
  foreach (reduction(max:maxres)) {
    res[] = b[] + (a[0,1] - a[0,-1])/(2.*y*Δ);
    foreach_dimension()
      res[] += (face_gradient_x (a, 0) -
		face_gradient_x (a, 1))/Δ;  
#endif // !TREE
    if (fabs (res[]) > maxres)
      maxres = fabs (res[]);
  }
  boundary (resl);
  return maxres;
}

The function axistream() takes as input the u velocity field (with x=z and y=r) and returns the corresponding axisymmetric stream function.

mgstats axistream (vector u, scalar ψ)
{
  scalar ω[];
  foreach() {
    ω[] = y*(u.y[1] - u.y[-1] - u.x[0,1] + u.x[0,-1])/(2.*Δ);
    ψ[] = 0.;
  }
  ψ[bottom] = dirichlet(0);
  boundary ({ψ});
  return mg_solve ({ψ}, {ω}, residual_psi, relax_psi, NULL, 0, NULL, 1);
}

See also