Longitudinal Flight Control of a Fixed-Wing Aircraft

Overview

• In this research, a dynamic model is established for a representative general-purpose fixed-wing aircraft based on flight performance parameters and fundamental aerodynamic parameters. The aircraft's inherent modes and flying qualities are analyzed to understand its characteristics and identify its limitations. A control system based on the C* index is designed for the aircraft, with subsequent performance analysis conducted to complete the control system and improve overall flight performance
  

  
  

• Aircraft mode and flight quality analysis are primarily based on the full-period dynamic model, with the source code available in allPeriod.m
  
  

• Flight control, root locus, stability margins, handling quality, and time-domain characteristics analysis are mainly based on the short-period model, with the source code available in shortPeriod.m and the Simulink model in shortPeriodSimulink.slx

1. Modeling

Parameter Determination

1. Typical Fixed-Wing Civil Aircraft Parameters

• Aircraft weight: $W = 12224 \ N​$
  

• Moments of inertia: $I_x=1420.9 \mathrm{~kg} \cdot \mathrm{~m}^2, I_y=4067.5 \mathrm{~kg} \cdot \mathrm{~m}^2, I_z=4786.0 \mathrm{~kg} \cdot \mathrm{~m}^2, I_{x z}=0 \mathrm{~kg} \cdot \mathrm{~m}^2$
  

• Wing area: $S=17.1 \mathrm{~m}^2$
  

• Mean aerodynamic chord: $c=1.74 \mathrm{~m}$
  

• Wingspan: $b=10.18 \mathrm{~m}$

2. Basic Flight Parameters

• Aircraft flying at sea level at $M a=0.158$, $C_{L *}=0.41, C_{D *}=0.05$
  

• Angle of attack derivatives: $C_{L \alpha}=4.44, C_{D \alpha}=0.33, C_{m \alpha}=-0.683​$
  

• Velocity derivatives: $C_{L V}=0.0, C_{D V}=0.0, C_{m V}=0.0$
  

• Dynamic derivatives: $C_{L \dot{\alpha}}=C_{z \dot{\alpha}}=0.0, C_{m \dot{\alpha}}=-4.36, C_{L q}=-C_{z q}=3.80, C_{m q}=-9.96$
  

• Control derivatives: $C_{L \delta_e}=-C_{z \delta_e}=0.355, C_{D \delta_e}=0.0, C_{m \delta_e}=-0.923$
  

• Assuming no variation in thrust with velocity, then $T_V = 0$
  

• Gravity acceleration: $g = 9.81 \ m/s^2$
  

• Atmospheric density: $\rho = 1.225 \ kg/m^3$

3. Additional Parameters

• Aircraft mass: $m = \frac{W}{g}$
  

• Reference velocity: $V_* = Ma \times 340 $
  

• Reference dynamic pressure: $q_* = \frac{\rho {V_*}^2}{2}$
  

• Assuming zero wind speed, level straight-line flight with no sideslip $\gamma_* = 0$ ，maintaining constant altitude

Model Formulation

The modeling is based on the simplified form of the longitudinal perturbation equations, neglecting the influence of altitude variations in perturbation motion and considering only the elevator control for longitudinal flight control. The system is represented in matrix form as:

$$\frac{\mathrm{d}}{\mathrm{d} t}\left[\begin{array}{c}\Delta V \\ \Delta \alpha \\ \Delta q \\ \Delta \theta\end{array}\right]=\left[\begin{array}{cccc}X_V & X_\alpha+g & 0 & -g \\ -Z_V & -Z_\alpha & 1 & 0 \\ \bar{M}_ V-\bar{M}_ {\dot{\alpha}} Z_V & \bar{M}_ \alpha-\bar{M}_ {\dot{\alpha}} Z_\alpha & \bar{M}_ q+\bar{M}_ {\dot{\alpha}} & 0 \\ 0 & 0 & 1 & 0\end{array}\right] \left[\begin{array}{c}\Delta V \\ \Delta \alpha \\ \Delta q \\ \Delta \theta\end{array}\right]+\left[\begin{array}{c}X_{\delta_{\mathrm{e}}} \\ -Z_{\delta_{\mathrm{e}}} \\ \bar{M}_ {\delta_{\mathrm{e}}}-\bar{M}_ {\dot{\alpha}} Z_{\delta_{\mathrm{e}}} \\ 0 \end{array}\right]\Delta \delta_{\mathrm{e}}$$
This formulation involves 11 longitudinal aerodynamic derivatives, given by:

