Laboratory Stand for Simplified Verification of Missile Dynamics for Didactic Purposes

The purpose of this study was to describe dynamics of missile based on bibliography, simplify the model to the easiest form, and verify it with the use of a designed laboratory stand. Next, a simple control system of the missile was designed. The stand is prepared for didactic purposes. A mathematical model was derived by applying Newton’s second law together with the earth’s coordinate system and the base coordinate system. Parameters of the actual rocket model were determined using the created laboratory stand. Synthesis of the rocket roll system was designed using Simulink PID Tuner application (MathWorks, Inc. Natick, Massachusetts, USA). A control system was based on a proportional-integral (PI) regulator. The designed control system was subjected to simulated tests in MATLAB Simulink (MathWorks, Inc. Natick, Massachusetts, USA).


Nomenclature
U, V, W, -linear velocities u, v, w -increments of velocities ϕ, Θ, Ψ -large angular displacements φ, θ, ψ -small angular displacements, P, Q, R -angular velocities p, q, r -increments of angular velocities p x , p y , p 2 -radius of turn around x, y, z axis ρ -total radius of turn H x , H y , H z -moments of momentum I x , I y , I 2 -moments of inertia J xz , J yz , J zx -area moments of inertia, F x , F y , F z -sums of forces acting on an object, F 0 -sums of forces in equilibrium, Fg x , Fg y , Fg z -forces of gravity a x , a y , a z -linear accelerations, ΔF -sum of external forces, out of equilibrium M, N, L -sums of moments acting on an object, M 0 -sum of moments inequilibrium ΔM -sum of external moments, out of equilibrium m -mass, v T -linear velocity, g -gravity acceleration, ω -rotational velocity v tan -velocity (tangential), r -rotational vector, U 0 , V 0 , W 0 -velocity vectors in equilibrium P 0 , Q 0 , R 0 -rotational velocities in equilibrium β -slip angle S 0 -reference surface S w -surface providing lift, S r -cross -section area of rocket b -span q dyn -dynamic pressure d -diameter R = ½ d -radius

Introduction
The purpose of this article is to show results of an engineer's thesis. In our project, algorithms are supposed to provide automatic control of specified roll speed of a given rocket. At the beginning of this study, we focused on bibliographical analysis of missile dynamics. The next step was to verify the applied mathematical model to the actual model. For this purpose, a laboratory stand was built, in which a flight of the rocket was simulated by providing a constant air flow, and rotation speed of the rocket and dynamic pressure of the air stream were examined. Obtained data were compared against a mathematical model. Finally, a control system was developed for the mathematical model of the rocket, and it was subjected to simulations tests. The first section focuses on obtaining equations for rocket's general motion. The second section uses equations of general motion to obtain equations in lateral motion and in the end transfer function of a rocket in roll mode. The third section depicts deriving the dimensionless coefficients used in transfer function of a rocket.

