# Numerical Study of Density Ratio Influence on Global Wave Shapes Before Impact

Stéphane Etienne, Yves-Marie Scolan, Laurent Brosset

► To cite this version:

Stéphane Etienne, Yves-Marie Scolan, Laurent Brosset. Numerical Study of Density Ratio Influence on Global Wave Shapes Before Impact. 37th International Conference on Ocean, Offshore and Artic Engineering (OMAE2018), Jun 2018, Madrid, Spain. 10.1115/OMAE2018-78624. http-01867902

HAL Id: hal-01867902

https://ensta-bretagne.hal.science/hal-01867902

Submitted on 16 Oct 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# NUMERICAL STUDY OF DENSITY RATIO INFLUENCE ON GLOBAL WAVE SHAPES BEFORE IMPACT

**Stéphane Etienne**

Ecole Polytechnique Montréal
Montréal, Québec
Canada

**Yves-Marie Scolan**

ENSTA-Bretagne
Brest
France

**Laurent Brosset**

GTT (Gaztransport & Technigaz)
Saint-Rémy-lès-Chevreuse
France

## ABSTRACT

The influence of the gas-to-liquid density ratio (DR) on the global wave shape before impact is studied through numerical simulations of the propagation of two different waves in a rectangular wave canal. Two different codes are used: the first one, named FSID, is a highly non-linear 2D bi-fluid potential code initially developed in the frame of SLOSHEL JIP (Kaminski et al. (2011)) to simulate incompressible inviscid free-surface flows without surface tension thanks to a desingularized technique and series of conformal mappings; the second one, named CADYF, is a bi-fluid high-fidelity front-tracking software developed by Ecole Polytechnique Montreal to simulate separated two-phase incompressible viscous flows with surface tension.

The first studied wave leads to a flip-through impact while the second one leads to a large gas-pocket impact. Each condition is studied with water and three different gases with increasing densities corresponding to DR=0.001, 0.003 and 0.005. The global wave shapes are compared a few tenths of second before the impact, before free surface instabilities triggered by the shearing gas flow have developed and also before any gas compressibility matters. Both codes give precisely the same global wave shapes.

Whatever the condition studied, it is shown that DR has an influence on these global wave shapes. The trends observed from the simulations are the same as those described in Karimi et al. (2016) obtained from sloshing model tests with Single Impact Waves (SIW) in a 2D tank with a low filling level. A small part of the mechanical energy of the liquid is progressively given to the gas. The larger the DR, the larger this transfer of energy from the liquid to the gas. This explains an increasing delay of the wave front for increasing DRs.

## CONTEXT AND OBJECTIVES

### Context

For any new project of a floating structure equipped with membrane Liquefied Natural Gas (LNG) tanks (LNG carriers, offshore LNG terminals such as FLNGs or FSRUs, LFS, LNG bunker ships) the dominant design loads of the LNG containment system derive from liquid impacts due to sloshing and are determined from sloshing model tests, usually at scale 1:40, performed with water. The motions

of the floating structure are calculated at scale 1 taking into account the coupling with the cargo motions. After adequately down-scaling these motions (as the gravity is the same at both scales, the time scale needs to be the square root of the geometrical scale), those motions are imposed by a six degree-of-freedom hexapod (Stewart platform) to the model tank. As obviously the fluids (liquid and gas) in the model tank cannot have all their properties relevantly scaled with regard to those of LNG and Natural Gas (NG), the question arose about a hierarchy between them or between the related dimensionless numbers in order to bias as less as possible the similarity between the flows at both scales. Assuming the liquid in the model tank is water, the first property of the gas to look at is its density, therefore the gas-to-liquid density ratio (DR).

From sloshing model tests performed with water and different ullage gases or vapor (Maillard and Brosset (2009), Ahn et al. (2012)), it has been shown that, statistically and for all level of probability to be considered, the heavier the ullage gas, the smaller the pressures. Based on this result, sloshing model tests for any project of LNG floating tank are now performed with a heavy gas made of a mixture of SF₆ and N₂ tuned in order to match the DR in real LNG tanks (close to 0.004).

Actually the reduction of the statistical pressures when using heavier gases should not be attributed only to the DR. Indeed all gas properties are involved during liquid impacts. For instance, gas compressibility is involved just before impact, when the gas cannot escape quickly enough to cope with the remaining space left by the advance of the wave, and during the impact when gas pockets are entrapped.

As the compressibility of a heavier gas tends to be larger and as any gas at small scale is far too stiff with regard to the relevantly scaled compressibility of NG, having a heavier gas in the model tank is more representative of the reality also from the gas compressibility point of view.

Nevertheless, as it is almost impossible to change the density of a gas keeping the other properties unchanged, the simplest way to look at the single influence of the DR, is to observe its influence before the other properties matter, therefore before any impact. Thus, the
comparison is limited to the wave shape for Single Impact Waves (SIW) starting from rest, before any start of the compression of the gas.

Such comparisons are described in Karimi et al. (2016) based on sloshing model tests with a 2D tank at scale 1:20 representing a transverse slice of the tank#2 of a LNG carrier filled at 20% of the tank height. Several SIW conditions obtained by short sway motions with two different liquids and many different ullage gases have been studied. The wave shapes before impact were recorded by a high speed camera and precisely compared when using the different liquids and gases before any compression of the gas. Two areas were distinguished for the wave shape comparison: (1) the wave front where the free surface keeps smooth and precisely repeatable, unaffected by any development of free surface instabilities, namely from the trough to the base of the crest. This area is referred to as the *global wave shape*; (2) the area around the crest where a plume of free surface instabilities with liquid ligaments or liquid sheets and droplets develop under the action of the shearing gas flow. The local flow highly varies in that area keeping the same signature (global characteristics of the instable structures) when repeating precisely the same global wave shape. This is the reason for the high variability of the local pressure measurements at the wall (see Frihat et al. (2016)).

Some results extracted from Karimi et al. (2016) are presented in Fig. 1.

<figure><img src="image_2.png" /><figcaption>**FIGURE 1.** Sloshing tests in a 2D tank filled at 20% of the tank height and excited with a SIW condition. Wave shape before impact at the same instant for increasing DRs from left to right. (b) and (c) correspond to the same DR with liquids of different densities.</figcaption></figure>

When repeating the same SIW excitation with two different liquids of different densities, water (1 000 kg/m³) or a solution of Sodium Poly-tungstate (PST) (1 800 kg/m³), but with respective ullage gases chosen in order to get the same DR, the *global wave shape* remains precisely the same (Fig. 1. (b) and (c)). Therefore, this shape is independent of the liquid density and only depends on DR.