$X_V = \frac{T_V \cos \left( \alpha_* + \varphi_T \right) }{m} - \frac{\left( C_{D V} + 2 C_{D*} \right)}{m V*}$

$X_\alpha=\frac{-T* \sin \left(\alpha*+\varphi_T\right)-C D \alpha q* S}{m}$

$X_{\delta_{\mathrm{e}}}=-C_{D \delta_{\mathrm{e}}} \frac{q * S}{m} $

$Z_V=\frac{T_V \sin \left(\alpha*+\varphi\right)}{m V*}+\frac{\left(C_{L V}+2 C_{L*}\right) q* S}{m V*^2}$

$Z_\alpha=\frac{\left(C_{D*}+C_{L \alpha}\right) q* S}{m V*}=\frac{T* \cos \left(\alpha* + \varphi_T\right)+C_{I \alpha} q* S}{m V*}$

$Z_{\delta_{\mathrm{e}}}=C_{L \delta_{\mathrm{e}}} \frac{q* S}{m V*}$

$\bar{M}_ V = \frac{ (C_{mV} + 2 C_{m*}) q * Sc}{2}$

$\bar{M}_ \alpha=C_{m \alpha} \frac{q* S c}{I_y}$

$\bar{M}_ {\dot{\alpha}}=C_{m \dot{\alpha}}\left(\frac{c}{2 V*}\right) \frac{q* S c}{I_y} $

$\bar{M}_ q=C_{m q}\left(\frac{c}{2 V*}\right) \frac{q* S c}{I_y}$

$\bar{M}_ {\delta_{\mathrm{e}}}=C_{m \delta_e} \frac{q* S c}{I_y} $

2. Inherent Modal and Flight Quality Analysis

Inherent Modal Analysis

1. Short-Period Mode (Approximate Model)

• Natural frequency: $\omega_{n, \ sp }=\sqrt{-\left(\bar{M}_ a+\bar{M}_ q Z_\alpha\right)} = 3.6138$
  

• Damping ratio: $\zeta_{\mathrm{sp}}=-\frac{\bar{M}_ q+\bar{M}_ {\dot{\alpha}}-Z_\alpha}{2 \omega_{\mathrm{n} . \mathrm{sp}}} = 0.6954$
  

• Period: $T_{sp} = \frac{2\pi}{\omega_{n, \ sp}} = 2.4194$
  

• Half decay time: $t_{1 / 2, sp}=-\frac{\ln 2}{\eta_{sp}} = 0.275$
  

• Oscillation cycles in half decay period: $N_{1 / 2, sp}=\frac{\ln 2 \sqrt{1-\xi^2_{sp}}}{2 \pi \xi_{sp}} = 0.114$

2. Phugoid Mode

• Natural frequency: $\omega_{\mathrm{n} , \ \mathrm{p}}= \sqrt{\eta^2_p + \omega_p^2} = 0.2137$
  

• Damping ratio: $\zeta_{\mathrm{p}}= \frac{-\eta_p}{w_n, \ p} = 0.0798$
  

• Period: $T_p = \frac{2\pi}{\omega_{n, \ p}} = 29.4923$
  

• Half decay time: $t_{1 / 2, p}=-\frac{\ln 2}{\eta_p} = 40.6338$
  

• Oscillation cycles in half decay period: $N_{1 / 2, p}=\frac{\ln 2 \sqrt{1-\xi^2_{p}}}{2 \pi \xi_{p}} = 1.379$

Inherent Flight Quality Analysis

1. Short-Period Mode

Firstly, perform a damping ratio analysis for the short-period mode

$0.35 < \zeta_{\mathrm{sp}} = 0.6954 < 1.3$

The damping ratio for the short-period mode meets the requirements

The flight quality index studied in this project, namely the $CAP$ and $C^*$ index, are specifically focused on the short-period mode

A servo actuator is integrated into the aircraft to establish a closed-loop system ( the short-period model used here will be introduced in Section 5 )

It can be obtained through calculation ( can be made using the index calculation formulas provided in Section 4 ) that $\zeta_{\mathrm{sp}} (0.1566)$, $CAP (136.2)$, $C^*$ index is significantly poor, and the actuator output oscillates

2. Phugoid Mode

