Nomenclature
-
$\alpha $
-
AoA (angle-of-attack)
-
$\beta $
-
sideslip angle
-
$\phi, \theta, \psi $
-
euler angles (roll, pitch, yaw)
- p, q, r
-
inertial angular rates about body’s axes (roll, pitch, yaw)
- l, m, n
-
aerodynamic moments (roll, pitch, yaw)
- L, Y, D
-
aerodynamic forces (lift, side-force, drag)
-
${\delta _e},{\delta _a}$
-
control deflections (elevator, aileron)
-
$\rho $
-
air density
- V
-
true air speed
-
${C_L},{C_Y},{C_D},{C_l},{C_m},{C_n}$
-
aerodynamic coefficients
1.0 Introduction
Active gust alleviation control technologies, designed to enhance ride comfort and mitigate hazardous vertical fluctuations, are currently being integrated into the implementation level of autopilots for large passenger aircraft during cruise phases [Reference Regan and Jutte1]. The typical control mechanism involves both alleviating gust loads and damping structural wing oscillations by utilising control surfaces that receive feedback from accelerometers embedded in both the fuselage and the wingtip.
However, the development of highly reliable stabilisation control for terrain-induced turbulence encountered by airplanes at low altitudes is still ongoing. This issue also applies to small aircraft operating in the near-ground layer, including general aviation and the recently developed urban air mobility [Reference Abdulrahim, Watkins, Segal, Marino and Sheridan2, Reference Watkins, Milbank, Loxton and Melbourne3]. Due to their increased sensitivity to turbulence [Reference Watkins, Thompson, Loxton and Abdulrahim4], small aircraft require robust automatic stabilisation control to effectively handle severe turbulence, which often lacks a mathematical model, such as strong winds passing through structures.
For wind gusts in free space at high altitudes, their typical spectral characteristics are well-known, as demonstrated by the Dryden model. This knowledge allows for the optimisation of automatic control systems and evaluation through simulations using artificially generated time series of wind disturbances [Reference Stevens, Lewis and Johnson5].
On the other hand, wind gusts occurring at low altitudes are challenging to model in advance due to their unsteady nature and variation depending on the terrain. The greater the uncertainty in the disturbance characteristics, the more conservative the control performance is achieved through fixed gain control.
To attain higher control performance, research has been conducted on actively utilising gust information. Fluctuations in the angle-of-attack (AoA) during a cruise can be considered as wind disturbances that induce unexpected aircraft motion. Consequently, methodologies that involve providing feedback from AoA sensors to the control system to counteract the disturbance effect just before pronounced oscillation has been studied for several decades [Reference Krag, Rohlf and Wuennenberg6, Reference Koenig and Hahn7].
However, a dilemma arises due to the response delay of the AoA sensor in turbulent air and the actuator delay, both of which decrease closed-loop stability margins when the control gain is increased to achieve higher performance. As a result, the effectiveness of directly counteracting the effects of wind disturbances with control surfaces is limited, and such a method has not yet been implemented in practice.
To address this problem, the authors have developed the flow-field-integrated flight control (FIFC) system, which can be easily installed in small aircraft as well. This method utilises in-flight pressure-field observations to actively control surface deflection, employing multiple static pressure taps distributed primarily on the wing surface of the aircraft.
Typically, pressure measurements on the aircraft surface are obtained through wind tunnel tests and computational fluid dynamics (CFD) analysis to understand flow-field properties, including pressure distribution and separation progression. In contrast to these measurements, which are focused on understanding aerodynamics, the proposed method employs pressure sensors to acquire more refined information about the atmospheric environment during flight than AoA sensors.
An auto-regressive moving average model of angular velocity, driven by the pressure field, is introduced, and the coefficient matrix in the model is updated using a Kalman filter. This equation automatically incorporates changes in aerodynamic and inertia properties during flight.
Flight control utilising in-flight pressure-field sensing has also been studied by several university research groups. Most of the research is experimental, utilising model airplanes and drawing inspiration from biological flight with turbulence-resistant performance.
One important issue is the optimal layout of discrete pressure observations, but currently, it remains primarily an experimental consideration. Researchers at RMIT University have published the results of experimental research, which include free flights conducted in a large wind tunnel facility. These studies have revealed the correlation between turbulence just before the wing and wing surface pressure fluctuations, offering insights into attitude control [Reference Mohamed, Watkins, Fisher, Marino, Massey and Clothier8, Reference Marino, Ravi and Watkins9, Reference Mohamed, Watkins, Clothier and Abdulrahim10].
Researchers at the University of Bristol present a method for constructing a control system using in-flight pressure-field sensing, with the expectation that AoA estimation and control would prevent flight instability during turbulence encounters. The proposed control system takes the form of a neural network [Reference Guerra-Laugan, Araujo-Estrada, Richards and Windsor11].
Additionally, researchers at the University of Central Florida demonstrate how to geometrically approximate the actual aerodynamic forces and moments acting on a wing using discrete pressure data. To maintain the accuracy of aerodynamic force estimation, it is necessary to position as many pressure taps as possible within a specified area. They further determine control commands based on the resulting equations of the angular velocity of motion [Reference Shen, Xu and Remeikas12, Reference Shen and Xu13, Reference Shen, Xu and Dickinson14, Reference Thompson, Xu and Dickinson15].
Other related studies have explored the incorporation of wind disturbance look-ahead information as a feed-forward loop [Reference Watkins, Mohamed, Abdulrahim, Marino, Clothier, Fisher, Wild and Ravi16, Reference Watkins, Abdulghani, Prudden, Tennent, Marino and Fisher17, Reference Giesseler, Kopf, Varutti, Faulwasser and Findeisen18, Reference Hamada, Kikuchi and Inokuchi19]. Doppler lidar and five-hole pitot tubes are utilised as sensors for detecting disturbances.
This paper presents a method for real-time updating of the equation of motion for angular rate using discrete pressure sensing and a Kalman filter type parameter estimation algorithm [Reference Ljung and Soderstrom20] (Chapter 4). Instead of directly estimating aerodynamic moments, the surface pressure sensing and the angular velocities from an onboard gyro sensor are directly linked through linear regression. This approach eliminates the need to estimate aerodynamic and inertia properties for constructing a mathematical model of the aircraft. The control law design is continuously updated in real time, providing an adaptive mechanism to account for uncertainties in aerodynamics and inertia properties.
The effectiveness of the FIFC is demonstrated through wind tunnel experiments and numerical simulations. In contrast to typical flight simulations that utilise predefined aerodynamic models, the computation of variations in the surface pressure field is necessary. A panel method for steady flow fields is employed in this study. Additionally, the accuracy of the simulation is verified through newly designed dynamic wind tunnel experiments.
Traditionally, the design and fabrication of the dynamic support mechanism and gust generator in a wind tunnel have been costly and involved complex mechanisms [Reference Owens, McConnell, Brandon and Hall21, Reference Owens, Brandon, Croom, Fremaux, Heim and Vicroy22, Reference Vicroy23, Reference Pattinson, Lowenberg and Goman24]. However, in this study, we devised a simple dynamic support system that can be constructed by students on a small budget in a university laboratory. We provide detailed information regarding the hardware and mathematical model of this dynamic support system.
The experimental airplane mounted on the support system is essentially a commercial model aircraft equipped with pressure sensors embedded in the wings. This airplane serves for both wind tunnel experiments and flight experiments.
The layout design of pressure taps on the aircraft surface is a crucial consideration as it directly impacts the ability to represent motion variables such as angular velocity and acceleration, as well as the resulting control performance. This paper presents an attempt to prioritise important pressure taps using LASSO regression [Reference Tibshirani30] as an initial step in discussing pressure observation network design issues. However, since validation through experimentation has not yet been conducted, it is discussed in the Appendix. The evaluation index used is the capability to approximate the time variability of angular velocity and acceleration.
If coupled with suitable guidance laws, the onboard program developed for wind tunnel experiments can be utilised directly for flight with minimal modifications. However, automatic flight experiments will be reported at a future date.
2.0 Experimental setup
2.1 Model aircraft
2.1.1 Airframe and sensors
An Opterra™ 2m flying wing manufactured by Horizon Hobby shown in Fig. 1 was used as the flying platform in the experiments and simulations.
The experimental airplane used in this study is a commercially available, tailless, remote-controlled model aircraft. This particular aircraft possesses minimal static stability, making it suitable for evaluating the effectiveness of active control performance against turbulence. Tailless airplanes also offer advantages in numerical analysis as the wing wake panels do not interfere with other components regardless of the aircraft’s attitude.
The aircraft specifications are presented in Table 1. It features a single control surface, called an Elevon, on each side, which works as both the elevator and aileron. The aircraft is lightweight and can be launched by hand.
The aircraft parameters
Opterra™ 2m flying wing.
As described in detail in the following section, static pressure taps are positioned on the surface, and MEMS pressure transducers are integrated inside the wing. A pitot tube is mounted in the nose section of the airframe to measure dynamic pressure and atmospheric pressure during flight. Flight motion data is captured using the XSENS MTi-670G Attitude Heading Reference System (AHRS). Conversely, in the enclosed environment of the wind tunnel, a magnetic field-enhanced AHRS, the XSENS MTi-630R, which does not rely on GPS information, was employed. For detailed specifications of the AHRS, see the XSENS website at https://www.movella.com/products/sensor-modules/xsens-mti-3-ahrs.
Define the coordinate system fixed to the fuselage as depicted in Fig. 2. In this coordinate system, the elevator is defined as positive for pitching up, and the aileron is defined as positive for rolling to the right.
Body-fixed coordinate system.
2.1.2 Pressure measurement system
Figure 3 illustrates the process of checking the feasibility of embedding a wired-connected pressure transducer for a sliced wing. Due to the limited space within the wing for wiring, except for areas occupied by the main spar and servo motors, the layout prioritised maximizing the number of measurement points. Figure 4 is a schematic diagram illustrating the placement of the multiple pressure sensors and the internal wiring within the wing. A 1-mm diameter brass pipe is flush-mounted with the wing surface and connected to a MEMS pressure transducer inside the wing through a urethane tube. The length of the tube is adjusted to be as short as possible to minimise measurement delay.
Tearing open the wing to examine its internal structure and assess potential space for electronic component wiring. Due to its weaker structural strength, the half-cut wing in the photo was not used in the experiment.
Total 32 pressure taps are symmetrically positioned on both the left and right wing surfaces. The left side displays the layout of the pressure taps, while the right side depicts the wiring inside the wing. Red and blue indicate the upper and lower pressure taps, respectively. Each pressure tap is numbered so that the number of symmetrically located taps adds up to 33. Differential pressure, obtained by subtracting the static pressure from the pitot tube, is measured between the wing surface pressure and the static pressure. The measured pressure values are collected on a sub-board located inside the wing and transmitted to the FCC (flight control computer) via signal lines. The pressure transducer and the pressure taps are connected using a urethane tube. The length of the tube is adjusted to minimise measurement delay, making it as short as possible.
The pressure transducer used is DLHR-L05D-E1BD-C-NAV8 manufactured by Amphenol All Sensors. The sampling rate of the sensor is 100 Hz. The final pressure taps arrangement is shown on the left side in Fig. 4, where red and blue indicate the measurement points on the upper and lower surfaces, respectively. The pressure taps are numbered so that they correspond to the symmetrically located pressure taps, and are defined so that the number of symmetrically located taps adds up to 33 (see Fig. 4). For example, the pressure taps symmetrically located at P5 is P28. Electromagnetic interference from the servo motor wiring was found and removed by installing ferrite cores in appropriate locations.
2.2 Dynamic support system
2.2.1 Motivation
Although the FIFC system uses aircraft surface pressure in real-time, it is difficult to virtually generate realistic pressure fluctuations in simulation, especially for an oscillating target in a gusty wind environment. Therefore a test environment that allows multi-degree-of-freedom motion in a wind tunnel is necessary to evaluate control performance. On the other hand, an outdoor flight test would provide a more accurate discussion based on the actual relationship between wing surface pressure fluctuations and aircraft motion. However, it is not appropriate to evaluate control performance by outdoor flight. This is because the atmospheric environment during actual flight is unknown, and it is not possible to create the same atmospheric environment again. For these reasons, wind tunnel experiments will be conducted using an easily craftable dynamic support system and a primitive but useful gust generator. Refer to Section 5.1.2 for the mechanism of gust generation.
2.2.2 Configuration of 5-DoF support system
This system is basically a rotatable seesaw. The configuration of the dynamic support system is shown in Fig. 5, and its installation in a wind tunnel cart is shown in Fig. 6. The dynamic support system comprises a steel base, a steel vertical Link 1, and a lightweight CFRP horizontal Link 2. The airplane is supported from below at the centre of gravity using a ping-pong ball joint. The airplane has three degrees of rotational freedom on the ball joint. The base of the support system consists of heavy steel plates and bearings to enable smooth yaw motion.
Configuration of dynamic wind tunnel testing model.
The picture of dynamic wind tunnel testing model. The safety fences are temporary installations set up in preparation, so they are removed during actual airflow testing.
The rotational motion of Link 1 around the vertical axis is called ‘link-yaw’, and the angle is denoted by
$\eta $
. The rotational motion of Link 2 around the fulcrum is called ‘link-pitch’ and the angle is denoted by
$\xi $
. See Fig. 7 for the definition of
$\xi $
and
$\eta $
.
Simplified link model. (Left) Nomenclature. (Right) Fixed coordinates of each link.
A counterweight is attached to the rear end of Link 2. This is done to ensure a balance between the lift force and the weight of the aircraft at the reference wind speed of 8m/s for the wind tunnel experiment. The actual cruising speed of the aircraft is approximately 20m/s. However, if a control failure occurs at the higher wind speed of 20m/s, the link angle rapidly reaches its limit, potentially causing damage to the support system. Therefore, the wind speed was set to 8m/s for this experiment. Consequently, a counterweight is placed at the end of Link 2 to achieve equilibrium between the lift force and the weight of the aircraft at a wind speed of 8m/s.
A simple gust generator using a polyethylene sheet was set upstream of the airplane (see Fig. 15 in Section 5.1.2). The polyethylene sheet was attached to a pole and moved manually up and down. The small flutter of the polyethylene sheet can generate a disturbance of about ±5 degrees in terms of angle-of-attack variation (see Fig. 16). In addition, the vertical movement of the pole can simulate abrupt changes in dynamic pressure.
2.2.3 The mathematical model of the dynamic support system
This chapter presents the mathematical model of the dynamic support system. As the support system has five degrees of freedom, directly deriving the equations of motion would be complex. To deal with this problem, the three degrees of freedom related to the aircraft’s attitude and the two degrees of freedom associated with the support system are considered separately.
Firstly, the three degrees of freedom of the aircraft’s attitude are represented by the following equations.
\begin{align}l = {I_{xx}}\dot p - {I_{xz}}\dot r + qr\!\left( {{I_{zz}} - {I_{yy}}} \right) - {I_{xz}}pq\end{align}
\begin{align}m = {I_{yy}}\dot q + rq\!\left( {{I_{xx}} - {I_{zz}}} \right) + {I_{xz}}\!\left( {{p^2} - {r^2}} \right)\end{align}
\begin{align}n = - {I_{xz}}\dot p + {I_{zz}}\dot r + pq\!\left( {{I_{yy}} - {I_{xx}}} \right) + {I_{xz}}qr\end{align}
Where
${I_{xx}}$
,
${I_{yy}}$
,
${I_{zz}}$
and
${I_{xz}}$
are the moment of inertia of the aircraft with respect to the body-fixed coordinate system.
Next, the equations of motion for the two degrees of freedom of the support system are derived using the Lagrangian equations. In this case, the aircraft model is treated as a mass point that generates the aerodynamic forces. Figure 7 shows the simplified link model and the link-fixed coordinate systems. In these link-fixed coordinate systems, the moment of inertia of each link is defined. The subscript 0 denotes the wind tunnel fixed coordinate system. The subscript 1 denotes the Link 1 fixed coordinate system, and the subscript 2 denotes the Link 2 fixed coordinate system.
If the aircraft, counterweight and link 2 are regarded as a single rigid body (hereinafter referred to as Link 2’), the moment of inertia of Link 2’ can be expressed as follows:
\begin{align}{I^{\prime}_{x{x_2}}} = {I_{x{x_2}}}\end{align}
\begin{align}{I^{\prime}_{y{y_2}}} = {I_{y{y_2}}} + {m_2}{d^2} + {m_a}l_a^2 + {m_c}l_c^2\end{align}
\begin{align}{I^{\prime}_{z{z_2}}} = {I_{z{z_2}}} + {m_2}{d^2} + {m_a}l_a^2 + {m_c}l_c^2\end{align}
The moments of inertia of Link 2 are denoted as
${I_{x{x_2}}}$
,
${I_{yy2}}$
, and
$Iz{z_2}$
. Note that since the aircraft and counterweight are considered as mass points, they do not contribute to the moment of inertia around the
${X_2}$
axis.
The kinetic energy of Link 1 is represented as
${T_1}$
, while the kinetic energy of Link 2 is denoted as
${T_{2{\rm{^{\prime}}}}}$
. These are expressed as follows:
\begin{align}{T_1} = \frac{1}{2}{I_{z{z_1}}}{\dot \eta ^2}\end{align}
\begin{align}T_2^{\rm{^{\prime}}} = \frac{1}{2}\left( {{{I^{\prime}}_{z{z_2}}}{{\dot \eta }^2}{\rm{co}}{{\rm{s}}^2}\xi + {{I^{\prime}}_{x{x_2}}}{{\dot \eta }^2}{\rm{si}}{{\rm{n}}^2}\xi + {{I^{\prime}}_{y{y_2}}}{{\dot \xi }^2}} \right)\end{align}
Next, the potential energy of this system will be derived. The potential energy of Link 1, Link 2, the aircraft, and the counterweight with respect to the wind tunnel ground can be expressed as follows:
\begin{align}{U_1} = {m_1}g\frac{{{l_1}}}{2}\end{align}
\begin{align}{U_2} = {m_2}g\!\left( {{l_1} + d\,{\rm{sin}}\,\xi } \right)\end{align}
\begin{align}{U_a} = {m_a}g\!\left( {{l_1} + {l_a}\,{\rm{sin}}\,\xi } \right)\end{align}
\begin{align}{U_c} = {m_c}g\!\left( {{l_1} - {l_c}\,{\rm{sin}}\,\xi } \right)\end{align}
the Lagrangian
${\mathcal{L}}$
is defined as follows:
\begin{align}{\mathcal{L}} = \left( {{T_1} + T_2^{\rm{^{\prime}}}} \right) - \left( {{U_1} + {U_2} + {U_a} + {U_c}} \right)\end{align}
In general, the Lagrange’s equation is expressed as follows:
\begin{align}\frac{d}{{dt}}\left( {\frac{{\partial {\mathcal{L}}}}{{\partial {{\dot q}_i}}}} \right) - \frac{{\partial {\mathcal{L}}}}{{\partial {q_i}}} + \frac{{\partial \bar D}}{{\partial {{\dot q}_i}}} = {Q_i}\!\left( {i = 1, \ldots, n} \right)\end{align}
Where
${q_i}$
and
${Q_i}$
are the generalised coordinates and the generalised forces respectively. In this case, generalised coordinates are
$\eta $
and
$\xi $
.
$\bar D$
is called the Rayleigh dissipation function which is to account for the frictional forces. To simplify the analysis, we assume that there is no friction acting on the support system, leading to
$\bar D = 0$
.
Let us denote the moments acting on Link 1 and Link 2 as
${M_\eta }$
and
${M_\xi }$
, respectively. By applying Equation (14), the following equations of motion can be derived:
\begin{align}{M_\eta } = \left( {{I_{z{z_1}}} + {{I^{\prime}}_{x{x_2}}}{\rm{si}}{{\rm{n}}^2}\,\xi + {{I^{\prime}}_{z{z_2}}}{\rm{co}}{{\rm{s}}^2}\xi } \right)\ddot \eta + \left( {{{I^{\prime}}_{x{x_2}}} - {{I^{\prime}}_{z{z_2}}}} \right)\dot \eta \dot \xi \,{\rm{sin}}\,2\xi \end{align}
\begin{align}{M_\xi } = {I^{\prime}_{y{y_2}}}\ddot \xi + \left( {{{I^{\prime}}_{z{z_2}}} - {{I^{\prime}}_{x{x_2}}}} \right){\dot \eta ^2}\,{\rm{sin}}\,\xi \,{\rm{cos}}\,\xi + \left( {{m_2}d + {m_a}{l_a} - {m_c}{l_c}} \right)g\,{\rm{cos}}\,\xi \end{align}
The moments exerted on Link 1 and Link 2 can be derived from the aerodynamic forces acting on the fuselage. Firstly, the aerodynamic forces acting on the aircraft, expressed in the body-fixed coordinate system, are transformed into the wind tunnel fixed coordinate system. This coordinate transformation can be achieved by multiplying the aerodynamic force vector by the following coordinate transformation matrix:
\begin{align}{^0}{{\textbf{F}}_a} = {}^{0}{{\textbf{T}}_B}^B{{\textbf{F}}_a}\end{align}
\begin{align}{{\textbf{F}}_a} = {{[^B}{X_a}\quad {}^B{Y_a}\quad {}^B{Z_a}]^{\rm{{\rm T}}}}\end{align}
\begin{align}{{\textbf{T}}_B} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}
{{\rm{cos}}\,\psi\, {\rm{cos}}\,\theta } & {{\rm{cos}}\,\psi\, {\rm{sin}}\,\phi\, {\rm{sin}}\,\theta - {\rm{cos}}\,\phi\, {\rm{sin}}\,\psi } & {{\rm{cos}}\,\phi\, {\rm{cos}}\,\psi\, {\rm{sin}}\,\theta + \,{\rm{sin}}\,\phi\, {\rm{sin}}\,\psi }\\[4pt]
{{\rm{cos}}\,\theta\, {\rm{sin}}\,\psi } & {\,{\rm{sin}}\,\phi\, {\rm{sin}}\,\psi\, {\rm{sin}}\,\theta + {\rm{cos}}\,\phi\, {\rm{cos}}\,\psi } & {{\rm{cos}}\,\phi\, {\rm{sin}}\,\psi\, {\rm{sin}}\,\theta - {\rm{cos}}\,\psi\, {\rm{sin}}\,\phi }\\[4pt]
{ - \,{\rm{sin}}\,\theta } & {{\rm{cos}}\,\theta\, {\rm{sin}}\,\phi } & {{\rm{cos}}\,\phi\, {\rm{cos}}\,\theta }\end{array}} \right]\end{align}
Where
${{\textbf{F}}_a}$
is the aerodynamic force vector,
${X_a}$
,
${Y_a}$
and
${Z_a}$
are the aerodynamic forces acting on the aircraft.
The subscript in the upper left corner denotes the coordinate system in which the vector is presented. The subscript
$B$
represents the body-fixed coordinate system (see Fig. 2), and the subscript
$0$
represents the wind tunnel fixed coordinate system (see Fig. 7). The
${^0}{{\textbf{T}}_B}$
is the transformation matrix that transforms the vector from the body-fixed coordinate system to the wind tunnel fixed coordinate system. Once the aerodynamic forces are expressed in the wind tunnel fixed coordinate system, it is easy to calculate the moments acting on each link.
The moments
${M_\eta }$
and
${M_\xi }$
are expressed as follows:
\begin{align}{M_\eta } = \left( {{ - ^0}{X_a}\,{\rm{sin}}\,\eta + {}^{0}{Y_a}\,{\rm{cos}}\,\eta } \right){l_a}\,{\rm{cos}}\,\xi \end{align}
\begin{align}{M_\xi } = { - ^0}{Z_a}{l_a}\,{\rm{cos}}\,\xi - {(^0}{X_a}\,{\rm{cos}}\,\eta + {}^{0}{Y_a}\,{\rm{sin}}\,\eta ){l_a}\,{\rm{sin}}\,\xi \end{align}
The previous equations (Equations (1), (2), (3), (15) and (16)) describe the equations of motion of the dynamic support system.
3.0 Flight simulation environment
3.1 Overview
This chapter describes the flight simulation configuration. We have developed a simulator that incorporates flow field calculations using a 3D panel method for steady flow. However, if the attitude and link angles change rapidly, involving unsteady flow, the control capability will not be sufficient, and the model support system’s declination will quickly saturate. Therefore, for simulating wind tunnel experiments, the flow field calculations are limited to a steady flow, and the steady panel method code is utilised. Two types of motion calculation blocks are available for interchange: one for dynamic wind tunnel experiments that include the support system, and the other for free flight simulations.
3.2 Mathematical model
Figure 8 shows the block diagram of the flight simulation environment. The panel method block receives inputs such as the angle-of-attack, angle of sideslip, elevator angle, and aileron angle and produces outputs including the pressure coefficient on the aircraft panel and the steady aerodynamic coefficients for forces and moments. The pressure coefficient is expressed by the following equation:
\begin{align}C_P^{\left( j \right)} = \frac{{{P^{\left( j \right)}} - {P_s}}}{{\frac{1}{2}\rho {V^2}}} = \frac{{{P^{\left( j \right)}} - {P_s}}}{{\bar q}}\end{align}
where
$\bar q\!\left( { = 1/2\rho {V^2}} \right)$
is the dynamic pressure and
${P_s}$
is the static pressure of the pitot tube. The panel method used in this simulation calculates the steady part of aerodynamic forces and moments, where the contribution of unsteady aerodynamic damping effects is not included. As a result, the aerodynamic coefficients are obtained as the sum of the steady aerodynamics from panel method and the damping effect.
\begin{align}{C_L} = {C_L}\!\left( {\alpha, \beta, {\delta _e},{\delta _a}} \right) + {C_{{L_q}}}\frac{{\bar cq}}{{2V}}\end{align}
\begin{align}{C_Y} = {C_Y}\!\left( {\alpha, \beta, {\delta _e},{\delta _a}} \right) + {C_{{Y_p}}}\frac{{bp}}{{2V}} + {C_{{Y_r}}}\frac{{br}}{{2V}}\end{align}
\begin{align}{C_l} = {C_l}\!\left( {\alpha, \beta, {\delta _e},{\delta _a}} \right) + {C_{{l_p}}}\frac{{bp}}{{2V}} + {C_{{l_r}}}\frac{{br}}{{2V}}\end{align}
\begin{align}{C_m} = {C_m}\!\left( {\alpha, \beta, {\delta _e},{\delta _a}} \right) + {C_{{m_q}}}\frac{{\bar cq}}{{2V}}\end{align}
\begin{align}{C_n} = {C_n}\!\left( {\alpha, \beta, {\delta _e},{\delta _a}} \right) + {C_{{n_p}}}\frac{{bp}}{{2V}} + {C_{{n_r}}}\frac{{br}}{{2V}}\end{align}
Simulation block diagram: The aerodynamic forces and moments are computed using a 3D panel method. A subset of the overall pressure field is utilised for feedback observation in the adaptive controller. The user has the option to select the EOM for either free flight or dynamic wind tunnel testing.
Among the dynamic stability derivatives, especially
${C_{{m_q}}}$
and
${C_{{n_r}}}$
are mainly due to the presence of a tail in usual aircraft configurations. Therefore, estimating these coefficients for tailless aircraft is difficult and involves relatively large uncertainties. In this simulation, the damping coefficients from the “Skywalker X8” which is a flying wing aircraft with a similar geometry [Reference Gryte, Hann, Alam, Roháč, Johansen and Fossen25] are used. However, since the coefficients are rough estimates, it would be desirable to assume a larger error bar and assess its impact in flight simulations. However, we will not address this issue any further, as it is outside the scope of the main purpose of this paper.
As for the drag coefficient, the simplest quadratic drag modeling is used as follows:
\begin{align}{C_D} = {C_{{D_{{\rm{min}}}}}} + \frac{{{C_L}^2}}{{\pi \cdot AR \cdot e}}\end{align}
The first term,
${C_{{D_{{\rm{min}}}}}}$
, represents the minimum drag coefficient mainly due to skin friction, which cannot be accounted for by the panel method targeting inviscid flow. On the other hand, the second term represents the lift-induced drag, which can be computed by the panel method. Here, the well-known quadratic drag formula is used. In Equation (28),
$e$
is Oswald efficiency and
$AR$
is the aspect ratio. In this study, it is set to
$e = 0.8$
. The Oswald span efficiency was calculated using Raymer’s formula [Reference Raymer26], based on
$AR$
and swept angle, and was roughly estimated to be 0.8, considering the effect of the wingtip. Since the Reynolds number differs between the dynamic wind tunnel test (8m/s) and outdoor flight (20m/s), the
${C_{{D_{{\rm{min}}}}}}$
of the former simulation should be set approximately 10% higher than that of the latter.
The six-component aerodynamic forces and moments are expressed by the following equation:
\begin{align}L = \bar qS{C_L}\end{align}
\begin{align}{Y_a} = \bar qS{C_Y}\end{align}
\begin{align}D = \bar qS{C_D}\end{align}
\begin{align}l = \bar qSb{C_l}\end{align}
\begin{align}m = \bar qS\bar c{C_m}\end{align}
\begin{align}n = \bar qSb{C_n}\end{align}
These terms are used to calculate the EOM (equation of motion).
\begin{align}m\!\left( {\dot u + qw - rv} \right) = {X_a} - mg\,\,{\rm{sin}}\,\theta \end{align}
\begin{align}m\!\left( {\dot v + ru - pw} \right) = {Y_a} - mg\,{\rm{cos}}\,\theta \,{\rm{sin}}\,\phi \end{align}
\begin{align}m\!\left( {\dot w + pv - qu} \right) = {Z_a} - mg\,{\rm{cos}}\,\theta \,{\rm{cos}}\,\phi \end{align}
\begin{align}{I_{xx}}\dot p - {I_{xz}}\dot r + qr\!\left( {{I_{zz}} - {I_{yy}}} \right) - {I_{xz}}pq = l\end{align}
\begin{align}{I_{yy}}\dot q + rq\!\left( {{I_{xx}} - {I_{zz}}} \right) + {I_{xz}}\!\left( {{p^2} - {r^2}} \right) = m\end{align}
\begin{align} - {I_{xz}}\dot p + {I_{zz}}\dot r + pq\!\left( {{I_{yy}} - {I_{xx}}} \right) + {I_{xz}}qr = n\end{align}
The subscript
$B$
which represents the body-fixed coordinate system is omitted. Since the block of panel method outputs lift and drag forces, they are converted to aerodynamic forces aligned with the body-fixed coordinate system using the following equation:
\begin{align}{X_a} = L \cdot{\rm{sin}}\,\alpha - D \cdot {\rm{cos}}\,\alpha \end{align}
\begin{align}{Z_a} = - L \cdot {\rm{cos}}\,\alpha - D \cdot{\rm{sin}}\,\alpha \end{align}
The angle-of-attack and the sideslip angle in no-wind conditions are expressed using the centre-of-gravity velocity of the aircraft,
$u,v,w$
, as follows:
\begin{align}\alpha = {\rm{ta}}{{\rm{n}}^{ - 1}}\frac{w}{u}\end{align}
\begin{align}\beta = {\rm{si}}{{\rm{n}}^{ - 1}}\frac{v}{V}\end{align}
\begin{align}V = \sqrt {{u^2} + {v^2} + {w^2}} \end{align}
When considering atmospheric disturbances, each component of the disturbance can be denoted as
${u_g},{v_g},{w_g}$
. The angle-of-attack
${\alpha _{{\rm{gust}}}}$
and the sideslip angle
${\beta _{{\rm{gust}}}}$
are obtained using the following equation:
\begin{align}{\alpha _{{\rm{gust}}}} = {\rm{ta}}{{\rm{n}}^{ - 1}}\frac{{w + {w_g}}}{{u + {u_g}}}\end{align}
\begin{align}{\beta _{{\rm{gust}}}} = {\rm{si}}{{\rm{n}}^{ - 1}}\frac{{v + {v_g}}}{V}\end{align}
The sensor block receives aircraft states as input and introduces white noise to the output. The sensors in the flight simulation are AHRS, the pitot tube, and wing surface pressure sensors. However, only the gyro sensor and wing surface pressure sensors are used for FIFC.
The actuator for control surfaces is represented by a first-order delay transfer function.
\begin{align}{\delta ^{{\rm{act}}}} = \frac{1}{{{T_a}s + 1}}{\delta ^{{\rm{cmd}}}}\end{align}
where
${\delta ^{{\rm{cmd}}}}$
is the control surface angle command,
${\delta ^{{\rm{act}}}}$
is the control surface angle moved by the actuator, and
${T_a}$
is the time constant of the actuator. In this simulation, it is set to
${T_a} = 0.1$
.
${\textbf{Remark}}\;\;:$
We can switch between two types of simulations; dynamic wind tunnel testing or free flight. The EOM of the dynamic wind tunnel testing is obtained in Section 2.2, where the aerodynamic forces
${X_a},{Y_a},{Z_a}$
are given by (30), (41) and (42). On the other hand, the EOM of free flight is given by the above Equations (35)–(42).
3.3 3D Geometry of model aircraft
The airframe geometry was obtained by scanning the actual model aircraft using a 3D surface profiling camera. The scanned data was then divided into rectangular panels. As the geometry data comprises a highly detailed mesh, the number of mesh grid points was reduced, and the smoothing process was applied so as to approximate the original curved surface. The resulting panel model, depicted in Fig. 9, consists of a total of 3,994 panels.
The Opterra panel model with aileron deflection of 10 degrees.
The left and right control surfaces can operate within a range of ±20 degrees, as illustrated in Fig. 10. The position of both control surfaces is determined by mixing of the elevator and the aileron commands, as follows:
\begin{align}{\delta _R} = {\delta _e} + {\delta _a}{\rm{\;\;\;\;s}}.{\rm{t}}. - 20\!\left[ {{\rm{deg}}} \right] \le {\delta _R} \le 20\!\left[ {{\rm{deg}}} \right]\end{align}
\begin{align}{\delta _L} = {\delta _e} - {\delta _a}{\rm{\;\;\;\;s}}.{\rm{t}}. - 20\!\left[ {{\rm{deg}}} \right] \le {\delta _L} \le 20\!\left[ {{\rm{deg}}} \right]\end{align}
where
${\delta _R}$
and
${\delta _L}$
denote the right elevon angle and the left elevon angle, respectively. The elevon angle is defined as positive when the trailing edge is lifted upwards.
The control surface deflection. (a)
$ - 20$
deg, (b)
$0$
deg, (c)
$ + 20$
deg.
3.4 Aerodynamic characteristics using 3D panel method
Aerodynamic characteristics were calculated using the panel method. Figure 11 shows the curves of the aerodynamic coefficients obtained from the angle-of-attack and the sideslip angle sweeps. Both the angle-of-attack and the sideslip angle varied within the range of
$ \pm 8$
degrees. These limits were chosen due to the limitations of the panel method, which does not account for phenomena such as separation and stall. As a result, the panel method may not provide accurate results when the angle of airflow incidence is large. Figure 12 shows the curves of the aerodynamic coefficients obtained from the elevator and aileron sweeps at
$\alpha = \beta = 0$
. Figure 13 shows the pressure field at 3 degrees of the angle-of-attack.
The curves of the aerodynamic coefficients obtained from the angle-of-attack and the sideslip angle sweeps.
The curves of the aerodynamic coefficients obtained from the elevator and the aileron sweep (
$\alpha = \beta = 0$
).
The surface pressure field at 3 degrees of the angle-of-attack calculated using the 3D panel method.
4.0 Flow-field integrated flight control
4.1 Control law derivation
This chapter presents a mathematical formulation for designing angular rate control laws that utilise real-time, in-flight pressure field sensing. As the surface pressure distribution generates an aerodynamic external force that induces rotational motion, it is reasonable to consider the ARMA model in the following form:
\begin{align}\left[ {\begin{array}{c}{{q_k}}\\[4pt]
{{p_k}}\\[4pt]
{{r_k}}\end{array}} \right] = {{\textbf{F}}_k}\!\left[ {\begin{array}{c}{{q_{k - 1}}}\\[4pt]
{{p_{k - 1}}}\\[4pt]
{{r_{k - 1}}}\end{array}} \right] + {{\textbf{G}}_k}\!\left[ {\begin{array}{c}{{q_{k - 2}}}\\[4pt]
{{p_{k - 2}}}\\[4pt]
{{r_{k - 2}}}\end{array}} \right] + {{\textbf{H}}_k}\!\left[ {\begin{array}{c}{P_k^{\left( 1 \right)}}\\[4pt]
{P_k^{\left( 2 \right)}}\\[4pt]
\vdots \\[4pt]
{P_k^{\left( N \right)}}\end{array}} \right]\end{align}
\begin{align}{{\textbf{F}}_k} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}
{{f_{11}}} & {{f_{12}}} & {{f_{13}}}\\[4pt]
{{f_{21}}} & {{f_{22}}} & {{f_{23}}}\\[4pt]
{{f_{31}}} & {{f_{32}}} & {{f_{33}}}\end{array}} \right]\end{align}
\begin{align}{{\textbf{G}}_k} = \left[ {\begin{array}{c@{\quad}c@{\quad}c}
{{g_{11}}} & {{g_{12}}} &{{g_{13}}}\\[4pt]
{{g_{21}}} & {{g_{22}}} & {{g_{23}}}\\[4pt]
{{g_{31}}} & {{g_{32}}} &{{g_{33}}}\end{array}} \right]\end{align}
\begin{align}{{\textbf{H}}_k} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}
{{h_{11}}} & {{h_{12}}} & \cdots &{{h_{1N}}}\\[4pt]
{{h_{21}}} & {{h_{22}}} & \cdots & {{h_{2N}}}\\[4pt]
{{h_{31}}} & {{h_{32}}} & \cdots & {{h_{3N}}}
\end{array}} \right]\end{align}
The matrices
${{\textbf{F}}_k}$
,
${{\textbf{G}}_k}$
,
${{\textbf{H}}_k}$
are the coefficient of the ARMA model to be estimated. The matrix
${{\textbf{H}}_k}$
is a parameter that represents how the surface pressure field contributes to the rotational motion. In other words, nominal aerodynamic forces, gust effects and inertia properties such as mass distribution and the position of the centre of gravity are all involved in the third term of the right side of the above equation. The index
$\left( N \right)$
of
$P_k^{\left( n \right)}$
is the identifier of pressure sensors. If the parameter matrices
${{\textbf{F}}_k},{{\textbf{G}}_k}$
, and
${{\textbf{H}}_k}$
are updated so that the Equation (51) holds at every time step, the angular velocity Equation (51) results in the latest model for control design.
Assuming that the coupling between longitudinal and lateral-directional motions is negligible, some parameters can be regarded as nearly zero:
\begin{align}{f_{12}} = {f_{13}} = {f_{21}} = {f_{31}} = {g_{12}} = {g_{13}} = {g_{21}} = {g_{31}} \approx 0\end{align}
This enables the design of the control system to be independent of longitudinal and lateral-directional motions. While this paper focuses solely on longitudinal control design, the same design procedure can be applied to lateral-directional motions as well.
Now, the ARMA model for pitch rate is written as follows:
\begin{align}{q_k} = {f_{11}}{q_{k - 1}} + {g_{11}}{q_{k - 2}} + \mathop {\mathop \sum \limits_{i = 1} }\limits^N {h_{1i}}P_k^{\left( i \right)}\end{align}
Considering the symmetry in the placement of pressure taps, the following assumption can be made regarding the coefficient
$h$
because the contribution of the symmetrically located pressures to the pitch rate is considered to be equal.
\begin{align}{h_{1i}} = {h_{1\!\left( {N + 1 - i} \right)}}{\rm{\;\;\;\;}}\left( {i = 1,2, \cdots, N/2} \right)\end{align}
This reduces the number of coefficients to be estimated to
$\left( {2 + N/2} \right)$
,
${f_{11}},{g_{11}},{h_{11}},{h_{12}}, \cdots, {h_{1N/2}}$
. Here, a classical parameter estimation scheme is applied using a linear Kalman filter. A column vector
$\boldsymbol\theta $
of parameters to be estimated is defined, and the observation equation is expressed as follows:
\begin{align}{q_k} = \phi _k^{\rm{{\rm T}}}{\boldsymbol\theta _k} + {w_{q,k}}\end{align}
\begin{align}{\boldsymbol\theta _k} = {\left[ {\begin{array}{*{20}{l}}{{f_{11}}} \quad {{g_{11}}}\quad {{h_{11}}} \quad{{h_{12}}}\quad \cdots \quad {{h_{1N/2}}}\end{array}} \right]^{\rm{{\rm T}}}}\end{align}
\begin{align}{\phi _k} = {\left[ {\begin{array}{*{20}{l}}
{{q_{k - 1}}} \quad{{q_{k - 2}}}\quad {P_k^{\left( 1 \right)} + P_k^{\left( N \right)}}\quad {P_k^{\left( 2 \right)} + P_k^{\left( {N - 1} \right)}} \quad\cdots \quad{P_k^{\left( {N/2} \right)} + P_k^{\left( {N/2 + 1} \right)}}\end{array}} \right]^{\rm{{\rm T}}}}\end{align}
The
${w_{q,k}}$
is assumed to be Gaussian noise with a given normal distribution
$N\!\left( {0,\sigma _{{w_q}}^2} \right)$
. Since the prior transition model of the coefficient vector
$\boldsymbol\theta $
cannot be given, we assume a Wiener model as follows:
\begin{align}{\theta _{k + 1}} = {\theta _k} + {{\boldsymbol{v}}_{q,k}}\end{align}
The
${{\textbf{v}}_{q,k}}$
is assumed to be Gaussian noise with a given normal distribution
$N\!\left( {\textbf{0},{{\boldsymbol{\Sigma }}_v}} \right)$
. The covariance
$\sigma _{{w_q}}^2$
and
${{\boldsymbol{\Sigma }}_v}$
are hyper parameters to be given appropriately. The standard Kalman filter is applied to the transition and observation Equations (61) and (58) to estimate
${\theta _k}$
.
\begin{align}{\textbf{P}}_k^ - = {{\textbf{P}}_{k - 1}} + {{\boldsymbol{\Sigma }}_v}\end{align}
\begin{align}{{\textbf{g}}_k} = \frac{{{\textbf{P}}_k^ - {\phi _k}}}{{\phi _k^{\rm{{\rm T}}}{\textbf{P}}_k^ - {\phi _k} + \sigma _{{w_q}}^2}}\end{align}
\begin{align}{\hat \theta _k} = {\hat \theta _{k - 1}} + {{\textbf{g}}_k}\!\left( {{q_k} - \phi _k^{\rm{{\rm T}}}{{\hat \theta }_{k - 1}}} \right)\end{align}
\begin{align}{{\textbf{P}}_k} = \left( {{\textbf{I}} - {{\textbf{g}}_k}\phi _k^{\rm{{\rm T}}}} \right){\textbf{P}}_k^ - \end{align}
where
${{\textbf{P}}_k}$
is the estimation error covariance matrix and
${{\textbf{g}}_k}$
is the Kalman gain. The estimated variables are indicated by the hat symbol. Using the estimated
${\hat \theta _k}$
, the angular velocity equation for control design is written as
\begin{align}{q_{k + 1}} = {\hat f_{11}}{q_k} + {\hat g_{11}}{q_{k - 1}} + \mathop {\mathop \sum \limits_{i = 1} }\limits^{N/2} {\hat h_{1i}}\!\left( {P_k^{\left( i \right)} + P_k^{\left( {N + 1 - i} \right)}} \right) + {b_{11}}{\rm{\Delta }}{\delta _{e,k}}\end{align}
where the effect of the elevator deflection
${\rm{\Delta }}{\delta _{e,k}}$
is added, and the coefficient
${b_{11}}$
indicates the elevator effectiveness. The third term on the right side of the Equation (66) is the external aerodynamic force approximated by the pressure measurement. The time history of the external force can be AR model representation:
\begin{align}{z_{m,k}} = {a_{11}}{z_{m,k - 1}} + {a_{12}}{z_{m,k - 2}}\end{align}
Since the coefficients
${a_{11}}$
and
${a_{12}}$
are unknown parameters, the Kalman filter is applied to obtain the estimated values
${\hat a_{11}}$
and
${\hat a_{12}}$
. Equations (66) and (67) are then summarised to form the augmented system as follows:
\begin{align}{\textbf{x}}_{k + 1}^{{\rm{lon}}} = {\textbf{A}}_k^{{\rm{lon}}}{\textbf{x}}_k^{{\rm{lon}}} + {{\textbf{B}}^{{\rm{lon}}}}{\rm{\Delta }}{\delta _{e,k}}\end{align}
\begin{align}{{\textbf{x}}_k^{{\rm{lon}}} = {{\left[ {{q_{k - 1}}\ {q_k}\ q_k^{{\rm{cmd}}}\ {z_{m,k - 1}}\ {z_{m,k}}} \right]}^{\rm{{\rm T}}}}}\end{align}
\begin{align}{\textbf{A}}_k^{{\rm{lon}}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}
0 & 1 & 0& 0& 0\\[4pt]
{{{\hat g}_{11}}} & {{{\hat f}_{11}}}& 0& 0 & 1\\[4pt]
0 & 0 & 1 & 0& 0\\[4pt]
0 & 0 & 0 & 0 & 1\\[4pt]
0 & 0 & 0 & {{{\hat a}_{12}}} & {{{\hat a}_{11}}}\end{array}} \right]\end{align}
\begin{align}{{\textbf{B}}^{{\rm{lon}}}} = {\left[ {\begin{array}{*{20}{l}}0 \quad {{b_{11}}} \quad 0\quad 0 \quad 0\end{array}} \right]^{\rm{{\rm T}}}}\end{align}
The increment of elevator angle
${\rm{\Delta }}{\delta _{e,k}}$
is determined such that
${q_k} \to q_k^{{\rm{cmd}}}$
, which can be obtained as the solution of a standard linear quadratic control problem that minimises the following cost function:
\begin{align}J_{{t_f}}^{{\rm{lon}}} = \mathop {\mathop \sum \limits_{k = {t_0}} }\limits^{{t_f}} \left[ {{{\bar q}_1}{{\left( {{q_k} - q_k^{{\rm{cmd}}}} \right)}^2} + {{\bar r}_1}{\rm{\Delta }}{\delta _{e,k}}^2} \right]\end{align}
\begin{align} = \mathop {\mathop \sum \limits_{k = {t_0}} }\limits^{{t_f}} \left( {{\textbf{x}}{{_k^{{\rm{lon}}}}^{\rm{{\rm T}}}}{{\textbf{Q}}^{{\rm{lon}}}}{\textbf{x}}_k^{{\rm{lon}}} + {{\bar r}_1}{\rm{\Delta }}{\delta _{e,k}}^2} \right)\end{align}
${\bar q_1}$
and
${\bar r_1}$
are positive weight coefficients that can be chosen by the designer. Since our problem assumes uncertain atmospheric turbulence, the future estimates of parameters may also experience significant variations. Therefore, it does not make sense to set long control intervals, only a single time step ahead is sufficient.
\begin{align}J_1^{{\rm{lon}}} = {\textbf{x}}{_{k + 1}^{{{\rm{lon}}}^{\rm{{\rm T}}}}}{{\textbf{Q}}^{{\rm{lon}}}}{\textbf{x}}_{k + 1}^{{\rm{lon}}} + {\bar r_1}{\rm{\Delta }}{\delta _{e,k}}^2\end{align}
The
${\rm{\Delta }}{\delta _{e,k}}$
that minimises this is obtained by setting
$\partial J_1^{{\rm{lon}}}/\partial {\rm{\Delta }}{\delta _{e,k}} = 0$
.
\begin{align}{\rm{\Delta }}{\delta _{e,k}} = \frac{{ - 1}}{{{{\textbf{B}}^{{\rm{lon}}}}^{\rm T}{{\textbf{Q}}^{{\rm{lon}}}}{{\textbf{B}}^{{\rm{lon}}}} + {{\bar r}_1}}}{{\textbf{B}}^{{\rm{lon}}}}^{\rm{{\rm T}}}{{\textbf{Q}}^{{\rm{lon}}}}{\textbf{A}}_k^{{\rm{lon}}}{\textbf{x}}_k^{{\rm{lon}}}\end{align}
This control input is the increment of the elevator angle, and the resulting elevator command is as follows:
\begin{align}\delta _{e,k + 1}^{{\rm{cmd}}} = \delta _{e,k}^{{\rm{cmd}}} + {\rm{\Delta }}{\delta _{e,k}}\end{align}
Note that anti-windup processing should be necessary to deal with the implementation issue associated with input saturation.
A comparison of experiment and simulation result.
5.0 Evaluation of FIFC
5.1 Dynamic wind tunnel testing
5.1.1 Link-pitch control
Wind tunnel tests were conducted in a large 6.5m × 5.5m low-speed wind tunnel at JAXA. As shown in Fig. 11, the directional stability of the tail-less aircraft is quite weak, so the link-yaw can quickly reach a critical point and damage the support system. Therefore, it was necessary to limit the degrees of freedom of the support system to two: the link-pitch angle
$\xi $
and Euler pitch angle
$\theta $
.
Considering the requirement for low computational cost for real-time processing by the onboard microcomputer, four pressure taps (P2, P5, P8 and P9) on the upper surface of the left wing and symmetrical four taps (P31, P28, P25 and P24) on the right wing were selected as the pressure sensor placement for control. In the chordwise pressure distribution of typical aerofoils, it is known that the pressure variation on the leading-edge side of the upper surface is highly sensitive to changes in angle-of-attack. Based on these observations, firstly, P2, P5 and P9 are selected to be evenly spaced along the leading-edge span direction. One additional point, P8, was selected from the trailing-edge side so that it was not aligned in the chord direction with the already selected points.
Initially, the control of the link-pitch angle was carried out without any wind disturbance. To achieve the desired link-pitch angle, the pitch rate command was generated based on the difference between the link-pitch angle command and the actual link-pitch angle, along with the link-pitch rate.
\begin{align}q_k^{{\rm{cmd}}} = {K_{p,\xi }}\left( {\xi _k^{{\rm{cmd}}} - {\xi _k}} \right) + {K_{d,\xi }}{\dot \xi _k}\end{align}
where
${K_{p,\xi }}$
represents the proportional gain of the link-pitch angle error, while
${K_{d,\xi }}$
represents the differential gain of the link-pitch angle. In this experiment,
${K_{p,\xi }} = 0.05$
and
${K_{d,\xi }} = - 0.4$
were used. The link-pitch angle command was set to
$\xi _k^{{\rm{cmd}}} = \pm 20\left[ {{\rm{deg}}} \right]$
.
Figure 14 presents the longitudinal states of the aircraft when the link angle command is given (
${\xi ^{{\rm{cmd}}}} = \pm 20\left[ {{\rm{deg}}} \right]$
, green line). The link-pitch angle, aircraft pitch rate, and aircraft pitch angle are sensor outputs, while the elevator deflection angle is the commanded value from the flight computer. The high-frequency oscillation of the pitch angular velocity is due to the characteristics of the onboard gyro. Observing the changes in pitch angle, it can be seen that the pitch angle adjusts in the direction of movement to follow the link angle command. Since this is a uniform flow in the test section, the angle-of-attack is almost the same as the pitch angle.
Next, we conducted a simulation of the dynamic wind tunnel test under the same conditions as the experiment. Figure 14 presents the simulation results shown in red. While some minor differences are observed, the experiment and simulation results display nearly identical behaviour. As no rapid aircraft movement was performed during the experiment to induce unsteady flow, the flow field calculation using the panel method can sufficiently reproduce the experiment accurately. Therefore, the validity of the simulation has been confirmed. Utilising this simulator, tuning of FIFC can be carried out within the simulation environment.
5.1.2 Control performance under gusty environment
A gust-generating sheet was hung upstream of the wind tunnel to confirm the transient performance and stability in turbulent flow. The polyethylene sheet was affixed to a pole and manually moved in an up-and-down motion (see Fig. 15). It can be pulled upward at a consistent rate to generate the same turbulence repeatedly.
Before starting the experiment, the time variation of the wind disturbance was measured at the nominal aircraft position using the ultrasonic anemometer. The results are shown in Fig. 16. It can be easily seen that dynamic pressure
$\frac{1}{2}\rho {U^2}$
dropped by 25
${\rm{\% }}$
and the equivalent angle-of-attack
${\rm{\Delta }}\alpha = \frac{W}{U}$
fluctuated in the range of ±5°.
The link pitch command was set to 0 degrees without any disturbance for the initial 15 seconds, and the gust generator was activated at approximately
$t = 155$
s. This procedure is shown in Fig. 17. The vertical green line in Fig. 18 indicates the time activating gust generator.
Gust generator installed upstream of the wind tunnel. The polyethylene sheet can be pulled upward at a consistent rate to generate the same turbulence repeatedly. (a) Wind speed at 0m/s. (b) Wind speed at 8m/s. Gusts are generated by small fluttering in the polyethylene sheet.
Comparison of wind speeds with and without gusts (measured by ultrasonic anemometer at the nominal nose position). The wind component parallel to the main flow is denoted by
$U$
, and the vertical component is denoted by
$W$
. The bottom graph shows the equivalent angle-of-attack variation, which varies within approximately
$ \pm $
5 degrees.
Snapshot of activating gust generator. (Frame 1) Wind speed is 8 m/s. Gust generator is not activated. (Frame 2) Gust generator is activated and polyethylene sheet rising. (Frame 3) The polyethylene sheet reaching to the top, the aircraft is descending due to the gust. (Frame 4) The aircraft return and keep the original position by the attitude control.
The time histories of the longitudinal state variables of the aircraft and the surface pressure fluctuations during gust generation. Pressure values are shown in red for the upper surface pressure and in blue for the lower surface pressure.
The longitudinal state variables of the aircraft are shown in the upper half of Fig. 18.
It can be observed that the pitch rate and link-pitch angle experience small oscillations during the gust encounter. Even in the gust environment, the aircraft can successfully track the pitch rate command and the resulting link-pitch angle maintains
$0$
degrees.
The lower half of Fig. 18 depicts the wing surface pressure variations of several pressure taps on the upper surface in red and on the lower surface in blue. Since the lift force is generated properly, the upper surface exhibits negative pressure values, while the lower surface displays positive pressure values. The pressure values indicate that gusts increase pressure fluctuations on both the upper and lower surfaces. As one would expect, the pressure on the leading edge is more sensitive to disturbances than the pressure on the trailing edge side. Compared to a single-point AoA sensor, multiple pressure sensor networks gather gust information from the wing surface, resulting in a satisfactory transient response of the link-pitch angle
$\xi $
.
Next, experiments were conducted under the same conditions, but with the control law switched to typical proportional-integral-derivative (PID) control. In this case, the elevator command is incremented as follows:
\begin{align}{\rm{\Delta }}{\delta _{e,k}} = {K_{P,q}}{e_k} + {K_{D,q}}{\dot e_q}\end{align}
\begin{align}{e_k} = q_k^{{\rm{cmd}}} - {q_k}\end{align}
\begin{align}{\dot e_q} = \frac{{{e_k} - {e_{k - 1}}}}{{{\rm{\Delta }}T}}\end{align}
In this experiment, the gains were adjusted to
${K_{P,q}} = 0.015$
and
${K_{D,q}} = - 0.0075$
. The gain coefficients were manually tuned in nonlinear simulations, observing the transient time-response to commands. During the tuning phase, a gain increase/decrease factor of ±6 dB and a maximum time delay of 0.2 seconds were independently introduced at the input port of the closed-loop to ensure a stability margin for the closed-loop control system.
${\rm{\Delta }}T = 0.05\!\left( {\rm{s}} \right)$
is the time step for control. The
$q_k^{{\rm{cmd}}}$
is generated using the same formula Equation (77) as in FIFC. The time responses of the FIFC method and PID control are compared and shown in Fig. 19. When gusts are applied, it is inevitable that the aircraft position deviates from the target position due to the sudden decrease in dynamic pressure and the resulting loss of lift. However, the extent of the initial fall between 5 and 10 seconds is kept smaller with FIFC than with PID control. Even during the settling phase after 10 seconds, the control by FIFC is able to suppress the steady-state deviation to a smaller level. These improvements in control performance are achieved through timely fine-tuning of adaptive control using in-flight pressure field sensing rather than fixed gains.
Comparison of time response when the FIFC method and PID control are applied (the upper figure shows the link pitch angle, and the lower figure shows the elevator displacement).
5.2 Multiple case simulation
In the previous section, the performance of the FIFC method was verified through dynamic wind tunnel experiments. However, the results represent a single case, and experiments for severe gusts cannot be conducted due to the limited operating range of the support system. In this chapter, the performance of the FIFC is evaluated by free-flight simulation for a large number of gust cases to confirm the reliability of the FIFC.
Firstly, in order to check the basic behaviour of the FIFC, set a stepwise pitch angle command
${\theta ^{{\rm{cmd}}}}$
. The pitch angle command is converted to the pitch rate command by appropriately selected proportional gain
${K_{P,\theta }}$
.
\begin{align}q_k^{{\rm{cmd}}} = {K_{P,\theta }}\!\left( {\theta _k^{{\rm{cmd}}} - {\theta _k}} \right)\end{align}
where
${K_{P,\theta }}$
is determined so as to achieve satisfactory transient time response for non-gust cases and sufficient gain and delay margin. In this simulation, the gain was adjusted to
${K_{P,\theta }} = 1.5$
.
Figure 20 shows the results of the simulation. The first 40 seconds are without wind disturbance, while the second 40 seconds simulate a severe gust environment. The actual vertical gusts applied are shown at the bottom of the figure. Although the lowest left and right graphs are exactly the same, they are for the convenience of comparing the upper and lower graphs. This wind disturbance is a real outdoor measured gust and was clipped data from Mauree et al. [Reference Mauree, Deschamps, Bequelin, Loesch and Scartezzini27, Reference Mauree, Lee, Coccolo and Scartezzini28]. Databases on wind measurement are publicly available in Ref. (29).
Figure 20(a) shows typical variables of the longitudinal motion, and Fig. 20(b) is the time variation of the internal variables of the FIFC, including the various parameter estimates and the diagonal elements of the estimation error covariance matrix obtained from the Kalman filter. The figure shows a distinctive feature of FIFC. That is, gust information is acquired from surface pressure in real-time, and it appears in parameter estimation, leading to adaptive control deflection. Command tracking performance is degraded compared to the undisturbed case, but stability is maintained even in the gusty environment.
The result of the pitch angle tracking simulation. (a) the vertical state variables of the aircraft. (b) the estimated parameters (Param. Est.) and the diagonal elements of the estimation error covariance (Error. Cov.) obtained from the Kalman filter. The vertical gust is acting from
$t = 40\!\left[ {\rm{s}} \right]$
.
5.2.1 FIFC vs fixed gain control in gust environment
A comparative evaluation was conducted between FIFC and PID control systems under the gust environment. Among atmospheric disturbances, vertical disturbances, in particular, lead to variations in the angle-of-attack. Vertical disturbances can result in sudden vertical accelerations and, in the worst case, lead to a stall.
In the field of gust alleviation control, the performance is often assessed by considering the aircraft’s vertical motion. Previous studies have proposed ride discomfort indices to evaluate passenger comfort in aircraft, which are derived from vertical acceleration and pitch rate. Based on this background, the present simulation evaluates the performance using multiple cases of vertical gust samples.
PID control system was used as a typical control system for comparison. The increment of elevator angle is obtained by PID control using the following equation:
\begin{align}{\rm{\Delta }}{\delta _e} = {K_{p,{\rm{lon}}}}\!\left( {\theta _k^{{\rm{cmd}}} - {\theta _k}} \right) + {K_{d,{\rm{lon}}}}{q_k}\end{align}
Since this is the elevator angle increment, the elevator angle is the sum of these values at each point in time. In this case, anti-windup processing is used to prevent saturation. The gains were determined following the same procedure as described in 5.1.2, considering both time-response and stability margin. Eventually, the gains were adjusted to
${K_{p,{\rm{lon}}}} = 0.03$
and
${K_{d,{\rm{lon}}}} = - 0.02$
.
The gust data used in this section is the same as that used in the previous section. This data consisted of long-duration gust recordings, from which appropriate segments of 40 seconds were extracted to create multiple sample cases. Each sample case was then classified into one of three gust intensities: Light, Moderate, or Severe, based on the variance value of the gust. For statistical evaluation, 1,500 samples were prepared for each gust intensity.
Figure 21 shows a subset of the classified gust samples, featuring 20 cases. The frequency response is shown in Fig. 22. Notably, ‘severe’ gusts exhibited fluctuations of approximately 3 degrees in terms of the angle-of-attack. The frequency response diagram shows that ‘severe’ gusts exhibit higher gain across the entire frequency range compared to the other gust intensities. Additionally, a line representing the Kolmogorov-5/3 rule is plotted for reference in the same figure.
The plot of gust samples for each gust strength. Using the actual measurement data [29], each sample case was then classified into one of three gust intensities: Light, Moderate, or Severe, based on the variance value of the gust.
The frequency response diagram of each gust strength. This graph shows that Severe gusts exhibit higher gain across the entire frequency range compared to the other gust intensities. The Kolmogorov-5/3 rule line is also plotted, indicating that each gust follows this rule.
The histogram illustrates the control cost at various gust strengths. The dash line represents the mean control cost of each control method. A comparison between the FIFC method and the PID control reveals that the distribution of control costs for the FIFC method is shifted towards the left side.
The flight simulation was conducted with the angle-of-attack not exceeding
$ \pm 8{\rm{\;\;deg}}$
, to ensure the accuracy of the panel method. The pitch angle command was set to
$3{\rm{\;\;deg}}$
. To evaluate the performance of the control system, the cost function is defined as the sum of the squares of the pitch angle errors.
\begin{align}\bar{\!J} = \mathop \sum \limits_{k = 1}^n {(\theta _k^{{\rm{cmd}}} - {\theta _k})^2}\end{align}
Figure 23 presents a histogram illustrating the control cost at various gust intensities. A comparison between the FIFC method and the PID control reveals that the distribution of control costs for the FIFC method is shifted towards the left side. The results demonstrate that FIFC exhibits a smaller standard deviation, indicating a more robust performance. The dashed line represents the mean control cost of each control method. This indicates that FIFC demonstrates superior performance compared to the conventional approach, regardless of turbulence intensity. However, it is worth noting that as the gust intensity increases, the average control cost of the FIFC also increases. Consequently, the performance difference between FIFC and the PID control diminishes. This trend contradicts the intended objective of FIFC and suggests the need for further investigation and research.
6.0 Conclusions
Based on the concept of integrating in-flight pressure sensing into inertial-sensor-based flight control, the feasibility of FIFC was verified through dynamic wind tunnel tests and flight simulations. The model airplane used in the tests was a modified version of a commercial product, making it cost-effective to conduct dynamic wind tunnel tests in a gust environment. The main advantage of the proposed FIFC over conventional flight control is that it updates the state-space equations of angular motions in-flight and optimises the control system according to the gust environment, without requiring a prior aerodynamic model of the aircraft or its mass characteristics.
The verification process of the effectiveness of FIFC is as follows: First, the time response of FIFC to step-wise link pitch angle commands was observed in an undisturbed wind tunnel. Furthermore, it was compared with simulation results based on the mathematical model of the dynamic support system, and it was found that they were in close agreement. This fact not only confirms the characteristics of FIFC but also demonstrates the high fidelity of the simulation. Next, a PID control system was designed using the aforementioned simulation, and the response characteristics of FIFC and PID control were compared in a wind tunnel experiment, revealing that FIFC exhibits superior transient response and steady-state characteristics. Since the gusty environment in the wind tunnel represents just a single case, the performance of FIFC and PID control was compared in simulations for a total of 1,500 gust cases (500 light, 500 moderate and 500 severe cases). The statistical results also indicated that FIFC has the potential for higher control performance and reliability compared to fixed gain control.
The placement of multiple pressure taps on the wing surface in the FIFC method is an interesting future research topic. Based on a few attempts, the influence of the pressure tap position is not surprisingly significant in terms of its contribution to control performance, as the parameters are adjusted to represent rotational motion using given pressure outputs. In the Appendix, preliminary considerations are also made regarding the optimal layout of the pressure sensor arrangement. Pressure observation locations are selected using sparse modeling, with the optimisation index being the ability to express angular velocity and acceleration as linear functions of the pressure field. Additional optimisation for sensor placement using gust samples of severe turbulence is currently under consideration.
Acknowledgements
We would like to express our gratitude for the support of the staff at the JAXA large low-speed wind tunnel (LWT1) and to Sagamido Corporation for their valuable advice in the design and fabrication of the dynamic support system.
Appendix A. Placement of pressure sensors
A.1 Overview
This appendix introduces an initial approach to the optimal layout design of pressure sensors. The surface pressure distribution itself governs the aerodynamic forces and moments that control the aircraft’s motion. From this perspective, one idea is to select a minimal number of pressure measurement points where the motion can be effectively approximated through linear regression. In this appendix, numerical examples will be presented using actual outdoor flight data to demonstrate the selection of measurement points. Due to the vast number of possible combinations, Least Absolute Shrinkage and Selection Operator (LASSO) regression is employed to sparsely select the pressure points on the model airplane that effectively describe the vertical motions, pitch angular velocity and vertical acceleration in a linear relationship.
A.2 LASSO regression
To optimise the configuration of pressure taps, LASSO regression can be utilised, which is a commonly used method for sparse modeling. Assume the following linear equation:
\begin{align}{y_k} = {\phi ^{\rm{{\rm T}}}}{{\textbf{x}}_k}\end{align}
where
${y_k}$
represents the target variable indicating a motion-related quantity such as angular velocity or acceleration,
${{\textbf{x}}_k} \in {{\rm{R}}^m}$
denotes the explanatory input consisting mainly of multiple measured pressures, and the subscript
$k$
denotes the time step. The vector
$\phi $
represents the coefficient to be estimated. Consider the dataset
${\textbf{y}}$
and
${\textbf{X}}$
, where
${\textbf{y}}$
consists of
$n$
time series of
${y_k}$
, and
${\textbf{X}}$
is an
$n \times m$
matrix consisting of
$m$
column vectors
${{\textbf{x}}_k}$
. A linear regression model is considered that expresses
${\textbf{y}}$
using the coefficient vector
$\phi $
and
${\textbf{X}}$
. Each element of the vector
$\phi $
represents the coefficient for each pressure measurement. Having fewer nonzero coefficients implies that the target variable can be explained with fewer pressure sensors. The problem of optimally placing as few pressure sensors as possible on the wing surface, for example, on approximately 4,000 panels as shown in Fig. 10, is of great importance. It is well-known that the sparse modeling problem can be reduced to the following optimisation problem [Reference Tibshirani30–Reference Brunton and Kutz32]:
\begin{align}\mathop {{\rm{min}}}\limits_\phi \frac{1}{2}\|{\textbf{y}} - {\textbf{X}}\phi \|_2^2 + \lambda \|\phi\|_1\end{align}
That is, the objective is to minimise the cost function by adding the
${l_1}$
norm regularisation term to the squared error between the estimate and the true value. This regression method is known as LASSO regression, which yields a sparse coefficient vector
$\phi $
. The parameter
$\lambda $
represents the weight coefficient of the regularisation term. Increasing this weight results in obtaining a sparser coefficient vector, but it also leads to a larger squared error between the estimated and true values. When applying LASSO to tasks such as image denoising or reconstruction, the number of unknown variables is typically greater than the given equations, resulting in an underdetermined system. In contrast, in this research, the number of unknown variables is smaller than the given equations, resulting in an overdetermined system.
A.3 Example 1: Pitch rate and pressure
Here, the pressure sensor placement is determined in terms of expressing the pitch rate using the surface pressure field. In this case,
${{\textbf{x}}_k}$
and
${y_k}$
are written as follows:
\begin{align}{{\textbf{x}}_k} = {\left[ {\begin{array}{*{20}{l}}
{{q_{k - 1}}} \quad {{q_{k - 2}}} \quad {P_k^{\left( 1 \right)} + P_k^{\left( {32} \right)}} \quad {P_k^{\left( 2 \right)} + P_k^{\left( {31} \right)}} \quad \cdots \quad {P_k^{\left( {16} \right)} + P_k^{\left( {17} \right)}}\end{array}} \right]^{\rm{{\rm T}}}}\end{align}
\begin{align}{y_k} = {q_k}\end{align}
Note that the variables of vector
${{\textbf{x}}_k}$
were standardised [Reference Tibshirani30]. LASSO is applied to obtain a sparse coefficient vector
$\phi $
. Unlike the recursive Kalman-Filter type formulation in Chapter 4, the coefficient vector
$\phi $
is determined as a fixed value for the entire time interval. By decreasing the weight coefficient
$\lambda $
of LASSO regression, the number of selected pressure taps gradually increases.
The results of LASSO regression are shown in Fig. A1.
The result of LASSO regression for the pitch rate model.
As the weight coefficient
$\lambda $
decreases, the number of selected pressure taps increases. While increasing the number of pressure taps reduces estimation error, the rate of decrease gradually flattens after a certain number of taps. Therefore, estimation can be sufficiently performed with only a few pressure taps. The second graph illustrates which coefficients of the explanatory variables become non-zero as the weight coefficient changes. This indicates that a pressure tap that was initially selected may no longer be selected in the middle of the process. The pitch rate is represented by green, the upper-pressure taps by red and the lower-pressure taps by blue. For the number and arrangement of pressure taps, refer to Fig. 4. It is intuitive that the pitch rate at 1 and 2 steps back are important to express the current state. The first pressure tap selected is the P12 located on the lower side, followed by P13 and the upper-pressure port P9. These pressure taps are located on the leading-edge side of the wing root. Once the pitch rate at time steps
$t - 1$
and
$t - 2$
is chosen, the mean squared error (MSE) decreases and eventually reaches convergence. This indicates that the pitch rate can be predominantly represented by the auto-regressive model. This case applies to stable flights in calm atmospheres; for flights in severely disturbed atmospheres, it is expected that the contribution of the pressure term will become significant.
The third graph represents the solution path commonly used in LASSO regression, showing the changes in the coefficient values
$\phi $
in a graph. The coefficients corresponding to unselected pressure taps have values of zero. Therefore, this graph also demonstrates the gradual increase in the number of non-zero coefficients.
The result of LASSO regression for the vertical acceleration model.
A.4 Example 2: Vertical acceleration and pressure
Next, modeling of the vertical acceleration will be considered, which is one of the aircraft state variables in vertical motion. In aircraft gust alleviation, the objective is to minimise changes in vertical acceleration. If certain pressure sensors on the wing surface can capture signs of large vertical accelerations, it may be possible to achieve innovative gust alleviation with minimal delay.
In this modeling, the standardised explanatory variables and the target variables are defined as follows:
\begin{align}{{\textbf{x}}_k} = {\left[ {\begin{array}{*{20}{l}}
{P_k^{\left( 1 \right)} + P_k^{\left( {32} \right)}} & {P_k^{\left( 2 \right)} + P_k^{\left( {31} \right)}} & \cdots & {P_k^{\left( {16} \right)} + P_k^{\left( {17} \right)}}\end{array}} \right]^{\rm{{\rm T}}}}\end{align}
\begin{align}{y_k} = {a_{z,k}}\end{align}
where
${a_{z,k}}$
is the vertical acceleration of the aircraft measured by AHRS at time step
$k$
.
The results of LASSO regression are shown in Fig. A2.
Similar to the validation for pitch rate, it can be observed that the number of selected pressure taps increases as the weight coefficients decrease. Both upper and lower surface pressure taps are initially selected simultaneously, indicating the importance of pressure taps on both the upper and lower surfaces for estimating vertical acceleration. Moreover, examining the arrangement of the selected pressure taps reveals that the leading-edge pressure taps are sequentially chosen before the trailing-edge pressure taps. This suggests that the leading-edge pressure taps capture pressure fluctuations that are correlated with vertical acceleration.
In the solution path, it can be observed that as the number of pressure taps used increases, the coefficient values decrease. This is because in cases where there are few pressure taps, the weight coefficients become large in order to express the target variable using fewer pressure taps. As the number of pressure taps increases, the contribution of each pressure tap to the target variable becomes smaller.
A.5 Brief discussion
In this appendix, the optimal layout design of pressure sensors for representing pitch angular velocity and vertical acceleration is discussed using LASSO regression. It is found that the important pressure sensors differ depending on the model used to express pitch angular velocity and vertical acceleration. These findings guide designers in determining the layout of pressure sensors, considering factors such as the number of pressure taps and spatial constraints on the aircraft surface.
In addition to representing motion variables, for the given FIFC structure, pressure sensor placement can also be optimised for typical severe turbulence, which leads to a new approach to gust alleviation control.
Open access
