An investigation of coupling ship/rotor flowfield using steady and unsteady rotor methods

ABSTRACT A coupling numerical methodology has been developed using the Navier–Stokes solver for ship/rotor flowfield simulation during the helicopter shipboard launch. The steady rotor model (SRM) based on the momentum source approach and unsteady rotor model (URM) based on the moving overset mesh method are respectively employed to account for rotor operations. Studies of coupling ship/rotor flowfield on two typical ships highlight complex interactions over the flight deck. It is shown that the rotor-shipboard vortex interaction and recirculation zone are more serious for small-sized-deck destroyers, while on the straight-through-deck ship the helicopter operations are sensitive to asymmetric distribution of lateral velocity caused by the island. Qualitative and quantitative comparisons in terms of vorticity and velocity between SRM and URM solvers have been conducted. The results demonstrate that both methods can capture complex interactions well. The velocity difference is within 1 m/s in most flowfield areas, except the rotor-wake dominant region. With the consideration of computational efficiency, the momentum source approach could be used to effectively conduct the study of helicopter shipboard operations.


Introduction
The launch and recovery of helicopters in a naval environment are extremely challenging and dangerous tasks. Due to the bluff-body flow separation, the ship airwake developed at the flight deck contains a low Mach number, unstable separated shear layer and recirculation zone, and seriously interacts with the rotor wake as the helicopter is in close proximity to the ship. The large velocity gradient in the flowfield disturbs the aerodynamic environment around the helicopter and results in a sudden change of local aerodynamic force, which might adversely influence the stability and controllability of the helicopter. Therefore, it is of great significance to carry out research into helicopter/ship coupling flowfield characteristics for landing flight safety.
Extant research (He, Liu, & Tan, 2014;Syms, 2004;Wang, Gao, & Liu, 2015;Zan, 2005) regarding isolated ship flowfields has been widely conducted in the last two decades. Early numerical simulations by Reddy, Toffoletto, and Jones (2000) using FLUENT and the Simple Frigate Shape (SFS) validated the presence of the horseshoe vortex over the flight deck. Zhang, Xu, and Ball (2009) used the modified ship geometry (SFS2) and the Cobalt solver to demonstrate that the change of windover-deck (WOD) has an obvious effect on flowfield CONTACT Yongjie Shi shiyongjie@nuaa.edu.cn behaviors. A parametric study using landing helicopter assault (LHA) was conducted by Polsky and Bruner (2000), who showed that the flowfield structure is insensitive to the Reynolds number within a certain range of wind speed. Comparisons of numerical methods on airwake simulation using Reynolds-averaged Navier-Stokes (RANS), Detached Eddy Simulation (DES) and hybrid RANS-LES (Large Eddy Simulation) methodologies have been performed in recent studies (Forrest & Owen, 2010;Lawson, Crozon, Dehaze, Steijl, & Barakos, 2012;Muijden, Boelens, Vorst, & Gooden, 2013;Thornber, Starr, & Drikakis, 2010). The results illustrate that the RANS solver is capable of predicting the fundamental flow characteristics of ship airwake, although there are some deficiencies in detail simulations of turbulent fluctuation. One-way coupling calculations (Forrest, Owen, & Padfield, 2009;Forrest, Owen, Padfield, & Hodge, 2012;Lee, Sezeruzol, Horn, & Long, 2005;Roper, Owen, Padfield, & Hodge, 2006) refer to simulations that only assume the effect of the ship wake on the helicopter, as commonly used in previous analyses of shipboard landing. However, a large discrepancy has been found in the helicopter control inputs compared with the experimental data, due to the rotor-on-ship effect being ignored. Recently, the moving overset mesh method was employed by Lee and Silva (2010) to investigate the interactional aerodynamics of ship deck and rotorcraft and to compute the peak and total load on a hangar door. Oruc, Horn, Polsky, Shipman, and Erwin (2015) used the efficient and simple momentum source approach (Rajagopalan & Mathur, 1989), which considers the time-averaged effect of the rotor across the disk to conduct the numerical calculation of rotor/ship flowfield for the development of a virtual dynamic interface simulation. In the work of Crozon, Steijl, and Barakos (2014), the two aforementioned methods were employed over the Canadian Patrol Frigate to demonstrate the importance of coupling effect on the wake and rotor inflow. Although the momentum source approach has been used in numerical studies of ship/rotor flowfield, whether this method can accurately predict various complex interactions that occur during the helicopter shipboard launch, and the specific difference in velocity distribution between two methods, is not clear.
In this paper, the numerical method of ship/helicopter coupling simulation has been developed with the help of previous studies (Shi, Xu, & Wei, 2016;Wang, Zhao, Xu, & Ye, 2013;Zhao, Xu, & Zhao, 2005), including the steady rotor model (SRM) based on the momentum source approach and the unsteady rotor model (URM) based on the moving overset mesh method. Qualitative analysis of complex interactions that arise from the ship/rotor wake in terms of vorticity and velocity for different wind angles and shipboard landing locations was performed for two typical ships: small-sized-deck destroyer and straight-through-deck ship. In addition, quantitative comparison of the velocity distribution and computational efficiency between the two methods was employed to demonstrate the scope of application of the momentum source approach in simulating the coupling ship/rotor flowfield.

