Loss of control in flight accident case study: icing-related tailplane stall

Abstract

In January 2017, a business jet flew in Norway on a short repositioning flight with two pilots onboard, no passengers or cargo. Initially, the take-off proceeded as normal but as the landing gear was retracted both pilots observed that the airspeed was rapidly approaching the flap limiting speed of 200kts. When the flaps were fully retracted at a height of approximately 2,100ft above ground level, the crew experienced a violent nose-down pitch motion. Control was regained at a height of approximately 170ft above ground level and, following the accident, data from the flight data recorder showed that the aircraft experienced –2.62G during the pitch upset. A tailplane stall due to icing was suspected; however, the flight data recorder, being limited to 36 parameters, was not able to confirm this. For expediency during the accident investigation process, a simplified, linear flight dynamics model was developed using Matlab/Simulink to assess static and dynamic stability for a range of tailplane efficiency factors to simulate the effects of tailplane icing.

Nomenclature

A

‘A’ matrix

AOA

angle-of-attack (deg)

B

‘B’ matrix

CFD

computational fluid dynamics

CG

centre of gravity

C_M

pitching moment coefficient

CRM

crew resource management

EU

European Union

FDR

flight data recorder

g

normal acceleration due to gravity (m/s)

h

height above ground level (feet)

ICAO

International Civil Aviation Organisation

kts

knots

LOC-I

loss of control in flight

L _W

lift due to main wing (N)

L _T

lift due to horizontal tailplane (N)

m

total aircraft mass (kg)

MIMO

multi-input multi-output system

MTOM

maximum take-off mass (kg)

NSIA

Norwegian Safety Investigation Authority

PF

pilot flying

PM

pilot monitoring

q

pitch rate (rad/s)

SISO

single-input and a single-output system

T

temperature (deg F)

TED

trailing edge down (deg)

TEU

trailing edge up (deg)

u

longitudinal speed in X-body axis (m/s)

${\bar V_h}$

tail volume ratio

w

vertical speed in Z-body axis (m/s)

${\delta _e}$

elevator deflection (deg)

$\varphi$

bank angle (deg)

$\theta$

pitch angle (deg)

$\eta$

horizontal tailplane efficiency

$\tau$

throttle input [0–1]

${\tau _e}$

angle-of-attack effectiveness of the elevator

1.0 Introduction

In January 2017, a business jet flew in Norway on a short repositioning flight with a crew of two pilots and no passengers or cargo on board. Initially, the take-off proceeded as normal but as the landing gear was retracted, the pilots observed that the airspeed rapidly approached the flap limiting speed of 200kts. As the flaps were retracted at a height above ground level of approximately 2,100ft, the crew experienced a violent nose-down pitch motion restrained only by their seat belts, as the aircraft started banking sharply to the left. Following the accident, a full investigation was conducted by the Norwegian Safety Investigation Authority, supported by industry and academic partners [1] in accordance with international regulation [2, 3]. It is likely that the commander (pilot flying) and first officer (pilot monitoring) experienced different levels of startle and/or surprise during the upset. Control was regained at a height of approximately 170ft above ground level. Data from the flight data recorder showed that the aircraft experienced –2.62G during the pitch down upset and +5.99G during the pull-out. A tailplane stall due to icing was suspected, however the flight data recorder – limited to 36 parameters – was not able to verify this. The aircraft involved – a Cessna 560 Encore – was manufactured in 2003 and is type is certified for two pilots, commander and co-pilot with a cabin capacity of seven passengers.

1.1 The flight

The crew, which consisted of a Commander and a First Officer, had flown from Bern, Switzerland to Gardermoen, Norway with a passenger on board. After disembarking the passenger at Gardermoen, the aircraft was scheduled to fly to its home base at Sandefjord Airport Torp with no passengers on board. The crew planned to make the ground stop as short as possible and, if the weather conditions permitted, they would avoid de-icing. During the ground stop at Gardermoen, only one engine was stopped while the First Officer completed an external inspection of the aircraft. He did not observe any ice or anything out of the ordinary on the areas of the aircraft that could be inspected. According to the crew, the snow did not accumulate on the wings before departure, they could only see melted water on the wing surfaces and therefore decided not to de-ice the aircraft. When the crew requested taxi clearance, they were assigned a different runway than expected. This change entailed a longer taxi time thus longer exposure to the prevailing weather conditions. The aircraft’s ground stop lasted approximately 15 minutes at an air outside temperature of 0°C. The taxiways and runway were covered with 3–6mm of slush and it was snowing when the aircraft took off. After flying from Switzerland for more than two hours in approximately minus 50°C, the aircraft’s surfaces (fuselage and wings) were more than likely chilled.

Initially, the take-off proceeded as normal as indicated by the FDR plot (Fig. 1) [1]. The landing gear was retracted and both pilots observed that the speed was rapidly approaching 200kt, the maximum allowable speed with flaps deployed. As the flaps were retracted, the crew experienced a violent nose-down movement and the pilots were “hanging by their seat belts”, while the aircraft started banking sharply to the left. Following the accident, data from the FDR showed that the aircraft at this moment experienced –2.62G. The Commander (PF) did not trust the instruments while the First Officer (PM) who did not experience ‘startle’, likely to have had a higher level of situation awareness quickly took control and started a pull-out from the dive. The aircraft descended below the cloud base, and even though it was dark, the pilots could see the ground level. The control was regained, and the aircraft levelled off 170ft above ground level. Analysis of the FDR showed the aircraft was overstressed to +5.99G during the pull-out. The crew made a MAYDAY call to Air Traffic Control. The engines were left in take-off position during the entire pull-out and the speed increased to 325kt. Once control was regained, the MAYDAY was cancelled, and the flight continued towards Torp where an approach and landing took place without further problems. The NSIA investigation did not reveal any technical malfunctions in the aircraft or it’s control systems.

Figure 1. Extract from flight data recorder output: X-axis represents count of subframes/time and Y axis (left) individual FDR parameters (colour coded).

