Numerical Simulation on the Impact of A Liquid Square on Rigid Plate and Liquid Layer

Fluids and structures impact is one of the common phenomena in nature, and it widely exists in engineering practice, including ship hydrodynamic slamming, wave impact on offshore platforms, plunging wave on coastal structures, emergency landing of aircrafts at sea as well as impact of ultra-cold droplets and ice lumps under aviation conditions. In this paper, a two dimensional (2-D) solver for Navier-Stokes equations is developed and applied in the numerical simulation of the impact on a rigid plate by a liquid square. The computational domain is discretized by Finite Volume Method (FVM). The Volume of Fluid (VOF) technique is used to track the free surface and the Piecewise-Linear Interface Construction (PLIC) is used for reconstruction. The Continuum Surface Force (CSF) model is used to account for the surface tension. The convective term and the diffusive term are upwind and centrally differenced respectively. The Inner Doubly Iterative Efficient Algorithm for Linked Equations (IDEAL) is used to decouple the pressure and velocity. Based on the proposed techniques, collapse of water column is simulated and convergence study is performed for the validation of the numerical solver. Then the impact of a free falling liquid body is simulated, and the effect of volume and initial height of the liquid body is analyzed. In addition, the impact on a plate with a liquid layer is also simulated to study the effect of falling height on a liquid floor.


Introduction
Recently, the collapse of a liquid body and the subsequent impact on structures has attracted much attention and extensive relevant simulations have been performed. Numerical investigation of the problem can reveal the complicated fluid mechanism, such as the two-phase flow, strong nonlinear interface motion and the splashing jet. Also, the problem has great importance in engineering practice, namely dam breaking. Usually, the liquid square is located on the bottom and the impact occurs on the side wall by the collapsing liquid front. In this paper, we turn to a more complicated case as the liquid square is free falling. Then, the liquid square hits the bottom with a velocity and the impact will be much stronger.
Tracking of the interface is critical when developing a numerical solver for the simulation. Currently, the relevant methods can generally be classified into Lagrangian (moving-grid) and Eulerian (fixed-grid) types. For violently deformed interface, the Lagrangian methods will cause additional efforts for the remeshing. Thus, Eulerian method is more preferred in the present problem. The VOF method is one of the most widely used Eulerian techniques for interface tracking (Meier, 2002). Pioneering work on VOF is derived from Noh and Woodward (1976). However, until SOLA-VOF was developed by Hirt and Nichols (1981), the method was not widely used. Nowadays, some commercial CFD codes, such as CFX or Fluent, have incorporated simple versions of VOF (Ubbink, 1997). Approximate interface locations are found from the technique called interface reconstruction according to the fluid phase information. There are generally five well-known types of the interface reconstruction techniques, which are Simple Line Interface Calculation (SLIC) by Noh and Woodward (1976), Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) by Ubbink (1997), Donor-Acceptor Method by Hirt and Nichols(1981), Flux Line-segment model for Advection and Interface Reconstruction (FLAIR) method by Ashgriz and Poo (1991) and PLIC method by Youngs (1982), among which the PLIC method has been proved to be accurate and stable enough for strong nonlinear flow (Jofre et al., 2010;Zhao and Li, 2003;Chen, 2015).
Surface tension plays an important role in the evolution of the free surface, especially for strong nonlinear liquid motion. Numerically, surface curvature must be determined from a discrete phase field, when significant surface tension is considered in the flow. As formulating SOLA-VOF (Hirt and Nichols, 1981), the surface tension was accounted by a second derivative of interface line. This method has been proved to be rather crude and inaccurate (Meier et al., 2000). The CSF method developed by Brackbill et al. (1992) was the first method which can meet both requirements of robustness and reasonable accuracy (Meier et al., 2002). By casting interfacial surface force into a corresponding volume force, the surface tension item can be included in the discretized Navier−Stokes equation in CSF. The CSF has been widely used by various researchers and has been included in commercial codes such as CFX-4 and Star-CD for its advantages (Meier et al., 2002). Kothe et al. (1996) tested the CSF model with 2-D and 3-D static drop tests, and obtained satisfying results.
The equations of motion for the two-fluid system are closely coupled, and this raises the requirement for special solving methods. Nowadays, many methods have been developed to solve the Navier−Stokes equations, among them, SIMPLE-family (Semi-Implicit Method for Pressure-Linked Equations) algorithms are most well-known (Patankar and Spalding, 1972;Chow, 1979). But the methods are still not satisfying, even very poor for some flow cases (Sun et al., 2008). IDEAL algorithm was then proposed by Sun et al. (2008) to enhance the convergence rate and stability of the iteration process. IDEAL can almost fully guarantee the coupling between velocity and pressure. The inner doubly iterative processes for the pressure equation in IDEAL can almost completely overcome the approximations in the SIMPLE algorithm. Sun et al. (2009) has also proved robustness and efficiency of the IDEAL on three-dimensional incompressible fluid flow and heat transfer problems by comparing IDEAL with the other three widely used algorithms, such as SIMPLER (SIMPLE Revised, Patankar, 1980), SIMPLEC (SIMPLE Consistent, Van Doormaal and Raithby, 1984) and PISO (Pressure Implicit with Splitting of Operator, Issa, 1986).
In this paper, the IDEAL algorithm, PLIC-VOF technique and CSF module are briefly introduced. An efficient and robust numerical solver is then developed based on the Navier−Stokes equations. Simulation for 2-D dam-break is performed firstly and the results are compared with those from previous literature for the validation of the proposed method. On the basis, the impact of a liquid square on a rigid plate is simulated. The free surface evolution and the pressure distribution on the plate are analyzed. The effects of the falling height and liquid volume on the impact are discussed. Furthermore, the impact on the plate with a liquid layer is investigated.

