Identification of Lateral-Directional Aerodynamic Parameters for Aircraft Based on a Wind Tunnel Virtual Flight Test

Abstract

In the early stages of aircraft design, a scaled model of the aircraft is installed in a wind tunnel for dynamic semi-free flight to approximate real flight, and the test data are then used to identify the aerodynamic parameters. However, the absence of the translational motion of the test model makes its flight dynamics different from those in free flight, and the effect of this difference on parameter identification needs to be investigated. To solve this problem, a 3-DOF wind tunnel virtual flight test device is built to fix the test model on a rotating mechanism, and the model is free to rotate in three axes through the deflection of the control surfaces. The flight dynamics equations for the wind tunnel virtual flight test are established and expressed as a decoupled form of the free flight force and the influence of the test support frame force on the model’s motions through linearization. The differences between wind tunnel virtual flight and free flight are analysed to develop a model for the identification of aerodynamic parameters. The selection of the lateral-directional excitation signal and the design method of its parameters are established based on the requirements for the identifiability of the aerodynamic derivatives, and a step-by-step method for the identification of aerodynamic force and moment derivatives is established. The aerodynamic parameter identification results of a blended wing body aircraft show that the identification method proposed in this paper can obtain results with high accuracy, and the response of the modified motion model is consistent with that of the free flight motion model.

1. Introduction

The identification of aerodynamic parameters and the study of the dynamics of aircraft based on actual test flights of real aircraft is costly and risky [1,2]. In the early stage of aircraft design, a wind tunnel virtual flight test is carried out to fly the test model dynamically in the wind tunnel, and then complete the identification of aerodynamic parameters and compare and correct with the traditional wind tunnel force measurement data [3,4].

At present, the most widely used wind tunnel virtual flight test is a dynamic flight test with three degrees of freedom [5]. As shown in Figure 1, in the wind tunnel virtual flight test, the test model is connected with the support rod through a three-degree-of-freedom rotating mechanism and installed in the wind tunnel test section, so that the linear displacement of the model is fixed, but it has three angular motion degrees of freedom. Open-loop and closed-loop control of the test model are achieved by driving the control surface directly or by using the flight control law to drive the control surface, respectively [6].

The traditional wind tunnel test and aerodynamic modelling technology are based on the linear superposition principle [7], as shown in Figure 2. Static wind tunnel tests are used to obtain constant aerodynamic data, small amplitude forced oscillation tests are used to obtain dynamic derivative data, and the large-amplitude forced oscillation test is used to obtain aerodynamic data at a large angle of attack, so that the aerodynamic model of the aircraft can be finally constructed [8].

Traditional wind tunnel force measurements are sometimes inaccurate. The traditional wind tunnel test technique decomposes the parameters of aircraft motion and attitude to study them, which separates the interaction and influence between aerodynamics and motion, and the test results obtained cannot fully reflect the aerodynamic characteristics during real flight. Particularly in the large angle of attack, the aircraft’s translation, rotation, and control surface deflection are coupled with each other, when the traditional wind tunnel measurement results often differ from the aerodynamic model of real flight. Data from wind tunnel virtual flight tests can be used to identify the aerodynamic parameters and to compare and correct them with traditional wind tunnel force measurements.

Traditional wind tunnel tests mainly consist of static force measurement tests or forced oscillation tests, where the aircraft cannot move autonomously. In the wind tunnel virtual flight test, the motion of the model is realised by manipulating the surface deflection, and the flight control computer can be installed in the aircraft to control the motion of the model through the closed-loop flight control law [9], so it can realise integrated aerodynamic/motion/control system research at the same time.

The wind tunnel virtual flight test fills the gap between the subscale model actual test flights and the traditional wind tunnel force measurement test [10] and can simulate flight conditions more realistically before the actual test flights, which has a wide background of applications in aerodynamic parameter identification, flight control law verification, and stability evaluation [11].

Research on the identification of aerodynamic parameters based on 3-DOF wind tunnel virtual flight tests has been widely performed [12]. Cranfield University conducted open-loop maneuvering tests for the BAe Hawk aircraft, which is fixed to a 3-DOF rotating mechanism with belly support, to identify the aerodynamic derivatives in the longitudinal- and lateral-directional linearization matrix of the aircraft [13,14]. However, it does not take into account the differences in the equations of motion between wind tunnel virtual flight and free flight due to displacement constraints [4], and this difference needs to be analysed to propose a correction method for the identification model to better simulate the motion of free flight.

In addition, the test model does not have side displacements, and therefore, the side force derivatives cannot be identified directly from the motion of the model. References [11,15] proposed installing a strain gauge balance in the test device to measure the aerodynamic force of the model in real time, but did not explain the specific identification method for side force derivatives.

In order to build the aerodynamic parameter identification model for wind tunnel virtual flight tests, a flight dynamics model is first needed. The common practice currently is to omit the aircraft’s equations of translational motion and only consider its rotational motion [5]. This approach does not take the forces of the support device into account in the equations of motion and is not conducive to analysing the differences between wind tunnel virtual flight and free flight dynamics of the test model.

In order to ensure the accurate identification of aerodynamic derivative, it is necessary to design the appropriate frequency range and amplitude of the excitation signal. If the frequency and amplitude of the signal are inappropriate, the lateral-directional modal characteristics of the aircraft cannot be fully excited, thus reducing the accuracy of the identification results [16]. The existing excitation signal design methods have two kinds: one is to establish a simulation model and calculate the time domain to analyse the response size of each aerodynamic force or moment component in the process of aircraft motion, so as to ensure that the aerodynamic force or moment magnitude caused by the aircraft motion variables (β, p, r, ϕ) and control variables (δ_a, δ_r) meet the identifiable requirements, to determine the frequency and amplitude range of the signal [17,18]; and the accuracy of the identification results of various aerodynamic derivatives can be guaranteed by designing complex orthogonal multi-sine optimization excitation signals [19,20]. These two methods require a large amount of calculation and complicated steps, and they are not suitable for the excitation signal design of aircraft with closed-loop flight control law, because the flight control law will have a feedback effect on the initial signal input, resulting in a large difference between the final control surface deflection and the initial input signal.

To solve the above problems, a method for identifying the lateral-directional aerodynamic parameters of an aircraft based on wind tunnel virtual flight tests is proposed in this paper. First, the flight dynamics equations for the wind tunnel virtual flight test are derived and decomposed into two parts, the free flight motion and the additional motion resulting from displacement constraints by linearizing the equations of motion. The differences in the lateral-directional motion characteristics between wind tunnel virtual flight and free flight are analysed, and the identification model of the aerodynamic parameters is established. Second, the Bode diagram analysis of the motion model is used to establish a method for the selection of the lateral-directional excitation signal and the design of its parameters to ensure the identification accuracy of each aerodynamic derivative. Third, for the characteristics of wind tunnel virtual flight tests where the aerodynamic forces of the model can be measured by a force balance, a step-by-step identification method of side force derivatives, yawing, and rolling moment derivatives is proposed. Finally, the lateral-directional aerodynamic derivatives are identified according to the test data, and the motion model of wind tunnel virtual flight is modified. The correctness of the method established in this paper is verified by comparing the response of the modified motion model with that of the free flight motion model.

By decoupling the equations of motion of the aircraft, this paper innovatively establishes a wind tunnel virtual flight motion model, considering the influence of the support device, and visually analyses the differences between it and the free flight motion model. The accurate free flight motion model can be obtained by correcting the wind tunnel virtual flight motion model obtained from the identification. Based on this idea, this paper proposes a step-by-step identification method of aerodynamic parameters based on wind tunnel virtual flight, an excitation signal design method for closed-loop control aircraft based on amplitude–frequency characteristic analysis, and a correction method of identification model.

2. Flight Dynamics Model for Wind Tunnel Virtual Flight Test

