Numerical Simulation of Folding Tail Aeroelasticity Based on the CFD/CSD Coupling Method

Abstract

This paper presents a CFD/CSD coupling method for aeroelastic simulation of folding tail morphing aircraft. The unsteady aerodynamic analysis is based on an in-house computational fluid dynamics (CFD) solver for the Euler equations, and emphasis is made on developing an efficient dynamic mesh method for the tail’s hybrid fold motion/elastic vibration deformation. The structural dynamic analysis is based on the computational structural dynamics (CSD) technique for solving the structural equation of motion in modal space. The aeroelastic coupling was achieved through successive iterations of CFD and CSD computations in the time domain. An adaptive multi-functional morphing aircraft allowing tail fold motion was selected to be studied. By using the developed method, aeroelastic simulation and mechanism analysis for fixed configurations at different folding angles and for variable configurations during the folding process were performed. The influence of folding rate on tail aeroelasticity and its influence mechanism were also analyzed.

1. Introduction

Compared with the traditional fixed configuration aircraft, the morphing aircraft can change its configuration in different mission stages so as to maximize the overall flight performance. In particular, the folding wing aircraft has attracted much attention in recent decades [1]. Taking the folding tail flying wing aircraft as an example, it has the abilities of high subsonic voyage, transonic maneuver, and supersonic cruise. When it needs to fly stealthily, the tail is folded to achieve good stealthy performance and high lift-to-drag ratio. When it needs to take off, land, or perform maneuvers, the tail is unfolded and can realize lateral-directional flight control through tail deflection. However, such aircraft is an active variable structure with lower stiffness and diverse materials, and the flow during the folding process is characterized by multi-scale unsteadiness and stronger nonlinearity. These may all lead to the presence of folding wing/tail aeroelastic problems.

In view of the fact that the aeroelastic tests of morphing aircrafts have the characteristics of complex testing technique, high cost, and long period, most of the current aeroelastic studies are carried out by theoretical and numerical means. All of them can be divided into two categories: one refers to those for configurations at fixed folding angles and the other one refers to those for time-varying configurations during the folding process. For the former one, the major focus is on the aeroelastic stability issues. For example, Liu et al. [2] presented a rapid aeroelastic analysis method for a folding wing aircraft by coupling the ROM (reduced order modeling) for nonlinear aerodynamics with the ROM for nonlinear structure. They used it to predict the onset of flutter at various folded configurations, and the results show that both the flutter dynamic pressure and the flutter frequency increase with the folding angle. Love et al. [3] presented a design optimization study for actuation systems of a folding wing aircraft in terms of folding wing angle, actuated system stiffness, and structural weight. The results can be applied in the design process and flight criteria specification for morphing aircraft structures. Snyder et al. [4] established a finite element model of a folding wing using MSC PATRAN/NASTRAN and analyzed the effects of the folding angle, actuator stiffness, and structural weight on the flutter characteristics. Dowell and Tang [5] modeled the folding wing structure based on the linear plate theory and divided it into three components: the fuselage, the inboard wing, and the outboard wing. By using a component modal analysis and a three-dimensional vortex lattice aerodynamic model, the structural dynamic behavior and aeroelastic stability of the folding wing are investigated.

When the wing is in fold motion, the structural dynamics change continuously, and the aeroelastic system is also varying with time, increasing the difficulty of aeroelastic analysis. Some researchers have been concerned with this issue, and the major focus is on the transient structural or aeroelastic analysis during the folding process. Zhao and Hu [6] developed a set of differential-algebraic equations that govern the time evolution of the folding wing based on the floating frame approach, and then used them to solve the transient responses of the wing during the morphing process. The equations can be integrated with the flow solver to carry out the transient aeroelastic analysis. Hu et al. [7] studied the aeroelastic characteristics of a folding wing during the folding process, where the structural model is built using the flexible multi-body dynamics approach and the aerodynamic model is built using the Kriging agent model technique. The results show that the folding and unfolding processes have opposite influences on the aeroelastic stability of the folding wing, and the influences become more significant at higher folding and unfolding rates. Reich and Bowman [8,9] developed a tool that combines a multi-body dynamic simulation package with an aerodynamic code for load computation and a Matlab-based flight control system, and they used it to simulate the flight of a folding wing vehicle. Xu et al. [10] constructed an aeroelastic simulation platform based on the secondary development of ADAMS software. This platform is applied to simulate the flight folding process, calculate the hinge moments of the folding wing, and investigate the influences of the folding rate and the aircraft’s center of gravity position on the results. Results show that for the fast folding process, the steady-state simulation errors of the hinge moments are substantially large, and a transient method is necessary. Dario et al. [11] described the morphing unmanned aerial vehicles (UAV) dynamics as a linear parameter-varying model and its simplification through the small convex hull approach. The ZAERO software is used to obtain aerodynamic coefficients at different folding angles and Mach numbers. A multi-loop controller is formulated for simulation of the entire folding process of the aeroelastic morphing UAV. Recently, Shen et al. [12] and Chen et al. [13] reported their studies on dynamics and analysis of tensegrity morphing airfoils, which could be a viable option for adaptive morphing wings.