For the phugoid mode, the analysis mainly focuses on the damping ratio

$\zeta_p = 0.0798 \ge 0.04$

The damping ratio for the long-period mode satisfies the requirements

3. The Research Objective

From above analysis, it can be concluded that the aircraft's inherent modal are well-designed, with good ability to resist disturbances and maintain balance. However, the $\zeta_{\mathrm{sp}}$, $CAP$, $C^*$ index and actuator output of the short-period mode is very poor and needs to be improved through the control system

Therefore, a control system will be designed specifically for the short-period mode to not only further optimize the $CAP$ index but also greatly improve the $C^*$ index, thereby enabling the aircraft to achieve excellent longitudinal flying qualities through the designed control system

4. $C^*$ Control Law Design

$C^* = n_n+\frac{V_{\mathrm{co}}}{g} q$

Here, $V_{co}$ represents the crossover speed, generally ranging from $120 \sim 132 m/s$，here selected as $V_{co} = 122 m/s$

The control law is designed for the error in $\Delta C^*$, the error is represented as:

$e(t) = \Delta C^{**} - \Delta C^*$

where $\Delta C^{**}​$ is the desired value of $\Delta C^*​$

PI control will be applied to the design of the $C^*$ control law, in a form similar to the following:

$u(t) = {K_P e(t) + K_I \int e(t) dt}$

The actuator model is defined as:

$G_\delta(s) = \frac{-1}{0.1s+1}$

5. Control System Design and Analysis

Control System Design

The dynamics model is simplified to the short-period model, focusing on elevator control for aircraft:

$$\left[\begin{array}{c}\Delta \dot{\alpha} \\ \Delta \dot{q}\end{array}\right]=\left[\begin{array}{cc}-Z_\alpha & 1 \\ \bar{M}_\alpha-\bar{M}_{\dot{\alpha}} Z_\alpha & \bar{M}_q+ \bar{M}_{\dot{\alpha}}\end{array}\right]\left[\begin{array}{c}\Delta \alpha \\ \Delta q\end{array}\right] + \left[\begin{array}{c}-Z_{\delta_e} \\ \bar{M}_{\delta_e}-\bar{M}_{\dot{\alpha}} Z_{\delta_e}\end{array}\right] \Delta \delta_e$$
The relationship between the change in normal load and the change in angle of attack is:

$\Delta n_n = \frac{C_{L\alpha} q S \Delta \alpha}{W}$

The output of the closed-loop system is:

$\Delta C^* = \Delta n_n +\frac{V_{\mathrm{co}}}{g} \Delta q$

The control system consists of a Stability Augmentation System (SAS) and a Control Augmentation System (CAS):

• S.A.S
  

  • Normal load feedback with an amplifier gain of $K_{n1}$
    

  • Pitch damper with an amplifier gain of $K_{q1}$ and a washout network time constant of $1$
    

• C.A.S
  

  • PI controller with a proportional coefficient $K_P$ and an integral coefficient $K_I$
    

  • Control amplifier gain $K_a$
    

  • Command feedforward compensator with an amplifier gain $K_{ff}$, which adjusts a zero in the system without altering the poles
    

Establish the following control system in Simulink:

The output $\Delta C^*$ of the control system is shown below:

The $C^*$ index converges stably to the input value

$\Delta q$ is depicted in the following graph:

The $\Delta q$ converges stably to $3.1^\circ/s$

$\Delta \delta_e​$ is depicted in the following graph:

$\Delta \delta_e$ converges stably to $-1.8^\circ$

Root Locus Analysis

Since feedforward compensator $K_{ff}$ does not affect closed-loop poles, it is temporarily neglected for root locus analysis, considering $K_{n1} = 0.05, K_{q1} = 0.01, K_P = 0.001, K_I = 0.04$

The open-loop transfer function of the system is:

$G_{ol}(s) = K^* \frac{s^3 + 43.817s^2 + 155.497s + 112.68}{s^5 + 16.0256s^4 + 79.6081s^3 + 203.0028s^2 + 137.325s}$

where $K^* = K_G^* K_H^* = 1.4836 K_a K_H^* = 1.4836 K_a$

Among them,

$K^*$ is the open-loop root locus gain

$K_G^*$ is the forward path root locus gain

$K_H^*$ is the feedback path root locus gain

$K_a$ is the control amplifier gain

Root locus plots are as follows:

From the root locus diagram, it can be seen that when $0 < K^* < 12.5$, the root locus remains in the left half of the $s$-plane, indicating that the system is stable

Stability Margin Analysis

Considering $K_a = 0.35, K_{ff} = 0.03$

$G_{ol}(s) = \frac{45.027(s+1) (s+0.4613) (s+2.817)}{(s+9.925) (s+0.9876) (s^2 + 5.113s + 14.01)}​$

Nyquist plots are shown below:

Bode plots are shown below:

The gain margin is $\inf$, and the phase margin is $105.1°$

Flight Quality Analysis

1. $\zeta_{sp}$ and $CAP$ Index

The closed-loop transfer function is determined to be:

$G_{cl}(s) = \frac{45.03s^3 + 192.6s^2 + 206.1s + 58.51}{s^5 + 16.03s^4 + 80.14s^3 + 225.8s^2 + 218.2s + 58.54}$

This results in a natural frequency for the short-period mode of $\omega_{n, sp} = 3.67$ and a damping ratio of $\zeta_{sp} = 0.61$

The Control Anticipation Parameter ($CAP$) is calculated as:

$CAP = \frac{\omega_{\mathrm{n}, \mathrm{sp}}^2}{(V \star / g) Z_\alpha} = 1.21$

This study primarily focuses on the control during the cruising phase of the aircraft, which is categorized under Type B flight phase. The short-period flight quality requirements for this phase are illustrated in the figure below:

It can be observed that the aircraft meets Level 1 flight quality requirements during the cruising phase (Moreover, the designed control system can also meet Level 1 flight quality under Type C flight phase)

2. $C^*$ Index

$C^* = C^(0) + \Delta C^$

Here, $C^(0)$ represents the initial value of $C^$ index and $C^*(0) = 1$ because the aircraft is in steady level flight

When the input is $1g$, the $C^*$ index is depicted in the following figure:

Here, $C^_\infty = C^(0) + \Delta C^{**} = 1 + 1 = 2$. This figure demonstrates that the $C^*$ index of the aircraft meets the envelope constraints

Time-Domain Characteristics Analysis

1. Dynamic Performance Indicators

The closed-loop poles are shown in the following figure:

The closed-loop system has five poles, including three real poles and one pair of complex conjugate poles:

$P_1 (-0.4405, 0)$

$P_2 (-0.9783, 0) $

$P_3 ( -2.254, 2.893)$

$P_4 ( -2.254, -2.893)$

$P_5 (-10.1, 0)$

There are three zeros in the system:

$Z_1 (-2.82, 0)$

$Z_2 (-1, 0)$

$Z_3 (-0.46, 0)$

The pole $P_1$ and zero $Z_3$ , as well as the pole $P_2$ and zero $Z_2$, are located too close to each other near the imaginary axis, which leads to mutual attenuation of their effects on the system response. The distance of pole $P_5​$ from the imaginary axis is relatively greater compared to $P_3, P_4​$. Therefore, the system response is primarily determined by the conjugate pole pair $P_3, P_4​$ and the zero $Z_1​$

The complex conjugate poles $P_3， P_4$ represent an underdamped system with natural frequency $w_n = 3.67$ and damping ratio $\zeta = 0.61$

The zero $Z_1$ mainly reduces peak time, accelerating system response, but increases overshoot

Response to Step Input:

• Rise Time:

  Rise time of the underdamped system: $t_{r1}=\frac{\pi-arccos\zeta}{\omega_d} = 0.77 \ s$
  

  Actual rise time: $t_r = 0.31 \ s$
  

• Peak Time:

  Peak time of the underdamped system: $t_{p1} = 1.09 \ s$
  

  Actual peak time: $t_p = 0.74 \ s$
  

• Overshoot:

  Overshoot of the underdamped system: $\sigma_1 \% = \mathrm{e}^{-\pi \zeta / \sqrt{1-\zeta^2}} \times 100 \% = 8.7 \%$
  

  Actual overshoot $\sigma \%= 18 \%$
  

• Settling Time:

  Settling time of the underdamped system: $t_{s1} = \frac{3.5}{\zeta \omega_n} = 1.55 \ s$
  

  Actual settling time: $t_s = 2.45 \ s$
  

2. Steady-State Performance Indicators

The open-loop system has a pole at the origin, categorizing it as a Type I system, which allows it to track a step input with zero steady-state error