Transonic flutter characteristics of an airfoil with morphing devices

Abstract

An investigation into transonic flutter characteristic of an airfoil conceived with the morphing leading and trailing edges has been carried out. Computational fluid dynamics (CFD) is used to calculate the unsteady aerodynamic force in transonic flow. An aerodynamic reduced order model (ROM) based on autoregressive model with exogenous input (ARX) is used in the numerical simulation. The flutter solution is determined by eigenvalue analysis at specific Mach number. The approach is validated by comparing the transonic flutter characteristics of the Isogai wing with relevant literatures before applied to a morphing airfoil. The study reveals that by employing the morphing trailing edge, the shock wave forms and shifts to the trailing edge at a lower Mach number, and aerodynamic force stabilization happens earlier. Meanwhile, the minimum flutter speed increases and transonic dip occurs at a lower Mach number. It is also noted that leading edge morphing has negligible effect on the appearance of the shock wave and transonic flutter. The mechanism of improving the transonic flutter characteristics by morphing technology is discussed by correlating shock wave location on airfoil surface, unsteady aerodynamics with flutter solution.

Introduction

Morphing wing technology has been developed to improve aerodynamic efficiency and flight performance by changing the wing shape adaptively during flight. ¹ It is regarded as one of the potential and feasible approach contributing to next generation green transport aircraft. ² However, the coupling between the aerodynamics and structures of a morphing wing may have more significant influence to the aeroelastic characteristics especially in transonic regime. The drop of transonic flutter speed, the so-called ‘transonic dip’, ³ , ⁴ is a major concern for a conventional transport aircraft. For a morphing wing, this aeroelastic phenomenon raises a new challenge to be addressed.

In recent years, many efforts have been made to assess the aeroelastic stability for different types of morphing concepts. The flutter mechanism of a span-morphing wing was analysed by Huang et al., ⁵ revealing the rigid-body modes especially pitch motion of the aircraft have a significant effect on flutter characteristics. The flutter and divergence characteristics of a large civil aircraft with morphing winglets and adaptive flap tabs were systemically examined to assess the robustness and safety of adoption in the real structure. ⁶ Hu et al. ⁷ paid more attention to the nonlinear aeroelastic characteristics of a folding wing with cubic connection stiffness in the quasi-steady condition and during the morphing process, and variable motion types, such as chaos motion, were observed in the numerical calculation. Flutter and divergent speeds significantly changed during the transition between take-off, climb, cruise and loiter phases of a fully morphing wing studied by Ünlüsoy and Yaman. ⁸ In particular, unsteady aerodynamic analysis of a morphing wing conducted by Kan et al. ⁹ and Xiang et al. ¹⁰ showed that the stall can be delayed by implementing a flexible periodical trailing-edge deflection. These results suggest that applying morphing technology on a conventional wing may have a favourable effect on the flutter characteristic. Flutter suppression by active control of morphing wing was also studied. Based on a low-fidelity aeroelastic model and cantilever uniform morphing wings, the feasibility of active flutter suppression by using morphing device was studied by Ajaj and Friswell. ¹¹ Investigation was also extended into the dynamic behaviour and stability of an axially morphing wing in supersonic airflow ¹² to demonstrate that a proposed morphing law is effective for flutter suppression.

From a practical point of view so far, morphing technology is more suitable for small unmanned vehicles at low flight speed. ¹³ Hence most of the study on aeroelastic characteristics of morphing wing has been focused in subsonic air flow. Few investigation has been found into the effect of morphing wing on the transonic flutter despite the transonic flutter is a critical concern for large commercial aircraft. ¹⁴ , ¹⁵

Generally speaking, transonic flutter is a nonlinear aeroelastic phenomenon due to shock wave acting on the wing. In transonic regime, the steady flow parameters vary with spatial position in the flow field around the wing that generates unsteady aerodynamic forces induced by the wing motion. To simplify the analysis, it is assumed that the parameters of flow field and the shock wave motion vary in a linear fashion with the wing motion of small perturbation in transonic steady air flow. This is usually called dynamically linear but statically nonlinear aerodynamics. ¹⁶

Wing transonic flutter analysis is a particularly difficult task but also an attractive academic topic for aeroelasticians due to the complexity in aerodynamics modelling. The coupled CFD/CSD (Computational Structural Dynamics) approach, or time marching approach based on CFD, has long been used to obtain the wing transonic flutter characteristic since the 1990s. ¹⁷ , ¹⁸ However, this approach is time-consuming and requires much computation cost even to perform only a single solution of wing transonic flutter. To overcome the disadvantage of the coupled CFD/CSD time marching method, ROM approaches,^19–21 based on CFD technique are developed to calculate the wing flutter speed with the assumption of dynamically linear aerodynamics. ¹⁶ Amongst these approaches, the system identification method ¹⁹ , ²² is a robust and effective technique to build the transonic aerodynamic ROM, which is used in the present study.

