A compressible two-phase flow model for pressure oscillations in air entrapments following green water impact events on ships

3.1.Algorithm

The governing equations are discretized by means of a finite volume method on a fixed Cartesian grid with staggered variables. The velocities are defined on cell faces while the density, the pressure and the curvature are defined in cell centers.

Using cell labelling, we distinguish between the liquid phase, the gas phase and representations of structures in the domain. Cells completely filled by structures are labelled B(ody) cells. The cells filled with air are labelled E(mpty) (the E-label name is a residue from when ComFLOW was one phase). Cells with some liquid, adjacent to E-cells are labelled S(urface) cells. All other cells are labelled F(luid) cells. This means that a F-cell never connects with an E-cell. Note that a F-cell is not necessarily completely filled with liquid.

Time integration of the momentum equation is implicit for the pressure, and explicit for the convective and diffusive terms of the momentum equation. The convective and diffusive terms are integrated in time with a second-order Adams–Bashforth scheme. The time discrete version of the continuity equation is

By substituting Eq. (6) in Eq. (5), a Poisson equation for the pressure is obtained

The density in the cell center at time level n is calculated by

The liquid fraction is indicated by Fs and solved using Eq. (4) by reconstructing the free surface with SLIC in every cell. The flux through a cell face is calculated as the velocity times the area of the cell face times the time step [20].

3.2.Viscosity

Wemmenhove [35], and also Plumerault [27], neglect the term (∇u)T in the viscous term, as is common to do. Wemmenhove et al. [36] do include the term in their mathematical description, but do not evaluate its effect. The effect of neglecting this term is evaluated for the 2D rising bubble case in Section 5.3. In matrix form the stresses become as follows. For readability the compressible term is left out of Eq. (6). For the final simulations in the Results section, the compressible term is included.

where τ indicates the shear stress. The boxed terms are added when the term (∇u)T is not neglected. As the viscosity is variable around the interface and ∇·u≠0 for air, the term ∇·(μ(∇u+∇uT)) is not equal to ∇·(μ∇u) as assumed by others [27,35].

The dynamic viscosities at the top and bottom of the staggered control volume in horizontal direction are needed to find the derivative, see Fig. 2. These are found by linear interpolation between the corner point viscosity values μn,w and μn,e, μs,w and μs,e, respectively. These corner point values between adjacent cell centers are found by harmonic averaging [27]. For computing the local average viscosity at a pressure point, Eq. (8) is used for μ(Fs,μl,μg), where μl is the dynamic viscosity of the liquid and μg of the gas.

Fig. 2.

Corner points of μ for staggered control volume.

4.1.Local height function

SLIC has flotsam and jetsam (small droplets disconnecting from the surface) as major drawback. As a remedy, SLIC is combined with a local height function, consisting of three cells around the S-cell in all axis directions [14]. Instead of updating the volume fraction of the S-cells separately, the height function is updated, after which the water is redistributed depending on its original surface orientation. With a local height function flotsam and jetsam are practically absent from simulations. We use the same height function to assess the local curvature for the application of surface tension.

4.2.Curvature

To implement surface tension, the curvature κ needs to be calculated in every center of a S-cell. To calculate the mean curvature, the local height function based on the surface orientation is used [20]. The grid-aligned free surface orientation is determined by rounding the gradient of the height function. When the free surface is oriented in x-direction (i.e. horizontally), the curvature for the free surface can be calculated from (see Fig. 3)

Fig. 3.

Notation for curvature κ.

4.3.Gravity-consistent density interpolation

Like the pressure, the density is defined at cell centers. In the discretization of the governing equations, the density is also needed at the cell faces. Several alternatives for calculating the density at cell faces are available. We employ a cell-weighted average of the adjacent-cell center values, see Fig. 4

Fig. 4.

Notation for cell-weighted averaging.

It is demonstrated by several authors [13,16,30,36] that this method leads to spurious velocities around the free surface. From the perspective of offshore applications where the gravity force dominates, spurious velocities are caused by an imbalance between gravity and the pressure gradient.

To balance these forces, both terms need to be discretized in the same way. The requirement ∇×(ρg)=0 can be found from the momentum equation. To solve in a way that meets the requirement, [36] came up with a gravity-consistent method (without reference to whether it is applicable for SLIC)

Fig. 5.

SLIC (dashed lines) of interface for density interpolation to cell faces.

