# Helicopter Rotor Flow Analysis Using Mapped Chebyshev Pseudospectral Method and Overset Mesh Topology

Authors: Dong Kyun Im; Seongim Choi
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Unsteady helicopter rotor flows are solved by a Chebyshev pseudospectral method with overset mesh topology which employs Chebyshev polynomials for solution approximation and a Chebyshev collocation operator to represent the time derivative term of the unsteady flow governing equations. Spatial derivative terms of the flux Jacobians are discretized implicitly while the Chebyshev spectral derivative term is treated in explicit form. Unlike the Fourier spectral method, collocation points of standard Chebyshev polynomials are not evenly distributed and heavily clustered near the extremities of the time interval, which makes the spectral derivative matrix ill-conditioned and deteriorates the stability and convergence of the flow solution. A conformal mapping of an arcsin function is applied to redistribute those points more evenly and thus to improve the numerical stability of the linear system. A parameter study on the condition number of the spectral derivative matrix with respect to the control parameters of the mapping function is also carried out. For the validation of the proposed method, both periodic and nonperiodic unsteady flow problems were solved with two-dimensional problems: an oscillating airfoil with a fixed frequency and a plunging airfoil with constant plunging speed without considering gravitational force. Computation results of the Chebyshev pseudospectral method showed excellent agreements with those of the time-marching computation. Subsequently, helicopter rotor flows in hovering and nonlifting forward flight are solved. Moving boundaries of the rotating rotor blades are efficiently managed by the overset mesh topology. As a set of subgrids are constructed only one time at the beginning of the solution procedure corresponding to the mapped Chebyshev collocation points, computation time for mesh interpolation of hole-cutting between background and near-body grids becomes drastically reduced when compared to the time-marching computation method where subgrid movement and the hole-cutting need to be carried out at each physical time step. The number of the collocation points was varied to investigate the sensitivity of the solution accuracy, computation time, and memory. Computation results are compared with those from the time-marching computation, the Fourier spectral method, and wind-tunnel experimental data. Solution accuracy and computational efficiency are concluded to be great with the Chebyshev pseudospectral method. Further applications to unsteady nonperiodic problems will be left for future work.


---

## Body

## 1. Introduction

Challenges in helicopter rotor flow analysis are several. Blade wake interactions need to be solved accurately requiring higher mesh resolution in the wake-trailing regions. Moving boundaries of the rotating blades need to be handled dynamically in time [1–5]. A typical local, finite difference method for the time integration of the flow governing equations of the URANS (unsteady Reynolds-averaged Navier-Stokes equations), such as BDF2 (2-step backward differentiation formula) or BDF3 (3-step backward differentiation formula), of the dual-time-stepping method [6, 7] requires a small time step for solution stability. Although the pseudo-time integration for the inner iteration and a multigrid technique can make the solution converge faster, the time-marching computation is still expensive even with highly parallel computation as it needs to fully resolve unsteady transient flows before flows reach a periodic steady-state. On the other hand, a spectral method for time integration has been proposed [8–16] as an alternative to a time-marching, dual-time stepping method with much reduced computational cost at equivalent solution accuracy. With the flow solution approximation by a discrete Fourier series, a time derivative term of the governing equations reduces to a spectral derivative term and removes time dependence of the flow solutions. The time-spectral derivative is the multiplication of the spectral derivative matrix and solution vectors at all time instances. The flow solutions advance in a pseudo-time domain with periodic steady-state remaining throughout the iterations until flows converge. Variants of the time-spectral method have been developed by various researchers. Hall et al. [8] originally developed a harmonic balance method (HBM) to solve internal flows of a two-dimensional compressor and use a discrete Fourier series for the solution approximation. The improved HBM by Thomas et al. [9, 10] also approximates flow variables by a discrete Fourier series and represents the time derivative term of the governing equations in Fourier coefficients but solves the FFT-transformed governing equations directly in the time domain by taking the inverse transform of the frequency domain governing equations. This approach is more convenient as the flow solutions are integrated directly in time with residuals checked simultaneously in the time domain. Based on the same principle, McMullen et al. [13, 14] developed a nonlinear frequency domain (NLFD) method which is slightly different from the HBM in the details of the implementation. Although the flow approximation by the Fourier series and the formulation of the governing equations in the frequency domain are the same as the HBM, the flow solutions are directly solved in the frequency domain, while the convergence is checked by the residuals reformulated in the time domain by taking the IFFT of the conservative flow variables. Gopinath and Jameson [12] proposed a time-spectral method which replaces the time derivative term of the governing equations through the spectral derivative matrix operator using Fourier spectral differentiation. Similar to the HBM developed by Thomas et al. [9, 10], it integrates the spectral derivative term of the governing equations directly in the time domain. Ekici et al. [11] solved helicopter rotor flows using the harmonic balance method with periodic boundary conditions both in time and in space, which makes the computation domain more compact than that of the time-spectral method and contributes to further decreasing computation time and memory in comparison to the time-marching computation. Choi et al. [17, 18] used the time-spectral CFD computation method for fluid-structure interaction analysis of the rotor flows in unsteady forward flight of the UH60-A helicopter. Adjoint sensitivity analysis has been carried out in the time-spectral form of the URANS equations to compute the sensitivity of aerodynamic performance metrics (drag, lift, torque, etc.) with respect to hundreds of shape design parameters. Despite the rigid blade assumption, aerodynamic shape optimization using time-spectral and adjoint sensitivity analysis was validated with the FSI analysis afterward with good performance and shown to be effective in solving unsteady (i.e., periodic steady) flows. A parameter study on the accuracy and computation efficiency of flow solutions with respect to varying number of harmonics has been carried out at various flight conditions. A Chebyshev pseudospectral method [19–22] uses Chebyshev polynomials for the solution approximation of the boundary-value problems, either periodic or nonperiodic, and the spectral derivative matrix is derived correspondingly to continuous differentiation. Standard Chebyshev polynomials are defined in the interval of [-1,1] and are a member of an orthogonal polynomial set. They can approximate function values with great precision and are applied practically to mathematics, physics, and engineering areas. Dinu et al. [19] applied the Chebyshev pseudospectral method to an aeroelastic problem and solved unsteady aerodynamic forces accurately. Niu et al. [20] solved unsteady cylinder flows with the Chebyshev pseudospectral method applied to the spatial derivative terms of the momentum and energy conservation equations. Flow solutions demonstrated great agreements with those of the time-marching computation which uses the multistage Runge-Kutta method. Khater et al. [21] applied the method to the spatial derivative term of the flow governing equations of the Burgers equations in various forms and excellent agreements were observed in comparison to the exact solutions. Im et al. [23] applied a Chebyshev collocation operator to the time derivative term of the URANS equations and employed a conformal mapping function to improve the numerical stability of the method. They applied the method to the oscillating airfoil with a prespecified frequency and compared the computation results with those of the dual-time stepping method and the harmonic balance method. In the current study, helicopter rotor flows in a hovering condition and forward level-flight are solved using the mapped Chebyshev pseudospectral method applied to the time integration of the URANS equations. To accurately solve the blade-vortex interaction and vortex wake flows trailing from the rotor blades, two key approaches are made: conformal mapping of an arcsin function and the overset mesh topology. The arcsin function is used for the conformal mapping of the uneven collocation points of the standard Chebyshev polynomials and yields points more evenly distributed over the physical time interval of [tinitial,tfinal]. Improved numerical stability is shown with respect to varying control parameters of the function. The application of the overset mesh topology is particularly advantageous with the mapped Chebyshev pseudospectral method, as the overset mesh construction and a hole-cutting procedure is carried out once and for all for subtime intervals corresponding to mapped collocation points of the Chebyshev polynomials. Cost for identifying hole-cutting topology between the near-body and background grids is reduced considerably when compared to the time-marching computation where it is carried out at each physical time step with the size often limited to a small value due to the numerical stability. The ratio of the number of subtime intervals between the time-spectral and the finite difference method is about (10) : O(102), and the corresponding cost saving is more than an order of magnitude. Despite some disadvantage of increase in computation memory per solution iteration due to the storage of meshes and flow solutions at all subtime intervals, the mapped Chebyshev pseudospectral method has shown cost effectiveness in the present study. In the current work, the overset mesh topology will be used for the mapped Chebyshev pseudospectral method to solve helicopter flows of Caradonna and Tung’s rotors in both hovering and unsteady nonlifting level flights. Aerodynamic results including section Cp distribution, thrust, and torque show excellent agreements in all computations and with wind-tunnel experimental data.


