Flight Stability of Rigid Wing Airborne Wind Energy Systems

: The ﬂight mechanics of rigid wing Airborne Wind Energy Systems (AWESs) is fundamentally different from the one of conventional aircrafts. The presence of the tether largely impacts the system dynamics, making the ﬂying craft to experience forces which can be an order of magnitude larger than those experienced by conventional aircrafts. Moreover, an AWES needs to deal with a sustained yet unpredictable wind, and the ensuing requirements for ﬂight maneuvers in order to achieve prescribed control and power production goals. A way to maximize energy capture while facing disturbances without requiring an excessive contribution from active control is that of suitably designing the AWES craft to feature good ﬂight dynamics characteristics. In this study, a baseline circular ﬂight path is considered, and a steady state condition is deﬁned by modeling all ﬂuctuating dynamic terms over the ﬂight loop as disturbances. In-ﬂight stability is studied by linearizing the equations of motion on this baseline trajectory. In populating a linearized dynamic model, analytical derivatives of external forces are computed by applying well-known aerodynamic theories, allowing for a fast formulation of the linearized problem and for a quantitative understanding of how design parameters inﬂuence stability. A complete eigenanalysis of an example tethered system is carried out, showing that a stable-by-design AWES can be obtained and how. With the help of the example, it is shown how conventional aircraft eigenmodes are modiﬁed for an AWES and new eigenmodes, typical of AWESs, are introduced and explained. The modeling approach presented in the paper sets the basis for a holistic design of AWES that will follow this work.


Introduction
Airborne Wind Energy (AWE) is the field of wind energy which aims at harvesting power from the wind through airborne systems. Airborne Wind Energy Systems (AWESs) can access wind resources at higher altitudes and they can obtain the same power output of conventional wind turbines with a drastic reduction of mass, which could lead to a cost decrease [1,2]. Many concepts are being developed and they are classified according to how the aerodynamic force, needed for the power production, is generated. A classification of the currently pursued AWE concepts, their flight operations and developers is given in [3]. The present paper is focusing on the development of AWESs which produce aerodynamic force by flying crosswind. Crosswind AWESs can generate power in two ways [3]: with on-board wind turbines (Fly-Gen AWESs) or with a generator placed at the ground station (Ground-Gen AWESs), which can be fixed or moving [3]. Fly-Gen AWESs produce power on board and transmit it to the ground station through electric wires embedded in the tether [4,5]. Ground-Gen AWESs with fixed ground station produce power cyclically. In the power generation phase, they produce power by pulling and unwinding the tether, which is connected to a generator. In the recovery phase, they fly back, spending some power, while the tether is wound back. Ground-Gen AWESs are being developed with soft kites [6,7] or rigid wing aircrafts [8][9][10].
To become a viable alternative to conventional wind turbines, AWESs flying crosswind need to display a sufficient reliability, thus allowing operations over long time frames. Therefore, design choices should keep robustness to external disturbances and failures as a premium goal. This in turn can be achieved through an increased knowledge of the dynamics of the flying craft, as a means to quantify how some characteristics of the AWES impact its performance, so as to allow suitably sizing the machine to meet all expectations.
The flight mechanics problem of AWESs is fundamentally different from the one of conventional aircraft for mainly three reasons: the presence of the wind, needed for the power generation (the average wind speed does not influence conventional aircraft flight stability analyses, as the inertial coordinate system can be considered to move with the average wind speed, but this is not the case for AWES), the presence of the tether, which transfers all forces acting on the airborne unit to the ground, and the need for maneuvers to keep the aircraft in the correct trajectory and attitude. It is therefore clear that conventional design techniques assuming straight and level flight can hardly be used for AWESs.
Makani Power [4], before the shutdown, started investigating the possibility of designing the flight path and aerodynamic characteristics of the AWES to achieve passive stability [11,12]. Stable AWESs would maintain the flight path with the least amount of control activity. Previous works [13][14][15] analyze the flight stability of kites not undertaking the crosswind motion for power generation. In this work, the kite crosswind motion needed for power generation is included in the analysis.
To achieve stability, the flight path should be properly chosen: a circular path makes the problem nearly axial-symmetric, and there exists one turning radius which ensures maximum tether force and potentially maximum power production [16]. As this turning radius maximizes the tether force, a stable-by-design AWES, if perturbed from this path, would tent to return to the operational point which maximizes tether force. This flight path is then selected for the present investigation.
The AWES aerodynamic design should help, as much as possible, the control system to achieve robust operations. An analytical modelling approach to study flight stability is presented in this text. The six d.o.f. equations of motion are linearized about a fictitious steady-state motion of the AWES in the circular path. No real steady-state can be achieved during power generation because of the continuous maneuvers of AWESs. Indeed, gravitational forces and aerodynamic forces related to the mean elevation angle act periodically on the kite. The fictitious steady-state motion, also called trim in this paper, is then computed by considering all the fluctuating terms as disturbances. The peculiarity of the selected flight path is that, if the fictitious steady-state is considered, centrifugal forces, induced by the constant turning maneuver, are balanced by the radial component of the force on the tether. In this way, lift is not used for the turn maneuver. The derivatives of external forces and moments are finally taken about the fictitious steady-state, to formulate the linearized problem.
Finally, the stability of the system can be studied by looking at the eigenvalues of the linearized system. The procedure to study stability is summarized in Figure 1.
Step A

Find trim solution
Step C

Formulate linearized problem
Step B

Compute external force derivatives
Step D

Study eigenvalues
The steady state is defined by considering all fluctuating terms as disturbances.  The linearized model allows for a simplified investigation of the influence of main wing geometry (position, area, aspect ratio, sweep, dihedral), main wing aerodynamics, control surfaces aerodynamics and geometry (area, aspect ratio, position), tether attachment position, tether mechanical and aerodynamic properties and mass properties of the AWES on flight stability. Also, since it does not primarily impact the dynamics of the flying craft, in this paper the power generation mechanism is not modeled. In this way, the model presented here can be used, with small modifications, to model both Ground-Gen and Fly-Gen AWESs.
As previously stated, the so-obtained model of the plant allows a good physical understanding of the effect of design choices on performance, and especially on stability. This in turn enables the formulations of design guidelines for the geometry and aerodynamics of rigid wing AWESs. Furthermore, the quality and characteristics of the model introduced in this paper bends itself to an adoption inside iterative optimal design tools, as well as for control design and tuning tasks.
An AWES designed to be stable, with the approach proposed in this work, has no guarantees to converge to the prescribed trajectory without control inputs. The formulation proposed here is meant to study the system dynamics of the AWES set to a state representative of its flight during the power generation loop. Therefore, the term stability in this paper refers to the ability of the system of returning to the fictitious steady-state if perturbed from it. As the fictitious steady-state is representative of the flight during the power generation loop, flight dynamic performances and control efforts over the loop are expected to benefit from a stable design.
The paper is organized as follows. Section 2 presents the non-dimensional linearized equation of motion derived for a generic AWES and Section 3 deals with the derivation of the external forces acting on the AWES. These two sections represent the mathematical basis for this paper and for subsequent work on the holistic AWES design and are therefore intended to be very comprehensive and detailed. Section 4 deals with the numerical implementation of the above mathematical model in MATLAB ® , while Section 5 presents a preliminary validation of this tool. Finally, Section 6 reports the results for different sets of configurations of increasing complexity and Section 7 closes the manuscript summarizing the main findings of this research.

