Fluid-structure interaction by the mixed SPH-FE method with application to aircraft ditching

This paper deals with numerical simulation of fluid-structure interaction as it occurs during aircraft ditching – an emergency condition where an aircraft is forced to land on water. The work is motivated by the requirement for aircraft manufactures to analyze ditching as part of the aircraft certification process requested by airworthiness authorities. The strong interaction of highly non-linear fluid flow phenomena and structural responses requires a coupled solution of this transient problem. Therefore, an approach coupling Smoothed Particle Hydrodynamics and the Finite Element method within the commercial, explicit software Virtual Performance Solutions has been pursued. In this paper, several innovative features are presented, which allow for accurate and efficient solution. Finally, exemplary numerical results are successfully compared to experimental data from a unique test campaign of guided ditching tests at quasi-full scale impact conditions. It may be concluded that through the application of state-of-the-art numerical techniques it has become possible to simulate the coupled fluidstructure interaction as occurring during ditching. Therefore, aircraft manufacturers may significantly benefit from numerical analysis for design and certification purposes.


INTRODUCTION
The present work is motivated by the objective to simulate the complex water entry of aircraft during an emergency landing on water.Predicting hydrodynamic loads and the associated structural response is of high interest for the analysis of aircraft ditching, but also in the maritime field where lightweight structures are increasingly being used.
Ditching must be analyzed by aircraft manufacturers in order to fulfil certification requirements [1,2].Current practices involve sub-scale experimental testing and comparison with aircraft that previously passed ditching certification.In addition, semi-analytical methods are used to simulate the global motion of the rigid body aircraft during the ditching.The structural capacity is subsequently evaluated by applying worst case pressure distributions that are derived from hydrodynamic loads acting on the rigid body; hence the fluid-structure interaction (FSI) is uncoupled.Occurring structural deformations, however, dynamically *Corresponding Author: E-mail: paul.groenenboom@esi-group.nl change the boundary conditions for the fluid during ditching, potentially affecting the hydrodynamic loading, and should therefore be taken into account within the analysis.
In order to simulate the fully coupled FSI, an explicit numerical simulation approach combining Smoothed Particle Hydrodynamics (SPH) for the fluid and Finite Element (FE) method for the structure is selected.The SPH method is particularly suited for this application as it naturally handles large fluid deformations and large displacements of both the interface with structures and free surfaces.The FE method is the most appropriate tool to model structures under dynamic loading.Both methods are Lagrangian which significantly facilitates their coupling.Fundamentals of the simulation approach are addressed in section 2.
Arising from the relatively large forward velocity of the aircraft during ditching and the duration of the event, large fluid domains are typically required.Accurate resolution of the fluid-structure interaction simulation, on the other hand, requires a locally fine fluid model which implies that the computational demands will become very high if not prohibitive.Additionally, the quality of the numerical solution has to be maintained at a high level.
In order to overcome some of the limitations, the SPH solver within the commercial explicit FE analysis code Virtual Performance Solutions (VPS, formerly known as PAM-CRASH) has been significantly enhanced over the last years.Novel features comprising an improved pressure evaluation, particle position regularization, translating periodic boundary conditions and the option for non-uniform particle distributions are introduced in section 3.
In section 4, improvements are applied to a model of guided ditching tests.These experiments, comprising oblique water impact of aeronautical structures at quasi-full scale impact conditions, were recently conducted by CNR-INSEAN in order to support the development and validation of numerical tools.Herein, exemplary numerical results are compared to experimental data to outline the capabilities of the numerical simulation, yet also its limitations are discussed.

COMPUTATIONAL APPROACH
The numerical results presented in this publication have been performed with the commercial software tool VPS [3].VPS is a general purpose finite element code with an explicit solver optimized for dynamic, strongly non-linear structural mechanics.The software contains finite element formulations for thin shells, solid elements, membranes and beams with material models allowing for plasticity and failure of metals, plastics, rubbers, foams and composites.Furthermore, VPS incorporates an SPH solver.The software also contains robust contact algorithms to enable various parts within a model to interact.Optimal numerical performance is approached using the Distributed Memory Parallel (DMP) version.