Computational Fluid Dynamics (CFD) solver
The in-house solver, Rotorcraft AeroDynamics and Aeroacoustics Solver (RADAS), is used to perform the numerical simulations of the coupling rotor/ship flowfield. The three-dimensional, unsteady RANS equations are employed as the governing equations, which can be written as: W is the vector of conserved variables, and F(W) and F v (W) are the convective and the viscous flux vectors, respectively. Rrepresents source term, and β is the logical variable for the rotor model. ρ, p and E are density, pressure and energy. q and q b are absolute and blade moving velocity, respectively. n represents vector, and i is the unit normal vector in the x, y and z direction.
The inviscid terms are computed using the secondorder upwind Roe scheme (Roe, 1981) and the viscous terms are computed using the second-order central difference. The dual-time stepping method is applied with the second-order lower-upper symmetric Gauss-Seidel scheme (Yoon & Jameson, 1988) to simulate the unsteady flow phenomenon. A two-equation k − ε turbulence model (Krishnamurty & Shyy, 1997) is used for the RANS closure. The Courant-Friedrichs-Lewy (CFL) number of 5 is set up in the solver. Validation of RADAS is demonstrated through a wide range of rotorcraft applications (Shi et al., 2016;Wang et al., 2013;Zhao et al., 2005); e.g. single rotor, dual-rotor interaction, and rotor/fuselage interaction.
To model the aerodynamic interaction stems from the helicopter rotor, two rotor models are employed. In Equation (1), β = 1 represents the SRM, and β = 0 is the URM. Details of the two models are provided in the following sections.

Unsteady rotor model
One solver adopted in this paper is an URM, which is based on moving overset mesh. As shown in Figure 1, the overset grid system contains several C-O or C-H type body-fitted grids for each blade and a large-scale fixed background grid. In the solution, each blade grid can move arbitrarily in the background grid to account for the complex blade motions, such as rotation, flap and pitch. The topological relationship between the blade nearbody grid and the background grid is established using the hole-cutting (Chiu & Meakin, 2013) and donor-cellidentifying process (Zhao et al., 2005), and the flowfield data between the two sets of grids are exchanged using a tri-linear interpolation scheme (Fan, Xu, & Shi, 2014).

