Next Article in Journal
Graph Convolutional Networks with Bidirectional Attention for Aspect-Based Sentiment Classification
Next Article in Special Issue
Robust Quadrotor Control through Reinforcement Learning with Disturbance Compensation
Previous Article in Journal
EMG Based Gesture Recognition Using the Unbiased Difference Power
Previous Article in Special Issue
Horizontal Wind Effect on the Aerodynamic Performance of Coaxial Tri-Rotor MAV
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Bifurcation Flight Dynamic Analysis of a Strake-Wing Micro Aerial Vehicle

by
Mirosław Nowakowski
1,*,
Krzysztof Sibilski
2,*,
Anna Sibilska-Mroziewicz
3 and
Andrzej Żyluk
1
1
Air Force Institute of Technology, 01-494 Warsaw, Poland
2
Faculty of Power and Aviation Engineering, Institute of Aviation Engineering and Applied Mechanics, Warsaw University of Technology, 00-665 Warsaw, Poland
3
Faculty of Mechatronics, Warsaw University of Technology, 02-525 Warsaw, Poland
*
Authors to whom correspondence should be addressed.
Appl. Sci. 2021, 11(4), 1524; https://doi.org/10.3390/app11041524
Submission received: 10 January 2021 / Revised: 31 January 2021 / Accepted: 2 February 2021 / Published: 8 February 2021
(This article belongs to the Special Issue Unmanned Aerial Vehicles)

Abstract

:
Non-linear phenomena are particularly important in -flight dynamics of micro-class unmanned aerial vehicles. Susceptibility to atmospheric turbulence and high manoeuvrability of such aircraft under critical flight conditions cover non-linear aerodynamics and inertia coupling. The theory of dynamical systems provides methodology for studying systems of non-linear ordinary differential equations. The bifurcation theory forms part of this theory and deals with stability changes leading to qualitatively different system responses. These changes are called bifurcations. There is a number of papers, the authors of which applied the bifurcation theory for analysing aircraft flight dynamics. This article analyses the dynamics of critical micro aerial vehicle flight regimes. The flight dynamics under such conditions is highly non-linear, therefore the bifurcation theory can be applied in the course of the analysis. The application of the theory of dynamical systems enabled predicting the nature of micro aerial vehicle motion instability caused by bifurcations and analysing the post-bifurcation microdrone motion. This article presents the application of bifurcation analysis, complemented with time-domain simulations, to understand the open-loop dynamics of strake-wing micro aerial vehicle model by identifying the attractors of the dynamic system that manages upset behaviour. A number of factors have been identified to cause potential critical states, including non-oscillating spirals and oscillatory spins. The analysis shows that these spirals and spins are connected in a one-parameter space and that due to improper operation of the autopilot on the spiral, it is possible to enter the oscillatory spin.

1. Introduction

Classic methods of testing dynamic aircraft stability enable analysing their dynamic properties in the framework of minor disturbances of a steady straight flight. However, Micro Aerial Vehicles (MAV) operating in open space, where instantaneous gusts of wind can exceed 10 m/s (which amounts for 25% of their cruising velocity) are exposed to sudden flight parameter changes, the angles of attack and slip, in particular. At strong gusts of wind, the disturbance velocity reaches values comparable to the flight velocity. This issue has been reviewed in, among others, the publications [1,2,3]. In average weather conditions, the changes of attack angle are sudden, and their amplitude amounts from +60° to −30°. The velocity changes from 20 to 130 km/h, and the altitude from 250 to 25 m [3]. This is why there was a need to analyse the dynamic properties of micro aerial vehicles over the entire operational range, including the range of near-critical and super-critical angles of attack. Due to the high non-linearity of the differential equations describing flight dynamics, well-developed and described methods of modal analysis cannot be applied in this case. There was a need to develop a new method for testing aircraft motion stability, which would enable studying its dynamic properties throughout the entire range of usable angles of attack and slip. In the second half of the 20th century, researchers suggested a continuation method based on the bifurcation theory that enabled analysing dynamic properties of an aircraft, over a wide range of flight states. Articles on this matter, published in the 1970s include: [4,5]. The cited research work concerned the analysis of a fast inertial barrel roll and spin. Numerous studies devoted to the bifurcation analysis of aircraft flight dynamics and aerodynamics appeared in the 1970s, 80s and 90s. These include the AGARD report from 1985 [6], and many papers on the issues of non-stationary aerodynamics and non-linear dynamics of flight in the perspective of the dynamical system theory and the bifurcation theory: [7,8,9,10,11,12,13,14].
The dynamical system theory enabled a global analysis of a system of strongly non-linear ordinary differential equations describing the state of a dynamical system, and in this aspect it is a generalization of the aforementioned dynamic stability analysis involving a steady straight flight of an aircraft. The first step in this method is the assessment of quasi-steady state stability. The quasi-steady state is determined by equating the derivatives of state vector to zero and solving a system of non-linear algebraic equations. The local stability of a dynamical system can be assessed based on the Hartman-Grobman [15] theorem. This stability is determined by the eigenvalues of a locally linearised equations of MAV flight dynamics. In our case, this system is linearized around the equilibrium position, [7,8,9,10,11,12,13,14]. If even one of eigenvalues has a positive real part, then the equilibrium position is unstable. It can be proved that if a linearized equation system is non-singular, then the steady state of a dynamical system is a continuous function of state parameters [15]. It can also be proved that steady states of differential equations describing an aircraft motion are continuous functions of rudder surface deflections [16]. Stability changes occur when at least one eigenvalue of a locally linearized differential equation system of aircraft motion changes its sign. The changes in the steady state stability lead to a qualitatively various system response and are called bifurcations. Stability limits can be determined through searching for eigenvalues with zero real zero part [17,18]. In the bifurcation theory, the stability of quasi-steady states is tested as a function of the so-called bifurcation parameters [15,19,20]. The usual assumption in the course of a flight dynamics bifurcation analysis is that such parameters are the rudder surface deflections within the range of change from the minimum to the maximum value of these angles. It enables obtaining an image of all steady states (or quasi-steady) of aircraft motion. The source literature often refers to this computation process as “global analysis” (cf. work [10,12,21,22,23,24,25,26,27,28,29]). This is why, this type of aircraft flight analysis will be called a “global analysis of equilibrium position stability”. Of course, classic states defining the stability of a steady and straight micro aerial vehicle flight are one of the points determined in the course of a global analysis of equilibrium position stability. In this perspective, this analysis is a generalization of the classical aircraft stability analysis.
This article presents the application of bifurcation analysis, complemented with time-domain simulations, to understand the open-loop dynamics of strake-wing micro aerial vehicle model by identifying the attractors of the dynamic system that manages upset behaviour. A number of factors have been identified to cause potential critical states, including non-oscillating spirals and oscillatory spins. The analysis shows that these spirals and spins are connected in a one-parameter space and that due to improper operation of the autopilot on the spiral, it is possible to enter the oscillatory spin.
The Cobra maneuver was also studied from the point of view of the bifurcation theory and the results of the time histories of the flight parameters during this maneuver were presented.

2. Bifurcation and Continuation Analysis

A bifurcation can be defined as a qualitative change in the system dynamics as a parameter is varied. In flight dynamics, a qualitative change is usually understood as a change in the stability of the aircraft. Bifurcation analysis and continuation methods are based on the principles of the dynamical system theory. The basis of the dynamical systems theory are described in many books (it can be mention here: [15,16,17,20,30,31]). These methods consist in finding and tracking solutions numerically, in a selected range of parameters, in order to generate bifurcation diagrams. Bifurcation diagrams allow to identify qualitative changes in the dynamic response of the system.

2.1. Bifurcation Theory and Bifurcation Types

