# Distributed Modeling and Control of Large-Scale Highly Flexible Solar-Powered UAV

Authors: Rui Wang; Zhou Zhou; Xiaoping Zhu; Xiaoping Xu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The modeling, stability, and control characteristics of a large scale highly flexible solar-powered UAV with distributed all-span multielevons were presented. A geometrically nonlinear intrinsic beam model was introduced to establish the structural/flight dynamics coupled equation of motion (EOM); based on it, the explicit decoupled linear flight dynamics and structural dynamics EOM were derived through mean axis theory. Undeformed, deformed, and flexible models were compared through trimming and modal analysis. Since the deformation of wing has increased the UAV’s moment of inertia about the pitch axis, the frequency of short period mode has obviously decreased for the deformed model. The serious coupling between short period mode and 1st bending mode also significantly influences the roots of short period mode of flexible model. So flexible model was the only one which is able to accurately estimate the flight dynamics behaviors and was selected as the later control model. Forty distributed elevons and LQG/LTR controller were employed to control the attitude and suppress the aeroelastic deformation of the UAV simultaneously. The dynamics performance, robustness, and simulation results show that they were suitable for large scale highly flexible solar-powered UAV.


---

## Body

## 1. Introduction

The wing of high aspect-ratio solar-powered UAV is highly flexible and the surface density of structure is much lower than other kinds of aircraft, which may lead to the huge scale of solar-powered UAV in order to carry sufficient payload. These characteristics will lead to significant geometrically nonlinear aeroelastic deformation during flight and strong coupling between aeroelasticity and flight dynamics. If these issues have not been taken into full account, large error will occur. One of the key recommendations from the flight mishap investigation of the most famous solar-powered UAV “Helios” was to “develop more advanced, multidisciplinary (structures, aeroelastic, aerodynamics, atmospheric, materials, propulsion, controls, etc) “time-domain” analysis methods appropriate to highly flexible, “morphing” vehicles” [1].In the modeling of “morphing” vehicles, fruitful works have been presented by Hodges, Cesnik, and their coworkers. The analysis method of Hodges et al. underlying NATASHA (Nonlinear Aeroelastic Trim And Stability of HALE Aircraft), a computer program that is based on geometrically exact, fully intrinsic beam equations and a finite-state induced flow model [2–4]. The University of Michigan’s Nonlinear Aeroelastic Simulation Toolbox (UM/NAST) is developed by Cesnik and his coworkers [5, 6], which includes nonlinear strain-based beam model, unsteady aerodynamics with simplified stall models, and a six-degrees-of-freedom flight dynamic formulation.In the study of flight dynamics of highly flexible UAV, [2] found that the behavior is distinctly different from that of a rigid aircraft while taking into account the flexibility effects. Because the low-frequency modes of structure are completely coupled to the flight dynamics mode, the phugoid as well as the short period mode is affected by wing flexibility significantly. Reference [3] furthered the study, the paper presented a theoretical basis for the flight dynamic response estimation of a highly flexible flying wing, multiple engines, and multiple control surfaces are taken into account. The results show that the phugoid mode may become unstable during large dihedral angles caused by large mass of payload pod. The pair of complex-conjugate short period roots merges to become two real roots. This is because the deformed “U”-shape leads to an order-of-magnitude increase in the pitch moment of inertia and so there is corresponding decrease in the frequency. Thus, this highly flexible aircraft does not show a classical short-period mode in its deformed state. Similar studies by Raghavan and Patil [4] and Su and Cesnik [5] confirmed this result. All the references paid attention to the nonlinear model and the results analysis, but the modeling of linear dynamics was ignored.In [7], an optimal Static Output Feedback controller was designed for flutter suppression and gust load alleviation. However, only an elevon on the wing tip was employed in the controller, and the performance was highly dependent on sensor selection and placement.In [8, 9], dynamic inversion control laws were presented for the trajectory control for Very Flexible Aircraft (VFA). However, in the use of dynamic inversion, there is an assumption that the stability and control derivatives remain fairly constant during the flight. This is not the case with highly flexible UAV.In [10], a reduced-order H∞ controller was employed to control a high speed civil transport aircraft, and the results implied that the reduced-order H∞ controller was robust in the presence of unstructured uncertainty arising due to the model reduction error. However, the reduced order controller is weakly robust when employed to control the full-order model with parametric uncertainty.LQG/LTR is another popular robust control method. The primary objective of LQG/LTR is to recover the robustness of LQG, an excellent full-state feedback robust controller. The procedure of LQG/LTR method is as follows: firstly, design a target feedback loop, which satisfies desired frequency performance, such as disturbance rejection and uncertainty deviation of the unmodeled dynamics, and then design a model based compensator to recover the properties of the target feedback loop [11]. LQG/LTR has been used widely in flight control, but only a few references are concentrated on the highly flexible aircraft.In [12], an adaptive LQG/LTR controller was employed to control a very flexible aircraft, but the dynamics model was greatly simplified to capture only the unstable phugoid mode, and there were only 3 traditional control effectors, throttle, elevator, and aileron, available for the controller.The large scale highly flexible UAV discussed in this paper is shown in Figure1. It is equipped with 8 propulsion system nacelles (including secondary batteries, electrical motor, and propeller) and 5 vertical tails. The span is 70 m and the mean aerodynamics chord is 2.44 m. In order to enhance the control ability, trailing-edge elevon spans the entire wing. Similar to the Helios [1], the elevon is divided into 40 distributed parts to release the notable aeroelastic deformation caused by large scale highly flexible wing.Figure 1
High aspect-ratio highly flexible solar-powered UAV.For the aeroelastic flight dynamics analysis of such a UAV, most papers ignored the aerodynamics force distribution characteristic along the span and assumed the center of gravity (CG), elastic center (EC), and aerodynamics center (AC) to be coincident for simplification. In this paper, a more practical UAV model is proposed, in which the AC and EC is separated, the 3D quasisteady aerodynamics force is considered, and aerodynamics derivatives in every wing sections are derived for accuracy. In order to get a more accurate characteristic of solar-powered UAV, the mass parameters derived from conceptual design are employed, shown in Table1.Table 1
The mass parameters of high aspect-ratio solar-powered UAV.
Coefficient of mass
CG (from leading edge)
Elastic center
Structure
0.3
35%
35%
Solar array
0.15
45%
/
Propulsion system
0.5
0%
/
Payload
0.05
(at the CG of UAV)
/In this paper, a nonlinear coupled structural and flight dynamics model based on geometrically nonlinear intrinsic beam was established first, and then its explicit linear dynamics model is derived for decoupled dynamics modal analysis by mean axis theory. At last, a LQG/LTR controller is designed with distributed 40 elevons for attitude control and aeroelastic deformation suppression for the large-scale highly flexible solar-powered UAV.


## 2. Modeling of Nonlinear Structural/Flight Dynamics

### 2.1. Structural Dynamics Model

The geometrically nonlinear structural dynamics and flight dynamics equations are based on the intrinsic beam equation derived by Hodges et al. shown as follows [3, 13]:(1)P-˙+Ω-~P-=F′+K-~F-+f-,H-˙+Ω-~H-+V-~P-=M′+K-~M-+e~2+γ-~F-+m-,γ-˙=V′+K-~V-+e~2+γ-~Ω-,κ-˙=Ω′+K-~Ω-,P^˙+Ω^~P^=f^,H^˙+Ω^~H^=m^.The equations above are the core of structural dynamics, but they still cannot be used for the analysis of aerodynamics and flight dynamics without the relationship between displacements and rotations. Hence, the kinematics motion must be established to obtain the relationship between strain and displacement. The position vector of the origin of the deformed beam frame (B-) in the root frame (R-) and their rotational relationship are [14]:(2)rBR=CRBγ+e2,CBR′=-κ~+k~CBR.The spatial derivatives can be replaced by the following equations, namely, central difference method [15]:(3)xi′=12hxi+1-xi-1,xi′′=1h2xi+1-2xi+xi-1,xi′′′=12h3xi+2-2xi+1+2xi-1-xi-2,xi′′′′=1h4xi+2-4xi+1+6xi-4xi-1+xi-2,where h is the distance between two points.The mean axis is applied to simplify and linearize the equation of motion. In the mean axis (M-) frame, the linear and angular momentums due to elastic deformation are zero at every instant, which eliminates the coupling between the rigid body modes and the structural mode and results in a reduced-order system with only six degrees-of-freedom.Since the equations in (1) are written in theB-frame, the equations of motion of rigid mode can be derived by premultiplying all terms by CMB to obtain components in theM- frame and integrating in it along the beam axis as follows [16]:(4)MtotV˙cgM+ΩMM×VcgM=∑CMBifnBi,(5)IcgMΩ˙MM+I˙cgMΩMM+ΩMM×IcgMΩMM=∑rM×CMBifnBi+∑CMBimiBi-∑CMBnIiBΩ˙iB-∑C˙MBnIiBΩiB,where the linear and angular velocities are defined as(6)VcgM=1Mtot∑μiCMBiVi,ΩMM=IcgM-1∑rBiMM×μiCMBiVi.The relationship of coordinate transformation is(7)CMB=CMgCgB.The total moment of inertia of the UAV is solved as follows:(8)IcgM=∑CMBiIiCMBiT+riTriE-ririTmi.If the relative position of the rigid body does not vary significantly,{I˙cg}M can be approximated to be zero. If the motion of the UAV is limited to longitudinal only, the last two terms at the right hand of (5) can also be equal to zero.


