henry.h

Go to the documentation of this file.

/** @file henry.h
 */
/**
# Advection/diffusion of a soluble tracer
 
We consider the transport and diffusion of a tracer \f$c\f$ with different
solubilities in the two-phases described by
[two-phase.h](/src/two-phase.h).
 
The diffusion coefficients in the two phases are \f$D_1\f$ and \f$D_2\f$ and
the jump in concentration at the interface is given by
\f[
c_1 = \alpha c_2
\f]
The advection/diffusion equation for \f$c\f$ can then be written
\f[
\partial_t c + \nabla\cdot(\mathbf{u} c) = 
   \nabla\cdot\left(D\nabla c - D \frac{c(\alpha - 1)}
                     {\alpha f + (1 - f)}\nabla f\right)
\f]
with \f$f\f$ the volume fraction and
\f[
D = \frac{D_1 D_2}{D_2 f + D_1 (1 - f)}
\f]
the harmonic mean of the diffusion coefficients (see [Haroun et al.,
2010](#haroun2010)).
 
The diffusion coefficients and solubility are attributes of each
tracer.
 
The *stracers* list of soluble tracers must be defined by the calling
code. */
 
attribute {
  double D1, D2, alpha;  //D1 for f = 1, D2 for f = 0
  scalar phi1, phi2; // private
}
 
extern scalar * stracers;
 
/**
## Defaults
 
On trees we need to ensure conservation of the tracer when
refining/coarsening. */
 
#if TREE
/** @brief Event: defaults (i = 0) */
void event_defaults (void) 
{
  for (int _s = 0; _s < 1; _s++) /* scalar in stracers */ {
#if EMBED
    s.refine = refine_embed_linear;
    set_prolongation (s, refine_embed_linear);
#else
    s.refine  = refine_linear;
#endif
    set_restriction (s, restriction_volume_average);
  }
}
#endif // TREE
 
/**
## Advection
 
To avoid numerical diffusion through the interface we use the [VOF
tracer transport scheme](/src/vof.h) for the temporary fields
\f$\phi_1\f$ and \f$\phi_2\f$, see section 3.2 of [Farsoiya et al.,
2021](#farsoiya2021). */
 
static scalar * phi_tracers = NULL;
 
/** @brief Event: vof (i++) */
void event_vof (void) 
{
  phi_tracers = f.tracers;
  for (int _c = 0; _c < 1; _c++) /* scalar in stracers */ {
    scalar phi1 = {0} /* new scalar */, phi2 = {0} /* new scalar */;
    c.phi1 = phi1, c.phi2 = phi2;
    scalar_clone (phi1, c);
    scalar_clone (phi2, c);
    phi2.inverse = true;
    
    f.tracers = list_append (f.tracers, phi1);
    f.tracers = list_append (f.tracers, phi2);
 
    /**
    \f$\phi_1\f$ and \f$\phi_2\f$ are computed from \f$c\f$ as
    \f[
    \phi_1 = c \frac{\alpha f}{\alpha f + (1 - f)}
    \f]
    \f[
    \phi_2 = c \frac{1 - f}{\alpha f + (1 - f)}
    \f]
    */
      
    for (int _i = 0; _i < _N; _i++) /* foreach */ {
      double a = c[]/(f[]*c.alpha + (1. - f[]));
      phi1[] = a*f[]*c.alpha;
      phi2[] = a*(1. - f[]);
    }
  }
}
 
/**
## Diffusion
 
We first define the relaxation and residual functions needed to solve
the implicit discrete system
\f[
\frac{c^{n + 1} - c^n}{\Delta t} = 
\nabla\cdot (D \nabla c^{n + 1} + \beta c^{n + 1})
\f]
see section 3.2 of [Farsoiya et al., 2021](#farsoiya2021).
 
Note that these functions are close to that in [poisson.h]() and
[diffusion.h]() but with the additional term \f$\nabla\cdot (\beta c^{n
+ 1})\f$. */
 
struct HDiffusion {
  vector D;
  vector beta;
};
 