The aircraft’s anti and de-icing systems on the wings and tailplane were switched on and the aircraft’s tailplane rubber de-icing “boots” were in automatic mode and inactive during the take-off and when the event occurred. NSIA’s assessment is that the systems were not suitable to remove this relevant type of ice and snow. This accident shows the significance of functioning crew resource management (CRM) in the cockpit when an unexpected and extreme flight situation occurs. In this instance, the First Officer’s situation awareness and initial pull-out saved the crew and aircraft.

1.2 Review of loss of control in flight accidents where icing was contributory/causal factor

Loss of control in flight (LOC-I) or the departure from controlled flight continues to threaten flight safety. It is the primary category of fatal accidents in both commercial aviation (CA) [4] and general aviation (GA) [5]. Business aviation (BA) is considered a sub-sector of GA and in the United States – the world’s largest single aviation market – 65% of flights are conducted for business or public services [6].

Within the CA and GA aviation sectors, weather related environmental factors are frequently cited as either causal or contributory factor(s) [7, 8] to LOC-I events. Icing is frequently cited as a contributory factor, particularly in BA and CA where aircraft are certified for flight in icing conditions [9].

In CA involving commercial jets and turboprops, during the period 1981 to 2010, there were 31 icing-related events resulting in 765 fatalities. Twelve occurred during take-offs with ice contaminated wings. Twenty upsets during up and away flight resulted in 243 fatalities, two of which involved tailplane stalls [10].

The effects of icing on main wing and tailplane aerodynamics and the relationship to stability and control is known but the ability to distinguish between the two events is less understood by pilots, hence regulatory authorities have published guides to raise pilot awareness [11, 12].

1.3 Previous methods for the analysis of the effects of icing on aircraft performance and flying qualities

Multiple methods have been developed to better understand how icing affects large and small aircraft performance and flying qualities. These methods have included ground and flight-based experimental methods, computational methods and flight dynamic modelling and simulation [13].

Broeren et al., developed a database of airplane aerodynamic characteristics with simulated ice accretion using a 12ft low speed wind tunnel and a 3.5% scale generic transport model for a range of configurations and angle-of-attack and sideslip ranges [14]. Ratvasky et al., investigated the effects of tailplane icing on aircraft performance and aerodynamics with a combination of icing wind tunnel testing, dry air aerodynamic wind tunnel testing, flight testing using a De Havilland Twin Otter and analytical code development [15].

Ratvasky and Ranaudo used flight test generated data to obtain the stability and control derivatives of the baseline airplane and artificially ‘iced’ the configuration using ice shapes for different thrust conditions. Parameter identification manoeuvres were performed to determine the derivatives in a range of flight conditions and configurations. The use of ice shapes reduced elevator effectiveness and longitudinal static stability decreased substantially as a result [16].

More recently, Cao et al., used an engineering prediction method for the estimation of stability and control derivatives using individual component computational fluid dynamic calculation and narrow strip theory. Predicted results compared favourably with actual flight test results generated using system identification methods [17].

Sibilski et al., used a classical model of an airplane to analyse the flight dynamics due to icing to describe an aircraft’s motion in support of air accident investigation. Equations of motion were defined in vector form and generalised aerodynamic force equations were developed. The effects of icing on aerofoil characteristics using experimental the results of wind tunnel testing [18] were used to estimate the change in lift coefficient and pitching moment due to icing on the main wing only. The analysis suggested that the deterioration of the aircraft’s flight properties is rapid and gives the pilot little advanced notice of subsequent events. Six degree of freedom computational models developed by Bragg et al. showed that large changes in airspeed, angle-of-attack and elevator angle as ice accretes for constant power and altitude [19].

Time and cost constraints for the accident investigation – the subject of this paper – forced the development of a simple method to enable parametric analysis of the effects of tailplane icing using a limited data set available from the flight data recorder, radar and weather reports.

1.4 Special investigations: modelling & simulation

A review of loss of control in flight accidents where icing was contributory/causal factor was conducted. Theory related to the effects of icing on main wing and tailplane aerodynamics was undertaken and the relationship to stability and control investigated. Due to the limited flight data parameter set, it was decided to use modelling to simulate the effects of tailplane icing and compare with available, known data. Several modelling methods were considered including computational fluid dynamics, scale model wind tunnel testing, scale model flight testing etc., but due to time and cost constraints associated with these methods, it was decided to develop a representative (similar but not exact) model using the AAA aircraft design software and desktop mathematical modelling and simulation software. A ‘generic business jet’ linear flight dynamics model was implemented in Matlab/Simulink using aircraft geometry, mass and balance, initial flight conditions from the recorded flight data and stability and control derivatives for a similar business jet aircraft (Cessna 550 – Fig. 2). The model used was two dimensional, isolated for the longitudinal (dominant) motions about the lateral axis only. The model was linearised about the trim condition obtained from the FDR data. The airframe was treated as a rigid body. Aircraft static and dynamic stability of the generic business jet was assessed for a range of tailplane efficiency factors to simulate the effects of tailplane icing.

Figure 2. Generic business jet.

2.0 Method

In order to assess the static and dynamic stability, modelling of the total aircraft pitching moment is needed and individual contributions of all major components are required, not simply tail lift and wing lift.

2.1 Static stability

The total aircraft pitching moment is given by Ref. 20:

(1) \begin{equation}{{\bf{C}}_{\bf{m}}} = {{\bf{C}}_{{{\bf{m}}_{\bf 0}}}} + {\rm{\;}}{{\bf{C}}_{{{\bf{m}}_{\boldsymbol{\alpha }}}}} + {\rm{\;}}{{\bf{C}}_{{{\bf{m}}_{{{\bf{i}}_{\bf{h}}}}}}} + {{\bf{C}}_{{{\bf{m}}_{{{\boldsymbol{\delta }}_{\bf{e}}}}}}}\end{equation}

where ${C_{{M_0}}}$ is the value of C _M for α = i _h = δ _e = 0

