# 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), $$\nabla \cdot (\epsilon \nabla {\varphi}^{n+1})=-{\rho}_{e}^{n+1}\phantom{\rule{1em}{0ex}}\mathrm{\text{and}}\phantom{\rule{1em}{0ex}}({\rho}_{e}^{n+1}-{\rho}_{e}^{n})=\Delta t\nabla \cdot (K\nabla {\varphi}^{n+1})$$ where ${\rho}_{e}$ is the charge density, and $\varphi $ 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, $$\nabla \cdot [(K\phantom{\rule{0.167em}{0ex}}\Delta t+\epsilon )\nabla {\varphi}^{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 φ[], 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 $-{\rho}_{e}^{n}$,

```
foreach()
rhs[] = - rhoe[]*cm[];
if (K.x.i) {
```

We compute the coefficients of the Laplacian: $(Kdt+\epsilon )$.

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

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

` mgphi = poisson (φ, rhs, f);`

Finally ${\rho}_{e}^{n+1}$ is computed as ${\rho}_{e}^{n+1}=-\nabla \cdot (\epsilon \nabla {\varphi}^{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, ε);
}
```