Abstract
Several new systematic methods for high fidelity and reliability calculation of static and single dynamic derivatives are proposed in this paper. Angle of attack step response is used to obtain static derivative directly; then translation acceleration dynamic derivative and rotary dynamic derivative can be calculated by employing the step response motion of rate of the angle of attack and unsteady motion of pitching angular velocity step response, respectively. Longitudinal stability derivative calculations of SACCON UCAV are taken as test cases for validation. Numerical results of all cases achieve good agreement with reference values or experiments data from wind tunnel, which indicate that the proposed methods can be considered as new tools in the process of design and production of advanced aircrafts for their high efficiency and precision.
1. Introduction
How to determine an aircraft’s stability characteristics at the boundary of flight envelope is one of the most complicated and important aspects in the process of advanced aircraft design, as well as the aerodynamic and flight dynamic analysis. Normally, the stability of the aircraft is categorized as two parts, that is, static stability and dynamic stability [1–3]. The former one relates to the initial response of the aircraft with disturbance, which represents the trend that the aircraft returns to its original balanced state under restoring moments. In contrast, the dynamic stability indicates the property that the vehicle recovers to its original position after disturbance with dynamic damping moments in terms of time response process.
Both the static and dynamic stabilities play vital roles in flight control system design and aerodynamic optimization. In 1911, Bryan [4] treated stability as aerodynamic derivative assuming that the aerodynamic forces and moments are functions of the instantaneous values of the disturbance velocities, control angles, and their rates. For slow motion at low angle of attack, velocities and control angles that can model the aerodynamic loads independently are often called the static derivatives, while their rates used to calculate the stability at higher angle of attack and rates are called dynamic derivatives. Traditionally, several methods have been studied to estimate derivatives in aircraft design. Among them, wind tunnel testing has a well-known physical realism and is able to acquire detailed static and dynamic derivatives. However, the corresponding applications are limited by blockage, scaling, and Reynolds number effects together with support interference issues that prevent the proper modeling of the full-scale vehicle behavior [5]. Flight testing method may be the most physically accurate way but it usually cannot be implemented in the early stage of aircraft design due to the prohibitively huge cost. With the development of computational fluid dynamics (CFD) solvers, it is possible to acquire stability derivatives with high fidelity and robust parallel CFD methods, which have become a promising way for advanced aircraft design.
Recently many research groups have studied a variety of techniques to compute stability and control derivatives based on CFD. Ronch et al. [6, 7] presented a time-domain solver to compute dynamic derivatives of full aircraft configurations. They investigated the longitudinal stability derivatives of the Standard Dynamic Model (SDM) and the Transonic Cruiser (TCR) by using the Fourier integral and linear regression mathematical model. The results of both methods agreed well with experimental data. However, these methods could not present more detailed information about the stability derivatives. Combined with Cartesian meshing, Murman [8] developed a nonlinear, reduced frequency approach to simulate the response with a forced oscillation to obtain dynamic stability derivatives. The method was rather efficient but not so accurate. Park and Green [9] used a steady way to quickly evaluate the derivatives with the modification of CFL3D code. It was with high efficiency but poor universality since not all the derivatives could be obtained. Based on arbitrary Lagrangian-Eulerian (ALE) formulation with a dynamically deforming mesh algorithm, Oktay and Akay [10] proposed an unsteady method to calculate the dynamic derivatives of the Basic Finner Missile (BFM) model. Predicted results showed good agreement with wind tunnel data. However, the mesh deformation would highly increase the cost time. Park et al. [11, 12] presented a dual-time stepping algorithm combined with a multigrid DADI method, which was also a traditional method to simulate the forced oscillation to identify stability derivatives, while it is accurate but inefficient. Green et al. [13] employed a low-order, time dependent panel method CFD code to analyze the stability derivatives of an F-16XL model, but the method was highly questionable for other configurations. Mader and Martins [14] modified dual-time method to calculate dynamic derivatives and their sensitivities to further optimize the aerodynamic shape of aircraft. The stability-constrained optimization demonstrated the viability of the method. However, it was still the initial step to optimize the dynamic characteristics of the aircraft only considering the simplest derivatives. Limache and Cliff [15] examined the analytic sensitivity method to compute dynamic derivatives, and they showed that the multiple derivatives could be obtained from a single solution, but it is not an easy job to extend this method to general application.
Up to now, although significant progress has been made for the stability derivative identification with CFD, it is inconvenient to use most of those existing numerical methods for the further analysis on flight dynamics. The reason is that only the poor accurate static derivatives and combined forms of dynamic derivatives can be calculated. Moreover, it is difficult to apply these complicated methods to general analysis to satisfy the requirement of advanced aircraft design.
In this article, several new systematic methods to calculate static and single dynamic derivatives are developed based on CFD technique. With aircraft forced to motion in step response form, the static derivative and each single dynamic derivative can be quickly simulated from the identification of unsteady aerodynamic forces and moments. Longitudinal numerical results of a fly wing configuration using the new methods are compared to reference data and experimental measurements. Different from conventional ways, not only the new methods are recommended for their high precision and computational efficiency, but also the high fidelity and reliability single dynamic derivatives could be obtained. In addition, they can be generalized to the lateral and directional stability derivatives simulation.
2. Numerical Method
We chose Ansys Fluent [16] for the calculation of the steady and unsteady flow fields in this article, which provides comprehensive modeling capabilities for a wide range of incompressible and compressible, laminar, and turbulent flow problems. Moreover, it has been widely used and tested in the field of aerospace.
Three-dimensional Navier-Stokes equations can be written aswhere indicates the control volume, denotes the vector of conserved variables ( is the total energy per unit mass), represents the density, and , , and are the components of velocity given by the Cartesian velocity vector . In addition, denotes the inviscid component of the flux vectors, while is the viscous flux vector, which contains terms of the heat flux and viscous forces exerted on the body. Re is the Reynolds number.
In this study, shear-stress transport (SST) model is used which is developed by Menter to effectively combine model in the near-wall region for its robust and accurate prediction and model in the far field with free-stream independence.
The rigid dynamic mesh technique is adopted to model flows when the whole domain is flexible against time in order to identify the stability derivatives for several unsteady motions. This model is a special case of general dynamic mesh motion wherein the nodes move rigidly in the whole zone. It can be conveniently used in cases of single freedom motions with only one body.
3. Dynamic Derivative Calculation Method
3.1. Traditional Methods
The static derivative is the slope of pitching moment coefficient, and its conventional estimation is obtained by interpolating the steady aerodynamic forces and moments, as shown in Figure 1. This method is simple but not accurate, especially when used to analyze the dynamic characteristics of aircraft. Furthermore, it cannot capture the interference of motion to static derivative.
The classical dynamic derivative simulation method is based on wind tunnel testing technique. By imposing a forced sinusoidal motion around the aircraft center of gravity, longitudinal dynamic derivative values from the time-history of the forces and moments are calculated based on the assumption that the aerodynamic coefficients are linear functions of the angle of attack , pitching angular velocity , and their rates and . The mathematical expression of the forces and moments in terms of these variables during the applied pitching can be obtained by a Taylor series expansion as
The higher order terms, and , are often discarded; then the unsteady moment coefficient can be described as
The harmonic motion in pitch defines the kinematic relations for the angle of attack, pitching angular velocity, and rate as
Substituting (4) into (3), the new equation is
The unsteady moment coefficient periodically varies with stable amplitude over several cycles, and the primary effects of unsteady characteristic can be omitted.
When and , (5) is, respectively, expressed aswhere indicates the nondimensional reduced frequency of the applied motion. Representing the pitching moment coefficient damping, in (6) is usually called out-of-phase derivative or combined dynamic derivative due to the fact that and are difficult to separate. While in (7) named as in-phase derivative is comprised of static derivative and rotational derivative .
As discussed above, time cost is rather big when this method is used to calculate the unsteady several cycles’ unsteady motion. Although balance method or time spectral method [19–21] have been used to alleviate computation cost by converting the unsteady calculation from time domain to frequency domain, only coupling terms of dynamic derivative can be obtained in these methods.
3.2. New Systematic Methods Based on Step Response to Calculate Stability Derivatives
Aiming at overcoming the disadvantage of conventional approaches, new systematic methods in this paper for calculating reliable single static and dynamic derivatives are mainly based on step response motion. Actually the concepts of step response of a system in a given initial state consist of time evolutions of its outputs when its control inputs are Heaviside step functions [22]. In electronic engineering and control theory, step response is the time behavior of the outputs when its inputs change from zero to one in a very short time, which can be simply represented in Figure 2. Drawing lesson from this concept and combined with CFD, we calculate the unsteady aerodynamic forces and moments by imposing step response motions to the aircraft and then directly obtain the static and dynamic derivatives through simple identifications. Here we take longitudinal cases as example to introduce the new methods in detail.
(1) Static Derivative Calculation. Static derivative is the moment coefficient about the store body pitch axis due to angle of attack. It indicates the effect of aerodynamic moment increment caused by unit change of angle of attack. The aircraft can be classified as statically unstable if , or stable if else. In order to simulate this derivative, a step response motion with constant additional angle of attack is applied to the aircraft; see Figures 3 and 4. Note that the additional angle of attack should be transformed to additional velocity perpendicular to free-stream velocity with following relation:
Then
Hence, this unsteady motion is only affected by , and the unsteady moment coefficient can be expanded as
It should be noticed that (10) has neglected the higher order terms of , which is validated to be accurate enough to simulate linear or weakly nonlinear cases at small angles of attack. The static derivative is readily solved with the unsteady and steady aerodynamic coefficients by
In summary, this method involves unsteady effect and reflects the physical realism compared to the conventional one.
(2) Acceleration Derivative Calculation. The translation acceleration derivative represents the finite time that aerodynamic loads at the tail lag changes in downwash converted downstream from the wing. Hence, it is also called the derivative of lag of wash. This derivative indicates the aerodynamic moment increment caused by the unit change of angle of attack rate . Since the form of is similar to the pitching angular velocity when the free-stream velocity is unchanged, the single derivative is difficult to be calculated using the conventional methods. Now we can use the step response motion to solve this problem.
Shown in Figure 5, we use a method to force the aircraft to move with constant angle of attack rate , which is called the step response motion of angle of attack rate. Furthermore, the additional angle of attack should be transformed to additional velocity perpendicular to free-stream velocity; see Figure 6.
When the constant rate of angle of attack is given, the instantaneous additional angle of attack at time can be defined as
Combine (12) with the relation
Then the additional velocity is written as
Since aircraft is approximately in the uniformly accelerative motion for small disturbance, the unsteady moment coefficient related to only and is expressed as
Equation (15) has discarded higher order terms of and based on small-disturbance theory, and this is a reasonable assumption which will be proved later. Similarly, the final additional angle of attack can be defined as the same form as that used for calculating , and now the transition acceleration derivative value is able to obtained by subtracting the instantaneous unsteady coefficient calculating from that of calculating at time , aswhere indicates the nondimensional angle of attack rate.
(3) Rotary Derivative Calculation. Rotary derivative is the angular velocity derivative, which indicates the change of aerodynamic coefficient due to the increment of pitching angular velocity. When the aircraft rotates around its center of gravity, the pitching moment, which is produced by additional aerodynamic forces of the body behind or in front of the center of gravity, always seems to impede the movement. That is reason why it is named as damping moment. Correspondingly, the derivative can also be called damping derivative.
Still using step response method, we can identify this single dynamic derivative by enforcing the aircraft pitching around the center of gravity with constant angular velocity , as shown in Figures 7 and 8, and the additional angle of attack is given by
As the velocity of free-stream is constant, we have
Then the unsteady aerodynamic coefficient related to all the variables , , and is described as
Higher order terms of and are legitimately discarded in (19) as previous. The pitching angular velocity and simulation time are defined as same values as that in the calculation procedure of ; then the final instantaneous additional angle of attack will be also identical in these two cases. As a result, the single derivative is obtained by subtracting the instantaneous unsteady coefficients of the two cases at time aswhere implies the nondimensional pitching angular velocity.
Now the new systematic methods are summarized as the flow diagram in Figure 9. Starting with steady case, all the single static and dynamic derivatives can be calculated following these steps.
4. CFD Modeling
4.1. Computational Model and Grid
A fly wing wind tunnel testing model, namely, a generic UCAV geometry called SACCON “Stability and Control CONfiguration,” is considered as test case to verify the methods. The model is used by an Applied Vehicle Technology Task Group (AVT-161) established by the NATO Research and Technology Organization (RTO) to assess the state-of-the-art in CFD methods for the prediction of static and dynamic stability and control characteristics of military vehicles in the air and sea domains. The wind tunnel tests of both static and dynamic aerodynamics were carried out in the Low Speed Wind Tunnel Braunschweig (DNW-NWB) [18].
The SACCON UCAV model has a lambda wing platform with a leading edge sweep angle of 53°, seen in Figure 10. The root chord is approximately 1 m and the wing span is 1.53 m, while the mean aerodynamic chord is 0.479 m and the reference wing area is 0.77 m2. The moment reference point (MRP) locates 0.6 m from the head of the wing, and the point of rotation (POR) for dynamic testing is 0.255414 m behind the MRP. More detailed description of the geometry can be seen in [18].
The computational mesh is generated by ICEM CFD with approximately 7 million cells (as shown in Figures 11 and 12). The extent of the boundary layer is 4.9 mm and the thickness of the first layer is 0.0036 mm. The largest value of does not exceed the value of 1.2. The mesh cells are clustered around the leading edges and the volume mesh has enhanced density around the aircraft to insure adequate resolution of the flow filed.
4.2. Static Aerodynamic Results
The static computations without a sting are taken under the same conditions with wind tunnel tests. The Mach number is 0.147 while the Reynolds number based on mean aerodynamic chord is . Comparisons of lift, drag, and pitching moment coefficient are, respectively, shown in Figures 13–15. The lift and drag coefficients have good agreement with the wind tunnel data or the reference data [23] throughout the angle of attack below 17 degrees, while both the lift coefficients calculation and reference data diverge when ranging up to larger angle of attack. Therefore, none of the simulations well match the test data.
The changing process of pitching moment coefficients are more complicated with the increment of the angle of attack. First of all, both the calculation and reference data have bad agreement with the experiment results when the angle of attack is below nearly 15 degrees, including the magnitude and slope of profiles. The possible reason may be that the sting and wind tunnel wall effects are not considered in the simulation. The predictions become better when the angle of attack increases further. All the methods capture the sharp reduction in nose-up pitching moment, with a sharp increase afterwards in positive moment between 17 and 20 degrees. This indicates that all the methods can well describe the basic vortex variations during the static process. Detailed analyses on the nonlinear variations of the pitching moment can be seen in [23–25] and we present the surface pressures and streamlines of several angles of attack in Figures 16(a)–16(d).
(a) 0 deg incidence
(b) 5 deg incidence
(c) 10 deg incidence
(d) 15 deg incidence
4.3. Static Derivative Simulation
According to step response method in the paper, an additional angle of attack is adopted to impose the aircraft to move unsteadily, and we set the additional velocity as . The unsteady aerodynamic coefficients are shown in Figure 17. Due to the increment of constant angle of attack, the unsteady aerodynamic loads arising from forced motions converge to a steady value after the decay of the initial transient.
Now we can calculate the static derivatives according to the new method shown in (11). With parameters shown in Table 1, it is only necessary to derive the additional angle of attack from the difference between the unsteady aerodynamic force of step motion and the steady force at initial angle of attack. When the static derivatives are obtained, the traditional method is also used to identify the static derivatives as comparison. With the difference of pitching moment coefficient corresponding to the two angles of attack before and after the initial angle of attack, the static derivative can be further calculated by deriving the difference of angles of attack. The steady data used in this method are shown in Figure 18, and the results can be seen in Table 2.
The static derivatives with new method show good agreement with reference data in Table 3, while conventional values have large deviation in comparison. It indicates that the proposed new method is more effective and precise. Moreover, if the interference factors such as sting and wall effects are considered in detail, this method may achieve better results.
4.4. Acceleration Derivative Simulation
Different from calculating static derivative, the translation acceleration derivative is obtained from step response motion of the angle of attack rate, where the additional angle of attack will be linear change with time. In this paper, the constant angle of attack rate is given as , and the instantaneous angle of attack is . At time , keeps the same with that in the calculation of . The variations of the unsteady aerodynamic moment coefficients are shown in Figure 19. We can see that as the time steps increase, the unsteady aerodynamic loads also exhibit an approximately linear relationship at these initial angles of attack.
Based on the unsteady aerodynamic forces calculated in Section 4.3 and this section, the single acceleration derivative can be obtained with the new method shown in (16). Since the difference between the two unsteady cases is only the effect of the rate of angle of attack at time , then the acceleration derivative is calculated by deriving with the difference of the two unsteady aerodynamic forces. Table 4 describes the specific calculation process.
Numerical results of in Table 4 are in good agreement with experimental measurements obtained at ten and fifteen degrees. A corresponding increment in the translation acceleration derivative is predicted when initial angle of attack increases, which can be related to the vortices creating at the wing leading edge and moving along the upper wing surface. Furthermore, the local pressure and aerodynamic forces will be altered during the dynamic process, as well as the longitudinal dynamic stabilities.
4.5. Rotary Derivative Simulation
The pitching angular velocity step response motion with is applied to identify the rotary derivative , and the instantaneous additional angle of attack is . When the time is , the angle of attack is 0.01 rad and its rate is . The other values are the same with those in the calculation of . The unsteady aerodynamic moment coefficients are shown in Figure 20. We see that, at these angles of attack, all the coefficients vary linearly as time increases, but this relationship may be broken in larger angle of attack.
When the calculating time is , the addition angles of attack are 0.01 rad for both of the two cases in Section 4.4 and this section. However, the difference is that the influence factor of the translational step motion in Section 4.4 is the rate of angle of attack , while the factor of rotational step motion in this section is . Thus, we can further calculate the single rotary derivative with (20) by deriving with the difference of aerodynamic forces of the two unsteady cases. The detailed process of calculation can be seen in Table 5.
Table 5 compares the simulation results of rotary derivative with wind tunnel data. Reasonable agreements are obtained between the two methods at ten and fifteen degrees with a maximum error 4.2%. The values of derivative indicate that the damping moment caused by pitching angular velocity will be reduced at large angle of attack.
4.6. Efficiency and Application of the Methods
The major computational cost for traditional methods is the computation of time-accurate simulations in response to periodic motions, which is avoided in the new systematic methods with simple step response motions. If we take the case of SACCON UCAV as an example, a typical calculation requires about 20 hours of CPU time on 24 processors to identify the combined dynamic derivatives, while only about 12 hours are needed to calculate the static and two single dynamic derivatives with the new methods.
For more details, the time-consuming comparison of static and dynamic derivatives calculation with traditional and new methods at a given angle of attack is shown in Figure 21. Note that both the steady and unsteady calculations adopt the second-order accuracy method. The traditional method can only achieve the combined dynamic derivative by using the sinusoidal motion, and at least three cycles should be simulated to ensure the periodic changes of unsteady aerodynamic forces. While the computational time is one second for both the three unsteady cases in new method with the same conditions as traditional methods.
If we check the total time to identify the static and dynamic derivatives, it takes about 2 and 20 hours for the calculation of static derivative and combined dynamic derivatives, respectively. In contrast, the new method can reduce the total time to about 12 hours with 4 hours decreasing in each unsteady case. It can be seen that the new method in this paper not only obtains more precise results of static and single dynamic derivatives but also greatly enhances the computational efficiency.
Furthermore, if the angle of attack in the paper is replaced with angle of sideslip, the methods can be extended to directional stability derivative simulation or further modified and used for the calculation of lateral stability derivatives if combined with the relation among angle of attack, sideslip, and rolling angle.
5. Conclusions
New systematic CFD methods based on step response motions are proposed in this study, which effectively and accurately estimate static and single dynamic stability derivatives. Beginning with the calculation of steady cases, the procedures of the methods employ angle of attack, rate of angle of attack, and pitching angular velocity step response to identify the static derivative and two single dynamic derivatives: translation acceleration and rotary derivatives.
Compared with conventional ways, those new methods can provide more detailed and precise static and single dynamic derivatives with lower computation power cost. Besides being applied to directional stability simulation by replacing the angle of attack with angle of sideslip, a modification with the relation among angle of attack, angle of sideslip, and rolling angle will extend the methods to lateral stability derivatives computation.
It is very convenient to use these methods for the production of advanced aircraft. However, the presented results are initial steps in the process to develop new ways for the calculation of static and single dynamic stability derivatives, which is still based on small-disturbance theory with linear or weakly nonlinear flow at medium-small angle of attacks. Future work will focus on the deep evaluation and modification of these methods at high angle of attack where more complex and nonlinear aerodynamic phenomena take place.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
The authors would like to acknowledge the support of National Natural Science Foundation of China (Grant no. 11672236), Aeronautical Science Foundation of China (Grant no. 2015ZA53013), and the Fundamental Research Funds for the Central Universities (Grant no. 3102016HKJJ004).