# A Modified Fifth-Order WENO-HLLC Riemann Solver for Transonic Viscous Flows around Helicopter Rotor in Hover

Authors: Li Xu; Fengfeng Zhao; Jingjing Cai
Journal: Shock and Vibration
 (2022)

---

## 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 onL1 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.


---

## Body

## 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 onL2 norm measurement [14]. The weighting relationship of the fifth-order WENO-JS scheme is ωk−ck=OΔx2. 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 ωk−ck=OΔx3 [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 usingL1 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 theL1 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 ωk−ck=OΔx3.The organization of the article is as follows. Section2 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:(1)ddt∭VWdV+∯∂VH−Hv⋅ndS+∭VGdV=0,where V is a control volume with ∂V as the boundary, and n is the unit outward normal of ∂V. W is the conserved variable vector, H is the convective flux vector, Hv is the viscous flux vector, and G caused by rotation can be treated as a source term, given as follows:(2)W=ρρuρvρwρE,H=ρq−qbρuq−qb+pIxρvq−qb+pIyρwq−qb+pIzρHq−qb+pqb,Hv=0τxxIx+τxyIy+τxzIzτxyIx+τyyIy+τyzIzτxzIx+τyzIy+τzzIzf5Ix+g5Iy+h5Iz,G=0ρΩ×qxρΩ×qyρΩ×qz0.ρ,p,E,H and q=uIx+vIy+wIz are the density, pressure, total energy, total enthalpy, and velocity vector of the fluid, respectively. qb is the velocity vector of the grid. Ω denotes the angular velocity vector of the rotating frame. τi,j,i,j=x,y,z are the mean viscous stresses. f5, g5 and h5 are given by the following:(3)f5=uτxx+vτxy+wτxz+k∂T∂x,g5=uτxy+vτyy+wτyz+k∂T∂y,h5=uτxz+vτyz+wτzz+k∂T∂z.Here,T and k 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 p∞, free stream temperature T∞ , and mean aerodynamic chord length l 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 M∞γ/Re∞.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) i,j,k with a structured grid in space, resulting in a system of ordinary differential equations in time:(4)ddtVi,j,kWi,j,k=RHSi,j,k.With residuals(5)RHSi,j,k=−∑fH−Hv⋅nΔSf+Vi,j,kGi,j,k.Wi,j,k and Gi,j,k are the volume-averaged of variables W and G, respectively. Subscript f is taken as i±1/2,j,k,i,j±1/2,k,i,j,k±1/2 which denote the indexes of the control surfaces and ΔS 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):(6)L+DD−1U+DΔWi,j,k=RHSi,j,k.Here,ΔW=Wn+1−Wn. L,D and U are as follows:(7)L=−Ai−1,j,k++Bi,j−1,k++Ci,j,k−1+,D=Vi,j,kΔt+σλA+λB+λC+γμρPrξ2+η2+ζ2I,U=Ai+1,j,k−+Bi,j+1,k−+Ci,j,k+1−,with(8)A±=12A±σλA+γμρPrξ2I,B±=12B±σλB+γμρPrη2I,C±=12C±σλC+γμρPrζ2I.A=∂H⋅S/∂W|ξ, B=∂H⋅S/∂W|η, C=∂H⋅S/∂W|ζ are Jacobian matrices on outward normal ξ, η, ζ of interface,Respectively.λA, λB and λC are their spectral radiuses, and for example,λAi,j,k=1/2q−qbi,j,k⋅Si−1/2,j,k+Si+1/2,j,k+ai,j,k=1/2Si−1/2,j,k+Si+1/2,j,k, where ai,j,k 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.Figure1 shows the structure of the HLLC Riemann solver. In the figure, WL and WR denote the initial left and right states. WL∗ and WR∗ denote the left and right states to the contact discontinuity. SL and SR are the fastest signal velocities, and S∗ 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.Figure 1
The structure of the HLLC Riemann solver.The HLLC numerical fluxesH^=H⋅n in the interface i+1/2,j,k are defined as follows:(9)H^i+1/2HLLC=H^L,SL>0,H^L∗,SL≤0<S∗,H^R∗,S∗≤0<SR,H^R,SR≤0,with the intermediate fluxes H^L∗ and H^R∗ still to be determined. The Rankine-Hugoniot conditions are imposed on each of the waves of velocity SL, S∗ and SR, resulting in the following:(10)H^L∗=H^L+SLWL∗−WL,(11)H^R∗=H^L∗+S∗WR∗−WL∗,(12)H^R∗=H^R+SRWR∗−WR.These three equations involve four unknown vectorsWL∗, WR∗, H^L∗, H^R∗. It can be seen from (10-12) that finding solutions for the two intermediate state vectors WL∗ and WR∗ 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:(13)pL∗=pR∗=p∗,qL∗=qR∗=q∗,where q=unx+vny+wnz.(14)S∗=qv∗−qb,with qb=ubnx+vbny+wbnz. From expressions (10), (12), (14), the pressure in the “∗” domain can be obtained as follows:(15)pL∗=pL+ρLSL−qL+qbS∗−qL+qb,pR∗=pR+ρRSR−qR+qbS∗−qR+qb.SincepL∗=pR∗, we have the following:(16)S∗=pR−pL+ρLSL−qL+qbqL−qb−ρRSR−qR+qbqR−qbρLSL−qL+qb−ρRSR−qR+qb.Algebraic manipulation of (10 and 12) and applying the values pL∗ and pR∗ from (15) gives the intermediate states as follows:(17)WK∗=ρK∗ρuK∗ρvK∗ρwK∗ρEK∗=ρKSK−qK−qbSK−S∗1ρuKρvKρwKρEK+1SK−S∗0p∗−pKnxp∗−pKnyp∗−pKnzp∗S∗−qb−pKqK,for K=L,R. The corresponding intermediate fluxes can be written as follows:(18)H^K∗=H^K+SKWK∗−WK.The left and right wave speeds can be simply estimated by(19)SL=minqL−qb−cL,q˜−qb−c˜,SR=maxqR−qb+cR,q˜−qb+c˜,where c 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 cellsIi=xi−1/2,xi+1/2,i=1,2,… with a uniform interval length Δx=Δxi=xi+1/2−xi−1/2. Considering that the right state can be obtained by symmetry, we only discuss the high-order reconstruction of the left state Wi+1/2L at the cell interface xi+1/2. u represents an arbitrary component of the vector W. The fifth-order WENO reconstruction is carried out on three stencils:(20)Ski=Ii+k−2,Ii+k−1,Ii+k,k=0,1,2.By interpolating the three nodesxi+k−2,xi+k−1,xi+k,k=0,1,2 on each stencil, three different reconstructions for the left state values ui+1/2L are produced as follows:(21)ui+1/2L,0=13ui−2−76ui−1+116ui,ui+1/2L,1=−16ui−1+56ui+13ui+1,ui+1/2L,2=13ui+56ui+1−16ui+2,which satisfy(22)ui+1/2L,k=ui+1/2+akΔx3+OΔx4,k=0,1,2.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 Si=S0i,S1i,S2i=Ii−2,Ii−1,Ii,Ii+1,Ii+2, the fifth-order central upwind interpolation scheme can be obtained as follows:(23)ui+1/2=∑k=02ckui+1/2L,k+OΔx5,with c0=0.1,c1=0.6,c2=0.3.Then,(24)∑k=02ckui+1/2L,k=ui+1/2+bΔx5+OΔx6.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:(25)ui+1/2L=∑k=02ωkui+1/2L,k,where weighting factors ωk,k=0,1,2 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]:(26)ωk=αk∑k=02αk.with(27)αk=ckβk+εp.The introduction of0<ε≪1 is to avoid the denominator being zero. Smoothness indicator βk, a function of ui+k−2,ui+k−1,ui+k, is used to measure the smoothness of the stencil Ski. WENO flux (25) can be expressed as follows:(28)ui+1/2L=∑k=02ckui+1/2L,k+∑k=02ωk−ckui+1/2L,k.Letω+ and ω− denote the weighting factors of ui+1/2L and ui−1/2L, respectively. Substituting (22) and (24) into (28), we obtain the following:(29)ui+1/2L=ui+1/2+b+Δx5+OΔx6+∑k=02ui+1/2+αk+Δx3+OΔx4ωk+−ck.In the same way,(30)ui−1/2L=ui−1/2+b−Δx5+OΔx6+∑k=02ui−1/2+αk−Δx3+OΔx4ωk−−ck.For the fifth-order WENO scheme, the expression(31)ui+1/2L−ui−1/2LΔx=ui+1/2−ui−1/2Δx+OΔx5,should be true, and then the sufficient condition is as follows:(32)ωk±−ck=OΔx3,k=0,1,2.


