src/henry.h

    Advection/diffusion of a soluble tracer

    We consider the transport and diffusion of a tracer c with different solubilities in the two-phases described by two-phase.h.

    The diffusion coefficients in the two phases are D_1 and D_2 and the jump in concentration at the interface is given by \displaystyle c_1 = \alpha c_2 The advection/diffusion equation for c can then be written \displaystyle \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) with f the volume fraction and \displaystyle D = \frac{D_1 D_2}{D_2 f + D_1 (1 - f)} the harmonic mean of the diffusion coefficients (see Haroun et al., 2010).

    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
event defaults (i = 0)
{
  for (scalar s in stracers) {
#if EMBED
    s.refine = s.prolongation = refine_embed_linear;
#else
    s.refine  = refine_linear;
#endif
    s.restriction = restriction_volume_average;
    s.dirty = true;
  }
}
#endif // TREE

    Advection

    To avoid numerical diffusion through the interface we use the VOF tracer transport scheme for the temporary fields \phi_1 and \phi_2, see section 3.2 of Farsoiya et al., 2021.

    static scalar * phi_tracers = NULL;

event vof (i++)
{
  phi_tracers = f.tracers;
  for (scalar c in stracers) {
    scalar phi1 = new scalar, phi2 = 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);

    \phi_1 and \phi_2 are computed from c as \displaystyle \phi_1 = c \frac{\alpha f}{\alpha f + (1 - f)} \displaystyle \phi_2 = c \frac{1 - f}{\alpha f + (1 - f)} 

        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 \displaystyle \frac{c^{n + 1} - c^n}{\Delta t} = \nabla\cdot (D \nabla c^{n + 1} + \beta c^{n + 1}) see section 3.2 of Farsoiya et al., 2021.

    Note that these functions are close to that in poisson.h and diffusion.h but with the additional term \nabla\cdot (\beta c^{n + 1}).

    struct HDiffusion {
  face vector D;
  face 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;
  face vector D = p->D, beta = p->beta;

  scalar c = a;
  foreach_level_or_leaf (l) {
    double n = - sq(Delta)*b[], d = cm[]/dt*sq(Delta);
    foreach_dimension() {  
      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;
  face vector D = p->D, beta = p->beta;
  double maxres = 0.;
#if TREE
  /* conservative coarse/fine discretisation (2nd order) */
  face vector g[];
  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[];
    foreach_dimension()
      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[];
    foreach_dimension()
      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;
}

event tracer_diffusion (i++)
{
  free (f.tracers);
  f.tracers = phi_tracers;
  for (scalar c in stracers) {

    The advected concentration is computed from \phi_1 and \phi_2 as \displaystyle c = \phi_1 + \phi_2 and these fields are then discarded.

        scalar phi1 = c.phi1, phi2 = c.phi2, r[];
    foreach() {
      c[] = phi1[] + phi2[];
      r[] = - cm[]*c[]/dt;
    }
    delete ({phi1, phi2});

    The diffusion equation for c is then solved using the multigrid solver and the residual and relaxation functions defined above.

        face vector D[], beta[];
    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

    [farsoiya2021]

    Palas Kumar Farsoiya, Stéphane Popinet, and Luc Deike. Bubble-mediated transfer of dilute gas in turbulence. Journal of Fluid Mechanics, 920(A34), June 2021. [ DOI | http | .pdf ]
    [haroun2010]

    Y Haroun, D Legendre, and L Raynal. Volume of fluid method for interfacial reactive mass transfer: application to stable liquid film. Chemical Engineering Science, 65(10):2896–2909, 2010. [ DOI ]

    Usage

    Tests