Simple multiple reference frame for high-order solution of hovering rotors with and without ground effect

In the present work, the aerodynamic performance of the Caradonna and Tung and S-76 in hover were investigated using a simpliﬁed concept of multiple reference frame (MRF) technique in the context of high-order Monotone Upstream Centred Scheme for Conservation Laws (MUSCL) cell-centred ﬁnite volume method. In the present methodology, the frame of reference is deﬁned at the solver level by a simple user input avoiding the use of mesh interface to handle the intersections between frames of reference. The calculations were made for both out-of-ground-effect (OGE) and in-ground-effect (IGE) cases and compared with experimental data in terms of pressure distribution, tip-vortex trajectory, vorticity contours and integrated thrust and torque. The predictions were obtained for several ground distances and collective pitch angle at tip Mach number of 0.6 and 0.892. © 2021 The Author(s). Published by Elsevier Masson SAS. This is an open access article under the CC BY license


Introduction
The helicopter is known by its wide versatility of flight condition with distinctive ability to take off and land from almost any place.Hovering near ground is one of the most important flight operations frequently encountered on urban mobility and emergency medical and rescue services.Even though this operation is found on a day-by-day routine, it involves a great risk factor, especially near coastal areas.On sand terrains such as beaches or deserts, the close proximity to the ground restrict the in-flight visibility, due to dust or sand on the air, this phenomenon is known as "brown-out".This interferes with the pilot's sense of spatial orientation and therefore increases the risk of accidents [1].On urban areas, the obstacle proximity is also another risk factor since the aerodynamic performance is heavily influenced by its surrounding.This is important not only on the design concept of the helicopter or drones but also on noise and safety of the surrounding.However, accurately predicting this vortex dominated flow-field and aerodynamics is a difficult task for most commercial Computational Fluid Dynamics (CFD) codes since higher-order resolution is almost requirement [2].
Until '80s, modelling hovering rotors was limited to analytical lifting line theory and prescribed wake models [3].At the time, the current state of CFD methods were insufficient in conserving the wake and computations of the rotor flowfield were based on Euler solvers formulated in a rotating reference frame using the Jameson's finite volume scheme with explicit Runge-Kutta time stepping [4] and Navier-Stokes coupled with an external wake model [5].The following decade was characterised by significant advances in computational power, numerical techniques and turbulence modelling.Srinivasan et al. [6] used a Reynolds Averaged Navier Stokes (RANS) solver using an implicit, upwind, finite difference scheme on a C-H mesh topology.The authors remarked the role of well meshed domains and need of dissipative numerical schemes to capture the tip vortices.Higher-order computations introduced later that decade demonstrated efficiency in capturing the vortex system up to several revolutions [7].However, even with this numerical advance, it still requiring a minimal mesh refinement in the path of the tip vortices, which is unknown beforehand.Traditionally, most of the CFD codes were mainly based on that Overset mesh methods to reproduce the blade movement [8,9].Since the meshes are mainly cartesian, implementing Adaptive Mesh Refinement (AMR) was a natural move to automatically refine and coarsen the mesh locally on those complicated vortex structures [10].The 00's was marked by the improvement of both higher-order schemes coupled with AMR techniques.Hariharan et al. [11] used a seventh-order Weighted Essentially Non-Oscillatory scheme (WENO) on a block-structured overset meshes to perform Out of Ground Effect (OGE) of hovering rotor.Postdam et al. [12] performed a AMR steady-state solution of the V-22 tiltrotor using the OVERFLOW code.The numerical dissipa-tion of the schemes demonstrates a direct impact on the prediction capabilities of blade-vortex interaction.Numerical simulations of rotors in-ground-effect (IGE) was only introduced on the following decade.Kutz et al. [13] carried out unsteady RANS simulations of a helicopter on hover and forward flight using low-order schemes.These computations were compared with experimental values showing the strong impact of the ground on the hover performance and the brown-out problem.Kalra et al. [14] studied the effect of different tip shapes on a hovering rotor near ground using a third-order Monotone Upstream-Centred Scheme for Conservation Laws (MUSCL) scheme on the blade mesh and fifth-order WENO resolution on the background mesh.Detailing comparisons regarding the tip vortex trajectory and core size were made against experimental results.Hariharan et al. [15] performed hover simulations of the Caradonna-Tung rotor at several radii above the ground using a hybrid overset mesh-based solver with fifth-order accuracy.More recently, Wang and Zhao [16] presented a design optimization by changing the twist and chord length to improve the aerodynamic performance of an unmanned helicopter.Hwang et al. [17] assessed the ground effect of the S-76 hovering rotor for three blades geometries using overset unstructured mesh approach and high-order WENO reconstruction in the off-body regions.Hu et al. [18] coupled CFD with a discrete element method to investigate the dust cloud and the motion of sediment particles with the ground effect.Stokkermans et al. [19] uses high-fidelity CFD to study aerodynamic interactions between rotors of a compound helicopter.Pasquali et al. [20] presented a numerical and experimental correlation to predict the rotor aerodynamic in the presence of close ground.A research effort has been made on using this dissipative numerical techniques to understand the wake-breakdown phenomena and braid secondary vortex [21,22] The use of high-order schemes on block structured blocks has promoted the development of overset meshes to assess rotary wings problems.This powerful combination established this technique as standard to handle the relative motion of the blade.Chen-Long et al. [23] assessed the capability of the Overset and other methods such as sliding-mesh and moving reference frame to perform a low-order simulation of the hovering Caradonna and Tung rotor.These methods were compared in terms of velocity and pressure field and differences were pointed in terms of accuracy, mesh and convergence time.It is highlighted that sliding-mesh and moving reference frame are slightly less accurate than Overset in capturing shock waves, but they require considerable lower computational cost.Steijl et al. [24] presented the third-order MUSCL sliding-mesh work on a block structured mesh to simulate a rotor-fuselage interaction.This technique has been successfully used in LES simulations of tidal-stream turbines, presenting great agreement with experimental power and thrust predictions [25].Ramirez et al., [26] presented a new technique to preserve the high-order stencil across the sliding interfaces.The moving reference frame simple which avoid the mesh movement by using a non-inertial mathematical formulation.It has been widely used in turbo-machinery [27] and rotary wings [4] applications.Zhong et al. [28] presented a block incomplete lower-upper decomposition on C-H structured inviscid and Navier-Stokes solver written in rotating reference frame using the Osher's approximate Riemann solver and MUSCL reconstruction.The findings demonstrated enhanced convergence acceleration with larger time steps.In a Multiple Reference Frame (MRF), the governing equations are switched between frames enabling to simulate two or more domains which may be either stationary or in motion.The connection between these regions is handled by a mesh interface created at the Computer-aided design (CAD) level.This surface is required exclusively for numerical purpose and has no physical meaning.It is mainly responsible for exchanging data between the set of equations for these regions.One of main drawbacks of the MRF approach, is the interface location definition [29].The position of the interface is relevant and when misplaced may lead to a time consuming procedure of changing the geometry and re-meshing.Moreover, meshing algorithms may struggle to maintain a good mesh quality when an interface is placed near the blade.In order to avoid this non-physical mesh interface, Remaki et al. [30] proposed a simplified technique to virtually split the domain.This approach was implemented o on the second-order cell-vertex based CFD code, SU 2 , and compared with experimental and commercial software for 2D and 3D cases.The algorithm was also compared with the classical MRF multi-zone mesh and it shows a perfect agreement between them.The steady-state flow assumption on MRF simulations is also very attractive for industrial applications and has shown to be accurate for several industrial problems such as rotor-stator in turbomachinery [31], fans [32], wind turbines [33] and stirred tanks [34].
In the present study, the virtual interface algorithm concept is extended to a high-order unstructured cell-centred solver.Solutions are obtained on hybrid (mixed-element) three dimensional unstructured meshes using the Spalart-Allmaras turbulence model with rotational correction.The governing equations are discretized using the absolute velocity formulation, the Harten-Lax-van-Leer-Contact (HLLC) Riemann solver [35], the spatial discretisation employs the MUSCL-MOGE variant limiter of Tsoutsanis [36].The solution is advanced in time is with a Lower-Upper Symmetric Gauss-Seidel (LU-SGS) implicit backward Euler time integration unless otherwise stated.The CFD solver labelled Unstructured Compressible Navier-Stokes 3D (UCNS3D) [37] which is validated, assessed and evaluated across a wide range of flow conditions [38][39][40][41][42][43][44][45][46][47][48][49].The proposed algorithm is used to simulate viscous flow around the Caradonna and Tung [50] and S-76 [51] rotor in hover with and without the ground effect.The Caradonna and Tung [50] case was chosen to validate the predicted results of pressure distribution and vortex trajectory and to investigate the thrust ratio and vortex system for three different ground distances.The size of the rotational reference frame and the impact of the mesh resolution on the wake breakdown are discussed.The OGE and IGE predicted results of the rotor aerodynamic performance of S-76 [51] rotor are also compared with experimental data in terms of collective blade pitch, thrust and torque.

