Free Surface Impact in a Biphasic SPH Simulation

Free Surface Impact in a Biphasic SPH Simulation G. OGER, P. FERRANT, B. ALESSANDRINI Fluid Mechanics Laboratory (CNRS UMR6598), Ecole Centrale de Nantes [email protected] Abstract In most of cases, free surface impacts may be correctly described without taking care of air effects. But depending on the shape and dimensions of the solid involved in the impact, these effects may have some non negligible influences on the impact features at the entry beginning, as it would be observed in the occurrence of air-cushioning. In this paper, through a two-dimensional numerical simulation of a body in free motion within a multi-phase flow using Smoothed Particles Hydrodynamics (SPH) method, these influences are taken into account and discussed. The application presented here deals with fluids having a high density difference, implying a need for some special care on the SPH formulation to adopt. Successive effects of air and water on the body are included, and numerical results are validated through some comparisons with experimental data, showing good agreements. Introduction The meshless and lagrangian properties of SPH allow to solve easily some problems involving one or several solids acting with very large imposed motions through a fluid flow. Recent works [3][5] showed the possibility for this method to be enhanced with a specific treatment to deal with free solids in a weakly coupled approach. This has been achieved through a procedure consisting in calculating forces on the solid from fluid particles pressure. We aim now at simulating cases where air effects play a crucial role in the global fluid (or fluid-solid) behavior, as it would be observed with air-cushion effects or air entrapment in breaking waves for instance. The standard SPH formulation [4] naturally allows multi-phase flow simulations with a weak density difference between the various fluids considered. But for high density differences, this formulation introduces some strong instabilities at the interface, which ends in contaminating the whole simulation. Some recent works [2] show that another formulation is needed to cope with this problem. The study presented here combines the treatment of free moving boundaries with an air-water implementation through the standard test case of a free falling wedge in the air, ended with its impact on the free surface. Successive effects of both air and water on this prism are attempted to be captured as accurately as possible. First a global presentation of SPH is provided, followed by practical computational considerations. Finally, the simulation of this prism in air and water is discussed, and comparisons with experimental data are provided. SPH method for compressible fluid flows SPH scheme is based on weighted interpolating points discretizing the fluid domain. Using a kernel function in a convolution, these points can be used to determine the gradients appearing in partial differential equations of conservation laws without any underlying mesh. For free surface flows, the system of equations we need to solve is the classical Navier-Stokes one. Since we always assume the fluid to be inviscid in the application presented here, shear stresses are neglected and this system reduces to Euler equations as follows. d~v dt dρ dt P = ~g − ~ ∇P ρ (1) ~v = −ρ.∇~   γ ρ = κ −1 ρ0 (2) (3) W (q =       |~r| )=C  h     2 3 − q 2 + 12 q 3 1 6 (2 0 − q)3 if 0 ≤ q < 1 if 1 ≤ q < 2 (4) else x This scheme is written in a Lagrangian way, so that particles follow the fluid movement as d~ v . One dt = ~ of the main SPH features consists in considering any fluid media as compressible, resulting in the use of an equation of state called Tait’s equation (3), linking the pressure to the density and thus allowing the above system of equations to be closed. Variables are convoluted using a kernel function throughout a set of interpolating points as follows Z f (~r) ≈ f (~x)W (~r − ~x, h)d~x (5) D where h is a spatial discretisation parameter and D a circular domain (in two dimensional context) defined with a radius equal to 2h. This approximation makes sense provided the kernel function W tends to a delta function in the limit h → 0, making the approximation (5) to become a strict equality. The kernel function used in this paper is defined by equation (4). Note that this function has a compact support since its value is zero everywhere out of D. Furthermore, the gradient of any function f can be approximated by extending (5) writing equation (6), which reduces to equation (7) thanks to an integration by parts and by neglecting the surface term using the compact support property of W. ~ (~r) ≈ ∇f Z ~ (~x)W (~r − ~x)d~x (6) ∇f ~ (~r) ≈ ∇f D Z ~ (~r − ~x)d~x f (~x)∇W (7) D In order to enhance the numerical efficiency, the formulae are symmetrized [4] following equations (8) and (9). Writing this finally leads to the formulation defined by equations (10) to (12), which seems to be the most commonly used by SPH practitioners. ~ ∇P ρ ~v ρ∇~ ~ = ∇  P ρ  P ~ + 2 ∇ρ ρ ~ v ) − ~v ∇ρ ~ = ∇(ρ~ (8) (9) dx~i dt d~ vi dt dρi dt = v~i = ~g − (10) X mj j = X Pi Pj + 2 ρ2i ρj ! ~ (~ ∇W ri − r~j , h) ~ (~ mj (~ vi − v~j ) ∇W ri − r~j , h) (11) (12) j In the context of multi-phase flows, the various fluids simulated with SPH can be seen as only one fluid with a sudden density jumps at interfaces. No specific treatment such as dynamic boundary condition is implemented at these interfaces. Furthermore SPH implicitly fulfills any kinematic boundary condition thanks to its Lagrangian feature. Note that SPH method described with equations (11) and (12) is not suitable with multi-phase flows computations except maybe for the case of small density ratios. In deed, for fluids having large density differences as it would be the case for a air-water computation for instance, some strong instabilities ~ in the symmetrization equations appear, leading to a poor simulation. This fact is due to the presence of ∇ρ (8) and (9). Actually, because of the density discontinuities at fluid interfaces, the density gradient becomes singular and its interpolation results in a sum of errors. ~ does not appear while A solution to overcome this problem is to change the global formulation so that ∇ρ interactions remain symmetrized [2]. This can be done writing equations (13) and (14), leading to the alternative formulation equations (15) to (17). ~ ∇P ~v ∇~ ~ + P ∇1 ~ = ∇P ~ v − ~v ∇1 ~ = ∇~ (13) (14) dx~i dt d~ vi dt dρi dt = v~i (15) = ~g − X j = ρi mj (Pi + Pj ) ~ ∇W (~ ri − r~j , h) ρi ρj X mj ~ (~ (~ vi − v~j ) ∇W ri − r~j , h) ρ j j (16) (17) ~ terms would The choice of the signs in (13) and (14) is not naı̈ve. Indeed, changing the signs of the ∇1 still be mathematically correct, but for strict consistency considerations, the choice for the SPH formulation of (2) imposes the formulation to adopt for (1), as can be found in [1]. Finally, this ordinary differential equation system is integrated in time by schemes such as Runge-Kutta, Leap-Frog or Predictor-Corrector. The case treated here concerns the capture of a free body dynamics in coupled interactions with two various fluids. For this we used the procedure described in [3] to extract forces due to fluid(s) on its boundaries. Test case of a free falling wedge impact Experimental device These technics is applied to the standard validation test case of a free-falling wedge impacting the free surface. The experimental device is described in figure 1. Some more details on this experiment can be found in [6]. Figure 1: Test section description At t=0 s, this free-falling wedge is dropped from 0.61 meters above the free surface with a five degrees clockwise initial heel angle and no initial velocity. This test section was instrumented with a pair of triaxial piezoresistive accelerometers dedicated to the measurement of acceleration time histories. This pair of accelerometers were separated transversely in order to measure the roll dynamics in addition to vertical dynamics. The free fall generates two whirlwinds above the prism, and the impact on water creates two free surface jets blowing the air vortices during their exit, and imposing a strong vertical deceleration as well as a transverse self-righting of the solid. Experiments are supposed to be realized so that the flow could be regarded as two-dimensional. Numerical simulation In the SPH simulation, the tank size has been chosen to ensure no interaction between the wedge and the sound wave generated by the impact. The nominal sound speed chosen for water is 80 m/s. Because of its nearly incompressible feature, the equation of state (3) used for water is computed with γ = 7, whereas knowing its natural adiabatic behavior γ = 1.4 for air. The hydrostatic configuration of air makes the pressure differs from 0 at the interface, this pressure having to be the same for the two fluid in this critical zone. The sound speed chosen for air is 150 m/s. This value is sufficient to preserve the weakly compressible aspect of air flow in the circumstance of M a < 0.1. Note that it would be possible to use the actual air sound speed, that is 340 m/s, but such a value would require some smaller time steps with a similar final result in the application presented here. ρ0 = 1000kg/m3 was used as the nominal density for water, and ρ0 = 1.29kg/m3 for air. A variable smoothing length computation has been achieved in order to concentrate the accuracy in the free fall area in the air, and in the impact area in water, limiting the total particle number needed and the associated CPU time. Finally 230,000 air and water particles were necessary to achieve this simulation. This test case has already been successfully simulated using a monophasic SPH computation [3][5]. However, the very beginning of vertical acceleration time history turned out to be rather subtle: the use of a quite crude simulation resolution gave the illusion of a very good capture of slamming force at the impact beginning, but this property disappeared while we refined our particle distribution, revealing a time difference of order 10−3 s between SPH and experimental results, and imposing the need for matching curves in time to make the comparison possible. We supposed this phenomena to be due to air influence, by the mean of air cushion effects that could not be captured by a monophasic calculation. We tried here to capture this effect. Figure 2 shows a comparison between SPH and experiments for vertical and angular acceleration time history. 13 5 12 10 9 8 ÿ g 7 6 SPH Experiment 4 θ̈ (103 deg.s−2 ) SPH Experiment 11 3 2 1 5 0 4 3 -1 2 1 -2 0 -1 0.35 0.375 0.4 t(s) 0.425 -3 0.35 0.375 0.4 0.425 t(s) Figure 2: Vertical and angular accelerations results As shown in the left plot of this figure, the vertical acceleration time history of the wedge is correctly predicted when compared with experimental data. Maximum load appearing at t=0.365 s is well captured in time as well as in amplitude, and the main trends of loads during water entry are well approached. The same comment can be made concerning the angular acceleration time history in the right plot, namely the three peaks are captured in amplitude with a good agreement with experimental data. Nevertheless, it should be underlined that a defect appear from t=0.375 s, where a stiffer slope than in the experiment is obtained in both vertical and angular time evolutions. It mainly leads to a time difference in the third peak of the angular acceleration. This behavior finds its origin in the presence of instabilities appearing at the right corner of the wedge during the exit of the jet, due to some high velocity gradients between the two phases in this sharp angle area. 1 0.5 Biphasic SPH Monophasic SPH ÿ g Experiment 0 -0.5 -1 0.34 0.345 0.35 0.355 0.36 t(s) Figure 3: Numerical air cushion effects In a first attempt this biphasic computation results in a globally correct approach of the main phenomena, and allows to focus on the wedge vertical acceleration evolution at the impact beginning. Right plot of figure 3 emphasizes the overpressure due to air (color pressure scaled) expulsion between the wedge and the free surface (black points) at instant t=0.35 s. Left plot of figure 3 shows a zoom of the left plot of figure 2 around this instant. It compares the experimental data with the monophasic and biphasic SPH results, for the same used smoothing length, at the water impact ignition. This plot clearly displays the air influence on numerical results. Nevertheless, it appears that experimental data displays a more progressive evolution than the SPH biphasic one, and that numerical air effects are overestimated. Note that the shape of the test section used in the experiment does not have such sharp knuckle angles, and may lead to a quite different aerodynamic behavior. Further work is still needed in this way. Conclusion The free fall in the air of a low deadrise angled wedge and its impact on free surface has been simulated using a SPH based method, the whole dynamics of this solid being obtained by pressure integration of both air and water pressure forces. Results concerning vertical and angular acceleration time histories have been compared with experimental data, showing good agreements in time synchronization as well as in amplitude. On the ground of this biphasic implementation, air cushion effects at the impact ignition was studied, emphasizing a vertical force that can not be captured by a monophasic SPH computation, but giving a less good agreement when compared with the experiment. Improvement concerning instability problems at media interfaces and the validation of air cushion loads capture need to be pursued. References [1] J. Bonet and T. S. L. Lok. Variational and momentum preservation aspects of smooth particle hydrodynamics formulations. Comput. Methods Appl. Engrg, 180:97–115, 1999. [2] A. Colagrossi and M. Landrini. Numerical simulation of interfacial flows by smoothed particle hydrodynamics. Journal of Computational Physics, 191:448–475, 2003. [3] M. Doring, G. Oger, P. Ferrant, and B. Alessandrini. Sph simulations of floating bodies in waves. In Proceedings of the 19th International Workshop on Water Waves and Floating Bodies. April 2004. [4] J. J. Monaghan. Smoothed particle hydrodynamics. Annu. Rev. Astron. Astrophys., 30:543–574, 1992. [5] G. Oger, M. Doring, L. A., P. Ferrant, and B. Alessandrini. Sph and finite volume simulations of a wedge water entry. In Proceedings of ISOPE’04 . May 2004. [6] L. Xu, A. W. Troesch, and R. Petterson. Asymmetric hydrodynamic impact and dynamic response of vessels. Proceedings of 17th International Conference on OFFSHORE MECHANICS and ARCTIC ENGINEERING, pages 98–320, 1998.