Most of the above studies used a linear aerodynamic model for aeroelastic analysis, and thus it may bring larger errors for nonlinear flows such as transonic flows and separated flows. In recent years, the computational fluid dynamics/computational structural dynamics (CFD/CSD) coupling method has attracted much attention due to its advantages of fewer assumptions, higher simulation fidelity, and stronger versatility. By coupling the unsteady flow and structural dynamics solvers in the time domain, this method achieves a more realistic aeroelastic simulation. So far, it has gained a lot of applications in fixed configuration aircraft aeroelasticity. For example, Liu et al. [14] developed a parallel integrated CFD/CSD simulation program for flutter prediction of an aeroelastic system, and the code has been used to calculate transonic flutter of a two-dimensional airfoil and a three-dimensional wing. Yang et al. [15] summarized their researches in the recent years on the CFD/CSD-based aeroelastic methods and partly engineering applications, and they analyzed new flutter phenomena and mechanisms. Guruswamy et al. [16] has presented a perspective on CFD/CSD-based time accurate computational aeroelasticity methods since the late 1970s. The authors have also developed an in-house CFD/CSD code for aeroelastic simulation in many engineering areas [17,18,19]. By contrast, the relevant studies on morphing aircrafts are rarely seen. It is also observed that nearly all existing folding aircraft aeroelastic studies concentrate on the folding wing configurations, but few on the folding tail configurations, lacking the understanding of the characteristics and mechanisms of folding tail aeroelasticity. Note that some recent works on the aerodynamic characteristics of tail morphing have been reported [20,21].

This paper presents a CFD/CSD coupling method for folding tail aeroelasticity, and the major focus is on developing an efficient dynamic mesh method for tail’s hybrid fold motion/elastic vibration deformation and designing a flowchart of the calculation process. Based on the developed method, we performed an aeroelastic simulation and mechanism analysis of a folding tail morphing aircraft, and emphasis was on the flutter characteristics of the tail at different folding angles, the response characteristics during the folding process, and the influence of folding rate on tail aeroelasticity and its influence mechanism.

2. CFD/CSD Coupling Method

2.1. Aerodynamic Model

2.1.1. Governing Equations

For folding tail morphing aircrafts, the folding time period is large as compared with the typical time scale of flexible vibration, and a long-running computation may be required for simulation of the entire folding process. For the sake of computational efficiency, the Euler equations are used as the governing equations for unsteady flows, though the present method can also apply to the Navier–Stokes equations. Using the arbitrary Lagrangian Eulerian approach, the three-dimensional Euler equations can be written in conservative, differential form as

(1)$\begin{array}{l}\frac{\partial \rho }{\partial t}+\nabla \cdot [\rho (v-v_g)]=0\\ \frac{\partial \rho v}{\partial t}+\nabla \cdot [\rho v(v-v_g)+pI]=0\\ \frac{\partial \rho E}{\partial t}+\nabla \cdot [\rho E(v-v_g)+pv]=0\end{array},$

where $\rho$, $p$, and $E$ denote the density, the pressure, and the total energy per unit mass, respectively; $\mathit{v}$ is a vector of fluid velocities; ${\mathit{v}}_g$ is a vector of grid velocities; and $\mathit{I}$ represents the identity matrix.

The multi-block structured finite volume method is used for discretization of the governing equations, and we have

(2)$\frac{\partial }{\partial t}{\displaystyle {\int }_{\Omega }W}d\Omega +{\displaystyle {\oint }_S{F}_{c}dS}=0,$

where $\Omega$ is a control volume with a surface element $dS$, and $W$ and $F_c$ are the vectors of conservative variables and convective fluxes, respectively.

2.1.2. Numerical Technique

The Roe scheme is employed to evaluate the convective fluxes. To achieve a higher accuracy, the third-order monotone upstream centered scheme for conservation laws scheme is used to reconstruct the flow variables at cell interfaces. The reconstruction is based on the characteristic flow variables. The Van Albada’s limiter is used to prevent the generation of oscillation and spurious solutions in shock regions.

The dual-time stepping method is used for simulation of the unsteady flow, which transforms the unsteady problem into the steady-state problem in pseudo time. Through this method, the time stepsize limitation due to the stability condition can be relaxed and many acceleration techniques originally designed for steady computations can still be applied. By approximating the time derivative in Equation (1) with a three-point backward difference and treating other terms implicitly, we have

(3)$\frac{3{\Omega }^{(n+1)}{W}^{(n+1)}-4{\Omega }^{(n)}{W}^{(n)}+{\Omega }^{(n-1)}{W}^{(n-1)}}{2\Delta t}=-R^{(n+1)},$

where $\Delta t$ is the physical time stepsize; $R$ is the vector of residuals; and the superscripts $n-1$, $n$, and $n+1$ represent the last, current, and next time levels, respectively.

Introducing the pseudo time $\tau$, so that

(4)${\Omega }^{(n+1)}\frac{d{W}^{*}}{d\tau }=-\left(R(W^{*})+\frac{3W^{*}}{2\Delta t}-Q^{*}\right),$

where $W^{*}$ and $R(W^{*})$ are the approximations to $W^{(n+1)}$ and $R^{(n+1)}$, respectively, and $Q^{*}$ is the source term.

