#### Abstract

A two-phase flow model is developed to study violent impact flow problem. The model governed by the Navier-Stokes equations with free surface boundary conditions is solved by a Constrained Interpolation Profile (CIP)-based high-order finite difference method on a fixed Cartesian grid system. The free surface is immersed in the computation domain and expressed by a one-fluid density function. An accurate Volume of Fluid (VOF)-type scheme, the Tangent of Hyperbola for Interface Capturing (THINC), is combined for the free surface treatment. Results of another two free surface capturing methods, the original VOF and CIP, are also presented for comparison. The validity and utility of the numerical model are demonstrated by applying it to two dam-break problems: a small-scale two-dimensional (2D) and three-dimensional (3D) full scale simulations and a large-scale 2D simulation. Main attention is paid to the water elevations and impact pressure, and the numerical results show relatively good agreement with available experimental measurements. It is shown that the present numerical model can give a satisfactory prediction for violent impact flow.

#### 1. Introduction

Violent water impact may occur in many hydrodynamics problems associated with important coastal and offshore engineering applications. Wave breaking in harbors, coastal areas, and offshore platforms, liquid sloshing in tank, and green water on decks are the most common examples of this class of problems. Accurate evaluation of such impact forces and possible structure responses is important for structure safety and disaster prevention. Water impacts are usually characterized by nonlinear phenomena, distorted free surfaces, and large amplitude structure responses for freely floating bodies, and their analysis is very complex. Analytical methods are only available for simple cases such as linear problems while laboratory tests are limited by high cost and technical limitations of the experimental facilities. As a result, there is an increasing interest in numerically simulating distorted free surfaces and their violent impacts on engineering structures.

Due to the possible similarity between green water impact and dam-break flow, dam-break flow theory has been widely used to study violent impact problems due to green water incidents [1–4]. Yilmaz et al. [3] developed a semianalytical solution for a dam-break flow to simulate green water on a deck. Kleefsman et al. [4] studied green water impact problem using a Navier-Stokes solver combined with a VOF model for free surface modeling. The studies modeled a dam-break flow to mimic the water flow on a deck without considering the body-wave interaction. In order to model violent impact problem, analytical solutions and laboratory tests have been used. When the assumption of hydrostatic pressure is made, analytical solutions can be obtained, known as the shallow water equations. Among the solutions, the widely simplest one for a frictionless dry flat bed is Ritter’s solution [5, 6]. Laboratory experiments have been performed to study the violent impact flow problems. Cox and Ortega [7] performed a small-scale laboratory experiment to quantify a transient wave overtopping a horizontal deck and a fixed deck. Ariyarathne et al. [8] studied the green water impact pressure due to plunging breaking waves impinging on a simplified, three-dimensional model structure in the laboratory.

For numerical study of violent impact flow problems, the nonlinear distorted free surface is one of its main difficulties. Among the available strategies to numerically construct an interface, the VOF method is one of the most popular methods in water-surface capturing, first introduced by Hirt and Nichols [9]. The advantages of the VOF method are its mass conservation and being easy to implement. Many improved VOF-typed schemes have been proposed, like PLIC-VOF [10], tangent of hyperbola for interface capturing (THINC) [11], THINC/WLIC [12] weighed line interface calculation (WLIC), and THINC/SW [13] scheme. Another difficulty is the simulation of multiphase flow. Thus we need a scheme for treating both liquid and gas fluids with large density ratios simultaneously in one program. Fully implicit solvers can handle this procedure, but the convergence of iteration in a highly distorted state is still a problem.

