Nonlinear F-16 Model Description
L. Sonneveldt Control & Simulation Division Faculty of Aerospace Engineering Delft University of Technology The Netherlands Version 1.0, March 2010
c Control & Simulation Division, Delft University of Technology Copyright All rights reserved L. Sonneveldt 1
Nomenclature State variables: V T α, β φ,θ,ψ q0 , q1 , q2 , q3 p, q, r xE , yE , zE pow
= = = = = = = =
total velocity, m/s angle of attack, rad angle of sideslip, rad Euler angles, rad quaternion components body-axis roll, pitch and yaw rates, rad/s aircraft position w.r.t. reference point, m power setting, %
= = = = =
throttle setting, 0 − 1 elevator deflection, rad aileron deflection, rad rudder deflection, rad leading edge flap deflection, rad
Control variables: δth δe δa δr δLEF
Additional parameters b c¯ C F T g1 , g2 , g3 H eng I x , I y , I z I xz ¯ M, ¯ N ¯ L, m M ps q¯ S T T s/b T w/b u,v,w ¯ Y¯ , Z ¯ X, xcg xcgr γ ρ ∗
= = = = = = = = = = = = = = = = = = = = = = =
reference wing span, m mean aerodynamic chord, m Aerodynamic coefficient of * thrust force, N gravity components, m/s2 engine angular momentum, kg.m2 /s moments of inertia kg.m2 product moment of inertia kg.m2 rolling, pitching and yawing moments, N m mass, kg Mach number static pressure, P a dynamic pressure, P a reference wing surface area, m2 Temperature, K Rotation matrix body-fixed to stability-axes reference frame Rotation matrix body-fixed to wind-axes reference frame body-axis velocity components, m/s axial, lateral and normal force components, N center of gravity location, m reference center of gravity location, m flight path angle, rad air density, kg/m3
2
1
Aircraft Dynamics
In this section the equations of motion for the F-16 model are derived, this derivation is based on [1, 2, 3]. A very thorough discussion on flight dynamics can be found in the course notes [4].
1.1
Reference Frames
The reference frames used in this thesis are the earth-fixed reference frame F E , used as the inertial frame and the vehicle carried local earth reference frame F O with its origin fixed in the center of gravity of the aircraft which is assumed to have the same orientation as F E ; the wind-axes reference frame F W , obtained from F O by three successive rotations of χ, γ and µ; the stability-axes reference frame F S , obtained from F W by a rotation of −β; and finally the body-fixed reference frame F B , obtained from F S by a rotation of α as is also indicated in Figure 1. The body-fixed reference frame F B can also be obtained directly from
Figure 1: Aircraft reference frames: earth-fixed reference frame F E , body-fixed reference frame F B , stability-axes reference frame F S and wind-axes reference frame F W F O by three successive rotations of yaw angle ψ, pitch angle θ and roll angle φ. All reference frames are right-handed and orthogonal. In the earth-fixed reference frame the zE -axis points to the center of the earth, the xE -axis points in some arbitrary direction, e.g. the north, and the yE -axis is perpendicular to the xE -axis.
3
The transformation matrices from F B to F S and from F B to F W are defined as T s/b
1.2
=
cos α 0 − sin α
0 1 0
sin α 0 cos α
,
T w/b
=
cos α cos β − cos α sin β − sin α
sin β cos β 0
sin α cos β − sin α sin β cos α
.
Aircraft Variables
A number of assumptions has to be made, before proceeding with the derivation of the equations of motion: 1. The aircraft is a rigid-body , which means that any two points on or within the airframe remain fixed with respect to each other. This assumption is quite valid for a small fighter aircraft. 2. The earth is flat and non-rotating and regarded as an inertial reference . This assumption is valid when dealing with control design of aircraft, but not when analyzing inertial guidance systems. 3. The mass is constant during the time interval over which the motion is considered , the fuel consumption is neglected during this time-interval. This assumption is necessary to apply Newton’s motion laws. 4. The mass distribution of the aircraft is symmetric relative to the X B OZ B plane, this implies that the products of inertia I yz and I xy are equal to zero. This assumption is valid for most aircraft. Under the above assumptions the motion of the aircraft has six degrees of freedom (rotation and translation in three dimensions). The aircraft dynamics can be described by its position, orientation, velocity and angular velocity over time. pE = (xE , yE , zE )T is the position vector expressed in an earth-fixed coordinate system. V is the velocity vector given by V = (u,v,w)T , where u is the longitudinal velocity, v the lateral velocity and w the normal velocity. The orientation vector is given by Φ = (φ,θ,ψ)T , where φ is the roll angle, θ the pitch angle and ψ the yaw angle, and the angular velocity vector is given by ω = ( p, q, r)T , where p,q and r are the roll, pitch and yaw angular velocities, respectively. Various components of the aircraft motions are illustrated in Figure 2. The relation between the attitude vector Φ and the angular velocity vector ω is given as 1 sin φ tan θ cos φ tan θ ˙ 0 cos φ − sin φ Φ= ω (1) sin φ cos φ 0 cos θ cos θ
Defining V T as the total velocity and using Figure 2, the following relations can be derived: V T α β
=
u + v + w 2
2
w u v = arcsin V T = arctan
4
2
(2)
Figure 2: Aircraft orientation angles φ, θ and ψ, aerodynamic angles α and β, and the angular rates p, q and r. All angles and rates are defined positive in the figure. Furthermore, when β = φ = 0, the flight path angle γ can be defined as γ = θ − α
1.3
(3)
Equations of Motion for a Rigid Body Aircraft
The equations of motion for the aircraft can be derived from Newton’s Second Law of motion, which states that the summation of all external forces acting on a body must be equal to the time rate of change of its momentum, and the summation of the external moments acting on a body must be equal to the time rate of change of its angular momentum. In the inertial, Earth-fixed reference frame F E , Newton’s Second Law can be expressed by two vector equations [4] F
=
M
=
d (mV) dt dH dt E
E
(4) (5)
where F represents the sum of all externally applied forces, m is the mass of the aircraft, M represents the sum of all applied torques and H is the angular momentum. 1.3.1
Force Equation
First, to further evaluate the force equation (4) it is necessary to obtain an expression for the time rate of change of the velocity vector with respect to earth. This process is complicated by the fact that the velocity vector may be rotating while it is changing in magnitude. Using the equation of Coriolis in appendix A of [1] results in F=
d (mV) dt
B
+ ω × mV ,
(6)
whereω is the total angular velocity of the aircraft with respect to the earth (inertial reference frame). Expressing the vectors as the sum of their components 5
with respect to the body-fixed reference frame F B gives V ω
= iu + jv + kw = i p + jq + kr
(7) (8)
where i, j and k are unit vectors along the aircraft’s xB , yB and zB axes, respectively. Expanding (6) using (7), (8) results in F x F y F z
= m(u˙ + qw − rv) = m(v˙ + ru − pw) = m(w˙ + pv − qu)
(9)
where the external forces F x , F y and F z depend on the weight vector W, the aerodynamic force vector R and the thrust vector E. It is assumed the thrust produced by the engine, F T , acts parallel to the aircraft’s xB -axis. Hence, E x E y
= F T = 0
E z
= 0
(10)
The components of W and R along the body-axes are W x W y W z
= −mg sin θ = mg sin φ cos θ = mg cos φ cos θ
(11)
¯ = X ¯ = Y ¯ = Z
(12)
and Rx Ry Rz
¯ , Y ¯ where g is the gravity constant. The size of the aerodynamic forces X ¯ is determined by the amount of air diverted by the aircraft in different and Z directions. The amount of air diverted by the aircraft mainly depends on the following factors: • the total velocity V T (or Mach number M ) and density of the airflow ρ, • the geometry of the aircraft: wing area S , wing span b and mean aerodynamic chord ¯c, • the orientation of the aircraft relative to the airflow: angle of attack α and side slip angle β, • the control surface deflections δ, • the angular rates p,q,r, 6
There are other variables such as the time derivatives of the aerodynamic angles that also play a role, but these effects are less prominent, since it is assumed that the aircraft is a rigid body. This motivates the standard way of modeling the aerodynamic force: ¯ = q¯SC XT (α,β,p,q,r,δ,...) X ¯ = q¯SC YT (α,β,p,q,r,δ,...) Y ¯ = q¯SC ZT (α,β,p,q,r,δ,...) Z
(13)
where q¯ = 12 ρV T 2 is the aerodynamic pressure. The air density ρ is calculated according to the International Standard Atmosphere (ISA) as given in Section B. The coefficients C XT , C Y T and C ZT are usually obtained from (virtual) wind tunnel data and flight tests. Combining equations (11) and (12) and the thrust components (10) with (9), results in the complete body-axes force equation: ¯ F T − mg sin θ X + ¯ + mg sin φ cos θ Y ¯ mg cos φ sin θ Z + 1.3.2
= m(u˙ + qw − rv) = m(v˙ + ru − pw)
(14)
= m(w˙ + pv − qu)
Moment Equation
To obtain the equations for angular motion, consider again Equation (5). The time rate of change of H is required and since H can change in magnitude and direction, (5) can be written as dH M= dt
B
+ω×H
(15)
In the body-fixed reference frame, under the rigid body and constant mass assumptions, the angular momentum H can be expressed as H = I ω
(16)
where, under the symmetrical aircraft assumption, the inertia matrix is defined as I x 0 −I xz 0 I 0 I = (17) y −I xz 0 I z
Expanding (15) using (16) results in M x M y M z
= pI ˙ x − rI ˙ xz + qr(I z − I y ) − pqI xz = qI ˙ y + pq(I x − I z ) + ( p2 − r2 )I xz = rI ˙ z − pI ˙ xz + pq(I y − I x ) + qrI xz .
(18)
The external moments M x , M y and M z are those due to aerodynamics and engine angular momentum. As a result the aerodynamic moments are 7
M x M y M z
¯ = L ¯ − rH eng = M ¯ + qH eng = N
(19)
¯ M ¯ and N ¯ are the aerodynamic moments and H eng is the engine where L, angular momentum. Note that the engine angular momentem is assumed to be parallel to the body x-axis of the aircraft. The aerodynamic moments can be expressed in a similar way as the aerodynamic forces in Equation (13): ¯ = q¯SbC lT (α,β,p,q,r,δ,...) L ¯ = q¯S ¯ M cC mT (α,β,p,q,r,δ,...) ¯ = q¯SbC nT (α,β,p,q,r,δ,...) N
(20)
Combining (18) and (19), the complete body-axis moment equation is formed as ¯ = pI L ˙ x − rI ˙ xz + qr(I z − I y ) − pqI xz ¯ − rH eng = qI M ˙ y + pq(I x − I z ) + ( p2 − r2 )I xz ¯ + qH eng = rI N ˙ z − pI ˙ xz + pq(I y − I x ) + qrI xz .
1.4 1.4.1
(21)
Gathering the Equations of Motion Euler Angles
The equations of motion derived in the previous sections are now collected and written as a system of twelve scalar first order differential equations. 1 ¯ (X + F T ) m 1 ¯ = pw − ru + g sin φ cos θ + Y m 1 ¯ = qu − pv + g cos φ cos θ + Z m
u˙ = rv − qw − g sin θ +
(22)
v˙
(23)
w˙
(24)
¯ + c4 (N ¯ + qH eng ) p˙ = (c1 r + c2 p)q + c3 L ¯ − rH eng ) q˙ = c5 pr − c6 ( p2 − r2 ) + c7 (M ¯ + c9 (N ¯ + qH eng ) r˙ = (c8 p − c2 r)q + c4 L
(25) (26) (27)
φ˙ = p + tan θ(q sin φ + r cos φ) θ˙ = q cos φ − r sin φ q sin φ + r cos φ ψ˙ = cos θ
(28) (29)
8
(30)
x˙ E
= + = + =
y˙ E z˙E
u cos ψ cos θ + v(cos ψ sin θ sin φ − sin ψ cos φ) w(cos ψ sin θ cos φ + sin ψ sin φ) u sin ψ cos θ + v(sin ψ sin θ sin φ + cos ψ cos φ) w(sin ψ sin θ cos φ − cos ψ sin φ) −u sin θ + v cos θ sin φ + w cos θ cos φ
(31) (32) (33)
where 2 Γc1 = (I y − I z )I z − I xz
Γc4 = I xz
Γc2 = (I x − I y + I z )I xz
c5 =
Γc3 = I z
c7 =
2 Γc8 = I x (I x − I y ) + I xz
I z −I x I y
c6 =
1 I y
I xz I y
Γc9 = I x
2 with Γ = I x I z − I xz .
1.4.2
Quaternions
The above equations of motion make use of Euler angle approach for the orientation model. The disadvantage of the Euler angle method is that the differential equations for p˙ and r˙ become singular when pitch angle θ passes through ± π2 . To avoid these singularities quaternions are used for the aircraft orientation presentation. A detailed explanation about quaternions and their properties can be found in [3]. With the quaternions presentation the aircraft system representation exists of 13 scalar first order differential equations: 1 ¯ X + F T + 2(q1 q3 − q0 q2 )g m 1 ¯ = pw − ru + Y + 2(q2 q3 + q0 q1 )g m 1 ¯ = qu − pv + Z + (q02 − q12 − q22 + q32 )g m
u˙ = rv − qw + v˙ w˙
(34) (35) (36)
¯ + c4 (N ¯ + qH eng ) p˙ = (c1 r + c2 p)q + c3 L ¯ − rH eng ) q˙ = c5 pr − c6 ( p2 − r2 ) + c7 (M ¯ + c9 (N ¯ + qH eng ) r˙ = (c8 p − c2 r)q + c4 L
q˙ q˙ ˙ = q˙
0
q
1 2
q˙3
x˙ E y˙ E z˙E
=
0 = 1 p 2 q
q02 + q12 − q22 − q32 2(q1 q2 + q0 q3 ) 2(q1 q3 − q0 q2 )
r
− p 0 −r q
−q −r r −q 0 p − p 0
2(q1 q2 − q0 q3 ) q02 − q12 + q22 − q32 2(q2 q3 + q0 q1 )
9
q q q
0 1 2
q3
(37) (38) (39)
2(q1 q3 + q0 q2 ) 2(q2 q3 − q0 q1 ) q02 − q12 − q22 + q32
(40)
u v w
(41)
where
q q q
0 1 2
q3
cos φ/2cos θ/2cos ψ/2 + sin φ/2sin θ/2sin ψ/2 = ± sin φ/2cos θ/2cos ψ/2 − cos φ/2sin θ/2sin ψ/2 cos φ/2sin θ/2cos ψ/2 + sin φ/2cos θ/2sin ψ/2 cos φ/2cos θ/2sin ψ/2 − sin φ/2sin θ/2cos ψ/2
Using (40) to describe the attitude dynamics means that the four differential equations are integrated as if all quaternion components were independent. Therefore, the normalization condition |q| = q02 + q12 + q22 + q32 = 1 and the derivative constraint q0 q˙0 + q1 q˙1 + q2 q˙2 + q3 q˙3 = 0 may not be satisfied after performing an integration step due to numerical round-off errors. After each integration step the constraint may be re-established by subtracting the discrepancy from the quaternion derivatives. The corrected quaternion dynamics are [6]
q˙ ′ = q˙ − δ q,
(42)
where δ = q0 q˙0 + q1 q˙1 + q2 q˙2 + q3 q˙3 . 1.4.3
Wind Axes Force Equations
For control design it is more convenient to transform the force equations (34)(36) to the wind-axes reference frame. Taking the derivative of Eq. (2) results in [3] 1
˙ T V
=
α˙
=
q − ( p cos α + r sin α)tan β −
β˙
=
p sin α − r cos α +
m
(−D + F T cos α cos β + mg1 )
1 mV T
1 (L + F T sin α − mg3 ) mV T cos β
(Y − F T cos α sin β + mg2 )
(43) (44) (45)
where the drag force D, the side force Y and the lift force L are defined as ¯ ¯ sin β − Z sin ¯ D = −X cos α cos β − Y α cos β ¯ ¯ cos β − Z sin ¯ Y = −X cos α sin β + Y α sin β ¯ ¯ L = X sin α − Z cos α and the gravity components as g1 g2 g3
2
= g (− cos α cos β sin θ + sin β sin φ cos θ + sin α cos β cos φ cos θ) = g (cos α sin β sin θ + cos β sin φ cos θ − sin α sin β cos φ cos θ) = g (sin α sin θ + cos α cos φ cos θ)
Control Variables and Engine Modeling
The F-16 model allows for control over thrust, elevator, ailerons, and rudder. The thrust is measured in Newtons. All deflections are defined positive in 10
the conventional way, i.e. positive thrust causes an increase in acceleration along the xB -axis, a positive elevator deflection results in a decrease in pitch rate, a positive aileron deflection gives a decrease in roll rate and a positive rudder deflection decreases the yaw rate. The F-16 also has a leading edge flap, which helps to fly the aircraft at high angles of attack. The deflection of the leading edge flap δLEF is not controlled directly by the pilot, but is governed by the following transfer function dependent on angle of attack α and static and dynamic pressures: δLEF = 1.38
2s + 7.25 q¯ α − 9.05 + 1.45 s + 7.25 ps
The differential elevator deflection, trailing edge flap, landing gear and speed brakes are not included in the model, since no data is publicly available. The control surfaces of the F-16 are driven by servo-controlled actuators to produce the deflections commanded by the flight control system. The actuators of the control surfaces are modeled as a first-order low-pass filters with certain gain and saturation limits in range and deflection rate. These limits can be found in Table 1. The gains of the actuators are 1/0.136 for the leading edge flap and 1/0.0495 for the other control surfaces. The maximum values and units for all control variables are given in Table 1. Table 1: The control input units and maximum values Control Elevator Ailerons Rudder Leading edge flap
units deg deg deg deg
MIN. -25 -21.5 -30 0
MAX. 25 21.5 30 25
rate limit ± 60 deg/s ± 80 deg/s ± 120 deg/s ± 25 deg/s
The Lockheed-Martin F-16A is powered by an after-burning turbofan jet engine, which is modeled taking into account throttle gearing and engine power level lag. The thrust response is modeled with a first order lag, where the lag time constant is a function of the current engine power level and the commanded power. The commanded power level to the throttle position is a linear relationship apart from a change in slope when the military power level is reached at 0.77 throttle setting [5]: P c (δth ) =
64.94δ
if δth ≤ 0.77 . 217.38δth − 117.38 if δth > 0.77 th
(46)
Note that the throttle position is limited to the range 0 ≤ δth ≤ 1. The derivative of the actual power level P a is given by [5] ˙a = P
1 (P c − P a ) , τ eng 11
(47)
where
P 60 P = 40P 5.0 1 = 5.0 τ c
c
c
1
if if if if
if if if if
∗
τ eng
eng
1 τ eng
∗
P c P c P c P c
1 ∗ τ eng
≥ 50 and ≤ 50 and < 50 and < 50 and P c P c P c P c
≥ 50 and ≤ 50 and < 50 and < 50 and
1.0 = 0.1 1.9 − 0.036(P − P ) c
P a P a P a P a
a
≥ 50 < 50 ≥ 50 < 50 P a P a P a P a
≥ 50 < 50 ≥ 50 < 50
if (P c − P a ) ≤ 25 if (P c − P a ) ≥ 50 . if 25 < (P c − P a ) < 50
The engine thrust data is available in a tabular form as a function of actual power, altitude and Mach number over the ranges 0 ≥ h ≥ 15240 m and 0 ≥ M ≥ 1 for idle, military and maximum power settings [5]. The thrust is computed as F T =
a + (T mil − T idle ) P 50 P a 50 T mil + (T max − T mil ) 50
T
idle
−
if P a < 50 . if P a ≥ 50
(48)
The engine angular momentum is assumed to be acting along the xB -axis with a constant value of 216.9 kg.m 2 /s.
3
Geometry and Aerodynamic Data
The relevant geometry data of the F-16A can be found in Table 2. The aerodynamic data of the F-16 model have been derived from low-speed static and dynamic (force oscillation) wind-tunnel tests conducted with sub-scale models in wind-tunnel facilities at the NASA Ames and Langley Research Centers [5]. The aerodynamic data in [5] is given in tabular form and is valid for the following subsonic flight envelope: • −20 ≤ α ≤ 90 degrees; • −30 ≤ β ≤ 30 degrees. Two examples of the aerodynamic data for the F-16 model can be found in Figure 3. The pitch moment coefficient C m and the C Z both depend on three variables: angle of attack, sideslip angle and elevator deflection. The various aerodynamic contributions to a given force or moment coefficient
12
0.2 2
0.1 0
1
−0.1 ) −−0.2 (
) − (
m −0.3 C
0
Z C
−1
−0.4 −2
−0.5 −0.6
20
−20
0
20
40
60
20
−3 −20
0
0
−20
80
0 20
40
60
80
beta (deg) alpha (deg)
(a)
C m
−20 beta (deg)
alpha (deg)
for
δe
=0
(b)
C Z
for
δe
=0
Figure 3: Two examples of the aerodynamic coefficient data for the F-16 obtained from wind-tunnel tests. are summed as follows. For the X-axis force coefficient C XT : C XT
δLEF 25 q¯ c δLEF C Xq (α) + δC XqLEF (α) 1 − 2V T 25
= C X (α,β,δe ) + δC XLEF 1 − +
where δC XLEF = C XLEF (α, β) − C X (α,β,δe = 0 o ) For the Y-axis force coefficient C YT : C YT
δ δC + δC 1 − δ 25 δ 25 20 δ δ rb δC + C (α) + δC (α) 1 − 30 2V 25 pb δ C (α) + δC (α) 1 −
= C Y (α, β) + δC YL EF 1 −
LEF
+
LEF
+ +
Yδ a
Y δr
2V T
Yδ a
a
LEF
r
T
Y p
Y r
Yr LEF
Yp LEF
LEF
LEF
25
where δC Y LEF δC Yδ a δC Y δa LEF
δC Yδ r
= C YL EF (α, β) − C Y (α, β) = C Yδ a (α, β) − C Y (α, β) = C Yδ a (α, β) − C Y LEF (α, β) − δC Yδ a LEF
= C Yδ r (α, β) − C Y (α, β)
13
For the Z-axis force coefficient C ZT : δLEF 25 q¯ c δLEF C Zq (α) + δC ZqLEF (α) 1 − 2V T 25
C ZT
= C Z (α,β,δe ) + δC ZLEF 1 −
+
where δC ZLEF = C ZLEF (α, β) − C Z (α,β,δe = 0 o ) For the rolling-moment coefficient C lT : C lT
= C l (α,β,δe ) + δC lLEF 1 −
δLEF 25
δ
1 − δ 25 20 δ δ rb δC + C (α) + δC (α) 1 − 30 2V 25 pb δ C (α) + δC (α) 1 − + δC (α)β δC
+
lδa
+
lδr
+
2V T
LEF
+ δC lδa
LEF
r
lr
T
lp
a
LEF
lrLEF
LEF
lpLEF
lβ
25
where δC lLEF δC lδa δC lδa LEF
δC lδr
= C lLEF (α, β) − C l (α,β,δe = 0 o) = C lδa (α, β) − C l (α,β,δe = 0 o ) = C lδa (α, β) − C lLEF (α, β) − δC lδa LEF
= C lδr (α, β) − C l (α,β,δe = 0 o )
For the pitching-moment coefficient C mT : C mT
= C m (α,β,δe ) + C ZT +
x
cgr
q¯ c C mq (α) + δC mqLEF 2V T
δ −x + δC 1− 25 δ (α) 1 − + δC (α) + δC cg
LEF
mLEF
LEF
25
m
mds (α, δe )
where δC mLEF = C mLEF (α, β) − C m (α,β,δe = 0 o ) For the yawing-moment coefficient C nT : C nT
δLEF c¯ − C Y T xcgr − xcg 25 b δLEF δa + δC nδa + δC nδa 1− LEF 25 20 δr rb δLEF + δC nδr + C nr (α) + δC nrLEF (α) 1 − 30 2V T 25 pb δLEF + C np (α) + δC npLEF (α) 1 − + δC nβ (α)β 2V T 25
= C n (α,β,δe ) + δC nLEF 1 −
14
where δC nLEF δC nδa δC nδaLEF δC nδr
4
= C nLEF (α, β) − C n (α,β,δe = 0 o ) = C nδa (α, β) − C n (α,β,δe = 0 o ) = C nδaLEF (α, β) − C nLEF (α, β) − δC nδa = C nδr (α, β) − C n (α,β,δe = 0 o )
MATLAB/Simulink script files
Descriptions are included in the files. Most files are based on or copied from the F-16 models by R. S. Russell [7] and Y. Huo. The S-functions F_16dyn.c
and F_16dynam.c
both contain the aircraft dynamics, they only differ in the use of quaternions vs. Euler angles. the runF16model.m
script file allows the user to trim the aircraft model for steady-state flight conditions. Other files: tgear.m trim_fun.m lofi_f16_aerodata.c hifi_f16_aerodata.c engine_model.c ISA_atmos.c mexndinterp.c
15
References [1] J. H. Blakelock. Automatic Control of Aircraft and Missiles . John Wiley & Sons, 2nd edition, 1991. [2] M. V. Cook. Flight Dynamics Principles, pages 11–29. Heinemann, 1997.
Butterworth-
[3] B. L. Lewis and F. L. Stevens. Aircraft Control and Simulation , pages 1– 54,110–115. John Wiley & Sons, 1992. [4] J. A. Mulder, W. H. J. J. van Staveren, J. C. van der Vaart, and E. de Weerdt. Flight dynamics, lecture notes ae3-302. Technical report, Delft University of Technology, 2006. [5] L. T. Nguyen, M. E. Ogburn, W. P. Gilbert, K. S. Kibler, P. W. Brown, and P. L. Deal. Simulator study of stall post-stall characteristics of a fighter airplane with relaxed longitudinal static stability. Technical report, NASA Langley Research Center, 1979. [6] K. Refson. Moldy’s User’s Manual . Department of Earth Sciences, May 2001. [7] R. S. Russell. Nonlinear f-16 simulations using simulink and matlab. Technical report, University of Minnesota, 2003.
16
A
F-16 Geometry Table 2: F-16 parameters. Parameter aircraft mass (kg) wing span (m) wing area (m2 ) mean aerodynamic chord (m) roll moment of inertia (kg .m2 ) pitch moment of inertia (kg.m2 ) yaw moment of inertia (kg.m2 ) product moment of inertia (kg.m2 ) product moment of inertia (kg.m2 ) product moment of inertia (kg.m2 ) c.g. location (m) reference c.g. location (m) engine angular momentum (kg.m 2 /s)
B
Symbol m b S ¯c I x I y I z I xz I xy I yz xcg xcgr H eng
Value 9295.44 9.144 27.87 3.45 12874.8 75673.6 85552.1 1331.4 0.0 0.0 0.3¯ c 0.35¯ c 216.9
ISA Atmospheric Model
For the atmospheric data an approximation of the International Standard Atmosphere (ISA) is used [4]. T = p
=
T + λh if T if p 1 + λ p e p RT γRT 0
(h=11000) 0
h T 0
−
(h=11000)
ρ = a =
−
h ≤ 11000 h > 11000 g0 Rλ
if h ≤ 11000
g
RT (h=11000)
(h−11000)
if h > 11000
where T 0 = 288.15 K is the temperature at sea level, p0 = 101325 N/m2 the pressure at sea level, R = 287.05 J/kg.K the gas constant of air, g0 = 9.80665 m/s2 the gravity constant at sea level, λ = dT/dh = −0.0065 K/m the temperature gradient and γ = 1.41 the isentropic expansion factor for air. Given the aircraft’s altitude (h in meters) it returns the current temperature ( T in Kelvin), the current air pressure ( p in N/m2 ), the current air density (ρ in kg/m3 ) and the speed of sound (a in m/s).
17