Applying the gravity-consistent averaging method with SLIC, it prevents spurious velocities in many, but not all circumstances. Especially near cells with volume fractions of 0.5, the gravity-consistent method gives large errors in combination with SLIC. This because of the lower accuracy of the free surface reconstruction compared to PLIC. As an example, the requirement ∇×(ρg)=0 rewritten in integral form ∮ρgdS=0 is calculated, assuming that both S-cells in Fig. 5 have a volume fraction of nearly 0.5 with the same free surface orientation. This is worked out with numbers in Table 1; it gives the sum of ρg over the dashed red lines, using a gravity vector of g=[−10,−10]T[m/s2] perpendicular to the free surface. The non-zero residue of the gravity-consistent method in Table 1 yields spurious velocities, whereas the cell-weighted method of density averaging does not. Note that besides the free surface configuration illustrated in Fig. 5, there are many other configurations that have non-zero residues, leading to spurious velocities.

The effect of the residuals in Table 1 is illustrated in Fig. 6 for a domain size of 1 [m]×1 [m], 30×30 cells, and the orientation of the gravity vector and free surface mentioned above. The maximum spurious velocity reached after 3 [s] using the gravity-consistent method is 0.4 [m/s] while using cell-weighted density method a maximum velocity of 1·10−9 [m/s] is obtained. It was found that gravity-consistent density averaging combined with SLIC produces spurious velocities, just like cell-weighted averaging. The errors for cell-weighted density averaging appear to be smaller and more similar for different free surface configurations than gravity-consistent density averaging. For this reason, combined with the fact that it requires less computational effort, we chose cell-weighted averaging together with SLIC for the remainder of this article.

Fig. 6.

Liquid fraction- and velocity field: (b) gravity-consistent averaging (max. 0.4 [m/s]) versus (c) cell-weighted averaging (max. 1·10−9 [m/s]).

4.4.Capillary forces

Because we chose an aggregated-fluid approach to keep computational costs in check, the capillary force is added to the momentum equation as a body force. Two options considered for the implementation of the body force are Continuum Surface Force (CSF) [4] and Sharp Surface Force (SSF) [13]. Of the two, SSF is formally more accurate, but CSF is less involved and has similar practical accuracy, because the error resulting from the imbalance between pressure and surface tension is dominated by how the curvature of the free surface is estimated [1].

The density in CSF is averaged between phases (ρ˜=12(ρl+ρg)) to reduce spurious velocities in high density ratio flows [4]. The delta function δΓ in F below Eq. (2) is equal to |∇Fs| and the surface tension coefficient σ is assumed constant. The force is discretized over a (momentum) control volume consisting of two half (continuity) cells where either cell can have a surface orientation with a certain curvature. It was demonstrated by [13] that a face-centered CSF implementation performs better than a cell-centered one. Our discretization of ∫V1ρnFσndV in x-direction becomes

Fig. 7.

Cell face value of κ when using SLIC.

5.1.2D standing viscous capillary wave

Fig. 8.

Setup for simulation of capillary wave.

Standing wave simulations can be used to asses the performance of the numerical method for free-surface waves. All important free-surface dynamics are included, while the domain is conveniently limited [34]. Standing capillary waves are driven by surface tension. We simulated them with zero gravity to verify the CSF model used for the representation of surface tension described in Section 4.4. The setup with a domain of 1 [m]×2 [m] is shown in Fig. 8. The density and the dynamic viscosity of the compressible air is ρg=1 [kg/m3] and μl=0.01 [kgm/s]. The density and viscosity of the liquid are varied in three simulations. The dispersion relation for a non-viscous capillary wave at zero gravity is given by [22].

Fig. 9.

Capillary standing wave for different ratios compared with analytical solutions (a) ρ¯ and μ¯=1,000; (b) ρ¯ and μ¯=100; (c) ρ¯ and μ¯=10.

5.2.2D planar oscillating rod

Another test case for the CSF model is an initially square 2D planar rod of liquid in gas where oscillations are generated by capillary forces. This case has a direct relation to an oscillating air pocket entrapped by a wave impact. Our numerical results are compared with the results of [31] who used ANSYS Fluent.

For the simulation, a square of liquid with an area of 4·10−4 [m2] is used. This 2D square should become, due to the capillary forces, a 2D circle with a diameter of approximately 2.26·10−2 [m]. The surface tension coefficient in the simulations equals σ=2.36·10−2 [N/m]. The domain size is 0.04 [m]×0.04 [m] with 40 equally spaced cells in both directions. A fixed time step of 1.0·10−3 [s] is used. The liquid density and gas density are unvarying and equal to ρl=790 [kg/m3] and ρg=1.2 [kg/m3], respectively. The dynamic viscosity of the gas and liquid phase is 1.0·10−3 [Pa·s]. These settings were also used by [31]. A major difference is that we used 1,600 cells, where they used 25,262 prism elements.

