Modelling and Simulation of Surface to Surface Missile General Platform

: In this paper, a general platform of a surface ‐ to ‐ surface missile model is presented. An accurate missile system model was made with consideration of the rigid body dynamics in space with six degrees of freedom and equation of motion starting from Newton second law to the non ‐ linear state space model of the missile. Computer simulations have been executed using the MATLAB / Simulink software. Verification of the system model has been made.


Introduction
Modelling of flight dynamics is limited to some constraints, such as picking a frame reference coordinate, and derivation of dominant equations of motion.In this derivation, the vehicle is treated as a rigid body when considering its rotational motion, which reduces the degrees of freedom from infinite to only six.The equations of motion can be described in two groups; kinematic equations and dynamic / kinetic equations [1].
The aerodynamic and propulsion force and moment vectors are obtained from Newton's laws of motion and gravity, as well as the aerothermodynamics principles [2].A MathWork's Simulink software toolbox; Aerospace Block-Set [1] has been used to construct the flight model including: 6 degrees-of-freedom equations of motion, linearization of aerodynamics, also a turbofan engine block, 2 nd and 3 rd order actuators, several standard atmosphere and earth gravity models, and finally the wind disturbance models [1].

Earth-Centred Inertial (ECI) Reference Frame <OXcYcZc>
This is a coordinate frame with its origin at the Earth's centre, which moves with the Earth, but with a fixed orientation relative to the stars [3].

Body-fixed Reference Frame <CXBYBZB>
This is a coordinate frame which is located at the vehicle centre of gravity with origin C as shown in Fig. 1.The frame moves with the flying vehicle and the X-axis is directed towards the longitudinal axis of the vehicle, the Y-axis is directed towards the lateral axis rightward from the rear and the Z-axis is directed downwards to complete the orthogonal tetrahedral frame Fig. 2, [3].

Wind-axes Reference Frame <CXwYwZw>
This is a coordinate frame that is located at the vehicle centre of gravity.The frame moves with the flying vehicle and the Xw-axis is directed against the free stream velocity vector.This frame results from rotating the body-fixed reference frame with an angle of attack α around the negative YB-axis and then from the rotation around the resulted Zaxis with side slip angle β, as shown in Fig. 2.

Firing-point Reference Frame <OXfYfZf>
The axes are fixed on the Earth and rotate with it while the origin is placed on the firing point.The OXf axis is taken in the direction of firing direction and the axis OZf tends to be directed downward vertically.Namely, it takes the direction of the gravity, and OYf is taken so that the coordinates establish the right-handed system [3], as shown in Fig. 3.

Fig. 3 Firing-point reference frame
It is necessary to choose an inertial reference frame (which is considered fixed) in the development of the vehicle mathematical model.Inertial frame is considered to be non-rotating and non-accelerating (but it may be moved with constant rectilinear velocity) relative to the average position of the "fixed" distant stars.If the Earth is assumed to be flat and non-rotating, the firing point coordinate system could be considered as the inertial coordinate system [3].

Transformation between Reference Frames
According to Euler, the orientation of one reference frame with respect to another is completely defined by a set of three successive virtual perpendicular rotations in a specific order.The orientation is then completely defined by the values of the angles of rations called Euler angles and the order of them is widely agreed in aeronautics to be yawing, pitching then rolling [3].

Firing-point and Body-fixed Frames
This transformation is used to get the components of any vector defined in the firingpoint reference frame into the body-fixed reference frame axes.The transformation matrix results from a sequence of rotations (ψ, θ, ϕ) used with a specific order to describe the instantaneous attitude according to the firing reference frame [3].FB Body f f cos cos cos sin sin cos sin sin sin cos cos cos sin sin sin sin cos sin sin cos sin cos sin cos cos sin sin cos cos

Wind-axes and Body-fixed Frames
This transformation is used to get the components of any vector defined in the windaxes reference frame (e.g.aerodynamic forces) into the body-fixed reference frame axes.The transformation matrix results from a sequence of basis rotations (β, α) used to describe the instantaneous orientation of velocity vector starting from the wind frame [3].
FB Body f f cos cos cos sin sin cos sin sin sin cos cos cos sin sin sin sin cos sin sin cos sin cos sin cos cos sin sin cos cos