### 2.2. Force and Moment Model

#### 2.2.1. Internal Forces and Moments Model

The wing of solar-powered UAV may undergo large deformations during normal flight, while the strain is small, exhibiting geometrically nonlinear behaviors. Therefore, the small-strain assumption is suitable for this problem. The internal forces and moments are linearly related to the strains and curvatures, shown in (9) [17]:(9)FM=RSSTT-1γκ.The wing is assumed to be a prismatic and isotropic beam with the shear center and tension axis coincident with the reference axis, by which the above equation can be derived as follows:(10)FxFyFz=GAxEAGAzγxγyγz,MxMyMz=EIxGJEIzκxκyκz,where G represents the shear modulus; E represents Young’s modulus; J represents the torsional constant; Ix and Iz represent the area moment of inertia along the x and z direction.


#### 2.2.2. Aerodynamics Model

In order to analyze the motion of the flexible UAV exactly, the aerodynamic force of vibrating-wings should be analyzed. Unsteady aerodynamics models that are useful for aeroelastic analysis are usually Theodorsen unsteady aerodynamic model and finite state induced flow model of Peters and Cao [18]. The Theodorsen model is derived in frequency domain, but it can be converted into time domain by rational function approximation. The finite state induced flow model is derived in time domain, and it can be used in a state-space model easily. These two types of aerodynamics are equivalent and can be derived from one another through Laplace and Fourier transforms.In this paper, the Theodorsen unsteady aerodynamics model is adopted, which has been used widely for airfoil or large aspect ratio wings, and its quasisteady form is relative to the traditional aerodynamics derivatives used in rigid aircraft.For the Theodorsen model, the wing section is described in Figure2, where the definitions of points are as follows: LE: leading edge, G: center of gravity of the wing section,  g: center of gravity of the UAV, A: aerodynamics center, M: middle point of chord, E: elastic center, the position in x direction of s-frame is a-b, C: control point of unsteady aerodynamics, and D: aerodynamics center of control surface.Figure 2
Illustration of wing section, coordinate systems, and variables.The definitions of variables are as follows:α: angle of attack,  b: length of 1/2 chord, u: deformation defined in B-frame, and s-frame: the origin is at LE, x along the chord and pointing to the trailing edge, z perpendicular to x and point up.The formulation of Theodorsen unsteady aerodynamics force and its moment to the elastic center are(11)CL=πbα˙V+h¨V2-ba-α¨V2+2πCkh˙V+α+b12-a-α˙V,CM,E=π2ba-h¨V2-b2V218+a-2α¨-12-a-bα˙V+πa-+12Ckh˙V+α+b12-a-α˙V=πbα˙V-14-12a-+h¨V212a--ba-α¨V218+a-212b+14+12a-2πCk·h˙V+α+b12-a-α˙V,where 2π is the theoretic lift slope of airfoil, the steady aerodynamics center is at 1/4 chord and unsteady aerodynamics center and its control point is at 3/4 chord. All these variables can be replaced by the actual airfoil data. C(k) is the Theodorsen function, produced by circulation. k=ωb/V is the reduced frequency. If k is low enough and C(k)=0, the flow will become quasisteady. As the flight dynamics and control of solar-powered UAV on low frequency is more concerned, quasisteady flow assumption is adopted for analysis in this paper.By applying the quasisteady flow assumption and ignoring the terms in high orders, (11) becomes(12)CL=CLαα+u˙zV+1-a-bα˙V,CM,E=CLα14+12a-α+u˙zV+CLα14a--12a-2bα˙V,where the airfoil data have been replaced by the actual ones.The physical meaning of the above equations is that the lift of wing section is produced by steady and quasisteady flow together, and the steady part is produced by the angle of attackα of the section and its increment u˙z/V0 caused by relative elastic deformation. The quasisteady part is caused by α˙.It should be noted that, in most cases, the angles of attack in every sections (denoted asαB) are not easy to be measured directly. Hence, the relationship between αB and incident angle at the CG of UAV expressed in aM-frame (notated as αM and βM) can be calculated as(13)-sin⁡αB=CaBB3,1=0,0,1·CaBaMCaMMCMB1,0,0T=cos⁡⁡αMsin⁡βMsin⁡φBM,x-sin⁡αMcos⁡⁡φBM,x·cos⁡⁡φBM,ycos⁡⁡φBM,z+cos⁡⁡βMsin⁡φBM,x·cos⁡⁡φBM,ysin⁡φBM,z-sin⁡αMsin⁡βMsin⁡φBM,x+cos⁡⁡αMcos⁡⁡φBM,x·sin⁡φBM,y,where φBM is the Euler angle from B-frame to M-frame and the subscripts x, y, and z represent the 3 components of Euler angle, which correspond to the roll, pitch, and yaw angle of the aircraft, respectively.If the wing rotates along thex axis only, (13) can be simplified as(14)sin⁡αB=-cos⁡⁡αMsin⁡βMsin⁡φBM,x+sin⁡αMcos⁡⁡φBM,x.The magnitude of drag is small, and its influence on the flight dynamics characteristics is limited, so the drag on every wing section can be divided from the total drag. The expression of total drag is(15)CD=CD0+ainduceCL2.The aerodynamics force expressed inaM-frame can be derived by the following axis conversion:(16)faaM=CaMaBfaaB,where(17)CaMaB=1000cos⁡⁡φBM,x-sin⁡φBM,x0sin⁡φBM,xcos⁡⁡φBM,x.Substituting (17) into (16), the nondimensional force coefficients at the aM-coordinate are derived as follows:(18)CDaM=CDaB,CCaM=CLaBsin⁡φBM,x,CLaM=CLaBcos⁡⁡φBM,x.Changing the reference point of moment into CG, the aerodynamics moment to CG of UAV is(19)maM=ma,0M+rM×CMaMfaaM.It should be noted that the method proposed in this section can be used to calculate either the aerodynamics force or the moment of vertical tails, although the incident angle and Euler angle may vary significantly.


#### 2.2.3. 3D Effect Correction of Aerodynamics

In the analysis of the aeroelastic flight dynamics of solar-powered UAV, most papers ignored the aerodynamics force distribution characteristics along the wing span and considered the flow condition to be same as the airfoils. Actually, the wing span is finite, which will make the flow of each section affect each other; that is to say, 3D effect must be considered for precision. In this paper, 3D effect is corrected based on lifting-line theory, and it applies well to wings that have aspect ratio greater than 5 and sweep angle less than 20 degrees.According to lifting-line theory, the lift of airfoil is determined by effective angle of attackαeff, which is formulated as follows:(20)αeff=αgeo-αind,where αgeo and αind represent geometric and induced angle of attack, respectively.According to the relationship of trailing vortex and circulation, the equation above becomes(21)αgeoy0=Γy0πV∞cy0+αL=0+14πV∞∫0ldΓ/dydyy0-y,where y0 is an arbitrary spanwise position on the wing and 2π refers to theoretical lift slope of airfoil, which can be substituted by true airfoil value in the actual computation. Since (21) is a nonlinear integral-differential equation, coordinate transformation method should be applied before solving it for simplification [19].After solving the circulation distributionΓ(y) from (21), the lift distribution can be obtained from the Kutta-Joukowski theorem:(22)Ly0=ρV∞Γy0.By repeating the above steps for different geometric angles of attack and taking a derivative, one can obtain the lift slope of each cross section of the wing considering 3D effect correction.


## 2.1. Structural Dynamics Model

The geometrically nonlinear structural dynamics and flight dynamics equations are based on the intrinsic beam equation derived by Hodges et al. shown as follows [3, 13]:(1)P-˙+Ω-~P-=F′+K-~F-+f-,H-˙+Ω-~H-+V-~P-=M′+K-~M-+e~2+γ-~F-+m-,γ-˙=V′+K-~V-+e~2+γ-~Ω-,κ-˙=Ω′+K-~Ω-,P^˙+Ω^~P^=f^,H^˙+Ω^~H^=m^.The equations above are the core of structural dynamics, but they still cannot be used for the analysis of aerodynamics and flight dynamics without the relationship between displacements and rotations. Hence, the kinematics motion must be established to obtain the relationship between strain and displacement. The position vector of the origin of the deformed beam frame (B-) in the root frame (R-) and their rotational relationship are [14]:(2)rBR=CRBγ+e2,CBR′=-κ~+k~CBR.The spatial derivatives can be replaced by the following equations, namely, central difference method [15]:(3)xi′=12hxi+1-xi-1,xi′′=1h2xi+1-2xi+xi-1,xi′′′=12h3xi+2-2xi+1+2xi-1-xi-2,xi′′′′=1h4xi+2-4xi+1+6xi-4xi-1+xi-2,where h is the distance between two points.The mean axis is applied to simplify and linearize the equation of motion. In the mean axis (M-) frame, the linear and angular momentums due to elastic deformation are zero at every instant, which eliminates the coupling between the rigid body modes and the structural mode and results in a reduced-order system with only six degrees-of-freedom.Since the equations in (1) are written in theB-frame, the equations of motion of rigid mode can be derived by premultiplying all terms by CMB to obtain components in theM- frame and integrating in it along the beam axis as follows [16]:(4)MtotV˙cgM+ΩMM×VcgM=∑CMBifnBi,(5)IcgMΩ˙MM+I˙cgMΩMM+ΩMM×IcgMΩMM=∑rM×CMBifnBi+∑CMBimiBi-∑CMBnIiBΩ˙iB-∑C˙MBnIiBΩiB,where the linear and angular velocities are defined as(6)VcgM=1Mtot∑μiCMBiVi,ΩMM=IcgM-1∑rBiMM×μiCMBiVi.The relationship of coordinate transformation is(7)CMB=CMgCgB.The total moment of inertia of the UAV is solved as follows:(8)IcgM=∑CMBiIiCMBiT+riTriE-ririTmi.If the relative position of the rigid body does not vary significantly,{I˙cg}M can be approximated to be zero. If the motion of the UAV is limited to longitudinal only, the last two terms at the right hand of (5) can also be equal to zero.


