Challenges of Fully-Coupled High-Fidelity Ditching Simulations

: An important element of the process of aircraft certiﬁcation is the demonstration of the crashworthiness of the structure in the event of an emergency landing on water, also referred to as ditching. Novel numerical simulation methods, that incorporate the interaction between ﬂuid and structure, open up a promising way to model ditching in full scale. This study focuses on two main issues of high-ﬁdelity ditching simulations: the development of a suitable ﬂuid-structure coupling framework and the generation of the structural model of the aircraft. The ﬁrst issue is addressed by implementing a partitioned coupling approach, which combines a ﬁnite volume hydrodynamic ﬂuid solver as well as a ﬁnite element structural solver. The developed framework is validated by means of two ditching-like experiments, which consider the drop test of a rigid cylinder and a deformable cylindrical shell. The results of the validation studies indicate that an alternative to the standard Dirichlet-Neumann partitioning approach is needed if a strong added-mass effect is present. For the full-scale simulation of aircraft ditching, structural models become more complex and have to account for damage. Due to its high localization, the damage creates large differences in model scale and usually entails severe non-linearities in the model. To address the issue of increasing computational effort for such models, the process of developing a multi-scale model for the simulation of the failure of fuselage frames is presented.


Introduction
Commercial transport aircraft often fly above open waters and, thus, have to prove the safe landing on water. This process is also referred to as ditching and is part of the certification specifications by the European Aviation Safety Agency (EASA) and the Federal Aviation Administration (FAA) (see [1,2]). Those regulations aim for minimizing the risk of injury for the occupants of the aircraft during the impact phase as well as allowing a safe evacuation in the subsequent floatation phase. This requires an accurate prediction of the structural behaviour of the aircraft during the impact phase as well the damage and its influence on the floatation capabilities.
A common practice to certify aircraft for ditching is to compare the new structural design with existing aircraft configurations [3]. As this approach impedes the development of novel, unconventional aircraft configurations, contemporary aircraft design relies on alternative methods to analyse ditching. Experimental investigations of ditching are limited to low horizontal speed (e.g., helicopters [4]) or are conducted on model scale [5,6] where they suffer from scaling effects, particularly with respect to the accurate representation of hydrodynamic effects [7]. Consequently, there is a motivation to employ numerical tools to simulate ditching. Those tools can be broadly multiscale structural model of the aircraft fuselage. The particular challenge is to accurately model structural failure while limiting the computational effort to a manageable amount. A summary along with conclusions are given in Section 4.

Numerical Simulation Framework
This study employs the coupling environment ifls, which has already been used to model various multi-physics problems (see e.g., [25][26][27]). Ifls follows a modular approach that allows to couple two or more black-box solvers in a partitioned way. The usage of the interpreted, high-level programming language Python facilitates a flexible and easy implementation of the desired coupling scheme. Time-critical operations are performed using compiled code as well as the NumPy library. Figure 1 shows the individual processing units of ifls. Simulation steering and solver communication are performed by the control code and the co-processes, respectively. Each of those units is run in a separate thread. A summary along with conclusions are given in Section 5.

