Numerical Seakeeping Analysis for a Floating Helicopter After Ditching in Water

A recent report of global airline fatal accidents elaborated in the Annual Safety Review 2021 [1] shows a decreasing tendency of fatal accidents involving large airplanes for passenger and cargo operations over the last 10 years. This decreasing tendency is also observed for fatal and non-fatal accidents in all helicopter operations. Currently, 75% of accidents and incidents are related to human factors [2]. A reason for such a decrease in fatalities is the requirement of an investigation of accidents and some incidents, which helps to understand the causes and prevent similar incidents in the future.

A type of accident involving aircraft (airplanes and helicopters) is the so-called ditching: the process of a planned and controlled emergency landing of an aircraft in water. National Transportation Safety Board [3] classifies this as “a planned event in which a flight crew knowingly makes a controlled emergency landing in water (excludes floatplane landings in normal water landing areas).” In other words, an aircraft designed to land on water, such as a floatplane or flying boat, is not defined as a ditching accident. As ditching is a process conducted within an emergency situation, it is not unusual to yield injuries and fatalities. A recent example of an accident with fatalities involving helicopter ditching occurred in 2018, while a helicopter tried to land on the East River [4]. The pilot was able to release his restraint underwater and successfully egress from the helicopter, but unfortunately, none of the passengers were able to egress, and they all drowned.

Generally, all helicopters must be designed to be compliant with the Certification Specification of different authorities (e.g., European Union Aviation Safety Agency (EASA), Federal Aviation Administration (FAA), Civil Aviation Authority (CAA)). As a guide for showing compliance, EASA has published the means of compliance (AMC) for Small Rotorcraft [5] for helicopters in offshore operations. These documents define requirements for helicopters to withstand impact loads during ditching and to remain dynamically stable (flotation stability) in irregular seas. European Union Aviation Safety Agency [6] sets specifications for ditching (as an impact or immersion process) and the flotation stability of helicopters, requiring a realistic scenario for experimental or numerical tests. Numerical investigations of flotation stability are recommended to be performed for irregular sea conditions.

So far, semi-empirical methods have been used to predict the ditching of aircraft. Increasing efforts have been made recently to investigate ditching with experimental tests and more complex numerical simulation methods. The project “Smart Aircraft in Emergency Situations” [7] dealt with the development of methods to obtain load distributions that occur when aircraft ditches. In the EU-funded project “Increased Safety & Robust Certification for ditching of Aircrafts & Helicopter” (SARAH) [8], numerical and experimental methods were used to investigate the ditching of airplanes and helicopters, but the flotation stability was not addressed.

After landing on water, it is important to ensure the stability of the helicopter in irregular waves to avoid capsizing and ensure the safety of the crew passengers to get out of the helicopter and wait for rescue from the emergency unit. As helicopter landings in water are rare, it makes the event study of capsizing in helicopters more scarce. Thus, the field of dynamic stability in maritime engineering is an option that can be further studied to be adapted and applied to the current problem. Dynamic stability involves an assessment of an object’s capacity to return to a stable position after some disturbance. For maritime engineering, stability is widely analyzed for vessels in rough seas. Manderbacka et al. [9] presented some insights into contemporary research on ship stability and identified the possible directions for future research. As it is pointed out, the computational fluid dynamics (CFD) is increasingly becoming a practical tool to assess a ship’s dynamics and stability. Neves [10] focused on the dynamic stability of ships, presenting some of the main developments in these fields over the years, with special emphasis on covering the 12th International Conference on the Stability of Ships and Ocean Vehicles (STAB 2015). It discussed that the main objective of the study of ship dynamics is to ensure safety against capsizing, but it is not a simple task mainly for two reasons: (1) for irregular random seaway, capsizing is a stochastic phenomenon and (2) capsizing in ships are rare events, so it makes dependent of several ship parameters, as a consequence of the very low probability levels involved.

Due to the engine that is typically positioned above the cabin, the center of gravity (CG) of helicopters is in a higher position. Additionally, compared to airplanes, which have a large wetted area due to the wings when landing in water, the helicopter tends to have a higher probability of capsizing than airplanes. For this reason, the emergency flotation system (EFS) for helicopters was developed. Reilly [11] developed an EFS for the Boeing CH-46 helicopter, which led to increased stability for negligible weight. The system was first successfully tested in model tests and then in the full-scale. Wilson and Tucker [12] developed an automated EFS in which the buoyancy devices are filled with compressed helium gas in an emergency. An overview of helicopter ditching and flotation systems can be found in Ref. [13]. Flotation systems were also tested concerning their effectiveness in various ditching scenarios [14].