To solve the above-mentioned difficulties, a constrained interpolation profile (CIP) method combined with a surface capturing scheme, THINC, is applied to the study. The Constrained Interpolation Profile (CIP) method was developed by Yabe et al. [14] to treat both compressible and incompressible fluids with large density ratios simultaneously in one program to simulate the interaction of gas with a liquid and/or solid, which is based on CIP method in Takewaki et al. [15], Takewaki and Yabe [16], Yabe and Takei [17], Yabe and Aoki [18], and Yabe et al. [19]. The CIP method is a compact upwind scheme with subcell resolution for the advection calculation. In the CIP method, both the advection function and its spatial derivatives are used to construct an interpolation approximation of high accuracy within one grid cell. Since the spatial derivatives are also employed, the interface profile inside the grid is retrieved, and the subcell resolution can be obtained. Furthermore, the CIP method can treat all the phases of matter from solid state through liquid state and from two-phase state to gas state without restriction on the time step from high-sound speed [20, 21]. The CIP-based model was introduced to obtain a robust flow solver of Navier-Stokes equations for dam-break flow problems by Hu and Kashiwagi [22]. In [22], a dam-break and oscillation experiments were performed to validate the 2D numerical model with the CIP method for the solution of advection equation of Navier-Stokes equation and also for the free surface treatment. After that, Hu and Kashiwagi [23] presented an enhanced model for nonlinear wave-body interactions, in which the THINC scheme was combined with the model to treat the violent free surface. Recently, Hu et al. [24, 25] extended the CIP model to 3D simulation of violent sloshing and conducted a 3D simulation of wave-body interactions using CIP method combined with THINC scheme. Zhao and Hu [26] presented an enhanced model to treat body motions due to extreme waves, in which the THINC/WLIC scheme was introduced for the free surface capturing. The above numerical results have shown that this model is able to deal with the free surface problems with reasonable accuracy when compared to experimental data.

The dam-break flow is widely used to examine the performance of various numerical techniques designed for simulating the surface interfacial and impact problems. The development of water flow along the floor after the sudden break of the dam has been a conventional target for numerical studies, but our research interest is in the second stage of the flow development, that is, the flow after impact on the vertical wall. In the second stage, overturning and breaking of the free surface as well as air entrapment are observed, and the computation of these complicated phenomena is a more challenging subject. Many models [27, 28] show departure from the experiments significantly for both impact pressure and free surface elevations after the impact. Park et al. [29] proposed a VOF level-set method for simulating two-phase flows and got satisfactory results.

The objective of the study is to examine the performance of the model based on the CIP method and THINC scheme for simulating the violent impact and the dam-break flow problems. We use two dam-break experiments: a small-scale experiment and a relatively large-scale experiment to validate the numerical model. In the small-scale case, both a 2D- and a 3D-CIP-based models are used for analyzing the water collapse flow and the pressure on the opposite wall to examine the difference in simulations in different dimensions. In the relative large-scale case, we analyze the effect of air cushion caused by the backward plunging water reentering the main forward flow, which plays an important role in the calculation of impact pressure and water height on the second stage of the flow.

This paper is organized as follows. In Section 2, we briefly introduce the CIP-based numerical model for two fluids (water/air). The flow solver and the treatment of free surface are briefly introduced. After that, in Section 3, we perform two representative numerical examples, compared with experimental results and other numerical simulations. The paper closes with a summary of our main conclusions in Section 4.

#### 2. A CIP-Based Model

##### 2.1. Governing Equations

We consider an unsteady, viscous (laminar), and incompressible flow. The governing equations are as follows: where . The last term on the right-hand side of (2) stands for the body pressure, such as gravity.

##### 2.2. The Fractional Step Approach

The time evolution of (2) is performed by a fractional step method in which the equation is divided into three calculation steps: an advection step and two nonadvection steps. In the advection calculation step, the CIP method is applied to solve the hyperbolic equation.

*Advection Phase:*

The advection phase computation of (3) is conducted by the A-type CIP scheme where the advection part can be calculated by a semi-Lagrangian procedure after interpolation function is determined. The principal of the CIP method is to use the grid point value and its spatial derivative (gradient) in two grid points to form a cubic interpolation function to approximate the profile. In the CIP method, each grid point has the gradient information of the physical values, as well as the interpolated values between a grid and the neighbor grid. Hence, the CIP has an ability to solve the advection phase with high accuracy and stability. Details of the CIP scheme can be found in the paper by Yabe et al. [14]. The nonadvection phase is divided into a diffusion part denoted by nonadvection phase (i), which includes a viscous term and a source term, and a state-related part to treat the velocity-pressure coupling, denoted by nonadvection phase (ii).