A morphing wing model with the leading edge and trailing edge devices from previous study ²³ , ²⁴ was employed in this paper. Since the transonic flutter is more sensitive to airfoil profile than structural property, the mass and stiffness distribution was kept constant for different morphing configurations in the present study. This setting enables us to address the effect of the aerodynamic shape of a morphing wing on the transonic flutter characteristics. The results indicate that the morphing trailing edge device can improve the transonic flutter characteristic, but the morphing leading edge has negligible influence. The mechanism of improving the transonic flutter characteristics by morphing technology is investigated by altering the angle of attack of the airfoil due to the aerodynamic similarity of them.

In this paper, the governing equations of motion and technical methods are formulated for flutter analysis of the morphing wing in the next section. The validation study of the methods using the classical Isogai wing is performed. Subsequently, the numerical results and discussion of the transonic flutter characteristic of the morphing wing are presented. Conclusions are drawn in the last section.

Aeroelastic airfoil model with morphing leading and trailing edges

The governing equation of a wing section motion

A typical wing section modelled as an aeroelastic airfoil system with plunging ( $h$ ) and pitching ( $\alpha$ ) motion is illustrated in Figure 1. The elastic axis of the airfoil (E point) is located at a distance of $ab$ after the mid-chord point; the gravity centre (G point) is located at $x_{\alpha }b$ after the elastic axis, where $b$ is the semi-chord length, $m$ is the mass per unit span, $S_{\alpha }=x_{\alpha }b$ is the first moment of inertia about the elastic axis, and $I_{\alpha }=m{r}_{\alpha }^2$ is the moment of inertia about the elastic axis. The bending stiffness and torsion stiffness of coefficient $K_h=m{\omega }_h^2$ and $K_{\alpha }=I_{\alpha }{\omega }_{\alpha }^2$ are modelled by springs attached to the elastic axis.

Figure 1.

A typical two-dimensional aeroelastic airfoil system.

The governing equations of motion of such a 2-D system were derived from the Lagrange equations according to Dowell et al., ¹⁶ written as