Some studies of flotation stabilities in airplanes can be observed in Ref. [15], which carried out the first extensive investigations on this topic. Ahlvin and Brown [16] described the first developments of criteria and minimum requirements for aircraft flotation stability. However, the topic of flotation stability for helicopter ditching is more scarce. Studies on this topic for maritime applications can be applied to the current problem, but carefully, especially because helicopters have different geometry than ships.

Flotation stability plays a crucial role in the safety of helicopters in water, preventing capsizing and ensuring the well-being of the crew during emergency landings or ditching maneuvers at sea. This study aims to investigate the seakeeping behavior of helicopters under realistic scenarios, considering a realistic helicopter geometry. CFD techniques, specifically the cell-centered finite volume method and boundary element methods (BEM), are employed to assess motion responses and load distribution. While an initial study was presented in Ref. [17], it lacked certain numerical recommended practices and a comprehensive analysis of irregular waves. This research addresses these gaps, ensuring reliable numerical results. Natural periods are estimated, and motion and force analysis in regular waves are conducted across different amplitudes and directions. Additionally, the helicopter’s behavior is simulated in irregular waves based on recommended guidelines, with a focus on peak movements such as roll motion and amplitude variations.

The studied geometry is a realistic representation of a helicopter, with realistic dimensions and mass distribution. This geometry is composed of the helicopter fuselage, skid, and the EFS, composed of four floaters already activated and the gas cylinder, shown in Fig. 1, along with the adopted notation and symbols.

Figure 2 shows the main dimensions of the helicopter. Note that the CG is not aligned with the symmetry plane of the fuselage, being slightly to the right, because of the weight of the cylinder used to fill the floaters.

The mass configuration is such that it represents one of the heaviest conditions on the envelope of flight. Table 1 presents some of the main characteristics of the helicopter, including some hydrostatic main particulars. The characteristic length $L$ is defined as the length of the submerged part of the helicopter in the hydrostatic condition.

The Regulations AMC to CS 27.801(e) and 27.802(c) [5] recommend testing the helicopter with irregular waves, representing a Northern North Sea wave climate, given by a JONSWAP wave spectrum with peak enhancement factor of $γ=3.3$⁠. Table 2 shows the significant wave heights $Hs$⁠, the mean wave periods $Tz$⁠, the peak wave periods $Tp$ (defined as $Tp=1.29Tz$⁠), and the significant steepness (defined as $Ss=2πHs/gTz2$⁠). All cases represent waves in sea state 6.

Two methods are used: CFD (by using the finite volume method) and the BEM. Panel methods are a good approach to obtaining the response amplitude operator (RAO) with low computational time, which makes it attractive. Combined with the statistics of the sea behavior, it allows to assess the short-term statistics of motions and loads. However, the helicopter geometry may not be suitable to be used with the traditional potential method approach. To investigate if the method could be feasibly adopted, results are compared with CFD coupled with the equation of motions in an implicit way, to understand the limitations of the potential method for this specific floating body. A similar comparison was also done in Ref. [18], but focused in Wind Energy Offshore Platform: it investigates the motion predictions of floating bodies in extreme waves, focusing on an offshore platform in a wind farm. It critically compares the accuracy, applicability, and computational costs of three different seakeeping models: a frequency-domain boundary element method, a partly nonlinear time-domain method, and a nonlinear viscous model based on the URANS equations.

CFD simulations are conducted based on Navier–Stokes equations, using the cell-centered finite volume method. From now on, the expression CFD is referred to this method in specific.

Another approach to simulate the behavior of a floating offshore body is to solve the potential flow problem using the radiation and diffraction theory. To solve it, the BEM is adopted. The fluid flow is assumed to have some simplifications to be modeled as linear potential flow, such as:

• The fluid is inviscid and incompressible, and the fluid flow is irrotational;

• The body has zero or very small forward speed;

• The incident regular wave is of small amplitude compared to its length (small slope); and

• The motions are to the first-order and hence must be of small amplitude.

Some conditions need to be satisfied to solve the boundary value problems, as described in the following:

Second-order hydrodynamic loads are described by Refs. [19,20]. These loads are proportional to the square of the wave amplitude and have frequencies equal to the sum and difference of pairs of incident wave frequencies. They can be expressed in three terms: mean drift forces, difference-frequency forces, and sum-frequency forces.