Governing equations
Governing equations for a two-dimensional incompressible flow in primitive form are the continuity equation and the Navier−Stokes equations (Greaves, 2006). The continuity equation states that, in any steady state process, the mass rate entering a system is equal to the one leaving the system. And the Navier−Stokes equations describe the force balance at a given point within a fluid. They can be written as: ∂ρ ∂t + ∇ · (ρV) = 0; (1) where is the fluid density; t is the time; V = (u, v, w) is the velocity vector and u, v, w are the velocity components in x, y, z directions respectively; is the mixed dynamic coefficient of viscosity; P is the pressure; g = (0, 0, −g) is the gravitational acceleration; is the surface tension and is the surface tension coefficient; is the local interface curvature; S = (S x , S y , S z ) is the source item (Chen, 2015). In Eq. (2), the second item of the left and the first item of the right are called the convective term and the diffusive term respectively.
In the discretization of the governing equations, Gauss' theorem is applied to the continuity equation. The discrete convective term and the discrete diffusive term are obtained by the first-order upwind scheme and the central differencing scheme respectively. Then the discrete governing equations can be written as: where, n is the number of faces of a control volume, is the face velocity vector, is the face area vector, A is the coefficient obtained by the first-order upwind scheme, the subscript refers to the location of the variables, with p refers to center of the control volume and f refers to its faces respectively (Chow,1979), is the time step, V p is the control volume, F f is the face flux, is the kinematic viscosity. The discrete Navier−Stokes equations are solved based on the FVM method and the staggered mesh, using the IDEAL method. In the staggered mesh, the velocities are calculated at the control volume faces and all other scalar variables are computed at the center. This type of mesh can avoid the checkerboard pressure field effectively.
The IDEAL algorithm applies multiple iterations in each step, which has greatly improved the consistency of the initial fields and completeness of the discretized pressure equation. The first inner iteration for the pressure equation aims to get initial pressure that approximates the final result of the current iteration enough. The second inner iteration is acted on the pressure equation, too, aiming to obtain the final velocity and pressure of the current iteration which can almost fully satisfy the governing equations. The inner doubly iterations can improve the two approximations in the SIMPLE algorithm greatly and make the IDEAL fully implicit (Sun et al., 2008). Sun et al. (2009) analyzed the performance of the IDEAL algorithm for solving three-dimensional incompressible fluid flow and heat transfer problems by systematic comparison with the other three most commonly used algorithms (SIMPLER, SIMPLEC and PISO). The main conclusions were as follows.
(1) Among the four algorithms, IDEAL algorithm was the most robust and effective algorithm.
(2) For the five problems studied, the IDEAL algorithm can converge at almost any time step.
(3) When each algorithm calculated at its own optimal time step, the computing time of IDEAL algorithm was much less than that of the other three algorithms.
With the remarkable advantages, the proposed IDEAL algorithm can be widely used in the calculation of three-dimensional incompressible fluid flow and heat transfer problems. Therefore, the IDEAL algorithm used in this paper has more advantages than other algorithms in terms of computational efficiency and accuracy, and can deal with the velocity-pressure coupling problem in Navier−Stokes equation more effectively.