When comparing the same SIW condition with increasing DRs (see Fig. 1. (a), (b) and (d)), the wave front (below the crest region) becomes more inclined backward with regard to the impacted wall as though the breaking process was progressively impeded by the heavier gas.

It can also be observed in Fig. 1., that the characteristics of the local wave shape around the crest are modified. As these characteristics are highly dependent on the global wave shape itself when keeping the DR unchanged, the direct local influence of DR can hardly be discriminated from the consequence of its influence on the global wave shape. However, this consequence is important to notice

because it leads to modifications of the statistical distributions of local pressures when repeating many times the same SIW.

## Objectives

In the present paper, focus will be made on the influence of DR on the global wave shapes of SIWs in a wave canal as studied by numerical simulations with two different codes before any compression of the gas and before large development of free surface instabilities. The codes are both incompressible bi-fluid solvers.

Four different SIWs have been studied with the two codes for five DRs. For the sake of brevity, only the results for two waves and for three DRs will be presented within this paper but illustrating the general trends. The wave profiles at different times just before impact and the time traces of the different components of the energy will be provided.

It will be shown that the difference between the results of the two codes is small enough with regard to the difference for two different studied DRs with the same code. Therefore any of them is adapted for studying the influence of DR on the global wave shape and their relevance is reinforced by each other. Furthermore the same qualitative trends as observed experimentally by Karimi et al. (2016) on the influence of DR on the global wave shape are also observed numerically. The evolution of the energy distribution between liquid and gas enables to understand the influence of DR on the global wave shape and confirms the explanations already given in Braeunig et al. (2009).

The two different codes and the different calculation cases used for this study are described respectively in the next two sections. The results and analysis are presented in the last two sections.

## THE TWO NUMERICAL MODELS

### FSID

FSID stands for Free-Surface IDentification. It was initially developed in the frame of SLOSHEL project (see Kaminski et al. 2011). It is further developed by the second author for various applications in Naval Architecture and Coastal Engineering where breaking waves are expected. FSID simulates highly nonlinear two-dimensional two-fluid flows separated by a single continuous interface in the frame of potential theory, assuming therefore that the flow in both fluids is irrotational and incompressible. The surface tension is also not taken into account into the model. In practice one fluid is a liquid and the other is a gas. The above restrictive assumptions are considered as well fulfilled during the generation and further propagation of a water wave in presence of any gas or of a LNG wave in presence of NG before any impact occurs and more precisely before any development of free surface instability due to the shearing gas flow, which happens generally before any compression of the gas.

The fluid domain D is a rectangle possibly restricted by chamfers or by different kinds of solid shapes adjacent to the walls, including a quarter of ellipse, lying on the floor beside one of the vertical walls. This enables to describe for instance a transverse section of a membrane LNG tank on a floating structure or a wave canal with a given bathymetry adjacent to the impacted wall favoring the breaking process. The fluid domain is split in two compact domains $D_1$ (for the liquid) and $D_g$ (for the gas) separated by a single continuous interface at any time.

The tank might be animated by three-degree-of-freedom forced motions in its plane. Assuming it represents a transverse section of a
tank of LNG carrier; the three degrees of freedom define sway, heave and roll motions. In that case, the velocity potential, defined for each fluid in a reference system attached to the earth, is decomposed in two parts: a relative velocity potential (relative velocities with regard to a reference system attached to the tank) following a Neumann condition at the walls and a complementary potential associated to the velocity induced by the solid motion of the reference system attached to the tank with regard to the fixed reference system. This complementary potential includes the Stokes-Joukowski potential (Joukowski (1885)) related to the roll-induced solid motion.

The problem could be directly solved with a standard panel method. In practice the desingularized technique described in Cao et al. (1991) and Tuck (1998) is used: for each fluid problem (liquid or gas), the different components of the potential are represented by a finite set of sources (Green function of Rankine type). These sources are located at a short distance from the interface outside the considered fluid domain. Special attention has been paid on the relevant choice of the desingularizing distance.

A succession of conformal mappings is applied so that the Green functions $G_f$ and $G_g$ attached respectively to the liquid domain and to the gas domain satisfy the boundary conditions on the walls of $D_f$ and $D_g$ respectively.

These two combined techniques (desingularized technique and series of conformal mappings) make the code fast enough to perform direct visualizations of the flow for any SIW generation on standard laptops.

The code is not suited for dealing with strong gas flow shearing the free surface. The simulation of a SIW stops as soon as the shearing gas flow tends to generate Kelvin-Helmholtz-type free surface instabilities. This generally happens soon before any wave impact against a wall and before any compression of the gas could actually happen.

Comparisons between wave shapes captured with a high speed camera during sloshing tests and wave shapes as calculated by FSID in the same conditions showed good agreement for one-, two- or three-degree-of-freedom forced motions when considering low filling levels with water and air in a rectangular tank as long as the basic assumptions of FSID are still valid, namely as long as the free surface remains continuous (Scolan and Brosset (2017)).

That is why GTT has been using FSID for years with different partners in order to quickly generate relevant inflow conditions before impact for further impact studies by using more sophisticated CFD codes taking into account physical properties of the fluids such as compressibility, viscosity, surface tension or phase change.

More details on the numerical model of FSID can be found in Scolan (2010), Scolan et al. (2016) and Scolan and Brosset (2017).

## CADYF

Developed by Ecole Polytechnique of Montreal (EPM), CADYF is a high-precision front-tracking solver of the Navier-Stokes equations simulating separated viscous two-phase flows with surface tension. Adaptivity in space (adaptive remeshing) and time (hp-adaptivity) enable to yield accurate predictions while keeping computational cost low.

The incompressible Navier-Stokes equations are solved for two Newtonian fluids (usually a liquid and a gas) on deforming domains using an Arbitrary Lagrangian-Eulerian (ALE) formulation (Hay et al. (2014)). At the interface between the two immiscible fluids, mass and

momentum are conserved. The kinematic condition follows from mass conservation and, without phase change, indicates that the fluid velocity is continuous across the interface for viscous fluids. The dynamic condition follows from the momentum conservation law and states that forces acting on the fluid at the interface are in equilibrium. The discontinuity of pressure at the interface corresponds to the surface tension. Slip or no-slip boundary conditions can be applied on the outer boundaries of computational domains.