(2) \begin{equation}{{\bf{C}}_{{{\bf{m}}_{\bf 0}}}} = {{\bf{C}}_{{{\bf{m}}_{{\bf{a}}{{\bf{c}}_{{\bf{wf}}}}}}}} + {\rm{\;}}{{\bf{C}}_{{{\bf{L}}_{{{\bf 0}_{{\bf{wf}}}}}}}}\!\left( {{\bf{h}} - {{\bf{h}}_{\bf 0}}} \right) + {{\bf{C}}_{{{\bf{L}}_{{{\boldsymbol{\alpha }}_{\bf{h}}}}}}}{\boldsymbol{\eta }}\frac{{{{\bf{S}}_{\bf{T}}}}}{{\bf{S}}}\!\left( {{{\bf{h}}_{\bf 0}} - {\bf{h}}} \right){{\boldsymbol{\varepsilon }}_{\bf 0}}\end{equation}

and ${C_{{M_\alpha }}}$ = , the change in total aircraft pitching moment coefficient with a change in angle-of-attack

(3) \begin{equation}{{\bf{C}}_{{{\bf{m}}_{\boldsymbol{\alpha }}}}} = {{\bf{C}}_{{{\bf{L}}_{{{\boldsymbol{\alpha }}_{{\bf{wf}}}}}}}}\!\left( {{\bf{h}} - {{\bf{h}}_{\bf 0}}} \right) - {{\bf{C}}_{{{\bf{L}}_{{{\boldsymbol{\alpha }}_{\bf{h}}}}}}}{\boldsymbol{\eta }}\frac{{{{\bf{S}}_{\bf{T}}}}}{{\bf{S}}}\!\left( {{{\bf{h}}_{\bf 0}} - {\bf{h}}} \right)\left( {{\bf 1} - \frac{{{\bf{d}{\boldsymbol\varepsilon }}}}{{{\bf{d}\boldsymbol{\alpha }}}}} \right)\end{equation}

and ${C_{{M_{{i_h}}}}}$ = , the change in total aircraft pitching moment coefficient with a change of horizontal tailplane angle, i_h for α = δ _e = 0

(4) \begin{equation}{{\bf{C}}_{{{\bf{m}}_{{{\bf{i}}_{\bf{h}}}}}}} = - {{\bf{C}}_{{{\bf{L}}_{{{\boldsymbol{\alpha }}_{\bf{h}}}}}}}{\boldsymbol{\eta }}\frac{{{{\bf{S}}_{\bf{T}}}}}{{\bf{S}}}\!\left( {{{\bf{h}}_{\bf 0}} - {\bf{h}}} \right) = - {{\bf{C}}_{{{\bf{L}}_{{{\boldsymbol{\alpha }}_{\bf{h}}}}}}}{\boldsymbol{\eta }}{\bar{\textbf{V}}_{\bf{h}}}\end{equation}

where ${{\bar{\textbf{V}}}_{\bf{h}}}$ , the tail volume ratio is given by:

(5) \begin{equation}{{\bar{\textbf{V}}}_{\bf{h}}} = \frac{{{{\bf{S}}_{\bf{T}}}}}{{\bf{S}}}\!\left( {{{\bf{h}}_{\bf 0}} - {\bf{h}}} \right)\end{equation}

and ${C_{{M_{{\delta _e}}}}}$ = , the change in total aircraft pitching moment coefficient with a change of elevator angle, δ _e for α = i _h = 0

(6) \begin{equation}{{\bf{C}}_{{{\bf{m}}_{\boldsymbol{\delta }}}_{\bf{e}}}} = - {{\bf{C}}_{{{\bf{L}}_{{{\boldsymbol{\alpha }}_{\bf{h}}}}}}}{\boldsymbol{\eta }}{{\bar{\textbf{V}}}_{\bf{h}}}{{\boldsymbol{\tau }}_{\bf{e}}}\end{equation}

and horizontal tailplane efficiency, η ranges from 0 to 1.

Thus, the contributions to the pitching moment fall into two major categories [21], contributions that vary with angle-of-attack and this that do not vary with angle-of-attack. The tailplane efficiency (η) can be used to model the degradation of lift due as a result of tailplane icing.

The total aircraft pitching moment about the aircraft’s centre of gravity [22] consists of contributions from wing, tail and fuselage (Fig. 3). Each contribution generates moments that vary with angle-of-attack and contributions that are independent of angle-of-attack (constant). A negative total pitching moment slope represents positive static stability – the aircraft returning to the trimmed flight condition following a disturbance (e.g., sudden change of tail lift).

Figure 3. Contributions to total aircraft pitching moment about the CG.

As ice builds up on the tail, it becomes less aerodynamically efficient, increased negative tail lift is needed to keep maintain the trimmed flight condition and prevent the aircraft tending to nose down and overspeed with flaps extended. This is achieved by increasing elevator trailing edge up to compensate (Fig. 4).

Figure 4. Forces and moments on an aircraft and effect of elevator deflection (adapted from Ref. 23).

When flaps are retracted (Fig. 5), the lift generated by the wing decreases and the point of lift moves forward. The aircraft tends to nose UP as the nose down pitching moment due to wing lift decreases. Increased trailing edge down (TED) elevator is required to compensate. At the tailplane, downwash angle decreases hence negative tail lift decreases. The aircraft tends to nose DOWN as the pitching moment due to negative tail lift decreases. This time, increased trailing edge UP (TEU) elevator is required to compensate. The resultant pitching moment about the CG is the net effect of these contributions. The effects due to tailplane are usually dominant (Fig. 3).

Figure 5. Forces and moments on an aircraft and effect of flap retraction (adapted from Ref. 23).

2.2 Dynamic stability

The longitudinal dynamic stability characteristics of the aircraft are dependent upon the airspeed, pitching moment variation with angle-of-attack, moment of inertia, lift to drag ratio and aerodynamic tail damping. The system may be exhibit positive, neutral or negative dynamic stability [24].

