I Hammani

S Marrone

A Colagrossi
email: andrea.colagrossi@cnr.it

G Oger

D Le Touzé

Detailed study on the extension of the δ-SPH model to multi-phase flow

Keywords: Smoothed Particle Hydrodynamics, multi-phase flows, water impact, air-cushion effects, δ-SPH model

, allows to improve the evaluation of the pressure field, removing numerical noise and improving also the particles spatial distribution.

The time stepping and the choice of the speeds of sound for the different phases are discussed, showing that this choice is driven not only by physical consideration but also by numerical constraints linked to the stability of the scheme. To this aim, comparisons are provided to a Riemann-SPH multi-phase model. In particular we show that the proposed δ-SPH multi-phase solver has a different stability region than its Riemann-SPH counterpart, allowing for bigger time steps for some density and speed of sound ratios, while in other conditions the Riemann-SPH multi-phase model is more convenient.

A series of validating tests are carried out over different benchmarks widely used in the SPH literature. As a final test-case the water entry of a corrugated panel involving the entrapment of an air cavity is considered to show how the proposed multi-phase δ-SPH method is able to accurately treat complex water impact events.

Introduction

Multi-phase flows can be very important for many engineering applications (mixing/separation devices, engines, propellers, etc.). The problem of simulating multi-phase flows has been the subject of many researches in the recent years due to the growing industrial demands, in the maritime, coastal and aeronautical engineering sectors. For the simulation of these flows, a multi-phase solver is required. However, because of the complex interfaces dynamics between the phases and of the numerical constraints, the derivation of robust models is not straightforward. In this context, the Smoothed Particle Hydrodynamics (SPH) method has proven to be a powerful method. Indeed, unlike standard mesh-based methods in SPH the fluid elementary volumes are followed in their Lagrangian motion. Therefore, the interface between two fluids will not be diffused and will remain sharply described, and no treatment is required to track them.

The present article revolves around the improvement of the multi-phase SPH formulation by Grenier et al. [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF][START_REF] Grenier | Viscous bubbly flows simulation with an interface sph model[END_REF], as it is suitable for simulating both interfacial and multiphase flows. Modifications have been made mainly in order to produce high quality pressure fields via the adaptation of a diffusive operator added in the continuity equation. The model is derived following two different ways leading to two sets of equations which present some differences, though ultimately they show little to no practical differences in the results. A series of validating tests were carried out over many well-known benchmarks, proving the versatile capabilities of the proposed multi-phase SPH model.

In 2003 Colagrossi and Landrini [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF] derived the first SPH multi-phase model able to treat density discontinuities with high density ratio in the fluid domain. In order to remove nonphysical high-frequencies from the pressure field, they used a periodic density filtering which, however, introduced numerical instabilities close to the free-surface for long time simulations. [START_REF] Español | Smoothed dissipative particle dynamics[END_REF] [START_REF] Español | Smoothed dissipative particle dynamics[END_REF] presented an SPH-DPD scheme where the particle volumes are evaluated through kernel summation; these volumes were called thermodynamic volumes, being noncoincident with the geometric ones. This approach was used by [START_REF] Hu | A multi-phase sph method for macroscopic and mesoscopic flows[END_REF] [START_REF] Hu | A multi-phase sph method for macroscopic and mesoscopic flows[END_REF] for another multi-phase SPH model. Conversely to [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF], in the latter the density field is directly evaluated through the ratio between the particle masses and their thermodynamic volumes without the use of a continuity equation integrated in time. However, this SPH model cannot be used if a free-surface is present since the thermodynamic volumes cannot be evaluated in its proximity because of the truncation of the kernel support. To circumvent this issue, [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF][START_REF] Grenier | Viscous bubbly flows simulation with an interface sph model[END_REF] proposed a multiphase model for interfacial and free-surface flows. They used a Volumetric Strain Rate equation for the time evolution of thermodynamic volumes, while the density field was evaluated by a Shepard correction.

The use of the Shepard kernel for the density field in [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF][START_REF] Grenier | Viscous bubbly flows simulation with an interface sph model[END_REF] is, however, not sufficient to completely eliminate the high frequency instabilities of the pressure field. In order to improve that model, in the present work we propose the use of the diffusive term introduced in the δ-SPH scheme by [START_REF] Antuono | Numerical diffusive terms in weakly-compressible sph schemes[END_REF]. Since in [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] the continuity equation is expressed for the thermodynamic volumes, the diffusive term has to be modeled in a different way with respect to [START_REF] Antuono | Numerical diffusive terms in weakly-compressible sph schemes[END_REF].

Another class of SPH models producing accurate pressure fields are the Riemann-SPH models, deriving from the early works of [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF] [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF], in an ALE (Arbitrary Lagrangian Eulerian) version, and [START_REF] Parshikov | Smoothed particle hydrodynamics using interparticle contact algorithms[END_REF] [START_REF] Parshikov | Smoothed particle hydrodynamics using interparticle contact algorithms[END_REF], in a Lagrangian one. In these Riemann-SPH models the standard SPH scheme is stabilized by solving Riemann problems at each pair interaction as in the Finite Volume Method for compressible flows. An extension to multi-phase flows of Vila's ALE Riemann-SPH [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF] was proposed by [START_REF] Leduc | A sph-ale 38 method to model multiphase flows with surface tension[END_REF] [START_REF] Leduc | A sph-ale 38 method to model multiphase flows with surface tension[END_REF] where mass fluxes were prevented at the interface to preserve its contact discontinuity nature, through the use of an approximate acoustic Riemann solver and an adequate choice of the interface velocity. This strategy was later extended by [START_REF] Guilcher | Simulation of liquid impacts with a two-phase parallel SPH model[END_REF] [START_REF] Guilcher | Simulation of liquid impacts with a two-phase parallel SPH model[END_REF] in a multi-phase Riemann-SPH model where renormalized SPH interpolators were used. In the present paper, we use a multi-phase Riemann-SPH model extended directly from the Lagrangian Riemann-SPH scheme of Parshikov and Medin [START_REF] Parshikov | Smoothed particle hydrodynamics using interparticle contact algorithms[END_REF] where no mass fluxes exist anywhere and an approximate acoustic Riemann solver is used everywhere. In this way the resulting multi-phase Riemann-SPH scheme is very close to the proposed (Lagrangian) multi-phase delta-SPH scheme but for the nature of the diffusion terms (approximate Riemann solver vs. δ-SPH terms). It is therefore a natural SPH model to compare with the proposed δ-SPH one. In particular, we show that the two SPH variants, δand Riemann-SPH, have different regions of stability when a multi-phase flow is considered. It is interesting that for some density and speed of sound ratios the δ-SPH model allows for bigger time steps, but in other conditions it is the Riemann-SPH model which is more convenient. Beside these differences in the stability of the two schemes, and CPU costs, we did not find substantial differences in terms of accuracy of the solution, both giving accurate pressure fields. However, this topic deserves further investigations which go beyond the present article. Note that the same kind of Riemann-SPH multi-phase model was used recently by Meng et al. [START_REF] Meng | A multiphase sph model based on roe's approximate riemann solver for hydraulic flows with complex interface[END_REF] but with a Roe's approximate Riemann solver instead of an acoustic one.

The present article is organized as follows:

• Section 2: the governing equations adopted in this work are briefly recalled.

• Section 3: the proposed numerical scheme is derived in two mathematically distinct frameworks. The first way is presented in section 3.1 where the SPH model is derived starting from the Volumetric Strain Rate equation. Then, in section 3.3 the SPH equations are derived using the direct kernel summation for the particle volumes. • Section 4: the Riemann-SPH model used for comparison is briefly presented.

• Section 5: time stepping and choice of the speeds of sound for the different phases are addressed. The comparison between δand Riemann-SPH is discussed. • Finally, in section 6 the proposed SPH model is tested against different testcases of increasing complexity. The advantages of the model are highlighted, especially when simulating interfacial flows with the presence of a free-surface.

Governing equations

In the present work we consider a domain Ω where different fluids or phases are present. The domain is confined by solid boundary Ω B and also by a free-surface Ω F . The solid boundary consists of an external frontier which contains Ω and also by solid bodies which can be inside the fluid domain and which can move across Ω F . Figure 1 show a sketch of the problem which can be handled by the SPH model proposed in this article.

The flow evolution is governed by the compressible Navier-Stokes Equation (NSE):

                 Dρ Dt = -ρ divu , Du Dt = - ∇p ρ + div ρ + g + F S De Dt = -p divu + : ½ ρ , Dr Dt = u, p = f (ρ) , (1)
where u is the flow velocity, r the position of the material points, p is the fluid pressure, ρ the fluid density, e the specific internal energy, the viscous stress tensor, ½ the rate of strain tensor, g is a generic specific body force and finally F S is the surface tension acting between different fluids/phases in Ω and on the free surface Ω F . Phase changes are not considered and we assume that all the fluids can be considered barotropic, hence, the pressure is linked only with the density field, ignoring the pressure dependency on the entropy. This condition is verified since in this work we are focused on flow with limited compressibility. Also thermal conductivity is neglected and only immiscible fluids are considered.

In this work the viscous and surface-tension models are not considered since the main interest is on free-surface flows and on water impacts with air-entrainment phenomena. Being the flow simulated inviscid, the numerical schemes need to be stabilized with diffusive terms both in the continuity and in the momentum equation.

Regarding the inclusion of viscous and surface-tension SPH models, similarly to what shown in [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] and in [START_REF] Grenier | Viscous bubbly flows simulation with an interface sph model[END_REF], they can be included in the present method without any restriction.

