# Prediction of Dynamic Stability Using Mapped Chebyshev Pseudospectral Method

Authors: Jae-Young Choi; Dong Kyun Im; Jangho Park; Seongim Choi
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

A mapped Chebyshev pseudospectral method is extended to solve three-dimensional unsteady flow problems. As the classical Chebyshev spectral approach can lead to numerical instabilities due to ill conditioning of the spectral matrix, the Chebyshev points are evenly redistributed over the domain by an inverse sine mapping function. The mapped Chebyshev pseudospectral method can be used as an alternative time-spectral approach that uses a Chebyshev collocation operator to approximate the time derivative terms in the unsteady flow governing equations, and the method can make general applications to both nonperiodic and periodic problems. In this study, the mapped Chebyshev pseudospectral method is employed to solve three-dimensional periodic problem to verify the spectral accuracy and computational efficiency with those of the Fourier pseudospectral method and the time-accurate method. The results show a good agreement with both of the Fourier pseudospectral method and the time-accurate method. The flow solutions also demonstrate a good agreement with the experimental data. Similar to the Fourier pseudospectral method, the mapped Chebyshev pseudospectral method approximates the unsteady flow solutions with a precise accuracy at a considerably effective computational cost compared to the conventional time-accurate method.


---

## Body

## 1. Introduction

The conventional time-marching approach and time-spectral method have been widely used to solve various types of unsteady flow problems that are governed by partial derivative equations. These methods have been developed to approximate the solutions of PDEs by different approaches. To estimate the time derivative term in the unsteady flow governing equations, the time-marching method employs a finite difference method, while time-spectral methods used in this study apply a discrete Fourier expansion and Chebyshev polynomials.In this study, two different basis functions are used for the spectral method in order to solve time-dependent problems. The first kind is based on a frequency domain method by casting a discrete Fourier expansions. Hall et al. [1] first proposed an application of the frequency-domain-based approach to a time-dependent nonlinear flow problem. The method utilizes the Fourier spectral method to the temporal discretization and approximates the time derivative terms in the unsteady transport equations. The transformation to frequency domain from time-domain turns the unsteady form of the governing equations into the time-independent equations by treating time derivative term as a source term, and this allows ones to solve unsteady state problems as steady-state manner at multiple time instance points [2]. The frequency domain-based time-spectral method can deliver very accurate approximation on periodic unsteady problems at an efficient computational cost while it satisfies Nyquist theorem. However, the computational efficiency will become rather slower with an increasing number of harmonics without knowledge of frequency information prior [3].If this is the case and the problem is nonperiodic, the governing equation can be represented with a series of Chebyshev polynomials. The Chebyshev polynomials are one class of the orthogonal Lagrangian polynomials, and the main features of the use of the Chebyshev polynomials are its general application capability to either periodic or nonperiodic problems compared to the frequency-based time-spectral method [3]. The classical Chebyshev polynomials can provide remarkable solution resolution near the boundaries [3], but this feature of the Chebyshev polynomials is not favorable for current unsteady problems for its high potential outcome of the numerical instability. To mitigate such numerical instability, an application of conformal mappings [4–9] to the traditional Chebyshev pseudospectral method is considered. Im et al. [2] have applied conformal mapping techniques to unsteady flow problems, and an inverse sine mapping function is employed for even distribution of the time instance points over the range for current problem.The mapped Chebyshev pseudospectral method can provide the solution accuracy as fine as that of time-accurate method at an effective efficiency in computational cost. Im et al. [10] solved two-dimensional airfoil flows under forced oscillation by the Euler and Navier-Stokes solvers with the mapped Chebyshev pseudospectral method. The results obtained very similar characteristics compared to those of the Fourier pseudospectral method in terms of both solution accuracy and computational efficiency.In this paper, a mapped Chebyshev pseudospectral method is expended to three-dimensional case and employed to approximate solutions of the time-dependent flow past both two- and three-dimensional bodies under a forced harmonic oscillation. The results are used to calculate dynamic derivatives of the oscillating bodies. The standard NACA 0012 and a tailless aircraft, the Lockheed Martin Tactical Aircraft Systems-Innovative Control Effector (ICE) configuration, [11–13] are subjected to the unsteady flow simulations. The solutions are obtained by computation of Euler and Navier-Stokes equations using time-accurate method and time-spectral methods.The details of Chebyshev pseudospectral method with a mapping function and the calculation method of dynamic derivatives are explored in Section2 and Section 3. The solutions from the present method are compared with those from the time-accurate method, the Fourier pseudospectral method, and experimental data [12] to validate the accuracy of the present method in Section 4. Finally, dynamic derivatives are obtained through investigations of the unsteady flow solutions from each of computational methods in Section 5.


## 2. Time-Spectral Methods

### 2.1. Fourier Pseudospectral Method

Frequency domain approach can provide very accurate solutions for the time-dependent transport equations of conservation of mass, momentum, and energy with a temporal periodicity at an efficient computational cost, compared to the traditional time-accurate method. The frequency-based time-spectral method uses a discrete Fourier expansions to estimate the temporal derivative terms in the unsteady flow governing equations.(1)It=∂U∂t+Ru=0,(2)U=ρρuρvρE,where U is the vector of conservative variables and R is the residual of the spatial discretization of the fluxes.The time spectral representation of the governing equation can be obtained by the approximation ofI, U, and R with discrete Fourier expansions truncated at NHth harmonic as follows.
(3)It≈I^0+∑n=1NHI^ancosωnt+I^bnsinωnt,(4)Ut≈U^0+∑n=1NHU^ancosωnt+U^bnsinωnt,(5)∂Ut∂t≈∑n=1NHωnU^ansinωnt−ωnU^bncosωnt,(6)Rt≈R^0+∑n=1NHR^ancosωnt+R^bnsinωnt,where T is the period and ω is the angular frequency. If (3), (4), (5), and (6) are compared to each other, one can find the relationships between Fourier coefficients, and the resulting equations can be represented as
(7)I^0=R^0=0,(8)I^an=ωnU^an+R^an=0,(9)I^bn=−ωnU^bn+R^bn=0.The vectors of the Fourier coefficients in (7), (8), and (9) have a size of NT=NH+1, and a general form of the system of the equations can be represented as follows:
(10)ωMU^+R^=0,where
(11)U^=U0Uc1Uc2⋮UcNHUs1Us2⋮UsNH,R^=R0Rc1Rc2⋮RcNHRs1Rs2⋮RsNH.However, a direct calculation of (10) is difficult to solve due to its highly nonlinearity that requires N3 order of operations [1]. The computational expenses increase rapidly with the number of harmonics. Hall et al. [1] suggested an alternative way to resolve such problems by casting inverse Fourier transformation to (10) as shown in
(12)U^=FU,R^=FR,(13)ωF−1MFU+F−1FR=0,where F is the Fourier transformation matrix and F−1 is inverse Fourier transformation matrix. The pseudotime derivative term is added to (13) to solve with the pseudo time-stepping, and the final form of the governing equations becomes
(14)∂U∂τ+ωDHBMU^+R^=0.The temporal derivative term is treated as a source term, and the equation becomes a time-independent function.


### 2.2. Chebyshev Pseudospectral Method