Given the trimmed flight condition and static stability of the generic business jet, dynamic analysis can be undertaken to consider the effects of small disturbances (perturbations) such as turbulence (external) or control inputs (internal). Using defined aircraft notation and axes state variable, control inputs and matrix/vector notations are defined. Systems with more than one input and more than one output are known as multi-input multi-output systems (MIMO). Systems that have only a single-input and a single-output are defined (SISO). The aircraft in longitudinal and pitching motion maybe modelled using a state-space model. The longitudinal aircraft dynamics are linearised about the setpoint (airspeed, pitch angle, pitch rate, elevator deflection and throttle setting) and can be written in state space form [22]:

(7) \begin{equation}{\dot{\textbf{X}}} = {{\bf{A}}_{{\bf{long}}}}{\bf{X}} + {{\bf{B}}_{{\bf{long}}}}{\bf{u}}\end{equation}

The elements of matrix A are stability derivatives describing the effect of state variables on forces and moments. The elements of matrix B are control derivatives representing the effects of elevator and throttle commands on the body referenced forces and moments.

2.3 Transfer functions, elevator to pitch

Using this matrix method, transfer functions (input-output relationships) are derived using desktop computing mathematical modelling tool Matlab to determine the relationship between Input: Elevator Deflection ( ${{\rm{\delta }}_e}$ ) and Output: Pitch Angle ( $\theta$ ).

2.4 Effects of flap retraction

The effects of flap retraction are simulated using a Simulink switching model. This type of model enables the dynamic analysis to account for changes in stability and control derivatives as a result of flap configuration changes. Given stability and control derivatives for the selected aircraft, transfer functions are estimated for different flap configurations and tailplane efficiencies.

3.0 Modelling and simulation results

The commercial AAA aircraft software design package (using applied theory described earlier) was used for modelling of static stability using a generic business jet in a range of conditions similar to the accident aircraft. Dynamic analysis was conducted using custom developed Matlab models and Simulink.

3.1 Static stability

For the generic business jet, using the given flight condition of V= 204 KTAS, H = 4,000ft pressure height, $\textrm{T} = 29^{\circ}\textrm{F}$ , aft CG/low MTOM with flap = 7 degrees, the variation of pitching moment with angle-of-attack was obtained for a range of horizontal tailplane efficiencies from 1.0 (100%) to 0.5 (50%) to simulate the effects of icing on the aircraft tailplane (Fig. 6).

Figure 6. Pitch stability – pitching moment vs angle-of-attack @204kts.

The results show that the pitching moment versus angle-of-attack gradient decreases as the horizontal tailplane efficiency decreases, hence aircraft static stability also decreases as horizontal tailplane efficiency decreases.

The results show that the pitching moment versus angle-of-attack gradient decreases as the horizontal tailplane efficiency decreases, hence aircraft static stability also decreases as horizontal tailplane efficiency decreases. This trend is also apparent by examining Equation (3) showing that the tail is associated with a negative (hence stabilising) contribution to the pitching moment coefficient, C _Mα.

The results also show that the angle-of-attack for trimmed flight conditions increases as tailplane efficiency decreases. Analysis of horizontal tailplane efficiency versus stick-fixed static margin suggests that static margin decreases as the horizontal tailplane efficiency decreases and that the aircraft is neutrally (statically) stable when the horizontal tailplane efficiency is approximately 0.68 for the generic business jet model in given flight conditions (Fig. 7).

Figure 7. Tailplane efficiency factor vs static margin.

Further analysis of elevator deflection versus horizontal tailplane efficiency (Fig. 8) suggests that increasing UP elevator (–ve elevator deflection according to convention) is required to maintain trimmed flight as the horizontal tailplane efficiency decreases. The range of elevator deflection for the generic business jet model was 20 degrees UP and 15 degrees DOWN. The results suggest that as horizontal tailplane efficiency decreases below approximately 0.2 (20%), there is insufficient UP elevator to maintain trimmed flight.

Figure 8. Elevator deflection vs tailplane efficiency.

3.2 Dynamic stability

Having established the trimmed flight condition and static stability of the generic business jet the dynamic analysis was undertaken to consider the effects of small disturbances (perturbations) such as turbulence (external) or control inputs (internal). Before conducting the dynamic analysis, the aircraft, the system of notation and axes were defined using the right-hand rule [25]. State variable, control inputs and matrix/vector notations were defined as follows:

(8) \begin{equation}\boldsymbol{m}{\dot{\boldsymbol{V}}_{\boldsymbol{I}}} = \boldsymbol{m}\!\left( {{{\dot{\boldsymbol{V}}}_{\boldsymbol{B}}} + {\boldsymbol\omega _{\boldsymbol{B}/\boldsymbol{I}}} \times {\boldsymbol{V}_{\boldsymbol{B}}}} \right) = {\boldsymbol{F}_{\boldsymbol{Aero}}} + {{\boldsymbol{F}}_{\boldsymbol{gravity}}} + {\boldsymbol{F}_{\boldsymbol{T}}}\end{equation}

(9) \begin{equation}\boldsymbol{I}{\dot{\boldsymbol\omega}_{\boldsymbol{I}}} = \boldsymbol{I}{\dot{\boldsymbol\omega}_{\boldsymbol{B}}} + {\boldsymbol\omega _{\boldsymbol{B}/\boldsymbol{I}}} \times \boldsymbol{I}{\boldsymbol\omega _{\boldsymbol{B}}} = {\boldsymbol{M}_{\boldsymbol{Aero}}} + {\boldsymbol{M}_{\boldsymbol{gravity}}}\end{equation}

where F _T is the throttle force that only acts along the X axis (other components on thrust due to misalignment were verified to be negligible). The throttle acceleration will be denoted $\tau = \dfrac{{{F_T}}}{m}$ longitudinal motion only affects the X and Z axes and the pitch rotation (Fig. 9), with v=r=p=Y=L=N=0:

(10) \begin{equation}\dot{\boldsymbol{U}} = - {\boldsymbol{g}}{\bf{sin\boldsymbol\Theta }} - {\boldsymbol{QW}} + \frac{{{\boldsymbol{X}} + {{\boldsymbol{F}}_{\boldsymbol{T}}}}}{\boldsymbol{m}}\end{equation}