After obtaining the solution of the above equation, Equation (2) is naturally satisfied and thus $W^{(n+1)}=W^{*}$. In order to solve Equation (4), a first-order implicit scheme is adopted, and its linearized formulation is given by (omitting the “asterisk” that denotes the unsteady term)

(5)$\left[\Omega \left(\frac{1}{\Delta \tau }+\frac{3}{2\Delta t}\right)I+\frac{\partial R}{\partial W}\right]\Delta W^{(n)}=-R^{(n)},$

where $\Delta W^{(n)}=W^{(n+1)}-W^{(n)}$.

The implicit lower-upper symmetric Gauss–Seidel scheme is employed to solve the above equations, from which the conservative variables are solved.

Three types of boundary conditions are used in this work, namely, the farfield boundary, the slip wall boundary, and the symmetry boundary. In order to improve the computational efficiency, the Message Passing Interface (MPI) technique is used for parallel programming. Within the framework of multi-block structured grids, its main idea is to split the whole computational domain into several sub-domains, and the computational task on each sub-domain is carried out by a CPU core. The interprocess communication is realized by the MPI.

2.1.3. Dynamic Mesh Method

The folding tail morphing aircraft aeroelasticity belongs to the moving boundary problem, which involves rigid deformation due to tail folding and elastic deformation due to tail vibration. Such a hybrid deformation puts forward higher requirements for mesh deforming capabilities. In this work, we adopted a hybrid dynamic mesh method, namely, the RBF-TFI method [22], which combines the strong deformation ability of the RBF (radial basis function) method and high computational efficiency of the TFI (transfinite interpolation) method.

The RBF method is a desirable approach for scattered data interpolation, which introduces an interpolation function to describe the deformation in the whole space. For any deformation component, the RBF interpolation function can be expressed as a weighted sum of basis functions. Here, the thin plate spline basis function [23] is used, which gives

(6)$\mathsf{\Delta }\mathit{x}={\displaystyle \sum _{n=1}^N{a}_n{‖x-x_n^{*}‖}_2^2\log \left({‖x-x_n^{*}‖}_2+\delta \right)},$

where $\mathsf{\Delta }\mathit{x}$ denotes the deformation at any grid point $x$ and its value at $N$ control points are known, $x_n^{*}$ is a vector of the coordinates at the n-th control point, $a_n$ is the unresolved coefficient, and $\delta$ is a small quantity.

Writing Equation (6) for all control points, $a_n$ is solved from the combined equations by the Gaussian Elimination method with complete pivoting. Then, the deformations at all grid points can be calculated.

The TFI method is one of the most widely used algebraic techniques for mesh generation, which is computationally very efficient. In three-dimensional cases, the grid inside the computational domain can be generated using the Boolean sum of the interpolation functions, i.e.,

(7)$\mathsf{\Delta }\mathit{x}=U+V+W-UV-VW-UW+UVW,$

where $U$, $V$, and $W$ are the univariate interpolation functions in each of the coordinate directions in the computational space, and the other terms are the tensor products. More details can be found in [22].

For clarity, the major steps of the RBF-TFI method are as follows:

• (1). Select some characteristic points at the surface boundaries (e.g., farfield, wall, and symmetry) as the control points needed by the RBF interpolation, whose deformations are known.

• (2). Use the RBF method to calculate the deformations at all edges of grid blocks according to the deformations at control points.

• (3). Use the two-dimensional TFI method to calculate the deformations at all domains of grid blocks according to the deformations at all edges.

• (4). Use the three-dimensional TFI method to calculate the deformations at all interior points of grid blocks according to the deformations at all domains.

The practical application shows that for a million-scale grid of morphing aircraft, it takes seconds to generate dynamic mesh by the RBF-TFI method, and the quality of the deformed mesh is comparable to the original undeformed one.

2.2. Structural Model

Assuming a linear elastic model, the structural deformation vector $\Delta r$ at any location $(x,y,z)$ and any time $t$ is described using the modal superposition method, i.e.,

(8)$\Delta r\left(x,y,z,t\right)={\displaystyle \sum _{i=1}^n{q}_i\left(t\right)h_i\left(x,y,z\right)},$

where $n$ is the number of used vibration modes, $h_i$ is the i-th mode shape vector at any location $(x,y,z)$, and $q_i$ is the i-th generalized coordinate. Note that in the Cartesian coordinate system, the values of $h_i$ change during the tail folding process.

Applying Lagrange’s equations, the structural equation of motion in modal space can be written as

(9)$M{q}_{tt}+G{q}_t+Kq=A,$

where $M$, $G$, and $K$ are the generalized mass matrix, the generalized damping matrix, and the generalized stiffness matrix, respectively, $A$ is the generalized aerodynamic force with its i-th component

(10)$A_i={\displaystyle ∯p(x,y,z,t)\cdot h_i(x,y,z)dS},$

where $p$ denotes the unsteady surface pressure obtained from the CFD computation.

Because of the mode shape orthogonality, after the mass normalization, $M$ and $K$ become diagonal matrices, which consist the following components

(11)$K_{ii}={\omega }_i^2{M}_{ii}, {\omega }_i=2\pi f_i$

Introducing the state variable $E$ to reduce the order of the structural equation of motion, that is