In our case bifurcation analysis is applied to autonomous dynamical system of general form:
x ˙ = f ( x , μ )
where x n is vector of n state variables, μ m is vector of m control parameters, f is nonlinear vector field governing system dynamics. The bifurcation for Equation (1) is every qualitative change of the phase portrait occurring upon the passage of parameter µ through a certain point µ0. Point µ0 is called the bifurcation point for Equation (1). By establishing the set of parameters μ = μ 0 and selecting the initial conditions x(t = 0) = x0, the system of differential Equation (1) can be integrated for the selected initial condition. In this way, one can study the evolution of vector field x over time. In order to thoroughly study the dynamics of the system, this process can be repeated an infinite number of times starting from different initial conditions and for an innumerable fixed combination of parameter values µ. This task is exhausting and tedious. The major problem in this exercise is the selection of initial conditions, which is a non-trivial task when dealing with systems that are nonlinear. An alternative and more efficient approach to the analysis of nonlinear dynamical systems described by Equation (1) is based on the asymptotic bifurcation and continuation method. The bifurcation and continuation method begins with the calculation of steady states of the equilibrium type of the Equation (1), which comes down to solving a system of nonlinear algebraic equations:
x ˙ = 0     f ( x , μ ) = 0
and computing the eigenvalues of the Jacobi matrix:
J = f x
in each equilibrium state. The numerical scheme for solving both problems is called the continuation algorithm. It is an algorithm of the predictor-corrector type. A continuation algorithm is used to solve the system of Equation (2) as a function of a single parameter of the system μ μ . In other words, the task comes down to determining the zeros of family f : n n , namely, to determine the solutions to stationary, non-linear algebraic Equation (2). The dimension of the µ parameter vector is called the bifurcation dimension. For a one dimensional case µµ∈ ℜ1.
The bifurcation theory concerning non-linear ordinary differential equations deals with a system of first-order differential Equation (1), which is a mathematical model of a dynamical system in an n-dimensional Euclidean space ℜn. If the system (1) has asymptotically stable stationary solutions x=0, then for all solutions x(0) belonging to this neighbourhood:
-
trajectory x(t) fulfils the condition: |x(t)| < ε for t > 0;
-
|x(t)| → 0 for t → ∞.
Solving the problem involves finding answers to the question of how a parameter change μμ locally affects the neighbourhood of point x = x0. Due to the fact that for all μ the Equation (2) is satisfied, and Equation (1) can be expressed as:
x ˙ = R x + σ ( x , μ )
where in R is a square characteristic matrix (Jacobian matrix) with elements provided by the equation:
[ R ]   i , j = f i ( x 0 , μ ) x j ,
and non-linear vector function σ fulfils the conditions:
σ ( x 0 , μ ) = 0 , σ i ( x 0 , μ ) x j = 0 .
and finally:
x ˙ = R x
The Hartman-Grobman theorem applies within the process of studying the stability of a stationary solutions to equation [15,16,19,20]. It states that if all eigenvalues of a Jacobian matrix R of a linearized system (1) lie in the left complex half-plane, i.e.,
Re ( λ j ) < 0 , i = 1 , 2 , , n
then there is a certain continuous, homomorphic transformation of variables reducing a locally non-linear system of Equation (1) to a linear system. This means that if a stationary solution to a linearized system of equations is asymptotically stable, then also the solution to a non-linear system is stable. The Hartman-Grobman theorem also indicates that every qualitative change in the nature of solutions to a system of non-linear equations describing a dynamical system is indicated by the appearance of zero real parts of the eigenvalues of a linearized system characteristic matrix R.
In mathematical terms, bifurcation of equilibrium positions takes place, when a eigenvalue of Jacobi matrix (3), (5) of system (1), estimated in the state of equilibrium, intersects the imaginary axis. Similarly, in the case of an oscillating solution, bifurcation occurs when the Floquet multiplier intersects the unit circle. The results shown in this article discuss five bifurcation types; they all have a codification of one, which means that they are encountered when a single continuation parameter changes.
Basic bifurcation types are [32]:
A.
A saddle-node bifurcation, also called a saddle limit point, occurs when the real eigenvalue of a Jacobian matrix (5) estimated in a state of equilibrium, intersects the imaginary axis. There is no equilibrium on one side of the bifurcation point (locally), whereas there are two equilibrium branches on the other (e.g., one stable and one unstable) (Figure 1).
B.
Hopf bifurcation occurs, when a complex pair of Jacobian (5) eigenvalues, assessed at equilibrium, intersect the imaginary axis. In this case, the equilibrium changes stability and a periodic orbit is formed, which can be stable or unstable (Figure 1).
C.
The limit point or periodic orbit fold bifurcation occur when the real Floquet multiplier intersects the unit circle at 1; as for the states of equilibrium, then there are no periodic orbits on one side of the bifurcation (locally), whereas there are two periodic orbits on the other (Figure 1).
D.
A period-doubling bifurcation occurs when the real Floquet multiplier intersects the unit circle at −1. The periodic orbit loses stability when a new period orbit appears with a period (approximately) twice as long (Figure 1).
E.
The Neimark-Sacker bifurcation or torus bifurcation appears, when the periodic orbit becomes unstable, namely, when a pair of complex Floquet multipliers intersects the unit circle and an additional oscillation frequency is introduced. The outcome is a torus dynamic, which can be periodic (blocked) or quasi-periodic (Figure 1).
To sum up, it can be concluded that:
The bifurcation theory is a set of mathematical results aimed at analysing and explaining unexpected modifications in the asymptotic behaviour of non-linear systems of differential equations, when their parameters slowly change.
Starting from the asymptotic state approximation, for a set value of bifurcation parameters, the computer code, in the course of a continuation process, determines curves for the solution of x(μ), which constitute a set of solutions to non-linear algebraic equations (2). In the case of each of the points of the solutions to system (2) and (7), the stability of solutions is determined:
{ Equilibrium   points :                                                                                      f ( x , μ ) = 0 Limit   points :                                                                                                 f ( x , μ ) = 0 and   at   least   one   of   eigenvalues :                                                                  λ = 0   Hopf   points   :                                                                                                 f ( x , μ ) = 0   and   at   least   Re ( λ i j ) = 0   of   pair   of   complex   eigenvluwes :             λ i j = ± 2 i π / T Periodic   orbits :                                                                                           x ( T ) = x ( 0 ) + 0 T f ( x , μ ) d t
The continuation process assumes that all functions of the system of Equation (2) are continuous and differentiable.
Summarising, the numerical continuation methods are a set of tools that enables obtaining information required to conduct bifurcation analyses. These methods utilize the predictor-corrector technique for finding, tracking and constructing equilibrium curves or periodic orbits of a differential Equation (1), as solutions of a properly defined system of algebraic Equation (2). Information on the stability is calculated based on analysing the J eigenvalues of Jacobian matrix (5) for equilibrium states, or based on Floquet multipliers (in the case of periodic orbits). Bifurcations can be detected, as well as tracked at various bifurcation parameter values.

2.2. Continuation Software

Numerical methods for solving bifurcation problems appeared relatively recently and complement the analytical achievements in this field. The following fundamental difficulties encountered when applying these methods can be distinguished:
-
numerical instabilities associated with calculations in close proximity to bifurcation points,
-
issues related to parametrization in close proximity to bifurcation points and limit points,
-
structures of bifurcation branches,
-
determination whether bifurcation actually takes place,
-
problems associated with the convergence of the Newton-Raphson method at singular points.
The system of Equation (2) can have numerous solutions. These can be isolated solutions and they may not exist at all. There is no theory, which could be used as a base to determine which case you are dealing with. Therefore, this issue is quite challenging, because when using numerical methods, the question whether all solutions have already been found still stands. Among the many methods, the ones applied the most when the f functions are smooth is the Newton-Raphson method.
At the time, there are several software packages available, which are designed to analyse non-linear dynamical systems, e.g., MatCont [33,34] or KRIT [35]. The most recognizable computer program designed for bifurcation analysis of non-linear dynamical systems is the one developed at the Canadian Concordia University, by a team led by Prof. Eusebius Doedel, the AOTO97 package—program developed in the FORTRAN language and AUTO2000—developed in the C language. The latest (as of 2021) version of these packages is the AUTO07P [36]. The description of subsequent versions of the AUTO package can be found in: [37,38,39]. AUTO07P was developed in the FORTRAN language, for the UNIX operating system environment. A team at the University of Pittsburgh in America, led by prof. Bard Ermentrout, developed XPPAUT [40], which is a version of the AUTO package compatible with the WINDOWS system. A comprehensive description of this package, together with a manual, can be found in the textbook by prof. Bard Ermentrout [41]. A MATLAB system toolbox—Dynamical System Toolbox [42] was also created based on the AUTO [43].
The best-known program designed for the bifurcation analysis of homogeneous ordinary differential equations is that developed by a team lead by Prof. Doedel from the Canadian Concordia University called AUTO. Differential equations describing MAV flight dynamics were analysed using an AUTO version implemented in MATLAB (Dynamical System Toolbox [42]).

3. Micro Aerial Vehicle Mathematical Model

3.1. Reference Frames

Coordinate systems rigidly associated with the object (MAV for example) or associated with the inflow of air streams, and their combinations are usually selected to represent motion equations in aviation. Figure 2 shows coordinate systems, used to derive aircraft motion equations. Details of deriving aircraft motion equations can be found in numerous textbooks and monographs (e.g., [44,45,46,47,48,49]), which is why will not discuss them in this work. Below you can find micro aerial vehicle motion equations expressed in a so-called-semi constrained coordinate system. This means that centre of mass motion equations are written within a velocity system of coordinates (Figure 2b), while the equations of MAV spherical motion relative to the centre of mass are written within a system related to MAV axes of inertia.

3.2. Equations of Motion

In the course of implementing a mathematical dynamics model of an aircraft into MATLAB Dynamical System Toolbox, it is more convenient to present aircraft centre of mass motion equations in a velocity coordinate system Oxayaza (Figure 2) [49]. Furthermore, the equations omitted the aircraft yaw (heading) angle Ψ and the aircraft centre of mass position relative to the system related to the Earth, since these values do not influence the dynamic properties of an aircraft [25]. This enabled the state vector dimension to be reduced to eight components:
x = [ V , α , β , P , Q , R , Θ , Φ ] T
where:
V  flight velocity,
α  angle of attack
β  slip angle,
P  banking angular velocity,
Q  pitching angular velocity,
R  yawing angular velocity,
Θ  pitch angle,
Ф  bank angle.
In the general case of bifurcation analysis of aircraft flight dynamics, the components of the vector of bifurcation parameters m may be the control surface deflection angles (aileron, elevator, rudder, thrust) and other parameters such as the position of the center of gravity. In this particular case of bifurcation analysis of strake-wing microdrone flight dynamics, the bifurcation parameters vector has following form:
μ = [ T , δ e , δ e l v ] T  
where:
T  propeller thrust
δe  angle of symmetrical elevon deflection
δelv  angle of asymmetrical elevon deflection,
The general form of the micro-airplane motion equations takes the form of autonomous differential Equation (1). Components of vector f have the form [44,45,46,47,48]:
f 1 = 1 m [ T cos α P X a ] g [ cos Θ sin Φ sin β ( sin Θ cos α cos Θ cos Φ sin α ) cos β ] f 2 = Q ( P cos α + R sin α ) tan β + 1 m V 0 cos β [ T sin α + P Z a m g ( sin Θ sin α + cos Θ cos Φ cos α ) ] f 3 = P sin α R cos α 1 m V [ T cos α m g ( sin Θ cos α cos Θ cos Φ sin α ) ] sin β + m g cos Θ cos Φ cos β P Y a }
f 4 = ( J X J Z J X J 2 X Z J X J Z ) Q R D + ( 1 J Y J X J Z ) J X Z P Q J X D + + q S b J X D { C l 0 ( α , β ) + C l R ¯ R ¯ b 2 V 0 + C l δ e l v δ e l v + [ C l P ¯ + ( 2 C l P ¯ β β + 2 C l P ¯ 2 P ¯ b 2 V 0 ) ] P ¯ b 2 V 0 } + L T J x D + + J X Z q S b J X J Z D { C n 0 ( α , β ) + C n R ¯ R ¯ b 2 V 0 + C n δ e l v δ e l v + [ C n P ¯ + ( 2 C n P ¯ β β + 2 C n P ¯ β P ¯ b 2 V 0 ) ] P ¯ b 2 V } + J x z N T J x J z D f 5 = q S c A J Y { C m 0 ( α , β ) + C m α ˙ 1 V 0 + C m Q ¯ Q ¯ c A 2 V 0 + C m δ e δ e } + M T J Y + J Z J X J Y Q P + J X Z J Y ( R 2 P 2 ) f 6 = ( J 2 X Z J X J Z J Y J X J Z ) P Q D + ( J Y J Z J X 1 ) J X Z Q R J Z D + + J X Z J X J Z { C l 0 ( α , β ) + C l R ¯ R ¯ b 2 V 0 + C l δ e l v δ e l v + [ C l P ¯ + ( 2 C l P ¯ β β + 2 C l P ¯ 2 P ¯ b 2 V 0 ) ] P ¯ b 2 V 0 } + J X Z L T J X J Z + + q S b J Z D { C n 0 ( α , β ) + C n R ¯ R ¯ b 2 V 0 + C n δ e l v δ e l v + [ C n P ¯ + ( 2 C n P ¯ β β + 2 C n P ¯ β P ¯ b 2 V 0 ) ] P ¯ b 2 V }
f 7 = P + Q sin Φ tan Θ + R cos Φ tan Θ f 8 = Q cos Φ R sin Φ
where:
D = 1 J 2 X Z J X J Z ,
q = 1 2 ρ V 0 2 is the dynamic pressure,
P ¯ = b P 2 V 0 ,   Q ¯ = c A Q 2 V 0 ,   R ¯ = b P 2 V 0 are dimensionless angular velocities.
The aerodynamic characteristics of the “Bee” MAV shown in Figure 3 were identified based on aerodynamic water tunnel testing. The static and dynamic aerodynamic loads were measured using a five-component aerodynamic balance. The testing was conducted over a wide range of angles of attack and slip. A wide description of these tests can be found in the works [49,50]. The identified aerodynamic characteristics presented in the work [Sibilski et al.] were shown in the form of graphs and were available in tabular form. Examples of aerodynamic derivative waveforms are shown in Figure 4.
The linear interpolation of aerodynamic data tables commonly used in simulation studies, due to the lack of derivative continuity, cannot be used in the case of continuation tests, since the continuation software of the AUTO type can misidentify bifurcation points. This is why a different method for interpolating aerodynamic data had to be used. Smooth, differentiable state parameter function were created as a result of interpolating aerodynamic characteristics. The pchip function interpolation in the MATLAB software was used for this purpose, while maintaining linear interpolation and a multi-variative orthogonal function. Block interpolation was also implemented using a 3rd-order spline function. Owing to the structure of the “Bee”, which only has elevons, the bifurcation parameters were the elevon deflection angles: δe and δelv.