## 2.2. Force and Moment Model

### 2.2.1. Internal Forces and Moments Model

The wing of solar-powered UAV may undergo large deformations during normal flight, while the strain is small, exhibiting geometrically nonlinear behaviors. Therefore, the small-strain assumption is suitable for this problem. The internal forces and moments are linearly related to the strains and curvatures, shown in (9) [17]:(9)FM=RSSTT-1γκ.The wing is assumed to be a prismatic and isotropic beam with the shear center and tension axis coincident with the reference axis, by which the above equation can be derived as follows:(10)FxFyFz=GAxEAGAzγxγyγz,MxMyMz=EIxGJEIzκxκyκz,where G represents the shear modulus; E represents Young’s modulus; J represents the torsional constant; Ix and Iz represent the area moment of inertia along the x and z direction.


### 2.2.2. Aerodynamics Model

In order to analyze the motion of the flexible UAV exactly, the aerodynamic force of vibrating-wings should be analyzed. Unsteady aerodynamics models that are useful for aeroelastic analysis are usually Theodorsen unsteady aerodynamic model and finite state induced flow model of Peters and Cao [18]. The Theodorsen model is derived in frequency domain, but it can be converted into time domain by rational function approximation. The finite state induced flow model is derived in time domain, and it can be used in a state-space model easily. These two types of aerodynamics are equivalent and can be derived from one another through Laplace and Fourier transforms.In this paper, the Theodorsen unsteady aerodynamics model is adopted, which has been used widely for airfoil or large aspect ratio wings, and its quasisteady form is relative to the traditional aerodynamics derivatives used in rigid aircraft.For the Theodorsen model, the wing section is described in Figure2, where the definitions of points are as follows: LE: leading edge, G: center of gravity of the wing section,  g: center of gravity of the UAV, A: aerodynamics center, M: middle point of chord, E: elastic center, the position in x direction of s-frame is a-b, C: control point of unsteady aerodynamics, and D: aerodynamics center of control surface.Figure 2
Illustration of wing section, coordinate systems, and variables.The definitions of variables are as follows:α: angle of attack,  b: length of 1/2 chord, u: deformation defined in B-frame, and s-frame: the origin is at LE, x along the chord and pointing to the trailing edge, z perpendicular to x and point up.The formulation of Theodorsen unsteady aerodynamics force and its moment to the elastic center are(11)CL=πbα˙V+h¨V2-ba-α¨V2+2πCkh˙V+α+b12-a-α˙V,CM,E=π2ba-h¨V2-b2V218+a-2α¨-12-a-bα˙V+πa-+12Ckh˙V+α+b12-a-α˙V=πbα˙V-14-12a-+h¨V212a--ba-α¨V218+a-212b+14+12a-2πCk·h˙V+α+b12-a-α˙V,where 2π is the theoretic lift slope of airfoil, the steady aerodynamics center is at 1/4 chord and unsteady aerodynamics center and its control point is at 3/4 chord. All these variables can be replaced by the actual airfoil data. C(k) is the Theodorsen function, produced by circulation. k=ωb/V is the reduced frequency. If k is low enough and C(k)=0, the flow will become quasisteady. As the flight dynamics and control of solar-powered UAV on low frequency is more concerned, quasisteady flow assumption is adopted for analysis in this paper.By applying the quasisteady flow assumption and ignoring the terms in high orders, (11) becomes(12)CL=CLαα+u˙zV+1-a-bα˙V,CM,E=CLα14+12a-α+u˙zV+CLα14a--12a-2bα˙V,where the airfoil data have been replaced by the actual ones.The physical meaning of the above equations is that the lift of wing section is produced by steady and quasisteady flow together, and the steady part is produced by the angle of attackα of the section and its increment u˙z/V0 caused by relative elastic deformation. The quasisteady part is caused by α˙.It should be noted that, in most cases, the angles of attack in every sections (denoted asαB) are not easy to be measured directly. Hence, the relationship between αB and incident angle at the CG of UAV expressed in aM-frame (notated as αM and βM) can be calculated as(13)-sin⁡αB=CaBB3,1=0,0,1·CaBaMCaMMCMB1,0,0T=cos⁡⁡αMsin⁡βMsin⁡φBM,x-sin⁡αMcos⁡⁡φBM,x·cos⁡⁡φBM,ycos⁡⁡φBM,z+cos⁡⁡βMsin⁡φBM,x·cos⁡⁡φBM,ysin⁡φBM,z-sin⁡αMsin⁡βMsin⁡φBM,x+cos⁡⁡αMcos⁡⁡φBM,x·sin⁡φBM,y,where φBM is the Euler angle from B-frame to M-frame and the subscripts x, y, and z represent the 3 components of Euler angle, which correspond to the roll, pitch, and yaw angle of the aircraft, respectively.If the wing rotates along thex axis only, (13) can be simplified as(14)sin⁡αB=-cos⁡⁡αMsin⁡βMsin⁡φBM,x+sin⁡αMcos⁡⁡φBM,x.The magnitude of drag is small, and its influence on the flight dynamics characteristics is limited, so the drag on every wing section can be divided from the total drag. The expression of total drag is(15)CD=CD0+ainduceCL2.The aerodynamics force expressed inaM-frame can be derived by the following axis conversion:(16)faaM=CaMaBfaaB,where(17)CaMaB=1000cos⁡⁡φBM,x-sin⁡φBM,x0sin⁡φBM,xcos⁡⁡φBM,x.Substituting (17) into (16), the nondimensional force coefficients at the aM-coordinate are derived as follows:(18)CDaM=CDaB,CCaM=CLaBsin⁡φBM,x,CLaM=CLaBcos⁡⁡φBM,x.Changing the reference point of moment into CG, the aerodynamics moment to CG of UAV is(19)maM=ma,0M+rM×CMaMfaaM.It should be noted that the method proposed in this section can be used to calculate either the aerodynamics force or the moment of vertical tails, although the incident angle and Euler angle may vary significantly.


### 2.2.3. 3D Effect Correction of Aerodynamics

In the analysis of the aeroelastic flight dynamics of solar-powered UAV, most papers ignored the aerodynamics force distribution characteristics along the wing span and considered the flow condition to be same as the airfoils. Actually, the wing span is finite, which will make the flow of each section affect each other; that is to say, 3D effect must be considered for precision. In this paper, 3D effect is corrected based on lifting-line theory, and it applies well to wings that have aspect ratio greater than 5 and sweep angle less than 20 degrees.According to lifting-line theory, the lift of airfoil is determined by effective angle of attackαeff, which is formulated as follows:(20)αeff=αgeo-αind,where αgeo and αind represent geometric and induced angle of attack, respectively.According to the relationship of trailing vortex and circulation, the equation above becomes(21)αgeoy0=Γy0πV∞cy0+αL=0+14πV∞∫0ldΓ/dydyy0-y,where y0 is an arbitrary spanwise position on the wing and 2π refers to theoretical lift slope of airfoil, which can be substituted by true airfoil value in the actual computation. Since (21) is a nonlinear integral-differential equation, coordinate transformation method should be applied before solving it for simplification [19].After solving the circulation distributionΓ(y) from (21), the lift distribution can be obtained from the Kutta-Joukowski theorem:(22)Ly0=ρV∞Γy0.By repeating the above steps for different geometric angles of attack and taking a derivative, one can obtain the lift slope of each cross section of the wing considering 3D effect correction.


## 2.2.1. Internal Forces and Moments Model

The wing of solar-powered UAV may undergo large deformations during normal flight, while the strain is small, exhibiting geometrically nonlinear behaviors. Therefore, the small-strain assumption is suitable for this problem. The internal forces and moments are linearly related to the strains and curvatures, shown in (9) [17]:(9)FM=RSSTT-1γκ.The wing is assumed to be a prismatic and isotropic beam with the shear center and tension axis coincident with the reference axis, by which the above equation can be derived as follows:(10)FxFyFz=GAxEAGAzγxγyγz,MxMyMz=EIxGJEIzκxκyκz,where G represents the shear modulus; E represents Young’s modulus; J represents the torsional constant; Ix and Iz represent the area moment of inertia along the x and z direction.


