Comparative Study of Coaxial Main Rotor Aerodynamics in the Hover with the Usage of Two Methods of Computational Fluid Dynamics

The paper is focused on numerical modeling of the aerodynamic characteristics of a full-scale coaxial main rotor in hover. The simulation was performed using two approaches of computational fluid dynamics (CFD): the original free wake model (FWM) developed by the authors and the unsteady Reynolds-averaged Navier–Stokes (URANS) equations method based on the Ansys Fluent software. The structure of the rotor vortex wake, flow images, vorticity and induced velocity fields, total and distributed aerodynamic characteristics of the rotor, including the rotor performance and figure of merit diagrams, have been analyzed. A comparison between FWM/URANS based calculations and calculated/experimental data by other authors has been performed. A satisfactory match of these data confirms reliability of used methods when modeling the aerodynamic characteristics of coaxial rotors. The FMW demonstrates a significant advantage in speed and resources intensity when calculating the total aerodynamic characteristics of a coaxial rotor. The URANS method makes it possible to model with significant accuracy the effects associated with the blade vortex interactions and pressure distribution over the blade surface. Finally, conclusions about most effectiveness of joint application of considered methods in solving complex problems of coaxial rotor aerodynamics has been made.


INTRODUCTION
The coaxial main rotor, consisting of a pair of rotors mounted on one axis and rotating in different directions, has been the feature of the Kamov company helicopters for more than 70 years. According to the coaxial main rotor scheme, almost the whole model line of Kamov helicopters was created, from the Ka-8 to the Ka-52. The absence of the power consumption on the rotation of the tail rotor and a number of other advantages of the coaxial main rotor allowed coaxial helicopters to hold their place in the global helicopter fleet. At the same time, the special features of the coaxial main rotor scheme in horizontal flight made it the basis for the advanced blade concept (ABC) technology first implemented by the Sikorsky Aircraft company in the S-69 high-speed helicopter demonstrator. Since then, the coaxial main rotor has become the basis for many high-speed helicopter projects, which has significantly increased interest in this aerodynamic scheme of the main rotor. Therefore, coaxial rotor aerodynamics has been the subject of a large number of researches and nowadays it does not lose its relevance.
A number of experimental works in wind tunnels (WT) and on special experimental stands were focused on the study of coaxial rotor aerodynamics (Harrington 1951;Nagashima and Nakanishi 1983;Coleman 1997). The experimental studies and flight tests have been conducted by TSAGI and Kamov company starting from the 1950s (Bourtsev et al. 2001;Petrosian 2004).
In recent years, numerical modeling of aerodynamics has taken an increasing place in the process of creating and improving aircraft. This is primarily due to the growing capabilities of computer technology, which allows implementing in practice mathematical models that adequately represent the complex physical processes of flow around various elements of aircraft. The use of numerical modeling in solving problems of rotorcraft aerodynamics makes it possible to supplement significantly and replace partially model experimental tests in WT and flight tests. Everything mentioned above applies to helicopters, since the calculation of the aerodynamic characteristics of the main rotor is one of the most difficult problems of aerodynamics, especially in flight modes, when the structure of the vortex wake behind the rotor is significantly nonlinear and nonstationary.
Currently, there are two main approaches in numerical modeling of helicopter rotor aerodynamics that belong to the class of computational fluid dynamics (CFD) methods. These are approaches based on various free vortex wake models and based on the solution of Navier-Stokes equations by the finite volume method (FVM) using various turbulence models.
Modern vortex models represent a vortex wake behind a rotor in a nonstationary and nonlinear setting. This allows to model accurately the field of induced velocities around the rotor and the interaction of vortex of the rotor wake, which is especially important in the study of coaxial rotors. Many recent researches study coaxial rotor aerodynamics on the basis of various vortex models (Yana and Rand 2012;Singh and Friedmann 2017;Fernandes 2017;Ignatkin et al. 2018;Tan et al. 2018;Kritsky et al. 2020). Also, it is possible to mention the work performed on the basis of the original vorticity transport model (VTM) (Kim and Brown 2010). It should be noted that, in general, vortex models are not too demanding on computing resources and are often able to solve applied problems of helicopter rotor aerodynamics using conventional personal computers.
The main problem of FVM based CFD methods is the simulation of turbulent flows. To date, there are no universal turbulence models that can cover all the turbulence scales. Currently, there are four main directions of solving the problem of modeling turbulent flows. Firstly, there is the direct numerical simulation (DNS), which is a method for solving unsteady Navier-Stokes equations without any turbulence models (Moin and Mahesh 1998). Secondly, it is necessary to mention the method for solving unsteady Reynolds-averaged Navier-Stokes (URANS) equations using algebraic or differential turbulence models (Menter 1993;Wilcox 1994;Spalart and Allmaras 1992). Next is the large eddy simulation (LES) method based on the solution of unsteady Navier-Stokes equations with modeling the effect of sub-grid scale vortex (Ferziger 2000). Finally, there is the method of detached eddy simulation (DES), which is a combination of two approaches mentioned above, where URANS is used in the external flow zone and LES is used in the flow separation zone with large eddies (Strelets 2001).
Finite volume method based CFD methods require powerful supercomputer clusters for solving problems of helicopter rotor aerodynamics. The URANS method is currently used more often, as it is the least resource-intensive. Now, the URANS method is successfully used both for calculating single (Ignatkin and Konstantinov 2012;Garipova et al. 2014;Ridhwan et al. 2017;Abalkin et al. 2020) and coaxial rotors (Deng et al. 2019;Dacheng et al. 2019;Kinzel et al. 2019). Less often, the DES method is used (Dehaeze et al. 2018).
The present work is focused on numerical modeling and comparative analysis of the aerodynamic characteristics of a coaxial main rotor based on two calculation methods: the free wake model (FWM), developed by the authors Ignatkin et al. (2009) and the URANS method, based on the Ansys Fluent software. This approach allows evaluating possibilities and perspectives of using both models to solve problems of a helicopter coaxial main rotor aerodynamics. The results obtained also serve as validation of the FWM.