Multi-phase δ-SPH model

In the present section an extension of the δ-SPH model to multi-phase flow is derived in two different ways relying on the two main branches of the multi-phase SPH literature, namely [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] and Espanol et al. (2003) [START_REF] Español | Smoothed dissipative particle dynamics[END_REF]. Both of them are based on the use of equations for the particle volumes evolution so as to be able to treat problems involving discontinuity of the density field. It is then proven that the two derivations, even presenting differences in the equations, yield very similar results.

In both formulations presented in this section, similarly to the standard SPH model, the mass of the generic i-th particle, m i , is evaluated by the initial condition at time t = 0:

m i = ρ i0 V i0 , (2)
where ρ i0 is the initial density evaluated through the pressure field and the equation of state of the specific fluid whereas the particle volumes V i0 are known from the particles positions inside the fluid domain at the initial time. Usually at t = 0 the particles are regularly positioned in the domain. Therefore, at least at the initial instant, the particle volumes coincide with the geometrical ones, and can be initialized as such. One can use also the particle packing technique described in [START_REF] Colagrossi | Particle packing algorithm for SPH schemes[END_REF]; in that case the volume can be assigned by the obtained tessellation. Differently from the densities and volumes of the particles, the masses m i remain constant during the time evolution guaranteeing the conservation of the total mass.

Derivation from Grenier et al. 2009

In [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] the continuity equation is expressed through the Volumetric Strain Rate equation:

D δV Dt = δV div(u) , (3)
where δV is the infinitesimal volume associated to the material point which moves with the flow velocity u and D/Dt is the Lagrangian derivative.

Similarly to the SPH model presented in [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF] [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF] the use of the continuity equation allows to simulate flow in presence of a free surface as originally described by [START_REF] Monaghan | Simulating Free Surface Flows with SPH[END_REF] [START_REF] Monaghan | Simulating Free Surface Flows with SPH[END_REF]. This is not the case when adopting the direct evaluation of the volume or the density from the particle positions as in [START_REF] Hu | A multi-phase sph method for macroscopic and mesoscopic flows[END_REF].