Rocket dynamics -rocket's roll mathematical model
The earth's and base coordinate systems were used (Figure 1) with the following assumptions [2]: 1. Axes OX and OZ lie in the plane of symmetry of a rocket, and J xy and J yz are equal to zero. 2. Rocket's mass stays constant through dynamic analysis. 3. Rocket has a rigid body, so that any two points remain fixed to each other. 4. The earth is an inertial reference and the atmosphere is fixed with respect to the earth.
Knowing kinematic relation between rocket body frame F B and F E frame, as well as gravity force coefficients and aerodynamic loads, general equations of motion can be derived, and these are shown in [2]: ( ) Along with equations of linear motion, the following equations are obtained: For synthesis of control algorithms, Eqs (1)-(6) are linearized. In a linearization process, it is assumed that rocket motion is decoupled to lateral and longitudinal motion. In this article, lateral motion is described. As in deriving equations for general motion of a rocket, the coordinate system associated with the object is used. If state disturbance is created during equilibrium, the X axis is displaced relative to the Z axis from the rocket velocity vector, and the slip angle β is formed. The positive direction of this angle is shown in Figure 2. Equations for lateral motion can be obtained from Eqs (1), (3), and (5). It is assumed that Q is negligibly small. Assuming a flight with a slip angle of zero (the flight speed is tangent to the flight path), then 0 , and U u =   At the beginning, the rocket is in an unaccelerated flight, so P 0 and R 0 are zero. Thus, P = p and R = r. With these assumptions, assuming small disturbance theory, we obtain the following equations: and presenting forces and moments in functions of , , , β ψ φ φ  , we can formulate the following equation: Plane perpendicular to OX Figure 1. Sketch of fixed and rocket axes.
If partial derivatives are considered linear at a given moment of the flight, then the differentials can be replaced by incremental changes. Thus, Eq. (10) is transformed into Eq. (11): From analysis of gravitational forces [2], we can formulate the following equations: Since in Eq. (11), the initial values of β, ψ, and ϕ are equal to zero, the differences D can be dropped, and then substituting Eq. (11) in Eq. (7) gives us Eq. (14): As U 0 ≌ U, we drop the index 0. By dividing Eq. (14) by S 0 and q dyn , substituting Eqs (12) and (13), and introducing dimensionless aerodynamic derivatives, we get the following equation: Eqs (1)-(3) can be developed and written as Eqs (16) Dividing by S 0 , q dyn , and b and using the record in the form of dimensionless coefficients, Eqs (16) and (17) take the following form: Eqs (18)-(20) are linearized equations describing the dynamics of the rocket in lateral motion. In analogy to the description of aircraft dynamics, approximations of individual types of motion can be introduced. Since the considerations in this paper concern the rotational movement of the rocket, this approximation will be applied. We assume that the rocket moves only with a rotational movement around its own axis. Then, the only equation necessary to describe the dynamics is the torque equation 2-30, in which it was assumed that ψ, β, and δ r take the zero value: Then, we can determine the transfer function of the rotational speed of the rocket in relation to the angle of deflection of the ailerons, using the following equation: The time constant is given by the following equation: The reference area, S 0 , is usually taken as the cross-sectional area of the missile. Cl p results from the change in the angle of attack of the aerodynamic surfaces caused by the rotational speed. Dropping the aileron causes rise in the angle of attack, and rising causes the opposite. These changes in angles of attack result in different lift and drag forces acting on aerodynamic surfaces. After the vectors of lift and drag forces along the X and Z axis of equilibrium are divided, it can be seen that the Z components generate a torque counteracting the rotational speed. Typically, the drag component is ignored and only the lift force is taken into account. This topic is more broadly defined in chapter 9 of the study by Perkins and Hage [3] and in NACA TR 635 of the study by Pearson [4].
To determine the value of Cl δ a , we need to know the appropriate parameters of the rocket, i.e., the ratio of the chord of the wing at the root to the chord of the wing at the tip λ, aileron effectiveness τ, and aspect ratio AR. Blakelock [2] gives a graphical relation between C l p , Cl δ a , and mentioned parameters, and these parameters were used in the article. Figure 3 shows the dimensions of the rocket completed as part of the work. It is characterized by the following parameters: Based on Figure 3, knowing that the length of the aileron is equal to the length of the wing: % b 2 =1, we read the value Clδ a τ ≈ 8 [2], with τ = 1, and so Cl δ a = 0.8. Then, knowing that λ ≈ 0.5 and AR ≌ 7, we get Cl p = -0.52 [2].