(12)$E={\left[q_1,q_2,\cdot \cdot \cdot ,q_n,q_{t1},q_{t2},\cdot \cdot \cdot ,q_{tn}\right]}^T,$

then Equation (9) can be rewritten as a system of linear equations

(13)$E_t=\left[\begin{array}{cc}\mathbf{0}& I\\ -M^{-1}K& -M^{-1}G\end{array}\right]E+\left[\begin{array}{c}\mathbf{0}\\ M^{-1}\end{array}\right]A,$

where $\mathbf{0}$ denotes the zero matrix and $I$ denotes the identity matrix. The second-order predictor–corrector scheme [24] is used to solve the above equations so as to achieve an efficient coupling between fluid and structure.

2.3. CFD/CSD Coupling for Folding Tail Morphing Aircraft

Because of different physical characteristics in the fluid and solid domains, it is usually not desirable to generate a point-to-point mesh at the fluid–structure interface, and the deformations at structural mesh need to be interpolated onto the aerodynamic mesh. The thin plate spline method is used for interpolation. Its main idea is to introduce a load vector corresponding to a mode shape vector to solve the differential equation based on the virtual work principle. Traditional aeroelasticity studies require a repetitive deformation interpolation from structural to aerodynamic mesh and a repetitive force interpolation from aerodynamic to structural mesh, and the latter one relies on additional conservation conditions. This brings greater complexity and inconvenience to the CFD/CSD computation. Instead, here, we introduced the so-called spatial mode, that is, after the mode shape on aerodynamic surface mesh is obtained by interpolation, it is then extended to the whole spatial mesh using the aforementioned dynamic mesh method. By doing so, only a one-time mode shape interpolation from structural to aerodynamic mesh is required, and the efficiency of aeroelastic computation is improved.

The computational sequence of the present CFD/CSD coupling method for aeroelastic simulation of folding tail morphing aircraft is given below:

• (1). Given the range and rate of tail fold motion, generate the initial computational mesh for the unfolded configuration.

• (2). Calculate the rigid deformation according to the instantaneous folding angle at time t and transform the mode shape matrix.

• (3). Calculate the elastic deformation according to the instantaneous aerodynamic force at time t for the current folded configuration.

• (4). Combine the folding and elastic deformations and update the aerodynamic mesh using the RBF-TFI method.

• (5). Calculate the unsteady flow field and obtain the instantaneous aerodynamic force at time t + ∆t.

• (6). Repeat the procedures (2)–(5) until the entire tail folding process is finished.

For more clarity, Figure 1 shows the flowchart of the present aeroelastic simulation of folding tail morphing aircraft.

2.4. Code Validation

Before using the present method to study the folding tail aeroelasticity, we first validate our code through an AGARD 445.6 wing case (see Figure 2). It is a standard numerical example for examining the aeroelastic computation method [25]. The wing is composed of NACA65A004 airfoils, and it has an aspect ratio of 1.65, a taper ratio of 0.66, a span of 0.762 m, a root chord of 0.559 m, and a quarter-chord swept angle of ${45}^{\circ }$. An H-H topology multi-block structured mesh is used, which has a total of 0.62 million cells. The physical time stepsize is set as $5\times {10}^{-5}\mathrm{s}$. A CFL number of 100 is used, and the number of pseudo time steps is chosen to make the residual reduce by two orders or reach a maximum value of 50.

The weakened model reported in [25] is used. The first four natural modes are considered for flutter computation, which are the first bending, the first torsion, the second bending, and the second torsion modes. See Figure 3 for their respective mode shapes and natural frequencies. Note that the modal displacements shown in the figure have been normalized by their respective maximum values.

The CFD/CSD coupled method is applied to predict the flutter characteristics at several Mach numbers (0.499, 0.678, 0.901, 0.960, 1.072, 1.141). Specifically, at each Mach number, the free-stream dynamic pressure q is increased until a neutrally stable state is reached, in which case the flutter velocity and flutter frequency are obtained. For comparison with experiment, the non-dimensional flutter speed index $V^{\ast }$ and the flutter-frequency ratio ${\omega }^{\ast }$ are introduced

(14)$\begin{array}{cc}{V}^{\ast }=\frac{U_{\infty }}{b{\omega }_{\alpha }\sqrt{\mu }},& \hspace{1em}\end{array}{\omega }^{\ast }=\frac{\omega }{{\omega }_{\alpha }},$

where $U_{\infty }$ is the free-stream velocity, $b$ is the wing root semichord, ${\omega }_{\alpha }$ is the natural frequency of the first torsion mode, and $\mu$ is the mass ratio.

Figure 4 shows the computed time histories of generalized displacements at near neutrally stable states for typical Mach numbers. From these figures, we can not only determine the flutter point, but also find that this weakened model is prone to classical bending–torsion flutter. At $Ma=0.499$, flutter exhibited an explosive behavior that caused a sudden growth of the generalized response despite only a small change in dynamic pressure, while at $Ma=0.96$ and $Ma=1.072$, milder flutter behaviors were observed.