As usually performed in the SPH framework, the divergence operator can be approximated as: div(u

) i = j (u j -u i ) • ∇W i j V j . (4)
where u i and u j are the flow velocity associated with the generic particles i and j, ∇W i j is the spatial derivative of kernel function, with respect the position of particle i, and V j the volume of the j-th particle. The kernel function W i j adopted in this work is the C2 Wendland function defined in [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree[END_REF]. The support of its domain is equal to 2h, where h is the noted as smoothing length. The W i j is a function of the relative distance between particle i-th an j-th, i.e. W i j = W(|r i -r j |; h).

From the above equations, the time evolution for the particle volumes can be derived as:

dV i dt = V i div(u) i . (5)
In order to integrate in time eq. ( 5) the initial particles volumes V i0 have to be initialized as described in the previous section. From these volumes also the radius of the particles ∆x i can be derived and the ratio h/∆x i is set equal to 2 (i.e. about fifty particles are present in the support of the kernel function W in a 2D framework).

The [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] model differs from other SPH models where a continuity equation is written in terms of density. Indeed, similarly to Espanol et al.

(2003) [START_REF] Español | Smoothed dissipative particle dynamics[END_REF], the particles volumes are evaluated before the density. It follows that the particles densities can be calculated directly through the ratio:

ρ i (t) = m i V i (t) , (6)
or through an interpolation formula. For example, in Grenier et al. (2009) [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF], after the time integration of eq. ( 5), a Shepard interpolation of the particles masses is used to evaluate the density field:

ρ i = j∈χ m j W S i j ; W S i j := W i j k∈χ W ik V k , ∀i ∈ χ , (7)
where χ indicates a generic fluid in the domain Ω (i.e. the above summations are computed only using particles belonging to the same phase). The function W S i j

is known in the literature as Shepard kernel and it satisfies by construction the normalization property:

j∈χ W S i j V j = 1; ∀i ∈ χ , (8)
for more details see, e.g., [START_REF] Fries | Classification and Overview of Meshfree Methods[END_REF]. Thanks to this relation the density discontinuities can be treated explicitly since, conversely to the kernel W i j , the summation of the Shepard kernel is not affected by the truncation of the kernel support. On the other hand, it is worth to note that the Shepard kernel W S requires the a priori knowledge of the particles volumes V k in eq. ( 7).

Once the density field is evaluated, the pressure is obtained through the use of a stiff equation of state:

p i = f χ (ρ i ) = ρ 0χ c 2 0χ γ χ ρ i ρ 0χ γ χ -1 + p b , ∀ i ∈ χ (9)
where γ χ is the polytropic coefficient of phase χ, p b is a background pressure and c 0χ is the speed of sound of the phase χ. In the present work p b is constant in time and in space within all the fluid domain. It is set at the beginning of the simulation in order to avoid negative pressure to appear in Ω during the time evolution and circumvent the inception of the tensile-instability (see, e.g., [START_REF] Le Touzé | A critical investigation of smoothed particle hydrodynamics applied to problems with free-surfaces[END_REF]). The latter is a classical problem for SPH models since it leads to spurious fragmentations of the fluid phases.

Concerning the choice of speed of sound, section 5 has been devoted to this topic. Indeed, beside physical aspects we show that also stability constrains are linked to this choice. It is worth to stress that state equation ( 9) is adopted for both gas and liquid phases, as the latter are treated as weakly-compressible media.

Following the work of Colagrossi et al. (2009) [START_REF] Colagrossi | Theoretical considerations on the free-surface role in the smoothed-particle-hydrodynamics model[END_REF], using the principle of virtual work the smoothed divergence operator defined by eq. ( 4) leads to the following smoothed pressure gradient:

∇p i = j (p i + p j )∇W i j V j . (10)
Since the Shepard kernel is used in eq. ( 7), in [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] also the smoothed velocity divergence and the smoothed pressure gradient operators were modified as:

                 div(u) i = j (u j -u i ) ∇W i j Γ i V j , ∇p i = j p i Γ i + p j Γ j ∇W i j V j Γ i := k W ik V k (11)
where Γ i is used as renormalized factor for taking into account possible kernel support truncation. However, conversely to the summation used in eq. ( 7) for the evaluation of the densities, the summation for the Γ i factor is extended also to particles of different fluids/phases. It follows that Γ i remains close to one and it plays a role only when the particles are close to the free surface. From a practical point of view we found that the effects of the Γ i terms are always negligible, their use does not improve the accuracy of the scheme and for this reason they will not be adopted in the rest of this work. A simple heuristic proof of this is given in section 6.2.

The scheme proposed by [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] is able to treat interfacial flows in the presence of a free-surface. However, as other weakly-compressible SPH schemes, it presents numerically high frequency oscillations on the pressure field.

In order to remove the genesis of these spurious acoustic waves the δ-SPH model by [START_REF] Antuono | Numerical diffusive terms in weakly-compressible sph schemes[END_REF] [1] can be adopted. The latter consists in adding a diffusive term in the continuity equation which is able to filter high-frequency numerical pressure noise. The pressure signals which are characterized by wave numbers with wave length smaller than the SPH kernel support are dissipated thanks to the use of this diffusive term.

It is useful here to remind the δ-SPH scheme in the single-phase context:

                           dρ i dt = -ρ i j (u j -u i ) • ∇W i j V j + δ h c 0 j D i j • ∇W i j V j ρ i du i dt = - j (p j + p i ) ∇W i j V j + ρ i f i + α h c 0 ρ 0 j π i j ∇W i j V j dr i dt = u i ; (12)
where c 0 is the artificial speed of sound, ρ 0 the density of the fluid at rest condition, and α and δ the two parameters for the artificial viscous terms. The δ parameter is always set equal to 0.1, while α ranges between 0.01 up to 0.1 depending by the problem at hand. For example for impulsive and violent impact problems α = O(0.1) is generally adopted in order to stabilize the scheme. Regarding the particle interaction terms D i j and π i j of the diffusive terms, they are defined as:

D i j := 2 (ρ i -ρ j ) - 1 2 (∇ L ρ i + ∇ L ρ j ) • r i j r i j |r i j | 2 , π i j := u i j • r i j |r i j | 2 , (13)
where r i j := (r i -r j ), u i j := (u i -u j ) and ∇ L ρ is the renormalized density gradient defined as:

∇ L ρ i := j (ρ j -ρ i ) L i ∇W i j V j , L i :=         - j r i j ⊗ ∇W i j V j         -1 (14)
(for more details see [START_REF] Antuono | Free-surface flows solved by means of SPH schemes with numerical diffusive terms[END_REF]).

Since the present method is based on a continuity equation for the particle volumes (see eq. ( 5)), the adaptation of the diffusion term proposed in [START_REF] Antuono | Numerical diffusive terms in weakly-compressible sph schemes[END_REF] is not trivial. An "intuitive" way to adapt the D i j term is to simply rewrite it in terms of particles volumes in the place of particle densities. However, this method proved to be unsuccessful in practice and we propose to rewrite the D i j term as follows:

D V i j := V i 2 1 - ρ j ρ i - 1 ρ i (∇ L ρ i + ∇ L ρ j ) • r i j r i j |r i j | 2 . ( 15
)
Even if it is added in the continuity equation written for the particles volumes and not for the density as in eq. ( 12), this new diffusive term is still depending on the density variations.

To sum up, the complete set of ODEs for the proposed model, written for a generic i th particle of generic phase χ, is:

                               dV i dt = V i j (u j -u i ) • ∇W i j V j + δ h c 0χ j∈χ D V i j • ∇W i j V j ρ i du i dt = - j (p j + p i ) ∇W i j V j + ρ i f i + α h c 0χ ρ 0χi j∈χ π i j ∇W i j V j dr i dt = u i , ρ i = j∈χ m j W S i j , p i = f χ (ρ i ) (16)
The parameters δ is always set to 0.1 independently of the simulation, and in the case of viscous flow the α-term is replaced by the real one, same as in singlephase δ-SPH [START_REF] Marrone | δ-SPH model for simulating violent impact flows[END_REF]. Note that the sums in the diffusive terms, as well as in the density equation, are made only over particles of the same phase as particle i. This is motivated by the fact that we do not want to alter the explicit treatment of the interface discontinuities by allowing a diffusion mechanism across the different phases. As discussed in [START_REF] Meringolo | A dynamic δ-SPH model: How to get rid of diffusive parameter tuning[END_REF] the δ-term alone is sufficient to stabilize the scheme, however the combination of the δ and α terms can improve the quality of the solution without introducing higher levels of numerical dissipation. This mechanism happens also in the Riemann-SPH scheme where numerical dissipation enters through the fluxes in both the continuity and the momentum equations (see Sec. 4). Nonetheless, in Riemann-SPH the control of the numerical diffusion is more complex since it is governed by the Riemann solver determining the fluxes and not simply by two parameters as in the present scheme.

Note that we chose to evaluate the density through eq. ( 7) mainly for coherence with the original work [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF]. Due to the presence of the diffusive term, in the above scheme ( 16) one could adopt either eq. ( 7) or ( 6) with negligible effects on the obtained solution.

Beside the use of the artificial diffusion terms, when treating multi-fluid flows the addition of a numerical surface tension is also required to stabilize the interfaces between the different phases, as explained in the following subsection.

Artificial surface-tension model

Simularly to what already underlined by Grenier et al. (2009) [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF], when using the models presented above to simulate multi-phase flows where surface tension effects are negligible, a non-physical inter-penetration of particles belonging to different phases may occur, leading to a numerical fragmentation of the fluid interfaces. In order to prevent this, a small repulsive force has been introduced in [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] modifying the pressure gradient as:

∇p i = j (p i + p j )∇W i j V j + χ j∈χ (|p i | + |p j |)∇W i j V j ( 17
)
where χ is a parameter ranging between 0.01 and 0.1. The second summation applies to particles that do not belong to the fluid of the i-th particle. This set of particle is noted by χ. The above choice naturally implies that on the free-surface the artificial pressure term of eq. ( 17) is null.

As shown in [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF], the last term of eq. ( 17) models a cohesion force, i.e. a surface tension, as also explained in [START_REF] Nugent | Liquid drops and surface tension with smoothed particle applied mechanics[END_REF].

In [START_REF] Szewc | Spurious interface fragmentation in multiphase sph[END_REF] [START_REF] Szewc | Spurious interface fragmentation in multiphase sph[END_REF] it is shown in an heuristic way how the spurious interface fragmentation in multiphase SPH can be controlled by the coefficient χ . However, it needs to be increased with a factor 1/h when increasing the spatial resolutions in order to maintain the effect of the second term of eq. ( 17).

Derivation from Español and Revenga (2003)

In the work of Español and Revenga [START_REF] Español | Smoothed dissipative particle dynamics[END_REF], the volume V i associated to a generic i th particle is written in terms of the SPH kernel function:

V i = 1 j W i j (18)
Therefore, V i is a direct function of the positions of the neighbouring particles of i. As highlighted in Español and Revenga [START_REF] Español | Smoothed dissipative particle dynamics[END_REF], V i does not coincide with the geometrical volume occupied by the particle, and in particular the sum of these volumes does not necessarily coincide with the total volume of the fluid domain, V T , i.e. i V i V T . This error decreases when increasing the ratio h/∆x if the kernel function is not truncated, which means that the formula ( 18) cannot be used if a free surface is present in the fluid domain. Since the summation of V i is not equal to V T , V i cannot be considered as the geometric volume of particle i, and for this reason in [START_REF] Español | Smoothed dissipative particle dynamics[END_REF] the quantity V i is called "thermodynamic volume", in the sense that it is directly linked to W i j which determines also the inter-particle forces in the momentum equation.

Once the particles volumes are evaluated by their positions, the density of particle i is calculated as:

ρ i (t) = m i V i (t) = m i j W i j (19)
In contrast to the classic SPH formula, ρ i = j m j W i j , used in the astrophysical context (see, e.g., [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF]), eq. ( 19) handles density discontinuities by preventing neighbouring particles masses from contributing to the evaluation of the particle density, especially in presence of other phases.

However, it falls short when a free-surface is present. Indeed, the assumption ( 18) yields wrong volumes near the free-surface due to the truncation of the kernel support, which in turn results in erroneous densities and pressures. Thus, in the present context of multi-phase with free-surface, we circumvent this issue by writing the time derivative of (18) as follows:

dV i dt = - j dW i j /dt ( j W i j ) 2 , dW i j dt = -(u j -u i ) • ∇W i j ( 20
)
where in the right equation the time dependency of the particle positions in the Kernel function has been considered. Eq. ( 20) yields the following expression:

dV i dt = V 2 i j (u j -u i ) • ∇W i j (21)
As discussed in section 3.1, the use of ( 21) implies an hypothesis on the initial conditions for the particles volumes. At the initial instant, the particle volumes coincide with the geometrical ones, and can be initialized as such, while this is not the case when Eq. ( 18) is directly used. Once again, since the volumes are evaluated by the time integration of eq. ( 21), the density can be also calculated through the Shepard interpolation in eq. ( 7) while the particle masses m i are set using the initial conditions and remain constant in time, ensuring mass conservation.

Following the work by Colagrossi et al. (2017) [START_REF] Colagrossi | Discussion of stokes' hypothesis through the smoothed particle hydrodynamics model[END_REF] the pressure gradient linked to the continuity eq. ( 21) is given by:

∇p i = j        p i V 2 i + p j V 2 j V i V j        ∇W i j V j . ( 22
)
Again this is derived using a variational principles on the energy equation of the flow.

To sum up, the complete set of ODEs for this second model, written for a generic i th particle of generic phase χ, reads as:

                               dV i dt = V 2 i j (u j -u i ) • ∇W i j + δ h c 0χ V i j∈χ D V i j • ∇W i j m i du i dt = - j (p j V 2 j + p i V 2 i ) ∇W i j + m i f i + α h c 0χ ρ 0χi j∈χ π i j 2 (V 2 i + V 2 j ) ∇W i j dr i dt = u i , ρ i = j∈χ m j W S i j , p i = f χ (ρ i ) (23)
where we maintained the same idea of adding the diffusive δ term to the volume equation in order to improve the quality of the pressure field. The system of eqs. ( 23) can be seen as an extension of the SPH model proposed by [START_REF] Hu | A multi-phase sph method for macroscopic and mesoscopic flows[END_REF] improving its pressure field evaluation and making it applicable also to free-surface flows.

Comparing eqs. [START_REF] Marrone | δ-SPH model for simulating violent impact flows[END_REF] with the corresponding ones ( 16) they differ in terms of the particle volumes averaging. However, since in this work the problems studied are characterized by compressible effects which are always limited, it follows that the variations of particle volumes within the kernel support remain quite small. As a consequence the differences in using the two formulations are in practice always negligible and they reduce when increasing the spatial resolution, (i.e. reducing the particle sizes ∆x i ). In Sec. 6.3.1 a case run by using either eqs. [START_REF] Hernquist | Treesph-a unification of sph with the hierarchical tree method[END_REF] or eqs. ( 23) is presented confirming their equivalence for incompressible problems and when the smoothing length is constant. Note that, for CPU efficiency, in the last test case of section 6.4 we used regularly-varying smoothing lengths far from the zone of interest of the problem, meaning that these two formulations there slightly differ.

More generally, the two formulations could present differences for problems with large compressibility of the fluids (problems not considered here). In such a condition h-variable models like the ones exposed in [START_REF] Hernquist | Treesph-a unification of sph with the hierarchical tree method[END_REF][START_REF] Richard | Variable smoothing lengths and energy conservation in smoothed particle hydrodynamics[END_REF] should be adopted. The latter models are more complicated and present different drawbacks. Note that to limit this effect a dedicated splitting/merging algorithm can be used to dynamically recover close-to-constant volumes even in case of very large compressibility as proposed in [START_REF] Pn Sun | Derivation and validation of δ-SPH model for simulating strongly-compressible multiphase flows[END_REF].

We underline that in the following part of the article the numerical results presented are obtained through the set of eqs. ( 16).

Multi-phase Riemann-SPH scheme used for comparison

As stated in introduction the multi-phase Riemann-SPH scheme we use here is directly derived from the Lagrangian Riemann-SPH model of Parshikov and Medin [START_REF] Parshikov | Smoothed particle hydrodynamics using interparticle contact algorithms[END_REF]. As for the δ-SPH scheme, in this model the particles masses m i remain constant throughout the simulations. In Riemann-SPH schemes a Riemann problem is solved at the interface between pair particles i and j located at rij = (r i + r j )/2, with an interface velocity ūij = (u i + u j )/2 and a normal direction n i j = ∇ i W i j /||∇ i W i j ||. The system of equations to solve then reads:

                           dV i dt = V i j 2 (u E -u i ) • ∇W i j V j ρ i du i dt = - j 2 p E ∇ i W i j V j + ρ i f i dr i dt = u i + δu i , ρ i = m i / V i (24)
where u E and p E are the Riemann problem solution at the interface between particles i and j. In the present scheme, it is obtained by the approximate acoustic Riemann solver of [START_REF] Murrone | A five equation reduced model for compressible two phase flow problems[END_REF] and using the piecewise linear MUSCL (Monotonic Upstream Scheme for Conservation Laws) method proposed by Van Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov's method[END_REF] to reconstruct the left and right states of the Riemann problem at the interface between i and j, with a limiter. This drastically lowers the numerical dissipation of the Riemann solver while increasing the scheme order. In this scheme the contact discontinuity at a physical interface is accounted for through the fact that the interface velocity is equal to the Riemann problem velocity at the interface ūij = ūE . Finally, in this scheme a particle shifting technique is used where the shifting law for δu i is the one proposed in [START_REF] Oger | Sph accuracy improvement through the combination of a quasi-lagrangian shifting transport velocity and consistent ale formalisms[END_REF].

5 Time stepping, choice of the sound speeds and numerical stability constraints

A fourth-order Runge-Kutta scheme is used to march in time the system of eqs. [START_REF] Hernquist | Treesph-a unification of sph with the hierarchical tree method[END_REF]. As described in [START_REF] Antuono | Numerical diffusive terms in weakly-compressible sph schemes[END_REF], in order to reduce CPU costs and improve the stability of the scheme, the diffusive δ-term is updated outside the sub-time steps of the Runge-Kutta scheme while it is "frozen" inside.

The time step of the simulations is set as:

∆t = min i (∆t c i ) , ∆t c i = K 2h c 0 χ ∀i ∈ χ (25)
where ∆t c i is the time step due to the acoustic constraints, K is the CFL factor, and 2h is the support of the C2 Wendland kernel (see, e.g., [START_REF] Antuono | Numerical diffusive terms in weakly-compressible sph schemes[END_REF]). The time step ∆t c i depends on the specific speed of sound of each phase χ. The K factor is set smaller than the maximum value K max which depends on the specific time scheme, the chosen kernel function and can be also problem dependent.

As stated in section 1, this work focuses on interfacial flows in the presence of a free-surface. More specifically, problems where a liquid X interacts with a gas pocket Y and a free-surface FS. The gas pocket Y in some conditions can be compressed by the liquid phase (for example in the water-entry problem of section 6.4). The compressibility of the gas is linked to its speed of sound c 0Y which cannot be chosen for computational convenience but has to be equal to the real one. In general, for cases where the compressibility is important for the gas pockets dynamics, the Euler dimensionless number should be preserved in the simulation. The latter is defined as:

Eu := P 0Y ρ X U 2 X , P 0Y := ρ 0Y c 2 0Y γ Y ( 26
)
where P 0Y is the ambient pressure of the gaseous phase and U X is the speed of the liquid acting against the gas pocket.

Conversely the liquid phase X is treated as a weakly-compressible medium, i.e. for computational convenience the speed of sound can be reduced with respect to the real one up to the limit:

c 0X ≥ 10 max         U max , ∆p max ρ 0X         (27)
where U max and ∆p max are the maximum expected velocity and pressure variation in the phase X . The condition [START_REF] Martin | Part iv. an experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF] guarantees the weakly-compressible regime, that is, fluid density variations for X shall remain smaller than 1% (see, e.g., [START_REF] Marrone | Prediction of energy losses in water impacts using incompressible and weakly compressible models[END_REF]). This allows the adoption of a lower value of the sound speed with respect to the the actual value which would lead to very small time steps and, therefore, high computational costs.