The set of equations is discretized in space using the finite-element method. A stabilized (PSPG-SUPG) mixed formulation is used resulting in a fully coupled solution procedure solved by a modified Newton method. The velocity and pressure variables are discretized using P1-P1 elements on adaptive grids generated by an adaptive remeshing algorithm coupled with an automatic mesh generation procedure. This allows for the simulation of extremely large interfacial deformations as those induced by the shearing gas flow before a wave impact.

The adaptive remeshing algorithm monitors the maximum mesh deformation $|E|$ from the pseudo-solid approach. A new calculation grid is automatically generated as soon as $max(|E|)$ becomes larger than a chosen threshold $E_{max}$. The mesh generation procedure is based on an advancing front technique for which local size of elements are defined from a field of mesh sizes stored over a background grid.

Time integration is performed by a hp-adaptive algorithm based on the Backward Differentiation Formulas (BDF) which delivers solutions of prescribed accuracy while optimizing computational efficiency. In practice, the adaptive time-stepping procedure automatically selects the integration step size and order to guarantee that the solution time error is below the user selected error tolerance $\epsilon$ (Hay et al. (2015)).

The key ingredient of the present numerical model is a front tracking approach in which interfaces are aligned with the mesh edges. It yields accurate predictions of interfacial flows by preventing any diffusion of interfaces and allowing for fine physical modeling at interfaces. The deformation of interfaces with time induces grid deformations that are naturally accounted for by the ALE formulation. The pseudo-solid approach is used to propagate the interface displacement with time throughout the computational domain. The jump conditions arising from the kinematic and dynamic conditions are implemented by using 1D zero-thickness interfacial elements.

In the interface tracking approach the interface must be a material surface in the normal direction. In the framework of the ALE formulation this means that, in the normal direction, the ALE velocity is equal to the fluid velocity at the interface. The choice of the ALE velocity in the tangential direction to the interface is arbitrary. Similarly to what has been done for free surface approaches, we set the time evolution of the pseudo solid displacement in the tangential direction so as to preserve the mesh regularity along the interface. In practice the node displacement along the interface is chosen such that the normalized curvilinear abscissa of interfacial nodes is kept constant.

The explicit representation of interfaces in the front tracking approach allows for a very accurate description and prediction of interfacial flows for a low computational time. This is its main advantage over front capturing (e.g. VOF, level-set) or meshless methods such as SPH. However, its main limitation is that it cannot directly simulate very large deformation of interfaces or change of topology of interfaces. This limitation can be largely alleviated by using mesh adaptation.
Comparisons between wave shapes captured with a high speed camera during sloshing tests and wave shapes as calculated by CADYF in the same conditions showed good agreement for forced motions of a rectangular tank with a low filling level of water and air as ullage gas, at least until the development of local perturbations on the interface (Hay et al. (2016)).

The ability of the code to capture adequately the growth of free surface instability such as Kelvin-Helmholtz's, at least for academic conditions is shown in Fortin et al. (2018).

# CALCULATION CASES

## History of the calculation cases

The generation and propagation of four unidirectional waves have been simulated in 2D by CADYF and FSID for five density ratios (20 cases for each code). These wave cases have been first proposed in Guilcher et al. (2013) and Guilcher et al. (2014) with water and air. They lead to a flip-through impact (denoted FTI) and three gas-pocket impacts with entrapped gas-pockets of different sizes (the impact with the largest gas pocket is denoted LGPI). For the sake of brevity, only the results corresponding to the waves leading to the FTI and to the LGPI are presented here with only three DRs (6 cases for each code). They are sufficient to illustrate the general conclusions.

In Guilcher et al. (2013) and Guilcher et al. (2014), the focus was made on the impact simulation by a SPH code solving the compressible Euler equations for both the liquid and the gas. The pressure maps $p(y, t)$ at the wall (y measures a distance along the wall, $t$ is the time and $p$ the pressure) were compared for the different studied conditions. The wave generation and propagation were already carried out by FSID but with a mono-fluid version at that time. The SPH simulations were initiated with the liquid velocity and pressure fields as calculated by FSID in the liquid prior to the impact time and with a gas at rest. This coupling strategy between an incompressible solver for the wave generation and a compressible solver for the impact simulation saved a large amount of computational time. It has been adopted by different authors to simulate either the FTI or the LGPI or both (Costes et al. (2013), Behruzi et al. (2016)).

In the present paper, the focus is made on the wave generation and propagation before the impact and before any compressibility influence in the gas in order to compare the wave profiles before impact for different DRs.

A complete CADYF simulation of the FTI (propagation and impact) with water and air has already been presented in Hay et al. (2016) showing good comparison with results from Guilcher et al. (2014) in terms of pressure map.

FSID simulations of the wave propagation leading to the LGPI have already been presented in Scolan and Brosset (2017) for water without gas or with gas for a DR=0.001.

In the following, the notations FTI and LGPI will refer both to the wave generation and to the impact itself.

## Definition of the calculation cases

A 20 m by 12 m rectangular tank is considered. At the left bottom corner defined as the origin O of a reference system (x, y), there is a local bathymetry made of a quarter of an ellipse whose center is O, the main half-dimensions of which are 18 m and 2.8 m. The tank is filled with water and a gas, both considered as incompressible.

The simulation is artificially started with a given initial interface

shape considering both fluids at rest. The shape of the initial interface is defined by:

$$y = h + a \tanh[r(x - x_p x_l)] \quad (1)$$

The tank geometry and an initial position of the interface corresponding to the FTI are presented in Fig. 2.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Geometry and initial condition for the FTI.</figcaption></figure>

The different parameters needed for the definition of the initial wave shape for FTI and LGPI are given in Table 1.

<table><caption>TABLE 1. Parameters defining the initial wave shape for FTI and LGPI according to (1).</caption><thead><tr><th></th><th>h (m)</th><th>a (m)</th><th>x<sub>p</sub></th><th>x<sub>l</sub> (m)</th><th>r (m<sup>-1</sup>)</th></tr></thead><tbody><tr><th>FTI</th><td>7.6</td><td>3.6</td><td>2.5</td><td>20.0</td><td>0.36</td></tr><tr><th>LGPI</th><td>7.6</td><td>3.6</td><td>3.1</td><td>20.0</td><td>0.44</td></tr></tbody></table>

The fluid material properties of water and of the three gases used in the simulations are given in Table 2.