With an absence of the frequency term, the Chevyshev pseudospectral method can be used for both periodic and nonperiodic problems. The Chebyshev pseudospectral method is based on a series of Chebyshev polynomials to approximate the solution as an arbitrary function. The Chebyshev pseudospectral method can exhibit general applications to unsteady flow problems, compared to the frequency pseudospectral method. However, the method itself inherently suffers from clustering of the collocation points near the boundaries. This uneven distribution of the collocation points can easily lead to numerical instabilities due to the ill conditioning of the spectral derivative matrix [10]. A need for conformal mapping arose to mitigate the numerical instabilities, and the conformal mapping functions have been investigated by a number of researchers. Im et al. [10] considered several types of mapping functions in his study of mapped Chebyshev pseudospectral method. He considered four different types of conformal mapping functions: inverse sine mapping function [9, 14], polynomial mapping function [4], tangent function [5–7], and hyperbolic sine function [8]. His study on conformal mapping showed the inverse sine function which represented the most favorable performance in terms of even distribution and convergence characteristics. The detail of the other mapping function study can be found in [9, 10, 14].


### 2.3. Inverse Sine Function

The standard Chebyshev pseudospectral method uses a cosine function to provide a distribution of collocation points over the interval−1,1.
(15)yi=cosiπN.In this study, inverse sine mapping function is used for even distribution of the collocation points as it showed the most favorable performance, compared to the other mapping functions [10]. The inverse mapping function, proposed by Kosloff and Tal-Ezer [9, 14], can be derived into either symmetric or nonsymmetric transformation forms over the interval −1,1, and the nonsymmetric transformation in (16) is considered to deliver more flexibility in the distribution of the collocation points over the domain.
(16)gy,α,β=arcsinαβ+α−βα+β,(17)α=cosiπN,β=cosiπN,where α and β are the parameters that control the distribution of the collocation points near the boundaries [10, 14]. In this study, the values of α and β are used below 0.995. The time derivative term of the Chebyshev collocation operator is used with the inverse sine mapping function and defined as
(18)dudt=dydtdudy=1g′dudy,g′=yi=Tmaxαβ1−αy1+βy.The final form of the mapped Chebyshev collocation operator becomes(19)dudt=1g′dudy=αβTmax1−αy1+βyDch.


### 2.4. Chebyshev Collocation Operator

Chebyshev pseudospectral method uses a series of orthogonal polynomials to approximate the state and flux variables in the governing equation at each collocation point. However, the collocation points are projected onto the interval of−1≤t≤1 by Chebyshev-Gauss-Lobatto point, and it results in the undistribution of collocation points clustered toward the extremes.An arbitrary and smooth continuous functionut can be represented with a series of orthogonal polynomials in the domain −1≤t≤1 as follows:
(20)ujt=∑n=0NanTntj=∑n=0NancosjπnN,forj=0,1,2,…,N,(21)an=2Ncn∑j=0N1cjujTntj=2Ncn∑j=0NujcjcosjπnN,forn=0,1,2,…,N,(22)ci=2,ifi=0,N,1,otherwise,where an in (21) is the Chebyshev coefficient of the function ut, and it can be defined in a vector form as
(23)a→=DTu→=2NT0t04T0t12…T0tN4T1t02T1t1…T1tN2⋮⋮⋱⋮TNt04TNt12…TNtN4u0u1⋮uN.The matrixDT in (23) includes the Chebyshev polynomials of the function, u→. The TN in the matrix represent the Nth degrees of Chebyshev polynomial. The differential of u→ can be represented in a vector form as
(24)u→′=D~Tb→=T0t0T1t0…TNt0T0t1T1t1…TNt1⋮⋮⋱⋮T0tNT1tN…TNtNb0b1⋮bN,where the elements of b→ in (24) are the Chebyshev coefficients and D~T corresponds to the Chebyshev polynomials of the u′t. The Chebyshev coefficient u→′ can be represented as
(25)b→=Ga→=01030…G0,N00408…G1,N00060…G2,N⋮⋮⋮⋮⋮⋱⋮00000…GN,Na0a1a2⋮aN,(26)ci=2,ifi=0,N,1,otherwise,Gij=2jci,ifi+jodd,0,otherwise,where the element Gi,j of the matrix G in (26) is defined by the value of i+j. If the value of i+j is an odd number, Gij becomes 2j/ci or turns into 0 otherwise. If the interval is transformed from the interval of −1≤t≤1 to the interval 0≤t≤Tmax, the differential can be computed by dividing each component of the Gij by Tmax/2. Finally, the time derivatives of an arbitrary function ut can be computed as
(27)u′=D~Tb→=D~TGa→=D~TGDTu→=Dchu→.Here,Dch in 27 represents a Chebyshev collocation operator with the matrix size of N+1,N+1.


### 2.5. Implementation

The Chebyshev pseudospectral method with mapping function is applied to the three dimensional unsteady Reynolds Average Navier-Stokes equations as follows:(28)∂Q→∂t+∂F→∂x+∂G→∂y+∂H→∂z=0→,(29)Q→=∑n=0NQnTntj,(30)F→=∑n=0NFnTntj,(31)G→=∑n=0NGnTntj,(32)H→=∑n=0NHnTntj,where Q represents the conservative variables, and F, G, and H denote the residual of viscous and inviscid flux vectors in the X, Y, and Z direction. These vectors of conservative and flux variables can be approximated with (29), (30), (31), and (32), and the time differential of the Q in the interval of 0≤t≤Tmax can be replaced with the mapped Chebyshev collocation operator derived in (28) as
(33)1g′DchQ→+∂F→∂x+∂G→∂y+∂H→∂z=0→,where Q→, F→, G→, and H→ have the size of N+1,1 matrix and each component is defined at the points which are redistributed by the arcsine mapping function.
(34)Q→=Q0Q1Q2⋮QN,F→=F0F1F2⋮FN,G→=G0G1G2⋮GN,H→=H0H1H2⋮HN.The pseudotime stepping is applied to (33).
(35)∂Q→∂τ+1g′DchQ→+∂F→∂x+∂G→∂y+∂H→∂z=0→.The diagonalized alternate direction implicit (D-ADI) technique [15] was applied to the Chebyshev pseudospectral method for the time integration. In this study, explicit discretization of the Chebyshev pseudospectral operator and implicit discretization of residual term are taken account into the governing equation with the pseudotime stepping to take advantages of the implicit method for the larger time step size and the advantage of the explicit method for diagonal dominance of the Jacobian matrix as shown in (36). The explicit application of the Chebyshev operator removes its influence on the nondiagonal elements of the matrix, and the source term can be directly added into the flux term at the pseudotime level [10].
(36)∂Q→∂τ+1g′DchQ→n+∂G→n+1∂x+∂H→n+1∂y=0→.The implicit integration of the spatial flux used the linearization technique as follows:(37)F→n+1=F→n+UΔQ→n,(38)G→n+1=G→n+VΔQ→n,(39)H→n+1=H→n+WΔQ→n.TheU, V, and W in (37), (38), and (39) represent the flux Jacobian matrices with the size of N+1,N+1. Throughout the replacements of flux terms in (36) with (37), (38), and (39), (36) becomes
(40)ΔQnΔτ+∂∂xUΔQ→n+∂∂yVΔQ→n+∂∂zWΔQ→n=−∂F→n∂x−∂G→n∂y−∂H→n∂z−1g′DchQ→n.The flux terms in (37), (38), and (39) are the function of the variables at the specific time instance over the domain of 0≤t≤Tmax, and only diagonal elements of the matrix U, V, and W are considered as below.
(41)U=∂U0∂Q00…00∂U1∂Q1…0⋮⋱⋮0…0∂GN∂QN,V=∂G0∂Q00…00∂G1∂Q1…0⋮⋱⋮0…0∂GN∂QN,W=∂H0∂Q00…00∂H1∂Q1…0⋮⋱⋮0…0∂HN∂QN.Each of the diagonal element is at different time instances and denoted with subscript zero throughN and dependent only on the flow variables in that time instance. The key point here is that the diagonal elements of the simplified diagonal Jacobian matrix are identical to those of the flux Jacobian matrix in the time-domain solution algorithm. Therefore, computation of additional terms in Jacobian matrices of V and W is not required. Finally, (28) can be represented as (42), employing a diagonally implicit method.
(42)I+δxVΔτ+δyWΔτΔQ→n=−∂U→n∂xΔτ−∂G→n∂YΔτ−∂H→n∂ZΔτ−Δτg′DchQ→n,(43)T00…00T1…0⋮⋮⋱⋮00…TNΔQ0ΔQ1⋮ΔQN=S0S1⋮SN,where
(44)Ti=I+δxUiiΔτ+δxViiΔτ+δyWiiΔτ,(45)Si=−∂U→in∂xΔτ−∂G→in∂yΔτ−∂H→in∂zΔτ−Δτg′DchQ→in.Here, the vectorS contains the Chebyshev spectral source terms and flux terms as shown in (42). The Tn in (45) can be computed by the method proposed by Pulliam and Chaussee [15]. Therefore, the final form of the governing equations is represented such that the mapped Chebyshev collocation operator is treated as a source term in the governing equation in a similar way to the harmonic source term of the frequency domain method [2, 16].