*Nonadvection Phase (i):**Nonadvection Phase (ii):*
An implicit scheme is used for the time integration of (5) as
Taking divergence of (6) and substituting by using (1), we obtain the following pressure equation:

This is a Poisson type equation for the pressure calculation. Equation (7) is assumed to be valid for liquid and gas phases. A solution of it provides the pressure distribution in the whole computation domain. In this treatment, the boundary condition for pressure at the interface between different phases is not required.

##### 2.3. Free Surface Capturing Method

In the dam-break simulation, two density functions and are defined to denote liquid and gas phases, respectively. One has . The free surface is treated as an immersed interface. The free surface boundary is distinguished by a density function , which is solved by the equation

After all of the density functions have been determined, any physical property , such as the density or the viscosity, can be calculated for each computation cell as follows:

The drawback of the averaging process of (9) is that the computational accuracy is reduced to first order in terms of cell size at the interfaces.

The THINC scheme proposed by Xiao et al. [11] is used for free surface capturing of incompressible flow. Some test examples indicate that the scheme has the features we need for our computations: mass conservation, a lack of oscillation, and smearing at the interface. Similar to the CIP scheme, the profile of inside an upwind computation cell is approximated by an interpolation function. Instead of using a cubic polynomial in the CIP scheme, the THINC scheme uses a hyperbolic tangent function, which is shown as follows: where , , and are parameters to control the quality of the numerical solution. Parameters and are used to avoid interface smearing and are determined as Parameter determines the steepness of the jump in the interpolation function varying from 0 to 1. Larger corresponds to sharper variation of . However, too large may result in numerical instabilities. In this paper is chosen according to many computation tests. Also a numerical flux is used in calculation of the density function; therefore, mass conservation could be satisfied. Figure 1 illustrates the concept of the one-dimension THINC scheme. Multidimensional computations are performed by a dimensional splitting method.

#### 3. Validation of Numerical Model

The dam-break flow is widely used as a benchmark test of violent impact flow on the deck to check the computation accuracy for a largely distorted free surface. Its simple initial and boundary conditions make the validation of computation and experiment quite straightforward. In this study, we reproduce two dam-break experiments numerically: one is in a small-scale tank performed by Hu and Kashiwagi [22], in which we also conduct a 3D calculation to test the ability of simulating complex distorted free surface flow with 3D air-water interfaces, and the other is a relative large-scale experiment performed at MARIN by Zhou et al. [30] to examine the applicability of the model in different scale cases.

##### 3.1. Small-Scale Experiment

###### 3.1.1. Problem Setup

This test case was originally performed for validation of a 2D-CIP-based model by Hu and Kashiwagi [22]. As this is a small-scale tank, the viscous effect at the tank walls might be an important factor. In this case, we use 2D- and 3D-CIP-based models to analyze the water collapse flow and the pressure on the opposite wall. The initial conditions are plotted in Figure 2, where point A denotes the pressure sensor point. In the simulation, variable grids are used. The grid points are concentrated near the floor and the right-hand wall. Grid 2 is used as that in Hu and Kashiwagi [22].

###### 3.1.2. 2D Computations

In the simulations, three interface capturing methods, VOF, CIP, and THINC, are used for the 2D computation, and the results for the free surface variations (s and 1.1 s) are depicted in Figure 3. The free surface is indicated by density function contours, with the three lines showing , 0.5, and 0.95. The middle line is . The distance between the two lines and 0.95 roughly represents the transient distance from liquid to gas. This thickness should be zero in a physical problem. As pointed out by Hu and Kashiwagi [22], the free surfaces computed by the THINC scheme are very compact, while those obtained by the VOF and CIP schemes are diffusive after the water collapses. For the VOF method, there appear flotsam and jetsam which are droplets disconnecting from the free surface and the surface shape is not good. It should be noted that the VOF method used in this paper is Simple Line Interface Calculation (SLIC), which is responsible for the result. Since there are many improved VOF methods, this drawback may have been overcome.

(a) s |

(b) s |