<table><caption>TABLE 2. Fluid material properties.</caption><thead><tr><th></th><th>Density<br>(kg.m<sup>-3</sup>)</th><th>Viscosity*<br>(Pa.s)</th><th>Surface tension**<br>(N.m<sup>-1</sup>)</th></tr></thead><tbody><tr><th>Water</th><td>1000</td><td>0.001</td><td>/</td></tr><tr><th>Gas1</th><td>1.0</td><td>1.75 10<sup>-5</sup></td><td>0.072</td></tr><tr><th>Gas2</th><td>3.0</td><td>1.75 10<sup>-5</sup></td><td>0.072</td></tr><tr><th>Gas3</th><td>5.0</td><td>1.75 10<sup>-5</sup></td><td>0.072</td></tr></tbody></table>

*Dynamic viscosity

**Surface tension at the interface between water and any studied gas

Only the density of the water and the density of the gases are used by FSID. The viscosity and tension surface are only used by CADYF. It will be shown later by comparison with FSID results that the low values of surface tension and viscosity used for CADYF simulations have no influence on the global wave shape. There would have an influence on the development of free surface instabilities triggered by the shearing gas flow just before impact. But studying those instabilities is not the purpose of the present paper.

For both codes a slip boundary conditions is applied on the outer boundaries of the computational domains.

# RESULTS

Let us define the different components of the mechanical energy into the different fluids:
* $E_p^l$ is the potential energy into the liquid,

* $E_k^l$ is the kinetic energy into the liquid,

* $E_p^g$ is the potential energy into the gas,

* $E_k^g$ is the kinetic energy into the gas,

* $E_m^l = E_p^l + E_k^l$ is the mechanical energy into the liquid,

* $E_m^g = E_p^g + E_k^g$ is the mechanical energy into the gas,

* $E_p = E_p^l + E_p^g$ is the total potential energy,

* $E_k = E_k^l + E_k^g$ is the total kinetic energy,

* $E_m = E_p + E_k = E_m^l + E_m^g = E_p^l + E_k^l + E_p^g + E_k^g$ is the total mechanical energy.

The potential energy for each fluid is considered as null when both fluids are at rest with a flat horizontal interface between them. Whatever the studied wave, as the wave generation started with a chosen initial non-horizontal wave shape with both fluids at rest, the initial total mechanical energy is the total initial potential energy:

$$E_m(0) = E_p(0) = E_p^l(0) + E_p^g(0) \quad (2)$$

The different components of the energy are denoted with a ~ when normalized by the total available energy $E_m(0)$. Thus, $\tilde{E}_m(0) = 1$.

Figure 3. and Fig. 4. show the time traces of the normalized energy components, $\tilde{E}_p^l$ (blue solid line), $\tilde{E}_k^l$ (blue dotted line), $\tilde{E}_p^g$ (green solid line), $\tilde{E}_k^g$ (green dotted line) and $\tilde{E}_m$ (black line) as calculated by CADYF for respectively the FTI and the LGPI and for DR=0.005. The corresponding curves as obtained by FSID are not provided on the graphs. The differences between the two codes could hardly be discriminated.

Whatever the wave, as the DR is small, the mechanical energy into the gas $E_m^g(t)$ remains a very small proportion of the total initial mechanical energy $E_m(0)$ during the wave propagation. Thus, the evolutions of $\tilde{E}_p^g(t)$ and $\tilde{E}_k^g(t)$ which are crucial to explain the influence of DR on the global wave shape before impact can hardly be distinguished on Fig. 3. and Fig. 4.

Therefore, we define the change of respectively liquid, gas and total mechanical energy by:

$$\Delta E_m^l(t) = E_m^l(t) - E_m^l(0) \quad (3)$$

$$\Delta E_m^g(t) = E_m^g(t) - E_m^g(0) \quad (4)$$

$$\Delta E_m(t) = E_m(t) - E_m(0) \quad (5)$$

Figure 5. and Fig. 8. show the time evolutions of $\Delta \tilde{E}_m^g$ (green lines), $\Delta \tilde{E}_m^l$ (blue lines) and $\Delta \tilde{E}_m$ (black lines) for DR=0.001, DR=0.003 and DR=0.005 for respectively the FTI and the LGPI. The three curves of the same color are differentiated by a thicker line for a larger DR. Fig. 5. provides the curves as calculated by CADYF (solid lines) and by FSID (dotted lines) while Fig. 8 provides the curves only as calculated by CADYF.

Figure 6. and Fig. 9. respectively show the global wave profiles at two different instants for the FTI (t=1.47 s and t=1.67 s) and for the LGPI (t=1.82 s and t=2.02 s) as calculated by CADYF for DR=0.001 (orange lines), DR=0.003 (red lines) and DR=0.005 (pink lines).

Figure 7. and Fig. 10. present a close up of the wave profiles in the area of the impacted wall for respectively the FTI and the LGPI at the same instants as for respectively Fig. 6. and Fig. 9. and for the three

DRs with the same code of colors as previously. As the close-up enables to discriminate them, the wave profiles as calculated both by CADYF (solid lines) and by FSID (dotted lines) have been displayed in the figures when available.

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Time traces of $\tilde{E}_p^l$ (blue solid line), $\tilde{E}_k^l$ (blue dotted line), $\tilde{E}_p^g$ (green solid line), $\tilde{E}_k^g$ (green dotted line) and $\tilde{E}_m$ (black line) for the FTI and DR=0.005, as calculated by CADYF.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Time traces of $\tilde{E}_p^l$ (blue solid line), $\tilde{E}_k^l$ (blue dotted line), $\tilde{E}_p^g$ (green solid line), $\tilde{E}_k^g$ (green dotted line) and $\tilde{E}_m$ (black line) for the LGPI and DR=0.005, as calculated by CADYF.</figcaption></figure>

## ANALYSIS

### Transfer between the components of the liquid mechanical energy

In Fig. 3. and Fig. 4., we observe that the liquid is progressively animated thanks to a transfer from its potential energy to its kinetic energy. The evolution of this transfer depends on the initial wave shape.
Flip-Through Impact

<figure><img src="image_6.png" /></figure>

FIGURE 5. Time traces of $\Delta\tilde{E}_m^g$ (green lines), $\Delta\tilde{E}_m^l$ (blue lines) and $\Delta\tilde{E}_m$ (black lines) for DR=0.001 (thin line), DR=0.003 (intermediate thickness), DR=0.005 (thick line) as calculated by CADYF (solid lines) and FSID (dotted lines) for the FTI.

<figure><img src="image_7.png" /></figure>

FIGURE 6. FTI: global wave profiles as obtained by CADYF for DR=0.001 (orange), DR=0.003 (red) and DR=0.005 (pink) at two instants t=1.47 s and t=1.67 s.

<figure><img src="image_8.png" /></figure>