## 2. Mapped Chebyshev Pseudospectral Method

A Chebyshev pseudospectral method uses orthogonal Lagrange polynomials to evaluate state and control parameters of a given function to approximate or reconstruct. Either Legendre-Gauss-Lobatto (LGL) or Chebyshev-Gauss-Lobatto (CGL) points are used to evaluate function derivatives. However, the LGL or CGL points are heavily clustered around the end points of the interval, and the corresponding spectral derivative matrix becomes ill-conditioned with a large condition number. If the original function shows strong nonlinearity around the center of the interval where only a few CGL points are available for the spectral interpolation, the prediction of function value becomes very inaccurate around the region. The uneven distribution of the collocation points also causes numerical instability of the Chebyshev pseudospectral method [24–28]. Conformal mapping based on the arcsin function is applied to the CGL points to reduce the condition number of the spectral derivative matrix and improve the stability.


### 2.1. Chebyshev Collocation Operator [29]

An arbitrary functionu(t) defined in time over the interval of [a,b] can be represented as Chebyshev polynomials with N+1 GCL points:(1)ujt=∑n=0NanTntj=∑n=0Nancos⁡jπnN,forj=0,1,2,…,N(2)an=2Ncn∑j=0N1cjujTntj=2Ncn∑j=0Nujcjcos⁡jπnN,forn=0,1,2,…,N(3)ci=2,ifi=0,N1,otherwise.A vector of polynomial coefficients is written as(4)a→=DTu→=2NT0t04T0t12…T0tN4T1t02T1t1…T1tN2⋮⋮⋱⋮TNt04TNt12…TNtN4u0u1⋮uN.A vector of function derivatives can be written similarly as(5)u→′=D~Tb→=T0t0T1t0…TNt0T0t1T1t1…TNt1⋮⋮⋱⋮T0tNT1tN…TNtNb0b1⋮bN;then the coefficient vector of the function derivative is represented as(6)b→=Ga→=01030…G0,N00408…G1,N00060…G2,N⋮⋮⋮⋮⋮⋱⋮00000…GN,Na0a1a2⋮aN,(7)Gij=2jci,ifi+jeven0,otherwise,ci=2,ifi=0,N1,otherwise.A variable of Gi,j in (7) is defined only when the value of (i+j) becomes an even number and is zero for the odd number. For the interval of [a,b], not [-1,1], the elements of Gi,j are divided by (b-a)/2. Therefore, the time derivative of the arbitrary function is summarized as(8)u′=D~Tb→=D~TGa→=D~TGDTu→=Dchu→,where Dch represents the Chebyshev collocation operator and is a matrix of size [N+1,N+1].


### 2.2. Conformal Mapping

The purpose of the conformal mapping is to mitigate the problem of ill-conditioning of the spectral derivative matrix. Collocation points or nodes in a given time interval are rewritten by the mapping function ofg(yi), where yi is the standard Chebyshev collocation points represented as a cosine mesh.(9)ti=gyiyi=cos⁡iπN,i=0,1,…,N.To have an interval bounded by [0,a], not by [-1,1], a linear transformation is applied as(10)ti=1Agyi-B,(11)A=1Tmaxgy0-gyN,B=gyN.Then, the function derivative of u′ is written by the chain rule,(12)dudt=dydtdudy=1g′dudy.For the conformal mapping function of g(y), the arcsin function is used which was originally suggested by Kosloff and Tal-Ezer [25, 26]. A distribution characteristic is controlled by the constants of α,β with the value of α close to one resulting in more even distribution whereas smaller values push the collocation points towards the two ends of the interval. Analytic form of g(y) is(13)gy=sin-12αβy+α+βα+β,α=cos⁡iπN,β=cos⁡jπN.Variables of A and B, shown in (11), are the functions of α, β, and time and can be written as below over the time interval of [a0,a1], where a0 is typically set to be zero for the interval of time.(14)A=-1Tmaxε-θ,B=ε,ε=sin-12αβ+α-βα+β,θ=sin-1-2αβ+α-βα+β.The interval size of g(y) is shown in Figure 1. The graph shows the distribution of Δg(y) of inverse sine function. The number of mapping points is 37, where 0<g(y)<Tmax. The distribution of the inverse sine mapping function is consistent and uniform and is less sensitive to the g(y) values. A derivative function of the mapping function, g′, is as follows:(15)g′yi=Tmaxαβ1-αy1+βy.It is a diagonal matrix of size (N+1) with zero off-diagonal terms. Therefore, the mapped Chebyshev collocation operator with the mapping function applied is represented as(16)dudt=1g′dudy=αβTmax1-αy1+βyDch.The size of the matrix shown in (16) is N+1×(N+1), and it represents the mapped Chebyshev pseudospectral collocation operator mapped to [0,a] by a tangent function. The value of Dch corresponds to the interval of [-1,1] as 1/g′ takes into account the [0,a1]. Table 1 shows the values of the condition number for the standard and the mapped Chebyshev spectral derivative matrices with the total number of collocation points increasing from 24 to 96 by a factor of 2, where the values of α and β are optimized for the lowest condition number for each case. Reduction of the condition number by an order of magnitude is shown in all cases in Table 1 with the application of the conformal mapping of the arcsin function shown in (16).Table 1
Comparison of condition number for the differential operator.
i,j N α , β Cond. number (Mapped Cheby.) Cond. number (Stand. Cheby) 2 50 0.9918 3.1899×1016 2.7738 × 1 0 17 2 60 0.9943 6.8668×1016 3.6066 × 1 0 16 2 70 0.9958 1.9928×1016 6.1467 × 1 0 16 2 80 0.9968 3.8054×1016 3.6405 × 1 0 17 2 90 0.9975 6.9001×1016 2.1414 × 1 0 17Figure 1
Distribution ofΔg(y) of inverse sine mapping function.


## 2.1. Chebyshev Collocation Operator [29]

An arbitrary functionu(t) defined in time over the interval of [a,b] can be represented as Chebyshev polynomials with N+1 GCL points:(1)ujt=∑n=0NanTntj=∑n=0Nancos⁡jπnN,forj=0,1,2,…,N(2)an=2Ncn∑j=0N1cjujTntj=2Ncn∑j=0Nujcjcos⁡jπnN,forn=0,1,2,…,N(3)ci=2,ifi=0,N1,otherwise.A vector of polynomial coefficients is written as(4)a→=DTu→=2NT0t04T0t12…T0tN4T1t02T1t1…T1tN2⋮⋮⋱⋮TNt04TNt12…TNtN4u0u1⋮uN.A vector of function derivatives can be written similarly as(5)u→′=D~Tb→=T0t0T1t0…TNt0T0t1T1t1…TNt1⋮⋮⋱⋮T0tNT1tN…TNtNb0b1⋮bN;then the coefficient vector of the function derivative is represented as(6)b→=Ga→=01030…G0,N00408…G1,N00060…G2,N⋮⋮⋮⋮⋮⋱⋮00000…GN,Na0a1a2⋮aN,(7)Gij=2jci,ifi+jeven0,otherwise,ci=2,ifi=0,N1,otherwise.A variable of Gi,j in (7) is defined only when the value of (i+j) becomes an even number and is zero for the odd number. For the interval of [a,b], not [-1,1], the elements of Gi,j are divided by (b-a)/2. Therefore, the time derivative of the arbitrary function is summarized as(8)u′=D~Tb→=D~TGa→=D~TGDTu→=Dchu→,where Dch represents the Chebyshev collocation operator and is a matrix of size [N+1,N+1].


## 2.2. Conformal Mapping

