src/ehd/implicit.h

Ohmic conduction

This function approximates implicitly the EHD equation set given by the electric Poisson equation and the ohmic conduction term of the charge density conservation (the advection term is computed elsewhere using a tracer advection scheme), (εϕn+1)=ρen+1and(ρen+1ρen)=Δt(Kϕn+1) where ρe is the charge density, and ϕ is the electric potential. K and ε are the conductivity and permittivity respectively.

Substituting the Poisson equation into the conservation equation gives the following time-implicit scheme, [(KΔt+ε)ϕn+1]=ρen

We need the Poisson solver and the timestep dt.

#include "poisson.h"
extern double dt;

The electric potential and the volume charge density are scalars while the permittivity and conductivity are face vectors. In mgphi we will store the statistics for the multigrid resolution of the electric poisson equation.

scalar φ[], rhoe[];
face vector ε[], K[];
mgstats mgphi;

event defaults (i = 0)
{

The restriction/refine attributes of the charge density are those of a tracer otherwise the conservation is not guaranteed.

#if TREE
  rhoe.restriction = restriction_volume_average;
  rhoe.refine = refine_linear;
#endif

By default the permittivity is unity and other quantities are zero.

  foreach_face()
    ε.x[] = K.x[] = fm.x[];
  boundary ((scalar *){ε, K});
}

event tracer_diffusion (i++)
{
  scalar rhs[];

The r.h.s of the poisson equation is set to ρen,

  foreach()
    rhs[] = - rhoe[]*cm[];

  if (K.x.i) {

We compute the coefficients of the Laplacian: (Kdt+ε).

    face vector f[];
    foreach_face()
      f.x[] = K.x[]*dt + ε.x[];
    boundary ((scalar *){f});

The poisson equation is solved to obtain ϕn+1.

    mgphi = poisson (φ, rhs, f);

Finally ρen+1 is computed as ρen+1=(εϕn+1).

#if TREE
    foreach_face()
      f.x[] = ε.x[]*(φ[] - φ[-1])/Δ;
    boundary_flux ({f});
    foreach() {
      rhoe[] = 0.;
      foreach_dimension()
	rhoe[] -= f.x[1] - f.x[];
      rhoe[] /= cm[]*Δ;
    }
#else // Cartesian
    foreach() {
      rhoe[] = 0.;
      foreach_dimension()
	rhoe[] -= (ε.x[1]*(φ[1] - φ[]) -
		   ε.x[]*(φ[] - φ[-1]));
      rhoe[] /= cm[]*sq(Δ);
    }
#endif
    boundary ({rhoe});
  }

In the absence of conductivity, the electric potential only varies if the electrical permittivity varies,

  else
    mgphi = poisson (φ, rhs, ε);
}

Usage

Tests