The second-order forces are represented using the quadratic transfer function (QTF) matrix. The QTF matrix quantifies the nonlinear interactions between different wave components. The QTF matrix has two main formulations: the Newman approximation and the full QTF matrix. The Newman approximation method requires far fewer data, making it less computationally intensive, as it only requires the diagonal QTF data defining the mean wave drift load. The Newman approximation then extrapolates the mean wave drift QTFs to approximate the off-diagonal QTF values, which allows the calculation of the slowly-varying part of the wave drift load. However, it is not appropriate for the sum frequency load. The full QTF matrix is generally preferable as it provides a more accurate representation of the wave loads: it uses the full QTF matrix, which includes both the diagonal and off-diagonal elements. The complete implementation of these formulations, including the computation of the QTF matrix, can be found in Refs. [21,22].

Initially, simulations are conducted in the frequency domain for various frequencies and wave encounter angles. The results of these simulations are then used to perform stability analysis of the hydrodynamic response. This analysis determines the initial position of the system for the subsequent time-domain analysis. Finally, it is possible to perform the time-domain analysis in irregular waves, which comprise several frequencies.

Simulating the hydrodynamics of offshore structures in waves usually considers the scenario of the semi-infinite domain. While the CFD needs to discretize the whole finite domain, BEM needs discretization mainly on the boundaries. This simplicity makes BEM promising for this type of hydrodynamic analysis, having reduced computational costs when compared to CFD, especially when assessing several sea state scenarios, or simulating for a considerable longer time.

One limitation of this model is that it does not consider the viscosity of the fluid, which may be important as the capsizing is associated with the rolling motion. Also, during the seakeeping calculation using BEM, it assumes large bodies, with a small KC number (Keulegan–Carpenter number, KC, indicates the relative importance of the drag forces over inertia forces for bluff objects in an oscillatory fluid flow. For small KC number, inertia is dominant).

This section outlines the numerical setup for each method, CFD and BEM, incorporating recommended practices derived from previous results in the maritime field [23–26].

To simulate the six degrees-of-freedom of the helicopter, two computational domains (regions) are proposed: one represents the ocean, and the other contains the helicopter. The dimensions of each region and the boundary conditions are shown in Fig. 3. Note that the helicopter region is covered by the overset boundary, which interpolates with the main domain mesh.

The overset technique is employed to simulate the motion of the helicopter. In this approach, the region containing the helicopter moves, while the region representing the ocean remains static. A key feature of the overset technique is the interpolation between these two overlapping meshes at each time-step. Figure 4 illustrates the overset technique in use: the top part shows the overlapping meshes after the helicopter mesh is moved based on its equation of motion. The bottom part depicts the mesh after interpolation in the overlap region: in regions away from the helicopter, the governing equations are solved in control volumes of the domain mesh; overlap regions use the helicopter mesh; and at the edge, an intermediate cell zone with interpolation between the two meshes is performed in both directions. An essential advantage of this method is that it captures the dynamics of the helicopter’s motion within the ocean environment, eliminating the need for constant remeshing.

The numerical treatment for velocity inlet boundary conditions is that velocity is specified, while the pressure gradient is zero. This boundary condition is applied to the top, bottom, left, right, and front surfaces of the domain region. For pressure outlet boundary conditions, velocity gradients are zero, while the pressure is defined. This boundary condition is applied in the back surface of the domain region. The wave model defines the volume of fractions on the boundary conditions, as well as velocities and pressures.

In both regions, hexahedral mesh topology is used. In the domain region, volume controls for mesh refinement are made in regions close to the helicopter, decreasing the density of elements according to the wake distance. The elements are refined in the waterline region to better represent the water–air interface. Volume control is used to refine the mesh within the possible trajectories of helicopter motions, maintaining similar element size to perform the overset interpolation.

The domain and helicopter meshes are presented in the bottom part of Fig. 3. The domain region has 6.90 M elements and the helicopter region, 1.02 M. Elements have sizes similar to the overset refine zone in the domain region, as a recommendation of Ref. [27] to interpolate these two meshes to generate an equivalent unified mesh. Recommended practices guidelines for ship CFD application recommend using no less than 40 grid points per wavelength on the free surface and, in irregular waves, using at least 20 grid points for the shortest wavelength [28]. This practice was taken into consideration to define the element size in the free surface.

The pressure and velocity fields are coupled in an implicit way using the semi-implicit pressure linked equations (SIMPLE) algorithm. The wave model imposes the wave elevation and wave kinematics on the inlet surface and the wave-induces pressures on the outlet surfaces, according to the pre-specified wave characteristics. The high-resolution interface capturing discretization scheme [29] is used to retain a sharp interface between water and air.