4. Methodology of Bifurcation Tests in Aircraft Flight Dynamics

MAV motion is described through a system of highly non-linear ordinary differential equations. For a classical model of a non-deformable micro aerial vehicle with movable control surfaces, motion equations are presented by relationships (1), (12), (13), and (14).
The dynamical system theory enables analysing solutions to a system of highly non-linear ordinary differential equations describing aircraft motion, depending on slow parameter changes (so-called bifurcation parameters). When analysing aircraft flight dynamics, it is assumed that the bifurcation parameters are the control vector components (i.e., deflection angles of the elevator, rudder, ailerons, thrust vector, etc.). The first stage in the analysis of a non-linear system of differential equations in the dynamical system theory is assessing the stability of steady states of a system of differential equations describing aircraft flight dynamics (1), (12), (13), and (14). The steady state is determined by equating the derivatives to zero and solving a system of algebraic Equation (2).
Given the experience from numerous research (based on the bifurcation analysis and continuation technique), the following, three-stage test diagram for a non-linear aircraft motion was formulated [8,9,10,11,12,13,14,18,21,22,23,24,25,26,27,28,34,48,49,51,52,53,54,55,56]
  • The first stage involves determining all parameters of a dynamical system. The fundamental task is to study all possible equilibrium states and periodic orbits, and the analysis of their local stability. This test should be very thorough. The outcome of the attempted test should be a determined global structure of the state space (e.g., phase portraits) of all discovered attractors (steady states and closed orbits). Approximated graphic representations of the calculations are crucial in this case, since they enable diagnosing the obtained results.
  • The second stage involves, based on information on the evolution of phase portraits together with parameter changes, predicting dynamical system behaviour. Next, based on the knowledge on the type of present bifurcations and the current position of system parameters relative to stable areas, further aircraft behaviour is predicted. Information on the range of parameter changes is also important for these analyses and predictions. Rapid parameter changes and higher differences between steady and transient states are also observed.
  • The third step involves a numerical simulation, which enables verification of the expected aircraft behaviour. Waveforms of transient system characteristics for significant state parameter changes upon a dynamical system parameter change are obtained.

5. Bifurcation Flight Dynamic Analysis of a Micro Aerial Vehicle