Equations of Motion
The vehicle equations of motion can be calculated from Newton's second law of motion.
The equations show that the summation of all forces applied externally on a body should be equal to the time rate of change of its momentum.The time rates of change take into consideration the inertial reference frame, under certain hypotheses: • The intended vehicle is a surface-surface tactical missile utilizing solid propellant engine for propulsion in the form of boost-sustainer motor and the missile is aerodynamically controlled using four all-movable fins in the form of "X"-configuration, as shown in Fig. 2. • The origin of the body axes system is located at the centre of gravity of the vehicle and fixed to the vehicle and it rotates with it.The axes are taken, as shown in Fig. 2, with XB forward from vehicle centre of gravity towards the vehicle nose tip, YB out to the right and ZB downward.• The axes XB and ZB lie in the plane of symmetry of the vehicle, and axes XB and YB lie in the other plane of symmetry.Hence, the cross moments of inertia Ixy, Ixz and Iyz have all zero value.The missile has the same mass distribution in the axes YB and ZB , so, the second moments of inertia of both axes are equal (Iyy = Izz).

• The equations of motion are composed of:
-Dynamic equations: they consist of vector force and moment equations.
-Kinematic equations: they consist of attitude and trajectory equations.• Although the vehicle mass changes considerably due to the propellant burning in the powered phase and it remains constant after engine shut-off, the amount of propellant consumed during certain integration time step in the powered region of the vehicle flight may be safely neglected [1].

Force Equations
The force equation is derived with respect to the body-fixed reference frame to get the acceleration with respect to body frame as shown in the next equations: From Eq. 3, the force equation can be written in the form:

Moment Equations
To get the equations of angular motion, it is necessary to consider that the moment affecting the vehicle is equal to the changing rate of its angular momentum HB, then: For fixed controls airframe, the overall moment MB acting on airframe with respect to body axes resulted from the moment of the resultant aerodynamic force due to the offset of the centre of pressure (the point of action of resultant aerodynamic force) from the centre of mass, in addition to the moment resulting from the thrust force due to the offset between the point of action of the thrust and the centre of mass.The moment equation [3,4] could be written as:

Attitude Equations
The R matrix describes the attitude of the vehicle as it changes with time.Its elements are functions of the three Euler angles and the body-axes rates (p, q, r) produce a corresponding set of attitude angular rates, depending on the particular instantaneous attitude.
The relation between Euler angles derivative and angular velocity of airframe [3, 5] could be written as follows: 1 tan sin tan cos 0 cos sin sin cos 0 cos cos Eq. 7 represents the kinematic relation between changing rate of Euler angles of the vehicle and body angular velocity.Eq. 7 has been commonly used for aircraft simulation, but it is not suitable for around-the-Earth flight, all-attitude flight, and simulation of spinning bodies.

Navigation Equations
The location of the vehicle mass centre relative to the Earth is given by the spherical polar coordinates R (geocentric radius), µ (longitude) and λ (latitude), as shown in Eq. 2. But, in case of calculation of position of vehicle mass centre relative to an inertial frame (firing-plane frame), the position with respect to the inertial frame is the integration of the velocity in that frame, i.e.:

Wind Equations
The aerodynamic forces on vehicle are created by its motion relative to the surrounding air.In case of existence of wind (atmosphere is not at rest), the air itself is considered in motion relative to the inertial reference frame.It is assumed that the local wind has north, east, and down velocity components Vwx, Vwy and Vwz, respectively, in the firing plane.Then the absolute velocity of the missile with respect to the firing-plane is the vector sum of the velocity of missile with respect to the air VR and the velocity of wind with respect to the firing-plane wf transformed to body-frame [2].The velocity of the vehicle centre-of-mass with respect to the air is then given by:

Missile Simulation
The core of the simulation is a non-linear model of the vehicle dynamics which includes twelve Ordinary non-linear Differential Equations (ODE's or state equations) which are Eqs 4, 6, 7, 8, and a large number of equations of output.This model can be broken down into a number of different modules which are almost independent from the kind of moving vehicle under consideration.

Vehicle characteristics and flight conditions
For calculating the dynamic trajectory [6], the airframe of the vehicle is specified with a hypothetical frame of cruciform shape, conical head and four all-movable tail fins of X-configuration, as shown in Fig. 4.

Thrust
A general form of the thrust equation for the rocket motor of solid propellant [7,8] is given by: ( )

T mV P P A
From Eq. 10, the thrust T0 at sea level can be written as: ( )

T mV P P A
where P0 is the atmospheric pressure at sea level.
The thrust TH at certain altitude H is given by: This equation yields:

T mV P P A P P A
Finally: ( ) Then, from the resulting data of rocket motor static test at sea level T0 and the pressure at certain altitude H, the thrust at that altitude TH can be calculated.

Atmosphere
According to ISA model, the air temperature T changes with altitude [1,9] Where T [K] is the temperature at altitude H [m] and T0 is the temperature at sea level and it equals 288.15 K.
Then the atmosphere layers will be divided into two types: a) Gradient layer in which the temperature changes with factor τd b) Isothermal layer in which the temperature remains constant (τd =0) 0.0065 [°/m] for from 0 to 11000 m 0 for from 11000 m to 20 000 m 0.001[°/m] for from 20 000 m to 32 000 m 0.0028 [°/m] for from 32 000 m to 47 000 m The equation of atmospheric pressure can be derived to be [1,9]: • For isothermal layer The density is given by: • Isothermal layer The sonic speed at altitude H is given by: 20.05 where χ is the exponent ratio and it is equal to 1.4 for air [4] and go is the gravitational acceleration at sea level

Aerodynamic Forces and Moments
The aerodynamic forces are defined in terms of dimensionless coefficients, the flight dynamic pressure and the reference area [5] as follows: The resultant aerodynamic force acts at the point called the centre of pressure and this point is usually not located at the centre of gravity of the vehicle.The aerodynamic moment definition is the same as the aerodynamic forces in terms of dimensionless coefficients, dynamic pressure, reference area and a characteristic length as follows: In case of zero aerodynamic roll, the angle of attack and sideslip angle can be defined in terms of the velocity components [4]

Aerodynamic Coefficients
For simplicity, the aerodynamic forces and moments acting on the aircraft are defined in terms of dimensionless aerodynamic coefficients Cx, Cy, Cz, Cl, Cm, and Cn and they are called total coefficients of moments and forces.They depend on the control surfaces deflection, the aerodynamic angles and the angular rates [3,10].

Construction of Simulation
The open-loop simulation is performed using Matlab / Simulink package.The block diagram is shown in Fig. 5.

Fig. 5 Block diagram of the vehicle
From Fig. 5, it can be concluded that: • The control actions are not taken into account.
• Neither the disturbance of the atmosphere is taken into account.
• The vehicle does not contain any sensors and is not guided.

OPENLOOP SIMULINK / MATLAB Model
This part discusses the Simulink / Matlab modelling of Missile system.The model is based on the derivation of the mathematical model that was listed before in the mathematical modelling part.The complete set of equations of motion is coded into Simulink blocks leading to a missile simulator that helps in the control process [11].

Missile Simulator
The missile simulator is the main system block and it includes various subsystems.It has 12 output signals representing the position and the orientation of the vehicle and their rates or velocities considering the body reference frame.Each internal subsystem of the missile simulator is illustrated.

Dynamics Block
The dynamic block is responsible for producing 12-output control signal that controls the equations of motion.It produces these signals according to the input of the mathematical modelling of the missile and its movements.The input to the block includes the forces on the body, body momentum, mass, mass of inertia, gravity and time respectively.

5.7.3
Aerodynamics Block Aerodynamics block is responsible for calculating the body forces and momentums for the dynamics block as an output to control the missile body according to the input signals; Velocity, Atmosphere density ρ, Acceleration, Height and Centre of gravity respectively.

5.7.4
Atmosphere Block The atmosphere block is responsible for calculating the ρ (atmosphere density), acceleration, gravity, height and the pressure at certain level as an output signal according to the position of the missile as an input.