(11) \begin{equation}\dot{\boldsymbol{W}} = - {\boldsymbol{g}}{\bf{cos\boldsymbol\Theta }} + {\boldsymbol{QU}} + \frac{{\boldsymbol{Z}}}{{\boldsymbol{m}}}\end{equation}

(12) \begin{equation}\dot{\boldsymbol{q}} = \frac{{\boldsymbol{M}}}{{{{\boldsymbol{I}}_{{\boldsymbol{yy}}}}}}\end{equation}

Figure 9. Aircraft trim condition at the time of the event (Adapted from Ref. 23).

Equations are linearised around a constant pitch $\theta$ ₀ and constant longitudinal velocity u ₀, with q ₀=w ₀=0 representing a zero rate of change in pitch and altitude.

The trim condition is given by:

(13) \begin{equation}{\boldsymbol{X}_{\textbf{0}}} - \boldsymbol{mg}{\bf{sin}}\!\left( {{\boldsymbol\theta _{\bf{0}}}} \right) = \bf{0}\end{equation}

(14) \begin{equation}{\boldsymbol{Z}_{\bf{0}}} + \boldsymbol{mg}{\bf{cos}}\!\left( {{\boldsymbol\theta _{\bf 0}}} \right) = {\bf 0}\end{equation}

(15) \begin{equation}{\boldsymbol{M}_{\bf 0}} = {\bf 0}\end{equation}

More precisely, U,W,Q and Θ are linearised around u ₀, w ₀, q ₀, $\theta_0$ such that:

\begin{equation*}\boldsymbol{U} = {\boldsymbol{u}_{\bf 0}} + {\bf{u}}\end{equation*}

\begin{equation*}{\bf{W}} = {\boldsymbol{w}_{\bf 0}} + {\bf{w}}\end{equation*}

\begin{equation*}{\bf{Q}} = {\boldsymbol{q}_{\bf 0}} + {\bf{q}}\end{equation*}

\begin{equation*}{\boldsymbol{\Theta }} = {\boldsymbol\theta _{\bf 0}} + \theta \end{equation*}

Θ is replaced by $\theta$ ₀+ $\theta$

By taking the differences and substituting Θ by $\theta$ ₀+ $\theta$ , we have:

(16) \begin{align}\dot{\boldsymbol{u}} & = \frac{\bf 1}{\boldsymbol{m}}\!\left( {\left( {\frac{{\boldsymbol\partial \boldsymbol{X} + \boldsymbol\partial {\boldsymbol{X}_{\boldsymbol{T}}}}}{{\boldsymbol\partial \boldsymbol{u}}}} \right)\boldsymbol{u} + \frac{{\boldsymbol\partial \boldsymbol{X}}}{{\boldsymbol\partial \boldsymbol{w}}}\boldsymbol{w} - \boldsymbol{mg}\boldsymbol\theta {\bf{cos}}{\boldsymbol\theta _{\bf 0}}} \right)\nonumber\\[4pt]\dot{\boldsymbol{w}} & = \frac{\bf 1}{\boldsymbol{m}}\!\left( {\frac{{\boldsymbol\partial \boldsymbol{Z}}}{{\boldsymbol\partial \boldsymbol{u}}}\boldsymbol{u} + \frac{{\boldsymbol\partial \boldsymbol{Z}}}{{\boldsymbol\partial \boldsymbol{w}}}\boldsymbol{w} + \frac{{\boldsymbol\partial \boldsymbol{Z}}}{{\boldsymbol\partial \boldsymbol{q}}}\boldsymbol{q} - \boldsymbol{mg}\boldsymbol\theta {\bf{sin}}{\boldsymbol\theta _{\bf 0}} + {\boldsymbol{u}_{\bf 0}}\boldsymbol{q}} \right)\nonumber\\[4pt]\dot{\boldsymbol{q}} & = \frac{\bf 1}{\boldsymbol{I}_{\boldsymbol{yy}}}\!\left( \frac{\boldsymbol\partial \boldsymbol{M}}{\boldsymbol\partial \boldsymbol{u}}\boldsymbol{u} + \frac{\boldsymbol\partial \boldsymbol{M}}{\boldsymbol\partial \boldsymbol{w}}\boldsymbol{w} + \frac{\boldsymbol\partial \boldsymbol{M}}{\boldsymbol\partial \dot{\boldsymbol{w}}} \dot{\boldsymbol{w}} + \frac{\boldsymbol\partial \boldsymbol{M}}{\boldsymbol\partial \boldsymbol{q}}\boldsymbol{q} \right)\\[4pt]\dot{\boldsymbol\theta} & = \boldsymbol{q}\nonumber\end{align}

3.3 MIMO (multiple input multiple output) aircraft longitudinal state space model

Systems with more than one input and more than one output are known as multi-input multi-output systems (MIMO). Systems that have only a single input and a single output are defined (SISO). The aircraft in longitudinal and pitching motion maybe modelled using a state-space model. The longitudinal aircraft dynamics, linearised about the setpoint (V= 204 KTAS, $\theta$ = 24 degrees pitch), can be written as a state space model:

(17) \begin{equation}{\dot{\textbf{X}}} = {{\bf{A}}_{{\bf{long}}}}{\bf{X}} + {{\bf{B}}_{{\bf{long}}}}{\bf{u}}\end{equation}

Matrix/vector notations in bold

A: State matrix, depends on stability derivatives and setpoint condition

B: Input matrix, depends on control derivatives and actuator layout

X: State vector, components are u, w, q, $\theta$

u: Control input vector, components are η, τ

where:

\begin{equation*}{\bf{X}} = \left[ {\begin{array}{*{20}{c}}u\\[2pt]w\\[2pt]q\\[2pt]\theta \end{array}} \right],{\bf{u}} = \left[ {\begin{array}{*{20}{c}}{{\delta _e}}\\[2pt]\tau \end{array}} \right]\end{equation*}

u: X axis speed, w : Z axis speed.

q: pitch rate, $\theta$ : pitch, ${{\rm{\delta }}_e}$ : elevator deflection, τ : throttle command.