#### 3.3.2. The Construction of the Improved Fifth-Order WENO Scheme

In [18], Ha et al. introduced a novel smoothness indicator by defining Ln,kuxi+1/2 to be as follows:(33)Ln,kuxi+1/2≔∑xl∈Skick,lnuxl,n=1,2,where the coefficients ck,ln:xl∈Ski,n=1,2 are obtained by solving the system(34)∑xl∈Skick,lnxl−xi+1/2mm!=Δxn,n=m,0,n≠m,m=0,1,2.Then,Ln,ku:=Ln,kuxi+1/2,n=1,2 are expressed as follows:(35)L1,ku=1−kui−2+k+2k−3ui−1+k+2−kui+k,L2,ku=ui−2+k−2ui−1+k+ui+k.Ln,ku follows that(36)Ln,ku=dnudxnxi+1/2Δxn+OΔx3,n=1,2.Based on the termsLn,ku:=Ln,kuxi+1/2n=1,2, one kind of smoothness indicator with L1 norm instead of L2 norm for the WENO-JS scheme was given by the following [18]:(37)βk=ξL1,ku+L2,ku,ξ∈0,1,k=0,1,2.By Taylor expanding for smoothness indicatorsβk, it follows that(38)β0=ξu′Δx−2324u′′′Δx3i+1/2+u″Δx2−32u′′′Δx3i+1/2+OΔx4,β1=ξu′Δx+124u′′′Δx3i+1/2+u″Δx2−12u′′′Δx3i+1/2+OΔx4,β2=ξu′Δx+124u′′′Δx3i+1/2+u″Δx2+12u′′′Δx3i+1/2+OΔx4.Thus,(39)βk=AΔx1+OΔx,ωk−ck=OΔx.The nonlinear weighting factors were given by the following:(40)ωkNS=αkNS∑l=02αlNS,αkNS=ck1+ςβk+ε2,where(41)ς=12β0−β22+gL1,1u2,gx=x31+x3.This fifth-order WENO scheme is WENO-NS [18].Because three stencils provide an unbalanced contribution to the evaluation pointxi+1/2 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 Ln,ku follow that(42)L1,0u=u′Δx−2324u′′′Δx3+u4Δx4−11991920u5Δx5+OΔx6,L1,1u=u′Δx+124u′′′Δx3+11920u5Δx5+OΔx6,L1,2u=u′Δx+124u′′′Δx3+11920u5Δx5+OΔx6.And(43)L2,0u=u″Δx2−32u′′′Δx3+2924u4Δx4−1116u5Δx5+OΔx6,L2,1u=u″Δx2−12u′′′Δx3+524u4Δx4−116u5Δx5+OΔx6,L2,2u=u″Δx2+12u′′′Δx3+524u4Δx4+116u5Δx5+OΔx6.Thus, we construct local smoothness indicators(44)βk=ξL1,ku2+L2,ku12,ξ∈0,1,k=0,1,2,so as to further reduce the weight of the stencil containing discontinuities. Their Taylor expansions can be written as follows:(45)β0=ξu′2Δx2+−2312ξu′u′′′+u″2Δx4+2ξu′u4−3u″u′′′Δx5+OΔx6,β1=ξu′2Δx2+112ξu′u′′′+u″2Δx4−u″u′′′Δx5+OΔx6,β2=ξu′2Δx2+112ξu′u′′′+u″2Δx4+u″u′′′Δx5+OΔx6.To ensure that three stencils provide balanced contribution to the evaluation pointxi+1/2, the smoothness indicators are further modified to the following:(46)β^0=β0,β^1=1+ηβ1,β^2=1−ηβ2,where the parameter η∈0,1 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:(47)γ=L2,2u−L2,0u−2L1,2u−L1,0u2.Its Taylor expansion can be written by the following:(48)γ=u4Δx4+OΔx52=u42Δx8+OΔx9.Then the new smoothness indicators are defined as follows:(49)β˜k=β^k+εβ^k+ε+γ,k=0,1,2.Furthermore, the new nonlinear weighting factors are given by the following:(50)ωk=αk∑l=02αk,αk=ckβ˜k,=ck1+γβ^k+ε.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ωk−ck=OΔx3.First, consider the convergence order when there are no extreme points. Letε=0, and then the local smoothness indicators can be written by the following:(51)β^k=AΔx21+OΔx2,with A=ξui+1/2′2, 1+ηξui+1/2′2 or 1−ηξui+1/2′2 from (45 and 46). Let B=u42 for (48), and then γ=BΔx81+OΔx. Thus,(52)αk=ck1+γβ^k,=ck1+BAΔx61+OΔx.From (50) and ∑k=02ck=1, the relation between the new weights ωk and the ideal weights ck is as follows:(53)ωk−ck=OΔx6,k=0,1,2,which satisfies the condition ωk−ck=OΔx3.Second, for the convergence order at an extreme point whereui+1/2′2=0 and ui+1/2″2≠0, the local smoothness indicators can be written by the following:(54)β^k=AΔx41+OΔx,with A=ui+1/2″2, 1+ηui+1/2″2 or 1−ηui+1/2″2. Similarly, it can be obtained that(55)αk=ck1+BAΔx41+OΔx.And then(56)ωk−ck=OΔx4,k=0,1,2,which also satisfies the condition ωk−ck=OΔx3.


