Velocity Potential Model

    We want to obtain a divergence-free velocity extrapolation from the field velocity with phase change. A possible approach, proposed in different works with small differences (Scapin et al. 2020, Malan et al. 2021, Palmore et al. 2019), consists in solving an additional Poisson equation that allows a velocity potential \phi to be computed: \displaystyle \nabla \cdot \left( \alpha \nabla \phi \right) = \dfrac{\dot{m}}{\Delta t} \left( \frac{1}{\rho_g} - \frac{1}{\rho_l} \right) \delta_\Gamma \\ The stefan velocity can be obtained from the velocity potential as: \displaystyle \mathbf{u}^S = -\Delta t \alpha \nabla \phi We then calculate the extended velocity by subtracting the stefan velocity from the field velocity. The resulting extended velocity field will be divergence-free by construction: \displaystyle \mathbf{u}^E = \mathbf{u} - \mathbf{u}^S


    We initialize the velocity potential field ps, the stefan velocity ufs, and the divergence-free extended velocity ufext.

    scalar ps[];
    face vector ufs[], ufext[];
    mgstats mgpsf;

    The volume expansion term is declared in evaporation.h.

    extern scalar stefanflow;

    Helper functions

    We define the function that performs the projection of the stefan velocity onto a field with divergence equal to the volume expansion term.

    mgstats project_sv (face vector uf, scalar p,
         (const) face vector alpha = unityf,
         double dt = 1.,
         int nrelax = 4)
      scalar div[];
        div[] = stefanflow[]/dt;
      mgstats mgp = poisson (ps, div, alpha,
          tolerance = TOLERANCE/sq(dt), nrelax = nrelax);
        ufs.x[] = -dt*alpha.x[]*face_gradient_x (ps, 0);
      return mgp;

    Boundary conditions

    It is important to impose for ps the same boundary conditions and the same tolerance used in the Poisson equation for the pressure p.

    ps[right] = neumann (neumann_pressure(ghost));
    ps[left]  = neumann (- neumann_pressure(0));
    #if AXI
    ufs.n[bottom] = 0.;
    ufs.t[bottom] = dirichlet(0);
    ps[top]    = neumann (neumann_pressure(ghost));
    #else // !AXI
    #  if dimension > 1
    ps[top]    = neumann (neumann_pressure(ghost));
    ps[bottom] = neumann (- neumann_pressure(0));
    #  endif
    #  if dimension > 2
    ps[front]  = neumann (neumann_pressure(ghost));
    ps[back]   = neumann (- neumann_pressure(0));
    #  endif
    #endif // !AXI

    Extended velocity

    We perform the projection of the stefan velocity by overloading the event end_timestep, defined in navier-stokes/centered.h.

    event end_timestep (i++, last)

    We set to zero the face-centered stefan velocity ufs and the velocity potential ps.

        ufs.x[] = 0.;
        ps[] = 0.;

    We solve the Poisson equation using the multigrid solver.

      mgpsf = project_sv (ufs, ps, alpha, dt, mgpsf.nrelax);

    We compute a divergence-free extended velocity by subtracting the stefan velocity from the field velocity.

        ufext.x[] = uf.x[] - ufs.x[];


    This approach works fine when the field velocity is larger than the velocity due to the phase change. An example is a falling droplet, or a droplet in forced convective conditions. Static droplets, evaporating in reduced gravity conditions, can be simulated using this method just for small vaporization rates or for small density ratio values. If these conditions are not met, the double pressure velocity coupling approach should be preferred.



    LC Malan, Arnaud G Malan, Stéphane Zaleski, and PG Rousseau. A geometric vof method for interface resolved phase change and conservative thermal energy advection. Journal of Computational Physics, 426:109920, 2021.


    Nicolò Scapin, Pedro Costa, and Luca Brandt. A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows. Journal of Computational Physics, 407:109251, 2020.


    John Palmore Jr and Olivier Desjardins. A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows. Journal of Computational Physics, 399:108954, 2019.