Figure 5 shows a comparison of flutter speed index and frequency ratio as a function of Mach number between the results from experiment and the present method. Also presented are those from some other inviscid computations, including Silva et al. [26] using the FUN3D code, Silva et al. [27] using the CFL3Dv6 code, and Lee and Batina [28] using the CFL3Dv1 code. A transonic dip around $Ma=0.96$ was observed in all the results. In the subsonic flow regime, the present results agreed well with the experimental data, while an overestimated flutter point was predicted in the high subsonic flow regime and an underestimated flutter point was predicted in the supersonic flow regime. Similar phenomena can also be found in other references. These results verify the effectiveness and accuracy of the developed CFD/CSD code for aeroelastic simulation.

3. Results and Discussions

3.1. Model Description

The model being studied in this work is an adaptive multi-functional morphing aircraft, which enables wing sweeping and tail folding. The main design parameters of aircraft include the length of 18.92 m, the wingspan of 14.5 m, and the tail span of 7.5 m. Here, the configurations with a fixed wing sweep angle of ${45}^{\circ }$ and a variable tail fold angle were considered. When the tail was unfolded, the aircraft had long-range cruise performance, and when the tail was folded, the aircraft was able to meet the requirements of lateral-directional flight control at subsonic and take-off/landing states. The schematic diagram for these two typical configurations is shown in Figure 6.

The fluid domain was set to make the distance to farfield boundary about 30 chord lengths and was meshed using an H-grid topology. The whole mesh consisted of 2.88 million cells, which was proven to be fine enough through a mesh-independence study. The physical time stepsize was set as $2\times {10}^{-5} \mathrm{s}$. A CFL number of 100 was used, and the number of pseudo time steps was chosen to make the residual reduce by two orders or reach a maximum value of 50.

For this morphing aircraft, the tail can fold between $0^{\circ }$ and ${60}^{\circ }$. As demonstrated above, the RBF-TFI method was used to automatically generate dynamic mesh at any folding angle $\delta$, and four typical meshes at different values of $\delta$ are presented in Figure 7 (tail surfaces highlighted with red).

Since the tail aeroelasticity is of our major concern, the fuselage and the wing were assumed to be rigid, and the tail was assumed to be elastic. The finite element model of the tail structure is displayed in Figure 8, which consists of beam, shell, and solid elements. A simple torsional spring stiffness was used to represent the potential flexibility in the actuation system.

Based on the structural modal analysis using PATRAN/NASTRAN 2012 software, the first six natural modes were considered. Their respective natural frequencies are given in Table 1. Two primary mode shapes are shown in Figure 9. Note that when the tail folds, the corresponding mode shape will change with the folding angle (i.e., rotate about the actuating shaft).

3.2. Aeroelastic Simulation and Mechanism Analysis

This section first studies the flutter characteristics at different fixed folding angles, then studies the aeroelastic response characteristics during the tail folding process, and finally analyzes the influence of the folding rate on the tail aeroelasticity and its influence mechanism.

A typical freestream flow condition near the flutter point was selected to be studied, which includes the Mach number $M{a}_{\infty }=0.8$, the angle of attack $\alpha =0^{\circ }$, and the dynamic pressure $q_{\infty }=40 \mathrm{kPa}$. Figure 10 shows the generalized displacements at this condition and at another condition with a lower dynamic pressure $q_{\infty }=38 \mathrm{kPa}$ for the unfolded configuration ($\delta =0^{\circ }$). It was observed that the tail may encounter the aeroelastic instability at $q_{\infty }=40 \mathrm{kPa}$. The potential flutter type is the bending–rotational flutter, which seems to be of a rather mild nature.

3.2.1. Flutter Characteristics at Fixed Folding Angles

We first computed the aeroelastic responses of the tail at fixed folding angles, which is equivalent to assuming a very slow folding process. Figure 11 shows the computed time histories of generalized displacements at several folding angles. It is seen that in all cases, the first bending and the rotational mode played a dominant role, while responses of other modes were much smaller. For a clear comparison, Figure 12 shows the responses for the two dominant modes, correspondingly. In cases of smaller folding angles ($\delta \le {15}^{\circ }$), the aeroelastic stability improved with the increase in $\delta$ and turned from an unstable state to a stable state at $\delta ={15}^{\circ }$, while it changed little in cases of higher folding angles. This phenomenon is thought to be related to the wing interaction. To be specific, at small $\delta$, the tail and the wing were almost on the same plane, and the tail was in a disturbed airflow environment, but at higher $\delta$, the tail was raised and it was in a “cleaner” airflow environment. At least for this case, the wing interference had negative effects on the aeroelastic stability of the tail. It was also observed that as the folding angle increased, the static deformation of the tail (relative to its unloaded position) decreased, and the response amplitude generally decreased. Consequently, it can be concluded that for a very slow folding process, the tail fold plays a role in reducing aerodynamic load and improving aeroelastic stability.

3.2.2. Response Characteristics during the Folding Process

