# 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

[malan2021geometric] |
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. |

[scapin2020volume] |
Nicolò Scapin, Pedro Costa, and Luca Brandt. A volume-of-fluid method for interface-resolved simulations of phase-changing two-fluid flows. |

[palmore2019volume] |
John Palmore Jr and Olivier Desjardins. A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows. |