The purpose of the conformal mapping is to mitigate the problem of ill-conditioning of the spectral derivative matrix. Collocation points or nodes in a given time interval are rewritten by the mapping function ofg(yi), where yi is the standard Chebyshev collocation points represented as a cosine mesh.(9)ti=gyiyi=cos⁡iπN,i=0,1,…,N.To have an interval bounded by [0,a], not by [-1,1], a linear transformation is applied as(10)ti=1Agyi-B,(11)A=1Tmaxgy0-gyN,B=gyN.Then, the function derivative of u′ is written by the chain rule,(12)dudt=dydtdudy=1g′dudy.For the conformal mapping function of g(y), the arcsin function is used which was originally suggested by Kosloff and Tal-Ezer [25, 26]. A distribution characteristic is controlled by the constants of α,β with the value of α close to one resulting in more even distribution whereas smaller values push the collocation points towards the two ends of the interval. Analytic form of g(y) is(13)gy=sin-12αβy+α+βα+β,α=cos⁡iπN,β=cos⁡jπN.Variables of A and B, shown in (11), are the functions of α, β, and time and can be written as below over the time interval of [a0,a1], where a0 is typically set to be zero for the interval of time.(14)A=-1Tmaxε-θ,B=ε,ε=sin-12αβ+α-βα+β,θ=sin-1-2αβ+α-βα+β.The interval size of g(y) is shown in Figure 1. The graph shows the distribution of Δg(y) of inverse sine function. The number of mapping points is 37, where 0<g(y)<Tmax. The distribution of the inverse sine mapping function is consistent and uniform and is less sensitive to the g(y) values. A derivative function of the mapping function, g′, is as follows:(15)g′yi=Tmaxαβ1-αy1+βy.It is a diagonal matrix of size (N+1) with zero off-diagonal terms. Therefore, the mapped Chebyshev collocation operator with the mapping function applied is represented as(16)dudt=1g′dudy=αβTmax1-αy1+βyDch.The size of the matrix shown in (16) is N+1×(N+1), and it represents the mapped Chebyshev pseudospectral collocation operator mapped to [0,a] by a tangent function. The value of Dch corresponds to the interval of [-1,1] as 1/g′ takes into account the [0,a1]. Table 1 shows the values of the condition number for the standard and the mapped Chebyshev spectral derivative matrices with the total number of collocation points increasing from 24 to 96 by a factor of 2, where the values of α and β are optimized for the lowest condition number for each case. Reduction of the condition number by an order of magnitude is shown in all cases in Table 1 with the application of the conformal mapping of the arcsin function shown in (16).Table 1
Comparison of condition number for the differential operator.
i,j N α , β Cond. number (Mapped Cheby.) Cond. number (Stand. Cheby) 2 50 0.9918 3.1899×1016 2.7738 × 1 0 17 2 60 0.9943 6.8668×1016 3.6066 × 1 0 16 2 70 0.9958 1.9928×1016 6.1467 × 1 0 16 2 80 0.9968 3.8054×1016 3.6405 × 1 0 17 2 90 0.9975 6.9001×1016 2.1414 × 1 0 17Figure 1
Distribution ofΔg(y) of inverse sine mapping function.


## 3. Implementation for Unsteady Reynolds-Averaged Navier-Stokes Equations

Unsteady Reynolds-averaged Navier-Stokes (URANS) equations in three dimensions are written as(17)∂Q∂t+∂F∂x+∂G∂y+∂H∂z=0,where Q represents conservative variables and the solutions of RANS equations. Flux vectors in x, y, and z directions are F, G, and H, respectively, and include both viscous and inviscid fluxes. To apply the spectral method, vectors of solutions and fluxes of Q, F, G, and H are approximated by Chebyshev polynomials with N+1 CGL points.(18)Q→jt=∑n=0NQ^chnTntj,forj=0,1,2,…,N,F→jt=∑n=0NF^chnTntj,forj=0,1,2,…,N,G→jt=∑n=0NG^chnTntj,forj=0,1,2,…,N,H→jt=∑n=0NH^chnTntj,forj=0,1,2,…,N.Unlike the Fourier spectral method which uses a discrete Fourier series with equally distributed collocation points, the cosine mesh points of the standard Chebyshev polynomials have uneven distribution and those are redistributed by the mapping function of g(y) as(19)tj=Tmaxε-θgyi-ε.After we replace the time derivative term in (17) by the mapped Chebyshev pseudospectral derivative operator, then (17) can be rewritten as(20)1g′DchQ→+∂F→∂x+∂G→∂y+∂H→∂z=0→,where Q→, F→, G→, and H→ are the vectors of size (N+1), and their elements correspond to mapped collocation points by the arcsin function.(21)Q→=Q0Q1Q2⋮Qn,F→=F0F1F2⋮Fn,G→=G0G1G2⋮Gn,H→=H0H1H2⋮Hn(22)∂Q→∂τ+1g′DchQ→+∂F→∂x+∂G→∂y+∂H→∂z=0→.


### 3.1. Partially Implicit Discretization with Chebyshev Pseudospectral Term

In this study, for stable time integration the diagonal ADI technique was applied to Chebyshev pseudospectral approach. Chebyshev pseudospectral source term is considered explicitly, so the influence of the source term in nondiagonal matrix is ignored. For this reason the source term is added in flux term at subtime level. Though the Chebyshev pseudospectral source term is used explicitly, Diagonal ADI technique that treats together with the flux term and the source term can give the implicit convergence characteristic. In other words, the convergence rate of the existing time domain solver is maintained. First, the combination of explicit source and implicit residual value is derived as(23)∂Q→∂τ+1g′DchQ→n+∂F→n+1∂x+∂G→n+1∂y+∂H→n+1∂z=0→.The flux terms of F→(n+1), G→(n+1), and H→(n+1) are linearized with flux Jacobian matrices of V, W, and Z, (24)F→n+1=F→n+VΔQ→n,G→n+1=G→n+WΔQ→n,H→n+1=H→n+ZΔQ→n,where V, W, and Z are the flux Jacobian matrices of size (N+1)×(N+1) with each element to be a 5 × 5 submatrix and can be expressed as (25)V=∂F0∂F0∂F0∂F1…∂F0∂FN∂F1∂F0∂F1∂F1…∂F1∂FN⋮⋱⋮∂FN∂F0…∂FN∂FN-1∂FN∂FN,W=∂G0∂Q0∂G0∂Q1…∂G0∂QN∂G1∂Q0∂G1∂Q1…∂G1∂QN⋮⋱⋮∂GN∂Q0…∂GN∂QN-1∂GN∂QN,Z=∂H0∂Q0∂H0∂Q1…∂H0∂QN∂H1∂Q0∂H1∂Q1…∂H1∂QN⋮⋱⋮∂HN∂Q0…∂HN∂QN-1∂HN∂QN.Substituting (24) in (23) can be expressed as(26)ΔQnΔτ+∂∂xVΔQ→n+∂∂yWΔQ→n=-∂G→n∂x-∂H→n∂y-1g′DchQ→n.The flux term in (23) is a function of variables of each subtime level, and (25) for each element except the diagonal term can be considered 0. Therefore, simplified implicit flux Jacobian matrix is represented as (27)V=∂F0∂F00…00∂F1∂F1…0⋮⋱⋮0…0∂FN∂FN,(28)W=∂G0∂Q00…00∂G1∂Q1…0⋮⋱⋮0…0∂GN∂QN,(29)Z=∂H0∂Q00…00∂H1∂Q1…0⋮⋱⋮0…0∂HN∂QN.The element of the flux Jacobians of V, W, and Z is a 5 × 5 submatrix and corresponds to the mapped collocation point. A noteworthy point here is that a simplified diagonal Jacobian matrix is the same as the flux Jacobian matrix in time domain solver. Therefore, the additional Jacobian matrix does not create in this method to be able to apply the Chebyshev pseudospectral method easily. With (29) substituted and multiplied by infinitesimally small time step Δτ, (26) is rewritten as(30)I+δxVΔτ+δyWΔτ+δzZΔτΔQ→n=-∂F→n∂xΔτ-∂G→n∂yΔτ-∂H→n∂zΔτ-Δτg′DchQ→n(31)T00…00T1…0⋮⋮⋱⋮00…TNΔQ0ΔQ1⋮ΔQN=S0S1⋮SN,where(32)Ti=I+δxViiΔτ+δyWiiΔτ+δzZiiΔτSi=-∂F→in∂xΔτ-∂G→in∂yΔτ-∂H→in∂zΔτ-Δτg′DchQ→in.The final form of the mapped Chebyshev pseudospectral method implemented in the URANS equations has a similar discretized structure of the steady form of RANS equations with Chebyshev derivative source terms added to the right hand side.