The above conditions where c 0Y is the real speed of sound of the gaseous phase and c 0X is a reduced speed of sound for the liquid, can lead to the counter-intuitive situation in which the air phase is modelled using a speed of sound larger than the liquid one. However, as explained in [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF][START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] and shown in section 6.4 the above choice remains a good approximation provided that the constraint ( 27) is respected.

In [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF] [10] a further constraint, this time on the ratio between the speeds of sound c 0X and c 0Y , was included and reads:

P 0Y = P 0X ⇒ c 0X c 0Y = γ X ρ 0Y γ Y ρ 0X (28)
The nature of this constraint was not justified by [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF]. In this work we show that it is linked to a stability constraint of the scheme and can be rewritten in a more complete fashion. In particular we found that varying the K factor we found further regions of stability which can be generally indicated as:

K < f c 0X c 0Y , η η := γ X ρ 0Y γ Y ρ 0X < 1 ( 29
)
where f is a monotonically non-increasing function of the variables (c 0X /c 0Y ) and η.

Numerical stability of the scheme through the Bagnold test

In order to derive the regions of stability given by the eq. ( 29) a 1-D Bagnold problem [START_REF] Bagnold | Interim report on wave-pressure research[END_REF] has been considered. In particular, the same configuration adopted in [START_REF] Guilcher | Simulation of liquid impacts with a two-phase parallel SPH model[END_REF] has been used for the test case: a liquid fluid patch is confined between two gas pockets under gravity forcing. In the initial condition the velocity is zero everywhere (see sketch in Fig. 2). For all the simulations the following parameters have been used: for the liquid γ X = 7, and for the gas γ Y = 1.4, δ = 0.1, α = 0.07, = 0.02 for both phases. This simple test case is used to study the conditions in which the scheme is stable. Several combinations of density ratios and speeds of sound have been tested.

In Fig. 3 the stability region is studied in the plane (c 0Y /c 0X , K). Three different ratios ρ 0X /ρ 0Y = 10, 100, 1000 have been considered which correspond to three η values 0.71, 0.22, 0.071 given by its definition in eq. ( 29). For this simple case the maximum K reached is K max = 1.4 much higher than the value 0.75 used in SPH simulations discussed in section 6. It is worth noting that in this 1D test case the advection of the particles is very small which is not the case for the other SPH simulations presented in the paper. The dots in Fig. 3 represent the maximum K value to obtain a stable SPH simulation for a given couple (c 0Y /c 0X , η). The unstable simulations, not shown here for the sake of clarity, appear as affected by numerical high frequency pressure oscillations inside the two gas cavities with increasing amplitude in time.

Using a set of forty-five simulations among the more than three hundred performed, the stability limits have been simplified using just two lines for each considered value of η, where the inclined lines present all the same steepness. A rule for the stability can then be established with the equation:

K < max K * -κ 1 log c 0X c 0Y 1 η , κ 2 η , η := γ X ρ 0Y γ Y ρ 0X < 1 ( 30
)
where K * = 1.13 is the obtained CFL number corresponding to the condition c 0X /c 0Y = η described in eq. ( 28). In [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF] the latter was needed to maintain the computation stable. The coefficients κ 1 = 0.82 and κ 2 = 1.20 do not depend on η while they could be dependent on the dissipation of the scheme, i.e. in our case on the parameters δ and α.

In most of the simulations performed in this article the ratio c 0X /c 0Y is always less than one, since we consider the X phase as a liquid and Y phase as a gas. As explained above, in such a condition, where c 0Y is higher than c 0X , it follows that the first condition of equation ( 30) dominates. The condition in eq. ( 28) should not be adopted since c 0X = η c 0Y could, in certain cases, not satisfy the weakly- compressible criterion [START_REF] Martin | Part iv. an experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF]. Instead, having fixed the adequate speed of the sound of the heavy phase through eq. ( 27), the K parameter should then be adjusted from eq. ( 30). As a consequence higher CPU costs are sometimes required.

In Fig. 4 the stability regions evaluated with the proposed δ-SPH scheme are compared with the ones given by the Riemann-SPH model described in section 4. The comparison shows that for the case c 0X /c 0Y = η the CFL number decreases to K = K + = 0.62 for the Riemann-SPH model, which is 1.8 lower than K * . On the other hand, the limits obtained for the Riemann-SPH model present a totally different shape than those for the δ-SPH model. Indeed, the CFL remains almost constant in the whole range of the ratio c 0X /c 0Y with just a small reduction of 0.1 when passing from c 0X /c 0Y < 1 to c 0X /c 0Y > 1. This behaviour is much more usual, the stability limit being only imposed by the CFL condition [START_REF] Marrone | Analysis of freesurface flows through energy considerations: Single-phase versus two-phase modeling[END_REF] for this model.

Anyway, from this plot comparing the two SPH models it is possible to conclude that, for a given density ratio ρ 0X /ρ 0Y , there are ratios of c 0X /c 0Y for which the Riemann-SPH model can be more CPU demanding than the δ-SPH one and viceversa. For example, considering the low density ratio ρ 0X /ρ 0Y ≤ 10 the CFL number for the δ-SPH model is always greater than the Riemann-SPH one. Conversely, this is not the case when increasing ρ 0X /ρ 0Y . For instance in the air-water case, in the condition c 0X /c 0Y > 0.15, the Riemann-SPH model remains stable with CFL numbers greater than for the δ-SPH one. The theoretical reasons of such kind of stability limits are still missing, especially for the proposed δ-SPH model exhibiting these peculiar stability slopes when c 0X /c 0Y < 1. This topic thus deserves further investigation, also in order to find possible solutions to go beyond the limits presented here.

Validation of the proposed model

The proposed model was validated through three classic test cases: a hydrostatic problem, an oscillating elliptic fluid patch and a dam-break. These benchmarks have the advantages of being easy to implement and of having available reference solutions, for which reasons they are widely used within the SPH community. In order to highlight the capabilities of the model, each simulation was run in singlephase and multi-phase configurations.

