src/test/cyl_axi.c

Charge relaxation in an axisymmetric insulated conducting column

A conducting rigid cylinder of radius Ro=0.1 is immersed in an insulating medium. Initially an uniform charge volume density is set in the cylinder (ρe(x,0)=0.5). As time passes the charge migrates from the bulk to the interface of the cylinder but the total charge in the cylinder section is preserved. Q(t)=Σρe(x,t)dΣ=Qo=πRo2ρe(x,0). The outer electric field reaches a steady-state.

A more detailed discussion of this simulation is given in Lopez-Herrera et al, 2011.

The charges, initially uniformly distributed, tend to accumulate at the interface due to electrostatic repulsion (charge relaxation).

Since the total charge at the column remains constant, the electric potential distribution at the outer medium will also remain constant.

If the macro NAVIER is set to 1 a true EHD problem is solved otherwise the problem is purely electrostatic. The ohmic conduction is calculated implicitly although the conservative explicit scheme could also be used.

#define NAVIER 1
int LEVEL;

#include "grid/multigrid.h"
#include "axi.h"
#if NAVIER
# include "navier-stokes/centered.h"
# include "ehd/implicit.h"
# include "ehd/stress.h"
#else
# include "ehd/implicit.h"
#endif
#include "fractions.h"

Far away from the conducting column the electric potential is zero.

φ[top]   = dirichlet(0);
#if NAVIER
p[top]     = dirichlet(0);
u.n[top]   = neumann(0);
#endif

#define β 3.
#define cond  3.
#define rhoini 0.5
#define R 0.1
#define A (R*R*rhoini/2.)
#define cylinder(y) (R - y)

scalar f[];

event init (t = 0) {

The conducting media will be defined through the scalar f. This will be useful also to define the electrical properties of the media and the initial distribution of charge density.

  face vector s[];
  vertex scalar ψ[];
  foreach_vertex ()
    ψ[] = cylinder(y);
  fractions (ψ, f, s);
  foreach()
    rhoe[] = rhoini*f[];
  boundary ({rhoe});

The electrical conductivity and permittivity are defined on faces. The face permittivity is constructed from the center values of the VOF scalar f by averaging. On the other hand, face conductivity is set using the reconstructed surface fractions s.

Note that the electrostatic properties ɛ and K, i.e electrical permittivity and conductivity, are multiplied by the face metric factors.

  foreach_face(){
    double ff = (f[] + f[-1])/2.;
    ε.x[] = (ff*β + (1. - ff))*fm.x[];
    K.x[] = cond*s.x[]*fm.x[];
  }
  boundary ((scalar *){ε, K});
}

This event checks the conservation of the total charge.

event chargesum (i++) {
  double Q = statsf(rhoe).sum;
  static double Q0;
  if (i == 0)
    Q0 = Q;
  else
    assert (fabs(Q - Q0) < 1e-7);
}

At the final instant, when the charge is fully relaxed on the interface, the radial electric field and the pressure distribution are output in order to compare with analytical results.

event epfield (t = 10) {

The name of the file with the results is formed with the concatenation of Er with the corresponding level.

  char name[80];
  sprintf (name, "Er-%d", LEVEL);
  FILE * fp = fopen (name, "w");
  scalar ee[];
  foreach() {
    double Ey = (φ[0,-1] - φ[0,1])/(2*Δ);
    ee[] = Ey - (y < R ? 0 : 0.5*sq(R)*rhoini/y);
#if NAVIER
    fprintf (fp, "%g %g %g %g\n", y, Ey, p[], rhoe[]);
#else
    fprintf (fp, "%g %g %g\n", y, Ey, rhoe[]);
#endif
  }
  fclose (fp);

We also log the error on the norm of the electric field.

  norm n = normf (ee);
  fprintf (stderr, "%d %g %g %g\n", LEVEL, n.avg, n.rms, n.max);
}

int main() {
  X0 = -0.5;
  L0 = 1.;
  DT = 1;
  TOLERANCE = 1e-7;
  for (LEVEL = 6; LEVEL <= 8; LEVEL++) {
    N = 1 << LEVEL;
    run();
  }
}
Radial electric field distribution as a function of grid refinement.

Radial electric field distribution as a function of grid refinement.

Pressure distribution for different grid refinements.

Pressure distribution for different grid refinements.

Convergence of error on the radial component of the electric field

Convergence of error on the radial component of the electric field

See also