The simulation results are presented in terms of diameter as a function of time, shape of the rod and the pressure. The change of the rod diameter over time is compared with the results of [31] in Fig. 10(a).

Fig. 10.

(a) Diameter and (b) average pressure of the bubble over time compared with [31].

The figure shows that the diameter oscillates towards the theoretical diameter and that the oscillations decrease over time as a result of both physical and numerical dissipation (the latter becoming less with grid refinement). The pressure in the bubble at equilibrium should be equal to 2σ/D=2.09 [Pa]. The pressure in Fig. 10(b) converges to a value somewhat higher than the analytical value due to a systematic error made for the curvature [29]. Renardy and Renardy [29] showed that the integral effect of the curvature converges to a value different than the analytical value. This is also confirmed by the results of [7] for an initially static 2D droplet. According to [15], in the inviscid limit, the angular frequency of oscillation for a 2D planar rod is given by [22] as

In Figs 11 and 12 the volume fraction and pressure over time are compared. Note that a different color scale is used than by [31], but these graphs are presented to demonstrate that our shape and our pressure maxima and minima match with their results, especially at the beginning. There is less of a match in Fig. 12(e). This is because the size of the oscillations in our method did not attenuate by the same amount as in [31] at the time of the snapshot; our method has less dissipation. The same conclusion is found from Fig. 10(b).

Fig. 11.

Volume fraction against time compared with [31] (left); (a) 0.01 [s], (b) 0.05 [s], (c) 0.09 [s], (d) 0.13 [s].

Fig. 12.

Pressure against time compared with [31] (left); (a) 0.49 [s], (b) 0.51 [s], (c) 0.67 [s], (d) 0.75 [s], (e) 0.79 [s].

5.3.2D rising bubble

The following test case is for the combination of buoyancy (gravity), viscosity and surface tension. Our results are compared with the benchmark for a rising bubble [17]. This benchmark was created due the absence of analytical solutions and used for quantitative comparison of incompressible interfacial flow codes. The initial fluid configuration of the 2D rising bubble test case is shown in Fig. 13.

Fig. 13.

Flow domain for 2D rising bubble.

In the benchmark, both water and air are incompressible. The density of the water and air are 1,000 [kg/m³] and 100 [kg/m³], respectively. Further parameter values are shown in Table 2.

Table 2

2D rising bubble parameter variations. Simulation 3 corresponds to the benchmark values

Test	μw[kgms]	μa[kgms]	σ[Nm]	Re [–]
Benchmark	10	1	24.5	35
	0.01	1·10−3	0	35·103
	0.01	1·10−3	24.5	35·103
	10	1	24.5	35

The spatial mean rise velocity vc is found from the simulations and compared with the benchmark. It is calculated as

Before making the comparison between our implementation, the results of [36] and the benchmark, we investigated the setup with parameter variations. These simulations are indicated in Table 2 with numbers ranging from 1 to 3. For these simulations, a grid of 40×120 cells was used. Figure 14 shows the spatial mean velocity of the rising bubble for the three simulations. Figure 15 shows the different rising bubble geometries.

Fig. 14.

Spatial mean velocity for all the rising bubble cases using the original method ComFLOW with a grid size of 40×120.

Fig. 15.

Snapshots of the 2D rising bubble in order of time for the three cases given in Table 2.

From the figures, we find that the evolution of the spatial mean velocity and the geometry of the bubble are highly dependent on surface tension and viscosity. Without surface tension (simulation 1), the rising velocity after the penetration of the jet is lower than with surface tension, because the bubble becomes wider as it rises. With a high surface tension for simulation 2, the bubble does not become as wide and does not slow down as much. The bubble reaches a higher maximum velocity and the bubble’s acceleration (after 2.2 [s]) occurs earlier in simulation 2 than in simulation 3 due to larger viscous stresses in simulation 3. Simulation 3 has the same parameters as the benchmark.

With the implementation in [35], we found by varying parameters, that we were never able to capture the benchmark’s spatial mean velocity when it is at maximum. After more careful consideration, it was concluded that the viscosity model was incomplete. Upon adding the boxed terms in Eq. (10) a better comparison with the benchmark was obtained. This is demonstrated in Fig. 16. With the same grid size of 80×240 and only the implementation of the missing viscous stress components, the difference in mean spatial velocity with the benchmark was reduced from 3.0% [35] to 0.3% (present implementation).

Fig. 16.