Figure 5 and Figure 6 show single-parameter bifurcation diagrams for equilibrium positions: angles of attack α and angular velocity of banking P for different values of the bifurcation parameter, namely, elevator deflection angles δe from a range of δe ∈ (−35°, 30°) for angles of attack equilibrium position, δe ∈ (−35°, 40°) for banking angular velocity equilibrium positions. Figure 5 shows enlarged bifurcation diagrams for elevator deflections in the range of δe ∈ (3°, 10°), corresponding to steady states at low and moderate angles of attack for α ∈ (0°, 54°), the second area corresponds to a range of high angles of attack α ∈ (65°, 83°). Various types of micro aerial vehicle dynamics, corresponding to seven flight regimes, can be classified in these areas. These regimes are marked with the letters A to G. Regime A means steady symmetrical flight equilibrium states, for angles of attack in the range of α ∈ (2°, 23°). Area B, for angles of attack in the range of α ∈ (15°, 50°), corresponds to various motion states with stable and unstable orbits. Area C, for angles of attack in the range of α ∈ (35°, 55°), corresponds to deep steady spirals. Area D corresponds to inverted spirals. Area E corresponds to steady spins. Areas F and G correspond to transient spins.
The aircraft equilibrium positions in each of these two areas were analysed depending on the bifurcation parameter, namely elevator deflection angles (symmetrical elevon deflections) δe and asymmetrical elevon deflections δelv. During continuation tests, following the methodology described in Section 4, it was assumed that the propeller thrust does not change, and elevon deflections are only symmetrical (or δelv = 0). Continuation analyses were conducted for positive and negative changes in the bifurcation parameter (symmetrical elevon deflections δe), starting with δe = 20°. The calculations were conducted for the elevator deflection angle range of change of −35° ≤ δe ≤ 40°, so that it was possible to detect almost all solutions corresponding to the quasi-steady flight states.
At the starting point of the continuation analysis (δe = 20°, α = 3°), the aircraft was conducting a steady symmetrical flight. As the elevator deflection angle changed, the elevator angle of attack initially increased, and the deflection angular velocity remained at zero. This dynamic regime was denoted with the letter A. Its corresponding range of angles of attack is 2° < α < 25° (Figure 5 and Figure 7). The range of steady states denoted with the letter A corresponds to a symmetrical flight of the MAV (angular velocity P = 0). When the elevator deflection angle decreases below δe≈10 and the angle of attack α ≈ 28 degrees, the micro aerial vehicle enters the dynamics regime denoted with the letter B, with both anti-symmetrical, as well as symmetrical motions. They can be associated with unstable Dutch roll, wing-rock oscillations, spiral motion instabilities, as well as unstable phugoids. Ant-symmetrical oscillations result from Hopf bifurcations in B. In this area, when the elevator deflection angle continues to decrease, reaching a value of δe ≈ 4° and α ≈ 34°, the spiral mode loses stability, which leads to the bifurcation of the stable and almost symmetrical solution. Two asymmetrical equilibrium position branches appear. At higher angles of attack, the spiral model becomes unstable along the almost symmetrical solution branch.
Development of equilibrium position asymmetry can be observed at continued reduction in the elevator deflection angle, and at δe ≈ 4°, the micro aerial vehicle enters the range of equilibrium states marked with C, which corresponds to spiral motions with a high amplitude of bank angles (Figure 7). Stable position branches corresponding to spiral motions exist on both sides of the straight line defining unstable horizontal flight equilibrium positions. Due to the aerodynamic load asymmetry, the micro aerial vehicle banks to the left wing (with a negative banking angular velocity P), since the equilibrium branch representing the breaking away is connected with A (Figure 6). The equilibrium branch representing rightward MAV banking is practically detached from branch A. The bifurcations therein lead to the appearance of unstable equilibrium positions, where P ≈ 0 [°/s]. However, it should be noted that the spiral motion direction is determined in practice by the nature of the disturbance or the transient dynamics state of a micro aerial vehicle. As shown in Figure 5 and Figure 6, equilibrium branches of dynamics regime C (within the physical range of elevator deflections), do not converge on the straight line corresponding to the horizontal flight conditions. The MAV will remain at one of two spiral branches only up to a value of δe = −35°. Branch C, shown on the bifurcation diagrams (Figure 5 and Figure 6), also represents undesirable steep spirals. Steep spirals can be deemed undesirable flight conditions, therefore, it is important to determine a recovery strategy. Based on a bifurcation diagram (Figure 5), a micro aerial vehicle can be easily recovered from such a flight state by reducing the angle of attack to a level below which no spiral branches are formed. This can be achieved by increasing the elevator deflection angle. The ability to find and determine branches, such as steep spirals, stresses the advantage of continuation and bifurcation analyses over classical linear methods. Although classical linear methods for analysing flight stability can identify stability changes, they are unable to find stable and unstable branches of quasi-steady flight states (Figure 8).
The slight lateral instability is also present at low angles of attack (α ≤ 2°), (branch D in Figure 5 and Figure 6). In the case of the bifurcation analysis and continuation test pattern in question, which assumes constant thrust, flight velocity reaches the highest and not always realistic values at low, negative angles of attack, which entails an increase in the negative angle of MAV pitch, until inverted flight conditions. Thus, although regime D can be deemed an inverted spiral, it is not representative for realistic flight conditions and will not be discussed in more detail (Figure 8).
Equilibrium positions of quasi-steady flight states also appear at high angles of attack 65° < α < 85° (bifurcation diagrams shown in Figure 5 and Figure 6). There are unstable equilibrium position branches for negative deflection angular velocities and angles of attack within a range of 80° < α < 85°, and elevator deflections angles in the range of −35° < δe < 10°. There is a certain area of stable equilibrium positions for positive values of deflection angular velocity. The differences arise from the asymmetry of MAV loads. Left- and right-handed spins differ in terms of the state of equilibrium of inertial and aerodynamic forces, which undoubtedly impacts the local stability of spin branches (just as in the case of lower angles of attack α). Stable equilibrium positions, corresponding to spins executed at a constant angular velocity, are represented by branch E of the bifurcation diagrams. It can be seen that there is also a slight area of stable spins near the elevator deflection angle of δ e 12 o . The waveform of flight parameters indicates that attractors are present within branch E, hence, this equilibrium branch does not play a significant role in the dynamics of disturbed spins [26,27].
For an angle of attack range of 65° < α < 85°, most equilibrium positions are unstable. There can be no steady spins on unstable branches, while one can expect solutions of deterministic chaos nature on these branches [13]. The spin dynamics is significantly impacted by two Hopf bifurcations present at δe = 23.4° and δe = −13.2° and a period-doubling bifurcation present at δe = 12.5°. All three bifurcations appear for angles of attack of α > 65°. In the case of an elevator deflection angle of δe = −13.2°, the periodic orbits created through the Hopf bifurcation are unstable throughout the entire range of bifurcation parameters. The period-doubling bifurcation appearing at δe = 12.5° can lead to a spin of chaotic nature [13]. The periodic orbit created through Hopf bifurcation at δe = 23.4° is more significant. It is initially unstable and then branches into a stable periodic orbit through the torus bifurcation at δe = −15.9°. The stable periodic orbit of branch F (Figure 5) is significant, since in this case, the average angle of attack is α ≈ 70°, with high values of deflection angular velocities. Therefore, it can be concluded that branch F corresponds to steep oscillation spin states. It can also be inferred that MAV oscillation spins occur at angles of attack of approximately 70°. Also, a second stable periodic orbit can be identified on branch G (Figure 6) within the range of high angles of attack. This solution is true for a closed curve that is not connected with other curves. It is the outcome of torus bifurcation present on branch F (for δe = −15.9°). Due to torus bifurcations, the periodic orbit becomes unstable, and MAV flight dynamics equation solutions are attracted by branch G periodic orbit. Solutions to micro aerial vehicle motion equations over branch G represent more rapid, oscillation spins, relative to branch F solutions [26,27].

6. Numerical Verification of Predicted MAV Behaviour

According to the bifurcation test methodology, the third step in continuation analysis involves numerical simulations, which enable verifying the predicted aircraft behaviour. Waveforms of transient system characteristics for significant state parameter changes upon a dynamical system parameter change are obtained. Examples of “wing rock” oscillation, unsteady spin and the “Cobra” manoeuvre simulations are shown below.

6.1. “Wing Rock” Oscillation Simulation

