1. Introduction
To lower carbon emissions, the demand for natural gas (NG) as well as the sea transport of liquefied natural gas (LNG) has been increasing notably. Meanwhile, renewable and clean energy, such as solar energy, wind energy and hydrogen, represent a growing share in the global energy production. Among these, although hydrogen is gaining more attention, the production and transportation of liquified hydrogen (LH2) will increase significantly in the future. It is known that ships transport LNG at extremely low temperatures in membrane tanks, which are vulnerable to local sloshing-induced high impact pressures. Thus, impact loads acting on the walls of tanks partially filled with LH2 or LNG need to be considered. Although numerical methods are used to investigate the impact pressure induced by sloshing, not all physical phenomena are always modelled. This includes, for example, phase changes of liquid to gas and vice versa, strong mixing of liquid and gas phases that generates gas bubbles in the liquid, compressibility of the mixture consisting of liquid and gas and a complex form of liquid–gas interfaces with splash formation. In general, existing numerical methods have not been able to simulate all these physical phenomena adequately.
Research on sloshing phenomena has been carried out both experimentally and numerically for decades. Abramson (Reference Abramson1966) was the first to discuss scale effects, based on the Pi theorem, to establish the range of the parameters in model tests in terms of the similarity to full-scale conditions. Since the 1970s, Faltinsen devoted his theoretical and numerical research on sloshing. Faltinsen (Reference Faltinsen1974, Reference Faltinsen1977, Reference Faltinsen1978) developed a two-dimensional nonlinear theory of sloshing in rectangular tanks by predicting internal waves in closed basins. He then implemented his theory to obtain numerical solutions of the free-surface motion inside moving bodies. Faltinsen et al. (Reference Faltinsen, Rognebakke, Lukovsky and Timokha2000), Faltinsen & Timokha (Reference Faltinsen and Timokha2001) and Faltinsen, Rognebakke & Timokha (Reference Faltinsen, Rognebakke and Timokha2007) further developed numerical schemes based on the potential flow assumption. Miles (Reference Miles1972, Reference Miles1976) applied a modal analysis to theoretically develop free-surface waves in closed basins. Bagnold (Reference Bagnold1939) derived a model to predict wave-impact pressures with gas entrapment by studying the compression of the gas pocket, based on a one-dimensional assumption to calculate impact pressures on neighbouring walls. Bredmose & Bullock (Reference Bredmose and Bullock2008) investigated the scaling of wave-impact pressures in trapped air pockets by extending Bagnold’s approach for two and three axisymmetric pockets and comparing their predictions with numerical simulations of two-dimensional wave impacts. Using a Riemann solver, Bredmose, Peregrine & Bullock (Reference Bredmose, Peregrine and Bullock2009) numerically studied violent breaking wave impacts and the effect of air. Their model for a homogeneous mixture of an incompressible liquid phase and an ideal gas enabled simulation of the behaviour of trapped air pockets and the associated propagation of pressure waves. Bogaert et al. (Reference Bogaert, Léonard, Brosset and Kaminsk2010), by performing sloshing experiments at different scales to investigate scale effects of wave impacts, found that their results compared favourably with results obtained from the one-dimensional Bagnold-type model. Bredmose, Bullock & Hogg (Reference Bredmose, Bullock and Hogg2015) continued their work on violent breaking wave impacts to numerically and analytically study scale effects and ventilation. For several kinds of impacts, they showed that the effect of ventilation reduces the maximum impact pressure, the maximum force and the pressure impulse. Braeunig et al. (Reference Braeunig, Brosset, Dias and Ghidaglia2010) introduced a one-dimensional analytical model of a piston to determine the impact pressures caused by the collapse of a trapped gas bubble. They showed that the qualitative time histories are similar to those from model experiments. However, the application of their model requires a precise knowledge of all physical and thermodynamic parameters, which is not always possible in an experiment. Ancellin, Ghidaglia & Brosset (Reference Ancellin, Ghidaglia and Brosset2012) presented an analytical model that is more robust than that of Braeunig et al. (Reference Braeunig, Brosset, Dias and Ghidaglia2010). However, its application requires a high safety margin, and it is primarily suitable only for events involving trapped gas bubbles. Antuono & Lugni (Reference Antuono and Lugni2018) developed a modal method to predict global forces and moments for wave-induced sloshing in rectangular tanks. Results from their modal method compared well with experimental data obtained from sloshing tests. They extended their model to better predict the most violent wave breaking cases by introducing a diffusive term into the continuity equation. Fan et al. (Reference Fan, Jain and van der Meer2024) studied the air cushioning below impacting flat- and wave-structured disks. They analysed the free surface deformation, slamming loads and the air flow escaping from below the disk prior to impact. They found that the air radially flowing outwards across the liquid surface causes a free-surface condition that influences the details of the subsequent impact dynamics.
Olsen & Hysing (Reference Olsen and Hysing1974), by carrying out model tests to investigate the effects of viscosity of the liquid, found that viscosity does not have a dominant influence on slamming pressures. Bass et al. (Reference Bass, Bowles, Trudell, Navickas, Peck, Yoshimura, Endo and Pots1985), after reviewing modelling criteria for scaled LNG sloshing experiments, researched mainly the influence of viscosity, density ratio, ullage pressure and compressibility. Their results proved that the effects of viscosity and density ratio are insignificant for sloshing-induced impact pressures, whereas the compressibility of ullage space is perceptible. Lyu et al. (Reference Lyu, el Moctar, Potthoff and Neugebauer2017) conducted a series of simulations on the same aspects of sloshing and arrived at conclusions similar to Bass’s experimental studies. Zou et al. (Reference Zou, Wang and Cai2015a
,
Reference Zou, Wang, Cai and Lib
) carried out experimental and numerical studies to investigate the effects of the liquid’s viscosity (varying from 0.00152 to 0.25cSt) and the gas’s compressibility caused by violent sloshing for a roll excited tank filled with liquid (at filling levels of 20 %, 55 % and 70 %). Baudry & Rousset (Reference Baudry and Rousset2017) measured the impact pressure and the wave elevation of the liquid with high viscosities and compared these pressures with those using water. They observed that the amplitudes of wave elevation noticeably reduce as viscosity increases for a highly viscous liquid. Unlike water, such a liquid prevents breaking waves and splashing. Using a large-eddy-simulation approach, Jin & Lin (Reference Jin and Lin2019) carried out their research by considering a wide range of viscosities (varying from
$10^{-6}$
to
$10^{-2}$
m
$^2$
s−1), which covered also the viscosity of LNG. They concluded that, as viscosity increases, the maximum pressure decreases and the rise time increases. Jin et al. (Reference Jin, Tang, Tang, Mi, Wu, Liu and Huang2020), by experimentally studying the effect of viscosity (varying from
$1.14\times 10^{-6}$
to
$1.1\times 10^{-3}$
m
$^2$
s−1) at intermediate liquid depths, concluded that viscosity significantly affects the frequency response of the impact pressure and that high viscosity reduces nonlinearity due to the associated viscous damping. With regard to viscous effects on sloshing, no consistent conclusion is available regarding the influence of viscosity on sloshing-induced impact pressures.
Another aspect is the influence of the density ratio between gas and liquid on sloshing and sloshing-induced impact pressures. Although Bass et al. (Reference Bass, Bowles, Trudell, Navickas, Peck, Yoshimura, Endo and Pots1985) found no significant influence of density ratio on sloshing, several further investigations in recent decades arrived at opposite conclusions. Experimental observations of Ahn et al. (Reference Ahn, Kim, Kim, Lee, Kim and Park2012), obtained in a two-dimensional model tank of 1 : 40 scale and a three-dimensional model tank of 1 : 50 scale, filled with a mixture of sulphur hexafluoride (SF6) gas and liquid nitrogen, demonstrated that the sloshing-induced impact pressure decreases as the density ratio increases. Karimi & Brosset (Reference Karimi and Brosset2014) and Karimi et al. (Reference Karimi, Brosset, Ghidaglia and Kaminski2015, Reference Karimi, Brosset, Ghidaglia and Kaminski2016, Reference Karimi, Brosset, Kaminski and Ghidaglia2017) conducted two-dimensional experimental investigations at scales of 1 : 20 and 1 : 40 on density ratios varying from 0.0002 to 0.006 with irregular excitations. Globally, no evidence exists showing that the tested density ratio affects the wave shapes far from the impact zones, while locally the density ratio has significant consequences for the induced pressure as it evidently changes the shape of the breaking waves. Maillard & Brosset (Reference Maillard and Brosset2009) carried out model tests to research the influence of the density ratio between ullage gas and the liquid matters in a rectangular tank undergoing forced harmonic surge motions. The density ratio varied from 0.00005 to 0.0058. As this covered the density ratio of 0.004 for NG and LNG, it proved to be important to predict impact pressures. Yung et al. (Reference Yung, Sandström, He and Minta2010) also conducted similar experiments and arrived at the same results. Lee, Ahn & Kim (Reference Lee, Ahn and Kim2021), by studying the influence of density ratio (varying from 0.0012 to 0.0057) in a two-dimensional rectangular tank, observed great differences between the magnitude of the impact pressure when density ratios differ. Recent research on the effects of gas–liquid density ratios arrived at the same conclusion in that the first impact pressure clearly decreases as the density ratio increases, although the influence of the density ratio on the impact pressure needs to be investigated systematically and under different excitations, especially under irregular excitations. Ezeta, Kimmoun & Brosset (Reference Ezeta, Kimmoun and Brosset2023) experimentally studied wave impacts generated by solitary waves and the influence of the associated bathymetry profile as a function of ullage pressure. The density ratio between gas and liquid phases was adapted by using a mixture of nitrogen and water vapour for the gas phase. For ullage pressures exceeding 1005 mbar, they found a ‘backward inclination’ of impacting waves, which leads to smaller gas pockets upon impact compared with cases of lower ullage pressures. They also found that, for cases of higher ullage pressures, mean impact pressures at the wave crest are reduced when the density ratio is increased because impact velocities decrease.
Further research gradually focused on the effects of phase transition, due to its relevance in the transportation of fluids close to their boiling points, as well as its potentially significant influence on impact pressures. This critical phenomenon is one of the most crucial aspects of sloshing research. The basic mechanics of cavitation and boiling are similar. Brennen (Reference Brennen1995) described the precise distinction between cavitation and boiling with respect to temperature and vapour pressure. Olsen & Hysing (Reference Olsen and Hysing1974) presented their experimental results with water at 20, 95 and 100
$^\circ$
C. However, the influence of phase transition seems to be limited. On the contrary, Maillard & Brosset (Reference Maillard and Brosset2009) observed that the thermodynamic property of the gas significantly affects impact pressures, especially when phase transition occurs. Based on a simple one-dimensional piston model and experimental results, Braeunig et al. (Reference Braeunig, Brosset, Dias and Ghidaglia2010) demonstrated that, statistically, the impact pressures reduce and the pressure oscillations disappear due to phase transitions between the liquid and its vapour, which generate gas-pocket impacts. Ancellin et al. (Reference Ancellin, Ghidaglia and Brosset2012) included phase change effects in the generalised Bagnold model, in which the conservation of the total energy of the liquid–vapour system was considered to describe a thermodynamic relaxation process from unbalanced conditions to the liquid/gas equilibrium. Their model explained the mitigation of impact pressures and the complete damping of gas-pocket oscillations. Ancellin, Brosset & Ghidaglia (Reference Ancellin, Brosset and Ghidaglia2016) then studied numerically the effects of phase change on wave-impact loads. Their simulations were based on the one-dimensional Bagnold-type piston impact and phase change rates were modelled according to the Hertz–Knudsen relation. They compared their results with those obtained from Bagnold’s law. Grotle et al. (Reference Grotle, Halse, Pedersen and Li2016), by conducting experiments with liquefied nitrogen, confirmed the remarkable influence of the liquid’s temperature lower than the saturation temperature on the pressure. Kim, Lee & Kim (Reference Kim, Lee and Kim2017) conducted a series of drop tests to study the effects of phase transition on the impact pressure with gas pockets. Srikanth (Reference Srikanth2019) modelled the cryogenic propellant phase change processes for long duration spaceflights. His numerical investigations of phase change processes, for example, were modelled based on the Hertz–Knudsen–Schrage equation. He also studied the influence of the accommodation coefficient, which controls the phase change rates, on the motions of cryogenic propellants. Lee et al. (Reference Lee, Ahn and Kim2021) contributed to this topic by carrying out an experimental study using the combinations of water and air and the liquid and gas phases of Novec. They observed liquefaction occurring on the wall of the tank (local phase transition) inside the gas pocket during sloshing impact. During the sloshing impact, the absolute pressure exceeded the vapour pressure. This resulted in the gas condensing and a significant reduction of the impact pressure. However, Lee et al. (Reference Lee, Ahn and Kim2021) also demonstrated that the influence of phase transition on the decrement of the impact pressure was less dominant in their experiments than the influence of the density ratio. Based on these studies, the influence of phase transition needs to be investigated individually without the influence of viscosity and the density ratio, regardless of the varying fluid properties during phase transition. Besides, it is difficult to quantify these changes. Ezeta et al. (Reference Ezeta, Palacios Muñiz, Fan, Kim, Couty, Brosset and van der Meer2025) studied phase transition processes during the large-scale wave impact of a boiling liquid. The liquid was in thermal equilibrium with its own vapour analogue to the transport of cryogenic liquids. For breaking wave impacts onto a wall, they reported wave-impact pressures for boiling liquids that were up to two orders of magnitude larger than those in comparable water–air mixtures. The authors traced these pressures back to the full collapse of the entrapped vapour pocket.
Li et al. (Reference Li, Shen, Zhu, Ren and Zhang2024) investigated liquid sloshing characteristics in tanks with arbitrary excitation, numerically considering the energy equation as well as phase change processes. They applied arbitrary excitation forces to the tank such that liquid sloshing and oscillation forces could be directly simulated without mesh motion. They validated the solver with respect to phase change processes and showed that the evaporation rate of LNG increased under higher tank motion amplitudes. Bellur et al. (Reference Bellur, Médici, Hermanson, Choi and Allen2018) investigated temperature distributions at a solid–fluid interface during phase change processes of cryogenic sloshing both experimentally and numerically. The authors developed a thermal model of heat transfer of a cryostat environment which can accurately capture steady and transient temperatures. Fu et al. (Reference Fu, Sunden, Chen and Huang2015) developed source terms to simulate phase change processes involved in self-pressurised cryogenic tanks under microgravity. The numerical method treated the interface between the phases using a combination of a level set approach with the volume of fluid method. Their results showed that small vapour regions changed the pressure rise during sloshing. Zuo, Wu & Huang (Reference Zuo, Wu and Huang2021) evaluated the validity of popular liquid–vapor phase change models to simulate cryogenic self-pressurisation processes. Pressures and temperatures were computed based on the different phase change models and compared with experimental data. The authors found the coefficients of the Lee and Tanasawa models for phase change to significantly influence the prediction of liquid nitrogen self-pressurisation. The energy jump model failed to accurately predict pressures during the self-pressurisation process with reasonable accuracy. Himeno et al. (Reference Himeno, Ohashi, Anii, Daichi, Sakuma, Watanabe, Inoue, Umemura, Negishi and Nonaka2018) numerically and experimentally investigated phase change processes and the pressure drop during violent sloshing of cryogenic fluid. They studied the non-isothermal sloshing behaviour of liquid nitrogen and liquid hydrogen and found violent sloshing impacts to amplify phase change processes. To investigate gas entrapment and scale effects during water entry of a falling body, Lugni et al. (Reference Lugni, Wang, Faltinsen, Bardazzi, Lucarelli and Duan2021) performed experiments to locally measure pressures inside the gas pocket. They found that the natural frequencies of cavity oscillation and the associated pressures can be scaled based on the Froude and Euler laws of similarity.
Until now, measuring flow velocities in a moving tank has rarely been attempted. Rognebakke (Reference Rognebakke2002) carried out experiments to measure elevations of a liquid–gas interface at three places inside a rectangular tank using capacitive wave probes. Faltinsen & Timokha (Reference Faltinsen and Timokha2009) described the use of wave probes to measure the changing free water surface. These studies can only be considered as a first step in the recording of velocities. Although videos may provide information about the temporal behaviour of the free water surface, local velocities in the flow field cannot be determined. Gavory (Reference Gavory2005) performed experiments to measure flow velocities in an LNG model tank using a particle image velocimetry (PIV) system. The associated laser and a camera were mounted on a three degrees-of-freedom moving platform, and the harmonic tank motions and flow velocities were averaged over 100 motion periods. Although he rather roughly described the qualitative measurement of the flow field, he concluded that a PIV system may open new ground for future research to investigate details of impact phenomena, to validate numerical predictions and to test loads acting on installations in tanks. Lyu et al. (Reference Lyu, el Moctar, Potthoff and Neugebauer2017) also employed a PIV system to measure time-dependent flow velocities. However, a distinctive feature of their method was that flow velocity measurements were time resolved rather than averaged over one or more periods. Furthermore, their experimental set-up was designed to counteract the effects of camera vibrations, which are inevitable during measurements. Their procedure and the existing equipment were also suitable for stereoscopic measurements, where velocities perpendicular to the light-cutting plane were also captured. This allowed measurement of the three-component velocity vector field.
Based on the literature review, it can be stated that extensive experimental and numerical studies dealt with the influence of fluid properties and phase transition on sloshing and sloshing-induced impact pressures. Some of the findings are contradictory. In addition, there are various mathematical models to predict sloshing phenomena. These include both potential flow methods and methods based on solving the Navier–Stokes equations. To date, Bagnold models have mainly been used to investigate the phase transition effects. Only a few publications documented the use of mass transfer models. The existing methods often model only parts of the physical phenomena relevant for sloshing.
Here, we present a numerical approach based on the conservation of mass, momentum and energy that takes into account compressibility, thermal effects, phase transition, viscous effects and surface tension. To model the multiphase flow, we solved a transport equation for the volume fraction (volume of fluid), which we extended to account for the compressibility of the liquid and gas phases. We investigated the effects of flow properties and phase transitions for single-impact wave with gas pocket. However, our numerical approach is not restricted to these kinds of impacts. To validate our numerical method, we relied on experiments performed in our group using water, air, SF6 and Novec 649. Section 2 describes our numerical method, used to simulate sloshing flows. Section 3 provides a description of the sloshing test cases, the test scenarios and the experimental set-up. Section 4 introduces our numerical set-up. Section 5 presents an extensive verification and validation study. Section 6 describes the results of our investigation into the effects of flow compressibility, viscosity, density ratio and gas density and phase transition on sloshing-induced pressures. Finally, § 7 summarises our main findings.
2. Numerical method
We used a homogeneous mixture approach that considers the liquid and gas phases as continua. It simulates the free-surface flow on a Eulerian grid. The properties of the effective fluid (mixture) consisting of liquid and gas phases vary according to the volume fraction of each phase in the mixture. The gas volume fraction is defined as
$\alpha _g=V_g/V$
with the volume of gas,
$V_g$
, in a control volume,
$V$
. For a two-phase flow, the liquid volume fraction is defined as
$\alpha _l=1-\alpha _g$
. Using the volume fraction variables, the density,
$\rho$
, and viscosity,
$\mu$
, of the mixture are calculated as follows:
\begin{equation} \rho = \alpha _g \rho _g + (1-\alpha _g) \rho _l,\quad \mu = \alpha _g \mu _g + (1-\alpha _g) \mu _l. \end{equation}
Index ‘g’ denotes properties of the gas phase; index ‘l’, properties of the liquid phase. Properties of the mixture are then used to calculate the flow properties for the mixture. The continuous flow of the homogeneous mixture is calculated based on satisfying the mass conservation equation and solving the Navier–Stokes equations for a Newtonian fluid. The mass conservation equation yields
\begin{equation} \frac {\partial \rho }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u} ) =0, \end{equation}
where
$t$
is time and
$\boldsymbol{u}$
is the flow velocity of the mixture. The momentum equation reads
\begin{align} \frac {\partial (\rho \boldsymbol{u})}{\partial t} + \boldsymbol {\nabla }\boldsymbol{\cdot }( \rho \boldsymbol{u} \boldsymbol{u} ) = \boldsymbol{\nabla }\boldsymbol{\cdot }\left( \mu \left( \boldsymbol{\nabla }\boldsymbol{u} + \left( \boldsymbol{\nabla }\boldsymbol{u}\right)^{\mathrm{T}}\right) \right) - \boldsymbol{\nabla }\left (p + \mu \frac 23 \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} \right ) + \rho \boldsymbol{g} + \sigma \kappa \boldsymbol{\nabla }\alpha _l , \end{align}
where
$p$
is the pressure. The left-hand side comprises the time derivative of momentum and the convection term. The first two terms on the right-hand side represent diffusion and pressure, respectively. Symbol
$\boldsymbol{g}$
is the vector of gravitational acceleration. The last term on the right-hand side accounts for forces owing to surface tension,
$\sigma$
, between liquid and gas phases, where
$\kappa$
is the curvature of the interface. As surface tension acts directly at the interface between the two phases, it cannot explicitly be applied in interface capturing methods, where the interface may be smeared over multiple control volumes. We used the continuum surface model of Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992) to calculate the surface tension force as a volume force acting on the control volume. The curvature of the interface is defined as follows:
\begin{equation} \kappa = - \boldsymbol{\nabla }\boldsymbol{\cdot }\frac {\boldsymbol {\nabla }\alpha _l}{|\boldsymbol {\nabla }\alpha _l|}, \end{equation}
where the gradient of the volume fraction of the liquid phase,
$\alpha _l$
, points into the normal direction of the interface. Surface tension vanishes in control volumes occupying only one phase (
$\boldsymbol{\nabla }\alpha _l=\boldsymbol{0}$
).
The volume fraction of the liquid phase,
$\alpha _l$
, is obtained by solving the following transport equation considering the phase transition and the compressibility of the liquid and gas phases (Yu et al. Reference Yu, Goldsworthy, Brandner and Garaniya2017; Yuan, Song & Liu Reference Yuan, Song and Liu2019):
\begin{align} \frac {\partial \alpha _l}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha _l \boldsymbol{u}) - \alpha _l ( \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}) + \boldsymbol{\nabla }\boldsymbol{\cdot }(\alpha _l(1-\alpha _l) \boldsymbol{u}_r) = \nonumber \\ \alpha _l(1-\alpha _l)\left (\frac {1}{\rho _g} \frac {\mathrm{d}\rho _g}{\mathrm{d}t} - \frac {1}{\rho _l} \frac {\mathrm{d}\rho _l}{\mathrm{d}t}\right ) + (\dot {m}_v + \dot {m}_c) \left [\frac {1}{\rho _l} - \alpha _l\left (\frac {1}{\rho _l} - \frac {1}{\rho _g}\right )\right ]\! , \end{align}
where
$\dot {m}_c$
and
$\dot {m}_v$
are the mass transfer rates per unit volume accounting for condensation and vaporisation, respectively. The relative velocity between the liquid and gas phases,
$\boldsymbol{u}_r$
is modelled according to Weller (Reference Weller2008) as follows:
\begin{equation} \boldsymbol{u}_r = c_\alpha |\boldsymbol{u}| \frac { \boldsymbol{\nabla }\alpha _l}{| \boldsymbol{\nabla }\alpha _l|} , \end{equation}
where
$c_\alpha$
is a coefficient to control the artificial compression of the interface. The mass transfer rates per unit volume,
$\dot {m}_c$
and
$\dot {m}_v$
, that account for phase change were determined via the approach of Phan et al. (Reference Phan, Nguyen, Duy, Kim and Park2022). This approach is based on the theory of critical phenomenon dynamics, as explained by Hohenberg & Halperin (Reference Hohenberg and Halperin1977) and Lee (Reference Lee1980). The theory assumes that the mass transfer rates depend on the difference between value of a flow quantity and its critical value. While Lee (Reference Lee1980) used the difference between the actual and critical temperatures at which the phase transition begins to calculate non-equilibrium phase change processes in a two-phase mixture, Phan et al. (Reference Phan, Nguyen, Duy, Kim and Park2022) used the difference between the absolute and vapour pressures. This approach is analogous to the Hertz–Knudsen equation and we used it to derive the source terms for condensation and vaporisation as follows:
\begin{eqnarray} \dot {m}_c = c_c \alpha _l (1 - \alpha _l) \rho _l \frac {\max ((p-p_v),0)}{p_v}, \end{eqnarray}
\begin{eqnarray} \dot {m}_v = c_v \alpha _l (1 + \alpha _{{nuc}} - \alpha _l) \rho _l \frac {\min ((p-p_v),0)}{p_v} , \end{eqnarray}
where
$p_v$
is the vapour pressure,
$c_c$
is the coefficient for condensation and
$c_v$
is the coefficient for vaporisation. Here,
$\alpha _{{nuc}}$
is the volume fraction of the non-condensable gas nuclei and enables the vaporisation of the liquid phase in the absence of a larger liquid–gas interface. Based on previous studies, a non-condensable gas content of
$\alpha _{{nuc}}$
=
$1\times 10^{-6}$
was assumed.
The system of equations contains the momentum equation (2.3), which is solved to obtain a prediction for the velocity. To close this system of equations, an additional equation for the pressure is derived by forming the divergence of the momentum equation and combining it with the mass conservation equation (Jasak Reference Jasak1996). A derivation of the pressure equation including phase transition can be found in Peters & el Moctar (Reference Peters and el Moctar2020). The pressure equation considering phase transition and the compressibility of the liquid and gas phases reads as follows:
\begin{equation} \left ( \frac {\alpha _l \psi _l}{\rho _l} \frac {\partial p}{\partial t} + \frac {\alpha _g\psi _g}{\rho _g} \frac {\partial p}{\partial t}\right ) - \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \frac {1}{a_p} \boldsymbol{\nabla }p \right ) = (\dot {m}_{v} + \dot {m}_{c})\left ( \frac {1}{\rho _l} - \frac {1}{\rho _g} \right ) - \boldsymbol{\nabla }\boldsymbol{\cdot }\left ( \frac {\boldsymbol{H(u)}}{a_p} \right )\! . \end{equation}
The coefficient of compressibility of the liquid phase,
$\psi _l$
, is calculated from the Tait equation of state; the coefficient for the gas phase,
$\psi _g$
, from the ideal gas law. Here,
$\boldsymbol{H(u)}$
is a matrix which contains the off-diagonal part of the coefficient matrix and all other source terms of the discretised momentum equation without the pressure gradient and
$a_p$
is the diagonal coefficient of the coefficient matrix of the discretised momentum equation (Jasak Reference Jasak1996; Peters & el Moctar Reference Peters and el Moctar2020).
We used the PIMPLE algorithm (Weller Reference Weller2005), which is a combination of the SIMPLE (semi-implicit method for pressure linked equations) (Caretto et al. Reference Caretto, Gosman, Patankar and Spalding1973) and the PISO (pressure-implicit with splitting operators) algorithms (Issa Reference Issa1986) to sequentially solve the system of equations in an iterative manner. The main steps of the algorithm are as follows: first, a prediction of the velocity is obtained by solving the momentum equation (2.3); this velocity field does not yet fulfil the mass conservation equation. Second, the pressure equation is solved to update the pressure field. Third, the velocity field is corrected using the new pressure. This process is repeated iteratively until the residuals for the flow quantities decrease below pre-defined thresholds.
To compute temperature changes in sloshing flows, the energy equation needs to be solved for the temperature,
$T$
. It reads as follows:
\begin{align} \frac {\partial \rho T}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u} T) - \boldsymbol{\nabla }\boldsymbol{\cdot }(\beta \boldsymbol{\nabla }T ) = \nonumber \\ \left ( \frac {\alpha _l}{C_{v,l}} + \frac {1-\alpha _l}{C_{v,g}} \right ) \left ( (\dot {m}_c + \dot {m}_v) H_v -\frac {\partial \rho k}{\partial t} - \boldsymbol{\nabla }\boldsymbol{\cdot }(\rho \boldsymbol{u} k ) - \boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u}p) \right )\! . \end{align}
Here,
$k=0.5|\boldsymbol{u}|^2$
and, when multiplied by the density, yields the kinetic energy,
$\beta = \alpha _g \beta _g + \alpha _l \beta _l$
is the thermal diffusivity of the homogeneous mixture and
$\beta _g$
and
$\beta _l$
are the thermal diffusivities of the gas phase and the liquid phase, respectively. Also,
$C_{v,g}$
is the specific isochoric heat capacity of the gas phase;
$C_{v,l}$
, of the liquid phase. On the right-hand side of (2.10), the term containing the latent heat of vaporisation,
$H_v$
, accounts for temperature changes due to phase change (Phan et al. Reference Phan, Nguyen, Duy, Kim and Park2022). This means that a certain amount of heat must be transported away from the liquid–gas interface to cause a phase change from gas to liquid (condensation), or transported to the interface to cause a phase change from liquid to gas (vaporisation). While condensation leads to an increase in the temperature of the surrounding fluid, vaporisation leads to a decrease in temperature.
The above governing equations were discretised using the finite volume method. The second-order approximation mid-point rule, central differencing scheme and linear interpolation scheme were used to approximate the area and volume integrals, differentials and to interpolate the values at the control volume (CV) faces using values at the CV centres, respectively. The gradients in the diffusive fluxes were approximated using Gauss’ algorithm. We used the implicit Euler method for time integration. The above described mathematical models were implemented in OpenFOAM-11. More details can be found in Greenshields (Reference Greenshields2023).
3. Test case description and experimental set-up
Extensive experiments were conducted in our group (Neugebauer, el Moctar & Potthoff Reference Neugebauer, el Moctar and Potthoff2014) using the international sloshing benchmark tank of Loysel et al. (Reference Loysel, Chollet, Gervaise, Brosset and De Seze2012). The experimental results of these investigations will be published in a separate paper. We will therefore limit ourselves to relevant experimental results that we used to validate our numerical method. Our tank’s length, width and height were 946, 118 and 670 mm, respectively. To enable the use of videos to record optical measurements, the tank walls were constructed of Plexiglas. The solid walls of 75 mm thickness significantly reduced the deformation of the tank walls at reduced pressures inside the tank. Figure 1 shows a photograph of the test tank and the arrangement of the pressure sensors of type Kulite XTM190 with a diameter of 3.8 mm. The test tank was fixed to the hexapod, and this hexapod was then excited to simulate the tank motions. Here, tank motions are given in a right-handed Cartesian coordinate system, which is shown in figures 2 (a) and 7 (a). Pressure sensors were always read simultaneously at a frequency of 50 kHz. To measure liquid and gas temperatures, one sensor was located near the tank bottom at an end face and another senor near the centre of the tank ceiling. Close to the centre of the ceiling, another sensor was located to also monitor the global absolute pressure in the tank. Temperatures and pressures in the laboratory itself were also recorded. The liquid–gas interface was recorded using a Phantom v9.1 high-speed camera with a resolution of 1600
$\times$
1200 pixels.
Top view of sloshing tank of length
$L_T=946$
mm, width
$B_T=118$
mm, height
$H_T=670$
mm and sensor arrangement, unit is mm. Filling level = 85 %.
Time-resolved PIV measurements were carried out. The high-speed camera captured an area of 200 mm × 200 mm, see figure 2. The Phantom v9.1 high-speed camera was focused on the centre plane of the tank to reduce wall effects. The field of view comprises approximately 1000
$\times$
1000 pixels, which corresponds to a spatial resolution of 0.2mm
$\times$
0.2 mm. A frame rate of 0.1 kHz was used for the PIV measurements. The minimal separation time between the two images of a double frame maybe as low as 10
$\mathrm{\unicode{x03BC} }$
s. As particles, polyamide spheres with a diameter of 100
$\mathrm{\unicode{x03BC} }$
m were used. The PIV calculation was done by a cross-correlation algorithm for an interrogation window size of 32
$\times$
32 pixels. The flow velocity was measured at four positions, see figure 2. All experimental test runs were repeated several times.
Side view of the tank with field of view area with positions of the velocity measurement points A (
$x=312$
mm,
$y=32$
mm), B (
$x=376$
mm,
$y=32$
mm), C (
$x=376$
mm,
$y=148$
mm) and D (
$x=466$
mm,
$y=148$
mm). Sketch adapted from Lyu et al. (Reference Lyu, el Moctar, Potthoff and Neugebauer2017).
To ensure a reliable repeatability and an accurate evaluation of individual impacts, tank motions for the single-impact wave (SIW) were specified to comply with the international sloshing benchmark tests. The repeatability of the characteristics of the time histories of the sloshing-induced pressure and the time instances of first impact pressure was favourable. The tests started with the liquid at rest and, thus, they were independent of previous events. Although these motions and the associated measurements lasted only a few seconds, the high-speed video recordings captured the responses during the entire experiment. Thus, it was possible to obtain detailed time series of the results. A wave characterised by an inclined profile contacting the tank wall and forming a large gas pocket was examined. The tank motions were generated by superimposing a harmonic function with a tangent hyperbolic function. This function acted as a ramp function that delayed the envelope of these motions’ amplitudes attaining their maxima. Equation (3.1) describes this tank motion. Tank motions were generated over a two-period time interval (Lyu et al. Reference Lyu, el Moctar, Potthoff and Neugebauer2017). Figure 3 plots a normalised time history of the tank’s motion.
\begin{equation} x(t) = \begin{cases} A \tanh \left ( \pi \frac {t^2}{T_{{surge}}} \right ) \sin \left ( 2\pi \frac {t}{T_{\textit{surge}}} \right )\!, \quad & t \leqslant T_{\textit{surge}}, \\[5pt] A \tanh \left ( \pi \frac {(2T_{\textit{surge}} - t)^2}{T_{\textit{surge}}} \right ) \sin \left ( 2\pi \frac {t}{T_{\textit{surge}}} \right )\!, \quad & t \gt T_{\textit{surge}} . \end{cases} \end{equation}
Here,
$x$
represents the surge motion of the tank,
$t$
is time,
$T_{\textit{surge}}=0.975$
s is the period of motion and
$A=32$
mm is its surge amplitude.
Time history of normalised tank motion see (3.1).
To investigate the effects of the density ratio between the liquid and gas phases on the impact pressure acting on the tank walls, the composition of the gas phase was varied. Figure 4 shows the test set-up. A vacuum pump and a gas tank containing SF6 were connected via valves and gas lines. The pump allowing the absolute pressure to be controlled was used to generate a target pressure in the tank. The air and SF6 were then mixed in a defined ratio. Figure 5 shows the arrangement of pressure sensor arrays at the top of the sloshing tank.
Experimental set-up for density ratio investigations.
Arrangement of pressure sensor arrays.
Safety concerns did not allow the use of deep-cold gases to investigate the effects of phase transition on sloshing-induced impact pressures. Besides, the cold temperatures would have affected the installed measuring sensors. Furthermore, boiling the water would have significantly reduced the density ratio, rather than heating the methane and hydrogen and, also, it would have had an undesirable influence on the measuring equipment and on the test tank itself. Therefore, for these tests we chose a fluid that boils already at slightly elevated temperatures and moderately reduced pressures, namely, 3 M Novec 649. This meant, however, that temperatures and pressures had to be regulated during these test runs. Figure 6 shows the experimental set-up. The tank was filled with liquid from a higher reservoir via the lower opening. The fluid was heated by increasing the room temperature in the laboratory, as this was sufficient to raise the temperature of the fluid inside the tank. The heat transfer into the tank, therefore, took place through the tank walls. The evaporating fluid was discharged via valves located at the tank’s ceiling. The gaseous fluid was liquefied in a condenser and fed back into the reservoir. The vacuum pump and the throttle valve, installed to obtain the required system pressure, were located behind the capacitor and, therefore, reacted exclusively with the condensate. Thus, air flowing through the throttle valve did not influence the composition of the gas phase in the tank. This was important to ensure that the desired density ratio was maintained.
Experimental set-up for phase transition investigations.
3.1. Test matrix
All tests were performed for a filling level of 85 %. To investigate the effects of the density ratio between liquid and gas phases on the sloshing-induced impact pressures, water was used for the liquid phase. When changing the density ratio, the old gas mixture was moved out of the tank by filling the tank with water. Subsequently, the tank was filled with either ambient air or SF6. Then the internal tank pressure was decreased by means of the vacuum pump and the tank pressure was adjusted by adding the second fluid until the desired mixing ratio at ambient pressure was obtained. Using air and SF6 made it possible to vary the gas density between 1.25 and 6.63 kg m−
$^3$
. Table 1 lists the density ratios we considered, comprising the ratio for gas and liquid methane at
$-162\,^\circ$
C, the ratio of air and water, the ratio for three gas mixtures and water and the ratio for SF6 and water. With air as the lightest gas, we obtained the smallest density ratio; with SF6 as the heaviest gas, the largest density ratio. For each of these combinations, the experimental runs were repeated 25 times to enable a statistical analysis of the measured data.
Investigated fluids and corresponding DRs.
We investigated the influence of phase transitions, i.e. the vaporisation and condensation, on the impact pressure acting on tank walls. For these runs, a liquid that vaporised at reasonable combinations of temperature and pressure and that attained the density ratio of gas and liquid NG at
$-162\,^\circ$
C and hydrogen at
$-253\,^\circ$
C was used. The fluid we selected was the coolant liquid Novec 649. Besides, this fluid is sufficiently safe and was, therefore, suitable for use with the experimental investigations. Table 2 lists the associated thermodynamic properties, comprising saturation pressure, gas and liquid densities, density ratio and kinematic viscosity of the liquid phase. To obtain a pure two-phase mixture in the tank and to cause a phase transition by small local pressure changes, we specified two boiling point temperatures of 20 and 30
$^\circ$
C.
Thermodynamic properties of saturated liquid and gas phases of Novec 649.
The circuitry installed to adjust the density ratio, described above, included a buffer system that was developed for the fluids to comply with the unsteady processes during the experiments. At the beginning of a test run, the tank was filled to its specified level (85 %). Two minutes before the tank started to move, the pressure in the system was lowered to the boiling point of the coolant liquid by means of a vacuum pump and held at this temperature by actuating a mechanical air leakage valve. Before each experiment, the heavier gas of the refrigerant displaced the air remaining in the tank. Over the time span of reduced pressure, the coolant liquid continued to disperse, condensed in the capacitor and was then stored in a reservoir. After each experiment, the ambient pressure was restored to prevent further vaporisation. To prepare for the next experiment, the coolant liquid was pumped back into the tank.
Table 2 lists the fluid properties at the investigated boiling points. Table 3 summarises the test scenarios, comprising the three combinations of fluid temperature and absolute internal tank pressure: 20
$^\circ$
C and 1.0 bar, 30
$^\circ$
C and 0.5 bar and 20
$^\circ$
C and 0.3 bar. At the ambient temperature and pressure, the liquid phase did not vaporise. By comparing the vaporisation state with the non-vaporisation state at a similar density ratio, we examined the influence of phase transition on the impact pressure by excluding the effects of the density ratio. These runs were also repeated 25 times.
Test scenarios for the phase transition investigations using Novec 649.
4. Numerical set-up
We employed the numerical method described in § 2. Figure 7 (a) shows the tank schematically with its principal dimensions and the computational domain of the fluid inside the tank. The tank translated horizontally according to motion equation (3.1) causing a single-impact wave (SIW). The tank’s motion generated a wave travelling between the tank’s sidewalls, which caused the formation of a gas pocket at the tank’s top. The inertia of the liquid flow towards the tank’s top wall caused the liquid to compress the gas pocket, thereby generating high impact pressures. This process of compressing a gas pocket by an approaching liquid mass is characterised as a Bagnold-type impact (Bagnold Reference Bagnold1939).
Computational domain (a) and numerical grid (b), top right corner marked by a red rectangle. The CV’s edge length in the refined top region was 1 mm.
This benchmark tank was very narrow. Therefore, we treated the flow in the tank as being quasi-two-dimensional (quasi-2-D) and, consequently, we performed 2-D simulations. We assumed the flow to be compressible. For simulations without phase transition, we considered the flow to be isothermal and neglected solving the energy equation (2.10). For simulations with phase transition, we solved the energy equation for the temperature and considered the effect of the latent heat of vaporisation during phase change processes. We specified no-slip velocity conditions for the side, top and bottom walls. Figure 7 (b) shows not only the numerical grid that idealised the computational domain inside the tank, but also the grid’s region (surrounded by a red rectangle) in the top right corner of the tank subject to sloshing impacts. As seen, the locally refined grid, with a CV edge length of 1 mm, at the liquid–gas interface extended up to the top of the tank. During the entire simulation time of a run, the liquid–gas interface was always located within the region of highest mesh resolution. We used an implicit Euler scheme for time integrations and second-order approximations for spatial discretisation. The average Courant number was kept below 0.03, and the time step was adapted accordingly.
The horizontally translating tank induced an SIW inside the tank that led to a single impact event at the roof of the tank. Figure 8 shows the associated liquid–gas interface shapes of the fluid, obtained at the four consecutive times of
$t$
=
$1.5$
,
$1.55$
,
$1.60$
and
$1.65$
s. The liquid–gas interface is seen moving upwards along the right tank wall with a small angle of inclination until it finally forms a gas pocket in the top right corner. As the gas pocket forms and is compressed by the inertia forces of the approaching liquid mass, the pressure around and inside of the gas pocket is increased, thereby generating a sloshing-induced impact in the top right corner of the tank. Without considering phase transition, we first performed computations for various density ratios between the liquid and gas phases listed in table 1. Then we considered the phase change processes of vaporisation and condensation using the mass transfer model described in § 2.
Motion of the liquid–gas interface leading to the impact in the top right corner of the tank.
5. Verification and validation
5.1. Discretisation study
To investigate the effect of time discretisation and the Courant number, we performed 2-D simulations of the sloshing flow for the tank filled with water and a mixture of SF6 and air for a density ratio of
$2.5\times 10^{-3}$
on the finest grid G3 using the four successively smaller time step sizes of
$\Delta t$
=
$7.00\times 10^{-5}$
,
$5.00\times 10^{-5}$
,
$3.54\times 10^{-5}$
and
$2.50\times 10^{-5}$
s. To comply with the experimental set-up described in § 3, a virtual sensor represented pressure sensor 3 from the experiments, located at the tank’s top as shown in figure 5 (b). Figure 9 plots the corresponding time series. As seen, the resulting pressure peaks between simulations using the largest and smallest time step sizes differed by only a small amount of approximately 1.2 %. Consequently, for consecutive oscillations of the enclosed gas pocket, numerical diffusion increased only slightly, leading to a stronger damping of pressure amplitudes at larger time step sizes.
Time histories of the pressure acting on sensor 3 for the tank filled with water and a mixture of SF6 and air for a density ratio of
$2.5\times 10^{-3}$
. Computations performed on the finest grid G3 with different time step sizes
$\Delta t$
.
Time histories of the pressure acting on sensor 3 for the tank filled with water and a mixture of SF6 and air for a density ratio of
$2.5\times 10^{-3}$
(without phase transition) computed on numerical grids G1, G2 and G3. CV denotes the control volume.
To quantify the discretisation errors, we performed 2-D simulations of the sloshing flow for the same case of the tank filled with water and a mixture of SF6 and air for a density ratio of
$2.5\times 10^{-3}$
on three consecutively refined numerical grids comprising
$4.7\times 10^{4}$
(G1),
$9.3\times 10^{4}$
(G2) and
$1.9\times 10^{5}$
(G3) CVs. We used a refinement ratio of
$\sqrt {2}$
between edge lengths of the consecutive grids in all directions and specified time step sizes for simulations on the different grids based on the time step study (see figure 9), such that the average Courant number was kept the same and did not exceed a value of 0.03. Taking into account spatial discretisation, for every numerical grid, the dimension of the virtual sensor was generated to be nearly equal to those of the physical sensors. We chose this sensor because it was mounted closest to the sidewall where the gas pocket was formed and because it best represented the impact area. Figure 10 plots time histories of the pressure acting on sensor 3, obtained on grids G1, G2 and G3. Table 4 lists the associated parameters identifying the three simulations performed on these three grids for the tank filled with water and air. As seen, the highest impact pressure was obtained on G3, i.e. on the grid with the highest spatial resolution. The differences in first impact pressures between consecutive grids became smaller with higher spatial resolution. The difference in the timing of the first impact pressure can be attributed to the variation in the shape of the liquid–gas interface prior to the closure of the gas pocket. The coarse grid failed to capture disturbances in the liquid–gas interface. As the grid resolution increased, the volume enclosed in the gas pocket decreased slightly and the frequency of pressure oscillations increased. This is generally consistent with the Minnaert formula (see (5.6)).
Results of the discretisation study.
We estimated the discretisation error
$\epsilon$
and the discretisation independent solution
$\varphi _\infty$
(corresponding to a solution obtained on an infinitesimally small spatial and temporal discretisation) using the constant Courant number approach of Oberhagemann (Reference Oberhagemann2017) and el Moctar, Schellin & Söding (Reference el Moctar, Schellin and Söding2021). Error
$\epsilon$
was estimated as follows:
\begin{equation} \epsilon \approx \delta = \varphi _i - \varphi _\infty = a_1 \varUpsilon _i + a_2 \varUpsilon _i^2 + \ldots, \end{equation}
where
$\varphi _i$
is the solution on grid
$i$
(e.g.
$i$
= 1, 2, 3 identify, respectively, the coarse, medium and fine grids),
$\varphi _\infty$
is the grid independent solution and
$a_1$
and
$a_2$
are constant coefficients. The non-dimensional refinement ratio
$\varUpsilon _i$
then reads as follows:
\begin{equation} \varUpsilon _i = \sqrt {0.25\left [ \left ( \frac {1}{r_x}\right )^{2(i-1)} + \left ( \frac {1}{r_y}\right )^{2(i-1)} + \left ( \frac {1}{r_z}\right )^{2(i-1)} + \left ( \frac {1}{r_t}\right )^{2(i-1)} \right ]}, \end{equation}
where
$r_x$
,
$r_y$
,
$r_z$
and
$r_t$
are the refinement factors in the
$x$
,
$y$
,
$z$
directions and at time
$t$
, respectively. In this study, the refinement factors were constant
\begin{equation} r_x=r_y=r_z=r_t=\sqrt {2}. \end{equation}
We used the least-squares formulation to minimise
\begin{equation} S(\varphi _\infty ,a_1) = \sqrt {\sum _{i=1}^n \left(\varphi _i - \left(\varphi _\infty + a_1 \varUpsilon _i^2\right)\right)^2}, \end{equation}
which assumes a second-order polynomial formulation of the discretisation error
$\epsilon$
. Here,
$n$
is the number of refinement levels. In this study
$n$
= 3.
Here, we define
$\varphi _1$
,
$\varphi _2$
and
$\varphi _3$
as the first impact pressure obtained from simulations on grids G1, G2 and G3. We then use these values to minimise (5.4) and to obtain an estimate for the grid independent solution, denoted
$\varphi _\infty$
.
To evaluate the convergence of the solution, we defined the convergence ratio
$R$
as
\begin{equation} R = \frac {\varphi _3 - \varphi _2}{\varphi _2 - \varphi _1}. \end{equation}
Monotonic convergence is obtained when
$0\lt R\lt 1$
, and divergence when
$R\gt 1$
.
To demonstrate the effect of the grid refinement ratio on pressure peaks, figure 11 plots the pressure peaks against the grid refinement ratio obtained from our simulations. The three red squares identify pressure peaks obtained on the three grids, with the grid refinement ratio decreasing from 1.0 for the coarsest grid to 0.5 for the finest grid. The black dot marks the extrapolated grid-independent pressure peak, which of course corresponds to a refinement ratio of zero, and the black curve represents the associated second-order polynomial. Table 4 lists the results from the discretisation study. Compared with the extrapolated grid independent peak pressure,
$\varphi _\infty$
, the peak pressure computed on the finest grid G3,
$\varphi _3$
, deviated by 3.8 %; on the medium grid G2,
$\varphi _2$
, by 9.4 %; and on the coarsest grid G1,
$\varphi _1$
, by 17.02 %. The numerical solution converged monotonically with a convergence ratio of
$R=0.725$
.
Extrapolated discretisation-independent sloshing-induced pressure peak vs. non-dimensional refinement ratio
$\varUpsilon$
for the tank filled with water and a mixture of SF6 and air for a density ratio of
$2.5\times 10^{-3}$
.
5.2. Comparison of 2-D and 3-D simulations
As described in § 4, we assumed the flow inside the tank to be two-dimensional because of its very small depth. To investigate 3-D effects, we performed computations of the 3-D flow inside the 3-D tank corresponding to the physical tests described in § 3. To ensure that the discretisations in the 2-D and 3-D tanks were the same, we extruded the grid in the third dimension. The edge length of the CVs in the third dimension was kept almost equal to the edge length in the other two dimensions, resulting in a total number of
$5.8\times 10^{6}$
almost cube-shaped CVs for the 3-D grid. Figure 12 plots comparative time histories of the pressure acting on sensor 3 mounted at the roof of the tank filled with water and air. As seen, the maximum pressures at the sensor during the initial gas pocket compression compared favourably. For the 2-D case, the higher amplitudes of pressure minima and maxima deviated only slightly from the 3-D case, and their frequencies corresponding to the oscillating gas pocket also agreed well with each other. The natural frequency of a gas bubble oscillating in an infinite liquid can be calculated according to Minnaert (Reference Minnaert1933) as
\begin{equation} f_M = \frac {1}{2\pi a} \sqrt {\frac {3 \gamma p_\infty }{\rho _l}} \, . \end{equation}
Here,
$a$
is the equilibrium radius of the bubble,
$p_\infty$
is the pressure in the liquid at infinity (1 bar) and
$\gamma$
= 1.4 is the exponent of adiabatic compression of air. We calculated
$a$
= 5.92 cm as a representative radius from the volume of the gas pocket (see table 5). This yields a Minnaert frequency of
$f_M$
=
$55.14$
Hz. From the pressure time history of our 3-D simulation (see figure 12), we obtained an oscillation frequency of
$f_{\textit{osc}}$
= 46.72 Hz. This value agrees fairly well with the value obtained from Minnaert’s equation (5.6).
Comparison of the time history of the pressure acting on sensor 3 computed on the 2-D grid G2 and 3-D grid consisting of
$5.8\times 10^{6}$
CVs for the tank filled with water and air.
Peak pressure, enclosed gas volume and mass at impact, oscillation frequency of pressure, vertical distance to sensor at impact and vertical velocity of liquid–gas interface before impact for different DRs.
5.3. Validation
To validate our numerical results, we first compared the time history of the computed and measured vertical velocity and sloshing-induced pressure acting on the roof of the tank filled with water and air. All computations were performed on grid G2. Figure 13 plots the measured and computed time histories of vertical velocity at the four sensor positions described in figure 2. The solid line identifies the vertical velocity from the computations; the squares, the vertical velocity from the measurements. As seen, computed values agree well with experiments, showing that the numerical method predicted the flow velocity in the field of view area at each time step fairly well. The differences between the computed and measured high-frequency oscillations after the impact were mostly due to the lower temporal resolution of the PIV measurements compared with the numerical simulations. Positions A and B were always wetted before the time instant 1.5 s, which enabled particle velocities to be monitored during this time of the test run. As sensor positions C and D were positioned higher up, they are wetted only during part of the time of the test run.
Computed (black line) and PIV-measured (red squares) time histories of the vertical velocity
$U_z$
at observed positions A, B, C and D (see figure 2) for the tank filled with water and air. Computations were performed on grid G2.
Depending on the type of impact, the measured values of impact loads may be scattered. This is particularly true for a flat impact, (Tödter et al. Reference Tödter, el Moctar, Neugebauer and Schellin2020). The scenario considered here is a gas pocket impact, for which the repeatability is good. Figure 14 plots the time history of the mean (black solid line) and standard deviation (purple) of the pressure acting on sensor 3 obtained from multiple repetitions of the experiments.
Time history of the mean (black solid line) and standard deviation (purple) of the pressure acting on sensor 3 (see figure 5) obtained from multiple repetitions of the experiments for the tank filled with water and air. The figure on the right shows an enlarged view of the figure on the left.
Figure 15 plots time histories of the pressure acting on sensor 3 as a mean value obtained in the experiment (red line) and as computed in the simulation (black line). The computed impact pressure peak and the time history of pressure during the initial compression of the gas pocket agree well with the measured values. After the first compression, differences in simulation and measured results become more apparent. After the impact, more gas bubbles are trapped in the liquid phase, and the mixing of the two phases is more pronounced. This changes the compressibility of the ventilated liquid phase. The homogeneous mixture approach may not have been able to accurately predict the ventilated flow. An Euler–Lagrange approach would have been a better choice, (Peters & el Moctar Reference Peters and el Moctar2020). Furthermore, the 3-D effects are more pronounced after the impact, and the pressure amplitudes in the 3-D flow are smaller than in the 2-D flow, see figure 12.
Time histories of the computed pressure (black line) and measured mean pressure (red line) acting on sensor 3. The tank was filled with water and air.
To validate our numerical results for different density ratios, we performed computations of the sloshing flow inside the tank filled with water and different gases, see table 1. Figure 16 plots the time histories of the computed (solid lines) and measured (squares) pressure acting on sensor 3 for density ratios (DRs) of
$2.5\times 10^{-3}$
and
$5.3\times 10^{-3}$
. For both DRs, the computed and measured first impact pressures agreed well. For the second pressure oscillation, computed and measured pressure amplitudes and frequencies agreed well for the DR of
$5.3\times 10^{-3}$
, whereas only a fair agreement was reached for the DR of
$2.5\times 10^{-3}$
. After the third oscillation, deviations between measured and computed frequency as well as pressure amplitude became more apparent. Note that the flows became chaotic after the first impact as more gas bubbles were trapped in the liquid.
Computed (solid line) and measured (squares) time histories of the pressure acting on sensor 3 for different DRs. Computations performed on the finest grid G3.
In addition to the comparison between our numerical and experimental results, we compared our computed impact pressure peaks with theoretical values. Bagnold (Reference Bagnold1939) derived an equation to describe the pressure on a 1-D wall induced by a compressed gas pocket enclosed between a liquid and the wall
\begin{equation} \left ( \frac {p_{\textit{max}}}{p_0} \right )^{\frac {\gamma }{\gamma - 1}} + (\gamma - 1) \left ( \frac {p_{\textit{max}}}{p_0} \right )^{-1/\gamma } - \gamma = \frac 12 \frac {\gamma - 1}{p_0} \rho _l u_0^2 (\alpha -1) . \end{equation}
Here,
$p_{\textit{max}}$
is the maximum impact pressure at the wall,
$p_0$
is the undisturbed pressure in the far field,
$u_0$
is the initial velocity of the liquid mass towards the wall and
$\alpha$
is the ratio of the length of the liquid mass compared with the length of the enclosed gas pocket. Here,
$\gamma =C_p/C_v$
is the exponent of adiabatic compression, calculated from the ratio of the isobaric heat capacity
$C_p$
and the isochoric heat capacity
$C_v$
. Bredmose & Bullock (Reference Bredmose and Bullock2008) extended this pressure law to 2-D cases with the following equations:
\begin{equation} \left ( \frac {p_{\textit{max}}}{p_0} \right )^{\frac {\gamma }{\gamma - 1}} + (\gamma - 1) \left ( \frac {p_{\textit{max}}}{p_0} \right )^{-1/\gamma } - \gamma = \frac {\gamma - 1}{p_0} \rho _l u_0^2 \ln {\alpha } . \end{equation}
As described in § 4, we performed 2-D simulations of the two-phase sloshing flow for different DRs. However, the closure of the gas pocket occurs mostly due to the vertical motion of the liquid–gas interface at a small angle relative to the tank’s top wall (see figure 8 at times
$t$
= 1.60 and 1.65 s). Figure 17 plots the normalised impact pressure
$p^*=(p-p_0)/p_0$
versus the scale factor
$S=\rho _lu_0^2\alpha /p_0$
computed using the Bagnold equation for one dimension (5.7), for two dimensions (5.8) and for
$\gamma =1.4$
(air) and
$\gamma =1.1$
(SF6), in comparison with results obtained using our numerical method for different DRs. The initial velocity of the liquid mass,
$u_0$
, was taken from table 5. In the Bagnold equations,
$\alpha$
is the ratio of the length of the liquid mass to the length of the enclosed gas pocket. This was calculated using the vertical distance of the liquid–gas interface to the top wall (as given in table 5) and the height of the liquid column in the tank. For
$\gamma =1.4$
,
$\alpha$
was calculated based on the air–water case (DR =
$1.25\times 10^{-3}$
), and for
$\gamma =1.1$
, it was calculated based on the SF6–water case (DR =
$6.63\times 10^{-3}$
). Normalised impact pressures,
$p^*$
, agreed well with the pressure law of a 2-D Bagnold-type impact, as given in Bredmose & Bullock (Reference Bredmose and Bullock2008).
Comparison of the normalised impact pressures using Bagnold equations (1-D, 2-D) for
$\gamma =1.4$
(air) and
$\gamma =1.1$
(SF6) in comparison with results obtained using our numerical method for five different DRs (see cases 1–5 in table 1). For each case, the normalised impact pressure,
$p^*$
, the liquid density,
$\rho _l$
, the velocity of the approaching liquid column,
$u_0$
, and the ratio of the height of the liquid mass compared with the height of the enclosed gas pocket,
$\alpha$
, varied. The initial tank pressure
$p_0=1$
bar was constant for all cases.
We computed the sloshing flow at the boiling point 1, namely, at an initial temperature of 20
$^\circ$
C, an initial pressure of 0.3 bar in the tank and a vapour pressure of
$p_v$
= 0.326 bar (see table 3). We estimated the vaporisation and condensation coefficients based on previous experiments for Novec 649 to be
$c_c$
=
$c_v$
= 3.25 s
$^{-1}$
. Figure 18 plots time histories of the pressure acting on sensor 3, obtained from our experimental measurements (red line and red squares) and from computed pressures, taking into account the phase transition (black line). Taking into account the phase transition, we found that the first impact pressure and the subsequent pressure oscillations were dampened. Our computed pressure time history compared favourably with the experimental measurements regarding the dampening of both the first impact pressure and the subsequent oscillatory pressures.
For the tank containing Novec 649, time histories of absolute pressures acting on sensor 3, obtained for boiling point 1 (20
$^\circ$
C and 0.3 bar,
$p_v$
= 0.326 bar). The red line and red squares identify experimental measurements; the black line, computed results with phase transition (
$c_c$
=
$c_v$
= 3.25 s
$^{-1}$
).
In addition to the above scenario, we computed the sloshing flow at boiling point 2 (initial temperature of 30
$^\circ$
C, initial pressure of 0.5 bar and vapour pressure
$p_v$
= 0.496 bar), taking into account the phase transition. Figure 19 plots time histories of the pressure acting on sensor 3, obtained from experimental measurements (red line and red squares), and from our simulation (black line). Our numerical results generally compared favourably with the experimental results. As in the previous case, the first impact pressure and the subsequent pressure oscillations were damped.
A detailed analysis of the effects of phase transition on the sloshing-induced impact pressure and the subsequent pressure oscillations is presented in § 6.4.
For the tank containing Novec 649, time histories of absolute pressures acting on sensor 3, obtained for boiling point 2 (30
$^\circ$
C and 0.5 bar,
$p_v$
= 0.496 bar). The red line and red squares identify experimental measurements; the black line, computed results with phase transition (
$c_c$
=
$c_v$
= 3.25 s
$^{-s}$
).
6. Results
Without considering phase transition, we first performed computations for various viscosities and various DRs between the liquid and gas phases listed in table 1. Then we considered the phase change processes of vaporisation and condensation using the mathematical model described in § 2.
6.1. Effects of compressibility
To investigate the influence of the flow compressibility on the sloshing-induced impact pressure involving gas pockets, we computed the sloshing flow inside the tank filled with water and air with a DR of
$1.25\times 10^{-3}$
. First, we assumed the air to be incompressible. Second, we considered the air to be compressible relating pressure, density and temperature according to the ideal gas law. Figure 20 plots the time histories of the pressure acting on sensor 3 for the incompressible air (black line) and compressible air (red line). The incompressibility assumption of the air led to an overestimation of the impact pressure peak compared with compressible air which absorbs a part of the energy of the impact by compression. As expected, the compression and expansion of the gas pocket cannot be observed for the incompressible air. This observation agrees well with former work of Lyu et al. (Reference Lyu, el Moctar, Potthoff and Neugebauer2017). It should be noted that these effects are mainly valid for sloshing-induced impact pressures involving gas pockets or strongly mixed liquid and gas.
Time histories of pressure acting on sensor 3 for tank filled with water and incompressible air (black) and compressible air (red).
6.2. Effects of viscosity
To investigate the effects of the liquid viscosity we computed the sloshing flow inside the tank filled with water and air (see table 1, case 1), for a liquid with zero viscosity, and a liquid with five times the viscosity of liquid water. The remaining fluid properties were identical for all three liquids. Figure 21 plots the time histories of the pressure acting on sensor 3 for a tank filled with air and inviscid liquid water (black line), liquid water (red line) and liquid with five times the water viscosity (blue line). Peak pressures between the case of inviscid water and water with five times the viscosity deviated by less than 3 %. Viscosity slightly damped the sloshing impact pressure oscillations.
Time histories of pressure acting on sensor 3 for a tank filled with air and inviscid liquid water (black line), liquid water (red line) and liquid with five times the viscosity of water (blue line).
6.3. Effects of density ratio
Based on the experiments, we computed the sloshing flow for different DRs between the gas phase (a mixture of air and SF6) and the liquid phase (water) as described in table 1. Figure 22 plots the sloshing-induced impact for the case of air and water. As the liquid–gas interface approached the top right corner of the tank a gas pocket was enclosed and compressed consecutively by the approaching liquid mass. During the gas compression process, the high pressure at the liquid–gas interface was transported towards the top right corner of the tank where the sensor array was located. At approximately 1.614 s, the highest pressures in the gas pocket were observed. Afterwards, the gas expanded again so that the pressure in the enclosed pocket decreased.
Sloshing-induced impact generates high pressure in top right corner as the gas pocket is closed and compressed by the mass of the approaching liquid.
Figure 23 plots the time histories of the pressure acting on sensor 3 for different DRs. For lower DRs between the gas phase and the liquid phase, first impact pressures were higher, which agreed well with the experimental measurements given in figure 16. The higher the DR, the smaller the pressure peak; see figure 24. These findings were reported in the literature (Bogaert et al. Reference Bogaert, Léonard, Brosset and Kaminsk2010; Lee et al. Reference Lee, Ahn and Kim2021). Pressure peaks occurred earlier and pressure amplitudes in the following oscillations were higher for the lower DRs. The frequency of oscillations was higher for the lower DRs. All these tendencies were consistently observed for all DRs and agreed well with the experimental measurements.
Time histories of pressure for different DRs.
Impact pressure peaks for different DRs.
In our simulations, the change of density of our gases, e.g. pure air, pure SF6 or a mixture of these two gases, was calculated based on the ideal gas law. To better understand the effects of DR on impact pressures, we analysed the motion of the liquid–gas interface prior to the impact event. Figure 25 plots the liquid–gas interface position at the time instance of the first impact pressure on sensor 3 for simulations of different DRs. The liquid–gas interface was reconstructed as an iso-surface for a liquid volume fraction of
$\alpha _l$
= 0.5. Note that
$x$
- and
$y$
-axes are scaled differently for better visualisation. For lower DRs, the volume of the enclosed gas pocket was smaller than for higher DRs and the interface was closer to the pressure sensors located in the top right corner of the tank.
We believe that the main reason for this observation is related to the inertia of the gas. Shortly before the liquid–gas interface reached the top of the tank to close the gas pocket, some of the gas escaped horizontally from the pocket shaped as a high-velocity jet. The velocity of this jet was higher for the lighter gases. The shear flow across the liquid–gas interface generated vortices, indicating Kelvin–Helmholtz instabilities, which caused deformations of the interface. At lower densities, and thus at higher velocities, the disturbances of the liquid–gas interface were more pronounced. Note that the liquid phase was kept the same (water) for all DRs.
Liquid–gas interface positions, reconstructed for
$\alpha _l$
= 0.5, in top right corner of the tank at the time instance of the first impact pressure acting on sensor 3 for different DRs. Same coordinate system as shown in figure 2. Note that x- and y-axes are scaled differently.
Figure 26 plots the liquid–gas interface for the time instances of 1.3 and 1.6 s, for the lowest DR, DR =
$1.25\times 10^{-3}$
, of pure air and water (black line) and for the highest DR, DR =
$6.63\times 10^{-3}$
, of pure SF6 and water (red line). Initially, before the interface approached the sensor in the top right corner of the tank, the location of the interface was almost the same for both DRs. As the liquid–gas interface moved upwards along the right wall at an increasing velocity and began to enclose the gas volume in the top right corner, for the lower DR, the interface moved faster and approached the wall further. Note that the air density (lower DR) is more than five times smaller than the density of SF6 (higher DR). The inertia of the gas for a smaller density magnitude is, therefore, also much smaller, so that more gas can escape the top right corner during formation of the enclosed gas pocket. Figure 27 shows the liquid volume fraction with constant sized velocity vectors for a case of SF6 and water (DR =
$6.63\times 10^{-3}$
) after 1.6 s of simulation time. Disturbances appear at the liquid–gas interface in the left part of the tank, which are more significant for the higher DRs. Vortex formation was observed along the interface, indicating the generation of Kelvin–Helmholtz instabilities, which may have induced these disturbances.
Positions of the liquid–gas interface for two DRs at time instances of 1.3 and 1.6 s as the interface approaches the top right corner of the tank before the impact. Same coordinate system as shown in figure 2. Note that x- and y-axes are scaled differently.
Vortex formation in the gas phase along the liquid–gas interface causing weak disturbances of the shape of the interface DR
$=6.63\times10^{-3}$
.
Figure 28 plots the position of the liquid–gas interface approaching at four different time instances leading up to the impact. For both DRs, each of the four lines represents the location of the liquid–gas interface at the same time. In case of the lower DR, as the interface approached the tank roof, it accelerated, separated itself into two main parts and formed multiple gas pockets at the tank’s roof.
Higher first impact pressures were observed in cases of lower DRs because the gas pocket was smaller and the liquid–gas interface was closer to the pressure sensor when the gas pocket closed. In these cases, the pressure wave generated at the liquid–gas interface had to travel a shorter distance to reach the pressure sensor on the top wall of the tank. Furthermore, a larger gas volume can absorb more energy of the impact during compression than a smaller gas volume. Table 5 lists the enclosed gas volume, the enclosed gas mass at the time instance of first impact pressure, the oscillation frequency derived from the time history of pressure, the vertical distance of the interface at the right wall from the sensor at the time instance of maximum pressure and the vertical velocity of the interface,
$u_0$
, between 1.5 and 1.6 s of simulation time, i.e. the time just before the impact. Figure 29 plots the enclosed gas volume and vertical distance of the liquid–gas interface to pressure sensor 3 at the time instance of first impact pressure for different DRs. Figure 30 plots the oscillation frequency of the enclosed gas volume and the vertical velocity of the liquid–gas interface before impact for different DRs. For the larger DRs, the enclosed gas mass was larger and the oscillation frequency of the gas pocket decreased. The vertical distance between the liquid–gas interface and the sensor was smaller for smaller DRs, which was a result of the lower gas inertia and, consequently, the higher velocity of the interface prior to the impact. Both the enclosed gas volume and the correlated vertical distance of the interface to the sensor at the tank increased with increasing DR. With increasing DR, the oscillation frequency and the impact velocity of the interface consistently decreased.
Positions of the liquid–gas interface for two different DRs of
$1.25\times 10^{-3}$
(black) and
$6.63\times 10^{-3}$
(red) at four different time instances as the interface forms the gas pocket and generates the highest impact pressure acting at the right corner of the tank. Same coordinate system as shown in figure 2. Note that x- and y-axes are scaled differently.
Enclosed gas volume and vertical distance of the liquid–gas interface (FS) to pressure sensor 3 at the time instance of the highest pressure versus DR.
Oscillation frequency and vertical velocity of the liquid–gas interface (FS) before impact versus DR.
To rule out the possibility that the effects of DR on the impact pressure, as described above, are not significantly influenced by the compressibility of the gas phase, we computed the sloshing flow for the smallest (
$1.25\times 10^{-3}$
) and the largest (
$6.63\times 10^{-3}$
) DRs, assuming the gas phase to be incompressible. Figure 31 shows the time history of the pressure for both DRs without compressibility effects. Whereas in the compressible flow case the pressure fell by 47.9 % from 16.92 kPa for DR =
$1.25\times 10^{-3}$
to 8.82 kPa for DR = 6.63–3, in the incompressible flow case, it fell by 33.5 % from 22.59 kPa for DR =
$1.25\times 10^{-3}$
to 15.03 kPa for DR =
$6.63\times 10^{-3}$
. This shows that the effects of the DR on the impact pressure peak are slightly influenced by the compressibility of the gas phase. Another way to investigate the influence of DR on impact pressures without affecting the compressibility is to use the same components for liquid and gas phases and to change the initial absolute pressure and, thereby, to change the density of the gas phase under isothermal conditions. The gas adapts its density under these conditions according to the pressure change based on the ideal gas law. However, in our investigations, we specified the density of the gas phase by changing the composition of the gas and, thus, changing its molar mass.
Time histories of pressure acting on sensor 3 for a tank filled with water and incompressible air with a DR of
$1.25\times 10^{-3}$
(black) and a tank filled with water and incompressible SF6 with a DR of
$6.63\times 10^{-3}$
(red).
To investigate the relationship between the closure of the gas pocket and the first impact pressure at sensor 3 for the lowest DR of DR =
$1.25\times 10^{-3}$
, figure 32 plots the liquid–gas interface at four consecutive time instances, each 2.0 ms apart, before impact and at the time of impact,
$t_{imp}$
, when the computed first impact pressure is reached. At the first four time instances, the gas pocket was not yet fully closed. While the motion of the liquid–gas interface along the right tank wall almost stagnated, the interface approached the tank’s top at a distance of approximately 0.200 m from the sidewall to close the gas pocket. The maximum pressure was computed about 1.0 ms after closure of the gas pocket.
Positions of the liquid–gas interface for a water–air mixture at consecutive time instances leading to the closure of the gas pocket and the associated first impact pressure at sensor 3. Time
$t_{\textit{imp}}$
marks the time instance of maximum pressure at the sensor location; times lower than
$t_{\textit{imp}}$
indicate times just before the impact. Same coordinate system as shown in figure 2. Note that x- and y-axes are scaled differently.
To investigate a more practically relevant fluid, we computed the sloshing flow for a mixture of LNG and NG, with DR =
$4.0\times 10^{-3}$
, which is close to the DR of a mixture of air and SF6 and water having DR =
$3.9\times 10^{-3}$
. The density of NG used was assumed to be 1.874 kg m−
$^3$
, which is 50 % larger than the density of air. The density of LNG was assumed to be 450 kg m−
$^3$
, which is approximately 55 % smaller than the density of water. Figure 33 plots the time histories of the pressure acting on sensor 3 for the mixtures of air and water (black line), and for NG and LNG (red line). The first impact pressure for the air–water mixture is significantly higher than for the NG–LNG mixture. This agrees well with our assumption that the density of the gas phase, leading to more or less inertia of the gas, has a significant effect on the impact pressure, in contrast to the DR itself.
Time histories of pressures acting on sensor 3, for a water–air mixture (black line) with DR =
$1.25\times 10^{-3}$
and a LNG–NG mixture (red line) with DR =
$4.0\times 10^{-3}$
.
Normalised time histories of normalised pressures acting on tank roof (sensor 3) for simulations with the same DR between gas and liquid phases; black line represents the simulation for 0.2 times the densities of air and water, red line represents the simulation for half the densities of air and water, blue line represents the simulation for air and water, purple line represents the simulation for twice the densities of air and water and orange line represents the simulation for five times the densities of air and water.
Based on our results, we assume that, not the DR itself, but the density magnitude, has a significant effect on the formation of the gas pocket and consequently on the induced first impact pressure. To demonstrate this, we computed the sloshing flow for five gas–liquid mixtures with adjusted densities of the gas and liquid phases by keeping DR constant at DR =
$1.25\times 10^{-3}$
. Here,
$\rho _{\textit{ref},g}$
refers to the reference density of air and
$\rho _{\textit{ref},l}$
refers to the reference density of liquid water. In four cases, the initial densities of the gas phase
$\rho _{0,g}$
and the liquid phase
$\rho _{0,l}$
were specified by multiplying the reference gas density
$\rho _{\textit{ref},g}$
and the liquid density
$\rho _{\textit{ref},l}$
, respectively, by a factor less than or greater than unity. Figure 34 plots the normalised time histories of the normalised pressures acting on the tank’s roof (sensor 3) for simulations with adjusted gas and liquid densities by keeping the same DR. Here,
$t_{\textit{imp}}$
is the time instance of maximum pressure acting on sensor 3. The time scale used for normalisation is
$T=\sqrt {\rho _{0,l}/p_0} L_{\textit{ref}}$
with the initial tank pressure
$p_0$
= 1 bar and the reference length
$L_{\textit{ref}}$
= 670 mm, which was chosen according to the tank height. The vertical velocity of the liquid–gas interface,
$u_0$
, was obtained by averaging over a time from 1.5 to 1.6 s. Even though
$u_0$
decreased, for simulations with higher densities, the maximum normalised impact pressures were smaller. This aligns with our previous assumption that gases of higher densities, owing to their increased inertia, escape the gas pocket slower.
To validate the scaling of pressures, we compared our numerical simulations for previously investigated cases (see figures 23, 33 and 34) with Bagnold equations for one dimension (5.7), and for two dimensions (5.8). Figure 35 plots normalised impact pressures
$p^*$
using Bagnold equations for one dimension (black lines), for two dimensions (red lines) and for
$\gamma =1.4$
(air; solid lines) and
$\gamma =1.1$
(SF6; dashed lines) in comparison with results obtained using our numerical method. Numerical results are presented for a case of air and water (black diamond), cases of different DRs with gaseous mixtures of air and SF6 and liquid water (blue crosses), cases of constant DR of DR =
$1.25\times 10^{-3}$
and different initial liquid and gas densities (green squares) and a case of LNG (purple ‘x’). For all cases, the normalised impact pressures agreed fairly with the Bagnold 2-D equation given by Bredmose & Bullock (Reference Bredmose and Bullock2008).
Comparison of the normalised impact pressures using Bagnold equations for one dimension (black lines), for two dimensions (red lines) and for
$\gamma =1.4$
(air; solid lines),
$\gamma =1.1$
(SF6; dashed lines) in comparison with results obtained using our numerical method. Numerical results are presented for a case of air and water (black diamond), cases of different DRs with gaseous mixtures of air and SF6 and liquid water (blue crosses) taken from figure 23, cases of constant DR of DR =
$1.25\times 10^{-3}$
and different initial liquid and gas densities (green squares) taken from figure 34 and a case of NG-LNG (purple ‘x’) taken from figure 33. For each case, the normalised impact pressure,
$p^*$
, the liquid density,
$\rho _l$
, the vertical velocity of the liquid–gas interface,
$u_0$
, and the ratio of the height of the liquid mass compared with the height of the enclosed gas pocket,
$\alpha$
, varied. The initial tank pressure
$p_0$
= 1 bar was constant for all cases.
6.4. Effects of phase transition
Under ambient conditions, the liquid in the tank did not boil. For the boiling point 1 and 2 scenarios, the DRs between gas and liquid phases were
$2.69\times 10^{-3}$
and
$4.08\times 10^{-3}$
, respectively. The DR for the boiling point 2 scenario corresponded closely to the DR of NG–LNG. The mass transfer rate source terms, denoted by
$\dot {m}_v$
for vaporisation and
$\dot {m}_c$
for condensation, depend linearly on their respective coefficients
$c_v$
and
$c_c$
(see (2.7) and (2.8)). Higher values of these coefficients increase the mass transfer rates, which means that the associated vaporisation and condensation processes become faster. We computed the sloshing flow at a temperature of 20
$^\circ$
C and an initial tank pressure of
$p_0$
= 0.3 bar (boiling point 1). At 20
$^\circ$
C, the vapour pressure of Novec 649 was
$p_v$
= 0.326 bar, which meant that the liquid immediately began to vaporise. On the other hand, as soon as the pressure rose above the vapour pressure, the vapour began to condensate.
Time histories of the absolute pressure acting on sensor 3 at the roof of the tank filled with Novec 649 (85 % liquid, 15 % gas) for boiling point 1 (20
$^\circ$
C and 0.3 bar,
$p_v$
= 0.326 bar) for different condensation rates and suppressed vaporisation (
$c_v$
=
$1\times 10^{-12}$
s
$^{-1}$
).
To separately investigate the influence of condensation and vaporisation, we suppressed the respective phase transition process by setting the corresponding coefficients,
$c_c$
or
$c_v$
, to a small value. Our preliminary investigations revealed that the phase transition process did not occur for values of
$c_c$
=
$c_v$
$\leqslant$
$1\times 10^{-12}$
s
$^{-1}$
, because the relaxation time was greater than
$1\times 10^{12}$
s. Figure 36 plots the time history of the sloshing-induced pressure acting on sensor 3 for various condensation coefficients
$c_c$
, where the vaporisation coefficient was set to
$c_v$
=
$1\times 10^{-12}$
s
$^{-1}$
to suppress vaporisation. As can be seen in figure 36, higher condensation rates (i.e. larger values of
$c_c$
) decreased the first impact pressures. This is because a condensation process began when the absolute pressure in the gas pocket rose above the vapour pressure of 0.326 bar. This increase in absolute pressure led to an increase in the condensation rate, causing some of the gas in the enclosed pocket to condense. This condensation process consumed part of the energy, resulting in a lower first impact pressure. Consequently, the volume of the gas pocket for subsequent oscillations after the first impact was reduced. In cases involving phase transition, the absolute pressure in the gas pocket only slightly exceeded the vapour pressure of
$p_v$
= 0.326 bar for
$c_c$
= 1.0 s
$^{-1}$
, resulting in minimal condensation. In the other cases with phase transition, the absolute pressure did not exceed the vapour pressure in subsequent oscillations after the first impact and, therefore, condensation did not occur. For
$c_c$
= 1.0 s
$^{-1}$
,
$c_c$
= 3.25 s
$^{-1}$
and
$c_c$
= 6.5 s
$^{-1}$
, the first impact pressure was approximately 10 %, 25 % and 38 % lower, respectively, than for the case of no phase transition. A higher condensation rate also led to an increased oscillation frequency because the volume of the gas pocket was reduced when the gas partially condensed during impact.
For the tank containing Novec 649, the top part of the tank showing the initial vaporisation of bubbles close to the large liquid–gas interface, obtained after 0.40 s, with
$c_c$
=
$c_v$
= 3.25 s
$^{-1}$
.
The initial vaporisation process may exert a considerable influence on the shape of the liquid–gas interface. Figure 37 shows the growth of vapour bubbles beneath the liquid–gas interface 0.4 s following the commencement of the tank motions. As the initial pressure in the tank was below the vapour pressure of Novec 649 at 20
$^\circ$
C (0.326 bar), the vaporisation process was initiated. As the generation of vapour increased, the absolute pressure in the gas phase also increased until it reached the vapour pressure. Because of the effect of buoyancy, the vapour bubbles rose upwards and diffused into the larger gas volume above the liquid–gas interface. In comparison with cases without phase transition, this process resulted in the disturbance of the liquid–gas interface and the influence on the formation of the gas pocket and the corresponding impact process.
Time histories of the pressure acting on sensor 3 at the roof of the tank filled with Novec 649 (85 % liquid, 15 % gas) for boiling point 1 (20
$^\circ$
C and 0.3 bar,
$p_v$
= 0.326 bar) for different vaporisation rates and suppressed condensation (
$c_c$
=
$1\times 10^{-12}$
s
$^{-1}$
).
To investigate the influence of vaporisation on the sloshing impact process, we varied the vaporisation coefficients and suppressed the condensation process. First, we computed the sloshing flow by considering vaporisation to occur up to 1.55 s. As the initial pressure in the tank was lower than the vapour pressure, vaporisation occurred until the absolute pressure in the tank reached the vapour pressure. Second, we computed the sloshing flow starting from the same initial conditions at 1.55 s (i.e. the same absolute pressure of 0.326 bar), using different vaporisation coefficients to study the influence of vaporisation on sloshing-induced impact pressure. Figure 38 shows the time history of the sloshing-induced pressure acting on sensor 3. As expected, the first impact pressure was unaffected by vaporisation. After the impact at approximately 1.625 s, the gas pocket expanded. This resulted in a pressure drop to below the vapour pressure, leading to liquid vaporisation and an increased gas pocket volume. The higher vaporisation coefficient led to an increased volume of the gas pocket. This caused a decrease in pressure amplitude after the first impact and dampened subsequent pressure oscillations. Furthermore, this higher vaporisation rate led to a lower oscillation frequency (after the first impact) because the gas pocket’s volume increased due to vaporisation.
7. Conclusions
We numerically investigated the influence of phase transition and fluid properties on sloshing-induced pressures. To simulate the sloshing flow, we solved the mass, momentum and energy conservation equations as well as a transport equation for the volume fraction to capture the interface between the liquid and gas phases. Both these phases were considered compressible. We extended the conservation equations as well as the transport equation to account for the flow’s compressibility and phase transition. The latter was represented via a mass transfer model. We verified the numerical method by quantifying spatial and temporal discretisation errors. The numerical solutions converged monotonically. To validate our numerical method, we compared our computed time histories of sloshing-induced pressures and velocities with experimental data and analytical values. Computed and PIV measured time histories of velocities compared favourably. For various DRs, based on a mixture of air and SF6 for the gas phase and liquid water, computed and sensor measured pressures agreed favourably. Furthermore, the computed and measured effects of phase transition of Novec 649 at different boiling points on sloshing-induced pressures also agreed fairly well. Finally, the computed sloshing-induced first impact pressures compared favourably also with Bagnold-based calculations.
The gas pocket type of impact pressures decreased considerably with the increasing DR between gas and liquid. For the highest DR of
$6.6\times 10^{-3}$
, the first impact pressure decreased by approximately 50 % compared with the case with the lowest DR of
$1.25\times 10^{-3}$
. The principal effect was largely related to the shape and motion of the liquid–gas interface. The higher DR reduced the velocity of the interface, the distance between interface and tank top and the enclosed gas volume at the time of impact. These effects were primarily due to the increased inertia of the gas with its higher density. The increased volume and mass of the gas pocket lowered the oscillation frequency of the gas pocket and decreased the pressure amplitudes after the impact. We concluded that, for this gas pocket kind of impact, the effects of the DR described above were mainly due to the magnitude of the gas density and the related gas inertia, which significantly influenced the liquid–gas interface and, consequently, the formation of the gas pocket. Thus, we assumed that the DR itself was not the relevant parameter. Furthermore, our results showed that the magnitude of gas density on impact pressure was slightly influenced by the compressibility of the gas.
Regarding the gas pocket impact, its phase transition significantly affected sloshing-induced pressures. During the impact phase, gas condensation occurred, consuming part of the kinetic energy and leading to a decrease of the first impact pressure. Furthermore, after the impact, the pressure amplitudes, induced by the dynamics of gas pocket, decreased. The effects of the vaporisation process occurred after the impact. Immediately after this impact, the gas pocket expanded and the pressure dropped below the vapour pressure, inducing vaporisation and damping the pressure oscillations.
Liquid viscosity hardly affected impact pressures. This was expected because the sloshing flow was pressure dominated.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10643.
Acknowledgements
We thank U. Lantermann for co-operation and fruitful discussions on numerical modelling involving phase transition and thermal effects.
Funding
We gratefully acknowledge the computing time granted by the Center for Computational Sciences and Simulation (CCSS) of the University of Duisburg-Essen and provided on the supercomputer magnitUDE (DFG grants INST 20867/209-1 FUGG, INST 20876/243-1 FUGG) at the Zentrum für Informations- und Mediendienste (ZIM). Numerical modelling and simulations were performed in the framework of a research project funded by the German Research Foundation (DFG) grant 533951202. The experiments were performed in the framework of a research project funded by the DFG grant 298747191.
Declaration of interests
The authors report no conflict of interest.
Open access