Spatial mean velocity compared with the benchmark of [17] with a grid resolution of 80×240.

To investigate how our method deals with two merging interfaces, the free surface and the air-water interface of the bubble, an additional simulation was performed with a similar air bubble interface configuration and a lowered free surface. It is shown in Fig. 17 for a numerical simulation how the rising air bubble protrudes through the free surface. Note that this event was not part of the benchmark.

Fig. 17.

Snapshots of the 2D rising bubble passing through the interface in order of time.

5.4.1D shock tube

By entrapping an air pocket between the water and the structure, the pocket is compressed and can have a cushioning effect on the peak pressure during a wave impact [5,27]. The modelling of the compressibility of the air is tested with the simulation of a shock wave. Note that in our case a non-conservative momentum equation is solved which results in diffused shock waves. The interest for our slamming applications, however, is not in the exact position of the shock, but rather on the associated pressure levels.

The simulation is based on [10], who derived an analytical solution for a 1D shock tube. The tube is simulated with unit length in two simulations. It is completely filled by gas (Fs=0), using 400 cells for one simulation and 600 cells for the other. On either side, the velocity and the gradient of the pressure are set to zero. The time step is unvarying and equal to 3.33·10−7 [s] and the specific ratio for air, γ=1.4. The initial values in the domain are

Figure 18 shows the initial configuration of the shock tube, divided in a driver section with the higher pressure and a driven section with a lower pressure. The figure also shows the relevant stages of the evolution of the pressure. When released, two propagating fronts are created, moving in opposite direction, the shock front and the rarefraction. The pressure immediately upstream of the shock is called the contact surface (p2). The Mach number associated with these two pressure levels can be found from

The pressure downstream of the shock after reflection from the domain wall has taken place (p3 in Fig. 18) can be calculated with

Fig. 18.

Initial condition and relevant stages of the pressure in a shocktube simulation.

The results of the numerical simulations are plotted in Fig. 19. The figure shows the pressure in the domain for different moments in time. The simulated shock front moves with an average speed of ≈560 [m/s] which corresponds with Ma=1.71. This is in agreement with the analytical results.

When using 400 grid cells, pressure values p2=324 [kPa] and p3=915 [kPa] are found. When using 600 grid cells, pressure levels p2=324.5 [kPa] and p3=908.7 [kPa] are found. Figure 19 shows wiggles near the shock front that originate from using central discretization of the pressure with an underresolved shock. As expected, the wiggles and the range in space over which they occur become smaller with increasing grid resolution.

Fig. 19.

(Reflected) shock front in terms of the pressure at different time levels (400 grid cells and 600 grid cells in the enlargement.

It is demonstrated that the simulation results for the 1D shock tube converge for a larger number of grid cells and that they converge to the analytical values. The wiggles observed near the shock front become smaller with an increased number of grid cells. They do not grow in time and are not expected to interfere with our interpretation of the pressure levels in wave slamming events with enclosed air pockets.

5.5.Dam-break experiment

The final comparison before moving on to our main result is for the onset of a wave impact event. The present implementation is validated against the experiments of [23] who focus on the evolution of the free surface in a dam-break event. A dam-break is a characteristic model for wave impact events.

The domain and initial condition for the experiment by [23] is shown in Fig. 20. The size of the domain is a=0.584 [m] and b=0.350 [m]. The size of the dam of water is l0=0.292 [m] and h0=0.146 [m]. The parameters for water and air are set to ρw=1·103 [kg/m3], ρg=1 [kg/m3], μw=1·10−3 [kg/ms] and μg=1·10−4 [kg/ms]. The gravitational constant is set to g=9.81 [m/s2] and the surface tension is equal to σ=7.2·10−2 [N/m].

Fig. 20.

Setup dam-break 2D case.

When the dam is released, the initial water level drops and a front propagates towards the opposite end of the domain. The free surface in [23] was measured along the left wall as an elevation h(t) and as a position of the front l(t) along the bottom. The simulation results are shown in Fig. 21, where h(t)/h0 is plotted against dimensionless time. The simulation results are compared to the experiment of [23] as well as more recent experimental and numerical results [11,19,21]. Three different grid resolutions, 20×12,40×24 and 80×48, were used. It is intriguing to observe that the simulation results converge away from the experimental results, i.e. the coarsest-grid simulation has the best agreement with the experiments. This is consistent with [19], but at present there is no explanation. The differences may be caused by not representing the friction between fluid and bottom well and by 3D effects [19,24].

Fig. 21.

Change of 2D dam over time in (a) height (h) and (b) length (l).