，*，，2

1.National Key Laboratory of Rotorcraft Aeromechanics，College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，P.R.China；2.Rotor Aerodynamic Key Laboratory，China Aerodynamic Research and Development Center，Mianyang 621000，P.R.China

（Received 1 October 2020；revised 15 September 2021；accepted 20 October 2021）

Abstract:The discontinuous Galerkin（DG）method is established and innovatively conducted on accurately simulating the evolution of blade-tip vortex and the aerodynamic characteristics of helicopter rotor.Firstly，the Reynolds-Averaged Navier-Stokes（RANS）equations in rotating reference frame are employed，and the embedded grid system is developed with the finite volume method（FVM）and the DG method conducted on the blade grid and background grid respectively.Besides，the Harten-Lax-Van Leer contact（HLLC）scheme with high-resolution and low-dissipation is employed for spatial discretization，and the explicit third-order Runge-Kutta scheme is used to accomplish the temporal discretization.Secondly，the aerodynamic characteristics and the evolution of blade-tip vortex for Caradonna-Tung rotor are simulated by the established CFD method，and the numerical results are in good agreement with experimental data，which well validates the accuracy of the DG method and shows the advantages of DG method on capturing the detailed blade-tip vortex compared with the FVM method.Finally，the evolution of tip vortex at different blade tip Mach numbers and collective pitches is discussed.

Key words：blade-tip vortex；discontinuous Galerkin（DG）method；Navier-Stokes（N-S）equations；aerodynamic characteristics；helicopter rotor

0 Introduction

Intense vortices form at blade tips with the rotation of the helicopter rotor，and as the rolling up and shedding of vortices，the flowfield of rotor is dominated by the vortex wake.Therefore，helicopters are manipulated in the environment of rotor vortex flowfield，and the blade-tip vortex shedding from the front blade might strike on rear blade which causes strong blade vortex interaction（BVI）.Moreover，the blade-tip vortex falling downward and backward will impact the surface of fuselage and tail rotor and other parts of helicopter，resulting in serious aerodynamic interactions［1］.These phenomena will seriously affect the aerodynamic performance of the rotor，and cause the reduction of aerodynamic performance of fuselage and tail rotor.So it has theoretical and practical value to simulate the evolution of blade-tip vortex，which will further help analyze the aerodynamic characteristics of helicopter.

A lot of studies on the evolution of blade-tip vortex have been carried out.Previously，most of the researches on blade-tip vortex were carried out through experiment.In 2000，Heineck et al.［2］acquired the velocity fields for wake ages ranging from 0° to 270° in hover by particle image velocimetry（PIV）technique，and the experiment represented a major step toward understanding the detailed structure of a rotor wake.Martin et al.［3］further tested the blade-tip vortex by PIV and laser Doppler velocimetry（LDV）equipment.It can be found by contrast that PIV is much faster than LDV and has the advantage of being a whole-plane measurement technique.However，PIV is more sensitive to calibration errors and can suffer from spatial resolution issues if large regions of the flow are to be interrogated［4］.In addition，the whole process of experiment has a long period and high cost including model design，rotor manufacturing，test preparation，measurement process and data process.Compared with the experiment，the CFD method can quickly and fully capture the rotor flowfield for blades with different shapes，considering the viscosity and compressibility of the air，which makes up for the deficiency of the experiment.

The previous CFD methods［5-7］，such as the small disturbance equation and the full potential equation，are restricted by the irrotational，inviscid and isentropic assumptions，which cannot obtain the information of the vortex field.With the improvement of computing capability，the numerical simulation method of rotor in hover based on Euler equations and Navie-Stokes equations has gradually matured.Srinivasan and Baeder［8］employed the N-S equations based on the Roe scheme without specifying any wake models to establish a CFD method for solving the flowfield of rotor in hover.The analyzed aerodynamic performance of rotor is in good agreement with the experimental data，but it is difficult to capture blade-tip vortex with large velocity gradient due to the coarse grids and excessive numerical dissipation.In order to reduce the numerical dissipation，some researchers have proposed high-order numerical methods to improve the ability to capture the tip vortex of the rotor.Zhao et al.［9］established a rotor flowfield code based on the Roe scheme and thirdorder monotone upstream-centred schemes for conservation law（MUSCL）scheme.Hariharan［10］introduced a high-order weighted essentially non-oscillatory（WENO）scheme to effectively capture the blade-tip vortex.After that，different high-order numerical methods were successively developed to simulate the flowfield［11］.Finite volume method（FVM）is mainly used for the flowfield simulation of rotor at present，and most of FV methods improve the numerical accuracy by extending stencils［12］，which is difficult to choose the stencils to computer the derivatives in practice for unstructured grids.Unlike the FV methods，where the derivatives are reconstructed using the mean values of neighboring cells，the discontinuous Galerkin（DG）method［13］solves the equations for the derivatives in a manner similar to the mean variables and provides a unified formulation to represent the numerical polynomial solutions in each element，which reduces the numerical dissipation in the simulation，and has a great advantage in simulating the evolution of the blade-tip vortex.