## 2.2.2. Aerodynamics Model

In order to analyze the motion of the flexible UAV exactly, the aerodynamic force of vibrating-wings should be analyzed. Unsteady aerodynamics models that are useful for aeroelastic analysis are usually Theodorsen unsteady aerodynamic model and finite state induced flow model of Peters and Cao [18]. The Theodorsen model is derived in frequency domain, but it can be converted into time domain by rational function approximation. The finite state induced flow model is derived in time domain, and it can be used in a state-space model easily. These two types of aerodynamics are equivalent and can be derived from one another through Laplace and Fourier transforms.In this paper, the Theodorsen unsteady aerodynamics model is adopted, which has been used widely for airfoil or large aspect ratio wings, and its quasisteady form is relative to the traditional aerodynamics derivatives used in rigid aircraft.For the Theodorsen model, the wing section is described in Figure2, where the definitions of points are as follows: LE: leading edge, G: center of gravity of the wing section,  g: center of gravity of the UAV, A: aerodynamics center, M: middle point of chord, E: elastic center, the position in x direction of s-frame is a-b, C: control point of unsteady aerodynamics, and D: aerodynamics center of control surface.Figure 2
Illustration of wing section, coordinate systems, and variables.The definitions of variables are as follows:α: angle of attack,  b: length of 1/2 chord, u: deformation defined in B-frame, and s-frame: the origin is at LE, x along the chord and pointing to the trailing edge, z perpendicular to x and point up.The formulation of Theodorsen unsteady aerodynamics force and its moment to the elastic center are(11)CL=πbα˙V+h¨V2-ba-α¨V2+2πCkh˙V+α+b12-a-α˙V,CM,E=π2ba-h¨V2-b2V218+a-2α¨-12-a-bα˙V+πa-+12Ckh˙V+α+b12-a-α˙V=πbα˙V-14-12a-+h¨V212a--ba-α¨V218+a-212b+14+12a-2πCk·h˙V+α+b12-a-α˙V,where 2π is the theoretic lift slope of airfoil, the steady aerodynamics center is at 1/4 chord and unsteady aerodynamics center and its control point is at 3/4 chord. All these variables can be replaced by the actual airfoil data. C(k) is the Theodorsen function, produced by circulation. k=ωb/V is the reduced frequency. If k is low enough and C(k)=0, the flow will become quasisteady. As the flight dynamics and control of solar-powered UAV on low frequency is more concerned, quasisteady flow assumption is adopted for analysis in this paper.By applying the quasisteady flow assumption and ignoring the terms in high orders, (11) becomes(12)CL=CLαα+u˙zV+1-a-bα˙V,CM,E=CLα14+12a-α+u˙zV+CLα14a--12a-2bα˙V,where the airfoil data have been replaced by the actual ones.The physical meaning of the above equations is that the lift of wing section is produced by steady and quasisteady flow together, and the steady part is produced by the angle of attackα of the section and its increment u˙z/V0 caused by relative elastic deformation. The quasisteady part is caused by α˙.It should be noted that, in most cases, the angles of attack in every sections (denoted asαB) are not easy to be measured directly. Hence, the relationship between αB and incident angle at the CG of UAV expressed in aM-frame (notated as αM and βM) can be calculated as(13)-sin⁡αB=CaBB3,1=0,0,1·CaBaMCaMMCMB1,0,0T=cos⁡⁡αMsin⁡βMsin⁡φBM,x-sin⁡αMcos⁡⁡φBM,x·cos⁡⁡φBM,ycos⁡⁡φBM,z+cos⁡⁡βMsin⁡φBM,x·cos⁡⁡φBM,ysin⁡φBM,z-sin⁡αMsin⁡βMsin⁡φBM,x+cos⁡⁡αMcos⁡⁡φBM,x·sin⁡φBM,y,where φBM is the Euler angle from B-frame to M-frame and the subscripts x, y, and z represent the 3 components of Euler angle, which correspond to the roll, pitch, and yaw angle of the aircraft, respectively.If the wing rotates along thex axis only, (13) can be simplified as(14)sin⁡αB=-cos⁡⁡αMsin⁡βMsin⁡φBM,x+sin⁡αMcos⁡⁡φBM,x.The magnitude of drag is small, and its influence on the flight dynamics characteristics is limited, so the drag on every wing section can be divided from the total drag. The expression of total drag is(15)CD=CD0+ainduceCL2.The aerodynamics force expressed inaM-frame can be derived by the following axis conversion:(16)faaM=CaMaBfaaB,where(17)CaMaB=1000cos⁡⁡φBM,x-sin⁡φBM,x0sin⁡φBM,xcos⁡⁡φBM,x.Substituting (17) into (16), the nondimensional force coefficients at the aM-coordinate are derived as follows:(18)CDaM=CDaB,CCaM=CLaBsin⁡φBM,x,CLaM=CLaBcos⁡⁡φBM,x.Changing the reference point of moment into CG, the aerodynamics moment to CG of UAV is(19)maM=ma,0M+rM×CMaMfaaM.It should be noted that the method proposed in this section can be used to calculate either the aerodynamics force or the moment of vertical tails, although the incident angle and Euler angle may vary significantly.


## 2.2.3. 3D Effect Correction of Aerodynamics

In the analysis of the aeroelastic flight dynamics of solar-powered UAV, most papers ignored the aerodynamics force distribution characteristics along the wing span and considered the flow condition to be same as the airfoils. Actually, the wing span is finite, which will make the flow of each section affect each other; that is to say, 3D effect must be considered for precision. In this paper, 3D effect is corrected based on lifting-line theory, and it applies well to wings that have aspect ratio greater than 5 and sweep angle less than 20 degrees.According to lifting-line theory, the lift of airfoil is determined by effective angle of attackαeff, which is formulated as follows:(20)αeff=αgeo-αind,where αgeo and αind represent geometric and induced angle of attack, respectively.According to the relationship of trailing vortex and circulation, the equation above becomes(21)αgeoy0=Γy0πV∞cy0+αL=0+14πV∞∫0ldΓ/dydyy0-y,where y0 is an arbitrary spanwise position on the wing and 2π refers to theoretical lift slope of airfoil, which can be substituted by true airfoil value in the actual computation. Since (21) is a nonlinear integral-differential equation, coordinate transformation method should be applied before solving it for simplification [19].After solving the circulation distributionΓ(y) from (21), the lift distribution can be obtained from the Kutta-Joukowski theorem:(22)Ly0=ρV∞Γy0.By repeating the above steps for different geometric angles of attack and taking a derivative, one can obtain the lift slope of each cross section of the wing considering 3D effect correction.


## 3. Model Linearization

### 3.1. Linearized Aerodynamics

The movement of wing section is separated into two parts: the rigid movement of UAV and the elastic movement of structure. Then, the lift coefficients applied on the wing section can be deployed in theaB-frame as follows:(23)CLα,CLq=0.5CLα+2CLα34-x-gs,CLp=-2CLαygsl,CLφy=CLα,CLu˙z=CLα1V,CLφ˙y=CLα1-a-c2V.It should be noted that the variables in the above equations are measured in theaB-frame and B-frame.According to the 2nd equation of (12), the aerodynamics moment to the elastic center (EC) is(24)Cmα,E=-CLαx-AEB≈CLα14+12a-,Cmq,E=0.5CLαx-CEB+2CLαx-Cgsx-AEB,Cmp,E=-2CLαygslx-AEB,Cmφy,E=Cmα,E,Cmu˙z,E=Cmα,E1V,Cmφ˙y,E=CLα14a--12a-2c2V.From the 2nd equation of (18), the difference of lift expressed in the aM-frame and aB-frame is only cos⁡⁡φBM,x. Hence, by replacing CLα in (23) with CLαcos⁡⁡φBM,x, the lift expressed in the aB-frame will be derived.Substituting the expression of lift ataM-frame into (19), the pitch moment about the center of gravity of the UAV can be obtained as follows:(25)Cmα,g=-CLαcos⁡⁡φMB,xx-As-x-gs,Cmq,g=CLαcos⁡⁡φMB,x0.5x-CgB-2x-CgB2,Cmp,g=-2CLαygslx-CgB,Cmφy,g=Cmα,g,Cmu˙z,g=Cmα,g1V,Cmφ˙y,g=-2CLαcos⁡⁡φMB,x12-x-gs2bV.


### 3.2. Linearized Flight Dynamics Equation