Governing equations
In order to avoid any relative motions between parts of the computational domain, the frame of reference is decomposed into two parts: rotational and stationary.This procedure is known as Multiple Reference Frame (MRF) and works by switching the observer view of the problem.By defining an axisymmetric subdomain with a rotating reference frame attached to the blade at a particular position, an unsteady problem can be reformulated as steady.Then, the governing equation in each subdomain is computed in a different manner.A useful mathematical approach to deal with more than one frame of reference is using the absolute velocity formulation, as it does not require special reference transformation at the interface between the subdomains.In this sense, the governing equations at each subdomain are computed with respect to its reference frame but the velocities are stored in the absolute frame.This allows a precise computation of the fluxes on non uniform meshes [4] and is widely used on turbomachinery applications.For each frame of reference the equations are recast and solved in terms of absolute variables, thereby the flow in the farfield is uniform but not the relative flow.The relative velocity components u r , v r and w r and absolute velocity components u, v and w are related by the rotational velocity components u , v and w in the Cartesian coordinates x, y and z as follows: ( For constant clockwise angular velocity in all components 1 , 2 and 3 and the radius, r, as the distance of the element from the axis origin, the rotational velocity yields ( The three-dimensional, compressible Reynolds Averaged Navier-Stokes (RANS) equations are formulated in the inertial reference frame in the conservative form can be written as: where U is the vector of the conserved mean flow variables and turbulence model variable, x denotes the coordinate of a point of the domain, F c and F v are the convective and viscous flux vectors, respectively, and S is the source term vector due to the inertia forces and rotational corrected version of the Spalart-Allmaras (SA) turbulence model [52].
Using the absolute velocity formulation the equation (3) in the rotating reference frame yields: where ρ is the density, p is the pressure, ν is the turbulent viscosity SA working variable.Assuming a calorically perfect gas, the total energy per unit volume is calculated by E = p/(γ − 1) + 0.5ρ(u 2 + v 2 + w 2 ), where γ = 1.4 is the ratio of specific heats for air at normal atmospheric condition.The kinetic laminar viscosity, ν l , is computed by the Sutherland's law relating the dynamic viscosity, (ν l = μ l /ρ) to the ideal gas temperature.The kinematic turbulent viscosity, ν t is approximated by the SA turbulence model, the implementation is comprehensively described in [53].
The domain decomposition of the frame zones can be performed at two stages of CFD simulation: CAD or CFD solver.The former splits the computational domain into two mesh blocks using a interface boundary condition to perform the variable exchange and interpolation.In this sense, these two parts are independently created.In the latter approach, the domain is virtually decomposed by a geometry definition.Thereby, the elements of each part are identified on the solver and then treated according to its frame.

Virtual multiple reference frame
This approach employs a shape geometry formulation on the solver algorithm to identify in which reference frame the element is located.The rotational zone can be defined as cylinder by providing two points tridimensional coordinates ( P 1 and P 2 ), radius (r) and rotational velocity (ω), as shown in Fig. 1.Using this information we can identify whether the element centre ( P e ) frame of reference; Let the rotating origin, P o , defined as the cylinder centre as ( P o = (P 1 + P 2 )/2) and the rotating axis as −−→ The distance between P e and P o is P o P e and its vector −−→ ) r e = abs(P o P e sin(θ )), (7) the following criteria define at which frame the element lies: Fig. 2 shows in detail a schematic representation of the element identification and reconstruction on the virtual interface framework.The element is treated in the rotational frame of reference (grey-filled) if its cell-centre (blue box) locates inside the virtual interface (red dashed line).Even though some elements have edges inside the virtual interface, element j for instance, they are not considered on the rotational frame, since their cell-centre is out of this zone.In order to improve the high-order reconstruction, the rotational velocity is computed at each Gaussian quadrature points (white box) taking into account its relative distance to the rotating axis ( −−→ P 1 P 2 ), not the cell-centre position (blue box).The main reason for this choice is the gain on precision when computing the velocity contravariant and fluxes, especially in large tetrahedral elements.Two additional computation tasks are performed on the calculation process to change the frame from inertial frame to rotational: inclusion of extra source term and volume flux correction It is worth mentioning that this virtual decomposition is easy to implement on open source codes such as UCNS3D and can be easily expanded for multi-rotor problems such as unmanned aerial vehicles, drones and main rotor-tail interaction.A similar technique is also implemented on the 2nd-order cell-vertex compressible solver, SU 2 [30].In terms of parallel computing, this procedure and allocation of all additional variables is applied in each processor.Hence, the virtual decomposition is decoupled with the Message Passing Interface.For more information regarding the domain decomposition and parallel computations the reader is referred to [43].

Numerical framework
The computational domain is discretized on three dimensional shaped elements of type hexahedra, tetrahedra, prisms and pyramids.The element has a volume |V i |.An ordinary differential equation is achieved by integrating the equation ( 3) over the unstructured mesh using the finite volume formulation, yields (9) where U i and S i is a volume average of the conserved variable and source term vectors, respectively, within the control volume of element i. Summation limits, N f correspond to the number of adjacent faces per element and N qp is the number of quadrature points employed to approximate the surface integrals.|A l | is the surface area of the corresponding face, and α correspond to the Gaussian points, with coordinate vector x α and weights, ω α , over the adjacent face.Both weight and distribution depend on the order of the Gaussian quadrature rule applied.Thereby, higher integration rule improves the flux approximation.The interface fluxes are computed using extrapolated values, which are obtained from a polynomial reconstruction from the element-averaged data.

Spatial discretisation
The reconstruction algorithm is based on the approach given by Tsoutsanis et al. [54] and involves the decomposing of each element into unit tetrahedra.This procedure transforms the system of equations from physical space to a reference space, this is performed to reduce the scaling effect inherent in stencil with different element size.The MUSCL method has it origins with the pioneering work of Van Leer [55][56][57] and Kolgan [58,59], Barth and Jepersen [60] extended the scheme by applying a limiter to reduce the slopes of high gradient regions and preserve monotonicity of the solution, so no new local extrema are created.
All polynomials, for all the faces and for all the quadrature points are then limited to prevent any spurious oscillations from contaminating the solution.This slope limiter requires the minimum and maximum values from the stencil formed by the considered cell i and the direct side neighbours.The MUSCL-MOGE variant employed in this work the minimum and maximum values are determined not only from the direct side neighbours but from the entire stencil neighbourhood as detailed in [36].

Fluxes
The Riemann problem is solved with the approximate Harten-Lax-van-Leer-Contact (HLLC) solver [35] for the convective fluxes.This is also used to deal with the convective part of turbulence transport equation.Taking into account the rotational velocity invariant, Û , on the wave speed, S, for right and left state and its reconstructed solution, the flux function yields: (10) where Ŵ * ± is calculated for the element itself or its neighbour "+".The approximations of the wave speeds, S ± and S * , are iteratively calculated as explained in [61].The rotational velocity invariant, Û , act as a flux correction term on the rotation zone to compute the relative velocity as it is on the absolute formulation in equation ( 9).This correction is performed only on cells in the rotating zone.Thereby, in each domain, the governing equations are expressed in their relative reference-frame.
For the evaluation of the viscous fluxes the unlimited k-exact least square reconstruction is used for the gradients and they are then averaged from two discontinuous states as detailed in [54,62].
For the gradients additionally penalty terms are included following the formulation of Gassner et al. [63] for suppressing odd-even decoupling modes in the numerical solutions [64], in the following manner: (12) where L int is the distance between the cell centres of adjacent cells, and α = 4/3 similarly to previous approaches [64,65].

Temporal discretisation
The implicit approach is much desired in steady state solutions of high speed flows, especially dealing with complex geometries in RANS problems which requires a small element to ensure a y + lower than 1.In order to accelerate the solution towards steady state, the time step is computed locally taking into account the absolute velocity and the minimum edge at each element.Thereby, the solution advances in time with the maximum permissible time step of the control volume.
Since the goal is a steady state flow, the time step is computed locally taking into account the rotational velocity and the minimum edge of each element.In such a technique, the solution advances in each control volume with maximum permissible time step and thereby accelerate the convergence.The method marches in time using the Lower-Upper Symmetric Gauss-Seidel (LU-SGS) and implicit backward Euler time integration schemes, as it presents advantages in terms of parallelisation and computational cost.The finite volume formulation shown in equation ( 9) is rewritten the semi-discrete form where R i is the right-hand side residual of the conservative equation, which convergence approach the machine precision.Employing the first-order backward Euler implicit time stepping scheme, equation (13) yields Linearising in time, equation (14) becomes where R i tends to be equal to zero and ∂U is the flux Jacobian, which contains the linearisation of convective and viscous flux vectors and the source terms.Thereby, equation (15) gives where I stands for the identity matrix.Once the linear system is solved, the solution at each element i is simply updated as The residual is linearised by a first order approximation of the numerical flux using the Rusanov flux function as (17) where the convective and viscous eigenvalues are written respectively where n ij is the normal vector to the element interface, V ij is the contravariant absolute velocity at the face j on the cell i and a ij is the speed of sound.
The linearisation of the Rusanov flux of equation (17) gives the block diagonal and off-diagonal elements as: The linear system of the form A X = B in equation ( 16) is solved by its factorisation into three parts given by the following equation: where U * is the intermediate state and the lower, diagonal and upper operators are written as in which ∂S ∂U is the contribution of the source Jacobian term and it yields: and it is added directly on the diagonal component to keep its dominance [66].The diagonal elements of the matrix are stored and inverted directly and the off-diagonal ones are calculated at every stage.No numerical test were performed to evaluate the contribution of this term on the diagonal dominance of the matrix, however its inclusion has shown desirable convergence properties [67].

Results
In order to provide a strong level of correlation between the prediction and experimental, the accuracy of the developed framework is evaluated on two distinct hovering rotorcraft flow cases:  Caradonna and Tung [68] and S-76 [51].This choice was mainly based on the geometry information available and flow properties.Both cases have experimental data without the fuselage.It is well know that fuselage has a non-negligible effect on the main rotor.Including the fuselage modelling significantly increases the complexity of the rotor analyses in terms of geometry fidelity and mesh generation.Second, this experimental dataset covers a great range flow regimes, in terms of Mach tip velocity, angle of attack and ground effect.For a full description of the geometry and details of the experimental test data, the reader should refer to [68] and S76.

Rotational domain size influence
After a consistent refinement study on the blade, the computational mesh of 2.5M element was selected to perform a study of the impact of the rotational subzone size on the solution.The specifications of the meshes are shown in Table 1.The computational domain has a cylindrical shape with 35 radii in each direction, in which the rotor is placed at the centre.An angular velocity of 265.9 rad/s is applied to rotational domain, and at the farfield at the outer surfaces.The no-slip adiabatic wall boundary conditions are applied to the blade.Quadrilateral elements were used on the blade wall due to its benefits on capturing the surface preserving its curvature near the leading edge.These elements grow at 1.2 ratio on layers normal to their base forming a dominant hexahedraltype elements.The first wall cell distance of 2 × 10 −6 m is set.
The volume elements on the rest of domain were generated by the advancing front mesh algorithm [69].It is important to state that no local refinement was applied on this mesh, and larger elements are expected away from the blade.Thereby as the rotational domain increases, it would be poorly defined and less resemble a revoluted surface, in this case a cylinder.Three different sizes were assessed in this analysis: Small, Baseline and Big, as shows Table 2.The main difference among them is the length and radius of the rotational domain.
The effect of the domain size on the computed pressure coefficient and velocity profile at the interface, Figs. 4 and 3 respectively.There is no significant difference related to the rotational frame size.However, this region must be treated carefully as it impacts the numerical stability of the solver.Regions far from the rotor are mainly populated with large tetrahedral elements, especially using advancing front and Delaunay meshing algorithms.Large growth ratio on elements at the virtual interface damage the stability of the iterative solver.We believe this numerical detriment is mainly due to the high stiffness of the equation at large jump in element size across neighbours.This impact not only the gradient calculation at this location but also the linear solver.Since it also detriment the diagonal dominance of linear system.Numerical tests using hexahedral dominants meshing algorithms, such as Voxel, indicate an improvement on the interface definition and the stability of the linear solver, but it significantly increases the number of cells.In terms of wake physics capturing, previous work concerning the rotational domain size, which uses a mesh interface to handle the reference frame regions, suggests keeping the domain large contributes to decreased numerical errors due to interpolation exchange between the domains [29].However, other numerical works using a VMRF technique on cell-vertex discretization show good agreement with both multi-zone MRF and experimental data [30], even in cases within small gap between rotor and stator.

Wake Mesh refinement
The rotorcraft physics is clearly recognised by its particular strong vortical structures which develop at the tip of the blade.At each rotation, the blade interacts with the vortex generated by its previous passage.This phenomenon, known as Blade-Vortex Interaction (BVI), is one of the most complex feature of hovering flows and highly sensitive to the numerical approach and spatial resolution [70].It also significantly affects the rotor aerodynamic performance.An accurate computation of the wake trajectory and strength is much desired and challenging in rotor flow analysis.Dealing with such a high diffusive physics requires a local mesh refinement along the tip vortex trajectory, which are not known before the simulation starts.Automatic mesh refinement techniques are much desired but require a complex element decomposition in the solver computation.Typically, second-order schemes would require 20 to 30 points across the core of the vortex to numerically   convect the vortices without being largely dissipated [15].Higherorders schemes, on other hand, are quite useful to capture these highly diffusive regions but still needing some degree of local refinement.
A refinement study on both blade and wake was performed using only low-order schemes to assess the multiple frame capability in capture the wake and then the mesh is selected to higher order computations.For all simulations, the rotational zone has set according to the baseline size showed on Table 2. Through a consistent refining process on the blade and wake, three mesh levels was generated ranging from 2.5 × 10 6 to 15.4 × 10 6 elements, as tabulated in Table 3.A local refinement is performed at the tip blade radial station (r/R = 1) at near wake, covering up to 160 chords downstream the rotor.Fig. 5 shows the velocity magnitude contour at the near wake on those meshes.It is clear the improvement on capturing the wake velocity magnitude decay as the element size decreases.It is also noted that the virtual interface has no significant effect on the wake profile.As result, the Fine mesh was down-selected to perform the higher-order computations.
Fig. 6 shows the effect of the numerical-scheme order on capturing the pressure profile at the tip blade, r/R = 0.96.The flow at this section is highly influenced by the tip vortex and characterised by a shock at x/c ≈ 0.25.The impact of the numerical dissipation on the pressure profile is clear.The low-order predictions struggle to capture sharp pressure drop and its recovery, presenting small spurious oscillations at 0.35 < x/c < 0.6.While on the higher-order solutions are in good agreement with the experiment through entire span section.Fig. 7 shows the predicted tip vortex trajectories of the computed solutions compared with the experiment hot-wires measure- ments [68] and prescribed wake models [3].In comparison to the experimental data, the vortex radial contraction, Fig. 7(a), is well predicted by computations up to the first blade passage.After this point the computed solutions overpredict the vortex contraction.This deviation may be associated the effect of the rotor hub which is neglected in the present computational model.The fourth-order scheme performs slightly better than others, however the effect on the numerical scheme on the radial contraction is negligible.Overall, the computed solutions of the vortex contraction agree with the prescribed Kocurek wake models [3].Fig. 7(b) it can be seen that the CFD results accurately predict the slow convection of the tip vortices seen in the vertical displacement ( Z /R) up to a full cycle.
Fig. 8 shows the helical wake structure of the tip vortex using the same level of iso-vorticity contours.It is clear the advantages of using higher-order schemes on convecting the vortex structure downstream the rotor.It improves the predictions of the wake sheet and vortex tip flow phenomena, preserving the helical vortex filaments up to several radii downstream the rotor.Looking at the second-order prediction, the wake system evolves downstream the rotor before being consumed by the ringvortex.It is noticeable how this scheme struggles to deal with the large amount of diffusion.This numerical nature contributes to the vortex system settle down and achieves this stable ring state structure as seen on the couple radii downstream the rotor.As the wake is transported downwards, it can be clearly seen the interaction primary and secondary structures, which stretch between the tip vortices making an S-shaped path also seen experimentally [71].The presence of the vortex-ring close to the rotor blade contributes to a larger induced flow and under predictions of thrust coefficient [15].On higher-order calculations, the ring vortex persists on the converged solution, but it is pushed further down.As we increase the scheme-order, it becomes more evident how the helical system convects down and breakdown toroidal vortex into smaller scales.It can be seen that the horseshoe vortices once present on low-order scheme strain into vortex tube as the numerical dissipation decreases.This behaviour is also seen in similar works regarding the influence of the local mesh refinement on the vortex wake breakdown [22].It is worth to highlight capturing in detail this wake breakdown is challenging for unstructured solver without an automatic refinement meshing algorithm.Until recently this secondary vortex production has been seen as overemphasized by high-order solutions and not believe to be an accurate physical phenomenon, but state-ofthe-art volumetric flow measurement has confirmed its through the Lagrangian volumetric PIV variant technique [71].The reason for these phenomena are not known at the time and research concerning computed wake breakdown physics still ongoing.For a detailed overview of the recent advancements the reader must refer to recent papers published at the AIAA Hover Workshop [21,72,73].
Fig. 9 shows in detail the iso-vortex surface of the fourth-order MUSCL simulation coloured by vorticity magnitude detailing the mesh at resolution at the interface.An interesting phenomenon is seen regarding the helicity of the tip vortex on the wake.It is not clear if helicity of the vortex system it is damaged by unstructured mesh stretching [21] or interface between the frames.In order to investigate the nature of this issue, an additional comparison with single rotating reference frame simulation was carried out using the same mesh on second-order MUSCL as can be seen in Fig. 10.It is noteworthy that there is no significant difference on the preservation of the tip vortex trajectory and structure between the MRF and the SRF simulations.However, it is interesting that the wake breakdown and ring vortex are captured only on the MRF.This issue may be related to a number of factors that still on development research such as turbulent energy transfer [74], turbulence modelling, stretched mesh [21] and local refinement [22].

In-ground effect 4.3.1. Caradonna-Tung rotor
In order to assess the efficiency of the Multiple Reference Frame near wall surfaces, the hovering Caradonna and Tung rotor was placed near ground at three different distances ( Z /R = 0.5, 1.0, 3.0) as shown in Table 4 and Fig. 11.The local refinement was applied on the wake region.Inviscid wall conditions was set at the ground in order to avoid an excessive number of elements.The rotational zone was kept as the same size as showed previously on the wake refinement section.A Courant-Friedricks-Lewy (CFL) number is fixed to 20, which is computed based on the smallest edge in the spatial domain.As in the previous solutions, the simulation running strategy is tailored for convergence acceleration; the second-order solution is used to initialise third-order simulations and subsequently the later is used for the fourth-order.Fig. 12 shows iso-vorticity magnitude surfaces at the three different height above the ground.At the blade root, a strong upwash flow is observed on the IGE computations, in contrast to the OGE.The vortex-ring once present not far from the rotor in OGE simulations, on IGE it spreads due to the ground presence.The wake impinges on the ground and rolls up due to the influence of the incoming flow.The ground acts as a barrier constrain the flow in the axial direction, then it expands radially as it reaches the ground.This contact formation of a well-defined ring-like vortex structure under the rotor is known as ground vortex.It is clear the improvement of the higher-order schemes on capturing the development of the ground vortex and its breakdown into the braid vortex near ground.The dissipative nature of these schemes provides in more detail the information about the formation of intense vortex structures due to the proximity to the ground, such as the horseshoe vortex mixing in many directions and scales.This complex wake system would also contribute to explain the dust transport mixing when the rotor is hovering/lading near dusty-ground terrains.Note that the ring-like structure expands in size as it reaches the ground.For the IGE calculations of Z /R equal to 0.5, 1.0 and 1.25 it touches the ground at 1.3, 1.11 and 1.07 radii, respectively.It should be highlighted that capturing these structures are challenging for unstructured meshes.It also should be considered that the free-slip condition applied at the ground surface may cause a significant effect on promoting this radial expansion.Results of other numerical works also shown that the ground vortex is at least three times as strong as the tip vortices shed by the rotor [75].Fig. 13 plots the near-ground to out-ground thrust ratio at the three rotor positions computed here and compares it with historical data, image-vortex theory and other numerical predictions (data collected from [15]).
Fig. 14 shows slice view vorticity and velocity magnitude contours for the Z /R = 0.5, respectively, detailing the mesh of the virtual interface.Through all IGE solutions, it is noted that this interface has no significant effect on the velocity magnitude of the flow as it remains continuous at the interface.Concerning the vorticity magnitude, there is small effect on the helicity of the blade tip vortex as it move from the rotational frame to the stationary frame.The tip blade vortex strength seems to be more easily dissipated on the stationary frame.As this strong core structure passes through the interface, it seems to be dissipated at faster decay.This behaviour may be associated to gradient-diffusion approximation of the turbulence model on the rotational and inertial frame.Theoretical analysis suggests differences due to the rapid turbulent energy transfer rate due to the rotation and its relation with the helicity, an axial energy flux correction in rotating turbulence modelling has been recently proposed by [74].
To investigate the effect of the ground on the rotor wake, the predicted tip vortex trajectories are compared in Fig. 15 in terms of vortex age.The computations of the fourth-order accurate MUSCL are presented for radial and vertical positions for a range of the azimuth angle from zero to 400 degrees.For the radial contraction, Fig. 15(a), the results of the OGE simulation show great agreement with prescribed wake-models of Kocurek-Tangler [3].Taking into    [76] and other numerical works [15].account the ground effect, the IGE and OGE are quite similar up to the first blade passage at 180 degrees, in which the tip vortex contract radially.After this point the slope for IGE simulations changes considerable, especially for Z /R = 0.5, as the tip vortices start to expand as it approaches the ground.In terms of vertical displacement (z/R) of the tip vortices seen in Fig. 15(b), apart from small fluctuations for Z /R = 0.5 there is no significant differences between OGE and IGE calculations.Up to the first passage a slow convection is seen in vertical displacement.With the presence of the other tip vortex from the previous blade, the downwash decay of the vortex system appears to be faster.The presence of the ground seems to have no effect when compared to the change in the downwards velocity caused by the blade.
The tip vortex core size is measured based on the same vorticity magnitude level for both IGE and OGE simulations.Fig. 16 shows the growth of the vortex core radius normalized by the chord (c = 0.1905 m).The tip vortex experiences a rapid growth of radius up to 30 deg.Then, it remains about the same size until the blade passage, where it suddenly contracts due to the BVI as shown in the detail.Overall, in terms of vortex core size there is no direct correlation with the ground distance, even though it experiences a considerable change in shape due induced flow, especially after 240 degrees age at Z /R = 0.5.

S-76
Among the experimental available, the S-76 rotorcraft is one of the main benchmark case for near ground and has been used as validation case for many numerical studies [17].The experimental study was conducted in a 1/4.71scaled hovering rotor.The 4-bladed rotor is based on the aerofoil profiles, with a linear twist of 10 deg along the blade's span.The rotor has the radius of 56.1 inches and 3.1 chord length.The rotor performance was measured for several collective pitch angle operating at Mach tip number of 0.55, 0.6 and 0.65 using five tip shapes.The ground effect was considered at three different height-to-radius (z/R) ground positions [51].This vast experimental dataset consists of many rotor scenarios combinations with isolated main rotor, fuselage and tail rotor.
In the present study, the simulations were performed for isolated main rotor using the swept-tapered tip blade at tip Mach number of 0.6 and collective pitch of 9. Two flight configurations are assessed here, out-ground and the near-ground (z/R = 0.75).
For this rotor, the OGE simulation was mainly used for a preliminary comparison for the IGE computations.Then, the IGE simulations, were extended for collective pitch angles of 7 and 11  degrees.In order to reduce the number of elements the flow is assumed as axisymmetric and then the periodic condition was used.For all simulations, the computational domain has a quarter of a cylindrical shape with radius of 15R and length of 17R.To apply periodicity that the symmetry plane, a cylinder shape hub extending from the whole axial direction with 5% radius rotor R. The ground was treated as free-slip wall to avoid an excessive number of elements at this location.The mesh refinement strategy is conducted in similar manner as previous case, the meshes are generated based on the number of points on the aerofoil section, blade  span and regions of interest.The mesh details of the meshes used on the OGE and IGE cases are tabulated in Table 5.For both cases 50 prismatic layers are placed at boundary layer, with the first element height of 1.0 × 10 −6 m.In terms of refinement, the main difference between those meshes is related at the wake.On the OGE mesh, the volumetric local refinement is applied only at the tip blade (r/R = 1) near the wake covering up to 0.7 rotor radii downstream the rotor.While on the IGE mesh refinement cover a  cone-shaped volumetric region between the rotor and the ground as shows Fig. 19.
The hover solutions are highly unsteady with challenging convergence.Having near-zero velocity makes this equation system numerically stiff, especially using more than one reference frame with additional strong source term on a selected part of the domain.Steady-state solutions were achieved using implicit Jabobi iteration algorithm using a CFL number fixed to 10.The convergence of the integrated load on the blade for the 9 degrees with IGE solutions, shown in Fig. 17, has seemingly converged within couple thousands iterations.The simulation running strategy is tailored for convergence acceleration; the second-order solution is used to initialise third-order simulations and subsequently the later is used for the fourth-order It can be seen that the computation in all orders has run long enough to the point that the relative difference between the previous time-step is unnoticeable.
In Fig. 18 the pressure coefficient along the chord is compared between IGE and OGE flow conditions for three blade sections.The difference of pressure distribution between IGE and OGE predictions is relatively small, except the pressure peaks at the leading edge, especially at low blade span station (r/R = 0.2).Other works using seventh-order WENO simulations under the same flow conditions present similar values [17].
To investigate the accuracy of the Multiple Reference Frame in terms of integrated load, the computed results for the pitch angle of 9 degrees are selected as tabulated in Table 6.The comparison is made between the rotor not only with and without the ground effect but also concerning the numerical spatial order.The comparison is made using three main parameters: trust coefficient (C t ), torque coefficient (C q ), and Figure of Merit (FoM).The results re also compared to the experiment [51] in order to measure the improvement in accuracy within the computational cost.For the OGE predictions, it can be seen that both C t and C q are underpredicted.A small improvement on the torque coefficient as we increase the order.As the FoM is an relation between C t and C q (F oM = C 3/2 t /( √ 2C q )), this improvement in accuracy with the or-  der can be misread if we look at the FoM solely.For the IGE results, the improvement with a more diffusive numerical scheme is more evident.While the thrust coefficient computed using low-order scheme under-predict the experimental value, the highest order over-predict it just a bit.The increase in the computed thrust due to the in ground effect is 14.6% while the experiment report approximately 10.4%.Overall, the predicted figure of merit compares well with the experiment for both cases.In Fig. 20 the predicted blade loading (C t ) and torque coefficient (C q ) for three collective blade pitch angle for IGE simulations are compared with the experimental data.Both coefficients are normalized with the rotor solidity (σ ).It is interesting that the computed integrated thrust is slightly over predicted high collective pitch angle and overshoot for the lowest degree.Overall, for all three collective pitch angles, the computed aerodynamic performance parameters C t and C Q agrees well with experimental data.

Conclusion
In the present paper, a MRF algorithm was applied in a highorder unstructured mixed mesh flow solver to simulate two helicopter rotors in hover with and without the ground effect.The domain decomposition into rotating and stationary frames was defined at the solver level, avoiding the use of an interface mesh surface defined at CAD.The numerical solutions were compared with available experiment data, such as vortex trajectory, thrust, torque and pressure coefficients.The flow behaviour of the vortex system for both IGE and OGE conditions were also analysed.The resolution of the vortex path and wake breakdown were considerably improved with increased scheme order.It is shown that the MRF algorithm coupled with a high-order spatial accuracy enhances the prediction of the aerodynamic performance of helicopter blades in hover with and without the ground effect.

Fig. 2 .
Fig. 2. Detail of VMRF interface.(For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Vortex age and trajectory of computed solutions in terms of vortex radial contraction (a) and downwash rate (b) against azimuth angle compared with experiment.

Fig. 8 .
Fig. 8. Iso-vortex surfaces coloured by vorticity magnitude contours comparing the effect of order resolution of vortex wake.

Fig. 11 .
Fig. 11.Mesh cut planes of IGE mesh for the Caradonna and Tung rotor.

Fig. 12 .
Fig. 12. Iso-vorticity surface at different ground distance on Caradonna rotor comparing the effect of the spatial accuracy order on the wake vortex.

Fig. 14 .
Fig. 14.Contour of vorticity (a) and velocity (b) magnitude detailed at the interface for Caradonna and Tung rotor at Z /r = 0.5.

Fig. 15 .
Fig. 15.Vortex age and trajectory of computed solutions in terms of vortex radial contraction (a), downwash rate (b) and against azimuth angle compared with experiment.

Fig. 16 .
Fig. 16.Size of the vortex core versus wake age (in degrees) for the Caradonna and Tung rotor.

Fig. 17 .
Fig. 17.Thrust force convergence for the IGE solution at pitch angle equal to 9 degrees.

Fig. 18 .
Fig. 18.Comparison of pressure coefficient distributions between IGE and OGE at collective pitch angle of 9 deg at different radial sections.

Fig. 19 .
Fig. 19.Schematic of the computational mesh and boundary conditions for IGE and OGE simulations.

Fig. 20 .
Fig. 20.Comparison of predicted and experimental aerodynamic performance parameters for S-76 rotor with ground effect against collective pitch angle: (a) blade loading coefficient C t /σ (where σ is rotor solidity), (b) torque coefficient C Q /σ .

Table 1
[68] specification for the rotational size influence simulations of the Caradonna and Tung[68]rotor flow.

Table 3
[68]ification of the meshes, employed for the simulations of the Caradonna and Tung[68]rotor in ground effect flow.

Table 4
Specification of the meshes, employed for the simulations of the Caradonna and Tung rotor in ground effect flow.

Table 5
Specification of the meshes employed for the simulations of S-76 rotor.

Table 6
Predicted aerodynamic performance parameters for S76 rotor at collective pitch 9 degrees.