It is easy to construct arbitrary order polynomials to improve the accuracy of numerical simulations and enhance the capture of flow details by employing the DG method.Besides，it is compact and each element is independent in DG method，so it can be highly parallelizable which could help to improve the efficiency of numerical method.Luo et al.［14］applied the DG method based on the Taylor series expansion to simulate the classical flow around the cylinder and the two-dimensional steady state airfoil flowfield，and the calculated results are in good agreement with the experimental data.Then，the ONERA M6 Wing［15］is numerically simulated，the pressure coefficient distribution is very close to the experimental values and more details of the flowfield are captured.At present，the DG method is mostly used in the numerical simulation of the two-dimensional flowfield and the three-dimensional fixedwing flowfield［16-18］，and has not yet simulated the flowfield of the rotating rotor.In order to improve the accuracy of capturing the tip vortex of the rotor，the numerical simulation of rotor flowfield based on DG method has good academic value.

In this paper，the DG method is established and innovatively conducted on accurately simulating the evolution of blade-tip vortex and the aerodynamic characteristics of helicopter rotor.Embedded grid technology is adopted and blade grid with good orthogonality is generated and appropriately refined based on Poisson equations.The finite volume method is used to numerically simulate the flowfield of the blade grid.In order to reduce the numerical dissipation of the background grid，the DG method is used to simulate the evolution of the blade-tip vortex.The blade-tip vortex in hover is numerically simulated using the established method and the parameterized calculation of Caradonna-Tung（C-T）rotor is carried out with different blade-tip Mach numbers and collective pitches.

1 Computational Methods

1.1 Embedded grid methodology

Embedded grid methodology［19］is adopted to discretize the calculation domain of rotor flowfield in which O—C—O type body-fitted grid is generated around rotor blade，and the cylindrical Cartesian grid is adopted for the background grid［9］.Fig.1 shows the the whole embedded grid system.

Fig.1 Embedded grid system around rotor

FVM is conducted on the blade grid.In order to reduce the numerical dissipation and capture more revolutions of tip vortex with minimal diffusion，the DG method is conducted on the background grid which is dominated by the vortex wake.Considering that the viscous effects on background gird have little influence on the revolution of rotor flowfield，the background grid away from the rotor surface is treated as inviscid in order to reduce the calculation consumption.

1.2 Governing equations

The Navier-Stokes equations in rotating reference frame are employed as governing equations for simulating the quasi-steady flowfield of rotor in hover，and the expression is

whereUis the conservation variable，considering the contravariant velocity－Ωof cell surface and original convective flux，the modified convective flux vector，the viscous flux vector，Qthe source term，Sthe boundary ofV，andnthe unit outward normal vector to the boundary，which can be written as

whereρ，p，Eandhdenote the density，pressure，specific total energy and the specific total enthalpy of the fluid，respectively.uiis the velocity of the flow in the coordinate direction，Ωthe rotating angular velocity of the rotor，τthe viscous stress tensor，uitthe rotating velocity of rotor in the coordinate direction，Θthe work of viscous stresses and heat conduction in the fluid，andδijthe Kronecker delta function.

Aimed at well simulating separated airflow over surface of rotor，one-equation Spalart-Allmaras（S-A）［20］model is used to simulate turbulent flows.

1.3 DG method

1.3.1 DG method on background grid

Euler equations are employed as governing equations on background grid and the governing equations are discretized using a DG formulation.Firstly，the following weak formulation is introduced，which is obtained by multiplying the above conservation law by a test functionW，integrating over the domainV，and performing an integration by parts