As final test case, a more challenging problem is addressed: the water impact of a corrugated panel with an entrapped air cavity. This problem presents also a complex free-surface dynamics and, therefore, is a good candidate for showing how the proposed SPH method is able to treat such kind of water impact flows.

For the simulations shown in the present section K has been set equal to 0.75, when not otherwise specified. Indeed with such a value the scheme has been found to be stable even for violent flow conditions.

Long-time evolution for the hydrostatic test-case

In this section the evolution of a tank partially filled with water at rest in hydrostatic pressure distribution is simulated. This simple test case was chosen in order to test the ability of the model in maintaining a stable free-surface in both single-and two-phase configurations. Indeed, new SPH formulations can induce non-physical particle re-settlement due to lack of consistency of the scheme in presence of boundaries and volume forces. For instance in [START_REF] Antuono | Numerical diffusive terms in weakly-compressible sph schemes[END_REF] it was shown that the use of diffusive terms in the continuity equation can induce instabilities at the free-surface if the latter is not properly taken into account. In the present study, for the singlephase configuration a two-dimensional tank is half-filled with water at rest. For the air-water simulation, the other half is filled with air, and the density ratio is 1000. The filling height is H and a downward gravity acceleration g is present. The computed single-phase and two-phase solutions up to t g/H = 62.64 are shown in Fig. 5. Similarly to [START_REF] Antuono | Numerical diffusive terms in weakly-compressible sph schemes[END_REF], in both solutions the free-surface is stable and the hydrostatic solution is well respected. However, in the two-phase configuration a non-physical separation occurs at the interface region between air and water. This problem was reported in [START_REF] Kruisbrink | An sph multi-fluid model based on quasi buoyancy for interface stabilization up to high density ratios and realistic wave speed ratios[END_REF], and a correction was proposed based on a buoyancy model via The kinetic energy of each single phase is made dimensionless using its initial potential energy.

the addition of a numerical term to the momentum equation. Here, we simply add a background pressure p b = 0.1 ρ w gH in order to resolve the issue, since we found that the interface separation is strongly linked to the appearance of negative pressure in the interface zone during the first time steps. The particles generally stay on the initial Cartesian grid, although a small perturbation is observed just underneath the free-surface and around the air-water interface.

Fig. 6 shows the normalized kinetic energy evolution up to t g/H = 300. At the beginning of the simulation a small amount of potential energy is converted in kinetic energy due to the relaxation of the particle Cartesian lattice (see, e.g., [START_REF] Colagrossi | Particle packing algorithm for SPH schemes[END_REF]).

Then the kinetic energy rapidly decreases towards zero for both single-and twophase simulations, showcasing the good stability and robustness of the model. The kinetic energy of the air-water simulation is two orders of magnitude smaller than its single-phase counterpart. This is due to the addition of the background pressure in the multi-phase case, which induces an higher level of numerical dissipation.

Oscillating drop: single/two-phase simulations

In the present section a 2D fluid patch evolving in a quadratic potential field Φ(x, y) = 1/2Ω 2 (x 2 + y 2 ) is considered. The fluid is inviscid and the radius of the disk is R. The velocity field is assumed to have the following form:

u = A(t)x v = -A(t)y (31)
where in the initial condition A(t = 0) is set equal to the parameter A 0 . The ratio Ω/A 0 is set equal to 1. Following the study by Monaghan and Rafiee [START_REF] Monaghan | A simple sph algorithm for multi-fluid flow with high density ratios[END_REF], under these conditions the free surface evolves periodically with an elliptical shape. In the present work the objective of this test-case is to validate the proposed SPH scheme monitoring the evolution of the mechanical energy, and by comparing our solutions of both single-and two-phase oscillating drops with their analytical counterparts.

Single-phase case

Fig. 7 presents the initial and deformed shapes of the fluid patch for three different time instants. The obtained free-surface configurations are compared with the analytical solution provided in [START_REF] Monaghan | A simple sph algorithm for multi-fluid flow with high density ratios[END_REF], resulting in a very good agreement. This is further illustrated in the left plot of Fig. 8, where the predicted evolution of the ellipse semi-axis a(t) is plotted against the analytical solution for 8 periods of oscillation.

A convergence study was made in order to evaluate the energy conservation properties of the scheme. Right plot of Fig. 8 plots the evolution of the normalized mechanical energy for three increasing discretization ratios R/∆x = 50, 100, 200. Similarly to Antuono et al. [START_REF] Antuono | Energy balance in the δ-SPH scheme[END_REF], for increasing numbers of particles, smaller dissipations of the mechanical energy occur due to the presence of the δ-term in the continuity equation, as its contribution goes to zero when the inter-particle distance goes to zero. For the finest resolution the dissipated energy after 10 oscillation periods is about 1% with respect to the initial one.

The present test-case is also used to study the influence of the integrals Γ in the Grenier's model which have been removed in the present scheme (see section 3.1). Left plot of Fig. 9 shows the time history of the pressure recorded at the center of the oscillating bubble evaluated with and without taking into account the integral Γ (see eq. ( 11)). The two solutions remain quite close, even if the scheme without the Γ factors seems to remove the acoustic pressure oscillations in a more effective way.

In left plot of Fig. 9 a pressure peak is observed for both solutions around t = 0.2T . This is due to the destruction of the particles tessellation which originate from the initial Cartesian distribution. This leads to a numerical pressure wave which is damped later in time. This phenomenon was reported in [START_REF] Colagrossi | Particle packing algorithm for SPH schemes[END_REF], where a particle packing algorithm was used to circumvent this issue. Also, in our computations we found that initializing the particles on a polar distribution helps relieving this issue, as shown in right plot of Fig. 9. Besides this, the initial particle configuration does not affect the solution in the later evolution.

Multi-phase case

Here we consider two initial concentric circular patches of fluid. The heavier fluid occupies the inner circular region of radius R/2, while the lighter fluid occupies the region between the outer circle of radius R and the inner circle. Both fluids are left to evolve in the quadratic potential Φ(x, y) = 1/2Ω 2 (x 2 + y 2 ), and we keep again the ratio Ω/A 0 = 1. The heavier fluid needs to be in the inner region, indeed if the phases are swapped the flow will be unstable under the action of the potential Φ and a Rayleigh-Taylor instability will take place.

This second benchmark is particularly interesting for the validation of the present model since the multi-fluid domain is confined by a free surface. In the initial condition the pressure is correctly initialized, and unlike [START_REF] Monaghan | A simple sph algorithm for multi-fluid flow with high density ratios[END_REF], no damping technique is needed to initialize the simulation. Fig. 10 illustrates the initial configuration of the problem, with a density ratio ρ inner /ρ outer = 1000.

Fig. 11 shows the evolution of the outer ellipse semi-axis a(t) for two different density ratios, ρ inner /ρ outer = 10 and 1000. As predicted by the analytical solution, the time evolution of the outer axis, a(t), is not affected by the density ratio ρ inner /ρ outer . This is not the case for the pressure at the center of the inner fluid which is plotted in Fig. 12. Once again, a small perturbation at the end of the first quarterperiod occurs due to the initial Cartesian distribution of the particles. Apart from that, the pressure predicted by current model follows closely the incompressible solution for both density ratios.

Fig. 13 shows the pressure in the inner and outer phases along the horizontal line y = 0 at the initial time and after 2 periods for density ratio ρ inner /ρ outer = 1000. For this case, c 0X = 15A 0 R and c 0Y = 212A 0 R. The pressure after 2 periods matches quite accurately the initial distribution as expected. At the interface, some particles, especially those from the lighter phase, deviate from the expected value. Beside these drawbacks in the local solution, the proposed model reproduces the analytical incompressible kinematics with a satisfactory accuracy and yields smooth pressure fields in both lighter and outer phases.

Dam Break flow: single-phase and air-water configurations

The third test case is a dam-break flow impacting a vertical wall. It is one of the most popular simulations in the SPH community, due to the existence of experimental results [START_REF] Martin | Part iv. an experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF][START_REF] Stansby | The initial stages of dam-break flow[END_REF][START_REF] Buckner | Green water on ship-type offshore structures[END_REF] for different setups (tank geometry, water length and height, dry/wet bed etc.). Moreover, it is a challenging test case because of the free-surface fragmentation, vorticity generation due to multiple breaking processes, occurrence of water impacts.

In the present study, the dam-break configuration follows the recent experiments made by Lobovsky et al. [START_REF] Lobovský | Experimental investigation of dynamic pressure loads during dam break[END_REF], where different impact pressure measurements on the vertical wall are provided. The probe locations in our simulations were chosen in accordance with the ones from their experiments.

Single-phase simulation

The initial single-phase problem is illustrated in Fig. 14. The discretization ratio is H/∆x = 200. Snapshots of the dam-break flow evolution in the single-phase configuration are shown in Fig. 15. The model clearly handles well the presence of the free-surface. Also, thanks to the addition of the diffusive terms in the present model, no spurious oscillations of pressure are observed. This is the main improvement upon [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF] as illustrated in Fig. 16, which highlights the differences in the resulting pressure fields obtained with the proposed model and the original Grenier et al. ( 2009) model [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF]. Finally, in Fig. 17 a snapshot of the pressure field at time t(g/H) 1/2 = 6.8 is depicted using the present model, i.e. eqs. ( 16) derived in Sec.

Two-phase simulation