## 3.1. The Finite Volume Method

The finite volume method is used to discretize (1) over each control volume (grid cell) i,j,k with a structured grid in space, resulting in a system of ordinary differential equations in time:(4)ddtVi,j,kWi,j,k=RHSi,j,k.With residuals(5)RHSi,j,k=−∑fH−Hv⋅nΔSf+Vi,j,kGi,j,k.Wi,j,k and Gi,j,k are the volume-averaged of variables W and G, respectively. Subscript f is taken as i±1/2,j,k,i,j±1/2,k,i,j,k±1/2 which denote the indexes of the control surfaces and ΔS 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):(6)L+DD−1U+DΔWi,j,k=RHSi,j,k.Here,ΔW=Wn+1−Wn. L,D and U are as follows:(7)L=−Ai−1,j,k++Bi,j−1,k++Ci,j,k−1+,D=Vi,j,kΔt+σλA+λB+λC+γμρPrξ2+η2+ζ2I,U=Ai+1,j,k−+Bi,j+1,k−+Ci,j,k+1−,with(8)A±=12A±σλA+γμρPrξ2I,B±=12B±σλB+γμρPrη2I,C±=12C±σλC+γμρPrζ2I.A=∂H⋅S/∂W|ξ, B=∂H⋅S/∂W|η, C=∂H⋅S/∂W|ζ are Jacobian matrices on outward normal ξ, η, ζ of interface,Respectively.λA, λB and λC are their spectral radiuses, and for example,λAi,j,k=1/2q−qbi,j,k⋅Si−1/2,j,k+Si+1/2,j,k+ai,j,k=1/2Si−1/2,j,k+Si+1/2,j,k, where ai,j,k 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.Figure1 shows the structure of the HLLC Riemann solver. In the figure, WL and WR denote the initial left and right states. WL∗ and WR∗ denote the left and right states to the contact discontinuity. SL and SR are the fastest signal velocities, and S∗ 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.Figure 1
The structure of the HLLC Riemann solver.The HLLC numerical fluxesH^=H⋅n in the interface i+1/2,j,k are defined as follows:(9)H^i+1/2HLLC=H^L,SL>0,H^L∗,SL≤0<S∗,H^R∗,S∗≤0<SR,H^R,SR≤0,with the intermediate fluxes H^L∗ and H^R∗ still to be determined. The Rankine-Hugoniot conditions are imposed on each of the waves of velocity SL, S∗ and SR, resulting in the following:(10)H^L∗=H^L+SLWL∗−WL,(11)H^R∗=H^L∗+S∗WR∗−WL∗,(12)H^R∗=H^R+SRWR∗−WR.These three equations involve four unknown vectorsWL∗, WR∗, H^L∗, H^R∗. It can be seen from (10-12) that finding solutions for the two intermediate state vectors WL∗ and WR∗ 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:(13)pL∗=pR∗=p∗,qL∗=qR∗=q∗,where q=unx+vny+wnz.(14)S∗=qv∗−qb,with qb=ubnx+vbny+wbnz. From expressions (10), (12), (14), the pressure in the “∗” domain can be obtained as follows:(15)pL∗=pL+ρLSL−qL+qbS∗−qL+qb,pR∗=pR+ρRSR−qR+qbS∗−qR+qb.SincepL∗=pR∗, we have the following:(16)S∗=pR−pL+ρLSL−qL+qbqL−qb−ρRSR−qR+qbqR−qbρLSL−qL+qb−ρRSR−qR+qb.Algebraic manipulation of (10 and 12) and applying the values pL∗ and pR∗ from (15) gives the intermediate states as follows:(17)WK∗=ρK∗ρuK∗ρvK∗ρwK∗ρEK∗=ρKSK−qK−qbSK−S∗1ρuKρvKρwKρEK+1SK−S∗0p∗−pKnxp∗−pKnyp∗−pKnzp∗S∗−qb−pKqK,for K=L,R. The corresponding intermediate fluxes can be written as follows:(18)H^K∗=H^K+SKWK∗−WK.The left and right wave speeds can be simply estimated by(19)SL=minqL−qb−cL,q˜−qb−c˜,SR=maxqR−qb+cR,q˜−qb+c˜,where c 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 cellsIi=xi−1/2,xi+1/2,i=1,2,… with a uniform interval length Δx=Δxi=xi+1/2−xi−1/2. Considering that the right state can be obtained by symmetry, we only discuss the high-order reconstruction of the left state Wi+1/2L at the cell interface xi+1/2. u represents an arbitrary component of the vector W. The fifth-order WENO reconstruction is carried out on three stencils:(20)Ski=Ii+k−2,Ii+k−1,Ii+k,k=0,1,2.By interpolating the three nodesxi+k−2,xi+k−1,xi+k,k=0,1,2 on each stencil, three different reconstructions for the left state values ui+1/2L are produced as follows:(21)ui+1/2L,0=13ui−2−76ui−1+116ui,ui+1/2L,1=−16ui−1+56ui+13ui+1,ui+1/2L,2=13ui+56ui+1−16ui+2,which satisfy(22)ui+1/2L,k=ui+1/2+akΔx3+OΔx4,k=0,1,2.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 Si=S0i,S1i,S2i=Ii−2,Ii−1,Ii,Ii+1,Ii+2, the fifth-order central upwind interpolation scheme can be obtained as follows:(23)ui+1/2=∑k=02ckui+1/2L,k+OΔx5,with c0=0.1,c1=0.6,c2=0.3.Then,(24)∑k=02ckui+1/2L,k=ui+1/2+bΔx5+OΔx6.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:(25)ui+1/2L=∑k=02ωkui+1/2L,k,where weighting factors ωk,k=0,1,2 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]:(26)ωk=αk∑k=02αk.with(27)αk=ckβk+εp.The introduction of0<ε≪1 is to avoid the denominator being zero. Smoothness indicator βk, a function of ui+k−2,ui+k−1,ui+k, is used to measure the smoothness of the stencil Ski. WENO flux (25) can be expressed as follows:(28)ui+1/2L=∑k=02ckui+1/2L,k+∑k=02ωk−ckui+1/2L,k.Letω+ and ω− denote the weighting factors of ui+1/2L and ui−1/2L, respectively. Substituting (22) and (24) into (28), we obtain the following:(29)ui+1/2L=ui+1/2+b+Δx5+OΔx6+∑k=02ui+1/2+αk+Δx3+OΔx4ωk+−ck.In the same way,(30)ui−1/2L=ui−1/2+b−Δx5+OΔx6+∑k=02ui−1/2+αk−Δx3+OΔx4ωk−−ck.For the fifth-order WENO scheme, the expression(31)ui+1/2L−ui−1/2LΔx=ui+1/2−ui−1/2Δx+OΔx5,should be true, and then the sufficient condition is as follows:(32)ωk±−ck=OΔx3,k=0,1,2.


