fixedbubblevelocity.c

Fixed Flux Bubble Expansion

Fixed Flux Bubble Expansion

Expansion of a bubble with a fixed vaporization rate. This test case was proposed by Tanguy et al., 2014 and it is interesting because it can be used to test the velocity profile with the Stefan flow, and how the velocity changes when the volume expansion term in the continuity equation is shifted and/or distributed.

We perform this test case using the Velocity Potential approach, and we compare the results both with the theoretical velocity profile and with the numerical results obtained by Tanguy et al., 2014. The same test is repeated using the Velocity Extrapolation approach, providing a better and smoother velocity profile for the advection of the volume fraction, at the cost of a higher computational time.

The animation shows the growth of the bubble in time.

Evolution of the interface position and the volume fraction field

Phase Change Setup

If the Velocity Potential approach is used, we divide the mass balance by the density of the gas phase, and we avoid the Stefan flow shifting procedure. If Velocity Extrapolation is used, we divide the mass balance by the liquid phase density, and we shift the expansion term toward the gas phase.

#if EXTRAPOLATION
# define INT_USE_UF
# define CONSISTENTPHASE2
# define SHIFT_TO_GAS
#else
# define NOSHIFTING
# define BYRHOGAS
#endif

Simulation Setup

We use the centered Navier-Stokes equations solver with the evaporation source term in the projection step. The evporation model is combined with the fixed flux mechanism, that imposes a constant vaporization flowrate.

#include "grid/multigrid.h"
#include "navier-stokes/centered-evaporation.h"
#if EXTRAPOLATION
# include "navier-stokes/velocity-extrapolation.h"
# define ufext ufext1
#else
# include "navier-stokes/velocity-potential.h"
#endif
#include "two-phase.h"
#include "tension.h"
#include "evaporation.h"
#include "fixedflux.h"
#include "view.h"

Boundary Conditions

Outflow boundary conditions are imposed on every domain boundary since the bubble is initially placed at the center of the domain. We set the BCs also for the velocity potential ps.

u.n[top] = neumann (0.);
u.t[top] = neumann (0.);
p[top] = dirichlet (0.);
#if EXTRAPOLATION
ps1[top] = dirichlet (0.);
ps2[top] = dirichlet (0.);
#else
ps[top] = dirichlet (0.);
#endif

u.n[right] = neumann (0.);
u.t[right] = neumann (0.);
p[right] = dirichlet (0.);
#if EXTRAPOLATION
ps1[right] = dirichlet (0.);
ps2[right] = dirichlet (0.);
#else
ps[right] = dirichlet (0.);
#endif

u.n[left] = neumann (0.);
u.t[left] = neumann (0.);
p[left] = dirichlet (0.);
#if EXTRAPOLATION
ps1[left] = dirichlet (0.);
ps2[left] = dirichlet (0.);
#else
ps[left] = dirichlet (0.);
#endif

u.n[bottom] = neumann (0.);
u.t[bottom] = neumann (0.);
p[bottom] = dirichlet (0.);
#if EXTRAPOLATION
ps1[bottom] = dirichlet (0.);
ps2[bottom] = dirichlet (0.);
#else
ps[bottom] = dirichlet (0.);
#endif

Model Data

We set the initial radius of the bubble, the value of the vaporization rate, and the maximum level of refinement.

int maxlevel = 7;
double D0 = 0.002, mEvapVal = -0.1;
double effective_radius0;


int main (void) {

We set the density and viscosity values.

  rho1 = 1000., rho2 = 1.;
  mu1 = 0.001, mu2 = 1.78e-5;

We set the surface tension coefficient.

  f.sigma = 0.07;

We change the dimension of the domain and we shift the origin.

  L0 = 0.008;
  origin (-0.5*L0, -0.5*L0);

We modify the number of layers used in velocity-extrapolation.h in order to extrapolate the gas phase velocity from the second layer of cells close to the interface. This is necessary if the expansion term is shifted toward the gas phase.

#if EXTRAPOLATION
  nl1 = 2;
#endif

We setup the grid and run the simulation.

  init_grid (1 << maxlevel);
  run();
}

#define circle(x,y,R)(sq(R) - sq(x) - sq(y))

We initialize the volume fraction field, and from that field we compute the numerical initial radius.

event init (i = 0) {
  fraction (f, -circle(x,y,0.5*D0));

  double gasvol = 0.;
  foreach(reduction(+:gasvol))
    gasvol += (1. - f[])*dv();
  effective_radius0 = sqrt (gasvol/pi);
}

We adapt the grid according to the velocity and the volume fraction fields.

#if TREE
event adapt (i++) {
  adapt_wavelet_leave_interface ({u.x,u.y}, {f},
      (double[]){1.e-3,1.e-3}, maxlevel, 1);
}
#endif

Post-Processing

The following lines of code are for post-processing purposes.

Exact Solution

We write a function that computes the analytic solution for the bubble radius.

double exact (double t) {
  return effective_radius0 - mEvapVal/rho2*t;
}

Output Files

We write the bubble radius from the numerical simulation, the exact radius and the relative error.