Another simulation of the same dam-break flow was run, this time in a multi-phase configuration. The adopted density ratio is ρ water /ρ air = 1000. The sound speeds are chosen as c 0X = 10 √ gH and c 0Y = 343 √ gH. The background pressure value here is P b = ρ w gH. Fig. 18 shows the dam-break evolution and the pressure field in the water phase. The latter is in very good agreement with its single-phase counterpart, as the pressure levels inside the water are almost identical in both cases up to the cavity closure, after which the air cushioning effects become relevant (see third plot of Fig. 18) and the two flow evolutions do not behave in the same manner any more. This visual agreement is confirmed by monitoring the pressure signal at probe P 1 where the maximum pressure on the wall is expected. As observed in Fig. 19 both models are in very good agreement up to around t(g/H) 1/2 = 6, which corresponds to the first cavity closure. Then the two-phase model predicts a pressure oscillation due to air-cushioning effect of the entrapped air bubble, similarly to what observed in [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF]. A different behaviour is observed in the single-phase simulation, in which this cavity is void. In this case the cavity collapses in a flat impact around t(g/H) 1/2 = 8.0, transferring some of the mechanical energy into internal energy, in the form of traveling acoustic waves (see also [START_REF] Marrone | Analysis of freesurface flows through energy considerations: Single-phase versus two-phase modeling[END_REF] for a detailed discussion about this phenomenon). In Fig. 19 the pressure recordings from [START_REF] Lobovský | Experimental investigation of dynamic pressure loads during dam break[END_REF] are also reported. A good agreement between numerical output and experimental data is achieved for the first pressure peak. The pressure peak evaluated by the SPH models are practically identically and in agreement with the maximum pressure value predicted by the potential flow theory (see, e.g., [START_REF] Marrone | δ-SPH model for simulating violent impact flows[END_REF] for details), whilst the maximum pressure recorded in the experiments appears higher of 32%. Further discrepancies between the numerical predictions and the experimental signal are also observed further in time. These differences are to be ascribed to relevant 3D effects observed in [START_REF] Lobovský | Experimental investigation of dynamic pressure loads during dam break[END_REF] and deserve supplementary investigations.

Fluid impact of a corrugated panel with trapped gas cavity

In order to test the ability of the proposed model to correctly compute the pressure in challenging problems such as the water impact, the problem of the water entry of a corrugated panel on a liquid free-surface with a trapped gas cavity is considered. Khabakhpasheva et al. provided a semi-analytical solution for this problem in [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF] derived from the Wagner theory. In that work they studied the initial stage of an incompressible liquid impact onto both rigid and elastic corrugated panels, accounting for a compressible gas pocket trapped between the corrugations (see sketch in figure 20). The corrugation shape function is described by the following function:

f (x) = h 1 -cos 4 πx 2 L (32)
where h and L are the corrugation's height and half-width respectively. This panel shape is noted in literature as Mark III panel which is a type of containment tank of Liquified Natural Gas (LNG) carriers (for more details see, e.g., [START_REF] Marrone | Challenges on the numerical prediction of slamming loads on lng tank insulation panels[END_REF]).

The distance between the two corrugations is 2L (see Fig. 20) whereas the total panel length is 2c. Similarly to [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF] the following dimensions are adopted in all the simulations: h = 3.6 cm, 2 L = 5.6 cm and 2L = 34 cm. At the initial time t = 0, the corrugations touch the free-surface at two points x = ±L and the fluids are at rest. The gas cavity is bounded by the panel and the free-surface and the panel penetrates the liquid surface with a purely vertical velocity U. In [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF], the gas in the cavity is modelled adopting the polytropic state equation ( 9) discussed in section 3.1. In that work the pressure inside the cavity is assumed as a function of the gas volume, which means that the pressure is considered homogeneous in the gas phase. The solution provided in [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF] is general and valid for any density ratio between the liquid and the gas phase. It is worth to note that in [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF] outside the gas cavity a free-surface is considered and the SPH simulations have been set using the same condition.

Since the focus is on high density ratios configurations, in the following test cases the water and air stand for the liquid and gas phases respectively, with a density ratio of 1000 and the polytropic coefficients are set to γ air = 1.4 and γ water = 7. In order to correctly take into account the compressibility effects in the air cavity the same value used in [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF] for the gas speed of sound is adopted, that is c 0air = 343 m/s. As for the water sound speed, formula ( 27) is used to guarantee the weaklycompressible hypothesis in the impact region, this being the common choice in SPH to approximate incompressible fluids. Based on the expected maximum pressure variation ∆P max = 1 atm in the cavity, the sound speed in water is obtained as c 0water = 10 ∆P max /ρ water = 100 m/s.

Note that, following the stability conditions described in section ( 5), the CFL is set equal to K = 0.2, with α = 0.07. In order to avoid the mixing of air and water, the artificial surface tension correction ( 17) is used, and the parameter χ is set to 0.02.

In the following, the dimensionless variables are denoted by the superscript *. The length L is taken as reference length scale of the problem, unless specified otherwise; U and h/U are chosen as, respectively, the velocity and time scales; the ratio P 0air = c 2 0air ρ 0air /γ air is taken as the pressure reference scale in the gas cavity.

Two discretizations were tested, L/∆x = 160, 640. These values are based on the work of Marrone et al. [START_REF] Marrone | Challenges on the numerical prediction of slamming loads on lng tank insulation panels[END_REF], where the same type of corrugated panel was tested with a dead-rise angle of 4 • , and in which these resolutions yielded satisfactory, convergent results. Note that the panel impact generates several acoustic waves which would be reflected off of the domain walls if the latter is not large enough, resulting in noisy pressure fields even in the air phase. Therefore, the domain is made sufficiently large and its dimensions were chosen simply by computing the distance which acoustic waves would travel during the impact at the adopted water sound speed. The variable-h scheme discussed in [START_REF] Oger | Two-dimensional sph simulations of wedge water entries[END_REF][START_REF] Oger | On distributed memory mpi-based parallelization of sph codes in massive hpc context[END_REF] was used since the focus is on the gas cavity and its vicinity. The particle size varies between the air cavity to the domain walls from fine to coarse respectively, with a magnification factor of 100 between the smallest and biggest resolutions. Since the problem is 2D symmetrical with respect to the line x * = 0, a symmetry boundary condition was used. The adopted discretization for the finest case is shown in Fig. 21.

Results and analysis: Kinematics

The focus here is on the case where the impact velocity and the non-dimensional panel length are set to U = 2 m/s and c * = 2. The obtained flow kinematics is studied by comparing the shapes of the inner and outer free-surface with the ones provided in [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF], and by tracking the time evolution of the inner and outer contact points, namely c * 1 and c * 2 .

Fig. 22 shows the evolution of the air/water flow through four snapshots at different instants. After the panel hits the water surface, its corrugations penetrate it, and the air cavity starts to compress. Shortly after, the extremities of the cavity start to creep under the corrugations and towards the outside of the panel. This air movement is accompanied with the water jet in the outer region freely moving towards the outside of the panel

The contact point is defined here as the abscissa of the intersection between the water front and the panel. In our simulations, this was implemented as follows: • First, a control zone is defined by setting a criterion of distance to the corrugation wall. In particular, only particles with minimal distance to the panel smaller than ∆x are considered in the control zone, in order to guarantee the presence of at least one particle (see Fig. 23).

• Then, at each iteration, and for each side of the corrugation, we select the water particle lying in the control zone, with smallest x-coordinate within the ones to the left of the corrugation, and respectively with highest x-coordinate within the ones to the right of the corrugation. The contact point abscissa c * 1 is then considered to be the x-coordinate of the left water particle selected, and respectively the x-coordinate of the right water particle selected for c * 2 .

In Fig. 24 the evolution of the SPH contact points c * 1 and c * 2 are plotted against the ones from [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF]. The general trend of the predicted outer contact point follows closely its semi-analytical counterpart. Also, the slope of both curves changes at around t * = 0.5, when the water front reaches the intersection between the corrugation and the horizontal part of the panel, where a change of convexity occurs. Finally, both curves stagnate at c * 2 = 2 since it corresponds to the extremity of the panel.

The time evolution of the contact points is subdivided in four stages:

Stage 1) for t * ≤ 0.45 is the compression stage of the entrapped cavity. At t * = 0.41 an increase of the inner contact point abscissa is registered, which signifies the creeping of the cavity towards the outside of the corrugation. Stage 2) it is characterized by the formation of a high speed water jet which sucks the cavity towards the outside of the corrugation. Stage 3) after the jet reaches the end of the panel, the suction effect becomes dialled down. As the panel continues to penetrate the water, the cavity resumes its compression in this stage. Stage 4) around t * = 0.85, the gas cavity begins its decompression, as better explained in the next subsection 6.4.2 where the dynamic part of the solution is investigated.

Note that the curves do not begin exactly at the same starting x-coordinate. This is a consequence of the adopted definition for the contact point. Indeed, as shown in Fig. 23, at the beginning of the simulation the particles identified as contact points do not lie exactly on the axis x * = 1, this difference decreasing as the resolution increases.

There are some differences between the SPH and semi-analytical slopes which are limited and reflect the different natures of the solutions being compared. As for the inner contact point c * 1 , it is strictly linked to the volume variations of the air cavity. In Khabakhpasheva et al. [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF] the air cavity was prevented from reaching the head of the corrugation by imposing the condition ċ1 ≤ 0; with this constraint the horizontal velocity of the inner contact point always points towards the center of the cavity, which means that the air will never flow under the corrugation. In the present simulations however, no condition on the inner contact point's velocity is enforced, so that the cavity is left to freely evolve in time.