FIGURE 7. FTI: close-up of the wave profiles near the impacted wall as calculated by CADYF (solid lines) and FSID (dotted lines).

Large Gas-Pocket Impact

<figure><img src="image_9.png" /></figure>

FIGURE 8. Time traces of $\Delta\tilde{E}_m^g$ (green lines), $\Delta\tilde{E}_m^l$ (blue lines) and $\Delta\tilde{E}_m$ (black lines) for DR=0.001 (thin line), DR=0.003 (intermediate thickness), DR=0.005 (thick line) as calculated by CADYF for the LGPI.

<figure><img src="image_10.png" /></figure>

FIGURE 9. LGPI: global wave profiles as obtained by CADYF for DR=0.001 (orange), DR=0.003 (red) and DR=0.005 (pink) at two instants t=1.82 s and t=2.02 s.

<figure><img src="image_11.png" /></figure>

FIGURE 10. LGPI: close-up of wave profiles near the impacted wall as calculated by CADYF (solid lines) and FSID for DR=0.001 (dotted line).
For the FTI the horizontal velocity of the wave front and the
vertical velocity of the wave trough are progressively increasing until
the impact leading to a progressive increase of the liquid kinetic
energy while the global profile of the free surface is flattening because
the amplitude of the wave front is progressively decreasing (Fig. 6.)
leading to a progressive decrease of the liquid potential energy almost
down to zero. A flip-through impact corresponds to an almost
complete transfer from potential to kinetic energy during the
generation of the wave.

For the progressive wave leading to a gas pocket impact (LGPI),
the liquid kinetic energy and the liquid potential energy reach
respectively a maximum and a minimum value. These values remain
almost constant as soon as the wave front is formed. This can be
understood by the fact that there are only small variations of the
horizontal velocity of the wave front and that the vertical velocity of
the trough is much smaller than the horizontal velocity of the front.
This leads to only a small decrease of the amplitude of the wave
through a slow progressive elevation of its trough until the impact
(Fig. 9.). Actually, it was checked that the smaller the size of the
entrapped gas pocket, the smaller the plateau of remaining potential
energy. Therefore, the transfer from potential energy to kinetic energy
is more and more complete when leading to smaller gas-pocket
impacts. This is important to notice because smaller gas pocket
impacts induce larger impact pressures.

Starting from rest, the gas will also be progressively animated
taking a part of the mechanical energy of the liquid. We will show
later that this amount of mechanical energy transferred from the liquid
to the gas depends directly on DR and explains the differences on the
global wave shape before impact for different DRs. For the time being,
due to the range of energy displayed for both waves, Fig. 3. and Fig. 4.
only show that the global share of mechanical energy taking by the gas
is very small with regard to that of the liquid. It should be noted that
the figures are given with the highest DR studied for which the share
of the mechanical energy taken by the gas is the largest.

The total mechanical energy as shown on Fig. 3. and Fig. 4. looks
very much constant for both waves. This can be checked more
carefully by looking at ΔĒ<sub>m</sub> in Fig. 5. and Fig. 8. For the FTI, ΔĒ<sub>m</sub>
remains lower than 10<sup>-5</sup> whatever the code and whatever the DR. For
the LGPI, Fig. 8. shows a quick drop of ΔĒ<sub>m</sub> with CADYF around
t=2 s whatever the DR. ΔĒ<sub>m</sub> remains lower than 10<sup>-5</sup> until the sudden
growth of kinetic energy into the gas around this time, whatever the
DR. The energy dissipation observed here is mainly a numerical
artefact: numerical energy dissipation is cumulated at each remeshing
operation mainly due to a too simple linear interpolation scheme used
in CADYF. The drop of mechanical energy observed on the figure
corresponds to the sudden increase of the number of remeshing
operations when the curvature of the interface becomes very high at
some locations due to the development of free surface instabilities.
The development of a higher order interpolation scheme will soon
solve this issue. Anyway, this issue starts at a time when the maximum
velocity of the gas flow in between the crest and the wall is already
very high (around 200 m/s obtained by CADYF at t=2 s) and for
which the incompressibility assumption is not any longer valid.

This means that the dissipation, either due to the actual viscosity
influence or to any numerical artefact remains very low as long as the
incompressibility assumption is valid.

**Comparison between CADYF and FSID results**

We can observe in Fig. 7. and Fig. 10. that the differences on the
global wave shape before impact as obtained by CADYF and FSID

simulations are hardly visible and much smaller than the differences
generated by the influence of DR that we wanted to capture.

FSID ensures the continuity of normal velocity (and pressure) but
there is a discontinuity of tangential velocity at the interface that is
calculated. The program stops the simulation as soon as this
discontinuity becomes too large and perturbations at the interface start
to grow. It would be possible to delay the time when the program stops
by artificially reducing the resolution at the interface (smaller number
of Lagrangian markers). This would inhibit the development of large
jumps of the tangential velocity at the interface but would be paid by
worse energy/mass conservation. The development of free surface
instabilities is inherent to liquid impacts. FSID is intended to capture
quickly and precisely the global wave shape until this first
development.

On the other hand, the last developments of CADYF with the
adaptive mesh refinement depending on the proximity to the free
surface and on its local curvature (Hay et al. (2016), Fortin et al.
(2018)) have been introduced especially to enable the capture of
instabilities at the interface or at least their initial stage of development
before fragmentation. In reality, even for a globally unidirectional
wave (2D global flow), the instabilities soon develop in 3D with the
apparition of ligaments, thin liquid sheets and droplets (see Fig. 1.).
Therefore, only 3D simulations with CADYF could, in the best case,
able a realistic capture of the free surface instabilities prior to a
liquid impact if carried out with extreme mesh refinements. As this is
not our purpose here, CADYF is only limited by a possible change of
the liquid domain topology which happens for instance when the tip of
a wave crest hits a wall. Therefore, for the studied cases, CADYF is
able to capture the global wave shape until the impact.

For the FTI, FSID was able to simulate the flow for the three DRs
until t=1.67 s, therefore until the start of the pressure rise at the wall.
For the LGPI, the last moment as simulated by FSID for DR=0.001
was at t=1.82 s, a few tenths of second before the impact. This
explains that only the wave profile obtained by FSID at t=1.82 s for
DR=0.001 is displayed in Fig. 9. This last simulated time might be
later for larger DRs. For instance it was possible to simulate the flow
until t=2.02 s for DR=0.003. At each time the comparison was
possible the global wave shapes obtained by both codes were almost
superimposed.

