src/test/ehd_axi_stress.c

    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, \displaystyle \frac {\partial_y (y \epsilon \partial_y \phi)}{y} + \partial_x(\epsilon \partial_x \phi) = - \rho_e with \rho_e the electric charge density. If we look for a solution of the form, \displaystyle \phi (x, y) = y^a \cos x \sin y the charge density must then be (for \epsilon = 1), \displaystyle \rho_e = -y^{-2 + a} \cos x \left( (1 + 2 a) y \cos y + (a^2 - 2 y^2) \sin y \right)

    The corresponding electrical stresses \mathbf{F} are calculated from the divergence of the Maxwell stress tensor, \displaystyle \mathbf{F} = \nabla \cdot \mathbf{M} Their components in the x- and y- directions are, \displaystyle \mathbf{F}|_y = \frac{\partial_y (y M_{yy})}{y} + \partial_x M_{yx} - \frac{M_{\theta \theta}}{y} and \displaystyle \mathbf{F}|_x = \frac{\partial_y (y \partial_y M_{xy})}{y} + \partial_x M_{xx} with \displaystyle M_{x x} = \epsilon (E_x^2 - E_y^2)/2 \displaystyle M_{y y} = \epsilon (E_y^2 - E_x^2)/2 \displaystyle M_{x y} = M_{y x} =\epsilon (E_y E_x) and \displaystyle M_{\theta \theta} = -\epsilon (E_y^2 + E_x^2)/2

    E_{x/y} = -\partial_{x/y} \phi 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 \phi and the electrical stresses are recovered.

    #include "grid/multigrid.h"
    #include "axi.h"
    #include "run.h"
    #include "poisson.h"
    
    mgstats mg;
    scalar phi[];
    face vector epsilon[];
    face vector a[], alpha[];

    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))
    
    phi[top] = dirichlet(0.);
    phi[bottom] = dirichlet(0.);
    
    int main()
    {
      L0 = 2.*pi;
      periodic (right);
      for (N = 16; N <= 128; N *= 2)
        run();
    }
    
    event init (i = 0) {
      scalar rhoe[], rhs[];
      foreach() {
        rhoe[] = RHOE(x,y);
        rhs[] = -rhoe[]*cm[];
        phi[] = 0.;
      }
    
      foreach_face() {
        epsilon.x[] = fm.x[];
        alpha.x[] = fm.x[];
        a.x[] = 0.;
      }
      boundary ((scalar *){rhoe, phi, epsilon, alpha, a});
      
      mg = poisson (phi, rhs, epsilon);
    }
    
    #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.

    ftitle(a,b) = sprintf("%.0f/x^{%4.2f}", exp(a), -b)
    fx(x)=a+b*x
    fit fx(x) 'log' u (log($1)):(log($4)) via a,b
    fx2(x)=a2+b2*x
    fit fx2(x) 'log' u (log($1)):(log($2)) via a2,b2
    fy(x)=c+d*x
    fit fy(x) 'log' u (log($1)):(log($7)) via c,d
    fy2(x)=c2+d2*x
    fit fy2(x) 'log' u (log($1)):(log($5)) via c2,d2
    set xlabel 'Resolution'
    set ylabel 'Error'
    set logscale
    set xrange [8:256]
    set cbrange [1:2]
    set xtics 8,2,256
    set grid ytics
    set yrange [:]
    plot 'log' u 1:4 t 'Fx:max', exp(fx(log(x))) t ftitle(a,b), \
         'log' u 1:2 t 'Fx:norm1', exp(fx2(log(x))) t ftitle(a2,b2), \
         'log' u 1:7 t 'Fy:max', exp(fy(log(x))) t ftitle(c,d),     \
         'log' u 1:5 t 'Fy:norm1', exp(fy2(log(x))) t ftitle(c2,d2)
    Convergence (script)

    Convergence (script)