## 3.1. Partially Implicit Discretization with Chebyshev Pseudospectral Term

In this study, for stable time integration the diagonal ADI technique was applied to Chebyshev pseudospectral approach. Chebyshev pseudospectral source term is considered explicitly, so the influence of the source term in nondiagonal matrix is ignored. For this reason the source term is added in flux term at subtime level. Though the Chebyshev pseudospectral source term is used explicitly, Diagonal ADI technique that treats together with the flux term and the source term can give the implicit convergence characteristic. In other words, the convergence rate of the existing time domain solver is maintained. First, the combination of explicit source and implicit residual value is derived as(23)∂Q→∂τ+1g′DchQ→n+∂F→n+1∂x+∂G→n+1∂y+∂H→n+1∂z=0→.The flux terms of F→(n+1), G→(n+1), and H→(n+1) are linearized with flux Jacobian matrices of V, W, and Z, (24)F→n+1=F→n+VΔQ→n,G→n+1=G→n+WΔQ→n,H→n+1=H→n+ZΔQ→n,where V, W, and Z are the flux Jacobian matrices of size (N+1)×(N+1) with each element to be a 5 × 5 submatrix and can be expressed as (25)V=∂F0∂F0∂F0∂F1…∂F0∂FN∂F1∂F0∂F1∂F1…∂F1∂FN⋮⋱⋮∂FN∂F0…∂FN∂FN-1∂FN∂FN,W=∂G0∂Q0∂G0∂Q1…∂G0∂QN∂G1∂Q0∂G1∂Q1…∂G1∂QN⋮⋱⋮∂GN∂Q0…∂GN∂QN-1∂GN∂QN,Z=∂H0∂Q0∂H0∂Q1…∂H0∂QN∂H1∂Q0∂H1∂Q1…∂H1∂QN⋮⋱⋮∂HN∂Q0…∂HN∂QN-1∂HN∂QN.Substituting (24) in (23) can be expressed as(26)ΔQnΔτ+∂∂xVΔQ→n+∂∂yWΔQ→n=-∂G→n∂x-∂H→n∂y-1g′DchQ→n.The flux term in (23) is a function of variables of each subtime level, and (25) for each element except the diagonal term can be considered 0. Therefore, simplified implicit flux Jacobian matrix is represented as (27)V=∂F0∂F00…00∂F1∂F1…0⋮⋱⋮0…0∂FN∂FN,(28)W=∂G0∂Q00…00∂G1∂Q1…0⋮⋱⋮0…0∂GN∂QN,(29)Z=∂H0∂Q00…00∂H1∂Q1…0⋮⋱⋮0…0∂HN∂QN.The element of the flux Jacobians of V, W, and Z is a 5 × 5 submatrix and corresponds to the mapped collocation point. A noteworthy point here is that a simplified diagonal Jacobian matrix is the same as the flux Jacobian matrix in time domain solver. Therefore, the additional Jacobian matrix does not create in this method to be able to apply the Chebyshev pseudospectral method easily. With (29) substituted and multiplied by infinitesimally small time step Δτ, (26) is rewritten as(30)I+δxVΔτ+δyWΔτ+δzZΔτΔQ→n=-∂F→n∂xΔτ-∂G→n∂yΔτ-∂H→n∂zΔτ-Δτg′DchQ→n(31)T00…00T1…0⋮⋮⋱⋮00…TNΔQ0ΔQ1⋮ΔQN=S0S1⋮SN,where(32)Ti=I+δxViiΔτ+δyWiiΔτ+δzZiiΔτSi=-∂F→in∂xΔτ-∂G→in∂yΔτ-∂H→in∂zΔτ-Δτg′DchQ→in.The final form of the mapped Chebyshev pseudospectral method implemented in the URANS equations has a similar discretized structure of the steady form of RANS equations with Chebyshev derivative source terms added to the right hand side.


## 4. Numerical Validation

### 4.1. Pitching Airfoil with a Given Frequency

Helicopter rotor blades are in pitching and flapping motions during the forward flight to overcome the imbalance of aerodynamic forces and moments between the advancing and the retreating side. Flows on the advancing side have higher speed due to the positive contribution from rotation speed, whereas flows on the retreating side experience lower speed due to negative contribution. This causes imbalance of the aerodynamic forces and moments, and to overcome the issue a pitch angle is lowered on the advancing side and increased on the retreating side. In earlier research on the helicopter rotor blade design, aerodynamic analyses of pitching airfoil have been carried out in various wind-tunnel facilities to investigate the associated flow characteristics. Experimental data from AGARD Report 702 [30] are used in the present study to validate the mapped Chebyshev pseudospectral method. A C-type mesh of size 309 × 81 is constructed and shown in Figure 2 to solve two-dimensional URANS equations. Mach number is 0.755 and a pitching motion is described by(33)αAOA=α0+αmsin⁡kcMinft∗,where the damping frequency is kc=0.1628∘, the mean angle of attack is α0=0.016∘, and the pitching amplitude is αm=2.51∘. t∗ is the nondimensional time and Reynolds number based on the chord length is 5.5 × 106. For the pseudo time integration, the mapped Chebyshev spectral derivative term, a DADI (diagonal alternate direction implicit), method is applied [15, 23, 31]. A Roe’s upwind flux-difference splitting (FDS) method [32] is employed for the flux term discretization with a third order MUSCL scheme to improve solution accuracy. A two-equation turbulence model of k-ω WD+ [33] is adopted in the flow solver to compute the turbulent eddy viscosity coefficient. The number of collocations for the mapped Chebyshev polynomial approximation varies from 4 to 15, with the values of α and β used in the arcsin mapping function to be 0.99. A convergence history of the l2 norm of the density residual is shown in Figure 3 over solution iterations. Results of the mapped Chebyshev pseudospectral computation with the varying number of collocation points or subtime intervals are plotted in different colors and symbols. The HBM computation with 9 harmonics, which correspond to 19 subtime intervals due to the characteristics of the conjugate symmetry for the real values, is carried out and the results are compared with those of the mapped Chebyshev pseudospectral computation. Similar convergence rates are obtained for the mapped Chebyshev pseudospectral computation regardless of the number of subtime intervals, although CPU time per solution iteration increases linearly to the number of subtime intervals. The similar convergence rate between the HBM and the mapped Chebyshev pseudospectral method is explained by the fact that both methods have the same spectral form of the unsteady governing equations with the source term of the spectral derivative operator added to the right hand side while the details of the spectral derivative matrix are different. Computation results of the unsteady lift, drag, and moment coefficients are plotted in Figure 4, respectively, with respect to the varying number of subtime intervals, and compared with those of the harmonic balance computation with 9 harmonics (i.e., 19 subtime intervals). For the mapped Chebyshev pseudospectral computation, solutions are plotted as symbols in different colors only at discrete collocation points corresponding to the specific subtime intervals. The HBM results are plotted in continuous solid black lines constructed through the Fourier spectral interpolation. With slight discrepancies in the lift and the drag coefficient in computations with 4 and 6 subtime intervals, the present method with subtime intervals more than 6 seems to converge with good agreements with those of the HBM computation.Figure 2
NACA0012 Grid (309 × 81).Figure 3
Comparison of convergence history.Figure 4
Comparison of lift, drag, and moment coefficients.
(a)
Lift coefficient (b)
Drag coefficient (c)
Moment coefficient


### 4.2. Plunging Airfoil