Coordinate Systems
Three coordinate systems are defined to derive the equations of motion. The ground coordinate system F G (Figure 2) is centered at the ground station and it is inertial. Its origin is denoted as G. Z G points upwind and X G toward the ground. F G can also be interpreted as the wind reference frame. As F G is treated as an inertial frame, the inertial force components caused by the rotation of the wind vector are neglected. A second coordinate system F R (Rotating) is defined such that it moves on a circumference of constant radius R 0 , rotating around its Z R axis: X R points along the tangential direction and Y R outward. X R and Y R define the rotor plane ( Figure 2). The origin of F R is named R and lies on the plane connecting the ground station, the circumference center and the tether attachment point on the kite. With this definition, F R does not have a constant rotational speed.
The rotation matrix R G→R , which describes a coordinate transformation from F G to F R , is defined by two sequential planar rotations, associated with the mean elevation angle β and the angular position Ψ of F R along the circumferential trajectory of the craft (see Appendix A for rotational matrix notation) Ground coordinate system F G and rotating coordinate system F R . The stability coordinate system F S (Figure 3) is centered at the tether attachment point and is moving and rotating with the kite. Its origin is denoted as S. Equations of motion are written in F S , which makes the formulation of the external forces compact. As F S is defined from the aircraft trajectory, and not from a reference frame attached to the kite, inertial and geometrical properties of the system need to be expressed in F S for each operational point (i.e., wind speed). This procedure is detailed in Section 4. The rotation matrix which describes a coordinate transformation from F R to F S is defined with three sequential planar rotations, around the third, second (rotated) and first (rotated) axes of the frame respectively In steady-state, F S and F R coincide. Thus, this definition of coordinate systems is particularly useful when studying perturbations about the steady-state, which is the case in the formulation adopted in this paper.
Rotating coordinate system F R and stability coordinate systems F S . On the left, a graphical representation of the three sequential planar rotations from F R to F S . As angles φ, θ and ψ are assumed to be small and R R RS is linearized about the condition φ = θ = ψ = 0, they can be interpreted as rotation angles about axes X R , Y R and Z R , which is shown on the right.
As the final aim of this work is to linearize the system about its steady-state, (φ, θ, ψ) are assumed to be small and R R RS is linearized, so that where 1 is the 3 by 3 identity matrix and the symbol (·) × applied to a vector represents its corresponding skew-symmetric tensor form. With these assumptions,

Position
To formulate the equations of motion in F S and evaluate the tether elastic force, the position of the ground station with respect to F S , expressed in F S (X S S→G ), is to be found. It can be expressed as the summation of X S S→R and X S R→G : where with R 0 the circular path radius, Z 0 the distance of the F G from the rotor plane AA (see Figure 2) and X R RS = [x rs , y rs , z rs ] T (see Figure 3).

Angular Velocity
The angular velocity of F S with respect to F G , expressed in F S takes the form

Relative Wind Speed
The wind velocity is along the Z G axis and has a negative direction. A constant and uniform wind field, parallel to the ground, is considered. Its definition by components in the stability coordinate system can be written as In the formulation of dynamic equilibrium, the fluctuating terms, functions of Ψ, are considered as disturbances. Therefore, V S w,Ψ can be decomposed in a non-fluctuating and a fluctuating term The fluctuating term V S w, f l needs to be considered when studying the system in the time domain. Furthermore, for systems with no elevation angle, the fluctuating terms are null. This case is representative of a setup where the ground station is set on top of a tower or other elevated anchor points. The contribution of the fluctuating terms is along the X S and Y S axes. On the X S axis, the relative wind is dominated by the kite motion, while on Y S axis the kite experiences the fluctuating terms as side-slip velocities.
By considering a generic position in F S , namely X S = [x, y, z] T , the undisturbed relative wind speed V S r in F S is defined as the the subtraction of the kite motion from the wind velocity, yielding where V ≡ [u, v, w] T = V S GS is the velocity of F S with respect to F G , expressed in F S . The trim state (fictitious steady-state) is indicated in this paper with the subscript 0, therefore yielding for the baseline condition Here u 0 is the kite velocity along X S at trim, r 0 = − u 0 R 0 (valid for a left-turning kite), ξ = x R 0 , η = y R 0 and G is the system glide ratio (comprehensive of the drag acting on the tether), defined as G is defined here as a velocity ratio and not the lift-to-drag ratio, as commonly done in aeronautics. However, for AWESs flying crosswind, the two definitions coincide.

Equations of Motion in F S
The stability coordinate system F S is centered on the tether connection, so the AWES center of mass can be elsewhere. The equations of motion for a non-barycentric moving coordinate system are where m is the kite mass, 1 is the 3-by-3 identity matrix, TCG is the position of the center of mass in F S , I is the inertia tensor by components in F S , F and T t are the external forces and moments.

Linearized Non-Dimensional Equation of Motion
To make the equations of motion non-dimensional, a diagonal matrix is introduced as where b is the reference wing span and c the reference chord. The unit force is defined as where ρ is the air density, A is the reference wing area and u 0 the AWES translational velocity (evaluated at the tether anchor point).
A non-dimensional form of the six-dimensional equation of motion is wherein and Now, in order to study the stability of the system, the equations of motion are linearized about the trim (fictitious steady-state) condition, yielding The second term in the latter form can be expressed as a function of the incremental velocity ∆u = u − u 0 as Upon substitution, the linearized equation of motion is therefore where A c = S −1 F cd × 0 − u × 0 m and ∆f = ∆X, ∆Ỹ, ∆Z, ∆L, ∆M, ∆Ñ T .