THE OBJECT OF STUDY AND ASSUMPTIONS
This paper considers the coaxial main rotor of Ka-226 helicopter (Vassiliyev et al. 2007;Bourtsev et al. 2007). The coaxial main rotor structurally consists of two rotors: upper rotor (UR) and lower rotor (LR) with the same radius located one above other on the same axis of rotation and rotating in different directions. The planes of rotation of the rotors are separated by a distance h. The blades of the rotor are rectangular in shape in the plan view and have a blade twist. The main parameters of the rotor (Bourtsev et al. 2007) are shown in Table 1. The aerodynamic characteristics of the rotor in hovering modes are determined for the range of blade pitch angles of the upper and lower rotor blades in the range of blade pitch angles θ = 6−16°. Blade pitch angles of the upper θ UR and lower θ LR rotor blades are selected so that upper and lower rotors are balanced by the values of rotor torque Q. The required blade pitch angles of the upper and lower rotors were previously calculated with the rapidly FWM method and then applied in the URANS method. This significantly reduced the calculation time on the base of high resourceintensive URANS approach.
For both calculation methods, the rotor blades were modeled to be absolutely rigid for bending and twisting. At hover, the blades of the upper and lower rotors rotate with fixed cone angles, depending on the operating mode of the rotor. Used models include this point. In the FWM based calculations, the blade cone angles were determined by the solution of the blade flapping motion equations. The calculation of the blade flapping motion was not implemented in the URANS method. On each of the calculation modes, the blades were initially installed with the cone angles defined on the base of FWM calculations.
The description of the FWM and URANS methods used for numerical study are given in the following sections.

