UFDC Home  myUFDC Home  Help 



Full Text  
CONTROL STRATEGIES FOR SUPERCAVITATING VEHICLES By ANUKUL GOEL A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2002 ACKNOWLEDGMENT S I would like to express my sincere gratitude to my committee chairman, Dr. Andrew Kurdila, for his invaluable guidance throughout the course of this project. I would also like to thank him for giving me this opportunity to work on such a fascinating project. I would also like to thank my committee cochair, Dr. Richard C. Lind, for his invaluable guidance and inspiration throughout the project. I would also like to show my sincere appreciation to Dr. John Dzielski, Dr. Nor man FitzCoy and Jammulamadaka Anand Kapardi for their valuable contributions to this project. I would also like to express my gratitude to all the members, past and present, of the Supercavitation Project. I would also like to thank the Office of Naval Research for the support of the research grant for the project. On a personal note, I would like to thank all my friends and family members whose support helped me to aim towards my goals. TABLE OF CONTENTS page ACKNOWLEDGMENTS . LIST OF TABLES . . . .. ii . . . . v i LIST OF FIGURES ............ ABSTRACT . . . . CHAPTER 1 INTRODUCTION . . . . . 1 1.1 Cavitation .... .. . . ... 1.2 Types of Supercavitating Projectiles . 1.3 Related Research ............. 1.4 Overview of this Thesis . ..... 2 CONFIGURATION OF VEHICLE . 2.1 Cavitator ........ 2 .2 F in s . . . . . 2.3 Maneuvering .............. 3 NONLINEAR EQUATIONS OF MOTION ..... 3.1 Kinematic Equations of Motion ................ 3.1.1 Orientation of the Torpedo ............... 3.1.2 Orientation of the Cavitator ............. 3.1.3 Orientation of Fins ................... 3.1.4 Angle of Attack and Sideslip ....... 3.1.5 Kinematic Equations .................. 3.2 Dynamic Equations of Motion ................ 3.2.1 Forces on Cavitator .. ................ 3.2.2 Forces on Fins . . . . . . 3.2.3 Gravitational Forces .. ............... 3.2.4 Equations of M otion .. ............... 4 LINEARIZATION ......... 4.1 Linearization ....... 2 4 5 7 . . 1 . . 2 . . 34 4.1.1 Need for Linearization ................. 34 4.1.2 Generic Form of Equations of Motion . . .. . 35 4.1.3 Small Disturbance Theory . . . . . 35 4.1.4 Stability and Control Derivatives . . .. . 37 4.2 State Space Representation . . . . . . 42 5 CONTROL DESIGN SETUP . . . . . 46 5.1 Openloop Performance . . . . . . . 47 5.2 ClosedLoop Problem . . . . . . . 50 5.3 Robustness of the Controller . . . . . . 52 5.3.1 G ain M argin . . . . . . . 52 5.3.2 Phase M argin . . . . . . . 52 5.3.3 Uncertainty In Parameters . . . . . 53 5.3.4 Controller Objective . . . . . . 53 6 LQR CONTROL . . . . . . 54 6.1 LQR Theory.......... ... . . ............ 54 6.2 Control Synthesis . . . . . . . 58 6.3 Nominal Closedloop Model . . . . . . 60 6.3.1 M odel . . . . . . . . 60 6.3.2 Linear Simulations . . . . . . 60 6.3.3 Gain and Phase Margins . . . . . .... 63 6.4 Perturbed Closedloop Model . . . . . . 64 6.4.1 M odel . . . . . . . . 66 6.4.2 Linear Simulations . . . . . . 69 6.4.3 Gain and Phase Margins . . . . . .... 72 7 NONLINEAR SIMULATIONS . . . . 74 7.1 Nonlinear Simulations for Nominal System . . . . 74 7.2 Nonlinear Simulations for Perturbed System . . . . 79 8 CONCLUSION . . . . . . 83 8.1 Sum m ary . . . . . . . . 83 8.2 Future Work.......... ... . . ............ 83 APPENDIX A REFERENCE FRAMES AND ROTATION MATRICES . 84 B NUMERICAL TECHNIQUES. . . . . 86 B.1 Interpolation of Force Data . . . . . . 86 B.I.1 Extrapolation Scheme .. . . . . . 86 B .1.2 C avitator . . . . . . . 87 B .1.3 F in s . . . . . . . . 87 B.2 Numerical Linearization .................. ........ 88 REFERENCES . . . . . . 90 BIOGRAPHICAL SKETCH . . . . . 91 LIST OF TABLES Table 5.1 5.2 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 B.1 B.2 page . . 46 . . 5 1 . . 64 . . 67 . . 67 . . 67 . . 68 . . 68 . . 68 . . 69 . . 70 Gain and Phase Margin for Perturbed Closedloop System: 20% error in clfin Grid For Experimental Cavitator Data .. ................. Grid For Experimental Fin Data .. ................... Control Param eters . . .. . . . Control Constraints . . . . . . Gain and Phase Margin with LQR Controller . ...... Percentage Variation in A Matrix due to 20% Variation in cc . Percentage Variation in B Matrix due to 10% Variation in c . Percentage Variation in A Matrix due to 20% Variation in cdc . Percentage Variation in B Matrix due to 20% Variation in cdc . Percentage Variation in A Matrix due to 20% Variation in clfin . Percentage Variation in B Matrix due to 20% Variation in clfin . Percentage Variation in A Matrix due to 20% Variation in cdfin . Percentage Variation in B Matrix due to 20% Variation in cdfin 1 LIST OF FIGURES Figure page 1.1 Tip Vortex Cavitation . . . . . . . 2 1.2 Formation of Cavity. . . ................... 3 1.3 Supercavitating Vehicle . . . . . . . 4 2.1 Supercavitating Vehicle . . . . . . . 7 2.2 Cavitator and Fins . . . . . . . 9 3.1 Bodyfixed and Inertial Frames .. . . . . . 11 3.2 Principle Planes of Symmetry for the Torpedo . . .. . 12 3.3 Euler Angles of Rotation ................... . . 12 3.4 Cavitator Reference Frame . . . . . . 13 3.5 Rudder and Fin Reference Frames . . . . . 14 3.6 Rudder 1 Fin Reference Frames . . . . . 16 3.7 Rudder 2 Fin Reference Frames . . . . . 16 3.8 Elevator 1 Fin Reference Frames . . . . . .... 17 3.9 Elevator 2 Fin Reference Frames . . . . . 17 3.10 Angle of Attack (a) and Sideslip () . . . . . 18 3.11 Cavitator: (a) Angle of Attack and Sideslip and (b) Hydrodynamic Forces 25 3.12 Fin G eom etry . . . . . . . . 28 5.1 Simulink Model for Open Loop Simulation . . . . 48 5.2 OpenLoop Response for Torpedo: w,p,q . . . . 48 5.3 Variation of Eigenvalues with Change in Velocity . . . 50 5.4 Loop G ain . . . . . . . . 53 6.1 Controller for Tracking when Plant has an Integrator. . . . 57 6.2 Controller for Tracking when Plant has no Integrator. . . . 57 6.3 Eigenvalues for the Closedloop System . . .. . 60 6.4 Pitch Command Tracking for Linear System : q,u . . . 61 6.5 Pitch Command Tracking for Linear System : c, c . .. 61 6.6 Pitch Command Tracking for Linear System : ei, e . . .. 62 6.7 Roll Command Tracking for Linear System : p,r . . . 62 6.8 Roll Command Tracking for Linear System : 6r1, 6 . . 62 6.9 Breakpoints for Calculating the LoopGain for a Tracking Controller . 63 6.10 Gain and Phase Margin: Longitudinal Outerloop . . . 64 6.11 Gain and Phase Margin: Longitudinal Innerloop . . . 65 6.12 Gain and Phase Margin: Longitudinal Allloop ..... . . 65 6.13 Gain and Phase Margin: Lateral Outerloop.. .. ....... .. 66 6.14 Gain and Phase Margin: Lateral Innerloop . . . . 66 6.15 Gain and Phase Margin: Lateral Allloop . . . . 69 6.16 Eigenvalues for the Perturbed Closedloop System: 20% Error in clfin 70 6.17 Pitch Command Tracking for Perturbed Linear System : q, u . . 71 6.18 Pitch Command Tracking for Perturbed Linear System : 8c, . 71 6.19 Pitch Command Tracking for Perturbed Linear System : 6ei, e . 71 6.20 Roll Command Tracking for Perturbed Linear System : p,r . . 72 6.21 Roll Command Tracking for Perturbed Linear System : 1, . 72 7.1 Complete Nonlinear Simulation with LQR Controller . . .. 75 7.2 Pitch Command Tracking : p, q .. . . . . . 76 7.3 Pitch Command Tracking: . . . . . 76 7.4 Pitch Command Tracking: e,,pe . . . . . 76 7.5 Pitch Command Tracking : 81, {x,y,z} Trajectory . . 77 7.6 Roll Command Tracking: p, q .. . . . . . 77 7.7 Roll Command Tracking: . . . . . 78 7.8 Roll Command Tracking: 8ei,e . . . . . 78 7.9 Roll Command Tracking: ,1,,1 . . . . . 78 7.10 Roll Command Tracking: {x,y,z} Trajectory . . . . 79 7.11 Roll & Pitch Command Tracking: p, q . . . . . 79 7.12 Roll & Pitch Command Tracking: .. .. .. ..... 80 7.13 Roll & Pitch Command Tracking: 1, . . . . 80 7.14 Roll & Pitch Command Tracking: ,1, 1 . . . . 80 7.15 Roll & Pitch Command Tracking: {x,y,z} Tracking .. . ...... 81 7.16 Response for 20% Variation inclfi: u,w . . .. . 81 7.17 Response for 20% Variation in clfi: p, q . . .. . 82 7.18 Response for 20% Variation in clfi,: ,1 . . . . 82 7.19 Response for 20% Variation in clfi,: {x,y,z} Trajectory . . 82 A.1 Rotation of Frames.. . . ................... 84 B.1 Shape Function for One Dimensional Quadratic Scheme . . 87 Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Master of Science CONTROL STRATEGIES FOR SUPERCAVITATING VEHICLES By Anukul Goel December 2002 Chair: Andrew J. Kurdila Cochair: Richard C. Lind Major Department: Mechanical and Aerospace Engineering Underwater travel is greatly limited by the speed that can be attained by the vehicles. Usually, the maximum speed achieved by underwater vehicles is about 40 m/s. Supercav itation can be viewed as a phenomenon that would help us to break the speed barrier in underwater vehicles. The idea is to make the vehicle surrounded by water vapor while it is traveling underwater. Thus, the vehicle actually travels in air and has very small skin friction drag. This allows it to attain high speeds underwater. Supercavitation is a phe nomenon which is observed in water. As the relative velocity of water with respect to the vehicle increases, the pressure decreases and subsequently it evaporates to form water vapor. Supercavitation has its drawbacks. It is really hard to control and maneuver a su percavitating vehicle. The first part of this work deals with modeling of a supercavitating torpedo. Nonlinear equations of motion are derived in detail. The latter part of the work deals with finding a controller to maneuver the torpedo. A controller is obtained by LQR synthesis. For the synthesis, it is assumed that the cavity is fixed and the torpedo is situated symmetrically in the cavity. The robustness analysis of the LQR controllers is carried out by calculating the gain and phase margins. Simulations of the response for a perturbed system also have been studied. CHAPTER 1 INTRODUCTION Achieving high speeds is an important issue for underwater vehicles. Even the common fastest underwater vehicles are restricted to travel at speeds around 40 ms1. The reason for this restriction is the drag induced by skin friction. When a body moves in a fluid, a layer of the fluid clings to the surface of the body and is dragged with it. This interaction causes high drag forces on the body and is commonly termed skin friction drag. The drag force in water, unlike air, is dominated by skin friction drag as compared to other sources such as pressure drag. Over the years, extensive research has been done to boost the speed of underwater vehicles. However, most of the studies were mainly focused on streamlining the bodies and improving their propulsion systems. Even though these have proven to give advancements in speed, there has not been a considerable reduction in skin friction drag. In the late 1970's, the Russian Navy was able to engineer a torpedo, the Shkval (squall) [1], that achieved a speed over 80 ms1. This phenomenal improvement in speed was made possible by supercavitation. The idea was to envelop the torpedo in a gas so that it has little contact with water. The Shkval was able to achieve a tremendous reduction in skin friction drag and exhibit high speed. 1.1 Cavitation As the speed of an underwater vehicles increases, i.e., as the relative velocity of water with respect to the vehicle increases, the pressure decreases. The speed can be increased to the point the pressure goes below the vapor pressure of water and then bubbles of water vapor are formed. This process is known as cavitation. Cavitation is mostly observed at sharp covers of a body where the speed can reach high magnitudes. A classic example for cavitation is at the tip of propellers, like the one shown in Figure 1.1. Since the propeller 2 is rotating at high speed, the relative velocity at the tips is high enough so that water at the tips vaporizes. Cavitation has been extensively researched, but remains a challenge for underwater vehicles. Figure 1.1 Tip Vortex Cavitation [2] Cavitation can be harmful in many cases. The cavitation region is usually turbulent. When the cavitation is not steady, the pressure drop is momentary and very quickly the cavitation bubble encounters a region of high pressure that forces the bubble to collapse like a tiny bomb. This collapse causes high levels of noise and also may cause considerable damage to the surface of the body. Figure 1.2 shows the various stages of formation of cavity. It shows a bulletlike body traveling in water at high speed. The cavitation starts as vapor bubbles near the region of high relative velocity. As the speed is further increased, the bubbles merge to form a large area of vapor. On further increase in speed, the whole of the body is covered in vapor. This stage is called the supercavitation. At this point, the condition is similar to traveling in air. The skin friction drag is tremendously reduced, and thus high speed can be attained in this phase. 1.2 Types of Supercavitating Projectiles The first version of the Russian Shkval was a crude supercavitating vehicle. It was unguided and had a small range of about 5 miles. Now, there are various advanced forms of supercavitating bodies. One class of supercavitating bodies, called RAMICS (Rapid PFigure 1.2 Formall tion of Cavity [3] Airborne Mine Clearance System), is a helicopterborne weapon that destroys surface and nearsurface marine mines by firing supercavitating rounds at them. These are small bullet like, flat nosed projectiles that are designed to travel stably through air and water. As the RAMICS can produce a cavitation bubble, they can travel at high speed underwater and can pierce the mines to destroy them. As they are fired from a distance in air, they need to be designed to travel in both phases. The RAMICS are better than conventional bullets, as conventional bullets are rapidly slowed down by drag in water. Another kind of a supercavitating projectile is a subsurface gun system using Adapt able HighSpeed Undersea Munitions (AHSUM). These are supercavitating "kinetickill" bullets, fired from guns fitted under submerged hull of submarines. These sonar directed bullets would be used to intercept undersea missiles. The RAMICS and AHSUM are uncontrolled small range supercavitating projectiles. The next level of supercavitating projectiles is larger torpedoes, with higher speeds and longer range. These vehicles are much more complex. They require the design of a launch station. They require detailed studies of hydrodynamics, acoustics, guidance and control, propulsion, etc. Even if the vehicle is designed to be an uncontrolled torpedo, it is still a challenge to produce and maintain a cavity around the vehicle. The cavity is usually produced by the nose of the vehicle, which is specially shaped for the purpose. The nose is called a cavitator. The cavitator may not be sufficient to produce the cavity. Thus air can be blown at the nose and various sections of the body so as to maintain and produce the IJmillcd t.aUl  r^Try Figure 1.2 Formation of Cavity [3] Airborne Mine Clearance System), is a helicopterbome weapon that destroys surface and nearsurface marine mines by firing supercavitating rounds at them. These are small bullet like, flat nosed projectiles that are designed to travel stably through air and water. As the RAMICS can produce a cavitation bubble, they can travel at high speed underwater and can pierce the mines to destroy them. As they are fired from a distance in air, they need to be designed to travel in both phases. The RAMICS are better than conventional bullets, as conventional bullets are rapidly slowed down by drag in water. Another kind of a supercavitating projectile is a subsurface gun system using Adapt able HighSpeed Undersea Munitions (AHSUM). These are supercavitating "kinetickill" bullets, fired from guns fitted under submerged hull of submarines. These sonar directed bullets would be used to intercept undersea missiles. The RAMICS and AHSUM are uncontrolled small range supercavitating projectiles. The next level of supercavitating projectiles is larger torpedoes, with higher speeds and longer range. These vehicles are much more complex. They require the design of a launch station. They require detailed studies of hydrodynamics, acoustics, guidance and control, propulsion, etc. Even if the vehicle is designed to be an uncontrolled torpedo, it is still a challenge to produce and maintain a cavity around the vehicle. The cavity is usually produced by the nose of the vehicle, which is specially shaped for the purpose. The nose is called a cavitator. The cavitator may not be sufficient to produce the cavity. Thus air can be blown at the nose and various sections of the body so as to maintain and produce the cavity. Figure 1.3 shows a supercavitating torpedo traveling underwater. It can be seen that the whole of its body is enveloped by a cavity. This is the kind of vehicle that has been studied in this work. Figure 1.3 Supercavitating Vehicle [1] 1.3 Related Research Research in the field of supercavitation has been going on from the early 1900's. But earlier research was aimed at reduction of cavitation so as to reduce noise and body damage. In the 90's the focus shifted to exploiting the effects of supercavitation. The work shown in May [4] has an extensive collection of experimental data for vari ation of forces on various supercavitating shapes. Coefficients of lift and drag are plotted with the variation of cavitation number for shapes like disks, cones, ogives and wedges. The work done in this research makes use of a CFD database provided in Fine [5]. This database has values for coefficients of lift and drag for conical cavitators, which are func tions of the half angle of the cone and the angle of attack. This database also has coefficients of lift and drag for wedges, which are functions of wetted surface of the wedge, angle of attack and sweepback angle (we discuss the definition of these geometric quantities such as the angle of attack and half angle later in this thesis). This information is useful to calculate the forces on fins of the torpedo. In late 90's a lot of research has been done on the dynamics of supercavitating vehicles. Work shown in Kulkami and Pratap [6] and Rand et al. [7] deals with studying dynamics of uncontrolled supercavitating projectiles. A dynamic model for RAMICS and AHSUM has been developed. It is shown that these projectiles rotate inside the cavity. This rotation leads to impacts between the tail of the projectile and the cavity wall. The frequency of the 5 impact increases with time. These projectiles are very short range and have a small time of flight, on the order of a few seconds. The work shown in Dzielski and Kurdila [8] focuses on the formulation of a control problem for a supercavitating torpedo. A dynamical model for a fin controlled torpedo has likewise been developed. The model also includes a formulation for the cavity. It is observed that the weight of the body forces it to skip inside the walls of the cavity and the vehicle is unstable. A control system is designed for the torpedo and results of closedloop simulations have been presented. 1.4 Overview of this Thesis This work aims at formulating the control design for a supercavitating torpedo. Equa tions of motion for the torpedo are derived, and linear control methodologies are applied for pitch and roll control of the torpedo. The main purpose of this work is to analyze these controllers and ascertain their robustness to various errors and uncertainties. Chapter 2 of this thesis briefly describes the configuration of the supercavitating torpedo for which the equations of motion have been developed. A detailed derivation of the equations of motion for the torpedo has been carried out in Chapter 3. Various reference frames have been used to obtain the kinematic equations of the torpedo. Dynamic equations are derived using Newton's Laws. Various forces experienced by different regions of the torpedo have been explained. Chapter 4 describes linearization of the equations of motion using small disturbance theory. Numerical methods are used for this purpose. It is observed that the linearization for a simple trim can be very complicated. Chapter 5 describes the control synthesis for the torpedo. Openloop dynamics are shown. The closedloop problem and various constraints on the torpedo have been stated. Chapter 6 formulates a Linear Quadratic Regulator (LQR) control design for the tor pedo. Controllers for pitch and roll rate control of the torpedo are obtained. The results for linear closedloop system and a perturbed liner system have been shown. 6 The results for pitch and roll rate control for the nonlinear closedloop system have been presented in Chapter 7. The simulations for a perturbed system model also have been presented. CHAPTER 2 CONFIGURATION OF VEHICLE Although supercavitation can be a very helpful phenomenon, it presents significant challenges in modeling and control of supercavitating vehicles. As a significant portion of the vehicle is located in the cavity, the control, guidance and stability of the vehicle has to be managed by very small regions in front and aft of the vehicle. Also generation of the cavity can be a significant problem. The major problems associated with the supercavitat ing vehicles may be summarized as: * generation and maintenance of cavity * balancing the weight of the vehicle * control and guidance * stability Figure 2.1 is a figure of a supercavitating torpedo that is modeled in this work. The main parts of the torpedo are the cavitator in the front and the four fins in the aft portion. The cav itator is used to generate and maintain the cavity. The Cavitator and the four fins together are also used for control and stability of the vehicle. 2.1 Cavitator The cavitator is the element that generates a cavity around the torpedo. Several types of cavitators, including cones, wedges and plates have been investigated [4]. This project will consider a conical cavitator as shown in Figure 2.1. The main parameter that defines Figure 2.1 Supercavitating Vehicle [8] Figure 2.1 Supercavitating Vehicle [8] a conical cavitator is its halfangle. The coefficients of lift and drag, as functions of half angle, for the cavitator have been generated using CFD and tabulated in [5]. The cavitator in this model has one degree of freedom defined by a rotation angle about an axis perpendicular to its axis of symmetry. The shape and location of the cavity is a nonlinear function of this cavitator deflection angle and half angle of the cone. This shape determines the position where the cavity intersects the body of the vehicle. Thus, the amount of wetted area of the vehicle is dependent on these angles, which in turn determines the efficiency of supercavitation achieved. As stated earlier, a large portion of the vehicle is in the cavity. The lift produced by the cavitator combined with the lift produced by the fins helps in balancing the weight of the body. 2.2 Fins Although the cavitator is capable of providing enough lift to sustain the body in water, it does not usually provide sufficient forces and moments to stabilize and control the ve hicle. For this purpose fins are required in the aft portion of the vehicle. The fins help in counteracting the moments produced by the cavitator and also providing some amount of lift to support the weight of the body. There are four fins placed symmetrically along the girth of the vehicle, near the tail. The fins also are the control surfaces, as they can provide differential forces, thus causing control moments. Two of the fins shown in Figure 2.2 are horizontal (placed parallel to the axis of rotation of cavitator), called elevators, are used to affect the longitudinal dynamics for the vehicle. The other two fins are called the rudders and are used to affect the lateral dynamics to the vehicle. The fins have two degrees of freedom, a sweepback angle and an angle of rotation. These angles will be explained further in Chapter 3. 2.3 Maneuvering The motion of the vehicle is controlled by all five control surfaces, the four fins and the cavitator. In a straight line motion the cavitator and the elevators balance the weight of the Figure 2.2 Cavitator and Fins vehicle. A propulsion force at the tail pushes the vehicle forward. The rudders usually have a zero deflection in such a case. The vehicle depends on a banktoturn strategy for making a turn, and cannot use tradi tional missile strategies such as skidtoturn. This dependence arises because the hull of the vehicle is incapable of providing a lift force. The fins are the main lift generating surfaces. A differential force from the fins can be used to generate a force towards the center of the turn. CHAPTER 3 NONLINEAR EQUATIONS OF MOTION The dynamics of the torpedo, whose configuration was discussed in Chapter 2, can be mathematically formulated. A complete derivation of the equations of motion is given in this chapter. Section 3.1 deals with the setup of reference frames and derivation of the kinematic equations. The dynamic equations of motion are derived in Section 3.2. 3.1 Kinematic Equations of Motion The definition of a suitable coordinate system and degrees of freedom is a precursor to any formulation of equations of motion. The derivation of the equations of motion of multibody systems is highly simplified by defining various reference frames and relations between them. Appendix A describes briefly the concept of reference frames and rotation matrices. These concepts will be used extensively in the derivation of equations of motion. The derivation of the equations of motion will be done in two parts. First, the kinematic equations will be derived. These include the formulation of the position vectors, velocities and accelerations of various points on the torpedo. Next, the dynamic equations will be derived. Finally, Newton's Laws yield the dynamic equations of motion. 3.1.1 Orientation of the Torpedo Six degrees of freedom (DOF) are required to describe the position and orientation of the torpedo when it is considered a rigid body. Of these, three DOFs give the location of a point fixed on the torpedo. The other 3 DOFs give the orientation of the torpedo in a fixed reference frame. The position of the torpedo in a reference frame can also be obtained by the integration of its velocity in that reference frame. The torpedo is assumed to be moving in an earthfixed reference frame E, centered at any conveniently chosen point and described by the basis vectors (el e^, 3). The earthfixed 3 e Figure 3.1 Bodyfixed and Inertial Frames reference frame is shown in Figure 3.1. The vector 63 points in the downward direction, i.e., the direction of the gravity. The vectors e\ and j2 are placed in the horizontal plane such that the basis vectors form a righthanded coordinate system. As shown in the figure, ei points to the right and ^2 points outside the plane of the paper. A bodyfixed frame B is defined by the basis vectors (bi, b2, b3) so as to simplify the derivation of the equations of motion. The frame B is centered at 0, the center of gravity of the torpedo, and moves with the torpedo. The reference frame B is shown in the Figure 3.1. It can be seen in Figure 3.2 that the torpedo has two planes of symmetry namely b1b2 and b1b3. The plane b1b3 is called the longitudinal plane and plane b1b2, the lateral plane. Transformation from E to B can be achieved by 3 rotations. Many such sets of rotations are possible. The rotation steps chosen here are the Euler angles of rotation, which are the yaw, pitch and roll. As there are three rotations to be performed, two intermediate reference frames with basis vectors (21,42,43) and (C1i,2,j3) will have to be introduced to perform the transformation. The rotations, to be performed in order, are listed below. 1. Rotate the frame E about 63 through a yaw angle, ?, to obtain the frame (21,42,23). 2. Rotate (h1,22,h3) about x2 through a pitch angle, 0, to obtain the frame (\1,Y2,Y3). 3. Rotate (i1,J2,j3) about Uf through a roll angle, (, to obtain the frame B. Figure 3.2 Principle Planes of Symmetry for the Torpedo 3.3 Figure 3.3 Euler Angles of Rotation Figure 3.3 shows the above rotations in order. Based on the above definition of rotations, the transformation matrix from E to B can be written as in equation 3.1. 1 0 0 CO 0 so CY SwT 0 C' SO 0 1 0 ST CY 0 SO Cc SO 0 CO 0 0 COCY COST SO CYTSSO COST SISOST + CYCc SICO CISOCY + SO ST CSSOST CS C(CCO 0 fe 0 j2 1 83 j ] e2 e3 bi b2 b3 (3.1) 3.1.2 Orientation of the Cavitator As described earlier, the cavitator has only one degree of freedom. It can rotate in the longitudinal plane about an axis parallel to the b2 axis. The orientation of the cavitatorfixed axes with respect to the body fixed axes is shown in Figure 3.4. A.) bi Ap e ACP b3 ^ \ A b3 Figure 3.4 Cavitator Reference Frame The deflection of the cavitator is given by an angle, 8c, which is positive for a positive cavitator rotation about the b2 direction. The geometric center of rotation of cavitator is denoted by P. CP is the center of pressure for the cavitator, which is at a distance Acp from P, along c1. From Figure 3.4, the rotation matrix from B to cavitator fixed frame C, can be written as in Equation 3.2. ei C8c 0 SSc bh c2 = 0 1 0 b2 (3.2) C3 SSc 0 CSc h3 k ^ k > 3.1.3 Orientation of Fins Figure 3.5 shows the orientation of the finfixed reference frames. For convenience, all the fin frames are indicated by basis vectors (11,12,~3). In text they will be represented as (f1,f2,f3) fn, where subscript fin corresponds to a particular fin. Rudder (1 fRu Rudder 22 FRONT VIEW f2 TOP VIEW Elevator l EAt Elevator 2 Figure 3.5 Rudder and Fin Reference Frames All the fins have 2 DOFs, namely the sweepback angle (Ofi,) and the fin rotation (fi,). These can be defined as Sweepback angle (Ofi,): This parameter is the angle of rotation of a fin about its f3 axis. The direction of rotation for positive sweepback is such that the fin is moved toward the rear portion of the torpedo. Due to this definition, the positive sweepback is along the negative 3 direction for all fins. Sweepback angle determines the amount of fin that is enveloped in the cavity. a, f * Fin Rotation (8,i,): This parameter is the angle of rotation of the fin about its f2 axis, in positive the f2 direction. Fin rotation determines the magnitude and direction of the forces acting on the fins. The order of rotation in the above case is important to obtain the correct equations. Sweepback has to be performed before fin rotation. An intermediate reference frame G with basis vectors (gI, 2,g 3) is introduced so as to simplify the derivation of rotation ma trix from B to the fin coordinates. Sweepback aligns the finfixed frames with the interme diate frame G and then a fin rotation puts the fin frame in actual orientation with the fins. As the second rotation is identical in all cases, the transformation matrix from frame G to fin frame Ffi, can be written as in Equation 3.3. l C8fin 0 S8fin gl f2 = 0 1 0 I 2 (3.3) f3 fin S8fn 0 Cafn1 g3 fin The rotation matrix for sweepback and the transformation matrices from B to Ffin frame for each of the fins can be derived easily, and are summarized below. * Rudder 1 Figure 3.6 shows the sweepback and fin rotation for rudder 1. The matrices for transformation from B to R1 can be derived as in Equation 3.4 and Equation 3.5. Si [ COR1 0 SOR1 I1 S2 = SORI 0 COR i 2 (3.4) g3 R1 0 1 0 J3 fCR1 0 SR1 C 0 SOR f i 2 0 1 0 SOR1 0 COR1 2 (3.5) S3 SSR1 0 CR1 0 1 0 3 Rudder 2 Figure 3.7 shows the sweepback and fin rotation for rudder 2. The matrices for transformation from B to R2 can be derived as in Equation 3.6 and Equation 3.7. g COR2 0 SOR2 ~j S2 = SOR2 0 COR2 2 (3.6) g3 1 0 b3 SCs8 0 S8R2 COR2 0 S 2 0 0 1 0 SSR2 0 COR2 2 (3.7) 1i3 R SSR2 0 CR2 0 1 0 ( 3 OR1 gz,f2 6R1 f 4 g3 f3 Figure 3.6 Rudder 1 Fin Reference Frames Figure 3.7 Rudder 2 Fin Reference Frames * Elevator 1 Figure 3.8 shows the sweepback and fin rotation for Elevator 1. The matrices for transformation from B to El can be derived as in Equation 3.8 and Equation 3.9. 0 b\ 0 b2 COEI SOE1 SOE1 COE1 0 0 0 b2 1 3 3 * Elevator 2 Figure 3.9 shows the sweepback and fin rotation for Elevator 2. The matrices for transformation from B to E2 can be derived as in Equation 3.10 and 61 3 g2 g3 El COE1 SOE1 0 CSEl 0 S8EI SOE1 COEI 0 S8E1 0 C E1 (3.8) (3.9) g2,f2 fl \ g3 f3 Figure 3.8 Elevator 1 Fin Reference Frames b,2 E2 3 31 E >E2 0E2 Figure 3.9 Elevator 2 Fin Reference Frames ; ________^^ yI Figure 3.9 Elevator 2 Fin Reference Frames Equation 3.11. SOE2 0 r COE2 0 b2 0 1 3 SSE2 COE2 0 SOE2 CSE2 0 SOE2 COE2 0 10 b2 1 I 3 Equations 3.2 to 3.11 will be used in later sections to transform the forces on fins and cavitator to the bodyfixed frame. g2 g83 E2 f E2 I j3 P 2 SCOE2 SOE2 0 C8E2 0 0 1 S8E2 O (3.10) (3.11) 3.1.4 Angle of Attack and Sideslip The bodyfixed reference frame has been defined, but the velocity of various points on the body of the torpedo is yet to be defined. The torpedo is considered as a rigid body. If the velocity of a certain point on a rigid body is known, the velocity at any other point on that body can be found by knowing the rotation matrices. Thus, V = ub1 + vb2 + wb3 will be taken as the velocity of CG of the torpedo. The velocity of other points on the torpedo can be defined subsequently. Two very useful parameters, angle of attack and angle of sideslip can be defined in conjunction with the orientation of the body axis with the velocity of CG. Figure 3.10 shows these parameters and their geometric interpretation. V ubi +( v)(b ,)+wb, V b3 Figure 3.10 Angle of Attack (a) and Sideslip (3) Angle of attack, a, can be defined as the angle between the projection of velocity V onto b2b3 plane and the b1 axis. Angle of attack is positive when the nose of the torpedo points above the velocity vector of the torpedo. As before, an intermediate frame F given by (, f12,13) can be described by rotation of the B frame by an angle a. Thus the rotation matrix from Fbody to B can be written. by Ca 0 Sa fi b2 = 0 1 0 f2 (3.12) b3 Sa 0 Ca f3 J body The Angle of sideslip, 3, is defined as the angle between the actual velocity V and the projection of V onto b2b3 plane. Again, a frame Gbody can be defined by rotation of Fbody by an angle equal to 3 in negative f3 direction, thus giving a rotation matrix as in Equation 3.13. Sg2 S= S C3 0 /2 (3.13) 3 body 1 3 bo Velocity of CG of the torpedo in the Gbody frame can be written as Vgl, where V is magnitude of V. It will be seen that drag and lift on the torpedo can be obtained in this frame. Thus a transformation from Gbody to B is important. It is given by Equation 3.14. b1 CPCaC C cSp S( ]i b2 = S3 C3 0 g2 (3.14) b3 C3Sc SS3 Ca J body3 Using Equation 3.14 V can be rewritten as in Equation 3.15. V = Vgl (3.15) = VCp3Ca b VS b2 VCS b3(3.15) H V W where V2 = V2 = 2 +2 +w2. From Figure 3.10, relations between the velocity compo nents and the angles of attack and sideslip can be derived. w tan = (3.16) u v sin 3 = (3.17) Though the matrix GbodyB in Equation 3.14 has been defined for the bodyfixed ref erence frame and velocity of CG of the torpedo, the equation is valid for any other rigid part of the body like the fins and cavitator. Thus, in case of a fin, the velocity V would correspond to a point (like the tip, center of pressure etc.) on that fin, and Gfin_B matrix would correspond to the finfixed reference frame. In this case the velocity of center of pressure of the fin will be used to define the above parameters. Thus, obtaining afi, and 3fin is a two step process: 1. Obtain the velocity of center of pressure of fin. VCPbody = Vcg +E oB X rcgCP (3.18) where Vcpbody is velocity of CP of fin in B frame, Vcg is the velocity of CG of the torpedo in E frame, E o is angular velocity of B in E, and rcgcp is position vector from CG to CPf,. Equation 3.18 can be rewritten as in Equation 3.19. { fin bl b2 b3 Vfpin, v + p q r (3.19) Wfin B w g Xfin Yin Zin where rcgCp = Xfinbl + Yinb + Zfinb3. 2. Transform the velocity (in E) of CP of fin from frame B to frame of corresponding fin. This transformation is obtained by using rotation matrices derived in Equations 3.3 to 3.11. UR1 VR1 WR1 UR2 VR2 il12 }: R2 C6R1 0 S6R1 C8R2 0 S8R2 SSR1 0 S8R2 0 C6R2 UEI CsE1 0 SE1I VE 1 0 1 0 WE1 l E SE1 0 C8E1 UE2 C6E2 0 S8E2 VE2 0 1 0 1, E2 S8E2 0 C8E2 The left hand terms in Equations 3.20 COR1 SOR1 0 COR2 SOR2 0 COEI SOEI 0 COE2 SOE2 0 0 SORI UR1 0 CORI VR1 1 0 WR1 B 0 SOR2 U R2 0 COR2 VR2 1 0 12 RB SE1 o0 UEI COE1 0 VE1 0 l1 WE1 B SOE2 0 UE2 COE2 0 VE2 0 1 lI B to 3.23 give the velocity components at the CP for corresponding fins, in that fin frame. These can be used in Equations 3.16 and 3.17 to find the angle of attack and sideslip for a particular fin. 3.1.5 Kinematic Equations Velocity of the CG of the torpedo has been defined in the previous section. That velocity defines the translational kinematics for the torpedo. Analogous to this quantity, angular velocity is required to define the rotational kinematics. The angular velocity of the hull has components p, q and r in the frame B. EoB Apb + qb2 + rB3 (3.24) (3.20) (3.21) (3.22) (3.23) As the transformation from E to B has already been defined in terms of rotational motions give by Euler angles, the angular rates can also be obtained by differentiation of Euler angles. Thus, another form of Equation 3.24 can be written as in Equation 3.25. EB = Ye3 + Ox2 + [1 (3.25) The rotation matrices from Equation 3.1 can be substituted into Equation 3.25 to obtain Equation 3.26. E oB = (6 SO )bi + (YC SO + OCD)&b2 + (COCO OS) b3 (3.26) Equations 3.24 and 3.26 can be equated to obtain Equation 3.27. p So 0 1 Y q = CSO COS 0 (3.27) r COCO SO 0 6 Equation 3.27 can be rewritten as in Equation 3.28. p = SOY q = TCOSS + OCO (3.28) r = TCOCO OSO To apply Newton's Laws, acceleration of the CG is required. The values of p, q, r will be the angular accelerations of torpedo in B. The translational acceleration can be calculated by time differentiation of V in Newtonian frame. A differentiation formula can be used to find the time derivative, in some frame, for a vector defined in some other related frame. S(v) = (v) +I o xv (3.29) dt I dt B where, subscript I denotes Newtonian (inertial) frame, and B is the bodyfixed frame. I'B is angular velocity of the body (or bodyfixed frame) in the I frame, v is the velocity in I frame, of the point where acceleration is desired. Using the formula the acceleration of CM of torpedo in E can be obtained. f b Z2 3 EaCM + p q r (3.30) w u v w u +qw vr = + ur pw (3.31) w + pv uq Similarly, the rotational acceleration will be required in the frame E. This can be written analogous to Equation 3.30. b, bi2 b3 EOI= + p q r r p q r (3.32) P = q r The position of torpedo can be found by integrating the velocity. Let (x,y, z) represent the coordinates of CG in frame E. Thus, the time derivative of these coordinates in E should represent the velocity components of the torpedo in E frame. When these time derivatives are transformed to bodyfixed frame, they would be equivalent to the velocity components in bodyfixed frame. x u y = v (3.33) Ez w Equation 3.1 can be substituted in Equation 3.33 to obtain Equation 3.34. x COCY COsT so u z C)SOCY +SST CS OSY SOCY CCOC w 3.2 Dynamic Equations of Motion Now that the accelerations of various parts of the torpedo are known, Newton's Laws can be used to derive the dynamic equations of motion. Newton's laws state that the rate of change of momentum is equal to the sum of external force applied on the body. Equation 3.35 can be obtained from Newton's laws by an assumption that the mass of the torpedo is constant. This assumption is valid for a short period of time. The equations are YF=P (3.35) = ma where P is the linear momentum, m is mass of the body, a is the acceleration and F is all the forces of the body. Using Equation 3.31 in Equation 3.35, Newton's Laws for the torpedo can be rewritten as in Equation 3.36. ui + qw vr m v+ur pw =F (3.36) w + pv uq Newton's laws can be extended to rotation. Equation 3.37 are the Newton's Laws for rotational motion. (3.37) =Ic + E (B x H where H (= IcMEoB) is the angular momentum, IcM is moment of inertia matrix of the body, a is the angular acceleration and M is all the moments on the body. It should be noted that above stated Newton's laws are only valid when the quantities are calculated in an inertial frame of reference. Thus, the quantities have been calculated in frame E. Using Equation 3.32, the Newton's Law for rotation can be written as in Equation 3.38. 1I 0 0 1p bI b2 b3 0 12 0 q + p q r =AM (3.38) 0 0 13 r Ilp I2q I3r To derive the equations, the forces on each of the parts will be calculated individually, and then expressed in bodyfixed reference frame, i.e., summation will be done in body reference frame. The rotation matrices derived in previous sections will be used extensively for this purpose. Various types of forces are experienced by a moving torpedo in water. These forces can be summarized as follows: * Hydrodynamic Forces: These are the forces exerted by the fluid on the torpedo. Thus they exist whenever the fluid is in contact with body. For supercavitating ve hicle, most of the body (hull) is inside the cavity. Only a portion of the fins and the cavitator are in contact with the fluid. Thus the lift and drag on cavitator and fins are only hydrodynamic forces. The coefficients of lift and drag for the fins and cavitator are functions of their angle of attack, immersion, sweepback angle, angle of rotation etc. and are tabulated in a database [5]. This database will be interpolated and ex trapolated to calculate the numerical values for the forces. Occasionally, a part of the hull comes in contact with water and might experience some forces. These forces are known as planing forces. It is assumed that the vehicle is centered in the cavity. Thus the planing forces are neglected in the summation of the hydrodynamics forces. FHydrodynamic = FR1 + FR2 + FE1 + FE2 + Fc (3.39) MHydrodynamic = MR1 +MR2 +ME1 +ME2 +Mc (3.40) Gravitational Forces: This is the gravity forces on the body. As the summation of moments will be taken about the center of gravity, the gravitational forces will not contribute to the summation on moments. They will appear only in summation of forces. Propulsive: The torpedo has a propulsion system, which is similar to that of rockets. The line of action of the propulsive force is assumed to be passing through center of gravity and along b1 axis. Thus this force will also contribute to the sum of forces, and not moments. The propulsive forces are provided by burning the fuel, but for simplicity it will be assumed that the mass of the torpedo remains constant. The total forces and moments are expressed in terms of these components. F = FHydrodynamic + FGrav + FProp (3.41) M = MHydrodynamic (3.42) 3.2.1 Forces on Cavitator Figure 3.11 shows the forces acting on the cavitator. Coefficient of lift (clc) and drag (cdc) for the cavitator are functions of angle of attack, ac, and halfangle, yi, of the cavi 2 tator. Halfangle, for a cone, is the angle made by axis of the cone with any line passing through the vertex and lying in the surface of the the cone. This parameter defines the main geometry of the conical cavitator. The values of cl/ and cdc are determined using CFD and tabulated in [5]. These values have been extrapolated to calculate lift and drag. Lc = pV Sccl (c, ) (3.43) 2 2 Dc = 2pV2Secdc(oc,7y) (3.44) where Sc is the crosssectional area of the cavitator base. Directions of lift (Lc) and drag (Dc) are as shown in Figure 3.11(b). These can be transformed to the body axes using 3.2 and 3.14 for the cavitator. b2 = CB(c) x GcavC(ac, 3c) x &2 (3.45) SC^2 M n Lc Dc g3 (a) (b) Figure 3.11 Cavitator: (a) Angle of Attack and Sideslip and (b) Hydrodynamic Forces Thus the forces on cavitator, in body frame, can be written. F JR Fc = Fcy Fc, B D c(ac, Yi) = 0 (3.46) Lc(ac,,7 ) 20 cav C8c 0 SSc Cp cCac CacXSc Sac Dc (ac 7) = 0 1 0 Sc CPc 0 0 sac 0 C6c CcSac ScSc Cuc Lc((c, 7) where Fc is a 3x1 matrix with each row corresponding to each direction in B basis. The forces are assumed to be acting at CP of the cavitator. Once the forces have been calculated, the moment about any point can be calculated. Mc = rcpcav X Fc (3.47) where rcPcav is the position vector from that point to CP of cavitator. It is assumed that the CP lies on b1 when cavitator deflection is 0, and its coordinates with respect to body origin O, in this case, are (Xc, 0, 0). Thus from Figure 3.4, the coordinates of CP can be written for the case when the cavitator is deflected. Xc +ACPC8c rcpcav = 0 (3.48) AcPSc body body The moments on the cavitator in bodyfixed can be obtained by substituting Equations 3.46 and 3.48 in Equation 3.47. C c 0 SSc M = [(Xc+ AcpCc) bi AcpS8c3] x 0 1 0 SSc 0 CSc (3.49) CC(,ca, Ccas Sa, D(ac, i) SIp cpc 0 0 CPcSoc( SacSPIc Cac Lc (oc,7i) 3.2.2 Forces on Fins Fin forces are also extrapolated from the CFD database [5], which gives the values of coefficients of lift (clfin) and drag (cdfin) for fins. These values are functions of angle of attack Cafi,, fin sweepback Ofi,, fin rotation 8fi,, fin immersion Ifi, and fin geometry. Figure 3.12 shows these parameters graphically, and they can be defined as follows: * Fin Geometry: The geometry parameters for fins are L and S, and wedge half angle (ha), as shown in Figure 3.12. These parameters are fixed according to the values given by the database, so as to calculate the forces accurately. Fin Immersion: As the fin is partially wetted by fluid, the wetted length can be rep resented by parameter So along fin Yaxis. The immersion Ifi, can be defined as the ratio of the wetted length of the fin to its true length. fin = (So/S)fin (3.50) Ifin determines the total hydrodynamic force on the fin. Fin Rotation (8fi,): As defined earlier, this is rotation about fin j2 axis. This deter mines the direction of the hydrodynamic force. Thus fin rotation is used as primary control for the torpedo. Fin Sweepback (Ofi,,): As defined earlier, this is rotation about fin f3 axis. Sweepback determines the wetted surface of the fin, thus the hydrodynamic force on the fin. Angle of Attack: Angle of attack for the fin is calculated as described in Figure 3.12 and section 3.1.4. The database gives clfin and cdfin as a function of (fi,, Ofin and Ifin, thus lift and drag on the fins can be calculated by the normalizing factor. 1 2 2 Lfin = 2 V SinCfin (3.51) Dfin = p2Sfincdfin (3.52) tions 3.3 to 3.11. Figure 3.12 Fin Geometry Where Sfi is the length of the fin as shown in Figure 3.12 These forces have directions as shown in Figure 3.12. Before substituting in Equation 3.39, these forces have to be transformed to bodyfixed reference frame. This process involves following two rotations: 1. Rotate the frame Ffin (which has Lfi, and Dfin along its basis vectors) to align it with the finfixed frame using Equation 3.14 and 2. Rotate the above obtained finfixed frame to obtain the bodyfixed frame using Equa tions 3.3 to 3.11. Thus the forces on the fins in bodyfixed frame axis can be obtained. * Rudder 1 FRl,x COR1 SOR1 0 CSR1 0 S8R1 FRIy 0 0 1 0 1 0 FRI,z B SOR1 COR1 0 SSR1 0 C6SR S CPRICR1 CaRISPR1 SaR1 DR1I _SR1 CfR1 0 0 CPR1SaR1 SaR1SPR1 CaR1 LR1 Rudder 2 FR2,x COR2 SOR2 0 F c8R2 0 S R2 {FR2, Y 0 0 1 0 1 0 FR2,z B SOR2 COR2 0 SSR2 0 C8R2 SCR2CaR2 Ca(R2SR2 S(R2 DR2 SCR2 CSR22 0 0 CIR2SaR2 SaR2SPR2 CUR2 LR2 J * Elevator 1 FEI,x COEI SOE1 0 CE1 0 S8Ei FE1, = SOE1 COE1 0 0 1 0 FEI,z B 0 0 1 SSEI 0 CSE CE1CEcoE1 CoESIE1 SE(1 DE 1 S3E1 CPEl 0 0 C3E1SaE1 SaE1SE1 CaE1 LE1 Elevator 2 FE2,x1 COE2 SOE2 0 CE2 0 S8E2 FE2,y = SE2 COE2 0 0 1 0 FE2z 0 0 1 S8E2 0 CSE2 (3.56) SCE2CaE2 C(E2SIE2 S(E2 DE2 SPE2 CPE2 0 0 CPE2S(E2 S(E2SPE2 C(E2 LE2 Equations 3.53 to 3.56 give the components of hydrodynamic forces on fins in bodyfixed frame. What remains is to find the moment of these forces about CG of body. The moments can be obtained in analogous to Equation 3.47. MIfin = rf_Cp X Ffi (3.57) In this case, the root of fins is fixed to body, and it can move thus moving the CP of fin relative to root. The position of CP from root is also know with respect to fin coordinates. rCroot =Xotb I+Yrb2+Zr t3 (3.58) oiotcp = Aj1 +A 2 (3.59) where (fi,2,/3) is finfixed coordinates for corresponding fin. Equations 3.58 and 3.59 can be added by using matrices given in Equations 3.3 to 3.11. Thus, the position vector from CG to CP of fins can be obtained. * Rudder 1 XR1 r1ot COR1 SOR1 0 YR1 + 0 0 1 ZR1 Z ot R1 COR1 0 SC8R1 0 SSR1 ~ Ay(3 0 1 0 Ayj S8R1 0 C8R1 0 R1 Rudder 2 SYR2 2 Yt 0 0 1 ZR2 B Z SR2 CR2 0 OL.2 (3.61) C8R2 0 SSR2 ( 0 1 0 Ayp SSR2 0 CsR2 0 R2 Elevator 1 YE I y( oot S E1I COEI 0 ZE 1 B Z oEo B O 0  0 E (3.62) C8E1 0 S8E A 0 1 0 AyCp SSE1 0 CSE1 0 El Elevator 2 XE2 Xr 20t C[ E2 SOE2 0 YE2 Y t + SOE2 COE2 0 ZE2 J R I ot I 0 1 B C8E2 0 8E2 (3.63) C8E2 0 CSE2 0 E Equations 3.60 to 3.63 give the position vector from CG to CP of the fins. These equa tions in conjunction with Equations 3.53 to 3.56, used in 3.57, gives the external moments on fins about the CG. b1 b2 b3 Mfi = Xfi. Yfin Zfin (3.64) Ffin,x Ffiny Ffin,z 3.2.3 Gravitational Forces For simplicity, mass (m) of the torpedo is assumed to be constant over time. This is not the case in reality but the approximation is reasonable for considering short time maneuvers. Acceleration due to gravity, g, is also assumed to be constant as torpedo moves in space. The direction of gravity is given by 63 axis. Thus, the gravitational force can be written as in Equation 3.65. Frav = mgd3 (3.65) Equation 3.65 can be reexpressed in body frame of reference using Equation 3.1. COCY COST SO 0 Fgrav = CYS4SO COSY SISOST + COC SICO 0 CISOCY +S ST CISOST SICY CODCO mg E (3.66) So = mg SOCO C(CO B 3.2.4 Equations of Motion Now that a mathematical formulation of forces on the torpedo has been achieved, these equations can be substituted into Equations 3.36 to 3.42 to obtain the dynamic equations of motion. Thus the force equations of motion can be summarized as in Equation 3.67. 3.2.4.1 Force Equations i + c qwvr Fprop FRx m r + u pWr = Fimmersion + 0 + FR y + w+ pv uq 0 FR1,z FR2,x FE1,x FE2,x SO FR2,y + FEly + FE2y +mg SOCO + (3.67) FR2,z FE1,z FE2,z COI C CSc 0 S81 CpcCac CacSpc Sc1 Dc(aey) 0 1 0 SPc CPc 0 0 Sc 0 C C CP ccSc S4cSPc Cac Lc(ac,7y) Some issues should be noted for Equation 3.67. The planing forces have not been included in the equations of motion. These forces are neglected by assumption that the vehicle is centered in the cavity. The propulsion force is constrained to be along negative b\ axis. 3.2.4.2 Moment Equations Il 0 0 p b b 2 b3 bI b2 b3 0 12 0 Y + p q r = XR YR1 ZR + 0 0 13 Ilp Iq I3r FRI,x FRIy FR1,z b1 b2 b3 b1 b2 b3 XR2 YR2 ZR2 + XE1 YE1 ZE1 + (3.68) FR2,x FR2, FR2,z FEI,x FEly FEI,z b1 b2 b3 b1 b2 b3 XE2 YE2 ZE2 + Xc + ApC8c 0 AcpSSc FE2,x FE2y FE2,z Fc,x Fcy Fc,z Some issues should be noted for Equation 3.68. * Some of the terms in Equation 3.68 are shown as a determinant. They need to be expanded and written as components in bodyfixed frame so as to equate the left hand and righthand terms. Moments due to gravitation do not appear because the moments are taken about CG. Again, the moments due to planing forces have been neglected. 3.2.4.3 Orientation Equations j } 0 C P S1 SSO CO r 3.2.4.4 Position Equations x COCY COsT so u y = CS~YSOC SCY SISOSY +CYCI SCO v (3.70) z CASOCY +SASY C(SDSY SOCY COCO w Equations 3.67 to 3.70 are a complete set of equations of motions for the supercavitating torpedo. CHAPTER 4 LINEARIZATION 4.1 Linearization 4.1.1 Need for Linearization The equations of motion for the torpedo are identical to airplane equations of motion but the forces terms on the righthand side of these equations are different. Thus, the linearization can be carried out similarly, as shown in Nelson [9]. The equations of motion, as in the case of a supercavitating torpedo, are represented by a set of firstorder differential equations. x f(x, u) (4.1) using f : 9)" 9" as a nonlinear function of a timevarying vector x E 9" and u E 9". For control design, the system dynamics are observed at some trim conditions by giving perturbations to states of the system at that trim. The dynamics associated with these perturbations are obtained by linearization. An advantage by linearization is that most of the control methodology is based on linear equations of motion. A controller is designed initially for the linear system and then tested for the actual nonlinear system. Yet, there are few disadvantages for this process * Linearized equations can predict the system performance only in a small range of op erations. The linearized equations change as the operating point of system changes, thus making it difficult for simulating true behavior of system. Information relating to nonlinearities like hysteresis, backlash, coulomb friction, dis continuities etc. may be lost by linearizing the system. A controller that is good for linearized system might have very poor performance for the nonlinear equations. An iterative process may be needed to find a controller that is good for nonlinear equations. 4.1.2 Generic Form of Equations of Motion The equations of motion in Equations 3.67 and 3.70 can be written in simplified form using sums of total forces and moments acting on the body. m(u + qw vr +gSO) = X m(v + ru pw gCOSc) = Y (4.2) m(w + pv qu gCOCO) = Z IP +qr(IIy) =L (4.3) Iyq+rp( ) =M (4.4) Ir + pq(Iy I) =N (4.5) 0 C0 0 = 0 CO SO q (4.6) { 1 SO CO r x COCY COS So u y = CTSISOSC SY SISOSY+CTCO SICO v (4.7) z COSOCY +SOSY CSOSY SOCY COC w E These equations of motions are coupled by the state vector, x, and are dependent on the control vector, u. x = {u, v, w,, p, q,r, T, 0, (, x,y, z} (4.8) U = {OR1, R2, OE1, OE2, 8R1, 2 R2, 8E1, 8E2, 6c,Fprop 4.1.3 Small Disturbance Theory The small disturbance theory will be used for linearization of equations of motion. According to the theory the linearization will be carried about an operating point (reference flight condition), i.e., the equations thus derived will be valid for the torpedo operating at and near that value of vector x. The operating point is chosen to correspond to the trim condition in Equation 4.9. xo = {uo, vo,wo,, qo, o,Yo, o, o,xo,yo,zo} (4.9) = {uo,0, 0 0, 0, 0, 0, 0, 0, 0,} This corresponds to straight and level flight with constant velocity. As the torpedo may be traveling at other flight conditions, the linearization at those conditions would be carried out numerically, which will be explained in later sections. A value of uo is found numerically, that satisfies the equations of motion for a given value of xo. Then a disturbance of Ax is given to the equations of motion from the reference condition thus changing the flight conditions to xo + Ax. Several assumptions are made to carry out the linearization: * The disturbances from reference flight condition are small. Thus the terms involv ing higher order of disturbances A will be neglected. Furthermore first order terms involving A will be approximated as in Equation 4.10. Sin(A) = (A) (4.10) Cos(A) = 1 The propulsive forces and mass are assumed to be constant. Planing and immersion forces are neglected for this analysis. The linearization procedure is resolved for the force equation in bf direction. This equation relates the force, X, to the states. m(u + qw ru+gSO) = X (4.11) Using the flight condition from Equation 4.9 in Equation 4.11 gives the value of force at the reference trim condition. mgSOo = Xo (4.12) The perturbation equation, i.e., the equation with flight condition disturbed by Ax can be obtained by substituting the perturbed flight condition into Equation 4.11. m[ (uo + Au) + (qo + Aq) (wo + Aw) (ro + Ar) (uo + Au) (4.13) +gS(o+ A)]= Xo + AX Equations 4.12 and 4.13 can be combined to give the linearized form of Equation 4.11 for straight and level flight condition. m(Au+gAOCOo) = AX (4.14) Proceeding in a similar way all other equations of motion can be linearized. The lin earized equations for straight level flight are shown in Equation 4.15 to Equation 4.18. 4.1.3.1 Force Equations m(Au +gAOCOo) = X m(Av + uoAr gAcCOo) = AY (4.15) m(Au uoAq+gAOSOo) = AZ 4.1.3.2 Moment Equations IAp = AL IyAq = AM (4.16) IAr = AN 4.1.3.3 Orientation Equations COo A = Aq (4.17) A = Ap +TOoAr 4.1.3.4 Position Equations Ax = SOn,, AO + COoAu + SOoAw Ay =Av (4.18) Az = CO,,n, ,A SO ,A/ + COoAw 4.1.4 Stability and Control Derivatives The variations in total force and moment are often difficult to compute.These variations in forces can be calculated by chain rule for derivatives. As stated in Equation 4.8, these are functions of state variables x and control variables u. Thus for example AX can be written by chain rule as in Equation 4.19. AX =X,,A +XAv +X,, Ai +XpAp +XqAq + XAr + XAT + XAO + X~A( + XpropAFrop (4.19) + XO eR1 R+ +XOR 2 + XOEO E1 + xO E2E +XaR6R1 +XaR28R2 +Xa6E1 +X26E22 +X&6cc where the subscripted X represents its partial derivative with respect to its subscript. X, = (4.20) Iux=x0 Each of these subscripted variables that have a subscript of state variable are known as stability derivatives and the ones with a control variable as subscript are known as a control derivative. There can be as many stability derivatives as there are forces and state and control variables. Many of these are negligible, depending on the reference flight condition. These dependencies are known usually by experimentation or numerical calculations. For example, the force, X, is observed to depend mainly on a subset of the state and control variable. Thus only the stability derivatives that correspond to these variables have to be retained in Equation 4.19, when straight and level flight is considered. X = funct(u,w, 6E1,6E2, OE1, E2, 6c,Fprop) (4.21) The next problem is calculating numerical values of these derivatives. In most cases it is easy to calculate these numerically or using a symbolic manipulation software. For some reference points, it is possible to do the calculation manually. The calculation ofX,, will be done manually for straight and level flight. X, = (FR,x + FR2,x + FE,x + FE2,x + F,x + Fprop,x) (4.22) Expressions for each of the terms in Equation 4.22 have been derived in Chapter 3. For example, the expression for the force on cavitator is shown in Equation 4.23. cp>cCac CacSc Sac Dc (ac, yi) Fc,x = C8c 0 SSc SIc CIc 0 0 (4.23) cp3cSac ScacSpc Cac Lc (Cc, 7) In Equation 4.23, ac, Pc, and thus Lc and Dc are the only functions of u. Thus the partial derivatives with respect to u can be obtained. Fc, au ScLcS13j~7 CCkCI3 ~~Ps cc~ SacSPL^+C(aXCC C(a, U au u a S3 aC 0 S(XCCaPCL I+C( aa Sa^_ 0 SSc CPcC c CScSPc Sac 1 0 Sic C3c 0 0 Cc CcSocs SacSPc Cac It can be seen that ea, aP, L and Dc are terms required to be calculated. These can be calculated from equations 3.16 and 3.17, which are restated in Equations 4.25 to Equation 4.27. tan(ac) =  Uc ye tan(3c) = c+ V; /2+g+w; (4.25) (4.26) (4.27) The velocity components in Equation 4.27 can be found using Equation 3.2. CIc 0 SSc SICcs LaP + CP3S o SP3Sa^CL +CPCCQL Dc (ac,) C6c 0 + 0 Lc(c c,7) c S8c Dc( 0c, l 0 Lkc( ocY)7C TH 1 ^r (4.24) uc C8c 0 Vc = 0 1 We S8c 0 Uc U Vc = V + Wc W )IB )B SSc uc 0 vc C8c w B b1 b2 b3 p q r Xc Yc Zc Now the velocity components can be obtained for the reference flight condition that is stated in Equation 4.9. uc C 0cUO vc g= 0 WJ SScuo ( ); c (4.30) The variation of ac can be obtained by differentiating Equation 4.25 at reference flight condition. sec2 (c)dac ucdwc wcduc U2 c ucdwc wcduc dac = CScdw SScduc d U = Uo Uo (4.31) (4.32) at (xo, uo) Similarly the variation of Pc can be obtained by differentiating Equation 4.26 at reference flight condition. S (VYdvc vcdVc) Vc V (4.33) dvc = at (xo,uo) Now, using Equations 4.28 and 4.29, variation of velocity components of cavitator with respect to u can be obtained. CSc au (4.34) (4.28) (4.29) Finally, combining Equations 4.32 to 4.34, the variations of ac and 3P with respect to u can be obtained at reference flight condition. (xc Cs8c dw S6Sc uc au uo au uo au C CSSSc S8cC8c (4.35) uo uo =0 apc 1 v, au Uo au (4.36) =0 Clearly, two of the derivatives that are required to calculate stability derivatives have been obtained. It was previously stated that lift and drag are calculated using the coefficient of lift and drag. 1 L, = pV2S cclc 2 (4.37) 1 Dc = pV2Sccd 2 These forces can be differentiated by chain rule the derivative would be like in Equation 4.38. aLC 1 Vc a\cc c1 u = pSc 2Vclu + Vc u (4.38) au 2 [ du duc ou c The Equation 4.38 requires two derivatives. One of the derivatives is calculated in Equation 4.35. The other derivative can be calculated using Equation 4.27. Vc a \U2+V2+W21 au au [ c ] w uc L v, + VC a +WC3 uc + au a+ (4.39) uV V1 Thus the derivative of the lift and drag forces can be obtained. 3L c pSScVcl, (4.40) au aDc u= pScVccdc (4.41) du Thus all the derivatives required to calculate righthand side of Equation 4.24 have been cal culated. All the terms on righthand side of the Equation 4.20 can be calculated in a similar manner. It is tedious to calculate the derivatives analytically in such a way. The complexity increases for other flight conditions. For practical purposes these derivatives are calculated using symbolic manipulation software like 'Mathematica' or by using numerical methods. The numerical methods used to calculate the derivatives have been described in Appendix B. 4.2 State Space Representation Equations 4.15 to 4.18 are a complete set of linearized equations of motion for the torpedo. They can be represented in a more convenient form known as the State Space Form. The state space equations are a set of firstorder differential equations. x = Ax +Bu y = Cx +Du x E 9n,u E Wy E 9m (4.42) A E n x 9n,B E n x p C E Z x 9~,D 9Z x W Equation 4.42 is a generalized form of state space representation for any system. Each of the terms in the equations has a particular importance for describing the dynamics of the system. * State Variable x: The state variables for a system are a set of variables, when known at time to and along with input u, are sufficient to determine the state of the system at any time t > to. All the states of the system need not be measurable. * Input Variable u: This is the control surface deflections. * Output Variable y: The output variables are the measured parameters. These may or may not be same as the state variables. The output variables are usually considered to be measurable but sometimes they are estimated. The matrices A, B, C and D may either be constant or timevarying functions. In the case of the supercavitating torpedo, the state vector is of size 12 (n) and the control vector is of size 10 (p). x= [A Aw Aq Ag Av Ap Ar A AT Ax Ay Az] SJ (4.43) u = c 6E1 8E2 6R1 8R2 OEl 0E2 OR1 0R2 Aprop Some of these controls may not be needed for some maneuvers. From the linearized equa tions it can be observed that the state variables are not coupled by the states {(,x,y,z}. These four states can be removed from the analysis for now. The system becomes a 8 state system. These states can be further divided into longitudinal and lateraldirectional dynam ics. The variables Au,Aw, Aq,AO correspond to longitudinal dynamics, which also means the dynamics in b1b3 plane. The variables Av, Ap, Ar, Ac correspond to lateral dynamics, which is the dynamics in flb2 plane. Sometimes the lateral and longitudinal equations can be decoupled. Thus if the torpedo is making a pull climb/descent to a certain depth, usually its dynamics can be represented by longitudinal state variables. The plant matrix A can be divided into four parts. A F= Along Acoupl (4.44) A A 1 (4.44) Acoup2 Alatd When A is divided as in equation 4.44, where each element is a 4 x 4 matrix, Along would correspond to longitudinal dynamics and Alatd would correspond to lateral dynamics. Acoupl and Acoup2 are coupling matrices. It is required that the coupling matrices become negligi ble for the equations to be decoupled. If these parts are not negligible, the equations cannot be decoupled, and a 8 state model will be required to be considered. From linearized 44 equations, the four parts of the A matrix for the torpedo can be written as in Equation 4.45 to Equation 4.48. x. qo + x xu Zx qo +  Mu Mw 0 0 o, oP Lv Lp Nv Np(IyIx)qo T 1 0 1 Along Alatd Acoupl w wo+ uo + L T/ m MZq Cy Cco Lr C xp m vo + Mp(IxIz)r iy 0 Po+ Lw Ix Nw 0 gCGo+x gCQoSo +  ME ly 0 uo + gCOoCoo + r. (IzIy)0q Le Ix Ix Nr N kI@0 TTo qoC(oSO, ,, ,,,, vo +x, x o m m ZI gS0oCeo + 0 Mr(IxIz)po Me Iy Iy Sco Scoqo Ccoro I gSOoSQco + Lq(IzIy)ro Lo Nq(IyIx)PO No IS IS SIoTO0 Scoqo +C4oro Similarly B is a 8 x 10 matrix, whose elements are just the control derivatives according to their locations in the matrix. The first 4 rows correspond to longitudinal dynamics and the last 4 correspond to lateral dynamics. X8c XFpropX m m Blong . 0 *.. 0 (4.49) 4x10 ro+ Po + Mv 0 Lu ro + Ix N, 0 (4.45) (4.46) (4.47) (4.48) YSc YFprop, m m Blatd= : (4.50) 0 *** 0 S4x10 Now the complete state space representation for the torpedo can be written as in Equation 4.51 which gives two sets of equations. The first set is the longitudinal equations and the second set is the lateraldirectional equations. T1 Xlong = Au Aw Aq AO xlatd = Av Ap Ar AI] U= 8c 8E1 6E2 8R1 8R2 OEl 0E2 OR1 OR2 A o (4.51) SXlong Along Acoupl Xlong + Blong X10 Xlatd Acoup2 Alatd Xlatd Blatd 8x1 L 8x8 8x1 L 8x10 CHAPTER 5 CONTROL DESIGN SETUP This chapter deals with the control design for the torpedo described in previous chap ters. Various parameters associated with the control are restated in Table 5.1. Table 5.1 Control Parameters Longitudinal Lateral State u, w, q, 0 v, p, r, 0, Y Control 6c, ei, 8e2 8rl, 8r2 It should be noted that Y has been included in the states though it was observed in the state matrices that all other variables are independent of Y. It will be seen later that the inclusion of Y in the feedback states helps in improvement of performance. Also, for longitudinal control, two elevators and the cavitator are required. Similarly for lateral direction, the rudders should be enough for control. There are various control methods, like linear quadratic regulator (LQR) synthesis, u synthesis etc., which can be used to design a controller. Each of these methods has advan tages and disadvantages. LQR method gives a constant gain controller which is based on minimization of a quadratic performance index and considers the problem of robustness only in terms of gain and phase margins. psynthesis deals with robustness with respect to a wide variety of uncertainties to minimize an infinitynorm matrix but the resulting con troller can be of high order. Regardless of complexity and robustness, each design method presents difficulties. The LQR method was chosen for controller synthesis for the torpedo. This method was chosen because 1. It is easy to implement in a nonlinear simulation. 2. The resulting robustness for the LQR controller was acceptable. 3. It is straightforward to vary some design parameters to achieve performance goals. This chapter explains various problems associated with the control synthesis and the system model used for synthesis of the controller. 5.1 Openloop Performance Initially the equations of motion for the torpedo have been linearized for straight and level flight at a forward velocity of 75 ms1 x= {75, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 (5.1) It is found that the cavitator and two elevators are sufficient to maintain the above trim. It is also assumed that the value of propulsive force required to maintain this trim is fixed at the required value. U = {Sc, 5el, 5e2, 8rl, 5r2,Fprop} (5.2) = {0.0067, 0.0106, 0.0106,0,0, 4.0023e + 03} where the deflections are given in radians, and Fprop is in Newtons. As these parameters are obtained numerically, it may not be a perfect trim. The system may have some non zero accelerations, and consequently may tend to deviate from straight and level flight. To check this, the openloop dynamics are simulated at this condition without any feedback. Figure 5.1 shows the Simulink model used for openloop simulation. The control surface deflections are fixed at their trim values for this simulation. The closedloop system is obtained using the equations of motion that were derived in Chapter 3. The state derivatives are then integrated to obtain the state at the next instant. Figure 5.2 shows the openloop response for the torpedo at this trim condition. It can be seen that the openloop system is unstable. The simulation is carried out at the trim that is shown in Equation 5.1, i.e., all the values shown in these figures have to be zero. The system is seen to have oscillation about nonzero states. The state and control matrices obtained for this trim condition are shown in Equations 5.3 to 5.6. There are 5 control variables, {Sc, 8el, 8e2, 8rl, r2}. The matrices corresponding S NL Equation Integrator L J Torpedo at Foi  simt Clock Time Figure 5.1 Simulink Model for Open Loop Simulation Figure 5.1 Simulink Model for Open Loop Simulation vIvvvv~vvwovvvv'fvvv I0 20 40time(s40 4 35 3 12 5 S2 01 5 1  80 100 0 040 80 100 0 20 40 60 times) Figure 5.2 OpenLoop Response for Torpedo: w,p, q Control Trim state r Plotting 20 15 5" 5 80 100 to the lateral dynamics are of dimension 5 because the state Y is included in the lateral dynamics. Thus the lateral states are now {v,p,r,1, Y}. 4.5204 1.5417 1.3110 9.8100 0.2616 15.7648 78.5888 0 Along = (5.3) 0.0000 1.2077 3.5614 0 0 0 1.0000 0 32.3010 69.0608 69.0608 0 0 406.0942 303.3736 303.3736 0 0 Blong= (5.4) 158.4675 45.1531 45.1531 0 0 0 0 0 0 0 12.0422 0.0002 71.6011 9.8100 0 0.1813 54.2281 0.3004 0 0 Alatd= 1.1437 0.0025 1.2528 0 0 (5.5) 0 1.0000 0 0 0 0 0 1.0000 0 0 0 0 0 366.60511 366.60511 0 14297.086 14297.086 17276.994 17276.994 Blatd= 0 1.4129523 1.4129523 54.561629 54.561629 (5.6) 0 0 0 0 0 0 0 0 0 0 The longitudinal eigenvalues are {21.1414,4.5137,1.8262,0.0178} and the lat eral eigenvalues are {0, 54.2289,0.0002,6.6472 + 7.2683,6.6472 7.2683}. The eigenvalues clearly show that the system is unstable. It can also be seen that the longi tudinal dynamics have no oscillatory modes. Figure 5.3 shows the variation of eigenvalues Longitudinal Egienvalues for u=60:5:90 Lateral Egienvalues for u=60:5:90 1 1 0.5 5 0 20 Rea10 long) 0 10 0 60 20 0 20 Reral (oneal (at) (a) Longitudinal (b) Lateral Figure 5.3 Variation of Eigenvalues with Change in Velocity for the torpedo for different velocities. State values are fixed except for forward velocity, which is varied from 60 ms1 to 90 ms1. The figures show that the variation is small and most of the eigenvalues stay in negative half of complex plane. 5.2 ClosedLoop Problem As stated earlier the control problem can be subdivided into various problems. Each can be solved to get a final controller. The ultimate goal of the controller design is to achieve a desired trajectory or reach a particular point with minimization of some perfor mance criteria. The achievement of this goal requires addressing maneuvering, trimming, guidance and navigation. This thesis will consider the basic problem of maneuvering. So the problem is to be able to track a certain pitch and roll command while maintaining cer tain performance criteria. The performance criteria that the controller is required to meet are: * Track a step command in pitch or roll rate of size up to 30 deg/s. * Maintain an overshoot less than 15%. * Have a rise time of less than 0.6 sec. * Have a steady state error of less than 5%. Besides meeting the above mentioned performance criteria, the controller is also required to stabilize the closedloop system. Table 5.2 Control Constraints Cavitator Deflection 15 < 8 c < +15 Cavitator Rate 25/s <, 8 < +25/s fins 60 < 8f < +600 Fin Rate 25/s < 8f < +250/s Thrust 0 < Fprop < 30000N Various bounds are placed on the control surface deflections and rates. These bounds are listed in Table 5.2. The bounds are included in the nonlinear simulations and it is required that there is no saturation of deflection or the rates at least for the commands having the rate 30 deg/s. An actuator model is included in the system so as to constrain the rates of the control surface motion. This actuator is realized as a low pass filter, Ac = 80 Addition of this filter synthesizes a controller that has slower control deflections. Let qcomm(t) be a function of time, defining the desired pitch rate profile. The aim of the controller is to find a control law u(t) that yields an achieved pitch rate, qachi(t), to minimize the optimization problem stated in Equation 5.7. find u(t) that minimizes (t) = qachi(t) qcomm(t) subject to Umin < U < Umax (5.7) Umin < ll < lmax x = Ax + Bu where, umin and Umax are lower and upper bounds on control deflections. The quantities Umin and Umax are lower and upper bounds on control deflection rates. The problem becomes a disturbance rejection problem, when the commanded variable is 0 for all time. This is an optimization problem,where the state space equation is an equality constraint and the control surface bounds are inequality constraints. 5.3 Robustness of the Controller A control system that is designed to accommodate the uncertainties in a mathematical model is called a robust control system. Usually the response of a model does not accurately match the response of the true system. A robust control system should give the desired performance not only during the simulations of the model, but for the true system also. Various parameters can be introduced in the model to simulate uncertainties. Random noise can be added to output signal to simulate measurement errors, the gains in signals can be changed to model uncertainty in response etc. Gain and phase margins are generally used to predict the robustness of a control system. Physically, these are just the factors by which the feedback gain can be increased and yet have a stable real system. A formal definition of these can be given by using a frequency analysis for a feedback system. 5.3.1 Gain Margin Figure 5.4 shows a typical feedback system involving a plant, P, and a controller, K. The gain for the system in dotted region is known as the loop gain. The loop gain is important because it determines stability. Errors in the predicted loop gain could cause errors in predicted stability. The gain margin is the largest factor by which this gain can be increased and still have a stable system. Physically, it means if the response of the torpedo for a given elevator input is higher by a factor of the gain margin, the torpedo is still stable. The gain margin is usually expressed in decibel (db) units, and can be easily obtained from the Bode plots for the system. The gain margin is a measurement of the magnitude on the Bode plot, at the point where the phase is 180. 5.3.2 Phase Margin Gain is a valid robustness criteria when the system has real eigenvalues. But usually the eigenvalues have imaginary components and thus the phase is also a concern. Phase margin is the measure of the maximum possible phase lag before the system becomes unstable. From the Bode plot, phase margin is the phase when the magnitude of the gain is zero. Figure 5.4 Loop Gain 5.3.3 Uncertainty In Parameters Another factor that can determine the robustness of a controller is its response to errors in known parameters. As stated earlier, the coefficients of lift and drag are calculated from a CFD database. The accuracy of the model depends on accuracy of this data. Robustness of a controller can be assessed by introducing errors in the data and checking how it performs. The following variations have been introduced in the system to check for performance of the system with intrinsic uncertainties: * 120% error in C1 of Cavitator. * 120% error in Cd of Cavitator. * 120% error in C1 of all the Fins. * 120% error in Cd of all the Fins. 5.3.4 Controller Objective In terms of robustness, the controller is required to meet various performance objec tives. These objective can be summarized as: * The closedloop system should have a gain margin of at least 6 dB. * The closedloop system should have a phase margins of at least 45 deg. CHAPTER 6 LQR CONTROL 6.1 LQR Theory The tracking problem, like the one given in Equation 5.7, can be solved by using a combination of feedback and feedforward control [10]. The Linear Quadratic Regulator (LQR) problem is to find an optimal feedback matrix K such that the statefeedback law u = Kx minimizes the linear quadratic cost function shown in Equation 6.1. J(u)= (xTQx + uTRu + 2xTNu)dt (6.1) 0 The basic idea of LQR control is to bring the state of the system close to zero. A linear system can be represented in the state space form as in Equations 6.2 and 6.3. The matrices A and B are the state and control matrices. The variable x represents the state vector, y is the output vector and u is the input vector. S= Ax+Bu (6.2) y=x (6.3) The LQR controller is realized by a constant gain matrix K, such that the feedback u = Kx makes x go to zero. By a modification to this law, the LQR method can also be used for tracking. The state vector x is of size n. x= {X1,X2,", XnT (6.4) Let the tracking problem be for the state xl to track a step command rl. The idea is to make (xl rl) go to zero using a LQR controller. The new control law can be chosen as in Equation 6.17. x1 rl X2 u= K (6.5) S x Equation 6.2 can be rewritten by substituting the new control law. x = Ax + Bu x1 r1 X2 (6.6) = Ax BK \ x, For simplicity, assume that there is only one control, u (this is different from velocity u). The controller K is of size n x 1 and it can be expanded in its elements. K =[kl,k2,... ,k,] (6.7) Equation 6.6 can be rewritten by substituting the K in its expanded form. xi ri X2 k = Ax B[kl, k2, ,kn] (6.8) \ xn =Ax BKx +Bklrl =(A BK)x +Bkirl It should be noted that the command rl is a step command. The steadystate dynamics of the system can be obtained from Equation 6.8. (o) = (A BK)x() +Bkirl (6.9) The error dynamics can be obtained by subtraction Equation 6.9 from Equation 6.8. k(t) 2() = (A BK) (x(t) x()) (6.10) e = (ABK)e (6.11) where e = (x(t) x(o)). So, the tracking problem is cast as a regulator problem. The new state vector is the steadystate error e, which is made zero using the regulator. Figure 6.1 shows the block diagram for this closedloop system. It is required that the closedloop system has an integrator somewhere so as to make the steadystate error go to 0 [10]. That is, e has to go to zero rather than e so as to achieve good tracking. Figure 6.2 shows the new configuration of a system that has no integrator and thus an integrator has to been included during design. Thus, the integral of the actual error has to be made to go to zero so as to achieve a good tracking. e= (rixl) (6.12) The state space equation for the system with this modification can be written. S=Ax+Bu e= ri xi (6.13) = rl Cx where xl = Cx. It can be seen that the error equation is similar to state equation. Thus e can be considered as another state, .i.e, the system now has n + 1 states with state vector x = {xl,x2, ,Xn,, T. So a new formulation of the state space equation can be written, A 0 B 0 X= x+ u+ ri (6.14) C 0 0 1 x=Ax+Bu+Ir (6.15) which is similar to Equation 6.8. The error dynamics of this system represent the form of state space equations, for which a LQR controller can be derived. The LQR controller K, will be a constant matrix of size n + 1 as the system now is of size n + 1. 57 K= [kl,k2,.. ,kn lk,+I] (6.16) Then, the new control law can be written as in Equation 6.17. u = K x = [k, k2, kn I kn+1 S(6.17) = [ki,k2, ,kn]x + [kn+l ] =Kx +ki which is represented in Figure 6.2. r +0 x =Ax+Brk X y=x K Figure 6.1 Controller for Tracking when Plant has an Integrator. k1 =Ax+Bu +0 f+ y=x K C  Figure 6.2 Controller for Tracking when Plant has no Integrator. 6.2 Control Synthesis The torpedo system does not have an integrator in the system. A tracking controller can be obtained from LQR method by the process described in Section 6.1. The controller is obtained for a trim state of straight and level flight at 75ms 1. The linearized dynamics are first separated into the longitudinal and lateral dynamics as given in Table 5.1. The controls used in longitudinal direction are the cavitator and 2 elevators. The controls in lateral direction are the 2 rudders. It is observed that for straight and level flight the longitudinal and lateral dynamics are practically decoupled. Once the state and control matrices have been obtained, the main variables that the LQR controller depends on are the weighting matrices Q, R and N. In this case the cross coupling matrix N is chosen to be 0. The matrices Q and R penalize the cost function for higher state and control values respectively. A higher value in Q matrix would cause a better track ing. A larger R would constrain the control surface deflection. An optimum combination of the matrices is obtained iteratively, so as to get good tracking with achievable control deflections. The matrices for longitudinal pitch rate tracking are given in Equation 6.18. Qlong = diag([0, 0, 0, 0, 10]) (6.18) Rlong = diag([5,4]) The first four numbers in the Qiong correspond to weightings on the four longitudinal states. They are chosen to be 0. We do not want to restrict the states from changing. This is especially important for weightings on q and 0. A weighting on these variables would restrict them from changing. The last number, 10, is weighting on the error. This is chosen to be high to penalize the tracking error. A higher value of weighting would give a better tracking, but it is seen that it would require very high control rates. The first number in the weighting matrix Rlong, 5, corresponds to cavitator weighting and the number 4 is for elevator weighting. Elevator weighting is chosen to be smaller so as to encourage the controller to use more elevator than the cavitator. This gives a more stable performance. The control matrices obtained for the longitudinal dynamics are given in Equations 6.19 and 6.20. 1.1182 k1 = (6.19) 0.9681 0.0000 0.0040 0.1016 0.0000 1.4195 1.1184 K = (6.20) 0.0000 0.0042 0.0995 0.0000 1.3981 1.1308 Similar process is involved in the design of the lateral controller. Initially the lateral controller is designed with only four state feedback, and Y is neglected in the feedback. In this case it is found that the torpedo has high sidewash and deviates considerably from the original path, even when the pitch angle is 0. To avoid this, Y is included in the feedback states. It is also observed that a penalty on the yaw motion causes the controller to com mand a very fast control surface motion. Also, a continuous correction of control surface deflection is required to prevent the yaw motion entirely. Thus an optimum combination of the weighting matrices is obtained that would prevent a very high yaw motion but would still have slow control surface motion. Qlatd = diag([0, 0, 0, 0, 0,.1]) (6.21) Rlatd = diag([1000, 1000]) The first 5 numbers correspond to 5 states and the last number is weighting for the error. The Rlatd is of dimension 2 as only the rudders are included in the synthesis. The weighting on the rudders is high as it is observed that the roll rate is very sensitive to the rudder deflection. The control matrices obtained for the lateral dynamics are given in the Equations 6.22 and 6.23. 0.0071 k, = (6.22) 0.0071 =0.0005 0.1253 0.0132 0.0019 0.0000 (6.23) K=10.3 (6.23) 0.0005 0.1254 0.0132 0.0026 0.0000 The feedback matrix K for lateral dynamics is of size 2 x 5, which is shown in Equation 6.23. 6.3 Nominal Closedloop Model 6.3.1 Model Figure 6.3 shows the eigenvalues for the closedloop longitudinal and lateral systems. It can be seen that both systems are stable as all the eigenvalues are in the left half of the complex plane. Also, each of the dynamics has one eigenvalue at the origin, which is introduced due the integrator in the system. 15 8 6 10 4 5 2 5 E 2 5 4 10 6 6 10o0 800 600 4P(lon o00 0 200 1o00 800 600 4 d00 0 200 (a) Longitudinal (b) Lateral Figure 6.3 Eigenvalues for the Closedloop System 6.3.2 Linear Simulations The response of the closedloop linear system has been shown in Figures 6.4 to 6.8. The simulations for lateral and longitudinal systems have been carried out separately as the linear system is decoupled into lateral and longitudinal. Figure 6.4 shows the tracking obtained for a 15 deg/s pitch rate command. The com mand is achieved in 0.17s with an overshoot of 3.95% and with no steadystate error. Fig ures 6.5 and 6.6 show the control surface deflections and rates required to achieve the com mands. Though there are some quick deflections, the rates are still under the constraints. 5 Time(s) 10 15 0 5 times) 10 Figure 6.4 Pitch Command Tracking for Linear System : q, u 5 times) 10 15 0 5 times) 10 Figure 6.5 Pitch Command Tracking for Linear System : 8c, 8c Figure 6.7 shows the roll rate tracking obtained for a 15 deg/s roll rate command. The command is achieved in 0.53s with no overshoot and a 0 steadystate error. The variation of the yaw rate is also shown in the figure and it can be seen that the yaw motion is coupled with the roll. At the end of the roll doublet, the torpedo has a nonzero yaw angle thus changing the direction of motion. The control surface deflections required for the roll rate command are shown in Figure 6.8. The rudder deflection is small for reasons explained in the next chapter. 0.01 0.02 0.03 0 20 I/) ao) 1 10 10 1" 0 S5 times) 10 15 0 5 times) 10 Figure 6.6 Pitch Command Tracking for Linear System : 8el, 8el Commanded  Achieved 5 10 15 ) 5 10 Time(s) times) Figure 6.7 Roll Command Tracking for Linear System : p,r U) 0) 0) t0) g! e n . 5 times) 10 15 0 5 times) 10 Figure 6.8 Roll Command Tracking for Linear System : 6r1, 6r rf rL r 6.3.3 Gain and Phase Margins The LQR tracking system shown in Figure 6.2 is obviously more complex than the system shown in Figure 5.4. Thus, the loop gain can be defined in many ways in this case. The block diagram can be broken at different points so as to simplify it to the form shown in Figure 5.4. Figure 6.2 is redrawn in Figure 6.9 which shows the possible breakpoints for this system. For understanding, the output of plant P is divided into two parts, one is the achieved value of the commanded variable (ra) and the other is remaining states of the plant P (x). The break points are numbered 1 to 3. The system can be broken at each of these points to give a loop gain. These gains will be named outerloop, innerloop and allloop gains respectively. rl + 1 x Ax + Bu ra  K +v 2 Figure 6.9 Breakpoints for Calculating the LoopGain for a Tracking Controller Gain and Phase margins for each of the above possible break points have been calcu lated for both the longitudinal and lateral controllers. Table 6.1 lists the gain and phase margins for the torpedo with LQR controller that was obtained in previous sections. All margins are quite high and meet the desired conditions of 6dB for gain and 45 deg for phase margin. Figures 6.10 to 6.15 show the corresponding bode plots for the data given in Table 6.1. Also, the lateral controller is unable to stabilize the unstable spiral mode. Thus the closedloop system is inherently unstable due to this pole and would consequently have negative gain margin. Numerous simulations show that the affect of spiral mode is negligible, Table 6.1 Gain and Phase Margin with LQR Controller Longitudinal Gain Margin(db) Phase Margin (deg) 1 21.056(at 47.498 rad/s) 64.846(at 9.0625 rad/s) 2 327.87(at 0 rad/s) 77.118(at 25.925 rad/s) 3 o 57.606(at 20.845 rad/s) Lateral Gain Margin(db) Phase Margin (deg) 1 22.964(at 0 rad/s) 2 250.51 (at 0 rad/s) 3 50.36 (at 0 rad/s) i.e., the time to double for the instability is considerably larger than the maneuvering time of the torpedo. So, the closedloop system model is reduced by removing the spiral mode from the model. The gain and phase margins in Table 6.1 are for this reducedorder system and reflect the robustness of the dominant dynamics. Bode Diagram Gm = 21.056 dB (at 47.498 rad/sec), Pm = 64.846 deg (at 9.0625 rad/sec) 0 S50 1  100 90 a135 P180 225 27 01 /sc102 3 10 10 (rad/sec) 10 10 Figure 6.10 Gain and Phase Margin: Longitudinal Outerloop 6.4 Perturbed Closedloop Model A perturbed system model is formed by adding an error to the values of coefficients of lift and drag for the fins and cavitator. New values of trim deflection are obtained for the perturbed model and thus a new set of A and B matrices is obtained. Tables 6.4 to 6.9 Bode Diagram Gm = 327.87 dB (at 0 rad/sec), Pm = 77.118 deg (at 25.925 rad/sec) 50 0 50 45 0 45 90_ 135 1804 13 100 10 (rad/sec) 102 10 Figure 6.11 Gain and Phase Margin: Longitudinal Innerloop Bode Diagram Gm = Inf, Pm = 57.606 deg (at 20.845 rad/sec) 100 0 100 90 1135 180 100 (rad/sec) 102 Figure 6.12 Gain and Phase Margin: Longitudinal Allloop show the percentage variation of the elements of A and B matrices for a 20% change in coefficients of lift and drag of cavitator and fins. Few elements in the state and control matrices change. In most cases, the change in elements of A and B matrices is a linear function of the change in a coefficient. For example, in Table 6.4 there are 8 terms that show a variation due to a 20% variation in coefficient of lift of the cavitator. The term A(3,1) shows a large variation but its numerical value is negligible. The term A(3,2) shows a 34% variation but this term is also small compared to other terms. Remaining terms in the matrix show very small variation. Some terms in the B matrix show a 20% variation. Bode Diagram Gm = 22.964 dB (at 0 rad/sec), Pm = Inf 20 30 840 50 180 S135 90 100 10 (rad/sec) 102 103 Figure 6.13 Gain and Phase Margin: Lateral Outerloop Bode Diagram Gm = 250.51 dB (at 0 rad/sec), Pm = Inf 20 30 140 50 90 145 0 Q45  90 M o 10 101 (rad/sec) 102 103 Figure 6.14 Gain and Phase Margin: Lateral Innerloop Thus some terms in controllability matrix change considerably. This would mean that for an error in these coefficients, the response would show some difference in control surface deflection. As it is observed that the closedloop system has good gain and phase margins, this effect on B matrix should not be of much concern. 6.4.1 Model Figure 6.16 shows the eigenvalues for the perturbed closedloop longitudinal and lateral systems. An error of 20% is included in the value of coefficient of lift for the fins. It can Table 6.2 Percentage Variation in A Matrix due to 20% Variation in clc u w q O v p r c u 0.46 0.62 w 5.52 6.86 1.58 q 1.05e5 34.8 13.58 v p r ~( ~~ Ty Table 6.3 Percentage Variation in B Matrix due to 10% Variation in cc 8c 8el 8e2 861 8r2 U u w 20 q 20 0 v p r (D Ty Table 6.4 Percentage Variation in A Matrix due to 20% Variation in cdc u w q O v p r (Y u 10 6 7.6 w 1.54 0.36 q 658 7.8 3 v 2 20 0.4 p 2 15.6 r 10 0.8 8.4 ~( ~ Ty Table 6.5 Percentage Variation in B Matrix due to 20% Variation in cdc 8c 8el 8e2 8r1 6r2 u 20 w q 0.12 0 v p r ~(~ Ty Table 6.6 Percentage Variation in A Matrix due to 20% Variation in clfi, u w q O v p r ( u 1.22 0.64 w 14.4 10.0 0.9 q le5 60 3 v 16.2 1.2 p 19 36.0 r 25.4 0.94 10.2 Table 6.7 Percentage Variation in B Matrix due to 20% Variation in clfi, 8c 8el 6e2 6r1 6r2 u w 20 20 q 20 20 O v 20 20 p 20 20 20 20 r 20 20 ~(~ Ty Bode Diagram Gm = 50.36 dB (at 0 rad/sec), Pm = Inf 50 60 C 70 80 180 S135 10 10 (rad/sec) 10 10 Figure 6.15 Gain and Phase Margin: Lateral Allloop Table 6.8 Percentage Variation in A Matrix due to 20% Variation in cdfi, u w q O v p r ( Y u 9.4 24.0 12.4 w 1.34 0.12 q 68 2.6 0.4 v 20 1.76 2.3 0.13 p 2.3 1.12 0.62 r 2.74 18.26 1.12 be seen that the longitudinal dynamics show some perturbation in the damping while the lateral system relatively unchanged. 6.4.2 Linear Simulations Figures 6.17 to 6.19 show the response of the perturbed system for a 15 deg/s pitch rate doublet command. The perturbation to the system is a 20% error in the value of the coefficient of lift of the fins. It can be seen that the performance criteria are met even in the case of the perturbed system. It can also be observed that the overshoot is increased for the perturbed system. The performance is achieved at the cost of small perturbations in Table 6.9 Percentage Variation in B Matrix due to 20% Variation in cdfi, 8c 8el 8e2 8rl 8r2 u 20 20 w q 1.36e2 v p r 20 20 12.6e3 2.6e3 _(____ ______ 15 8 6 10 4 5 2 E 2 5 4 10 6 1100 800 600 400 200 0 1400 800 600 400 200 0 200 Real Real (a) Longitudinal (b) Lateral Figure 6.16 Eigenvalues for the Perturbed Closedloop System: 20% Error in clfi, other states of the system. As the control effectiveness of the control surfaces is changed, the amount of control surface deflection is also changed by a constant factor. Figures 6.20 to 6.21 show the response of the perturbed system for a 15 deg/s roll rate doublet command. The perturbation to the system is a 20% error in the value of the coefficient of lift of the fins. It can be seen that the performance criteria are met even in case of perturbed system. In this case also, it can be observed that there is a perturbation in other states.  No Uncertainty +20% in cl 20% in clf /  No Uncertainty V + +20% in cl ...... 20% in cl 5 times) 10 times) Figure 6.17 Pitch Command Tracking for Perturbed Linear System : q, u o0 05 40 5 10 15 0 5 10 times) times) Figure 6.18 Pitch Command Tracking for Perturbed Linear System : 8c, 8c 5 times) 10 times) 20 U) ) 10 0 1 0 10 S20 15 0  No Uncertainty +20% in cl ...... 20% in clf 5time(s) 10 Figure 6.19 Pitch Command Tracking for Perturbed Linear System: 6ei, 5e1 I 74 E "73.5 :3 5 Time(s) 10 Times) 72.5 S 720 15 0 1 5 1 S05 " 0 05 20 15 10 o 5 0 " o 5 <10  No Uncertainty _ +20% in clf 20% in clf Figure 6.20 Roll Command Tracking for Perturbed Linear System : p,r  No Uncertainty 002 ___ +20% In clf 0.1 ... 20% In clf 003 0.15 0 5 10 15 0 5 10 15 times) times) Figure 6.21 Roll Command Tracking for Perturbed Linear System : 8rl, 6r 6.4.3 Gain and Phase Margins Table 6.10 lists the gain and phase margins for the perturbed closedloop system. The perturbed system also has good gain and phase margins. Comparing the values with Table 6.1, it can be seen that there are small changes in the values except for the lateral allloop. The last value is increased to o showing an improvement for the perturbed system. From the analysis of the perturbed closedloop system it can be said that the linear model is robust to various uncertainties in the system. 73 Table 6.10 Gain and Phase Margin for Perturbed Closedloop System: 20% error in clfin Longitudinal Gain Margin(db) Phase Margin (deg) 1 21.193(at 49.599 rad/s) 68.981(at 8.7966 rad/s) 2 320.33(at 0 rad/s) 77.605(at 25.925 rad/s) 3 Io 60.305(at 22.552 rad/s) Lateral Gain Margin(db) Phase Margin (deg) 1 24.391(at 0 rad/s) o 2 278.23 (at 0 rad/s) C 3 (at 0 rad/s) Co CHAPTER 7 NONLINEAR SIMULATIONS The controller for longitudinal and lateral dynamics have been obtained separately. That is, the longitudinal controller is to achieve a required pitch rate and the lateral controller achieves a given roll rate. Once the controllers have been found for linear systems, they are employed with the nonlinear torpedo model and the performance is checked using numeri cal simulation. Figure 7.1 shows the complete nonlinear simulation model for the torpedo with the LQR controller. The nonlinearity in the model is provided by both the aerodynamics and control surface constraints. These constraints, given in Table 5.2, are not directly included in the linear model. So it is important to find their effects on the nonlinear simulation. It should be noted that there is a constraint on thrust, but thrust is assumed to be constant with time. The thrust constraint is used to find a trim value for various trajectories. 7.1 Nonlinear Simulations for Nominal System The simulations have been carried out for various commands with the nominal sys tem. The 'Nominal system' is the nonlinear torpedo model assuming that it is completely accurate. No uncertainty is included in the model. The simulations show good tracking response while meeting all the performance criteria. The response of the vehicle to a longitudinal command is simulated and shown in Fig ures 7.2 to 7.5. These figures show the response for a pitch rate doublet of 15 deg/s. The rise time for the pitch rate command of 15 deg/s is 0.18s and there is an overshoot of 11.53%. The steadystate error is .8%. The controller is able to command pitch rates as high as 30 deg/s. It is observed that the vehicle motion is confined to longitudinal plane Derivative r2dl Figure 7.1: Complete Nonlinear Simulation with LQR Controller x10 20 10 0 0 0 5 times) 10 15 0 5 times) 10 15 Figure 7.2 Pitch Command Tracking : p, q 1 20 0.5 10I 0) 0 00 0.5  20 10 30 0 5 times) 10 15 0 5 times) 10 15 Figure 7.3 Pitch Command Tracking : c, 8c only. This shows that the controller allows pure longitudinal motion to be uncoupled from the lateral motion. 2 20 15 1.5 10 5 Cn z 0 0.5 16 5 10 0 15 0 5 times) 10 15 0 5 times) 10 15 Figure 7.4 Pitch Command Tracking : e,1, e1 2 107 300 15 200 120 starting point 0100 800 600 5 10 15 0 0 4 times) 10o 20 0 x(m) Figure 7.5 Pitch Command Tracking : 8ir, {x,y,z} Trajectory The response of the vehicle to a lateral, roll rate, command is shown in Figures 7.6 to 7.10. A roll rate command of 15 deg/s is achieved in .52s with an overshoot of 0% and a steadystate error of 0.09%. The controller is able to command a roll rate motion of as high as 50 deg/s before a saturation of control surface rate is reached. It is observed that there is some longitudinal motion in this case. This longitudinal motion has been reduced by inclusion of Y in the feedback states to the controller. It can be seen that the rudder deflection for a roll rate command is small. This is expected as the terms corresponding to roll rate from rudder are an order of 3 times larger than the terms corresponding to pitch rate from elevators. It is assumed that the control surface deflection is achievable. 20 0.5 120 00 SI 0.5 0 i1.5 10 2 20 5 times) 10 15 2 5 times) 10 15 Figure 7.6 Roll Command Tracking: p, q Figures 7.11 to 7.15 show the response of the torpedo for a combined roll and pitch rate command, similar to a windup turn. In this case, the torpedo is given a 5 deg/s roll 5 times) 10 15 5 times) 10 Figure 7.7 Roll Command Tracking: c, 8c 2 V 2.5 1 31 5 times) 10 15 0 5 times) 10 15 Figure 7.8 Roll Command Tracking: 8el, 8el 0.5 5 time(s)10 15 times) Figure 7.9 Roll Command Tracking: 8,1, 6,1 0 to 0.4 0.3 0.2 0.01 0.02 0.03 * 1000 500 Starting Point 0 1000 1500 500 1000 500 y(m) 0 0 x(m) Figure 7.10 Roll Command Tracking: {x,y,z} Trajectory rate command from 2 to 12 seconds, 5 deg/s pitch rate command from 12 to 22 seconds, 5 deg/s pitch rate command from 22 to 32 seconds and then 5 deg/s roll rate command from 32 to 42 seconds. As the vehicle motion is a little different from the actual trim, it can be seen the vehicle has considerable sidewash. Despite this sidewash, the controllers give a good tracking performance. The rise times for roll and pitch commands are .5s and 0.22s respectively. The overshoot for the same are 0% and 0% respectively. 4 I 4 (n i(n 0 0, __2 '2 4 4 6o 10 20 30 40 50 0 10 20 30 40 50 times) times) Figure 7.11 Roll & Pitch Command Tracking: p, q 7.2 Nonlinear Simulations for Perturbed System The performance of the controllers is studied using the simulation with a perturbed system model. An error is assumed in the values of various coefficients and a correction factor is added. Response of the closedloop nonlinear system is not much affected by the variations in coefficients of lift and drag. It is observed that the controller commands the 2i 0    10 20 30 40 50 times) Figure 7.12 Roll & Pitch Command Tracking: 8c, 6c 6 4 (n 2 0 2 40 10 20 30 40 50 times) Figure 7.13 Roll & Pitch Command Tracking: e,1, e, 0.01 0.04 0 10 20 tme(30 40 50 0 times) 10 20time(s 0 40 50 times) Figure 7.14 Roll & Pitch Command Tracking: 8,r, 6 2000 1000 Starting Point 0 1000> 2000 1000 1500 1000 0 500 y(m) 1000 0 x(m) Figure 7.15 Roll & Pitch Command Tracking: {x,y, z} Tracking 757 02 7576 ____ 01   t0 755 756 / 014 . 753 E 302 75 2 .75. __ 20% In cl. 0. 4 1 __ 20% In cl 72 in 04 20% in ci 20% in cl20% in cl 74 051 074 1 0 5 10 15 20 time time Figure 7.16 Response for 20% Variation in clfi,: u,w system to a new trim state which is also a straight and level flight, with change in speed and control deflections. After that, the system follows a pitch or roll command as well as before. Figures 7.16 to 7.19 show the response for one such case. In this case a roll doublet is commanded to the system, and there is an error of 20% in the value of clfi,. It can be clearly seen that the vehicle has gone to another trim state and then it follows the command equally well. There is almost no change in the trajectory of the vehicle. The control surface deflections are similar with a constant offset. Such response has also been checked for other cases. The affect of error is similar in all cases. By above analysis it can be said that the LQR controller designed for the torpedo for pitch and roll rate tracking commands is fairly stable and can be expected to achieve good performance for the real torpedo.  No Uncertainty 20% in cl .....20% in cl 0) 5 1.5 2 2.5 0 3 20 0 Figure 7.17 Response for 20% Variation in clfin: p, q 15 20 Figure 7.18 Response for 20% Variation in clfin: 8c, 861 No Uncertainty 20% inclf 1000 20% in clf E 500 N 0 1000 500 y(m) 1000 500 n x(m) Figure 7.19 Response for 20% Variation in clfin: {x,y,z} Trajectory No Uncertainty 20% in clf 20% in clf 15 5 10 time CHAPTER 8 CONCLUSION 8.1 Summary A dynamical model for a supercavitating vehicle has been obtained. The vehicle is found to be openloop unstable, and a controller for stabilizing the pitch and roll rate motion has been obtained. The LQR controller shows good tracking performance for the vehicle and all the control objectives are met. The controller is also found to be robust to errors in cavity prediction and velocity changes. This robustness is further demonstrated by the fact that the closedloop system has high gain and phase margins. 8.2 Future Work The dynamical analysis of the vehicle has been derived with an assumption that the cavity shape is fixed. The openloop cavity dynamics need to be modeled and included in the synthesis. Robust control methodologies like psynthesis can be applied to obtain a more robust controller. A robust control design could include the uncertainties in the model during the control synthesis. The LQR controllers obtained are typically known as the 'innerloop' controllers. An outerloop controller is also needed for guidance and navigation. The idea is that the outer loop controller can be modeled for tracking the trajectory in space, based on the closedloop dynamics of the innerloop model. APPENDIX A REFERENCE FRAMES AND ROTATION MATRICES X3 Y3 YX2 X1 ,X2 Figure A. 1 Rotation of Frames Figure A.1 shows two frames X{xix2X3} and Y{yly2y3}. Y is rotated from X by an angle 0 about xaxis. Thus the basis vectors of frame Y can be written in terms of basis vectors of X frame. y2 = X2Cos(O) +X3sin(0) (A.1) y3 = x2sin(O) +X3cos(O) This relation can also be expressed in terms of matrices. y1 1 0 0 x1 Y2 = 0 cos(O) sin(O) x2 (A.2) y3 0 sin(O) cos(O) x3 This was a case of simple rotation. The matrix above in square brackets is known as the rotation matrix from X to Y and is represented as X_Y The rotation matrix can be generalized for a case when the two reference frames are arbitrarily oriented. 85 (yl,xi) (y, x2) (y1, x3) XJ = (y2,X1) (y2,X2) (y2,X3) (A.3) (Y3,X1) (Y3,X2) (Y3,X3) where (,) means the dot product of the two vectors. Thus, Y = XY*X (A.4) APPENDIX B NUMERICAL TECHNIQUES B.1 Interpolation of Force Data This section describes the numerical technique used to obtain the values of coefficients of lift and drag for cavitator and fins. B.1.1 Extrapolation Scheme For a better result, a quadratic interpolation/extrapolation scheme is used. Thus 3 points would be required to obtain an interpolated or extrapolated data value. Figure B.1 shows the shape functions used for one dimensional interpolation. Say, points {xi1, i,xxi+1} are used to find the value of a function f at point x. The value of f(x) would be given by a parameter a and the three shape function N1, N2 and N3. N = 1 + a + or (xi x ii) (Xi Xi1) (xi+ xi1 )2 (xi+ 1 1)2 a2 (B. 1) (xi xi) (xi xi+) (i xi)_ (xi ,Yi+) N3 (xi1 xi) (x+i1 xi+) j (xi xi+I) (xi xi+I) where, the value of a shape function can be obtained by finding the value of a. x xi1 (B.2) a=  (B.2) xi+1 xi aE[0,1] for x [xi1,xi+1] a < for x < xil a > 1 for x > xi+l Thus a E [0, 1] for interpolation and it is greater than 1 or less than 0 for extrapolation. N2 N\ / \ 3 0 t 0.5 (D I O / \ 0.5 I + 0 0.2 0.4 0.6 0.8 1 Alpha Figure B. Shape Function for One Dimensional Quadratic Scheme f(x) = N1 *f(xi) + N2 f(xi) +N3 *f(xi+l) (B.3) This method can be extended for 2D and 3D as in case for cavitator and fins respectively. B.1.2 Cavitator The coefficients of lift (clc) and drag (cdc) for the cavitator are functions of half angle (ha) of cavitator cone and angle of attack for cavitator (Xc). The CFD data [5] is avail able for combination of points given in Table B.1. Equation B.3 can be extended for 2D cavitator. 3 3 f(xc, ha) = I N(1, i)N(2,j)f(oc (i),ha(j)) (B.4) i=1j=1 ac (i) Value of ac at ith node ha (j) Value of ha at jth node N(1, i) ith Shape function for ac N(2,j) jth Shape function for ha B.1.3 Fins The coefficients of lift (clfin) and drag (cdfin) for the fins are functions of angle of attack (af) for fin, immersion (Sf) and sweepback angle (Of). The CFD data is available for combination of points given in Table B.2. Equation B.3 can be extended for 3D fin. Table B.1 Grid For Experimental Cavitator Data Half Angle (ha deg) I {15, 30, 45, 60, 75, 90} Angle of Attack (ac deg) {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15} Table B.2 Grid For Experimental Fin Data Immersion (Sf) {0.1,0.3,0.5,0.7,0.9} Sweepback (Of deg) {0,15,30,45,60,70} Angle of Attack (af deg) {0,1,2,3,4,5,6,7,8,9,10,12,15} 1. S>0&S<0.1&0>0 1. S > 0 & S < 0. 1 & 0 > 0 2. S>0.1 & S< 0.3 & > 30 3. S>0.3 &S< 0.5 & > 45 4. S> 0.5 &S< 0.7 & 0 > 60 Data not Available for <0.9 5. S > 0.7 & S < 0.9 & 0 > 70 6. S>0.9&S< 1 &0>0 7. a<0 8. a> 15 f(Sf,, O ), = I I IN(1, i)N(2, j)N(3,k)f(Sf(i), O(j),ayf(k)) i=lj=1 k=l S(i) Of(jW af (k) N(1, i) N(2,j) N(3,j) (B.5) Value of Sf at ith node Value of Of at jth node Value of af at kth node ith Shape function for Sf jth Shape function for f/ kth Shape function for af B.2 Numerical Linearization Numerical linearization can be done by the 'linmod' command in the Matlab Simulink toolbox. This can also be done by noting that, the terms in the A and B matrices are the derivatives of state rates with respect to states and controls. For example, suppose xo and uo represent the state and control values at trim. It should be noted that xo is a 9 x 1 (excluding the positions {x,y, z}) vector and uo is 5 x 1 (cavitator and four fins). xo = uo wo qo Oo vo Po ro (o To uo = {c0 e,10 e,20 ,10 ,r20 (B.6) The equations of motion are of the form as in equation B.7. x = f(x,u) (B.7) where the function f is a vector function having 9 outputs and there are 9 states. The code for nonlinear equation of motion, takes x and u as inputs and give the value of x as output. Let e define a very small change. Now the element A (i, j) can be calculated as in equation B.8. S f (xo + (j), uo)i f(xo, U)i A(i,j) = 1 < i,j < 9 (B.8) where, 7(j) means a matrix of size xo with all zeros except jth element, which is equal to e, and f. represents the ith element of vector f. An element B(i, j) also can be obtained in a similar way. f (xo, uo +e (j)) f/(xo, uo)i B(i,j) = 1 f where, 7(j) means a matrix of size uo with all zeros except jth element, which is equal to F. 