### 2.6. Computational Efficiency and Solution Accuracy

The main advantage of the time-spectral method is the spectral accuracy and its convergence rate. The time-spectral method can achieve an order of accuracy equivalent to that of the time-accurate method at a considerably effective computational cost. Although the time-spectral solutions can be corrupted by aliasing errors if the number of the time sampling points is less than doubles of the Nyquist frequency or3/2−rule [17], one can find a number of sampling points with which the aliasing error is not significant to the time-spectral solution accuracy through frequency spectrum analysis.A precise quantification of the aliasing errors is important to find an optimal number of sampling points. It is ideal if the number of sampling points is chosen by the Nyquist theorem or3/2−rule since it can deliver exactly the same order of accuracy compared to that of the time-accurate method. However, the computational cost highly depends on the number of sampling points, and this will be lagging off the convergence rate of the time-spectral method with a larger number of sampling points. The total computational cost will exceed that of the time-accurate method. Thus, the trade-off between the solution accuracy and the computational cost should be considered.The number of time instance points is determined by the magnitude of the aliasing error. The aliasing error is quantified by calculating the norm of the discrepancy of the amplitude in the frequency contents between the time-accurate and time-spectral methods as follows.(46)L2,aliasing=∑i=1NCi,TA−Ci,TS2N,where Ci,TA is the amplitude at each wave number in frequency contents of the time-accurate method and Ci,TS is that of the time-spectral method. In this study, the time-spectral solutions are obtained with the number of time instance points to find an appropriate number that maintains the spectral accuracy at an effective computational efficiency.


## 2.1. Fourier Pseudospectral Method

Frequency domain approach can provide very accurate solutions for the time-dependent transport equations of conservation of mass, momentum, and energy with a temporal periodicity at an efficient computational cost, compared to the traditional time-accurate method. The frequency-based time-spectral method uses a discrete Fourier expansions to estimate the temporal derivative terms in the unsteady flow governing equations.(1)It=∂U∂t+Ru=0,(2)U=ρρuρvρE,where U is the vector of conservative variables and R is the residual of the spatial discretization of the fluxes.The time spectral representation of the governing equation can be obtained by the approximation ofI, U, and R with discrete Fourier expansions truncated at NHth harmonic as follows.
(3)It≈I^0+∑n=1NHI^ancosωnt+I^bnsinωnt,(4)Ut≈U^0+∑n=1NHU^ancosωnt+U^bnsinωnt,(5)∂Ut∂t≈∑n=1NHωnU^ansinωnt−ωnU^bncosωnt,(6)Rt≈R^0+∑n=1NHR^ancosωnt+R^bnsinωnt,where T is the period and ω is the angular frequency. If (3), (4), (5), and (6) are compared to each other, one can find the relationships between Fourier coefficients, and the resulting equations can be represented as
(7)I^0=R^0=0,(8)I^an=ωnU^an+R^an=0,(9)I^bn=−ωnU^bn+R^bn=0.The vectors of the Fourier coefficients in (7), (8), and (9) have a size of NT=NH+1, and a general form of the system of the equations can be represented as follows:
(10)ωMU^+R^=0,where
(11)U^=U0Uc1Uc2⋮UcNHUs1Us2⋮UsNH,R^=R0Rc1Rc2⋮RcNHRs1Rs2⋮RsNH.However, a direct calculation of (10) is difficult to solve due to its highly nonlinearity that requires N3 order of operations [1]. The computational expenses increase rapidly with the number of harmonics. Hall et al. [1] suggested an alternative way to resolve such problems by casting inverse Fourier transformation to (10) as shown in
(12)U^=FU,R^=FR,(13)ωF−1MFU+F−1FR=0,where F is the Fourier transformation matrix and F−1 is inverse Fourier transformation matrix. The pseudotime derivative term is added to (13) to solve with the pseudo time-stepping, and the final form of the governing equations becomes
(14)∂U∂τ+ωDHBMU^+R^=0.The temporal derivative term is treated as a source term, and the equation becomes a time-independent function.


## 2.2. Chebyshev Pseudospectral Method

With an absence of the frequency term, the Chevyshev pseudospectral method can be used for both periodic and nonperiodic problems. The Chebyshev pseudospectral method is based on a series of Chebyshev polynomials to approximate the solution as an arbitrary function. The Chebyshev pseudospectral method can exhibit general applications to unsteady flow problems, compared to the frequency pseudospectral method. However, the method itself inherently suffers from clustering of the collocation points near the boundaries. This uneven distribution of the collocation points can easily lead to numerical instabilities due to the ill conditioning of the spectral derivative matrix [10]. A need for conformal mapping arose to mitigate the numerical instabilities, and the conformal mapping functions have been investigated by a number of researchers. Im et al. [10] considered several types of mapping functions in his study of mapped Chebyshev pseudospectral method. He considered four different types of conformal mapping functions: inverse sine mapping function [9, 14], polynomial mapping function [4], tangent function [5–7], and hyperbolic sine function [8]. His study on conformal mapping showed the inverse sine function which represented the most favorable performance in terms of even distribution and convergence characteristics. The detail of the other mapping function study can be found in [9, 10, 14].


## 2.3. Inverse Sine Function