Numerical Simulation Framework
This study employs the coupling environment ifls, which has already been used to model various multi-physics problems (see e.g., [25][26][27]). Ifls follows a modular approach that allows to couple two or more black-box solvers in a partitioned way. The usage of the interpreted, high-level programming language Python facilitates a flexible and easy implementation of the desired coupling scheme. Time-critical operations are performed using compiled code as well as the NumPy library. Figure 1 shows the individual processing units of ifls. Simulation steering and solver communication are performed by the control code and the co-processes, respectively. Each of those units is run in a separate thread.  [25][26][27]).
The control code implements the selected partitioned coupling scheme and sends control commands to the co-processes. Those commands provide direct control over the respective solvers and employ a predefined set of instructions that includes, for example, solving a time step or defining and obtaining boundary conditions. Additional instructions for sending and receiving coupling quantities enable access to the coupling data. NumPy is used to store those data and perform arithmetic operations, thus circumventing the speed limitations of Python. If coupling data is exchanged between different surface meshes, the send and receive instructions implicitly invoke appropriate interpolation routines, which are implemented in converter objects and employ methods from the visualisation toolkit (VTK) (The visualization toolkit is Open-Source software, which can be downloaded from https://www.vtk.org.) As the interpolation is computationally more expensive, all interpolation filters are implemented in C++. Another advantage of the usage of VTK is that it facilitates the easy visualization of the simulation results.
Each co-process implements the instruction set used by the control code where basic features common to all co-processes, for example, inter-thread communication and quantity array management, are separately implemented in a generic parent class. Communication between a coprocess and the corresponding solver can be achieved in different ways. In this case, a socket-based approach is employed, which has superior performance compared to a simple file-based data exchange. Further information on the coupling environment ifls as well as its underlying principles can be found in Reference [28].
For the simulation of the fluid domain, the finite-volume solver interDyMFoam, which is part of the OpenFOAM 5 toolbox, is employed. It solves the incompressible Navier-Stokes equations and is able to model multiphase flow by means of the Volume of Fluid (VOF) method. The fluid equations are solved in the Arbitrary Lagrangian-Eulerian (ALE) form to account for mesh movement in the convective term. The communication between ifls and OpenFOAM is achieved by using a socket- Coupling environment ifls (see [25][26][27]).
The control code implements the selected partitioned coupling scheme and sends control commands to the co-processes. Those commands provide direct control over the respective solvers and employ a predefined set of instructions that includes, for example, solving a time step or defining and obtaining boundary conditions. Additional instructions for sending and receiving coupling quantities enable access to the coupling data. NumPy is used to store those data and perform arithmetic operations, thus circumventing the speed limitations of Python. If coupling data is exchanged between different surface meshes, the send and receive instructions implicitly invoke appropriate interpolation routines, which are implemented in converter objects and employ methods from the visualisation toolkit (VTK) (The visualization toolkit is Open-Source software, which can be downloaded from https://www.vtk.org.) As the interpolation is computationally more expensive, all interpolation filters are implemented in C++. Another advantage of the usage of VTK is that it facilitates the easy visualization of the simulation results.
Each co-process implements the instruction set used by the control code where basic features common to all co-processes, for example, inter-thread communication and quantity array management, are separately implemented in a generic parent class. Communication between a co-process and the corresponding solver can be achieved in different ways. In this case, a socket-based approach is employed, which has superior performance compared to a simple file-based data exchange. Further information on the coupling environment ifls as well as its underlying principles can be found in Reference [28].
For the simulation of the fluid domain, the finite-volume solver interDyMFoam, which is part of the OpenFOAM 5 toolbox, is employed. It solves the incompressible Navier-Stokes equations and is able to model multiphase flow by means of the Volume of Fluid (VOF) method. The fluid equations are solved in the Arbitrary Lagrangian-Eulerian (ALE) form to account for mesh movement in the convective term. The communication between ifls and OpenFOAM is achieved by using a socket-based data exchange as presented in Reference [29]. The structure is simulated with the implicit Abaqus solver assuming linear elastic material behaviour.
The general partitioned FSI problem can be formulated as follows. The computational domain Ω is divided into the structural domain Ω S and the fluid domain Ω F including their respective boundaries ∂Ω S and ∂Ω F . Both domains share the common boundary Γ i def = ∂Ω S ∩ ∂Ω F , referred to as the FSI interface, on which the coupling values are exchanged. On the fluid side of the interface, the Dirichlet boundary condition is applied, where v f and v s represent the fluid and structural velocity vector. Equation (1) corresponds to a no-slip velocity boundary condition in the fluid domain. Conversely, the fluid forces are applied as a Neumann boundary condition to the structural domain: The quantities σ f and σ s represent the fluid and structural stress tensor, respectively, whereas n denotes the unit normal vector.
Ditching exhibits typical characteristics of a strongly coupled FSI problem. For that reason, the implicit fixed-point iteration procedure summarized in Figure 2 is employed. The variables u Γ and v Γ contain the nodal displacements and velocities at the interface Γ i , whereas p Γ represents the nodal forces. The operator S describes the solution in the corresponding computational domain. The mapping of the coupling values at the FSI interface from fluid to structure and vice versa is expressed by T Γ,fs and T Γ,sf . The subscripts s and f denote the structural and fluid domain. The index n represents the current time increment, whereas i stands for the current iteration.
is applied, where and represent the fluid and structural velocity vector. Equation (1) corresponds to a no-slip velocity boundary condition in the fluid domain.
Conversely, the fluid forces are applied as a Neumann boundary condition to the structural domain: , on Γ . ( The quantities and represent the fluid and structural stress tensor, respectively, whereas denotes the unit normal vector.
Ditching exhibits typical characteristics of a strongly coupled FSI problem. For that reason, the implicit fixed-point iteration procedure summarized in Figure 2   The fixed-point scheme is usually suitable for strongly coupled problems, although it requires several equilibrium iterations per time step. To accelerate the convergence, the Aitken relaxation method is used to calculate the relaxation factor [30]. Convergence is checked at the end of each iteration by the following displacement-based criterion: where the constant value ensures that the denominator does not become zero. The fixed-point scheme is usually suitable for strongly coupled problems, although it requires several equilibrium iterations per time step. To accelerate the convergence, the Aitken relaxation method is used to calculate the relaxation factor ω [30]. Convergence is checked at the end of each iteration by the following displacement-based criterion: where the constant value C ensures that the denominator does not become zero.

Drop Test of Rigid Cylinder
The model problem of the drop test of a rigid cylinder is based on an experimental study by Greenhow and Lin [31]. The problem set-up as well as the boundaries of the computational domains are shown in Figure 3a. Boundary Γ f ,I represents the free atmosphere, whereas the boundaries Γ f ,II-IV are fixed walls. Table 1 summarizes the velocity and pressure boundary conditions which are valid in the fluid domain. Moreover, the interface conditions described in Section 2 are applied on the cylinder surface. The diameter of the cylinder is set to d = 0.11 m. The mass is chosen in such a way that the cylinder is neutrally-buoyant. The impact velocity amounts to v entry = 2.955 m/s.

Drop Test of Rigid Cylinder
The model problem of the drop test of a rigid cylinder is based on an experimental study by Greenhow and Lin [31]. The problem set-up as well as the boundaries of the computational domains are shown in Figure 3a. Boundary , represents the free atmosphere, whereas the boundaries , are fixed walls. Table 1 summarizes the velocity and pressure boundary conditions which are valid in the fluid domain. Moreover, the interface conditions described in Section 2 are applied on the cylinder surface. The diameter of the cylinder is set to 0.11 m. The mass is chosen in such a way that the cylinder is neutrally-buoyant. The impact velocity amounts to | | 2.955 m/s.  The solid is meshed with rigid surface elements where point masses are used to model the mass of the cylinder. The fluid domain is spatially discretized by a two-dimensional block-structured mesh with an O-Grid around the cylinder. A section of the initial fluid mesh mesh02 is shown in Figure 3b. This mesh consists of 22k cells and is constructed from 90 elements along the outer edges , and , to 120 elements along , and , . The mesh density increases in closer proximity to the cylinder. Additionally, two h-refinement steps are performed. The resulting meshes are labelled mesh03 and mesh04 and contain 89k and 355k cells, respectively. For the simulation of the rigid validation case, the medium fine mesh03 is used.
Fluid and structural simulation are coupled employing the fixed-point iteration scheme  The solid is meshed with rigid surface elements where point masses are used to model the mass of the cylinder. The fluid domain is spatially discretized by a two-dimensional block-structured mesh with an O-Grid around the cylinder. A section of the initial fluid mesh mesh02 is shown in Figure 3b. This mesh consists of 22k cells and is constructed from 90 elements along the outer edges Γ f ,I and Γ f ,III to 120 elements along Γ f ,II and Γ f ,IV . The mesh density increases in closer proximity to the cylinder. Additionally, two h-refinement steps are performed. The resulting meshes are labelled mesh03 and mesh04 and contain 89k and 355k cells, respectively. For the simulation of the rigid validation case, the medium fine mesh03 is used.
Fluid and structural simulation are coupled employing the fixed-point iteration scheme described in Section 2. The displacement residual is set to r Γ,s,max = 10 −7 m, while a constant time step size of 2.5 × 10 −5 s is chosen. In the structural simulation, the time integration is done by an implicit Newmark-beta method, whereas the fluid solver uses the implicit Euler method. Furthermore, the volume of fluid method in combination with the so-called Multi-Dimensional Limiter for Explicit Solution (MULES) interface compression scheme is employed to model the free water surface. For a detailed description of that scheme, see for example [32]. Figure 4b depicts the depth of penetration over the time and shows good agreement between simulation and experiment. However, the illustration of the different stages of water impact in Figure 5 reveals that the free water surface is not perfectly sharp at the beginning of the simulation. This indicates a strong dependence of the MULES method on the mesh geometry, which was for example also found in Reference [33]. As the simulation results match the experimental results well, it can be assumed that this effect does not have a significant influence on the overall hydrodynamic loads. Limiter for Explicit Solution (MULES) interface compression scheme is employed to model the free water surface. For a detailed description of that scheme, see for example [32]. Figure 4b depicts the depth of penetration over the time and shows good agreement between simulation and experiment. However, the illustration of the different stages of water impact in Figure 5 reveals that the free water surface is not perfectly sharp at the beginning of the simulation. This indicates a strong dependence of the MULES method on the mesh geometry, which was for example also found in Reference [33]. As the simulation results match the experimental results well, it can be assumed that this effect does not have a significant influence on the overall hydrodynamic loads.

Drop Test of Deformable Cylinder
The elastic validation case is based on [34] by Arai and Miyauchi. In their study, the impact of cylindrical shells on water is both experimentally and numerically investigated. The parameters of the second validation case are summarized in Table 2. Limiter for Explicit Solution (MULES) interface compression scheme is employed to model the free water surface. For a detailed description of that scheme, see for example [32]. Figure 4b depicts the depth of penetration over the time and shows good agreement between simulation and experiment. However, the illustration of the different stages of water impact in Figure 5 reveals that the free water surface is not perfectly sharp at the beginning of the simulation. This indicates a strong dependence of the MULES method on the mesh geometry, which was for example also found in Reference [33]. As the simulation results match the experimental results well, it can be assumed that this effect does not have a significant influence on the overall hydrodynamic loads.

Drop Test of Deformable Cylinder
The elastic validation case is based on [34] by Arai and Miyauchi. In their study, the impact of cylindrical shells on water is both experimentally and numerically investigated. The parameters of the second validation case are summarized in Table 2.

Drop Test of Deformable Cylinder
The elastic validation case is based on [34] by Arai and Miyauchi. In their study, the impact of cylindrical shells on water is both experimentally and numerically investigated. The parameters of the second validation case are summarized in Table 2. As proposed in Reference [35], the friction of the sliding system in the experiment taken into account by assuming that the water entry velocity is reduced by 10% to the value shown in Table 2. Furthermore, following a simplified approach proposed by Ionina and Korobkin in Reference [36], additional masses due to cables and connectors are considered through an increase of the shell density by 11%, that is assuming those masses are uniformly distributed. In this context, ρ s,1 represents the density without additional mass, whereas ρ s,2 stands for the cylindrical shell with additional mass.
Despite different dimensions, the discretization of the fluid model is equivalent to the one described in Section 3.1.1, except that all three fluid mesh configurations are investigated in the deformable validation case. The time step is set to a constant value of 1 × 10 −5 s. The time integrators used in the fluid and structural subsystem are the same as for the rigid validation case. The structural mesh shown in profile in Figure 6 consists of elastic shell elements with linear form functions.  Poisson's ratio 0.34 As proposed in Reference [35], the friction of the sliding system in the experiment taken into account by assuming that the water entry velocity is reduced by 10% to the value shown in Table 2. Furthermore, following a simplified approach proposed by Ionina and Korobkin in Reference [36], additional masses due to cables and connectors are considered through an increase of the shell density by 11%, that is assuming those masses are uniformly distributed. In this context, , represents the density without additional mass, whereas , stands for the cylindrical shell with additional mass.
Despite different dimensions, the discretization of the fluid model is equivalent to the one described in Section 3.1.1, except that all three fluid mesh configurations are investigated in the deformable validation case. The time step is set to a constant value of 1 10 s. The time integrators used in the fluid and structural subsystem are the same as for the rigid validation case. The structural mesh shown in profile in Figure 6 consists of elastic shell elements with linear form functions. It represents the coarsest resolution and is used in combination with mesh02. For the finer configurations mesh03 and mesh04, h-refinement steps are similarly applied to the structural mesh. Again, the implicit coupling scheme described in Section 2 is used. In contrast to the rigid case, the displacement residual, which is set to , , 10 m, is complemented by a maximum number of 30 iterations per time increment.  It represents the coarsest resolution and is used in combination with mesh02. For the finer configurations mesh03 and mesh04, h-refinement steps are similarly applied to the structural mesh. Again, the implicit coupling scheme described in Section 2 is used. In contrast to the rigid case, the Aerospace 2019, 6, 10 8 of 15 displacement residual, which is set to r Γ,s,max = 10 −10 m, is complemented by a maximum number of 30 iterations per time increment. Figure 7a compares the experimental strain response at the bottom of the cylinder with the numerical values of all three mesh configurations for a structural density of ρ s = ρ s,2 . The solution is already sufficiently converged for mesh03. For all mesh configurations, the results reveal an overprediction of the magnitude of extreme values. For the two global extreme values, the difference amounts to 13% on average for mesh02 and 20% for mesh04. One reason for this difference could be the insufficient resolution of the experimental data, which does not reproduce higher frequent oscillations in the strain response and may prevent an accurate representation of extreme values. An additional observation is that the dynamic response of the cylinder is slightly faster in the simulation. oscillations in the strain response and may prevent an accurate representation of extreme values. An additional observation is that the dynamic response of the cylinder is slightly faster in the simulation. Similar deviations were found by Sun and Faltinsen in Reference [35] by their two-dimensional boundary element method (BEM). The corresponding results are shown in Figure 7b. The phase shift is similar to the one found in this study. However, the difference in the magnitude of extreme values is smaller. As explained by Sun and Faltinsen, the extreme values significantly depend on the force calculation at the beginning of the simulation. If those forces are calculated by the more accurate Wagner's method (see [37]), the magnitudes of the extreme values become larger.
In summary, two factors can lead to the differences observed between simulation and experiment. On the one hand, the resolution of the experimental data is relatively low. This could lead to inaccurate extreme values. On the other hand, certain assumptions are not fulfilled in reality, for instance, additional masses are not uniformly distributed within the cylinder. This may alter the dynamic behaviour of the structure. Additionally, Figure 7b compares the strain response for the different structural densities , and , . It becomes apparent that negligence of the additional mass leads to a less accurate approximation of the experimental data.
To assess the computational efficiency of the iterative coupling method, the number of iterations required each time increment during the moment of impact is compared for all mesh configurations in Figure 8. Although some increments reach the maximum number of 30 iterations, the average number of iterations per increment over the whole simulation is lower. It amounts to 24.5 for mesh02, 22.9 for mesh03 and 17.9 for mesh04. Anyhow, this number is still too large to directly employ such a coupling scheme for the high-fidelity simulation of a realistic aircraft structure.  Figure 7b. The phase shift is similar to the one found in this study. However, the difference in the magnitude of extreme values is smaller. As explained by Sun and Faltinsen, the extreme values significantly depend on the force calculation at the beginning of the simulation. If those forces are calculated by the more accurate Wagner's method (see [37]), the magnitudes of the extreme values become larger.
In summary, two factors can lead to the differences observed between simulation and experiment. On the one hand, the resolution of the experimental data is relatively low. This could lead to inaccurate extreme values. On the other hand, certain assumptions are not fulfilled in reality, for instance, additional masses are not uniformly distributed within the cylinder. This may alter the dynamic behaviour of the structure.
Additionally, Figure 7b compares the strain response for the different structural densities ρ s,1 and ρ s,2 . It becomes apparent that negligence of the additional mass leads to a less accurate approximation of the experimental data.
To assess the computational efficiency of the iterative coupling method, the number of iterations required each time increment during the moment of impact is compared for all mesh configurations in Figure 8. Although some increments reach the maximum number of 30 iterations, the average number of iterations per increment over the whole simulation is lower. It amounts to 24.5 for mesh02, 22.9 for mesh03 and 17.9 for mesh04. Anyhow, this number is still too large to directly employ such a coupling scheme for the high-fidelity simulation of a realistic aircraft structure.  For that reason, alternative coupling procedures are needed. Also, as the modelling of structural damage is a major component of high-fidelity ditching simulations, explicit time integrations schemes should be used in the structural simulation. A simple staggered coupling scheme, which is exchanging the coupling data only once per time step, seems to be a viable solution. However, this scheme is not stable if standard Dirichlet-Neumann boundary conditions are employed and, thus, requires alternative partitioning approaches.

Global FE-Model
Preliminary ditching simulations with a rigid aircraft similar to the A350 show that an approximately 10 m long section of the rear fuselage is subjected to water loads during the impact phase of a ditching event. The relevant section is shown in Figure 9. Though only a small part of the fuselage is subjected to the hydrodynamic impact loads, the whole aircraft needs to be structurally modelled to capture the overall dynamics during impact. If structural damage occurs, it is expected to be limited to the green region. To capture the overall aircraft dynamics outside the green region, a global FEM model (GFEM) based on the XRF-1, an Airbus large transport aircraft of the eXternal Research Forum (XRF), is used. Stringers and frames are modelled as rods and beams on a skin modelled as shells, whereas masses are defined as point masses which are assigned to whole fuselage sections. The total number of degrees of freedom in the model amounts to approximately 300k.

Fuselage Crash Simulations
The majority of the available methods to simulate fuselage damage under crash conditions primarily address barrel drop tests on solid ground. Two distinct numerical approaches can be found in the literature.
The first approach is to discretize the airframe as a lumped mass-spring-damper system. This leads to small models with only several hundred degrees of freedom. Examples of such an implementation are DRI-KRASH solver described in Reference [13] or a method by Weiß proposed in For that reason, alternative coupling procedures are needed. Also, as the modelling of structural damage is a major component of high-fidelity ditching simulations, explicit time integrations schemes should be used in the structural simulation. A simple staggered coupling scheme, which is exchanging the coupling data only once per time step, seems to be a viable solution. However, this scheme is not stable if standard Dirichlet-Neumann boundary conditions are employed and, thus, requires alternative partitioning approaches.

Global FE-Model
Preliminary ditching simulations with a rigid aircraft similar to the A350 show that an approximately 10 m long section of the rear fuselage is subjected to water loads during the impact phase of a ditching event. The relevant section is shown in Figure 9. Though only a small part of the fuselage is subjected to the hydrodynamic impact loads, the whole aircraft needs to be structurally modelled to capture the overall dynamics during impact. If structural damage occurs, it is expected to be limited to the green region.  For that reason, alternative coupling procedures are needed. Also, as the modelling of structural damage is a major component of high-fidelity ditching simulations, explicit time integrations schemes should be used in the structural simulation. A simple staggered coupling scheme, which is exchanging the coupling data only once per time step, seems to be a viable solution. However, this scheme is not stable if standard Dirichlet-Neumann boundary conditions are employed and, thus, requires alternative partitioning approaches.

Global FE-Model
Preliminary ditching simulations with a rigid aircraft similar to the A350 show that an approximately 10 m long section of the rear fuselage is subjected to water loads during the impact phase of a ditching event. The relevant section is shown in Figure 9. Though only a small part of the fuselage is subjected to the hydrodynamic impact loads, the whole aircraft needs to be structurally modelled to capture the overall dynamics during impact. If structural damage occurs, it is expected to be limited to the green region. To capture the overall aircraft dynamics outside the green region, a global FEM model (GFEM) based on the XRF-1, an Airbus large transport aircraft of the eXternal Research Forum (XRF), is used. Stringers and frames are modelled as rods and beams on a skin modelled as shells, whereas masses are defined as point masses which are assigned to whole fuselage sections. The total number of degrees of freedom in the model amounts to approximately 300k.

Fuselage Crash Simulations
The majority of the available methods to simulate fuselage damage under crash conditions primarily address barrel drop tests on solid ground. Two distinct numerical approaches can be found in the literature.
The first approach is to discretize the airframe as a lumped mass-spring-damper system. This leads to small models with only several hundred degrees of freedom. Examples of such an implementation are DRI-KRASH solver described in Reference [13] or a method by Weiß proposed in To capture the overall aircraft dynamics outside the green region, a global FEM model (GFEM) based on the XRF-1, an Airbus large transport aircraft of the eXternal Research Forum (XRF), is used. Stringers and frames are modelled as rods and beams on a skin modelled as shells, whereas masses are defined as point masses which are assigned to whole fuselage sections. The total number of degrees of freedom in the model amounts to approximately 300k.

Fuselage Crash Simulations
The majority of the available methods to simulate fuselage damage under crash conditions primarily address barrel drop tests on solid ground. Two distinct numerical approaches can be found in the literature.
The first approach is to discretize the airframe as a lumped mass-spring-damper system. This leads to small models with only several hundred degrees of freedom. Examples of such an implementation are DRI-KRASH solver described in Reference [13] or a method by Weiß proposed in Reference [38]. Obvious advantages of those approaches are their simplicity of application and reduced computational effort. This allows examining many different configurations and crash scenarios. However, the parameters needed for those models have to be obtained experimentally or from experience. Furthermore, the results need to be carefully reviewed as the range of validity cannot be reliably predicted.
The other approach is the simulation on a high-fidelity level using explicit FE solvers [39]. Although this approach tends to be more accurate, experimental verification still remains an important requirement. Another major disadvantage is the high computational cost. Besides the model size, the computational effort is dependent on the physical time to be simulated as well as the refinement level in the model. Especially the latter can cause problems as the maximum time step size is directly linked to the size of the smallest element in the model. A second disadvantage is that the required degree of granularity in the model needs to be determined. In this context, the term granularity not only refers to the element size but also to the level of structural detail, for instance, if discrete rivets with failure criterions are modelled. For barrel drop tests, suitable degrees of granularity were determined by comparisons between simulations and experiments. However, the question if those results are applicable to the simulation of ditching events still remains unanswered.
If a direct simulation of damage with the explicit FE method is employed, the structural models need to be fine enough to account for relevant damage effects. To illustrate this, the simulation of a breaking frame loaded by a bending moment is shown in Figure 10b. As frames serve as circumferential stiffeners, their failure has the biggest impact on the overall stiffness of the fuselage [38]. The damage is highly localized and the damage region is magnitudes of order smaller compared to the overall aircraft dimensions. Reference [38]. Obvious advantages of those approaches are their simplicity of application and reduced computational effort. This allows examining many different configurations and crash scenarios. However, the parameters needed for those models have to be obtained experimentally or from experience. Furthermore, the results need to be carefully reviewed as the range of validity cannot be reliably predicted. The other approach is the simulation on a high-fidelity level using explicit FE solvers [39]. Although this approach tends to be more accurate, experimental verification still remains an important requirement. Another major disadvantage is the high computational cost. Besides the model size, the computational effort is dependent on the physical time to be simulated as well as the refinement level in the model. Especially the latter can cause problems as the maximum time step size is directly linked to the size of the smallest element in the model. A second disadvantage is that the required degree of granularity in the model needs to be determined. In this context, the term granularity not only refers to the element size but also to the level of structural detail, for instance, if discrete rivets with failure criterions are modelled. For barrel drop tests, suitable degrees of granularity were determined by comparisons between simulations and experiments. However, the question if those results are applicable to the simulation of ditching events still remains unanswered.
If a direct simulation of damage with the explicit FE method is employed, the structural models need to be fine enough to account for relevant damage effects. To illustrate this, the simulation of a breaking frame loaded by a bending moment is shown in Figure 10b. As frames serve as circumferential stiffeners, their failure has the biggest impact on the overall stiffness of the fuselage [38]. The damage is highly localized and the damage region is magnitudes of order smaller compared to the overall aircraft dimensions. Due to complex damage mechanisms, that are triggered by stability buckling in the inner flange or by plasticity effects, respectively, a very fine mesh is needed to directly simulate the damage. Consequently, a mesh size study is conducted for the case depicted in Figure 10b. The results are summarized in Figure 11, showing the bending angle over the applied momentum for different element edge lengths. The results indicate that very fine meshes are needed to get a reliable ultimate momentum. In addition, the post-failure stiffness is very dependent on the mesh size, as many damage effects are very small in scale. Due to complex damage mechanisms, that are triggered by stability buckling in the inner flange or by plasticity effects, respectively, a very fine mesh is needed to directly simulate the damage. Consequently, a mesh size study is conducted for the case depicted in Figure 10b. The results are summarized in Figure 11, showing the bending angle over the applied momentum for different element edge lengths. The results indicate that very fine meshes are needed to get a reliable ultimate momentum. In addition, the post-failure stiffness is very dependent on the mesh size, as many damage effects are very small in scale. Aerospace 2018, 5, x FOR PEER REVIEW 11 of 15 Aerospace 2018, 5, x; doi: FOR PEER REVIEW www.mdpi.com/journal/aerospace Figure 11. Results of the model in Figure 10b with different element side lengths.
The complete ditching event is expected to have a duration of several seconds and, as the computation effort grows with the simulated time, an FE model that accounts for structural damage within the region relevant for ditching by using a sufficiently fine mesh while simultaneously capturing the whole aircraft is not feasible. Therefore, the authors suggest using a multiscale model for the structural simulation.

Multiscale Model for Ditching Simulation
To model the section of the aircraft relevant for ditching, a combination of the two methods outlined in Section 3.2.2 is used. The green section shown in Figure 9 is modelled using an FE discretization with an element size between 15 and 20 mm. The water loads are applied to the fuselage skin. To allow the correct simulation of the load paths within the model, stringers, clips and frames are modelled as shells. The resulting model has approximately 6M degrees of freedom and, therefore, can be handled by commercial explicit FE solvers. The blue region in Figure 9 resembles the GFEMmodel. Both parts of the model are connected by multi-point-constraints.
In addition to the refined mesh, a better approximation of the mass distribution is needed as inertia forces are important for the loading state during highly dynamic events. The structural masses of the refined region are directly incorporated by material densities while passenger and cargo masses are separately modelled. Afterwards, the difference in mass to the original GFEM section is evenly smeared over all elements of the model. In this way, the original mass as well as the overall mass properties of the aircraft, which are essential for the computation of the correct trajectory during ditching, are preserved. Although the moments of inertia in the green section are not preserved exactly, their influence on the overall moment of inertia is relatively small.

Modelling Damage
At this stage, the structural model is still too coarse to accurately account for damage effects as discussed in Section 3.2.2. Most parts of the structure are subjected to characteristic loading states. The stringers in combination with the skin are for example primarily loaded by longitudinal stresses introduced by bending of the fuselage, whereas the frames keep the fuselage cross sections in shape and are mainly loaded by bending moments and hoop stresses. The clips transfer the loads from the skin to the frame and are therefore subjected to shear stresses. These observations can be used to introduce reduced order models that can capture damage effects using a relatively coarse FE discretization. This approach is called Kinematics Model and was described by Waimer et al. in Reference [13]. For instance, the damage of a frame can be modelled as a non-linear hinge element with the centre of rotation located at the skin. The hinge is closed for uncritical loads and starts to open as soon as the frame starts to break. The moment to opening angle curve is obtained either from experimental results or, as done in his case, by using the results of a detailed FE simulation as shown for example in Figure 11. Due to the highly localized damage, its overall effect on the structural behaviour can be captured by such approaches with sufficient accuracy.

Moment/kNm
Angle/deg Figure 11. Results of the model in Figure 10b with different element side lengths.
The complete ditching event is expected to have a duration of several seconds and, as the computation effort grows with the simulated time, an FE model that accounts for structural damage within the region relevant for ditching by using a sufficiently fine mesh while simultaneously capturing the whole aircraft is not feasible. Therefore, the authors suggest using a multiscale model for the structural simulation.

Multiscale Model for Ditching Simulation
To model the section of the aircraft relevant for ditching, a combination of the two methods outlined in Section 3.2.2 is used. The green section shown in Figure 9 is modelled using an FE discretization with an element size between 15 and 20 mm. The water loads are applied to the fuselage skin. To allow the correct simulation of the load paths within the model, stringers, clips and frames are modelled as shells. The resulting model has approximately 6M degrees of freedom and, therefore, can be handled by commercial explicit FE solvers. The blue region in Figure 9 resembles the GFEM-model. Both parts of the model are connected by multi-point-constraints.
In addition to the refined mesh, a better approximation of the mass distribution is needed as inertia forces are important for the loading state during highly dynamic events. The structural masses of the refined region are directly incorporated by material densities while passenger and cargo masses are separately modelled. Afterwards, the difference in mass to the original GFEM section is evenly smeared over all elements of the model. In this way, the original mass as well as the overall mass properties of the aircraft, which are essential for the computation of the correct trajectory during ditching, are preserved. Although the moments of inertia in the green section are not preserved exactly, their influence on the overall moment of inertia is relatively small.

Modelling Damage
At this stage, the structural model is still too coarse to accurately account for damage effects as discussed in Section 3.2.2. Most parts of the structure are subjected to characteristic loading states. The stringers in combination with the skin are for example primarily loaded by longitudinal stresses introduced by bending of the fuselage, whereas the frames keep the fuselage cross sections in shape and are mainly loaded by bending moments and hoop stresses. The clips transfer the loads from the skin to the frame and are therefore subjected to shear stresses. These observations can be used to introduce reduced order models that can capture damage effects using a relatively coarse FE discretization. This approach is called Kinematics Model and was described by Waimer et al. in Reference [13]. For instance, the damage of a frame can be modelled as a non-linear hinge element with the centre of rotation located at the skin. The hinge is closed for uncritical loads and starts to open as soon as the frame starts to break. The moment to opening angle curve is obtained either from experimental results or, as done in his case, by using the results of a detailed FE simulation as shown for example in Figure 11. Due to the highly localized damage, its overall effect on the structural behaviour can be captured by such approaches with sufficient accuracy.

Creating the Multiscale Model
To create a reduced order damage model without having access to full scale ditching experiments, it is necessary to identify locations of the structure where damage might occur as well as the characteristic loading states of the corresponding parts. In addition to the frames, the vertical struts are considered as well as they have a significant influence on the damage characteristics [38]. For the creation of the multiscale model, the iterative process shown in Figure 12 is proposed. To create a reduced order damage model without having access to full scale ditching experiments, it is necessary to identify locations of the structure where damage might occur as well as the characteristic loading states of the corresponding parts. In addition to the frames, the vertical struts are considered as well as they have a significant influence on the damage characteristics [38]. For the creation of the multiscale model, the iterative process shown in Figure 12 is proposed. The initial FE model is loaded by water loads obtained from a ditching simulation of a rigid structure with similar outer geometry. Highly stressed regions are identified and the characteristic loading states are extracted. For the identified regions, FE damage models with a high level of granularity are created and loaded with their respective characteristic loading states. Based on the results of the detailed simulations, appropriate surrogate models, that are able to reproduce the kinematics of the damage, are chosen and provided with the calculated load to deformation curve. The surrogate models are then introduced into the fuselage FE model. The procedure is repeated as the occurrence of damage may change the load path of the structure. If no further critical locations are found the iteration process is stopped.

Conclusions and Outlook
In this study, a fluid-structure coupling framework suitable for ditching simulations is presented and validated. The partitioning of the computational domain is done by the standard Dirichlet-Neumann approach. It is shown that this framework works well for two-dimensional cases with rigid and deformable structures if an implicit fixed-point iteration scheme is used. However, in the deformable case, this scheme proves to be inefficient for the iterative coupling scheme and even becomes unstable if an explicit structural solver is used together with a simple staggered coupling approach. Those results indicate that carefully designed coupling procedures need to be developed to increase computational efficiency and facilitate the usage of explicit structural solvers in the simulation of ditching problems. In this context, a simple staggered coupling scheme would be favourable but certainly requires partitioning approaches other than standard Dirichlet-Neumann boundary conditions. Possible options include operator-splitting techniques (e.g., [40]), mixed boundary conditions (e.g., [41,42]) or the weak introduction of coupling conditions (e.g., [43]). For full-scale ditching simulations, which include the failure of the structure, the time scales of fluid and structural simulation are expected to differ significantly. Hence, sub-cycling of the structural solver will likely be necessary to keep the computational effort at a manageable amount.
Furthermore, as soon as horizontal speed is introduced and mesh deformation becomes an issue, the overset grid technique is a viable option to allow large movements of the fluid-structure boundary while at the same time preserving the quality of the fluid mesh [44].
The second key point of this study pertains the development of a multiscale structural model of the aircraft fuselage. It is shown that a structural model that on the one hand can capture the overall dynamics of the aircraft and, on the other hand, directly simulate structural damage is not feasible. Due to severe nonlinearities in the structure, explicit time integration is needed. However, simply refining damage relevant regions would lead to impracticable short time steps and very poor overall computational efficiency. A method is presented how the highly localized damage behaviour of  The initial FE model is loaded by water loads obtained from a ditching simulation of a rigid structure with similar outer geometry. Highly stressed regions are identified and the characteristic loading states are extracted. For the identified regions, FE damage models with a high level of granularity are created and loaded with their respective characteristic loading states. Based on the results of the detailed simulations, appropriate surrogate models, that are able to reproduce the kinematics of the damage, are chosen and provided with the calculated load to deformation curve. The surrogate models are then introduced into the fuselage FE model. The procedure is repeated as the occurrence of damage may change the load path of the structure. If no further critical locations are found the iteration process is stopped.

Conclusions and Outlook
In this study, a fluid-structure coupling framework suitable for ditching simulations is presented and validated. The partitioning of the computational domain is done by the standard Dirichlet-Neumann approach. It is shown that this framework works well for two-dimensional cases with rigid and deformable structures if an implicit fixed-point iteration scheme is used. However, in the deformable case, this scheme proves to be inefficient for the iterative coupling scheme and even becomes unstable if an explicit structural solver is used together with a simple staggered coupling approach. Those results indicate that carefully designed coupling procedures need to be developed to increase computational efficiency and facilitate the usage of explicit structural solvers in the simulation of ditching problems. In this context, a simple staggered coupling scheme would be favourable but certainly requires partitioning approaches other than standard Dirichlet-Neumann boundary conditions. Possible options include operator-splitting techniques (e.g., [40]), mixed boundary conditions (e.g., [41,42]) or the weak introduction of coupling conditions (e.g., [43]). For full-scale ditching simulations, which include the failure of the structure, the time scales of fluid and structural simulation are expected to differ significantly. Hence, sub-cycling of the structural solver will likely be necessary to keep the computational effort at a manageable amount.
Furthermore, as soon as horizontal speed is introduced and mesh deformation becomes an issue, the overset grid technique is a viable option to allow large movements of the fluid-structure boundary while at the same time preserving the quality of the fluid mesh [44].
The second key point of this study pertains the development of a multiscale structural model of the aircraft fuselage. It is shown that a structural model that on the one hand can capture the overall dynamics of the aircraft and, on the other hand, directly simulate structural damage is not feasible. Due to severe nonlinearities in the structure, explicit time integration is needed. However, simply refining damage relevant regions would lead to impracticable short time steps and very poor overall computational efficiency. A method is presented how the highly localized damage behaviour of aeronautical structures can be incorporated in full-size aircraft models while retaining a model size that can be handled by today's computers. This is done by reduced order models that account for the correct kinematics of structural damage and the combination of a GFEM model with a locally refined section of interest. An iterative process is proposed how relevant regions, where reduced order damage models should be placed, can be identified.