The blended wing body (BWB) aircraft studied in this article refers to a new layout aircraft that eliminates the tail and uses a V-tail, instead of horizontal and vertical tails [21], as shown in Figure 1. This configuration can effectively improve the lift-to-drag ratio and fuel efficiency of the aircraft [22]. The V-tail replaces the horizontal and vertical tails, the same direction deflection of the V-tail is used as the elevator δ_e during pitch control, and the reverse deflection is used as the rudder δ_r during yaw control, thus decoupling the operation between longitudinal and directional. Define the V-tail control surface downwards deviation as positive, upwards deviation as negative, and V-tail surface angles of the left and right sides as δ_V−L and δ_V−R, respectively. The equivalent elevator and rudder angles are defined as δ_e = 0.5 × (δ_V−L + δ_V−R) and δ_r = 0.5 × (δ_V−L − δ_V−R).

If the aerodynamic derivatives over the full envelope of the aircraft are to be identified through a wind tunnel virtual flight test, the attitude angle of the test model needs to be precisely controlled, so that it varies within the desired range. The open-loop manipulation does not accurately control the attitude of the model during wind tunnel tests, and the model may deviate from the original equilibrium point after the excitation signal input. Therefore, it is necessary to design control augmentation laws for the BWB test model, carry out control surface excitation for the closed-loop controlled model, and complete the study of aerodynamic parameter identification for wind tunnel virtual flight tests.

The support system is shown in Figure 3, and the connection point of the 3-DOF rotation mechanism coincides with the centre of mass of the model. The support mechanism allows the aircraft to rotate freely in 3-DOF around the centre of mass with no translational motion. The range of permitted rotation is ±45° in pitch and roll and ±180° in yaw, and it can rotate to the maximum angle at the same time. The sequence of rotations follows the conventional definition of Euler angles.

The dimensional information of the test model is shown in Figure 4. The test model is flown dynamically in a wind tunnel with airborne sensors mounted on the model to measure aircraft motion parameters and using the test measurement data for aerodynamic parameter identification. In real flight tests, the aerodynamic forces on an aircraft can be indirectly derived by measuring the triaxial acceleration at the center of mass of the aircraft. While in the wind tunnel virtual flight test, a strain gauge balance is installed at the connection between the three-axis rotating mechanism and the centre of mass of the aircraft, allowing the aerodynamic forces of the model and the forces exerted on the model by the support devices to be collected in real time. The data on the longitudinal and lateral-directional motion that can be directly measured during wind tunnel virtual flight tests are shown in

Table 1. Measurable data from wind tunnel virtual flight test.

Parameters	Description	Instruments
ϕ, θ, ψ	Roll angle, pitch angle, yaw angle	Inertial measurement unit
p, q, r	Roll rate, pitch rate, yaw rate	Inertial measurement unit
α, β	Angle of attack and sideslip	Numerical solution
δa, δr	Aileron, rudder deflection	Rotary encoder
Fx, Fy, Fz	Force of support device in the x-, y-, and z-axis directions (body coordinate system)	Strain gauge balance

.

2.1. Flight Dynamics Modelling

Generally, the aircraft is symmetrical about the Ox_bz_b plane under the body coordinate system, i.e., I_xy = I_yz = 0. The lateral-directional rotation dynamics and kinematic equations of the test model around the centre of mass are the same as free flight, as shown in Equation (1) [23].

$$\{\begin{array}{l}L\hspace{0.17em}=\hspace{0.17em}{I}_x\dot{p}-I_{xz}\left(\dot{r}+pq\right)-\left(I_y-I_z\right)qr\\ N=I_z\dot{r}-I_{xz}\left(\dot{p}-qr\right)-\left(I_x-I_y\right)pq\\ \dot{\varphi }=p+q\sin \varphi \tan \theta +r\cos \varphi \tan \theta \\ \dot{\psi }=q\sin \varphi \sec \theta +r\cos \varphi \sec \theta \end{array}$$

(1)

where I_x, I_y, and I_z are the inertia of the three axes; I_xz is the cross inertia; p, q, and r are the roll rate, pitch rate and yaw rate, respectively; and L and N are the rolling and yawing moments of the aircraft, respectively.

The lateral-directional dynamics equation describing the centroid of an aircraft in the body coordinate system is shown in Equation (2) [24]:

$$m\left(\frac{dv}{dt}+ru-pw\right)=F_{y-b}=0$$

(2)

where m is the mass of the model; u, v, and w are the components of airspeed in the body coordinate system; and F_y−b is the combined external force on the y-axis of the aircraft under the body coordinate system. Since the model has no translational motion, the combined external force F_y−b = 0.

In the wind tunnel virtual flight test, the external forces on the model include gravity, the triaxial aerodynamic forces, and the force exerted by the support rod on the model. Project all external forces under the body coordinate system, as shown in Equation (3).

$$\left[\begin{array}{c}\backslash \\ F_{y-b}\\ \backslash \end{array}\right]=\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]+L_{bg}·\left[\begin{array}{c}0\\ 0\\ mg\end{array}\right]+\left[\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right]$$

(3)

where L_bg is the rotation matrix from the ground coordinate system to the body coordinate system; F_x, F_y, and F_z are the forces of the support rod on the model in the x-, y-, and z-axis directions, respectively; and X, Y, and Z are the aerodynamic forces on the model in the body coordinate system, respectively.

Substituting the combined external force F_y−b of the model in the body coordinate system derived from Equation (3) into Equation (2), the centre-of-mass dynamics equation of the model is obtained as shown in Equation (4).

$$m\left(\frac{dv}{dt}+ru-pw\right)=Y+mg\sin \varphi \cos \theta +F_y$$

(4)

The angle of attack and sideslip, which are important variables in wind tunnel tests, are usually measured by vane sensors. A wind vane sensor needs to be installed at the nose of the aircraft, which may lead to a decrease in the accuracy of the airflow angle measurement, due to the interference of the fuselage. In real flight test, it is necessary to install multiple sensors on the fuselage to measure the flow angle of the aircraft [25]. However, in the wind tunnel virtual flight test, they can be directly solved by the three-axis attitude angles ϕ, θ, ψ. Due to the high accuracy of the attitude angle measurement, the angle of attack and sideslip are also solved with high accuracy. According to the definition of the angle of attack and sideslip, their calculation formulas are as follows [26]:

$$\{\begin{array}{l}\alpha ={\tan }^{-1}(\frac{w}{u})\\ \beta ={\sin }^{-1}(\frac{v}{V})\end{array}$$

(5)

As shown in Figure 5, define the ground coordinate system with the origin O_g at the centre of mass of the model, with Ox_g pointing in the direction of wind velocity, Oz_g vertically downwards, and Oy_g perpendicular to the Ox_gz_g plane to the right. In the wind tunnel test, assuming that the direction of the incoming flow is constant, the x-axis of the air coordinate system, ground coordinate system, and flight path coordinate system all coincide. Thus, u, v, and w can be solved using the attitude angle and the coordinate system transformation matrix, which, in turn, directly solves for the airflow angles α and β, as shown in Equation (6).

$$\{\begin{array}{l}\left[\begin{array}{l}u\\ v\\ w\end{array}\right]=L_{bg}\left[\begin{array}{l}V\\ 0\\ 0\end{array}\right]=\left[\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]\left[\begin{array}{l}V\\ 0\\ 0\end{array}\right]\\ a_{11}=\cos \theta \cos \psi \\ a_{21}=\sin \varphi \sin \theta \cos \psi -\cos \varphi \sin \psi \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ a_{31}=\cos \varphi \sin \theta \cos \psi +\sin \varphi \sin \psi \hspace{0.17em}\\ \alpha =\arctan (w/u)=\arctan (a_{31}/a_{11})\\ \beta =\arcsin (v/V)=\arcsin (a_{21}/V)\end{array}$$

(6)

2.2. Linearization and Decoupling for Equations of Motion

The motion response of the wind tunnel virtual flight for the test model can be decoupled as the superposition of two parts of motion, namely (1) the motion of unconstrained free flight, and (2) the effect on the motion of the model by the force of the support device on the model. Through linearization, the lateral-directional motion equation can be described as the decoupling form of the effects of free flight aerodynamic force, gravity, and supporting force on the model motion.

In the linearization of the equations of motion, the amount of side force, rolling, and yawing moment perturbations, relative to the equilibrium flight motion, can be expressed as linear relationships, with ΔY, ΔL, and ΔN expressed as Equation (7): [24].

$$\{\begin{array}{l}\Delta Y=Y_{\beta }\Delta \beta +Y_r\Delta r+Y_p\Delta p+Y_{\delta r}\Delta \delta r\\ \Delta L=L_{\beta }\Delta \beta +L_p\Delta p+L_r\Delta r+L_{\delta a}\Delta \delta a+L_{\delta r}\Delta \delta r\\ \Delta N=N_{\beta }\Delta \beta +N_p\Delta p+N_r\Delta r+N_{\delta r}\Delta \delta r\end{array}$$

(7)

As the displacement of the test model is constrained, the combined force of gravity, aerodynamic forces, and the forces acting on the model by the support device are zero. Making the combined force F_y−b of the model in Equation (3) be 0, the expression of the force acting on the model by the support device F_y can be obtained:

$$F_y=-Y-mg\sin \varphi \cos \theta$$

(8)

We linearized Equation (8) with a base state parameter of ϕ_* = 0°, and the result of the linearization is shown in Equation (9).

$$\begin{array}{l}\Delta F_y=-\Delta Y-mg\cos {\theta }_{*}\Delta \varphi \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=-Y_{\beta }\Delta \beta -Y_p\Delta p-Y_r\Delta r-mg\cos {\theta }_{*}\Delta \varphi -Y_{\delta r}\Delta \delta r\end{array}$$

(9)

Equation (9) shows that ΔF_y can be expressed in the form of a linear superposition of aerodynamic forces and gravity. The linear expression for the force ΔF_y acting on the model by the support device can be obtained from Equation (9) as:

$$\Delta F_y=F_{y-}{}_{\beta }\Delta \beta +F_{y-}{}_p\Delta p+F_{y-}{}_r\Delta r+F_{y-}{}_{\varphi }\Delta \varphi +F_{y-}{}_{\delta r}\Delta \delta r$$

(10)

where F_y−β, F_y−p, F_y−r, F_y−ϕ, and F_y−δr denote the derivatives of the side force generated by the support device in the body coordinate system, with respect to the sideslip angle, roll rate, yaw rate, roll angle, and rudder deflection, respectively; the specific expressions are shown in Equation (11).

$$\{\begin{array}{l}{F}_{y-}{}_{\beta }=-Y_{\beta }\hspace{1em}\hspace{1em}{F}_{y-}{}_{\varphi }=-mg\cos {\theta }_{*}\\ F_{y-}{}_p=-Y_p\hspace{1em}\hspace{1em}{F}_{y-}{}_{\delta r}=-Y_{\delta r}\\ F_{y-}{}_r=-Y_r\end{array}$$

(11)

The lateral-directional equations of motion for the wind tunnel virtual flight shown in Equations (1) and (4) are linearised, and the expressions ΔY, ΔL, ΔN, and ΔF_y, shown in Equations (7) and (10), are substituted into these linearised equations of motion to obtain the linearization results for the test model at the base state of V_* and α_*, as shown in Equation (12):

$$\begin{array}{l}\left[\begin{array}{c}\Delta \dot{\beta }\\ \Delta \dot{p}\\ \Delta \dot{r}\\ \Delta \dot{\varphi }\end{array}\right]=\underset{}{\underbrace{\left[\begin{array}{cccc}{\overline{Y}}_{\beta }& {\alpha }_{*}+{\overline{Y}}_p& {\overline{Y}}_r-1& g\cos {\theta }_{*}/V_{*}\\ {\overline{L}}_{\beta }& {\overline{L}}_p& {\overline{L}}_r& 0\\ {\overline{N}}_{\beta }& {\overline{N}}_p& {\overline{N}}_r& 0\\ 0& 1& \tan {\theta }_{*}& 0\end{array}\right]}}\left[\begin{array}{c}\Delta \beta \\ \Delta p\\ \Delta r\\ \Delta \varphi \end{array}\right]+\underset{}{\underbrace{\left[\begin{array}{cc}0& {\overline{Y}}_{{\delta }_r}\\ {\overline{L}}_{{\delta }_a}& {\overline{L}}_{{\delta }_r}\\ {\overline{N}}_{{\delta }_a}& {\overline{N}}_{{\delta }_r}\\ 0& 0\end{array}\right]}}\left[\begin{array}{c}{\delta }_a\\ {\delta }_r\end{array}\right]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{A}_1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{B}_1\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\underset{}{\underbrace{\left[\begin{array}{cccc}-{\overline{Y}}_{\beta }& -{\overline{Y}}_p& -{\overline{Y}}_r& -g\cos {\theta }_{*}/V_{*}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]}}\left[\begin{array}{c}\Delta \beta \\ \Delta p\\ \Delta r\\ \Delta \varphi \end{array}\right]+\underset{}{\underbrace{\left[\begin{array}{cc}0& -{\overline{Y}}_{{\delta }_r}\\ 0& 0\\ 0& 0\\ 0& 0\end{array}\right]}}\left[\begin{array}{c}{\delta }_a\\ {\delta }_r\end{array}\right]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{A}_2\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{B}_2\end{array}$$

(12)

where A₁ and B₁ represent the free flight lateral-directional stability matrix and the control matrix of the control surface, respectively, and A₂ represents the matrix of the additional motion of the model caused by the support force. In addition, the deflection of the V-tail control surface generates side force and the constraints of the support device have an impact on the control matrix, as shown in B₂. The expressions for each lateral-directional dynamic derivative in Equation (12) are shown in

Table 2. Lateral-directional dynamic derivatives.

Side Force Derivatives	Yawing Moment Derivatives	Rolling Moment Derivatives
$\stackrel{-}{Y}$β = $\frac{Y_{\beta }}{m{V}_{*}}$ = $\frac{C_{y\beta }{\stackrel{-}{q}}_{*}S}{m{V}_{*}}$	$N_{\beta }$ = $C_{n\beta }{\stackrel{-}{q}}_{*}S$b	$L_{\beta }$ = $C_{l\beta }{\stackrel{-}{q}}_{*}S$b
$\stackrel{-}{Y}$p = $\frac{Y_p}{m{V}_{*}}$ = $\frac{C_{yp}{\stackrel{-}{q}}_{*}S}{m{V}_{*}}\frac{b}{2V_{*}}$	$N_p$ = $C_{np}{\stackrel{-}{q}}_{*}S$b$\frac{b}{2V_{*}}$	$L_p$ = $C_{lp}{\stackrel{-}{q}}_{*}S$b$\frac{b}{2V_{*}}$
$\stackrel{-}{Y}$r = $\frac{Y_r}{m{V}_{*}}$ = $\frac{C_{yr}{\stackrel{-}{q}}_{*}S}{m{V}_{*}}\frac{b}{2V_{*}}$	$N_r$ = $C_{nr}{\stackrel{-}{q}}_{*}S$b$\frac{b}{2V_{*}}$	$L_r$ = $C_{lr}{\stackrel{-}{q}}_{*}S$b$\frac{b}{2V_{*}}$
$\stackrel{-}{Y}$δa = $\frac{Y\delta a}{m{V}_{*}}$ = $\frac{{C_{y\delta a}\stackrel{-}{q}}_{*}S}{m{V}_{*}}$	$N_{\delta a}$ = $C_{n\delta a}{\stackrel{-}{q}}_{*}S$b	$L_{\delta a}$ = $C_{l\delta a}{\stackrel{-}{q}}_{*}S$b
$\stackrel{-}{Y}$δr = $\frac{Y\delta r}{m{V}_{*}}$ = $\frac{{C_{y\delta r}\stackrel{-}{q}}_{*}S}{m{V}_{*}}$	$N_{\delta r}$ = $C_{n\delta r}{\stackrel{-}{q}}_{*}S$b	$L_{\delta r}$ = $C_{l\delta r}{\stackrel{-}{q}}_{*}S$b
$\stackrel{-}{L}$I = $\frac{L_i+(I_{xz}/I_z)N_i}{I_x-I_{xz}^2/I_z}$, $\stackrel{-}{N}$i = $\frac{N_i+(I_{xz}/I_z)L_i}{I_z-I_{xz}^2/I_x}$, i∈(β, p, r, δa, δr)

.

To simulate the free flight motion of the aircraft through the results of aerodynamic parameter identification from the wind tunnel virtual flight test, the effects of the A₂ and B₂ matrices in Equation (12) need to be eliminated. For ease of description, Equation (12) is rewritten into the form shown in Equation (13). The mathematical expressions in the box are the differences in the lateral-directional equations of motion between the wind tunnel virtual flight and the free flight.

$$\begin{array}{l}\left[\begin{array}{c}\Delta \dot{\beta }\\ \Delta \dot{p}\\ \Delta \dot{r}\\ \Delta \dot{\varphi }\end{array}\right]=\underset{}{\underbrace{\left[\begin{array}{cccc}{\overline{Y}}_{\beta }+\boxed{-{\overline{Y}}_{\beta }}& {\alpha }_{*}+{\overline{Y}}_p+\boxed{-{\overline{Y}}_p}& {\overline{Y}}_r-1+\boxed{-{\overline{Y}}_r}& \frac{g\cos {\theta }_{*}}{V_{*}}+\boxed{-\frac{g\cos {\theta }_{*}}{V_{*}}}\\ {\overline{L}}_{\beta }& {\overline{L}}_p& {\overline{L}}_r& 0\\ {\overline{N}}_{\beta }& {\overline{N}}_p& {\overline{N}}_r& 0\\ 0& 1& \tan {\theta }_{*}& 0\end{array}\right]}}\left[\begin{array}{c}\Delta \beta \\ \Delta p\\ \Delta r\\ \Delta \varphi \end{array}\right]+\underset{}{\underbrace{\left[\begin{array}{cc}0& {\overline{Y}}_{{\delta }_r}+\boxed{-{\overline{Y}}_{{\delta }_r}}\\ {\overline{L}}_{{\delta }_a}& {\overline{L}}_{{\delta }_r}\\ {\overline{N}}_{{\delta }_a}& {\overline{N}}_{{\delta }_r}\\ 0& 0\end{array}\right]}}\left[\begin{array}{c}\Delta {\delta }_a\\ \Delta {\delta }_r\end{array}\right]\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}A\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}B\end{array}$$

(13)

The lateral-directional aerodynamic parameter identification model for the wind tunnel virtual flight test can be established according to Equation (13). As the displacement of the model is constrained, it is not possible to identify the side force derivatives of the model directly from the motion of the model; therefore, a step-by-step identification method for the side force derivatives, rolling, and yawing moment derivatives needs to be established.

3. Design Method for Excitation Signals Based on Frequency Domain Analysis

In wind tunnel virtual flight tests, displacement constraints lead to changes in the dynamics of the model, and appropriate excitation signals need to be designed to fully excite the aircraft’s dynamic characteristics, thus ensuring the accuracy of the identification of the lateral-directional aerodynamic derivatives. This paper establishes a method for the selection of excitation signal types and their parameter design based on frequency domain analysis.

3.1. Selection of Excitation Signal Type

The traditional design method of the excitation signal is to solve for the optimal excitation signal by minimizing the deviation or covariance of the aerodynamic parameter identification results [16]. For open-loop controlled aircraft, the above excitation signal design methods are applicable. However, for aircraft with a control augmentation law, since the above signal design method does not consider the feedback effect of the flight control law (as shown in Figure 6), the flight control law may have a significant effect on the excitation signal [4], so it cannot guarantee that the actual maneuvering signal of the aircraft can adequately excite the motion characteristics of the aircraft.

This paper proposes a method for excitation signal type selection based on frequency domain analysis. The most commonly used excitation signals include dipole square wave, 3211-multipole square wave, sine wave, and frequency sweep signals. The dipole square and sine wave signals have a single frequency range and are relatively unsuitable as excitation signals. The 3211-multipole square wave and frequency sweep signals are loaded on the control surface of the test model, and the data before and after the action of the excitation signals are analysed in the time domain and frequency domain to select the signal type that is less influenced by the flight control law.

Taking directional manipulation as an example, by inputting 3211 and frequency sweep excitation signals to the rudder, the waveform changes of the signals are observed after the feedback effect of the flight control law, and the typical results are shown in Figure 7.

As seen from Figure 7, the flight control law has a large effect on the 3211-multipole square wave signal, with a large irregular change in the waveform, while it has a small effect on the frequency sweep signal, with a slight change in frequency. The frequency spectrum of the time-domain changes of the above two types of excitation signals is analysed by fast Fourier transform, and the frequency band changes of the signal energy distribution are observed, as shown in Figure 8.

From Figure 8, it can be seen that after the feedback influence of the flight control law, the spectral shape of the 3211 signal changes significantly and the frequency range of the energy distribution of the signal becomes narrower, while the spectral shape of the frequency sweep signal changes to a relatively small extent and the frequency range of the energy distribution of the signal does not change significantly.

According to the above analysis, when a 3211-multipole square wave is input into the rudder of an aircraft with a closed-loop flight control law, the final rudder signal will have a large degree of irregular change, compared with the initial input excitation signal. However, if the sweep signal is input, the time domain waveform of the signal does not change much, and the frequency band width and position of the signal energy distribution do not change significantly. In addition, the long duration of the sweep signal allows for a continuous input signal to the control surface to excite the desired motion pattern of the aircraft. These laws also apply to the excitation signals of the ailerons. In summary, it is recommended to select the frequency sweep signal as the excitation signal for lateral-directional aerodynamic parameter identification.

3.2. Design of Excitation Signal Parameters

The design parameters of the frequency sweep signal include the low frequency limit ω_l, the high frequency limit ω_h, and the signal amplitude |A|. In this paper, a method of excitation signal parameter design based on frequency domain analysis is proposed to design the frequency band and amplitude magnitude of the excitation signal. Bode diagram analysis of the equations of motion of the aircraft is carried out to observe the amplitude response of each parameter to be identified with the change in excitation signal frequency. The frequency range of the excitation signal is determined by ensuring that the amplitude response of the aerodynamic forces or moments caused by each motion variable (β, p, r, ϕ), as well as the manipulation variables (δ_a, δ_r), is sufficiently large. Finally, the amplitude of the excitation signal is adjusted, so that the signal has high energy in the designed frequency band.

Equation (14) is the observation equation for the lateral-directional motion of the aircraft, which contains the observation equations for the side force, rolling moment, and yawing moment. The Bode diagram analysis of the lateral-directional equations of motion shown in Equation (13) is carried out to observe the amplitude response of the component of side force, rolling moment, and yawing moment in Equation (14), with a change of rudder and aileron frequency.

$$\{\begin{array}{l}\Delta Y=Y_{\beta }\Delta \beta +Y_r\Delta r+Y_{\delta r}\Delta {\delta }_r\\ \Delta \dot{p}={\overline{L}}_{\beta }\Delta \beta +{\overline{L}}_p\Delta p+{\overline{L}}_r\Delta r+{\overline{L}}_{{\delta }_a}\Delta {\delta }_a+{\overline{L}}_{{\delta }_r}\Delta {\delta }_r\\ \Delta \dot{r}={\overline{N}}_{\beta }\Delta \beta +{\overline{N}}_p\Delta p+{\overline{N}}_r\Delta r+{\overline{N}}_{{\delta }_r}\Delta {\delta }_r\end{array}$$

(14)

Taking directional manipulation as an example, the magnitude response curves $\left|{\stackrel{-}{N}}_{\beta }\Delta \beta \left(\omega \right)/\delta r\right|$, $\left|{\stackrel{-}{N}}_p\Delta p\left(\omega \right)/\delta r\right|$, $\left|{\stackrel{-}{N}}_r\Delta r\left(\omega \right)/\delta r\right|$, and $\left|{\stackrel{-}{N}}_{\delta r}\Delta \delta r\left(\omega \right)/\delta r\right|$ of ${\stackrel{-}{N}}_{\beta }$Δβ, ${\stackrel{-}{N}}_p$Δp, ${\stackrel{-}{N}}_r$Δr and ${\stackrel{-}{N}}_{\delta r}\Delta \delta r$ with the change in rudder frequency are drawn according to the yawing moment observation equation shown in Equation (14). The amplitude response curves $\left|{\stackrel{-}{L}}_{\beta }\Delta \beta \left(\omega \right)/\delta r\right|$, $\left|{\stackrel{-}{L}}_p\Delta p\left(\omega \right)/\delta r\right|$, $\left|{\stackrel{-}{L}}_r\Delta r\left(\omega \right)/\delta r\right|,$ and $\left|{\stackrel{-}{L}}_{\delta r}\Delta \delta r\left(\omega \right)/\delta r\right|$ of ${\stackrel{-}{L}}_{\beta }\Delta \beta$, ${\stackrel{-}{L}}_p\Delta p$, ${\stackrel{-}{L}}_r\Delta r$, and ${\stackrel{-}{L}}_{\delta r}\Delta \delta r$ with the change in rudder frequency are drawn according to the rolling moment equation, drawing the amplitude response curves |Y_βΔβ(ω)/δ_r|, |Y_rΔr(ω)/δ_r|, and |Y_δrΔδ_r(ω)/δ_r| of Y_βΔβ, Y_rΔr, and Y_δrΔδ_r with the variation of rudder frequency according to the side force equation.

Figure 9 shows the amplitude response curves of each yawing moment component with the change in rudder frequency. In the same way, the amplitude–frequency characteristic curves of the rolling moment and side force components with the frequency of the excitation signal can be obtained. The frequency corresponding to the peak of the resonance is approximately equal to the Dutch roll mode frequency of the aircraft. When the excitation signal frequency is lower or higher than a certain range, the response amplitude of each yawing moment component decreases, and the identification accuracy of all yawing moment derivatives decreases.

To ensure the identifiability of each yawing moment derivative, the magnitude-frequency response magnitudes of the yawing moment components ${\stackrel{-}{N}}_{\beta }\Delta \beta$, ${\stackrel{-}{N}}_r\Delta r$, ${\stackrel{-}{N}}_p\Delta p$, and ${\stackrel{-}{N}}_{\delta r}\Delta \delta r$ caused by the sideslip change Δβ, yaw rate change Δr, roll rate change Δp, and rudder deflection Δδ_r should not differ significantly. The time domain curves of each yawing moment component, with different frequencies of rudder signal inputs, are given in Figure 10.

In Figure 10a, when the excitation signal frequency is too low, the response of the angular acceleration of the yaw angle Δ$\dot{r}$ is almost zero. In this case, the static stability moment generated by ${\stackrel{-}{N}}_{\beta }\Delta \beta$ and the manipulation moment generated by ${\stackrel{-}{N}}_{\delta r}\Delta \delta r$ almost cancel each other out, resulting in a small yawing moment, which makes it difficult to excite the yaw rate response; therefore, it is not possible to accurately identify all the yawing moment derivatives. In Figure 10c, when the excitation signal frequency is too high, the aircraft response is much slower than the manipulation speed, resulting in small changes in the yaw rate, roll rate, and sideslip, leading to small response amplitudes in ${\stackrel{-}{N}}_{\beta }\Delta \beta$, ${\stackrel{-}{N}}_p$Δp, and ${\stackrel{-}{N}}_r\Delta r$. Therefore, all yawing moment derivatives cannot be accurately identified. In Figure 10b, when the excitation signal frequency is appropriate, the amplitude response of all yawing moment components are on the same order of magnitude, and the identification results will be accurate. In summary, to ensure the identification accuracy, the excitation signal must be in a certain frequency range.

In the same way, for the rolling moment observation equation shown in Equation (14), the magnitude response magnitudes of the rolling moment components ${\stackrel{-}{L}}_{\beta }\Delta \beta$, ${\stackrel{-}{L}}_r\Delta r$, ${\stackrel{-}{L}}_p\Delta p$, and ${\stackrel{-}{L}}_{\delta r}\Delta \delta r$ caused by the sideslip, yaw rate, roll rate, and rudder surface deflection should not differ significantly. The magnitude response magnitudes of each side force component Y_βΔβ, Y_rΔr, and Y_δrΔδ_r in the side force observation equation should not differ significantly.

In page 44 of reference [27], it is stated that, when the component of a certain aerodynamic force or moment occupies at least 10% of the maximum aerodynamic force or moment component, the aerodynamic derivative corresponding to the component is considered to be identifiable. To further improve the accuracy of the identification results, a more stringent requirement is used, i.e., when the component of an aerodynamic force or moment accounts for at least 10% of the full aircraft aerodynamic force or moment, the aerodynamic derivative corresponding to that component is considered identifiable. Therefore, the sum of the response amplitude of each aerodynamic force and moment component at all frequencies is calculated in the Bode diagram, and the value corresponding to the total response amplitude 1/10 is calculated, as shown in Equation (15).

$$\{\begin{array}{l}{\left|A\right|}_{Cn-{\omega }_i}=\frac{1}{10}{\displaystyle \sum _{{\omega }_i}\left({\left|A\right|}_{{\overline{N}}_{\beta }\Delta \beta }+{\left|A\right|}_{{\overline{N}}_p\Delta p}+{\left|A\right|}_{{\overline{N}}_r\Delta r}+{\left|A\right|}_{{\overline{N}}_{{\delta }_r}\Delta \delta r}\right)}\\ {\left|A\right|}_{Cl-{\omega }_i}=\frac{1}{10}{\displaystyle \sum _{{\omega }_i}\left({\left|A\right|}_{{\overline{L}}_{\beta }\Delta \beta }+{\left|A\right|}_{{\overline{L}}_p\Delta p}+{\left|A\right|}_{{\overline{L}}_r\Delta r}+{\left|A\right|}_{{\overline{L}}_{{\delta }_r}\Delta \delta r}\right)}\\ {\left|A\right|}_{CY-{\omega }_i}=\frac{1}{10}{\displaystyle \sum _{{\omega }_i}\left({\left|A\right|}_{Y_{\beta }\Delta \beta }+{\left|A\right|}_{Y_{\delta r}\Delta \delta r}+{\left|A\right|}_{Y_r\Delta r}\right)}\end{array}$$

(15)

where ${\left|A\right|}_{{\stackrel{-}{N}}_{\beta }\Delta \beta }$, ${\left|A\right|}_{{\stackrel{-}{N}}_p\Delta p}$, ${\left|A\right|}_{{\stackrel{-}{N}}_r\Delta r}$, and ${\left|A\right|}_{{\stackrel{-}{N}}_{\delta r}\Delta \delta r}$ are the response amplitudes of each yawing moment component when the rudder signal frequency is ω_i; ${\left|A\right|}_{{\stackrel{-}{L}}_{\beta }\Delta \beta }$, ${\left|A\right|}_{{\stackrel{-}{L}}_p\Delta p}$, ${\left|A\right|}_{{\stackrel{-}{L}}_r\Delta r}$, and ${\left|A\right|}_{{\stackrel{-}{L}}_{\delta r}\Delta \delta r}$ are the response amplitudes of each component of the rolling moment; and ${\left|A\right|}_{Y_{\beta }\Delta \beta }$, ${\left|A\right|}_{Y_{\delta r}\Delta \delta r}$, and ${\left|A\right|}_{Y_r\Delta r}$ are the response amplitudes of each side force component. ${\left|A\right|}_{Cn-{\omega }_i}$, ${\left|A\right|}_{Cl-{\omega }_i}$, and ${\left|A\right|}_{CY-{\omega }_i}$ are 1/10 of the total response amplitudes of the yawing moment, rolling moment, and side force at a rudder signal frequency of ω_i, respectively.

According to Equation (15), 1/10 of the total response amplitudes of the yawing moment, rolling moment and side force at different rudder frequencies can be calculated and connected to a curve. The dashed line in Figure 9 is the demarcation line corresponding to 1/10 of the total response amplitude of the yawing moment, which is defined as the identifiable boundary of yawing moment derivatives corresponding to different rudder frequencies. Using the same method, identifiable boundaries for the rolling moment, as well as the side force derivatives, can also be obtained.

When the response amplitude of a force or moment corresponding to an aerodynamic derivative is above the identifiable boundary, it means that the derivative is identifiable at this frequency; otherwise, it is not. When the frequency of the rudder signal is in the range ω_δr·Cn$\in$[ω₁~ω₂] in Figure 9, i.e., the amplitude responses of all yawing moment components are above the identifiable boundary, all yawing moment derivatives are identifiable. Therefore, the frequency band of the rudder excitation signal determined by the yawing moment derivatives is in the range of ω_δr·Cn. The same method can be used to obtain the frequency band range of the rudder excitation signal determined by rolling moment, as well as side force derivatives, such as ω_δr·Cl and ω_δr·CY. To ensure that all lateral-directional aerodynamic derivatives are identifiable, it is necessary to take the intersection of ω_δr·Cn, ω_δr·Cl, and ω_δr·CY to obtain the frequency band of the rudder excitation signal ω_δr.

$${\omega }_{\delta r}={\omega }_{\delta r\cdot Cn}\cap {\omega }_{\delta r\cdot C_Y}\cap {\omega }_{\delta r\cdot C_l}$$

(16)

In the same way, the frequency band of the aileron excitation signal ω_δa can be obtained.

$${\omega }_{\delta a}={\omega }_{\delta a\cdot Cn}\cap {\omega }_{\delta a\cdot C_Y}\cap {\omega }_{\delta a\cdot C_l}$$

(17)

where ω_δa·Cn, ω_δa·Cl, and ω_δa·CY are the frequency band ranges of the aileron signals determined by the yawing moment, rolling moment, and side force derivatives, respectively.

The amplitude of the aileron and rudder signals need to be within the appropriate range. A large excitation signal amplitude will lead to a large change in the aircraft sideslip angle, thus introducing nonlinear aerodynamic effects; if the signal amplitude is too small, it will not be easy to excite the aircraft’s lateral-directional motion mode. Therefore, first, it is necessary to ensure that the variation in the aircraft sideslip angle is approximately ±2° to ensure that the aerodynamic derivatives are in the linear segment. Then, the amplitude of the excitation signal is determined according to the requirement of the change in sideslip angle under the consideration of the feedback effect of the flight control law.

4. Step-by-Step Identification Method for the Aerodynamic Parameters

4.1. Identification Method of Side Force Derivatives

In the wind tunnel virtual flight test, the position of the test model is fixed, and the combined external force of the whole aircraft is zero, so the side force derivatives of the model cannot be directly identified by the motion of the model. However, the test model has a strain gauge balance installed at the centre of mass to measure the side force of the whole aircraft in real time, so the measured value of the side force Y_m can be used as an observed variable to identify side force derivatives. The identification model of side force derivatives is shown in Equation (18).

$$Y_m=Y_{\beta }\Delta \beta +Y_p\Delta p+Y_r\Delta r+Y_{\delta r}\Delta {\delta }_r$$

(18)

where Y_m is the side force of the whole test model directly measured by the force balance; Δβ, Δp, Δr, and Δδr are the amount of change in sideslip angle, roll rate, yaw rate, and V-tail control surface, which can be measured directly by the sensors. The parameter to be identified is Θ = [Y_β, Y_p, Y_r, Y_δr].

Since the side force measured by the force balance is the actual aerodynamic force on the model without the force of the support device, the identification results based on the measured data are the side force derivatives of free flight.

For the side force observation equation shown in Equation (18), the side force derivatives can be directly identified using the least squares method. The general form of least squares identification is shown in Equation (19) [28]:

$$\{\begin{array}{l}y=X\Theta \\ z=X\Theta +v\end{array}$$

(19)

where Θ is a vector composed of the aerodynamic derivatives to be identified; y is the aerodynamic force of the theoretical output from the test model, i.e., the side force; z is the side force actually measured by the force balance; X is the matrix of independent variables composed of parameters, such as the motion variables (β, p, r, ϕ) and manipulation variables (δ_a, δ_r) of the aircraft; and v is the measurement noise matrix. In the wind tunnel test, the measurement accuracy of the strain balance is high, and v can be regarded as white noise with a mean value of 0, as shown in Equation (20):

$$E(v)=0\hspace{1em}\hspace{1em}E(v{v}^T)={\sigma }^{2}I$$

(20)

The least squares solutions of the aerodynamic lift and drag derivatives to be identified are:

$$\widehat{\Theta }={(X^{T}X)}^{-1}{X}^{T}z$$

(21)

Y_β, Y_p, Y_r, and Y_δr are identified by the least square method and substituted into Table 2 to solve the aerodynamic derivatives C_yβ, C_yp, C_yr, and C_yδr.

4.2. Identification Method of Rolling and Yawing Moment Derivatives

Selecting β, p, r, and ϕ as the state variables for the lateral-directional aerodynamic parameter identification, the state equation of the lateral-directional aerodynamic parameter identification model for wind tunnel virtual flight is shown in Equation (22) [16].

$$\{\begin{array}{l}\dot{p}=(c_{1}r+c_{2}p)q+c_{3}L+c_{4}N\\ \dot{r}=(c_{5}p-c_{2}r)q+c_{4}L+c_{6}N\\ \dot{\varphi }=p+\tan \theta \cdot (q\cdot \sin \varphi +r\cdot \cos \varphi )\\ \dot{\beta }=-r\cos \alpha +p\sin \alpha \end{array}$$

(22)

In Equation (22), the expressions of rolling and yawing moments L and N are shown in Equation (7), and the expressions of c₁, c₂, c₃, c₄, c₅, and c₆ are shown in Equation (23). The lateral-directional moment derivatives to be identified are Θ = [L_β, L_p, L_r, L_δa, L_δr, N_β, N_p, N_r, N_δr].

$$\{\begin{array}{l}{c}_1=\frac{(I_y-I_z)I_z-I_{xz}{}^2}{I_x{I}_z-I_{xz}{}^2}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_4=\frac{I_{xz}}{I_x{I}_z-I_{xz}{}^2}\\ c_2=\frac{(I_x-I_y+I_z)I_{xz}}{I_x{I}_z-I_{xz}{}^2}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_5=\frac{I_x(I_x-I_y)-I_{xz}{}^2}{I_x{I}_z-I_{xz}{}^2}\\ c_3=\frac{I_z}{I_x{I}_z-I_{xz}{}^2}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_6=\frac{I_x}{I_x{I}_z-I_{xz}{}^2}\end{array}$$

(23)

p_m, r_m, β_m, and ϕ_m measured by the experiment are chosen as the observation variables, so the observation equation is shown in Equation (24).

$$\{\begin{array}{l}{p}_m=p+v_p\\ r_m=r+v_r\\ {\beta }_m=\beta +v_{\beta }\\ {\varphi }_m=\varphi +v_{\varphi }\end{array}$$

(24)

where p, r, β, and ϕ are theoretical outputs and v_p, v_r, v_β, and v_ϕ are measurement noise. Aerodynamic derivatives of the test model can be obtained by identifying aerodynamic parameters of the nonlinear motion model shown in Equations (22)–(24). The most widely used method is the output error method based on maximum likelihood estimation [29]. General nonlinear dynamic equations can be expressed in the form shown in Equation (25):

$$\{\begin{array}{l}\dot{x}(t)\hspace{0.17em}\hspace{0.17em}=f[x(t),u(t),\Theta ]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x(t_0)=x_0\\ y(t)\hspace{0.17em}\hspace{0.17em}=g[x(t),u(t),\Theta ]\\ z(t_k)=y(t_k)+v(t_k)\end{array}$$

(25)

where x is the state variable (p, r, β, ϕ); u is the input variable (δ_a, δ_a); y is the output variable (p, r, β, ϕ); z is the observed variable (p_m, r_m, β_m, ϕ_m); f and g are the general nonlinear functions; and Θ is the parameter to be identified. Assume that the measurement noise v is Gaussian white noise with zero mean and a covariance matrix of R.

The output error method based on the maximum likelihood method is used to estimate the aerodynamic derivatives to be identified in Equations (22)–(24). The optimization function and iterative algorithm are shown in Equations (26) and (27).

$$J(\Theta ,R)=\frac{1}{2}{\displaystyle \sum _{k=1}^N{[z(t_k)-y(t_k)]}^T{R}^{-1}[z(t_k)-y(t_k)]}+\frac{N}{2}In\left|R\right|$$

(26)

Starting from the specified initial value Θ₀, the Gauss–Newton solution algorithm shown in Equation (27) is used to iterate to find the optimal parameters to be identified [29].

$$\{\begin{array}{l}{\Theta }_{i+1}={\Theta }_i+\Delta \Theta \\ \Delta \Theta =-\frac{{\displaystyle \sum _{k=1}^N{\left[\frac{\partial y(t_k)}{\partial \Theta }\right]}^T{R}^{-1}\left[z(t_k)-y(t_k)\right]}}{{\displaystyle \sum _{k=1}^N{\left[\frac{\partial y(t_k)}{\partial \Theta }\right]}^T{R}^{-1}\left[\frac{\partial y(t_k)}{\partial \Theta }\right]}}\end{array}$$

(27)

5. Verification of Identification Results Based on the Wind Tunnel Virtual Flight Test

In conclusion, the lateral-directional aerodynamic parameter identification method of the test model based on wind tunnel virtual flight is formed. The specific identification process is shown in Figure 11. First, the lateral-directional flight dynamics equations of wind tunnel virtual flight are established, and the equations of motion are linearised and decoupled to analyse the differences between the lateral-directional motion equations of wind tunnel virtual flight and free flight, so as to establish the identification model. Second, the amplitude–frequency characteristics for the equations of motion are analysed to determine the type of excitation signal, and the frequency range and amplitude of the aileron and rudder signals based on the identifiability requirements of the aerodynamic derivatives. Third, the least squares method is used to identify the side force derivatives of the wind tunnel virtual flight test, which can measure the aerodynamic force of the whole aircraft in real time through the force balance. Based on the equations of motion for wind tunnel virtual flight, the identification of rolling and yawing moment derivatives is completed by maximum likelihood estimation. Finally, the verification of the identification results is completed.

5.1. Facility for Wind Tunnel Virtual Flight Tests

The wind tunnel virtual flight test was carried out in a FL-10 wind tunnel at the AVIC Aerodynamics Research Institute [30]. As shown in Figure 12, the test model is installed in the closed test section of the wind tunnel, the size of which is 8 m × 6 m, and the maximum wind speed of the wind tunnel is 110 m/s. The test model scaling ratio k = 1/9 and the design parameters of the test model can be obtained according to the similarity criterion between the full-size aircraft and the scaled-model [31,32], as shown in

Table 3. Design parameters of the test model.

Parameters	Proportions	Full-Size Aircraft	Test Model
Wing span b (m)	1/9	36	4.00
Mean aerodynamic chord c (m)	1/9	10.41	1.16
Wing area S (m2)	(1/9)2	241	2.98
Mass m (kg)	(1/9)3	49,149	67.42
Pitch moment of inertia Iy (kg·m2)	(1/9)5	4,044,856	68.50
Yaw moment of inertia Iz (kg·m2)	(1/9)5	5,166,787	87.50
Roll moment of inertia Ix (kg·m2)	(1/9)5	1,210,504	20.50
Product of inertia Ixz (kg·m2)	(1/9)5	171,242	2.90

.

The flight test control and measurement system are shown in Figure 13. The aircraft model is equipped with attitude measurement sensors, wind vane sensors, rotary encoders, flight control computers, and servos. The IMU attitude sensor is used to measure the three-axis attitude angle and angular velocity of the model; the wind vane sensor is used to measure the angle of attack and sideslip angle of the model; the rotary encoder is used to measure the deflection angle of the rudder surface; and servos drive the deflection of the control surfaces. The flight control computer collects data from the on-board sensors and runs control law algorithms to control the attitude and altitude of the test model through the control surfaces. In addition, the ground station is developed and designed based on the LabVIEW software platform, which can display the data during the wind tunnel virtual flight test in real time.

The lateral-directional control law structure of the model is shown in Figure 14. The direction of the roll axis adopts the roll angle command configuration ϕ_cmd, and the direction of the yaw axis is controlled by the sideslip angle command β_cmd.

The speeds V = 30 m/s and α = 6° are chosen as the initial trim flight state for the wind tunnel test, and the lateral-directional control law parameters of the test model are K_p = 0.25, K_P−ϕ = 2.18, K_I−ϕ = 0.14, K_β = 0.63, K_r = 0.37, K_P−β = 0.10, K_I−β = 0.50, and K_ari = 0.33.

5.2. Design of the Excitation Signal

5.2.1. Determination of the Rudder and Aileron Frequency Range

According to the analysis in Section 3.1, the frequency sweep signal is selected as the input signal for the rudder, as well as the aileron. The lateral-directional motion of the model is shown in Equations (18) and (22) are established based on the prior aerodynamic data of the test model. The amplitude–frequency response curves of the side force, yawing moment, and rolling moment components varying with the frequency of the rudder and aileron are shown in Figure 15 and Figure 16, respectively.

The Dutch roll mode frequency of the test model at V = 30 m/s and α = 6° flight state is approximately 2.5 rad/s. The resonance peaks of the amplitude–frequency characteristic curves of the aerodynamic parameters in Figure 15 and Figure 16 are approximately equal to the Dutch roll mode frequency. The frequency band ranges of the rudder excitation signals determined by the yawing moment, rolling moment, and side force aerodynamic derivatives determined by the method presented in Section 3 are ω_δr·Cn = [1.1, 5.3] rad/s, ω_δr·Cl = [0.7, 8.7] rad/s, and ω_δr·CY = [0.8, 12.0] rad/s, respectively. Taking the intersection of ω_δr·Cn, ω_δr·Cl, and ω_δr·CY, the band range of the rudder excitation signal is obtained as ω_δr = ω_δr·Cn$\cap$ω_δr·CY$\cap$ω_δr·Cl = [1.1, 5.3] rad/s. In the same way, the band range of the aileron excitation signal can be obtained as ω_δa = ω_δa·Cn$\cap$ω_δa·CY$\cap$ω_δa·Cl = [1.2, 5.2] rad/s.

In summary, to ensure the identifiability of the lateral-directional aerodynamic derivatives of the test model, the optimal frequency range for the rudder, as well as the aileron frequency sweep excitation signals, is ω_δr = [1.1, 5.3] rad/s and ω_δa = [1.2, 5.2] rad/s, approximately between 0.5 and 2 times the Dutch roll modal frequency.

5.2.2. Determination of the Signal Amplitude

The frequency range of the aileron and rudder frequency sweep signal is selected as ω_δa = ω_δr = 1.5~5 rad/s in this experiment. To make the sideslip of the test model change within the range of ±2°, the amplitude of the aileron excitation signal is set to 6°, and that of the rudder excitation signal is set to 4°. The lateral-directional excitation signal is shown in Figure 17.

5.3. Identification Results of Aerodynamic Derivatives

5.3.1. Identification of Side Force Derivatives

The wind tunnel virtual flight test is carried out under the conditions of V = 30 m/s, α = 6°, β = 0°, and δ_e = −5.6°. The outer ring of the flight control law is the instruction of 0° sideslip angle and 0° roll angle. The excitation signal is directly input on the V-tail control surface and aileron on both sides, and the side force Y_m of the test model is measured in real time through the force balance.

The side force Y_m, sideslip angle Δβ_m, yaw rate Δr_m, roll rate Δp_m, and V-tail control surface deflection angle Δδ_r of the model measured in the test are substituted into Equation (19), and the side force derivatives of the model are identified by the least squares method, shown in Equation (21). The identification results are shown in

Table 4. Identification results of the lateral-directional aerodynamic parameters.

Parameters		Traditional Wind Tunnel Measurements	$\mathbf{Identification} \mathbf{Results} \widehat{\mathsf{\Theta }}$	Deviation (%)	Standard Deviation $\frac{\mathbf{s}\left(\widehat{\mathsf{\Theta }}\right)}{\left|\widehat{\mathsf{\Theta }}\right|}\%$
Side force derivatives	Cyβ	−0.400	−0.391	−2.25	2.60
	Cyδr	0.121	0.111	−8.26	6.55
	Cyr	0.719	0.689	−4.17	5.89
Rolling moment derivatives	Clβ	−0.083	−0.090	8.43	3.55
	Clp	−0.085	−0.080	−5.88	4.36
	Clr	0.229	0.237	3.49	1.59
	Clδa	−0.031	−0.029	−6.45	6.47
	Clδr	0.020	0.021	5.00	7.23
Yawing moment derivatives	Cnβ	0.073	0.070	−4.11	1.23
	Cnp	−0.0102	−0.010	−1.96	8.25
	Cnr	−0.066	−0.063	−4.55	4.05
	Cnδr	−0.062	−0.060	−3.23	3.55

, with the first column of the table showing the aerodynamic derivatives measured in the traditional wind tunnel test. Comparing the identification results with the conventional wind tunnel measurements, it can be seen that the identification results of side force derivatives are close to the reference values, and the deviations are within 10%.

The identification results of the side force derivatives are substituted into Equation (18) to establish the side force observation model of the wind tunnel virtual flight test. Given the excitation signals used in the wind tunnel virtual flight test, the simulation results of the side force observation model, established based on the identification results, are compared with the side force data measured by the test, as shown in Figure 18. The simulation results of the identification model are in high agreement with the data measured by the wind tunnel virtual flight test, and the maximum errors of the side force curves are within 5%.

5.3.2. Identification of Rolling and Yawing Moment Derivatives

The measured angle of attack, sideslip angle, three-axis attitude angle, and three-axis angular rate are taken into the maximum likelihood solution equations shown in Equations (25) and (26) to identify the rolling and yawing moment derivatives of the model, and the identification results are shown in Table 4.

It is to be highlighted that the results for side-force came from least squares (i.e., equation-error method), whereas the results for yawing and rolling moment were obtained via output error method (using the maximum likelihood principle).

According to Table 4, the deviation between the identification results of the yawing and rolling moment derivatives and the conventional wind tunnel measurements is within 10%. The identification results of the rolling and yawing moment derivatives are substituted into Equation (22) to obtain the lateral-directional motion model for the wind tunnel virtual flight test. Given the rudder and aileron excitation signals used in the test, the simulation results of the identified model are compared with the corresponding test data, as shown in Figure 19. The motion response of the identified model is basically consistent with the data of the wind tunnel virtual flight test, and the maximum deviation of all motion parameters is within 5%.

Goodness of fit (GOF) is used to evaluate the consistency between the lateral-directional motion response of the identified model and the wind tunnel virtual flight test data. The GOF formula is shown in Equation (28), where z_i is the experimental measurement, $\stackrel{-}{z}$ is the mean value of the observed data, and y_i is the theoretical output; a GOF value close to 1 indicates that the two sets of data fit well [33].

$$GO{F}_i=1-\frac{{(z_i-y_i)}^T(z_i-y_i)}{{(z_i-\overline{z})}^T(z_i-\overline{z})}$$

(28)

The GOF between the lateral-directional motion response of the identification model and the wind tunnel virtual flight test is shown in

Table 5. The goodness of fit between the identification model and the test data.

GOF	p	r	β	ϕ
Identification model and wind tunnel virtual flight test	0.98	0.99	0.97	0.98

. The GOF of all motion variables are greater than 0.95. It is proven that the lateral-directional motion response of the identified model is consistent with the wind tunnel virtual flight test data.

In summary, the step-by-step identification method of lateral-directional aerodynamic derivatives introduced in this paper can obtain accurate identification results, and the deviation of all the aerodynamic derivatives from the conventional wind tunnel measurements is within 10%.

5.4. Verification of the Modified Identification Model

By substituting the aerodynamic derivatives identified in Table 4 into the A₁ and B₁ matrices in Equation (12) and eliminating the A₂ and B₂ matrices generated by the influence of the support device, a modified lateral-directional motion model for wind tunnel virtual flight can be constructed, as shown in Equation (29).

$$\left[\begin{array}{c}\dot{\beta }\\ \dot{p}\\ \dot{r}\\ \dot{\varphi }\end{array}\right]=\left[\begin{array}{cccc}-0.33& 0.12& -0.96& 0.32\\ -32.60& -5.51& 1.51& 0\\ 4.22& -0.32& -0.42& 0\\ 0& 1& 0.11& 0\end{array}\right]\left[\begin{array}{c}\Delta \beta \\ \Delta p\\ \Delta r\\ \Delta \varphi \end{array}\right]+\left[\begin{array}{cc}0& 0.10\\ -19.48& 7.04\\ -0.17& -4.46\\ 0& 0\end{array}\right]\left[\begin{array}{c}{\delta }_a\\ {\delta }_r\end{array}\right]$$

(29)

V = 30 m/s and α = 6° are selected as the steady flight state of the test model. The 3211-multipole square waves are input to the V-tail control surface, while the other control surfaces are kept fixed. Comparing the simulation results of the identified model of wind tunnel virtual flight before and after the correction with the simulation data of the free flight motion model, the results are shown in Figure 20. The motion models for both wind tunnel virtual flight and free flight are linearised motion models. It should be clarified that the free flight motion model is based on the results of traditional wind tunnel test force measurements, and the correction of the wind tunnel virtual flight motion model is obtained by eliminating the influence of the test rig.

The goodness-of-fit (GOF) is shown in

Table 6. Goodness of fit of the wind tunnel virtual flight test with free flight.

GOF	p	r	β	ϕ
Flight dynamics model before correction	0.84	0.82	0.77	0.83
Flight dynamics model after correction	0.97	0.99	0.98	0.97

, comparing the responses of the identified model before and after the correction with the free flight simulation results.

As shown in Figure 20, the response of the modified identification model for wind tunnel virtual flight is in high agreement with the simulation results of the lateral-directional motion of free flight. A small phase shift between the modified wind tunnel virtual flight motion response and the free flight motion response can be seen in Figure 20, but we believe that this small difference can be neglected. Since the motion models of both wind tunnel virtual flight and free flight are linearised motion models, the inaccuracy of the linearised motion models may lead to small differences in the responses of the two. The GOF of each motion data of the identified model before the correction is less than 0.85, while the GOF of each motion data of the model after the correction is greater than 0.95.

In summary, compared with the original wind tunnel virtual flight motion model, the modified motion of the model matches better with the motion of the free flight model, indicating that the modified model can more accurately characterise the lateral-directional motion of the aircraft in free flight.

6. Conclusions

(1)

The lateral-directional flight dynamics equations of the wind tunnel virtual flight test model are established. By linearizing the equations of motion to describe the wind tunnel virtual flight test as a decoupled form of the effects of free flight aerodynamic forces and support forces on the model motion, the differences between the lateral-directional dynamics of wind tunnel virtual flight and free flight can be analysed more intuitively, thus establishing a model for the identification of aerodynamic parameters.

(2)

Based on the identifiable requirements of aerodynamic derivatives, the amplitude–frequency characteristics of the equations of motion are analysed to establish the type selection and parameter design method of the lateral-directional excitation signal. The frequency of the lateral-directional excitation signal needs to be between 0.5–2 times the frequency of the Dutch roll mode. Therefore, a suitable actuator needs to be selected for different flight states, so that the deflection rate of the control surface is fast enough. In addition, to identify the aerodynamic model at high angle of attack or high sideslip, it is necessary to design appropriate flight control law to ensure the stability of the test model.

(3)

In this paper, a step-by-step method for the identification of side force, rolling, and yawing moment derivatives is established. The identification of the side force derivatives can be achieved by measuring the aerodynamic force of the test model with a force balance. The differences between the identification results of the aerodynamic derivatives and the conventional wind tunnel measurements are within 10%. The lateral-directional motion response of the identified model is basically consistent with the wind tunnel virtual flight test data, and the GOF of all motion variables are greater than 0.95.

(4)

The modified identification model can well-simulate the lateral-directional motion of the conventional wind tunnel test model, and the goodness-of-fit is higher than 0.95, which verifies the correctness of the proposed method.

The proposed method of identifying the lateral-directional aerodynamic parameters of aircraft based on wind tunnel virtual flight tests is simple in principle, and the identification results are highly accurate. The wind tunnel virtual flight test enables the accurate identification of the lateral-directional aerodynamic model at an early stage of aircraft design, which can be compared and corrected with the force measurement results of conventional wind tunnel tests. The proposed method has outstanding advantages in shortening the development cycle of aircraft and reducing risks and costs, but the application of the proposed method to the parameter identification of rotorcraft or spacecraft needs further study.