The Siemens star-ccm+ 2021 double-precision, a cell-center finite volume method software, is used for CFD simulations, including the mesh generator and the post-processing tools.

To discretize the helicopter into panels, the blades and the skid were not included. Taking it into account, the mass of the helicopter in the BEM simulation was defined as slightly lighter, due to lack of buoyancy from the skid volume. Figure 5 shows the mesh topology of the helicopter: 37,516 elements were used and the panels were conveniently positioned to match with the water surface at the hydrostatic condition. A standard wave grid resolution (⁠$81×51$⁠), also shown at the top of Fig. 5, was used, with a grid size factor of 2.

To solve the radiation and diffraction theory using BEM, the software ansys aqwa is adopted. A boundary integration approach is employed in aqwa to solve the fluid velocity potential governed by the above control conditions. In this approach, the frequency-domain pulsating Green’s function in finite depth water is introduced, which obeys the same linear free surface boundary condition, seabed condition, a and far-field radiation conditions. The water depth was set to a very high number (1000 m), so it can be considered as a deep water.

For the motion analysis of a ship in waves, most hydrodynamic forces on the hull can be computed using potential theory. However, roll damping is notably influenced by viscosity. Consequently, potential theory calculations tend to overestimate roll amplitude during resonance and lack accuracy [30]. Some corrections have been applied successfully to study the ship behavior, such as including a viscous damping coefficient to the roll motion. However, for the helicopter case, with a different geometry from a typical vessel, it should be further investigated the necessity to apply some corrections. The effect of viscosity in a non-vessel geometry is also assessed in Ref. [18]. Hence, it has opted to use numerical methods to estimate the natural periods of the roll, pitch, and heave, instead of using classical theoretical formulas typically applied for ships. To compare the effect of neglecting the viscous effect, these natural periods were obtained in CFD with and without viscous effects, and compared with BEM results.

For the CFD simulations, the helicopter was perturbed from the hydrostatic equilibrium in one direction, rotating it around the CG in the $x$⁠, $y$ directions, or sunk in the $z$ direction. Simulations were conducted in calm waters. In each motion, three different initial conditions were imposed: for the roll and pitch motions, an initial offset of equilibrium in 5 deg, 10 deg, and 15 deg; and the heave, $Δz=−0.20,−0.40$⁠, and $−0.60$ m. Also, in each case, a condition with inviscid and viscous flow models was taken into consideration, totalizing 18 simulations. Results for the roll, pitch, and heave motions are shown in Fig. 6.

All results have a period represented by more than 100 time-steps, following the recommendation of Ref. [28]. The estimation of natural periods is made by observing the peak time of each oscillation. The mean of the first five oscillation periods was taken into account. Table 3 summarizes the obtained natural periods for each model in CFD simulations.

Results for all three motions show a good agreement between the viscous and inviscid models, even for the roll motion. For vessels, it is common to apply viscous damping correction because the rolling motion generally does not produce much of a disturbance to the fluid. For slender (ship-like) bodies, the wave-making contribution is so small that the resistance of frictional drag on the hull surface plays an important role [31]. A helicopter with four floaters does not have a slender body like typical hulls, so the viscous contribution is not as pronounced as in ships. The results suggest that the viscous damping coefficient is not essential for BEM simulations due to the minor role that viscosity plays in these floating bodies. For the same reason in CFD simulations, the inviscid model is adopted instead of viscous one.

It is possible to estimate the natural periods using BEM by scanning a range of frequencies and solving the eigenvalue problem based on added mass and damping values, discarding modes that the frequency does not match the scanning frequency and retaining the matching modes as actual physical modes [22]. Three modes were found, one being associated with the roll motion, another with the pitch motion, and the third one, a combination of rolling, pitch, and heave motions. The average of each simulation in CFD and the obtained modes in BEM are shown in Table 4. BEM and CFD presented similar natural periods for the roll and pitch motions.

The RAO of the floating helicopter is obtained by using the hydrodynamic radiation/diffraction model of the bem software. To compute the RAO, 30 wave periods were simulated, from 0.8 s to 15 s, corresponding to a wavelength of $0.1L$ to $37L$⁠, a larger range than the minimum recommended by Ref. [32] (⁠$0.5Lpp$ to $2.0Lpp$⁠, in which $Lpp$ is the length between perpendiculars). The encounter waves vary from $μ=0deg$ to 180 deg, in 30 deg steep. Results of the roll, pitch, and heave RAO are shown in Fig. 7.