SPH METHOD
The SPH method is an interpolation method in which each particle describes a fixed amount of material in a Lagrangian reference frame.Despite its appearance as a collection of points, SPH is aimed at solving continuum mechanics of fluids and solid materials.Since the method does not require a mesh, SPH allows simple handling of the motion and topology changes of material surfaces.Evaluation of relevant field variables such as density and velocity at the location of a selected particle is conducted by interpolation over all neighbor particles in a region of influence as illustrated in Figure 1.The size of the sphere (or circle in 2D) of this region, is characterized by a smoothing length h.
Unless mentioned otherwise, the SPH method discussed here is basically similar to that of Monaghan [4,5].Within the SPH method, the derivative of the density for particle i may be solved from the continuity equation: (1) where the sum extends over all neighbor particles and ∇W is the gradient of the smoothing kernel evaluated at the distance between particles i and j.Optionally, renormalization of the density field may be performed [6].The updated velocity may be obtained from the momentum equation: (2) in which an artificial viscosity is defined by: (3) with: (4) The values α and β determine the strength of the artificial viscosity required to suppress numerical shocks.Values which are as low as possible should be used to avoid the flow becoming too viscous.The corrective constant η is used in (4) to avoid creating a singularity when particles are approaching one another.The choice for the smoothing kernel is a quintic Wendland kernel (5) which is compact within 2h.
(5) In Equation 5, q is the ratio of smoothing length h and distance between neighbors r, and θ is a normalization constant depending on the problem dimension.The flow will be assumed to be nearly incompressible implying that the pressure field is obtained from an equation of state (EOS) model.In this work the Murnaghan model, which is also known as Tait model, is used: (6) with background pressure p 0 , speed of sound c 0 , density ρ 0 , and adiabatic coefficient γ.The background pressure must be zero in this application to avoid free surface dispersion.Additionally, a cut-off pressure limiting the minimum pressure may be included as a simple approximation to cavitation.
Particles are assumed to interact mutually only if they are sufficiently close to each other; this being established by a nearest neighbor search algorithm that, for efficiency reasons, does not need to be repeated each computational cycle.Second-order accurate leap-frog time stepping is used for the explicit time integration of the above rate equations.The numerical time step is set at a fraction of the well-known CFL criterion based on the current size and sound speed of all particles and finite elements present in the model.
Symmetry with respect to reflection in one or more of the Cartesian planes by means of ghost particles may be employed to reduce the size of some of the models, if appropriate.

COUPLING WITH STRUCTURES
Interaction between particles, representing a fluid, and finite elements, representing moving or deformable structures, may be modeled by one of the sliding interface contact algorithms available in VPS.Such algorithms prevent penetration between selected structures, with sliding allowed in most cases.These sliding interfaces are based on the well-known penalty formulation, where geometrical interpenetrations between so-called slave nodes and adjacent master faces are penalized by counteracting forces that are proportional to the penetration depth.The contact algorithm will automatically detect when a particle (slave) penetrates a segment (master) of the outer surface of the finite element model of the structure.The contact thickness indicates the distance away from a contact face where physical contact is established.For SPH particles as slaves, the contact thickness should be representative of the particle spacing, possibly augmented by one-half of the shell thickness of the structure.This type of contact has been validated by simulating the vertical motion of floating bodies.It has been found that the correct floating body position is reached with a resolution much smaller than the particle size when the thickness defined for the contact equals one half the particle spacing and the artificial viscosity coefficients are taken significantly smaller than the values normally applied for shocks.In that case the upward force is also correct.This contact may also be used to define rigid boundary conditions limiting the SPH domain.A tied contact may be used in case sliding does not need to be considered, such as when an artificial interface between particles and volumetric elements is used to represent a contiguous volume of fluid.Coupling of the SPH-FE through the contact interface has been validated for a range of cases [7].Industrial applications include sloshing [8] and the impact of aeronautical structures on water [1,9].