FREE WAKE MODEL
The FWM of the rotor (Ignatkin et al. 2009) developed at the Helicopter Design Department of the Moscow Aviation Institute is based on the lifting line theory and the blade element theory. Each blade element is modeled by an attached vortex segment located on a quarter of the blade element chord c ( Fig. 1) with the control point in the center of the segment. For each time step Δt a quadrangular contour consisting of vortex segments with a constant circulation Γ, equal to the attached vortex circulation, descends from the blade element. The circulation of the attached vortex changes along the blade radius and depends on the blade azimuthal position. The system of vortex contours creates a free vortex wake in the form of a grid of longitudinal and transverse vortex segments (Fig. 1). The vortex wake grid is deformed at each calculated step under the influence of external and induced velocity fields.
where the coefficients of the lift force C L and the drag force C D of the blade element are determined at the found angle of attack α and the total flow velocity W based on airfoil steady test data in wind tunnel. Determining the attached vortex circulation is implemented through an iterative method.
A key part of the considered model is the calculation of induced velocities from the vortex segment, which is an element of the free vortex wake behind the rotor blades. Experiments prove (Leishman 2006) that the distribution of induced velocity close to the vortex axis is different from the analytical expressions of Biot Savart (where induced velocity value on the axis is equal to infinity) or Rankine (with a linear law of velocity inside the vortex core). Formulas that describe more accurately the velocity distribution inside the vortex core are also proposed by Scully and Vatistas (Leishman 2006). The described model uses an approach that allows to determine the vorticity field from the vortex segment and the induced velocity from this vorticity field on the basis of an exact solution of the vortex source diffusion (Ignatkin et al. 2009).
Vorticity of the vortex segment of length L, taking into account its diffusion, is determined by Eq. 2 (Ignatkin et al. 2009): where Γ is the circulation, ν is the coefficient of kinematic viscosity, r 1 is the distance from the center of the vortex segment with length L, t is the diffusion time (Fig. 1). Note that by integrating Eq. 2 for an infinite vortex line, it is possible to obtain the analytical solution proposed earlier by Lamb and Oseen, this model is known as a Lamb-Oseen vortex (Leishman 2006). Knowing the vorticity field, the induced velocity of the vortex segment length L can be determined by Eq. 3 (Ignatkin et al. 2009): where D is the area where the vorticity field is calculated, dτ is the elementary three-dimensional volume and r 2 is the distance from the velocity calculation point to the point in the vorticity field (Fig. 1). Equation 2 does not have an analytical solution in finite limits. The numerical solution of Eq. 3 will lead to a large of computing resources costs. In order to reduce the calculation time, pre-calculated tables of induced velocities from a normalized vortex segment are used, depending on the length of the segment L, the degree of its diffusion (time t) and the distance from the segment to the point.
In the model showed in Fig. 1, rotor blade can make flapping motions relative to horizontal hinge.
To determine the angle of the blade flapping β at the current calculation step i, the Taylor series expansion up to the second derivative is used (Eq. 4): The angular blade flapping velocity is also determined through Taylor series expansion up to the second derivative. The angular acceleration of the blade flapping is defined as Eq. 5: β i =β i-1 +β̇i -1 •∆t+0.5•β̈i -1 •∆t 2 (4) β̈•I HH =M HH (5) where I HH is the moment of inertia of the blade relative to horizontal hinge and M HH is the moment from external forces (aerodynamic force, gravity and centrifugal force). The FWM developed by the authors has been validated in the article (Ignatkin et al. 2009) by comparisons of calculation data with experimental data and calculations based on other models. The model provides calculation of aerodynamic characteristics for both single and multi-rotor configurations, taking into account the aerodynamic interference.
The scheme of the calculation model of the coaxial rotor is shown in Fig. 2  The time step used in the calculations has been corresponding to the rotation of the rotor at an azimuthal angle of 6°. Calculations of each of the operating modes of the rotor were performed for the number of revolutions of the rotor equal to n = 30. The time spent on calculating a single regime was about 24 h. A high-performance personal computer equipped with a processor with 8 computing cores (16 threads) was used.

UNSTEADY REYNOLDS-AVERAGED NAVIER-STOKES METHOD
Numerical modeling based on the FVM was performed in the Ansys Fluent software. In the calculations, the blades were modeled with a rigid mounting to the rotor hub. The blade cone angles for each rotor thrust values were obtained on the basis of an FWM.
The URANS equations with the k-w SST turbulence model (Menter 1993) were used for creating a viscous turbulent flow of compressible gas.
The computational mesh had a block-structured topology (Fig. 3) and contained about thirty million hexahedral cells. The computational volume was divided into several cylindrical areas, each of which had its own mesh discretization characteristics (Fig. 3b). Close to the upper and lower rotors there were located the rotation areas of the upper and lower rotors (Fig. 3b) with the calculated mesh of the smallest scale. Under the rotor areas, there was a cylindrical area containing a mesh for the far vortex wake (Fig. 3b). The environment area was situated around these three areas. This area had the roughest computational mesh. The upper border of this region was located at a distance of 2.5 R from the rotor and the lower and lateral border at a distance of 4 R (Figs. 3a and 3b). A sliding interface condition was set at the boundaries of the rotor areas. Due to the periodic symmetry of the rotor, each computational mesh area was initially created for a sector of 120° (Fig. 3с) and then combined into a full cylindrical area (Fig. 3b). Representation of the computational volume as separate areas with different scales of the computational grid allowed to achieve a good balance between the accuracy of modeling and the number of cells. This is especially important because of extreme resource intensity of the method used. Figure 4 shows the computational mesh of the coaxial rotor. The rotor blades were divided into 20 equal sections (chambers) along the radius (Fig. 4a). The blade surface grid was constructed in such a way that the first node was located in the area of the viscous velocity profile Y + ≤ 1 (Fig. 4b). Numerical simulation of the rotor was carried out in an unsteady setting: taking into account changes in the flow parameters over time. The initial turbulence parameters were selected based on the conditions of the average intensity of the developed turbulent flow. The relative turbulent viscosity was assumed µ t /µ = 5. The value of the turbulent intensity was assumed 1%.
Preliminary test calculations were performed to study the mesh convergence. There were some meshes used with a different number of cells created by sequential, reducing the size of cells and increasing their number. The results of test calculations showed that the final mesh used in this work is sufficient for reliable calculation of aerodynamic characteristics of the coaxial main rotor under study.
The calculation of one mode for the number of revolutions of the rotor n = 25 took about 120 h with powerful supercomputer cluster. Figure 5 shows the performance diagrams of the coaxial rotor in hover с ТΣ = f(с QΣ ), obtained on the basis of two methods used. Performance diagram contains the sums of time-average values of thrust (с ТΣ = с ТUR + с ТLR ) and torque (с QΣ = с QUR + с QLR ) coefficients of the upper and lower rotors. There is a good agreement between the results of calculations using FWM and URANS.  Figure 6 shows the dependence of the rotor efficiency in hover figure of merit (FoM) on the values of с ТΣ /σ (where σ is the rotor solidity) obtained from calculations for both models, as well as experimental data (Vassiliyev et al. 2007). The FoM is the ratio of the ideal power required to hover to the actual power required (Leishman 2006). On the basis of the momentum theory (Leishman 2006), the values of FoM can be expressed through rotor thrust and torque coefficients as FoM = с ТΣ 3/2 /(2·с QΣ ). There is a good agreement between the results of calculations, both among themselves and with the experimental data. Figure 7 shows the dependence of the thrust coefficients of the upper and lower rotors on the values of с ТΣ /σ. There is a good agreement between the calculation data obtained on the basis of both methods used. The ratio between the thrust of the upper and lower rotors is in the range of 17−25% depending on the value of с ТΣ /σ, which corresponds to the experimental data (Bourtsev et al. 2001;Petrosian 2004).