Throttle history was recorded within flight data but elevator history was not.

Using matrix notation, the state transition is given by:

(18) \begin{equation}{\bf{A}} = \left[ {\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}{{\boldsymbol{X}_{\boldsymbol{Tu}}} + {\boldsymbol{X}_{\boldsymbol{u}}}} & {{\boldsymbol{X}_{\boldsymbol{w}}}} & {{{\bf{X}}_{\boldsymbol{q}}} - \boldsymbol{m}{\boldsymbol{w}_{\bf 0}}}& { - \boldsymbol{g}{\bf{cos}}\!\left( {{\boldsymbol\theta _{\bf 0}}} \right)}\\[5pt]{\dfrac{{{\boldsymbol{Z}_{\boldsymbol{u}}}}}{\boldsymbol{f}}}& {\dfrac{{{{\bf{Z}}_{\boldsymbol{w}}}}}{\boldsymbol{f}}} &{{\boldsymbol{Z}_{\boldsymbol{q}}} + \dfrac{{{\boldsymbol{u}_{\bf 0}}}}{\boldsymbol{f}}}& { - \dfrac{{\boldsymbol{g}{\bf{sin}}\!\left( {{\boldsymbol\theta _{\bf 0}}} \right)}}{\boldsymbol{f}}}\\[13pt]{{{\bf{M}}_{{\bf{Tu}}}} + {{\bf{M}}_{\boldsymbol{u}}} + {{\bf{M}}_{{\dot{\textbf{w}}}}}\dfrac{{{{\bf{Z}}_{\boldsymbol{u}}}}}{\boldsymbol{f}}{\rm{\;}}}& {{{\bf{M}}_{\boldsymbol{w}}} + \dfrac{{{{\bf{M}}_{{\bf{T\alpha }}}}}}{{{\boldsymbol{u}_{\bf 0}}}} + {{\bf{M}}_{\dot{\boldsymbol{w}}}}\dfrac{{{{\bf{Z}}_{\boldsymbol{w}}}}}{\boldsymbol{f}}}& {{{\bf{M}}_{\dot{\boldsymbol{w}}}}\dfrac{{{{\bf{Z}}_{\boldsymbol{w}}}}}{{\bf{f}}} {{\bf{M}}_{\boldsymbol{q}}} + {{\bf{M}}_{\dot{\boldsymbol{w}}}}\dfrac{{{{\bf{Z}}_{\boldsymbol{q}}} + {{\bf{u}}_{\bf 0}}}}{\boldsymbol{f}}}& { - {{\bf{g}}{\dot{\textbf{w}}}}\dfrac{{\bf{sin}\! \left( {{\boldsymbol\theta _{\bf 0}}} \right)}}{\boldsymbol{f}}}\\[13pt]{\bf 0} &{\bf 0} &{\bf 1} &{\bf 0}\end{array}} \right]\end{equation}

where ${\boldsymbol{f}} = {\bf 1} - {{\boldsymbol{Z}}_{\dot{\boldsymbol{w}}}}$

All stability derivatives are in dimensional form and represent the partial derivatives of a force or moment with respect to a state variable, divided by mass in the case of forces or inertia in the case of moments, for example:

\begin{equation*}{{\boldsymbol{X}}_{\boldsymbol{u}}} = \frac{\bf 1}{\boldsymbol{m}}\frac{{\boldsymbol\partial \boldsymbol{X}}}{{\boldsymbol\partial \boldsymbol{u}}},{\rm{\;}}{{\boldsymbol{M}}_\boldsymbol{q}} = \frac{\bf 1}{{{\boldsymbol{I}_{\boldsymbol{yy}}}}}\frac{{\boldsymbol\partial \boldsymbol{M}}}{{\boldsymbol\partial \boldsymbol{q}}}, \textrm{etc}.\end{equation*}

The dimensional stability derivatives were derived from the non-dimensional stability derivatives of the aircraft using the equations defined by Roskam [26].

The non-dimensional derivatives were lift, drag, pitching moment coefficients ${C_{L1}},{C_{D1}},{C_{M1}}$ respectively, defined at speed setpoint u₀ and pitch setpoint $\theta$ ₀, and ${C_{L0}},{C_{D0}},{C_{M0}}$ that are constants depending on aircraft geometry.

1. The throttle coefficients: ${C_{TX0}},{C_{MT0}}$

2. Speed dependent coefficients: $\;{C_{Lu}},{C_{Du}},{C_{Mu}},{C_{TXu}},{C_{MTu}}$

3. Force coefficients defined with respect to angle-of-attack and pitch rate: $\;{C_{D\alpha }},{C_{L\alpha }},{C_{L\dot \alpha }},{C_{Lq}}$

4. Moment coefficients defined with respect to angle-of-attack and pitch rate: $\;{C_{M\alpha }},{C_{M\dot \alpha }},{C_{Mq}},{C_{MT\alpha }}$

5. These were derived using aircraft design software for stability and control analysis [27].

The control matrix is given by:

(19) \begin{equation}{\bf{B}} = \left[ {\begin{array}{c@{\quad}c}{{{\boldsymbol{X}}_{{\boldsymbol\delta _{\boldsymbol{e}}}}}} & {{{\boldsymbol{X}}_{\boldsymbol\tau} }}\\[9pt]{\dfrac{{{{\boldsymbol{Z}}_{{\boldsymbol\delta _{\boldsymbol{e}}}}}}}{{\boldsymbol{f}}}}& {\bf 0}\\[13pt]{{{\bf{M}}_{\boldsymbol\eta} } + {{\bf{M}}_{\dot{\boldsymbol{w}}}}\dfrac{{{{\boldsymbol{Z}}_{{\boldsymbol\delta _{\boldsymbol{e}}}}}}}{{\boldsymbol{f}}}} & {\bf 0}\\[13pt]{\bf 0} & {\bf 0}\end{array}} \right]\end{equation}

where the dimensional control derivatives were obtained from the non-dimensional control derivatives. ${C_{D{\delta _e}}},{C_{L{\delta _e}}},{C_{M{\delta _e}}}$ of the aircraft at pitch and speed setpoint.