The earlier triggering of instabilities for breaking waves than for
non-breaking waves (slosh waves including flip-through) has already
been observed during wave impact tests in flumes and during sloshing
model tests for SIWs. The growth of the perturbations depends on the
shearing gas flow at the interface. For slosh waves, the free surface
might remain smooth until the end of the run-up of the jet along the
vertical wall. For breaking waves, this shearing gas flow is associated
to strong vortices located in a gas layer close to the free surface. These
vortices are well captured by CADYF as can be seen in Fig. 11 for the
LGPI and DR=0.001 at t=1.82 s. This vorticity layer remains thin from
the trough to the base of the crest where the free surface remains
smooth. It becomes larger at the crest level and on the shoulder of the
wave where the instabilities start to develop. The vortices cannot be
captured by FSID as it assumes an irrotational flow in both fluids.
With such an assumption, the larger the curvature of the interface, the
stronger the shearing gas flow and thus, the stronger and earlier the
development of instabilities.
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Vorticity field in the gas and in the liquid as calculated by CADYF for LGPI and DR=0.001 at t=1.82 s.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 12. Velocity magnitude in the gas and in the liquid as calculated by CADYF for LGPI and DR=0.001 at t=1.82 s. Vmax(liquid)=9.6 m/s, Vmax(gas)=32.2 m/s</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 13. Velocity vectors in the gas (green) and in the liquid (red) as calculated by FSID for LGPI and DR=0.001 at t=1.82 s. Vmax(liquid)=9.6 m/s, Vmax(gas)=22.2 m/s</figcaption></figure>

The velocity fields as calculated by CADYF and FSID for the LGPI and DR=0.001 at t=1.82 s are respectively represented in Fig. 12. and Fig. 13. The velocity field as obtained by CADYF is represented by a color plot of the velocity magnitude in Fig. 12. whereas the velocity field as obtained by FSID is represented by vectors in Fig. 13. Indeed, because of the adaptive mesh refinement algorithm used by CADYF, the mesh is refined close to the free surface and especially in the area of large vortices. For this reason a plot of the velocity vectors as obtained directly by CADYF would not be much legible.

The liquid flows as simulated by both codes are pretty much the same with a maximum velocity magnitude of 9.6 m/s in both Fig. 12. and Fig. 13. On the other hand the gas flows as simulated by both codes is different in the vorticity layer near the free surface. The maximum velocity magnitude obtained by CADYF is 32.2 m/s in Fig. 12. whereas that obtained by FSID is 22.2 m/s in Fig. 13.

These local discrepancies between the gas flows as calculated by the two codes explain the small discrepancies between the two estimations of $\Delta\vec{E}_m^g$ and thus between the two estimations of $\Delta\vec{E}_m^l$ that can be observed in Fig. 5. Anyway these discrepancies remain small with regard to the differences brought by the influence of DR. The influence of DR can therefore be studied in the same way through any set of curves obtained either by CADYF or FSID.

## Influence of DR on the global wave shape before impact

Whatever the wave studied, FTI or LGPI, the DR has an influence on its global wave shape before impact as can be observed on respectively Fig. 7. and Fig. 10. For increasing DRs, the following trends are observed:

1. at any time the wave front is delayed and the delay is increasing during the wave propagation,

2. at any time this delay is getting larger from the trough (with a vertical motion along the wall) to the crest (with an almost horizontal velocity) with the exception of the final stage of the flip-through when the trough velocity is becoming as large as the crest velocity,

3. the wave front becomes more inclined backward with regard to the impacted wall as though the breaking process was impeded,

4. the crest of the LGPI gets a rounder shape (smaller curvature).

Actually the third point is the direct consequence of the second point. The three first observations correspond to those described in Karimi et al. (2016) from sloshing model tests with water and gases for different DRs and recalled in the introduction. The fourth point does really make sense only when the shape of the wave crest has already a small curvature for a small DR. This can happen for breaking waves with thick crests. Only in that case, the free surface remains sufficiently smooth to enable the characterization of a global shape. In the other cases with a sharp crest, the global shape of the crest is hidden by a plume of droplets and small structures. However, this trend has been obtained with the two codes each time FSID was able to simulate a late development of the breaking wave.

The global delay of the wave front which is increasing during the propagation of the wave can be explained by the transfer of mechanical energy from the liquid to the gas as shown in Fig. 5. and Fig. 8. respectively for the FTI and the LGPI. Indeed, whatever the studied wave and the DR considered, the gain of mechanical energy of the gas corresponds to an equivalent loss from the liquid. This transfer
is increasing over time starting from zero at rest. The rate of this growth is almost constant for the FTI for which it has already been noticed that the kinetic energy of the liquid is increasing until the end (Fig. 3.), while there is a strong reduction of this rate for the LGPI (before the final surge of the escaping gas around t=2 s) for which a plateau of liquid kinetic energy is reached soon before the impact (Fig. 4). This tends to show that the transfer of energy from the liquid to the gas is especially a transfer of kinetic energy. The most important observation is that the larger the DR, the larger the transfer of mechanical energy (and likely kinetic energy) from the liquid to the gas. Therefore, the liquid is globally and progressively slowed down for increasing DRs which explains the observed global delay.

This is not a totally trivial result. Indeed, in the definition of the kinetic energy, the velocity which is counted to the square has more influence than the density. One could have imagined a reduction of the gas velocities for a larger gas density leading to a global reduction of the gas kinetic energy. As the opposite is observed, it means that in the range of gas density studied, the reduction of the gas velocities obtained with a larger gas density is small enough to be more than compensated by the increase of gas density and leads to an increase of kinetic energy.

Figure 14. shows the time evolutions of the velocity magnitudes as calculated by FSID on both sides of the interface along its curvilinear abscissa $s$ in the area of the impacted wall and in the last instants of the calculations for the FTI and for DR=0.001. The curvilinear abscissa of the interface starts from the impacted wall at any time ($s(t)=0$ at wall). The velocity magnitudes in the gas and in the liquid are respectively represented in green and red.

<figure><img src="image_15.png" /><figcaption>FIGURE 14. Velocity magnitudes (m/s) on both sides of the interface along its curvilinear abscissa $s$ (starting at the impacted wall for $s=0$) at different instants, for the FTI with DR=0.001, as calculated by FSID. Velocity in the gas and in the liquid respectively represented in green and red.</figcaption></figure>