PRESSURE CORRECTION
One of the outstanding problems of the weakly compressible SPH solution is the rather large variation in pressure, both in time and space.The option of density renormalization (using the Shepard filter) has proven to be one method that allows reducing these variations, in particular for temporal variation, but the remaining variations in space are still quite large.
The mass diffusion option proposed by Molteni and Colagrossi [10] provides a significant reduction of these pressure variations.This formulation was later extended with a viscous term usually referred to as the δ-SPH scheme [11,12] and has been demonstrated to provide superior pressure results.
A slightly different formulation may be obtained by inspection of the time-discretized version of the momentum equation resulting in the following expression for the time rate of the density: (7) The most suitable value of parameter ξ could be determined from numerical tests.In case the pressure is related to the density and the time step is determined by the CFL criterion for the particle, and some density variations may be ignored, the correction will be similar to the original correction according to the mass diffusion proposed by Molteni and Colagrossi [10].
Furthermore, equation ( 7) may be demonstrated to be equivalent to the first-order approximation to the Riemann problem for fluid flow [13] (Rusanov flux), where the free parameter ξ is set to 0.5 and this will be the value adopted as the default for flow simulations.The pressure correction yields a significantly smoother pressure distribution than the reference SPH simulation but without significant effects on the kinematics of interfaces and on the velocity distribution.The total fluid volume (or average density) is also not modified.The additional computational cost is negligible while in practice run times decrease slightly due to reduced numerical noise in the density field.The successive reduction of the pressure oscillations in the water beneath an impacting plate by these improved SPH options is shown in the pressure contours and the time histories in Figure 2 for a typical, two-dimensional water entry test case [14].