The 2D impact pressure using different free surface capturing methods is shown in Figure 4 together with the experimental data by Hu and Kashiwagi [22]. The first peak pressure occurs at the impact instant, while the second peak pressure is induced when the overturning water hits the underlying water and causes a jet splash-up. The three free surface capturing methods can well estimate the impact instant, the second pressure peak instant, and the duration of the first peak pressure. However, roughness and overestimation of pressure are detected in numerical results at the later stage of the simulation. The difference between the numerical results and the experimental data may be caused by the three-dimensional water-air interactions.

###### 3.1.3. 3D Computations

The full scale 3D computations are run using the same scale as the experiment [22] to minimize errors due to scale effects. Figure 5 shows a snapshot of the simulated three-dimensional free surfaces by VOF (a) and THINC schemes (b). In the VOF result, there are many small droplets, close to the free surface, which are due to the reconstruction and displacement of the free surface. Evident significant improvement of free surface by THINC is found from the figure. There are almost no flotsam and jetsam appearing. The free surface is much smoother and the number of droplets is much smaller than that in the original VOF method.

**(a)**

**(b)**

The simulation of 3D dam breaking has been run on three different grids to investigate the behavior under grid refinement. The number of grid points used is , , and . All three grids are focused towards the bottom of the tank and the opposite wall. Figure 5 shows the pressure at the opposite wall. Expected small differences can be found as the grid resolution increases. The three grids predict the pressure equally well, compared with the experimental data. The simulation of 3D dam breaking has been run using three different free surface capturing methods. The middle grid () is adopted. Figure 6 also displays the predicted pressure using different free surface capturing methods. All three free surface capturing methods can predict the first pressure peak well. But its value is overestimated numerically. For the second pressure peak, the THINC scheme can predict it well. Both the CIP scheme and the VOF method underpredict the second pressure peak both in time and value. Beyond the second pressure peak, the presented code overpredicts the pressure. But the same trench is shown as the experimental results and the 3D results are more in line with the experimental results than the 2D results as shown in Figure 4.

**(a)**

**(b)**

Figure 7 displays the evolution of the 3D water collapse and interaction with the opposite wall. During the initial stage of the simulation, the flow remains almost 2D when running up the vertical wall and overturning back to the water tank. After s, however, the splash flow produced by the overturning water mass quickly develops into a 3D breaking wave pattern with the presence of small water droplets and trapped air bubbles. It is clearly seen that the CIP-based model combining the THINC free surface capturing scheme is effective enough to resolve distorted free surface flow with complex 3D air-water interfaces.

##### 3.2. Large-Scale Experiment

###### 3.2.1. Problem Setup

In this section, a relative large-scale experimental test performed at MARIN by Zhou et al. [30] is adopted to evaluate the performance of the proposed numerical method. We analyse the free surface with three surface capturing methods and compare our numerical results of water elevation and green water impact pressure with the experimental results. A set-up overview of this problem is illustrated in Figure 8.

The experimental water tank is 3.22 m in length, 1.0 m in width, and 2.0 m in height. The initial water column is placed behind the flap with the size of 1.2 m in length, 1.0 m in width, and 0.6 m in height. For comparison, we use the water heights measured by two standard capacitive wave gauges. The two probe points, H1 and H2, locate at 2.725 m and 2.228 m away from the left wall of the tank, respectively. The pressure history measured by a circular pressure transducer of 0.09 m diameter at a point P2, 0.16 m above the bottom on the right wall of the tank is presented for comparison.

The effects of both the grid spacing and the time step are carefully investigated, and the values shown in what follows are considered reasonable choices. In the computation, a variable grid is used, in which the grid points are concentrated near the floor and the right wall. The grid number is , and the minimum grid spacing and maximum grid spacing are 2 cm and 6 cm. The total simulation time length is s, and the time step is set to s.

###### 3.2.2. Free Surface Profile

