# Electrohydrodynamic stresses

The EHD force density, ${\mathbf{\text{f}}}_{e}$, can be computed as the divergence of the Maxwell stress tensor $\mathbf{\text{M}}$, ${M}_{ij}=\varepsilon \left({E}_{i}{E}_{j}-\frac{{E}^{2}}{2}{\delta }_{ij}\right)$ where ${E}_{i}$ is the $i$-component of the electric field, $\mathbf{\text{E}}=-\nabla \varphi$ and ${\delta }_{ij}$ is the Kronecker delta.

We need to add the corresponding acceleration to the Navier–Stokes solver.

If the acceleration vector a (defined by the Navier–Stokes solver) is constant, we make it variable.

``````event defaults (i = 0) {
if (is_constant (a.x))
a = new face vector;
}``````

We overload the acceleration event of the Navier–Stokes solver to add the electrohydrodynamics acceleration term.

``````event acceleration (i++) {
assert (dimension <= 2); // not 3D yet
vector f[];
foreach_dimension() {
face vector Mx[];
foreach_face(x)
Mx.x[] = ε.x[]/2.*(sq(φ[] - φ[-1,0]) -
sq(φ[0,1] - φ[0,-1] +
φ[-1,1] - φ[-1,-1])/16.)/sq(Δ);
foreach_face(y)
Mx.y[] = ε.y[]*(φ[] - φ[0,-1])*
(φ[1,0] - φ[-1,0] +
φ[1,-1] - φ[-1,-1])/sq(2.*Δ);
boundary_flux ({Mx});``````

The electric force is the divergence of the Maxwell stress tensor $\mathbf{\text{M}}$.

``````    foreach()
f.x[] = (Mx.x[1,0] - Mx.x[] + Mx.y[0,1] - Mx.y[])/(Δ*cm[]);
}``````

If axisymmetric cylindrical coordinates are used an additional term must be added.

``````#if AXI
foreach()
f.y[] += (sq(φ[1,0] - φ[-1,0]) +
sq(φ[0,1] - φ[0,-1]))/(8.*cm[]*sq(Δ))
*(ε.x[]/fm.x[] + ε.y[]/fm.y[] +
ε.x[1,0]/fm.x[1,0] + ε.y[0,1]/fm.y[0,1])/4.;
#endif
boundary ((scalar *){f});``````

To get the acceleration from the force we need to multiply by the specific volume $\alpha$.

``````  face vector av = a;
foreach_face()
av.x[] += α.x[]/fm.x[]*(f.x[] + f.x[-1])/2.;
}``````