Hovering rotor solutions by high-order methods on unstructured grids

Article history: Received 5 August 2019 Received in revised form 2 December 2019 Accepted 16 December 2019 Available online 19 December 2019


Introduction
In the not-so-distant past, academic codes were the backbone of high-order method development, these codes generally lacked the robustness and efficiency of low-order methods (≤ 2). Recently, commercial software developers have introduced third-order spatial discretisation schemes within their numerical arsenal; this is an affirmation that there is a commercial demand of higher order schemes as the maturity of the technology has reach wide adoption level. The scope of these methods is to implicitly increase the resolution of a CFD solution with a high-order approximation rather than refining the mesh; a labour intensive and generally difficult-to-automate task. These methods are applied successfully for CFD applications requiring low dissipation levels and enhanced discontinuities capturing abilities e.g. acoustics applications, turbulence modelling, transonic/supersonic flows and rotating wings amongst others. Helicopter aerodynamics encompass several challenges that high-order methods can address such as shock capturing and vortex dominated flow resolution, where the inherent low numerical dissipation properties of these methods complements the computed solutions.
The earliest methods for modelling rotorcraft aerodynamics were based on the Prandtl's lifting line theory, only in the 1980's advances of Euler based methods reached the point where simulating of an isolated rotor was feasible. Agarwal and Deese [1] investigated the flow around the Caradonna and Tung rotor [2] solv-* Corresponding author. ing the Euler's equations formulated in a rotating reference frame using the Jameson's finite volume scheme with explicit Runge-Kutta time stepping [3]. The work of Caradonna and Tung in 1981 [2], encompassed experimental measurements of a hovering rotor through a wide envelop of flow conditions, their findings provided the CFD community a tangible database for code validation. With an external wake model to compute the inflow velocities and wake effects Wake and Sankar [4] presented a finite difference Navier-Stokes approach using a hybrid implementation of the alternating direction implicit time marching scheme to reach a steady state solution of the ONERA hovering rotor [5].
Computational advances during the following decade, enabled developments of new Euler and Navier-Stokes methods in capturing the rotor wake without external models e.g. vortex confinement techniques [6]. Srinivasan et al. [7] simulated a hovering rotor by a Reynolds Averaged Navier Stokes (RANS) equations model using an implicit, upwind, finite difference scheme on a C-H grid topology; the authors highlighted the importance of well meshed domains for the farfield wake capture. At that time, the current state of CFD methods were insufficient in conserving the wake, since the numerical schemes needed to be dissipative for stability, producing a faster rate of diffusion of the sheets and vortices [8]. Thereby, great research efforts are made on the development of high-order methods and mesh adaptation techniques.
The initial step of a high-order framework involves the formulation of the governing equations; either in integral or differential form or a hybrid combination of theses, to enable the discretisation in time and space by a high-order approximation scheme. There are various flavours of formulating the Navier-Stokes equations. One of the most promising ones is the Discontinuous Galerkin (DG) approach [9] that inherents characteristics from the finite volume [10], and finite element method. A recent overview of high-order methods is presented by Z.J. Wang et al. [11], the paper summarises the collaborative work of several successful International Workshop on High-Order CFD Methods. Another popular formulation is the Flux Reconstruction (FR) technique [12], an extensive comparison of PyFR, an open-source solver based on a high-order FR implementation, with a commercial solver, is presented by Vermeira et al. [13] demonstrating promising efficiencies of high-order systems.
For fluid dynamics problems the finite volume method provides several benefits; flexibility, the ease of formulating transport equation systems, conservation properties on arbitrary-defined spatial elements, extensibility to multi-physics systems and robust numerical frameworks including schemes, solvers and limiters for hyperbolic systems. An example of these is the Weighted Essentially Non-Oscillatory scheme (WENO) and Monotone Upstream-Centred Scheme for Conservation Laws (MUSCL). On unstructured meshes the MUSCL scheme was first introduced by Barth and Jepersen [14] presenting a limiter which acts on the linear reconstruction and it is free of oscillations at strong gradients. Due to its robustness, this method gained popularity over the years and many authors contributed for high-order flux limiters on unstructured 3D grids [15][16][17][18]. WENO schemes [19] are an extension of the Essentially Non-Oscillatory (ENO) schemes first proposed by Harten [20]. Both approaches use an adaptive local stencil to enhance the resolution in smooth region preserving the non-oscillatory behaviour near the discontinuity, but while ENO select the smoothest stencil, WENO employs a weighted combination of all the candidates resulting in smoother solutions with better convergence property. Comparing to the classic ENO, such a improvement increase the stability and accuracy on unstructured grids [10]. Zhang and Shu [21] extended this approach for third-order on tetrahedral grids. Dumbser [22] presented a hybrid method up to sixth-order accurate for compressible Navier Stokes equations. Tsoutsanis et al. [23,24] improved the applicability of WENO schemes extending it for hybrid unstructured meshes. This scheme provides an accurate numerical framework on turbulent shock-wave boundary layer interaction on Implicit Large Eddy Simulations (ILES) [25] and for the RANS equations [26] and has been recently optimised for massively parallel computations [27].
Hariharan et al. [28] used a seventh-order overset WENO scheme on arbitrary block structured moving grids to capture the blade tip vortices of a hovering rotor, and their work outlines the clear benefits of the WENO scheme on capturing the wake turbulent structure and the challenges of overlap grids. Postdam et al. [29] performed steady state simulation of an isolated, half-span, and full-span V-22 tiltrotor hover configurations using the OVER-FLOW code and overset method. The numerical dissipation of the schemes demonstrates a direct impact on the prediction capabilities of blade-vortex interaction. Zhong et al. [30] 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. Costes et al. [31] study the extensive overview of the ONERA experience using CFD techniques for rotorcraft applications. They also emphasise the importance of developing high-order schemes on hybrid unstructured meshes. Hariharan et al. [32] review consistently the historical progress of the hover simulations including benchmark cases, grid aspects, turbulence models and results expected. The GOAHEAD project conducted one of the most prominent research over the recent years in terms of state-of-the-art on experimental and numerical studies. It conducted in a wind tunnel an extensive measurements data of a fully instrumented complete 4-bladed helicopter model, made of the NH-90 fuselage with the 7AD rotor. This database was used to compare with CFD simulations [33][34][35]. Recently Li et al. [36] presented an integral Arbitrary Lagrangian Eulerian (ALE) formulation of the unsteady RANS equations on the rotating reference frame using Roe's difference finite splitting scheme and the Spalart-Allmaras turbulence model on structured grids. The solver utilised a third-order MUSCL and fifth-order WENO schemes, unsteady computations of the blade cyclic operation were performed with reasonable accuracy by the ALE on non inertial formulation. Jimenez-Garcia and Barakos [37] describe the formulation fourth-order MUSCL scheme with correction terms for rotating flows. The proposed method was validated against two-and threedimensional test cases, the UH-60A [38], demonstrating reasonable CPU and memory overhead. Hwang et al. [39] 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.
This paper introduces the work related to the implementation of high-order MUSCL and WENO spatial discretisation schemes, up to fourth-order of accuracy, on a rotating reference frame, based on the k-exact finite volume formulation of the compressible Reynolds Averaged Navier-Stokes equations. Turbulence is accounted through the Spalart-Allmaras turbulence model; solutions are obtained on hybrid (mixed-element) three dimensional unstructured grids. The governing equations are formulated in a non inertial reference frame, the Harten-Lax-van-Leer-Contact (HLLC) Riemann solver [40], the spatial discretisation employs the MUSCL-MOGE variant limiter of Tsoutsanis [18] and the WENO implementation of Tsoutsanis et al. [23,24], 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) [41] has been the focus of research work by the authors for the last decade. The numerical framework is validated, assessed and evaluated across a wide range of flow conditions: inviscid flows [23], time dependent laminar, transitional and turbulent flows [42] and high-Reynolds number flows for external aerodynamics [43]. We study the implementation of the schemes on three hovering rotor configurations, the sensitivity of the mesh on the solution is consistently analysed, the accuracy and the computational cost. The rotor configuration includes the Caradonna and Tung [2], PSP [44] and UH-60A [38]. To the author's knowledge, this is first time that high-order solutions are obtained for rotorcraft aerodynamics on mixed-element unstructured grids. The paper is structured as follows: The governing equations on the rotating reference frame are presented in section 2. An overview of the numerical framework is detailed in section 3. Numerical solutions are presented in the following section, where we discussed in depth the findings. The conclusions of the present work are outlined in the last section.