RESULTS AND DISCUSSION
A feature of a coaxial main rotor is the mutual aerodynamic interference between upper and lower rotor. This is associated with the periodical passage of the blades over each other. As a result, there are instantaneous pulsations of aerodynamic loads, for example rotor thrust. Such pulsations in particular were analyzed experimentally (Deng et al. 2019) and by calculation (Singh and Friedmann 2017).    Figure 9 shows a fragment of this dependence for a single turn of the rotor. Here, in addition to the data obtained on the basis of the FWM, the results of calculations based on the URANS method are also presented. It can be seen that the results obtained on the basis of the FWM and the URANS method show a good match.  In Fig. 10a, the shape of the vortex wake in two projections is demonstrated. For ease of analysis, not all vortex wake grids are shown, but only blade tip vortices. They are shown in blue color for the upper rotor and in brown color for the lower one. There is clearly a compression of the vortex wake behind the rotor, which is noted in all works considering the rotor wake in the hover (Kim and Brown 2010;Yana and Rand 2012;Ridhwan et al. 2017). It is also seen that the helix pitch of the tip vortex descending from the upper rotor is much greater than that of the lower rotor vortex. This feature of the vortex wake of a coaxial rotor is described in a number of experimental (Nagashima and Nakanishi 1983) and computational works (Kim and Brown 2010;Yana and Rand 2012;Fernandes 2017). Figure 10b shows a visualization of the flow around the rotor made using airflow lines plotted in the OXY plane (Fig. 2). In accordance with the relative position of the upper and lower rotors blades relative to each other at the moment of calculation on the right side, the OXY plane passes through the blade axes and on the left between the blades. This circumstance affects the observed asymmetry of the flow regarding the rotor axis. Figure 10b also shows diagram of the vertical component of induced velocity at various distances from the rotor. The observed features of induced velocity distribution cause the differences noted above in the rate of descent of the blade tip vortices of the upper and lower rotors. As follows from the diagrams shown in Fig. 10b, the blade tip vortices of the upper rotor go down in a much more intensive induced flow. As shown in Fig. 10a, airflow lines show that due to the compression of the upper rotor flow, the coaxial rotor has a larger actual disk area than that of an equivalent conventional rotor with the same blade radius. A well-known result of this feature is high value of the efficiency (FoM values) of the coaxial rotor in hover (Ignatkin et al. 2018). Figure 11 shows the results of visualization of the vortex wake structure behind the coaxial rotor in hover for с ТΣ = 0.01, obtained on the basis of the URANS method. Visualization is performed using isosurfaces of vorticity. In Fig. 11a, the isosurface for vorticity value of ω = 50 s -1 is shown in isometry. Figure 11b shows a plane section with vorticity contours. The locations of vortex cores and their diffusion during the time are clearly visible. High vorticity value areas indicate discrete vortex bundles that exist under the rotor at a distance approximately equal to the blade radius and then slowly blur forming a continuous area of vorticity. In general, the presented visualization data by URANS modeling are in good agreement with the data of FWM calculation. Figure 12 shows diagram of calculated blade tip vortices trajectories of the upper and lower rotor for the values of с ТΣ ≈ 0.01 obtained on the basis of the FWM and URANS method. Experimental data (Nagashima and Nakanishi 1983) and calculated data (Kim and Brown 2010) are also presented. Here, ψ is the blade tip vortex azimuthal coordinate, y/R is the non-dimensional vertical coordinate describing the tip vortex descend and r/R is the coordinate describing the tip vortex radial compression. These results coincide with visualization data of the coaxial rotor vortex wake shape (Figs. 10 and 11). It can be seen that, according to the authors cited (Nagashima and Nakanishi 1983;Kim and Brown 2010), the speed of descent of the upper rotor blade tip vortices is significantly higher than from the lower one. At the same time, the results of calculations using the FWM and the URANS method are in good agreement with each other, as well as with experiment (Nagashima and Nakanishi 1983) and calculations (Kim and Brown 2010). The best match is obtained by comparing of the flow compression values r/R (ψ). For the upper rotor it was r/R ≈ 0.72 (FWM) and r/R ≈ 0.74 (URANS). For the lower rotor it was r/R ≈ 0.87 (for both methods). Free wake model calculations of the vortices trajectory of the lower rotor blade coincide with the data of calculations (Kim and Brown 2010) and experiments (Nagashima and Nakanishi 1983) a bit worse than calculations based on URANS.
In general, as well as for the data on visualization of the rotor flow, the results of calculations of the aerodynamic characteristics of the coaxial rotor, made on the basis of both used calculation methods, show a good match and agree with the data of other authors. This indicates the reliability of the methods used in modeling of coaxial rotor aerodynamics.
The high performance of the FWM developed and applied by the authors demonstrates the possibility of its application for a large volume of calculations based on conventional personal computers. The results obtained in this case, taking into account the model assumptions, allow to solve a significant range of problems in the aerodynamics of the coaxial rotor. The exceptions are individual tasks, such as tasks related to the calculation of aerodynamic characteristics of rotors with complex blade shapes and tasks related to the need of calculating the pressure distribution over the blade surface, including problems of aeroacoustics.
Modeling of helicopter rotor aerodynamics with the URANS based FVM method requires significant computational resources and a large amount of time to create and calculate grids. The use of FVM is very effective in solving problems related to the need for accurate modeling of the spatial flow around blades and the blade-vortex interaction (BVI). For example, problems connected with the influence of complex-shaped blade tips on the rotor aerodynamic characteristics, problems of studying the pressure distribution on the blade surface, etc. This method is also a powerful tool for obtaining data to refine vortex models.  Thus, each of the used approaches, both the FWM and the URANS based FVM method currently have their own niche. Therefore, their combination will be most effective in solving complex problems of helicopter rotor aerodynamics.