The different behaviour can be observed in Fig. 24 for c * 1 at t * 0.45. As expected, the condition ċ1 ≤ 0 strongly influences the kinematics of the air cavity. Indeed, the constrained contact point does not deviate much from its original position whereas, without constraints, the air cavity starts to move towards the outer region under the corrugation, as seen in Fig. 22. Moreover, up to t * = 0.25 both models predict the same evolution of the air cavity.

Results and analysis: Dynamics

In this subsection, the pressure evolution inside the gas cavity is analyzed and compared to its semi-analytical counterpart. First, the simulation results using coarse (L/dx = 160) and fine (L/dx = 640) discretization levels are plotted in Fig. 25. The SPH solution can be viewed as a succession of 4 stages as done for the discussion of the solution of the kinematic (see subsection 6.4.1):

Stage 1) A compression of the gas up until t * = 0.45 during which the pressure grows linearly akin to an elastic reaction, due to the homogeneous water load around it as shown in Fig. 26. Stage 2) during the high speed water jet motion, the cavity is sucked towards the outside of the corrugation inducing a relaxation in the gas cavity during which the pressure keeps decreasing. Stage 3) After the jet reaches the end of the panel, the suction effect becomes dialled down and is counterbalanced by the water pressure from under the gas cavity. Therefore, as the panel continues to penetrate the water, the cavity resumes its compression in this stage. Stage 4) Around t * = 0.85, the pressure reaches its second peak then the gas cavity begins its decompression as the loads balance out during this final stage.

Interestingly, both SPH and semi-analytical solutions predict these four phases. The major difference is in Phase 2 where the SPH pressure decreases but the semianalytical one stagnates. In [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF], the pressure is directly linked to volume through the following equation:

P(t) = [V(t)] -γ -1 (33)
and the volume itself is strongly dependent of the inner contact point evolution c 1 (t). The speed ċ1 is computed at each time step and is set to zero if the computed value is positive. The reasoning behind this condition is extensively covered in [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF]. Thus the contact point is prohibited from ever moving towards the outside of the corrugation, which means that the cavity volume either decreases or stagnates. Consequently, following (33), the semi-analytical pressure can only stagnate or increase. Nevertheless, both models predict the elastic compression in Phase 1, a second pressure rise in Phase 3, and a final decompression in Phase 4, although at different times/pressure levels due to the differences in constraints between the SPH and semi-analytical frameworks. The simulation has been performed also using the Riemann-SPH model (see Sec. 4) and in Fig. 27 the comparison of the time evolution of the pressure inside the cavity for the δand Riemann-SPH models is shown, confirming that the two SPH models provide results with comparable accuracy. As a final study, the effect of the panel length on the cavity pressure is investigated. The numerical and geometrical configurations are kept unchanged, however this time the panel length c * = 3. The simulation was run with the fine discretization ratio L/∆x = 640. Fig. 28 compares the evolutions of the cavity pressures for both panel lengths c * = 2 and c * = 3. Fig. 28 regroups the semi-analytical and computed pressure solutions for both panel lengths c * = 2 and 3, from which some general conclusions can be made. First, in accordance with [START_REF] Khabakhpasheva | Fluid impact onto a corrugated panel with trapped gas cavity[END_REF], the shorter the panel, the lower the pressure inside the cavity. Also, right away it can be seen that both pressure curves follow the same general trend, in the sense that they both predict the aforementioned 4 phases. The curves are superposed during Phases 1 and 2 which is expected since the water does not reach c * yet. The solution with the larger panel agrees the conclusions made for the shorter panel, meaning that the SPH pressure follows closely the semi-analytical one in terms of compression + relaxation in Phase 3 and 4, albeit at different pressure levels since the air cavity is not constrained within the SPH model, while it is set to move only towards the inside within the semi-analytical model.

Conclusion

In this paper, we proposed a multi-fluid weakly-compressible SPH model in which the δ-SPH diffusive terms were extended to multi-phase flow, and adapted to the volumes equation of the Grenier et al. model [START_REF] Grenier | An hamiltonian interface sph formulation for multi-fluid and free surface flows[END_REF][START_REF] Grenier | Viscous bubbly flows simulation with an interface sph model[END_REF] in order to eliminate the spurious oscillations of pressure. Two methods of deriving the model were presented. It was found that despite the different theoretical basis of the two derived variants no substantial differences were observed on the results.

The time-stepping and stability limits in terms of CFL number of the proposed model were studied in detail heuristically on the 1D Bagnold test case. It permitted to draw a stability map depending on the sound speed and density ratios of the two phases considered in the problem. The obtained complex stability region was compared to the one of a Riemann-SPH multi-phase model whose formulation was also given in this article. As the proposed δ-SPH multi-phase model, this multiphase Riemann-SPH model can handle both the presence of fluid interfaces and a free-surface, and provides similar accuracy in terms of pressure field quality. The time stepping and the choice of the speeds of sound for the different phases were discussed, highlighting how this choice is driven by physical consideration and numerical constraints linked to the stability of the scheme. Concerning the latter, it was shown that the two SPH variants, δand Riemann-SPH, have different regions of stability when a multi-phase flow is considered. For some density and speed of sound ratios the δ-SPH model allows for larger time steps, but in other conditions the contrary holds.

Different test-cases using high-density ratios were then run for validating the proposed δ-SPH multi-phase model. First, it was checked that this model correctly maintains the hydrostatic solution both in single-and two-phase configurations.

The second test case simulated single-and two-phase oscillating drops, with a focus on the assessment of the quality of the pressure fields in both light and heavy phases. The kinematics and pressure fields were found to be in very good agreement with their analytical counterparts.

The third test case concerned a water impact generated by a dam-break flow.

As for the previous test-cases the simulations were conducted with single-phase and then with air-water. Pressure signals on the impact wall were compared to the experimental data of Lobovsky et al. [START_REF] Lobovský | Experimental investigation of dynamic pressure loads during dam break[END_REF]. A close match was found on both configurations (single and two-phases) in the early stage of the water impact. Then, air bubble oscillation effects were observed in the pressure signals in the air-water simulation. During this second impact stage discrepancies with the experimental data were observed and linked to the 2D framework adopted in this work.

As a final benchmark a more challenging test case was studied, in which the presence of a free-surface is combined with the presence of an entrapped air cavity generating a complex water impact flow. It involved the water entry of a corrugated MarkIII panel. The corrugations are small rigid structures under the flat panel, and their shape is given by an analytical function. The SPH results were compared with a semi-analytical approach based on the Wagner theory. The water flow and pressure distribution in the air cavity were carefully analysed, challenging the validity of the assumptions made in the semi-analytical Wagner solution, and studying the influence of the panel length. It was especially found that the trapped air pocket could flow out of the corrugation, which was prevented in the semi-analytical solution. Despite this, good agreement was obtained with this solution. It was also confirmed that both the proposed model and the Riemann-SPH model were predicting very similar flow dynamics and pressure. Therefore, these promising results make this simulation a solid benchmarking case for the validation of other schemes.

Future works could revolve around the extension and deeper understanding of the stability regions of the scheme, on a larger set of canonical tests and, if possible, in comparison to an analytic prediction of the stability criteria. Further, the introduction of a particle shifting technique (see, e.g., [START_REF] Pn Sun | A consistent approach to particle shifting in the δ-plus-sph model[END_REF]) in the proposed scheme could lead to further improvement of the numerical solutions.

Fig. 1 .

 1 Fig. 1. Sketch of the fluid domain. The domain Ω is confined by a solid surface, different phases a, b are present with an empty region, confined by the free surface ∂Ω F . Moving solid bodies can be also present inside the domain and Ω can be also confined by a solid surfaces. All the above solid frontiers are are indicated with ∂Ω B .

Fig. 2 .

 2 Fig. 2. Sketch of the 1D Bagnold test. A patch of fluid X is confined by two gas pockets Y. At t = 0 the fluid velocity is zero and the fluid is subjected to the gravity force.

Fig. 3 .

 3 Fig.3. Region of stability for a two-phase Bagnold problem. Each bullet represents the maximum CFL number needed for a given density ratio ρ 0X /ρ 0Y and speed of sound ratio c 0X /c 0Y , X and Y being the heavier and lighter phase, respectively. The solid lines represent the regression lines delimiting the stability region for each η value, η being the parameter defined in formula[START_REF] Meringolo | A dynamic δ-SPH model: How to get rid of diffusive parameter tuning[END_REF]. The yellow bullets in the graph represent the points where c 0X /c 0Y = η and K = K * .

Fig. 4 .

 4 Fig. 4. Two-phase Bagnold problem: comparison of the stability regions of the δ-SPH and of the Riemann-SPH. Each bullets represent the maximum CFL number needed for a given density ratio ρ 0X /ρ 0Y and speed of sound ratio c 0X /c 0Y , X and Y being the heavier and lighter phase, respectively. The yellow bullets in the graph represent the points where c 0X /c 0Y = η and K = K * while the green bullets are the points where c 0X /c 0Y = η and K = K + .

Fig. 5 .

 5 Fig. 5. Hydrostatic solution: particle configuration and pressure distribution in single-phase (left) and two-phase (right) at t g/H = 62.64.

Fig. 6 .

 6 Fig.6. Evolution of the kinetic energy for the single-phase and air-water hydrostatic cases. The kinetic energy of each single phase is made dimensionless using its initial potential energy.

