**State-Space Modeling of a Rocket for Optimal Control System Design**

Aliyu Bhar Kisabo and Aliyu Funmilayo Adebimpe

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/intechopen.82292

#### Abstract

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.

Keywords: rocket, six degrees of freedom (6DoFs), state space, trimming, linearization

## 1. Introduction

Over the years, a number of authors [1–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.

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 [5].

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 [6]. 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 [7] 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 �<sup>1</sup> (also called the g-force or mass-specific force) has units of acceleration or ms . So, it is not actually a force at all but a type of acceleration:

$$\begin{aligned} \dot{u} &= \frac{F\_{A\_{x\_b}} + F\_{P\_{x\_b}} + F\_{\mathcal{S}\_{x\_b}}}{m} - (qw - rv), m/s^2\\ \dot{v} &= \frac{F\_{A\_{y\_b}} + F\_{P\_{y\_b}} + F\_{\mathcal{S}\_{y\_b}}}{m} - (ru - pvv), m/s^2\\ \dot{w} &= \frac{F\_{A\_{z\_b}} + F\_{P\_{z\_b}} + F\_{\mathcal{S}\_{z\_b}}}{m} - (pv - qu), m/s^2 \end{aligned} \tag{1}$$

where

FAxb , FAyb , FAzb ¼ components of aerodynamic force vector FA expressed in the body.

coordinate system, N.

Fg ¼ components of gravitational force vector Fg expressed in body coordinate xb , Fgyb , Fgzb system, N.

, Fpyb , Fp Fp zb ¼ components of thrust vector Fp expressed in the body coordinate system, N. xb

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.

u, \_ v,\_ <sup>w</sup>\_ <sup>¼</sup> components of linear acceleration expressed in body coordinate system, m/s<sup>2</sup> .

The aerodynamic forces have the following components:

$$\begin{aligned} F\_{A\_{\overline{\nu}}} &= -0.5 \rho V\_M^2 \mathcal{C}\_A \mathcal{S} \\ F\_{A\_{\overline{\nu}}} &= 0.5 \rho V\_M^2 \mathcal{C}\_{N\_{\overline{\nu}}} \mathcal{S} \\ F\_{A\_{\overline{\nu}}} &= 0.5 \rho V\_M^2 \mathcal{C}\_{N\_{\overline{\nu}}} \mathcal{S} \end{aligned} \tag{2}$$

where

CA = aerodynamic axial force coefficient, dimensionless.

CN = aerodynamic normal force coefficient, dimensionless.

CNy = coefficient corresponding to component of normal force on yb-axis

$$\mathbf{C}\_{N\_y} = \mathbf{C}\_N \frac{-\upsilon}{\sqrt{\upsilon^2 + w^2}}, \text{ dimensionless}$$

CNz = coefficient corresponding to component of normal force on zb-axis

$$\mathbf{C}\_{N\_{\varepsilon}} = \mathbf{C}\_{N} \frac{1}{\sqrt{\upsilon^{2} + w^{2}}}, \text{ dimensionless}$$

<sup>2</sup> S = aerodynamic reference area, m .

VM = magnitude of velocity vector of the center of mass of the rocket, m/s.

r = atmospheric density, kg/m<sup>3</sup> .

The propulsive forces in (1), as depicted in (Figure 1), are computed [8] as follows:

$$\begin{aligned} F\_{p\_{z\_b}} &= F\_p \cos \gamma\_2 \cos \gamma\_{1'} N \\ F\_{p\_{y\_b}} &= F\_p \cos \gamma\_2 \sin \gamma\_{1'} N \\ F\_{p\_{z\_b}} &= -F\_p \sin \gamma\_{2'} N \end{aligned} \tag{3}$$

where

Fp = magnitude of total instantaneous thrust force vector, N.

Fpxb, Fpyb, Fpzb = components of thrust vector F<sup>p</sup> expressed in body coordinate system, N. γ<sup>1</sup> <sup>=</sup> angle measured from xb-axis to projection of thrust vector Fp on xb yb-plane, rad (deg).

Figure 1. Propulsion force from the nozzle of a rocket.

γ<sup>2</sup> = 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, F<sup>p</sup> will be computed as given in (3). Here we assume that F<sup>p</sup> 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

$$F\_p = F\_{pref} + \left(p\_{ref} - p\_a\right) A\_{e\prime} N \tag{4}$$

where

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:

$$\begin{aligned} F\_{\mathbb{S}\_{\mathbb{S}\_{\varepsilon}}} &= 0, N \\ F\_{\mathbb{S}\_{\mathbb{S}\_{\varepsilon}}} &= 0, N \\ F\_{\mathbb{S}\_{\mathbb{S}\_{\varepsilon}}} &= m \mathbb{g}, N \end{aligned} \tag{5}$$

where

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

$$\mathbf{g} = \mathbf{g}\_0 \left[ \frac{\mathbf{R}\_e^2}{\left(\mathbf{R}\_e + h\right)^2} \right] \text{, m/s}^2 \tag{6}$$

where

g = acceleration due to gravity, m/s2 .

g0 = acceleration due to gravity at earth surface (nominally 9.8 m/s<sup>2</sup> ), m/s<sup>2</sup> .

h = altitude above sea level, m.

Re = radius of the earth, m.

The gravitational force expressed in body coordinates is computed by multiplying (5) by matrix (7):

$$\begin{aligned} \left[ \mathbf{T}\_{\epsilon/\mathsf{b}} \right] = \begin{bmatrix} c\theta c\psi & s\phi s\theta c\psi - c\phi s\psi & c\phi s\theta c\psi + s\phi s\psi \\ c\theta s\psi & s\phi s\theta s\psi + c\phi c\psi & c\phi s\theta c\psi - s\phi c\psi \\ -s\theta & s\phi c\theta & c\phi c\theta \end{bmatrix} . \text{dimensionless} \tag{7} \end{aligned} \tag{8}$$

where

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

$$
\boldsymbol{\upsilon}\_{\boldsymbol{e}} = \left[ \boldsymbol{T}\_{\boldsymbol{e}/\boldsymbol{b}} \right] \boldsymbol{\upsilon}\_{\boldsymbol{b}} \tag{8}
$$

Hence, considering (5) we can write the following:

$$
\begin{bmatrix} F\_{\mathcal{S}\_{\mathfrak{z}\_b}} \\ F\_{\mathcal{S}\_{\mathfrak{z}\_b}} \\ F\_{\mathcal{S}\_{\mathfrak{z}\_b}} \end{bmatrix} = \begin{bmatrix} T\_{b/\epsilon} \end{bmatrix} \begin{bmatrix} F\_{\mathcal{S}\_{\mathfrak{z}\_\ell}} \\ F\_{\mathcal{S}\_{\mathfrak{z}\_\ell}} \\ F\_{\mathcal{S}\_{\mathfrak{z}\_\ell}} \end{bmatrix} / N \tag{9}
$$

The terms Fgxb, Fgyb, and Fgzb are the components of the gravitational force substituted into (1). The mass in (1) is given below:

$$m = m\_0 - \frac{1}{I\_{sp}} \int\_0^t F\_{p\_{\text{ref}}} dt \, \text{kg} \tag{10}$$

where

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

$$\dot{\eta} = \frac{L\_A + L\_p - qr(I\_z - I\_y)}{I\_x}, rad/s^2 \left(\text{deg}/s^2\right)$$

$$\dot{\eta} = \frac{M\_A + M\_p - rp(I\_x - I\_z)}{I\_y}, rad/s^2 \left(\text{deg}/s^2\right) \tag{11}$$

$$\dot{r} = \frac{N\_A + N\_p - pq(I\_y - I\_x)}{I\_z}, rad/s^2 \left(\text{deg}/s^2\right)$$

where

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-m<sup>2</sup> .

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).

p,\_ q,\_ r = \_ components of angular acceleration ω\_ expressed in body coordinate (roll, pitch, and yaw, respectively), rad/s<sup>2</sup> (deg/s<sup>2</sup> ).

$$L\_A = 0.5 \rho V\_M^2 \text{C}\_l \text{S}d, Nm$$

$$M\_A = 0.5 \rho V\_M^2 \text{C}\_m \text{S}d, Nm \tag{12}$$

$$N\_A = 0.5 \rho V\_M^2 \text{C}\_n \text{S}d, Nm$$

where

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

$$\begin{aligned} \mathbf{C}\_{l} &= \mathbf{C}\_{l}\boldsymbol{\delta}\_{r} + \frac{d}{2V\_{M}} \left( \mathbf{C}\_{l} \boldsymbol{p} \right) \\ \mathbf{C}\_{m} &= \mathbf{C}\_{m\_{\eta f}} - \mathbf{C}\_{N\_{z}} \frac{\mathbf{x}\_{cm} - \mathbf{x}\_{ref}}{d} + \frac{d}{2V\_{M}} \left( \mathbf{C}\_{m\_{\eta}} + \mathbf{C}\_{m\_{u}} \right) \boldsymbol{q} \text{ dimensionless} \\ \mathbf{C}\_{n} &= \mathbf{C}\_{n\_{\eta f}} + \mathbf{C}\_{N\_{y}} \frac{\mathbf{x}\_{cm} - \mathbf{x}\_{ref}}{d} + \frac{d}{2V\_{M}} \left( \mathbf{C}\_{n\_{r}} + \mathbf{C}\_{n\_{\theta}} \right) \boldsymbol{r} \end{aligned} \tag{13}$$

where

Clp = roll damping derivative relative to roll rate p, rad�<sup>1</sup> (deg�<sup>1</sup> ).

Cl<sup>δ</sup> = slope of curve formed by roll moment coefficient C1 versus control surface deflection, rad�<sup>1</sup> (deg�<sup>1</sup> ).

Cmref = pitching moment coefficient about reference moment station, dimensionless.

Cmq = pitch damping derivatives relative to pitch rate q, rad�<sup>1</sup> (deg�<sup>1</sup> ).

Cm \_ <sup>α</sup>\_ α (slope of curve formed by C<sup>α</sup> verses α), rad�<sup>1</sup> (deg�<sup>1</sup> ). = pitch damping derivative relative to angle of attack rate

CNy = coefficient corresponding to component of normal force on yb-axis, dimensionless.

CNz = coefficient corresponding to component of normal force on zb-axis, dimensionless.

\_ <sup>r</sup>\_ r, rad�<sup>1</sup> (deg�<sup>1</sup> ). Cn = yaw damping derivative relative to yaw rate

Cnref = yawing moment coefficient about reference moment station, dimensionless.

Cn <sup>=</sup> yaw damping derivative relative to angle of sideslip rate \_ \_ <sup>β</sup>, rad-l (deg-1). <sup>β</sup>

xcm = instantaneous distance from rocket nose to center of mass, m.

xref = distance from rocket nose to reference moment station, m.

δ<sup>r</sup> = 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

$$\begin{aligned} \mathsf{C}\_{m\_{\text{ref}}} &= \mathsf{C}\_{m\_{\text{a}}} \alpha + \mathsf{C}\_{m\_{\text{b}}} \delta\_{\eta} \\ \mathsf{C}\_{n\_{\text{rf}}} &= \mathsf{C}\_{n\_{\text{f}}} \beta + \mathsf{C}\_{n\_{\text{b}}} \delta\_{\zeta} \end{aligned}, \text{dimensionless} \tag{14}$$

where

Cmref = pitching moment coefficient about reference moment station (this is the static value normally measured in the wind tunnel.), dimensionless.

Cm<sup>α</sup> = slope of curve formed by pitch moment coefficient. Cm versus angle of attack, α rad�<sup>1</sup> (deg�<sup>1</sup> ) slope of tune formed.

Cm<sup>δ</sup> = slope of tune formed by pitch moment coefficient Cm versus control surface deflection for pitch, δ<sup>p</sup> rad�<sup>1</sup> (deg�<sup>1</sup> ).

Cn<sup>β</sup> = slope of curve formed by yawing moment coefficient Cn versus angle of sideslip β, rad�<sup>1</sup> (deg�<sup>1</sup> ).

Cn<sup>δ</sup> = slope of curve formed by yaw moment coefficient Cn versus effective control surface deflection for yaw, δ<sup>y</sup> rad�<sup>1</sup> (deg�<sup>1</sup> ).

Figure 2. Fin control and designation for control.

α = angle of attack, rad (deg).

β = angle of sideslip, rad (deg).

δη = δ<sup>1</sup> = δ<sup>3</sup> effective control surface deflection causing pitching moment, rad (deg).

δζ = δ<sup>4</sup> = δ<sup>2</sup> = 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

$$\begin{aligned} \alpha &= \operatorname{Tan}^{-1}\left(\frac{w}{\mu}\right) \text{or } \alpha\_t = \tan^{-1}\left[\sqrt{\left(\frac{v^2 + w^2}{\mu}\right)}\right] \\ \beta &= \sin^{-1}\left(\frac{v}{V\_M}\right) \\ \phi &= \tan^{-1}\left(\frac{v}{w}\right) \end{aligned} \tag{15}$$

where ∅ is aerodynamic roll angle, rad (deg), and α<sup>t</sup> 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® [9] and Flexible Structures Simulator (FSS) [10]. 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.

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.,

$$\begin{aligned} \dot{\phi} &= p + (q \sin \phi + r \cos \phi) \tan \theta, rad/s(\text{deg}/s) \\\\ \dot{\theta} &= q \cos \phi - r \sin \phi, rad/s(\text{deg}/s) \\\\ \dot{\psi} &= \frac{q \sin \phi + r \cos \phi}{\cos \theta}, rad/s(\text{deg}/s) \end{aligned} \tag{16}$$

where

θ = 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.


Table 1. Aerodynamic coefficients as a function of angle of attack.

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:

$$\dot{u} = \frac{F\_{A\_{\rm b}} + F\_{P\_{\rm b}} + F\_{S\_{\rm b}}}{m} - (qv - rv), m/s^2$$

$$\dot{v} = \frac{F\_{A\_{\rm b}} + F\_{P\_{\rm b}} + F\_{S\_{\rm b}}}{m} - (ru - pw), m/s^2$$

$$\dot{w} = \frac{F\_{A\_{\rm a}} + F\_{P\_{\rm b}} + F\_{S\_{\rm a}}}{m} - (pv - qu), m/s^2$$

$$\dot{p} = \frac{L\_A + L\_p - qr(I\_z - I\_y)}{I\_x}, rad/s^2 \text{(deg/s}^2\text{)}$$

$$\dot{q} = \frac{M\_A + M\_p - rp(I\_x - I\_z)}{I\_y}, rad/s^2 \text{(deg/s}^2\text{)}\tag{17}$$

$$\dot{r} = \frac{N\_A + N\_p - pq(I\_y - I\_x)}{I\_z}, rad/s^2 \text{(deg/s}^2\text{)}$$

Figure 3. Six-degree-of-freedom motion of a rocket.

$$\begin{aligned} \dot{\phi} &= p + (q \sin \phi + r \cos \phi \quad \tan \theta, r \text{ad}/s(\text{deg}/s)) \\\\ \dot{\theta} &= q \cos \phi - r \sin \phi, r \text{ad}/s(\text{deg}/s) \\\\ \dot{\psi} &= \frac{q \sin \phi + r \cos \phi}{\cos \theta}, r \text{ad}/s(\text{deg}/s) \end{aligned}$$

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.

#### 2.4. Trimming

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 [11]. The set of conditions are the rocket's accelerations. The variables associated with the trim problems can be divided into three categories:


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

$$
\boldsymbol{\rho}\_{\rm la} = \begin{bmatrix} \dot{\boldsymbol{v}} & \dot{\boldsymbol{p}} & \dot{\boldsymbol{r}} \ \boldsymbol{\beta} \end{bmatrix}^{\rm T}. \tag{18}
$$

$$
\sigma\_{lo} = \begin{bmatrix} \dot{u} & \dot{w} & \dot{q} & a \end{bmatrix}^T. \tag{19}
$$

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:

$$\mathbf{c}\_{\text{lat}} = \begin{bmatrix} \delta\_{\text{a}} & \delta\_{\text{r}} & \phi & \psi \end{bmatrix}^{\text{T}},\tag{20}$$

$$\mathbf{c}\_{\rm lb} = \begin{bmatrix} \delta\_e & \boldsymbol{\pi} & \boldsymbol{\theta} \end{bmatrix}^T,\tag{21}$$

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):

