fs

Free-Surface Boundary Condition

At the time of imposing a Free Surface Boundary Condition, the methodology to be applied depends on whether Monophase or Multiphase simulations are considered.

Monophase simulations

One main feature of SPH is the posibility of carrying out Monophase simulations even in the presence of a Multiphase separation Free-Surface, dramatically reducing the computational cost. In fact, there are a lot of water-air Mulltiphase practical applications, which has been modelled by Monophasic simulations, documented on the literature.

Unfortunatelly, as Colagrossi et al. already demonstrated, consistenty issues are inexorably arising when this approach is considered, unless the fields involved does not vanish at the Free-Surface. To illustrate that, let's depart again form the situation depicted in the following picture:

Let's start considering we have information enough to compute the integrals along both phases, then we can assert that:

$$\langle \nabla f(\mathbf{x}) \rangle = \int_{\bar{\Omega}} f(\mathbf{y}) \nabla W(\mathbf{y} - \mathbf{x}) d\mathbf{y} + \int_{\Omega^{*}} f(\mathbf{y}) \nabla W(\mathbf{y} - \mathbf{x}) d\mathbf{y} = \nabla f(\mathbf{x}) + \mathcal{O}(h).$$

However, we are considering the possibility of neglecting the volume integral along the second phase, $\Omega^{*}$, such that we can write:

$$\int_{\bar{\Omega}} f(\mathbf{y}) \nabla W(\mathbf{y} - \mathbf{x}) d\mathbf{y} = \nabla f(\mathbf{x}) - \int_{\Omega^{*}} f(\mathbf{y}) \nabla W(\mathbf{y} - \mathbf{x}) d\mathbf{y} + \mathcal{O}(h),$$

with

$$\nabla W(\mathbf{y} - \mathbf{x}) d\mathbf{y} = \mathcal{O}\left(\frac{1}{h}\right).$$

Indeed, unless $f(\mathbf{y}) = \mathcal{O}(h)$ or smaller, the differential operator computed in this way will become fully inconsistent, getting bigger errors while the resolution is increased.

To workaround this situation 2 approaches can be considered. Please, notice that in the best case errors of $\mathcal{O}(1)$ will be obtained.

Enforcing small values along the Free-Surface

The first, and most obvious one, is enforcing the field get small values at the free-surface. That's the case of the pressure in CFD applications. Effectively, a null reference pressure can be condidered, $p_0 = 0 \mathrm{Pa}$, because close to the Free-Surface the pressure field takes the form $p(\mathbf{y}) = p_0 + \mathcal{O}(h)$

Symmetrizing the differentail operator

When the previous approach cannot be applied, a second approach is develop taking advantage of the following identity:

$$\int_{\bar{\Omega}} f(\mathbf{x}) \nabla W(\mathbf{y} - \mathbf{x}) d\mathbf{y} + \int_{\Omega^{*}} f(\mathbf{x}) \nabla W(\mathbf{y} - \mathbf{x}) d\mathbf{y} = 0,$$

which allows conveniently rewriting the trucated SPH differential operator,

$$\int_{\bar{\Omega}} \left(f(\mathbf{y}) - f(\mathbf{x})\right) \nabla W(\mathbf{y} - \mathbf{x}) d\mathbf{y} = \nabla f(\mathbf{x}) - \int_{\Omega^{*}} \left(f(\mathbf{y}) - f(\mathbf{x})\right) \nabla W(\mathbf{y} - \mathbf{x}) d\mathbf{y} + \mathcal{O}(h),$$

where this time $\left(f(\mathbf{y}) - f(\mathbf{x})\right) = \mathcal{O}(h)$, provided that the field is a regular enough function.

Multiphase simulations

When all the involved phases plays a main role in the physical problem, or either the truncation errors are unacceptable, then Multiphase simulations can be considered. As has been commented above, the main drawback of Multiphase simulations is the larger computational costs.

In principle, Multiphase simulations do not differs at all from Monophase simulations, consistently extending the integrals to both phases volumes. The most tipical exception is the case of 2 phases with a big density ratio, in which density diffusion should be taken into account. In such a case, a variable kernel length, $h$, should be considered to avoid dramatically large numbers of neighbours, requiring a special formulation as described by Zisis (2017).