stress.h

Go to the documentation of this file.

/** @file stress.h
 */
/**
# Electrohydrodynamic stresses
 
The EHD force density, \f$\mathbf{f}_e\f$, can be computed as the
divergence of the Maxwell stress tensor \f$\mathbf{M}\f$,
\f[
M_{ij} = \varepsilon (E_i E_j - \frac{E^2}{2}\delta_{ij})
\f]
where \f$E_i\f$ is the \f$i\f$-component of the electric field,
\f$\mathbf{E}=-\nabla \phi\f$ and \f$\delta_{ij}\f$ is the Kronecker delta.
 
We need to add the corresponding acceleration to the 
[Navier--Stokes solver](/src/navier-stokes/centered.h).
 
If the acceleration vector *a* (defined by the Navier--Stokes solver)
is constant, we make it variable. */
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) {
  if (is_constant (a.x))
    a = {0} /* new vector */;
}
 
/**
We overload the 
[acceleration event](/src/navier-stokes/centered.h#acceleration-term) 
of the Navier--Stokes solver to add the electrohydrodynamics acceleration
term. */
 
/** @brief Event: acceleration (i++) */
void event_acceleration (void) {
  assert (dimension <= 2); // not 3D yet
  vector f[];
  for (int _d = 0; _d < dimension; _d++) {
    vector Mx[];
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      Mx.x[] = epsilon.x[]/2.*(sq(phi[] - phi[-1,0]) - 
                               sq(phi[0,1] - phi[0,-1] + 
                                  phi[-1,1] - phi[-1,-1])/16.)/sq(Delta);
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      Mx.y[] = epsilon.y[]*(phi[] - phi[0,-1])*
      (phi[1,0] - phi[-1,0] + 
       phi[1,-1] - phi[-1,-1])/sq(2.*Delta);
 
    /**
    The electric force is the divergence of the Maxwell stress tensor
    \f$\mathbf{M}\f$. */
 
    for (int _i = 0; _i < _N; _i++) /* foreach */
      f.x[] = (Mx.x[1,0] - Mx.x[] + Mx.y[0,1] - Mx.y[])/(Delta*cm[]);
  }
 
  /**
  If [axisymmetric cylindrical coordinates](/src/axi.h) are used an 
  additional term must be added. */
 
#if AXI
  for (int _i = 0; _i < _N; _i++) /* foreach */
    f.y[] += (sq(phi[1,0] - phi[-1,0]) +
        sq(phi[0,1] - phi[0,-1]))/(8.*cm[]*sq(Delta))
    *(epsilon.x[]/fm.x[] + epsilon.y[]/fm.y[] +
      epsilon.x[1,0]/fm.x[1,0] + epsilon.y[0,1]/fm.y[0,1])/4.;
#endif
 
  /**
  To get the acceleration from the force we need to multiply by the
  specific volume \f$\alpha\f$. */
 
  vector av = a;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    av.x[] += alpha.x[]/fm.x[]*(f.x[] + f.x[-1])/2.;
}