$$
\sigma\_{i+1} - \sigma\_i = f\_i[\mathbf{c}\_{i+1} - \mathbf{c}\_i]. \tag{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:

$$\begin{aligned} \left.I\_{i} = \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}}\right|\_{c\_{(\bar{\mathbf{\boldsymbol{\theta}})}}} &= \frac{\begin{bmatrix} \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} \\ \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} \\ \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} \\ \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\theta}}}}{\partial \mathbf{\bar{\boldsymbol{\theta}}}} & \frac{\partial \mathbf{\bar{\boldsymbol{\$$

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.

#### 2.5. Linearization

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 firstorder differential equations and an algebraic output equation as [12]:

$$\begin{aligned} \dot{\mathbf{x}}\_1 &= f\_1(\mathbf{x}\_1, \dots, \mathbf{x}\_n, u) \\ \dot{\mathbf{x}}\_2 &= f\_2(\mathbf{x}\_1, \dots, \mathbf{x}\_n, u) \\ &\dots \\ \dot{\mathbf{x}}\_n &= f\_n(\mathbf{x}\_1, \dots, \mathbf{x}\_n, u) \\ y &= h(\mathbf{x}\_1, \dots, \mathbf{x}\_n, u) \end{aligned} \tag{25}$$

The column vector x = [x1,…xn] <sup>T</sup> is called the state of the system. The scalars u and y are called the control input and the system output, respectively, denoting

$$f(\mathbf{x}, \boldsymbol{\mu}) = \begin{bmatrix} f\_1(\mathbf{x}\_1, \dots, \mathbf{x}\_n, \boldsymbol{\mu}) \\ f\_2(\mathbf{x}\_1, \dots, \mathbf{x}\_n, \boldsymbol{\mu}) \\ \vdots \\ f\_n(\mathbf{x}\_1, \dots, \mathbf{x}\_n, \boldsymbol{\mu}) \end{bmatrix}. \tag{26}$$

Concisely, (26) is written as

$$\begin{aligned} \dot{\mathbf{x}} &= f(\mathbf{x}, \boldsymbol{\mu}), \\ \mathbf{y} &= h(\mathbf{x}, \boldsymbol{\mu}). \end{aligned} \tag{27}$$

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

$$\begin{aligned} \dot{\mathbf{x}} &= A\mathbf{x} + Bu \\ \mathbf{y} &= \mathbf{C}\mathbf{x} + Du. \end{aligned} \tag{28}$$

where A is nxn, B is nx1, C is 1xn, 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 <sup>x</sup><sup>∗</sup> = [x1…<sup>x</sup> ] T <sup>n</sup> obtained when <sup>∗</sup> <sup>∗</sup> <sup>u</sup> <sup>=</sup> <sup>u</sup> , noting that <sup>Δ</sup><sup>x</sup> <sup>¼</sup> <sup>x</sup> <sup>¼</sup> x , we define <sup>a</sup> coordinate transformation as follows:

$$
\Delta \mathbf{x} = \begin{bmatrix} \Delta \mathbf{x}\_1 \\ \vdots \\ \Delta \mathbf{x}\_n \end{bmatrix} = \begin{bmatrix} \mathbf{x}\_1 - \mathbf{x}\_1^\* \\ \vdots \\ \mathbf{x}\_n - \mathbf{x}\_n^\* \end{bmatrix}.
$$

<sup>∗</sup> <sup>∗</sup> Further, denoting <sup>Δ</sup><sup>u</sup> <sup>¼</sup> <sup>u</sup> � <sup>u</sup> , <sup>Δ</sup><sup>y</sup> <sup>¼</sup> <sup>y</sup> � h x<sup>ð</sup> <sup>∗</sup>; <sup>u</sup> <sup>Þ</sup>. The new coordinates <sup>Δ</sup>x, <sup>Δ</sup>u, and <sup>Δ</sup><sup>y</sup> 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 [12] state is given by

$$\begin{aligned} A\dot{\mathbf{x}} &= A\Delta \mathbf{x} + B\Delta u\\ \Delta y &= \mathbf{C}\Delta \mathbf{x} + D\Delta u \end{aligned} \tag{29}$$

where

� � � � � � � � � � 2 ∂f <sup>1</sup> <sup>∗</sup> <sup>∗</sup> <sup>∗</sup> ∂f <sup>1</sup> <sup>∗</sup> <sup>∗</sup> <sup>∗</sup> 3 � � 6 <sup>x</sup>1; …; xn; <sup>u</sup> <sup>⋯</sup> <sup>x</sup>1;…; <sup>x</sup> ; <sup>u</sup> <sup>∂</sup><sup>x</sup> <sup>n</sup> <sup>7</sup> <sup>∂</sup><sup>f</sup> <sup>6</sup> <sup>1</sup> <sup>∂</sup>xn <sup>7</sup> <sup>A</sup> ¼ ¼ <sup>6</sup> ⋯ ⋯ ⋯ <sup>7</sup>, <sup>∂</sup><sup>u</sup> <sup>x</sup><sup>∗</sup>,u<sup>∗</sup> <sup>4</sup> <sup>∂</sup><sup>f</sup> � � <sup>∂</sup><sup>f</sup> � � <sup>5</sup> <sup>n</sup> <sup>∗</sup> ∗ ∗ <sup>n</sup> <sup>∗</sup> ∗ ∗ <sup>x</sup>1; …; <sup>x</sup> ; <sup>u</sup> <sup>⋯</sup> <sup>x</sup>1;…; <sup>x</sup> ; <sup>u</sup> <sup>n</sup> <sup>n</sup> <sup>∂</sup>x<sup>1</sup> <sup>∂</sup>xn <sup>2</sup> <sup>3</sup> <sup>∂</sup><sup>f</sup> <sup>1</sup> � � <sup>∗</sup> ∗ ∗ <sup>x</sup>1;…; <sup>x</sup> ; <sup>u</sup> � � 6 <sup>n</sup> 7 ∂u B ¼ ¼ 6 6 ⋮ 7 7, <sup>∂</sup><sup>f</sup> <sup>6</sup> <sup>7</sup> <sup>∂</sup><sup>u</sup> <sup>x</sup><sup>∗</sup>,u<sup>∗</sup> <sup>4</sup> <sup>∂</sup><sup>f</sup> � � <sup>5</sup> <sup>n</sup> <sup>∗</sup> <sup>∗</sup> <sup>∗</sup> <sup>x</sup>1;…; <sup>x</sup> ; <sup>u</sup> <sup>n</sup> <sup>∂</sup><sup>u</sup> <sup>∂</sup><sup>h</sup> <sup>∂</sup><sup>h</sup> � � <sup>∂</sup><sup>h</sup> � � <sup>∂</sup><sup>h</sup> ∗ ∗ ∗ ∗ ∗ ∗ <sup>C</sup> <sup>¼</sup> <sup>¼</sup> <sup>x</sup>1; …; <sup>x</sup> ; <sup>u</sup> <sup>⋯</sup> <sup>x</sup>1;…; <sup>x</sup> ; <sup>u</sup> , D <sup>¼</sup> : <sup>n</sup> <sup>n</sup> <sup>∂</sup><sup>x</sup> <sup>x</sup><sup>∗</sup>,u<sup>∗</sup> <sup>∂</sup>x<sup>1</sup> <sup>∂</sup>xn <sup>∂</sup><sup>u</sup> ,u<sup>∗</sup> <sup>x</sup><sup>∗</sup>

� � Note that the matrices A, B, C, and D have constant coefficients in that all partial derivatives <sup>∗</sup> ∗ ∗ are evaluated at the numerical values x1;…; x ; u . <sup>n</sup>

#### 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)):

$$\begin{aligned} \dot{\upsilon} &= \frac{F\_{A\_{y\_1}} + F\_{P\_{y\_1}} + F\_{S\_{y\_1}}}{m} - (ru - pw), m/s^2 \\\\ \dot{p} &= \frac{L\_A + L\_p - qr\left(I\_z - I\_y\right)}{I\_x}, rad/s^2 \left(\text{deg}/s^2\right) \end{aligned}$$

$$\begin{aligned} \dot{r} &= \frac{N\_A + N\_p - pq\left(I\_y - I\_x\right)}{I\_z}, rad/s^2 \left(\text{deg}/s^2\right) \\\\ \dot{\phi} &= p + \left(q\sin\phi + r\cos\phi\right)\tan\theta, rad/s(\text{deg}/s) \\\\ \dot{\psi} &= \frac{q\sin\phi + r\cos\phi}{\cos\theta}, rad/s(\text{deg}/s) \end{aligned} \tag{30}$$

For a completely computed aerodynamic coefficients and their derivatives, (30) can be expressed in state-space form [13] as

$$
\begin{bmatrix}
\dot{\boldsymbol{v}} \\
\dot{p} \\
\dot{r} \\
\dot{\phi} \\
\dot{\psi}
\end{bmatrix} = \begin{bmatrix}
\mathcal{Y}\_{\boldsymbol{v}} & \mathcal{Y}\_{\boldsymbol{p}} & \mathcal{Y}\_{\boldsymbol{r}} & \mathcal{Y}\_{\boldsymbol{\phi}} & \mathcal{Y}\_{\boldsymbol{\psi}} \\
l\_{\boldsymbol{v}} & l\_{\boldsymbol{p}} & l\_{\boldsymbol{r}} & l\_{\boldsymbol{\phi}} & l\_{\boldsymbol{\psi}} \\
n\_{\boldsymbol{v}} & n\_{\boldsymbol{p}} & n\_{\boldsymbol{r}} & n\_{\boldsymbol{\phi}} & n\_{\boldsymbol{\psi}} \\
0 & \mathbf{1} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
0 & \mathbf{0} & \mathbf{1} & \mathbf{0} & \mathbf{0}
\end{bmatrix} \begin{bmatrix}
\boldsymbol{v} \\
\boldsymbol{p} \\
\boldsymbol{r} \\
\boldsymbol{\phi} \\
\boldsymbol{\psi}
\end{bmatrix} + \begin{bmatrix}
\mathcal{Y}\_{\boldsymbol{\xi}} & \mathcal{Y}\_{\boldsymbol{\zeta}} \\
l\_{\boldsymbol{\xi}} & l\_{\boldsymbol{\zeta}} \\
n\_{\boldsymbol{\xi}} & n\_{\boldsymbol{\zeta}} \\
\mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0}
\end{bmatrix} \boldsymbol{\zeta}.\tag{31}
$$

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 (M<sup>x</sup> and Mz), and two of the Euler angles from the 6DoF equations as given in (17):

$$\begin{aligned} \dot{u} &= \frac{F\_{A\_{\eta\_b}} + F\_{P\_{z\_b}} + F\_{\mathcal{S}\_{z\_b}}}{m} - (qv - rv), m/s^2\\ \dot{w} &= \frac{F\_{A\_{\eta\_b}} + F\_{P\_{z\_b}} + F\_{\mathcal{S}\_{z\_b}}}{m} - (pv - qu), m/s^2\\ \dot{q} &= \frac{M\_A + M\_p - rp(I\_x - I\_z)}{I\_y}, rad/s^2(\text{deg}/s^2) \end{aligned} \tag{32}$$
 
$$\dot{\theta} = q \cos \phi - r \sin \phi, rad/s(\text{deg}/s)$$

Just as with (31), (32) can also be re-expressed in state space as

$$
\begin{bmatrix}
\dot{u} \\
\dot{w} \\
\dot{q} \\
\dot{\theta}
\end{bmatrix} = \begin{bmatrix}
\mathbf{x}\_{\mu} & \mathbf{x}\_{w} & \mathbf{x}\_{q} & \mathbf{x}\_{\theta} \\
\mathbf{z}\_{\mu} & \mathbf{z}\_{w} & \mathbf{z}\_{q} & \mathbf{z}\_{\theta} \\
m\_{u} & m\_{w} & m\_{q} & m\_{\theta} \\
0 & 0 & 1 & 0
\end{bmatrix} \begin{bmatrix}
u \\ w \\ q \\ \theta
\end{bmatrix} + \begin{bmatrix}
\mathbf{x}\_{\eta} & \mathbf{x}\_{\tau} \\
\mathbf{z}\_{\eta} & \mathbf{z}\_{\tau} \\
m\_{\eta} & m\_{\tau} \\
0 & 0
\end{bmatrix} \begin{bmatrix}
\eta \\ \tau
\end{bmatrix} \tag{33}
$$

In MATLAB® the linmod [14] 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:


It can be shown that a typical trimmed and linearized model of the pitch plane (longitudinal dynamics) for a rocket [15] 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):

$$
\begin{bmatrix} \dot{\theta} \\ \dot{q} \\ \dot{w} \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 \\ 14.7805 & 0 & 0.01958 \\ -100.858 & 1 & -0.1256 \end{bmatrix} \begin{bmatrix} \theta \\ q \\ w \end{bmatrix} + \begin{bmatrix} 0 \\ 3.4858 \\ 20.42 \end{bmatrix} \delta\_\eta + \begin{bmatrix} 0 \\ 14.7805 \\ -94.8557 \end{bmatrix} a\_w$$
 
$$y = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} \theta \\ q \\ w \end{bmatrix}$$

where

$$A = \begin{bmatrix} 0 & 1 & 0 \\ 14.7805 & 0 & 0.01958 \\ -100.858 & 1 & -0.1256 \end{bmatrix}, B = \begin{bmatrix} 0 \\ 3.4858 \\ 20.42 \end{bmatrix}, G = \begin{bmatrix} 0 \\ 14.7805 \\ -94.8557 \end{bmatrix} \\ \text{C} = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}, D = 0.1$$