A typical phenomenon encountered when flying at large angles of attack is the so-called “wing rock”. It is a self-excited phenomenon, which occurs at subcritical and supercritical angles of attack. Despite “wing-rock” generally not being dangerous, it is advisable to examine it more thoroughly. Studying the dynamics of this motion was attempted in numerous research work (for example: [10,55,56,57,58,59,60]. All of the quoted work adopted quasi-static aerodynamic aircraft characteristics or presented different identification methods. However, due to the oscillatory changes in the angle of attack occurring at high, near-critical and supercritical angles of attack, the occurrence of phenomena associated with flow non-stationarity, including hysteresis, is obvious. “Wing rock” oscillation simulations were conducted based on aerodynamic derivatives obtained during water tunnel testing and identified through a pulse function.
Examples of simulation results are shown in Figure 9. The motion disturbance involved elevator displacement from 8° to 7.8°. This caused MAV displacement into an area of unstable steady states (Hopf bifurcation). The outcomes were gradually increasing oscillations, turning into a limit cycle, with a period of about 0.2 s.
Characteristic oscillations appeared as a result of the micro aerial vehicle moving into the area of unstable steady states. They primarily involve rocking of the MAV from wing to wing. That rocking has an amplitude of approximately 39° (for the banking angle), and 6° for the angle of attack and slip. The flight velocity is practically constant. Whereas the period of oscillations corresponding to MAV pitch angle change is equal to the period of phugoid motions. These oscillations are of limit cycle vibration character.

6.2. MAV Spin Simulation

The G and F branches of quasi-steady states (Figure 6 and Figure 7) correspond to states of steady spins. The bifurcation analysis enables “control matching”, which allows to recover from critical flight states. A sketch of a “bifurcation control matching” is shown in Figure 10. The equilibrium surface of banking angular velocities P and the bifurcation graph on the deflection plane of the elevator δe and elevons δelv show critical autorotation areas (excluding the spin and a rapid inertia barrel). Recovery trajectories encounter an “apex catastrophe”. Skipping through this singularity allows to “smoothly” achieve the desired point of zero banking rate.
Figure 11 show the results of a spin recovery simulation. The method of spin recovery was matched using the method shown in Figure 10 and involved displacement of control surface positions beyond the area of unstable equilibrium states and limit cycles. The analysis of simulation results shows the effectiveness of the “bifurcation control matching” method.

6.3. “Cobra” Manoeuvre Simulation

“Cobra” is one of aerobatic figures. The manoeuvre was first executed at the turn of the 1950s and 60s, by the pilots of the Swedish Air Force (Svenska Flygvapnet) flying J-35 Draken fighters. It has a Swedish name of “kort parad”—“short parade” and was part of standard short-range manoeuvring combat training for Swedish fighter pilots [48,49]. In the 1980s, the manoeuvre was executed by OKB Sukhoi test pilot, Igor Volk, during Su-27 spin tests. The Cobra was first demonstrated in public by Viktor Pugachev at the Le Bourget air show in 1989, on a MiG-29 fighter. This manoeuvre can be executed with an aircraft exhibiting excellent manoeuvrability and low thrust load or equipped with thrust vector control engines. Besides the spectacular impression it makes at air shows, enabling to demonstrate the manoeuvrability and turning abilities of modern combat aircraft, this manoeuvre, used in direct air combat, is primarily aimed at forcing the chasing foe to overtake through rapid deceleration, thus providing the chased aircraft with a convenient position to open fire. The effect of sudden deceleration is achieved owing to rapidly increasing aircraft drag resulting from vertical positioning of the aircraft body in the upward direction (perpendicular to the previous flight direction).
The “Cobra” is conducted at supercritical angles of attack, significantly exceeding the value under normal operating conditions. Figure 12 shows a diagram of the “Cobra”. This manoeuvre is divided into three basic stages:
  • Transition from horizontal flight into the phase of increasing the aircraft pitch angle, due to very rapid increase in the elevator deflection angle (sudden pulling of the control stick to a maximum), while simultaneously throttling the engine or engines,
  • manoeuvre phase in which, as a result of such action by the pilot, the aircraft nose rapidly rises up, until it reaches a very high angle of attack (even up to 120°),
  • the exit phase, which involves increasing the thrust and releasing the control stick, leading to the aircraft rapidly increasing its pitch angle, while simultaneously accelerating and returning to horizontal flight, with a minor altitude loss.
The “Bee” MAV has a stake-wing and its dynamics is similar to the dynamics of modern high-manoeuvring combat aircraft. Based on in-flight tests, it was concluded that it was able to execute a “Cobra” manoeuvre (Figure 12) [61]. The “Cobra” manoeuvre simulations were conducted based on aerodynamic and mass data of the “Bee” micro aerial vehicle. Simple analyses show that executing the “Cobra” manoeuvre will be possible without changing the flight altitude, when the vertical projection of the sum of external forces acting on the MAV in the course of the manoeuvre should be equal to zero. However, due to the fact that the computational thrust force value was initially negative, it was assumed that the MAV starts the manoeuvre with the engine off (zero thrust). It was also assumed that, in the course of executing a manoeuvre, the thrust depends on the difference in the aircraft’s angles of pitch Θ and attack α, and the flight velocity V [62,63].
Due to the fact that aircraft equipped with strake wings are characterized by wing-rock instability and are not spirally stable, it was impossible to obtain a fully symmetrical flight parameter waveform. This is associated with the occurrence of Hopf bifurcation and torus creation bifurcation on G, E and F branches (Figure 5 and Figure 6). The “Cobra” manoeuvre is initiated with a sudden downward displacement of elevons (δe = −35°). This causes a rapid transition of the MAV through C, G, E and F branches of steady flight states (Figure 5 and Figure 6).
Figure 13 and Figure 14 show the results of a digital simulation of the “Cobra” manoeuvre, taking into account the occurrence of limit cycles. Based on the analysis of these graphs, it can be concluded that all flight parameters are significantly changed during a Cobra figure, increasing their values. Using the terminology of the Dynamical System Theory, it can be said that the “Cobra” is unstable due to the presence of Hopf bifurcation (at an elevator displacement angle of δe = −23.1°) and torus bifurcation (for δe = −16.5°). Figure 14 shows Poincarè maps for selected state parameters. It can be concluded when the non-stationarity (hysteresis) of aerodynamic coefficients were taken into account, MAV motion equation solution irregularities of quasi-harmonic nature were obtained. The digital simulation took into account the fact that at high angles of attack, the aircraft is spirally unstable (branch C, Figure 5 and Figure 6) and has a tendency to wing-rock oscillations (branch B, Figure 6 and Figure 7). In the course of a manoeuvre, the aircraft attack and pitch angles rapidly increase, reaching maximum values of approximately 84° (for an attack angle of α) and 100° (for a pitch angle of Θ). The MAV bank and yaw angle waveforms (Figure 13) indicate that these angles increase over time. This is associated with the present area of spiral instability (branch C in bifurcation diagrams—Figure 5 and Figure 6). The development of wing-rock oscillations is also visible. The period of these oscillations varies at approx. 0.2 s. It should be noted that the amplitude and frequency of these oscillations are irregular (quasi-harmonic).

7. Conclusions

The article presents application examples of the theory of dynamical systems relative to studying the flight dynamic specifics of a micro aerial vehicle, constructed as a fixed-wing aircraft with a strake-wing. Such microdrons are characterized by high manoeuvrability and can fly with high, supercritical angles of attack. The bifurcation analysis shown in the article enabled identifying some of the numerous factors impacting the behaviour of a strake-wing MAV. More specifically, this analysis allowed to discover a number of stable attractors, associated with disturbances to the states of equilibrium. Compared with purely simulational studies, which require long-term computations leading to the disappearance of motion history transitions on weak attractors, creating bifurcation diagrams is much less computation-demanding. Furthermore, bifurcation diagrams show the fundamental structure of a dynamical system, hence they suggest, where and when the time domain simulation should be conducted in order to better explain the behaviour of a dynamical system (in this case, a micro aerial vehicle). Time domain simulations were presented as supplementary to bifurcation diagrams, in order to gain more clarity in terms of the nature of various attractor dynamics regimes. It was shown that an MAV was susceptible to steep spiral motion disturbances, induced by loss of stability in this flight range. Bifurcation analysis was also able to identify that elevator displacement aimed at recovery from a steep spiral glide can lead to an oscillating spin. Bifurcation diagram analysis indicates that the correct reaction is to restore elevator position to the value corresponding to balancing conditions in straight flight.
The approach adopted in this article can be extended in several ways. Firstly, MAV behaviour in an open control loop can be further tested through expanding the parameter space in question; this can cover additional control signals, as well as combining the centre of gravity and various damage scenarios, which can be added to the model. Another potentially beneficial extension is the use of bifurcation methods to analyse micro aerial vehicle models expanded with closed control loops, including bifurcation matching that enables recovery from critical flight states.

Author Contributions

Conceptualization, K.S. and A.Ż.; methodology, M.N.; software, A.S.-M.; validation, K.S. and M.N.; formal analysis, K.S., and A.Ż.; writing—original draft preparation, M.N. and K.S.; writing—review and editing, K.S.; visualization, A.S-M.; supervision, M.N., and A.Ż.; project administration, K.S.; funding acquisition, M.N., and A.Ż. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially supported by Air Force Institute of Technology by subsidy from financial resources for the maintenance and development of research potential, and the statutory funds of the Faculty Mechatronics of the Warsaw University of Technology granted for 2020/2021.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

Clbody axis rolling moment coefficient
Clprolling moment coefficient derivative with respect to to rolling rate
Cmbody axis pitching moment coefficient
Cpitching moment coefficient derivative with respect to angle of attack
Cmqpitching moment coefficient derivative with respect to pitching angular rate
Cnbody axis yawing moment coefficient
Cnryawing moment coefficient derivative with respect to to yawing rate
freduced frequency of model oscillation in water tunnel
fvector of generic nonlinear function
fi, i=1,…8components of f vector describing microdrone flight dynamics
gacceleration of gravity
JJacobi matrix
JX, JY, JZ, JXZmoments of inertia of microdrone
LTbody axis banking moment due to propulsion
mmass of microdrone
MTbody axis pitching moment due to propulsion
NTbody axis yawing moment due to propulsion
V, V0flight velocity
Pbody axia roll (banking) rate
PXax wind axis aerodynamic force component
PYay wind axis aerodynamic force component
PZaz wind axis aerodynamic force component
qdynamic pressure
Qbody axis pitching rate
Rbody axis yawing rate
Rstate matrix (Jacobian of linearised aircraft equations of motion)
Swing area
Tpropeller thrust
αangle of attack
βslip angle
δeangle of symmetrical elevon deflection
δelvangle of asymmetrical elevon deflection
Θpitch angle
λeigenvalue
μvector of bifurcation parameters (in this case microdrone control vector)
μsingle bifurcation parameter
ρair density
Фbank (roll) angle
Ψyaw (heading) angle
( ) . = d d t time derivative
( ) _ dimensionless quantity
MAVmicro aerial vehicle, microdrone, micro aircraft

References

  1. Abdulrahim, M.; Watkins, S.; Segal, R.; Sheridan, J. Dynamic sensitivity to atmospheric turbulence of Fixed-Wing UAV with varying configuration. J. Aircr. 2010, 47, 1873–1883. [Google Scholar] [CrossRef]
  2. Wróblewski, W.; Sibilski, K.; Garbowski, M.; Żyluk, A. The gust resistant MAV—Aerodynamic measurements, performance analysis, and flight tests (AIAA2015—1684 CP). In Proceedings of the AIAA SciTech Forum and AIAA Atmospheric Flight Mechanics Conference, Kissimmee, FL, USA, 5–9 January 2015. [Google Scholar] [CrossRef]
  3. Kowalski, M.; Sibilski, K.; Żyluk, A. Studies and tests of micro aerial vehicle during flight. J. KONES 2015, 22, 155–162. [Google Scholar] [CrossRef]
  4. Adams, W.W. SPINEQ: A Program for Determining Aircraft Equilibrium Spin Characteristics Including Stability, 1979, NASA TM 78759; The National Aeronautics and Space Administration: Washington, DC, USA, 1978. Available online: https://ntrs.nasa.gov/api/citations/19790002903/downloads/19790002903.pdf (accessed on 5 January 2021).
  5. Schy, A.A.; Hannaht, M.F. Prediction of Jump Phenomena in Roll-Coupled Maneuvers of Airplanes. J. Aircr. 1977, 14, 375–382. [Google Scholar] [CrossRef]
  6. Roberts, L.; Hamel, P.; Orlik-Ruckeman, K.J. (Eds.) AGARD CP-386 Unsteady Aerodynamics-Fundamental and Applications to Aircraft Dynamics; North Atlantic Treaty Organization Advisory Group for Aerospace Research and Development: Neuilly-sur-Seine, France, 1985; ISBN 92-835-0382-1. Available online: https://apps.dtic.mil/dtic/tr/fulltext/u2/a165045.pdf (accessed on 3 January 2021).
  7. Tobak, M.; Schiff, L.B. On the Formulation of the Aerodynamic Characteristics in Aircraft Dynamics, NASA TR-R-456; The National Aeronautics and Space Administration: Washington, DC, USA, 1976. Available online: https://ntrs.nasa.gov/api/citations/19760007994/downloads/19760007994.pdf (accessed on 5 January 2021).
  8. Tobak, M.; Schiff, L.B. The Role of Time-History Effects in the Formulation of the Aerodynamics of Aircraft Dynamics, NASA TM 78471; The National Aeronautics and Space Administration: Washington, DC, USA, 1978. Available online: https://ntrs.nasa.gov/citations/19780011113 (accessed on 5 January 2021).
  9. Hui, W.; Tobak, M. Bifurcation analysis of nonlinear stability of aircraft at high angles of attack, (AIAA 82-244 CP). In Proceedings of the 20th AIAA Aerospace Sciences Meeting, Orlando, FL, USA, 11–14 January 1982. [Google Scholar] [CrossRef]
  10. Carroll, J.V.; Mehra, R.K. Bifurcation analysis of non-linear aircraft dynamics. J. Guid. Contro Dyn. 1982, 5, 529–536. [Google Scholar] [CrossRef]
  11. Tobak, M.; Ünal, A. Bifurcation in Unsteady Aerodynamics, NASA TM 8316; The National Aeronautics and Space Administration: Washington, DC, USA, 1986. Available online: https://ntrs.nasa.gov/api/citations/19870002264/downloads/19870002264.pdf (accessed on 5 January 2021).
  12. Guicheteau, P. Bifurcation theory in flight dynamics and application to a real combat aircraft (ICAS-90-5.10.4 CP). In Proceedings of the 17th ICAS Congress, Stockholm, Sweden, 9–17 September 1990; Available online: https://www.icas.org/ICAS_ARCHIVE/ICAS1990/ICAS-90-5.10.4.pdf (accessed on 5 January 2021).
  13. Guicheteau, P. Bifurcation theory: A tool for nonlinear flight dynamics. Phil. Trans. R. Soc. Lond. A 1998, 356, 2181–2201. [Google Scholar] [CrossRef]
  14. Mehra, R.; Prasanth, R. Bifurcation and limit cycle analysis of nonlinear pilot induced oscillations (AIAA 98-4249 CP). In Proceedings of the 23rd AIAA Atmospheric Flight Mechanics Conference, Boston, MA, USA, 10–12 August 1998. [Google Scholar] [CrossRef]
  15. Wiggins, S. Introduction to Applied Nonlinear Dynamical Systems and Chaos, 2nd ed.; Springer: New York, NY, USA, 2003; ISBN 0-387-00177-8. Available online: https://www.springer.com/gp/book/9780387001777 (accessed on 5 January 2021).
  16. Tobak, M. On the Use of the Indicial Function Concept in the Analysis of Unsteady Motions of Wings and Wing-Tail Combinations; NACA Report 1188; The National Aeronautics and Space Administration: Washington, DC, USA, 1954. Available online: https://digital.library.unt.edu/ark:/67531/metadc65696/ (accessed on 5 January 2021).
  17. Kuznetsov, Y.A. Elements of Applied Bifurcation Theory; Springer: New York, NY, USA, 1998; ISBN 0-387-98382-1. Available online: https://www.springer.com/gp/book/9780387219066 (accessed on 5 January 2021).
  18. Ioos, G.; Joseph, D. Elementary Stability and Bifurcation Theory, 2nd ed.; Springer: New York, NY, USA, 2002; ISBN 978-1-4612-7020-2. Available online: https://link.springer.com/book/10.1007/978-1-4612-1140-2 (accessed on 5 January 2021).
  19. Keller, H.B.; Langford, W.F. Iterations, perturbations and multiplicities for nonlinear bifurcation problems. Arch. Ration. Mech. Anal. 1972, 48, 83–108. [Google Scholar] [CrossRef]
  20. Keller, H.B. Lecture Notes on Numerical Methods in Bifurcation Problems; Springer: New York, NY, USA, 1987; ISBN 3-540-20228-5. [Google Scholar]
  21. Abed, E.H.; Lee, H. Nonlinear Stabilization of High Angle-of-Attack Flight Dynamics using Bifurcation Control. In Proceedings of the 1990 American Control Conference, San Diego, CA, USA, 23–25 May 1990; pp. 2235–2238. [Google Scholar] [CrossRef] [Green Version]
  22. Avanzini, G.; de Matteis, G. Bifurcation analysis of a highly augmented aircraft model. J. Guid. Control Dyn. 1997, 20, 754–759. [Google Scholar] [CrossRef]
  23. Charles, G.; Lowenberg, M.; Stoten, D.; Wang, X.; di Bernardo, M. Aircraft Flight Dynamics Analysis and Controller Design Using Bifurcation Tailoring (AIAA-2002-4751 CP). In Proceedings of the AIAA Guidance, Navigation, and Control Conference and Exhibit, Monterey, CA, USA, 5–8 August 2002. [Google Scholar] [CrossRef]
  24. Goman, M.G.; Khramtsovsky, A.V. Application of continuation and bifurcation methods to the design of control systems. Philos. Trans. R. Soc. Lond. A 1998, 356, 2277–2294. [Google Scholar] [CrossRef]
  25. Goman, M.G.; Zagainov, G.I.; Khramtsovsky, A.V. Application of bifurcation methods to nonlinear flight dynamics problems. Prog. Aerosp. Sci. 1997, 33, 9–10, 539–586. [Google Scholar] [CrossRef]
  26. Gill, S.J.; Lowenberg, M.H.; Neild, S.A.; Krauskopf, B.; Puyou, G.; Coetzee, E. Upset Dynamics of an Airliner Model: A Nonlinear Bifurcation Analysis. J. Aircr. 2013, 50, 1832–1842. [Google Scholar] [CrossRef] [Green Version]
  27. Gill, S.J.; Lowenberg, M.H.; Neild, S.A.; Crespo, L.G.; Krauskopf, B.; Puyou, G. Nonlinear Dynamics of Aircraft Controller Characteristics Outside the Standard Flight Envelope, J. Guid. Control Dyn. 2015, 38, 2301–2308. [Google Scholar] [CrossRef] [Green Version]
  28. Eaton, A.J.; Howcroft, C.; Coetzee, E.B.; Neild, S.A.; Lowenberg, M.H.; Cooper, J.E. Numerical Continuation of Limit Cycle Oscillations and Bifurcations in High-Aspect-Ratio Wings. Aerospace 2018, 5, 78. [Google Scholar] [CrossRef] [Green Version]
  29. Angiulli, G.; Calcagno, S.; De Carlo, D.; Laganá, F.; Versaci, M. Second-Order Parabolic Equation to Model, Analyze, and Forecast Thermal-Stress Distribution in Aircraft Plate Attack Wing–Fuselage. Mathematics 2020, 8, 6. [Google Scholar] [CrossRef] [Green Version]
  30. Awrejcewicz, J.; Pyryev, J.; Kudra, G.; Olejnik, P. Mathematical and Numerical Methods of Bifurcation and Chaotic Dynamics Analysis of Mechanical Systems with Friction and Impact; Publishing House of the Lodz University of Technology: Łódź, Poland, 2006; ISBN 83-7283-173-4. (In Polish) [Google Scholar]
  31. Guckenheimer, J.; Holmes, P. Non-Linear Oscillators, Dynamical Systems, and Bifurcations of Vector Fields; Springer: New York, NY, USA, 2002; ISBN 978-1-4612-7020-1. Available online: https://link.springer.com/book/10.1007/978-1-4612-1140-2 (accessed on 5 January 2021).
  32. Magnitskii, N.A.; Sidorov, S.V. New Methods for Chaotic Dynamics; World Scientific Series on Nonlinear Science Series A: Singapore, 2006; Volume 58, ISBN 978-981-256-817-5. [Google Scholar] [CrossRef]
  33. Available online: https://sourceforge.net/projects/matcont/files/matcont/ (accessed on 5 January 2021).
  34. Dhooge, A.W.; Govaerts, W.; Kuznetsov Yu, A. Matconta: Matlab Package for Numerical Bifurcation Analysis of ODEs. ACM Trans. Math. Softw. 2003, 29, 141–164. [Google Scholar] [CrossRef]
  35. Goman, M.G.; Khramtsovsky, A.V. Computational framework for investigation of aircraft nonlinear dynamics. Adv. Eng. Softw. 2008, 39, 167–177. [Google Scholar] [CrossRef] [Green Version]
  36. Available online: http://indy.cs.concordia.ca/auto/ (accessed on 5 January 2021).
  37. Doedel, E.J.; Fairgrieve, T.F.; Champneys, A.R.; Sandstede, B.; Kuznetsov, Y.A.; Wang, X. Auto97: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont); Technical Report for Concordia University: Montreal, QC, Canada, 1998; Available online: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.44.9955 (accessed on 5 January 2021).
  38. Doedel, E.; Keller, H.B.; Kernevez, J.P. Numerical analysis and control of bifurcation problems. Int. J. Bifurc. Chaos 1991, 1, 493–520. [Google Scholar] [CrossRef]
  39. Doedel, E.J.; Oldeman, B.E. AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Equations; Technical Report for Concordia University: Montreal, QC, Canada, 2009; Available online: https://www.macs.hw.ac.uk/~gabriel/auto07/auto.html (accessed on 5 January 2021).
  40. Available online: http://www.math.pitt.edu/~bard/xpp/xpp.html (accessed on 5 January 2021).
  41. Ermentrout, B. Simulating, Analyzing, and Animating Dynamical Systems. A Guide to XPPAUT for Researchers and Students; SIAM: Philadephia, PA, USA, 2002; ISBN 978-0-89871-819-5. Available online: https://epubs.siam.org/doi/abs/10.1137/1.9780898718195 (accessed on 5 January 2021).
  42. Available online: https://www.mathworks.com/matlabcentral/fileexchange/32210-dynamical-systems-toolbox (accessed on 5 January 2021).
  43. Coetzee, E.B.; Krauskopf, B.; Lowenberg, M.H. The Dynamical Systems Toolbox: Integrating AUTO into MATLAB. In Proceedings of the 16th U.S. National Congress of Theoretical and Applied Mechanics, State College, PA, USA, 27 June–2 July 2010. [Google Scholar]
  44. Etkin, B.; Reid, L.D. Dynamics of Atmospheric Flight, 3rd ed.; John Willey & Sons Inc.: Hoboken, NJ, USA, 1996; ISBN 0-471-03418-5. [Google Scholar]
  45. Pamadi, B.N. Performance, Stability, Dynamics, and Control of Airplanes, 2nd ed.; AIAA: Reston, VA, USA, 2003; ISBN 978-1-62410-274-5. [Google Scholar] [CrossRef]
  46. Schmidt, L. Introduction to Flight Dynamics; AIAA, Ed.; Series AIAA: Reston, VA, USA, 1998; ISBN 1-56347-226-0. [Google Scholar] [CrossRef]
  47. Zipfel, P.H. Modeling and Simulation of Aerospace Vehicle Dynamics; AIAA, Ed.; Series AIAA: Reston, VA, USA, 2003; ISBN 978-1-62410-250-9. [Google Scholar] [CrossRef]
  48. Sibilski, K. Modeling and Simulation of Flying Vehicles Dynamics; MH Publishing House: Warsaw, Poland, 2004; ISBN 83-906620-8-6. [Google Scholar]
  49. Sibilski, K.; Lasek, M.; Sibilska-Mroziewicz, A.; Garbowski, M. Dynamics of Flight of Fixed Wings Micro Aerial Vehicle; Warsaw University of Technology Publishing House: Warsaw, Poland, 2020; ISBN 978-83-8156-124-2. [Google Scholar]
  50. Sibilski, K.; Nowakowski, M.; Rykaczewski, D.; Szczepaniak, P.; Żyluk, A.; Sibilska-Mroziewicz, A.; Garbowski, M.; Wróblewski, W. Identification of Fixed-Wing Micro Aerial Vehicle Aerodynamic Derivatives from Dynamic Water Tunnel Tests. Aerospace 2020, 7, 116. [Google Scholar] [CrossRef]
  51. Abramov, N.; Goman, M.; Khrabrov, A. Aircraft dynamics at high incidence flight with account of unsteady aerodynamic effects. In Proceedings of the AIAA Meeting Papers, AIAA Aymospheric Flight Mrchanics Conference and Exhibit, AIAA 2004-5274 CP, Providence, RI, USA, 16–19 August 2004. [Google Scholar] [CrossRef]
  52. Abramov, A.N.; Goman, M.; Khrabrov, A.; Kolesnikov, E.; Fucke, L.; Soemarwoto, B.; Smaili, H. Pushing Ahead—SUPRA Airplane Model for Upset Recovery. In Proceedings of the AIAA Modeling and Simulation Technologies Conference, (AIAA 2012-4631 CP), Minneapolis, MN, USA, 13–16 August 2012. [Google Scholar] [CrossRef]
  53. Pauck, S.; Jacobus Engelbrecht, J. Bifurcation Analysis of the Generic Transport Model with a view to Upset Recovery. In Proceedings of the AIAA Meeting Papers, AIAA Atmospheric Flight Mechanics Conference, AIAA 2012-4646, Minneapolis, MN, USA, 13–16 August 2012. [Google Scholar] [CrossRef]
  54. Cunis, T.; Condomines, J.-P.; Burlion, L.; la Cour-Harbo, A. Dynamic Stability Analysis of Aircraft Flight in Deep Stall. J. Aircr. 2020, 57, 143–155. [Google Scholar] [CrossRef]
  55. Jahnke, C.C.; Culick, F.E.C. Application of Bifurcation Theory to the High-Angle-of-Attack Dynamics of the F-14. J. Aircr. 1993, 31, 26–34. [Google Scholar] [CrossRef] [Green Version]
  56. Jahnke, C.C. On the Roll-Coupling Instabilities of High-Performance Aircraft. Phil. Trans. R. Soc. Lond. A 1998, 356, 2223–2239. [Google Scholar] [CrossRef]
  57. Dul, F.; Lichota, P.; Rusowicz, A. Generalized Linear Quadratic Control for Full Tracking Problem in Aviation. Sensors 2020, 20, 2955. [Google Scholar] [CrossRef] [PubMed]
  58. Lichota, P. Multi-Axis Inputs for Identification of a Reconfigurable Fixed-Wing UAV. Aerospace 2020, 7, 113. [Google Scholar] [CrossRef]
  59. Liebst, B.S. The dynamics, prediction, and control of wing rock in high-performance aircraft. Phil. Trans. R. Soc. Lond. A 1998, 356, 2257–2276. [Google Scholar] [CrossRef]
  60. Pietrucha, J. Modern Techniques for Active Modification of the Aircraft Dynamic Behaviour. J. Theor. Appl. Mech. 2000, 38, 132–156. Available online: http://www.ptmts.org.pl/jtam/index.php/jtam/article/view/v38n1p131 (accessed on 5 January 2021).
  61. Galiński, C.; Mieloszyk, J. Results of the Gust resistant MAV Programme. In Proceedings of the 28th International Congress of the International Council of the Aeronautical Sciences, Brisbane, Australia, 23–28 September 2012; Paper ICAS 2012-3.1.1. Available online: http://www.icas.org/ICAS_ARCHIVE/ICAS2012/PAPERS/186.PDF (accessed on 5 January 2021).
  62. Dżygadło, Z.; Kowaleczko, G.; Sibliski, K. Method of Control of a Straked Wing Aircraft for Cobra Manoeuvres. In Proceedings of the 20th Congress of International Council of Aeronautical Sciences, ICAS’96, Sorrento, Italy, 8–13 September 1996; Paper ICAS-96-3. 7.4. pp. 1566–1573, ISBN 1 56347-219-8. Available online: https://www.icas.org/ICAS_ARCHIVE/ICAS1996/ICAS-96-3.7.4.pdf (accessed on 5 January 2021).
  63. Sibilski, K. An Agile Aircraft Non-Linear Dynamics by Continuation Methods and Bifurcation Theory. In Proceedings of the 22nd Congress of International Council of Aeronautical Sciences, ICAS’2000, Harrogate, UK, 27 August–1 September 2000; Paper ICAS-2000-3. 11.2. ISBN 0 9533991 2 5. Available online: https://www.icas.org/ICAS_ARCHIVE/ICAS2000/ABSTRACTS/ICA3112.HTM (accessed on 5 January 2021).