Assume that the domainVis subdivided into a collection of non-overlapping elementsΩe.Then，the broken Sobolev spaceconsists of discontinuous vector-values polynomial functions of degreepis as follows

wheremis the dimension of conservative state vector，and the Sobolev space can be written as

whereαdenotes a multi-index anddthe dimension of space.After that，the following semi-discrete form by applying the weak formulation on each elementΩecan be obtained as

whereUhandWhdenote the finite element approximations to the analytical solutionUand the test functionW，respectively.Assume thatBis the basis of polynomial function of degreep，then the system ofNequations can be achieved as

whereNdenotes the dimension of the polynomial space.A Taylor series expansion at the centroid of the cell is adopted to represent the numerical polynomials solutions in each element，and the polynomial solutions can be expressed as

whereB1=1，，denotes the mean value in the cell，(xc，yc，zc) the centroid of the cell，and（Δx，Δy，Δz）the cell size.

In this paper，the Harten-Lax-Van Leer contact（HLLC）［21］approximate Riemann solver is adopted for the inviscid flux spatial discretization，which is applicable to both structured and unstructured grids and has low storage requirements.Gaussian integration quadrature formula is employed in discontinuous finite element calculations since they present the best accuracy for a given number of points.The HLLC flux is defined by

whereUlandUrare the solution polynomials at the left and right states of the cell interface.

1.3.2 Time integration

In this paper，the explicit method is employed in the temporal discretization of the semi-discrete system.Here，the explicit three-stage third-order TVD Runge-Kutta scheme［22］is adopted

whereMdenotes the mass matrix and it has the block diagonal structure.R(U)is the residual vector obtained by flux calculation.

1.3.3 Information transmission

The flowfield information of the contributing element of background grid needs to be transferred to the corresponding outer boundary grid of the blade.Because the DG method used in the background grid obtains high-order accuracy by polynomial interpolation in the element，which can calculate the conservative variables accurately at a certain point，it is more accurate when transmitting information to the blade grid.As shown in Fig.2，elementAis an outer boundary grid of the blade，and the center point a of the element falls on elementBof the background grid.On the background gird，the mean value and the partial derivative of the conservative variables are used to calculate the actual value at pointawhich is transferred to elementAof the blade grid.The specific formula is as follows

Fig.2 Schematic diagram of value transfer from background grid to blade grid

2 Numerical Examples

2.1 Validation of CFD method

The classic C-T rotor［23］is taken as validation case for simulation of aerodynamic characteristics of rotor.The rotor has two rectangular blades with aspect ratio of 6.NACA0012 airfoil is adopted in the rotor blade without twist.The hover condition with a collective pitch of 8°，and blade-tip Mach numberMatip=0.439 is conducted for numerical simulation.The number of grids around the blade is 213×49×110，and the number of background grids is 121×150×120.

Fig.3 shows the calculated and experimental pressure coefficient（Cp）distribution of several sections.In Fig.3，x/cis the chord position of the blade，andr/Ris the radial position of the blade section.It can be seen that the pressure coefficient distribution calculated by the method used in this paper is in good agreement with the experimental data in each section of the rotor blade，which shows the effectiveness of the method in this paper.

Fig.3 Comparison of pressure coefficient distribution at different sections of C-T rotor(Matip=0.439)

In order to capture the evolution process of blade-tip vortex better，the effects of special bladetip shape and twist distribution should be excluded.Therefore，the C-T model rotor is chosen to be simulated.The Jameson scheme，the fifth-order Roe-WENO scheme and the method established in this paper are used to simulate the hover condition respectively.The background grid is properly refined in order to capture more revolutions of the blade-tip vortex，and the number of grids around the blade is 259×49×180，and the number of background grids is 151×253×160.Table 1 shows the rotor thrust coefficient calculated by three different methods respectively.It can be seen that thrust coefficient calculated by the method adopted in this paper is more comparable to the experimental data，which illustrates the accuracy of the calculation of the aerodynamic performance of the rotor，and the error of thrust coefficient is less than 4%.

Table 1 Aerodynamic performance comparison between different methods