PLIC-VOF method
VOF method which was first reported by Nichols and Hirt (1975) is formulated based on the idea of a so-called volume fraction F. It is a scalar function, defined as the integral of a fluid's characteristic function in a control volume. Value of 0 indicates the presence of one fluid and 1 indicates the other fluid, when there is interface in the cell volume, (Seifollahi et al., 2008). Transport equation of F stands for transportation laws of volume fraction, which can be written as: The fluid properties of control volume are calculated in a form of a volume average of the properties of the two fluids. Thus, the dynamic viscosity coefficient and density are defined as (Greaves, 2006): where the subscripts w and a refer to water and air, respectively. Approximate interface locations are reconstructed from the phase information stored in F. To obtain a sharp interface, special, non-diffusive and geometry-based techniques are necessary to advect F in an appropriate way. Earlier versions use rather simple, crude reconstructions. Recently, PLIC-VOF was introduced, which approximates the interface by a straight line within each computational cell (Meier et al., 2002). Vector n is defined as the gradient of the interface. To obtain it numerically, we have (Wang, 2010) δx δy where, and are the grid sizes in the x and y directions, respectively, the subscripts i and j refer to locations of the cell or face of a control volume. Once the normalized unit vector n is calculated, a straight line is positioned perpendicular to it in such a way that it matches with the value of F in the cell. With the orientation of interface, twenty different cases may occur. Four of them have or to be zero. After symmetries and flips, the other sixteen can be simplified to four, as shown in Fig. 1 (Seifollahi et al., 2008).
In the paper, a typical example is selected to compare different VOF methods, and the calculation results of the changes of the free surface in the shear velocity field are given, as shown in Fig. 2.
It can be seen from the results that in the case of large deformation of the free surface caused by the shear velocity field, the FLAIR method gets very poor results, while the TVD method can get better results, but the boundary width of the free surface obtained after reverse shear is too large, and the outer isoline can no longer remain circular. In this example, the PLIC method can obtain more reasonable and accurate results in the shear velocity field, thus, the PLIC-VOF method is chosen in the paper.
2.3 CSF model Kothe et al. (1996) tested the CSF model with 2-D and 3-D static drop tests, and the results showed that the original CSF method can accurately simulate the surface tension effect in strongly nonlinear liquid motion and obtain very satisfying results. CSF converts the interfacial surface force into a corresponding volume force, which can then be included in the discrete Navier-Stokes equations. The surface tension can be explained by the Young−Laplace equation (Meier et al., 2000) σ κ κ where, is the surface tension coefficient, is the local interface curvature. In the CSF, it is postulated that the neighboring eight grids contain all the essential information concerning the local curvature, for the convenience to establish the polynomial approximation for using the nine volume fractions (Meier et al., 2000).