Three interface capturing methods, VOF, CIP, and THINC, are also used in this simulation, and the results for the free-surface variation are compared in Figure 9. Similarly, we can tell that the free surfaces computed by the THINC scheme are very compact; the thickness of the computed free surface is two to three times the size of the cell. On the other hand, those obtained by the VOF and CIP schemes are diffusive after the overturning water hits the free surface. For instance, the free surface by THINC is basically one line and all the time except , and are time and initial water height, resp.), while we can see noticeably three lines from that by CIP, especially at and when there is extremely distorted deformation of the free surface. For the VOF method, there appear flotsam and jetsam which are droplets disconnecting from the free surface and the surface shape is not good such as the position of m at and the tip of the overturning wave at . Also the VOF method used in this paper is SLIC.

**(a)**

**(b)**

**(c)**

The conservation of water mass in the computations for different interface-capturing schemes is checked. Figure 10 shows the variation of the total water volume in the tank over the computation time, with VOF “- ·,” CIP “- -,” and THINC “-*-”. Variation of up to −1.0% is found for the CIP scheme, while the variations in total mass for the VOF and THINC scheme are below 0.01%. The conservation is perfect. Among the three free surface capturing methods, VOF, CIP, and THINC, the VOF and THINC scheme act good in terms of mass conservation; while in terms of interface diffusion, the THINC scheme show perfect performance with very compact surface profile and less flotsam. In the following sections, the THINC scheme will be applied as the free surface capturing method.

Figure 11 shows the propagation in time of the water-front toe after the dam breaks. The numerical result (—) is compared with the experimental result (*) and the analytical shallow water solution of Ritter (- ·). The latter is not applicable to the initial instants of the phenomenon because the shallow-water conditions are not verified, so that the water-front location would be largely overestimated with respect to all the reported solutions. Therefore, the Ritter solution has been shifted laterally just to show the tendency of the numerical simulations to approach the analytical shallow-water result after a suitable time interval. We can see the numerical result asymptotically approaches the Ritter solution as time increases. Also the numerical result shows agreement with the experiment before and a reasonable deviation from the experiment after . However, the experiment has a lower progressive velocity than the numerical result. It may be due to the imperfect initial conditions or bottom roughness in the experiments and some physical effects not considered in the numerical model.

###### 3.2.3. Water Height and Impact Pressure

In Figure 12, we compare the water heights () at H1 and H2 points, respectively. In the experiments, standard capacitive wave gauges have been used which are sensitive to the wetted portion of the wire. Hence, the numerical values are deduced from the simulations by taking the water level and deducting the height of the (possibly present) entrapped air cavity. The initial evolution, and , is characterized by the sudden rise of the water level due to the transition from dry-deck to wet-deck conditions. The water-front shape determines the actual growth rate of . Therefore, the differences detected between numerical and experimental data are reasonably due to details of the initial conditions in the experiments and bottom roughness effects. As time passes, and , a smaller growth rate of the water level is observed, which corresponds to nearly flat interface above the wave gauges. The agreement between measurements and numerical data is still satisfactory. A second steep increase of is then observed due to the water overturning which gives an additional contribution to the water height. The overturning water from the right wall reaches at location H1 at and later , recorded by the gauge H2. The computed result of H2 agrees well with the experiments, while the numerical prediction slightly overpredicts the water height of H1. Later on, the experiment and the present simulations differ largely to some extent, especially at gauge H1. The reason for this may be that the water is entrapped with air during overturning and breaking. It makes the definition and measurement of water level difficult. And the fact that gauge H1 locates near the first breaking point of the overturning wave leads to the discrepancy between the experiment and computation. We make a try to analyze the effect of entrapped air on water level at gauge H1. Figure 13 displays the evaluation of the amount of entrapped air “*” and the water height including the contribution of entrapped air “-”. The calculated water height accords well with the experiment measurement before the overturning wave tip reconnects with the underlying water. It can be noticed that the treatment of considering the contribution of the entrapped air overestimates the water height at . After that, the predicted water height including the contribution of entrapped air agrees well with the experimental data. We can see that (1) before when the contributed height of air cushion is zero, the simulation water height agrees well with the experiment; (2) when the water begins to turn back until the wave tip hits the underlying water (around ), the simulation water height including air cushion is a little higher than the experimental measurement. This is because the definition of water level is complicated during this process. If we subtract this contribution from the total water height, the water height may accord with the experiment result; (3) after that, the overturned water meets with the underlying water (around ), and the water height including the contribution of entrapped air corresponds well to the experimental data. Still the limited information about the experiment does not allow for a better discussion of such comparison.