Linearized Dynamics: Explicit Formulation
Thirteen incremental state variables are necessary to fully describe the dynamics of the tethered AWES system, namely three translational speed components, three rotation rate components, a circumferential position, three spatial and three attitude components. This yields the following state vector: X = ∆u, ∆v, ∆w, ∆p, ∆q, ∆r, ∆Ψ, ∆x rs , ∆y rs , ∆z rs , ∆φ, ∆θ, ∆ψ T , where ∆Ψ describes the perturbed angular position with respect to the steady state (baseline) position. It can be noted that, as the problem has been formulated to be axialsymmetric, no dependence on this state variable is expected. To express kinematic relationship between state variables, the perturbed translational and angular velocity equations are considered. The translational velocity can be written as where The linearized form of the previous equation writes where (6)) and ω R× GR,0 X R GR = [u 0 , 0, 0] T = V 0 (X R GR is defined in Equation (5)). In matrix form The linearized angular velocity equation in matrix form (from Equation (6)) can be written as   ∆φ ∆θ ∆ψ Based on the equations introduced up to now, it is possible to write the linearized dynamics of the system in the form In this expression, the mass matrix M writes with 1 6×6 is the 6 by 6 identity matrix, R 0x = [R 0 , 0, 0] T andk = [0, 0, 1] T . Furthermore, the coefficient matrix A eq appearing on the r.h.s. yields The external force matrix A f is instead used to describe the linear variation of the external forces with respect to a variation of the state space variables. Its components are derived in Section 3. As a further step in dealing with dynamic equations, the system of linearized equations can be made non-dimensional. This is achieved by normalizing the state variables according to dimensional groups typically found in flight mechanics literature. In the process, a scaling matrix S x is defined as where in particular L 0 t is the tether length at rest. The array of non-dimensional state variables is then obtained as In a similar fashion, the non-dimensional mass, coefficient and derivatives matrices arẽ Therefore, the non-dimensional linearized equations of motion can be written as It can be noted at this level that the dynamic equilibrium is described by 12 scalar equations and 13 variables. To solve this, F S is taken such that it is in the same radial direction of F R . Therefore ∆x rs = ∆ẋ rs = 0. Thus, the column relative to this contribution can be removed from the system, which becomes balanced in 12 variables. The modelling choice of having F S aligned with F R comes from the need of having these two coordinates systems the closest possible, as the position of F S is linearized with respect to F R .

External Forces Derivatives
Considering the linearized dynamics introduced in Equation (32), in this section the contribution of the external forces are described. They will populate the matrix of derivativesÃ f . Since it is not clear which derivatives will influence the dynamics of AWESs most, all derivatives are analytically derived.

Modeling the Wing Contributions to Aerodynamics
The analytical formulation presented here is meant to model the physics of the AWES in a general way. However, a few assumptions are made to ease the derivation without imposing excessive constraints.
The main wing is assumed to have an elliptic planform with no twist. Wing dihedral and sweep are considered to be small. The induced velocities generated by the lifting surfaces are modelled with Prandtl's lifting line theory, which assumes a straight translatory motion. The helicoidal wake structure is then modelled as straight, which is a good approximation [17] for the problem investigated in this paper. The squared wing span is considered to be small compared to the circular path radius squared, or in analytic terms, b R 0 2 1.

Wing Coordinate System
When designing the main wing for stability, dihedral angle Γ and sweep angle Λ are expected to largely influence the results. Therefore, the aerodynamic derivatives shall consider these angles. F Wr , a reference system attached to the right wing, is defined to be with the Y Wr along the direction defined by the points along the quarter chord line, pointing to the tip of the right wing, as in Figure 4. The rotation matrix between F S and F Wr is defined by applying first a rotation around X S of the dihedral angle Γ, followed by a rotation around Z S of the sweep angle Λ. This yields where sweep and dihedral angles are assumed to be small and DS = [Γ, 0, Λ] T . Note that the dihedral angle Γ is positive with the right wing pointing down, and the sweep angle Λ is positive with the right wing pointing backwards.  The coordinate system of the left wing F Wl can be defined in similar manner. The rotational matrix from F S to F Wl is The rotational matrices of the two wings can be unified as follows where y |y| = y |y| R 0 R 0 = η |η| is the signum function of the wing span coordinate, over the full left and right wing span.
A given position y w along the wing axis Y W with respect to the tether attachment can be expressed in F S (X S T→W ) as where X S T→W 0 is the position of the wing root with respect to the tether attachment. By looking at the second row of Equation (36), for small dihedral and sweep angles, it is found that the coordinates along Y S and Y W are equivalent (TW y = y w ).

Relative Wind Speed in F W
Aerodynamic forces are computed in the wing coordinate system. Thus, the relative wind speed in the wing coordinate system F W is where the relative velocity V S r and its components are defined in Equation (9). Note that (x,y,z) in Equation (9) are still the coordinates in F S . Therefore, V W r is the relative wind velocity in F W given a point in the F S .
For the crossflow principle, only the components on the X W and Z W axes produce aerodynamic forces. The relative velocity in the wing coordinate system, needed for the evaluation of aerodynamic forces, does not take into account the component along the Y W axis. Taking these features into account, the modulus of the relative speed squared is written as the summation of the first and third rows of Equation (37) squared where terms proportional to the sweep and dihedral squared are neglected (small higher order terms). Note that the two first terms are not dependent on the dihedral and sweep, so they are the only surviving results for a straight wing. The inflow angle γ w , assumed to be small, is defined as the angle between the relative velocity and the X W axis measured in the (X W , Z W ) plane (see Figure 5 for a straight wing). Angle γ w is used to project lift and drag in F W ; it can be written as the ratio between the third and the first row of Equation (37) where the first term is for a straight wing and the odd functions take into account wing geometry. Terms proportional to the sweep and dihedral squared are neglected (small higher order terms). Figure 5. Velocity triangle and aerodynamic forces for the main wing.