Governing equations in rotating reference frame
The compressible Navier-Stokes equations are formulated in a blade-fixed reference frame with steady angular velocity. The equations are cast and solved in terms of absolute variables, thereby the flow in the farfield is uniform but the relative flow is not. This approach allows a precise computation of the fluxes on non uniform grids [1]. 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 compressible Navier-Stokes equations with the one-transportequation turbulence model of Spalart-Allmaras [45] for three dimensional flow in a coordinate system attached to the blades can be written in conservative form 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 SA turbulence model [46]: 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 law relating the dynamic viscosity, (ν l = μ l /ρ) to the ideal gas temperature: where S is the Sutherland temperature (S = 110.4 K) and the ref- where ν t is the eddy viscosity, δ ij is Kronecker delta and the subscripts i, j, k are the Cartesian components of x = x, y, z. The work performed by the viscous stress and heat conduction, i , can be written: where, Pr is the Prandtl number and Pr t is the turbulent Prandlt number. The kinematic turbulent viscosity, ν t is approximated by the SA turbulence model, the implementation is comprehensively described in [26].

Numerical framework
The computational domain is discretised 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 Eq. (3) over the unstructured mesh using the finite volume formulation, yields (8) 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 depends 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 process is based on the k-exact reconstruction and aims to build a high-order polynomial p i (x, y, z) of a specified order for each element V i which present the same cell average quantity U i yielding, The reconstruction algorithm is based on the approach given by Tsoutsanis et al. [25] 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 reconstructed polynomial at the transformed element, V i , is expanded over local polynomial basis functions labelled as λ k (ξ, η, ζ ), given by: (10) where ξ, η, ζ are the coordinates in the transformed reference frame, a k are the degrees of freedom. K relates to the order of the polynomial r by K = 1 6 (r + 1)(r + 2)(r + 3) − 1. The decomposition into tetrahedrals is only performed to improve the condition number of the reconstruction matrices, rather than computing the inter-element quantities such as fluxes, due to the advantage in using prismatic or hexahedral elements for the boundary layer for computing the gradients with high-accuracy [47]. To compute a k degrees of freedom at least k elements are required in the stencil.
However, using the minimum (M ≡ k) weakens the solution of the linear system and its matrix may become ill-conditioned. Thereby, M = 2K makes the method more robust and in the present study the stencil based compact algorithm is used for the central stencils and the type 3 for the directional stencils as introduced in the recent work of Tsoutsanis [48]. By using a least square approach, the system of equation (10) allows the solution for the unknown degrees of freedom a k . The QR decomposition solves the final form of the linear system. Limiting functions are used to deal with discontinuities solutions, such as shock waves, maintaining the numerical stability by suppressing non physical oscillations.