Next, we computed the aeroelastic responses of the tail during the folding process following the flowchart shown in Figure 1. The folding range was between $0^{\circ }$ and ${60}^{\circ }$, and the folding rate was taken as a typical value of $5^{\circ }/\mathrm{s}$. Figure 13 shows the time histories of generalized displacements for the two dominant modes. Compared with the previous results for fixed configurations, the aeroelastic responses during the folding process exhibited nonlinear behaviors. Initially, at smaller $\delta$, the aeroelastic stability was rather poor (see above), and the response curve showed a divergence trend until near $t=2 \mathrm{s}$ ($\delta ={10}^{\circ }$). After that, the amplitude of response started to decay, and the decay rate at between $t=2.8 \mathrm{s}$ ($\delta ={14}^{\circ }$) and $t=4.4 \mathrm{s}$ ($\delta ={22}^{\circ }$) was slightly lower than at the other time. Finally, the generalized displacements of both two modes converged to a certain value. All the above results show the influence of the folding process on the tail aeroelastic characteristics.

3.2.3. Influence of Folding Rate and Influence Mechanism

To further investigate the tail aeroelasticity during the folding process, several folding rates were considered, i.e., $5^{\circ }/\mathrm{s}$, ${10}^{\circ }/\mathrm{s}$, ${20}^{\circ }/\mathrm{s}$, and ${40}^{\circ }/\mathrm{s}$. Figure 14 and Figure 15 show the variations of generalized displacements with time and with folding angle for different folding rates, respectively. It is seen that for smaller folding rates ($5^{\circ }/\mathrm{s}$ and ${10}^{\circ }/\mathrm{s}$), the aeroelastic response of the folding tail exhibited a type of envelope, that is, the response increased first and then decreased. This was because the smaller the folding rate, the longer the time history of the tail experiencing a range of smaller folding angles, and the more serious the divergence of response. When the folding angle is further increased, the response will rapidly converge.

On the other hand, for larger folding rates (${20}^{\circ }/\mathrm{s}$ and ${40}^{\circ }/\mathrm{s}$), the time variation of response was relatively small, and the aeroelastic characteristics during this fast folding process differed significantly from those for fixed configurations. This is because the growth of amplitude in the case of smaller folding angles had not been fully developed before the tail was out of flutter instability. When the folding angle was increased to a larger value, the aeroelastic stability of the tail continuously improved, and the responses for different folding rates converged rapidly to the same value. Therefore, with the increase in the folding rate, the folding process had a greater impact on the tail aeroelasticity, and overall, this dynamic process played a role in improving the aeroelastic stability.

4. Conclusions

In this paper, a numerical study on folding tail aeroelasticity of the morphing aircraft based on the CFD/CSD coupled method is presented. First, we introduce the numerical algorithms for unsteady aerodynamic analysis, structural dynamic analysis, and aeroelastic coupling in the time domain. In particular, an efficient RBF-TFI dynamic mesh method was used for describing the tail’s hybrid fold motion/elastic vibration deformation. Then, a benchmark flutter test case was used to validate the developed CFD/CSD code. Finally, aeroelastic simulation and analysis of a typical folding tail morphing aircraft were conducted, including the flutter characteristics at different fixed folding angles, the aeroelastic response characteristics during the tail folding process, and the influence of folding rate on tail aeroelasticity and its influence mechanism. The main findings of this work are summarized below:

• (1). For a very slow folding (quasi-static) process, the increase in tail folding angle reduced the aerodynamic load and improves the aeroelastic stability. The wing interference had negative effects on the aeroelastic stability, at least for this case.

• (2). The dynamic process of tail folding had a significant influence on the aeroelastic characteristics, and the responses exhibited nonlinear behaviors.

• (3). The influence of the tail folding process became greater as the folding rate increased, and generally it had a positive effect on the aeroelastic stability, at least for this case.

Author Contributions

Conceptualization, J.W. and T.G.; Methodology, D.Z. and W.L.; Validation, D.Z.; Investigation, D.Z. and W.L.; Resources, J.W. and B.L.; Data curation, B.L.; Writing—original draft, D.Z.; Writing—review & editing, H.X.; Visualization, H.G.; Supervision, T.G. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data is contained within the article.

Acknowledgments

This work is also partially supported by High Performance Computing Platform of Nanjing University of Aeronautics and Astronautics.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Flowchart of aeroelastic simulation of folding tail morphing aircraft based on the CFD/CSD coupling method.

Figure 2. Computational mesh of the AGARD 445.6 wing model.

Figure 3. Mode shapes of the AGARD 445.6 wing model. (a) Mode 1, first bending (f = 9.599 Hz). (b) Mode 2, first torsion (f = 38.165 Hz). (c) Mode 3, second bending (f = 48.348 Hz). (d) Mode 4, second torsion (f = 91.545 Hz).

Figure 4. Time histories of generalized displacements at near neutrally stable states for typical Mach numbers. (a) [Forumla omitted. See PDF.]. (b) [Forumla omitted. See PDF.]. (c) [Forumla omitted. See PDF.]. (d) [Forumla omitted. See PDF.]. (e) [Forumla omitted. See PDF.]. (f) [Forumla omitted. See PDF.].

Figure 5. Comparison of flutter speed index and frequency ratio as a function of Mach number between different results. (a) Flutter speed index. (b) Frequency ratio.

Figure 6. Schematic diagram of the morphing aircraft for different configurations. (a) Unfolded configuration. (b) Folded configuration.