**(a)**

**(b)**

The evolution of the pressure field and free surface for the dam-break flow is presented in Figure 14. Following the figures of the pressure field, we can get the variation of the pressure on the vertical wall. As expected, there are two peaks of the pressure corresponding to times when the water hits the vertical wall and when the overturning water hits the underlying water, respectively. At , a plunging breaker with a large amount of air entrapped in the water impinges on the underlying layer of water, which generates great pressure against the vertical wall and at the breaking point of the bottom ( m). Subsequently, a second plunging breaker forms from the splashed water around and breaks again around , and the zone of concentrated pressure moves backward to m. During this process, the impact pressure on the right wall remains relatively large. After that, the pressure becomes smaller as the flow momentum weakens.

In Figure 15(a), we show the time history of impact pressure measured and computed at P2 on the right vertical wall. Also Figures 15(b) and 15(c) list the reported comparisons by methods of VOF [27] and SPH [28]. Although there is a phase lag of the numerical computation, especially the first peak, which may be caused by the shape of the water front due to the initial set of boundary condition and the effect of compressibility of the air which is not considered in this model, the numerical model predicts two peaks with no appreciable difference in amplitude. In more detail, the numerical results show some oscillations and a little higher than the pressures measured experimentally at P2. This phenomenon could be associated with the drawback of 2D model caused by fast circulatory flow around the entrapped air. Also, since the spatial and temporal pressure gradients are high, an exact pressure measurement at P2 cannot be determined experimentally as pressure sensors need an area to sense the impact pressure (in this experiment this area is a circle with diameter of 90 mm). For the pressure behavior after , the present simulation provides a more slightly improved agreement with the measurements than other numerical results reported in the literature, as shown in Figures 15(b) and 15(c). We can see that the tendency of the impact pressure can be predicted with acceptable accuracy.

**(a)**

**(b)**

**(c)**

#### 4. Concluding Remarks

In this paper, we present a numerical model adopting the CIP method as the base flow solver to deal with multiphase flow problems with complicated free surface deformation. The THINC scheme is adopted in the interface capturing calculation and compared with other two capturing methods, VOF and CIP, from many aspects. This model is used to simulate two violent impact flow problems due to 2D and 3D dam-break flows and validated with experimental results.

Comparisons among the three free surface capturing methods in the two experiments, VOF, CIP and THINC, show that the VOF method (SLIC) has a good performance of mass conservation but the surface profile is not good, and the CIP method causes mass loss and has thick free surface, while the THINC scheme works much better than the other two schemes in terms of mass conservation and the suppression of interface smearing.

In the small-scale experiment, we perform both 2D and 3D computations and get that the complex distorted free surface with 3D breaking wave pattern can be simulated well. In the large-scale experiment, water elevations at two points H1 and H2 in the wave tank are calculated numerically and compared with experimental measurements. Though there is discrepancy at H1 when the water is overturning from the vertical wall, which is mainly caused by the entrapped air, the overall tendency at H2 is pretty good.

The water impact pressure on the vertical wall is numerically computed in two experiments. Though there is some difference, especially in the larger-scale case, the tendency of the impact pressure can be predicted with acceptable accuracy by this model.

It is proved that the CIP-based model combining the THINC free surface capturing scheme is effective in resolving violent free surface flow problem and the 3D calculation can give more details of the violent flow problem. For more accurate computations, we need more numerical model investigation (like 3D simulation or compressible flow consideration) and supplementary experimental study.

#### Acknowledgments

This work is jointly supported by the National Natural Science Foundation of China (51209184 and 51279186), Fundamental Research Funds for the Central Universities (2012QNA4020), Educational Commission of Zhejiang Province of China (Y201225713), Key Laboratory of Water-Sediment Sciences and Water Disaster Prevention of Hunan Province (Grant no. 2013SS03), and Scientific Research Foundation of Shandong University of Science and Technology for Recruited Talents.