The standard Chebyshev pseudospectral method uses a cosine function to provide a distribution of collocation points over the interval−1,1.
(15)yi=cosiπN.In this study, inverse sine mapping function is used for even distribution of the collocation points as it showed the most favorable performance, compared to the other mapping functions [10]. The inverse mapping function, proposed by Kosloff and Tal-Ezer [9, 14], can be derived into either symmetric or nonsymmetric transformation forms over the interval −1,1, and the nonsymmetric transformation in (16) is considered to deliver more flexibility in the distribution of the collocation points over the domain.
(16)gy,α,β=arcsinαβ+α−βα+β,(17)α=cosiπN,β=cosiπN,where α and β are the parameters that control the distribution of the collocation points near the boundaries [10, 14]. In this study, the values of α and β are used below 0.995. The time derivative term of the Chebyshev collocation operator is used with the inverse sine mapping function and defined as
(18)dudt=dydtdudy=1g′dudy,g′=yi=Tmaxαβ1−αy1+βy.The final form of the mapped Chebyshev collocation operator becomes(19)dudt=1g′dudy=αβTmax1−αy1+βyDch.


## 2.4. Chebyshev Collocation Operator

Chebyshev pseudospectral method uses a series of orthogonal polynomials to approximate the state and flux variables in the governing equation at each collocation point. However, the collocation points are projected onto the interval of−1≤t≤1 by Chebyshev-Gauss-Lobatto point, and it results in the undistribution of collocation points clustered toward the extremes.An arbitrary and smooth continuous functionut can be represented with a series of orthogonal polynomials in the domain −1≤t≤1 as follows:
(20)ujt=∑n=0NanTntj=∑n=0NancosjπnN,forj=0,1,2,…,N,(21)an=2Ncn∑j=0N1cjujTntj=2Ncn∑j=0NujcjcosjπnN,forn=0,1,2,…,N,(22)ci=2,ifi=0,N,1,otherwise,where an in (21) is the Chebyshev coefficient of the function ut, and it can be defined in a vector form as
(23)a→=DTu→=2NT0t04T0t12…T0tN4T1t02T1t1…T1tN2⋮⋮⋱⋮TNt04TNt12…TNtN4u0u1⋮uN.The matrixDT in (23) includes the Chebyshev polynomials of the function, u→. The TN in the matrix represent the Nth degrees of Chebyshev polynomial. The differential of u→ can be represented in a vector form as
(24)u→′=D~Tb→=T0t0T1t0…TNt0T0t1T1t1…TNt1⋮⋮⋱⋮T0tNT1tN…TNtNb0b1⋮bN,where the elements of b→ in (24) are the Chebyshev coefficients and D~T corresponds to the Chebyshev polynomials of the u′t. The Chebyshev coefficient u→′ can be represented as
(25)b→=Ga→=01030…G0,N00408…G1,N00060…G2,N⋮⋮⋮⋮⋮⋱⋮00000…GN,Na0a1a2⋮aN,(26)ci=2,ifi=0,N,1,otherwise,Gij=2jci,ifi+jodd,0,otherwise,where the element Gi,j of the matrix G in (26) is defined by the value of i+j. If the value of i+j is an odd number, Gij becomes 2j/ci or turns into 0 otherwise. If the interval is transformed from the interval of −1≤t≤1 to the interval 0≤t≤Tmax, the differential can be computed by dividing each component of the Gij by Tmax/2. Finally, the time derivatives of an arbitrary function ut can be computed as
(27)u′=D~Tb→=D~TGa→=D~TGDTu→=Dchu→.Here,Dch in 27 represents a Chebyshev collocation operator with the matrix size of N+1,N+1.


## 2.5. Implementation

The Chebyshev pseudospectral method with mapping function is applied to the three dimensional unsteady Reynolds Average Navier-Stokes equations as follows:(28)∂Q→∂t+∂F→∂x+∂G→∂y+∂H→∂z=0→,(29)Q→=∑n=0NQnTntj,(30)F→=∑n=0NFnTntj,(31)G→=∑n=0NGnTntj,(32)H→=∑n=0NHnTntj,where Q represents the conservative variables, and F, G, and H denote the residual of viscous and inviscid flux vectors in the X, Y, and Z direction. These vectors of conservative and flux variables can be approximated with (29), (30), (31), and (32), and the time differential of the Q in the interval of 0≤t≤Tmax can be replaced with the mapped Chebyshev collocation operator derived in (28) as
(33)1g′DchQ→+∂F→∂x+∂G→∂y+∂H→∂z=0→,where Q→, F→, G→, and H→ have the size of N+1,1 matrix and each component is defined at the points which are redistributed by the arcsine mapping function.
(34)Q→=Q0Q1Q2⋮QN,F→=F0F1F2⋮FN,G→=G0G1G2⋮GN,H→=H0H1H2⋮HN.The pseudotime stepping is applied to (33).
(35)∂Q→∂τ+1g′DchQ→+∂F→∂x+∂G→∂y+∂H→∂z=0→.The diagonalized alternate direction implicit (D-ADI) technique [15] was applied to the Chebyshev pseudospectral method for the time integration. In this study, explicit discretization of the Chebyshev pseudospectral operator and implicit discretization of residual term are taken account into the governing equation with the pseudotime stepping to take advantages of the implicit method for the larger time step size and the advantage of the explicit method for diagonal dominance of the Jacobian matrix as shown in (36). The explicit application of the Chebyshev operator removes its influence on the nondiagonal elements of the matrix, and the source term can be directly added into the flux term at the pseudotime level [10].
(36)∂Q→∂τ+1g′DchQ→n+∂G→n+1∂x+∂H→n+1∂y=0→.The implicit integration of the spatial flux used the linearization technique as follows:(37)F→n+1=F→n+UΔQ→n,(38)G→n+1=G→n+VΔQ→n,(39)H→n+1=H→n+WΔQ→n.TheU, V, and W in (37), (38), and (39) represent the flux Jacobian matrices with the size of N+1,N+1. Throughout the replacements of flux terms in (36) with (37), (38), and (39), (36) becomes
(40)ΔQnΔτ+∂∂xUΔQ→n+∂∂yVΔQ→n+∂∂zWΔQ→n=−∂F→n∂x−∂G→n∂y−∂H→n∂z−1g′DchQ→n.The flux terms in (37), (38), and (39) are the function of the variables at the specific time instance over the domain of 0≤t≤Tmax, and only diagonal elements of the matrix U, V, and W are considered as below.
(41)U=∂U0∂Q00…00∂U1∂Q1…0⋮⋱⋮0…0∂GN∂QN,V=∂G0∂Q00…00∂G1∂Q1…0⋮⋱⋮0…0∂GN∂QN,W=∂H0∂Q00…00∂H1∂Q1…0⋮⋱⋮0…0∂HN∂QN.Each of the diagonal element is at different time instances and denoted with subscript zero throughN and dependent only on the flow variables in that time instance. The key point here is that the diagonal elements of the simplified diagonal Jacobian matrix are identical to those of the flux Jacobian matrix in the time-domain solution algorithm. Therefore, computation of additional terms in Jacobian matrices of V and W is not required. Finally, (28) can be represented as (42), employing a diagonally implicit method.
(42)I+δxVΔτ+δyWΔτΔQ→n=−∂U→n∂xΔτ−∂G→n∂YΔτ−∂H→n∂ZΔτ−Δτg′DchQ→n,(43)T00…00T1…0⋮⋮⋱⋮00…TNΔQ0ΔQ1⋮ΔQN=S0S1⋮SN,where
(44)Ti=I+δxUiiΔτ+δxViiΔτ+δyWiiΔτ,(45)Si=−∂U→in∂xΔτ−∂G→in∂yΔτ−∂H→in∂zΔτ−Δτg′DchQ→in.Here, the vectorS contains the Chebyshev spectral source terms and flux terms as shown in (42). The Tn in (45) can be computed by the method proposed by Pulliam and Chaussee [15]. Therefore, the final form of the governing equations is represented such that the mapped Chebyshev collocation operator is treated as a source term in the governing equation in a similar way to the harmonic source term of the frequency domain method [2, 16].