static void h_relax (scalar * al, scalar * bl, int l, void * data)
{
  scalar a = al[0], b = bl[0];
  struct HDiffusion * p = (struct HDiffusion *) data;
  vector D = p->D, beta = p->beta;
 
  scalar c = a;
  for (int _i = 0; _i < l; _i++) {
    double n = - sq(Delta)*b[], d = cm[]/dt*sq(Delta);
    for (int _d = 0; _d < dimension; _d++) {  
      n += D.x[1]*a[1] + D.x[]*a[-1] +
  Delta*(beta.x[1]*a[1] - beta.x[]*a[-1])/2.;
      d += D.x[1] + D.x[] - Delta*(beta.x[1] - beta.x[])/2.;
    }
    c[] = n/d;
  }
}
 
static double h_residual (scalar * al, scalar * bl, scalar * resl, void * data)
{
  scalar a = al[0], b = bl[0], res = resl[0];
  struct HDiffusion * p = (struct HDiffusion *) data;
  vector D = p->D, beta = p->beta;
  double maxres = 0.;
#if TREE
  /* conservative coarse/fine discretisation (2nd order) */
  vector g[];
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    g.x[] = D.x[]*face_gradient_x (a, 0) + beta.x[]*face_value (a, 0);
  foreach (reduction(max:maxres)) {
    res[] = b[] + cm[]/dt*a[];
    for (int _d = 0; _d < dimension; _d++)
      res[] -= (g.x[1] - g.x[])/Delta;
    if (fabs (res[]) > maxres)
      maxres = fabs (res[]);
  }
#else // !TREE
  /* "naive" discretisation (only 1st order on trees) */
  foreach (reduction(max:maxres)) {
    res[] = b[] + cm[]/dt*a[];
    for (int _d = 0; _d < dimension; _d++)
      res[] -= (D.x[1]*face_gradient_x (a, 1) -
    D.x[0]*face_gradient_x (a, 0) +
    beta.x[1]*face_value (a, 1) -
    beta.x[0]*face_value (a, 0))/Delta;     
    if (fabs (res[]) > maxres)
      maxres = fabs (res[]);
  }
#endif // !TREE    
  return maxres;
}
 
/** @brief Event: tracer_diffusion (i++) */
void event_tracer_diffusion (void) 
{
  free (f.tracers);
  f.tracers = phi_tracers;
  for (int _c = 0; _c < 1; _c++) /* scalar in stracers */ {
 
    /**
    The advected concentration is computed from \f$\phi_1\f$ and \f$\phi_2\f$ as
    \f[
    c = \phi_1 + \phi_2
    \f]
    and these fields are then discarded. */
    
    scalar phi1 = c.phi1, phi2 = c.phi2, r[];
    for (int _i = 0; _i < _N; _i++) /* foreach */ {
      c[] = phi1[] + phi2[];
      r[] = - cm[]*c[]/dt;
    }
    delete ({phi1, phi2});
 
    /**
    The diffusion equation for \f$c\f$ is then solved using the multigrid
    solver and the residual and relaxation functions defined above. */
 
    vector D[], beta[];
    for (int _i = 0; _i < _N; _i++) /* foreach_face */ {
      double ff = (f[] + f[-1])/2.;
      D.x[] = fm.x[]*c.D1*c.D2/(c.D1*(1. - ff) + ff*c.D2);
      beta.x[] = - D.x[]*(c.alpha - 1.)/
  (ff*c.alpha + (1. - ff))*(f[] - f[-1])/Delta;
    }
  
    restriction ({D, beta, cm});
    struct HDiffusion q;
    q.D = D;
    q.beta = beta;
    mg_solve ({c}, {r}, h_residual, h_relax, &q);
  }
}
 
/**
## References
 
~~~bib
@article{haroun2010,
  title = {Volume of fluid method for interfacial reactive mass transfer: 
           application to stable liquid film},
  author = {Haroun, Y and Legendre, D and Raynal, L},
  journal = {Chemical Engineering Science},
  volume = {65},
  number = {10},
  pages = {2896--2909},
  year = {2010},
  doi = {10.1016/j.ces.2010.01.012}
}
 
@hal{farsoiya2021, hal-03227997}
~~~
*/