#### 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.

### 4. Conclusion

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.

## Author details

Aliyu Bhar Kisabo\* and Aliyu Funmilayo Adebimpe

\*Address all correspondence to: aliyukisabo@yahoo.com

Centre for Space Transport and Propulsion (CSTP), Epe, Lagos-State, Nigeria

## References


**Chapter 3**

Provisional chapter

**Discrete Element Modeling of a Projectile Impacting**

DOI: 10.5772/intechopen.75550

From theoretical standpoint, it is difficult to analytically build a general theory and physical principles that critically describe the mechanical behaviour of granular systems. There are many substantial gaps in understanding the mechanical principles that govern these particulate systems. In this chapter, based on a two-dimensional soft particle discrete element method (DEM), a numerical approach is developed to investigate the vertical penetration of a non-rotating and rotating projectile into a granular system. The model outcomes reveal that there is a linear proportion between the projectile's impact velocity and its penetration downward displacement. Moreover, depending on the rotation direction, there is a significant deviation of the x-coordinate of the final stopping point of a rotating projectile from that of its original impact point. For negative angular velocities, a deviation to the right occurs while a left deviation has been recorded for positive angular

Keywords: discrete element method, granular bed, rotating projectile, penetration,

A granular material is any collection of many macroscopic discrete solids, whose typical size ranges from micrometres to centimetres such as sand, coal, sugar, and rice. It is obvious that these materials cannot be characterised as gas, fluid, or solid. But, they have hybrid state between these three phases depending on the average energy of the individual grains inside the granular system. Piles, snow avalanches, and planetary rings (interstellar dust) are, respec-

> © 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and eproduction in any medium, provided the original work is properly cited.

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use,

distribution, and reproduction in any medium, provided the original work is properly cited.

tively, considered as solid, fluid, and gas phases for granular materials [1].

Discrete Element Modeling of a Projectile Impacting

**and Penetrating into Granular Systems**

and Penetrating into Granular Systems

Additional information is available at the end of the chapter

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/intechopen.75550

Waseem Ghazi Alshanti

Waseem Ghazi Alshanti

Abstract

velocities.

1. Introduction

vertical impact
