Electrohydrodynamic stresses

The EHD force density, fe, can be computed as the divergence of the Maxwell stress tensor M, Mij=ɛ(EiEjE22δij) where Ei is the i-component of the electric field, E=ϕ and δ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[];
      Mx.x[] = ε.x[]/2.*(sq(φ[] - φ[-1,0]) - 
                               sq(φ[0,1] - φ[0,-1] + 
                                  φ[-1,1] - φ[-1,-1])/16.)/sq(Δ);
      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 M.

      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
    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.;
  boundary ((scalar *){f});

To get the acceleration from the force we need to multiply by the specific volume α.

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