## 2.6. Computational Efficiency and Solution Accuracy

The main advantage of the time-spectral method is the spectral accuracy and its convergence rate. The time-spectral method can achieve an order of accuracy equivalent to that of the time-accurate method at a considerably effective computational cost. Although the time-spectral solutions can be corrupted by aliasing errors if the number of the time sampling points is less than doubles of the Nyquist frequency or3/2−rule [17], one can find a number of sampling points with which the aliasing error is not significant to the time-spectral solution accuracy through frequency spectrum analysis.A precise quantification of the aliasing errors is important to find an optimal number of sampling points. It is ideal if the number of sampling points is chosen by the Nyquist theorem or3/2−rule since it can deliver exactly the same order of accuracy compared to that of the time-accurate method. However, the computational cost highly depends on the number of sampling points, and this will be lagging off the convergence rate of the time-spectral method with a larger number of sampling points. The total computational cost will exceed that of the time-accurate method. Thus, the trade-off between the solution accuracy and the computational cost should be considered.The number of time instance points is determined by the magnitude of the aliasing error. The aliasing error is quantified by calculating the norm of the discrepancy of the amplitude in the frequency contents between the time-accurate and time-spectral methods as follows.(46)L2,aliasing=∑i=1NCi,TA−Ci,TS2N,where Ci,TA is the amplitude at each wave number in frequency contents of the time-accurate method and Ci,TS is that of the time-spectral method. In this study, the time-spectral solutions are obtained with the number of time instance points to find an appropriate number that maintains the spectral accuracy at an effective computational efficiency.


## 3. Dynamic Derivative Analysis

In a previous study on stability analysis on tailless aircraft [18], stability and control behavior of a tailless aircraft are investigated for Innovative Control Effectors (ICE) configuration developed by Lockheed Martin. The system matrix of the ICE configuration is developed as follows:
(47)u̇ẇq̇α̇=η∂CA∂u+μu0CA0η∂CA∂w+μw0CA0η∂CA∂q−w0η∂CA∂α̇+gcosθ0η∂CN∂u+μu0CN0η∂CN∂w+μw0CN0η∂CN∂q−w0η∂CN∂α̇+gcosθ0η∂CA∂u+μu0CA0η∂Cm∂w+μw0Cm0η∂Cm∂q−w0η∂Cm∂α̇+gcosθ00010ΔuΔwΔqΔα,where ∂CA/∂u, ∂CN/∂u, ∂Cm/∂u, ∂CA/∂w, ∂CN/∂w, and ∂Cm/∂w are the static derivatives and ∂CA/∂q, ∂CN/∂q, and ∂Cm/∂q are the steady dynamic derivatives. These derivative can be easily found using steady CFD computations and finite difference method with rotating flow scheme [19]. However, the unsteady dynamic derivative terms, ∂CA/∂α̇, ∂CN/∂α̇, and ∂Cm/∂α̇, require time-marching unsteady CFD computations.The dynamic derivatives can be approximated based on the unsteady flow solutions, in which a forced sinusoidal oscillation is imposed at the center of gravity of aircraft [20]. From one cycle of unsteadiness during the harmonic oscillation, the mathematical model for the increment in aerodynamic forces and moments can be expressed as follows with an assumption of a linear relationship between the aerodynamic properties and flight states [21].
(48)ΔCk=∂Ck∂αΔα+c2V∞∂Ck∂α̇α̇+c2V∞∂Ck∂qq,where
(49)Δαt=αAsinωt,(50)Δα̇t=q=ωαAcosωt.Herein,αA is the amplitude of angle of attack and ω is an oscillation frequency. Equation (48) can be rewritten with (49) and (50) based on an assumption α̇t=q as shown in (51) [22].
(51)ΔCαt=αAkCkα̇+Ckqcosωt+αACkαsinωt,(52)ΔCαt=CCcosωt+CSsinωt.Herein,Ckα and Ckα̇+Ckq̇ represent stiffness (in-phase) and damping (out-of-phase) components of dynamic derivatives, respectively [21]. The determination of the dynamic derivatives can be achieved by comparison of the first harmonic flow solutions from the time-accurate and time-spectral methods. Prior to the comparisons, it was assured that the unsteady flow solutions from each method are fully converged. The solutions from each method are transformed into the frequency domain through discrete Fourier transform as shown below.
(53)CC0CC1⋮CCNHCS0CS1⋮CSNH=cos2π0×1/NT2cos2π0×2/NT2…cos2π0×NT/NT2cos2π1×1NTcos2π1×2NT…cos2π1×NTNT⋮⋮⋱⋮cos2πNH×1NTcosNHπN×2NT…cos2πNH×NTNTsin2π0×1NTsin2π0×2/NT2…sin2π0×NTNTsin2π1×1NTsin2π1×2NT…sin2π1×NTNT⋮⋮⋱⋮sin2πNH×1NTsinNHπN×2NT…sin2πNH×NTNTX0X1⋮XNT,where the vector of X′s and C′s represent the flow solutions and the Fourier coefficients of the corresponding flow solutions in the frequency domain, respectively. Here, the first Fourier coefficients, CC1 and SS1, are compared with the coefficients in (52).


## 4. Validation of Chebyshev Pseudospectral Method and Dynamic Derivatives

### 4.1. NACA 0012

A harmonic oscillation is imposed on the quarter chord of NACA 0012 airfoil to validate the accuracy of Chebyshev pseudospectral method. The computation of Navier-Stokes equations is performed with reduced frequencyk=0.0814, initial angle of attack α0=0.016, and pitching amplitude αA=2.51. The airfoil oscillates in a range of −2.494<α<2.526. The unsteady flow solutions are approximated at Mach number M=0.755 using the time-accurate method and the time-spectral method. A C-type structured mesh is employed for the NACA 0012 airfoil with 4096 points and 4260 elements as shown in Figure 1. The Chebyshev collocation operator with the mapping function is applied to the time differential term along with D-ADI method for time integration and periodic boundary conditions. The spatial numerical flux in each direction is discretized using Roe’s FDS with a third order MUSCL scheme to achieve higher order accuracy, and turbulent eddy viscosity is computed from the two-equation k−ω Wilcox-Durbin + turbulence model [23].Figure 1
Structured mesh for NACA 0012 airfoil.Unsteady flow solutions of oscillating NACA0012 are approximated by the time-accurate method (TA), frequency pseudospectral method (FPM), and Chebyshev pseudospectral method (CPM). In computation with the time-spectral methods, the spectral accuracy can be contaminated by aliasing errors due to a lack of number of sampling points. In this study, the aliasing errors are represented by the norm of the discrepancy in the amplitude at every wave number between the time-accurate method and the time-spectral method results. The results are plotted in Figure2(a). The aliasing error is quantified for a different number of the time instance points, and Figure 2(b) shows that the aliasing errors are reduced with an increasing number of the time instance points. This indicates the aliasing errors become less effective on the solution reconstruction as the number of the time instance points increases, and one can find an optimum number, in terms of computational cost, which provides the spectral accuracy fairly equivalent to the order of accuracy of the original wave form. As shown in Figure 2(b), the aliasing error converges at the time instance points of 15, and the time-spectral methods give a precise solution accuracy with respect to the time-accurate solutions. However, the computational time linearly increases with a number of collocation points over the interval as shown in Figure 3.Figure 2
Frequency analysis for the time-accurate method and the time-spectral method results.
(a)
Frequency contents of the mapped Chebyshev pseudospectral method for NACA0012 under pitching motion
(b)
L2 norms of difference in amplitude of frequency contents between the time-accurate method and the mapped Chebyshev pseudospectral method for NACA 0012 under pitching motionFigure 3
Computational time for NACA 0012 under pitching motion.The advantages in the computational efficiency vanish when the number of the time instance points is more than 50. In this study, the number of the time instance points of 15 is chosen as an optimal for the Fourier- and mapped Chebyshev pseudospectral method, and it gives twice faster computational efficiency compared to the time-accurate method.Normal force and pitching moment coefficients are compared with those from time-accurate method and the frequency pseudospectral method as shown in Figure4. The solid line in black represents the time-accurate method, and dotted and dashed lines in blue correspond to the present method, while the solid line in red represents the frequency pseudospectral method. The experimental data are shown in black circles [11]. With a small number of collocation points, the results show a large deviation from those of the time-accurate method and the frequency pseudospectral method results with 15 time instance points. However, the results of the present method are shown in a good agreement with the solutions from the other methods as the time instance points increases. With 7 or more number of time instance points, the results from Chebyshev pseudospectral method converge to the time-accurate method results. The time-spectral method results are also shown in reasonably good agreement with the experimental data.Figure 4
Normal force and pitching moment hysteresis of NACA 0012.
(a)
Normal force coefficient
(b)
Pitching moment coefficientDynamic derivatives are obtained for oscillating NACA 0012 from the unsteady flow solutions from the time-accurate method, the frequency pseudospectral method, and the mapped Chebyshev pseudospectral method. The number of collocation points varied from 5 to 15 for the Chebyshev pseudospectral method, and the frequency pseudospectral method is used with 15 time instance points. The results are summarized in Tables1 and 2. The dynamic derivatives from normal coefficients are presented in Table 1. With a fewer number of collocation points, the values of the dynamic derivatives show some discrepancy compared to the time-accurate method results. However, they are converging to the results from the time-accurate method and the frequency pseudospectral method as the number of collocation points increases.Table 1
Dynamic derivatives for NACA 0012 airfoil from normal force.
Stiffness (in-phase)
Damping (out-of-phase)
Time-accurate
−7.16
28.1
FPM (15)
−7.18
29.1
CPM (4)
−8.16
22.3
CPM (6)
−7.22
29.0
CPM (8)
−7.20
29.1
CPM (14)
−7.17
29.7Table 2
Dynamic derivatives for NACA 0012 airfoil from pitching moment.
Stiffness (in-phase)
Damping (out-of-phase)
Time-accurate
0.12
−3.96
FPM (15)
0.11
−3.79
CPM (4)
0.13
−3.09
CPM (6)
0.11
−3.60
CPM (8)
0.11
−3.75
CPM (14)
0.11
−3.82


### 4.2. ICE Configuration

The mapped Chebyshev pseudospectal method is applied to an investigation of the oscillating three-dimensional body. The ICE configuration is used to evaluate the accuracy of the present method. The unsteady flow solutions are obtained by computation of Euler equations with the mapped Chebyshev- and the Fourier collocation operator. A structured mesh is employed with 5.2 million elements within 9 blocks as shown in Figure5. A forced harmonic oscillation is imposed on the aircraft center of gravity with reduced frequency of k=0.1242, with an initial angle of attack of α0=0.0, and the amplitude of 4.16 and 8.33 degrees. The Mach number is given as 0.0266 which is the same as the wind tunnel test conditions.Figure 5
Structured mesh for ICE configuration.The experiment was performed with the Reynolds number of0.574×106 in Subsonic Air Force Research Laboratory Vertical Wind Tunnel (AFRL) at Wright-Patterson AFB, Ohio [12]. The ICE configuration was tested under various pitching rate conditions during the longitudinal stability test as presented in Table 3. The magnitude of amplitude varied from 4.16 and 8.33 degrees to match pitch rate qc/2V of 0.009 and 0.018 at 0 angle of attack [12].Table 3
Wind tunnel test conditions for forced oscillation.
Properties
Values
Freestream velocity
8.84 m/s
Reynolds number
0.574 × 106/m
Amplitude
4.16 and 8.33 deg.
Reduced frequency
0.1242
Moment reference
38% MACThe number of the time instance points varied to investigate the spectral accuracy of the time-spectral methods. The frequency contents of the unsteady solutions from the Fourier pseudospectral method and the mapped Chebyshev pseudospectral method are represented in Figure6(a). Similar to the pitching NACA 0012 case in the previous section, the discrepancy between the time-accurate method and the mapped Chebyshev pseudospectral method decreases as the number of the time instance points increases. As shown in Figure 6(b), the minimum number of the time instance points should be chosen more than five to minimize effects of the aliasing error on the solution accuracy. The reconstructed solutions from the present method are shown in Figure 7 through dotted and dashed lines in blue. The Fourier pseudospectral method solutions are represented in red solid line, and the time-accurate method is shown in black solid line. The reconstructed solutions from the Fourier pseudospectral method and the mapped pseudospectral method are shown in a good agreement with that of time-accurate method.Figure 6
Frequency analysis for the time-accurate method and the time-spectral method results.
(a)
Frequency contents for ICE configuration
(b)
L2 norms of difference in amplitude of frequency contents between the time-accurate method and the mapped Chebyshev pseudospectral method for ICE configuration under pitching motionFigure 7
Normal force and pitching moment hysteresis of ICE configuration.
(a)
Normal force coefficient
(b)
Pitching moment coefficientThe computational cost is also investigated and compared between the time-accurate and the time-spectral methods as shown in Figure8. The total computational time for the time-accurate method is about 333 hours. On the other hand, it linearly increases with the number of the time instance points for the mapped Chebyshev pseudospectral method and the Fourier pseudospectral method. The computational cost becomes the same order of magnitude when the number of time instance points is more than 25. From the frequency content analysis, the time-spectral methods with the time instance points more than eleven can provide the order of solution accuracy equivalent to that of the time-accurate method. This indicates that the mapped Chebyshev pseudospectral method can compute the solution at half cost of the time-accurate method for the current application on ICE geometry.Figure 8
Computational time for ICE configuration under pitching motion.To validate the time-spectral methods for the three-dimensional flow, the unsteady flow solutions are compared with experimental data as shown in Figure9. The circle symbols in black and solid line with diamond symbols in red indicate the experimental data and the Fourier pseudospectral method results, respectively. The solid line in blue corresponds to the results from the present method. Figure 9 shows a good agreement in normal force coefficients with the experimental data. The pitching moment, however, presents a large deviation at high pitching rate. The issues can be found in an inaccurate prediction of the flow separation caused by artificial dissipation term in Euler solver. This can become a significant concern for the calculation of the complex vortex structures on the delta wing, which are formed from the nose and the leading edge region. As shown in Figure 10(b), the pressure contours are dynamically changing as the angle of attack, but it is hard to find an evidence of the interactions between the vortex structure and the surface of the body, which would come out as flow separations and reattachments as it goes downstream.Figure 9
Comparison of normal force and pitching moment coefficient with experimental data.
(a)
Normal force coefficient
(b)
Pitching moment coefficientFigure 10
Surface pressure contours of ICE configuration at different collocation points.
(a)
Collocation points
(b)
Surface pressure contourThe calculation of dynamic derivatives is extended to three-dimensional case. Unsteady flow solutions of oscillating ICE configuration are obtained using the Fourier pseudospectral method with 11 time instance points and the mapped Chebyshev pseudospectral method with 5 to 11 time instance points. As shown in Tables4 and 5, the discrepancy in dynamic derivatives between the current method and time-accuracy method is decreasing with a larger number of the time instance points. The out-of-phase and in-phase components of dynamic derivatives show a good agreement with the Fourier pseudospectral method results with 11 time instance points and the time-accurate method.Table 4
Dynamic derivatives for ICE configuration from normal force.
Stiffness (in-phase)
Damping (out-of-phase)
Time-accurate
−5.65
−2.87
FPM
−5.79
−2.92
CPM (5)
−5.39
−2.80
CPM (7)
−5.79
−2.85
CPM (11)
−5.78
−2.94Table 5
Dynamic derivatives for ICE configuration from pitching moment.
Stiffness (in-phase)
Damping (out-of-phase)
Time-accurate
1.98
1.02
FPM (11)
2.12
1.03
CPM (5)
2.04
0.99
CPM (7)
2.02
1.01
CPM (11)
2.11
1.04