MUSCL
The MUSCL method has it origins with the pioneering work of Van Leer [49][50][51], Barth and Jepersen [14] 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. The scheme can be written as: (11) where U ij,α is the extrapolated reconstructed solution at face j, and at quadrature point α, U i is the value for the conserved variable of element i, and (ξ a , η a , ζ a ) are the coordinates of the quadrature point at the ij face. All polynomials, for all the faces and for all the quadrature points are then limited by the limiter φ i which is valid for cell i 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: where l = 1, .., L; L is the local numbering for the cell i and it's direct-side neighbours, where l = 1 corresponds to the cell i.
However in 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 [18]. Therefore the minimum and maximum values from the all the elements in the stencil including the considered cell i are considered: where m = 1, .., M; M is the local numbering for the cell i and it's stencil elements, where m = 1 corresponds to the cell i. The limiter seeks the minimum value of the slope limiter for all the quadrature points α that satisfy the following conditions: Where φ il,α corresponds to the slope limiter value at side l and quadrature point α for cell i. Then, the limiting function is applied, composed by three different states according to the difference of the unlimited reconstructed value U il,α at the quadrature points of the considered element U (i,l,α) , the minimum U min i and maximum U max i values from the entire stencil, yielding: However in the case of the MUSCL-MOGE variant the limiting function used is a new function min(1, y) which replaces the minimum function of Barth and Jepersen [14] to achieve higher order of accuracy in smooth regions of the flow, and the reader is referred to [18] for further details.

WENO
In WENO scheme the non-linearity is preserved by combining a number of reconstructed polynomials arising from various stencils, which are weighted according to the smoothness of its solution.
The key aspect is that the coefficient of the reconstruction polynomial are dependent only on the mesh and on the particular stencil considered, but not on the solution. In general, the polynomials are defined as where m s is the total number of WENO stencils. Applying the scheme on Eq. (10) it becomes p m (ξ, η, ζ ) = K k=1 a m k λ k (ξ, η, ζ ). (17) As the sum of all weights is unity, we obtain where ã k are the reconstructed degrees of freedom and the nonlinear weights, ω m is defined as where d m is the linear weight, I m smoothness indicator and is a small number to prevent division by zero. The typical values for is 10 −6 and b = 4, and for the linear weights a large value is assigned to the central stencil d 1 = 10000 and a value of d 2,3,...,s = 1 to the remaining directional stencils. The main idea is to define for each stencil that generate a polynomial approximation for the cell a smoothness indicator, which yields where β is a multi-index, r is the polynomial's order and D is the derivative operator. This smoothness function is a quadratic function of the degrees of freedom and can be expressed as a universal mesh-independent oscillation indicator matrix. Further details concerning the WENO characteristic reconstruction can be found in [23,25]. In addition, the pressure and density remain positive for each Gaussian quadrature points (α) through the reconstruction process by the positivity condition of Harten et al. [52]:

Fluxes
The Riemann problem is solved with the approximate Harten-Lax-van-Leer-Contact (HLLC) solver [40] 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: W * ± is calculated for the element itself or its neighbour "+". The approximation of the wave speeds, S ± and S * , are iteratively calculated as explained in [53]. 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 [25,54]. For the gradients additionally penalty terms are included following the formulation of Gassner et al. [55] for suppressing odd-even decoupling modes in the numerical solutions [56], in the following manner: where L int is the distance between the cell centres of adjacent cells, and α = 4/3 similarly to previous approaches [56,57].