PARTICLE REGULARIZATION
During a flow simulation an initially regular particle distribution may become irregular by creation of small void regions or by clumping.In such cases the accuracy of the SPH method may deteriorate rapidly.Particle regularization algorithms (PRA's) may be devised to reduce clumping by redistributing the particles at regular intervals throughout the analysis.Such an algorithm is expected to mitigate numerical effects such as the irregular pressure distribution for the standard SPH.Since the tension instability is related to the occurrence of irregular particle distributions, it is expected that a suitable PRA will alleviate this problem.If particles are defined to be in a different position than that obtained by direct integration of the individual particle kinematics, they no longer represent the same amount of material and a correction is required [15].This correction is not accounted for in the usual particle displacement correction methods (XSPH).
A proper redistribution algorithm should ensure that: 1.A more regular particle distribution is obtained within a limited number of calculation steps.2. The effect of the correction of a locally 'poor' distribution should remain local.3. The method should be numerically stable and not require excessive amounts of CPU. 4. The exterior surface should remain close to the original one.5.The updated mass, velocity, density and other properties should be defined properly.6.The total mass, momentum and energy should be conserved accurately.7. The modifications of the properties should be close to zero for initially smooth distributions.
Since the particle positions are redefined according to a particular type of PRA, the individual particle masses will no longer be conserved.With an appropriate interpolation method for the density and particle volume, the mass of any particle affected by the PRA of choice is corrected.
Various methods to displace the particles other than only by the flow field have appeared in the literature [15].One of the methods is the color gradient smoothing (CGS) which is based on an algorithm that has been proposed by Lind et al. [16] for incompressible SPH.There does not appear to be any reason as to why this could not be applied to weaklycompressible SPH as well.Lind et al. observed that the proposed shifting algorithm enhances both stability and accuracy of the SPH simulation for fluid flow.The basic concept of the shifting algorithm is to provide an additional displacement to the particles according to Fick's law for diffusion: (8) where is the flux, C the particle number density or color and Dʹ a diffusion coefficient.
The diffusion coefficient is in this case a numerical parameter for which the CFL stability criterion should be valid.Hence: (9) where h is the smoothing length, α a constant smaller than one and Δt the time step.Using the above relations, the displacement correction may be written as As also mentioned by Lind et al. [16], particles near the boundary should be exempt from this displacement correction since there will be a strong color gradient in the normal direction except when there would be several arrays of ghost particles.
Additional methods that may improve the particle distribution during the flow simulations have been tested but will not be discussed here.

INITIAL PARTICLE DISTRIBUTIONS
For many SPH simulations the possibility to start with a distribution of particles of nonuniform size may be of significant advantage since it allows for having small particles where needed and bigger particles where allowed.For that reason the option to allow particles to either fill a given domain or redistribute themselves according to a non-uniform spatial distribution for the particle size during the particle filling process is of great importance.
The Weighted Voronoi Tessellation (WVT) technique proposed in [17] is to start from a more or less arbitrary particle configuration and iterate towards a converged distribution.At each step, a displacement is imposed on particle i by: (11) in which λ is a strength parameter, h the smoothing length and r ij the neighbor distance.Diehl et al. [17] found optimal convergence when f decreases inversely proportional to the distance squared, so they propose: (12) In the above equations h ij is the average smoothing length and ε a small positive number to avoid divergence in case of particles being very close to each other.The original WVT algorithm was proposed in the context of astrophysics where boundaries are usually not of concern.For most other SPH applications, boundaries are highly relevant which led to extensions of the WVT algorithm.
The first extension of the WVT algorithm is the inclusion of boundaries either by a rectangular box or by finite elements.The second extension allows for an arbitrary boundary shape to be discretized into a set of plane segments, triangular or quadrilateral, that can be made to fit curved surfaces as accurate as required.This FE distribution may also be used to define the boundaries for subsequent SPH simulations.For the EWVT iteration, the displacements as defined by Eqns (11)(12) are limited to yield positions that are not closer to the boundary than a given search distance.It was observed that a superior distribution could be obtained by replacing the box by defining a suitable distribution of ghost particles on the opposite side of the wall.To do this properly, it is required defining these particles such that they resemble the reflection of the interior particles (adjacent to the wall locations) in the final configuration, which is a-priori unknown.For particles of non-uniform size it will not be possible to devise a proper distribution of ghost particles.Hence, an algorithm has been implemented to do this automatically.
A further extension of the WVT algorithm is the possibility to let particles become larger while they are displaced.For the standard WVT, the particles are not displaced any further when they have reached a spacing equal to 2h to their neighbors.By allowing the smoothing length to increase in proportion to the increasing particle diameter, the distribution will expand until the boundaries have been reached.The particle volume will increase according to h.This option allows filling domains of arbitrary shape and size starting from an initial particle distribution in the interior that can be kept very simple.A proper algorithm for increasing the size of the particles during expansion is by increasing h of all particles by the ratio of the particle density number (color) averaged over all particles for the updated configuration with that at the previous iteration step.
Although in a recent implementation using the WVT algorithm from Diehl at al. [17], the resulting particle distribution appeared to be quite reasonable [18], it was observed that the impact forces for water entry simulations were not always correct in case the particle size was not uniform, which was found to be related to a non-uniform particle number density within the domain.Hence, it was proposed to consider the color gradient smoothing from the particle regularization feature as discussed above (cf.section 3.2) as alternative.The CGS method is constructed in such a manner that it reduces the gradient of the color by displacing the particles and is, just as the WVT algorithm, a purely geometrical approach.The CGS method has been implemented as an option within the EWVT approach and is used as a subsequent, final step to generate stable, non-uniform distributions with a correct particle volume distribution.

PERIODIC BOUNDARIES
A method for reducing the CPU demands for flow simulations in extended or infinitely long domains is that of periodic boundaries.This feature allows particles leaving a user-defined rectangular domain to be entered at the other end with the same properties and is useful for a variety of flow studies.This type of boundary conditions has recently been extended to allow opposing boundaries to be translated according to the motion of a user-defined node.With this extension, the particle domain can follow the motion of a moving structure such as a ditching aircraft without introducing additional velocities for the particles themselves [19].A special feature that is beneficial for ditching simulation is the option to re-enter particles from the wake of the structure at undisturbed conditions but since particles on opposite sides of the periodic boundaries are considered as neighbors which influence each other, some disturbance to the reentering particles cannot be avoided.For that purpose, the translating periodic boundary condition may be combined with a damping zone feature [15] in which particles are returned to their initial conditions and positions (except for the direction in which they are translated).
Figure 3 illustrates how the combined features of the periodic boundary condition and the new damping zone may be used to significantly reduce the size of the domain.Many particles have left the domain on the 'downstream' side but despite the huge amount of spray-like distortion and the reduction of free surface level directly behind the structure, no disturbance is observable in the particles that have re-entered at the front face.With the domain translating at the same instantaneous horizontal velocity as the structure (linked by one reference node), the distance from the plate to the domain boundaries remains constant throughout the simulation.The test campaign comprises panels of different thickness, material and curvature which are accelerated to quasi-full-scale horizontal velocities similar to real aircraft ditching before they impact the water.All panels measure 1000 × 500 mm and are attached to a trolley that moves guided during the complete experiment in order to provide well-defined boundary conditions, which is beneficial for comparison with numerical results.The total mass of the structure impacting the water is about 832 to 840 kg depending on the test case.Key parameters varied include horizontal impact velocity (30-46 m/s), pitch angle (4-10 degrees), curvature (flat, convex, concave), material (aluminum, composite) and thickness (0.8-15 mm) of the panels, which allows for essential validation of numerical models.Measurements include up to 44 channels for accelerations, forces, strains and pressures.For more details on test campaign, test facility and initial results the reader is referred to [20].

NUMERICAL MODEL
The numerical modeling of the guided ditching test follows the SPH-FE approach presented in section 2 and combines several innovative features introduced in section 3. The following explanations will be made in reference to the model shown in Figure 4.
The structural model comprises the guide structure in coarse discretization in order to apply proper boundary conditions to the trolley motion, after experimental evidence has shown a vertical displacement of the guide near the impact location, which reduces the effective vertical impact velocity depending on the loading and on time.The guide model has been validated based on comparison with the displacement and velocity time histories of the experiment.An initial velocity vector depending on the specific test case is assigned to the trolley assembly.Trolley and panel are modeled in detail using classic FE shell and beam elements following typical modeling strategies used in crashworthiness analysis of aeronautical structures.Standard elastic-plastic material models with respective material properties for steel and aluminum AL2024-T3/T351 are applied.The connection of the panel to the L-frame is established using point link elements, which represent a purely elastic joint connection.Forces are extracted from the numerical model using beam elements, which are accordingly connected to the trolley to represent the load cells.Evaluation of strain results simply uses the element surface strain components as they are computed at the outer integration point, which is located on the surface location.Local pressure evaluation on the moving structure is accomplished using gauge particles as described in [14].
The fluid domain is modeled using translating periodic boundary conditions (cf.section 3.4), which allow significantly reducing the amount of particles required for the solution.In order to further reduce the computational time, the far field of the fluid domain is modeled using FE solids which are computationally cheaper compared to SPH particles.The two fluid domains interact based on a classic penalty contact algorithm.This coupled modeling approach has been used and validated for fluid-structure interaction problems [2,15].Compatibility with the translating particle domain was achieved by freezing the nodes of the solid elements located outside the translating domain.The minimum size of the water domain was established based on parameter studies varying depth, width and length.It was found that the SPH domain may be rather shallow whereas the rest of the pool depth could be discretized by FE solids without any effect on simulation results.The width of the SPH domain, however, showed a more important effect on simulation results and should be chosen wide enough to allow for an unconstrained free surface motion.
Modeling the fluid domain purely with particles, making use of non-uniform initial particle distributions as described in section 3.3, has shown to allow for large savings in O (10) in the amount of particles compared to a classic uniform particle distribution.Therefore, the computational effort as well as the amount of data (model as well as result files) could be reduced drastically.However, the available nearest neighbor search algorithm in VPS, a cell-linked list algorithm, defines its search cells based on the maximum smoothing length of all particles, which is related to particle volume as V ~ h d where d is the dimension of the problem.Particles that differ significantly in size consequently cause the solution to become less efficient, which adversely affects the efficiency gain due to this modeling approach.Possibilities to overcome this limitation will be investigated in the future.

NUMERICAL RESULTS
This section gives an insight into the numerical results achieved with the developed model presented in the previous section.
In general, qualitative comparison of numerical results with experimental underwater images shows similar behavior in terms of bow wave formation and shape.The fluid motion shows a stronger three-dimensional flow when increasing the pitch angle, because the outflow is less constrained for larger pitch angles.Also the curvature of the panels significantly changes the fluid behavior.Most interestingly, the dynamic structural response of flexible panels causes a time-dependent change of the jet root shape which starts as a nearly straight line and becomes highly three-dimensional when the structural deformation grows.These observations emphasize the necessity for the selected three-dimensional model.Quantitative analysis of force-and strain-time histories revealed very good agreement between simulation and experiment (see Figures 5 and 6 for exemplary results).Force-time curves show a slightly larger initial gradient, but generally match the experimental measurements for both test cases shown.The sharp drop of the vertical force towards the end refers to that time instant when the leading edge of the panel enters the water and the trolley box structure impacts the water.The one large oscillation measured in the experiment at 6 degrees pitch angle at approximately t = 15-25 ms is attributed to oscillations of the guide structure rather than a hydrodynamic effect.Apparently larger oscillations in the numerical result for the 10 degrees test case were found to be related to the particle disorder that is stronger for larger pitch angles.These oscillations can be reduced when applying particle regularization algorithms presented in section 3.2, which help to maintain more regular particle spacing.Another means to reduce oscillations in the numerical model is to refine the particle domain; this however comes at a large increase of computational requirements.
When comparing strain results, maximum strain values are in very good agreement with experimental measurements.However, it was found that there is a time lead for numerical strains.This discrepancy in time is larger for gauges located further ahead.It is currently assumed that the time lead is due to cushioning effects caused by the compressed air layer on top of the water in the experiment, which is generally known to have an important effect when angles between structure and water are small.In the experiment, the air may be compressed and thus push down the free surface, which delays the impact.The numerical model does not include air due to the large additional computational cost associated; therefore it cannot capture this phenomenon.To some extent this time lead may also be caused by the slightly overestimated force gradient at the beginning of the simulations, which results in a faster advance of the pressure peak zone that travels with the jet root along the panel.This argument is further supported by experimental evidence that revealed the simultaneous occurrence of strain and pressure peak values in time and space.
Deformable cases show different phenomena regarding the global hydrodynamic force acting on the structure.Under the presence of deformation, forces generally increase over the entire duration of the impact.Moreover, there is a distinct force peak just prior to the leading edge immersion, which is caused by the water exiting the pocket formed by the deformed structure.This becomes obvious when looking at the contour plot of the simulation in Figure 7, which clearly shows this mechanism.First, the water piles up underneath the structure while deformation increases (large Z-coordinate).Next, the water is pushed out of the pocket (reduced Z-coordinate), which results in larger constraints and thus increases the force.
Comparison of numerical and experimental force-time histories in Figure 8 shows excellent agreement in terms of force magnitudes.Also the simulations capture the aforementioned force peak.However, similar to the findings for the quasi-rigid panels presented above, the timing in the simulation is too early, which is believed to be caused by the missing effect of the air in the model.
It becomes significantly more difficult to capture the strain response of the thin panels.Figure 9 compares strain-time curves of the test case using a panel with 3 mm thickness impacting at 40 m/s horizontal velocity and 6 degrees pitch angle.Despite the time lead of the simulation mentioned, the simulation portrays very good qualitative shapes as well as strain magnitudes for all gauge positions.Strain results along the boundaries of the panel were found to be highly depending on the exact location of the element used to extract the strain as there are large spatial strain gradients in perpendicular direction to the boundaries.This mesh dependency is even stronger for test cases using the 0.8 mm thick panels, where there is a significant amount of bending and plasticity in the vicinity of the boundary, which makes it very challenging to predict the exact structural response.

CONCLUSIONS AND OUTLOOK
In this paper, a coupled SPH-FE approach to model fluid-structure interaction during water impact at large forward velocity has been presented.Several innovative modeling and simulation features were discussed: pressure correction, particle regularization, translating periodic boundary conditions, and non-uniform particle distributions.These features allow for accurate and efficient simulation of the regarded water impact as well as similar applications such as sloshing and slamming of waves on floating structures.It has been demonstrated that the approach is suitable to model fluid-structure interaction as it occurs during water impact at realistic aircraft ditching conditions.Global forces as well as strain responses of structures with different levels of flexibility have been used to validate the adequacy of the model to represent realistic ditching phenomena based on comparison with experimental data.In particular, the structural response of highly deformable structures has been simulated successfully.The presented approach is however reaching its limits when simulating very thin structures, which showed a large mesh dependency.Also the effect of air has been discussed; it is considered responsible for a small discrepancy in timing between experiment and simulation, yet its effects on the structural response within the investigated range of test conditions are minor.
In conclusion, the present work highlighted the importance of structural flexibility on the hydrodynamic loads acting.Consequently, coupled simulations should become the standard for the analysis of structural capacity under ditching loads.With the modeling and simulation techniques available, it becomes possible for aircraft manufacturers to significantly benefit from numerical analysis, which can reduce the substantial investments in time and budget necessary for experimental testing.Future work will be directed at further analyses of the mechanisms involved in this type of FSI, especially for other test conditions available from the unique campaign of guided ditching tests.Further code developments should focus on high performance computing capabilities in order to further reduce the solution times associated with such detailed simulation models.The ultimate goal is numerical simulation of ditching scenarios of full-scale aircraft models, which requires the possibility to employ finite element models of deformable aircraft and to include additional physical phenomena involved such as air-cushioning, suction effects and aerodynamics.The first requirement may be fulfilled by the use of increasingly complex aircraft models from rigid bodies to deformable models of the fuselage and other relevant parts in the structure that may come into contact with the water [21].Inclusion of air-cushioning effects is possible by definition of a region of SPH particles representing the air beneath the approaching structure, as has already been demonstrated for free-fall lifeboats [22].Suction due to the adhesion of water to the structure can be included in a simple manner and has been demonstrated to significantly influence the kinematics of the aircraft [23].Aerodynamic lift and drag effects may effectively be incorporated by phenomenological forces and moments acting on the aircraft [23].Although the computational demands for these complex phenomena remain high, it is expected that accurate simulation of aircraft ditching will soon become feasible.

Figure 1 :
Figure 1: Two-dimensional representation of the region of influence of particle i.The shading illustrates the strength of the smoothing kernel and its range

Figure 2 :
Figure 2: Pressure field at t = 30 ms in 2D flat plate impact test case using no correction (top), Shepard filtering with cycle frequency f = 20 Hz (center), and Rusanov flux with strength 0.5 (bottom)

4 .
GUIDED DITCHING SIMULATION 4.1.DESCRIPTION OF THE EXPERIMENTS The development of novel modeling and simulation features requires verification and validation.In order to provide experimental data for validation, an extensive test campaign of guided ditching tests has been conducted by CNR-INSEAN as part of the EU-funded FP7research project SMAES (SMart Aircraft in Emergency Situations, 2011-2014).

Figure 3 :
Figure 3: Exemplary impact simulation demonstrating the usage of translating periodic boundaries and the damping zone.Color represents vertical displacement to emphasize the wake behind the structure, which does not reenter at the front face

Figure 4 :Figure 5 :
Figure 4: Overview of numerical model of guided ditching test with zoom on trolley structural model

Figure 6 :
Figure 6: Comparison of strain-time histories in x-and y-direction at rear centerline position for quasi-rigid test case at 40 m/s horizontal impact velocity and pitch angle of 6 degrees (identifier 1122)

Figure 7 :
Figure 7: Volume cut of three contour plot instances of test cases at 6 degrees pitch angle, 40 m/s horizontal impact velocity and panel thickness of 0.8 mm (identifier 3122).Note that parts of the model are hidden for clarity of the illustration

Figure 8 :
Figure 8: Comparison of force-time histories for test cases at 6 degrees pitch angle, 40 m/s horizontal impact velocity and panel thickness of 3 mm (left, identifier 2122) and 0.8 mm (right, identifier 3122)

Figure 9 :
Figure 9: Comparison of strain-time histories in x-and y-direction along centerline for 3 mm case at 40 m/s impact velocity and pitch angle of 6 degrees (identifier 2122)