Cylindrical SPH Simulations of Water Entry

eprints@whiterose.ac.uk https://eprints.whiterose.ac.uk/ Reuse This article is distributed under the terms of the Creative Commons Attribution (CC BY) licence. This licence allows you to distribute, remix, tweak, and build upon the work, even commercially, as long as you credit the authors for the original work. More information and the full terms of the licence here: https://creativecommons.org/licenses/


INTRODUCTION
The Smoothed Particle Hydrodynamics (SPH) method is a robust numerical modeling technique without the use of fixed mesh grids. It originated from the 3D astrophysical applications [1] and its 2D form in a Cartesian coordinate has been widely used in the hydrodynamic calculations [2], as well as a variety of industry fluid flows [3][4][5][6]. Although fully 3D SPH simulations can disclose many interesting flow features, such as the large deformation of water surface, multi-interface and moving boundary [7,8], the CPU expense is quite demanding for a large flow system where numerous SPH particles have to be used. In spite of some rapid progresses in the computing acceleration technique, i.e. GPU and parallel algorithm [9,10], a practical 3D SPH simulation in engineering scale is still a challenging issue.
Quite a few interesting hydrodynamic problems involve the axisymmetric features, such as those happening during the water entry of a cylindrical plate and a circular sphere, or the collapse of an axisymmetric column of water and granular materials. Under a Cartesian coordinate system these must be computed by using truly 3D models, although some 2D simulations have been documented in special case studies [11]. The original concept of axisymmetrical cylindrical SPH was proposed by Coleman and Bicknell [12] and Stellingwerf [13], and later by Petschek and Libersky [14] for the hypervelocity impact and long rod penetration. Omang et al. [15] developed a new kernel function for spherical and cylindrical coordinates based on the fundamental interpolation theory of SPH for various shock tube tests. By this formulation they studied the problem of two colliding spherical shocks using a variable SPH particle resolution. Ming et al. [16] further expanded the study of Omang et al. [15] and investigated the effect of charge parameters on an underwater contact explosion. In their study, a dynamic boundary particle approach was proposed to improve the pressure computations near the symmetric axis.
One common feature of the above-mentioned approaches is that they are established based on a modified kernel function under the specific cylindrical coordinate system. Therefore, both the complexity of the numerical algorithm and the CPU cost increase significantly compared with the counterpart of Cartesian-coordinate SPH. To overcome this drawback, Lee and Cho [17] simply normalized the SPH kernel function and transformed the SPH formulations under Cartesian coordinate into axisymmetric space. This practice showed quite promising results on an application of a cylindrical rod impacting on rigid wall. Besides, a sensitivity study indicated that the numerical results are sensitive to the values of artificial viscosity and kernel smoothing length.
On the other hand, another type of cylindrical SPH approaches, which is based on the modification of SPH kernel approximation rather than the kernel function itself, was proposed by Yang FE-18-1095: Shao 4 et al. [18] for an explosively-driven metallic tube and Baeta-Neves and Ferreira [19] for a jet formation induced by the cylindrical-shaped charge. Following benchmark study of Brookshaw [20], both of them used a cylindrical plane density to replace with the original volume density in the Lagrangian form of SPH continuity and momentum equations. Besides, Seo and Min [21] modified both the SPH kernel function and kernel approximation for the study of elasto-plastic contact during the low-velocity impact process. Gong and Liu [22] also tentatively explored the idea of axisymmetrical SPH for simulating the vertical water entry of a circular disk under cylindrical coordinate. The latest work in this field could be attributed to Tavakkol et al. [23], in which a general curvilinear SPH framework was developed. This technique has been shown very efficient for the problems in which one dimension is much larger than the other, or with higher gradient of the parameters in only one dimension.
In this study, we aim to develop a simple and effective cylindrical SPH method for axisymmetric hydrodynamic applications. The numerical model will be based on the Lagrangian governing equations represented under a cylindrical coordinate system. The standard SPH kernel functions will be used, but the kernel approximations will be modified, which is similar to the study of Yang et al. [18]. This treatment should lead to a more concise numerical algorithm, and meanwhile, save considerable CPU expense in large hydrodynamic computations. In addition, further improvement will also be made on the density re-initialization of the particles for the purpose of maintaining numerical stability and accuracy.
The paper is structured as follows. Section 2 presents the cylindrical SPH governing equations and solution procedures, with the treatment of initial and boundary conditions. Section 3 is the model validation through a benchmark cylindrical water column collapse. Comparisons with a fully 3D SPH computation under the Cartesian coordinates are also made to demonstrate the advantage of the cylindrical SPH approach in view of its CPU cost and numerical error. Section 4 is the practical model application to water entry problems with a focus on the fluid-structure interactions, through the use of a cylindrical disc and circular sphere. The convergence of numerical scheme is investigated by comparing the SPH simulations with different particle resolutions in the spatial domain.

Governing Equations
By assuming a compressible fluid, the Lagrangian governing equations in a cylindrical coordinate can be written, in the form of continuity and momentum equations, as where is the density of fluid, t is the time, v is the flow velocity, p is the pressure, g is the gravitational acceleration, and r, z and are the coordinates in the radial, axial and circumferential directions, respectively.
For axisymmetrical problem with and , the governing equations (1) -(4) are simplified as

SPH Formulations in Cylindrical Coordinate
Following the study of Brookshaw [20], the fundamental SPH interpolations under axisymmetric cylindrical coordinates can be represented as (8) where is any physical variable, is the SPH kernel function, h is the smoothing length, and is the particle position vector represented in a cylindrical coordinate.
The discretized form of equation (8) by using SPH particles is represented as (9) where m is the particle mass, j is the neighboring particle, N is the number of neighbors, and is defined as the 2D plane density under a cylindrical coordinate. The first derivative of equation (9) is calculated as (10) By using relevant SPH principles, the second term on RHS of equation (5) can be written as (11) where , , , and are defined.
Therefore, the SPH continuity equation under a cylindrical coordinate reads On the other hand, the RHS of equation (6) can be written in the SPH discretized form as Further considering (14) with combining equations (13) and (14) to obtain (15) Similar derivations can also be applied to the RHS of equation (7) to reach (16) Therefore, the SPH momentum equations under a cylindrical coordinate become To stabilize SPH computations, it is common practice to add an artificial viscosity to the momentum equations (17) and (18). Here the damping term as proposed by Colagrossi and Landrini [24] is used where and are defined, and are the artificial sound speed, and are the velocity vector. The artificial viscosity coefficient  in Eq. (19) plays an important role in the free surface stresses [25].
Then the final form of SPH momentum equations with the viscosity term is as follows Here the improved Gaussian kernel function [24] is used for the proposed cylindrical SPH as (23) where is defined.
The pressure of fluids is computed by the Equation of State, which follows the weakly compressible SPH (WCSPH) approach [2] as (24) where is an exponential constant that normally takes value around 7.0 in hydrodynamic simulations, is the reference density with 1000 kg/m 3 for water, and the coefficient takes 700000.
Similar to standard WCSPH approach, the present cylindrical SPH scheme also uses XSPH to smooth out the particle fluctuations to stabilize the computation. The following correction is imposed after the velocity is solved from the governing equations (25) where  is a smoothing constant in the range of 0.0 ~ 1.0, depending on specific applications. In this paper, the value recommended in [24] is simply adopted, i.e.  = 0.5.

Density Re-initialization for Cylindrical SPH
Following equation (9), we can have By referring to the principles of density re-initialization under a Cartesian coordinate, the following representation can be derived from the division of above two equations as (28) This should be applied every 20 time steps in the computation, as recommended by [24]. Based on the previous computational experiences, the choice of time step seems not to influence this number.
It has been found that this correction technique is able to effectively reduce the density errors arising from the inconsistency of particle mass and volume, and thus significantly improve the SPH numerical accuracy. Similar treatment has not been explicitly documented in the previous cylindrical SPH studies.

Free Surfaces and Solid Boundaries
Due to its mesh-free nature, the free surface boundary conditions in SPH can be automatically satisfied and there is no special treatment required on the free surface. The present SPH model is a single-phase one without considering the influence of the air.
The solid walls are treated by the dynamic SPH boundary approach of Gomez-Gesteira and Dalrymple [26]. That is to say, one layer of the particle is placed on the solid boundary and additional three layers of the dummy particle are placed outside of the solid boundary to maintain the kernel compactness. The number of layers depends on the support length of SPH kernel function, which makes sure that enough particles can be searched near the solid wall boundary. This is independent of the computational particle size. The continuity equation (12) and pressure equation (24) are both solved on the wall and dummy particles. Therefore, the non-penetration boundary condition is satisfied on the solid walls.

MODEL VALIDATION THROUGH A CYLINDRICAL DAM-BREAK FLOW
The dam break flow caused by the collapse of a cylindrical column of fluids is widely found in the engineering field, such as those documented in [11]. During this process, the movement of underwater flow jet, the collapse of submarine medium and formation of the granular material cavity all demonstrate a cylindrical shape. Usually a fully 3D numerical simulation is required to provide accurate flow information under a Cartesian coordinate, although 2D models sometimes provide a reasonable approximation along certain particular planes. By using the proposed cylindrical SPH modeling approach, such a problem can be simplified to 2D but with sufficient accuracy. Meanwhile, the computational expense is also greatly reduced. In this section, we revisit the cylindrical dam-break flow experiment carried out by Martin and Moyce [27] for the purpose of validation of the proposed model.
In SPH model setup, a cylindrical column of water with base radius of is used to study the subsequent collapse and flooding process. The height of water column is , where represents the ratio of dam height to base radius. To compare with the laboratory experiment [27], two parameter settings of the water column, and 2, are considered. For analysis the length scale is normalized by base radius and the time scale is normalized by .
To study the convergence of numerical results, three different particle sizes are used in the initial setup. The particle spacing is given by , where is the number of particle in the radial direction.
Three alternative values of 40, 80 and 160 are used for the comparison of the results. The computational time step is fixed at The choice of time step was based on the Courant sound-speed stability criterion by following [24]. The particles are uniformly placed in 2D r-z plane with r = 0 being the central axis at the beginning of computation. So a true 3D flow scenario can be visualized through rotating the numerical results around axis r = 0 by an angle of 360°.
To compare the accuracy and efficiency of 2D cylindrical SPH with a 3D Cartesian SPH model, the latter computation is also carried out by using the particle numbers of 20 and 40 in the radial direction. Accordingly, a relatively larger time step of is used. To save CPU time, 3D SPH simulations are carried out only in the first quandrant of computational domain thanks to the axisymmetrical feature of the problem. That is to say, only a quarter of the water column located within x > 0 and y > 0 is simulated and the symmetric boundary conditions are imposed at x = 0 and y = 0 planes. In principle, the proposed "symmetric boundary conditions" are based on the standard SPH mirror particle boundary approach, which was initially developed to enforce the correct solid boundary conditions. The difference is that the symmetric boundary here ensures the fluid particles to move following the tangential direction of the symmetric plane and thus complement the motion of fluid particles on another side. Similar treatment has also been done by Liu et al. [28] for modeling the 3D tidal bores. Fig. 1 shows the initial particle configuration in 2D cylindrical SPH and 3D SPH computations for the height-to-base ratio of and radial particle number of , in both side and top views. Table 1 summarizes the CPU efficiency between 2D cylindrical SPH and 3D Cartesian SPH under the different height-to-base ratios of and the different particle numbers of m in the radial direction. The comparisons are made in terms of the total number of time step, total CPU time and CPU time per time step, respectively. It is clearly shown that in the 2D cylindrical SPH computations, the particle numbers in double those in for the same particle resolutions of . As a result, the CPU time per time step increases accordingly in the same rate. Under an identical , the doubling of m leads to an increase of the particle number by four times, and therefore, the same increase in CPU time per time step. On the other hand, the particle number doubles in 3D SPH simulations following the same increase in . Similarly, when m doubles the number of particles increases by eight times in a 3D situation. Although the total particle numbers are similar in 3D SPH with and in 2D with , their CPU expenses differ as much as 20 times. The reason is that the neighbouring particle search is carried out in the order of O(N 3 ) in a 3D SPH, while it is in the order of O(N 2 ) in 2D computations (where N is the total particle number). This could justify the superiority of proposed 2D cylindrical SPH for solving practical 3D axisymmetric problems.
To validate SPH computations, the time-dependent leading edge of the dam break flow relative to the axis of symmetry is shown in Figure 2, compared with the experimental data [27] for the two different height-to-base ratios. For each ratio, the different particle numbers are used in the 2D cylindrical and 3D Cartesian SPH computations, to demonstrate the convergence of numerical results. It is shown that the 2D cylindrical SPH computations can closely match with the experimental data when the particle number in radial direction reaches m = 160, while the same accuracy can be achieved with m = 40 in 3D SPH simulations. This behavior is observed for both and 2. Here it should be pointed out that a full 3D simulation can achieve good accuracy with less particle numbers, which is due to that more robust governing equations are used. In SPH point of view, the particles are more concentrated with their neighbors and therefore, sufficient particles can always be found in the computational process. On the other hand, the 2D cylindrical SPH, although it is equivalent to 3D model in principle, the situations could become different in the numerical aspect. In the initial setup, the particles closer to the symmetric axis possess a small volume. During the motion, these particles could move to a spacious space but still with their initially fixed volume. Then this should cause some numerical errors especially for a long-time SPH simulation. With regard to the 3D SPH results, for the case of , the total particle number is only half of the case . So it has a more strict requirement on the spatial resolution, leading to the obvious effect that results are much Besides, it clearly shows the pressure evolutions from being highly non-hydrostatic at the early stage of dam break to being almost quasi-static at the later stage. Very smooth and stable pressure fields have been obtained with negligible pressure noises. This implies the accuracy of pressure solution schemes in the present 2D cylindrical SPH approach.
As mentioned earlier, a density re-initialization technique is imposed to reset the density every FE-18-1095: Shao 13 20 time steps. This aims to eliminate the accumulation of density errors increasing with the computational time. According to [24], the mass conservation in SPH can be naturally satisfied because each particle has a fixed mass and if the number of particles is constant. However, by using the evolution equation of mass the numerical consistency cannot be enforced exactly between the mass, density and the occupied area. To alleviate this issue, the density field should be periodically re-initialized. To demonstrate the effectiveness of this correction technique, Figure 4  instability. This finding is also similar to that documented by Colagrossi and Landrini [24] in their study of two-phase dam-break flows.

MODEL APPLICATION IN FLUID-STRUCTURE INTERACTIONS
In this section, two benchmark water entry problems are investigated by using the proposed 2D cylindrical SPH modelling approach.

Water Entry of a Cylindrical Disc
The problem of solid body entering a free surface of water has received considerable attentions in the field. According to Glasheen and McMahon [29], as basilisk lizards and shore birds run along the water surface, they can support their body weight by slapping and stroking into the water with their feet, thereby producing a ventilated air cavity under the water. They can pull their feet out of the water before the air cavity closure. The relevant foot motions actually explore the hydrodynamic process of low-speed water entry. Due to its practical significance, early works most concentrated on the water entry with a relatively high falling speed, in which the influence of gravity and hydrostatic pressure can be ignored for the simplicity of problem. Then the flow conditions were mainly dominated by dynamic pressure. However, for the water entry with a low speed, the gravitational force cannot be ignored and the influence of hydrostatic pressure, which dominates over the dynamic counterpart, must be accounted for. This leads to a bigger change in the pressure condition during the downward motion of falling object.
In this section, the numerical simulation of a cylindrical disc entering the water is carried out, following the laboratory experiment of Glasheen and McMahon [29]. Two different radices of the disc are studied, i.e. and , respectively. The computational domain is shown in Figure   5, where the initial water depth is 0.4 m and the radius of cylindrical water flume is 0.2 m. The SPH computations are executed for the different combinations of disc radius and Froude number defined by , where is the initial entry velocity of disc when in contact with the still water surface. In the experiment [29], the falling velocity of object was kept nearly the same at .
This was achieved by carefully selecting the weight of disc to match with the flow resistance when the disc reaches its terminal falling velocity. As a result, the disc was subjected to zero acceleration during the later stage of entry process.
In SPH simulation the particles were initially arranged on a staggered grid system with the particle spacing 0.002 m. Then there are 19900 inner fluid particles and 2212 solid boundary particles, including both the disc and tank walls. A constant time step of is adopted. The computation stops when the air cavity region behind falling disc is completely sealed off. This used 220000 ~ 240000 time steps on a standard desktop computer (CPU AMD Athlon64 3000+) and cost around 24 CPU hours. The 2D cylindrical SPH results will be provided on one particular plane under cylindrical coordinate, and 3D scenario can be visualized by rotating the results fully around its central axis. Figure 6 shows the instantaneous closure of air cavity region behind the entry disc, for different disc radiuses R and Froude number r. It clearly discloses the influence of Fr on the closure shape of cavity region. Since Fr number represents the counterbalance between gravitational and inertia forces, it should generate different effects on the fluid particle motion. It is shown from Figure 6  To quantitatively validate the SPH computations, Figure 8 shows the normalized disc entry depth against the square-root of Froude number Fr 1/2 . Two different radiuses of disc are studied, i.e. R = 0.01 and 0.02 m, respectively. Meanwhile, the experimental data of [29] are also shown for a comparison. is defined as the depth of disc at the moment of cavity closure. It is obvious from the figure that the computed disc depth is linearly proportional to Fr 1/2 , which is consistent with the experimental results. Although the mass of disc used in experiment had been chosen to balance the fluid drag force during entry, the velocity of disc still decreases somewhat slightly at the initial impact.
In comparison, a constant vertical velocity is imposed on the falling disc in SPH computation, therefore the computed relative depth is always slightly higher than the experimental measurement.
Besides, for the two radiuses of disc, three different SPH particle sizes are used to check the convergence of model, i.e. 0.004 m, 0.002 m and 0.001 m, respectively. This is also shown in Figure 8 and it shows that the refined SPH computations are closer to the experimental data, and the difference between two refined SPH runs is smaller than that between two rough runs. Figure 9 (a) and (b) provide the particle pressure snapshot at certain time during the disc entry, for the study on the influence of density re-initialization on numerical stability. The positive performance of density correction scheme is fully demonstrated by the smooth and noise-free pressure patterns as shown in Figure 9 (a), while a very chaotic pressure field is found in Figure 9 (b) due to the accumulation of density error.

Water Entry of a Circular Sphere
In this section, we apply the model to a more general situation of circular sphere entry into the water to demonstrate the capacity of 2D cylindrical SPH. The high-performance simulation of a snooker ball impacting on free surface by Maruzewski et al. [30] indicated that 3D water entry process is very complex in terms of the flow evolution. Although this is geometrically simple, the sphere impact could constitute a difficult problem in various engineering fields.
Here the accuracy of cylindrical SPH is further examined by comparing the numerical results with experimental data of Aristoff et al. [31]. Rather than studying the standard time history of sphere displacement, we focus on more challenging hydrodynamic features such as the shape of air cavity, and the pinch-off time and depth. According to [31], the pinch-off time is defined as the time duration that the center of sphere spends to travel from the free surface to pinch-off point. These interesting issues have been investigated with the coupled Eulerian-Lagrangian method [32,33], but not in SPH field probably due to the requirement on high-resolution particles in a 3D simulation.
A schematic setup of the computational domain, including boundary and initial particle distributions, is shown in Figure 10. The radius of the falling sphere is with an initial impact velocity of 2.17 m/s that is normal to the water surface. The radius of the cylindrical tank is 0.14 m with a water depth of 0.3 m. The initial particle spacing is taken to be and a constant time step of is used. The artificial viscosity coefficient  in Eq. (19) is taken to be 0.003.
The original state is described by the uniform particle distribution within 2D r-z plane as shown in the figure, with pressure field being hydrostatic. The next cylindrical SPH simulations are carried out for the sphere with four different densities, as documented by Aristoff et al. [31]. The simulation stops at the time when the cavity region behind falling sphere is completely sealed off. Here it should be mentioned that a coupled SPH model has been used in the computation of the sphere motion. Since only one degree of the vertical motion is considered, the coupled model is relatively straightforward. In principle, the computed pressures are integrated over the surface of the sphere and then used to calculate its acceleration, which is similar to Gong et al. [34].
The computed particle snapshot and pressure field during the sphere entry are shown in Figure   11, for the four different sphere materials. Meanwhile, the experimental photos of Aristoff et al. [31] are shown for a comparison. It is shown that the cavity region and its evolution demonstrate a good agreement between the experimental photos and SPH results for all cases, in spite of a slight delay in the numerical pinch-off time. Besides, the vertical movement of sphere with time is also accurately reproduced. The figure shows a relatively larger error in case of the smaller density ratio as compared with the larger one. This could be attributed to the counterbalance between the gravitational force of the sphere and the hydrodynamic force arising from the neighboring fluids. For a low density sphere, these two forces are nearly equivalent and thus even a slight fluctuation in the pressure field could generate a larger impact on the sphere motion, while this pressure effect becomes much less substantial if the self-weight of the sphere dominates the kinematic motion.
The results show that the pressure field changes from the initial hydrostatic state when the sphere starts going downwards. Then the pressure wave quickly propagates radially downwards inside the tank until reaching the solid walls and reflecting back. Larger density sphere creates more substantial dynamic pressure waves, which is accompanied by a large water splash on the free surface.
The advantage of SPH computations has been clearly demonstrated by its tracking of the free surface deformation. Moreover, the smooth and noise-free pressure distribution supports the reliability of numerical solution scheme. Figure 11 also shows that the trajectory of sphere and the cavity shape near pinch-off point differ with the sphere density. With an increase in density, both the pinch-off depth and the depth of sphere at pinch-off increase accordingly. The pinch-off time becomes longer as well. As soon as the sphere enters into the water, an axisymmetric cavity region forms below the surface with a splashing curtain above.
With the further sphere descending, it transfers momentum to the fluids by driving it radially outwards.
Then this inertial expansion of the fluid is counterbalanced by the hydrostatic pressure, which eventually reverses the direction of radial flow and causes the collapse of cavity. The process of cavity collapse continues until to a moment of cavity pinch-off, when the initial cavity area is split into two smaller ones. Generally speaking, the above cavity evolution features computed by proposed cylindrical SPH model are highly consistent with the experimental studies of [31].
For a more quantitative validation of the numerical model, the time histories of the sphere entry depth are shown in Figure 12, for the four different sphere densities. This clearly demonstrates the accuracy of the present 2D cylindrical SPH computations for all the test conditions. Besides, to gain a more thorough understanding on the characteristics of water cavity during the deceleration of sphere, four normalized parameters as recommended in [31] are also compared between the SPH results, and the experimental data and theoretical solutions in [31]. The comparisons were carried out only for one specific Fr number that was used in the SPH simulations. Generally speaking, there is a satisfactory agreement among three sets of the data across all the sphere densities, but in some cases it can be observed that the SPH results seem to be closer to the theoretical values.

CONCLUSIONS
The paper proposed a cylindrical SPH numeral model for axisymmetric hydrodynamic applications. One distinct feature of the model is that it transforms a practical 3D axisymmetric problem into a 2D plane space, therefore significantly improves the computational efficiency. Three benchmark hydrodynamic tests are carried out, including the collapse of a cylindrical column of water and water entries of the cylindrical disc and circular sphere. Compared with documented experimental data, and also through the numerical sensitivity analysis, it was found that the proposed 2D cylindrical SPH computation can achieve almost the same accuracy as a full 3D SPH, but meanwhile the CPU expenses were greatly reduced.
In the future study, a two-phase cylindrical SPH model should be developed so as to address the effect of compressible air during the fluid-structure interactions. Besides, an explicit pressure and force comparisons would be important for more robust validations of the numerical model.  Comparison between experimental photos [29] and SPH simulations on sphere entry for different sphere densities Figure 12 The time history of sphere entry depth for the different sphere densities Figure 13 The characteristics of water entry cavity for different sphere densities. The symbols denote the dependence on Fr 1/2 of the normalized (a) pinch-off depth z pinch ; (b) pinch-off time t pinch ; (c) sphere depth at the pinch-off z(t pinch ); and (d) ratio of the pinch-off depth to depth at the pinch-off z pinch /z(t pinch )  Comparison between experimental photos [29] and SPH simulations on sphere entry for different sphere densities Figure 12 The time history of sphere entry depth for the different sphere densities Figure 13 The characteristics of water entry cavity for different sphere densities. The symbols denote the dependence on Fr 1/2 of the normalized (a) pinch-off depth z pinch ; (b) pinch-off time t pinch ; (c) sphere depth at the pinch-off z(t pinch ); and (d) ratio of the pinch-off depth to depth at the pinch-off z pinch /z(t pinch )