## 4.1. NACA 0012

A harmonic oscillation is imposed on the quarter chord of NACA 0012 airfoil to validate the accuracy of Chebyshev pseudospectral method. The computation of Navier-Stokes equations is performed with reduced frequencyk=0.0814, initial angle of attack α0=0.016, and pitching amplitude αA=2.51. The airfoil oscillates in a range of −2.494<α<2.526. The unsteady flow solutions are approximated at Mach number M=0.755 using the time-accurate method and the time-spectral method. A C-type structured mesh is employed for the NACA 0012 airfoil with 4096 points and 4260 elements as shown in Figure 1. The Chebyshev collocation operator with the mapping function is applied to the time differential term along with D-ADI method for time integration and periodic boundary conditions. The spatial numerical flux in each direction is discretized using Roe’s FDS with a third order MUSCL scheme to achieve higher order accuracy, and turbulent eddy viscosity is computed from the two-equation k−ω Wilcox-Durbin + turbulence model [23].Figure 1
Structured mesh for NACA 0012 airfoil.Unsteady flow solutions of oscillating NACA0012 are approximated by the time-accurate method (TA), frequency pseudospectral method (FPM), and Chebyshev pseudospectral method (CPM). In computation with the time-spectral methods, the spectral accuracy can be contaminated by aliasing errors due to a lack of number of sampling points. In this study, the aliasing errors are represented by the norm of the discrepancy in the amplitude at every wave number between the time-accurate method and the time-spectral method results. The results are plotted in Figure2(a). The aliasing error is quantified for a different number of the time instance points, and Figure 2(b) shows that the aliasing errors are reduced with an increasing number of the time instance points. This indicates the aliasing errors become less effective on the solution reconstruction as the number of the time instance points increases, and one can find an optimum number, in terms of computational cost, which provides the spectral accuracy fairly equivalent to the order of accuracy of the original wave form. As shown in Figure 2(b), the aliasing error converges at the time instance points of 15, and the time-spectral methods give a precise solution accuracy with respect to the time-accurate solutions. However, the computational time linearly increases with a number of collocation points over the interval as shown in Figure 3.Figure 2
Frequency analysis for the time-accurate method and the time-spectral method results.
(a)
Frequency contents of the mapped Chebyshev pseudospectral method for NACA0012 under pitching motion
(b)
L2 norms of difference in amplitude of frequency contents between the time-accurate method and the mapped Chebyshev pseudospectral method for NACA 0012 under pitching motionFigure 3
Computational time for NACA 0012 under pitching motion.The advantages in the computational efficiency vanish when the number of the time instance points is more than 50. In this study, the number of the time instance points of 15 is chosen as an optimal for the Fourier- and mapped Chebyshev pseudospectral method, and it gives twice faster computational efficiency compared to the time-accurate method.Normal force and pitching moment coefficients are compared with those from time-accurate method and the frequency pseudospectral method as shown in Figure4. The solid line in black represents the time-accurate method, and dotted and dashed lines in blue correspond to the present method, while the solid line in red represents the frequency pseudospectral method. The experimental data are shown in black circles [11]. With a small number of collocation points, the results show a large deviation from those of the time-accurate method and the frequency pseudospectral method results with 15 time instance points. However, the results of the present method are shown in a good agreement with the solutions from the other methods as the time instance points increases. With 7 or more number of time instance points, the results from Chebyshev pseudospectral method converge to the time-accurate method results. The time-spectral method results are also shown in reasonably good agreement with the experimental data.Figure 4
Normal force and pitching moment hysteresis of NACA 0012.
(a)
Normal force coefficient
(b)
Pitching moment coefficientDynamic derivatives are obtained for oscillating NACA 0012 from the unsteady flow solutions from the time-accurate method, the frequency pseudospectral method, and the mapped Chebyshev pseudospectral method. The number of collocation points varied from 5 to 15 for the Chebyshev pseudospectral method, and the frequency pseudospectral method is used with 15 time instance points. The results are summarized in Tables1 and 2. The dynamic derivatives from normal coefficients are presented in Table 1. With a fewer number of collocation points, the values of the dynamic derivatives show some discrepancy compared to the time-accurate method results. However, they are converging to the results from the time-accurate method and the frequency pseudospectral method as the number of collocation points increases.Table 1
Dynamic derivatives for NACA 0012 airfoil from normal force.
Stiffness (in-phase)
Damping (out-of-phase)
Time-accurate
−7.16
28.1
FPM (15)
−7.18
29.1
CPM (4)
−8.16
22.3
CPM (6)
−7.22
29.0
CPM (8)
−7.20
29.1
CPM (14)
−7.17
29.7Table 2
Dynamic derivatives for NACA 0012 airfoil from pitching moment.
Stiffness (in-phase)
Damping (out-of-phase)
Time-accurate
0.12
−3.96
FPM (15)
0.11
−3.79
CPM (4)
0.13
−3.09
CPM (6)
0.11
−3.60
CPM (8)
0.11
−3.75
CPM (14)
0.11
−3.82


## 4.2. ICE Configuration