One of the biggest advantages of the Chebyshev pseudospectral method is the capability to solve nonperiodic unsteady problems if the positions at the two end points of the time interval are known a priori. The Fourier spectral method, or the nonlinear frequency domain method, uses the Fourier transform which requires the specification of the motion frequency a priori with the periodic condition imposed at the ends of the time interval; however, the Chebyshev pseudospectral method has flexibility in having different values at the two boundary points of the time interval and is applicable to nonperiodic problems. Thus, the second validation problem is chosen as a plunging airfoil which moves downwards at a constant plunging speed prespecified in (34), where external gravity force is ignored. Freestream Mach number is 0.8 in a horizontal direction with the initial angle of attack of the airfoil to be zero. With the end point of the time interval arbitrarily set, Chebyshev collocation points are distributed based on the plunging speed and distance (Figure 5).(34)y˙=-P1·ωp,(35)ωp=κU∞Chord,(36)P1=0.0436·Chord,(37)κ=0.4.Moreover, it is reported in the reference paper [23] that unsteady flows around an airfoil plunging at constant speed with nonzero freestream velocity are equivalent to steady flows of an stationary airfoil fixed at a certain angle of attack with the same freestream velocity, which is explained graphically in Figures 6(a) and 6(b). An overset mesh topology for the plunging airfoil is constructed to facilitate the dynamic motion of the airfoil with a near-body grid of 257 × 45 and a background grid of 301 × 151 as shown in Figure 7.Figure 5
Chebyshev points in plunging airfoil.Figure 6
An airfoil plunging at constant speedy˙ and a stationary airfoil at the AOA of α.
(a)
Present method (b)
Time accurate methodFigure 7
An overset mesh topology for the plunging airfoil.Figure8(a) shows the schematics of the conventional time-marching method with the overset grid technique, which uses a local, finite difference-based approach. The overset mesh technique consists of two major steps, hole-cutting and donor cell identification. Wall boundaries are used for the hole-cutting and the cutting boundaries are constituted with a Cartesian mesh. The cutting surfaces are approximately made up with a collection of cubes. A zone of interface scheme and a hole-map algorithm are used for the chimera hole-cutting. For the donor cell identification, a stencil walk and gradient searching method is applied [34]. The stencil walk finds the nearest point to a given interpolation point by comparing the distance between the interpolation point and the donor candidate cells. The gradient search method checks whether the interpolation point actually lies inside the cube that is made with the eight donor cells or not [35]. For the time-marching method, time integration is carried out by a finite difference method, and therefore the application of the body movement, subgrid movement, hole-cutting and donor cell identification, and solving the RANS equations per dual-time step are made iteratively over the entire interval until flow solutions converge. An initial transient flow state needs to be fully resolved before flow reaches a periodic steady-state. Therefore, the overset grid technique needs to cut surfaces with a zone of interface scheme and hole-map algorithm and to interpolate with donor cell identification every physical time step. On the other hand, the mapped Chebyshev pseudospectral method constructs and stores a set of subgrid meshes, hole-cutting, and donor cell information corresponding to the mapped CGL points once at the beginning of the solution process, and the information is not modified throughout the iterations until flow solutions converge in a pseudo-time domain. This fact is greatly favorable in efficiency by reducing computation time considerably, which is more so when the size of grid becomes larger. Although it should be noted that computation memory increases accordingly from storing the meshes and flow solutions at all subgrid intervals, the overall computation efficiency has been shown higher in the Chebyshev pseudospectral computation for this plunging airfoil problem and the helicopter rotor flow analysis discussed in Section 5; Figures 9 and 10 show the contours of pressure and Mach number, respectively, around the airfoils. Contours shown in Figures 9(a) and 10(a) are computed for the plunging airfoil by the mapped Chebyshev pseudospectral method, whereas those shown in Figures 9(b) and 10(b) are solved by the steady RANS computation for the stationary airfoil with the angle of attack fixed at 8∘. A single block, C-type structured mesh is used for the steady RANS computation. Comparisons shown in Figures 9 and 10 show excellent agreements between the two computations, although slightly nonsmooth contours are observed in the Chebyshev pseudospectral computation especially in the overlapping mesh regions between the near-body and the background grid. Figure 11 shows a direct comparison of the surface pressure coefficients with good agreements.Figure 8
Comparison of the computation procedures between the time-marching and the Chebyshev pseudospectral method.
(a)
Conventional finite difference-based time integration (b)
Chebyshev pseudospectral methodFigure 9
Pressure contours.
(a)
Results by the Chebyshev spectral method for the plunging airfoil (b)
Results by the steady computation for the stationary airfoil withα=0Figure 10
Mach number contours.
(a)
Results by the Chebyshev spectral method for the plunging airfoil (b)
Results by the steady computation for the stationary airfoil withα=0Figure 11
Comparison of the surface pressure coefficients.


## 4.1. Pitching Airfoil with a Given Frequency

Helicopter rotor blades are in pitching and flapping motions during the forward flight to overcome the imbalance of aerodynamic forces and moments between the advancing and the retreating side. Flows on the advancing side have higher speed due to the positive contribution from rotation speed, whereas flows on the retreating side experience lower speed due to negative contribution. This causes imbalance of the aerodynamic forces and moments, and to overcome the issue a pitch angle is lowered on the advancing side and increased on the retreating side. In earlier research on the helicopter rotor blade design, aerodynamic analyses of pitching airfoil have been carried out in various wind-tunnel facilities to investigate the associated flow characteristics. Experimental data from AGARD Report 702 [30] are used in the present study to validate the mapped Chebyshev pseudospectral method. A C-type mesh of size 309 × 81 is constructed and shown in Figure 2 to solve two-dimensional URANS equations. Mach number is 0.755 and a pitching motion is described by(33)αAOA=α0+αmsin⁡kcMinft∗,where the damping frequency is kc=0.1628∘, the mean angle of attack is α0=0.016∘, and the pitching amplitude is αm=2.51∘. t∗ is the nondimensional time and Reynolds number based on the chord length is 5.5 × 106. For the pseudo time integration, the mapped Chebyshev spectral derivative term, a DADI (diagonal alternate direction implicit), method is applied [15, 23, 31]. A Roe’s upwind flux-difference splitting (FDS) method [32] is employed for the flux term discretization with a third order MUSCL scheme to improve solution accuracy. A two-equation turbulence model of k-ω WD+ [33] is adopted in the flow solver to compute the turbulent eddy viscosity coefficient. The number of collocations for the mapped Chebyshev polynomial approximation varies from 4 to 15, with the values of α and β used in the arcsin mapping function to be 0.99. A convergence history of the l2 norm of the density residual is shown in Figure 3 over solution iterations. Results of the mapped Chebyshev pseudospectral computation with the varying number of collocation points or subtime intervals are plotted in different colors and symbols. The HBM computation with 9 harmonics, which correspond to 19 subtime intervals due to the characteristics of the conjugate symmetry for the real values, is carried out and the results are compared with those of the mapped Chebyshev pseudospectral computation. Similar convergence rates are obtained for the mapped Chebyshev pseudospectral computation regardless of the number of subtime intervals, although CPU time per solution iteration increases linearly to the number of subtime intervals. The similar convergence rate between the HBM and the mapped Chebyshev pseudospectral method is explained by the fact that both methods have the same spectral form of the unsteady governing equations with the source term of the spectral derivative operator added to the right hand side while the details of the spectral derivative matrix are different. Computation results of the unsteady lift, drag, and moment coefficients are plotted in Figure 4, respectively, with respect to the varying number of subtime intervals, and compared with those of the harmonic balance computation with 9 harmonics (i.e., 19 subtime intervals). For the mapped Chebyshev pseudospectral computation, solutions are plotted as symbols in different colors only at discrete collocation points corresponding to the specific subtime intervals. The HBM results are plotted in continuous solid black lines constructed through the Fourier spectral interpolation. With slight discrepancies in the lift and the drag coefficient in computations with 4 and 6 subtime intervals, the present method with subtime intervals more than 6 seems to converge with good agreements with those of the HBM computation.Figure 2
NACA0012 Grid (309 × 81).Figure 3
Comparison of convergence history.Figure 4
Comparison of lift, drag, and moment coefficients.
(a)
Lift coefficient (b)
Drag coefficient (c)
Moment coefficient