It can be found that the form of (4) and (5) is similar to that of flight dynamics equations of normal rigid vehicles. Hence, they can also be linearized by small disturbance method and written into two sets of equations: longitudinal and lateral-directional. Different from the normal rigid vehicle, the external force such as aerodynamics, thrust, and gravity are related to rigid motion and elastic motion for elastic model. The linear equation of flight dynamics for elastic UAV is (26)x˙r=Arxr+Brur+Arexe,where xr is the independent variable of rigid motion, with xr=[v,α,q,θ]T in longitudinal equations; xe=[u˙xT,u˙zT,φ˙yT,uxT,uzT,φyT]T is the independent variable of elastic motion; ur is the control variable of rigid motion, ur=[δe,1⋯δe,40,δt]T; that is, all the elevons can be controlled independently, but all throttles act conformably again for simplification; Ar is the same as conventional rigid vehicles as shown in [20]; and Br is similar to the rigid model also, but the column is equal to the length of ur; Are is the influence matrix of elastic motion to rigid motion, for longitudinal motion:(27)Are=01×n01×n01×n01×n01×n01×n01×nZu˙z,1⋯Zu˙z,nZφ˙y,1⋯Zφ˙y,n01×n01×nZφy,1⋯Zφy,n01×nMu˙z,1⋯Mu˙z,nMφ˙y,1⋯Mφ˙y,n01×n01×nMφy,1⋯Mφy,n01×n01×n01×n01×n01×n01×n.The details ofAre can be found in [21].


### 3.3. Linearized Structural Dynamics Equation

Ignoring nonlinear coupling terms in (1), and assuming that the translational strain is small and the external moment and the moment of inertia of x and z direction are equal to zero, (1) can be reduced into(28)μV˙x=EIzκz′′+fx,μV˙z=-EIxκx′′+fz,IyyΩ˙y=EIyκy′+mz,κ˙x=Vz′′,κ˙y=Ωy′,κ˙z=-Vx′,μV˙y=fy.The equations above can be rearranged into a compact form as follows:(29)μx¨+EIzx′′′′=fx,μz¨+EIxz′′′′=fz,Iyyφ¨Bg,y-GJzφBg,y′′=my,μV˙y=fy.Since the focus of this paper is the flight dynamics characteristics of the total UAV, the structural dynamics characteristics are used to derive the increment of aerodynamics force caused by aeroelasticity only and many assumptions can be introduced to linearize the equations above for simplification.Equation (29) can be simplified to the form as in (30) by rewriting in M-frame and assuming that the structural elastic vibrations are small perturbations near the static trim state and the elastic motion in x direction is ignorable. In the case of geometry, it is assumed that the M-frame is parallel to R-frame and the sweep angle of wing is zero. Besides, in mean axis system, the elastic motion is influenced by rigid motion through aerodynamics force only:(30)μu¨xM,x+EIzuxM,x′′′′=cos⁡⁡φBM,zfB,x,μu¨zM,z+EIxuzM,z′′′′=cos⁡⁡φBM,xfB,z,Iyyφ¨BM,y-GJzφBM,y′′=mB,y,where φBM,y represents the angle between x-axis of B-frame and the oxy plane of M-frame. All variables in (30) are in incremental form and the notation Δ has been neglected for simplification.Replacing the spatial derivatives by (3), and knowing the fact that the wing tip is free, the boundary conditions can be defined that the bending moments and shear forces at the first and last nodes are zero for the 2nd equation in (30); that is [15],(31)uz,1′′=0,uz,n′′=0,uz,1′′′=0,(32)uz,n′′′=0.Thus, the structural dynamics equation can be derived as follows: (33)Muzu¨z+Kuzuz=Quz,where(34)uz=u1⋯unT,Muz=μ1⋱μn,Kuz=2-42-25-411-46-41⋱⋱⋱⋱⋱1-46-411-45-22-42∗1Δl4diag⁡EIz,i.For the 2nd equation in (30), it is the same as (33); that is,(35)Muxu¨x+Kuxux=Qux.For the 3rd equation in (30), the boundary conditions are [15](36)φBM,y,1′=0,φBM,y,n′=0.Then, the structural dynamics equation becomes(37)Mφyφ¨y+Kφyφy=Qφy,where(38)φy=φBM,y,1⋯φBM,y,nT,Mφy=I1⋱In,Kφy=--221-21⋱⋱⋱1-212-2∗1Δl2diag⁡GJi.In order to further simplify the linear structural dynamics equations, (33), (35), and (37) can be expanded in the modal coordinate system, and the final form can be derived as follows:(39)x˙ne=Anexne+Anerxr+Bneur,where the independent variables are(40)xne=x˙nxn,xe=μxn,xe=uxT,uzT,φyTT,where xn is the modal displacements of structural elastic motion, which is the symmetric modes for longitudinal motion; μ is the transfer matrix from modal coordinate to physical coordinate.For longitudinal motion, the coefficient matrices in (39) are(41)Ane=Mn-1Qn1Mn-1Qn2-KnI0,Aner=Mn-1Qn00,Bne=Mn-1Qn30,Mn=μTMμ,Kn=μTKμ,Qni=μTQi,Qn0=μTQ0,Qn1=μTQ1μ,Qn2=μTQ2μ,Qn3=μTQ3μ,Q0=Xv,1⋮Xv,nXα,1⋮Xα,n0n×10n×1Zv,1⋮Zv,nZα,1⋮Zα,nZq,1⋮Zq,n0n×1Mv,1⋮Mv,nMα,1⋮Mα,nMq,1⋮Mq,n0n×1,Q1=Xu˙x,1⋯Xu˙x,nXu˙z,1⋯Xu˙z,n0n×1Zu˙x,1⋯Zu˙x,nZu˙z,1⋯Zu˙z,nZφ˙y,1⋯Zφ˙y,n0n×1Mu˙z,1⋯Mu˙z,nMφ˙y,1⋯Mφ˙y,n,Q2=01×n01×nXφy,1⋯Xφy,n01×n01×nZφy,1⋯Zφy,n01×n01×nMφy,1⋯Mφy,n,Q3=Xδe,1⋯Xδe,nXδtZδe,1⋯Zδe,nZδtMδe,1⋯Mδe,nMδt,where Zu˙z is related to some derivatives (e.g., CLu˙z). The exact expressions are the same as the derivatives of general flight dynamics equations. More details about them can be referred to [21].


### 3.4. Total Linearized Equation

Combining (26) and (39), the final total linear equations can be derived as follows:(42)x˙r=Arxr+Arnexne+Brur,x˙ne=Anerxr+Anexne+Bneur.The independent variable is(43)x=xrT,xneTTand the state matrix is (44)A=ArArneAnerAne.And the control matrix is:(45)B=BrBne.The eigenvalues of the above state matrix are the modal eigenvalues of the elastic UAV considering the coupling of aeroelastic and flight dynamics:(46)x˙=Ax+Bu.


## 3.1. Linearized Aerodynamics

The movement of wing section is separated into two parts: the rigid movement of UAV and the elastic movement of structure. Then, the lift coefficients applied on the wing section can be deployed in theaB-frame as follows:(23)CLα,CLq=0.5CLα+2CLα34-x-gs,CLp=-2CLαygsl,CLφy=CLα,CLu˙z=CLα1V,CLφ˙y=CLα1-a-c2V.It should be noted that the variables in the above equations are measured in theaB-frame and B-frame.According to the 2nd equation of (12), the aerodynamics moment to the elastic center (EC) is(24)Cmα,E=-CLαx-AEB≈CLα14+12a-,Cmq,E=0.5CLαx-CEB+2CLαx-Cgsx-AEB,Cmp,E=-2CLαygslx-AEB,Cmφy,E=Cmα,E,Cmu˙z,E=Cmα,E1V,Cmφ˙y,E=CLα14a--12a-2c2V.From the 2nd equation of (18), the difference of lift expressed in the aM-frame and aB-frame is only cos⁡⁡φBM,x. Hence, by replacing CLα in (23) with CLαcos⁡⁡φBM,x, the lift expressed in the aB-frame will be derived.Substituting the expression of lift ataM-frame into (19), the pitch moment about the center of gravity of the UAV can be obtained as follows:(25)Cmα,g=-CLαcos⁡⁡φMB,xx-As-x-gs,Cmq,g=CLαcos⁡⁡φMB,x0.5x-CgB-2x-CgB2,Cmp,g=-2CLαygslx-CgB,Cmφy,g=Cmα,g,Cmu˙z,g=Cmα,g1V,Cmφ˙y,g=-2CLαcos⁡⁡φMB,x12-x-gs2bV.


## 3.2. Linearized Flight Dynamics Equation

It can be found that the form of (4) and (5) is similar to that of flight dynamics equations of normal rigid vehicles. Hence, they can also be linearized by small disturbance method and written into two sets of equations: longitudinal and lateral-directional. Different from the normal rigid vehicle, the external force such as aerodynamics, thrust, and gravity are related to rigid motion and elastic motion for elastic model. The linear equation of flight dynamics for elastic UAV is (26)x˙r=Arxr+Brur+Arexe,where xr is the independent variable of rigid motion, with xr=[v,α,q,θ]T in longitudinal equations; xe=[u˙xT,u˙zT,φ˙yT,uxT,uzT,φyT]T is the independent variable of elastic motion; ur is the control variable of rigid motion, ur=[δe,1⋯δe,40,δt]T; that is, all the elevons can be controlled independently, but all throttles act conformably again for simplification; Ar is the same as conventional rigid vehicles as shown in [20]; and Br is similar to the rigid model also, but the column is equal to the length of ur; Are is the influence matrix of elastic motion to rigid motion, for longitudinal motion:(27)Are=01×n01×n01×n01×n01×n01×n01×nZu˙z,1⋯Zu˙z,nZφ˙y,1⋯Zφ˙y,n01×n01×nZφy,1⋯Zφy,n01×nMu˙z,1⋯Mu˙z,nMφ˙y,1⋯Mφ˙y,n01×n01×nMφy,1⋯Mφy,n01×n01×n01×n01×n01×n01×n.The details ofAre can be found in [21].