5.7.5
Propulsion Block The propulsion block is responsible for calculating the thrust force needed for the dynamic block according to the calculated pressure.The missile parameters like the centre of gravity and mass of inertia are fed as inputs to the block.

Results of Simulation
Fig. 6 shows the results of simulation performed on Matlab / Simulink package.

Flight Path with Different Launching Angles
Changing of the launching angle will directly affect the range and the summit of the missile trajectory due to the change of the initial pitch angle θ0 which is one of the initial conditions [12].The previous results of simulation are calculated at launching angle 60˚.
In order to study the launching angle effect on various parameters, simulation is performed at different launching angles and then the trajectories of uncontrolled missile are drawn at the angles 40˚, 45˚, 50˚, 55˚, 60˚, 65˚, 70˚, 80˚, as shown in Fig. 7.The results clarify that the range increases with the launching angle till 65˚ and then decreases with angle increasing.It is concluded that the maximum range is obtained at launching angle θ0 = 60˚ -65˚, but increasing of launching angle increases the summit and it is preferred to keep the missile in the dense atmosphere in order to attain a quiet aerodynamic control force created by fins.Thus, it was suitable to choose the launching angle to be 60˚.Based on the results, it can be concluded that: • The total flight time of missile is 180.2 s.
• The max. altitude reached by the missile is 35.5 km.
• The total distance travelled by the missile on a flat Earth along inertial frame Xaxis is 115.5 km.• The powered phase is considered to be terminated at the beginning of a velocity decrease due to motor burnout.The beginning of the decrease occurs at the time (t = 13.05 s).So, the powered phase will be from the time (t = 0 to 13.05 s) and the unpowered phase will be from the time (t = 13.05 s till the end of flight).• From Fig. 6a, which describes centre of gravity and centre of pressure change with time, the missile will remain statically stable during the whole flight time since the centre pressure remains behind centre of gravity.• In Fig. 6, the response of the angle of attack with respect to the time mainly depends on the undamped natural frequency, damping ratio of short period mode (ωsp, ζsp), and slightly on undamped natural frequency and damping ratio of phugoid mode (ωph, ζph) [13].
• The undamped natural frequency of short period mode ωsp depends mainly on pitching moment because of the angle of attack derivative Mα.The damping ratio ζsp depends on pitching moment due to the pitch rate derivative Mq and normal force due to the angle of attack derivative Zα [3].• The frequency of the angle of attack response increases rapidly in the powered phase due to the rapid increase of pitching moment because of angle of attack derivative Mα in the same phase.Then the frequency decreases at the beginning of the unpowered phase due to the decrease of the same derivative.Then it returns to the increase as the dynamic pressure increases and then the derivative increases as well.• The damping ratio of the angle of attack response increases rapidly in the powered phase due to the rapid increase of both pitching moments because of the pitch rate derivative Mq and normal force because of angle of attack derivative Zα in the same phase.Then, the damping ratio decreases at the beginning of the unpowered phase due to decrease of the both same derivative.Then, it returns to the increase as the dynamic pressure increases and then the derivatives increase, too.

Model Verification
The model has been verified with real recorded telemetry data with the Russian missile type "9k52 Luna".Additionally, it has been verified with trajectory / projectiles equations, to hit a target at altitude y and range x when fired from point (0,0) considering the initial speed v so the required angle(s) of launch are: ( ) The square roots of Eq. 24 refer to the two potential launching angles, provided that they are not imaginary.This formula is used to find the desired launching angle for model verification such that the desired altitude reached by the missile is 35.5 km and the distance travelled by the missile along the flat earth is 115.5 km.By substitution we will find out that the launching angle is 60˚.

Conclusion
In this paper, a general platform of a surface-to-surface missile model is presented, focusing on discussing the system dynamics and Simulink model.An accurate missile system model was made with consideration of the rigid body dynamics in space with six degrees of freedom and equation of motion starting from Newton second law to the non-linear state space model of the missile.Computer simulations have been executed using the MATLAB / Simulink software.Verification of the system model has been made.The system response shows a good result after testing under various input signals.

Fig. 4
Fig. 4 Model configuration and main dimensions

Fig. 6a Fig
Fig. 6a Results of open-loop simulation -Angle of attack response and static stability