stress.h

Electrohydrodynamic stresses

Electrohydrodynamic stresses

The EHD force density, \mathbf{f}_e, can be computed as the divergence of the Maxwell stress tensor \mathbf{M}, \displaystyle M_{ij} = \varepsilon (E_i E_j - \frac{E^2}{2}\delta_{ij}) where E_i is the i-component of the electric field, \mathbf{E}=-\nabla \phi 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;
}

The logical place of these macros is in common.h as others. However, meanwhile…

#ifndef EMBED
#define face_avg_gradient_t1_x(a,i) \
  ((a[1,i-1] + a[1,i] - a[-1,i-1] - a[-1,i])/(4.*Delta)) 
#define face_avg_gradient_t2_x(a,i) \
  ((a[1,0,i-1] + a[1,0,i] - a[-1,0,i-1] - a[-1,0,i])/(4.*Delta)) 

#define face_avg_gradient_t1_y(a,i) \
  ((a[i-1,1] + a[i,1] - a[i-1,-1] - a[i,-1])/(4.*Delta)) 
#define face_avg_gradient_t2_y(a,i) \
  ((a[0,1,i-1] + a[0,1,i] - a[0,-1,i-1] - a[0,-1,i])/(4.*Delta))

#define face_avg_gradient_t1_z(a,i) \
  ((a[i-1,0,1] + a[i,0,1] - a[i-1,0,-1] - a[i,0,-1])/(4.*Delta)) 
#define face_avg_gradient_t2_z(a,i) \
  ((a[0,i-1,1] + a[0,i,1] - a[i-1,0,-1] - a[0,i,-1])/(4.*Delta)) 
#endif

If COEF is true the bilinear interpolation is used not only on the embed segment but for the faces of the mixed cells.

#if EMBED

#define COEF 1

#if COEF

static void electric_force_mixed (Point point, scalar cs, scalar phi,
				  face vector perm, coord * force)
{
  assert(cs[] > 0. && cs[] < 1.); //Only mixed cells

  double area = 0.;
  coord o, n, of, E, np;
  pseudo_t A1;
  pseudo_v A2;
  

  area = coef_bilinear_corrected_embed_gradient (point, phi, &A1, &A2, &o, &n);
  foreach_dimension()
    np.x = sign(-n.x);

We compute first the momentum flux through the embed segment.

  double En = 0., E2 = 0., fa = 0., perma = 0.;
  foreach_dimension() {
    E.x = A1.x.x/Delta;
    En += E.x*n.x;
    E2 +=  sq(E.x);
    perma += perm.x[] + perm.x[1];
    fa += fm.x[] + fm.x[1];
  }
  perma /= fa;
  
  foreach_dimension()
    force->x = perma*(En*E.x-E2/2.*n.x)*area;

  for(int ii = 0; ii < 2; ii++)
    foreach_dimension () {
    
      of.x = -o.x + (2*ii-1)/2.;
      of.y = -o.y + np.y*(1.-fs.x[ii])/2.;
    
      foreach_dimension()
	E.x = (A1.x.x + A1.x.y*of.y
#if dimension > 2
	       + A1.x.z*of.z
#endif
	       + 2.*A2.x*of.x)/Delta;

      force->x += (2.*ii-1)*perm.x[ii]/2.*(sq(E.x)
#if dimension > 2
					   -sq(E.z)
#endif
					   -sq(E.y));
      force->y += (2.*ii-1)*perm.x[ii]*E.x*E.y;
#if dimension > 2
      force->z += (2.*ii-1)*perm.x[ii]*E.x*E.z;
#endif
    }
}
#endif
#endif

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

event acceleration (i++) {
  vector f[];
  foreach_dimension() {
    face vector Mx[];
    foreach_face(x)
      Mx.x[] = epsilon.x[]/2.*(sq(face_gradient_x (phi,0))  
                               - sq(face_avg_gradient_t1_y(phi, 0))
#if dimension > 2			       
                               - sq(face_avg_gradient_t1_z(phi, 0))
#endif
			       );
    foreach_face(y)
      Mx.y[] = epsilon.y[]*face_gradient_y (phi, 0)
      *face_avg_gradient_t1_x (phi, 0);

#if dimension > 2
    foreach_face(z)
      Mx.z[] = epsilon.z[]*(face_gradient_z (phi, 0))*
                          *face_avg_gradient_t1_x (phi, 0);
#endif

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

    foreach() {
      double d = 0.;
      foreach_dimension()
	d += Mx.x[1,0] - Mx.x[];
      f.x[] = d/Delta;
    }
  }

If case of embed solid the electrical stresses on the solid must be taken into account.

#if EMBED || AXI
  foreach() {
#if EMBED
    if(cs[] > 0 && cs[] < 1.) {
#if COEF      
      coord F;
      electric_force_mixed (point, cs, phi, epsilon, &F);
      foreach_dimension() 
      	f.x[] = F.x/Delta;
#else
      double area = 0.;
      coord grad, n;
      area = bilinear_corrected_embed_gradient (point, phi, &grad, &n);
      double E2 = 0., En = 0., epsilona = 0., fa = 0.;
      foreach_dimension() {
	En += grad.x*n.x;
	E2 += sq(grad.x);
	epsilona += epsilon.x[] + epsilon.x[1];
	fa  += fm.x[] + fm.x[1];
      }
      epsilona /= fa;
      foreach_dimension()
      	f.x[] += epsilona*(En*grad.x-E2/2.*n.x)*area/Delta;
#endif
    }
#endif
  
#if AXI
    double E2 = 0., epsilona = 0., fa = 0.;
    foreach_dimension() {
      E2 += sq(center_gradient(phi));
      epsilona += epsilon.x[] + epsilon.x[1];
      fa += fm.x[] + fm.x[1];
    }
    epsilona /= fa;
    f.y[] += epsilona*E2/2.;
#endif
  }
#endif

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

  face vector av = a;
#if EMBED  
  foreach_face()
    av.x[] += (cs[] > 0. && cs[-1] > 0.) ?
    (f.x[]/rho[] + f.x[-1]/rho[-1])/2. :
    (f.x[] + f.x[-1])/(rho[] + rho[-1] + SEPS);
#else  
  foreach_face()
    av.x[] += (f.x[]/rho[] + f.x[-1]/rho[-1])/2.;
#endif
}