### 3.3.2. The Construction of the Improved Fifth-Order WENO Scheme

In [18], Ha et al. introduced a novel smoothness indicator by defining Ln,kuxi+1/2 to be as follows:(33)Ln,kuxi+1/2≔∑xl∈Skick,lnuxl,n=1,2,where the coefficients ck,ln:xl∈Ski,n=1,2 are obtained by solving the system(34)∑xl∈Skick,lnxl−xi+1/2mm!=Δxn,n=m,0,n≠m,m=0,1,2.Then,Ln,ku:=Ln,kuxi+1/2,n=1,2 are expressed as follows:(35)L1,ku=1−kui−2+k+2k−3ui−1+k+2−kui+k,L2,ku=ui−2+k−2ui−1+k+ui+k.Ln,ku follows that(36)Ln,ku=dnudxnxi+1/2Δxn+OΔx3,n=1,2.Based on the termsLn,ku:=Ln,kuxi+1/2n=1,2, one kind of smoothness indicator with L1 norm instead of L2 norm for the WENO-JS scheme was given by the following [18]:(37)βk=ξL1,ku+L2,ku,ξ∈0,1,k=0,1,2.By Taylor expanding for smoothness indicatorsβk, it follows that(38)β0=ξu′Δx−2324u′′′Δx3i+1/2+u″Δx2−32u′′′Δx3i+1/2+OΔx4,β1=ξu′Δx+124u′′′Δx3i+1/2+u″Δx2−12u′′′Δx3i+1/2+OΔx4,β2=ξu′Δx+124u′′′Δx3i+1/2+u″Δx2+12u′′′Δx3i+1/2+OΔx4.Thus,(39)βk=AΔx1+OΔx,ωk−ck=OΔx.The nonlinear weighting factors were given by the following:(40)ωkNS=αkNS∑l=02αlNS,αkNS=ck1+ςβk+ε2,where(41)ς=12β0−β22+gL1,1u2,gx=x31+x3.This fifth-order WENO scheme is WENO-NS [18].Because three stencils provide an unbalanced contribution to the evaluation pointxi+1/2 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 Ln,ku follow that(42)L1,0u=u′Δx−2324u′′′Δx3+u4Δx4−11991920u5Δx5+OΔx6,L1,1u=u′Δx+124u′′′Δx3+11920u5Δx5+OΔx6,L1,2u=u′Δx+124u′′′Δx3+11920u5Δx5+OΔx6.And(43)L2,0u=u″Δx2−32u′′′Δx3+2924u4Δx4−1116u5Δx5+OΔx6,L2,1u=u″Δx2−12u′′′Δx3+524u4Δx4−116u5Δx5+OΔx6,L2,2u=u″Δx2+12u′′′Δx3+524u4Δx4+116u5Δx5+OΔx6.Thus, we construct local smoothness indicators(44)βk=ξL1,ku2+L2,ku12,ξ∈0,1,k=0,1,2,so as to further reduce the weight of the stencil containing discontinuities. Their Taylor expansions can be written as follows:(45)β0=ξu′2Δx2+−2312ξu′u′′′+u″2Δx4+2ξu′u4−3u″u′′′Δx5+OΔx6,β1=ξu′2Δx2+112ξu′u′′′+u″2Δx4−u″u′′′Δx5+OΔx6,β2=ξu′2Δx2+112ξu′u′′′+u″2Δx4+u″u′′′Δx5+OΔx6.To ensure that three stencils provide balanced contribution to the evaluation pointxi+1/2, the smoothness indicators are further modified to the following:(46)β^0=β0,β^1=1+ηβ1,β^2=1−ηβ2,where the parameter η∈0,1 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:(47)γ=L2,2u−L2,0u−2L1,2u−L1,0u2.Its Taylor expansion can be written by the following:(48)γ=u4Δx4+OΔx52=u42Δx8+OΔx9.Then the new smoothness indicators are defined as follows:(49)β˜k=β^k+εβ^k+ε+γ,k=0,1,2.Furthermore, the new nonlinear weighting factors are given by the following:(50)ωk=αk∑l=02αk,αk=ckβ˜k,=ck1+γβ^k+ε.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ωk−ck=OΔx3.First, consider the convergence order when there are no extreme points. Letε=0, and then the local smoothness indicators can be written by the following:(51)β^k=AΔx21+OΔx2,with A=ξui+1/2′2, 1+ηξui+1/2′2 or 1−ηξui+1/2′2 from (45 and 46). Let B=u42 for (48), and then γ=BΔx81+OΔx. Thus,(52)αk=ck1+γβ^k,=ck1+BAΔx61+OΔx.From (50) and ∑k=02ck=1, the relation between the new weights ωk and the ideal weights ck is as follows:(53)ωk−ck=OΔx6,k=0,1,2,which satisfies the condition ωk−ck=OΔx3.Second, for the convergence order at an extreme point whereui+1/2′2=0 and ui+1/2″2≠0, the local smoothness indicators can be written by the following:(54)β^k=AΔx41+OΔx,with A=ui+1/2″2, 1+ηui+1/2″2 or 1−ηui+1/2″2. Similarly, it can be obtained that(55)αk=ck1+BAΔx41+OΔx.And then(56)ωk−ck=OΔx4,k=0,1,2,which also satisfies the condition ωk−ck=OΔx3.


