This paper focuses on the numerical modelling of wave impact events under air entrapment and aeration effects. The underlying flow model treats the dispersed water wave as a compressible mixture of air and water with homogeneous material properties. The corresponding mathematical equations are based on a multiphase flow model which builds on the conservation laws of mass, momentum and energy as well as the gas-phase volume fraction advection equation. A high-order finite volume scheme based on monotone upstream-centred schemes for conservation law reconstruction is used to discretize the integral form of the governing equations. The numerical flux across a mesh cell face is estimated by means of the HLLC approximate Riemann solver. A third-order total variation diminishing Runge–Kutta scheme is adopted to obtain a time-accurate solution. The present model provides an effective way to deal with the compressibility of air and water–air mixtures. Several test cases have been calculated using the present approach, including a gravity-induced liquid piston, free drop of a water column in a closed tank, water–air shock tubes, slamming of a flat plate into still pure and aerated water and a plunging wave impact at a vertical wall. The obtained results agree well with experiments, exact solutions and other numerical computations. This demonstrates the potential of the current method to tackle more general wave–air–structure interaction problems.
Prediction of wave loading is a key aspect of offshore and coastal engineering. Violent water waves may cause severe damage to offshore platforms, breakwaters and sea walls, etc. Such hazards have been observed frequently in the past in many countries. Recent examples are the winter storms that occurred in the United Kingdom from December 2013 to February 2014, during which huge waves destroyed railway lines and coastal sea defences . These disasters worldwide have resulted in billions of pounds of losses and threatened thousands of lives. Scientific investigations from engineering, environmental and other perspectives are a necessity to understand and mitigate these naturally occurring events.
To gain a deeper insight into these problems from the viewpoint of the hydrodynamics, great efforts have been made to study wave impacts on breakwaters, sea walls and liquid storage tanks, etc., through carefully controlled experiments [2,3], field measurements  and theoretical analysis [5,6] in the past several decades. The extreme impulsive pressures recorded in violent wave impact events can be tens or even hundreds of times those of impacts induced by ordinary non-breaking waves . Intentions to classify the impact events into different types can be found in the work of Schmidt et al. , Oumeraci & Partenscky  and Kirkgoz . Here, we follow the ideas of Lugni et al.  to divide these into three modes including (i) impact of an incipient breaking wave, (ii) impact of a broken wave with an air pocket and (iii) impact of a broken wave with water–air mixing (see pages 8 and 9 of  for more details). The second impact mode of an overturning wave with an air cavity (also named a plunging wave or plunging breaker) is of particular interest in this study.
Although air pockets and bubbles form in plunging breakers, the influence of air on waves was traditionally ignored owing to its small density (e.g. 1.225 kg m−3 at 15°C and atmospheric pressure) compared with water (around 1000 kg m−3 for fresh water and 1025 kg m−3 for seawater). However, laboratory and field observations [2,3,5,10,11] disclose that air may play an important role in the impact process. During the transition of a plunging wave, an air pocket (or pockets) may be trapped in the body of the wave and compressed by the water mass; thus, a portion of the wave energy will be transferred to the pocket. Once the wavefront impacts the surface of the structure, the air pocket starts expanding to release the stored energy. The strongest pressure peak in the form of a ‘cathedral-roof’ shape and subsequent pressure oscillations will be experienced by the structure. This distinct phenomenon has been discovered in experiments [2,10,11]. In addition to the entrapped air pocket, pure water might also be aerated with many small bubbles through biological production, capillary entrapment, white capping and wave breaking . These bubbles generally persist for many wave periods especially in seawater. The peak pressure and impact duration are strongly influenced by trapped air pockets and entrained air bubbles. This might be expected to be closely related to the compressibility of the air and water–air mixture [12–14]. Furthermore, negative pressures essentially gauge values below the atmospheric pressure have been recorded in field measurements  and laboratory experiments [2,3,8]. Bullock et al.  stated that negative pressures have the potential to induce large seaward forces resulting in the removal of blocks from masonry structures. Crawford  pointed out that such large forces could produce a sufficient overturning moment to cause overall failure of important structures such as breakwaters. Additionally, Lugni et al.  indicated that even small pressure fluctuations might induce local flow cavitation for conditions close to the cavitation threshold. Therefore, it is necessary to be aware of these issues and the need to include all relevant physics when theoretical analysis, experiments or numerical computations are used to investigate wave impact problems. However, we note that very little attention has been paid to the possibility of wave impact-induced flow cavitation following Lugni et al.'s work.
Compared with experimental investigations of plunging wave impacts, which have made significant progress regarding measurement of peak pressures and forces [2,3,10,11], numerical simulations are not yet adequate to fulfil industrial and academic requirements owing to the extreme complexity of these problems. Not to mention the challenge of resolving the free surface, which might overturn, break and experience further strong deformations, the compressibility of air and the water–air mixture and possible flow cavitation make the problem much more difficult. Traditional numerical wave tanks (NWTs) developed in hydrodynamics are mostly based on single-fluid incompressible potential flow theory [15–19]. Because air is not explicitly considered in the mathematical model, computation of entrapped air pockets and/or entrained air bubbles in waves cannot directly be achieved with these single-fluid NWTs. Two-fluid NWTs based on the incompressible Navier–Stokes equations have been proposed to simulate both the liquid and gas phases for violent wave breaking problems [20–25]. However, these treat both the water and air as incompressible fluids, which means the density of each fluid remains constant throughout the process. Unfortunately, compressibility effects in the air pocket and water–air mixture cannot be handled properly by these models, nor, importantly, cavitation (change of phase) effects.
More recently, researchers have started to explore the importance of the compressibility of air and the water–air mixture for wave impact problems. Peregrine et al.  and Bredmose et al.  proposed a weakly compressible flow model combined with a single-phase potential flow solver to compute wave impact events. A fully conservative flow model based on the compressible Euler equations was adopted in the impact zone to describe the water wave with entrapped air pocket or wave with entrained air bubbles. In the energy equation, only the compressibility of the gas phase was included without considering compressibility in the liquid phase. Systematic numerical analysis of wave impacts on vertical walls was conducted, and promising results were presented. However, model tests of a one-dimensional shock wave passing through a water–air interface exhibited strong non-physical pressure oscillations at the material interface. These oscillations are a well-known numerical artefact when using a conservative variable scheme and measures should be taken to preclude them. Plumerault et al.  also carried out studies of aerated-water wave problems. They described the flow with a three-fluid model for gas, liquid and gas–liquid components in the flow. Energy conservation was not enforced explicitly with the corresponding equation removed from the equation set in their mathematical model. Instead, a pressure relaxation method  was employed to solve the three-fluid model. Strong non-physical oscillations also arose at the material interface in their water–air shock tube results . They did not present results for wave impacts at structures like vertical walls, but gave solutions for a deep-water breaking Stokes wave in the incompressible limit. We observe that these works still have deficiencies in treating the material interface between water and air.
The objective of this work is to develop an appropriate numerical model for violent water wave impact problems. A fundamental requirement is that the model be able to deal with compressibility effects of air pockets and a water–air mixture and to produce physical solutions with no or very low order spurious oscillations. To achieve the stated objective, a quasi-conservative volume-fraction-based compressible two-phase flow model, which includes the advection of a volume-fraction function and conservation laws for mass, momentum and energy, is presented in this paper. This model can deal with dispersed-phase flows, so that it is more capable for water–air mixtures than other two-fluid models that can only simulate separated-phase flows [28,29].
The remainder of the paper is organized as follows. The compressible flow model for water waves with air cavity and aeration is presented in §2. The numerical method, which uses a finite volume approach, third-order monotone upstream-centred schemes for conservation laws (MUSCL) reconstruction and the HLLC approximate Riemann solver to compute the convective fluxes in the governing equations, is presented in §3. Numerical results for a gravity-induced liquid piston, free drop of a water column in a closed tank, water–air shock tubes, slamming of a flat plate into still pure and aerated water and a plunging wave impact at a vertical wall are presented in §4. Final conclusions are drawn in §5.
2. Flow model for dispersed water waves
In hydrodynamics, ordinary (non-breaking) water waves are conventionally considered as separated-phase flows, in which the free-surface separates the gas and liquid components completely as shown in figure 1a. However, for breaking waves as mentioned in §1, the water may be entrapped with air pockets or entrained air bubbles when the crest starts to curl. At the same time, water droplets may also be ejected into the air as spray as illustrated in figure 1b. In this case, the flow is of dispersed gas and liquid phases. In order to construct the mathematical model for wave impact problems, we adopt the homogeneous equilibrium approach to make the following assumptions:
The bubbly fluid is assumed to be a homogeneous mixture of air and water.
Each component obeys the conservation laws of mass, momentum and energy.
The mixture obeys the conservation laws of mass, momentum and energy.
The temperature T, pressure p and velocity V of all the phases and components are identical.
These assumptions are based on the belief that differences in the thermodynamic and mechanical variables will promote momentum, energy and mass transfer between the phases rapidly enough, so that equilibrium is reached [30,31]. The equilibrium model is usually considered an appropriate approach to treat free-surface flows .
(a) Mathematical model
For each individual fluid component i (i=1 for air, i=2 for water), its basic material properties can be described as follows:
Density ρi=mi/Ωi, where mi is mass and Ωi volume.
Pressure p=p1=p2, both components have the same pressure.
Velocity V=V1=V2, both components have the same velocity.
Internal energy .
Kinetic energy .
Total energy .
To determine the internal energy of fluid component i, an appropriate equation of state should be adopted. Because the ideal gas equation of state is not suitable for liquid, the stiffened-gas equation of state is used for both fluid components in this study. Therefore, the internal energy of component i can be described as 2.1 and the total energy is 2.2 where γi is a polytropic constant and pc,i is a pressure constant. The parameters Γi and Πi are defined as 2.3 Additionally, the speed of sound for each component can be calculated as 2.4
In order to describe the material properties of the homogeneous mixture, we introduce the volume fraction function α1 for air, defined as 2.5 Accordingly, we have α2=1−α1 for water. Based on these values, the material properties of the water–air mixture can be expressed by the following:
in which N=2. Substituting equation (2.2) into the formulation of mixture total energy, the pressure can be computed as 2.6 The speed of sound for the bubbly water–air mixture can be estimated by Wood's formula  2.7
The mathematical model used here for the flow of the water–air mixture consists of the mass, momentum and energy conservation laws for the mixture. A conservation law of mass for each component is also included. In particular, gravitational effects should be considered and included for water wave problems. Consequently, the underlying conservative part of the flow model can be expressed in the following form 2.8 in which is the vector of conservative variables, is the flux function, are the source terms and these are defined as 2.9 where u, v and w are the velocity components along the x, y and z axes, g is the gravitational acceleration and h is the enthalpy given by 2.10 In addition to the conservative part, the advection of volume fraction function Dα1/Dt also needs to be considered. Here, we adopt Kapila et al.'s  one-dimensional advection equation 2.11 extended to three dimensions 2.12 where K is a function of the volume fraction and sound speed given by 2.13 Equation (2.11) is derived from the pressure equilibrium assumption, and its right-hand side term assures that the material derivatives of the phase entropy are zero in the absence of shock waves. If we neglect the right-hand side of equation (2.11), then this is a standard transport equation for α1 as pointed out by Murrone & Guillard .
The overall flow model includes equation (2.8) and (2.13), and we write it in the following form 2.14a and 2.14b The underlying reason, we do not choose a fully conservative or primitive model, but a quasi-conservative model, is due to the following factors. Fully conservative flow models have the pre-mentioned difficulties at material interfaces where non-physical oscillations inevitably occur even for first-order schemes  owing to a non-physical pressure update [29,37] or negative volume fraction  during numerical computations. Although a primitive variable flow model can avoid these oscillations, difficulties arise when resolving strong shock waves to maintain the correct shock speed . For complicated problems consisting of both material interfaces and shock waves, combining the fully conservative and primitive variable model formulations has previously been found to be an effective strategy. However, this method is quite intricate as switching is required between the two models for the different regions . Quasi-conservative flow models, which combine the conservation laws with a non-conservative scalar (volume fraction or other material property) advection equation, have proved proficient and much simpler in the past [29,36].
If gravitational effects are excluded, and only the x-direction is considered, equation (2.14) reduces to a five-equation system which can be named the five-equation reduced model  or Kapila model . For the one-dimensional five-equation reduced model, Murrone & Guillard  proved that it can be derived from the two-pressure and two-velocity Baer–Nunziato equations in the limit of zero relaxation time. They also proved that the reduced system has five real eigenvalues (three eigenvalues are equal) and is strictly hyperbolic with five linearly independent eigenvectors. Important information about other mathematical properties of the one-dimensional system, which include the structure of the waves, expressions for the Riemann invariants and the existence of a mathematical entropy, can also be found in their work.
The three-dimensional two-pressure two-velocity Baer–Nunziato model has 11 equations in total . The model (2.14) is computationally less expensive as it deals with only seven equations.
3. Numerical method
(a) Treatment of the advection equation
System (2.14) has a volume fraction transport equation which is not conservative. As indicated by Johnsen & Colonius , directly using the non-conservative formula presents difficulties when dealing with the material interface owing to an inconsistency between the wave speeds (shock wave, rarefaction and contact discontinuity) and the velocity vector V. They advise transforming the advection equation into a conservative formulation to overcome this obstacle for one-dimensional multi-fluid problems. Based on their work, here, we move forward to reconstruct the volume fraction equation appropriately for three-dimensional multi-phase flow problems.
Introducing a vector f defined as 3.1 then the divergence of f is given by 3.2a 3.2b 3.2c Consequently, we can obtain 3.3 Substituting equation (3.3) into equation (2.12) to obtain the quasi-conservative form 3.4 then replacing the non-conservative equation in system (2.14) by equation (3.4), the integral formulation of the overall system can now be expressed as 3.5 in which Ω represents the flow domain, ∂Ω is its boundary, the vectors U, F and G are given by 3.6 where (nx,ny,nz) is a unit normal vector across the boundary ∂Ω, q=unx+vny+wnz and Nq=[0,0,0,nx,ny,nz,q]T.
(b) Spatial discretization
For three-dimensional problems, hexahedrons are used to fill the domain Ω. Without loss of generality, if Ω is a cuboid, we may simply partition it using a Cartesian mesh with I×J×K cells. Within the mesh, any cell can be indicated by the subscripts i, j and k. As shown in figure 2, a hexahedral volume cell is closed by six quadrilateral faces. For this cell, six neighbouring cells are selected to form a computational stencil.
Supposing the mesh does not vary with time, the discrete form of equation (3.5) for a mesh cell mi,j,k can be written as 3.7 where the subscript m stands for the mesh cell itself and the subscript f represents a mesh cell face. Introducing the residual defined as 3.8 equation (3.7) can be written as 3.9
As flow discontinuities (shock waves and material interfaces, etc.) may be expected to occur in multi-phase multi-component compressible flows, attention should be paid to the computation of surface flux terms. In this work, the surface flux term F across any mesh cell face f is evaluated by an approximate Riemann solution based on the numerical flux function 3.10 where the symbols + and − indicate the left and right sides of the face considering a normal vector nf. The conservative variables at neighbouring cell centres may be assumed piecewise constant and assigned directly to and respectively, and this very simple but diffusive treatment is first-order accurate in space. To improve the accuracy, the solution data are formally reconstructed using a third-order MUSCL scheme  based on the primitive variables W=(α1,ρ1,ρ2,u,v,w,p)T, written as 3.11a and 3.11b where the subscripts L and R represent the left and right neighbouring mesh cells, respectively, 3.12 in which the subscript LL indicates the left neighbour of the left cell, and RR represents the right neighbour of the right cell. The parameter κ provides different options of upwind or centred schemes 3.13 In this work, is adopted. To prohibit spurious oscillations introduced by the high-order interpolation, a slope limiter function ϕ is used to modify the reconstruction procedure as follows 3.14a and 3.14b In this study, we apply van Albada's limiter defined as 3.15 where 3.16 The reconstructed conservative variables at the left and right sides of a mesh cell face can be easily obtained as 3.17 The numerical flux term represented by equation (3.10) is calculated by employing the HLLC approximate Riemann solver (ARS) for the homogeneous mixture.
(c) The HLLC Riemann solver
The HLLC ARS originates from the work of Harten et al.  and Toro et al. . Its applications for separated multi-phase multi-component compressible flows can be found in the work of Hu et al.  and Johnsen & Colonius  etc. However, to the authors' best knowledge, no work has been reported thus far that solves dispersed multi-phase multi-component compressible flows governed by the one-dimensional Kapila model or three-dimensional system (2.14) with the HLLC ARS. The numerical flux term represented by equation (3.10) is calculated by the HLLC ARS defined as 3.18 in which the middle left and right states are evaluated by 3.19 the intermediate wave speed can be calculated by 3.20 and the intermediate pressure may be estimated as 3.21 The left and right state wave speeds are computed by 3.22a and 3.22b where is an averaged velocity component evaluated by the component Roe averages . The averaged speed of sound is calculated by the formulation proposed by Hu et al.  3.23 More detail of equation (3.23) can be found in reference .
(d) Temporal discretization
In this work, a third-order total variation diminishing Runge–Kutta scheme is applied to update the numerical solution from time level n to n+1 3.24a 3.24b 3.24c For a mesh cell m, the partial differential operator L is defined as 3.25
4. Results and discussion
(a) Gravity-induced liquid piston
This problem provides a fundamental test of the capability of a numerical method to deal with the compressibility of fluids. Figure 3a shows the set-up of a liquid piston. A closed tube of 15 m length is filled with two air pockets separated by a section of water in the middle. The density of air is ρ1=1 kg m−3, and the density of water is ρ2=1000 kg m−3. Initially, the velocity field is still, and the pressure is p=1 bar. Under gravity, the water in the tube will firstly drop, the lower air pocket will be compressed, and the upper air cavity will expand. As the pressure in the lower air pocket increases, the water will decelerate, but continue compressing the lower air pocket until its velocity reduces to zero. Then, the lower air pocket will expand and push the water upwards. If friction effects are ignored and the tube is adiabatic, the water body will keep oscillating like a piston. This problem is very close to the classical Bagnold piston problem . It has been used to benchmark numerical methods for fluid sloshing in a tank by the liquid sloshing community at the ISOPE 2010 conference .
To solve this problem numerically, we equally distribute 150 mesh cells along the y-direction in the tube. The volume fraction level of the gas phase is set to 1−10−5 in the two air pockets and it is set to 10−5 in the water section. During the computations, density and pressure, etc., are recorded. In particular, we place a numerical pressure gauge at the bottom of the tube to facilitate a direct comparison with Guilcher et al.'s calculated results obtained independently  and shown with the present results in figure 3b, where good agreement is noted. A relatively long-time computation of five oscillation cycles is then carried out. Timeseries of pressure and density at the two ends of the tube are shown in figure 4. There is no obvious decay of the oscillation amplitudes for these parameters. This implies that the conservation laws are well satisfied by the present method.
(b) Free drop of a water column in a closed two-dimensional tank
Figure 5 shows the configuration for this problem, which is a benchmark test proposed by the liquid sloshing community at the ISOPE 2010 conference . A rectangular water column (ρ2=1000 kg m−3) with width 10 m and height 8 m is initially at rest in the closed tank and air (ρ1=1 kg m−3) fills the remainder of the tank. The initial pressure is p=1 bar in the tank. Under gravity, the water column will drop and impact upon the bottom of the tank at around t = 0.64–0.65 s. The impact pressure at this moment is of particular interest to ship structural engineers for this type of problem, because it is fundamentally a key issue for the safety of liquefied natural gas carriers.
To obtain an accurate prediction of the pressure peak, use of a fine mesh is recommended by other researchers [45,48]. Therefore, we equally distribute 1500 mesh cells in the vertical direction and 200 mesh cells in the horizontal direction. Figure 6 gives two snapshots of the pressure field in the tank just before impact and at the impact time t=0.65 s. It is clearly shown that the highest pressure occurs at the bottom centre of the tank and the distribution of pressure is symmetric about the central section (y-axis). In figure 7, we present several snapshots of the volume fraction contours in the liquid tank. An interesting finding is that a small amount of air is trapped between the water body and the bottom surface of the tank. This body of air undergoes not only compression owing to the liquid impact (at t=0.65 s), but also expansion when the gas phase pressure exceeds the environmental liquid phase pressure. At t=0.75 s, the portion of trapped air pocket appears to be a very thin layer, then has a cylindrical shape at t=1 s and a half-cylindrical shape at t=1.2 s.
The time history of the absolute pressure at the bottom centre of the tank is plotted in figure 8. Guilcher et al.'s  results computed independently using an SPH code are represented by the red curve with ‘+’ symbol. Braeunig et al.'s computation  is represented by the blue dashed line. We use a black curve to illustrate the results obtained by the present compressible method. The green curve indicates the result obtained on a coarse mesh (200×700 cells) with the present method. We have also used the interFoam module from the open-source software OpenFOAM 2.1.1, which is based on an incompressible two-fluid finite volume method, to compute this problem on the same mesh and present these results with a green line. Comparing the impact time, the present compressible results agree well with Guilcher et al.'s solution and Braeunig et al.'s work (t=0.65 s), whereas OpenFOAM gives a slightly early prediction at t=0.64 s. The incompressible OpenFOAM solver also gives the lowest pressure peak at 8.7 bar. The others produce much higher predictions at over 20 bar. After the peak, the pressure begins to fall. The minimum pressure following the peak obtained by OpenFOAM is about 2 bar; Braeunig et al. produced a non-physical negative pressure with some oscillations ; Guilcher et al.  gave a value of 1 bar; the present compressible method obtains much lower, but positive values for the minimum pressure after the pressure peak. Only the present method permits the fluid to expand sufficiently far to be in tension.
(c) Water–air shock tube
As pointed out by Peregrine et al. , compression waves may form in the aerated water region after waves impact on a structure as the water–air mixture is a highly compressible fluid . Therefore, it is necessary to assess the capability of a numerical method to handle the propagation of shock waves through pure and aerated water in addition to resolving free surfaces. In computational gas dynamics, shock tubes are frequently used to test different types of numerical methods. Here, we consider a group of shock tube problems involving water–air mixtures. The shock tube initially has an imaginary membrane in the middle, which separates it into left and right parts filled with fluids at different states. In total, we report here on two of these shock tube tests, for which the initial conditions, stopping times and number of mesh cells used are listed in table 1. Gravitational effects are not included for these problems.
Computed results for the two cases are shown in figures 9 and 10. A close examination of these results shows reasonable agreement with exact solutions and independent numerical predictions. This verifies that the present method resolves wave speeds correctly and does not produce spurious non-physical oscillations near shocks or material interfaces.
In particular in figure 10, it can be seen, in contrast, that Plumerault et al.'s solution exhibits strong numerical oscillations around the shock and material interface, and the present method has superior performance. A material interface, which is a small step in the volume fraction α and density ρ, appears near the middle of the tube (x=0.5). Apart from the material interface, both the numerical results underpredict the volume fraction and overpredict the density in the region 0.5≤x≤0.6.
(d) Flat plate impact on pure and aerated water
The numerical method is used here to solve a water-entry problem of a rigid flat plate. The water can be pure or entrained with air bubbles, and pure water is firstly considered. Figure 11a shows the computational set-up corresponding to Verhagen's experiments for pure water impact problems . The width of the plate is 0.4 m and it drops with a fixed constant velocity of v=2.8 m s−1. The impact pressure is of particular interest, because it is an important parameter for the safety of structures slamming into water.
As suggested by Yang & Qiu , we also use a uniform mesh with 800×700 cells to discretize the computational domain. Yang and Qiu did not state the initial distance from the bottom surface of the plate to the free surface in their work, we set it to 0.1 m in this work. The present method is constructed under the Eulerian frame with a fixed finite volume mesh. To solve this problem, we fix the flat plate in the mesh and cause the water to move upward with the constant velocity v=2.8 m s−1. This strategy is appropriate for a short-period impact process .
A quantitative comparison of the gauge pressure at the bottom centre of the flat plate is made and shown in figure 11b. The present result agrees well with Yang and Qiu's independent computations with the cubic interpolated pseudo-particle (CIP) method  and Verhagen's experiments . After the peak, the absolute pressure also drops below atmospheric pressure. Verhagen's measurement seems to be not very convincing at around t = 0.0075–0.009 s, when absolute negative values appear in the graph. But generally, the numerical computations and laboratory experiments give almost the same pressure evolution history and peak value before t=0.007 s.
We move forward to use the numerical method to investigate the effects of entrained air bubbles on slamming pressures. The set-up for the aerated water slamming problem is similar to the pure water case as shown in figure 11a. We choose a flat plate of width L=0.25 m and an impact velocity of 4 m s−1. The initial distance from the bottom surface of the plate to the free surface is 0.1 m. The aeration levels in the water are 0% (pure water), 1, 2 and 5%, respectively. We use the same mesh as before and record the pressures on the bottom surface of the plate at three positions: (i) edge of the plate x=0L, (ii) quarter of the plate x=0.25L and (iii) centre of the plate x=0.5L.
Figure 12 shows the computed slamming pressures at these locations. It can easily be seen that the impact pressures are not equally distributed on the bottom surface of the plate. The pressures rise from the edge of the plate towards the centre. The entrained air bubbles in the water reduce the impact pressures at all three gauge positions. At the plate centre, 1% aeration reduces the pressure peak from 12 to 9 bar, and 5% aeration produces a two-third reduction of the peak value to 4 bar. Entrained air bubbles also extend the rise time in the impact pressure time history compared with pure water. The overall duration of the impact loading is also seen to be extended by aeration. Similar phenomena can be seen at the edge and quarter locations on the plate.
Figure 13 shows the pressure wave propagation within the fluid for 2 and 5% aeration levels (only half of the domain is presented due to the symmetry). Before the impact pressures have reached their peak values, the air and water beneath the plate have already been compressed (see the first and second columns of the figure). A compression wave is observed to propagate downwards and outwards. At the edge of the plate, a rarefaction wave is observed as expected. The speed of sound in the 2% aerated water case is higher than the 5% aerated water case. As a consequence, the compression wave propagates faster in 2% aerated water than in the 5% case. At t=0.026 s, the compression wave established in the 2% aerated water has reached the tank floor and has been reflected through a single Mach reflection. In the 5% aerated water case, the compression wave has not yet reached the tank floor.
(e) Plunging wave impact at a vertical wall
Here, we perform an exploratory calculation to establish the viability and promise of the method for violent wave impact simulations involving air pockets and aeration. Figure 14 shows the set-up for a plunging wave impact event. The length of the wave tank is 3 m and the height is 0.8 m. A step of 0.2 m height is placed in the bottom right part of the wave tank starting at x=1.75 m with a 45° slope to cause the wave to steepen and break. A piston-type wave maker is placed at the left boundary of the domain to generate waves. The still water depth is d=0.3 m. The whole NWT is divided into two subdomains. A two-fluid NWT based on the incompressible Navier–Stokes equations developed in our previous work  is used to deal with the left subdomain. This incompressible solver, which is named AMAZON-SC, adopts an interface-capturing method to treat the free surface as a discontinuity in density. We use the present compressible flow model (2.14) to handle the right subdomain, where an air pocket will be trapped or enclosed by the water body. The dashed line in figure 14 indicates the interface between the incompressible and compressible flow solvers. Buffer zones are used near the interface to exchange flow information between the two solvers. Within the buffer zones, one or two layers of mesh cells for each component solver as required are placed on the opposite side of the interface. The flow information including density, velocity and pressure at these mesh cells is obtained from the companion solver domain through interpolation. More details of the coupling of component flow solvers will be reported separately in future work. A background uniform Cartesian mesh is used to overlay the flow domain, and the basic mesh cell is a square with size of h=0.01 m. Solid boundaries not aligned with the Cartesian mesh in the left subdomain are treated using the cut-cell method [24,53].
Before computing the plunging breaker impact problem, we first conduct a simple test to generate a solitary wave using the incompressible solver. The solid step is removed from the NWT, and the right boundary of the domain is treated as an open boundary. The amplitude of the solitary wave is a=0.09 m. The wave is generated by prescribing the paddle movement according to Rayleigh's solitary wave theory (see ). Figure 15 shows the water surface elevations at x=1 m and x=2 m. The wave crest takes 0.52 s to travel between these two locations. Obviously, the phase speed of the solitary wave is equal to c*=1.92 m s−1. The theoretical phase speed for solitary waves can be calculated as . We obtain c=1.95 m s−1 for this test case, so the relative error between the computed and theoretical wave phase speeds is less than 1.5%.
The integrated numerical wave tank is now used to solve the plunging wave impact problem. The paddle is used again to generate a solitary wave with height a=0.2 m. In addition to the integrated NWT, we also use the established stand-alone in-house incompressible two-fluid NWT AMAZON-SC to solve this problem for the purposes of comparison. According to field measurements and laboratory experiments, the first pressure spike in the form of a ‘church-roof’ shape is a key to the safety to structures. Therefore, we focus on this phase of the impact in the current discussion. Computation of the subsequent wave evolution will be considered in future work.
Figure 16 shows the profiles of the free surface at different times. The red solid line represents the results with AMAZON-SC, the stand-alone incompressible solver and the blue dashed line indicates the solution obtained with the compressible solver (i.e. the integrated NWT). At t=2.13 s, the two solvers give almost the same profiles. We note that an obvious discrepancy of the free surface profiles appears at t=2.15 s. Although the wave crests are almost the same distance away from the vertical wall in the horizontal direction, the wave crest obtained by the incompressible solver is higher than the compressible solver. The wave trough obtained by the compressible solver moves upward along the wall faster than the incompressible solver. The free surface beneath the wave crest obtained by the compressible solver is closer to the vertical wall than the incompressible solver. At t=2.16 s, we note that the wave crest obtained by the incompressible solver moves upward significantly higher than the compressible solver, and this trend continues to t=2.17 s when the waves almost impact the wall. The water continuously moves upward after the wave impacts the wall. From the figures, it is not difficult to observe that the trapped air pocket predicted by the compressible solver is much smaller than the incompressible solver. It would seem that compressibility effects play an important role in changing the shape of the air pocket and the free surface. The incompressible assumption appears to lead an overestimate of the volume of the air pocket for this type of problem.
In figure 17, we present several snapshots of the pressure distribution in the wave field at different times. The first row illustrates the results with the stand-alone incompressible solver AMAZON-SC and the second row corresponds to computations with the compressible solver. For the compressible solver, we can clearly see that the pressure in the air pocket increases dramatically, and the pressure rise travels downstream along the vertical wall and tank bottom. At t=1.95 s, a second compression (pressure increase) in the air pocket is captured by the compressible solver. These phenomena are much less apparent in the predictions with the incompressible solver.
Quantitative comparison of pressures is made and presented in figure 18. We gauge the pressures in the air pocket. For this problem, the size of the air pocket is relatively large as its diameter is around 10 cm. A significant amount of the wave energy is stored in the pocket through compression. Consequently, the air pressure rises significantly to about 114 000 Pa as shown in the left part of the figure. Under the constant density assumption, the incompressible solver cannot deal with compressibility effects and only predicts a peak pressure of 104 000 Pa. The rise time of the pressure spike for the incompressible solver is almost 10 times that of the compressible solver. After the first peak, the air expands to a low pressure. The compressible solver predicts a fall to around 9500 Pa, whereas the incompressible solver predicts about 101 000 Pa.
We show the pressure distribution along the vertical wall at different times computed by the compressible solver in figure 19. The pressure in the region (0.3≤y≤0.4) is strongly influenced by the trapped air pocket. The air pocket is continuously compressed from t=2.156 to 2.167 s. It is noted that at t=2.184 s the air pocket undergoes expansion and the pressure reduces accordingly to subatmospheric values. Thus, this local region is experiencing seaward (suction) force. These numerical findings confirm Bullock et al.'s work  that negative gauge pressures indeed occur in violent wave impact events and the resultant seaward force has the potential to cause the removal of blocks from masonry structures.
In their laboratory experiments of overturning wave breaking on structures, Lugni et al. observed that after the strongest first impact pressure peak, the pressure decreases to a value lower than atmospheric pressure and a subsequent second pressure peak is observed much lower than the first one (see figs 4 and 5 of ). The numerical findings produced by the present compressible solver of a steep pressure spike followed by a negative gauge pressure and subsequent lower second pressure peak, etc., agree qualitatively well with Lugni et al.'s experiments , although the wave conditions are not exactly the same.
A compressible multi-phase flow model that improves the representation of the flow physics for violent wave impact problems involving an air cavity and aeration is presented. The model is based on the conservation laws of mass, momentum and energy in addition to a non-conservative advection equation for volume fraction. Detailed derivation of the model is provided. A high-order finite volume scheme based on MUSCL reconstruction and an HLLC approximate Riemann solver for compressible water–air mixtures is used to discretize the integral form of the mathematical equations.
The numerical method is assessed through a series of test cases including gravity-induced liquid piston motion and the free drop of a two-dimensional water column; water–air mixture shock tubes; the slamming of a flat plate into still pure water and aerated water and an exploratory calculation of a plunging wave impact at a vertical wall. The obtained results agree well with experiments, exact solutions and other numerical computations. These results verify the present method and clearly show its advantages over other numerical methods based on single-fluid or two-fluid incompressible flow models for complex violent wave impact problems.
The present method treats the material interface as a fluid flow discontinuity, which is captured by the HLLC Riemann solver. In large flow gradient regions, the slope limiter automatically reduces the accuracy of numerical method to lower orders to stabilize computations. The resulting dissipation will diffuse the interface across several mesh cells. Front-tracking, anti-diffusion or other techniques are planned to be implemented to sharpen the interface.
Our future work will concentrate on a detailed wave impact analysis under a range of wave conditions. We will also consider the computation of wave evolution after the impact. Further work on the coupling of flow solvers and a massively parallel implementation of the current method with extensions to floating-body and fluid structure interaction problems are planned in subsequent phases of work to be reported in future publications.
The authors acknowledge with gratitude financial support from the Engineering and Physical Sciences Research Council (EPSRC), UK under grant nos. EP/J010197/1, EP/J012793/1 and EP/K037889/1.
- Received July 17, 2014.
- Accepted September 4, 2014.
- © 2014 The Author(s) Published by the Royal Society. All rights reserved.