The elements of matrix A are stability derivatives describing the effect of state variables on forces and moments. The elements of matrix B are control derivatives representing the effects of elevator and throttle commands on the body referenced forces and moments.

Note: The B1 type model requires computation of flap aerodynamic parameters, which are not available. The simulations were generated by computing the A and B for different tailplane stall efficiency settings.

3.4 Transfer functions, elevator to pitch

Aircraft parameters (Appendix A, Table A1) were used define elements in respective matrices. Using this matrix method, transfer functions were derived using Matlab [28] to determine the relationship between Input: Elevator Deflection ( ${{\rm{\delta }}_e}$ ) and Output: Pitch Angle ( $\theta$ ). Transfer function coefficients were derived from expressions for the aircraft stability derivatives [26].

For the given flight condition (Appendix A, Table A1):

The elevator to pitch transfer function, denoted G(s), with 100% horizontal tailplane efficiency was found to be:

(20) \begin{equation}{\boldsymbol{G}_{\boldsymbol{elevator} \to \boldsymbol{pitch}{\rm{\;}}\left( {\bf{100}{\rm{\% \;}}{\bf{tail}}{\rm{\;}}{\bf{efficiency}}} \right)}}\!\left( \boldsymbol{s} \right) = \frac{{ - \bf{17.39}{\boldsymbol{s}^{\bf 2}} - {\bf{52.13}}\boldsymbol{s} - \bf{1.128}}}{{{\boldsymbol{s}^{\bf 4}} + \bf{5.521}{\boldsymbol{s}^{\bf 3}} + \bf{12.5}{\boldsymbol{s}^{\bf 2}} + \bf{0.3387}\boldsymbol{s} + \bf{0.2689}}}\end{equation}

where the negative numerator coefficients are consistent with a pitch down response to elevator down deflections and the positive coefficients are consistent with a stable response.

Decreasing the horizontal tailplane efficiency factor was found to significantly change the intrinsic characteristics of the transfer function. For example, a 50% tailplane efficiency factor results in the following transfer function:

(21) \begin{equation}{\boldsymbol{G}_{\boldsymbol{elevator} \to \boldsymbol{pitch}{\rm{\;}}\left( {50{\rm{\% \;}}{\bf{tail}}{\rm{\;}}{\bf{efficiency}}} \right)}}\!\left( \boldsymbol{s} \right) = \frac{{\bf{8.736}{\boldsymbol{s}^{\bf 2}} - \bf{26.17}\boldsymbol{s} - \bf{0.6315}}}{{{\boldsymbol{s}^{\bf 4}} + \bf{4.226}{\boldsymbol{s}^{\bf 3}} - \bf{0.03886}{\boldsymbol{s}^{\bf 2}} + \bf{0.0396}\boldsymbol{s} + \bf{0.1132}}}\end{equation}

The negative coefficient (–0.03886) in the denominator for Equation (15) is linked to an unstable pole, and the pitch response to elevator commands is therefore unstable when horizontal tailplane efficiency is reduced to 50%. Lower values of the tailplane efficiency factor were found to further increase instability.

Figure 10. Short period oscillation (SPO) with $100\%$ horizontal tailplane efficiency.

Figure 11. Long period oscillation (LPO) with $100\%$ horizontal tailplane efficiency.

3.5 Short period oscillation (SPO) and long period oscillation (LPO)

For model validation purposes, the longitudinal dynamics were analysed by excitation of the short period (Fig. 10) and long period (Fig. 11) modes using the eigenvectors of A associated with each mode to specify initial conditions as follows (i is the mode number):

(22) \begin{equation}{{\bf{x}}_{\bf{0}}} = {{\boldsymbol{E}}_{\boldsymbol{i}}}{{\boldsymbol{T}}_{\boldsymbol{i}}} + {\boldsymbol{E}}_{\boldsymbol{i}}^{\rm{*}}{\boldsymbol{T}}_{\boldsymbol{i}}^{\rm{*}}\end{equation}

where E _i and T _i are respectively the eigenvalue and eigenvector associated with the mode i, with i=1 for the phugoid mode and i=2 for the short period mode.

The variables simulated were relative changes (variations with respect to a trim condition). Zero value indicates the initial trimmed flight condition and all angles are in radians.

The results show that for a horizontal tailplane efficiency of 100%, the aircraft is statically and dynamically stable with heavy and moderate damping for the SPO and LPO, respectively. In addition, the model was independently verified by comparing the modal characteristics of the full state space model to the reduced order models of the two modes.

When the efficiency is reduced to 80%, the response for the LPO is also stable although with less damping of oscillations than the 100% efficiency case (Fig. 12). The results of the dynamic stability are in agreement with those of the static stability analysis.

Figure 12. Long period oscillation (LPO) with 80% horizontal tailplane efficiency.

Instability was found to start from approximately 50% tailplane efficiency factor, which was in line with the static stability analysis.

3.6 The effects of flap retraction

The effects of flap retraction were simulated using a Simulink switching model. This model enables the dynamic analysis to account for changes in stability and control derivatives as a result of flap configuration changes. Stability and control derivatives were determined for the generic business jet model using the commercial aircraft design software as before. Transfer functions were estimated for four different conditions using the switching model (Table 1).

Table 1. Switching conditions

Figure 13. Pitch response with tailplane stall then flap retraction.

For flap retraction following a tailplane stall, initially, the tailplane efficiency is 100% and the flap angle is 7 degrees (Fig. 13). At t=4 seconds, a tailplane stall is applied by a step reduction of 50% to horizontal tailplane efficiency, which destabilises the system response. The aircraft pitches down a further -30 degrees within 2 seconds. At t= 6 seconds, the retraction of flaps from 7 degrees to 0 degrees helps initially, but insufficiently to stabilise the response with only 50% tailplane efficiency. Note that the pitch response in Fig. 13 effectively represents an asymptotic pitch plot. The true pitch is limited by the upper and lower bounds of this asymptotic plot at the model transition times (t=4 seconds and t=6 seconds) and approaches this asymptotic plot at other times. This is a limitation of using a switched linear system model with a Matlab/Simulink implementation in this instance. The model switches between the discrete stability and control derivatives for each condition (tailplane efficiency and flap angle) resulting in discontinuities. This leads to discrete change in the elevator to pitch transfer function parameters. The asymptotic plot however still provides a good understanding of the pitch variations around the transition times, with an increasingly accurate description of the pitch variations and dynamics at other times.