Demonstration of the spatial discretisation methods on 2D Vortex Evolution
The 2D vortex evolution test problem introduced by Balsara and Shu [59] is used to demonstrate the accuracy of the employed numerical framework. In this test problem an isentropic vortex propagates at supersonic Mach number at 45 • across the domain, which is ideal for assessing the multi-dimensional discretisation algorithms used. The computational domain is given 10] with periodic boundary conditions applied on the sides. The unperturbed domain has an initial condition (ρ, u, v, p) = (1, 1, 1, 1), where temperature and density are defined as T = p/ρ, and S = p/r γ the adiabatic gas constant γ = 1.4 and the vortex perturbations are given by: The vortex strength = 5 and adiabatic gas constant γ = 1.4. The e L 2 and the e L ∞ errors are computed as follows: where U c x, t f and U e x, t f are the computed and exact solutions at the end of the simulation t f . The exact solution U e x, t f coinciding with the initial condition at t 0 . One type of hybrid unstructured mesh is employed consisting of quadrilateral and triangular elements for this test problem of 16, 32, 64, and 128 edges per side as shown in Fig. 1, and the simulation is run for a time of t f = 10 with a third-order Strong Stability Preserving Runge-Kutta method [60]. As expected all the schemes achieve convergence rates close to their theoretical ones as seen in Table 1, and the MUSCL3, and MUSCL4 being at least two times and three times faster than the corresponding WENO schemes respectively. The MOGE formulation of the MUSCL limiter version employed in this study exhibits high-order of accuracy at significant reduced computational cost compared to WENO schemes.

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 Eq. (8) 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, Eq. (28) yields Linearising in time, Eq. (29) becomes where R i tends to be equal to zero and ∂R n i ∂U is the flux Jacobian, which contains the linearisation of convective and viscous flux vectors and the source terms. Thereby, Eq. (30) gives where I stands for the identity matrix. Once the linear system is solved, the solution at each element i is simply updated as U n+1 The residual is linearised by a first order approximation of the numerical flux using the Rusanov flux function as 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 Eq. (32) gives the block diagonal and off-diagonal elements as: The linear system of the form A X = B in Eq. (31) 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 Jacobian of the non-inertial source term which is added directly on the diagonal component to keep its dominance [58]. The diagonal elements of the matrix are stored and inverted directly and the off-diagonal ones are calculated at every stage.

