sandbox/ecipriano/test/interfaceregression.c

    Interface Regression Velocity

    This module test the capability of the interface regression model implemented in evaporation.h. The interface gradients are computed at cell centeres, and so is the vaporization rate. Therefore, interpolations are needed to transform the vaporization rate in an interface regression velocity, defined on cell faces, and that allows the interface to be transported by a modified velocity. This velocity is not conservative, but the “amount” of non conservation reflects the amount of reference phase to be added or removed by the phase change phenomena.

    PROS: This approach works well also when the velocity field is zero everywhere, and the only contribution to the interface motion is the phase change (i.e., isolated droplet with zero surface tension). While these conditions are ideal, they can be exploited to decouple the effect of phase change from other physical phenomena. Apart from these numerical conditions, this method limits the possibility of under- and over-shoots, that easily arise when applying the phase change term in the VOF transport equation as an explicit source. In that case, the use of a redistribution algorithm is crucial.

    CONS: The interface regression velocity is distributed in a manner that guarantees the correct consumption of the reference phase (according to the material balance). This does not mean that the divergence of the velocity used is equal to the vaporization rate, localized at the interface. Therefore, it cannot be used for the transport of tracers associated with the vof field, using an advection equation in conservative form (e.g., conserved tracers in the all-mach solver). The transport of tracers in non-conservative form, as currently implemented in vof.h, works fine.

    Redistribution Procedure

    The vaporization rate, in each interfacial cell, is distributed among the faces in contact with that interfacial cell. There are two possible cases:

    1. If the face connects an interfacial cell with a pure cell, the interface regression velocity on that face is taken from the interfacial cell, weighted by the corresponding normal component.

    2. If the face connects two interfacial cells, the interface regression velocity is computed from a linear interpolation between the two consecutive vaporization rates, weighted by the interface normal.

    Interface Regression Velocity Distribution

    Interface Regression Velocity Distribution

    Phase Change Setup

    We suppress the expansion term in the continuity equation, it is not relevant for this test. Therefore, we do not need a method to compute the extended velocity. For this reason, after the declaration of the field \mathbf{u}_f we set the extended velocity to be equal to this one-field face velocity.

    #define DIFFUSIVE

#include "grid/multigrid.h"
#include "navier-stokes/centered-evaporation.h"
#define ufext uf
#include "two-phase.h"
#include "evaporation.h"
#include "fixedflux.h"
#include "view.h"

    Simulation Setup

    We set the value of the vaporization rate per unit of interface surface. We also declare the index of the simulation case (we run 3 different cases).

    double mEvapVal = -0.02;
int sim = 0;

int main (void) {
  origin (-0.5, -0.5);
  DT = 1.e-2;
  for (sim=0; sim<3; sim++) {
    init_grid (1 << 6);
    run();
  }
}

#define circle(x,y,R)(sq(R) - sq(x) - sq(y))
#define iplane(x,y)(x-y+1.e-5)
#define ellipse(x,y,R)(1 - sq(x/(1.2*R)) - sq(y/(0.8*R)))

event init (i = 0) {
  switch (sim) {
    case 0: fraction (f, circle(x,y,0.23)); break;
    case 1: fraction (f, iplane(x,y)); break;
    case 2: fraction (f, ellipse(x,y,0.3)); break;
  }
}

    Post-Processing

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

    Compute Divergence

    We compute the divergence of the modified velocity field, used by the vof advection. The interface regression velocity is computed in the vof event. Therefore, we compute the divergence in the tracer_advection event, which is executed right after vof.

    scalar divu[];

event tracer_advection (i++) {
  foreach() {
    divu[] = 0.;
    foreach_dimension()
      divu[] += ((uf.x[1] - vpc.x[1]) - (uf.x[] - vpc.x[]));
    divu[] /= Delta;
  }
}

    Time-Derivative of the Volume Fraction

    We want to check that, even if we distribute the vaporization rate, we obtain a variation of the volume fraction which is coherent with the cell-centered vaporization rate:

    \displaystyle \dfrac{\partial f}{\partial t} = \dfrac{f^{t+1} - f^t}{\Delta t} = \dfrac{\dot{m}}{\rho_l}

    To compare the two quantities, we compute the volume integrals of the vaporization rate and of the vof fraction time derivative.

    First, we store the old time before the vof advection.

    scalar fold[];

event vof (i++) {
  foreach()
    fold[] = f[];
}

    Then, we compute the volume integrals.

    scalar mEvapEff[], mEvapVol[];

event tracer_advection (i++) {
  foreach() {
    mEvapEff[] = (f[] - fold[])/dt;
    mEvapEff[] = fabs (mEvapEff[]) > F_ERR ? mEvapEff[] : 0.;
    vofrecon vr = vof_reconstruction (point, f);
    mEvapVol[] = mEvapTot[]*vr.dirac;
  }
  double totmEvapEff = rho1*statsf(mEvapEff).sum;
  double totmEvapVol = rho1*statsf(mEvapVol).sum;

  fprintf (stderr, "%g %g %g\n", t, totmEvapEff, totmEvapVol);
}

event closefile (t = end,last) {
  fprintf (stderr, "\n\n");
}

    Movie

    We write the animation with the divergence of the velocity field and the interface facets.

    void write_movie (char * name) {
  clear();
  draw_vof ("f", lw = 3);
  squares ("divu", spread=-1);
  cells();
  save (name);
}

event movie (t += 0.1; t <= 10) {
  switch (sim) {
    case 0: write_movie ("case1.mp4"); break;
    case 1: write_movie ("case2.mp4"); break;
    case 2: write_movie ("case3.mp4"); break;
  }
}

    Results

    The animations show the divergence of the velocity that transports the volume fraction.

    Sphere

    Plane

    Ellipse

    The following plot shows the comparison between the variation of liquid volume according to the vaporization rate, and the actual variation from the interpolation presented in this module.

    reset
set xlabel "time [s]"
set ylabel "Vaporization Rate [kg/s]"
set size square
set key bottom right
set grid

plot "log" index 0 u 1:2 w l lw 2 t "Time Derivative Sphere", \
     "log" index 0 u 1:3 w l lw 2 t "Vaporization Rate Sphere", \
     "log" index 1 u 1:2 w l lw 2 t "Time Derivative Plane", \
     "log" index 1 u 1:3 w l lw 2 t "Vaporization Rate Plane", \
     "log" index 2 u 1:2 w l lw 2 t "Time Derivative Ellipse", \
     "log" index 2 u 1:3 w l lw 2 t "Vaporization Rate Ellipse"
    Comparison between time derivative and vaporization rate (script)

    Comparison between time derivative and vaporization rate (script)