Convergence study
In the convergence study, four meshes as listed in Table 1 are used. The total number of grids is referred as n. The left, right, upper and bottom boundaries of the computational domain are all no-slip walls. The initial pressure on wet surface depends on the static pressure purely by gravity, and that on dry surface is defined as 0. The densities of water and air are 999.04 kg/m 3 and 1.226 kg/m 3 , respectively. Dynamic viscosity coefficient of water is 1.1379 kg/m 3 and for air is 1.7825 N•s/m 3 . The gravitational acceleration g is 9.8 m/s 2 .
At initial time, a liquid square with width d = 0.7 m and length l = 0.7 m is confined above the plate. Falling-height h  between the lower surface of the liquid square and the bottom boundary is chosen to be 0.1 m. The impact is simulated in a tank with the dimension of 2 m×1 m, as shown in Fig. 3. The origin is defined as the lower-left corner.
Interface evolutions of the liquid square under different grids numbers are shown in Fig. 4. Owing to the internal assumption of VOF method, there is inevitably a thin air layer just before the exact contact of the liquid square and the plate, which will lead to scattered cavities in the liquid and kind of divergence of the numerical simulation. This is almost inevitable and cannot be eliminated in all the numerical simulation as strong impact occurs. Besides, the effects of the cavities and drops on the impact force are almost negligible. Thus, from practical perspective, the interface converges very well with respect to the grids number.

Validation by data from previous literature
For further validation of the present numerical method, a two-dimensional dam-break model is simulated and the results are verified by previous simulation. The model by Col-agrossi and Landrini (2003) is adopted, as shown in Fig. 5. A water column is confined by a plate in a tank with the dimension of 3.22 m×1.8 m. When the plate is lifted at t = 0 s, the water column with initial height of 0.6 m and length of 1.2 m flows freely. The model is simulated based on the orthogonal mesh, and the number of cells is 161×90. The noslip wall condition is imposed on all the rigid boundaries. Fig.6 compares the free surface evolutions. There is an upward jet in Fig. 6a when t = 1.19 s. The upward jet starts to overturn away from the wall when gravitational effect counteracts the upward velocity, causing a downward liquid tongue as shown in Fig. 6b. Some air is trapped by the liquid tongue when it retouches the water surface, and a cavity presents in Fig. 6c. A secondary jet happens in Fig. 6d caused by the impact between the downward liquid tongue and the water surface. It can be seen from Fig. 6 that, both the present simulation and that by Colagrossi and Landrini (2003) can demonstrate the free surface evolution very well.

Results and discussions
The present method is validated by convergence study and comparison with previous literature in Section 3. On this basis, the effects of the falling height h, length of the square l and the liquid layer are discussed in this section.

Effect of falling height on the impact
To observe the differences when the liquid square falls from different heights, the free surface evolution, pressure distribution and the impact force are investigated as the square falls from 0.1 m, 0.2 m, 0.3 m and 0.4 m. The numerical model with h = 0.1 m is shown in Fig. 7. The liquid body is defined to be 0.4 m×0.4 m. The computational domain has its length to be 3 m and height to be 0.6 m, 0.7 m, 0.8 m and 0.9 m, respectively.
Free surface evolutions for different h are shown in Fig. 8, in which the time t is replaced by a dimensionless time

366
GUO Chun-yu et al. China Ocean Eng., 2020, Vol. 34, No. 3, P. 362-373 Fig.8a is captured at the moment when the liquid just starts to touch the plate as = 0. It can be seen at this moment, there is no obvious difference between the four interfaces. As the dimensionless time τ τ τ proceeds to 0.5, as shown in Fig. 8b, the central part of the interface is lower and the liquid fronts is wider for larger falling heights, which indicates the stronger collapse of the liquid square. As proceeds to 1, the liquid fronts as h= 0.3 and 0.4 m have already touched the side walls and upward jets have already been formed, while those as h = 0.1 and 0.2 m move much slower, as shown in Fig. 8c. As proceeds to 2, the upward jets as h = 0.2, 0.3 and 0.4 m touch the top wall of the computational domain, and small liquid tongues are formed falling against the side walls.
The pressure distributions on the bottom plate corresponding to time instants and the falling heights in Fig. 8 are shown in Fig. 9, in which pressure P is replaced by a dimensionless pressure . Initially as = 0, it can be  GUO Chun-yu et al. China Ocean Eng., 2020, Vol. 34, No. 3, P. 362-373 367 τ seen that the higher falling heights induce stronger pressure on the wetted bottom wall, especially within the central part of the liquid square. The liquid body's falling makes the surrounding air rotate violently. This causes the low-pressure areas on the plate at the just spreading corners of the liquid body. As goes on to 0.5 in Fig. 9b, the liquid square collapses and spreads out, the pressure peak decreases and the pressure distribution becomes more uniform. Because of the different falling heights, the moment when the liquid touches the side walls is different. As shown in Fig. 9b, the liquid front corresponding to h = 0.3 and 0.4 m has touched the side walls and induced a pressure upsurge due to a sec-τ ondary impact. Also, it can be seen that the pressure peaks corresponding to different heights start to fall. In Fig. 9b, the pressure peaks as h = 0.3 and 0.4 m have already fallen below that as h = 0.1 m. As = 1 in Fig. 9c, adverse pressure peaks occurred at the side walls caused by the secondary impact can be found compared with those in Fig. 9a.
The pressure continues to decrease as the time proceeds and turns to negative as h = 0.3 and 0.4 m, as shown in Fig. 9d.