## 3.3. Linearized Structural Dynamics Equation

Ignoring nonlinear coupling terms in (1), and assuming that the translational strain is small and the external moment and the moment of inertia of x and z direction are equal to zero, (1) can be reduced into(28)μV˙x=EIzκz′′+fx,μV˙z=-EIxκx′′+fz,IyyΩ˙y=EIyκy′+mz,κ˙x=Vz′′,κ˙y=Ωy′,κ˙z=-Vx′,μV˙y=fy.The equations above can be rearranged into a compact form as follows:(29)μx¨+EIzx′′′′=fx,μz¨+EIxz′′′′=fz,Iyyφ¨Bg,y-GJzφBg,y′′=my,μV˙y=fy.Since the focus of this paper is the flight dynamics characteristics of the total UAV, the structural dynamics characteristics are used to derive the increment of aerodynamics force caused by aeroelasticity only and many assumptions can be introduced to linearize the equations above for simplification.Equation (29) can be simplified to the form as in (30) by rewriting in M-frame and assuming that the structural elastic vibrations are small perturbations near the static trim state and the elastic motion in x direction is ignorable. In the case of geometry, it is assumed that the M-frame is parallel to R-frame and the sweep angle of wing is zero. Besides, in mean axis system, the elastic motion is influenced by rigid motion through aerodynamics force only:(30)μu¨xM,x+EIzuxM,x′′′′=cos⁡⁡φBM,zfB,x,μu¨zM,z+EIxuzM,z′′′′=cos⁡⁡φBM,xfB,z,Iyyφ¨BM,y-GJzφBM,y′′=mB,y,where φBM,y represents the angle between x-axis of B-frame and the oxy plane of M-frame. All variables in (30) are in incremental form and the notation Δ has been neglected for simplification.Replacing the spatial derivatives by (3), and knowing the fact that the wing tip is free, the boundary conditions can be defined that the bending moments and shear forces at the first and last nodes are zero for the 2nd equation in (30); that is [15],(31)uz,1′′=0,uz,n′′=0,uz,1′′′=0,(32)uz,n′′′=0.Thus, the structural dynamics equation can be derived as follows: (33)Muzu¨z+Kuzuz=Quz,where(34)uz=u1⋯unT,Muz=μ1⋱μn,Kuz=2-42-25-411-46-41⋱⋱⋱⋱⋱1-46-411-45-22-42∗1Δl4diag⁡EIz,i.For the 2nd equation in (30), it is the same as (33); that is,(35)Muxu¨x+Kuxux=Qux.For the 3rd equation in (30), the boundary conditions are [15](36)φBM,y,1′=0,φBM,y,n′=0.Then, the structural dynamics equation becomes(37)Mφyφ¨y+Kφyφy=Qφy,where(38)φy=φBM,y,1⋯φBM,y,nT,Mφy=I1⋱In,Kφy=--221-21⋱⋱⋱1-212-2∗1Δl2diag⁡GJi.In order to further simplify the linear structural dynamics equations, (33), (35), and (37) can be expanded in the modal coordinate system, and the final form can be derived as follows:(39)x˙ne=Anexne+Anerxr+Bneur,where the independent variables are(40)xne=x˙nxn,xe=μxn,xe=uxT,uzT,φyTT,where xn is the modal displacements of structural elastic motion, which is the symmetric modes for longitudinal motion; μ is the transfer matrix from modal coordinate to physical coordinate.For longitudinal motion, the coefficient matrices in (39) are(41)Ane=Mn-1Qn1Mn-1Qn2-KnI0,Aner=Mn-1Qn00,Bne=Mn-1Qn30,Mn=μTMμ,Kn=μTKμ,Qni=μTQi,Qn0=μTQ0,Qn1=μTQ1μ,Qn2=μTQ2μ,Qn3=μTQ3μ,Q0=Xv,1⋮Xv,nXα,1⋮Xα,n0n×10n×1Zv,1⋮Zv,nZα,1⋮Zα,nZq,1⋮Zq,n0n×1Mv,1⋮Mv,nMα,1⋮Mα,nMq,1⋮Mq,n0n×1,Q1=Xu˙x,1⋯Xu˙x,nXu˙z,1⋯Xu˙z,n0n×1Zu˙x,1⋯Zu˙x,nZu˙z,1⋯Zu˙z,nZφ˙y,1⋯Zφ˙y,n0n×1Mu˙z,1⋯Mu˙z,nMφ˙y,1⋯Mφ˙y,n,Q2=01×n01×nXφy,1⋯Xφy,n01×n01×nZφy,1⋯Zφy,n01×n01×nMφy,1⋯Mφy,n,Q3=Xδe,1⋯Xδe,nXδtZδe,1⋯Zδe,nZδtMδe,1⋯Mδe,nMδt,where Zu˙z is related to some derivatives (e.g., CLu˙z). The exact expressions are the same as the derivatives of general flight dynamics equations. More details about them can be referred to [21].


## 3.4. Total Linearized Equation

Combining (26) and (39), the final total linear equations can be derived as follows:(42)x˙r=Arxr+Arnexne+Brur,x˙ne=Anerxr+Anexne+Bneur.The independent variable is(43)x=xrT,xneTTand the state matrix is (44)A=ArArneAnerAne.And the control matrix is:(45)B=BrBne.The eigenvalues of the above state matrix are the modal eigenvalues of the elastic UAV considering the coupling of aeroelastic and flight dynamics:(46)x˙=Ax+Bu.


## 4. Model Verification

### 4.1. Geometrically Nonlinear Deformation

In order to validate the effectiveness of the geometrically nonlinear structural model and central difference method proposed in this paper, a cantilever beam with bending rigidityEI=1 and length l = 1 m is employed. While applying increasing moment M at the tip of the beam, the deformation of the beam is shown in Figure 3. It can be seen that while M=2π, the angular displacement is 2π, the tip coincides with the root and the beam is deformed into a circle, which is consistent with the analytical results [16]. The effectiveness of structural model is validated.Figure 3
Beam deformation with increasing tip moment.


### 4.2. Aerodynamics Model

Aerodynamics model is the most complex and important external force model in the problem of flight dynamics and control, so it must be checked carefully before further studies. The lift coefficient distribution along the wing span obtained from two methods are compared in Figure4; one is derived from the Thin-Strip lifting-line theory proposed in Section 3.1 in this paper, and the other is the Vortex-Lattice method from AVL software. It can be seen that the precision of Thin-Strip lifting-line theory model can satisfy the application requirements of this paper.Figure 4
Comparison of spanwise lift coefficient distribution.


### 4.3. Linear Structural Model

In order to investigate the error of different modal truncation, 4 linear structural models of the wing are compared with the geometrically nonlinear one: truncating at the 1st, 2nd, 4th, and 6th bending mode, respectively, and then applying a flapwise impulse force to the tip of the wing. The time domain responses of the wing tip of different models are shown in Figure5.Figure 5
Time domain responses of the wing tip of different truncated linear wing models.It can be seen that, if the deformation is not large enough, for the flapwise deformation, if we employ the first 6 bending modes to constitute the linear structural model, the responses of the wing tip are coincident with the nonlinear one well. Similarly, it can be found from simulation that the first 4 twist modes and 2 chordwise modes are precise enough. So, the first 6 flapwise bending, 4 twist and 2 chordwise bending modes are adopted to constitute the linear structural model of the wing in the following of this paper.


## 4.1. Geometrically Nonlinear Deformation

In order to validate the effectiveness of the geometrically nonlinear structural model and central difference method proposed in this paper, a cantilever beam with bending rigidityEI=1 and length l = 1 m is employed. While applying increasing moment M at the tip of the beam, the deformation of the beam is shown in Figure 3. It can be seen that while M=2π, the angular displacement is 2π, the tip coincides with the root and the beam is deformed into a circle, which is consistent with the analytical results [16]. The effectiveness of structural model is validated.Figure 3
Beam deformation with increasing tip moment.


## 4.2. Aerodynamics Model

Aerodynamics model is the most complex and important external force model in the problem of flight dynamics and control, so it must be checked carefully before further studies. The lift coefficient distribution along the wing span obtained from two methods are compared in Figure4; one is derived from the Thin-Strip lifting-line theory proposed in Section 3.1 in this paper, and the other is the Vortex-Lattice method from AVL software. It can be seen that the precision of Thin-Strip lifting-line theory model can satisfy the application requirements of this paper.Figure 4
Comparison of spanwise lift coefficient distribution.


## 4.3. Linear Structural Model