Caradonna-Tung rotor
The Caradonna and Tung [2] experiment is conducted in a controlled wind-tunnel facility at NASA Ames Research Center, where a two-bladed rotor is used. The blade is based on an untwisted and untapped NACA 0012 aerofoil with a radius of 1.143m and an aspect ratio of 6. Two main aspects of the flow around the blade were measured: the surface pressure distribution by pressure tubes in five radial sections for both upper and lower sides of the blade; and also the tip vortex strength and trajectory captured by hot-wire probes mounted underneath the rotor. A wide range of blade tip Mach numbers (M tip ) and collective pitch angles (θ c ) were assessed. The developed schemes and the nature of the solver is tailored for compressible flows, the case with a blade tip Mach number of M tip = 0.89 and pitch angles of θ c = 8 • is considered for the computational assessment. The corresponding operating conditions are considered fairly challenging both in terms of the experiment but also with respect to prediction abilities of CFD solvers.
The computational domain has a cylindrical shape with 40m radius and 80m length, in which the helicopter blade is placed at the center. An angular velocity of 265.9 rad/s is applied to the whole domain, and at the farfield surface, characteristic inlet/outlet Riemann boundary conditions are applied, in which the choice of inlet/outlet conditions is made on the extrapolated value of the contravariant velocity at the interface of the domain boundaries. No-slip adiabatic wall boundary conditions are applied to the blade.
Five mixed-element unstructured grids are generated to perform computations and study the discretisation impact combined with two numerical approaches MUSCL and WENO of various reconstruction orders. The specifications of the grids are shown in Table 2.
The blade surfaces are meshed with quadrilateral elements, as they capture efficiently the surface by preserving its curvature in the vicinity of the leading edge. The meshed surfaces of the blade are projected normal to their base and the inflation layers are grown. Hexahedra-type elements are generated to capture the boundary layer, they are aligned with the flow and maintain orthogonality in the near-wall region. The first wall cell distance of 10 −5 m is set by using the 1/7 power law suggested in [61] to ensure a y + < 1. The rest of the domain is composed by isotropic hexahedra and tetrahedra elements automatically generated by the advancing front mesh algorithm [62]. Through a consistent refining process on the blade span and upper/lower surfaces, the grid levels range from 0.83 × 10 6 to 6 × 10 6 elements. An additional fifth grid is constructed specifically to assess the numerical accuracy in capturing the tip vortex trajectory. The grids are shown in Fig. 2, where in red the cross sectional mesh is shown near the blade tip, in blue the surface cell distribution on the blades for the four consecutively refined meshes and in black the wake refinement on the fifth grid compared with the coarse grid. For the refined wake grid, the coarse grid is used as baseline and a volumetric local refinement is performed at the tip blade (r/R = 1) near the wake covering up to three rotor radii downstream the rotor.
Computations are performed by employing the two classes of spatial discretisation namely the MUSCL-TVD and the WENO. Highorder approximations are computed of third-and fourth-order for both scheme classes and compared with a standard second-order approach. We evaluate the effect of each scheme and discretisation order on all meshes by computing the integrated blade load of thrust coefficient C t given by where T is the thrust, R is the blade radius and U tip is the tip velocity. The computed thrust coefficients are shown in Fig. 3, where each scheme is shown with a different colour. At first glance, it is reasonable to believe that grid independent solutions could be obtained with all numerical approaches for the two most refined grids, particularly the WENO predictions demonstrating clear asymptotic trend. These predictions on the grids with the lowest number of elements, suggest a C t much closer to the grid independent value compared with their MUSCL counterparts. The fourth-order MUSCL solution exhibit opposite behaviour in terms of grid convergence monotonicity compared with the other computed thrust coefficients, and overall has the highest discrepancy compared with the experiment. This is also mirrored in the coefficient of pressure plots shown in Fig. 4 at five radial cross sections, this scheme is the only one producing spurious oscillation of pressure, particularly in the vicinity of the shock near the blade tip (r/R = 0.8). This erratic behaviour is partly attributed to the pronounced difference between the reconstruction stencil and the elements used for bounding the solution for a fourth-order MUSCL scheme, as it has been reported in [18], and unless the bounds of the entire reconstruction stencil are considered these types of oscillations can occur. Furthermore, from the C p sectional plots of Fig. 4 the flow behaviour can be depicted, the flow mid-way the blade is smooth, and the computed profiles are in good agreement with the experiment. At r/R = 0.5 it can be observed that higher than secondorder predictions capture the suction peak more precisely, particularly the WENO ones. The same trend is observed at the following section, r/R = 0.68, where the rotational velocity is not as high and the flow is still smooth. It is noteworthy that MUSCL4 present small spurious oscillations at 0.35 < x/c < 0.6. At r/R = 0.8 the flow is characterised by a shock at x/c ≈ 0.25 and the solutions fairly agree with the experiment. It is worth mentioning that even on the coarsest grid the schemes are able to predict the shock position with reasonable level of accuracy, particularly with the weighted reconstruction. The MUSCL predicts a more smeared shock, especially the fourth-order which slightly oscillates at x/c ≈ 0.4 due to the bounds difference as reported in [18]. The influence to the tip vortex is clear on the last two sections, r/R = 0.89 − 0.96, in which all schemes struggle to capture the pressure drop and its recovery. Overall, the computed solution fairly describes the pressure curve over the span sections. Nonetheless, small inaccuracy is expected in numerical solutions and may be due to deficiencies related to the turbulence model employed.
The wake and its particular strong vortical structures which develop at the tip of the blade are the most complex feature that characterises the flow over hovering rotors. For every rotation, the blade interacts with the vortex generated by its previous passage. This phenomena, known as Blade-Vortex Interaction (BVI), is sensitive to the numerical approach and spatial resolution [31]. In this sense, an accurate computation of the wake trajectory and strength is much desired and challenging in rotor flow analysis. Thereby, high-order schemes can greatly improve the accuracy of wake vortex estimations. The wake-refined grid, shown in Figure 3, is employed to assess the vortex capturing capabilities of the developed numerical schemes. Fig. 5 shows the predicted tip vortex trajectories of the computed solutions compared with the experiment hot-wires measurements [2]. The plots present the axial and The computed solutions slightly overpredict the vortex contraction, this is also observed in [36]. 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, but the effect on the numerical scheme on the radial contraction is negligible. In Fig. 5 (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 the advancing blade at 180 • . At higher degrees, the downwash rate increases due to the passage of the retreating blade and it is clear the advantage of the low dissipation nature of the high-order schemes on predicting this phenomenon. Looking at the slope of second-and third-order predictions, it is evident that the numerical dissipation can be major obstacle for accurately predicting the vortex's path. It is worth to draw attention to the precision of the fourth-order approaches on both reconstruction schemes which agree with experimental measurements. In this sense, the lower dissipation of highorder scheme helps preserve the vortex core strength as shown in terms of vorticity magnitude shown in Fig. 5 (c), whilst the vorticity levels are significantly lower for the second-order solutions at θ = 200. Fourth-order solutions are able to predict BVI, shown with the small bump just after first blade passage at ≈ θ = 180 • . Overall, it can be seen that both reconstruction schemes with respect to their approximation order behave similarly in terms of vortex contraction, convection and strength. Fig. 6 shows the vorticity magnitude contours at the cutting plane along the rotor center and blade. Looking at the vortex struc-ture convecting downstream the rotor it is noticeable the advantages of the numerical scheme as we increase its order: the tip vortex is captured for two full rotations on both third-order solutions and up to three from the fourth-order. The vortex core it is clear distinguishable on high-order scheme while it is more dissipated on second-order. It is also noted that fourth-order contours are not completely symmetrical as it is more efficient in computing the unsteadiness and instability of the flow. Fig. 7 shows the helical wake structure of the tip vortex using iso-surfaces for Q − criterion = 1000. It is evident that high-order method can improve the prediction of the wake sheet and vortex tip flow phenomena, better preserving the helical vortex filaments that trails from each blade tip, up to several rotations. It is worth noting that the wake resolution is affected by the mesh shape and refinement [32] and capturing a clear vortex core trajectory is challenging for many codes, particularly when unstructured grids/solvers are employed.  5. Vortex age and trajectory of computed solutions in terms of vortex radial contraction (a) and downwash rate (b) against azimuth angle compared with experiment [2] and published CFD data [36] and vorticity magnitude (c). Numerical solutions are computed on the wake-refined grid.

PSP rotor
The flow over the PSP (Pressure Sensitive Paint) rotor is simulated, Watkins et al. [63] carried out extensive experimental activities on this configuration. Compared with the Caradonna and Tung rotor case, the PSP entails a much more complex blade geometry and topology, a property that is in tune with the flexibility of the unstructured grid framework. The four-bladed rotor is made of three aerofoil profiles with a linear twist angle along its blade span: RC(4)-12 aerofoil at r/R < 0.65, RC(4)-10 aerofoil at 0.7 < r/R < 0.8 and RC(6)-08 aerofoil at 0.85 < r/R < 1. The last section having a 60% tapered and 30 • sweep angle. The main attributes of the blade are summarised in the Table 3. Further details concerning the aerofoil coordinates and blade radial twist and chord distribution are presented by Noonan [64] and Watkins et al. [63], respectively.
The spatial domain is composed of by three boundary conditions: viscous wall, rotational periodicity and farfield, depicted in Fig. 8. The non-slip wall condition is applied on the blade's surface, rotational periodic conditions permit efficient computational resourcing by modelling the rotor in an axisymmetric fashion. A grid refinement strategy is conducted, the grids are generated based on the number of points on the aerofoil section, blade span and regions of interest, details on the grids are tabulated in Table 4.
The number of elements on the blade span and aerofoil are chosen carefully to maintain good geometric quality metrics such as skewness, aspect ratio and growth rate. Fig. 8 depicts the blade surface mesh, consisted of quad-dominant mesh elements which contributes for the generation of hexahedrals on the boundary layer. To get a full resolution of the boundary layer, 60 nodes are positioned on the surface and orthogonally projected with a 1.2 grown rate, the height of the first element is positioned at 1.6 × 10 −6 to insure y + < 1. Three meshes are generated ranging from 2.64 to 15.2 millions elements, mainly composed by hexahedra and tetrahedra.
Steady state solutions are achieved with the implicit LU-SGS algorithm, which is extensively employed for high-Reynolds flows on unstructured grids as it is free of matrix inversion and has moderate memory requirements at reasonable convergence rate [65]. A Courant-Friedricks-Lewy (CFL) number is fixed to 20, which is estimated based on the smallest edge in the spatial domain. The simulation running strategy is tailored for convergence acceleration; the second-order solution is used to initialise third-order simulations, for both scheme classes, and subsequently the later is  used for the fourth-order ones, with its respective scheme class. Convergence acceleration procedures and techniques can also be extended to the spatial grid, with multigrid methods for example; however, the sequential procedure is strictly applied to the same mesh, to minimise the sources of computational uncertainty. The convergence decay of each scheme on the medium mesh is shown in Fig. 9; the log l 2 norms for the mass conservation equation are plotted against the accumulated iteration number. Both scheme types have similar convergence speed and level, with the MUSCL one slightly better, the restart on the higher order scheme results in a spike, followed with a sudden drop. The WENO achieves lower levels with a smoother trend, suggesting lower convergence levels in fewer iterations. Fourth-order MUSCL solution exhibits fluctuations similar to the previous case, this could be attributed to the combination of low dissipation and inefficiency of slope-limiters to dismiss spurious oscillations. A possible remedy to circumvent this behaviour of the MUSCL approach is with an artificial dissipation technique as recently presented by Jimenez-Garcia and Barakos [66] on rotating flows.
The pressure coefficient C p is compared with experimental data at two tip radial stations (r/R = 0.93 − 0.99) at M tip = 0.585 on the medium mesh illustrated in Fig. 10. Wind tunnel measurements consists of C p by two techniques: kulite transducers and pressure sensitive paint. CFD predicts the overall C p curve at a reasonable level of agreement compared with both experimental techniques. An undershoot of pressure is observed by simulations in the vicinity of the leading edge stagnation point, this trend is confirmed also by Jimenez and Barakos [67] high-order solutions by the HMB solver. The agreement is improved moving towards the pressure recovery zone; it is worth noting that the predicted pressure on the blade, appear to be independent of the scheme type and order and mesh. This is somewhat expected for this case, as the flow is subsonic and compared with the previous test-case the flow is relatively smoother without strong discontinuities. A property of rotorcrafts that assessed the aerodynamic performance and its efficiency is the Figure of Merit (FoM) yielding (38) where C Q is the power coefficient and C t the thrust coefficient; FoM evaluates the amount of generated thrust for given power as the ratio of the ideal power required to hover over the actual power required. Overmeyer et al. [68] conducted an experiment measuring the FoM of this case, the model employed included fuselage and rotor. The computed FoM against C t /σ is plotted for all grids and schemes in Fig. 11 compared with Overmeyer et al. [68] measured data points.
The θ c = 12 • case is used as the selecting point for grid and numerical schemes order, as it corresponds to the one of the highest thrust tested and thereby the most challenging for CFD methods. Lower thrust points were computed using only the fine grid and second and third-order MUSCL to demonstrate the validity of the calculations across C t /σ range. At high thrust, it is clear the range of dispersion concerning numerical method and space discretisation. Even though the precision requirements are not well defined in the industry, aerodynamicist define a numerical prediction to be robust when the FoM ranges less than ±0.02 from experimental data [71]. Using this recommendation, only the refined meshes achieve accurate results. Even though the agreement improves  with increasing the discretisation order, all solutions obtained on the coarse grid cannot reach the required accuracy in agreement with the experiment. However, the medium (green symbols) and fine (solid red symbols) meshes perform consistently better on both reconstruction schemes and even MUSCL2 ( ) reaches high precision level, particularly on the fine mesh (red ). The weighted fourth-order scheme ( ) demonstrates to have the smallest deviation from the measured data. It is worth to mention that even though the computed thrust and FoM showcase good agreement with experimental data, discrepancies related to fuselage presence and unsteadiness must be taken into account to assess numerical uncertainty. According to Rohit et al. [69], an isolated rotor configuration, such as the present case, predicts lower FoM values than when one is installed with the fuselage, this is mainly attributed to the effect of upwash vortices induced by the fuselage. Thereby, this may be one reason that most of simulations underpredict FoM when compared against wind tunnel data. As matter of numerical comparison, the present work presents similar results and trends with CFD solutions found in the public domain [67,69,70].
The computational footprint differs for the class of scheme and of course the discretisation order. The attributed cost of each method and the obtained accuracy will determine the overall efficiency of the method. Table 5 shows the computational time in work units per implicit iteration for the simulation of the PSP rotor flow. The work unit is the non-dimensional time required to compute one implicit iteration of the MUSCL2 method at each grid. A similar approach evaluates the cost associated with the num- ber of elements, in this case the MUSCL2 results using the coarsest mesh is taken as reference time and expressed as CPU time. The relative error associated to each simulation is given in percentage of deviation from the FoM experimental measurements (run 156) [72] under the same operating conditions. For a detailed overview of the parallelisation performance of the UCNS3D code, the reader can refer to [27]. All solutions are performed on the Delta, HPC facility, at Cranfield University. The system consists of two intel E5-2620 v4 processors with 16 cores per node and 128GB of shared memory.
The cost of WENO scheme is noticeable, as the whole process of computing the reconstruction polynomials arising from several stencils, computing the oscillations indicators and applying a nonlinear combination of them is significantly more expensive than a scheme (MUSCL) that uses one reconstruction polynomial per element. In general, at same grid level its cost is around three and seven times higher than third-and fourth-order, respectively. This price pays out in terms of accuracy, as WENO results are more accurate. On other hand, MUSCL presents a much cheaper procedure costing up 3 times as the lowest order on the same mesh and delivering a reasonable error level.
It can be seen that the computational cost dramatically escalates with the number of elements. Overall, this scalability factor is about the same ratio of elements between the meshes. In this sense, we can compare solutions with the similar precision using different grids. Looking for margin of error lower than 3%, WENO3 on the medium mesh presents the most efficient one since it consumes half of the computational time of the fine mesh using the same scheme. Increasing this margin up to 5%, with a small difference in accuracy MUSCL3 on the medium mesh presents the most

UH-60A rotor
The UH-60A "Blackhawk" rotorcraft is part of an extensive program carried out by NASA which includes a large experimental database measured in-flight level conditions. Quantities acquired include aerodynamic pressures, structural loads and rotor forces and moments [38]. The UH-60A database is extensively used for code validation for fluid [32,36] and structural solvers [73][74][75]. The configuration entails a 4-bladed rotor based on SC1095 and SC1094R8 aerofoil profiles, with a non linear twist along the blade's span [76]. The blade has an aspect ratio of 15.3 and a sweep angle of 20 • from r/R = 0.93 to its blade tip. The rotor operates at M tip = 0.628 with a collective pitch angle of 6 • , 8 • and The computational domain has a quarter of a cylindrical shape with radius of 20m and length of 120m. The rotor is placed 40m from the top surface with the same boundary conditions set-up with the PSP case. The grid is also refined in the same manner by increasing the number of points in both aerofoils and blade span. In this case, 50 prismatic layers were placed at boundary layer, with the first element height of 1.9 × 10 −6 m. Three meshes are generated ranging from 4.89 to 15.69 millions elements shown in Table 6. (See Fig. 12.) The simulations are performed by following the same strategy as with the PSP rotor flow configuration. The converged secondorder solution kick-starts the third-and fourth-order solutions. Fig. 13 illustrates the log of l 2 norm of momentum as function of the accumulated iterations for UH-60A rotor on the coarse mesh. Both scheme types on third-order approximation manage lower residual levels compared with the second-order, that plateaus  higher at ≈ 10 −4 . The fourth-order one, inheriting higher dispersion characteristics, asymptotically convergence earlier; it's worth noting that the implicit solver with the WENO spatial scheme appears to have more desirable convergence properties than its MUSCL counterpart, in terms of convergence's speed and threshold level.
The pressure coefficient on three radial stations (r/R = 0.400, 0.775 and 0.945) is extracted from the three-dimensional computed solutions, Fig. 14 depicts the C p profiles compared with the experiment [38]. In general, there is good agreement between CFD and measurements and as with the PSP case no significant variation between the computed pressure profiles and experimental data. We further evaluate the efficacy of the numerical schemes by analysing the computed thrust and power coefficient in terms of FoM. The coefficients are represented by the FoM and compared to four flight test program data summarised by Shinoda et al. [77]. The full scale test was conducted by Nagata et al. [78] at the NASA Ames 80x120 feet wind tunnel. The second data set refers to Airloads Program [79] that was conducted by US Army/NASA and included both shaft and total engine power measurements. The last two data sets used a model-scale at two test facilities: the Sikorksy Model Hover Test Facility [38] and Duits-Nederlandse Windtunnel (DNW) in Netherlands [80]. There is experimental scatter mainly due to the differences in operating conditions, local temperature magnitude, transducer locating tolerance and other factors related to ambient conditions. These and other issues, such as scalability and facility influences still not fully understood yet for both experimental [77,81,82] and numerical works [83].
The computed FoM is shown in Fig. 15 and compared with the aforementioned wind tunnel data points. Both numerical and experimental data demonstrates the challenges in capturing this  ulations employed both MUSCL and WENO schemes on three grid refinements and three discretisation orders, the CFD data is compared with four experimental data: Nagata et al. [78] Airloads Program [79] and by Lorber et al. [38,80] and published CFD data Zhao et al. [84] and Schmitz et al. [71]. complex unsteady flow physics at high thrust levels. Even consistent measurements which took part in highly instrumented wind tunnel tests significantly differ one from another. This uncertainty between the existing data sets were highlighted by Shinoda et al. [77] who pointed out the FoM ranging up to 0.05 under the same wind tunnel facility. As stated by Schmitz et al. [71] this poses a dilemma in the numerical validation process once the experimental discrepancy extrapolates the accuracy recommended by aerodynamicists of FoM varying not more than 0.02 from the dataset.
Taking into account these uncertainties in measurements and the lack of a fuselage on the numerical simulation, the computed performance fairly agrees with those experimental data set. At high thrust, a clear range of accuracy is seen in terms of mesh and reconstruction order. Therefore, precise solution is achieved by either refining the mesh or increasing the scheme order. Taking into comparison the medium and fine meshes, 7M and 15M elements, respectively, overall it seems that increasing the order of the reconstruction method is more effective than increasing the grid resolution, since the higher-order scheme better capture the physical complexity of the hovering rotor. At low thrust point (θ c = 6 • and 8 • ), the computed integrated load using the fine grid and MUSCL2 and WENO3 matches with the experimental data showing the validity of the calculations across the C t /σ range.

Conclusion
We study, analyse and evaluate the computational performance in terms of accuracy and cost of the UCNS3D solver on three hovering rotor flow configurations. Numerical solutions are obtained by the finite volume, k-exact, high-order (up to fourth order) discretisation approach by solving the 3D compressible RANS equations on mixed-element unstructured grids. The high-order numerical framework has proven to considerably enhance the accuracy of the predictions on given grids, with desirable convergence properties, particularly by the WENO approach. Pressure coefficients are compared with corresponding experimental data points where high-order solutions better approximate discontinuities and strong gradients. In terms of integrated load predictions the most accurate scheme, fourth-order WENO, demonstrates the smallest discrepancies compared with experiment but with the highest registered cost. Assuming a level of accuracy of FoM, the trade-off in terms of cost and accuracy demonstrates that the third-order approach would be have the best of both cost and accuracy. Resolution of the vortex path can be considerably improved with increased scheme order, up to three full rotations where the standard second-order scheme can achieve only one, on the same grid. Currently, research is conducted on improving the limiting functions of the MUSCL approach, as well as extending the single rotating frame to multiple including a sliding mesh approach to enable the undertaking of complete helicopter (including fuselage) simulations and comprehensive turbo-machinery flows.

Declaration of competing interest
The authors declare no competing interests.