Convergence of axisymmetric viscous terms

We wish to test the accuracy of the discretisation of viscous terms in axisymmetric coordinates. To do so, we “manufacture” a solution to the pure viscous diffusion problem t(yu)=x(yτxx)+y(yτxy)+sx t(yv)=x(yτyx)+y(yτyy)+sy2vy with τxx=2xu τxy=τyx=(yu+xv) τyy=2yv and where sx and sy are source terms to be determined. We look for solutions of the form u(x,y)=v(x,y)=yacosxsiny The viscous stress tensor is then (for μ=1) τxx=2yasinxsiny τxy=τyx=cosx(yacosy+aya1siny)yasinxsiny τyy=2cosx(yacosy+aya1siny) and the viscosity equation gives sx=ya[3ycosxsiny(2a+1)cosxcosya2y1cosxsiny+ysinxcosy+(a+1)sinxsiny] sy=ya[ysinxcosy+asinxsiny+3ycosxsiny2(2a+1)cosxcosy+2(1a2)y1cosxsiny] The game is then to solve the viscosity equation with the forcing terms defined above and recover the (stationary) solution for the velocity field.

We use the axisymmetric metric, the viscous solver and the standard time loop.

#include "grid/multigrid.h"
#include "axi.h"
#include "viscosity.h"
#include "run.h"

We impose the correct boundary conditions for u (the default boundary conditions for v are already correct).

vector u[];

u.t[top] = dirichlet(0);
u.t[bottom] = dirichlet(0);

int main()
  L0 = 2.*π;
  periodic (right);
  for (N = 16; N <= 128; N *= 2)

We solve only for the case a=1. The case a=0 is ill-posed because the forcing terms diverge near the axis of symmetry.

double a = 1.;
mgstats mg;

event init (i = 0) {
    u.x[] = u.y[] = 0.;
  boundary ((scalar *){u});

event integration (i++)
  double dt = 1.;

This is the time integration loop. We first add the forcing term…

  foreach() {
    u.x[] += dt*pow(y, a - 1.)*(3.*y*cos(x)*sin(y) - (2.*a + 1.)*cos(x)*cos(y)
				- sq(a)*cos(x)*sin(y)/y	
				+ y*sin(x)*cos(y) + (a + 1.)*sin(x)*sin(y));
    u.y[] += dt*pow(y, a - 1.)*(y*sin(x)*cos(y) + a*sin(x)*sin(y)
				+ 3.*y*cos(x)*sin(y)
				- 2.*(2.*a + 1.)*cos(x)*cos(y) +
				2.*(1. - sq(a))*cos(x)*sin(y)/y);
  boundary ((scalar *){u});

…and then solve for viscosity.

  mg = viscosity (u, fm, cm, dt);

event error (i = 20)

Twenty timesteps are enough to converge toward the stationary solution.

  scalar e[];
  foreach() {
    printf ("%g %g %g %g %g\n", x, y, u.x[], u.y[], pow(y,a)*cos(x)*sin(y));
    e[] = u.x[] - pow(y,a)*cos(x)*sin(y);
  norm n = normf (e);
  fprintf (stderr, "%d %g %g %g %d %d\n",
	   N, n.avg, n.rms, n.max, mg.i, mg.nrelax);

As expected we get second-order convergence.