## 4.2. Plunging Airfoil

One of the biggest advantages of the Chebyshev pseudospectral method is the capability to solve nonperiodic unsteady problems if the positions at the two end points of the time interval are known a priori. The Fourier spectral method, or the nonlinear frequency domain method, uses the Fourier transform which requires the specification of the motion frequency a priori with the periodic condition imposed at the ends of the time interval; however, the Chebyshev pseudospectral method has flexibility in having different values at the two boundary points of the time interval and is applicable to nonperiodic problems. Thus, the second validation problem is chosen as a plunging airfoil which moves downwards at a constant plunging speed prespecified in (34), where external gravity force is ignored. Freestream Mach number is 0.8 in a horizontal direction with the initial angle of attack of the airfoil to be zero. With the end point of the time interval arbitrarily set, Chebyshev collocation points are distributed based on the plunging speed and distance (Figure 5).(34)y˙=-P1·ωp,(35)ωp=κU∞Chord,(36)P1=0.0436·Chord,(37)κ=0.4.Moreover, it is reported in the reference paper [23] that unsteady flows around an airfoil plunging at constant speed with nonzero freestream velocity are equivalent to steady flows of an stationary airfoil fixed at a certain angle of attack with the same freestream velocity, which is explained graphically in Figures 6(a) and 6(b). An overset mesh topology for the plunging airfoil is constructed to facilitate the dynamic motion of the airfoil with a near-body grid of 257 × 45 and a background grid of 301 × 151 as shown in Figure 7.Figure 5
Chebyshev points in plunging airfoil.Figure 6
An airfoil plunging at constant speedy˙ and a stationary airfoil at the AOA of α.
(a)
Present method (b)
Time accurate methodFigure 7
An overset mesh topology for the plunging airfoil.Figure8(a) shows the schematics of the conventional time-marching method with the overset grid technique, which uses a local, finite difference-based approach. The overset mesh technique consists of two major steps, hole-cutting and donor cell identification. Wall boundaries are used for the hole-cutting and the cutting boundaries are constituted with a Cartesian mesh. The cutting surfaces are approximately made up with a collection of cubes. A zone of interface scheme and a hole-map algorithm are used for the chimera hole-cutting. For the donor cell identification, a stencil walk and gradient searching method is applied [34]. The stencil walk finds the nearest point to a given interpolation point by comparing the distance between the interpolation point and the donor candidate cells. The gradient search method checks whether the interpolation point actually lies inside the cube that is made with the eight donor cells or not [35]. For the time-marching method, time integration is carried out by a finite difference method, and therefore the application of the body movement, subgrid movement, hole-cutting and donor cell identification, and solving the RANS equations per dual-time step are made iteratively over the entire interval until flow solutions converge. An initial transient flow state needs to be fully resolved before flow reaches a periodic steady-state. Therefore, the overset grid technique needs to cut surfaces with a zone of interface scheme and hole-map algorithm and to interpolate with donor cell identification every physical time step. On the other hand, the mapped Chebyshev pseudospectral method constructs and stores a set of subgrid meshes, hole-cutting, and donor cell information corresponding to the mapped CGL points once at the beginning of the solution process, and the information is not modified throughout the iterations until flow solutions converge in a pseudo-time domain. This fact is greatly favorable in efficiency by reducing computation time considerably, which is more so when the size of grid becomes larger. Although it should be noted that computation memory increases accordingly from storing the meshes and flow solutions at all subgrid intervals, the overall computation efficiency has been shown higher in the Chebyshev pseudospectral computation for this plunging airfoil problem and the helicopter rotor flow analysis discussed in Section 5; Figures 9 and 10 show the contours of pressure and Mach number, respectively, around the airfoils. Contours shown in Figures 9(a) and 10(a) are computed for the plunging airfoil by the mapped Chebyshev pseudospectral method, whereas those shown in Figures 9(b) and 10(b) are solved by the steady RANS computation for the stationary airfoil with the angle of attack fixed at 8∘. A single block, C-type structured mesh is used for the steady RANS computation. Comparisons shown in Figures 9 and 10 show excellent agreements between the two computations, although slightly nonsmooth contours are observed in the Chebyshev pseudospectral computation especially in the overlapping mesh regions between the near-body and the background grid. Figure 11 shows a direct comparison of the surface pressure coefficients with good agreements.Figure 8
Comparison of the computation procedures between the time-marching and the Chebyshev pseudospectral method.
(a)
Conventional finite difference-based time integration (b)
Chebyshev pseudospectral methodFigure 9
Pressure contours.
(a)
Results by the Chebyshev spectral method for the plunging airfoil (b)
Results by the steady computation for the stationary airfoil withα=0Figure 10
Mach number contours.
(a)
Results by the Chebyshev spectral method for the plunging airfoil (b)
Results by the steady computation for the stationary airfoil withα=0Figure 11
Comparison of the surface pressure coefficients.


## 5. Helicopter Rotor Flow Analysis

### 5.1. Hovering Condition

A geometry of the Caradonna and Tung’s rotor is popularly used to validate the hovering condition of the helicopter rotor flow analysis [36]. The rotor system has two blades which have a rectangular-shape planform with no twist and taper in the span direction, and their cross section is NACA0012 airfoil. The span length is 3.75 ft (1.143 m) with an aspect ratio of 7. A collective pitch angle of 8∘ is applied to the blades at all azimuth angles with Mach numbers of 0.439 and 0.877. All computations are carried out in parallel in a Linux cluster using a total of 56 CPUs of i7 Core processors with 24 GB RAM per node, which is effective in handling the increased memory requirement for the mapped Chebyshev pseudospectral method with meshes and solutions at all subtime intervals. Figure 12 shows the multiblock, overset mesh topology used for the Caradonna and Tung’s rotor. A near-body mesh is H-type with 155×239×239 points and the background mesh has the grid of 97×106×217. A total of 11 Chebyshev collocation points, corresponding to 11 subtime intervals due to the periodicity of the rotor motion, are used for the analysis of the hovering flight. Computations of the time-marching method, the HBM, and the mapped Chebyshev pseudospectral method are carried out independently and resultant aerodynamic forces and moments are compared with the wind-tunnel experimental data. For the HBM and the mapped Chebyshev pseudospectral method, solution reconstruction is carried out at the postprocessing stage based on the flow solutions at a finite number of subtime intervals or time instances. A graphical schematic of the flow reconstruction is illustrated in Figure 13 for the helicopter rotor flow analysis. Flow solutions available only at the discrete and unevenly distributed mapped Chebyshev collocation points are interpolated at an arbitrary time instance in the time interval using the following relations. (38)Qjt=∑n=0NanTntj=∑n=0Nancos⁡jπMchn,forj=0,1,2,…,Mch,(39)an=2Ncn∑n=0N1cjQjTntj=2Ncn∑n=0N1cjQjcos⁡jπN,forn=0,1,2,…,N,(40)ci=2,ifi=0,N1,otherwise.A value of M represents the number of points on the interpolant function which will be constructed from the N+1 number of collocation points or subtime intervals. Therefore, Qk is the value from the new approximation through the reconstruction, whereas Qj represents the value at the original subtime intervals. Eqns. (38) and (39) are the same as (1) and (2); however, it should be noted that the values of M and N are different. Figures 14 and 15 show the plots of the surface pressure distribution on the blade at Mach numbers of 0.439 and 0.877, respectively. Spanwise distribution of surface pressure coefficients is plotted at different spanwise locations of 0.5, 0.68, 0.80, 0.89, and 0.96 of the blade radius. Results of the mapped Chebyshev pseudospectral method, HBM, and the time-marching computation show excellent agreements with each other at all spanwise sections, although slight discrepancies are observed in comparison to the wind-tunnel experiments around the region of maximum pressure suction. A total of 19 collocation points are employed for the HBM and 19 points for the Chebyshev pseudospectral method and reconstructed with 180 number of interpolation points. Comparison of thrust and torque, CT and CQ, as well as the figure of merit (FM) is shown in Table 2 where the computations are carried with the Mach number of 0.439. All computation results match very well with those of the wind-tunnel experiments [36]. Surface pressure distributions on two blades at initial nine mapped Chebyshev collocation points or time instances are shown in Figure 16(a). Those can be used for the reconstruction by the interpolation given by (38) and Figure 16(b) shows the blades and surface pressure distribution reconstructed at 36 subtime intervals or time instances.Table 2
Aerodynamic coefficients at Mach number of 0.439.
CT CQ F.M. Experiment 0.0046 Time-accurate Method 0.0051 0.00049 0.53 Harmonic Balance Method 0.0051 0.00049 0.53 Present (Chebyshev spectral Method) 0.0051 0.00049 0.53Figure 12
A multiblock, overset mesh topology around Caradonna and Tung’s rotor blade with the background grid of97×106×217 and the near-body grid of 155×239×239.Figure 13
Schematic for solution reconstruction using the mapped Chebyshev points.Figure 14
Surface pressure coefficients at the radial locations of 0.5, 0.68, 0.8, 0.89, and 0.96 at Mach number of 0.439.
(a)
r / R = 0.50 (b)
r / R = 0.68 (c)
r / R = 0.80 (d)
r / R = 0.89 (e)
r / R = 0.96Figure 15
Surface pressure coefficients at the radial locations of 0.5, 0.68, 0.8, 0.89, and 0.96 at Mach number of 0.877.
(a)
r / R = 0.50 (b)
r / R = 0.68 (c)
r / R = 0.80 (d)
r / R = 0.89 (e)
r / R = 0.96Figure 16
Contour plots of the blade surface pressure at Mach = 0.877.
(a)
Blades and surface pressure contours at 9 mapped Chebyshev collocation points (b)
Blades and surface pressure contours reconstructed at 36 collocation points