## 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 cellsIi=xi−1/2,xi+1/2,i=1,2,… with a uniform interval length Δx=Δxi=xi+1/2−xi−1/2. Considering that the right state can be obtained by symmetry, we only discuss the high-order reconstruction of the left state Wi+1/2L at the cell interface xi+1/2. u represents an arbitrary component of the vector W. The fifth-order WENO reconstruction is carried out on three stencils:(20)Ski=Ii+k−2,Ii+k−1,Ii+k,k=0,1,2.By interpolating the three nodesxi+k−2,xi+k−1,xi+k,k=0,1,2 on each stencil, three different reconstructions for the left state values ui+1/2L are produced as follows:(21)ui+1/2L,0=13ui−2−76ui−1+116ui,ui+1/2L,1=−16ui−1+56ui+13ui+1,ui+1/2L,2=13ui+56ui+1−16ui+2,which satisfy(22)ui+1/2L,k=ui+1/2+akΔx3+OΔx4,k=0,1,2.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 Si=S0i,S1i,S2i=Ii−2,Ii−1,Ii,Ii+1,Ii+2, the fifth-order central upwind interpolation scheme can be obtained as follows:(23)ui+1/2=∑k=02ckui+1/2L,k+OΔx5,with c0=0.1,c1=0.6,c2=0.3.Then,(24)∑k=02ckui+1/2L,k=ui+1/2+bΔx5+OΔx6.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:(25)ui+1/2L=∑k=02ωkui+1/2L,k,where weighting factors ωk,k=0,1,2 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]:(26)ωk=αk∑k=02αk.with(27)αk=ckβk+εp.The introduction of0<ε≪1 is to avoid the denominator being zero. Smoothness indicator βk, a function of ui+k−2,ui+k−1,ui+k, is used to measure the smoothness of the stencil Ski. WENO flux (25) can be expressed as follows:(28)ui+1/2L=∑k=02ckui+1/2L,k+∑k=02ωk−ckui+1/2L,k.Letω+ and ω− denote the weighting factors of ui+1/2L and ui−1/2L, respectively. Substituting (22) and (24) into (28), we obtain the following:(29)ui+1/2L=ui+1/2+b+Δx5+OΔx6+∑k=02ui+1/2+αk+Δx3+OΔx4ωk+−ck.In the same way,(30)ui−1/2L=ui−1/2+b−Δx5+OΔx6+∑k=02ui−1/2+αk−Δx3+OΔx4ωk−−ck.For the fifth-order WENO scheme, the expression(31)ui+1/2L−ui−1/2LΔx=ui+1/2−ui−1/2Δx+OΔx5,should be true, and then the sufficient condition is as follows:(32)ωk±−ck=OΔx3,k=0,1,2.


## 3.3.2. The Construction of the Improved Fifth-Order WENO Scheme

In [18], Ha et al. introduced a novel smoothness indicator by defining Ln,kuxi+1/2 to be as follows:(33)Ln,kuxi+1/2≔∑xl∈Skick,lnuxl,n=1,2,where the coefficients ck,ln:xl∈Ski,n=1,2 are obtained by solving the system(34)∑xl∈Skick,lnxl−xi+1/2mm!=Δxn,n=m,0,n≠m,m=0,1,2.Then,Ln,ku:=Ln,kuxi+1/2,n=1,2 are expressed as follows:(35)L1,ku=1−kui−2+k+2k−3ui−1+k+2−kui+k,L2,ku=ui−2+k−2ui−1+k+ui+k.Ln,ku follows that(36)Ln,ku=dnudxnxi+1/2Δxn+OΔx3,n=1,2.Based on the termsLn,ku:=Ln,kuxi+1/2n=1,2, one kind of smoothness indicator with L1 norm instead of L2 norm for the WENO-JS scheme was given by the following [18]:(37)βk=ξL1,ku+L2,ku,ξ∈0,1,k=0,1,2.By Taylor expanding for smoothness indicatorsβk, it follows that(38)β0=ξu′Δx−2324u′′′Δx3i+1/2+u″Δx2−32u′′′Δx3i+1/2+OΔx4,β1=ξu′Δx+124u′′′Δx3i+1/2+u″Δx2−12u′′′Δx3i+1/2+OΔx4,β2=ξu′Δx+124u′′′Δx3i+1/2+u″Δx2+12u′′′Δx3i+1/2+OΔx4.Thus,(39)βk=AΔx1+OΔx,ωk−ck=OΔx.The nonlinear weighting factors were given by the following:(40)ωkNS=αkNS∑l=02αlNS,αkNS=ck1+ςβk+ε2,where(41)ς=12β0−β22+gL1,1u2,gx=x31+x3.This fifth-order WENO scheme is WENO-NS [18].Because three stencils provide an unbalanced contribution to the evaluation pointxi+1/2 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 Ln,ku follow that(42)L1,0u=u′Δx−2324u′′′Δx3+u4Δx4−11991920u5Δx5+OΔx6,L1,1u=u′Δx+124u′′′Δx3+11920u5Δx5+OΔx6,L1,2u=u′Δx+124u′′′Δx3+11920u5Δx5+OΔx6.And(43)L2,0u=u″Δx2−32u′′′Δx3+2924u4Δx4−1116u5Δx5+OΔx6,L2,1u=u″Δx2−12u′′′Δx3+524u4Δx4−116u5Δx5+OΔx6,L2,2u=u″Δx2+12u′′′Δx3+524u4Δx4+116u5Δx5+OΔx6.Thus, we construct local smoothness indicators(44)βk=ξL1,ku2+L2,ku12,ξ∈0,1,k=0,1,2,so as to further reduce the weight of the stencil containing discontinuities. Their Taylor expansions can be written as follows:(45)β0=ξu′2Δx2+−2312ξu′u′′′+u″2Δx4+2ξu′u4−3u″u′′′Δx5+OΔx6,β1=ξu′2Δx2+112ξu′u′′′+u″2Δx4−u″u′′′Δx5+OΔx6,β2=ξu′2Δx2+112ξu′u′′′+u″2Δx4+u″u′′′Δx5+OΔx6.To ensure that three stencils provide balanced contribution to the evaluation pointxi+1/2, the smoothness indicators are further modified to the following:(46)β^0=β0,β^1=1+ηβ1,β^2=1−ηβ2,where the parameter η∈0,1 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:(47)γ=L2,2u−L2,0u−2L1,2u−L1,0u2.Its Taylor expansion can be written by the following:(48)γ=u4Δx4+OΔx52=u42Δx8+OΔx9.Then the new smoothness indicators are defined as follows:(49)β˜k=β^k+εβ^k+ε+γ,k=0,1,2.Furthermore, the new nonlinear weighting factors are given by the following:(50)ωk=αk∑l=02αk,αk=ckβ˜k,=ck1+γβ^k+ε.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ωk−ck=OΔx3.First, consider the convergence order when there are no extreme points. Letε=0, and then the local smoothness indicators can be written by the following:(51)β^k=AΔx21+OΔx2,with A=ξui+1/2′2, 1+ηξui+1/2′2 or 1−ηξui+1/2′2 from (45 and 46). Let B=u42 for (48), and then γ=BΔx81+OΔx. Thus,(52)αk=ck1+γβ^k,=ck1+BAΔx61+OΔx.From (50) and ∑k=02ck=1, the relation between the new weights ωk and the ideal weights ck is as follows:(53)ωk−ck=OΔx6,k=0,1,2,which satisfies the condition ωk−ck=OΔx3.Second, for the convergence order at an extreme point whereui+1/2′2=0 and ui+1/2″2≠0, the local smoothness indicators can be written by the following:(54)β^k=AΔx41+OΔx,with A=ui+1/2″2, 1+ηui+1/2″2 or 1−ηui+1/2″2. Similarly, it can be obtained that(55)αk=ck1+BAΔx41+OΔx.And then(56)ωk−ck=OΔx4,k=0,1,2,which also satisfies the condition ωk−ck=OΔx3.