τ τ
History of the impact force R integrated from pressure P on the plate is shown in Fig. 10. As discussed before, higher drop-height will induce greater initial impact, as shown in Fig. 10 at = 0. As proceeds until the liquid fronts touch  the side wall, the impact forces keep falling except that as h = 0.1 m. At this height, the initial impact velocity can be calculated as . The velocity is not large enough to generate a sudden strong impact. Moreover, at this height, the initial impact is even smaller than the static gravitational force. This indicates that the gravitational effect needs time to accumulate due to fluidity, which is very different from that of a solid body. Thus, the impact force as h = 0.1 m reaches its peak as = 0.5, while those for higher liquid square reach their peak at the initial time. When = 1, the liquid fronts as h = 0.2, 0.3 and 0.4 m have collided with the side wall, so different extent of force increases are induced. For h = 0.1 m, the total pressure increase happens at = 1.5 due to its slow initial impact velocity.

Effect of length on the impact
In this section, the effect of the liquid volume is analyzed by changing its length l. The computational domain is similar to that in Fig. 7. The falling height h is kept constant as 0.1 m, and l is varied from 0.4 m to 0.7 m. Fig. 11 presents the free surface evolution at different time instants. It is interesting to find that the fronts of the liquid squares under different lengths move in a very similar manner. The main difference between the four cases exists in the central part of the square, where the gravitational force is more dominant. τ τ Fig. 12 presents the time history of the pressure on the bottom plate. As been stated, at the initial stage of the impact, the air cavity will lead to the irregularity of the pressure, which is shown in Fig. 12a. Also, it is interesting to find that as = 0, the pressure peak by l = 0.4 m is stronger than those by 0.5 m, 0.6 m, and 0.7 m. As = 0.5, the peak value by l = 0.7 m is larger than the others. For pressure distributions of the plate, they have almost similar developments, except the general increment caused by the increment of the volume. Actually, when h = 0.1 m the liquid gravity rather than the impact force is dominant.
The impact force R on the plate by the liquid square with different volumes is presented in Fig. 13. When = 0, the square with least volume actually induces the strongest force on the plate. This coincides with the pressure distribution in Fig. 12a. Then as the gravitational effect gradually accumulates, the impact force increases. This stage corresponds to the first rising of the impact forces. As the vertical momentum is translated horizontally, the forces on the plate decrease. Finally, there will be an increase as the liquid fronts hit the side walls. Between the three stages, there are two inflection points. The first one shows at = 0.5 as l = 0.4 m, = 0.6 as l = 0.5 m, = 0.9 as l = 0.6 m and = 1 as l = 0.7 m. The second one shows at = 1.5 for all the four cases.