In order to investigate the error of different modal truncation, 4 linear structural models of the wing are compared with the geometrically nonlinear one: truncating at the 1st, 2nd, 4th, and 6th bending mode, respectively, and then applying a flapwise impulse force to the tip of the wing. The time domain responses of the wing tip of different models are shown in Figure5.Figure 5
Time domain responses of the wing tip of different truncated linear wing models.It can be seen that, if the deformation is not large enough, for the flapwise deformation, if we employ the first 6 bending modes to constitute the linear structural model, the responses of the wing tip are coincident with the nonlinear one well. Similarly, it can be found from simulation that the first 4 twist modes and 2 chordwise modes are precise enough. So, the first 6 flapwise bending, 4 twist and 2 chordwise bending modes are adopted to constitute the linear structural model of the wing in the following of this paper.


## 5. Flight Dynamics of Highly Flexible UAV

### 5.1. Trimming

Two kinetics models are employed for trimming: the “Undeformed” model represents the common rigid flight dynamics model; the “Deformed” model corresponding to (1) represents the flexible flight dynamics model.In level flight, assuming that all the 40 distributed elevons deflect conformably for trimming for simplification, the highly flexible UAV will deform under the aeroelastic effect. In fact, the deformed UAV in trimming is shown in Figure6, are the flapwise position and pitch angle of the UAV are shown in Figures 7 and 8.Figure 6
Deformed UAV in trimming.Figure 7
Flapwise position of UAV of different calculation model.Figure 8
Pitch angle of UAV of different calculation model.It can be seen that the maximum wing deformation is 4.58 m at the wing tip, and the pitch angle at wing tip is larger than the root for about 3°. Larger pitch angle brings larger lift, and larger lift at wing tip will bring larger deformation. So, all the distributed elevons deflecting conformably is not the optimal control case while considering aeroelasticity.


### 5.2. Modal Analysis

The “Underformed” and “Deformed” models are corresponding to conventional rigid vehicles, and the linear “Flexible” model is based on (42). The longitudinal modes of different calculation models are shown in Table 2.Table 2
Roots of longitudinal mode of different calculation models.
Short period
Phugoid
1st bending
Undeformed
-
6.710
±
3.925
i
-
0.046
±
1.003
i
/
Deformed
-
5.985
±
3.100
i
-
0.00
2
±
1.063
i
/
Flexible
-
0.394
±
1.707
i
0.00
2
±
1.0
31
i
-
3.819
±
1.459
iThe results showed that the difference of these models is huge in longitudinal mode. Since the deformation of wing has increased the UAV’s moment of inertia about the pitch axis, the frequency of short period mode has obviously decreased for the deformed model. The serious coupling between short period mode and 1st bending mode also significantly influences the roots of short period mode, which then results in the phugoid changing from stable to unstable. These tendencies are similar to the results presented by Hodges et al. Through these analyses, it is clear that the flexible model is the only one which is able to accurately estimate the longitudinal flight dynamics behaviors.


## 5.1. Trimming

Two kinetics models are employed for trimming: the “Undeformed” model represents the common rigid flight dynamics model; the “Deformed” model corresponding to (1) represents the flexible flight dynamics model.In level flight, assuming that all the 40 distributed elevons deflect conformably for trimming for simplification, the highly flexible UAV will deform under the aeroelastic effect. In fact, the deformed UAV in trimming is shown in Figure6, are the flapwise position and pitch angle of the UAV are shown in Figures 7 and 8.Figure 6
Deformed UAV in trimming.Figure 7
Flapwise position of UAV of different calculation model.Figure 8
Pitch angle of UAV of different calculation model.It can be seen that the maximum wing deformation is 4.58 m at the wing tip, and the pitch angle at wing tip is larger than the root for about 3°. Larger pitch angle brings larger lift, and larger lift at wing tip will bring larger deformation. So, all the distributed elevons deflecting conformably is not the optimal control case while considering aeroelasticity.


## 5.2. Modal Analysis

The “Underformed” and “Deformed” models are corresponding to conventional rigid vehicles, and the linear “Flexible” model is based on (42). The longitudinal modes of different calculation models are shown in Table 2.Table 2
Roots of longitudinal mode of different calculation models.
Short period
Phugoid
1st bending
Undeformed
-
6.710
±
3.925
i
-
0.046
±
1.003
i
/
Deformed
-
5.985
±
3.100
i
-
0.00
2
±
1.063
i
/
Flexible
-
0.394
±
1.707
i
0.00
2
±
1.0
31
i
-
3.819
±
1.459
iThe results showed that the difference of these models is huge in longitudinal mode. Since the deformation of wing has increased the UAV’s moment of inertia about the pitch axis, the frequency of short period mode has obviously decreased for the deformed model. The serious coupling between short period mode and 1st bending mode also significantly influences the roots of short period mode, which then results in the phugoid changing from stable to unstable. These tendencies are similar to the results presented by Hodges et al. Through these analyses, it is clear that the flexible model is the only one which is able to accurately estimate the longitudinal flight dynamics behaviors.


## 6. Distributed LQG/LTR Controller Design

### 6.1. Statement of LQG/LTR Control

For a standard state space system described in (46), the design procedure of LQG/LTR controller in the output port is as follows:(1)
Design a LQR regulatorF, which ensures that the singular value of open loop system FG(s) on a certain frequency satisfies the robustness requirement.(2)
Design a target loop: select appropriateΞ and Θ, so as to the singular value of full-state feedback loop transfer function Gc(s)G(s) on the considering frequency close to that of FG(s) adequately.Where, the controlled plantG(s) is(47)Gs=CsI-A-1Band the derived LQG/LTR controller plant Gc(s) is(48)Gcs=FsI-A+BF+LC-1L.The methodology of LQR is searching a full state feedback controlleru=-Fx to minimize the quadratic cost function:(49)J=min⁡⁡∫0∞xTQx+uTRudtsubject to the system dynamics described in (46).The expression ofF can be calculated by(50)F=R-1BTP,where P is the solution of the following algebraic Riccati equation:(51)ATP+PA-PBR-1BTP+Q=0.The solution of L in (48) is similar to F; one just needs to replace the matrix A,B,Q,R in (49)–(51) to AT,CT,Ξ+ρBBT,Θ, respectively, where Ξ+ρBBT represents the noise intensity; different noise intensity can be adjusted by ρ; Θ represents the sensor noise intensity. A well-designed LQG/LTR controller should be that, when ρ→∞, the singular value of Gc(s)G(s) on the considering frequency close to that of FG(s) adequately. But it should be noted that larger ρ will derive larger gain of Gc(s).


### 6.2. Designing LQG/LTR Pitch Angle Controller

Before designing LQG/LTR control law, one shall augment the plant to satisfy the specific performance requirement. For example, adding integrators will assure the zero steady-state errors to step inputs. So, in order to control the pitch angle precisely, the integrator of pitch angleθi and velocity Vi should be included in the control plant. Then, the control plant in (46) becomes(52)x˙θ˙iV˙i=A04×104×1000101×2n00100001×2n00xθiVi+B00u.The critical factor for a practical controller is that not all the state variables are easy to be measured in engineering. So only a part of variables can be used to implement the output feedback control. Because the LTR technology can recover the stable margin of LQR, which is known as an excellent state feedback robust controller, the robustness of LQG/LTR output feedback controller can be guaranteed.Because the flapwise bending motion of the wing is easy to be measured, souz together with V, α, q, θ, θi, Vi are selected to be the measured variables for the controller. Then, we can apply the LQG/LTR method described above in (52).If we selectρ=1×106, the singular values of pitch angle loop of FG(s) denoted by LQR and Gc(s)G(s) denoted by LTR are compared in Figure 9.Figure 9
Singular values of pitch angle loop.It can be seen from Figure9 that, after the loop transfer recovery, the two singular values coincide perfectly on low frequency band, which indicates that the response of the close loop system to common control such as step input is quick and the tracking error is small even with model deviation. On high frequency band, the singular value of LTR decreases more quickly than that of LQR, which will provide excellent robustness to the high frequency disturbance such as unmodeled high frequency structural dynamics and unsteady aerodynamics. So, the LTR satisfies the required performance well.


### 6.3. Flight Simulation

There are 40 distributed elevons in the large scale highly flexible solar-powered UAV. The details about control actuator are omitted temporally for brief while designing the LQG/LTR controller. But in the stage of flight simulation, two kinds of actuators are employed to verify the robustness of the controller. The model of the elevons actuator is(53)δea=10.05s+1u.The model of the throttles is(54)δt=10.5s+1u.Their bandwidths are 20 and 2 rad/s, respectively.Applying the LQG/LTR controller proposed above to the linear structural and flight dynamics model described in (52), when given a desired pitch angle θc=5°, the longitudinal response of solar aircraft is shown in Figures 10–14. As shown in Figure 10, the adjust time of pitch angle is about 3 seconds; because of the use of integrator, steady-state error is eliminated. Seen from Figure 11, the variation of angle of attack in the whole process does not exceed 1.4°, so the additional load it brings is very small, indicating that the overall movement of the highly-flexible UAV is relatively stable. Seen from Figure 12, the largest deflection angle of elevons is less than 20 degrees, and it occurs in the initial phase only, and when approaching a steady state, the deflections are no more than 5 degrees, far less than the limitation. It must be noted that the deflection angle of the 40 elevons are all different, because in addition to providing pitch angle control, distributed elevons also need to suppress the elastic deformation movement of the highly-flexible wing in the LQG/LTR controller. Seen from Figure 13, the maximum deformation of the wing occurs in the tip and the maximum magnitude of wing deformation is only about 0.09 m, far less than the deformation of 4.58 m in trimming, reminding that all the elevons deflect conformably in trimming. Seen from Figure 14, the max rate of twist angle is only 0.26 rad/s, corresponding to the reduced frequency of 0.0144. In a word, with distributed elevons, the use of LQG/LTR control method allows faster attitude control of the large scale highly flexible UAV while good aeroelastic deformation suppression is achieved simultaneously.Figure 10
Response of pitch angle.Figure 11
Response of angle of attack.Figure 12
Response of elevons at different stations.Figure 13
Response of vertical displacements of wing at different stations.Figure 14
Response of rate of twist angle of wing at different stations.In addition, reminding that the linear structural model coincide with the nonlinear one perfectly on the condition of small deformation (shown in Figure5). Since the LQG/LTR controller has restrained the structural deformation and incident angle in a relative small range, it can be deduced that the simulation results in this section are reliable.


## 6.1. Statement of LQG/LTR Control

For a standard state space system described in (46), the design procedure of LQG/LTR controller in the output port is as follows:(1)
Design a LQR regulatorF, which ensures that the singular value of open loop system FG(s) on a certain frequency satisfies the robustness requirement.(2)
Design a target loop: select appropriateΞ and Θ, so as to the singular value of full-state feedback loop transfer function Gc(s)G(s) on the considering frequency close to that of FG(s) adequately.Where, the controlled plantG(s) is(47)Gs=CsI-A-1Band the derived LQG/LTR controller plant Gc(s) is(48)Gcs=FsI-A+BF+LC-1L.The methodology of LQR is searching a full state feedback controlleru=-Fx to minimize the quadratic cost function:(49)J=min⁡⁡∫0∞xTQx+uTRudtsubject to the system dynamics described in (46).The expression ofF can be calculated by(50)F=R-1BTP,where P is the solution of the following algebraic Riccati equation:(51)ATP+PA-PBR-1BTP+Q=0.The solution of L in (48) is similar to F; one just needs to replace the matrix A,B,Q,R in (49)–(51) to AT,CT,Ξ+ρBBT,Θ, respectively, where Ξ+ρBBT represents the noise intensity; different noise intensity can be adjusted by ρ; Θ represents the sensor noise intensity. A well-designed LQG/LTR controller should be that, when ρ→∞, the singular value of Gc(s)G(s) on the considering frequency close to that of FG(s) adequately. But it should be noted that larger ρ will derive larger gain of Gc(s).


## 6.2. Designing LQG/LTR Pitch Angle Controller

Before designing LQG/LTR control law, one shall augment the plant to satisfy the specific performance requirement. For example, adding integrators will assure the zero steady-state errors to step inputs. So, in order to control the pitch angle precisely, the integrator of pitch angleθi and velocity Vi should be included in the control plant. Then, the control plant in (46) becomes(52)x˙θ˙iV˙i=A04×104×1000101×2n00100001×2n00xθiVi+B00u.The critical factor for a practical controller is that not all the state variables are easy to be measured in engineering. So only a part of variables can be used to implement the output feedback control. Because the LTR technology can recover the stable margin of LQR, which is known as an excellent state feedback robust controller, the robustness of LQG/LTR output feedback controller can be guaranteed.Because the flapwise bending motion of the wing is easy to be measured, souz together with V, α, q, θ, θi, Vi are selected to be the measured variables for the controller. Then, we can apply the LQG/LTR method described above in (52).If we selectρ=1×106, the singular values of pitch angle loop of FG(s) denoted by LQR and Gc(s)G(s) denoted by LTR are compared in Figure 9.Figure 9
Singular values of pitch angle loop.It can be seen from Figure9 that, after the loop transfer recovery, the two singular values coincide perfectly on low frequency band, which indicates that the response of the close loop system to common control such as step input is quick and the tracking error is small even with model deviation. On high frequency band, the singular value of LTR decreases more quickly than that of LQR, which will provide excellent robustness to the high frequency disturbance such as unmodeled high frequency structural dynamics and unsteady aerodynamics. So, the LTR satisfies the required performance well.


## 6.3. Flight Simulation

There are 40 distributed elevons in the large scale highly flexible solar-powered UAV. The details about control actuator are omitted temporally for brief while designing the LQG/LTR controller. But in the stage of flight simulation, two kinds of actuators are employed to verify the robustness of the controller. The model of the elevons actuator is(53)δea=10.05s+1u.The model of the throttles is(54)δt=10.5s+1u.Their bandwidths are 20 and 2 rad/s, respectively.Applying the LQG/LTR controller proposed above to the linear structural and flight dynamics model described in (52), when given a desired pitch angle θc=5°, the longitudinal response of solar aircraft is shown in Figures 10–14. As shown in Figure 10, the adjust time of pitch angle is about 3 seconds; because of the use of integrator, steady-state error is eliminated. Seen from Figure 11, the variation of angle of attack in the whole process does not exceed 1.4°, so the additional load it brings is very small, indicating that the overall movement of the highly-flexible UAV is relatively stable. Seen from Figure 12, the largest deflection angle of elevons is less than 20 degrees, and it occurs in the initial phase only, and when approaching a steady state, the deflections are no more than 5 degrees, far less than the limitation. It must be noted that the deflection angle of the 40 elevons are all different, because in addition to providing pitch angle control, distributed elevons also need to suppress the elastic deformation movement of the highly-flexible wing in the LQG/LTR controller. Seen from Figure 13, the maximum deformation of the wing occurs in the tip and the maximum magnitude of wing deformation is only about 0.09 m, far less than the deformation of 4.58 m in trimming, reminding that all the elevons deflect conformably in trimming. Seen from Figure 14, the max rate of twist angle is only 0.26 rad/s, corresponding to the reduced frequency of 0.0144. In a word, with distributed elevons, the use of LQG/LTR control method allows faster attitude control of the large scale highly flexible UAV while good aeroelastic deformation suppression is achieved simultaneously.Figure 10
Response of pitch angle.Figure 11
Response of angle of attack.Figure 12
Response of elevons at different stations.Figure 13
Response of vertical displacements of wing at different stations.Figure 14
Response of rate of twist angle of wing at different stations.In addition, reminding that the linear structural model coincide with the nonlinear one perfectly on the condition of small deformation (shown in Figure5). Since the LQG/LTR controller has restrained the structural deformation and incident angle in a relative small range, it can be deduced that the simulation results in this section are reliable.


## 7. Conclusions

In this paper, the modeling, stability, and control characteristics of a highly flexible large scale solar-powered UAV with distributed all-span multi elevons were presented.A geometrically nonlinear intrinsic beam model is introduced to establish the equation of motion of flexible solar-powered UAV. Based on it, the decoupled linear flight dynamics and structural dynamics equations of motions are derived in mean axis frame. Steady aerodynamics model is derived by theoretical airfoil data and corrected for 3D effect by lifting-line theory, quasisteady aerodynamics model is derived by Theodorsen unsteady aerodynamics model, and the aerodynamics derivatives of every wing sections are derived and transferred to different frames for the multirole in linear equations of flight dynamics or structural dynamics.The nonlinear coupled structural and flight dynamics model was employed to trim the large scale flexible UAV, under the assumption that all the 40 distributed elevons deflect conformably. The maximum aeroelastic defection of wingtip is 4.58 m in flapwise, about 13% to the halfspan. The pitch angle at wing tip is larger than the root for about 3°. Larger pitch angle brings larger lift and leads to larger deformation. So all the distributed elevons deflecting conformably are not the optimal control case while considering aeroelastic effects.Three different linear dynamics models, undeformed, deformed, and elastic, were compared for longitudinal mode analysis. The results showed that the difference of three models is huge. Since the deformation of wing has increased the UAV’s moment of inertia about the pitch axis, the frequency of short period mode has obviously decreased for the deformed model. The serious coupling between short period mode and 1st bending mode also significantly influences the roots of short period mode, which then results in the phugoid changing from stable to unstable. So, it is clear that the flexible model is the only one which is able to accurately estimate the longitudinal flight dynamics behaviors and control the large scale highly flexible UAV exactly.Forty distributed all-span elevons were employed to control the attitude and suppress the aeroelastic deformation. LQG/LTR methodology was employed to accomplish the dynamics performance and robustness requirement. The flight simulation showed that the designed LQG/LTR controller can provide the adjust time of pitch angle in 3 seconds without steady-state error, and the maximum magnitude of wing deformation is only about 0.09 m, far less than the deformation of 4.58 m in trimming. So, LQG/LTR control method with distributed elevons is suitable for large scale highly flexible solar-powered UAV.



---

*Source: 195390-2015-09-16.xml*