Abstract
This paper presents a modified fifth-order WENO-HLLC Riemann solver for the transonic compressible viscous flows around the helicopter rotor in hover. The HLLC approximate Riemann solver is proposed to discrete the convection term involving grid velocity of the Navier-Stokes equations. In order to solve the interface flow accurately, a modified fifth-order WENO scheme is presented by designing the smoothness indicators based on norm measurement. The improved WENO scheme can provide the optimal approximation order even at critical points. Numerical accuracy and robustness are validated by several benchmark inviscid flow problems. Then the numerical properties of the WENO-HLLC solver in conjunction with the implicit LU-SGS time integration method with high efficiency are further validated by simulating transonic viscous flows over RAE2822 airfoil and ONERA-M6 wing. The results show that the accuracy of calculating shock, discontinuity, and the vortex is significantly improved. Finally, the method is developed to compute the transonic vortex flow around the helicopter rotor with a domain discretized by overset grids. The results indicate that the proposed method is very robust and effective in acquiring high resolution for vortex wake.
1. Introduction
Rotor is the key component of a helicopter, which provides the lift to overcome the weight of the helicopter and the thrust during forward flight and the control force of helicopter balance, flight, and lifting. The aerodynamic performance of the rotor largely determines the flight quality, stability, maneuverability, vibration, and noise level of the helicopter. Obviously, in helicopter technology, rotor aerodynamics plays an important role. However, the existence of a tip vortex makes the flow field of the rotor extremely complex [1, 2]. In the field of helicopter aerodynamics and computational fluid dynamics (CFD), the accurate prediction for rotor flow field with vortex wake is still one of the most complex and challenging problems [3].
The rotor always works in or near the tail vortex system with complex distortion. Especially in hover and low speed forward flight, the relative velocity of the fluid at each spanwise position of the blade is different. As a result, the lift is mainly concentrated near the tip. And the lift must be zero at the tip, so the spanwise lift gradient near the tip is very large, resulting in a strong blade tip vortex. The strong tip vortex forms a downward spiral vortex system near the rotor disk, which dominates the whole flow field and seriously changes the effective angle of attack of the blade and then affects the aerodynamic performance of the rotor. In addition, the tip vortex has an important influence on the aeroelasticity and noise performance of the rotor. There is a phenomenon of blade vortex interference (BVI) in the rotor flow field [4]. In some flight conditions, there will be strong BVI. BVI is not only the main source of helicopter noise but also an important factor of helicopter vibration.
The important feature of flows over helicopter rotor, as mentioned above, is that the tip vortex plays a leading role in the flow field, and the movement and trajectory change of the vortex wake from the blade has an important impact on the aerodynamic characteristics of the rotor. In numerical calculation, the numerical dissipation caused by the truncation error of the numerical method would lead to distortion or failure to capture vortex wake. Therefore, accurate capture of vortex wake, shock, and discontinuity has always been a great challenge for computational mathematics and CFD. In the paper, a low dissipation numerical method will be developed to calculate the transonic vortex field of the hovering rotor.
Considering the computational cost and the uncertainty of most Godunov-type methods in computational details, the research of numerical methods for compressible flows mostly focuses on the construction of approximate Riemann solutions. In the approximate Riemann solver, HLLC (Harten-Lax-van Leer-Contact) is one of the simplest equilibrium solvers that can keep the isolated shock, discontinuity, and shear wave accurately in theory [5]. HLLC Riemann solver is developed on the basis of the HLL scheme. The central idea of the HLL scheme is to assume that the Riemann solution consists of three constant states separated by two waves [6]. The assumption of a two-wave configuration is correct only for hyperbolic systems of two equations but not for the larger systems such as Euler equations. Toro, Spruce, and Speares presented the HLLC scheme by remedying intermediate waves in the HLL scheme, whereby the missing contact and shear waves in Euler equations are restored. For the intermediate region of the Riemann solution fan, HLLC involves two-star states instead of one-star states in the HLL. Literature [7] showed that the HLLC solver can solve the problems of isolated shock, discontinuity, and shear wave, and initial positive density and pressure are maintained in the Riemann solution fan. HLLC has been widely concerned and developed in many fields such as multiphase flows, compressible flows, incompressible fields, magnetohydrodynamic, and so on [8–13]. In the paper, the HLLC Riemann solver in combination with an overlapping grid system is developed to calculate three dimensional compressible viscous transonic flow field of the hovering rotor.
To increase the accuracy of first-order HLLC scheme, an improved fifth-order weighted essentially nonoscillatory (WENO) scheme is proposed to reconstruct the left and right states of the Riemann solution. Jiang and Shu developed the classic WENO scheme known as the WENO-JS scheme, whose smoothness indicators were constructed based on norm measurement [14]. The weighting relationship of the fifth-order WENO-JS scheme is . Based on the fifth-order WENO-JS scheme, Borges, Carmona, Costa, and Don constructed a global smoothness indicator on the whole stencil including five cells and presented the fifth-order WENO-Z scheme by combining the global smoothness indicator and the local smoothness indicators to build new weightings of satisfying [15]. The WENO-Z scheme also provides a new idea of constructing smoothness and weighting factors [16, 17].
Fifth-order WENO-NS with first-order local smoothness indicators was introduced by using norm measurement [18]. Due to the unbalanced contribution of three stencils for WENO-NS, the WENO-P scheme was presented by introducing a parameter into local smoothness indicators [19]. On the basis of the WENO-P scheme, the MWENO-P scheme was developed by constructing a global smoothness indicator [20]. The construction of these WENO schemes all adopts the way of combining local smoothness indicators with a global smoothness indicator. These schemes could achieve high-order accuracy in smooth areas and guarantee a certain accuracy at extreme points.
Inspired by the above-reconstructing ideas, the paper applies the norm measurement to construct local smoothness indicators with second-order accuracy. Meanwhile, a high-order global smoothness indicator without the first to the third derivative is proposed. Modified weighting factors are reconstructed to satisfy .
The organization of the article is as follows. Section 2 is dedicated to the governing equations for the compressible viscous flows around the hovering rotor. Section 3 presents the numerical method involving the HLLC Riemann solver and the improved WENO scheme. The numerical results for all cases are conducted in Section 4. The conclusions are outlined in the last section.
2. Governing Equations
For the hovering rotor, the flow field is rotationally symmetric. Therefore, the flow field can be transformed into a steady state when the coordinate system is established on the rotor. The integral formulas of Reynold Averaged Navier-Stokes (RANS) equations based on the noninertial rotating coordinate system can be written as follows:where is a control volume with as the boundary, and is the unit outward normal of . is the conserved variable vector, is the convective flux vector, is the viscous flux vector, and caused by rotation can be treated as a source term, given as follows: and are the density, pressure, total energy, total enthalpy, and velocity vector of the fluid, respectively. is the velocity vector of the grid. denotes the angular velocity vector of the rotating frame. are the mean viscous stresses. , and are given by the following:
Here, and denote the temperature and thermal conductivity coefficient. Aimed at well simulating separated flow over the surface of rotorcraft, the one equation Spalart-Allmaras (S-A) turbulence model [21, 22] is employed to calculate the turbulence viscosity.
Different dimensionless systems will get different dimensionless equations. In the paper, free stream density , free stream pressure , free stream temperature , and mean aerodynamic chord length of rotor or wing are chosen as the references for the density, pressure, temperature, and length. The advantage of this dimensionless system is that the dimensionless RANS equations have the same form as (1), except that the viscosity coefficient and thermal conductivity coefficient are multiplied by .
When the source term and grid velocity are both zero, the governing equations become the mathematics model of flow around a fixed wing.
3. Numerical Formulation
3.1. The Finite Volume Method
The finite volume method is used to discretize (1) over each control volume (grid cell) with a structured grid in space, resulting in a system of ordinary differential equations in time:
With residuals and are the volume-averaged of variables and , respectively. Subscript is taken as which denote the indexes of the control surfaces and is the surface area.
The viscous fluxes are approximated based on the second-order central finite difference method. The numerical method for the inviscid fluxes (i.e., convection fluxes) is described in detail in the following subsection.
For temporal discretization, considering that the computational grid is large and the evolution of rotor wake takes a long time, we adopt the implicit LU-SGS time integration method [23] to solve (4):
Here, . and are as follows:with
, , are Jacobian matrices on outward normal , , of interface,
Respectively. , and are their spectral radiuses, and for example,, where is sonic speed. is a constant greater than 1 and taken as 1.05 here.
The above LU-SGS scheme is unconditionally stable. This implicit method with high efficiency completely avoids solving a huge linear system of equations in the whole computational domain, which effectively reduces computational complexity.
3.2. HLLC Riemann Solver
Since the numerical viscosity of the method for solving the convection term is too large, the rotor wake could dissipate or disappear in the transport process, resulting in the distortion of the tip vortex. It is of great significance to find a numerical method that can accurately calculate the complex flow field characterized by vortices in the field of rotor aerodynamics. In the paper, the HLLC approximate Riemann solver [5, 7] is proposed to solve the convection flux. As the governing equation of rotor flow is established in the noninertial rotating coordinate system fixed to the rotor, it is necessary to modify the formula of the HLLC Riemann solver.
Figure 1 shows the structure of the HLLC Riemann solver. In the figure, and denote the initial left and right states. and denote the left and right states to the contact discontinuity. and are the fastest signal velocities, and is the speed of a middle wave that can be regarded as a contact discontinuity. Thus HLLC scheme has the advantage of capturing contact discontinuity. Furthermore, the HLLC scheme does not need entropy correction to avoid instabilities.
The HLLC numerical fluxes in the interface are defined as follows:with the intermediate fluxes and still to be determined. The Rankine-Hugoniot conditions are imposed on each of the waves of velocity , and , resulting in the following:
These three equations involve four unknown vectors , , , . It can be seen from (10-12) that finding solutions for the two intermediate state vectors and is sufficient. In order to solve the algebraic problem, some extra conditions need to be imposed. The pressure and normal velocity of the interface are as follows:where .with . From expressions (10), (12), (14), the pressure in the “∗” domain can be obtained as follows:
Since , we have the following:
Algebraic manipulation of (10 and 12) and applying the values and from (15) gives the intermediate states as follows:for . The corresponding intermediate fluxes can be written as follows:
The left and right wave speeds can be simply estimated bywhere is the sonic speed and the superscript tilde denotes the Roe average [24].
3.3. The Fifth-Order WENO Reconstruction
3.3.1. General Frame of the Fifth-Order WENO Scheme
The reconstruction is performed dimension by dimension, so the one-dimensional space domain is considered without loss of generality. The calculation region is divided into uniform grid cells with a uniform interval length . Considering that the right state can be obtained by symmetry, we only discuss the high-order reconstruction of the left state at the cell interface . represents an arbitrary component of the vector . The fifth-order WENO reconstruction is carried out on three stencils:
By interpolating the three nodes on each stencil, three different reconstructions for the left state values are produced as follows:which satisfy
Numerical oscillation could happen in the region containing discontinuities. To avoid the problem as much as possible, an essentially nonoscillatory (ENO) scheme is to select the smoothest one among all the stencils [25]. Obviously, the three interpolation stencils provide the same accurate flow field information in smooth regions, but the ENO scheme finally only chooses one and has third-order accuracy, whereas, on the three stencils over five cells , the fifth-order central upwind interpolation scheme can be obtained as follows:with .
Then,
To achieve fifth-order accuracy and maintain the property of essentially nonoscillatory, WENO scheme [26] takes a convex combination of the reconstructions on the three stencils as follows:where weighting factors are nonlinear functions instead of constants. It is required that the weight of the discontinuous stencil should be as small as possible, and the accuracy at the extreme points should not be reduced. The weighting factors are defined as follows [14]:with
The introduction of is to avoid the denominator being zero. Smoothness indicator , a function of , is used to measure the smoothness of the stencil . WENO flux (25) can be expressed as follows:
Let and denote the weighting factors of and , respectively. Substituting (22) and (24) into (28), we obtain the following:
In the same way,
For the fifth-order WENO scheme, the expressionshould be true, and then the sufficient condition is as follows:
3.3.2. The Construction of the Improved Fifth-Order WENO Scheme
In [18], Ha et al. introduced a novel smoothness indicator by defining to be as follows:where the coefficients are obtained by solving the system
Then, are expressed as follows:
follows that
Based on the terms , one kind of smoothness indicator with norm instead of norm for the WENO-JS scheme was given by the following [18]:
By Taylor expanding for smoothness indicators , it follows that
Thus,
The nonlinear weighting factors were given by the following:where
This fifth-order WENO scheme is WENO-NS [18].
Because three stencils provide an unbalanced contribution to the evaluation point as shown in (38), based on WENO-NS, Kim, Ha, and Yoon [19] proposed a WENO-P scheme by introducing a parameter into the formulation of local smoothness indicators. On the basis of the WENO-P scheme, a fifth-order WENO scheme named as MWENO-P [20] was developed by constructing a global smoothness indicator. In the paper, by analyzing the smoothness indicators of the WENO-NS scheme, we will construct more optimal ones. The Taylor expansions of follow that
And
Thus, we construct local smoothness indicatorsso as to further reduce the weight of the stencil containing discontinuities. Their Taylor expansions can be written as follows:
To ensure that three stencils provide balanced contribution to the evaluation point , the smoothness indicators are further modified to the following:where the parameter is used to carry out the trimming.
In this work, the main idea of constructing a global smoothness indicator is to make the condition (32) hold even at the extreme value. Thus, we take a new parameter as follows:
Its Taylor expansion can be written by the following:
Then the new smoothness indicators are defined as follows:
Furthermore, the new nonlinear weighting factors are given by the following:
For the convenience of description, we refer to this WENO scheme as WENO-NZ scheme here.
3.3.3. Convergence Order of the New Fifth-Order WENO Scheme
As the discussion above, the sufficient condition for fifth-order convergence of the WENO scheme is .
First, consider the convergence order when there are no extreme points. Let , and then the local smoothness indicators can be written by the following:with , or from (45 and 46). Let for (48), and then . Thus,
From (50) and , the relation between the new weights and the ideal weights is as follows:which satisfies the condition .
Second, for the convergence order at an extreme point where and , the local smoothness indicators can be written by the following:with , or . Similarly, it can be obtained that
And thenwhich also satisfies the condition .
4. Results and Discussion
The accuracy and the robustness of the fifth-order WENO-NZ scheme are validated by computing a one-dimensional scalar convection equation, one-dimensional and two-dimensional inviscid compressible flows. For these inviscid problems, the third-order TVD Runge-Kutta scheme [27] is utilized for temporal discretization. Then, the numerical properties of the proposed WENO-HLLC method are studied for steady transonic viscous compressible flows over the RAE2822 airfoil and ONERA-M6 wing. To improve computational efficiency for the Navier-stokes equations and the large grid and the realization of full evolution of complex flows, the time integration is carried out by using the implicit LU-SGS method with high efficiency. Moreover, the method is applied to the steady transonic vortex flow around the helicopter rotor with a domain discretized by overset grids. Results are validated against experimental data and the results of the classical WENO-JS and WENO-Z schemes based on the HLLC solver. We take for WENO-JS and for WENO-Z and WENO-NZ. We have chosen for (45) and for (46) through a variety of numerical tests.
4.1. One-Dimensional Scalar Convection Problem
We check if the proposed method could achieve the desired optimal order by computing several typical initial value problems of one-dimensional scalar convection equation, which can be described by the following:
The numerical result will be compared with the exact solution . Errors are measured by and norms
Case 1. To test the accuracy and convergence order of the proposed method, we first consider the initial conditionThe solution is arbitrary smooth and has two critical points. The temporal accuracy and spatial accuracy are third-order and fifth-order in the numerical algorithm, respectively, so the time step is taken as . The periodic boundary condition is adopted and the results at are discussed.
Tables 1 and 2 contain the and errors and the convergence order for the WENO-JS, WENO-Z, and WENO-NZ schemes. It can be observed that with the grid refinement from 20 to 320, the errors gradually decrease, and whether error or error of both WENO-NZ and WENO-Z schemes reach about that is smaller than that of WENO-JS scheme, the WENO-NZ scheme approaches the desired convergence order.
Case 2. For the one-dimensional scalar convection equation, consider the initial conditionwhich can be regarded as a small disturbance added to Case 1. This problem still has two critical points where the first and third derivatives vanish, but the second derivative is nonzero. The computation conditions are the same as those of Case 1. As shown in Tables 3 and 4, for WENO-JS, WENO-Z, and WENO-NZ schemes, the and errors and the convergence order are given. It can be seen that WENO-Z and WENO-NZ gradually achieve the desired convergence order as the grid number ranges from 20 to 320, and their errors are smaller than that of WENO-JS. The error of Case 2 is larger than that of Case 1 under the same calculation circumstance, which indicates that the error will increase with the increase of flow complexity.
Case 3. This initial valve problem involves a jump discontinuity, which is described by the initial conditionThe number of the uniform grid cells is 200, and the Courant number is set to be . For the WENO-JS, WENO-Z, and WENO-NZ schemes, their approximate solutions at are plotted in Figure 2 against the exact solution. By observation, the WENO-NZ scheme performs better than the other two schemes near the discontinuity.
4.2. One-Dimensional Sod’s Shock Tube Problem
The governing equation of Sod’s problem is a one-dimensional compressible Euler equation with the following initial conditions of Riemann type
The solution contains a left rarefaction wave, a right traveling contract contact wave, and a right shock wave. The calculation interval is uniformly divided into 200 cells, and the Courant number is taken as . Transmissive boundary conditions are taken. Figure 3 shows the density profiles at for three fifth-order WENO schemes and the reference solution. The local enlarged plots show that WENO-NZ can capture shock better than the other two WENO schemes.
4.3. One-Dimensional Shock Entropy Wave Interaction Problem
The one-dimensional shock entropy wave interaction problem of Shu-Osher [28] is specified by the initial conditions
The problem describes the interaction of a right moving Mach 3 shock with a sine entropy wave. Figure 4 illustrates the density results at with a mesh size of N = 200 against the reference data regarded as an “exact” solution that is calculated by WENO-JS with a fine grid N = 4000. All the schemes can capture the compressed wavelike structures behind shock well, whereas the result of WENO-NZ is much less dissipative than those of WENO-JS and WENO-Z schemes.
4.4. Two-Dimensional Vortex Transport Problem
Since the transportation problem of two-dimensional isentropic vortices in free flow [29] has an analytical solution and no dissipation and dispersion errors, the computational performance could be evaluated by the capability of capturing vortex shape and vortex strength. The calculation domain is . The initial flow is given bywith free stream density, free stream pressure, and free stream velocities in two directions set as , , , , respectively. The distance from the point to the vortex core is expressed as . represents the dimensionless vortex strength.
Based on periodic boundary conditions, flow fields are simulated by using equispaced Cartesian grid , and , respectively. We observe the results after two periods of vortex propagation. Figure 5 gives error orders for densities on WENO-JS, WENO-Z, and WENO-NZ schemes. It can be seen that the error of WENO-NZ is the smallest, followed by the one of WENO-Z. The pressure contours of the isentropic vortex convection are illustrated in Figure 6 from which we can observe that when grids are sparse, especially for the grids, WENO-NZ has the highest resolution for vortex among the three methods, and WENO-Z has less numerical dissipation than WENO-JS. With the increase of grid number, the calculation results become better. There are almost no differences between the numerical results and the analytical solution for the grids.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
4.5. Steady Transonic Flow over RAE2822 Airfoil
The steady transonic flow over RAE2822 airfoil is solved using WENO-JS, WENO-Z, and WENO-NZ schemes. A C-type mesh of 201 × 61 points is used to discretize the domain, with outer boundary 10 chord lengths away. For the boundary layer, the grid spacing at the airfoil surface is about 10−6 times the chord length. The angle of attack is , the freestream Mach number is 0.731 [30], and the Reynolds number based on the chord is 6.5 × 106. No-slip condition is applied on the airfoil surface, and characteristic boundary conditions are enforced on the outer boundaries. The symmetric boundary condition is used at the cutting line along the trailing edge of the airfoil.
The density residual based on norm is shown in Figure 7. When the iteration is less than 8000 steps, the residuals of the three WENO schemes converge to five orders of magnitude, which indicates that they have good convergence performance.
Figure 8 shows the pressure coefficients at the airfoil surface, where the calculated results from three schemes are validated against experimental data [31]. It can be seen from Figure 8(b) that the WENO-NZ scheme can capture the shock better than WENO-Z and WENO-JS schemes. Considering that there is little difference among the three schemes, only pressure contours and streamlines for the WENO-NZ scheme are illustrated in Figure 9. The supersonic flow region on the upper surface and the shock are clearly visible. The above calculation results also show that the WENO-NZ scheme is validated for the viscous compressible flow over a domain discretized by a curvilinear grid.
(a)
(b)
4.6. Steady Transonic Flow around ONERA-M6 Wing
The fifth-order WENO schemes are validated for a classical three-dimensional flow problem through solving the transonic flow over the ONERA-M6 wing. As shown in Figure 10, the domain is discretized by a C-H type mesh with 161 points in the wrap-around direction, 51 points in the normal direction, and 55 points in the spanwise direction. The angle of attack is , the Reynolds number based on the mean aerodynamic chord is 11.72 × 106, and the freestream Mach number is 0.8395. No-slip wall boundary conditions are applied on the wing surface. Characteristic-based freestream boundary conditions are used at far-field boundaries. The symmetric boundary condition is enforced on the cutting surface of the trailing edge of the wing. The wing is assumed to be symmetrical about the root.
(a)
(b)
Figure 11 shows the density residual history curves for the WENO-JS, WENO-Z, and WENO-NZ schemes. The solution is iteratively updated until the residual converges to four orders of magnitude, and at this time, the solution has reached a steady state. In this process, the WENO-Z scheme with the fastest convergence speed iterates 5215 steps, followed in turn by the WENO-NZ scheme with 5890 iterative steps and the WENO-JS scheme with 6005 iterative steps.
Figure 12 shows the pressure coefficients on the wing surface at six spanwise locations, including , , , , , obtained by the three WENO schemes and experimental data [32]. The results for the three WENO schemes show a good agreement with experimental data. There is no significant difference among the three methods. But it can be observed from the local enlarged plots of the spanwise locations and , as shown in Figures 12(g)-12(h), that WENO-NZ captures the shock better than WENO-Z and WENO-JS. Figure 13 illustrates the pressure contours and the streamlines based on the WENO-NZ scheme. From Figures 12 and 13(a), two shocks can be seen clearly on the upper surface on the wing root side, and they get closer and closer with the approach to the wing tip. Meanwhile, Figure 13(b) shows the flow separations on the upper surface on the wing tip side.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
4.7. Transonic Flow over Caradonna-Tung Rotor in Hover
The improved WENO scheme is finally applied to solve compressible viscous flow around the hovering helicopter rotor, and just like the above examples, the new method is compared with WENO-JS and WENO-Z schemes. The classical Caradonna-Tung rotor with two blades is chosen [33]. Blade tip Mach number is , pitch angle is, and Reynolds number is 3.93 × 106. The coordinate system is built on the rotor, and then the flow is translated into a steady state. Because of the symmetry of the flow field, only one blade is considered. If the computational grid is also symmetrical, periodic boundary conditions can be applied directly on the plane of symmetry, which avoids the numerical dissipation caused by interpolation. In addition, the tip vortex and vortex wake generated from the interaction between rotor and fluid always hover below the rotor and drag far away. Therefore, overset grids composed of a C-H rotor grid and cylindrical background grid are adopted to solve the complicated flows. The rotor grid includes 199 points in the wrap-around direction, 49 points in the normal direction, and 61 points in the spanwise direction. The background grid includes 101 points in the circumferential direction, 111 points in the radial direction, and 161 points in axial direction. Figure 14 shows the sketch picture and the real grid is finer. Obviously, the overset grids can avoid distorting the mesh too much, and the background grid is not only conducive to the implementation of periodic boundary conditions on the symmetric plane, but also conducive to the capture of vortex wake.
(a)
(b)
4.7.1. Boundary Conditions
No-slip wall boundary conditions are enforced on the blade surface. Wall condition is applied to the radial inner surface of the background grid, where no-penetration conditions are used since the domain for the background grid is assumed to be inviscid. In addition, artificial boundary conditions, periodic boundary conditions, and modified far-field boundary conditions are also involved.
(1) Artificial Boundary Conditions. The outer boundary of the rotor grid is called the artificial outer boundary. Because the rotor grid is completely embedded in the background grid, the information on the artificial outer boundary is obtained by interpolating the information of the contributing cells on the background grid. For the background grid, a hole domain shown in Figure 15 is dug out. A layer of cells close to the hole area constitutes the artificial inner boundary, also called the hole boundary, whose information is obtained by interpolating the contributing cell information on the rotor grid. In order to avoid the overlap zone appearing in the place where the flow field changes sharply, the artificial inner boundary is far away from the rotor surface. Considering the characteristics of WENO interpolation, in order to effectively transfer the flow field information between rotor grid and background gird, we use three-layer artificial inner boundaries and artificial outer boundaries, so it is necessary to search three-layer contributing cells. The flow field information exchange between grids is realized by the trilinear interpolation method.
(a)
(b)
(c)
(2) Periodic Boundary Conditions. For the background grid, Figure 16 shows a schematic diagram of the periodic boundary. Considering the symmetry of the flow field, only flow around one blade is computed. The periodic boundary cells are rotationally symmetric, so the calculation method of realizing periodic boundary conditions is as follows:where subscripts 1 and 2 represent the positions of symmetric cells on the periodic boundaries.
(a)
(b)
(3) Modified Far-Field Boundary Conditions. Characteristic-based freestream boundary conditions are enforced at far-field boundaries. However, because the wake vortex system will disappear far below the plane of rotation, the far-field velocity is modified as follows according to the momentum theory [34]. In the circle with radius on the lower surface of the background grid, the velocity which direction is vertical down, and in the other part of the far-field boundaries, which direction points to the center of the rotational plane.
4.7.2. Computational Results
Figure 17 shows the pressure coefficients on the blade surface at four spanwise locations for the solution obtained by the WENO-J, WENO-Z, and WENO-NZ schemes and experimental data. The numerical solutions of the three schemes agree well with the experimental data. Meanwhile, it can be observed from Figure 17(c) that the WENO-NZ scheme captures shock better than the other two WENO schemes.
(a)
(b)
(c)
(d)
Figure 18 illustrates the pressure contours and streamlines with three WENO schemes on the blade surface near the blade tip. Flow separation and reattachment occur on the upper surface, and the separation region calculated by three WENO schemes is right in theory.
(a)
(b)
(c)
(d)
(e)
(f)
Figure 19 shows the vorticity contours on two axial sections. As can be observed from the predicted trajectories of the blade tip vortex for a wake age up to 630°, the computing accuracy on calculating tip vortices is significantly improved with fifth-order WENO and HLLC solver. Additionally, compared with WENO-JS and WENO-Z schemes, WENO-NZ captures vortex wake with a much higher resolution.
(a)
(b)
(c)
(d)
(e)
(f)
5. Conclusion
This paper presents a high-precision robust Riemann solver that can simulate compressible viscous flows of helicopter rotors successfully. In order to capture the complex flows with shock, discontinuity, and vortex, inviscid fluxes are discretized by HLLC approximate Riemann solver, which does not need entropy correction, and an improved version of fifth-order WENO scheme by devising smoothness indicators based on norm measurement is proposed to reconstruct the states on the Riemann interface. Based on the third-order TVD Runge-Kutta scheme for temporal discretization, the proposed method is verified by computing a variety of classic cases, including one-dimensional scalar convection problem, one-dimensional and two-dimensional inviscid compressible flows.
The results shows that the improved WENO scheme has higher accuracy and higher capture capability for discontinuity, shock and vortex than conventional WENO-JS and WENO-Z schemes. Based on Navier-Stokes equations, the numerical properties of the WENO-HLLC solver in conjunction with the implicit LU-SGS time integration method that can accelerate convergence are further validated by simulating steady transonic compressible viscous flows over RAE2822 airfoil and ONERA-M6 wing. Finally, the transonic vortex flow around the helicopter rotor in hover has been studied with the proposed method based on a domain discretized by overset grids. For the pressure coefficients on the blade surface, the numerical solutions are in good agreement with experimental data. The computing accuracy on calculating tip vortices is significantly improved with the fifth-order WENO and HLLC solver, and along with comparisons to WENO-JS and WENO-Z schemes, the proposed WENO scheme can capture vortical wake with higher resolution. It should be pointed out that although the number of grids is not too large, the long wake age up to is still clearly captured with the proposed method.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge the support of the National Natural Science Foundation of China (No. 11502141).