International Journal of Aerospace Engineering

Volume 2018, Article ID 2508153, 14 pages

https://doi.org/10.1155/2018/2508153

## Prediction of Dynamic Stability Using Mapped Chebyshev Pseudospectral Method

^{1}Aerospace and Ocean Engineering Department, Virginia Tech, Blacksburg, VA 24061, USA^{2}Department of Mechanical Design Engineering, Youngsan University, Yangsan, Republic of Korea

Correspondence should be addressed to Seongim Choi; ude.tv@1iohcs

Received 22 March 2018; Accepted 20 June 2018; Published 1 August 2018

Academic Editor: Mahmut Reyhanoglu

Copyright © 2018 Jae-Young Choi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

#### 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 Section 2 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. where is the vector of conservative variables and is the residual of the spatial discretization of the fluxes.

The time spectral representation of the governing equation can be obtained by the approximation of , , and with discrete Fourier expansions truncated at th harmonic as follows. where 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

The vectors of the Fourier coefficients in (7), (8), and (9) have a size of , and a general form of the system of the equations can be represented as follows: where

However, a direct calculation of (10) is difficult to solve due to its highly nonlinearity that requires 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 where is the Fourier transformation matrix and 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

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 .

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 , and the nonsymmetric transformation in (16) is considered to deliver more flexibility in the distribution of the collocation points over the domain. 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

The final form of the mapped Chebyshev collocation operator becomes

##### 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 by Chebyshev-Gauss-Lobatto point, and it results in the undistribution of collocation points clustered toward the extremes.

An arbitrary and smooth continuous function can be represented with a series of orthogonal polynomials in the domain as follows: where in (21) is the Chebyshev coefficient of the function , and it can be defined in a vector form as

The matrix in (23) includes the Chebyshev polynomials of the function, . The in the matrix represent the degrees of Chebyshev polynomial. The differential of can be represented in a vector form as where the elements of in (24) are the Chebyshev coefficients and corresponds to the Chebyshev polynomials of the . The Chebyshev coefficient can be represented as where the element of the matrix in (26) is defined by the value of . If the value of is an odd number, becomes or turns into 0 otherwise. If the interval is transformed from the interval of to the interval , the differential can be computed by dividing each component of the by . Finally, the time derivatives of an arbitrary function can be computed as

Here, in 27 represents a Chebyshev collocation operator with the matrix size of .

##### 2.5. Implementation

The Chebyshev pseudospectral method with mapping function is applied to the three dimensional unsteady Reynolds Average Navier-Stokes equations as follows: where represents the conservative variables, and , , and denote the residual of viscous and inviscid flux vectors in the , , and direction. These vectors of conservative and flux variables can be approximated with (29), (30), (31), and (32), and the time differential of the in the interval of can be replaced with the mapped Chebyshev collocation operator derived in (28) as where , , , and have the size of matrix and each component is defined at the points which are redistributed by the arcsine mapping function.

The pseudotime stepping is applied to (33).

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].

The implicit integration of the spatial flux used the linearization technique as follows:

The , , and in (37), (38), and (39) represent the flux Jacobian matrices with the size of . Throughout the replacements of flux terms in (36) with (37), (38), and (39), (36) becomes

The flux terms in (37), (38), and (39) are the function of the variables at the specific time instance over the domain of , and only diagonal elements of the matrix , , and are considered as below.

Each of the diagonal element is at different time instances and denoted with subscript zero through 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 and is not required. Finally, (28) can be represented as (42), employing a diagonally implicit method. where

Here, the vector contains the Chebyshev spectral source terms and flux terms as shown in (42). The 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 or [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 or 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. where is the amplitude at each wave number in frequency contents of the time-accurate method and 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: where , , , , , and are the static derivatives and , , and 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, , , and , 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]. where

Herein, 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 as shown in (51) [22].

Herein, and 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. where the vector of and represent the flow solutions and the Fourier coefficients of the corresponding flow solutions in the frequency domain, respectively. Here, the first Fourier coefficients, and , 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 frequency , initial angle of attack , and pitching amplitude . The airfoil oscillates in a range of . The unsteady flow solutions are approximated at Mach number 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 Wilcox-Durbin + turbulence model [23].