Laboratory Stand
The laboratory stand was built in order to verify mathematical model with respect to the real rocket model. Figure 4 shows the layout of the proposed laboratory stand. The factor responsible for angular movement is lift force on canards generated by dynamic pressure due to longitudinal movement. In the laboratory stand, a centrifugal fan was used, which provided air flow. To power the fan, a modified ATX power supply was used with a voltage supply of 3.3 V, 5 V, and 12 V and ground. A simplified rocket model was built, as shown in Figure 5.
The rocket model was placed on a brass rod allowing it to rotate. Figure 5 shows a modified manometer, which receives total pressure from the front; and the back of the manometer receives static pressure. As a result, dynamic pressure is derived by difference in fluid height.
To measure angular speed, opto-isolator KTIR0711S was used, wherein data from the sensor were processed by a mi-   crocontroller ATmega8 and sent to the USB port via a USB-UART converter. Measurement was carried out in the following way: to specified pin of the microcontroller, signal from opto-isolator was connected, which functioned as a sensor. In the back of our rocket, a marker was placed, at which our opto-isolator reacts ( Figure 6). The microcontroller, after receiving a voltage supply, interrupts the looping code and stores actual timer values. Afterward, the controller waits for the next event, at which a second timer value is stored. After ascertaining the values in-dicated by the two timers, the number of markers, and the resonator's frequency, we can calculate the rocket's roll speed.
Received data were displayed and saved via RealTerm terminal program.

Tests and Results
Canards were set contrary at an angle of 30° in relation to the longitudinal axis, and then, the fan was connected to a voltage supply of 5 V. In this way, the value of dynamic pressure obtained was 77.  Figure 7, we can see given values of rocket motion from derived data, compared to speed of a rocket model. The further step was to calculate friction coefficient between the rocket and laboratory stand. Mass on a string was attached to the rocket and then pulled upward at a given height. Time of the mass falling on the ground was measured. Given that rocket made exactly two turns, it was compared to time how this mass would fall on the ground without friction. Assuming constant angular acceleration value, in both conditions, we have compared value of acceleration without friction and with friction and roll data was multiplied, which can be seen on the graph on violet line. We can see satisfying correlation between a mathematical model and data from our laboratory stand.

Synthesis of Rocket Roll Control System
The final step was to design rocket roll control for the mathematical model of the rocket. We will apply the system shown in graph 13. Transfer function between rocket's roll speed and  canard's angle of inclination was derived at the beginning of this paper and is rewritten here as follows: Assuming servo dynamics, we can focus on choosing a controller for the system shown in Figure 8.
For the synthesis of the control system, Simulink's "PID Tuner" application was used. A proportional integral controller was used. Parameters of the controller were adjusted so that the following statements hold true: -Settling time was <0.1[s], -There was no overshoot/undershoot, -System was stable.
After adjusting the proportional gain (P ≌ 0.09) and integrating gain (I ≌ 1.23) of the controller, the system was commanded to maintain an angular speed of 4.9 [ rot sec ]. Response of the system is shown in Figure 9.

Summary and Conclusions
The objectives of the study were to determine the mathematical model of the rocket and to synthesize the algorithms of angular motion of the rocket. These goals have been fulfilled. The work was started by deriving the full mathematical model of the rocket and then determining the equations of dynamics in angular motion. Based on the derived equations, the dynamics of the rocket in angular motion were determined as the first-order inertial.
Another task was to verify the mathematical model's compatibility with the real object. The actual model of the rocket was designed and built, and its dynamic characteristics were determined on the basis of previously derived equations. In order to verify the relation between rocket's roll speed and canard's angle of inclination, a laboratory stand was created to measure angular speed. The stand can measure the dynamic pressure of the air flow by means of a differential pressure gauge and angular speed using a pulse counter and an optical isolator. The speed of the air flow is relatively low (11,26 [ m s ] -small Reynold's value), which is the laboratory stand maximum, when basic amateur missiles achieve speeds of 100 [ m s ]. This is a serious limitation, and the next step is the use of the presented model in a wind tunnel. By comparing the speed data of the mathematical model with the angular data from the measurements, considering the friction error, it was found that the assumed mathematical model of the rocket satisfies the dynamic properties of the rocket.  A rocket rotation control system has been designed. The regulator parameters have been selected using the Simulink PID Tuner application. Specific values of overshoot and settling time have been assumed, and on this basis, proportional and integrating gains have been chosen. The results made available in this article constitute preliminary work and can be developed in several directions. A further objective is to implement a designed control system in a missile using different autopilot systems and put it to tests in a real environment [8,9] and different missile configurations [10]. Another task is implementation of different control algorithms. The presented stand can also be used for hardware in loop simulation [8].