sandbox/ecipriano/src/navier-stokes/velocity-potential.h

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

Fields

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.

trace
mgstats project_sv (face vector uf, scalar p,
     (const) face vector alpha = unityf,
     double dt = 1.,
     int nrelax = 4)
{
  scalar div[];
  foreach()
    div[] = stefanflow[]/dt;

  mgstats mgp = poisson (ps, div, alpha,
      tolerance = TOLERANCE/sq(dt), nrelax = nrelax);

  foreach_face()
    ufs.x[] = -dt*alpha.x[]*face_gradient_x (ps, 0);
  boundary((scalar*){ufs});

  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.

  foreach_face()
    ufs.x[] = 0.;

  foreach()
    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.

  foreach_face()
    ufext.x[] = uf.x[] - ufs.x[];
}

Notes

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.

References