event output_data (i++) {
  char name[80];
  sprintf (name, "OutputData-%d", maxlevel);
  static FILE * fp = fopen (name, "w");

  double gasvol = 0.;
  foreach(reduction(+:gasvol))
    gasvol += (1. - f[])*dv();
  double effective_radius = sqrt (gasvol/pi);
  double relerr = fabs (exact(t) - effective_radius) / exact(t);
  fprintf (fp, "%g %g %g %g\n", t, effective_radius, exact(t), relerr);
  fflush (fp);
}

Velocity Profile

We write on a file the velocity profile at a specific time step.

event profiles (t = 0.005,last) {
  char name[80];
  sprintf (name, "Profiles-%d", maxlevel);
  FILE * fp = fopen (name, "w");

We create an array with the mass fraction profiles for each processor.

  Array * arruf = array_new();
  Array * arrus = array_new();
  for (double x = 0.; x < 0.5*L0; x += 0.05*L0/(1 << maxlevel)) {
    double val_uf = interpolate (uf.x,  x, 0.);
#if EXTRAPOLATION
    double val_us = 0.;
#else
    double val_us = interpolate (ufs.x, x, 0.);
#endif
    val_uf = (val_uf == nodata) ? 0. : val_uf;
    val_us = (val_us == nodata) ? 0. : val_us;
    array_append (arruf, &val_uf, sizeof(double));
    array_append (arrus, &val_us, sizeof(double));
  }
  double * val_ufp = (double *)arruf->p;
  double * val_usp = (double *)arrus->p;

We sum each element of the arrays in every processor.

  @if _MPI
  int size = arruf->len/sizeof(double);
  MPI_Allreduce (MPI_IN_PLACE, val_ufp, size, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  MPI_Allreduce (MPI_IN_PLACE, val_usp, size, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
  @endif

The master node writes the profiles on a file.

  if (pid() == 0) {
    int count = 0;
    for (double x = 0.; x < 0.5*L0; x += 0.05*L0/(1 << maxlevel)) {
      fprintf (fp, "%g %g %g\n", x, val_ufp[count], val_usp[count]);
      count++;
    }
    fprintf (fp, "\n\n");
    fflush (fp);
  }
  fclose (fp);
  array_free (arruf);
  array_free (arrus);
}

Movie

We write the animation with the evolution of the volume fraction field and the gas-liquid interface.

event movie (t += 0.0001; t <= 0.01) {
  clear();
  draw_vof ("f", lw = 1.5);
  squares ("f", min = 0., max = 1., linear = true);
  save ("movie.mp4");

#if EXTRAPOLATION
  clear();
  draw_vof ("f", lw = 1.5);
  squares ("uext1.x", spread = -1, linear = true);
  vectors ("uext1.x", scale = 1.e-3, lc = {1.,1.,1.});
  save ("movie-uext1.mp4");

  clear();
  draw_vof ("f", lw = 1.5);
  squares ("uext2.x", spread = -1, linear = true);
  save ("movie-uext2.mp4");
#endif
}

Results

We compare the velocity profile with the theretical solution and with the results obtained by Tanguy et al., 2014 (Fig. 5 (b), (c)).

reset
set xlabel "radius [m]"
set ylabel "Velocity [m/s]"
set key top right
set size square
set grid

plot "../data/tanguy-fixedbubblevelocity-theoretical.csv" u 1:2 w p pt 6 t "Theoretical", \
     "../data/tanguy-fixedbubblevelocity-b.csv" u 1:2 w p pt 6 t "Tanguy et al., 2014 b", \
     "../data/tanguy-fixedbubblevelocity-c.csv" u 1:2 w p pt 6 t "Tanguy et al., 2014 c", \
     "Profiles-7" u 1:2 w l lw 1.2 t "Field Velocity (potential)", \
     "../fixedbubblevelocity-extrapolation/Profiles-7" u 1:2 w l lw 1.2 t "Field Velocity (extrapolation)"

Radial velocity profile (script)

reset
set xlabel "t [s]"
set ylabel "Bubble Radius [m]"
set key top left
set size square
set grid

plot "OutputData-7" u 1:2 w l lw 2 t "LEVEL 7 potential", \
     "../fixedbubblevelocity-extrapolation/OutputData-7" u 1:2 w l lw 2 t "LEVEL 7 extrapolation", \
     "OutputData-7" every 20 u 1:3 w p pt 8 ps 1.5 t "Analytical"

Evolution of the radius in time (script)

Extrapolated Velocity Fields

The following animations show the evolution of the extrapolated velocities ufext1.x and ufext2.x. The liquid phase velocity is used for the transport of the volume fraction. We notice that the extrapolation procedure leads to a velocity field which is smooth and continuous across the interface. The extrapolation procedure is expensive, and for this reason it is performed just for a few steps, without the need to extrapolate the velocity field over the whole domain. We just need a few cells across the interface.

Evolution of the liquid extrapolated velocity (x component)

Evolution of the gas extrapolated velocity (x component)

References

[tanguy2014benchmarks]

Sébastien Tanguy, Michaël Sagan, Benjamin Lalanne, Frédéric Couderc, and Catherine Colin. Benchmarks and numerical methods for the simulation of boiling flows. Journal of Computational Physics, 264:1–22, 2014.