## 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 ε=10−6 for WENO-JS and ε=10−40 for WENO-Z and WENO-NZ. We have chosen ξ=0.2 for (45) and η=0.05 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:(57)∂u∂t+∂u∂x=0,−1≤x≤1,t>0,ux,0=ϕx.The numerical result will be compared with the exact solutionux,t=ϕx−t. Errors are measured by L1 and L∞ norms(58)L1=1N∑i=1N|ui−uexact,L∞=maxui−uexact,i=1,2,…,N.Case 1.
To test the accuracy and convergence order of the proposed method, we first consider the initial condition(59)ux,0=sinπx.
The 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Δt=Δx5/3. The periodic boundary condition is adopted and the results at t=2 are discussed.
Tables1 and 2 contain the L∞ and L1 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 L∞ error or L1 error of both WENO-NZ and WENO-Z schemes reach about 10−10 that is smaller than that of WENO-JS scheme, the WENO-NZ scheme approaches the desired convergence order.Table 1
L∞ errors and convergence order on Case 1: 1D scalar convection equation with initial condition (59).
NWENO-JSWENO-ZWENO-NZ202.5305E − 03 (—)3.2296E − 04 (—)3.1546E − 04 (—)408.8367E − 05 (4.84)1.0099E − 05 (5.00)1.0087E − 05 (4.97)802.8331E − 06 (4.96)3.1463E − 07 (5.00)3.1460E − 07 (5.00)1608.6156E − 08 (5.04)9.8084E − 08 (5.00)9.8083E − 09 (5.00)3202.6479E − 09 (5.02)3.0602E − 10 (5.00)3.0599E − 10 (5.00)Table 2
L1 errors and convergence order on Case 1: 1D scalar convection equation with initial condition (59).
NWENO-JSWENO-ZWENO-NZ201.4400E − 03 (—)2.0581E − 04 (—)2.0525E − 04 (—)404.4848E − 05 (5.01)6.4236E − 06 (5.00)6.4235E − 06 (5.00)801.3971E − 06 (5.01)2.0024E − 07 (5.00)2.0024E − 07 (5.00)1604.3612E − 08 (5.00)6.2437E − 09 (5.00)6.2437E − 09 (5.00)3201.3623E − 09 (5.00)1.9485E − 10 (5.00)1.9485E − 10 (5.00)Case 2.
For the one-dimensional scalar convection equation, consider the initial condition(60)ux,0=sinπx−1πsinπx,which 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 L∞ and L1 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.Table 3
L∞ errors and convergence order on Case 2: 1D scalar convection equation with initial condition (60).
NWENO-JSWENO-ZWENO-NZ201.4354E − 02 (—)3.4610E − 03 (—)5.6629E − 03 (—)401.0940E − 03 (3.71)2.2165E − 04 (3.96)2.0057E − 04 (4.82)808.9344E − 05 (3.61)6.5027E − 06 (5.09)6.3539E − 06 (4.98)1607.4868E − 06 (3.58)2.0489E − 07 (4.99)1.9868E − 07 (5.00)3204.0681E − 07 (4.20)6.4129E − 09 (5.00)6.2030E − 09 (5.00)Table 4
L1 errors and convergence order on Case 2: 1D scalar convection equation with initial condition (60).
NWENO-JSWENO-ZWENO-NZ204.9978E − 03 (—)1.5236E − 03 (—)2.0560E − 03 (—)403.5965E − 04 (3.80)6.4738E − 05 (4.56)6.8655E − 05 (4.90)801.6932E − 05 (4.41)2.2272E − 06 (4.86)2.2127E − 06 (4.96)1607.4194E − 07 (4.51)7.1964E − 08 (4.95)7.0557E − 08 (4.97)3202.5412E − 08 (4.87)2.2719E − 09 (4.99)2.2079E − 09 (5.00)Case 3.
This initial valve problem involves a jump discontinuity, which is described by the initial condition(61)ux,0=−sinπx−12x3,−1≤x<0,−sinπx−12x3+1,0≤x≤1.
The number of the uniform grid cells is 200, and the Courant number is set to beCFL=0.5. For the WENO-JS, WENO-Z, and WENO-NZ schemes, their approximate solutions at t=8.0 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.Figure 2
Numerical solution of Case3 at t=8.0 with grids of N=200 (1D scalar convection equation with initial condition (61)) based on WENO-JS, WENO-Z, and WENO-NZ schemes.


### 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(62)ρ,u,pL=1,0,1,ρ,u,pR=0.125,0,0.1.The solution contains a left rarefaction wave, a right traveling contract contact wave, and a right shock wave. The calculation interval0,1 is uniformly divided into 200 cells, and the Courant number is taken as CFL=0.5. Transmissive boundary conditions are taken. Figure 3 shows the density profiles at t=2 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.Figure 3
Density distribution of Sod’s problem att=2 computed by WENO-JS, WENO-Z, and WENO-NZ schemes with grids of N=200.


### 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(63)ρ,u,p=3.857143,2.629369,313,−5≤x<−4,1+0.2sin5x,0,1,−4≤x≤5.The problem describes the interaction of a right moving Mach 3 shock with a sine entropy wave. Figure4 illustrates the density results at t=1.8 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.Figure 4
Numerical solutions of shock entropy wave interaction att=1.8 computed by WENO-JS, WENO-Z, and WENO-NZ schemes with grids of N=200.