Liquid square impact on a plate with a liquid layer
In this section, we present a more complicated case as the square hits a 0.1 m deep liquid floor. In this case, strong liquid splash will appear as the square collapses. Owing to the existence of numerical simulation error, the splashing at the left and right liquid fronts may interfere with each other, which will ultimately affect the simulation result. Thus, only half of the domain is simulated to eliminate the asymmetric splashing at the left and right liquid fronts. It is reasonable Fig. 10. Impact force of the plate. Fig. 11. Free surface evolution with l varying at different time instants. GUO Chun-yu et al. China Ocean Eng., 2020, Vol. 34, No. 3, P. 362-373 to redefine h to be the distance between the liquid body and the water surface, as shown in Fig. 14 At the initial time when = 0, there is a liquid embossment at lower right corner of the liquid square. The embossment grows into splashing jet when = 0.25. The liquid front detaches from the splashing jet, causing the little droplets. Seen from Figs. 15b and 15c, the liquid tongue is elongated and deformed. The little droplets start to touch the side wall when = 0.5. When = 1, the splashing jet has developed into a liquid ribbon. The liquid stacks up at the domain center and forms the shape 'sea elephant'. Due to the gravity, the liquid ribbon shakes its head forwardly when = 1.5 and a small amount of air is trapped. When = 2, the liquid ribbon flaps to the side wall. The water surface on the right is represented by two distinguishable interfaces and much more air is trapped than that in the previous moment. τ Pressure distributions of the bottom plate are shown in Fig. 16. By comparing Fig. 16a with Fig. 9a, it is interesting to find that the pressure on the plate is much stronger if there is a liquid layer compared with no liquid layer. Actually, the square will not collapse as easily as it does without the layer due to the frictional effect, and the vertical impact will not be quickly released horizontally. As time proceeds to = 0.5, as shown in Fig. 16b, the gaps between the four pressure peaks are greatly weakened. The horizontal pressure distribution becomes irregular for h = 0.3 and 0.4 m. This is reasonable because higher falling height will cause more violent formation of the water head and the ribbon like splashing jet. As the convex liquid head propagates to the right, local pressure peak will occur accordingly, as shown in Fig. 16c and Fig. 16d.
The impact force on the bottom plate is shown in     Fig. 17. It is found that the greatest impact does not occur at the initial time = 0 for smaller falling heights, but it does occur at = 0 for larger falling heights. By comparing Fig. 17 with Fig. 10, it is found that the initial impact is much stronger if there is a liquid layer due to the frictional effect. When = 1, the impact force as h = 0.1 m has be-come the largest among the four heights, while h = 0.4 m, the smallest. The impact forces keep decreasing as h = 0.3 and 0.4 m, while the forces increase first and then keep decreasing as h = 0.1 and 0.2 m. The collapse propagates much slowly due to the effect of the liquid layer. Thus, the secondary upsurge is not shown in Fig. 17, as it is in Fig. 10.  GUO Chun-yu et al. China Ocean Eng., 2020, Vol. 34, No. 3, P. 362-373 371

Conclusions
Slamming is a drastically relative and interactive motion between structures and fluids in a short time. Slamming will produce a huge impact load, and its destructive power is quite amazing. Accidents caused by wave impact on ships and marine structures often lead to structural deformation, or even structural damage and fracture, such as the waist fracture of hull, resulting in a fatal blow. Therefore, this paper is of great reference value to solve the problems of fluid impact in engineering applications. In the present paper, a liquid square body impact has been simulated using a numerical solver based on Navier−Stokes equations. The computational domain is discretized by staggered mesh and structured rectangular elements. The velocity and pressure coupling is attained by the IDEAL algorithm. And PLIC-VOF method is implemented to track the interface evolution in the impact flow. Successful applications of the model to the liquid square impact of different cell numbers and classical dam-break problems show the convergence and accuracy of the model. The numerical solver is then used to investigate the impact of a free falling liquid body on a rigid plate. The effects of falling height and liquid volume are discussed. With the increase of falling height, the speed at which the liquid front touches the side walls and produces a secondary impact also increases, resulting in a larger pressure peak and impact force, accompanied by a more complex law of motion than at smaller falling heights. It is interesting to find that the increase of the liquid volume by length can hardly affect the motion of the collapsing front. A more complicated case is studied, in which a liquid body impacts on a water surface with a velocity. It is found that the impact is largely increased if there is a liquid layer.