sandbox/farsoiya/phase_change/epstein-plesset.c
Epstein-Plesset test
We simulate a static bubble which shrinks due to diffusion of gas.
The time evolution of the interface is given by Epstein and Plesset, 1950 as, \displaystyle \frac{d R}{d t}=-\mathcal{D} M_w\frac{\alpha c_s - c_i}{\rho}\left\{\frac{1}{R}+\frac{1}{\sqrt{\pi \mathcal{D} t}}\right\} and concentration,
\displaystyle c = c_0 + (\alpha c_s - c_0) \frac{R(t)}{r} \operatorname{erfc}\left( \frac{r - R(t)}{\sqrt{4 \mathcal{D} t}} \right)
|
|
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 15})
plt.figure()
t,rep,cep = np.loadtxt('../ref-epstein-plesset',delimiter='\t',unpack=True);
plt.plot(t/4, rep,'k',label='Epstein-Plesset')
ts, rb, cb = np.loadtxt('conc-plesset-spherical',delimiter=' ',unpack=True)
plt.plot(ts/4,rb,'k--',label='Basilisk');
# plt.ylim(0.97,1.007)
plt.legend();
plt.xlabel(r'$t\; \mathscr{D}_l/d_0^2$')
plt.ylabel(r'$r(t)/r(0)$')
plt.tight_layout()
plt.savefig('p001cbt.png')
plt.figure()
plt.plot(t/4, cep,'k',label='Analytical')
plt.plot(ts/4,cb,'k--',label='Basilisk');
# plt.ylim(0.97,1.007)
plt.legend();
plt.xlabel(r'$t\; \mathscr{D}_l/d_0^2$')
plt.ylabel(r'$c_l/c_{b0}$')
plt.tight_layout()
plt.savefig('p001clt.png')
plt.figure()
plt.savefig(' ')
References
[farsoiya2021] |
P. K. Farsoiya, Q. Magdelaine, A. Antkowiak, S. Popinet, and L. Deike. Bubble mediated single component gas transfer in homogeneous isotropic turbulence. Journal of Fluid Mechanics, 2021. submitted. |
[epstein1950stability] |
Paul S Epstein and Milton S Plesset. On the stability of gas bubbles in liquid-gas solutions. The Journal of Chemical Physics, 18(11):1505–1509, 1950. |
#define R0 1.
#define L 10. // size of the box
#define MIN_LEVEL 5
#define LEVEL 10
#define MAX_LEVEL 10
#define dR_refine (2.*L0/(1 << LEVEL))
#define F_ERR 1e-10
#define T_END 600
#define DT_MAX 10.
#define DELTA_T 1. // for measurements and videos
#include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "tension.h"
#include "phase-change.h"
#include "reduced.h"
#define D_V 1.
#define cinf 0
scalar c[], *stracers = {c};
Thanks to symmetry, we only simulate quarter bubble
c[right] = dirichlet(cinf);
c[top] = dirichlet(cinf);
In the main function of the program, we set the domain geometry to be ten times larger than the bubble.
int main() {
size (L);
origin (0., 0.);
N = 1 << LEVEL;
init_grid (N);
f.sigma = 0;
rho1 = 1.;
rho2 = 1.;
mu1 = 1./10.;
mu2 = mu1/20.;
G.x = 0.;
Z.x = 0.;
c.inverse = false; // false: bubble, true: drop
c.alpha = 0.8; // solubility
c.D = D_V; // Diffusivity in outside fluid
c.mw = 0.001; //Molecular weight
c.tr_eq = 1.;
run();
}
The initial position of the interface is defined with this function:
#define circle(x, y, R) (sq(R) - sq(x) - sq(y))
Before the first step, we initialize the concentration field (after having refined the grid around the future interface): c_s in the bubble and c_\infty in the liquid.
event init (i = 0) {
#if TREE
refine (level < MAX_LEVEL && circle(x, y, (R0 - dR_refine)) < 0.
&& circle(x, y, (R0 + dR_refine)) > 0.);
#endif
fraction (f, -circle(x, y, R0));
//~ fraction (f, y - L0*0.5);
foreach() {
u.x[] = 0.;
c[] = f[]*cinf + (1. - f[])*c.tr_eq;
}
foreach_face()
uf.x[] = 0.;
boundary({ u, uf,c});
}
#if TREE
event adapt (i++) {
adapt_wavelet ({f, c,u}, (double[]){1e-3, 1e-3, 1e-1, 1e-1,1e-1,1e-1},
minlevel = MIN_LEVEL, maxlevel = MAX_LEVEL);
}
#endif
Post-processings and videos
We now juste have to write post-processing events to save tthe effective radius of the bubble.
event outputs (t = 0.; t += DELTA_T; t <= T_END) {
double effective_radius;
scalar fc[];
foreach()
fc[] = 1. - f[];
effective_radius = pow(3.*statsf(fc).sum, 1./3.);
fflush(stderr);
char name[80];
sprintf (name, "conc-plesset-spherical");
static FILE * fp = fopen (name, "w");
fprintf (fp, "%.17g %.17g %.17g\n", t, effective_radius, interpolate(c, (1. + 0.2)*cos(M_PI/4), (1. + 0.2)*sin(M_PI/4) , 0) );
fflush (fp);
output_ppm (f, file = "f.mp4", box = {{0,0},{3,3}},
linear = true, min = 0, max = 1);
output_ppm (c, file = "c.mp4", box = {{0,0},{3,3}},
linear = true, min = 0, max = c.tr_eq);
}