CONCLUSIONS
On the basis of two CFD methods, the original FWM model developed at Moscow Aviation Institute and the URANS method, a comparative study of the aerodynamic characteristics of the coaxial main rotor of Ka-226 helicopter in hovering modes was performed.
The results of visualization the vortex wake forms, rotor flow structures, blade tip vortices trajectories, as well as total and distributed aerodynamic characteristics of the rotor, including the rotor performance and FoM diagrams, have been obtained and analyzed.
The results are compared with each other and with the available experimental and computational results by other authors. A good agreement of the results was reached on the basis of both approaches, as well as their satisfactory coincidence with the data of other authors.
The obtained values of the rotor efficiency (FoM) in hovering mode are up to 0.77. Based on the analysis of the tip vortex trajectory, the flow compression under the rotors was determined: for the upper rotor, it was in the range of r/R ≈ 0.72−0.74 and for the lower rotor about r/R ≈ 0.87.
The performed comparative studies confirm the reliability of the applied models for modeling the aerodynamic characteristics of a coaxial rotor in hovering mode. Both methods used are currently promising and have their own niche in solving practical problems of coaxial propeller aerodynamics.
An analysis of the time spent on creating computational models and grids, the capabilities of both methods, as well as required computational resources, allowed to conclude that it is advisable to apply jointly the approaches considered in solving complex problems of coaxial rotor aerodynamics. An example is applying the FWM in this work for preliminary calculation of coaxial rotor parameters, such as the required blade pitch angles and the corresponding blade cone angles of the upper and lower rotors. This approach has significantly reduced the time and computational resources required for the URANS method.

DATA AVAILABILITY STATEMENT
All data sets were generated or analyzed in the current study.