src/test/axi_rising_bubble.c
Soluble gas diffusing from a rising bubble
This is the example discussed in section 3.3.2 of Farsoiya et al., 2021.
The four cases illustrated below are considered. The left-half of each figure is the vertical velocity component and the right-half the tracer concentration when close to the stationary regime.
|
|
|
|
The Sherwood number characterises the gas diffusion from the bubble to the liquid. The computed Sherwood number is compared to the theory of Levich, 1962.
alpha_c = 1.0/30
r0 = 1.
sherwood(D) = (ci = $2/$5, \
co = $3/$6, \
area = $7, \
dtr = (tr - $2)/0.01, \
tr = $2, \
Sh = dtr/area/(ci*alpha_c - co)*2.*r0/D)
# Levich formula to calculate transfer rate using terminal velocity
levich(ut, D) = 2.*sqrt(D*ut/2./r0/pi)*2.*r0/D
# get the rise velocity averaged for the last steps
array ut[4]
array et[4] = [ 140, 250, 250, 150 ]
do for [i in "1 2 3 4"] {
stats 'log' index 'case '.i u 1:4 every ::et[1] nooutput
ut[i] = STATS_mean_y
}
array D[4] = [ 0.4472, 0.5029, 0.17416, 0.08944 ]
tr = 0.
set key center right
set xlabel 't U/d_0'
set ylabel 'Sh'
set xrange [0.01:4]
set yrange [0:14]
plot for [i = 1:4] \
'log' index 'case '.i u ($1*ut[i]/2./r0):(sherwood(D[i])) \
w l t 'Axi Case '.i lt i, \
for [i = 1:4] levich(ut[i], D[i]) w l t '' lt i dt 2
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 ] |
[levich1962] |
Veniamin Grigorʹevich Levich. Physicochemical hydrodynamics. Prentice-Hall Inc., 1962. |
#include "axi.h"
#include "navier-stokes/centered.h"
#include "two-phase.h"
#include "navier-stokes/conserving.h"
#include "tension.h"
#include "reduced.h"
#include "henry.h"
#include "view.h"
scalar c[], * stracers = {c};
double bubble_radius = 1.;
double box_size = 20.;
double conc_liq1 = 0, conc_gas1 = 1.;
double end_time[4] = {7.5,3,3,2};
int MAXLEVEL, dcase = 1;
int main (int argc, char **argv)
{
size (box_size);
MAXLEVEL = 9;
#if !TREE
N = 1 << MAXLEVEL;
#endif
rho1 = 1.;
rho2 = 0.01;
c.alpha = 1./30.;
TOLERANCE = 1e-4 [*];
for (dcase = 1; dcase <= 4; dcase++) {
switch (dcase) {
case 1:
c.D1 = 0.4472;
f.sigma = 10.0;
G.x = - 2.5;
break;
case 2:
c.D1 = 0.5029;
f.sigma = 10.0;
G.x = - 7.8125;
break;
case 3:
c.D1 = 0.17416;
f.sigma = 1.0;
G.x = - 10.0;
break;
case 4:
c.D1 = 0.08944;
f.sigma = 10.0;
G.x = - 7.8125;
break;
default:
fprintf (stderr, "Error: must specify case\n");
exit (1);
}
mu1 = c.D1;
mu2 = c.D1/20.;
c.D2 = c.D1*100.;
fprintf (stderr, "\n\n# case %d\n", dcase);
run();
}
}
event init (t = 0)
{
#if TREE
refine (sq(2.*bubble_radius) - sq(x - box_size*0.2) - sq(y) > 0 &&
level < MAXLEVEL);
#endif
fraction (f, - (sq(bubble_radius) - sq(x - box_size*0.2) - sq(y)));
foreach()
c[] = conc_liq1*f[] + conc_gas1*(1. - f[]);
}
#if TREE
event adapt (i++) {
adapt_wavelet ({f, c, u}, (double[]){0.01,0.01,0.01,0.01,0.01},
maxlevel = MAXLEVEL);
}
#endif
event extract (t = 0; t += 0.01; t <= end_time[dcase-1])
{
double yb = 0., vb = 0., vbx = 0., area = 0., ci = 0., co = 0.;
foreach (reduction(+:yb) reduction(+:vb) reduction(+:vbx)
reduction(+:ci) reduction(+:co)
reduction(+:area)) {
double dvb = (1. - f[])*dv();
vb += dvb; // volume of the bubble
yb += x*dvb; // location of the bubble
vbx += u.x[]*dvb; // bubble velocity
ci += c[]*(1. - f[])/(f[]*c.alpha + (1. - f[]))*dv();
co += c[]*f[]*c.alpha/(f[]*c.alpha + (1. - f[]))*dv();
if (f[] > 1e-6 && f[] < 1. - 1e-6) {
coord n = interface_normal (point, f), p;
double alpha = plane_alpha (f[], n);
// area of the bubble interface
area += y*pow(Delta, dimension - 1)*plane_area_center (n, alpha, &p);
}
}
if (i == 0)
fprintf (stderr, "t ci co vbx vb vbo area dt\n");
fprintf (stderr,"%g %g %g %g %g %g %g %g\n",
t,
ci*2.*pi, co*2.*pi,
vbx/vb, 2.*pi*vb, 2.*pi*statsf(f).sum, 2.*pi*area,
dt);
}
event pictures (t = end)
{
char name[80];
#if 0
sprintf (name, "dump-%d", dcase);
dump (name);
#endif
double ty[] = { - 0.6, - 0.5, - 0.5, - 0.5};
view (fov = 9, quat = {0.707,0.707,0,0},
ty = ty[dcase - 1],
width = 400, height = 800);
squares ("c", spread = -1, linear = true, map = cool_warm);
draw_vof ("f");
mirror ({0,1}) {
squares ("u.x", spread = -1, linear = true, map = cool_warm);
draw_vof ("f");
}
sprintf (name, "final-%d.png", dcase);
save (name);
}