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), \displaystyle \nabla \cdot( \epsilon \nabla \phi^{n+1}) = -\rho_e^{n+1} \quad \mathrm{and} \quad (\rho_e^{n+1}-\rho_e^n) = \Delta t\nabla \cdot (K\nabla \phi^{n+1}) where \rho_e is the charge density, and \phi is the electric potential. K and \epsilon are the conductivity and permittivity respectively.

    Substituting the Poisson equation into the conservation equation gives the following time-implicit scheme, \displaystyle \nabla \cdot [(K \, \Delta t + \epsilon) \nabla \phi^{n+1}] = -\rho_e^{n}

    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 phi[], rhoe[];
face vector epsilon[], 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()
    epsilon.x[] = K.x[] = fm.x[];
}

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

    The r.h.s of the poisson equation is set to -\rho_e^n,

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

  if (K.x.i) {

    We compute the coefficients of the Laplacian: (K dt +\epsilon).

        face vector f[];
    foreach_face()
      f.x[] = K.x[]*dt + epsilon.x[];

    The poisson equation is solved to obtain \phi^{n+1}.

        mgphi = poisson (phi, rhs, f);

    Finally \rho_e^{n+1} is computed as \rho_e^{n+1} = -\nabla \cdot (\epsilon \nabla \phi^{n+1}).

    #if TREE
    foreach_face()
      f.x[] = epsilon.x[]*(phi[] - phi[-1])/Delta;
    foreach() {
      rhoe[] = 0.;
      foreach_dimension()
	rhoe[] -= f.x[1] - f.x[];
      rhoe[] /= cm[]*Delta;
    }
#else // Cartesian
    foreach() {
      rhoe[] = 0.;
      foreach_dimension()
	rhoe[] -= (epsilon.x[1]*(phi[1] - phi[]) -
		   epsilon.x[]*(phi[] - phi[-1]));
      rhoe[] /= cm[]*sq(Delta);
    }
#endif
  }

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

      else
    mgphi = poisson (phi, rhs, epsilon);
}

    Usage

    Tests