Results follow what was expected based on the natural periods obtained in Sec. 5. For the roll motion, it shows a peak of RAO close to the wave period shown in Table 4, 2.15 s, for every encounter wave angle $μ$⁠. As expected, the RAO is higher as close as the encounter wave angle is to 90 deg. For the pitch motion, it shows a peak of RAO also close to the natural period of 1.76 s. In this case, higher RAOs are associated with an encounter angle close to 0 deg or 180 deg. For heave motion, it does not have any values over 1.0. It is possible to observe small peaks in regions close to the natural period of 1.33 s. However, different from the other two analyzed RAOs, these peaks have a small impact on increasing the RAO. Indeed, for higher wave periods, the RAO is close to 1, regardless of the encounter wave angle.

The truncation of this polynomial “should be chosen so that the expected order of grid convergence can be replicated” [33]. Therefore, for second-order approximations, at least a number of grid representation $ng$ equals three are required. When using more than the minimum required number of grid representations, the coefficients $a$ are obtained from the least-squares fitting.

The refinement factor $Υ$ between each consecutive mesh level is around $2$⁠, as recommended by ITTC. Figure 8 shows the adopted meshes for the BEM simulations, and Fig. 9, for the CFD simulations. As shown, all three meshes present similar topologies in different refinements.

The analyzed condition is the helicopter in regular head waves, $μ=0deg$⁠, height of $H=2.894$ m, and period of $T=7.9$ s. The peak force in the $x$-direction $Fx$ is adopted as the quantity of interest to estimate the discretization error and uncertainty for CFD and BEM.

Figure 10 shows the results for the BEM and CFD. The results of extrapolation for the first and second-order convergences $S1$ and $S2$⁠, as described in Eqs. (17) and (18), respectively, are also shown. Table 5 shows the value for BEM simulations of elements for each mesh, time-steps, $Fx$ forces, and the uncertainties. Table 6 shows the results for the CFD simulations.

In both BEM and CFD cases, a low uncertainty is noted across all mesh levels. For the reference mesh, an uncertainty of $Uϕ$ of 0.59% for BEM (first order) and 2.1% for CFD (second order) were observed. Even when considering the order with higher uncertainties, the percentage does not exceed 6%. The results from both methods suggest that the adopted meshes and the associated time-steps are suitable for the proposed analysis of the helicopter in waves. Hence, further analysis in the following sections will proceed using the obtained combination of meshes and time-steps.

This section aims to evaluate the RAO by testing four main approaches: the RAO with BEM in the frequency-domain (the resulted values in Sec. 6); the response in the time-domain using BEM; and CFD. The helicopter was submitted to regular waves with a wave period of $T=7.9$ s, the mean wave period of the highest significant wave height condition, as described in Sec. 2.2. Also, it is considerably far from natural periods obtained in Sec. 5, which is around 2 s.

As it was intended to compare with the frequency-domain RAO, only the degrees-of-freedom of roll, pitch, and heave were allowed. Surge, sway, and yaw motions were blocked. In this way, it prevents the body from yawing, for example, which could compromise a fair comparison in the same situation for a long time running between the frequency domain and time domain.

The simulation was conducted by analyzing at least ten oscillations. Considering the wave period of $T=7.9$ s, simulations in both BEM and CFD time were 94.8 s, as the first two oscillations were ignored. The mean amplitude of the response, by observing the peak and valley values, was taken into account to determine the roll, pitch, or heave amplitudes.

In the following plots, BEM: frequency domain denotes the values obtained in Sec. 6, thus the slope of these lines indicates the RAO. BEM: time domain is the RAO obtained with the time-domain approach.

Figure 11 shows the comparison for the heave motion. Five different head wave amplitudes were tested, from $A=1.4$ m (⁠$H=2.8$ m) to $A=3.6$ m (⁠$H=7.2$ m).

Note that all points converge to the same slope, showing good linearity between the heave and wave amplitude. Indeed, as the natural period of heave is considerably smaller than the adopted wave period, a unit value of RAO was expected and confirmed in all results.

Figure 12 shows the comparison in pitch motion. The same five amplitudes were tested in head waves.

For the pitch motion, the frequency-domain RAO shows a value lower than the time-domain. When compared to CFD, the frequency-based RAO is closer than the actual response. However, all methods result in very close responses, which allows claiming that these differences can be associated with some inherent modeling errors.

Figure 13 shows the comparison of roll motion. In this case, BEM: frequency domain and BEM: time domain results do not lie close to each other.

