Aerodynamic coefficients as a function of angle of attack.
This chapter is the first of two others that will follow (a three-chapter series). Here we present the derivation of the mathematical model for a rocket’s autopilots in state space. The basic equations defining the airframe dynamics of a typical six degrees of freedom (6DoFs) are nonlinear and coupled. Separation of these nonlinear coupled dynamics is presented in this chapter to isolate the lateral dynamics from the longitudinal dynamics. Also, the need to determine aerodynamic coefficients and their derivative components is brought to light here. This is the crux of the equation. Methods of obtaining such coefficients and their derivatives in a sequential form are also put forward. After the aerodynamic coefficients and their derivatives are obtained, the next step is to trim and linearize the decoupled nonlinear 6DoFs. In a novel way, we presented the linearization of the decoupled 6DoF equations in a generalized form. This should provide a lucid and easy way to implement trim and linearization in a computer program. The longitudinal model of a rocket presented in this chapter will serve as the main mathematical model in two other chapters that follow in this book.
- six degrees of freedom (6DoFs)
- state space
Over the years, a number of authors [1, 2, 3, 4] have considered modeling rocket/missile dynamics for atmospheric flights. In the majority of the published work on these mathematical models, trimming and locally linearization are done without detailed explanation to the variables in the decoupled airframe dynamics. As such, the easy-to-write computer programs to facilitate this process numerically have been impeded.
With the advent of fast processors and numerical software like MATLAB®, Maple, Python, etc., it is now possible to take a complex nonlinear 6DoF equation like that of a rocket and run a program that can trim and linearize it with ease. This has made the field of control system design to grow at an exponential rate .
It is a known fact that the mathematical models are developed with their use in mind. This means before delving into the realization of a model, one should be well informed of the purpose for which that mathematical model will serve . Especially, in the field of control system design, a mathematical model in transfer function might not be ideal for optimal control design. However, problem-solving environments (PSEs) like MATLAB®/Simulink come with built-in functions capable of transforming, say state-space model, to transfer functions. One should bear in mind that this is not without a “cost.”
2. Mathematical model
A nontrivial part of any control problem is modeling the process. The objective is to obtain the simplest mathematical description that adequately predicts the response of the physical system to all inputs. For a rigid dynamic body, its motion can be described in translational, rotational, and angular inclinations at all times.
2.1. Translational motion
An accelerometer is often used to measure force on a dynamic body. For a rocket in motion, these forces  are represented as given in (1). Actually, this is a measure of the specific force, i.e., the nongravitational force per unit mass in x,y,z-directions, respectively. The specific force (also called the g-force or mass-specific force) has units of acceleration or ms−1. So, it is not actually a force at all but a type of acceleration:
components of aerodynamic force vector FA expressed in the body.
coordinate system, N.
components of gravitational force vector Fg expressed in body coordinate system, N.
components of thrust vector Fp expressed in the body coordinate system, N.
m = instantaneous rocket mass, kg.
p, q, r = components of angular rate vector ω expressed in body coordinate system (roll, pitch, and yaw, respectively), rad/s.
u, v, w = components of absolute linear velocity vector Vx expressed in the body coordinate system, m/s.
components of linear acceleration expressed in body coordinate system, m/s2.
The aerodynamic forces have the following components:
CA = aerodynamic axial force coefficient, dimensionless.
CN = aerodynamic normal force coefficient, dimensionless.
CNy = coefficient corresponding to component of normal force on yb-axis
CNz = coefficient corresponding to component of normal force on zb-axis
S = aerodynamic reference area, m2.
VM = magnitude of velocity vector of the center of mass of the rocket, m/s.
ρ = atmospheric density, kg/m3.
Fp = magnitude of total instantaneous thrust force vector, N.
Fpxb, Fpyb, Fpzb = components of thrust vector Fp expressed in body coordinate system, N.
γ1 = angle measured from xb-axis to projection of thrust vector Fp on xb yb-plane, rad (deg).
γ2 = angle measured from projection of thrust vector Fp on xb yb-plane to the thrust vector Fp, rad.
lp = the distance from the aerodynamic center to center of mass, m.
l = the distance from the center of mass to nozzle, m.
If the rocket was designed for thrust vector control via gimbaling of the nozzle, Fp will be computed as given in (3). Here we assume that Fp is acting along the line of symmetry of the rocket because the nozzle is fixed (fin control). Hence, the magnitude of the thrust force Fp is calculated by
Ae = rocket nozzle exit area, m2.
Fpref = reference thrust force, N.
Pa = ambient atmospheric pressure, Pa.
Pref = reference ambient pressure, Pa.
The gravitational forces in (1) are computed as follows:
Fgxe, Fgye, Fgze = components of gravitational force vector Fg expressed in earth coordinate system, N.
g = acceleration due to gravity, m/s2.
m = instantaneous mass of rocket, kg.
The dependence of the acceleration due to gravity on the altitude of the rocket is given by
g = acceleration due to gravity, m/s2.
g0 = acceleration due to gravity at earth surface (nominally 9.8 m/s2), m/s2.
h = altitude above sea level, m.
Re = radius of the earth, m.
c = cosine function (cθ = cos θ), dimensionless.
s = sine function (sθ = sin θ), dimensionless.
θ = the Euler angle of rotation in elevation (pitch) of body frame relative to earth frame, rad (deg).
φ = the Euler angle of rotation in roll of body frame relative to earth frame, rad (deg).
ψ = the Euler angle of rotation in azimuth (heading) of body frame relative to earth frame, rad (deg).
[Te/b] = transformation matrix from earth to body.
A vector v expressed in the body coordinate system can be transformed to the earth coordinate system by the matrix equation
Hence, considering (5) we can write the following:
The terms Fgxb, Fgyb, and Fgzb are the components of the gravitational force substituted into (1).
The mass in (1) is given below:
Fpref = reference thrust force, N.
Isp = specific impulse of propellant, Ns/kg.
m0 = rocket mass at time zero (i.e., at the time of launch), kg.
t = simulated time, s.
2.2. Rotational motion
A gyroscope or gyro is a device that measures the angular acceleration or rotational motion of a dynamic body. On a rocket, this rotational motion can be described as
LA, MA, NA = components of aerodynamic moment vector MA expressed in body coordinate system (roll, pitch, and yaw, respectively), Nm.
Lp, Mp, Np = components of propulsion moment vector Mp expressed in body coordinate system (roll, pitch, and yaw, respectively), Nm.
Ix, Iy, Iz = components of inertia (diagonal elements of inertia matrix when products of inertia are zero), kg-m2.
p, q, r (P, Q, R) = components of angular rate vector ω expressed in body coordinate system (roll, pitch, and yaw, respectively), rad/s (deg/s).
= components of angular acceleration expressed in body coordinate (roll, pitch, and yaw, respectively), rad/s2 (deg/s2).
Cl = aerodynamic roll moment coefficient about center of mass, dimensionless.
Cm = aerodynamic pitch moment coefficient about center of mass, dimensionless.
Cn = aerodynamic yaw moment coefficient about center of mass, dimensionless.
d = aerodynamic reference length of body, m.
The aerodynamic moment coefficients are of the form
Clp = roll damping derivative relative to roll rate p, rad−1 (deg−1).
Clδ = slope of curve formed by roll moment coefficient C1 versus control surface deflection, rad−1(deg−1).
Cmref = pitching moment coefficient about reference moment station, dimensionless.
Cmq = pitch damping derivatives relative to pitch rate q, rad−1(deg−1).
= pitch damping derivative relative to angle of attack rate (slope of curve formed by Cα verses ), rad−1 (deg−1).
CNy = coefficient corresponding to component of normal force on yb-axis, dimensionless.
CNz = coefficient corresponding to component of normal force on zb-axis, dimensionless.
= yaw damping derivative relative to yaw rate , rad−1(deg−1).
Cnref = yawing moment coefficient about reference moment station, dimensionless.
= yaw damping derivative relative to angle of sideslip rate , rad-l (deg-1).
xcm = instantaneous distance from rocket nose to center of mass, m.
xref = distance from rocket nose to reference moment station, m.
δr = effective control surface deflection causing rolling moment, rad (deg).
Considering that the rocket we intend to control is via fin deflection, fins on the rocket will be designated as shown in Figure 2.
Hence, the following moment coefficients are also given as
Cmref = pitching moment coefficient about reference moment station (this is the static value normally measured in the wind tunnel.), dimensionless.
= slope of curve formed by pitch moment coefficient. Cm versus angle of attack, rad−1 (deg−1) slope of tune formed.
= slope of tune formed by pitch moment coefficient Cm versus control surface deflection for pitch, δp rad−1 (deg−1).
= slope of curve formed by yawing moment coefficient Cn versus angle of sideslip β, rad−1 (deg−1).
= slope of curve formed by yaw moment coefficient Cn versus effective control surface deflection for yaw, δy rad−1 (deg−1).
= angle of attack, rad (deg).
β = angle of sideslip, rad (deg).
δη = δ1 = δ3 effective control surface deflection causing pitching moment, rad (deg).
δζ = δ4 = δ2 = effective control surface deflection causing yawing moment, rad (deg).
The angle of attack, angle of sideslip, and roll angle required for the realization of the aerodynamic coefficients are
where ∅ is aerodynamic roll angle, rad (deg), and αt is the total angle of attack.
Table 1 gives a list of the aerodynamic coefficients that must be obtained for every rocket design before a model can be realized. There exist numerical and semi-numerical means of obtaining such coefficients. Examples of software that can do semiempirical computation of such coefficients and their derivatives are Missile Digital DATCOM®  and Flexible Structures Simulator (FSS) . Finally, a wind tunnel test is expected to validate and update all coefficients and their derivatives before the system engineer delves in the control design.
|1||CN||Normal force coefficient (body axis)|
|2||CL||Lift coefficient (wind axis)|
|3||CM||Pitching moment coefficient (body axis)|
|4||Xcp||Center of pressure in calibers from moment reference center|
|5||CA||Axial force coefficient (body axis)|
|6||CD||Drag coefficient (wind axis)|
|7||CY||Side force coefficient (body axis)|
|8||Cn||Yawing moment coefficient (body axis)|
|9||Cl||Rolling moment coefficient (body axis)|
|10||CNα||Normal force coefficient derivative with angle of attack|
|11||CMα||Pitching moment coefficient derivative with angle of attack|
|12||CYβ||Side force coefficient derivative with sideslip angle|
|13||Cnβ||Yawing moment coefficient derivative with sideslip angle (body axis)|
|14||Clβ||Rolling moment coefficient derivative with sideslip angle (body axis)|
|15||CMq||Pitching moment coefficient derivative with pitch rate|
|16||CNq||Normal force coefficient derivative with pitch rate|
|17||CAq||Axial force coefficient derivative with pitch rate|
|18||Pitching moment derivative with rate of change of angle of attack|
|19||Normal force derivative with rate of change of angle of attack|
|20||Clp||Rolling moment coefficient derivative with roll rate|
|21||Cnp||Yawing moment coefficient derivative with roll rate|
|22||CYp||Side force coefficient derivative with roll rate|
|23||Clr||Rolling moment coefficient derivative with yaw rate|
|24||Cnr||Yawing moment coefficient derivative with yaw rate|
|25||CYr||Side force coefficient derivative with yaw rate|
The third and final component to fully describe the motion of a rocket is its angular inclination or attitude. We chose the Euler angles to describe the attitude of the rocket.
2.3. The Euler angles
Missile attitude is required for a number of simulation functions including the calculation of angle of attack, seeker gimbal angles, fuze look angles, and warhead spray pattern. In simulations with six degrees of freedom, the missile attitude is calculated directly by integrating the set of equations that define the Euler angle rates, i.e.,
= the Euler angle rotation in elevation (pitch angle), rad (deg).
= the Euler angle rotation in roll (roll angle), rad (deg).
= rates of change of the Euler angles in heading, pitch, and roll, respectively, rad/s (deg/s).
P, q, r = components of angular rate vector w expressed in body coordinate system (roll, pitch, and yaw, respectively), rad/s.
The three heading angle of can be measured directly using a magnetometer. Note also that the Euler angles in (16) can all be derived by integrating gyroscopic measurements. As such we might not need an instrument that will measure it directly. Nevertheless, a magnetometer can be added to the instrumentation on board to measure heading.
Combining (1), (11), and (16) gives a complete six-degree-of-freedom equation of motion for a rocket in flight as shown in Figure 3. This could be written together as given in (17). In today’s modern aerospace industry, a single device like the MPU6050, a MEM-based integrated chip, can be used to give numerical values for state variables of (17) on any dynamic body. For control system design, the rocket system as described in (17) needs to be separated into the two planes (decouple); these are called the lateral (la) and longitudinal (lo) dynamic equations of motion:
Since the control system we intend to design is in a family of linear controllers and (17) is a nonlinear system of differential equations, trimming and linearization must be done.
To trim or find equilibrium values requires a good knowledge of advance computational techniques. A trim point, also known as an equilibrium point, is a point in the parameter space of a dynamic system at which the system is in a steady state. The trim problem for a rocket can be described as finding a set of suitable input values to satisfy a set of conditions. Hence, a trim point involves setting of its controls that causes the rocket to fly straight and level in the longitudinal plane. The suitable input values are the control surface deflections, the thrust setting, and the rocket attitude . The set of conditions are the rocket’s accelerations. The variables associated with the trim problems can be divided into three categories:
Flight condition variables
The objective variables need to be driven toward the specified values, often zero (i.e., steady flight with zero sideslip). The objective parameters are combined in the objective vectors ola (state vectors) and olo, for the lateral and longitudinal missile dynamics as
The sideslip angle is also included, since for most cases, there are multiple solutions to the trim problem, each with a different sideslip angle. In the desired solution, the sideslip angle should be zero. In that case, the drag is at a minimum. The control parameters are adjusted in order to drive the objective parameters to their specified values. Together, they form the control vectors cla and clo, described in (22). The control variables (input variables) describe the trimmed pilot control and the aircraft attitude:
Finally, the 12 states of the 6DoF equation of motion must be initialized with the initial state conditions. In MATLAB®, the trim command is used to find equilibrium points. The object of trimming is to bring the forces and moments acting on the rocket into a state of equilibrium. That is the condition when the axial, normal, and side forces and the roll, pitch, and yaw moments are all zero. Mathematically, trimming combines implicit and explicit Jacobian approaches. The Jacobian trim approach is the preferred method for rigid rocket. The Jacobian method is robust, since trim convergence is likely to occur even with rough estimates of the Jacobian and a rough first guess. Note that in general, any optimization routine could be used to solve the trim problem, as long as it is robust enough.
The Jacobian approach is based on the assumption that the change of the objective vector is linearly related to the change in the control input, which is shown in (22):
In (22), Ji is the Jacobian matrix evaluated near control input ci. Its entries are first-order partial derivatives and represent the effect of changes in each control input on acceleration.
Note that changes in lateral plane induce changes in the longitudinal plane and vice versa; thus, we can write (23) the Jacobian for the lateral dynamics and (24) for the longitudinal or pitch dynamics:
If the rocket is assumed to be at equilibrium, or trim condition, then the equations of motion can be linearized, and the 6DoF equation of motion can be resolved into their lateral and longitudinal states.
The system of first-order nonlinear differential equations of the rocket as presented by (17) is said to be in state variable form if its mathematical model is described by a system of n first-order differential equations and an algebraic output equation as :
The column vector x = [x1,…xn]T is called the state of the system. The scalars u and y are called the control input and the system output, respectively, denoting
Concisely, (26) is written as
where f and h are nonlinear functions of x and u; then, we say that the system is nonlinear. To linearize (26), we desire it to become
where A is n x n, B is n x 1, C is 1 x n, and D is all scalar in MATLAB®/Simulink.
One reason for approximating the nonlinear system (26) by a linear model of the form (28) is that, by so doing, one can apply rather simple and systematic linear control design techniques. Given the nonlinear system (26) and an equilibrium or trimmed points x∗ = […]T obtained when u = u∗, noting that we define a coordinate transformation as follows:
Further, denoting, . The new coordinates and represent the variations of x, u, and y from their equilibrium values. You have to think of these as a new state, new control input, and new output, respectively. The linearization of (26) is at x∗ in which the equilibrium or trim  state is given by
Note that the matrices A, B, C, and D have constant coefficients in that all partial derivatives are evaluated at the numerical values .
2.6. Lateral dynamics of a rocket
Equations of motion in the lateral plane are described by (30). Note that (30) comprises of one of the force equations (Fy), one of the momentum equations (My), and two of the Euler angles from the 6DoF equations (decoupled from (17)):
All the variables in A matrix of (31) are the lateral dimensionless aerodynamic stability derivatives with respect to the system state vectors. The variables in B matrix of (31) are the lateral dimensionless aerodynamic control derivatives with respect to the designated control surfaces.
2.7. Longitudinal dynamics of a rocket
The longitudinal dynamics in motion is also called the pitch plane. Equations describing the motion of a rocket in this plane can be described as given in (32). Note that (32) comprises two of the force equations (Fx and Fz), two of the momentum equations (Mx and Mz), and two of the Euler angles from the 6DoF equations as given in (17):
In MATLAB® the linmod  command is used to invoke linearization. The assumption made for decoupling the linear model is that the cross-coupling effects between the two modes are negligible. These assumptions are:
The rocket is designed with conventional control surfaces that do not give significant cross-coupling control between the lateral and longitudinal modes.
The rocket is symmetrical about the xz plane in which the inertia cross-coupling in (xy and xz planes) result to cross-coupling between the lateral and longitudinal modes is minimum.
It can be shown that a typical trimmed and linearized model of the pitch plane (longitudinal dynamics) for a rocket  is given as presented in (34). Notice that compared to (33), the velocity in x-direction is not included in (34). This is basically due to the fact that in this pitch plane, translational motion for a rocket is predominantly in the z-direction (velocity w):
3. Discussion of result
From (34), it can be seen that a three state variable models have been realized in state space. This implies that modern observer like the Kalman filter can be designed to estimate and predict the trajectory of such rocket dynamics. This mathematical model also can be used to design all the control algorithms that fall in the class of modern (optimal theory) control. Particularly, this model is important in the realization of the longitudinal autopilot system of a rocket.
Mathematical models are the bedrock of almost all scientific activities. Here we were able to define the nonlinear airframe dynamics of a rocket in 6DoF. We went further to decouple the 6DoF equations of motion for the rocket and presented forms in which the decoupled 6DoF is easily trimmed and linearized with a computer program like MATLAB®. The process presented here can be repeated for any size of the rocket and aircrafts/unmanned aerial vehicle (UAV). Note that if the aerospace vehicle being model is not a rocket, and a state-space model is needed, all the procedures outlined in this chapter will be the same. The only changes that will be accommodated will come from the numerical values of the aerodynamic coefficients and their derivatives.