3.7 Tailplane efficiency factor effects

At 100% horizontal tailplane efficiency, a large 20deg elevator down during 1 second is needed to initiate a large magnitude but stable phugoid response, which is only shown during the first 20 seconds for comparison purposes (Fig. 14). At 50% efficiency, a 1 second 2-degree elevator down input destabilises the system with growing oscillations and at 20% efficiency, a 1-degree elevator down is sufficient to produce very fast divergence without oscillations.

Figure 14. Pitch response with varying tailplane efficiency and elevator inputs.

Load factor changes with 1 second elevator down commands for varying horizontal tailplane efficiencies were also presented (Fig. 15). With 100% efficiency, load factor remains very close to 1 with a small elevator command, as expected. With 50% tailplane efficiency, load factor is reduced but changes are small. Major changes to load factor are however obtained with 20% tailplane efficiency and negative G is quickly reached in this case. This is believed to be closer to the conditions during the flight incident.

Figure 15. Load factor with varying tailplane efficiency and elevator inputs.

4.0 Discussion of results

The model input conditions were defined using approximated data obtained from the flight data recorder prior to the upset leading to loss of control in flight.

The results show that the gradient of total aircraft pitching moment versus angle-of-attack decreases as the horizontal tailplane efficiency decreases, hence aircraft static stability also decreases as ice builds up simulated by reduction in horizontal tailplane efficiency. Similarly, the stick-fixed static margin also decreases as horizontal tailplane efficiency decreases further confirming a reduction in static stability as ice builds up.

As the tailplane becomes less aerodynamically efficient, the aircraft will tend to nose down, increasing airspeed – as observed during the event when the aircraft approached flap limiting speed. Additional negative tail lift is needed to compensate by using increased up elevator, however at very low tailplane efficiency levels, the aircraft may run out of elevator authority thus the aircraft may continue to pitch downwards as observed during the event.

Baseline tests show that with 100% tailplane efficiency, the aircraft is statically and dynamically stable for the given flight condition. As tailplane efficiency is reduced to around 50%, the aircraft is statically and dynamically unstable.

The use of the switching model to simulate flap retraction following a tailplane stall, shows that although the aircraft pitches nose up, it continues to nose down at a rapid rate (approximately 20 degrees per second) since it is dynamically unstable.

The overall effects of tailplane efficiency on the aircraft dynamic response during the modelling and simulation are simulated by small elevator (down) perturbation, replicating possibly external disturbances such as turbulence. However, during the actual event, the rapid pitch down was observed shortly after the full retraction of flaps. Retraction of flaps reduces downwash at the tail, thus reducing the effective angle-of-attack seen locally at the tailplane and reduces the negative tail lift. It is this action that may have triggered the resultant nose down pitching motion resulting in a rapid descent rate given that the aircraft was statically and dynamically unstable at the time.

The results show that static stability decreases as tailplane efficiency decreases, simulating the onset of tailplane icing and that increasing UP elevator (-VE) is required to compensate. At low tailplane efficiency the aircraft pitch response to elevator commands becomes unstable. Flap retraction initially helps, but not sufficiently to stabilise the response with low tailplane efficiency and a tailplane stall destabilises the system response.

The build-up of ice resulted in the aircraft becoming neutral or negatively stable due the reduction in effectiveness of the horizontal tailplane. In this condition, the aircraft may diverge and either nose up or down following a disturbance. In the case of this aircraft, the retraction of the flaps (a disturbance) resulted in the aircraft nosing down (increased +ve stability) causing the aircraft to diverge nose down and not nose up.

5.0 Conclusions

The aim of this study was to provide further insight into loss of control inflight (LoC-I)/upset events in icing conditions. The main objective was to identify the probable characteristics of a LOC-I/upset event due to tailplane icing for a generic business jet. The lack of available stability, control and aerodynamic data for a specific aircraft make/model resulted in a generic business jet model being used for all analyses. Therefore, it has not been possible to replicate exact aircraft dynamics as evidenced by FDR data using modelling and simulation techniques. Flight data analysis and weather reports were used to determine flight conditions to be assessed, static and dynamic stability was assessed using established flight dynamics theory and modelling.

The modelling and what-if trends analysis does however illustrate similar trends to the recorded flight data, particularly in the case of a severe tailplane stall. The degradation/severity of tailplane aerodynamic characteristics due to icing was simulated using an assumed reduction in tailplane efficiency factor and classical theory supported by the commercial AAA aircraft design software package.

The results are applicable only for short time periods after a given disturbance since a linearised flight model was used about a trimmed flight condition, no pilot control inputs were available (e.g., yoke pitch/roll, rudder) and no external (environmental) disturbance data were available (e.g., turbulence). The results demonstrate that the generic business jet aircraft used in the analysis is statically and dynamically stable when horizontal tailplane efficiency is high. When horizontal tailplane efficiency is reduced (simulating a tailplane stall), the aircraft is statically and dynamically unstable, smaller and shorter elevator commands produce large pitch responses and negative G may be quickly reached within a short time period.

As with most incidents and accidents there are multiple contributing factors. NSIA determined that a probable explanation for the aircraft’s sudden dive is that the tailplane stalled as a result of icing caused by contamination from slush and spray from the runway and/or from falling sleet and snow. An infrequent event but one nonetheless from which lessons may be learnt in the interests of improving flight safety. For example, tailplane stall due to icing is a real threat; a tailplane stall due to icing is a well-known cause for accidents during icing conditions, the incident highlights the importance of de-icing before take-off in icing conditions and the incident underscores the deficiency of the pneumatic boot de-icing systems.