Comparing the time domain, the results match with the CFD results for values up to $A=1.447$ m (⁠$H=0.894$ m). Results from $A=1.736$ and 2.084 m (⁠$H=3.472$ and 4.168 m) have a considerable discrepancy between CFD and BEM. Figure 14 shows the CFD results with the waterline in four different encounter wave heights and at the moment with the highest roll angle. Note that for results equal and lower to $H=2.894$ m, the waterline covers partially the floaters. When the wave height increases, water covers the left floaters.

The greater the waterline height on the left floaters, the larger the submerged volume of the floaters and subsequently, the greater the buoyancy forces. However, once the floaters are fully submerged, the buoyancy on the float remains constant regardless of the waterline height. Consequently, the roll moment induced by the buoyancy of the floaters also remains unchanged. With this limitation on the roll moment, the helicopter tends to roll in the direction of the left floaters. This phenomenon is accurately modeled using CFD. In the case of BEM, even when accounting for Froude–Krylov effects (forces introduced by the unsteady pressure field generated by undisturbed waves) in the time-domain results, it becomes clear that such considerations are insufficient to adequately represent the abrupt change in submerged body dynamics during roll rotation.

To further observe this difference in the roll motion between CFD and BEM, first the force in the $x$-direction is analyzed in the condition of head waves (⁠$μ=0deg$⁠), shown in Fig. 15. Five amplitudes were analyzed, as made before for the heave and pitch RAO plots.

A small difference is observed between BEM and CFD. Although this difference could lead to some positioning differences in the long term, the percentage difference between CFD and BEM is small.

Figure 16 shows the force in the $y$-direction, when submitted to regular beam waves, $μ=90deg$⁠. The same input wave conditions as observed for the roll motion in Fig. 13 are analyzed.

No differences were observed between the BEM results of the time domain. However, a notable fact is that CFD and BEM agree well for wave heights up to 2 m. Then, BEM starts to underpredict the force in the $y$-direction, when compared to the CFD result. This result is connected to the previously discussed point that the BEM method was underpredicting the Froude–Kyrlov forces to correctly model the roll moment, which led to small rolling and small $y$-forces.

To further investigate it, first it is analyzed the condition in which the methods have similar results. Figure 17 shows the pressure distribution in the maximum rolling angle, for $H=2.010$ m and $μ=90deg$⁠. All methods use the same color bar scale.

Indeed, the forces and amplitudes of the rolling angles from BEM and CFD are close. It is possible to observe the similarities between the pressure distribution on the bottom part of the helicopter: in all cases, it is observed higher pressure value on the right side, especially on the back of the forward-right floater and on the front of the backward-right floater. Similar pressure values are also found in the other floaters as well as closer values in the bottom of the fuselage. Therefore, with the matching of the results between the methods, BEM can be a viable option to get the pressure distribution for moderate conditions, in which the wave height is not too high.

Since the peak roll angle is different for each method in the $H=4.168$ m case, Fig. 18 shows the pressure distribution at zero roll angle with a positive angular velocity (rising rolling angle). In this case, it is possible to note a clear difference between BEM and CFD: in the first, there is a higher pressure in the bottom of the right floaters whereas, in the CFD, the left floaters are the ones that have the highest pressures. While the CFD presents the acceleration occurring in the orientation to push the left floaters up, the BEM forces still present the rolling acceleration to push the right floaters to the bottom.

Figure 19 shows the pressure distribution for the maximum roll angle observed for regular waves at $H=4.168$ m. Previously in Fig. 14 it was shown that, for this wave height, the water covers the floaters in CFD simulation. The front view of the pressure distribution confirms these results by showing high values of the pressure in the top of the left floaters, as a consequence of being submerged by the wave. On the other hand, results using BEM show that panels above the water level in the hydrostatic condition underestimated the pressure, a consequence of the low rolling angle condition, which did not submerge the floaters.

Due to this reason, it is not recommended to use BEM for results higher than $H=2.894$ m. For results below the critical angle of 3.64 deg, results match with CFD. However, it should be taken into consideration the limitations of BEM when simulating rougher wave conditions.

In conclusion, it was observed that the BEM approach gives good results when compared to CFD approach for certain ranges. A linear relationship was observed between the wave amplitude and pitch, heave, and partially the roll, for the analyzed wave period and heights. Roll motion higher than 3.64 deg should be avoided using BEM, as it disagrees considerably with the CFD results. If the BEM simulation is conducted below the limits here tested, it should lead to providing similar motion results as with CFD.