$$\left\{\begin{array}{l}m\ddot{h}+S_{\alpha }\ddot{\alpha }+K_{h}h=-L\\ S_{\alpha }\ddot{h}+I_{\alpha }\ddot{\alpha }+K_{\alpha }\alpha =M_{e\alpha }\end{array}\right.$$

(1)

where $L=1/2\rho V^2{c}_l$ and $M_{e\alpha }=1/2\rho V^2{c}_m$ are the aerodynamic lift and moment about the elastic axis, respectively. $c_l$ is the lift coefficient, $c_m$ is the aerodynamic moment coefficient, and $\rho$ is the air density.

By introducing non-dimensional time $\tau ={\omega }_{\alpha }t$ and mass ratio $\mu =m/\pi \rho b^2$ , equation (1) can be rewritten as

$$\left\{\begin{array}{l}h″+x_{\alpha }\alpha ″+{\left(\frac{{\omega }_h}{{\omega }_{\alpha }}\right)}^{2}h=\frac{U^2}{\pi \mu }\left(-c_l\right)\\ x_{\alpha }h″+r_{\alpha }^2\alpha ″+r_{\alpha }^2\alpha =\frac{U^2}{\pi \mu }\left(2c_m\right)\end{array}\right.$$

(2)

where ${ \left(\right)}^{″}=d^2 \left(\right)/d{\tau }^2$ , $U=V/b{\omega }_{\alpha }$ is the non-dimensional air speed. Subsequently, the governing equation can be written in matrix form,

$$\mathit{M}{\mathit{\xi }}^{″}+\mathit{K}\mathit{\xi }=\frac{U^2}{\pi \mu }{\mathit{f}}_{\mathit{a}}$$

(3)

where $\mathit{M}=\left[\begin{array}{cc}1& x_{\alpha }\\ x_{\alpha }& r_{\alpha }^2\end{array}\right]$ is the mass matrix, $\mathit{K}=\left[\begin{array}{cc}{({\omega }_h/{\omega }_{\alpha })}^2& 0\\ 0& 1\end{array}\right]$ is the stiffness matrix. For this aeroelastic system, $\mathit{\xi }={\left\{\begin{array}{cc}h/b& \alpha \end{array}\right\}}^T$ and ${\mathit{f}}_{\mathit{\alpha }}={\left\{-\begin{array}{cc}{c}_l& {2c}_m\end{array}\right\}}^T$ serve as the generalized displacements and the generalized aerodynamic forces, respectively. The generalized aerodynamic forces corresponding to the generalized displacements in transonic air flow can be obtained from full CFD simulation or an aerodynamic ROM.

By defining the structural state vector ${\mathit{x}}_{\mathit{s}}={\left\{\begin{array}{cc}\mathit{\xi }& {\mathit{\xi }}^{\prime }\end{array}\right\}}^T$ , the structural motion equation of the aeroelastic system can be written as

$$\left\{\begin{array}{l}{\mathit{x}}_{\mathit{s}}^{\prime }={\mathit{A}}_{\mathit{s}}{\mathit{x}}_{\mathit{s}}+\frac{U^2}{\pi \mu }{\mathit{B}}_{\mathit{s}}{\mathit{f}}_{\mathit{a}}\\ \mathit{\xi }={\mathit{C}}_{\mathit{s}}{\mathit{x}}_{\mathit{s}}+{\mathit{D}}_{\mathit{s}}{\mathit{f}}_{\mathit{a}}\end{array}\right.$$

(4)

where ${\left(\mathrm{ }\right)}^{\mathrm{\prime }}=d\left(\right)/d\tau$ , ${\mathit{A}}_{\mathit{s}}=\left[\begin{array}{cc}0& \mathit{I}\\ -{\mathit{M}}^{-1}\mathit{K}& 0\end{array}\right]$ , ${\mathit{B}}_{\mathit{s}}=\left[\begin{array}{c}0\\ {\mathit{M}}^{-1}\end{array}\right]$ , ${\mathit{C}}_{\mathit{s}}=\left[\begin{array}{c}\mathit{I}\\ 0\end{array}\right]$ , and ${\mathit{D}}_{\mathit{s}}=\left[0\right]$ .

Aerodynamic ROM and flutter equation

Amongst the numerous methods, the system identification method is an effective and efficient technique to establish aerodynamic ROM. Following the suggestion from Cowan et al. ¹⁹ and Zhang and Ye, ²² the ARX model is used to establish the ROM of transonic aerodynamics. The time domain equation for Multi Input and Multi Output (MIMO) ARX model can be described as

$${\mathit{f}}_{\mathit{a}}\left(k\right)=\sum _{i=1}^{na}{\mathit{A}}_i{\mathit{f}}_{\mathit{a}}\left(k-i\right)+\sum _{i=0}^{nb-1}{\mathit{B}}_i\mathit{\xi }\left(k-i\right)$$

(5)

One advantage of ARX model is that the system response at any time step ${\mathit{f}}_{\mathit{a}}\left(k\right)$ is a linear combination of past inputs $\mathit{\xi }(k-i)$ and outputs ${\mathit{f}}_{\mathit{a}}(k-i)$ , so that this model is easy to establish the ROM mathematically. With an assumed model order consisting of $na$ past outputs and $nb$ inputs, the only task is to identify the constant coefficient matrices ${\mathit{A}}_i$ and ${\mathit{B}}_i$ .

In the present study, a so-called ‘3211’ signal developed by Cowan et al. ¹⁹ is utilized as the input of the CFD solve due to its broad frequency spectra. Then the least squares method is adopted to fit the time history of the output of the CFD solver, i.e., ${\mathit{f}}_{\mathit{a}}$ , to carry out the unknown coefficient matrices in equation (5).

One challenging problem with using the ARX model to build aerodynamic ROM is how to determine the order of the model. In theory, the order of aerodynamic ROM could vary at different flow conditions, for instance, Mach number, angle of attack even airfoil profile. In the present study, identifying the model order is treated as a minimization problem at specific flow condition, which can be written in a general form

$$\underset{\begin{array}{c}\mathit{na}=1\text{,}2,\dots , 12\\ \mathit{nb}=1\text{,}2,\dots , 12\end{array}}{\min }:J \mathrm{where}\mathrm{ } J=w\left|c_{l,\mathit{CFD}}-c_{l,\mathit{ROM}}\right|+(1-w)\left|c_{m,\mathit{CFD}}-c_{m,\mathit{ROM}}\right|$$

(6)

where $w$ is the weight factor, chosen from 0.2 to 0.4 herein.

Thus, the Genetic Algorithm (GA) can be implemented to search the most appropriate order of ARX. And with well-determined orders and corresponding coefficient matrix, the discrete-time ARX model can be transformed into the continuous-time form through Tustin approximation, which can be described in state-space form as

$$\left\{\begin{array}{l}{\mathit{x}}_{\mathit{a}}^{\prime }={\mathit{A}}_{\mathit{a}}{\mathit{x}}_{\mathit{a}}+{\mathit{B}}_{\mathit{a}}\mathit{\xi }\\ {\mathit{f}}_{\mathit{a}}={\mathit{C}}_{\mathit{a}}{\mathit{x}}_{\mathit{a}}+{\mathit{D}}_{\mathit{a}}\mathit{\xi }\end{array}\right.$$

(7)

where the aerodynamic state vector is

$${\mathit{x}}_{\mathit{a}}={\left\{\begin{array}{ccc}{\mathit{f}}_{\mathit{a}}(\tau -\mathrm{\Delta }\tau )& \dots & {\mathit{f}}_{\mathit{a}}(\tau -na\mathrm{\Delta }\tau )\\ \mathit{\xi }(\tau -\mathrm{\Delta }\tau )& \dots & \mathit{\xi }(\tau -(nb-1\left)\mathrm{\Delta }\tau \right)\end{array}\right\}}^T$$

and ${\mathit{A}}_{\mathit{a}}$ , ${\mathit{B}}_{\mathit{a}}$ , ${\mathit{C}}_{\mathit{a}}$ , ${\mathit{D}}_{\mathit{a}}$ are the coefficient matrices for the state-space model of aerodynamic ROM.

Then, introducing $\mathit{x} = {\left\{\begin{array}{cc}{\mathit{x}}_{\mathit{s}}^T& {\mathit{x}}_{\mathit{a}}^T\end{array}\right\}}^T$ and coupling the structural motion equation (4) with the aerodynamic ROM equation (7), the governing equation for the aeroelastic system in state-space form can be obtained

$${\mathit{x}}^{\prime }=\mathit{Ax}=\left[\begin{array}{cc}{\mathit{A}}_{\mathit{s}}+\frac{U^2}{\pi \mu }{\mathit{B}}_{\mathit{s}}{\mathit{D}}_{\mathit{a}}{\mathit{C}}_{\mathit{s}}& {\frac{U^2}{\pi \mu }\mathit{B}}_{\mathit{s}}{\mathit{C}}_{\mathit{a}}\\ {\mathit{B}}_{\mathit{a}}{\mathit{C}}_{\mathit{s}}& {\mathit{A}}_{\mathit{a}}\end{array}\right]\mathit{x}$$

(8)

Then, the critical non-dimensional flutter speed $U$ and flutter frequency ratio $\omega /{\omega }_{\alpha }$ can be obtained by conventional stability analysis, i.e., solving the eigenvalue of $\mathit{A}$ in equation (8) at different air speeds.

The morphing wing mechanism

The original morphing technology was developed in the DARPA Smart Wing project. ²⁵ , ²⁶ A new morphing mechanism called Eccentric Beam Actuation Mechanism (EBAM) for both leading and trailing edges has been designed, analysed and tested in previous research ²³ , ²⁴ ,^27–29 as shown in Figure 2(a). By rotating the EBAM at one end using an actuator, the wing shape is forced to deform depending on the EBAM rotation angle as illustrated in Figure 2(b) with the curvature of the chord line shown in Figure 2(c) to fulfill a specified morphing shape. The EBAM is connected to the skin through discs and stringers, which provide a pass-way to transfer the actuating force to the skin.

Figure 2.

A morphing wing: (a) leading and trailing edge device, (b) airfoil shape, (c) deformed chord line.

Taking the morphing trailing edge for example, the bending angle $\theta \left(x\right)$ varies from 0 at the rear spar to ${\theta }_{Te}$ at the trailing edge, as displayed in Figure 2(b). In the same way, ${\theta }_{Le}$ is defined as the deflection angle at the leading edge. Note that ${\theta }_{Le}$ and ${\theta }_{Te}$ are used to measure the deformation of morphing wing in the following sections. As presented in Figure 2(c), the horizontal and vertical displacements of the ith point along the morphing wing chord-line can be obtained from the geometric relationship described in equation (9). The position of the deformed upper and lower skin follows the original thickness distribution along the morphed chord-line. As for the upper and lower skins, the displacements are taken at the point with the same x-coordinate on chord-line in the current study.

$$\left\{\begin{array}{l}{u}_{i+1}=u_i+|x_{i+1}-x_i|\sin \left(\frac{{\theta }_{Te}}{L_m}{x}_i\right)\mathrm{t}\mathrm{an}\left(\frac{{\theta }_{Te}}{L_m}{x}_i\right)\\ v_{i+1}=v_i+|x_{i+1}-x_i|\sin \left(\frac{{\theta }_{Te}}{L_m}{x}_i\right)\end{array}\right.$$

(9)

where $u_i$ and $v_i$ are the horizontal and vertical displacement of the ith point with x-coordinate of $x_i$ , and $L_m$ is the length of morphing device. Note that $u_1=0$ and $v_1=0$ .

Isogai wing model

ANSYS Fluent was used to model the aerodynamics based on the fluid governing equations, in which the Euler and Navier-Stokes equations were solved, respectively, by using the pressure-based coupled algorithm in the current paper. In Fluent, a control-volume-based technique is employed to convert the general scalar transport equation to an algebraic equation, which is solved by using a point implicit (Gauss-Seidel) linear equation solver in conjunction with an algebraic multigrid (AMG) method. To deal with viscous flow problems, a classical one-equation turbulence model, the S-A model, is used in current simulations. For spatial discretization, the second-order upwind scheme is implemented to interpolate the convection terms. In terms of temporal discretization, a technique called bounded second order implicit time integration is employed for real-time advancement.

A Radial Basis Functions (RBF) interpolation for large mesh deformation ³⁰ is implemented to enhance the capability of mesh deformation in ANSYS Fluent via user-defined function (UDF). The RBF interpolation $s\left(\mathit{x}\right)$ , representing the displacement of the CFD mesh, ³⁰ , ³¹ can be expressed by a sum of base functions

$$s\left(\mathit{x}\right)=\sum _{i=1}^{n_{bd}}{\mathit{\alpha }}_i\varphi \left(\left|\mathit{x}-{\mathit{x}}_{bi}\right|\right)$$

(10)

where ${\mathit{x}}_{bi}$ is the centre for RBF, describing the displacement of boundary nodes, and $n_{bd}$ is the number of boundary node, | | is the norm biasing.

$$\left|\mathit{x}-{\mathit{x}}_{\mathit{i}}\right|=\sqrt{{(x-x_i)}^2+{(y-y_i)}^2}/R$$

(11)

where $R$ is the support radius, and 10 $c$ is adopted herein. $\varphi$ is the basis function, and Wendland's C2 function is implemented in the current investigation

$$\varphi \left(\left|\mathit{x}\right|\right)=\left\{\begin{array}{c}{\left(1-\left|\mathit{x}\right|\right)}^4(4\left|\mathit{x}\right|+1), \left|\mathit{x}\right|<1\\ 0, \left|\mathit{x}\right|\ge 1\end{array}\right.$$

(12)

When the motion of nodes on the boundary, i.e., ${\mathit{d}}_{\mathit{b}}$ , is specified, coefficients $\mathit{\alpha }={\left\{\begin{array}{ccc}{\alpha }_1& \dots & {\alpha }_{n_{bd}}\end{array}\right\}}^T$ can be obtained by the inverse of equation (10)

$$\mathit{\alpha }={\mathit{d}}_{\mathit{b}}{\mathit{M}}_{b,b}^{-1}$$

(13)

where ${\mathit{M}}_{b,b}$ is a $n_{bd}\times n_{bd}$ matrix containing the evaluations of the basis function ${\varphi }_{b_i,b_j}=\varphi \left({\mathit{x}}_{{\mathit{b}}_{\mathit{i}}}-{\mathit{x}}_{{\mathit{b}}_{\mathit{j}}}\right)$ . The displacement of all remaining nodes can be determined by equation (10).

The Isogai Wing model is of NACA 64A010 airfoil with the following parameters ³ , ⁴

$$a=-2.0, x_{\alpha }=1.8, r_{\alpha }=1.8655, \frac{{\omega }_h}{{\omega }_{\alpha }}=1.0,$$

$$\mu =60.$$

For a specified Mach number, ‘3211’ signal, as shown in Figure 3(a), is used as the inputs of CFD solver, and the corresponding coefficients of lift and moment can be obtained. Subsequently, the aerodynamic ROM is established for this Mach number. Typical time history of the aerodynamic coefficients at Mach 0.8 from CFD simulation and the aerodynamic ROM is shown in Figure 3(b).

Figure 3.

Comparison of aerodynamic force coefficient from CFD solver and ARX ROM for Isogai wing model at Mach 0.8 (Euler solution): (a) ‘3211’ signal, (b) lift and aerodynamic moment coefficients.

Based on the aerodynamic ROM, the root locus for different Mach numbers are used to determinate the critical flutter speed and flutter frequency. After depicting these results versus Mach numbers, the resulting flutter boundary in transonic regime is shown in Figure 4.

Figure 4.

Comparison of flutter boundary of Isogai wing model: flutter speed (a) and flutter frequency (b) versus Mach number from inviscid solution; flutter speed (c) and flutter frequency (d) versus Mach number from viscous solution.

For the Euler solutions as shown in Figure 4(a) and (b), the results are in good agreement with those obtained by transonic flutter solution in the frequency domain ³² and time marching solutions. ¹⁷ , ³³ , ³⁴ The results show that the flutter speeds in transonic regime are lower than in the subsonic regime and there are multiple values of flutter speed occurring between Mach 0.85 and 0.9, which forms the so-called S shape flutter boundary. ³⁵

Only a few studies of Isogai wing model in viscous flow are available in the existing literature. Aerodynamic ROM method are utilized to obtain the flutter solution at different Mach numbers as shown in Figure 4(c) and (d). Compared with the existing time marching solutions, ³³ , ³⁴ reasonably good agreements are achieved demonstrating the feasibility of the present methods. It is also found that the so-called S shape flutter boundary observed in the Euler computational results disappears for Navier-Stoke solutions.

Transonic flutter of a morphing wing

Flutter boundary of Isogai wing with morphing devices

The study was then extended to the Isogai wing model with morphing leading edge (LE) and trailing edge (TE) to investigate the transonic flutter characteristics. To focus on the effect of the aerodynamic profile variation of the morphing wing, the mass and stiffness distribution of the morphing wing keeps in the same values as the original Isogai wing.

The interface location of the LE and TE with the wing box is set at $0.3c$ and $0.6c$ chord of the morphing wing, respectively. Figure 5(a) shows the deformed airfoil profile of morphing TE and LE with different deflection angles forming a smooth aerodynamic surface. The aerodynamic pressure distributions were obtained at different Mach numbers based on Euler equations for the original and the morphing TE as shown in Figure 5(b) and (c). Smooth pressure distribution was also obtained for the morphing LE case despite some fluctuation is observed near the LE as shown in Figure 5(d).

Figure 5.

(a) Airfoil profile for morphing TE and LE and pressure distributions for (b) wing without morphing, (c) morphing TE ${\theta }_{te}=5°$ , (d) morphing LE ${\theta }_{le}=5°$ for Euler solution.

The effect of morphing TE and LE on the transonic characteristic of the wing can be considered via the aerodynamic ROM as well. By setting the morphing TE or LE at a specific deflection angle with zero mean angle of attack, the generalized aerodynamic forces due to ‘3211’ signal for the deformed airfoil shape can be obtained from CFD solver. The aerodynamic ROM can then be established at different Mach numbers. Subsequently, the process for searching the critical flutter condition over a range of flight speed was executed to determine the transonic flutter characteristics of the morphing wing.

Figure 6 shows the transonic flutter boundary of the wing with morphing LE and TE at different deflection angles based on Euler and Navier-Stokes equations. For the Euler results, for $5°$ morphing TE, the flutter speed drops with respect to the original Isogai wing between Mach 0.725 and 0.81, whilst flutter speed increases above Mach 0.82. Conversely, the morphing LE at $2°$ and $5°$ has a beneficial effect on the flutter characteristic in a small Mach number range below 0.85, but exhibits an adverse effect above Mach number 0.85. It is evident that the effect of morphing LE on transonic flutter is much less than the TE.

Figure 6.

Comparison of flutter boundary of Isogai wing with the morphing wing based on: (a) Euler equation, (b) Navier-Stokes equations.

From Figure 6(a), it is apparent that the minimum non-dimensional flutter speed is about 5.45 for the morphing TE model with deflection angle of $5°$ in the transonic region. Theoretically the transonic dip can be further increased to 5.53 with a proper morphing law, i.e. deploying morphing TE from Mach 0.81, increasing by 36.9% compared to original Isogai wing. The results indicate that employing a morphing TE is feasible and effective to improve the transonic flutter characteristic.

Similar effects of morphing LE and TE on the flutter speed are observed for Navier-Stokes solutions as shown in Figure 6(b). Ideally, the minimum flutter speed in transonic regime can be increased from 4.82 of the original model to 8.41 with $5°$ morphing TE. And the effect of morphing LE on bottom of flutter boundary is much less than that of morphing TE.

From linear aerodynamics, ³ , ⁴ the flutter speed generally decreases in transonic regime below Mach 1, and increases significantly above Mach 1. The current results, as shown in Figure 6, present the same trend for both original and morphing wings. To the authors’ knowledge, the exact Mach number, in which the flutter speed starts to increase in the transonic regime, has not yet been fully understood. For that reason, it is difficult to give a general conclusion on flutter characteristics over transonic wings.

According Bendiksen, ³⁶ , ³⁷ the Mach number freeze phenomenon, or transonic stabilization law, plays an important role in the wing transonic flutter. Flutter speed may increase significantly when the freestream Mach number reaches the freeze Mach number. It is defined as the freestream Mach number at which shock waves on both upper and lower surfaces reach the TE of an airfoil. ³⁶ , ³⁷ And these phenomena have been also observed in our previous numerical cases. ³⁸

In order to give physical-based explanations of morphing wing on transonic flutter, the relationship between the shock wave locations on airfoil surface, aerodynamic coefficient and the flutter characteristic is addressed.

To evaluate the unsteady characteristic of the transonic aerodynamics, the aerodynamic influence coefficient ²⁰ is employed in the present study. The input of CFD, i.e. the pitching motion of the aeroelastic airfoil, is taken as a sinusoid function

$$\alpha \left(t\right)= {\alpha }_0\mathrm{s}\mathrm{in}\left(\omega t\right)$$

(14)

where ${\alpha }_0$ denotes the amplitude and $\omega$ is the frequency. The corresponding output of CFD, taking aerodynamic moment coefficient as an example, is periodic which can be written as

$$c_m\left(t\right)={\left(c_m\right)}_c\cos \left(\omega t\right)+{\left(c_m\right)}_s\sin \left(\omega t\right)={\left(c_m\right)}_{cs}\sin \left(\omega t+\phi \right)$$

(15)

As for the coefficients ${\left(c_m\right)}_c$ and ${\left(c_m\right)}_s$ can be obtained by fitting the time history of aerodynamic moment coefficient to the above equation. The aerodynamic influence coefficient can be expressed as

$$c_{m\alpha }=\frac{{\left(c_m\right)}_{cs}}{{\alpha }_0}{e}^{i\phi }=\frac{{\left(c_m\right)}_{cs}}{{\alpha }_0}\left[\cos \left(\phi \right)+i\sin \left(\phi \right)\right]=\mathfrak{R}\left(c_{m\alpha }\right)+i\mathfrak{I}\left(c_{m\alpha }\right)$$

(16)

Then the aerodynamic influence coefficient of different reduced frequencies ( $k=\omega b/V$ ) at different Mach numbers are calculated based on Euler equations, as displayed in Figure 7.

Figure 7.

Aerodynamic influence coefficient due to pitching oscillation around $a=-2$ as a function of Mach number and reduced frequency ( $k=\omega b/V$ ) for: (a) NACA 64A010 airfoil, (b) NACA 64A010 with morphing trailing edge deflection $5^{°}$ .

From Figure 5(b), according to the Euler solution, it can be observed that the freeze Mach number is around 0.9 for the original NACA 64A010 airfoil according to the location of shock wave. When Mach number is greater than freeze Mach number, the unsteady aerodynamic moment coefficient become stable as clearly observed in Figure 7(a). Meanwhile, a strong stabilization of the flutter boundary occurs as shown in Figure 6(a), which coincides with the observation in Ref. 36.

Following the Mach number freeze phenomenon for the undeformed airfoil, the shock wave for the morphed airfoil with a TE deflection $5^{°}$ initially appears near Mach 0.75, and became stronger on the upper surface at Mach 0.775 and subsequently moved backward towards TE at Mach 0.85 as shown in Figure 5(c). Two interesting things are noted. Firstly, when a strong shock wave forms on the upper surface of the airfoil, the flutter speed stops decreasing and starts to increase with increasing Mach number. Secondly, when the shock wave reaches the TE, another freeze region emerges at a lower Mach number around 0.85 as shown in Figure 7(b). Consequently, the unsteady aerodynamic force over the airfoil stabilizes and can lead to an increase of the flutter speed, which is consistent with the flutter boundary shown in Figure 6(a).

In summary, a possible explanation in improving the transonic flutter characteristic is that employing the morphing TE makes shock wave and Mach number freeze region occur at an early Mach number.

Appearance of shock wave and transonic flutter characteristic

From previous results, it can be seen that the shock wave location plays a vital role in the transonic flutter characteristics. A more general approach to alter the appearance of a shock wave is to vary the angle of attack (AoA) of the airfoil. To study the effect of morphing devices on the transonic flutter characteristics, the analysis was conducted at different AoA in this section. Special attention was paid to the relationship of the appearance of the strong shock wave on the upper surface of the airfoil and the transonic flutter boundary. It should be noted that Euler equations are applied in this section.

Similar to a morphing wing, the AoA effect on the wing transonic characteristic can be considered via the establishment of aerodynamic ROMs. The AoA of the Isogai wing model was set to a specific value in CFD solver when calculating aerodynamic coefficients due to ‘3211’ signals. Then the wing transonic flutter characteristics at this AoA was obtained by eigenvalue analysis.

Since Euler equations are adopted, the AoA is restricted to a small range from $0°$ to $4°$ in the current study. The resulting pressure coefficient distributions at different AoAs are presented in Figure 8, which shows the Mach number when the shock wave on the upper surface approaches the TE. At AoA =  $3°$ , for instance, a local supersonic flow appears on the upper surface at Mach number near 0.7, and a strong shock wave comes apparent at Mach 0.8, and moves backward to TE around Mach 0.84.

Figure 8.

Pressure distribution across transonic region at angle of attack of (a) $1°$ , (b) $2°$ , (c) $3°$ , (d) $4°$ .

Figure 9 shows the flutter boundary of the original Isogai wing model at different AoAs. It is noted that the bottom of the flutter boundary moves up as AoA increases. For small AoAs between $1°$ and $2°$ , the AoA effect on flutter speed at the transonic dip is unapparent. This might be caused by a strong shock wave occurring on the upper surface around Mach 0.8. When the shock wave moves close to TE, however, the Mach number remains around Mach 0.9 similar to the case of AoA =  $0°$ . When the AoA is further increased, remarkable change of transonic flutter boundary takes place. Taking AoA =  $3°$ for example, the downward trend of flutter speed starts to reverse at Mach 0.775, which is consistent with the appearance of a strong shock wave on the upper surface. When the shock wave moves close to the TE, the flutter speed increases significantly around Mach 0.84 showing a big difference from the AoA= $0°$ case.

Figure 9.

Flutter boundary of Isogai wing model at different angle of attack.

Yates et al. ³⁹ also noted that the growth of the supersonic area on the wing surface causes the aerodynamic center of the wing move backward, and the flutter speed increases with the AoA. Moreover, the numerical simulation by Edwards et al. ⁴⁰ indicates that the flutter speed rises when the shock wave forms on the upper surface. The results obtained in the current study are consistent with the above observations.

Another approach to evaluating the effect of the shock wave on the flutter boundary can be made by altering the AoA at the same Mach number. As the AoA increases at Mach 0.8, the flutter speed decreases slightly at small AoAs and then recovers, as shown in Figure 10(a). When the AoA is larger than 2°, a strong shock wave emerges on the upper surface as shown in Figure 10(b). Meanwhile, the downwards trend of flutter speed stops and starts to grow gradually.

Figure 10.

Effect of AoA on the flutter speed for Isogai wing model at Mach number 0.8: (a) The flutter speed at different AoA, (b) Steady pressure coefficient distribution at different AoA.

As shown in Figure 11, the scenario at Mach 0.84 is slightly different from Mach 0.8 when a strong shock wave has already been established at AoA= $0°$ , and the shock wave approaches to TE at AoA= $3°$ . The flutter speed increases with the AoA from $0°$ as shown in Figure 11(a). This is consistent with the trend of flutter speed over AoA at Mach 0.8 after a strong shock emerges. When the shock wave moves to TE at AoA= $3°$ , significant improvement of flutter speed is observed.

Figure 11.

Effect of AoA on the Isogai wing flutter speed at Mach number 0.84: (a) flutter speed at different AoA, (b) steady pressure coefficient distribution at different AoA, real part of aeroelastic eigenvalue versus non-dimensional air speed at (c) AoA =  $0.5°$ , (d) AoA =  $2°$ , (e) AoA =  $3°$ , (f) AoA =  $3.5°$ .

It is noted that the flutter boundary varied with the AoA and formed an S shape as shown in Figures 10(a) and 11(a). This was also observed from the conventional flutter boundary corresponding to Mach number in Isogai wing model section. When the AoA is increased at a fixed Mach number, the flow velocity on the upper surface of the airfoil is accelerated. A similar effect on the flow field exists by increasing the Mach number at a fixed AoA. When the AoA is between $0°$ and $1°$ in Mach 0.84, multiple flutter speeds take place as shown in Figure 11(c). After AoA is increased slightly above $1.5°$ , only a single flutter speed is detected due to the change of aerodynamic damping caused by the varying AoA. Thus, an isolated stable region is observed in the high speed area at low AoAs in Figure 11(a). When the AoA is further increased to greater than a specific AoA $3°$ , a jump of the flutter speed is exhibited. This is caused by the transition of ‘hump mode’ as shown in Figure 11(e) and (f).

In summary, the above results reveal that by varying Mach number, AoA or morphing TE deflection angle, the appearance of a strong shock wave on the upper surface of the airfoil leads to a reverse of the downward trend of flutter speed. When the shock wave on the upper surface reaches to the TE, the flutter speed starts to increase significantly. Thus, the mechanism of improving the transonic flutter characteristics by morphing TE is due to the shock wave emerging and moving to the TE at an early stage with lower Mach number. It is also note that the morphing LE has negligible effect on the shock wave and the flutter boundary in the transonic dip region.

Time domain transonic flutter analysis of a morphing wing

Due to the fact that morphing change requires a period of time to be implemented in structure, the aerodynamic configuration is being altered during this process, which may lead to unexpected influences on the aeroelastic behavior. For that reason, a time domain analysis considering TE morphing process is performed to assess these effects, which is conducted by the time marching method based on CFD.

In the following simulations, the morphing in the TE starts at $\tau =100$ with a constant morphing speed $\frac{d{\beta }_{TE}}{d\tau }$ varying from $0.05°$ to $0.25°$ , and the final deflection angle ( ${\beta }_{TE}$ ) is set to $5°$ . The resulting aeroelastic responses of the system are presented in Figure 12. It is observed that the morphing TE has a stable effect on the aeroelastic system, since the divergent motion turns to convergent when the morphing TE is employed. In addition to this, it was also found that the morphing speed has negligible influence on the stability of aerodynamic system, and only affects the transition period of the motion shifting from divergence to convergence.

Figure 12.

Aeroelastic response of morphing wing with different morphing speeds: (a) plunging motion and (b) pitching motion at U  =  4.5 of Mach 0.85, (c) plunging motion and (d) pitching motion at U  =  4.6 of Mach 0.875. ( ${\beta }_{TE}=5°$ starts at $\tau =100$ ).

Also observed in Figure 12, the static aeroelastic deflections of pitching and plunging motions converge to non-zero positions when the TE morphing is deployed. In particular, a negative mean AoA occurs because of the negative pitching aerodynamic moment induced by employing the morphing TE. From previous section, it is known that the AoA can influence the flutter characteristic of the aeroelastic system. Thus, the flutter speed obtained from the time marching method in this section is a result of both mean AoA and morphing TE. Note the results in this section are different from Figure 6 where only deformed shape is considered and zero mean AoA is assumed.

The flutter speed boundary obtained by time marching approach is displayed in Figure 13, and the flutter speeds for the original Isogai wing and morphing wing with constant mean AoA are presented. It is observed that the minimum flutter speed of flutter boundary is improved by 16.8% to 4.72 with respect to the original Isogai wing. During the current analysis, it was noted that the mean AoA induced by morphing TE is about ${-2}^{°}$ from Mach 0.8 to 0.85. So the flutter boundary of the morphed wing at mean AoA =  ${-2}^{°}$ is calculated, and the minimum flutter speed agrees well with that of time marching analysis, as shown in Figure 12. Comparing to the results in Figure 6, the minimum flutter speed improvement in this section is less, and the effective Mach region to increase flutter speed narrows, which is caused by the change of mean AoA induced by morphing TE.

Figure 13.

Flutter boundary of Isogai wing with the morphing trailing edge from time marching method.

Conclusion

The transonic flutter characteristic of a 2D airfoil system with morphing trailing and leading edges is studied. The Fluent software based on Euler and Navier-Stokes equations is employed with an RBF interpolation model to improve the capability of mesh deformation in the CFD model via UDF. An ARX model is implemented to establish aerodynamic ROM with the GA used to identify the model order. The critical flutter speed and frequency can be obtained by solving the eigenvalue equations at specific flow conditions. The flutter boundary can be subsequently obtained with a type of S shape flutter characteristic corresponding to AoA and Mach number for Euler solutions.

The relationship between shock wave locations on transonic flutter characteristics for the present model shows that the transonic flutter boundary depends on the forming and location of the shock wave on the upper surface of the airfoil. Corresponding to Mach number, the flutter speed will revise the downward trend and increase significantly following a strong shock wave taking place and moving close to the TE.

It is observed that the morphing TE provides an effective way to improve the transonic characteristics. For the presente wing model, the morphing TE leads to an increase of the transonic flutter speed by 36.9% based on Euler equation or 74.5% based on Navier-Stokes equation. The study reveals the a possible mechanism of the improvement of transonic flutter characteristics of morphing wing is that it makes shock wave and Mach number freeze region occur at a lower Mach number. However, the morphing LE has negligible effect on the transonic flutter characteristics because of its weak influence to the shock wave taking place near the TE.

Time domain analysis show that that the morphing speed has negligible influence on the stability of aerodynamic airfoil, and the minimum flutter speed of flutter boundary is improved by 16.8% with respect to the original wing. The effective Mach region for employing morphing TE to improve transonic flutter characteristic narrows, which is caused by the change of mean AoA.