Steady rotor model
To reduce the complexity of the simulation of the unsteady rotor flow, a steady rotor method based on the momentum source approach is used. In this method, the solid rotary blades are replaced by an infinitely thin actuator disc. Moreover, this method has been used in many studies on rotor flowfield and rotor/fuselage interaction (Chuiton, 2004;Lee & Kwon, 2002;O'Brien & Smith, 2013).
The unsteady airloads of the blade are equaled to the time-averaged momentum source term and added to the right side of the governing equation. Therefore, the source term, R, in Equation (1) can be written as The aerodynamic force, dF, which consists of lift and drag, at the representative blade element on the rotor is estimated by the blade element momentum theory in the blade coordinate system. The resultant lift, dL, and drag, dD, are where V r is the relative velocity at blade section r, b is blade chord and C l , C d represent the lift and drag coefficients, respectively. By applying Newton's third law, the instantaneous force that acts on the flow from the blade is represented as −dF. Then, the mean force required in the momentum source approach can be obtained through averaging the time-varying result over one revolution. Accounting for the effect of the number of blades, the expression of momentum source term, S, in Equation (2) can be written as where, φ is the azimuthal angle and N is the number of blades.

Numerical method validation
There is no available experiment data of coupling rotor/ship flow, thus the simulations of isolated ship and rotor/fuselage interaction are performed to demonstrate the validation of the developed method.
(1) Isolated SFS ship simulation The CFD solver is first performed using modification of a simple frigate ship (SFS2). The wind tunnel experiment was conducted at the National Research Council of Canada (Zhang et al., 2009) for the 1:100 SFS2 model. The contour map of velocity at four vertical slices along the X-axis -i.e. x/L = 0.25, 0.5, 0.75 and 1.0 -are shown in Figure 2, where the predicted and test data are plotted separately. With the blocking effect of the hangar, a large-scale recirculation zone is formed over the ship deck, which results in a reduction of horizontal velocity. In these planes, the velocity reaches minimum value at the central symmetry line, and then increases gradually outward until it returns to the far field value. As shown, the typical characteristic of flow over ship is well captured by the CFD solver.
A grid independence analysis is then carried out with three levels of meshes. The test simulations were performed for 0 degrees WOD at a freestream speed of 15 m/s, and the grids were denoted: coarse grid (2.7 × 10 6 cells), baseline grid (4.3 × 10 6 cells) and refined grid (8.5 × 10 6 cells). The near-wall boundary conditions (y+) in three calculations adopted the same setting. Figure 3 presents the horizontal velocity along a lateral line 10.66 m above the flight deck, with wind tunnel experimental data also provided for reference. The horizontal velocity magnitudes were normalized by freestream velocity. It can be clearly seen that the three sets of results are similar, despite the different grid cell sizes between them, and the increase in the number of grid cells did not provide significant improvement as compared with experimental data. Figure 4 shows a comparison of the vorticity in the longitudinal symmetry plane for the coarse grid, baseline grid and refined grid. As shown, three sets of numerical results were in good agreement. (2) Rotor/fuselage interaction simulation A numerical case of a simple rotor/fuselage model, which was designed and experimented at the Georgia Institute of Technology (Park, Nam, & Kwon, 2003), is then carried out to investigate the difference between SRM and URM for simulation of rotor-caused aerodynamic interaction. The coupling rotor/fuselage model operates at a low speed flight case, with Ma tip = 0.295, Ma = 0.0295, and rotor shaft angle of -6.0 degrees (positive tilted forward). Figure 5 shows the contour map of vorticity on a vertical slice crossing over the central axis of the fuselage. The pressure contours are superposed on the fuselage. As shown, the formation and evolution of the rotor tip vortex, as well as the collision with the fuselage, are captured by URM. For the steady solution using SRM, the pressure perturbation on the fuselage surface induced by rotor wake, which stems from the time-averaged airload, is also observed in Figure 5(b). Figure 6, depicting a comparison of the predicted time-averaged pressure distribution with test data at the top and side lines on the fuselage, indicates that the SRM is suitable for solving complex aerodynamic interference problems of rotor machines.

Numerical setup
The 1/2-scale simplified landing platform dock-17 (LPD-17) and full-scale LHA ship are chosen to conduct the numerical simulation as they represent the typical destroyer geometry and the straight-through deck structure, respectively ( Figure 7). The 'Dolphin' helicopter rotor is employed, and its basic parameters are shown in Table 1. The numerical simulations were conducted at several WOD conditions and rotor positions, with 15 m/s of wind. The effect on the fuselage in the helicopter/ship coupling flowfield was analyzed in a recent paper by our research group (Su, Shi, Xu, & Zong, 2017). The results show that the fuselage created a blockage of the rotor induced downwash; however, it barely influenced the main flow characteristics of helicopter/ship coupling flowfield -e.g. the recirculation zone and vortex-vortex interaction -due to its narrow configuration. Thus, the effect of fuselage is ignored in present study.
A rectangular computational domain (8L × 6L × 6L, where L is the ship length) was employed, as shown in Figure 8, and the sea and ship surface were specified as no-slip walls. For the numerical simulation of ship airwake, the near-wall boundary condition (y+) is usually within the range of 50 ∼ 500 (Forrest & Owen, 2010). In the calculation, y+ is set as 50 and the first wall unit value is set as 2 mm in order to meet the requirement of standard wall function in turbulence model calculation.  Figure 9 shows three cases of the grid for the LPD-17 ship: (a) isolated ship, (b) actuator disk and ship, (c) shipboard rotor and ship. These grids are generated by commercial software (ICEM). The isolated ship grid (a) was refined in the vicinity of the ship, resulting in 4.1 million cells. The actuator disk was located in the refined grid area above the flight deck (b), with the isolated actuator disk grid containing 1.3 million cells. In case (c), the shipboard rotor grid was embedded in the ship grid and the overall grid size was 6.7 million cells. The information exchanges between shipboard rotor grid and ship grid used the same method as that introduced in Section 2.2. The grid setting on the LHA ship is similar to that on the LPD-17 ship, and the refined grid for SRM and URM are 7.5 and 8.8 million cells, respectively.
The solution convergence was determined by the rotor thrust coefficient and the residual. When the residual of the average value of fluxes is less than 10e-5 and the variation of the absolute value of thrust coefficient is less than 1%, the calculation can be considered as   a convergent result. For the coupling rotor/ship simulations, the solution convergence can be satisfied after 5000 ∼ 7000 iterations in most cases. As a result, to ensure the robustness of the calculations, 6,000 iterations are usually used for SRM solvers for the convergence. For the unsteady simulation, 30 revolutions were necessary to obtain the fully developed flowfield. Considering that the frequency of rotor rotation is larger than that of the shedding vortex of the ship airwake, to ensure the stability of the calculation, the time step is determined according to the period of rotor rotation. The analysis of time sensitivity in the rotor flowfield research has been investigated in our previous studies (Shi et al., 2016;Wang et al., 2013;Zhao et al., 2005). Thus, the rotor was resolved using 360 time steps per revolution, equal to 0.00047619 seconds per time-step.

LPD-17
Figure 10 compares iso-vorticity maps between the isolated ship and the coupling rotor/ship flowfield, where the rotor operates at the condition of Ma tip = 0.643, θ = 12.82°. When the air flows over the ship, it separates at the corner of the hangar because of its bluff body configuration, and the separated airflow then continues to move downstream and reattaches at the deck, which forms a large recirculation region behind the hangar. Due to the low pressure in the recirculation region, the airflow at two sides of the ship deck moves upward and two contra-rotating vortices are formed. The rotor/ship coupling case shown in snapshots (b) and (c) in Figure  10 illustrates that the flowfield above the landing deck changes greatly as the helicopter approaches the ship. The presence of the rotor stimulates a violent downdraft at the back of the hangar, which makes the shedding vortex separated from the top of the superstructure quickly fall down to the deck surface, while the airflows at two sides of the ship deck are pushed further outward to the edge of the ship deck. Although the URM solver can obtain a more accurate coupling regarding the ship/rotor flowfield, it is at the cost of a great number of computing resources. Table 2 summarizes the requirement of computation time for the two methods. All sample calculations were performed using 24 processors in IBM X440 workstation (Xeon E5-2650 2.2 GHz, 96 GB RAM). As shown, a typical simulation requires approximately 10 hours and 140 hours for SRM and URM, respectively. This indicates that the momentum source approach dramatically reduces the computational time for ship/rotor simulation.
Detailed comparisons between the two methods are then investigated. Figure 11 shows contour maps of vorticity at the longitudinal symmetric plane (x/L = 0). The streamlines colored by velocity are also shown in the figure. When the helicopter approaches the ship deck, part of the rotor is involved in the recirculation zone.  With the effect of the rotor's sucking and blowing, the vertical velocity component of the flowfield behind the hangar is increased, and the airflow is further downward; thereby, the reattachment points at the deck move forward and result in the reduced size of the recirculation zone. Comparing Figures 11 (a) and (b), there are some differences in the intensity of the vorticity in certain flowfields between SRM and URM, which could be explained by the fact that the unsteady simulations using the URM solver present the results at the moments the four blades are located at a 45°, 135°, 225°and 315°azimuth angle, respectively. However, it can be clearly seen that the spatial morphology and angle of deflection of the flowfield using SRM display very similar characteristics to the results for URM. Figures 12 and 13 show the contours of vorticity and velocity distribution in four different slices along the x direction -i.e. p1, q1, r1, and s1, at 0°WOD. The strong concentrated vortex shed from the advancing and    retreating sides of the rotor interact with the shipboard vortex at both sides of the deck and result in the strong vortex-vortex interaction due to the same rotation direction of two kinds of vortex, which is more serious for a destroyer with a small-sized deck.
In order to explore the change in coupling flowfield caused by the variable height of the rotor, Figure  14 presents the contours of velocity distribution when the rotor is 5 m above the flight deck. As the rotor gets closer to the flight deck, it further involves in the recirculation zone, and the influence of the airwake onto the rotor is deemed stronger, which may considerably alter the aerodynamic loading of the rotor system. In addition, due to the reduction of the space below the rotor, the airwake over the flight deck is seriously squeezed by the rotor wake and the ground effect is expected to play a more important role in this case. Generally, the momentum source method shows good   A quantitative comparative study is then conducted to further illustrate the specific differences of the two methods. The monitoring points are distributed in three planes parallel to the ship deck. The schematic of monitoring distributions is shown in Figure 15.
Comparisons of the time-averaged velocity at three planes are shown in Figures 16-18 respectively. Because of the symmetrical distribution of velocity at headwind case, the monitoring points on two longitudinal lines (y/L = 0 and y/L = −0.5) are selected for comparison. The mean velocity of URM is obtained by averaging the instantaneous results over the last two revolutions. Figure 16 shows the mean velocity components of the monitoring points above the rotor (at a height of 12 m). Disturbance of the y-component velocity mainly arises from the induced effect of the rotor wake, so the velocity magnitude is fairly minor. For the x-component velocity, the magnitude on the line (y/L = 0) is less than that on the line (y/L = −0.5), particularly at the locations close to the hangar. This is can be explained by the fact that the ship tower stands in the longitudinal symmetry plane and leads to a blockage of the airflow, which causes a reduction of x-component velocity on the line (y/L = 0). The z-velocity consists of the rotor-induced downwash and the vertical component of the separated airflow, and is affected by the airflow separated from the ship tower in this case; thus, the maximal difference of zvelocity magnitude between two solvers appears on the line (y/L = 0). Figures 17 and 18 show the mean velocity components at the two monitoring planes below the rotor. The rotor location is also marked in the figures for illustration purposes. Compared to the preceding case, the velocity components at these points fluctuate significantly since they are exposed to the strong rotor-induced downwash. The discrepancy in the x-component and z-component magnitude on the line (y/L = 0) is evidently larger than that on the line (y/L = −0.5), while a different trend for y-component velocity is exhibited. This results from the fact that the rotor wake in the longitudinal symmetry plane, y/L = 0, has a more powerful effect on the  longitudinal and vertical velocity, while the effect of wake is canceled out due to the symmetrical distribution in the y-direction. Overall, the discrepancy of velocity magnitude is less than 1 m/s in most areas; only in the rotorwake dominant region, especially at x = −16, can the maximal discrepancy reach approximately 3 m/s. Figure 19, which maps the iso-vorticity for the coupling ship/rotor flowfield at 15°WOD, illustrates the distinct features compared to the 0°WOD condition. The airflow at the right side of flight deck deflects leftwards and couples with the rotor wake, resulting in the stronger vortex-vortex interaction. A low-pressure area is subsequently formed under the windward side of the rotor, making the rotor wake fall down toward the right edge of the flight deck, which can be seen on slice r1 as a plane crossing over the rotor hub center ( Figure 20); while at the left side of flight deck, the shipboard vortex is blown away from the deck by the crosswind and can barely interact with the rotor-induced downwash. It can be deduced that the rotor airloads are mainly affected by the crosswind, rather than the recirculation region in the large WOD conditions.   Figure 21 shows comparisons of the mean velocity components at a height of 12 m. The velocity field is asymmetric in the crosswind condition, but the variation tendency of velocity components shows the similar characteristics with the headwind case. In addition, the separation point of the airflow at the ship tower deflects leftwards; thus, the x-component magnitude gradually decreases from the starboard (y/L = 0.5) to the port (y/L = −0.5); and in the y direction, an offset of 4 m/s in velocity magnitude can be seen in the figure, which is roughly equal to the magnitude of the y-component of the free stream (V ∞ sin 15 • = 3.88). The z-component is greatly affected by the mutual interaction of the shipboard vortex and rotor wake, which causes the positive induced velocity on the left side of the line (y/L = 0.5). For quantitative comparison, the discrepancy of velocity components between SRM and URM solvers is within 1 m/s.

LHA
Following the detailed comparative analysis with the LPD-17, the solver is applied to the LHA ship, which has the typical structure of a straight-through deck. Figure 22 shows the contours of iso-vorticity at 0°W OD, where the rotor operates at the condition of Ma tip = 0.643, θ = 12.82°. As the free stream flows over the ship, the vortex shed from the ship bow passes through the flight deck, and subsequently breaks up before it can influence the rotor. Unlike the destroyer with a small-sized deck, the airflow separated from the corner of the island is mainly responsible for helicopter      operations, which can be seen in Figure 23 as a longitudinal plane crossing over the coupling ship/rotor flowfield. Some differences associated with the distribution of the shedding vortex are also found in the vorticity contours, especially in the region behind the rotor. This is because the simulation using URM estimates the unsteady results, and the instantaneous flowfield at a certain time is presented in the paper. Furthermore, it is noted that the straight-through deck of LHA enables multiple helicopters to take off simultaneously, which may bring about a serious problem wherein the rotor wake of each helicopter disturbs the rotor aerodynamic environment of others. Figure 24 shows the contours of velocity distribution in four different slices along the x direction, i.e. p2, q2, r2, and s2 (in sequence from ship bow to ship stern), respectively. An obvious reduction of inflow velocity at the left edge of the island can be seen on slices p2-r2, and the asymmetric velocity distribution might cause an imbalance of the rotor loading. For slice s2, it appears that the more intense aerodynamic influences on the velocity distribution below the rotor are captured by the URM solver.
Finally, the results of vorticity at 15°WOD are shown in Figures 25 and 26. Due to the effect of crosswind, the vortex shed from the ship bow is partially blown away, and the landing region is surrounded by a large amount of airflow separated from the left and upper edge of the island. The contours of velocity magnitude plotted in Figure 27 illustrate that a weak recirculation zone is developed at the corner of the island, resulting in an increased velocity gradient in a lateral direction. In all, the SRM and URM solvers both capture the main characteristics of the coupling flowfield well.

Conclusions
Ship/rotor coupling simulations with steady and unsteady rotor models have been conducted over two typical ships.
Quantitative and qualitative studies on vorticity and velocity between two methods have been performed to analyze the characteristics of coupling flowfield for different wind angles and shipboard landing locations. The following conclusions can be drawn: (1) The calculation results demonstrate that both the momentum source approach and moving overset mesh method can predict the complex interactions over the flight deck. Though the unsteady nature of coupling flowfield is well captured by the moving overset mesh method, this is at the cost of a great number of computing resources. For a typical case, the calculation time requires approximately 10 hours and 140 hours for SRM and URM, respectively.
(2) The quantitative analysis of velocity distribution between SRM and URM solvers shows the specific difference is within 1 m/s in most areas of the flowfield, except for the rotor-wake dominant region, indicating that the momentum source approach is capable of the analysis of helicopter shipboard operations and that the flow details around the rotor are not required. (3) For the typical destroyer, the longitudinal velocity perturbation that arises from the recirculation zone is mainly responsible for helicopter operations at 0°WOD, and the ground effect then increases gradually as the rotor gets closer to the flight deck. In crosswind conditions, the mutual interaction between windward shipboard airwake and rotor wake results in the increase of lateral velocity fluctuation, which adversely affects the helicopter's stability. (4) For the straight-through-deck ship, due to the presence of the island, the separated airflow and ground effect will result in the asymmetric distribution of lateral velocity, which is more serious in crosswinds because of the increase of lateral velocity gradient. In addition, the complex interaction of multiple helicopters taking off simultaneously over the straight-through deck is expected as a considerable problem, and will be considered as part of future work.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work was supported by National Natural Science Foundation of China: [Grant Number 11302103].