### 5.2. Forward Flight

The same Caradonna and Tung rotors are solved at a nonlifting forward flight condition. Blades are not in a pitching or flapping motion and a collective pitch angle is set to be zero. Only forward speed is considered with corresponding tip Mach number of 0.5 and the advance ratio of 0.2. The CPU time is shown in Table3 with respect to the number of Chebyshev collocation points from 11 to 23 and is normalized by the computational time applied with one Chebyshev collocation point. The CPU time is observed with linearly increasing, which can be inferred from the diagonal form of the flux Jacobian matrices shown in (27)~(29) and (31). Table 3 also shows overset grid module time which is defined as the time for hole-cutting and donor cell identification with respect to the number of Chebyshev collocation points. For the time-marching computation for the unsteady analysis with this overset grid technique, the grid module time is about 3.4 seconds per physical time step which means 1° azimuth angle, which is inferred to be more than 20 minutes for one revolution which means 360° azimuth angle. The rotor flow usually converges after five revolutions and it only takes 1.7 hours for overset grid module time. However, by the present mapped Chebyshev pseudospectral method, the grid module time for the hole-cutting and mesh interpolation needs is needed only once at the beginning of the solution procedure and takes less than 2 minutes even with 23 Chebyshev points. Figure 17 shows the plots of the pressure coefficients at the 0.89 radial location of the blade when the azimuth angle varies from 30 deg to 180 deg by 30 deg. Results of the mapped Chebyshev pseudospectral method are computed with a total of 11, 13, and 15 collocation points and compared with the wind-tunnel data [37]. Strong shock appears in the advancing side of the rotor, approximately from 60 deg to 150 deg of the azimuth angle, demonstrating solution accuracy higher with more Chebyshev points. Computation results with a total of 11 points predict relatively well the shock structure in terms of the location and magnitude. Figure 18 shows the contours of surface pressure distribution on the blades at the mapped Chebyshev points which vary from 7, 9, to 11. Figure 18(d) shows the contours of surface pressure distribution on the blades reconstructed from the 11 points to 18 points. We checked the wall clock time with increasing interval numbers of HBM and present method in Figure 19. For time accurate method the dual-time stepping is applied with 360 intervals for one cycle and 20 subiteration for each interval, and three cycles for the periodic convergence of the unsteady flow. The computational time of present method is increasing linearly with respect to the number of intervals as like HBM, and Figure 19 shows present method is as efficient as HBM below 19 intervals.Table 3
Nondimensionalized time with respect to the varying number of Chebyshev collocation points and corresponding grid module time for mesh interpolation and hole cutting.
Chebyshev point CPU timeTn/T1 Grid module time (sec.) 11 7.21 37 15 15.02 51 19 19.81 63 23 25.48 76Figure 17
Pressure coefficients at the 0.89 radial location during the forward flight with an azimuth angle varying from 30 deg to 180 deg by 30 deg.
(a)
ψ = 30 deg (b)
ψ = 60 deg (c)
ψ = 90 deg (d)
ψ = 120 deg (e)
ψ = 150 deg (f)
ψ = 180 degFigure 18
Surface pressure distribution on the blades during the forward flight at the mapped Chebyshev points varying from seven, nine, and eleven in total. Surface pressure distribution on the blades at the reconstructed at 11 mapped Chebyshev points.
(a)
7 Chebyshev points (b)
9 Chebyshev points (c)
11 Chebyshev points (d)
Reconstructed from the 11 points to eighteenFigure 19
Wall clock time (in seconds) along with the number of interval points.


## 5.1. Hovering Condition

A geometry of the Caradonna and Tung’s rotor is popularly used to validate the hovering condition of the helicopter rotor flow analysis [36]. The rotor system has two blades which have a rectangular-shape planform with no twist and taper in the span direction, and their cross section is NACA0012 airfoil. The span length is 3.75 ft (1.143 m) with an aspect ratio of 7. A collective pitch angle of 8∘ is applied to the blades at all azimuth angles with Mach numbers of 0.439 and 0.877. All computations are carried out in parallel in a Linux cluster using a total of 56 CPUs of i7 Core processors with 24 GB RAM per node, which is effective in handling the increased memory requirement for the mapped Chebyshev pseudospectral method with meshes and solutions at all subtime intervals. Figure 12 shows the multiblock, overset mesh topology used for the Caradonna and Tung’s rotor. A near-body mesh is H-type with 155×239×239 points and the background mesh has the grid of 97×106×217. A total of 11 Chebyshev collocation points, corresponding to 11 subtime intervals due to the periodicity of the rotor motion, are used for the analysis of the hovering flight. Computations of the time-marching method, the HBM, and the mapped Chebyshev pseudospectral method are carried out independently and resultant aerodynamic forces and moments are compared with the wind-tunnel experimental data. For the HBM and the mapped Chebyshev pseudospectral method, solution reconstruction is carried out at the postprocessing stage based on the flow solutions at a finite number of subtime intervals or time instances. A graphical schematic of the flow reconstruction is illustrated in Figure 13 for the helicopter rotor flow analysis. Flow solutions available only at the discrete and unevenly distributed mapped Chebyshev collocation points are interpolated at an arbitrary time instance in the time interval using the following relations. (38)Qjt=∑n=0NanTntj=∑n=0Nancos⁡jπMchn,forj=0,1,2,…,Mch,(39)an=2Ncn∑n=0N1cjQjTntj=2Ncn∑n=0N1cjQjcos⁡jπN,forn=0,1,2,…,N,(40)ci=2,ifi=0,N1,otherwise.A value of M represents the number of points on the interpolant function which will be constructed from the N+1 number of collocation points or subtime intervals. Therefore, Qk is the value from the new approximation through the reconstruction, whereas Qj represents the value at the original subtime intervals. Eqns. (38) and (39) are the same as (1) and (2); however, it should be noted that the values of M and N are different. Figures 14 and 15 show the plots of the surface pressure distribution on the blade at Mach numbers of 0.439 and 0.877, respectively. Spanwise distribution of surface pressure coefficients is plotted at different spanwise locations of 0.5, 0.68, 0.80, 0.89, and 0.96 of the blade radius. Results of the mapped Chebyshev pseudospectral method, HBM, and the time-marching computation show excellent agreements with each other at all spanwise sections, although slight discrepancies are observed in comparison to the wind-tunnel experiments around the region of maximum pressure suction. A total of 19 collocation points are employed for the HBM and 19 points for the Chebyshev pseudospectral method and reconstructed with 180 number of interpolation points. Comparison of thrust and torque, CT and CQ, as well as the figure of merit (FM) is shown in Table 2 where the computations are carried with the Mach number of 0.439. All computation results match very well with those of the wind-tunnel experiments [36]. Surface pressure distributions on two blades at initial nine mapped Chebyshev collocation points or time instances are shown in Figure 16(a). Those can be used for the reconstruction by the interpolation given by (38) and Figure 16(b) shows the blades and surface pressure distribution reconstructed at 36 subtime intervals or time instances.Table 2
Aerodynamic coefficients at Mach number of 0.439.
CT CQ F.M. Experiment 0.0046 Time-accurate Method 0.0051 0.00049 0.53 Harmonic Balance Method 0.0051 0.00049 0.53 Present (Chebyshev spectral Method) 0.0051 0.00049 0.53Figure 12
A multiblock, overset mesh topology around Caradonna and Tung’s rotor blade with the background grid of97×106×217 and the near-body grid of 155×239×239.Figure 13
Schematic for solution reconstruction using the mapped Chebyshev points.Figure 14
Surface pressure coefficients at the radial locations of 0.5, 0.68, 0.8, 0.89, and 0.96 at Mach number of 0.439.
(a)
r / R = 0.50 (b)
r / R = 0.68 (c)
r / R = 0.80 (d)
r / R = 0.89 (e)
r / R = 0.96Figure 15
Surface pressure coefficients at the radial locations of 0.5, 0.68, 0.8, 0.89, and 0.96 at Mach number of 0.877.
(a)
r / R = 0.50 (b)
r / R = 0.68 (c)
r / R = 0.80 (d)
r / R = 0.89 (e)
r / R = 0.96Figure 16
Contour plots of the blade surface pressure at Mach = 0.877.
(a)
Blades and surface pressure contours at 9 mapped Chebyshev collocation points (b)
Blades and surface pressure contours reconstructed at 36 collocation points


