A weakly compressible SPH method for violent multi-phase flows with high density ratio

TL;DR

This paper introduces a novel weakly compressible SPH method for simulating violent multi-phase flows with high density ratios, improving efficiency and realism by using a low-dissipation Riemann solver and transport-velocity formulation.

Contribution

It develops a multi-phase WCSPH approach that enhances computational efficiency and eliminates unnatural voids and phase separation in high-density-ratio flows.

Findings

Demonstrates improved robustness in violent water-air flow simulations.

Achieves better accuracy compared to previous methods with similar sound speed.

Significantly increases computational efficiency for multi-phase flow modeling.

Abstract

The weakly compressible SPH (WCSPH) method is known suffering from low computational efficiency, or unnatural voids and unrealistic phase separation when it is applied to simulate highly violent multi-phase flows with high density ratio, such as that between water and air. In this paper, to remedy these issues, we propose a multi-phase WCSPH method based on a low-dissipation Riemann solver and the transport-velocity formulation. The two-phase Riemann problem is first constructed to handle the pairwise interaction between fluid particles, then modified for the fluid-wall interaction to impose the solid wall boundary condition. Since the method uses the same artificial speed of sound for both heavy and light phases, the computational efficiency increases greatly. Furthermore, due to the transport-velocity formulation employed for the light phase and application of the two-phase Riemann problem, the unnatural voids and unrealistic phase separation are effectively eliminated. The method is validated with several 2- and 3D cases involving violent water-air flows. The results have been compared with existing experimental data, previous numerical and analytical solutions, where the proposed method demonstrates good robustness, improved or comparable accuracy, respectively, comparing to previous methods with same choice of sound speed or those with much less computational efficiency.

Equations

$$\frac{d\rho }{dt}=-\rho \nabla \cdot \mathbf{v},$$

$$\frac{d\mathbf{v}}{dt}=-\frac{1}{\rho }\nabla p+\mathbf{g},$$

$$p=c^2(\rho -{\rho }_0),$$

$${\rho }_i=m_i\sum _j{W}_{ij},$$

$$\frac{d{\rho }_i}{dt}={\rho }_i\sum _j\frac{m_j}{{\rho }_j}{\mathbf{v}}_{ij}\cdot {\nabla }_i{W}_{ij}=2{\rho }_i\sum _j\frac{m_j}{{\rho }_j}({\mathbf{v}}_i-{\overline{\mathbf{u}}}_{ij})\cdot {\nabla }_i{W}_{ij},$$

$$\frac{d{\mathbf{v}}_i}{dt}=-\sum _j{m}_j\left(\frac{p_i+p_j}{{\rho }_i{\rho }_j}\right){\nabla }_i{W}_{ij}=-2\sum _j{m}_j\frac{{\overline{p}}_{ij}}{{\rho }_i{\rho }_j}{\nabla }_i{W}_{ij},$$

$$\{\begin{array}{l}({\rho }_L,U_L,p_L)=({\rho }_i,{\mathbf{v}}_i\cdot {\mathbf{e}}_{ij},p_i),\\ ({\rho }_R,U_R,p_R)=({\rho }_j,{\mathbf{v}}_j\cdot {\mathbf{e}}_{ij},p_j),\end{array}$$

$$\{\begin{array}{l}{U}^{\ast }=\overline{U}+\frac{p_L-p_R}{c({\rho }_L+{\rho }_R)},\\ p^{\ast }=\overline{P}+\frac{{\rho }_L{\rho }_R\beta (U_L-U_R)}{{\rho }_L+{\rho }_R},\end{array}$$

$$\frac{\text{d}{\rho }_i}{\text{d}t}=2{\rho }_i\sum _j\frac{m_j}{{\rho }_j}({\mathbf{v}}_i-{\mathbf{v}}^{\ast })\cdot {\nabla }_i{W}_{ij},$$

$$\frac{\text{d}{\mathbf{v}}_i}{\text{d}t}=-2\sum _j{m}_j\frac{p^{\ast }}{{\rho }_i{\rho }_j}{\nabla }_i{W}_{ij},$$

$${\tilde{\mathbf{v}}}_i(t+\delta t)={\mathbf{v}}_i(t)+\delta t\left(\frac{d{\mathbf{v}}_i}{dt}-2p_b\sum _j{m}_j\frac{1}{{\rho }_i{\rho }_j}{\nabla }_i{W}_{ij}\right).$$

$$\frac{d{\mathbf{v}}_i}{dt}=-2\sum _j{m}_j\frac{p^{\ast }}{{\rho }_i{\rho }_j}{\nabla }_i{W}_{ij}+2\sum _j{m}_j\frac{{\overline{\mathbf{A}}}_{ij}}{{\rho }_i{\rho }_j}\cdot {\nabla }_i{W}_{ij}+{\mathbf{g}}_i,$$

$$p^{\ast }=\frac{{\rho }_f{p}_w+{\rho }_w{p}_f}{{\rho }_f+{\rho }_w},$$

$$p_w=\frac{{\sum }_f\frac{p_f}{{\rho }_f}{W}_{wf}+\mathbf{g}{\sum }_f{\mathbf{r}}_{wf}{W}_{wf}}{{\sum }_f\frac{W_{wf}}{{\rho }_f}}.$$

$$\mathbf{v}_{i}^{n+\frac{1}{2}}=\mathbf{v}_{i}^{n}+\frac{\delta t}{2}\big{(}\frac{d\mathbf{v}_{i}}{dt}\big{)}^{n},$$

$${\mathbf{r}}_i^{n+1}={\mathbf{r}}_i^n+\delta t{\mathbf{v}}_i^{n+\frac{1}{2}}.$$

$$\rho_{i}^{n+1}=\rho_{i}^{n}+\delta t\big{(}\frac{d\rho_{i}}{dt}\big{)}^{n+1},$$

$$\mathbf{v}_{i}^{n+1}=\mathbf{v}_{i}^{n+\frac{1}{2}}+\frac{\delta t}{2}\big{(}\frac{d\mathbf{v}_{i}}{dt}\big{)}^{n+1}.$$

$$\delta t\leq 0.25\big{(}\frac{h}{c+U_{max}}\big{)},$$

$$\delta t\le 0.25\sqrt{\frac{h}{\left|\mathbf{g}\right|}}.$$

$$P=\frac{\sum p_f{W}_{sf}{V}_f}{\sum W_{sf}{V}_f+ϵ},$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Multi-phase flow phenomena, where different fluids are separated by phase interfaces, are ubiquitous in many natural and engineering problems. Due to the discontinuous properties and physical interactions across the dynamically evolving phase interface, the study of multi-phase flow problems is a challenging topic in computational fluid dynamics (CFD). Compared to the Eulerian mesh-based methods [1, 2, 3, 4], smoothed particle hydrodynamics (SPH), as a fully Lagrangian and mesh-free method, tracks the phase interfaces with no need for additional interface tracking or capturing algorithms, which are often suffered from serious numerical errors and instabilities. This advantage nominates SPH as a well suited candidate for the simulation of incompressible multi-phase flows and has attracted a great deal of attention in recent decades [5, 6, 7, 8, 9].

The most widely used SPH method for incompressible multi-phase flows is the weakly compressible SPH (WCSPH) method, which was first proposed for simulating single-phase flows with and without free surfaces [10, 11]. In the WCSPH method, the fluid pressure is related to density using an artificial equation of state allowing a small compressibility. When the flow is violent, e.g. inviscid or with high Reynolds number, the standard WCSPH formulation can produce spurious pressure oscillations, thereby may lead to numerical instabilities. To alleviate this issue, several methods, including artificial viscosity [12], density filtering [5] and those based on Riemann solver [13, 14, 15], have been proposed.

As the issue of numerical instabilities become even more serious in the simulation of violent multi-phase flows, especially for those with large density ratio, such as that of water and air, Colagorossi and Landrini [5] further proposed several combined modifications, including Mean Least Square (MLS) interpolation to frequently reinitialize the density field, gradient discretization based on volume (not density) to remedy the density discontinuity issues at the interface, a cohesion term for air to artificial surface-tension effects and XSPH formulation, which corrects the velocity term used in the density and position evolution equations, to regularize particle distribution of the heavy phase. These modifications to stabilize computations not only increase the complexity of the original simple WCSPH formulation, but also their effectiveness relies on a computationally very expensive approach, i.e. assigning the light phase with an artificial speed of sound typically 10 times larger than the heavy phase to obtain a reasonable pressure profile across the interface. Following Colagorossi and Landrini, Mokos et al. [16] have carried out high resolution simulations of violent water-air flows. The simulation results exhibit non-physical voids and phase separation. Such artifacts have also been observed in Gong et al. [17]. To avoid unnatural voids in such simulations a regularization method based on particle shifting [18] was introduced [19]. A further improvement for simulating high Reynolds number flows is also proposed by applying a density diffusive terms [15] for smoother pressure and density fields [20]. However, the very large speed of sound for the light phase is still required to produce acceptable simulation results.

Different from Colagorossi and Landrini [5], Chen et al. [21] proposed another set of combined modifications, including pressure gradient operators imposing continuously induced particle acceleration across the interface [22], a positive background pressure to avoid negative pressure and a similar density reinitialization to obtain a smooth pressure field. Since the same speed of sounds are used for the dense and light fluids, the computational efficiency increases greatly. A considerable drawback of this method is the excessive numerical dissipation, which is able to stabilize simulations but also leads to large discrepancies between their and the reference simulation results, such as those in Refs. [5, 16], especially when the density ratio is high. Zheng and Chen [23] proposed a further modification with a first-order density reinitialization and a local implementation of the artificial viscosity. This modification demonstrates enhanced accuracy, such as shorter delays of the dynamics, due to less numerical dissipation. However, as shown in their only violent water-air flow case, the particles still undergo spurious fragmentation, particularly in the vicinity of phase interface.

In the present study, we propose a simple WCSPH method to simulate highly violent multi-phase, typically liquid-gas, flows with wave breaking and impact events. Our method is based on two modifications of the original WCSPH formulation: one is a two-phase extension from the low-dissipation Riemann solver [24] to realize particle interactions and to stabilize the interface interaction; the other is applying the transport-velocity formulation [25] in the light phase to eliminate non-physical voids or fragmentations. As a result of the two-phase Riemann solver, the light phase experiences the heavy phase like moving wall boundary, while the heavy phase undergoes a free-surface-like flow with variable free-surface pressure. Since we use the same value for the artificial speed of sound in both light and heavy phases, the computational efficiency is much higher than those in Refs. [5, 8, 20]. The low-dissipation nature of the proposed method accurately captures the violent impacts and breaking events, without suffering from excessive numerical dissipation [21] or spurious fragmentation [23]. Note that, we do not use other stabilizing methods to tackle the density discontinuity at the phase interface [26, 7, 27], nor the artificial repulsive pressure force [8, 27]. Furthermore, density reinitialization schemes, deemed to be necessary in the literature [5, 21, 17, 23], are not required. In order to test the accuracy and robustness of the presented method, several 2D and 3D benchmark tests involving highly violent water-air flows with free surfaces are carried out and the results are validated against existing experimental data, previous numerical simulations and analytical solutions. The paper is organized as follows: Section 2 details the proposed multi-phase scheme. The numerical results are then presented and discussed in Section 3. Finally, the concluding remarks of the present study are given in Section 4.

2 Numerical method

For an inviscid and immiscible two-phase, typically liquid-gas, flow with large density ratio, i.e. $\rho_{l}/\rho_{g}\gg 1$, where $\rho_{l}$ and $\rho_{g}$ represent the densities of heavy and light phases, the mass and momentum conservation equations can be written respectively as

[TABLE]

[TABLE]

where $d/dt$ is the material or Lagrangian derivative, $\rho$ is the density, $\mathbf{v}$ is the velocity, $p$ is the pressure and $\mathbf{g}$ is the gravitational acceleration. To close the system, pressure is estimated from density via an artificial equation of state, within the weakly compressible regime. Here, we use a simple linear equation for both the heavy and light phases

[TABLE]

where $c$ and $\rho_{0}$ are the speed of sound and the initial reference density. Here, we assume that the speed of sound is constant and we set $c=10U_{max}$, where $U_{max}$ denotes the maximum anticipated velocity of the flow. Note that, different from Refs. [5, 8, 16], we use the same speed of sound for both heavy and light phases as in Refs. [21, 23].

In the WCSPH formulation, there are two different formulations to implement mass conservation, viz. particle summation and continuity equation [28, 29]. The former calculates density through a summation over all the neighboring particles

[TABLE]

where $m_{i}$ denotes the mass of particle $i$ and the smoothing kernel function $W(\left|\mathbf{r}_{ij}\right|,h)$ is simply substituted by $W_{ij}$, with $\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}$ and $h$ being the displacement vector between particle $i$ and $j$ and the smoothing length, respectively. The latter updates particle density by discretizing the continuity equation

[TABLE]

where $\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$ is the relative velocity and $\mathbf{\overline{u}}_{ij}=(\mathbf{v}_{i}+\mathbf{v}_{j})/2$ the average velocity between particle $i$ and $j$.

Following [28, 29], without taking the artificial viscosity into account, the momentum conservation equation can be discretized as

[TABLE]

with $\overline{{p}}_{ij}=(p_{i}+p_{j})/2$ being the average pressure between particle $i$ and $j$.

2.1 WCSPH based on a multi-phase Riemann solver

For an SPH methods based on Riemann solver [13, 14], an inter-particle Riemann problem is constructed along the unit vector $\mathbf{e}_{ij}=-\frac{\mathbf{r}_{ij}}{\left|\mathbf{r}_{ij}\right|}$ as shown in Fig. 1(a), with the initial left and right states reconstructed from particle $i$ and $j$ by piecewise constant approximation, respectively, as follows

[TABLE]

where subscripts $L$ and $R$ denote left and right states, respectively, and the discontinuity is assumed at $\mathbf{\overline{r}}_{ij}=(\mathbf{r}_{i}+\mathbf{r}_{j})/2$.

The solution of the Riemann problem can be shown using $x-t$ diagrams [30], as in Fig. 1(b). Three waves are present: $W_{L}$,$W^{*}$ and $W_{R}$ respectively at the left, middle and right. The left and right waves are either a shock or a rarefaction, while the middle wave is always a contact discontinuity separating the two intermediate states, namely, ($\rho^{*}_{L},U^{*}_{L},p^{*}_{L}$) and ($\rho^{*}_{R},U^{*}_{R},p^{*}_{R}$). The intermediate states are regarded to satisfy the interface condition $p^{*}=p^{*}_{L}=p^{*}_{R}$ and $U^{*}=U^{*}_{L}=U^{*}_{R}$. For multi-phase flows with the same speed of sound across the interface, following Ref. [31], the intermediate velocity and pressure can be approximated as

[TABLE]

where $\overline{U}=(\rho_{L}U_{L}+\rho_{R}U_{R})/(\rho_{L}+\rho_{R})$ and $\overline{P}=(\rho_{L}p_{R}+\rho_{R}p_{L})/(\rho_{L}+\rho_{R})$. In Eq. (8), each solution for velocity and pressure consists of two terms. The first term is a density-weighted average and the second term is corresponding to numerical dissipation. Due to the assumption of weak compressibility, the magnitudes of the second term (numerical dissipation) are much smaller than that of the first term. Note that, following Zhang et al. [24], a limiter $\beta=\mathrm{min}[3\mathrm{max}(U_{L}-U_{R},0),c]$ is applied for $p^{*}$ in Eq. (8) to reduce the numerical dissipation. Also note that, for a single phase interaction, Eq. (8) recovers the solution of the low-dissipation Riemann solver [24]. For two-phase particle interactions with large density ratio, Eq. (8) gives approximately the intermediate velocity from the heavy phase, which indicates that the light phase experiences the heavy phase as a moving wall boundary, and intermediate pressure from the light phase, which indicates that the heavy phase undergoes a free-surface-like flow with variable free-surface pressure [31]. Having the intermediate values determined, the mass and momentum conservation equations, i.e. Eqs. (5) and (6), can be rewritten as

[TABLE]

[TABLE]

where $\mathbf{v}^{*}=U^{*}\mathbf{e}_{ij}+(\overline{\mathbf{v}}_{ij}-\overline{U}\mathbf{e}_{ij})$. Here, for consistency, one can assign $\overline{\mathbf{v}}_{ij}=(\rho_{i}\mathbf{v}_{i}+\rho_{j}\mathbf{v}_{j})/(\rho_{i}+\rho_{j})$, which is a density-weighted average velocity between particle $i$ and $j$.

2.2 Transport velocity formulation

Due to the tensile instability, WCSPH results in particle clumping and unnatural void regions in fluid dynamics problems [32]. In particular, negative pressure gives rise to particle clustering. Such artifacts become even more serious in multi-phase SPH simulations of water-air flows [5, 19]. Colagrossi and Landrini [5] applied a background pressure to the momentum conservation equation to avoid the tensile instability, while Mokos et al. [19] additionally proposed a particle shifting method to prevent particle clumping.

Here, to prevent the unnatural void regions observed in the light phase, we apply the transport-velocity formulation [33, 34] on the particles of light phase, but not on those of the heavy phase because the pressure on latter generally is kept positive due to the gravity. In transport-velocity formulation, the particle advection velocity $\tilde{\mathbf{v}}$, which is used to advect the particles of light phase, is obtained by using a constant background pressure $p_{b}$ at every time step of size $\delta t$

[TABLE]

Note that, at the interface, water particles are also included in the summation of Eq. (11). With the transport-velocity formulation, the momentum conservation equation for light phase is discretized as

[TABLE]

where $\overline{\mathbf{A}}_{ij}=(\mathbf{A}_{i}+\mathbf{A}_{j})/2$. Here, $\mathbf{A}_{i}=\rho_{i}\mathbf{v}_{i}(\tilde{\mathbf{v}}_{i}-\mathbf{v}_{i})$ is an extra stress tensor as a consequence of advection velocity $\tilde{\mathbf{v}}_{i}$. Note that at the phase interface $\mathbf{A}_{j}$ vanishes for the particles of heavy phase because the transport-velocity formulation is not applied, i.e. $\tilde{\mathbf{v}}_{j}=\mathbf{v}_{j}$. Also note that, following transport-velocity formulation, we employ the density summation form of Eq. (4) for the light phase.

2.3 Wall boundary condition

In the present study, we use fixed dummy particles to impose the solid wall condition as proposed in Ref. [35]. To realize the fluid-wall interactions a Riemann problem is constructed between particles of fluids and wall dummy particles, as for fluid-fluid interactions (see Section 2.1). However, the intermediate pressure value is obtained as

[TABLE]

where subscripts $f$ and $w$ denote fluid and wall, respectively, to decrease the wall-induced numerical dissipation. Note that the intermediate velocity value $U^{*}$ is still obtained via Eq. (8), as for fluid-fluid interactions. Similar to Ref. [35], the pressure of wall dummy particles is calculated by the summation over all contributions of the neighboring fluid particles as

[TABLE]

Note that, by introducing $\rho_{f}$ into the above equation the contribution of heavy phase particles in $p_{w}$ at a triple point, where water, air and solid particles meet, is vanishing. Again, following Ref. [35], the density of wall dummy particles is obtained form pressure via Eq. (3).

2.4 Time integration

To integrate the equations of motion in time, the kick-drift-kick [28, 33] scheme is employed. The first half-step velocity is obtained as

[TABLE]

from which we obtain the position of particles at the next time step

[TABLE]

With the updated flow states, the density of heavy-phase particles at the new time step is then calculated as

[TABLE]

where the time increment of density is approximated via Eq. (5). Note that, the density of light phase particles at the new time step is obtained using the density summation equation, i.e. Eq. (4). Having new densities calculated, the pressure of particles and the time increment of velocity, can be obtained. Finally, the velocity of particles is updated for the new time step as

[TABLE]

For numerical stability, the time-step size is limited by the CFL condition

[TABLE]

and the body force condition

[TABLE]

3 Numerical results

In this section several test cases involving two-phase dam break and sloshing are considered to validate the proposed method for modeling violent flows with high density ratio. The 5th-order Wendland kernel [36] with a smoothing length of $h=1.3dx$, where $dx$ is the initial particle spacing, and a support radius of $2.6dx$ are used. In all cases, the maximum velocity is estimated as $U_{max}=2\sqrt{gH}$, where $H$ is the initial water depth, according to the shallow-water theory [37], for setting the speed of sound $c_{0}$. Here, we set the background pressure $p_{b}=4\rho_{0}c^{2}_{0}$ for the transport-velocity formulation in all of the test cases.

3.1 Two-phase dam break

We consider a two-phase dam break case, which has been experimentally studied by Zhou et al. [38]. Owing to its complex flow involving wave impact and breaking events, dam break is widely used to validate SPH-based methods [5, 35, 39, 16, 21]. Note that, in the two-phase model, the cushioning effect of entrapped air pockets during violent hydrodynamics can influence the impact loads.

The schematic of this problem is illustrated in Fig. 2.

A water column of height $H=1$ and width of $2H$ is located in a tank of length $5.366H$, and is surrounded by air. We consider an inviscid flow, where the density of water and air are set to $\rho_{w}=1$ and $\rho_{a}=0.001$, respectively. Using $g=1$ in negative $y$-direction as the gravity, all quantities correspond to their non-dimensional variables. Note that for simplicity, water with zero initial pressure is released immediately when the computation starts instead of being released from a gate holding the pressure-relaxed water as in the experimental setup. To record the pressure signals, a sensor $P$ is located on the downstream wall at $y/H=0.19$ from the bottom wall, which is slightly different from the experiments. Note that this adjustment gives better agreement between SPH and experimental results due to several uncertainties in the measurements [40]. In the simulation, the measured pressure is obtained by averaging the values from particles within a support radius $2.6dx$ with

[TABLE]

where subscripts $s$ and $f$ denote solid wall and fluids, respectively, $V_{f}=m_{f}/\rho_{f}$ is the volume of a fluid particle and $\epsilon=1.0\times 10^{-15}$ is a small positive number to avoid dividing by zero. Furthermore, two height sensors $h_{1}$ and $h_{2}$ are located on the bottom bed at $x_{1}/H=0.825$ and $x_{2}/H=1.653$ from the right wall, respectively, to monitor the time evolution of the water level.

Fig. 3 illustrates several snapshots during the time evolution.

It can be observed that the main features of the flow, viz. high roll-up along the downstream wall, a large reflected jet and a water surface breaking due to the re-entry of the backward wave, are well captured by the present method, as in previous single-fluid numerical simulations (see e.g. Refs. [41, 24, 42]). Moreover, the water-air interface is sharply maintained as in previous two-phase simulations (see e.g. Refs. [5, 43]). Note that, due to the transport-velocity formulation applied in air phase and the application of two-phase Riemann problem, no obvious unnatural voids or unrealistic phase separation appears.

To evaluate the accuracy of the captured phase interface, as shown in Fig. 4, the present results are compared with those obtained by the boundary element method (BEM) in Ref. [44], the level-set method in Refs.[45, 5], and previous two-phase SPH methods, namely, Refs. [5], [21] and [23]. All the three previous SPH simulations are with a spatial resolution of $H/dx=60$, which is comparable to the spatial resolution of our simulation, i.e. $H/dx=80$.

As it can be observed, a good agreement is achieved between the predicted interface location by the presented method herein, and the BEM and level-set predictions for the smooth part of the solutions. However, due to the remarkable physical uncertainties of the breaking waves, there exist discrepancies. On the other hand, as it is also observed in Ref. [44], the BEM method leads to excessive numerical dissipation and its results show considerable time delay.

The following comparison is noteworthy between the present method and the two classes of the previous SPH simulations, viz. Refs. [5] and [21, 23]: one can observe that the present method accurately predicts a sharply maintained phase interface (Fig. 4(c)), while is computationally much more efficient than the method presented in Ref. [5], due to algorithm simplicity and the same speed of sound for both the heavy and light phases, while there exist no unnatural void region as observed in Refs. [21, 23] (see Figs. 4(a) and (b)). Albeit the methods presented in Refs. [21, 23] are computationally efficient (also due to the same speed of sound for both the heavy and light phases), they suffer from spurious fragmentation at the phase interface, as well as non-physical void regions. As can be seen in Figs. 4(a) and (b), the excessive numerical dissipation exhibited in Ref. [21] (using a global artificial viscosity) is, to some degree, alleviated by the method presented in Ref. [23] (using a local artificial viscosity), however, the phase interface undergoes more severe fragmentations in the latter (see the figures in their Section 4.5).

Time history of the recorded impact pressure signals on the vertical downstream wall are compared against experimental data in Fig. 5. As expected, in a weakly compressible framework, the pressure profile exhibits high frequency oscillations, which are mitigated with increased spatial resolution. With the fairly low resolution in mind, a good agreement in the overall behavior of the pressure is noted.

The predicted water wave front with increasing particle resolution is compared against the experimental data provided by Refs. [46] and [47] in Fig. 6, as well as the analytical solution from shallow water theory [37]. The results are also compared with a previous single-phase SPH simulation [24]. Due to the consideration of air, compared to the previous single-phase SPH results, the present two-phase simulation predicts slightly slower wave front at time instants after $t\sqrt{g/H}\approx 2.0$, which suggest a closer agreement with experimental data. However, as observed in several previous works [5, 41, 48], a slower front wave propagation is observed in experiments, which is attributed to several uncertainties, viz. accuracy of the experiments, turbulent wall boundary layer and wall roughness. The asymptotic convergence of the predicted water wave front to the theoretical shallow water solution with increased spatial resolution can also be noted on a long time scale. Note that the shallow water assumption does not hold at initial time instants of the problem.

Fig. 7 plots the time history of the recorded water level signals at sensors $h_{1}$ and $h_{2}$ in comparison with experimental data and single-phase SPH predictions.

As mentioned in Ref. [5], since time evolution of the water level is determined by the wave front shape, initial conditions and the bottom bed roughness are influential. The differences between simulation results and the experiments are thus reasonable due to the present inviscid model and different initial conditions. Similar discrepancies are also observed in previous numerical studies [5, 41]. Due to the slightly faster wave front predicted by the inviscid simulation (see Fig. 6), higher water levels are recorded. However, the agreement between present numerical results and the experiments is satisfactory. In particular, the time evolution of water level on a longer time scale ($t\sqrt{g/H}\approx$ 7–10) is evidently in good agreement with experiments in higher spatial resolutions. Note that, the presently predicted $h_{1}$ signals and those obtained by the single-phase SPH method of Zhang et al. [24] are in a fair agreement. Nevertheless, the comparison for $h_{2}$ reveals the effect of air particles, as the peak value of $h_{2}$ does not experience an abrupt jump in the present results and is thus in a closer agreement with the measurements.

3.2 Two-phase liquid sloshing

Motion of a liquid which partially fills a tank due to the periodic motion of the latter is known as liquid sloshing. If the motion frequency is close to the natural frequency of the liquid due to gravity wave, resonant condition happens, which may cause enormous impact loads on the structures. Due to the highly non-linear phenomena in sloshing, this benchmark case has been used to validate many single-phase SPH methods (see e.g. [49, 50, 42]), as well as few two-phase counterparts [19]. In this section we consider a two-phase liquid sloshing case, which has been experimentally studied by Rafiee et al. [49].

The schematic of this problem is illustrated in Fig. 8.

A tank of height $D=0.9$ m and width of $L=1.3$ m is partially filled with water to 20% of the height, i.e. $H=0.18$ m, while the remainder is filled by air. We consider an inviscid flow, where the density of water and air are set to $\rho_{w}=1000~{}\mathrm{kg.m^{-3}}$ and $\rho_{a}=1~{}\mathrm{kg.m^{-3}}$, respectively. The gravity is $g=9.81~{}\mathrm{m.s^{-2}}$ in negative $y$-direction. The tank motion is defined by a sinusoidal excitation of $x=A_{0}sin(2.0f_{0}\pi t)$, where $A_{0}=0.1~{}m$ and $f_{0}=0.496~{}s^{-1}$ are the amplitude and frequency, respectively. $f_{0}$ is close to the natural frequency of water, thereby resulting in resonant condition. To record the pressure signals, three sensors $P_{1}$, $P_{2}$ and $P_{3}$ are located on the right wall at 0.165 m, 0.18 m and 0.195 m from the bottom, respectively. The same method is utilized to record pressure signals as explained in Section 3.1.

Fig. 9 illustrates several snapshots during the time evolution. It can be observed that the highly non-linear features of the flow, viz. strong impacts on the walls and largely deformed free surfaces due to re-entry of the backward waves, are well captured by the present method, as shown in the experimental and single-phase simulation results in Rafiee et al. [49].

Moreover, the phase interface is sharply maintained between air and water, however, air particles are trapped in water to form air bubbles due to highly violent hydrodynamic events. Again, evidently no unnatural voids or unrealistic phase separation appears. Fig. 10 shows the zoom-in view of the water surface breaking, entrapped air pockets and the pressure field at two different time instants from Fig. 9.

As it can be observed, the pressure field is quite smooth, while differences in its gradient at the phase interface are notable. Such near-interface pressure distribution is similar to the results obtained in Ref. [21] (their Fig. 19) where a pressure-gradient discontinuity is imposed at phase interface [22]. Note that, the air pockets entrapped during the impact events are well captured, being in agreement with experimental observations of Rafiee et al. [49].

As shown in Fig. 11, the time history of the recorded impact pressure signals on the right vertical wall are compared against the experimental data [49] and the simulation results obtained by a single-phase WCSPH method [24].

A good agreement in the overall behavior of the pressure is noted in the prediction of both time and value, with lower peak pressure values compared to those from experiments. Furthermore, in contrast to the single-phase SPH simulation, the peak pressure values predicted by the present method are much lower, implying the cushioning effect of entrapped air pockets, which is excluded in single-phase SPH simulations and can lead to highly overestimated peak pressure values [50].

3.3 3D two-phase dam break

In order to demonstrate the capabilities of the presented method to simulate large scale problems in 3D, the parallel computational framework proposed in Ref. [51] is adopted to carry out 3D simulations of the two-phase dam break problem. The configuration is identical to the case in Section 3.1 with the consideration of a third dimension of length $H$. To discretize the computational domain, particles are initially located on a regular grid with a particle spacing of $dx=H/40$.

Fig. 12 illustrates several snapshots during the time evolution of the free surface. Similar to the presented results in Section 3.1, the main features of the dam break problem are well captured also here.

Evidently, as shown in Fig. 12, the pressure field is quite smooth, while large gradient is observed near the phase interface. Again, the entrapped air pockets are well captured during the impact events and the phase interface is sharply maintained with no spurious fragmentation, and no particle penetration is observed.

4 Conclusions

In this paper, we have presented a multi-phase SPH method for the simulation of highly violent flows with large density ratio. With the particle interactions being handled by a multi-phase Riemann solver, the scheme is simple and robust. The method employs the same values for the artificial speed of sound in both heavy and light phases. This advantage, allows for larger time steps, which greatly improves the computational performance. The proposed method is validated against existing experimental data, previous numerical studies and analytical solutions for several water-air flow problems. It is observed that stable solutions with sharp interfaces are obtained without suffering from non-physical void regions which is observed in previous SPH simulations, especially in those using the same choice of artificial sound speeds. It is also observed that the accuracy of the present method is comparable to the previous less efficient methods using a much larger speed of sound in the light phase. Furthermore, the present method replicates the cushioning effect of entrapped air pockets, leading to more reasonable pressure peaks. Albeit the present work focuses on benchmark test cases, the method is sufficiently generic to be employed in more complex scientific and industrial problems, which is the main scope of future work.

Acknowledgments

The authors gratefully acknowledge the financial support by Deutsche Forschungsgemeinschaft (DFG HU1572/10-1, DFG HU1527/12-1) for the present work.