Figure 1. Unit circle—diagram showing basic periodic orbit bifurcations—periodic orbit bifurcations. 1. The actual eigenvalue exceeds +1. Periodic limit points appear in this case. 2. The actual eigenvalue exceeds −1. Period doubling bifurcation occurs in this case. In the proximity of this point, a stable periodic orbit with a period of T becomes unstable and a new stable periodic orbit, with a period of 2T appears. This type of stability loss leads to chaotic motions. 3. Two conjugated eigenvalues leave the unit circle. Stable or unstable trajectories surround an unstable bifurcation orbit [32].
Figure 1. Unit circle—diagram showing basic periodic orbit bifurcations—periodic orbit bifurcations. 1. The actual eigenvalue exceeds +1. Periodic limit points appear in this case. 2. The actual eigenvalue exceeds −1. Period doubling bifurcation occurs in this case. In the proximity of this point, a stable periodic orbit with a period of T becomes unstable and a new stable periodic orbit, with a period of 2T appears. This type of stability loss leads to chaotic motions. 3. Two conjugated eigenvalues leave the unit circle. Stable or unstable trajectories surround an unstable bifurcation orbit [32].
Applsci 11 01524 g001
Figure 2. Coordinate systems with written micro aerial vehicle motion equations; (a) systems related to axes of inertia (so-called “aircraft system”); (b) coordinate system related to flow (so-called “velocity system” or “wind axis system of co-ordinates”).
Figure 2. Coordinate systems with written micro aerial vehicle motion equations; (a) systems related to axes of inertia (so-called “aircraft system”); (b) coordinate system related to flow (so-called “velocity system” or “wind axis system of co-ordinates”).
Applsci 11 01524 g002
Figure 3. Strake–wing MAV “Pszczoła” (Bee) [49,50].
Figure 3. Strake–wing MAV “Pszczoła” (Bee) [49,50].
Applsci 11 01524 g003
Figure 4. Examples of aerodynamic derivative waveforms as function of angle of attack α, and reduced frequency f [50].
Figure 4. Examples of aerodynamic derivative waveforms as function of angle of attack α, and reduced frequency f [50].
Applsci 11 01524 g004
Figure 5. Bifurcation diagram. Quasi-steady states for various elevator deflection angles; α(δe).
Figure 5. Bifurcation diagram. Quasi-steady states for various elevator deflection angles; α(δe).
Applsci 11 01524 g005
Figure 6. Bifurcation diagram. Quasi-steady states for various elevator deflection angles; P(δe).
Figure 6. Bifurcation diagram. Quasi-steady states for various elevator deflection angles; P(δe).
Applsci 11 01524 g006
Figure 7. Bifurcation diagram. Quasi-steady states for various elevator deflection angles; α(δe).
Figure 7. Bifurcation diagram. Quasi-steady states for various elevator deflection angles; α(δe).
Applsci 11 01524 g007
Figure 8. Bifurcation diagram. Quasi-steady states for various elevator deflection angles; V(δe).
Figure 8. Bifurcation diagram. Quasi-steady states for various elevator deflection angles; V(δe).
Applsci 11 01524 g008
Figure 9. Simulation of Wing-rock oscillations; waveforms of selected flight parameters.
Figure 9. Simulation of Wing-rock oscillations; waveforms of selected flight parameters.
Applsci 11 01524 g009
Figure 10. Diagram of bifurcation control matching, which enables recovery from a spin (skipping through the so-called “cusp catastrophe”) [25] with the waveform of elevon deflection angle changes.
Figure 10. Diagram of bifurcation control matching, which enables recovery from a spin (skipping through the so-called “cusp catastrophe”) [25] with the waveform of elevon deflection angle changes.
Applsci 11 01524 g010
Figure 11. Spin simulation. Angle waveforms for selected flight parameters.
Figure 11. Spin simulation. Angle waveforms for selected flight parameters.
Applsci 11 01524 g011
Figure 12. “Cobra” manoeuvre phases (top), “Bee” MAV during a “Cobra” manoeuvre, time-lapse photos (bottom) [61].
Figure 12. “Cobra” manoeuvre phases (top), “Bee” MAV during a “Cobra” manoeuvre, time-lapse photos (bottom) [61].
Applsci 11 01524 g012
Figure 13. “Cobra” manoeuvre simulation. Time waveforms of selected flight parameters.
Figure 13. “Cobra” manoeuvre simulation. Time waveforms of selected flight parameters.
Applsci 11 01524 g013
Figure 14. “Cobra” manoeuvre. Poincare maps.
Figure 14. “Cobra” manoeuvre. Poincare maps.
Applsci 11 01524 g014
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Nowakowski, M.; Sibilski, K.; Sibilska-Mroziewicz, A.; Żyluk, A. Bifurcation Flight Dynamic Analysis of a Strake-Wing Micro Aerial Vehicle. Appl. Sci. 2021, 11, 1524. https://doi.org/10.3390/app11041524

AMA Style

Nowakowski M, Sibilski K, Sibilska-Mroziewicz A, Żyluk A. Bifurcation Flight Dynamic Analysis of a Strake-Wing Micro Aerial Vehicle. Applied Sciences. 2021; 11(4):1524. https://doi.org/10.3390/app11041524

Chicago/Turabian Style

Nowakowski, Mirosław, Krzysztof Sibilski, Anna Sibilska-Mroziewicz, and Andrzej Żyluk. 2021. "Bifurcation Flight Dynamic Analysis of a Strake-Wing Micro Aerial Vehicle" Applied Sciences 11, no. 4: 1524. https://doi.org/10.3390/app11041524

APA Style

Nowakowski, M., Sibilski, K., Sibilska-Mroziewicz, A., & Żyluk, A. (2021). Bifurcation Flight Dynamic Analysis of a Strake-Wing Micro Aerial Vehicle. Applied Sciences, 11(4), 1524. https://doi.org/10.3390/app11041524

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop