Pressure evaluation at arbitrary locations in SPH water impact simulations

961 III International Conference on Particle-based Methods – Fundamentals and Applications PARTICLES 2013 M. Bischoff, E. O˜nate, D.R.J. Owen, E. Ramm & P. Wriggers (Eds) PRESSURE EVALUATION AT ARBITRARY LOCATIONS IN SPH WATER IMPACT SIMULATIONS MARTIN SIEMANN1AND PAUL GROENENBOOM2 1German Aerospace Center Pfaffenwaldring 38-40 70569 Stuttgart, Germany [email protected] – http://www.dlr.de 2ESI Group Netherlands Rotterdamseweg 183 C 2629 HD Delft, The Netherlands [email protected] – http://www.esi-group.com Key words: Smoothed Particle Hydrodynamics, Fluid-Structure Interaction, Pressure evaluation, Fixed-Wing Aircraft Ditching, Hydrodynamic Phenomena Abstract. This paper reports the application of the SPH method to the simulation of structures impacting on water with the focus on pressure assessment at arbitrary positions. The study is motivated by the importance of correct and reliable pressure evaluation in aircraft ditching simulations. The presented approach refers to recent developments in the SPH solver of VPS/PAM-CRASH which now allows using dummy particles to probe pressures. Sensitivity studies towards SPH parameters on pressure results are conducted. Pressures near contact interfaces are also compared to the local contact force divided by the surface area. Furthermore, pressure correction methods referred to as Shepard filtering and Rusanov flux, the lowest order approximation of a Riemann solver, are applied and tested. A comparison with results from flat plate ditching experiments published by Smiley [1] is done. Initial validation is based upon results of an experimental campaign of rigid wedge impacts conducted by Battley et al. [2]. Recommendations for best practice are derived from the study and the application to aircraft ditching is discussed. 1 INTRODUCTION The prediction of global and local structural loads is of fundamental importance in water impact problems, e. g. aircraft ditching. These loads significantly differ from those in a crash on solid ground [3]. One key value affecting the structural response is the pressure acting along the structure. It may affect the global kinematic response of an aircraft, therefore leading to catastrophic failure of the aircraft accompanied by fatalities. 1 Pressure evaluation at arbitrary locations in SPH water impact simulations 962 Martin Siemann and Paul Groenenboom Ditching refers to an aircraft emergency situation which ends with the planned impact on water. A ditching event is typically described in four consecutive phases: approach, impact, landing and floatation. However, the present work considers exclusively the impact phase (figure 1). The presence of a relatively high forward velocity in fixed-wing aircraft ditching affects the pressure distribution and the interrelated hydrodynamic effects acting on the fuselage. Consequently, the pressure distribution highly influences structural loads and, as a result, the global aircraft kinematics which may determine the survivability of such an emergency situation. The analysis of ditching impact is therefore necessary to satisfy the airworthiness regulations in certification of novel aircraft. Readers are referred to publications of Toso [3], Ben´ıtez et al. [4], Climent et al. [5], and Lindenau et al. [6] for a more detailed view on all ditching phases. Figure 1: Impact phase of fixed-wing aircraft ditching with emphasis on the pressure distribution along the rear fuselage due to hydrodynamic phenomena, e. g. cavitation, ventilation, overpressure, suction, air entrapment, and air cushioning. vx vz overpressure suction In the present work, pressure assessment at arbitrary locations along moving structures is studied using coupled SPH-FE simulation models. Attention is put on improving the well-known discrepancy of numerical noise in SPH pressure results. Results of this work aim to allow studying the above mentioned hydrodynamic phenomena in more detail. This may contribute to the development of improved fluid and fluid-structure interaction models for aircraft ditching analysis as currently these are not satisfactory mastered by existing simulation tools. 2 FLUID-STRUCTURE INTERACTION: COUPLED SPH-FE MODELS The presented numerical simulations are based on fully coupled SPH and explicit FE simulation models which—among the variety of numerical simulation tools used to assess fluid-structure interaction—offer a convenient analysis tool. The commercial explicit Finite Element software VPS/PAM-CRASH (ESI Group) with Lagrangian formulation is used. The structural models are composed of rigid FE shell elements. Their movement is prescribed by an initial velocity vector and it is constrained by limiting the degrees of freedom to model the guidance of the structure according to the respective experiments. 2 963 Martin Siemann and Paul Groenenboom The fluid domain is modeled using the weakly compressible SPH method. Particles are initially arranged on a cubic lattice configuration with particle spacing s(orthogonally spaced) and are bounded in a box of rigid shell elements or hydrodynamic solid elements. The free surface does not necessitate boundary conditions. For two-dimensional cases, inplane displacement is constrained. The Wendland kernel function (5th degree class 2) with a radius of influence of twice the smoothing length his applied. This kernel function was shown to be superior over the renormalized Gaussian kernel in free-surface flow problems by Maci`a et al. [7]. The numerical stability is treated by using the standard MonaghanGingold artificial viscosity term [8]. The constitutive equation (1) relating fluid density ρto pressure prefers to the Murnaghan (Tait) equation of state. It defines the pressure p(ρ) as p(ρ)=p0+c2 0ρ0 γ  B ρ ρ0γ −1(1) wherein p0is the reference pressure, c0is the speed of sound in the fluid at the state ρ=ρ0,Bis the bulk modulus, ρ/ρ0is the ratio of current over initial mass density and γis the adiabatic exponent of the fluid. The Murnaghan equation of state allows representing a fluid with artificially increased compressibility. This approach is feasible for fluid-structure interaction problems where flow velocities uremain well below the corresponding speed of sound and hence compressibility effects are insignificant [9, 10]. Satisfying the equation c0≥10 max(u), this criteria may be expressed as B≥100 ρ0max(u)2 γ.(2) Assuming typical aircraft ditching conditions, the flow velocities are much lower than the true speed of sound in water. This condition allows using a reduced speed of sound which (automatically) increases the critical time step of the SPH particles (governing the simulation’s critical time step) and, hence, decreases the run time of the simulation. For the studies in this work, the bulk modulus is chosen depending on the corresponding maximum flow velocity in the regarded test cases. As compressibility is essential for impact phenomena, the influence of the bulk modulus on pressure results was verified and found to be negligible within the investigated range of maximum flow velocities. Coupling the SPH and the FE model is achieved by using a node-to-segment penalty contact formulation between the particles and the impacting FE structure. Contrary to results published by Aquelet and Souli [11], within this work, contact damping (stiffness proportional) was found to not influence the numerical noise in pressure results. Further information regarding the penalty contact algorithm of VPS/PAM-CRASH may be found in [12]. 3 964 Martin Siemann and Paul Groenenboom 3 SPH PRESSURE CORRECTION METHODS The standard weakly compressible SPH method is well known to give poor pressure distributions in terms of high-frequency oscillations (numerical noise) in time and space. This deficiency is counteracted by pressure correction methods which aim to yield a more regular pressure distribution. However, the effect of correction should be local and not smooth global flow characteristics. Conservation of mass, momentum, and energy should be maintained or changes in the conservation should remain very small. Furthermore, correction methods must be numerically stable and may not require considerably higher amounts of computational power. Established pressure correction methods like density re-initialization by Shepard filtering and Rusanov flux were recently implemented in the SPH solver of VPS/PAM-CRASH. In this section, the fundaments of the respective correction methods are presented and their superiority on pressure results is demonstrated in figure 2. 3.1 Density Re-Initialization using Shepard Filtering One method which may reduce numerical noise in SPH pressure field computation is referred to as the Shepard filter. This density re-initialization method was derived from an interpolation technique initially published by Shepard [13]. As the SPH pressure calculation is based on an equation of state including the fluid density, re-initializing the density by Shepard filtering directly influences the pressure distribution. The implemented SPH notation for the modified density ρireads ρi= N  j mj Wij N k mk ρkWjk (3) where mjand ρjare mass and respective mass density of particle jand Wij is the kernel function. The density field is periodically re-initialized at a user-defined cycle frequency fwith recommended values of 20 [14]. In general, numerical noise originating from the weakly compressible SPH solution may be reduced and the pressure field therefore may become much smoother. The additional computational cost of this correction method is negligible. 3.2 Rusanov Flux: A Non-Conservative Riemann Solution The Rusanov flux is an efficient and robust, but also more diffusive, numerical scheme to solve Riemann problems. This approximation of a Riemann solver achieves only first order accuracy compared to second order accuracy of the Riemann flux. The used formulation was proposed by Parshikov et al. [15,16] and later by Cha and Whitworth [17] as Godunov Particle Hydrodynamics (GPH). In 1D, the continuity equation is modified to dρ dt =−2 N  j mj(u∗ ij −ui)∇iWij (4) 4 965 Martin Siemann and Paul Groenenboom in which the intermediate velocity u∗ ij is the acoustic solution (5). u∗ ij =ρiciuR j+ρjcjuR i+Pi−Pj ρici+ρjcj (5) Assuming that variations in density and sound speed remain small, and that the time step obeys the Courant criterion, the continuity equation becomes ρn+1 −ρn ∆t= N  j mj(un i−un j)·∇iWij −2 N  j mj ρj (Pn i−Pn j)rij ·∇iWij (r2 ij +χh2)∆t(6) with the strength parameter in the order of one-half [23]. Further, the solution of the acoustic approximation for the velocities adds an expression similar to the regular artificial viscosity in the momentum equation. As previously shown by Groenenboom [23], this correction method does not significantly influence forces and kinematics of the impacting structure. The additional computational cost is negligible while in practice run times decrease slightly due to reduced numerical noise in the density field. 0 5 10 15 20 −4 −3 −2 −1 0 0 5 10 15 20 −4 −3 −2 −1 0 z/b 0 5 10 15 20 −4 −3 −2 −1 0 x/b Figure 2: Pressure field at t= 30 ms in 2D NACA flat plate test case using no correction (top), Shepard filtering with cycle frequency f= 20 (center), and Rusanov flux with strength =0.5 (bottom). 5 966 Martin Siemann and Paul Groenenboom 4 PRESSURE EVALUATION METHODS 4.1 SPH Pressure Gauge Particles VPS/PAM-CRASH now allows using dummy SPH particles to probe pressures. Main advantages are that pressure information may be obtained at any desired location (fixed or moving in space) and, since the pressure is averaged over a number of nearby particles, it will suffer less from the strong spatial fluctuations typical for SPH pressures. Pressure gauge particles are defined as a special type of particle and may be described as passive particles which probe the properties of nearby regular particles but do not contribute to the evaluation of their SPH properties. To facilitate pressure assessment, gauge particles may be attached to moving structures or be put at arbitrary locations. Refer to figure 3 for a detailed view. The gauge smoothing length hgis the main numerical parameter to be investigated and calibrated. It determines the amount of regular particles contributing to the pressure evaluation of the gauge and will clearly influence the stability in case of rather small values. Furthermore, it is expected that pressures assessed with gauge particles will be subject to smoothing in case the gauge smoothing length is chosen too large. In order to generalize findings, the gauge smoothing length is chosen proportional to that of the regular particles which is defined as gauge size ratio Π = hg/h. 4.2 Pressure via Contact Normal Force over Contact Area Another, to date frequently used method assessing pressures in numerical simulations is given by calculating the contact normal force over contact area. Here, pressure time histories depend on the ratio of regular particles per shell element and, furthermore, on the size of the contact area Acont. The latter may cause averaging over too many shell elements (too large area) which may smear out pressure peaks (local information) from the pressure time history. Figure 3 provides a schematic illustration. Figure 3: Detailed view on pressure evaluation methods: 2Dview showing position and size of gauge particle P(left), and pressure assessment via contact normal force over contact area in 3Dcases (right). hg P shell plane plate surface fluid particles vx vz shell plane Acont vx vz 6 967 Martin Siemann and Paul Groenenboom 5 TEST CASES 5.1 Two-Dimensional Rigid Wedge Vertical Impact This first test case was taken from Battley et al. [2], who performed motion-controlled vertical impact experiments using a rigid wedge of 10◦deadrise angle impacting with constant velocity of 3 m/s (servo-hydraulic controlled). Experimental pressure results are available for three positions along the center line between keel and chine which has a total length of 600 mm. Further information about the instrumentation, the data acquisition, and the servo-hydraulic slam testing system may be found in [2]. The numerical model is symmetrical with respect to the vertical axis. It consists of a wedge structural model and a two-dimensional water domain filled with particles with smoothing length of h=2mm. Pressure gauge particles with a smoothing length of hg=4mm (Π = 2) are positioned on the actual surface location of the wedge. Numerical pressure results using the previously described pressure gauge particles are validated against experimental data in figure 4. Numerical pressure time histories are CFC1000-filtered. Comparison between experimental and numerical pressure time histories shows good correlation of peak values. However, after peaks diminish, the numerical residual pressures are higher compared to the experimental ones. It is believed that the general overestimation of pressures in the simulation is caused by the two-dimensional nature of the numerical model. Nevertheless, it has been demonstrated that the pressure gauge particles allow for reasonable pressure results for this test case. 0 5 10 15 20 25 30 35 0 50 100 150 200 250 300 Time tin ms Pressure pin kPa exp sim Figure 4: Comparison of experimental and numerical pressure results: Time histories at positions P1–P3 and pressure contour plot at t= 20 ms. P1 P2 P3 0 20 40 60 80 100 [kPa] symmetry wedge, β= 10◦,v z=3 m s 7 968 Martin Siemann and Paul Groenenboom 5.2 Two-Dimensional NACA Flat Plate Ditching Due to lack of experimental data to validate the novel pressure assessment method using gauge particles under ditching conditions, flat plate ditching experiments published by Smiley [1] in 1951 are selected as a test case for this numerical parameter study. Nevertheless, initial velocities are smaller compared to fixed-wing aircraft ditching characteristics. The selected test configuration (run 4 in [1]) consists of a flat plate (1524 mm x 304.8mm x 18.288 mm, pitch angle α=6 ◦), which is connected to a trolley moving along a guidance structure during the entire experiment with a prescribed motion in terms of initial velocities in horizontal and vertical direction (vx= 13.35 m/s and vz=1.77 m/s). 18 flush mounted pressure gauges with 12.70 mm diameter and one bellow-type pressure gauge with 6.35 mm diameter were used. Data was acquired at a rate of 1 kHz. However, only eleven time instances are available for the selected case. As this is insufficient for a comparison of pressure time histories, the parameter study compares numerical pressure results to experimental maximum pressures along the center line which vary between 300 kPa and 400 kPa. One pressure gauge measured a higher maximum pressure of 427 kPa which is believed to be related to its smaller size (diameter 6.35 mm). The numerical model consists of a flat plate structural model and a two-dimensional water domain (see figure 2). Multiple pressure gauge particles are attached to the structure at the actual surface location of the plate. Within an extensive parameter study, gauge size ratios ranging between Π = 1 – 10 were investigated using different global finesse of b/s = 10 – 100 which refers to the particle spacing snormalized on the width bof the plate. Below results refer to a particle spacing of s=6.1mm (corresponding to b/s = 50) which gives reasonable results. Pressure data are compared at ξ=0.4(ξ∈[0,1] is the relative local plate coordinate in longitudinal direction measured from the trailing edge). Numerical pressure results assessed by pressure gauge particles are compared to maximum pressures from the respective experiments. Similar to observations in the experiments, numerical results show a sharp and immediate rise of the pressure during impact. The pressure peak travels with the ditching front along the plate until its full submersion where pressures become significantly lower. Pressure results of this study are analyzed in figure 5. Regarded pressure correction methods influence the oscillations which are smallest when using the Rusanov flux. Filtering of pressure signals is a critical issue as it may significantly alter peak pressure values. Hence the influence of filtering was studied for different gauge size ratios revealing that filtering reduces pressure peak values especially for ratios below Π = 3. Furthermore, there is a significant influence of the gauge size ratio on the pressure time histories. With increasing gauge size ratio, pressure peaks as well as pressure gradients and associated oscillations are reduced (smoothing effect), but the later residual pressures remain at similar level of magnitude. Based on the presented study, it is recommended to chose the gauge size ratio as small as possible but no less than Π = 2. 8 969 Martin Siemann and Paul Groenenboom 10 20 30 40 50 60 70 80 0 100 200 300 400 500 Time tin ms Pressure pin kPa uncorrected Shepard f = 20 Rusanov  =0.5 10 20 30 40 50 60 70 80 0 100 200 300 400 500 Time tin ms Pressure pin kPa Π=1.0 Π=1.5 Π=2.0 Π=2.5 Π=3.0 12345678910 0 100 200 300 400 500 600 700 800 900 Gauge Size Ratio Π in − Peak Pressure pmax in kPa CFC1000 CFC600 CFC180 CFC60 10 20 30 40 50 60 70 80 0 100 200 300 400 500 Time tin ms Pressure pin kPa Π=4.0 Π=5.0 Π=6.0 Π=8.0 Π = 10.0 Figure 5: Numerical pressure time histories at position ξ=0.4 comparing the influence of pressure correction methods (top left), gauge size ratio Π (top and bottom right), and filtering (bottom left). 5.3 Three-Dimensional Guided Ditching Test The novel pressure assessment method using gauge particles was applied to the threedimensional numerical model of the guided ditching tests to be carried out within the EC-funded research project SMAES (SMart Aircraft in Emergency Situations) [4,18–20]. The guided ditching simulation model is shown on the right hand side of figure 6. For the selected load case, the initial velocities are vx= 50.0m/s in horizontal and vz=1.5m/s in vertical direction. The flat plate measures 1000 mm x 500 mm x 15 mm and is assumed to be rigid. Throughout the simulation, the pitch angle of 10◦remains constant due to the guided motion. Fluid particles have a smoothing length of h= 10 mm and gauge particles of hg= 20 mm (Π = 2) respectively. The contact surface used for the pressure evaluation has an area of Acont = 1600 mm2. Results are compared at the position η=0 (center line regarding the lateral direction) and ξ=0.25. 9