Figure 7. Computational meshes at different tail folding angles. (a) [Forumla omitted. See PDF.]. (b) [Forumla omitted. See PDF.]. (c) [Forumla omitted. See PDF.]. (d) [Forumla omitted. See PDF.]. (e) [Forumla omitted. See PDF.].

Figure 7. Computational meshes at different tail folding angles. (a) [Forumla omitted. See PDF.]. (b) [Forumla omitted. See PDF.]. (c) [Forumla omitted. See PDF.]. (d) [Forumla omitted. See PDF.]. (e) [Forumla omitted. See PDF.].

Figure 8. Finite element model of the tail structure.

Figure 9. Two primary mode shapes of the tail structure. (a) Mode 1, first bending mode. (b) Mode 2, rotational mode.

Figure 10. Time histories of generalized displacements at [Forumla omitted. See PDF.] for the unfolded configuration ([Forumla omitted. See PDF.]). (a) [Forumla omitted. See PDF.]. (b) [Forumla omitted. See PDF.].

Figure 11. Time histories of generalized displacements at different fixed folding angles. (a) [Forumla omitted. See PDF.]. (b) [Forumla omitted. See PDF.]. (c) [Forumla omitted. See PDF.]. (d) [Forumla omitted. See PDF.]. (e) [Forumla omitted. See PDF.].

Figure 11. Time histories of generalized displacements at different fixed folding angles. (a) [Forumla omitted. See PDF.]. (b) [Forumla omitted. See PDF.]. (c) [Forumla omitted. See PDF.]. (d) [Forumla omitted. See PDF.]. (e) [Forumla omitted. See PDF.].

Figure 12. Time histories of generalized displacements for the two dominant modes. (a) First bending mode. (b) Rotational mode.

Figure 13. Time histories of generalized displacements during the tail folding process. (a) First bending mode. (b) Rotational mode.

Figure 14. Variations of generalized displacements with time for different folding rates during the tail folding process. (a) First bending mode. (b) Rotational mode.

Figure 15. Variations of generalized displacements with folding angle for different folding rates during the tail folding process. (a) First bending mode. (b) Rotational mode.

Table 1

Natural frequencies of the tail structure.

Mode Order	1	2	3	4	5	6
Natural frequency (Hz)	27.595	48.453	196.89	273.51	327.72	483.72

References

1. Rodriguez, A. Morphing Aircraft Technology Survey. Proceedings of the 45th AIAA Aerospace Sciences Meeting and Exhibit; Reno, NV, USA, 8–11 January 2007.

2. Liu, D.; Chen, P.-C.; Zhang, Z.; Wang, Z.; Yang, S.; Lee, D.-H.; Mignolet, M.; Liu, F.; Lindsley, N.; Beran, P. Continuous Dynamic Simulation for Morphing Wing Aeroelasticity. Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference; Palm Springs, CA, USA, 4–7 May 2009.

3. Love, M.; Zink, S.; Stroud, R.; Bye, D.; Chase, C. Impact of Actuation Concepts on Morphing Aircraft Structures. Proceedings of the 45th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics & Materials Conference; Palm Springs, CA, USA, 19–22 April 2004.

