Convergence of axisymmetric EHD stresses

We wish to test the accuracy of the EHD stresses in case of axisymmetric coordinates. To do so, we create an “ad hoc” solution for the electric poisson equation, y(yεyϕ)y+x(εxϕ)=ρe with ρe the electric charge density. If we look for a solution of the form, ϕ(x,y)=yacosxsiny the charge density must then be (for ε=1), ρe=y2+acosx((1+2a)ycosy+(a22y2)siny)

The corresponding electrical stresses F are calculated from the divergence of the Maxwell stress tensor, F=M Their components in the x- and y- directions are, Fy=y(yMyy)y+xMyxMθθy and Fx=y(yyMxy)y+xMxx with Mxx=ε(Ex2Ey2)/2 Myy=ε(Ey2Ex2)/2 Mxy=Myx=ε(EyEx) and Mθθ=ε(Ey2+Ex2)/2

Ex/y=x/yϕ are the axial and radial components of the electric field.

In this test we will set the charge density to check how the electric potential ϕ and the electrical stresses are recovered.

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

mgstats mg;
scalar φ[];
face vector ε[];
face vector a[], α[];

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.

#define PHI(x,y) y*cos(x)*sin(y)
#define RHOE(x,y) cos(x)*(-3*cos(y) + (2*y-1/y)*sin(y))

#define FY(x,y) (y == 0. ? 0. : sq(cos(x))/y*(sq(sin(y)) + \
					      y/2*(y + 5*y*cos(2*y)	\
						   - 2*(-2+sq(y))*sin(2*y))))

#define FX(x,y) cos(x)*sin(x)*sin(y)*(-3*y*cos(y)+(2*sq(y)-1)*sin(y))

φ[top] = dirichlet(0.);
φ[bottom] = dirichlet(0.);

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

event init (i = 0) {
  scalar rhoe[], rhs[];
  foreach() {
    rhoe[] = RHOE(x,y);
    rhs[] = -rhoe[]*cm[];
    φ[] = 0.;

  foreach_face() {
    ε.x[] = fm.x[];
    α.x[] = fm.x[];
    a.x[] = 0.;
  boundary ((scalar *){rhoe, φ, ε, α, a});
  mg = poisson (φ, rhs, ε);

#include "ehd/stress.h"

event plot_err (i = 0)
  face vector err[];

  foreach_face (x) {
    err.x[] = a.x[] - FX(x,y);
    printf ("x: %g %g %g %g %g  \n", x, y, a.x[], FX(x,y), err.x[]);
  foreach_face (y) {
    err.y[] = a.y[] - FY(x,y);
    printf ("y: %g %g %g %g %g \n", x, y, a.y[], FY(x,y),err.y[]);
  norm nx = normf (err.x), ny = normf (err.y);
  fprintf (stderr, "%d %g %g %g %g %g %g\n", N,
	   nx.avg, nx.rms, nx.max,
	   ny.avg, ny.rms, ny.max);

As expected we get second-order convergence.