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

c_1 = \alpha c_2

The advection/diffusion equation for $c$ can then be written

\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

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 = 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 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

\phi_1 = c \frac{\alpha f}{\alpha f + (1 - f)}

\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

\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

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

~~~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} ~~~

Basilisk CFD

Navier–Stokes

Saint-Venant / Shallow Water

Multilayer

Compressible

Advection & VOF

Interfacial Forces

Two-Phase

Elliptic Solvers

Electrohydrodynamics

Viscoelasticity

Grid System

Core

Coordinates & Metrics

Input / Output