A short-term statistical analysis is conducted using BEM and CFD to study the motion of the helicopter based on real case scenarios. The adopted irregular wave scenarios were described in Sec. 2.2, choosing three different significant wave heights, $Hs=4$⁠, 5, and 6 m. The regulation recommends 5 min study after ditching in sea conditions: “the probability of capsizing in a 5-minute exposure is acceptably low in order to allow the occupants to leave the rotorcraft and enter life rafts” [5]. However, to ensure a better statistical representation, both BEM and CFD were simulated twice the recommendation, 600 s. Given the potential for different wave periods, two approaches of QTFs to capture the potential influence of second-order wave forces on the helicopter’s response were studied. The Newman approximation and the full QTF matrix were chosen for calculating the QTF matrix, balancing accuracy with computational efficiency.

As there is no connection to fix the structure, such as a mooring system, the $x$-movement (wave direction) is deactivated in BEM simulations. For the CFD simulation, the domain region moves to ensure that the helicopter region stays inside the designated zone to allow proper hole-cutting for the overset mesh method.

The yaw orientation is shown in Fig. 20, in which the left-hand side shows the results with CFD and the right-hand side, with BEM, by also comparing the effects of the full QTF matrix calculation. Plots are split between methods for better visualization. It should be noted, however, the difference of scales in the $y$-axis of each plot.

Results with CFD show that, after some time, all conditions oscillate around 90 deg orientation, indicating that the beam waves put the helicopter in a “stable” position. A notable difference between each wave condition is the time to reach the yaw orientation of 90 deg, going earlier as the higher the significant wave height is. However, given enough time, each of the analyzed cases reaches the beam wave orientation.

Results with BEM also demonstrate a tendency to oscillate around the 90 deg orientation. This is also valid for the $Hs=6$ m condition, as it oscillates around 450 deg (360 deg + 90 deg), making a complete round before oscillating around the beam wave condition. The full QTF matrix calculation yields the same results as with the Newman calculation for the cases of 4 and 5 m, being barely possible to visually distinguish each curve on the plot. The only notable different condition is at $Hs=6$ m, and after 200 s. Thus, results suggest that the more accurate calculation of the QTF matrix to provide a more accurate second-order effect has no significant influence on the yaw angle for this test case, except for some more extreme cases, in which high wave amplitudes make a difference in the drift forces, resulting in a different yaw orientation in the long term run. In all BEM cases, the yaw angle reached more than 300 deg before decreasing and stabilizing in the beam wave condition. The difference of forces in the $y$-direction for higher wave amplitudes, as shown before in Fig. 16, helps to understand the different behavior of the yaw curve differences: compared to BEM, the curve from CFD results has a more noisy shape, with higher frequency.

Figure 21 shows the roll angle distribution. As discussed in Sec. 7, BEM presented some difficulties in simulating roll angles higher than 3.64 deg (for $H=2.894$ m and $T=7.9$ s). This is also observed in the short-term statistical results with irregular waves: the distribution of the roll angles is smaller when compared to the CFD. Regardless of that, results show similar results for different forms of QTF matrix calculation, especially from $Hs=4$ and 5 m. For $Hs=6$ m, there is a small difference of values, observed by the curves between 0 deg and 10 deg. Indeed, the yaw results shown before in Fig. 20 indicate a similar behavior of $Hs$ = 4 and 5 m, and a slightly different for $Hs=6$ m. Between CFD results, it shows a close similarity of the results between $Hs=5$ and 6 m.

Observing the figure, there is a 48.75% distribution of the roll motion being between −3.64 deg and 3.64 deg with $Hs=4$ m condition; and 39.51%, in the same range for $Hs=6$ m. Thus, the considerable distribution of roll motion above 3.64 deg discourages the use of BEM, as discussed previously in Sec. 7.2 and Fig. 14, that for higher roll angle the method does not model correctly the effects of the submersion of the floaters.

To provide a deeper analysis of the roll angle, Fig. 22 shows the exceedance rate of the absolute roll angles. The peak (and valley, since it is considering the absolute value of roll angle) of the roll motion along the time series were computed and divided by the total simulation time, 600 s.

It is observed a discrepancy in exceedance rate between BEM and CFD, which was also expected, based on the previous distribution of roll angle. It shows that the BEM underpredicts the exceedance rate of the roll angle for all the analyzed test cases. Therefore, for the analysis in this condition of sea state 6, with the analyzed helicopter, using BEM can not be recommended because it underpredicts the roll motion, important for the analysis of green water and capsizing events. Based on the CFD results, the exceedance rate of the maximum roll peak was 23.1 deg, six times per hour (⁠$Hs=4$ m); 24.3 deg, six times per hour (⁠$Hs=5$ m); and 29.0 deg, six times per hour (⁠$Hs=6$ m).