4. Snyder, M.P.; Sanders, B.; Eastep, F.E.; Frank, G.J. Vibration and Flutter Characteristics of a Folding Wing. J. Aircr.; 2009; 46, pp. 791-799. [DOI: https://dx.doi.org/10.2514/1.34685]

5. Tang, D.; Dowell, E.H. Theoretical and Experimental Aeroelastic Study for Folding Wing Structures. J. Aircr.; 2008; 45, pp. 1136-1147. [DOI: https://dx.doi.org/10.2514/1.32754]

6. Zhao, Y.; Hu, H. Prediction of transient responses of a folding wing during the morphing process. Aerosp. Sci. Technol.; 2013; 24, pp. 89-94. [DOI: https://dx.doi.org/10.1016/j.ast.2011.09.001]

7. Hu, W.; Yang, Z.; Gu, Y. Aeroelastic study for folding wing during the morphing process. J. Sound Vib.; 2016; 365, pp. 216-229. [DOI: https://dx.doi.org/10.1016/j.jsv.2015.11.043]

8. Reich, G.; Bowman, J.; Sanders, B.; Frank, G. Development of an Integrated Aeroelastic Multi-Body Morphing Simulation Tool. Proceedings of the 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference; Newport, RI, USA, 1–4 May 2006.

9. Bowman, J.; Reich, G.; Sanders, B.; Frank, G. Simulation Tool for Analyzing Complex Shape-Changing Mechanisms in Aircraft. Proceedings of the AIAA Modeling and Simulation Technologies Conference and Exhibit; Keystone, CO, USA, 21–24 August 2006.

10. Xu, H.; Han, J.; Yun, H.; Chen, X. Calculation of the Hinge Moments of a Folding Wing Aircraft during the Flight-Folding Process. Int. J. Aerosp. Eng.; 2019; 9362629. [DOI: https://dx.doi.org/10.1155/2019/9362629]

11. Baldelli, D.H.; Lee, D.-H.; Pena, R.S.S.; Cannon, B. Modeling and Control of an Aeroelastic Morphing Vehicle. J. Guid. Control Dyn.; 2008; 31, pp. 1687-1699. [DOI: https://dx.doi.org/10.2514/1.35445]

12. Shen, Y.; Chen, M.; Skelton, R.E. Markov data-based reference tracking control to tensegrity morphing airfoils. Eng. Struct.; 2023; 291, 116430. [DOI: https://dx.doi.org/10.1016/j.engstruct.2023.116430]

13. Chen, M.; Shen, Y.; Skelton, R.E. Model-Based and Markov Data-Based Linearized Tensegrity Dynamics and Analysis of Morphing Airfoils. Proceedings of the AIAA SCITECH 2024 Forum 2024; Orlando, FL, USA, 8–12 January 2024.

14. Liu, F.; Cai, J.; Zhu, Y.; Tsai, H.M.; Wong, A.S.F. Calculation of Wing Flutter by a Coupled Fluid-Structure Method. J. Aircr.; 2001; 38, pp. 334-342. [DOI: https://dx.doi.org/10.2514/2.2766]

15. Yang, G.; Zheng, G.; Li, G. Computational methods and engineering applications of static/dynamic aeroelasticity based on CFD/CSD coupling solution. Sci. China Technol. Sci.; 2012; 55, pp. 2453-2461. [DOI: https://dx.doi.org/10.1007/s11431-012-4935-1]

16. Guruswamy, G.P. Perspective on CFD/CSD-Based Computational Aeroelasticity during 1977–2020. J. Aerosp. Eng.; 2021; 34, 06021005. [DOI: https://dx.doi.org/10.1061/(ASCE)AS.1943-5525.0001329]

17. Guo, T.; Lu, Z.; Tang, D.; Wang, T.; Dong, L. A CFD/CSD model for aeroelastic calculations of large-scale wind turbines. Sci. China Technol. Sci.; 2013; 56, pp. 205-211. [DOI: https://dx.doi.org/10.1007/s11431-012-5028-x]

18. Lv, B.; Lu, Z.; Guo, T.; Tang, D.; Yu, L.; Guo, H. Investigation of winglet on the transonic flutter characteristics for a wind tunnel test model CHNT-1. Aerosp. Sci. Technol.; 2019; 86, pp. 430-437. [DOI: https://dx.doi.org/10.1016/j.ast.2019.01.025]

19. Zhou, D.; Lu, Z.; Guo, T.; Chen, G. Aeroelastic prediction and analysis for a transonic fan rotor with the “hot” blade shape. Chin. J. Aeronaut.; 2021; 34, pp. 50-61. [DOI: https://dx.doi.org/10.1016/j.cja.2020.10.018]

20. Ma, X.; Su, J.; Zhong, S.; Huang, Y.; Zhang, Y. Study of aerodynamic and stealthy performance for a multifunctional morphing tail. Acta Aerodyn. Sin.; 2020; 38, pp. 896-900.

21. Jiang, Q.-F.; Chen, S.-S.; Yang, H.; Zhang, Z.-G.; Gao, Z.-H. Unsteady aerodynamic characteristics of a morphing tail configuration. Phys. Fluids; 2024; 36, 036133. [DOI: https://dx.doi.org/10.1063/5.0199736]

22. Ding, L.; Guo, T.; Lu, Z. An Efficient Dynamic Mesh Generation Method for Complex Multi-Block Structured Grid. Adv. Appl. Math. Mech.; 2015; 6, pp. 120-134. [DOI: https://dx.doi.org/10.4208/aamm.2013.m199]

23. Harder, R.L.; Desmarais, R.N. Interpolation using surface splines. J. Aircr.; 1972; 9, pp. 189-191. [DOI: https://dx.doi.org/10.2514/3.44330]

24. Zhang, W.; Jiang, Y.; Ye, Z. Two better loosely coupled solution algorithms of CFD based aeroelastic simulation. Eng. Appl. Comput. Fluid Mech.; 2007; 1, pp. 253-262. [DOI: https://dx.doi.org/10.1080/19942060.2007.11015197]

25. Yates, E.C., Jr. AGARD Standard Aeroelastic Configurations for Dynamic Response I-Wing 445.6, Proceedings of the 61st Meeting of the Structures and Materials Panel, Oberammergau, Germany, 8–13 September 1985; Advisory Group for Aerospace Research and Development: Neuilly sur Seine, France, 1988.

26. Silva, W.A.; Chwalowski, P.; Perry, B., III. Evaluation of linear, inviscid, viscous, and reduced-order modelling aeroelastic solutions of the AGARD 445.6 wing using root locus analysis. Int. J. Comput. Fluid Dyn.; 2014; 28, pp. 122-139. [DOI: https://dx.doi.org/10.1080/10618562.2014.922179]

27. Silva, W.A.; Bartels, R.E. Development of reduced-order models for aeroelastic analysis and flutter prediction using the CFL3Dv6.0 code. J. Fluids Struct.; 2004; 19, pp. 729-745. [DOI: https://dx.doi.org/10.1016/j.jfluidstructs.2004.03.004]

28. Lee-Rausch, E.M.; Batina, J.T. Wing flutter boundary prediction using unsteady Euler aerodynamic method. J. Aircr.; 1995; 32, pp. 416-422. [DOI: https://dx.doi.org/10.2514/3.46732]