Fig. 7 .

 7 Fig. 7. Evolution of the single-phase oscillating drop at different times. The dashed lines represent the analytical solution at the maximum, the intermediate and the minimum elongations of the vertical abscissa.

Fig. 8 .

 8 Fig. 8. Left: Comparison between the predicted evolution of the semi-axis a(t) and the analytical solution. Right: Time history of the normalized mechanical energy variation for different discretizations.

Fig. 9 .

 9 Fig. 9. Time history of the pressure at the bubble center r = 0. Left: computed with and without taking into account the integral Γ. Right: computed with initial polar/Cartesian particle distributions.

Fig. 10 .

 10 Fig. 10. Evolution of the pressure field for the concentric circular bubble problem at different times. The dashed lines represent the analytical solution at the maximum, the intermediate and the minimum elongations of the vertical abscissa.

Fig. 11 .

 11 Fig. 11. Outer ellipse semi-axis time evolution for a density ratio of 10 (left) and 1000 (right). The δ-SPH solution is compared with the analytical solution.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Evolution of the pressure at the center of the inner bubble for density ratios of 10 and 1000.

Fig. 14 .

 14 Fig. 14. Initial configuration of the single-phase dam-break problem.

Fig. 15 .

 15 Fig. 15. Snapshots of the single-phase dam-break flow at different times. Colours are representative of the dimensionless pressure P/ρ w gH from 0 (blue) to 1 (red).

Fig. 16 .

 16 Fig. 16. Contours of the pressure for the dam-break flow at t(g/H) 1/2 = 9, using the original Grenier et al. model (left) and the proposed model (right).

Fig. 17 .

 17 Fig. 17. Contours of the pressure for the dam-break flow at t(g/H) 1/2 = 6.8, using the proposed model (eqs. (16), top) and its variant (eqs. (23), bottom).

  3.1 from the Grenier et al. (2009) model, and using its variant, i.e. eqs. (23) derived in Sec. 3.3 starting from Español and Revenga (2003). The two solutions are very close to each other confirming the equivalence of the two formulations, as underlined at the end of section 3.3.

Fig. 18 .

 18 Fig. 18. Snapshots of the multi-phase dam-break flow evolution at different times. Colours are representative of the dimensionless pressure (P-P b )/ρgH from 0 (blue) to 1 (red) while the solid gray line is the air-water interface.

Fig. 19 .

 19 Fig. 19. Time history of the pressure at the probe P 1 . Comparison between the current single/two-phase results and a typical dambreak event from Lobovsky et al.

Fig. 20 .

 20 Fig. 20. Sketch of the problem Fluid impact of a corrugated panel with trapped gas cavity corrugation panel.

Fig. 21 .Fig. 22 .

 2122 Fig. 21. Left: Multi-resolution discretization of the domain. Right: Zoom on the gas cavity vicinity.

Fig. 23 .

 23 Fig. 23. Definition of the contact point after the first iteration.

Fig. 24 .

 24 Fig. 24. Evolution of the SPH/semi-analytical internal and external contact points c * 1 and c * 2

Fig. 25 .

 25 Fig. 25. Comparison with the semi-analytical solution of the pressure inside the cavity, for a panel length c * = 2. The pressure curve is divided into 4 stages.

Fig. 26 .

 26 Fig. 26. Pressure field during the impact stages for four different time instants t * = 0.014, 0.42, 0.55, 0.90. The panel length is c * = 2.

Fig. 27 .

 27 Fig. 27. Comparison between δ-SPH and Riemann-SPH against the semi-analytical solution of the pressure inside the cavity, for a panel length c * = 2. The pressure curve is divided into 4 stages.

Fig. 28 .

 28 Fig. 28. Semi-analytical and computed pressure profiles in the cavity, for panel lengths of c * = 2 and 3.

Acknowledgments

The research leading to these results has received funding from the European Union's Horizon 2020 Research and Innovation Programme (Grant No. 724139). The SPH simulations performed under the present research have been obtained using the SPH-Flow solver, software developed within a collaborative consortium composed of Ecole Centrale de Nantes, NextFlow Software company and CNR-INM. This work was also partially supported by the SLOWD project (which received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 815044). The proposed SPH multi-phase model will be used in the context of violent sloshing wing dynamics within the framework of this project.


• 2 Barrier CMP processes performances were investigated: POR process B was compared to process B1

• CMP performances were compared by quantifying bonding yield on 3 different test masks with different design rules:

• Heterogeneous Mask 1 with pitch and density variations • Mask II, with fixed density and pitch variation • Mask III, with density variation and fixed pitch

• While we will not give precise process conditions, we will talk about a method to quantify bonding yield

What we will talk about

Qualitative Estimation of Bonding Performance

| 9

Qualitative Estimation of Bonding Performance  HR-SAM Images obtained for C1 Cu CMP Process show that white pixels ratio is significantly decreased and, also, their radial distribution on the wafer is significantly improved  Hybrid bonding quality is highly ameliorated by changing from C Cu CMP process to C1

Cu CMP process and CMP process signature is evidenced  From mask design perspective, at die scale, white pixels ratio seems reduced for C1 Cu CMP process but is difficult to conclude at this stage Qualitative Estimation of CMP Process Impact on Bonding Quality Barrier CMP  Concerning barrier process conditions, comparing B to B1 is difficult, once again: HR-SAM images show comparable white pixels ratio and distribution so choosing the most performant Barrier CMP process in terms for bonding quality is not obvious

  . BALAN, Y.EXBRAYAT, D.SCEVOLA • While BEOL CMP is a mature process, CMP for 3D hybrid bonding is much more challenging: Cu pad size is bigger while topography constraints are much more aggressive in order to avoid critical bonding voids • Quantifying CMP process performance is done by estimating final topography at different scales, with different metrology techniques • Correlation between topography amplitude and Scanning Acoustic Microscopy (SAM) images is done in order to improve bonding yield Intro | 6 CAMP2021 | V. BALAN, Y.EXBRAYAT, D.SCEVOLA • While BEOL CMP is a mature process, CMP for 3D hybrid bonding is much more challenging: Cu pad size is bigger while topography constraints are much more aggressive in order to avoid critical bonding voids • Quantifying CMP process performance is done by estimating final topography at different scales, with different metrology techniques • Correlation between topography amplitude and Scanning Acoustic Microscopy (SAM) images is done in order to improve bonding yield • In this work, we developed a image analysis method in order to have a direct quantification of the bonding yield from SAM images and propose best CMP process Intro | 7 CAMP2021 | V. BALAN, Y.EXBRAYAT, D.SCEVOLA • Improve a CMP process for Hybrid Bonding  Quantify Bonding Performance • CMP Process performances are investigated through High-Resolution Scanning Acoustic Microscopy HR-SAM • 3 Cu CMP processes performances were studied : standard, POR process C was compared to new consumables processes C1 & C2

•

  HR-SAM image of bonded wafers after POR CMP • Mask I: pitch & density variation | 10 Qualitative Estimation of Bonding Performance • Process signature impact both at • Wafer level  radial non-uniformity can be seen on the extracted band from the wafer image, indicating CMP process non-uniformity: bonding quality non-uniformity is quantified by white, unbonded zones, giving information on CMP process performance • Die-Level  we can clearly see that some white squares are repetitive  corresponding Cu arrays with a certain density and/or pitch are not bonded, giving information on mask design performance NU% | 11 Qualitative Estimation of CMP Process Impact on Bonding Quality Cu CMP (1)



  Changing for another Cu CMP process conditions, Cu C2, the HR-SAM Image show comparable white pixels ratio and distribution: the hybrid bonding quality seems equivalent  choosing the most performant CMP process is difficult in this case



  *  Surface Analysis Script  Radius Band Extraction  Image Thresholding  White Pixels Extraction  White / Black Pixels Ratio Unbonded %Area  CMP Process Performance * Rasband, W.S., ImageJ, U.S. National Institutes of Health, Bethesda, Maryland, USA, https://imagej.nih.gov/ij/, 1997-2018. 3/ FIJI  Surface Analysis Script  Die Extraction  Image Thresholding  White Pixels Extraction  White / Black Pixels Ratio Unbonded %Area Mask Design Impact 1/ High Resolution SAM Image Acquisition 4/ Unbonded Arrays Identification using die layout file  Design Rules While C1 Cu CMP process do show significant reduction of Unbonded %Area, meaning better bonding performance, C2 Cu CMP process lead to slight increase of white, unbonded pixel ratio: bonding performance of C2 process is worse Changing the Barrier CMP Process, from B to B1, helps to reduce Unbonded %Area  HR-SAM Image Analysis  Detect CMP process performance  Choose best process for bonding Similar Unbonded %Area at die scale for the 3 Cu CMP processes  White pixel distribution highlight repetitive unbonded Cu arrays  Mask design importance in terms of pitch & density Changing Barrier CMP process from B to B1 show significant white pixels reduction Barrier CMP have local impact, at Die scale, allowing to increase bonding process window in terms of density and pitch Center Die POR

Die Scale Quantification of Bonding Quality Unbonded Cu Arrays Identification  Based on mask design file (left-side), unbonded Cu arrays were identified and counted on HR-SAM Image (middle)  Graph on right-side shows that, for these specific conditions, a large number of Cu-arrays (white bar on graph) were not bonded:

 ~ 56% unbonded arrays (white bar)  ~38 bonded arrays (black bar)  ~6% mixed (grey bar) : these arrays, on limit between bonded/unbonded are important to focus on for CMP process optimization