### 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 −5,5×−5,5. The initial flow is given by(64)ρ=ρ∞1−γ−12γb2u∞2e1−γ21/γ−1,p=p∞1−γ−12γb2u∞2e1−γ2γ/γ−1,u=u∞−bu∞y−yce1−γ2/2,v=v∞+bu∞x−xce1−γ2/2,,with free stream density, free stream pressure, and free stream velocities in two directions set as ρ∞=1, p∞=1, u∞=0.2, v∞=0, respectively. The distance r from the point x,y to the vortex core xc,yc is expressed as r=x−xc2+y−yc2. b=5 represents the dimensionless vortex strength.Based on periodic boundary conditions, flow fields are simulated by using equispaced Cartesian grid32×32,64×64, and 128×128, respectively. We observe the results after two periods of vortex propagation. Figure 5 gives L2 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 32×32 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 128×128 grids.Figure 5
Bar graph ofL2 error orders for densities about WENO-JS, WENO-Z, and WENO-NZ schemes.Figure 6
Pressure contours for isentropic vortex convection after two periods.
(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 2.510, the freestream Mach number M∞ 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 onL2 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 7
Convergence history for the RAE2822 airfoil.Figure8 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.Figure 8
Pressure coefficients at the airfoil surface about the transonic flow for RAE2822 airfoil. (a) Full view. (b) Locally enlarged view.
(a)(b)Figure 9
Pressure contours and streamlines with WENO-NZ scheme for RAE2822 airfoil.


### 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 Figure10, 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 3.060, the Reynolds number based on the mean aerodynamic chord is 11.72 × 106, and the freestream Mach number M∞ 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.Figure 10
Mesh for ONERA-M6 wing. (a) Full view. (b) Local mesh around the wing.
(a)(b)Figure11 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 11
Convergence history for ONERA-M6 wing.Figure12 shows the pressure coefficients on the wing surface at six spanwise locations, including Z/C=0.20, Z/C=0.44, Z/C=0.65, Z/C=0.80, Z/C=0.90, Z/C=0.95 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 Z/C=0.90 and Z/C=0.95, 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.Figure 12
Pressure coefficients at the wing surface about the transonic flow over the ONERA-M6 wing.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 13
Pressure contours and streamlines with WENO-NZ scheme on the upper surface for ONERA-M6 wing. (a) Pressure contours. (b) Streamlines.
(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 Mtip=0.877, pitch angle is, 80 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.Figure 14
Overset grid sketch for Caradonna-Tung rotor. (a) Rotor grid. (b) Overset grids with rotor grid and background grid.
(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.Figure 15
Schematic diagram of artificial inner boundary on the background grid. (a) Hole domain. (b) A section of holeomain. (c) Hole boundary sketch.
(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:(65)ρ2=ρ1,p2=p1,u2v2w2=−10001000−1u1v1w1,where subscripts 1 and 2 represent the positions of symmetric cells on the periodic boundaries.Figure 16
Schematic diagram of periodic boundaries. (a) Schematic diagram of aerial view for two blades. (b) Periodic boundary.
(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 R/2 on the lower surface of the background grid, the velocity q∞=2MtipCT/2 which direction is vertical down, and in the other part of the far-field boundaries, q∞=1/4MtipR/r2CT/2 which direction points to the center of the rotational plane.


#### 4.7.2. Computational Results

Figure17 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.Figure 17
Pressure coefficients at the blade surface about the transonic flow over the Caradonna-Tung rotor.
(a)(b)(c)(d)Figure18 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.Figure 18
Pressure contours and streamlines with three WENO schemes on the blade surface for the Caradonna-Tung rotor. (a) WENO-JS. (b) WENO-Z. (c) WENO-NZ. (d) WENO-JS. (e) WENO-Z. (f) WENO-NZ.
(a)(b)(c)(d)(e)(f)Figure19 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.Figure 19
Vorticity contours with three WENO schemes for Caradonna-Tung rotor. (a) WENO-JS. (b) WENO-Z. (c) WENO-NZ. (d) WENO-JS. (e) WENO-Z. (f) WENO-NZ.
(a)(b)(c)(d)(e)(f)


## 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:(57)∂u∂t+∂u∂x=0,−1≤x≤1,t>0,ux,0=ϕx.The numerical result will be compared with the exact solutionux,t=ϕx−t. Errors are measured by L1 and L∞ norms(58)L1=1N∑i=1N|ui−uexact,L∞=maxui−uexact,i=1,2,…,N.Case 1.
To test the accuracy and convergence order of the proposed method, we first consider the initial condition(59)ux,0=sinπx.
The 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Δt=Δx5/3. The periodic boundary condition is adopted and the results at t=2 are discussed.
Tables1 and 2 contain the L∞ and L1 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 L∞ error or L1 error of both WENO-NZ and WENO-Z schemes reach about 10−10 that is smaller than that of WENO-JS scheme, the WENO-NZ scheme approaches the desired convergence order.Table 1
L∞ errors and convergence order on Case 1: 1D scalar convection equation with initial condition (59).
NWENO-JSWENO-ZWENO-NZ202.5305E − 03 (—)3.2296E − 04 (—)3.1546E − 04 (—)408.8367E − 05 (4.84)1.0099E − 05 (5.00)1.0087E − 05 (4.97)802.8331E − 06 (4.96)3.1463E − 07 (5.00)3.1460E − 07 (5.00)1608.6156E − 08 (5.04)9.8084E − 08 (5.00)9.8083E − 09 (5.00)3202.6479E − 09 (5.02)3.0602E − 10 (5.00)3.0599E − 10 (5.00)Table 2
L1 errors and convergence order on Case 1: 1D scalar convection equation with initial condition (59).
NWENO-JSWENO-ZWENO-NZ201.4400E − 03 (—)2.0581E − 04 (—)2.0525E − 04 (—)404.4848E − 05 (5.01)6.4236E − 06 (5.00)6.4235E − 06 (5.00)801.3971E − 06 (5.01)2.0024E − 07 (5.00)2.0024E − 07 (5.00)1604.3612E − 08 (5.00)6.2437E − 09 (5.00)6.2437E − 09 (5.00)3201.3623E − 09 (5.00)1.9485E − 10 (5.00)1.9485E − 10 (5.00)Case 2.
For the one-dimensional scalar convection equation, consider the initial condition(60)ux,0=sinπx−1πsinπx,which 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 L∞ and L1 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.Table 3
L∞ errors and convergence order on Case 2: 1D scalar convection equation with initial condition (60).
NWENO-JSWENO-ZWENO-NZ201.4354E − 02 (—)3.4610E − 03 (—)5.6629E − 03 (—)401.0940E − 03 (3.71)2.2165E − 04 (3.96)2.0057E − 04 (4.82)808.9344E − 05 (3.61)6.5027E − 06 (5.09)6.3539E − 06 (4.98)1607.4868E − 06 (3.58)2.0489E − 07 (4.99)1.9868E − 07 (5.00)3204.0681E − 07 (4.20)6.4129E − 09 (5.00)6.2030E − 09 (5.00)Table 4
L1 errors and convergence order on Case 2: 1D scalar convection equation with initial condition (60).
NWENO-JSWENO-ZWENO-NZ204.9978E − 03 (—)1.5236E − 03 (—)2.0560E − 03 (—)403.5965E − 04 (3.80)6.4738E − 05 (4.56)6.8655E − 05 (4.90)801.6932E − 05 (4.41)2.2272E − 06 (4.86)2.2127E − 06 (4.96)1607.4194E − 07 (4.51)7.1964E − 08 (4.95)7.0557E − 08 (4.97)3202.5412E − 08 (4.87)2.2719E − 09 (4.99)2.2079E − 09 (5.00)Case 3.
This initial valve problem involves a jump discontinuity, which is described by the initial condition(61)ux,0=−sinπx−12x3,−1≤x<0,−sinπx−12x3+1,0≤x≤1.
The number of the uniform grid cells is 200, and the Courant number is set to beCFL=0.5. For the WENO-JS, WENO-Z, and WENO-NZ schemes, their approximate solutions at t=8.0 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.Figure 2
Numerical solution of Case3 at t=8.0 with grids of N=200 (1D scalar convection equation with initial condition (61)) based on WENO-JS, WENO-Z, and WENO-NZ schemes.


## 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(62)ρ,u,pL=1,0,1,ρ,u,pR=0.125,0,0.1.The solution contains a left rarefaction wave, a right traveling contract contact wave, and a right shock wave. The calculation interval0,1 is uniformly divided into 200 cells, and the Courant number is taken as CFL=0.5. Transmissive boundary conditions are taken. Figure 3 shows the density profiles at t=2 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.Figure 3
Density distribution of Sod’s problem att=2 computed by WENO-JS, WENO-Z, and WENO-NZ schemes with grids of N=200.


## 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(63)ρ,u,p=3.857143,2.629369,313,−5≤x<−4,1+0.2sin5x,0,1,−4≤x≤5.The problem describes the interaction of a right moving Mach 3 shock with a sine entropy wave. Figure4 illustrates the density results at t=1.8 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.Figure 4
Numerical solutions of shock entropy wave interaction att=1.8 computed by WENO-JS, WENO-Z, and WENO-NZ schemes with grids of N=200.


## 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 −5,5×−5,5. The initial flow is given by(64)ρ=ρ∞1−γ−12γb2u∞2e1−γ21/γ−1,p=p∞1−γ−12γb2u∞2e1−γ2γ/γ−1,u=u∞−bu∞y−yce1−γ2/2,v=v∞+bu∞x−xce1−γ2/2,,with free stream density, free stream pressure, and free stream velocities in two directions set as ρ∞=1, p∞=1, u∞=0.2, v∞=0, respectively. The distance r from the point x,y to the vortex core xc,yc is expressed as r=x−xc2+y−yc2. b=5 represents the dimensionless vortex strength.Based on periodic boundary conditions, flow fields are simulated by using equispaced Cartesian grid32×32,64×64, and 128×128, respectively. We observe the results after two periods of vortex propagation. Figure 5 gives L2 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 32×32 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 128×128 grids.Figure 5
Bar graph ofL2 error orders for densities about WENO-JS, WENO-Z, and WENO-NZ schemes.Figure 6
Pressure contours for isentropic vortex convection after two periods.
(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 2.510, the freestream Mach number M∞ 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 onL2 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 7
Convergence history for the RAE2822 airfoil.Figure8 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.Figure 8
Pressure coefficients at the airfoil surface about the transonic flow for RAE2822 airfoil. (a) Full view. (b) Locally enlarged view.
(a)(b)Figure 9
Pressure contours and streamlines with WENO-NZ scheme for RAE2822 airfoil.


## 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 Figure10, 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 3.060, the Reynolds number based on the mean aerodynamic chord is 11.72 × 106, and the freestream Mach number M∞ 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.Figure 10
Mesh for ONERA-M6 wing. (a) Full view. (b) Local mesh around the wing.
(a)(b)Figure11 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 11
Convergence history for ONERA-M6 wing.Figure12 shows the pressure coefficients on the wing surface at six spanwise locations, including Z/C=0.20, Z/C=0.44, Z/C=0.65, Z/C=0.80, Z/C=0.90, Z/C=0.95 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 Z/C=0.90 and Z/C=0.95, 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.Figure 12
Pressure coefficients at the wing surface about the transonic flow over the ONERA-M6 wing.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 13
Pressure contours and streamlines with WENO-NZ scheme on the upper surface for ONERA-M6 wing. (a) Pressure contours. (b) Streamlines.
(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 Mtip=0.877, pitch angle is, 80 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.Figure 14
Overset grid sketch for Caradonna-Tung rotor. (a) Rotor grid. (b) Overset grids with rotor grid and background grid.
(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.Figure 15
Schematic diagram of artificial inner boundary on the background grid. (a) Hole domain. (b) A section of holeomain. (c) Hole boundary sketch.
(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:(65)ρ2=ρ1,p2=p1,u2v2w2=−10001000−1u1v1w1,where subscripts 1 and 2 represent the positions of symmetric cells on the periodic boundaries.Figure 16
Schematic diagram of periodic boundaries. (a) Schematic diagram of aerial view for two blades. (b) Periodic boundary.
(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 R/2 on the lower surface of the background grid, the velocity q∞=2MtipCT/2 which direction is vertical down, and in the other part of the far-field boundaries, q∞=1/4MtipR/r2CT/2 which direction points to the center of the rotational plane.


### 4.7.2. Computational Results

Figure17 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.Figure 17
Pressure coefficients at the blade surface about the transonic flow over the Caradonna-Tung rotor.
(a)(b)(c)(d)Figure18 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.Figure 18
Pressure contours and streamlines with three WENO schemes on the blade surface for the Caradonna-Tung rotor. (a) WENO-JS. (b) WENO-Z. (c) WENO-NZ. (d) WENO-JS. (e) WENO-Z. (f) WENO-NZ.
(a)(b)(c)(d)(e)(f)Figure19 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.Figure 19
Vorticity contours with three WENO schemes for Caradonna-Tung rotor. (a) WENO-JS. (b) WENO-Z. (c) WENO-NZ. (d) WENO-JS. (e) WENO-Z. (f) WENO-NZ.
(a)(b)(c)(d)(e)(f)


## 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.Figure 15
Schematic diagram of artificial inner boundary on the background grid. (a) Hole domain. (b) A section of holeomain. (c) Hole boundary sketch.
(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:(65)ρ2=ρ1,p2=p1,u2v2w2=−10001000−1u1v1w1,where subscripts 1 and 2 represent the positions of symmetric cells on the periodic boundaries.Figure 16
Schematic diagram of periodic boundaries. (a) Schematic diagram of aerial view for two blades. (b) Periodic boundary.
(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 R/2 on the lower surface of the background grid, the velocity q∞=2MtipCT/2 which direction is vertical down, and in the other part of the far-field boundaries, q∞=1/4MtipR/r2CT/2 which direction points to the center of the rotational plane.


## 4.7.2. Computational Results

Figure17 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.Figure 17
Pressure coefficients at the blade surface about the transonic flow over the Caradonna-Tung rotor.
(a)(b)(c)(d)Figure18 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.Figure 18
Pressure contours and streamlines with three WENO schemes on the blade surface for the Caradonna-Tung rotor. (a) WENO-JS. (b) WENO-Z. (c) WENO-NZ. (d) WENO-JS. (e) WENO-Z. (f) WENO-NZ.
(a)(b)(c)(d)(e)(f)Figure19 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.Figure 19
Vorticity contours with three WENO schemes for Caradonna-Tung rotor. (a) WENO-JS. (b) WENO-Z. (c) WENO-NZ. (d) WENO-JS. (e) WENO-Z. (f) WENO-NZ.
(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 to6300 is still clearly captured with the proposed method.



---

*Source: 2180845-2022-01-19.xml*