The vorticity iso-surfaces of C-T rotor at the collective pitch angle of 8° are shown in Fig.4.It can be seen that the blade-tip vortex simulated by Jameson-Schmidt-Turkel（JST）scheme is not complete，and the vortex sheds about 360° until it dissipates completely.The wake age captured by fifth-order Roe-WENO scheme can reach 720°.On the contrary，the blade-tip vortex simulated by the method established in this paper sustains to more than 800°before it dissipates completely.Fig.5 shows the vorticity iso-surfaces on the section simulated by the JST scheme and the DG method，respectively.The same grids are adopted in three different methods，which indicates the ability of Discontinuous Galerkin method used in this paper to reduce the numerical dissipation.It can be found that the method established in this paper could simulate the vortex flowfield of the rotor more accurately through the comparison.

Fig.4 Comparisons of vorticity iso-surfaces among different methods

Fig.5 Comparison of vorticity iso-surfaces

Figs.6（a，b）show the predicted trajectories of the blade-tip vortex in vertical and radial directions respectively，wherey/R，r/Rdenote the vertical position and the radial position，andφis the wake age.It can be seen that the vertical and radial displacements of blade-tip vortex simulated by the method established in this paper are in good agreement with the experimental data for a wake age up to 360°.

Fig.6 Blade-tip vortex location of C-T rotor

Fig.7 shows the vorticity iso-surfaces on two typical circumferential sections where the wake age between the sections differs by 90°.It can be seen from Fig.6 and Fig.7 that the blade-tip vortex will be contracted to the root of the blade when it develops downstream along the vertical direction and the contraction degree of the blade-tip vortex will be slowed down when wake age is greater than 360°，then，the tip vortex begins to expand outward.With the increase of wake age，the vorticity value of blade-tip vortex continues to decrease，resulting from physical dissipation and numerical dissipation.When the blade-tip vortex develops downward to a certain extent，the vertical speed of vortex-core motion increases.This sudden increase in downward convection is the result of the downwash from the blade and the evolving tip vortex on the first passage of the vortex.The dissipation of vortex continues to increase until the blade-tip vortex dissipates completely.

Fig.7 Vorticity iso-surfaces of different sections

2.2 Analysis on evolution of blade-tip vortex

In order to study the evolution of the blade-tip vortex of the rotor，the parameterized calculation of C-T rotor is carried out at different blade tip Mach numbers and collective pitches.Mach number ranges from 0.439 to 0.794，and collective pitch of rotor ranges from 5°to 9°.

Fig.8 shows the vertical positiony/Rand the radial positionr/Rof the C-T rotor tip vortex core with a collective pitch of 8° and blade tip Mach numbers are 0.439，0.526，0.612，0.727 and 0.794，respectively.It can be seen from Fig.8（a）that the larger the Mach number is，the slower the vertical descent speed of the vortex will be.It can be seen from Fig.8（b）that the larger the Mach number is，the larger the wake age at which the tip vortex begins to expand outward is，and its contraction speed is relatively slow.

Fig.8 Blade-tip vortex location at different blade-tip Mach numbers

Fig.9 shows the vertical position and radial position of the C-T rotor tip vortex core with the blade tip Mach number is 0.794 and the collective pitches are 5°，6°，7°，8°，and 9°，respectively.In Fig.9（a），the greater the collective pitch is，the stronger the downwash induced by the rotor is，and the faster the vertical descent speed is.It can be seen from Fig.9（b）that the collective pitch has little influence on the radial motion of the blade-tip vortex.Similarly，the wake age reaching the contraction limit is little affected by the collective pitch of rotor.

Fig.9 Blade-tip vortex location at different collective pitches

3 Conclusions

A numerical simulation method based on the DG method has been developed for the prediction of rotor flowfield in hover in this paper.The aerodynamic performance and evolution of the blade-tip vortex of C-T rotor are simulated and analyzed by the established method，and the following conclusions can be obtained：

（1）The pressure coefficient distribution calculated by the method established in this paper is in good agreement with the experimental data，and the error of the thrust coefficient is less than 4%.

（2）The comparison of the simulated results indicates that the DG method has the characteristic of lower numerical dissipation compared with the JST scheme and the Roe-WENO scheme，and it could predict the characteristics of blade-tip vortex more accurately.

（3）There is a contraction limit during the downward development of blade-tip vortex.After the limited position，the tip vortex begins to expand outward and the speed of vortex-core motion increases.The dissipation of vortex continues to increase until the blade-tip vortex dissipates completely.