There is a sudden rise of the liquid velocity magnitude close to the wall at the end of the calculation. This corresponds to the phenomenon of convergence of the wave front and the wave trough in a small area around the wall for a flip though impact. The magnitude of the gas velocity presents a bump along the curvilinear abscissa of the interface at any time which forms a wavelet on the green surface. Outside the area of the wavelet, the magnitude of the gas and liquid velocities are close. The maximum of the wavelet is moving toward the wall and its amplitude is progressively increasing. This maximum of the gas velocity magnitude is actually located at the tip of the wave front at

any time where the curvature of the interface is maximal. This is also clearly the place where the difference between the gas and liquid velocity magnitude is maximal. As the normal velocities are the same on both sides of the interface, the gap between the two surfaces is due to the tangential velocity of the gas. In the top right corners of the two surfaces corresponding to the end of the interface very close to the wall and at the end of the calculation, the green surface is progressively getting closer to the red surface. Indeed, this area represents the wave trough when the run-up is suddenly accelerating. Both the liquid and gas velocities are locally vertical due to the slip boundary condition at the vertical wall. Their magnitudes are the same on both sides of the interface due to the continuity of the normal velocity.

The same kind of graphs have been plotted and compared for the different DRs. All the graphs present the same characteristics as that of Fig. 14. Reduction rates of liquid and gas velocities are observed when comparing a higher DR to a smaller DR. Those rates depend on time and location. They are significant only in the areas where the velocities are important. Fig. 15 shows the velocity magnitudes on both sides of the interface along its curvilinear abscissa $s$ in the area of the impacted wall and for the five studied DRs and at the last instant of the calculations ($t=1.67$ s). The velocities in the liquid obtained with no gas (mono-fluid version of FSID) have been added as a reference.

<figure><img src="image_16.png" /><figcaption>FIGURE 15. Velocity magnitudes (m/s) as calculated by FSID on both sides of the interface along its curvilinear abscissa $s$ (starting at the impacted wall for $s=0$) at $t=1.67$ s, for the FTI with DR=0.001 (green), 0.002 (blue), 0.003 (pink), 0.004 (cyan), 0.005 (brown). Velocity in the gas and in the liquid respectively represented in thin and thick lines of the same color for the same DR. Velocity magnitude in the liquid for mono-fluid version of FSID in red.</figcaption></figure>

The reduction of the gas velocity is significant (27% between DR=0.005 and DR=0.001) but localized around the tip of the wave front which explains that globally, on the whole gas domain the change of mechanical energy into the gas is still able to increase. The reduction rate of the liquid velocity is larger at the tip of the wave front, where the liquid velocities are maximal, at any time and is increasing with time (this cannot be seen in Figure 15 which is only at $t=1.67$ s but has been checked). It becomes even larger at the extremity of the interface (the wave trough) during the start of the final run-up. This reduction rate of the liquid velocity magnitude at the wave trough observed in Fig. 15 at the final time is around 42% when comparing
DR=0.005 to DR=0.001. This explains the larger delay at the tip of the front than at the trough during the wave propagation and then the quick delay at the trough during the final stage of the flip-through when the trough velocity is becoming large. The same kind of analysis can be made from the LGPI results to explain the difference of behavior at the trough and at the crest but is not provided for the sake of brevity. Actually, even assuming a uniform reduction rate of the liquid velocity would be enough to explain a larger delay for a higher DR in the area of larger liquid velocity.

The difference of behavior between the trough and the crest and also at the end of the flip-through development suggests that the transfer of kinetic energy between the liquid and the gas is performed locally. During the wave propagation, a large share of the gas kinetic energy is located around the tip of the wave front. The corresponding kinetic energy retrieved from the liquid seems to be taken in the same area. In the late development of the FTI, the vertical velocity of the trough becomes such that the area of the trough becomes an important contributor to the gas kinetic energy. The corresponding kinetic energy retrieved from the liquid seems to be taken in the same area.

## CONCLUSIONS

Two codes have been used to simulate the generation and propagation of waves in a flume tank for different gas-to-liquid density ratios (DR). The wave generation depends only on the initial condition which is artificial but simple to implement. It consists in giving a shape to the free surface keeping the liquid and gas at rest. Thus the fluids are given an initial potential energy in an instable situation leading to the animation of the two fluids by a transfer from potential to kinetic energy.

The first code, FSID, solves the incompressible Euler equation for the liquid and gas in the frame of the potential theory. The second code solves the incompressible Navier-Stokes equations for the liquid and gas with a high-fidelity front-tracking approach with adaptive mesh refinement and taking into account surface tension.

The results concerning two different waves and three DRs have been shown but more calculations have been performed and the conclusions are general. One of the two waves leads to a flip-through kind of impact (FTI). The other leads two a large gas-pocket impact (LGPI). Both waves have already been studied in the literature with either a focus on the propagation phase or on the impact itself.

The objective of the present study was especially to compare the wave profiles for the different DRs before the impact and even before any other property of the fluids, as the gas compressibility, could matter.

The two codes gave very precisely the same wave profiles in the different conditions at any instant as long as FSID can deal with the shearing gas flow at the interface. The differences between the wave profiles are much smaller than the differences due to the influence of DR.

Whatever the condition, there are small differences between the simulations of the two codes on the velocity field when the wave front is close to the wall due to different assumptions on the continuity at the interface. These differences leads thus to small discrepancies on the kinetic energy into the gas and into the liquid when the gas flow becomes strong in between the wave and the wall. Nevertheless, here also these differences are small with regard to those induced by the change of DR considered. Therefore both codes are relevant for our objectives and all following conclusions could have been derived with

any of them.

Changing the DR modifies the wave profile before the impact and therefore modifies the nature of the impact as it had already been observed in Braeunig et al. (2009) from simple liquid impact simulations with a compressible bi-fluid solver. The modifications of the shape that have been observed by this numerical study match qualitatively well those described in Karimi et al. (2016) derived from high speed videos during 2D sloshing model tests for single impact waves (SIW).

When comparing the same wave with two different DRs, the wave front generated with the larger DR is progressively delayed with regard to the other. This delay is larger at the crest than at the trough as long as the trough velocity amplitude remains low with regard to the crest velocity amplitude. This leads to a backward inclination of the wave with regard to the vertical wall which is larger for larger DRs. This is the general case except for flip-through impacts for which the velocity of the trough increases suddenly at the end.

This delay and its spatial distribution can be explained by a transfer of mechanical energy from the liquid to the gas. Any gain of kinetic energy from the gas is taken locally from the liquid which therefore slightly slows down.

Nevertheless, increasing progressively the DR leads to a progressive increase of the gas kinetic energy which explains that the larger the DR, the larger the delay of the wave front at any time. This is not a totally trivial result: it means that in the range of DRs studied, the slight reduction of gas velocity for a larger DR is more than compensated by the larger density in the contribution to the gas kinetic energy.