The mapped Chebyshev pseudospectal method is applied to an investigation of the oscillating three-dimensional body. The ICE configuration is used to evaluate the accuracy of the present method. The unsteady flow solutions are obtained by computation of Euler equations with the mapped Chebyshev- and the Fourier collocation operator. A structured mesh is employed with 5.2 million elements within 9 blocks as shown in Figure5. A forced harmonic oscillation is imposed on the aircraft center of gravity with reduced frequency of k=0.1242, with an initial angle of attack of α0=0.0, and the amplitude of 4.16 and 8.33 degrees. The Mach number is given as 0.0266 which is the same as the wind tunnel test conditions.Figure 5
Structured mesh for ICE configuration.The experiment was performed with the Reynolds number of0.574×106 in Subsonic Air Force Research Laboratory Vertical Wind Tunnel (AFRL) at Wright-Patterson AFB, Ohio [12]. The ICE configuration was tested under various pitching rate conditions during the longitudinal stability test as presented in Table 3. The magnitude of amplitude varied from 4.16 and 8.33 degrees to match pitch rate qc/2V of 0.009 and 0.018 at 0 angle of attack [12].Table 3
Wind tunnel test conditions for forced oscillation.
Properties
Values
Freestream velocity
8.84 m/s
Reynolds number
0.574 × 106/m
Amplitude
4.16 and 8.33 deg.
Reduced frequency
0.1242
Moment reference
38% MACThe number of the time instance points varied to investigate the spectral accuracy of the time-spectral methods. The frequency contents of the unsteady solutions from the Fourier pseudospectral method and the mapped Chebyshev pseudospectral method are represented in Figure6(a). Similar to the pitching NACA 0012 case in the previous section, the discrepancy between the time-accurate method and the mapped Chebyshev pseudospectral method decreases as the number of the time instance points increases. As shown in Figure 6(b), the minimum number of the time instance points should be chosen more than five to minimize effects of the aliasing error on the solution accuracy. The reconstructed solutions from the present method are shown in Figure 7 through dotted and dashed lines in blue. The Fourier pseudospectral method solutions are represented in red solid line, and the time-accurate method is shown in black solid line. The reconstructed solutions from the Fourier pseudospectral method and the mapped pseudospectral method are shown in a good agreement with that of time-accurate method.Figure 6
Frequency analysis for the time-accurate method and the time-spectral method results.
(a)
Frequency contents for ICE configuration
(b)
L2 norms of difference in amplitude of frequency contents between the time-accurate method and the mapped Chebyshev pseudospectral method for ICE configuration under pitching motionFigure 7
Normal force and pitching moment hysteresis of ICE configuration.
(a)
Normal force coefficient
(b)
Pitching moment coefficientThe computational cost is also investigated and compared between the time-accurate and the time-spectral methods as shown in Figure8. The total computational time for the time-accurate method is about 333 hours. On the other hand, it linearly increases with the number of the time instance points for the mapped Chebyshev pseudospectral method and the Fourier pseudospectral method. The computational cost becomes the same order of magnitude when the number of time instance points is more than 25. From the frequency content analysis, the time-spectral methods with the time instance points more than eleven can provide the order of solution accuracy equivalent to that of the time-accurate method. This indicates that the mapped Chebyshev pseudospectral method can compute the solution at half cost of the time-accurate method for the current application on ICE geometry.Figure 8
Computational time for ICE configuration under pitching motion.To validate the time-spectral methods for the three-dimensional flow, the unsteady flow solutions are compared with experimental data as shown in Figure9. The circle symbols in black and solid line with diamond symbols in red indicate the experimental data and the Fourier pseudospectral method results, respectively. The solid line in blue corresponds to the results from the present method. Figure 9 shows a good agreement in normal force coefficients with the experimental data. The pitching moment, however, presents a large deviation at high pitching rate. The issues can be found in an inaccurate prediction of the flow separation caused by artificial dissipation term in Euler solver. This can become a significant concern for the calculation of the complex vortex structures on the delta wing, which are formed from the nose and the leading edge region. As shown in Figure 10(b), the pressure contours are dynamically changing as the angle of attack, but it is hard to find an evidence of the interactions between the vortex structure and the surface of the body, which would come out as flow separations and reattachments as it goes downstream.Figure 9
Comparison of normal force and pitching moment coefficient with experimental data.
(a)
Normal force coefficient
(b)
Pitching moment coefficientFigure 10
Surface pressure contours of ICE configuration at different collocation points.
(a)
Collocation points
(b)
Surface pressure contourThe calculation of dynamic derivatives is extended to three-dimensional case. Unsteady flow solutions of oscillating ICE configuration are obtained using the Fourier pseudospectral method with 11 time instance points and the mapped Chebyshev pseudospectral method with 5 to 11 time instance points. As shown in Tables4 and 5, the discrepancy in dynamic derivatives between the current method and time-accuracy method is decreasing with a larger number of the time instance points. The out-of-phase and in-phase components of dynamic derivatives show a good agreement with the Fourier pseudospectral method results with 11 time instance points and the time-accurate method.Table 4
Dynamic derivatives for ICE configuration from normal force.
Stiffness (in-phase)
Damping (out-of-phase)
Time-accurate
−5.65
−2.87
FPM
−5.79
−2.92
CPM (5)
−5.39
−2.80
CPM (7)
−5.79
−2.85
CPM (11)
−5.78
−2.94Table 5
Dynamic derivatives for ICE configuration from pitching moment.
Stiffness (in-phase)
Damping (out-of-phase)
Time-accurate
1.98
1.02
FPM (11)
2.12
1.03
CPM (5)
2.04
0.99
CPM (7)
2.02
1.01
CPM (11)
2.11
1.04


## 5. Conclusions

A mapped Chebyshev pseudospectral method is extended for three-dimensional unsteady flow problems. While the spectral method can find unsteady flow solutions with an accuracy equivalent to the time-accurate method, the frequency-based pseudospectral approach has been limited to periodic problems. However, the Chebyshev pseudospectral method can be an alternative time-spectral approach with its application to both periodic and nonperiodic problems. In light of such potential, the mapped Chebyshev pseudospectral method is employed to solve three-dimensional unsteady periodic flow in order to validate its spectral solution accuracy and computational efficiency with those of Fourier pseudospectral method and the traditional time-accurate method. The unsteady flow solutions of the flow past NACA 0012 and ICE configuration under the forced pitching oscillation are obtained by the time-accurate method, the Fourier pseudospectral method, and the mapped Chebyshev pseudospectral method. The results from the present method are compared with the experimental data for validations as well. The mapped Chebyshev pseudospectral method found unsteady flow solutions in a good agreement with the other methods for both two- and three-dimensional cases. The present method also provided a good agreement in normal force coefficient but showed a discrepancy in some degree in pitching moment coefficient, compared to the experimental data. Finally, dynamic derivatives are obtained based on the unsteady flow solutions from the time-accurate, the Fourier pseudospectral method, and the mapped Chebyshev pseudospectral method. With 11 time instance points, the values of dynamic derivatives based on the mapped Chebyshev pseudospectral method results are obtained similar to those of the time-accurate method and the Fourier pseudospectral method at the same number of time instance points. The mapped Chebyshev pseudospectral method showed computational cost equivalent to the Fourier pseudospectral method, and about 2.5 and 2 times faster computational efficiency for NACA 0012 and ICE configuration case, respectively, compared to that of time-accurate method. In general, the mapped Chebyshev pseudospectral method approximated the unsteady flow solutions for the periodic oscillating body with a reasonable accuracy compared to the experimental data and provided a good agreement with those from the traditional time-accurate method and the frequency pseudospectral as well.



---

*Source: 2508153-2018-08-01.xml*