## 5.2. Forward Flight

The same Caradonna and Tung rotors are solved at a nonlifting forward flight condition. Blades are not in a pitching or flapping motion and a collective pitch angle is set to be zero. Only forward speed is considered with corresponding tip Mach number of 0.5 and the advance ratio of 0.2. The CPU time is shown in Table3 with respect to the number of Chebyshev collocation points from 11 to 23 and is normalized by the computational time applied with one Chebyshev collocation point. The CPU time is observed with linearly increasing, which can be inferred from the diagonal form of the flux Jacobian matrices shown in (27)~(29) and (31). Table 3 also shows overset grid module time which is defined as the time for hole-cutting and donor cell identification with respect to the number of Chebyshev collocation points. For the time-marching computation for the unsteady analysis with this overset grid technique, the grid module time is about 3.4 seconds per physical time step which means 1° azimuth angle, which is inferred to be more than 20 minutes for one revolution which means 360° azimuth angle. The rotor flow usually converges after five revolutions and it only takes 1.7 hours for overset grid module time. However, by the present mapped Chebyshev pseudospectral method, the grid module time for the hole-cutting and mesh interpolation needs is needed only once at the beginning of the solution procedure and takes less than 2 minutes even with 23 Chebyshev points. Figure 17 shows the plots of the pressure coefficients at the 0.89 radial location of the blade when the azimuth angle varies from 30 deg to 180 deg by 30 deg. Results of the mapped Chebyshev pseudospectral method are computed with a total of 11, 13, and 15 collocation points and compared with the wind-tunnel data [37]. Strong shock appears in the advancing side of the rotor, approximately from 60 deg to 150 deg of the azimuth angle, demonstrating solution accuracy higher with more Chebyshev points. Computation results with a total of 11 points predict relatively well the shock structure in terms of the location and magnitude. Figure 18 shows the contours of surface pressure distribution on the blades at the mapped Chebyshev points which vary from 7, 9, to 11. Figure 18(d) shows the contours of surface pressure distribution on the blades reconstructed from the 11 points to 18 points. We checked the wall clock time with increasing interval numbers of HBM and present method in Figure 19. For time accurate method the dual-time stepping is applied with 360 intervals for one cycle and 20 subiteration for each interval, and three cycles for the periodic convergence of the unsteady flow. The computational time of present method is increasing linearly with respect to the number of intervals as like HBM, and Figure 19 shows present method is as efficient as HBM below 19 intervals.Table 3
Nondimensionalized time with respect to the varying number of Chebyshev collocation points and corresponding grid module time for mesh interpolation and hole cutting.
Chebyshev point CPU timeTn/T1 Grid module time (sec.) 11 7.21 37 15 15.02 51 19 19.81 63 23 25.48 76Figure 17
Pressure coefficients at the 0.89 radial location during the forward flight with an azimuth angle varying from 30 deg to 180 deg by 30 deg.
(a)
ψ = 30 deg (b)
ψ = 60 deg (c)
ψ = 90 deg (d)
ψ = 120 deg (e)
ψ = 150 deg (f)
ψ = 180 degFigure 18
Surface pressure distribution on the blades during the forward flight at the mapped Chebyshev points varying from seven, nine, and eleven in total. Surface pressure distribution on the blades at the reconstructed at 11 mapped Chebyshev points.
(a)
7 Chebyshev points (b)
9 Chebyshev points (c)
11 Chebyshev points (d)
Reconstructed from the 11 points to eighteenFigure 19
Wall clock time (in seconds) along with the number of interval points.


## 6. Conclusions

The mapped Chebyshev pseudospectral method was developed and implemented for URANS solvers with the overset mesh topology. It uses Chebyshev polynomials for solution approximation with the CGL points for the collocation nodes. The time derivative term of the governing equations is replaced by mapped Chebyshev spectral derivative term and discretized in explicit form for time integration in a pseudo time domain. Flux terms are discretized in implicit form and this different discretization strategy makes the integration of the PDE governing equations particularly efficient in both time and space. The governing equations in discrete form are similar to those of the steady RANS equations. Unlike a discrete Fourier series, Chebyshev polynomials have uneven distribution of the collocation points, which often leads to the ill-conditioned spectral derivative matrix and deteriorates the numerical stability of the linear system. An arcsin conformal mapping function is applied for redistributing the Chebyshev points more evenly with optimum value of the control parameters ofα and β to be 0.99. For the validation of the present method, airfoils in oscillation and plunging motions were solved separately to investigate the solution accuracy of the method in periodic and nonperiodic flows. The number of subtime intervals was also varied to demonstrate the sensitivity of the solution accuracy to the number of Chebyshev points. Convergence rate did not change to the varying number of subtime intervals, although computation time and memory to compute one pseudo-time step increase linearly. Only a small number of Chebyshev points were needed to accurately predict the airfoil solution problems. Subsequently, the Chebyshev pseudospectral method was applied to the helicopter rotor flow problems in hovering and nonlifting forward flight conditions using the Caradonna and Tung rotor. As a set of subgrids are constructed once at the beginning of the solution procedure only at given mapped Chebyshev points along with the hole-cutting method for solution interpolation, computation time is saved dramatically in comparison to the time-marching computation which carries out the mesh movement with the hole-cutting at each physical time step which is often limited to less than 10 due to the numerical stability limit. The time saving was on the orders of magnitude. Aerodynamic performance of thrust and torque as well as surface pressure coefficients at several radial locations around azimuth angles shows excellent agreements when compared to those of the time-marching and the HBM method. All computation results show good agreements with wind-tunnel experimental data for both flight conditions. The number of subtime interval more than 11 is needed for the solution accuracy of the forward flight, while the corresponding computation time increases linearly with respect to the number of the subtime intervals. Overall, the present approach of the mapped Chebyshev pseudospectral method is efficient in solving both periodic and nonperiodic problems at greatly reduced cost in comparison to the time-marching computation method, which becomes more effective with the overset topology.



---

*Source: 9745862-2018-05-15.xml*