Trim Condition
In trim conditions, the relative velocity, using Equation (38), and the inflow angle, using Equation (39), take the forms where the assumption η 2 = y R 0 2 1 has been invoked and terms proportional to sweep and dihedral angle neglected. Both the relative velocity squared and inflow angle are linear functions of η.
Considering strip theory [18], the aerodynamic lift coefficient in F W at a given y w position of the wing is computed as where C Lα can be computed from the airfoil lift coefficient slope C lα which is a good approximation for an elliptic wing with no twist and ∂α ∂y = ∂γ w ∂y ≈ − 1 R 0 G . A similar procedure can be applied to the drag coefficient, such that C D (y) = C D (y = 0) − C Dα η G , where C D (y = 0) and C Dα can be found with lifting line theory No dependence of the lift coefficient on the side slip angle is considered. From now on, C L (y = 0) will be indicated as C L and C D (y = 0) as C D .
Given this modeling, the non-dimensional aerodynamic forces given by the main wing can be written as where the relative wind velocity squared V r (y) 2 (Equation (38)) and the aerodynamic coef- (35)) brings them into F S , where the integration happens. The term C L (y)γ w (y) in the aerodynamic coefficient matrix models the propulsive lift. The term X S TW × indicates the skew symmetric matrix of the application location of the aerodynamic forces and it is needed to compute moments. The matrix S −1 (Equation (13)) brings moments into forces, the unit force F (Equation (14)) makes the equation non dimensional and c(y) is the chord.
To get to a closed form solution of the integral, the odd functions of the wing span coordinate in the integral variable are to be separated from those which are even or constant. The first part of the integral function of Equation (44) can be written using Equation (36) and The very last term, including effects proportional to dihedral and sweep squared, can be neglected. Also the second last term can be left out, since the aerodynamic loads are only in the first and third axes [0, With these two considerations, Equation (45) can be written highlighting the dependence on the wing span non-dimensional coordinate η where P and D are even functions of η. For readability, the following matrices, derived in Appendix B, are defined The aerodynamic coefficient matrix at trim C 0 , introduced in Equation (44), is where γ w0 (y) is given in Equation (40) and C L (y) in Equation (41). The definition of C R0 and C η is useful when performing derivatives with respect to a generic state variable @ (Equation (30)) where the derivative of the angle of attack with respect to a generic variable is equal to the derivative of the inflow angle with respect to the same variable ∂α ∂@ = ∂γ w ∂@ . The three new aerodynamic coefficient matrices, defined to highlight the dependence of the aerodynamic coefficient matrix at trim C 0 (Equation (48)) and of its derivative ∂C ∂@ (Equation (49)) with respect to η are Deleting the odd functions of y, which provide a null contribution along the wingspan, the following form for non-dimensional force and moment at trim is obtained which can be written explicitly as Here y η is the y coordinate of the right wing centroid y η = I y The expressions of the non dimensional aerodynamic forceX (which should not be confused with the non-dimensional state variablesX) andZ correspond to the classical results for AWES modeled as a point mass flying in a straight crosswind motion. This is in accord with the findings on the flight path introduced in [16] and used here. Conversely,Ỹ appears just because of dihedral and sweep and the difference in relative wind speed between the inner and outer wing. Indeed, loads in the two wings are different and so is their projection on the Y S axis. The momentsL andÑ appear because of the difference in relative wind speed between an inner and outer wing. It is also noteworthy thatM basically does not depend on dihedral nor sweep.

Procedure for the Evaluation of Derivatives
The derivative of forces and moments with respect to a generic variable @, using Equations (40), (48) and (49), can be formalized as follows where ∂γ w ∂@ and ∂V 2 r ∂@ in the latter expression can be found by applying recursively the chain rule. For example, derivatives with respect top can be written as The systematical application of the chain rule for the other derivatives is shown in Appendix C.

Straight Wing Derivatives
Equations (38) and (39) express the relative velocity squared and the inflow angle for a wing profile. The contributions of the straight wing, dihedral angle and sweep can be analyzed individually as they appear in different terms. Consequently, the derivatives of γ w and V 2 r with respect to U , V and W for the straight wing can be provided separately, and are given in Table 1. The chain rule, following the procedure given in Appendix C, can be applied to find the analytic expression of such derivatives.
To present one example of the methodology for deriving expressions as in Table 2, the derivation of ∂f w ∂ũ is here displayed. Following the procedure introduced in Section 3.1.4, the two terms necessary for the evaluating the derivatives (Equation (53)) are where ∂γ w ∂U and ∂V 2 r ∂U are found in Table 1 and ∂U ∂u = 1 is given in Table A1. Performing an integral along the wingspan, in order to find the aerodynamic derivative with respect toũ, Equation (53), considering Equation (55), takes the form Table 2 summarizes all the non-null derivatives for the case of a straight wing (the passages leading to these forms are provided in Appendix D).

Derivatives Due to Dihedral and Sweep Angle
After showing how to write the derivatives of the aerodynamic force and moments for the case of a straight wing, the terms depending on sweep and dihedral are considered. The derivatives of γ w and V 2 r with respect to U , V and W are given in Table 3, for the chain rule application. Table 3. Derivatives of γ w and V 2 r with respect to U , V and W due to dihedral and sweep angle.
In Table 4 all the non-null derivatives due to sweep and dihedral angle are reported. The passages leading to the results presented here are given in Appendixes E and F.

Control Derivatives: Ailerons
Equation (52) showed that in steady-state the linear variation of the relative wind and of the inflow angle are generating a moment L. To compensate for that, a control strategy to trim the kite and set this moment to zero could be based on the use of ailerons, as shown in Figure 6. For this reason, it is especially interesting to model the effect of ailerons on aerodynamics. This can be done again via strip theory, thus making use of most of the formulation already introduced. The corresponding control derivative, expressing the sensitivity of aerodynamic forces and moments to control variables, can be studied as follows. By hypothesis, a relatively small sensitivity to aileron control can be considered, since the size of ailerons (hence their control authority) is generally limited compared to that of the wing. Furthermore, derivatives and trim forces are computed neglecting the dihedral and sweep angle for simplicity and it is assumed that outside the wing sections where ailerons are deflected, the span-wise aerodynamic loading remains unaffected.
Assuming a symmetric displacement of the right and left ailerons, the force associated to their deflection is where ∆C = [∆C L γ w − ∆C D , 0, −∆C L ] T is the aerodynamic coefficient matrix and the function f (y) takes a value of 1 along the portion of the wingspan between the aileron extremities Contributions to lift and drag coefficients by the ailerons are modelled as deltas with respect to the main wing lift and drag. Typically, ∆C L and ∆C D are provided as function of the ailerons deflections. Therefore, the corresponding changes in the values of coefficients at trim are Consequently, the derivative of the aerodynamic coefficient matrix can be written (according to the same procedure introduced for other force and moment derivatives) as where it is assumed that ∆C L and ∆C D do not depend on the angle of attack. The force vector due to the ailerons at trim is therefore where the new terms appearing in the latter expression are P a = S −1 1; [TW x , 0, TW z ] × , Table 5 the derivatives due to the ailerons are reported. Note that ∂γ w ∂@ and ∂V 2 r ∂@ are the same as for the straight wing case (Table 1).
It can be observed that the moment around the X S axis due to ailerons can be computed asL Therefore, the necessary ∆C L to setL to zero, neglecting the contribution from sweep and dihedral, can be estimated as

Aerodynamics of the Horizontal Tail
For the sake of simplicity, the horizontal tail surface is lumped at the symmetry plane of the aircraft. Thus, only the relative wind velocity and the inflow angle at the mid-airfoil of the horizontal tail are to be considered. The inflow angle is the summation of the relative wind speed (Equation (9)) and of the downwash of the main wing where ε d is evaluated using the formulation in [19], which models the downwash in the aft tail considering also the contribution of the main wing sweep, yielding Here, for an elliptic wing κ v = 1, κ b = π 4 , and κ p takes into account the aft tail position with respect to the main wing. Finally κ s models the effect of the main wing sweep on the downwash. Quantities C L,w and AR w refer to the main wing. Based on the model just introduced, the variation in the downwash angle for a change in the angle of attack of the main wing can be computed as The non-dimensional force given by the horizontal tail can be written as where A h is the horizontal tail area and C L and C D its aerodynamic coefficients. The derivative of the force with respect to a generic variable @ is where The derivative of the aerodynamic coefficient matrix C h with respect to a generic variable @ is The derivative of the inflow angle γ h for the horizontal tail have two contributions, according to the model in Equation (64) where f ε (@) is clearly a function of the derivation variable. In particular, In Table 6, the derivatives with respect to the pertinent states of the AWES are provided. Furthermore, note that derivatives given in Table 1 can be used also for the horizontal tail.

Aerodynamics of the Vertical Tail
The vertical tail is meant to give lateral stability and to trim the aircraft over the circular path. As the tether attachment position is not located in the center of mass in general, the centrifugal forces generate a moment around the Z S axis in steady state. The vertical tail should then set the overall moment to zero in steady state, as shown in Figure 7. Figure 7. Main forces acting along the Y S axis. The vertical tail is stabilizing the effect of the centrifugal forces.
Similar to the horizontal tail, the vertical tail is modelled as a lumped lifting surface, located on the aircraft symmetry plane (i.e., TV y = 0) Considering the angle θ 0 , defined as the angle between the vertical tail quarter-chord axis and the Z S axis, the relative wind velocity perpendicular to the vertical tail sections can be written based on its components as The inflow angle on the vertical tail γ v is defined as The force given by the vertical tail can be written in vector form as where A v is the vertical tail area and C L and C D the aerodynamic coefficients of the vertical tail. Correspondingly, the derivative of the force vector with respect to a generic variable @ can be written as where T and the derivative of the aerodynamic coefficient matrix ∂C v ∂@ can be written as follows In Table 7 the derivatives of γ v and V 2 r with respect to U , V, W are reported, while in Table 8 the contributions to the force and moment derivatives with respect to the AWES state variables are provided for the vertical tail.

Tether Reaction Force
The tether is considered as a linear spring. Since the equations of motions have been written in a coordinate system centered in tether attachment, the tether is not generating any moment.
Treating the tether reaction force is possible in a linear framework again recurring to derivatives. The starting point of the derivation process is the writing of the reaction force components in the adopted body system, making use of a standard structural modeling for cables.
The non-dimensional elastic force yields by components where E is the tether Young Modulus, A t the tether sectional area, L t and ε are the tether length and the tether strain respectively, defined respectively as wherein L 0 t is the tether length at rest and L t is the current tether length. The reference tether length L t0 and the strain at trim ε 0 are respectively Similar to what has been done for the aerodynamic modeling of the horizontal and vertical tail, the derivatives of the tether force (Equation (77)), considering Equations (78) and (79), with respect to a generic variable @ take the general form , and in turn Φ is the opening angle of the ideal cone swept along the flight path, as shown in Figure 2. The derivatives necessary to estimate the elastic tether force derivatives (Equation (80)) are given in Table 9.
In addition to the elastic force, the tether introduces also a drag component, which can be modelled as a concentrated force at the tether attachment [16], yielding Here C ⊥ is the tether drag coefficient, d t the tether diameter and A the reference wing area. The tether mass and its dynamics are neglected for the sake of simplicity. In Table 10, the derivatives for the tether force are provided.

Modeling the Effect of Gravity Force
Since the center of mass is not located at the center of the coordinate system, gravity generates moments. This is generally dissimilar to the most typical treatment of the dynamic equations for aircraft, which are typically put in barycentric body components.
For AWESs, and considering the reference at the anchor point of the tethering system, the non-dimensional force due to gravity can be written as follows wherein R SGî can be decomposed in a fluctuating and a steady term, yielding Again, a derivative of the steady component of the gravitational force with respect to a generic variable @ can be written as where ∂R SR ∂@ is given in Table 9. In Table 11, the corresponding derivatives of gravitational force are provided.

Computational Implementation of the Methodology
The linearized dynamics of the system, introduced in Sections 2 and 3, are suitable for a quick parameterized analysis of a novel AWES. An implementation has been carried out in MATLAB ® environment, aimed at the eigenanalysis and forward-in-time simulation of a system with an assigned geometry. The suite has been named LT-GliDe (Linearized Tethered Gliding system Dynamics), and it is a module of an under-development multidisciplinary design and optimization framework for rigid wing AWES, named T-GliDe (Tethered Gliding system Design).
For a given AWES, the code computes the trim condition, evaluates the external force derivatives and extracts the eigenvalues of the linearized system (Equation (32)). The AWES main wing, horizontal and vertical tails and tether properties are described by the variables listed in Tables 12 and 13. The shortest number of variables is considered, to keep the problem as simple as possible, and to be able to analyze the influence of each property of the AWES configuration on in-flight stability performance. Table 12. Inputs needed for the derivatives evaluation and the eigenproblem formulation in LT-GliDe (Linearized Tethered Gliding system Dynamics) (first part).

Main Wing Horizontal Tail Vertical Tail Tether
The relative positioning of the lifting surfaces, as well as the inertia matrix of the craft, are typically given in a reference system attached to the AWES, which does not depend on the operational regime and might be centered in principle in any point of the aircraft. This coordinate system is called here body coordinate system F B . In Table 13, data of the lifting surfaces, and similarly the elements of the inertia matrix, are assigned in this coordinate system. To express the geometrical and inertial quantities given in F B into F S , where the equations have been written, θ BS is defined as the angle around the Y B axis which describes a rotation from F B to F S . A graphical representation of F B and F S is given in Figure 8, highlighting two possible flight conditions. Finally, the wind velocity and the elevation angle describe the operational regime.   Off-diagonal inertia θ SB Rotation from F B to F S Once the inputs are defined in LT-GliDe, the trim condition is computed. In particular, the turning radius R 0 , the kite velocity u 0 , the tether strain and the lifting properties of the control surfaces describe the trim condition. Algorithm 1 is used to find the trim condition for an AWES flying the circular path.
Algorithm 1: Algorithm for the trim evaluation.
The code is such that five different cases of increasing complexity can be analyzed. In Table 14 they are summarized. Results from the five cases are introduced, so that complexity is incrementally increased and they can be compared. Case a models the motion of an aircraft in steady and leveled flight without tether (this is basically a standard aircraft in steady, horizontal flight). Case b models the motion of a tethered AWES in a straight motion, with an imposed velocity. Case c models the motion of a tethered AWES in a straight crosswind motion. Case d models the motion of a tethered AWES in a circular motion, with an imposed velocity. Case e models the motion of a tethered AWES in a circular crosswind motion.
Once the trimmed condition has been computed, the linearized dynamics are populated and the eigenproblem can be studied. The eigenvectors of the linearized system give information about its eigendynamics, while the eigenvalues provide a crucial information on the stability of the system.

Eigenproblem Validation
To validate the problem formulation, an ultralight glider named Zefiro, designed internally at the Department of Aerospace Science and Technology of the Politecnico di Milano [20], has been considered. This aircraft has been selected because it is expected to have similar characteristics to AWESs, in term of mass properties and flight dynamics for straight and leveled flight, as per available internal studies and reports. The geometry and mass properties necessary for this study are given in Table 15. Zefiro features −5 • of dihedral angle and 5 • of sweep. Table 15. Geometry and mass properties of the ultralight aircraft Zefiro. The tether attachment position has been arbitrarily picked.

Main Wing Horizontal Tail Vertical Tail Geometry Mass
Once all the derivatives of external forces are evaluated and the linearized problem populated (Equation (32)), the eigenvalues can be computed and the system stability studied. To check the correct implementation of the equations of motion, the eigenfrequencies of short period, phugoid and Dutch roll are compared with values predicted by the models presented in [21], finding a good agreement. In Table 16 the eigenfrequencies and damping ratios computed with the different methods are shown for a velocity of 80 m/s.
In the Reduced LT-GliDe model, the computation of the eigenvalues is performed by considering the same state variables as the approximated analytical formulation of Stevens [21]. The values found in this case are indeed more similar to the values found with the analytical formulations.

Aerodynamic Derivatives Validation
To validate the aerodynamic derivatives analytically derived as shown in the previous part of the paper, the computed values have been compared with the outcome of TORNADO [22], a vortex lattice method.
In order to make TORNADO model the same aerodynamic problem, the MATLAB ® optimization function fmincon, with sqp algorithm, is used to impose the same lift coefficients in the lifting surfaces. In TORNADO, real airfoil data are used to model the wing and tail. Only the lifting surfaces are modelled and the fuselage is left out. The inflow is set to consider both atmospheric wind speed and aircraft motion. To be consistent with the analytic formulation, the drag coefficient at zero lift C D0 is included in the drag evaluation. The external forces derivatives are finally obtained by finite difference. In the optimization function, the pitch angle of each lifting surface (main wing and tails) and the ailerons deflection are included as design variables. The square of the difference between the lift coefficients found in LT-GliDe with Algorithm 1 and the lifting coefficients found in TORNADO by pitching the wings is included as objective function. Once the optimization converges (i.e., the two software model the same trim condition), the aerodynamic derivatives can be compared and used to study the eigenvalue problem.
The validation of the derivatives has been performed considering one element a time (e.g. only the main wing or only the vertical tail) and the full aircraft. In this way, the derivatives can be checked individually.
Moreover, the cases introduced in Table 14 can be analyzed in a sequential way, such that different terms in the derivatives can be studied. For example, ∂f w ∂ũ (Equation (56)) can be understood and expressed for the five cases as follow for a straight motion with wind (case c) for a circular motion with wind (case e) In Table 17, the eigenvalues for Zefiro in the circular motion case with wind at 8 m/s (case e), obtained with TORNADO and with LT-GliDe, are presented for different values of sweep and dihedral angles. The eigenvalues are considered close enough, so that the analytical models for the derivatives in LT-GliDe can be considered a good approximation of the vortex lattice solution of TORNADO. Table 17. Eigenvalues found with the analytical approach (LT-GliDe) and with the vortex lattice method (VLM) in the circular motion case with wind (case e) for different values of sweep and dihedral angle. Both the analytical approach introduced in this paper and the numerical one based on TORNADO can be used to estimate the aerodynamic derivatives of the aircraft. Once the aerodynamic derivatives are estimated (with one of the two methods), derivatives related to the tether and gravity, computed analytically, can be added to formulate the linearize problem and study the eigenproblem associated to the AWES dynamics. The advantage of the analytical approach is the negligible computational cost and the possibility of understanding analytically the relationship between design parameters and stability, which may turn particularly useful in a preliminary design phase where the configuration and geometry of the flying craft are still not frozen. The advantage of a numerical approach is the possibility to model any geometry. The two approaches are therefore to be seen as complementary: the analytical approach will be useful for physical understanding, preliminary design and optimization, while the numerical approach for the detailed aerodynamic design.

Results
The results introduced here are meant to show the potential of the proposed approach. A detailed analysis of a AWES designed for power generation is left to future works. In Tables 18 and 19, the eigenvalues for the five cases are reported and linked to each other. The straight and circular motion cases are analyzed in the following sections.

Case a
This case models the straight and leveled flight of a conventional aircraft. No novelty is introduced in this case, as the tether is not present in the model and the aircraft is flying straight. Conventional longitudinal and lateral modes can be found. A lift coefficient of the main wing of 0.09 is found at a velocity of 80 m/s. The aircraft is longitudinally stable, with a static margin of 20.3 %. For a conventional aircraft, eight eigenvalues are found. Typically, the longitudinal motion can be described by four state variables (ũ,w,q, θ), which correspond to four eigenvalues: two pair of complex conjugate eigenvalues describing the so-called short period and the phugoid modes. The lateral-directional motion can be described by four state variables (ṽ,p,r, φ), and it is typically characterized by a pair of complex conjugate (Dutch roll) and two real eigenvalues (spiral and roll subsidence mode).

Case b and c
For these cases, the tether is added and the mean elevation angle is set to zero, such that gravity is not influencing the results. The considered tether has a diameter of d t = 1 cm, the total length at rest is L 0 t = 400 m, the drag coefficient C ⊥ = 0.8 and the Young modulus E = 110 GPa. The linear spring constant associated to the tether is κ = A t E/L 0 t = 21 kN/m. The lift coefficient of the main wing is set to 0.9. Since the tether is part of the modelling, two new eigenvalues are expected, related to the position state variablesỹ rs andz rs . The motion is straight, thus the longitudinal and lateral motions can be studied separately.
In Table 18, the eigenvalues of case b and c are reported and associated with the eigenvalues of the airplane mode case. Eigenvalues and eigenmodes trends are consistent for the two cases, so that they are analyzed here together.  Figure 9 shows the longitudinal eigevalues for an increasing tether stiffness for case b. For E → 0, the tether deformation ε goes to infinity to produce enough force on the tether to balance the aerodynamic forces. As in this case the elastic effects are negligible, the tether force is a constant force applied on the tether attachment. This force could be understood as a gravitational force applied on the tether attachment and not on the center of mass.  As the elasticity grows, the short period increases the imaginary part. By performing the same plot for κ= (A t E)/L 0 t → ∞, it is found that the short period has a vertical asymptote corresponding to the eigenmode of the system spring-mass with the elasticity from the tether and the mass from the kite. The phugoid increases its real part till approximately κ ≈ 4 kN/m and later decreases it. A new real longitudinal eigenmode starts from zero and moves in the real axis in the negative direction. In this paper, this eigenmode is called Loyd mode, after Miles L. Loyd, who derived the first power equation of AWESs [23]. The Loyd eigenvalues collapse at very low tether stiffness to a value, suggesting that this mode is not influenced by the tether elasticity.
The Loyd mode involves mainly the longitudinal velocityũ. For AWESs, the longitudinal velocity in equilibrium is given by the balance of the projection of lift along the X S axis and drag. Whenũ is perturbed, the kite tends to go back to the equilibrium point. A graphical representation of this mode is given in Figure 10. In the figure, the AWES velocity is decreased of ∆u, such that the inflow angle γ w is increased with respect to the steady state value (γ w > γ w ). Lift is perpendicular to the relative wind speed and a positive force ∆X along the X S axis is generated. This force acts to restore the equilibrium condition. The eigenmode can be approximated by a simple first order differential equation in the general formũ =Ã(1, 1)ũ, such that the associated eigenvalue is where the aerodynamic coefficients are related to the full flying craft. For the case with no wind (case b), G → ∞, then For the straight motion case with wind (cases c), G = This approximation suggests that for a well designed AWES this eigenvalue is always real and negative, as the term inside the brackets in Equations (86) and (87) is positive.
To better characterize the longitudinal modes, the eigenvalues are tracked as function of the tether attachment position along the X S axis. Figure 11 shows the longitudinal eigenvalues for an increasing distance between the tether attachment and the center of mass. When TCG x = 0 m, the two points coincide, while when TCG x = 0.2 m (nominal value considered in the paper), the tether attachment is placed 20 cm aft. Interestingly, the phugoid for high TCG x collapses to two real eigenmodes, one of which moves along the real positive axis. The Loyd eigenvalue is basically fixed at low TCG x and it drastically decreases its real part approximately when the phugoid becomes unstable.  Figure 12 shows two lateral eigenvalues as functions of the tether stiffness. The lateral mode related to Dutch roll and to roll subsidence are basically not influenced by the tether stiffness, then they are left out of the plot. The spiral mode increases its real part until matching another lateral eigenvalue on the real axis. After the two coalesce, they become complex conjugates for further values of the parameter. The coalescence happens at really low stiffness values (κ ≈ 0 kN/m), and the eigenvalues collapse rapidly to a location, suggesting that this mode is not influenced by the tether elasticity. The eigenvalues for κ > 0 kN/m have the imaginary part much larger than the real part, suggesting that this mode is lightly damped. This new oscillating stable mode is dominated by a linear lateral motion (involvingỹ rs ), and small perturbations in φ and ψ and it is called here pendulum mode. Its natural frequency is ω n = 0.47 rad/s (case c).
By considering a pendulum composed by the tether and a kite subject to aerodynamic force, the oscillation frequency is ω pend = T mL t0 , where T is the tether force, m is the kite mass and L t0 is the tether length. The frequencies found in this simplified way is ω pend = 0.49 rad/s (case c), close to what has been found in the LT-GliDe analysis. By modifying the dihedral angle, while keeping the main wing aerodynamic center fixed to isolate the effect of the change, to Γ = 0 • , the pendulum mode becomes unstable: λ p = 0.04 ± 0.47. If it is modified to Γ = −10 • , the real part becomes more and more negative: λ p = −0.19 ± 0.42. Dihedral angle then largely influences the real part of the pendulum mode and should be considered to achieve stability. The spiral mode, for conventional aircrafts, depends on the dihedral angle: so does the pendulum mode for AWESs, which arises from it ( Figure 12).

Case d
The circular motion is now considered. In this case, the inflow varies along the main wing and ailerons are used to trim the aircraft. In particular, it is assumed that ailerons deflections do not influence the section drag coefficient (i.e., ∆C D = 0). Longitudinal and lateral modes cannot be separated and they shall be studied together. However, by looking at the numerical values in Table 19, short period, Dutch roll and roll subsidence have consistent values with the straight motion cases.  Figure 13 shows all the non-zero eigenvalues (apart from the roll subsidence) for the circular motion without wind (case d).
The pendulum mode is represented by a pair of complex conjugate eigenvalues starting from the origin. For κ > 1 kN/m, all points have the same value, suggesting that this mode is not influenced by the tether elasticity. The spiral mode eigenvalue moves along the real axis. For extremely low κ it moves towards the negative direction, and for κ ≈ 2 kN/m it starts increasing in its real part again. It is basically not influenced by the tether stiffness. The Loyd mode decreases its real part as in the straight motion case and from low tether stiffness its eigenvalues are fixed to the same value.

Case e
In case e the wind velocity is considered and the kite speed is computed from the force balance. Similar to the previous cases, in Figure 14 the eigenvalues are tracked for an increasing tether stiffness. Only the modes typical of AWESs are shown in figure. Interestingly, Loyd and spiral modes are coalescing and merging into a pair of complex conjugate eigenvalues. The corresponding mode is called positional mode in the paper. The positional mode is not influenced by tether stiffness, as for low values of rigidity the eigenvalues have the same value.
In Figure 15, the eigenvalues are tracked as functions of the tether attachment position. The Loyd and the spiral mode for an increasing distance between the center of mass and tether collapse to a pair of complex conjugate eigenvalues (positional mode). This new oscillating mode is stable and highly damped. By increasing again TCG x , the positional mode splits again in two real eigenvalues. The phugoid is also largely influenced by the tether attachment position. For TCG x = 0 m, this mode is unstable, while when moving the tether attachment aft, its real part gets negative.
The pendulum mode eigenvalue displays instead an opposite trend. While moving the tether attachment aft, its real part becomes positive and finally the two pairs of complex conjugate eigenvalues become two real and positive eigenvalues. For the straight motion case (Figure 11), for high TCG x the phugoid splits into two real eigenvalues, while for the circular motion case the pendulum splits instead. This suggests that, when designing the AWES for stability, the circular motion case with wind shall be considered to properly characterize the eigendynamics of the system: even if some of the eigenvalues for straight and circular motion are similar for some specific cases, trends over the full design space might change radically.
This analysis highlights the importance of carefully choosing the tether attachment position when designing for stability.

Conclusions
This paper introduces a methodology to study the flight stability of rigid wing Airborne Wind Energy Systems (AWESs). A fictitious steady state motion of the kite in a circular path is found by considering all fluctuating terms during the loop, such as gravity and aerodynamic forces related to the mean elevation angle, as disturbances. The selected circular path is characterized by a turning radius which makes the force on the tether to be maximized. Once the steady state is found, all derivatives of external forces are computed to formulate the linearized dynamics, and the flight stability can be studied from the ensuing model. As the fictitious steady state is representative of the flight over the loop, flight dynamic performances and control efforts are expected to benefit from a stable design.
The steady state forces and the aerodynamic derivatives are computed using an analytical approach. The computational cost of the linearized problem evaluation is negligible compared to other methods, which make this tool suitable for a preliminary design phase. Moreover, the analytical formulation allows for an intuitive and quantitative understanding of how stability is influenced by design parameters. Indeed, the modelling introduced in this work takes into account all the main design features of AWESs: main wing geometry (position, area, aspect ratio, sweep, dihedral), main wing aerodynamics, control surfaces aerodynamics and geometry (area, aspect ratio, position), tether attachment position, tether mechanical and aerodynamic properties and mass properties of the kite. Moreover, the non-uniform distribution of the inflow along the wing span, due to the rotational velocity of the kite, is also taken into account.
The presence of the tether introduces three new state variables compared to conventional aircraft flight mechanics. However, the variable describing the kite position in the circular path does not influence the linearized dynamics because the problem has been formulated to be axial-symmetric.
The model is implemented in MATLAB ® and the suite has been named LT-GliDe (Linearized Tethered Gliding system Dynamics). LT-GliDe is a module of an under-development multidisciplinary design and optimization framework for rigid wing AWESs, named T-GliDe (Tethered Gliding system Design). Five cases of increasing complexity can be studied: AWES operating in aircraft mode (without tether), AWES operating in tethered straight motion without wind, AWES operating in tethered straight crosswind motion, AWES operating in tethered circular motion without wind, AWES operating in tethered circular crosswind motion.
The aerodynamic analytical derivatives have been validated using vortex lattice method (TORNADO) and an internally-designed glider, called Zefiro, is used to validate the problem formulation for straight and leveled flight. Zefiro is then used as an example to show the features of this modelling approach. A tether is attached to Zefiro in a position such that all eigenvalues have a negative real part for the tethered circular crosswind motion case. The five cases are finally analyzed.
In straight tethered motion, longitudinal and lateral motion can be decoupled. Ten non-null eigenvalues are found, two more than a conventional aircraft in steady and leveled flight. It has been shown how the tether modifies the conventional eigenmodes and the physics of two new modes is explained. The Loyd mode is a mode typical of AWESs, which brings the kite back to the equilibrium longitudinal velocity after a perturbation. The pendulum mode arises by the merge of a new mode with the spiral mode and has the frequency of a pendulum characterized by the force applied on the tether. Its real part, defining the mode stability, is largely influenced by the main wing dihedral angle.
In circular motion, longitudinal and lateral equations are coupled and so are the eigenmodes. However, short period, Dutch roll and roll subsidence eigenvalues have consistent values with the straight motion cases. Eleven non-null eigenvalues are found in this case, three more than for conventional aircrafts. Two are due to the presence of the tether and one to the steady maneuver. Also in this case, the pendulum mode can be found, as well as the Loyd mode in the case without wind. For the crosswind circular motion case, the Loyd mode merges with the spiral mode to bear a pair of complex conjugate eigenvalues (positional mode).
It is investigated how the tether attachment position influences stability. By moving the tether attachment aft with respect to the center of mass, the phugoid becomes stable and the pendulum mode unstable. The tether attachment location should then be carefully designed to achieve stability of all modes. Flight mechanics should generally be studied with the circular crosswind motion model, as the eigendynamics for this case can be radically different from the more simplified other cases.
As LT-GliDe is suitable for preliminary design and optimization of rigid wing AWESs, it is being included in a system design framework (T-GliDe) and optimal designs will be analyzed in future works. It is expected that an aircraft designed with this approach will reduce control action and the need for control authority. Aero-elastic time-marching simulations are envisaged to study the eigenmodes and the controllability with higher fidelity codes.

Conflicts of Interest:
The authors declare no conflict of interest.

Nomenclature
The following symbols are used in this manuscript:  The rotational matrix, expressed in F A , which rotates from F A to F B is denoted with R A A→B = R A AB . Note that R A AB = R B AB . The rotational matrix R A AB can be used to rotate a vector a A , expressed in F A into a vector b A , expressed in the same coordinate system Since this operation is a rotation: b B = a A . It follows that To change the basis of a rotation tensor T, rotational matrices are to be applied To compose rotations, rotational matrices expressed in the same basis need to be applied in a sequential order To express rotations as a composition of planar rotations, the basis of the rotation matrix need to be expressed in the basis where the rotation happens. The rotation matrix R A A→C can be therefore be written as The following planar rotational matrices are defined (A7)

Appendix B. Arm Matrices Derivations
The following integrals are defined to have the dimensions of areas As the analytical derivatives are compared with the results from vortex lattice method, attention should be given to use the same locations for the velocities calculation. Indeed, TORNADO computes aerodynamic forces using the velocity at 3/4 chord (collocation points) and apply them at 1/4 chord to compute moments.
To be consistent with the vortex lattice method and take velocities at the collocation points, the following geometrical integrals are defined The even arm matrix P is defined as function of P 1 and P 2 : The even arm integral matrices are derived as follows Considering x and z as the collocation points position at 3/4 chord (X S TW (y) − c(y)/2, where X S TW is given in Equation (36)), the even arm matrices needed for the estimation of the derivatives with respect to angular rates are Similarly, P η,x and P η,z can be derived:

. Odd Arm Matrices
The odd arm matrix D is defined as function of D 1 and D 2 .
The odd arm integral matrices for the derivatives computations are derived as follows Considering x and z as the collocation points position at 3/4 chord, the odd arm matrices needed for the estimation of the derivatives with respect to angular rates are Similarly, D ηη,x and D ηη,z can be derived (A27 )  Table A1. Derivatives of U , V, W (Equation (9)) with respect to dimensional state variables.

Appendix D. Straight Wing Derivatives
Appendix D. 1