High-speed water impacts of flat plates in different ditching configuration through a Riemann-ALE SPH model

The violent water entry of flat plates is investigated using a Riemann-arbitrary Eulerian-Lagrangian (ALE) smoothed particle hydrodynamics (SPH) model. The test conditions are of interest for problems related to aircraft and helicopter emergency landing in water. Three main parameters are considered: the horizontal velocity, the approach angle (i.e., vertical to horizontal velocity ratio) and the pitch angle, α. Regarding the latter, small angles are considered in this study. As described in the theoretical work by Zhao and Faltinsen (1993), for small α a very thin, high-speed jet of water is formed, and the time-spatial gradients of the pressure field are extremely high. These test conditions are very challenging for numerical solvers. In the present study an enhanced SPH model is firstly tested on a purely vertical impact with deadrise angle α = 4°. An in-depth validation against analytical solutions and experimental results is carried out, highlighting the several critical aspects of the numerical modelling of this kind of flow, especially when pressure peaks are to be captured. A discussion on the main difficulties when comparing to model scale experiments is also provided. Then, the more realistic case of a plate with both horizontal and vertical velocity components is discussed and compared to ditching experiments recently carried out at CNR-INSEAN. In the latter case both 2-D and 3-D simulations are considered and the importance of 3-D effects on the pressure peak is discussed for α = 4° and α = 10°.


I. Introduction
The problem of the high-speed water entry is classically of interest in the naval field as far as slamming loads on ships are concerned. In this context several theoretical solution have been derived for simplified conditions and a large literature of experimental data is available. Water-entry problems are also very important in the aircraft ditching, that is, the emergency landing on water. The response of the vehicle to this kind of water impact is critical in terms of safety of the passengers and certifications issued by airworthiness authorities includes the success of the airframe in ditching tests. In this context few high-fidelity numerical methods have been developed so far (examples can be found in [1], [2]).
The SPH method has already shown promising results for the simulation of violent water impacts thanks to its accuracy and easiness in following the free-surface deformations (see e.g. [3]- [5]). In the present work an in-depth study and validation of the SPH model is provided for 2D water entries of flat panels with small deadrise angle. Both purely vertical and oblique impact velocity with high horizontal velocity component are studied (sections III and IV). These conditions are of interest for, respectively, helicopter and airplane ditching situations. To this aim the numerical outcome will be compared to experimental measurements and analytical solutions when available. The influence of the air phase is also addressed for the oblique water entry through a multiphase SPH model. In order to accurately resolve the high and localized pressure peaks developed at the impact a Riemannbased SPH solver is used within an Arbitrary Eulerian-Lagrangian (ALE) framework. The choice of the numerical parameters to be adopted, as e.g. the liquid compressibility, is critically and extensively discussed on the base of physical considerations peculiar of water-impact flows.

II. Adopted SPH scheme
In the present work Euler equations for compressible fluids are solved. Indeed, since the Reynolds number of the flow is quite high and only the impact stage is simulated (short timerange regime) viscous effects can be considered negligible. The weakly-compressible model is adopted; the fluid is, therefore, assumed to be barotropic and a classical stiffened state equation is used: where ρ 0 and p 0 are constant, c 0 is the speed of sound, and γ is a dimensionless parameter greater than 1 (in all of the following examples γ = 7 is used). When considering violent free-surface flows, the proper identification of the reference velocity U re f is crucial, as is discussed in the following sections. Considering the Mach number Ma = U re f /c 0 , the constraint Ma < 0.1 is is enforced to make compressibility effects negligible. Additionally, during violent impact events (i.e., flat impacts) the acoustic pressure p = ρ u c 0 can be reached, and in this case the pressure peak intensity becomes proportional to 1/Ma. On the other hand, in such a condition an incompressible constraint can induce singularities on the pressure field (see [6] for a discussion on the difference between these two models in impact situations). This is linked to the fact that for this kind of impacts the presence of the air phase is generally crucial and the singlephase approach can lead to incorrect pressure evaluations under the incompressible/weakly-compressible hypothesis (for a deeper discussion see also [4]). Being aware of these limits of the single-phase model, the results obtained in this paper have been produced considering a possible Mach dependency.
In the present work the Riemann-based solver described in [7] is adopted. In that work the ALE formalism is used allowing for maintaining a regular particle spatial distribution and smooth pressure fields while preserving the whole scheme conservation and consistency of the classical SPH scheme. The introduction of the Riemann-based solver in the SPH scheme leads to an increased stability and robustness of the scheme with respect to the standard SPH formulation. The formalism proposed by [8], Thanks to the introduction of Riemann-solvers the fluxes between particles are upwind oriented and the resulting scheme is characterized by good stability properties. The discrete Euler equations are written as follows: where ρ E , P E and v E are the solutions of the Riemann problem at the interface r i j = (r i + r j )/2, between particles i and j. The particle transport velocity v 0 is obtained as the summation of the particle velocity plus a small perturbation which helps preserving a regular particle distribution (details about the adopted model can be found in [7]).

III. Vertical water entry of a flat panel
In this first section the vertical impact of a flat plate is numerically investigated and validated. Specifically, results of the 2D single-phase simulations are described and compared to the experimental data from a wet drop test performed in [9]. In that work a flat panel (panel length L equal to 64 cm) impacting with a deadrise angle of 4°and a vertical impact velocity U of 6.0 m/s is studied. Measures of pressures at several positions along the plate are taken, allowing for a detailed control of the pressure peak repeatability and for possible 3D effects.
As described in the theoretical work by [10], for these small deadrise angles a very thin, high-speed jet of water is formed, and the time-spatial gradients of the pressure field are extremely high. This makes the test conditions very demanding for numerical solvers. From the potential flow theory by [10], the jet thickness at model scale is about 0.1 mm. It is worth noting that the theory in [10] is formulated for symmetric wedge impacts whereas in the present case the impact of a inclined single plate is considered. More details about the reference solution to be adopted are given in section III-B.
On the base of this theoretical data, the 2D simulation has been conducted using a very high spatial resolution, corresponding to a particle size Δx = 15.6 µm; by referring the latter to the panel length L, the ratio L/Δx is equal to 41, 000. The whole tank depth and width are, respectively, 3 m and 6 m. In order to manage such a small particle size a variable-h technique has been used. Specifically, the particle size gradually changes with a maximum magnification factor of 3,200 between the most refined region and the lowest resolution one (see figure 1). The total number of particles is about 3 million.
For small dead-rise angles water compressibility cannot be neglected when calculating the speed of the water jet U jet (see e.g. [11]- [13]). In the experimental conditions the speed of sound in water is Under such a condition an estimate of the jet speed is U jet ∼ 42U ∼ 250m/s. As noted in [11], this very thin, high-speed jet disintegrates at some distance from the root due to interaction with both the surrounding air and the surface of the body. In this case surface tension plays an important role on the jet evolution. This local and complex physics is not taken into account in the present numerical method but it is not supposed to play a relevant role on the pressure distribution.

A. Selection of the speed of sound for the SPH model
Before starting the SPH simulations the speed of sound needs to be specified. As mentioned in section II, the numerical Mach number usually adopted within the SPH simulations is Ma = U/c 0 = 0.1, which guarantees the weakly-compressible regime (i.e., compressibility plays a negligible role). Nonetheless, for this kind of impact the reference speed for the Mach Fig. 1. Vertical impact of a flat panel. Colors are representative of the SPH velocity intensity. Pressure probe number 12, which is the gauge used for the comparison with the experimental data, is also depicted. number can not be the impact velocity U. Indeed, considering that the intersection point between the horizontal undisturbed free-surface and the wedge surface has a speed equal to U inters = U/sin(α) � 14 U. According to Wagner theory the water jet formed during the impact has a speed higher than U inters . Therefore, if one chooses the speed of sound using the wedge speed U, i.e., c 0 = 10 U, the jet would not form at all, its speed being in the supersonic regime. This is a clear example where the weakly-compressible rule Ma ≤ 0.1 needs to be enforced in a proper way, considering the specific problem at hand.
In the present case the reference speed should be the water jet speed U jet which, however, is an unknown of the problem. Using the theory in [11], the estimate U jet = 40U can be used. Note that, using the latter constraint, the speed of sound in water would result even higher than the real one, c 0 = 400U versus c � 0 = 247U. This means that for the adopted model scale water compressibility effects are not negligible, at least inside the jet region. However, considering that in the jet zone the pressure is close to the ambient pressure, water compressibility effects should not play a relevant role on the local impact loads.
In order to satisfy the weakly-compressible assumption, at least in the impact region, the reference velocity used is U re f = � P max /ρ, where P max is an unknown of the problem and needs to be estimated. Then, an ex-post facto verification of the SPH simulations is required. Using the Wagner theory (which is valid for small deadrise angles as far as the air presence is negligible) the maximum pressure predicted is: corresponding to about 91 bar in our experiment. The water-hammer pressure for this impact is ρ c � 0 U = 88 bar which is the maximum pressure level that can be physically reached in the experiments. The P max predicted by potential flow theory is higher than the acoustic pressure, and this is a further indication that water compressibility cannot be neglected for this problem. Thus, considering P max ∼ 80 bar, the reference velocity is Taking this into account, a good compromise for the SPH speed of sound can be obtained by adopting a reference velocity U re f = 10U which implies a speed of sound c 0 = 100 U (i.e., Mach number Ma = U/c 0 = 0.01). The resulting time step is equal to Δt = 1.5 Δx/c 0 � 0.04 µs which is 250 times smaller than the sampling rate of the pressure probes used in the experiments (i.e., the SPH sample frequency is 25 MHz versus the 100 kHz of the experimental pressure probes). This aspect is expected to influence the observed pressure peak which, for the considered configuration, needs a high time/space resolution to be captured. Further, in order to verify the appropriateness of such a choice, also a simulation using c � 0 = 247U has been run, the results will be shown at the end of section III-B.  B. Comparison between SPH results, analytical solution and experimental data Figure 2 shows the pressure field predicted by the SPH for the flat panel impact. In the same plot the free surface deformation evaluated by the potential flow theory by [10] is reported. Left plot of figure 3 shows an enlarged view of the flow velocity predicted by the SPH in the area of the highest pressure levels. A thin water jet is formed with a thickness of about 0.1 mm corresponding to about ten particles, thus justifying the high L/Δx ratio needed to properly solve such a flow.
In the right plot of the same figure an enlarged view of the pressure field is shown (in this case the displayed pressure range is enlarged too). From this plot it is seen that the highpressure region is limited to an area of 1 mm 2 with a pressure peak of about 60 bar. In the same figure the size of the probe used in the experiments is depicted. Clearly, the pressure sensor has a size much larger than the pressure bulb formed below the water-jet root. This aspect will be discussed in the following section.
In figure 4 the time histories of the pressure measured at probe P 12 are shown. This probe is positioned 23.6 cm from the left-hand edge of the panel (see figure 1). The SPH solution is compared with the experiments and with the pressure peak predicted by Wagner theory. The SPH prediction is between the two reference data sets, and is characterized by highfrequency components due to the fragmented jet of water that pass over the numerical pressure probe. The latter has a dimension of 0.125 mm (which is the size of the SPH kernel support) and the SPH sampling rate is 25 MHz. Both the SPH and the analytical predictions overestimate the experimental data for which the maximum pressure recorded is 34 bar (in the SPH solution the pressure peak reaches 70 bar whereas the analytical prediction is 91 bar). It has to be noted that, even though the Wagner theory is referred to the case of a symmetric wedge entry, according to [14] and [15] the value of the pressure peak should be substantially the same of the case of an oblique flat panel impact (it will be shown in section III-C that this approximation can be quite rough for the present case). Conversely, regarding the entire pressure signal, it is not possible to compare to classical analytical solutions, such as in [10], since they are all formulated for the symmetric wedge entry case.
As mentioned above, air entrainment is expected to play a minor role for deadrise angles greater than 3° [10], [15]. Notwithstanding that, most of the experimental measurements available in the literature for angles close to 4°exhibit pressure peaks much smaller than the one predicted by potential theory (cf., [15]- [18]) and the values measured in the present study are in fair agreement with previous experiments by [15], [16]. As for possible 3D effects, these have been checked by comparing the pressure values on gauges aligned at the same distance from the piercing edge. For the probes positioned in the most central region no relevant differences have been observed. However, the maximum pressure impact measured in the experiment for small dead-rise angles can be also affected by: 1) Changes of the body velocity during the impact stage 2) Rotations of the body during the impact stage 3) Deformations of the wedge surface 4) Sampling rate of the pressure signals 5) Size of the pressure gauge Thanks to the experimental setup adopted the first three points can be neglected, while the last two can play an important role. In order to take into account the experimental sampling rate, the SPH signal has first been filtered using a moving average filter (MAF) reducing the numerical sampling rate from 25 MHz to 100 kHz. The result is illustrated in figure 5. The reduction of the SPH peak due to this filtering procedure is not enough to get a good agreement with the experimental data.
As a further step the SPH pressure has been measured integrating on a circular area equal to the size of the experimental pressure probes and then filtered at 100 kHz.   This result is shown in figure 6. In this case the SPH output is much closer to the pressure recorded in the experiment.
Summarizing, according to the analytical solution very narrow pressure peaks are expected for this kind of impact. Using pressure gauges with a size of few millimeters and a sampling rate of 100 kHz it is not possible to record such a localized event. Having in mind such limits, through SPH it is possible to get predictions close to the experimental data if the pressures are integrated over the experimental gauge area, even when using a speed of sound c 0 smaller than the real one c � 0 . In figure 7 the experimental pressure time histories recorded   fig. 9. The corresponding pressure peak 1/2ρ U 2 P , the actual average pressure peak measured at the probes P max and the related pressure coefficient C P (4) are also reported.
on a sequence of pressure probes is reported. Even if the experimental pressure peaks present some fluctuations, very similar pressure evolutions are recorded with an almost constant time shift at each probe. The SPH results show quite good repeatability of the pressure peaks along the wedge surface as well (see figure 8). In this regard, it is worth comparing the propagation velocity of the pressure peak along the plate, U P . Indeed, in water entry flows the value of the pressure peak should correlate to U 2 P (see e.g. [19], [20]) as: In figure 9 the calculated values of U P for both SPH and experimental results are reported for several probes. Only the probes far from the panel edges have been considered to avoid influences of either the initial impact stage or the final stage, for which the flow self-similarity is not applicable. In the same figure the values of U P obtained from a simulation adopting the real speed of sound c � 0 = 247U (Ma=0.004) are also reported. The difference between the simulation with Ma=0.01 and Ma=0.004 is very small (about 2% of the average value), confirming that within the weakly-compressible regime the Mach number effect is limited. For the sake of completeness in the same plot also the analytical value of U P (valid for a wedge entry problem) is reported, i.e.: In Table I The average values of U P are reported. For the simulation at Ma=0.01 the average value of U P is about 117 m/s corresponding to a pressure coefficient C P = 1.03 which is close to the expected value (eq. 4). On the other hand, since the value of U P predicted by the SPH is smaller than the analytical one, maximum pressure of the SPH cannot be equal to the one predicted by the Wagner theory as shown above. This difference is essentially due to the fact that Wagner theory is referred to a symmetric wedge entry, while the present case is asymmetric (see section III-C).
When considering the experimental pressure probes an average value of U P equal to 104 m/s is obtained (with a standard deviation σ = 3.1). This propagation velocity of the pressure peak corresponds to a pressure coefficient C P equal to 0.56. This confirms that the measurement system adopted is not able to record the real pressure peaks which are therefore underestimated as shown in this section.

C. Comparison between symmetric and asymmetric water entry
In the last subsection it has been shown that the pressure peak predicted by SPH is consistent with the water impact theory. Indeed, the calculated peak propagation velocity U P and the pressure peak P max give a pressure coefficient C P close to unity (see eq. 4). Nonetheless, it still remains a significant discrepancy between SPH results and the analytical solution for this impact angle. In order to investigate the source of this incongruity a further simulation has been performed.
The symmetric problem has been set by retaining the initial particle configuration of the asymmetric case and closing the fluid domain at the left edge of the panel with a vertical wall (see figure 10). The simulation has been run with Ma=0.01.
In figure 11 the peak propagation velocity for both symmetric and asymmetric problems are shown. It is clear that in the symmetric configuration the SPH solution is now much closer to the potential flow prediction and, because of this increase in the value of U P , a larger pressure peak is expected on the panel surface. The recorded pressure at 0.16 m from the left-hand edge of the panel is shown in figure 12 for both symmetric and asymmetric SPH solutions.
Consistently with the observed value of U P , the pressure peak of the symmetric solution is about 88 bar. In the same figures the analytical solutions from [21] and [10] are also reported. Evidently, the symmetric SPH solution is now much closer to the analytical prediction and the remaining difference can be mainly attributed to fluid compressibility. This is in agreement with the analytical work in [22]. In that work it is shown that, generally, asymmetric impacts induce smaller pressure peaks with respect to the symmetric case. This effect is emphasized when the deadrise angle is smaller (in [22] the smallest angle considered is 10°).

IV. Water entry with high horizontal speed at 4°pitch angle
In this section the ditching problem including a large horizontal velocity component is studied, comparing the SPH outcome to model test experiments. Given the higher complexity of the problem, in this section also possible influences of the air phase are investigated, since its role can not be a priori neglected.   [21] (dash-dot line) and by [10] (dash dot-dot line) are also reported.

A. Description of the experimental data
Guided ditching impact experiments [19], [23] were performed in the CNR-INSEAN towing tank, which is 470 m long, 13.5 m wide and 6.5 m deep. The dimension of the flat plate were 500 mm by 1000 mm (Fig. 13). Pressures at 18 points are measured through Kulite XTL123B pressure transducers. The sampling rate of the latter is 200 kS .s −1 while their dimension is 3.8 mm. The horizontal and vertical velocities at the impact are respectively 40 m.s −1 and −1.5 m.s −1 . In the experiments, these velocities are assumed to remain constant during the whole ditching impact.

B. Numerical methodology
In this second test case series, simulations have been conducted using Adaptive Particle Refinement (APR) technique [24]. This technique allows keeping a high spatial resolution around the flat plate during the simulation as the refinement areas move at the same speed of the plate. In order to avoid the pressure filtering described in section III, the numerical adopted pressure sensor size is the same as in the experiments (i.e 3.8 mm).
For both single-phase and two-phases models, simulations were performed with a nominal sound speed c water 0 = 1480 m/s and a nominal density ρ water 0 = 1000 kg/m 3 for water. For the two-phases approach, the nominal sound speed of air was c air 0 = 343 m/s and the nominal density was ρ air 0 = 1 kg/m 3 . More details about the two-phases model adopted can be found in [25], [26]. Comparisons with experiments are established in terms of pressure coefficient, which is defined as : where U and V refer respectively to the horizontal and vertical velocities.

C. Flat panel impact : single-phase approach
The numerical tank is 20 m long and 6 m deep. Simulations involving 6 refinement levels (Δx min = 1.5625 mm) and 8 refinement levels (Δx min = 0.390625 mm) were performed ( Fig. 14). At initialization, the total number of particles for simulation involving 6 and 8 refinement levels are respectively about 67, 000 and 110, 000. Figure 15 shows the pressure coefficient field predicted by the single-phase SPH simulation in the area of the highest pressure levels. It is worth nothing that the field appears very regular despite the contour lines cross several refinement interfaces. In Figure 16, SPH pressure time histories for two different probes are compared with the experimental data from [19], [23]. The origin of the time axis is based on P4 pressure peak. A good agreement between the finest SPH solution and the experimental data is observed in terms of loading process and peak maxima. The phase shift in the pressure peak is due to the different peak propagation velocity which can be attributed in this case to 3D effects.

D. Flat panel impact : two-phases approach
Simulations involving 6 refinement levels (Δx min = 1.5625 mm) and 8 refinement levels (Δx min = 0.390625 mm) were performed also for the two-phase simulations (Fig.  17). At initialization, the total number of particles are about 125, 000 and 300, 000 for the simulations involving, respectively, 6 and 8 refinement levels. In order to take into account air-cushion effects, the bottom left corner of the flat plate is located 3 mm above the interface at initial time. In Figure 19, SPH pressure time histories for two different probes are compared with the experimental data in [19], [23]. The origin of the time axis is based on the pressure peak registered at P4. Since differences between single-phase and two-phases approaches are small, one can conclude that the influence of air is negligible for this deadrise angle.