Section 6 has shown that beam waves result in higher RAO of roll motion, compared to head waves. Figure 20 has shown that in both methods after a certain time, the helicopter tends to align perpendicularly (⁠$μ=90deg$⁠) with the waves. Therefore, these high roll angles are expected to occur when the helicopter is in the beam wave orientation. To analyze it, Fig. 23 shows a scatter distribution between the peak absolute roll angle and its respective yaw angle. BEM and CFD were plotted separately for better visualization. In this case, the $y$-axis scales are the same in both methods.

The results confirm that higher peak roll angles are only found when the yaw angles are close to 90 deg. There are small points between 0 deg and 90 deg for the CFD cases because the helicopter tends to reach 90 deg so, during the yaw rotation, there is no much exposure of the helicopter between these orientations. In both scenarios, a noticeable pattern occurs, indicating a limiting relationship between the yaw angle and its peak roll angle. This limit appears to exhibit a linear correlation, where higher yaw angles reach to increased maximum peak roll angles. Notably, the BEM plot demonstrates a shallower slope, indicating smaller peak roll angles as the yaw angle rises. Since the helicopter moves to the 90 deg, even when landing in $μ$ = 0 deg, and the yaw orientation reaches higher rolling motions, retarding the rotation is encouraged to keep the helicopter as much as possible aligned with the waves.

The computational time required for the BEM simulation is considerably lower than that for the CFD, which is a great appeal. The 600 s BEM simulation was run on a 6-core workstation and completed in hours, whereas the 600 s CFD simulation was run on a high-performance computing system with 4 nodes and 24 cores per node and completed in weeks. However, results show a flaw in using BEM for moderate wave heights, especially in beam waves. Therefore, as the interest is to assess the motion behavior in severe wave scenarios, using BEM is not promising to analyze the short-term statistic of the flotation stability of helicopters.

This study investigated the seakeeping behavior of helicopters under realistic scenarios, based on the Northern North Sea wave climate. CFD techniques, specifically the cell-centered finite volume method and BEM, were employed to assess motion responses and load distribution. Along this work, some conclusions could be inferred:

• Recommended practices in simulations were adopted to ensure the reliability of the numerical results;

• Based on the decay tests with roll, pitch, and heave motions, the inviscid model in CFD produced similar results when compared to the model with viscosity;

• For the typical size of a helicopter, a natural period of around 2 s was observed for the roll and pitch motions; and 1.3 s, for heave motion. Therefore, RAO responses around these wave frequencies have the highest values;

• Analysis with regular waves with a period of $T=7.9$ s shows a good linearity behavior between the incident wave amplitude and the response of pitch and heave. For roll, linearity was observed up to $ϕ=3.64deg$⁠, which is too small to study capsizing;

• In all irregular wave conditions in sea state 6, it was observed that after a certain time, the helicopter tends to align perpendicularly (⁠$μ=90deg$⁠) with the waves;

• As expected, the higher the significant wave height, the higher the peak roll angles. These peak roll angles occur when the helicopter is aligned perpetually;

• Based on the CFD results, the exceedance rate of the maximum roll peak was 23 deg, six times per hour (⁠$Hs=4$ m); 24 deg, six times per hour (⁠$Hs=5$ m); and 29 deg, six times per hour (⁠$Hs=6$ m).

Although BEM is considerably less computationally demanding, unfortunately, it can not ensure the accuracy of the results on applications that are usually intended. For most realistic scenarios, CFD is more suitable than BEM for the helicopter ditching case. BEM is not appropriate for capsizing analysis, being only suitable for specific moderate wave conditions.

For future works, it is planned to use CFD to assess the capsizing in other scenarios of the helicopter, such as the water inside the cabin or a failure in one of the floaters to get properly activated, affecting the buoyancy of the EFS. It is intended to validate the results in the future with an experimental campaign.

The authors thank Airbus Helicopters for providing the initial geometry and characteristics of the helicopter that were used to inspire the geometry adopted in this work.

This study is part of the research EvoS-BraWa, which is supported by the Bavarian Ministry of Economic Affairs, Regional Development and Energy (StMWi) within the Bavarian Aviation Research Programme (BayLu25), under Grant ROB-2-3410.20-04-10-36/BLU-2109-0036.

The authors gratefully acknowledge the computing time granted by the Center for Computational Sciences and Simulation (CCSS) of the Univeristy of Duisburg-Essen and provided on the supercomputer magnitUDE (DFG Grants INST 20876/209-1 FUGG, INST 20876/243-1 FUGG) at the Zentrum für Informations- und Mediendienste (ZIM).

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.