Changing the DR modifies the wave shape before impact. The trends that have been listed above could as well magnify the impact pressures or mitigate them when increasing the DR depending on the initial reference wave shape chosen for a breaking wave impact. They are not totally sufficient to explain why, statistically, the higher the DR, the smaller the pressure are.

## REFERENCES

Ahn, Y., Kim, S. Y., Kim, K. H., Lee, S. W., Kim, Y. and Park, J. J. (2012). "Study on the Effect of Density Ratio of Liquid and Gas in Sloshing Experiment". In Proc. of the 22nd Int. Offshore and Polar Eng. Conf. (ISOPE), Rhodes, Greece. Vol. 3, pages 311-317.

Behruzi, P, Gaulke, D, Haacke, D and Brosset, L (2017). "Modeling of impact waves in LNG ship tanks", Int. Journal of Offshore and Polar Eng., Vol. 27, n°1, p 18-26.

Braeunig, J.-P., Brosset, L., Dias, F. and Ghidaglia, J.-M. (2009). "Phenomenological study of liquid impacts through 2D compressible two-fluid numerical simulations". In Proc. of the 19th Int. Offshore and Polar Eng. Conf. (ISOPE), Osaka, Japan. vol. 3, pages 21-29.

Cao, Y., Schultz, W.W. and Beck, R. F. (1991). "A three-dimensional desingularized boundary integral method for potential problems", Int. Journal of Numerical Methods in Fluids, Vol. 11 pp. 785-803.

Costes, J., Dias, F., Ghidaglia, J.-M. and Mrabet, A. (2013). "Simulation of Breaking Wave Impacts on a Flat Rigid Wall by a 2D Parallel Finite Volume Solver with Two Compressible Fluids and an Advanced Free Surface Reconstruction," in Proc. of the 23rd Int. Offshore and Polar Eng. Conf. (ISOPE), Anchorage, AK, USA, Vol 3, pp – 243-250.
Fortin, S., Etienne, S., Béguin, C., Pelletier, D. and Brosset, L. (2018). "Numerical Study of the Influence of Weber and Reynolds Numbers on the Development of Kelvin-Helmholtz Instability", *28th (2018) Int. Offshore and Polar Eng. Conf.* (ISOPE), Sapporo, Japan.

Frihat, M., Karimi, M.R., Brosset, L. and Ghidaglia, J.-M. (2016). "Variability of impact pressures induced by sloshing investigated through the concept of "singularization"", in Proc. of the 26th (2016) *Int. Offshore and Polar Eng. Conf.* (ISOPE), Rhodes, Greece, Vol. 3, p901-914.

Guilcher, P.-M., Couty, N., Brosset, L. and Le Touzé, D. (2013). "Simulations of breaking wave impacts on a rigid wall at two different scales with a two-phase fluid compressible SPH model", *Int. Journal of Offshore and Polar Eng.* Vol. 23, n° 4, 241-253.

Guilcher, P.-M., Jus, Y., Couty, N., Brosset, L., Scolan, Y.-M. and Le Touzé, D. (2014). "2D Simulations of breaking wave impacts on a flat rigid wall – Part 1: influence of the wave shape", in Proc. of the 24th (2014) *Int. Offshore and Polar Eng. Conf.* (ISOPE), Busan, South-Korea, Vol. 3, pp 232-245.

Hay, A., Yu, K. R., Etienne, S., Garon, A. and Pelletier, D. (2014). "High-order temporal accuracy for 3D finite-element ALE flow simulations". *Computers and Fluids*, Vol. 100, pp 204–217.

Hay, A., Etienne, S., Pelletier, D., and Garon, A. (2015). "hp-Adaptive time integration based on the BDF for viscous flows". *Journal of Computational Physics*, Vol. 291, pp 151–176.

Hay, A., Etienne, S., Pelletier, D. and Brosset, L. (2016). "Accurate prediction of sloshing waves in tanks by an adaptive two-fluid front-tracking approach", in Proc. of the 26th (2016) *Int. Offshore and Polar Eng. Conf.* (ISOPE), Rhodes, Greece, Vol. 3, p822-831.

Joukowski, N. (1885). "On motions of a rigid body with cavity filled by homogeneous liquid", *Journal of Russian Physical-Mathematical Society*, Vol. XVI, pp. 30-85.

Kaminski, M.L., Bogaert, H., Brosset, L. (2011), "Full and large scale wave impact tests for a better understanding of sloshing - Results of the Sloshel project", in Proc. of the 30th (2011) *Int. Conf. on Ocean, Offshore and Arctic Eng.* (OMAE), Rotterdam, The Netherlands.

Karimi, M.-R., Brosset, L., Ghidaglia, J.-M. and Kaminski, M.-L. (2016). "Effect of ullage gas on sloshing – Part II: local effects of gas-liquid density ratio". *European Journal of Mechanics B/Fluids*, Vol. 57, pp 82-100, DOI: 10.1016/j.euromechflu.2015.11.011.

Maillard, S. and Brosset, L. (2009). "Influence of density ratio between liquid and gas on sloshing model test results". *International Journal of Offshore and Polar Engineering*, 19(4): pp 271 - 279.

Scolan, Y.-M. (2010). "Some aspects of the flip-through phenomenon: A numerical study based on the desingularized technique". *Journal of Fluids and Structures*, Volume 26, Issue 6, pp 918–953.

Scolan, Y.-M., Karimi, M. R., Dias, F., Ghidaglia, J.-M. and Costes, J. (2014). "Highly nonlinear wave in tank with small density ratio", *29th Int. Workshop on Water Waves and Floating Bodies*, Osaka, Japan, April 2014.

Scolan, Y.-M., Hay, A. and Brosset, L. (2016). "Some aspects of high kinematics in breaking waves due to sloshing". *31st Int. Workshop on Water Waves and Floating Bodies*, Plymouth, USA, April 2016.

Scolan, Y.-M. and Brosset, L. (2017). "Numerical simulation of highly nonlinear sloshing in a tank due to forced motion", *International Journal of Offshore and Polar Engineering*, 27(1):11-17 March 2017, DOI: 10.17736/ijope.2017.jc678.

Tuck, E. O. (1998). "Solution of Nonlinear Free-Surface Problems by Boundary and Desingularised Integral Equation Techniques", Invited Lecture, *Computational Techniques and Applications: CTAC'97* Eds J. Noye, M. Teubner, A. Gill, World Scientific, Singapore, 11-26.