UFDC Home  myUFDC Home  Help 



Full Text  
ROBUST CONTROL OF SUPi i': / AITATING VEHICLES IN THE~ PRESii .:E OF DYINAMICl AND UNCERTAIN CAVITY By ANUKUL GOEL A DIIi'.i iT'; iON PRESENTED TO THIE GRADUATE SCHOOL OF THIE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUI~il i NTS FOR THE~ DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA~ ACKNO W1V i ii: 1' il ii `TS I would like to : : my sincere gratitude to my committee chairman, Dr. A~ndrew Kurdila, for hi s invaluable guidance throughout the course of this proj ect. I would also liket to thank him for giving me this opportunity to work on such a i : :::.:i T::. proj ect. I would also like to thank my committee cochair, Dr. Richard C. Lind, for his invaluable: guidance and inspiration throughout the i : : i I would like to thank Dr. Normnan FitzCoy, Dr. Brian Mann and Dr. Haniph Latchman for : on m~y committee. I would also like to showv my sincere .:~.:..: ..i: to D~r. John IDzielsk~i and Jamnmulamnadakia Anand Kapardi for their valuable contributions to this i 1 w .iIvould also like to .... m~y gratitude to all the members, past and present, of the Supercavitation r:l I would also like to thank the Office of Naval Research for the support of the research grant for the :.:; 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 ......i LIST OFTABLES . vi LIST OF FIGURES . ix ABSTRACT ......x CHAPTER 1 INTRODUCTION . 1 1.1 Cavitation . . 1 1.2 Aspects of Supercavitation . 3 1.3 Related Research . . 4 1.4 Overview of This Dissertation . 5 2 NONLINEAR EQUATIONS OF MOTION . 7 2.1 Configuration of the Vehicle . 7 2.1.1 Body . 7 2. 1.2 Cavitator . 8 2.1.3 Fins . 9 2. 1.4 Maneuvering . 9 2.2 Kinematic Equations of Motion . . 10 2.2.1 Orientation of the Torpedo . . 10 2.2.2 Orientation of the Cavitator . . 12 2.2.3 Orientation of Fins . . 13 2.2.4 Angle of Attack and Sideslip . . 16 2.2.5 Kinematic Equations . . 20 2.3 Dynamic Equations of Motion . . 22 2.3.1 Forces on Cavitator. . . 24 2.3.2 Forces on Fins . . 27 2.3.3 Gravitational Forces . . 30 2.3.4 Equations of Motion . . 31 2.3.4.1 Force equations . . 31 2.3.4.2 Moment equations . . 32 2.3.4.3 Orientation equations . . 32 2.3.4.4 Navigation equations . . 32 3 CAVITY AND PLANING DYNAMICS 3.1 MunzerReichardt Model . . 35 3.2 Logvinovich Cavity Model . . . 35 3.2.1 Logvinovich Theory of Independent Expansion . . 35 3.2.2 Cavity Centerline . . 36 3.3 Planing Model. . . 38 3.4 Planing Force Equations . . 39 3.5 Planing Kinematics . . 42 3.5.1 Calculation of Immersion . . 42 3.5.2 Method of Calculation of C O.... . . 42 4 LINEARIZATION . 47 4.1 Linearization . . 47 4.1.1 Need for Linearization . . 47 4. 1.2 Generic Form of Equations of Motion . . 48 4.1.3 Small Dishractuerbac Theory .. .. . . 48 4.1.3.1 Force equations . . 50 4.1.3.2 Moment equations . . 50 4.1.3.3 Orientation equations . . 50 4.1.3.4 Position equations . . 50 4. 1.4 Stability and Control Derivatives . . 50 4.2 State Space Representation . . 55 5 CONTROL DESIGN SETUP . 59 5.1 OpenLoop Performance for the Fixed Cavity Model . . 60 5.2 ClosedLoop Problem . . 63 5.3 Robustness of the Controller . . 64 5.3.1 Gain margin . . 65 5.3.2 Phase margin . . 65 5.3.3 Uncertainty in parameters . . 65 5.3.4 Controller objective. . . 66 5.3.5 pu analysis: 66 6 LQR CONTROL . 67 6.1 LQRTheory.. .. .. ................... .. 67 6.2 LQR Control for Fixed Cavity Model: 71 6.2.1 Control Synthesis. . . 71 6.2.2 Nominal Closedloop Model . . 73 6.2.2.1 Model . . 73 6.2.2.2 Simulations . . 73 6.2.2.3 Gain and phase margins . . 76 6.2.3 Perturbed Closedloop Model . . 77 6.2.3.1 6.2.3.2 6.2.3.3 Model Simulations Gain and phase margins 7 p/lHo SYNTHESIS CONTROL 7.1 Uncertainty 7.2 Synthesis Model. . 7.3 Control Objective and Constraints 7.4 p Controller for Fixed Cavity Model 7.4.1 Longitudinal controller . 7.4.2 Lateral controller . 8 HOMING CONTROL 8.1 Homing Control using Proportional Navigation . 8.2 Constant Missile and Target Velocity . . 8.3 Constant Velocity Target and Accelerating Missile . 8.4 Accelerating Target and Missile .... 8.5 Yaw Controller .... ..... .. 8.5.1 Yaw control using the LQR controllers .. 8.5.2 Yaw control using the p/~Ho controller . 9 CONCLUSION ... 9.1 Summary .... 9.2 Future Work . . . . 99 100 . . 101 102 104 . 104 . 105 107 107 107 APPENDIX A REFERENCE FRAMES AND ROTATION MATRICES .. B NUMERICAL TECHNIQUES ... B.1 Interpolation of Force Data . . B.1.1 Extrapolation scheme . . B.1.2 Cavitator .... B.1.3 Fins .. .. . B.2 Numerical Linearization .... C TESTBED DATA ... C.1 Fin Parameters for the Test Bed and FinData ..... C.2 Cavitator Parameters for the Test Bed and CavitatorData REFERENCES ... BIOGRAPHICAL SKETCH ... .109 111 ... .111 ... .111 .... .. .12 .........13 .... .. .13 115 .... .. .15 . .. .. .. .16 119 S 122 LIST OF TABLES Table page 5.1 Control Parameters . . 59 5.2 Control Constraints . . 64 6.1 Gain and Phase Margin with LQR Controller . . 77 6.2 Percentage Variation in A Matrix due to 20% Variation in clc . . 78 6.3 Percentage Variation in B Matrix due to 10% Variation in clc . . 79 6.4 Gain and Phase Margin for Perturbed Closedloop System: 20% error in clfi, 81 7.1 Gain and Phase Margins for the Longitudinal Controller . . 93 7.2 Gain and Phase Margins for the Lateral Controller . . 96 B.1 Grid For Experimental Cavitator Data . . 112 B.2 Grid For Experimental Fin Data . . 113 C.1 The Fin Variables Range Set. . . 116 LIST OF FIGURES Figure page Supercavitating Vehicle. Supercavitating Vehicle. Vehicle Body is Divided Into 4 Sections. Cavitator and Fins . BodyFixed and Inertial Frames Principal Planes of Symmetry for the Torpedo. Euler Angles of Rotation Cavitator Reference Frame Rudder and Fin Reference Frames Rudder 1 Fin Reference Frames Rudder 2 Fin Reference Frames Elevator 1 Fin Reference Frames. Elevator 2 Fin Reference Frames. Angle of Attack (u) and Sideslip (P) Cavitator: (a) Angle of Attack and Sideslip and I Fin Geometry Vehicle Undergoing TailSlap. CrossSection of Cavity. Planing Section The Angle . Velocity of Transom Relative to Frame B . Calculation of a with Assumption that the Cavit 1.1 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 2.15 3.1 3.2 3.3 3.4 3.5 3.6 (b) Hydrodynamic Forces y Surface is Planar . 5.1 Simulink Model for Open Loop Simulation 5.2 5.3 5.4 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 7.1 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.10 7.11 OpenLoop Response for Torpedo: w, p, q Variation of Eigenvalues with Change in Velocity Loop Gain. Controller for Tracking when Plant has an Integrator. Controller for Tracking when Plant has no Integrator. Eigenvalues for the ClosedLoop System . Pitch Command Tracking : q. Pitch Command Tracking : Sc. 6c Pitch Command Tracking : Se, Bel  Roll Command Tracking: p . Roll Command Tracking: Sc. c  Roll Command Tracking: Se, Bel . Breakpoints for Calculating the LoopGain for a Tracking Controller. Eigenvalues for the Perturbed ClosedLoop System: 20% Error in clfi, Response for 20% Variation in clfi,: p, q . Response for 20% Variation in clfi,: 6c,6el Calculation of Uncertainty Linear Fractional Representation. Synthesis Model for pu Controller. Pitch Angle Tracking for the Nominal Plant Cavitator and Elevator Deflections Frequency Response of the Pitch Angle Controller Comparison of Response for the Nominal and Perturbed Plants. Roll Angle Tracking for the Nominal Plant . Cavitator and Elevator Deflections Frequency Response of the Roll Angle Controller. Comparison of Response for the Nominal and Perturbed Plants. 8.1 Collision Triangle . . 99 8.2 Various Possibilities of Sequence of Commands to Obtain Yaw Rate Control 105 8.3 Sequence of Commands to Obtain Yaw Rate Control . . 106 8.4 Yaw Angle Response of the Vehicle for the Commands Shown in Figure 8.3 106 A.1 Rotation of Frames . . 109 B.1 Shape Function for One Dimensional Quadratic Scheme . . 112 C.1 Hydrodynamic Forces Acting at the Center of Pressure of Fin. . 116 C.2 Fin Cl and Cd Test Bed Data: Sweep Angle = 0. A) Fin Cl; B) Fin Cd .. 117 C.3 Fin Cl and Cd Test Bed Data: Sweep Angle = 15. A) Fin Cl; B) Fin Cd 117 C.4 Fin Cl and Cd Test Bed Data: Sweep Angle = 30. A) Fin Cl; B) Fin Cd 117 C.5 Fin Cl and Cd Test Bed Data: Sweep Angle = 45. A) Fin Cl; B) Fin Cd 117 C.6 Fin Cl and Cd Test Bed Data: Sweep Angle = 60. A) Fin Cl; B) Fin Cd 118 C.7 Fin Cl and Cd Test Bed Data: Sweep Angle = 70. A) Fin Cl; B) Fin Cd 118 Abstract of D~issertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of ii ::1 i: .: i i:y ROBUST CONTROL OF SUPli iAV ITATING VEHICLES IN THE~ PRESii .:E OF DYINAMICl AND UNCERTAIN CAVITY B~y Anukul Goel August I Chair: Andrew J. Kurdila Cochair: Richard C. Lind ..I~~~ ..I=1.0...i Mechanical and Aerospace Engineering Ulnderwvater travel is greatly limited by the speed that can be attained by the vehicles. L : :II :i, the maximnum speed achieved by underwater vehicles is about 40 mn/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 allowi~s it to attain high speeds underwater. ":: : : l.:IT :: is a phe nomenon which is observed in water. As the relative velocity of water with ** i to the vehicle increases, the pressure decreases and .::1 : :: ::11 it ; :I to form water vapor. '::i cavitationn 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 :i :7.1: torpedo. Nonlinear; .:::l : of motion are derived in i. i : i The latter part of the work deals with finding innerloop controllers to maneuver the torpedo. A controller is obtained by ILQR synthesis for pitch and roll rate control. Robust controllers are obtained by i:. ;! synthesis for tracking pitch and roll angle commands. 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 r::IThe innerloop controllers are used for guidance and navigation of the vehicle. It is observed that the pitch and roll angle controllers can be used for yawv rate control of the vehicle. This is applied to the horning guidance problem. A simplified case is solved for horning guidance. 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 corners of a body where the speed can reach high magnitudes. A classic example for cavitation is at the tip of propellers. Since the propeller 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. 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. The cavity behavior is generally discussed in terms of a dimensionless coefficient called the cavitation number o. This is calculated from Equation 1.1. P~ Pc (1/2)pwYo V3. where p, is the density of fluid and V. is the speed of the fluid relative to the body. The pressure in the cavity pc may be that of vapor pv or vapor plus air. This number relates the difference in the pressure between the cavity and its surroundings to the dynamic pressure (q, = (1/2)pwV3~) of the freestream. When cavitation occurs, the cavitation number will equal the negative of the lowest value of the coefficient of pressure. This value is referred to the critical cavitation number. This is the value above which cavitation starts to occur in the fluid. Supercavitation is a higher stage of cavitation when the body is entirely surrounded by the cavitation bubble. 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. Figure 1.1 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.1 Supercavitating Vehicle [1] 1.2 Aspects of Supercavitation Supercavitation is an extreme version of cavitation in which a single bubble is formed which envelops the moving object completely. The cavity takes on the shape necessary to conserve the constant pressure condition on its boundary. The shape is determined by factors like the body creating it, the cavity pressure and the force of gravity. With slen der axisymmetric bodies, supercavities take shape of extended ellipsoids beginning at the forebody and trailing behind, with the length dependent on the speed of the body [2]. Only small areas at the nose and on the afterbody remain in contact with the liquid. Expressions for maximum diameter and length of the cavity as functions of cavitator geometry, drag coefficient and cavitation number are available. The cavitation number and drag coefficient then determine the cavity geometry. The supercavity may consist com pletely of vapor (natural cavity), or have partial noncondensible gas composition (ven tilated cavity), which is injected from the nose. The determination of cavitation number and hence cavity structure in ventilated cavitation is more complicated than for a natural supercavity [3]. A high speed supercavitating projectile, while moving in the forward direction, may rotate inside the cavity. This rotation leads to a series of impacts between the proj ectile tail and the cavity wall, called the tailslap phenomenon. The impacts affect the trajectory as well as the stability of motion of the projectile. Despite the impacts with the cavity wall, the projectile, in certain cases, nearly follows a straight line path. May [4] has a detailed analysis from experimental data for traj ectories taken by unpowered proj ectiles upon water entry. Also the patterns of tailslap for such proj ectiles have been discussed. The wetted forces determine the vehicle speed and vehicle stability. Because the wet ted contact area is small, unwanted changes in the wetted contact area arising from local cavity breakdown may result in large destabilizing forces and moments for a supercavi tating vehicle. The wetted forces must be specified and understood to determine vehicle maneuverability. 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 in the later chapters). This information is useful to calculate the forces on fins of the torpedo. In late 90's a lot of research was done on the dynamics of supercavitating vehicles. Work shown in Kulkarni 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 proj ectiles rotate inside the cavity. This rotation leads to impacts between the tail of the proj ectile and the cavity wall. The frequency of the impact increases with time. These proj ectiles 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 Dissertation This work aims at formulating a control design for a supercavitating torpedo. Equations of motion for the torpedo are derived. Linear control methodologies, LQR and pusynthesis, have been applied to obtain various robust innerloop controllers for the torpedo. Chapter 2 briefly describes the configuration of the supercavitating torpedo. The loca tion and dimension of the control surfaces have been discussed. A detailed derivation of the equations of motion for the torpedo has been carried out in the later sections of the Chapter 2. Various reference frames have been used to obtain the kinematic equations of the tor pedo. Dynamic equations are derived using Newton's Laws. Various forces experienced by different regions of the torpedo have been explained. Chapter 3 describes the two cavity models that have been used in simulations. The cav ity contour equations and the cavity kinematics are presented. The hydrodynamic planing force acts on the hull of the vehicle and generally destabilizes the motion of the vehicle. The planing force equations and the planing kinematics are presented in the in the second half of this chapter. Chapter 4 describes linearization of the equations of motion using small disturbance theory. It is observed that the linearization, even for a simple trim, straightlevel flight, can be very complicated. Thus, numerical methods are used for this purpose. 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. Various robust ness analysis methods show that the controller is robust to numerical errors in coefficients of lift and drag. Chapter 7 describes robust control design for the torpedo using pu synthesis. Controllers for pitch and roll angle have been obtained, that are robust to various parametric and mul tiplicative uncertainties. Controllers for pitch and roll rate are also obtained that are robust to uncertainties in the cavitator. One of the missions of a missile is target homing. Design of homing controller for a supercavitating torpedo is discussed in the Chapter 8. The summary of the current work is presented in the concluding chapter of the disser tation, Chapter 9. A discussion of the future work is discussed in the next section of the thesis. Appendices A and B describe some of the notations and methods used for derivation and simulation of the equations of motion. A description of the CFD database [5] used to obtain numerical values for the code has been presented in the Appendix C. CHAPTER 2 NONLINEAR EQUATIONS OF MOTION 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 have 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 cavitator 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. Figure 2.1 Supercavitating Vehicle [8] 2.1 Configuration of the Vehicle 2.1.1 Body For simulation purposes the body of the torpedo has been designed with 4 separate sections. The first is cylindrical in shape with a very small radius compared to the main body. The cavitator pivots at the front of this section and forms the nose of the vehicle. The second is a section of a cone with a much larger radius than the cavitator, the third a Figure 2.2 Vehicle Body is Divided Into 4 Sections. cylindrical part making up most of the length of the vehicle. The last portion of the vehicle is the conical section accounting for the nozzle. These geometric components are depicted in Figure 2.2. 2.1.2 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.2. The main parameter that defines 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 Fine [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 are 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. An effort has also been made to design a 2DOF cavitator. 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.1.3 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.3 are horizontal (placed parallel to the axis of rotation of cavitator), called elevators and 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 in later sections of this chapter. Figure 2.3 Cavitator and Fins 2.1.4 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 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. 2.2 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. 2.2.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 (edi, e^2, e^3). The earthfixed reference frame is shown in Figure 2.4. The vector ^3 points in the downward direction, i.e., the direction of the gravity. The vectors e01 and e^2 are placed in the horizontal plane such that the basis vectors form a righthanded coordinate system. As shown in the figure, el points to the right and e^2 points outside the plane of the paper. A bodyfixed frame B is defined by the basis vectors (jby, j2, 3 ~) so as to simplify the derivation of the equations of motion. The frame B is centered at O, the center of gravity of the torpedo, and moves with the torpedo. The reference frame B is shown in the Figure 2.4. It can be seen in Figure 2.5 that the torpedo has two planes of symmetry namely fil2 and fil3. The plane $183 iS called the longitudinal plane and plane $182, the lateral plane. Figure 2.5 Principal Planes of Symmetry for the Torpedo Figure 2.6 Euler Angles of Rotation i: ex e3s Figure 2.4 BodyFixed and Inertial Frames 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 (1, 2, 3) and (fi1, 2, 3) 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 ^3 through a yaw angle, Y, to obtain the frame (1, 2,23). 2. Rotate (1, 2, 3) about 2 through a pitch angle, 8, to obtain the frame (fi, f2, 3) 3. Rotate (ftl,f2, 3) about ft through a roll angle, 4, to obtain the frame B. Figure 2.6 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 2.1i. CO 0 SO 0 10  SO 0 CO CGSY SAS8SY +t CYCO CAS8SY CY SO S 0 YSf ~Y CY 0 SG SOC8 CCOC8 1 0 0 0 A O 0 SO CA COCY CYSSOS CASY CASOCY +t SOSY E _B e82 e^3 (2.1) 2.2.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 cavitator fixed axes with respect to the body fixed axes is shown in Figure 2.7. The deflection of the cavitator is given by an angle, Sc, 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 St. 1 e3 dI e^3 b3 Figure 2.7 Cavitator Reference Frame From Figure 2.7, the rotation matrix from B to cavitator fixed frame C can be written as in Equation 2.2. di C6c 0 S6c bi(22 do = 0 1 0 b22 cl8c3 S6c 0 C6c b 2.2.3 Orientation of Fins Figure 2.8 shows the orientation of the finfixed reference frames. For convenience, all the fin frames are indicated by basis vectors (fi,}2}). In text they will be represented as (fi,f2,f3) ps, where subscript fin corresponds to a particular fin. All the fins have 2 DOFs, namely the sweepback angle (8fi,) and the fin rotation (6ip,). These can be defined as *Sweepback angle (6 ps,): 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 f3 direction for all fins. Sweepback angle determines the amount of fin that is enveloped in the cavity. Ruddei 1) Ruddel 2) FRONT VIEW Elevator 2 ) TOP VIEW Figure 2.8 Rudder and Fin Reference Frames Fin Rotation (6r;,): 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 (gl,Ag2,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 F piz can be written as in Equation 2.3. gzfi f It. t)R1 ~3 A Figure 2.9 Rudder 1 Fin Reference Frames (2.3) 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 2.9 shows the sweepback and fin rotation for rudder 1. The matrices for transformation from B to R1 can be derived as in Equation 2.4 and Equation 2.5. ~1 2 ceR1 S6eR1 0 C6R1 S6R1 0 S6eR1 ceR1 0 (2.4) (2.5) CiR1 0 SSR1 S6R1 ceR1 0 *Rudder 2 Figure 2.10 shows the sweepback and fin rotation for rudder 2. The matrices for transformation from B to R2 can be derived as in Equation 2.6 and Equation 2.7. (2.6) (2.7) ~1 2 eR 3se  R2 SeR2 ceR2 0 ftCS~fin 0 S~fin g1 1 0 g2 0 CS fin g3 fi #2 I 1 3 R1 }2 I fl 73 R1 1 2 $3 SSR1 0 C6R 1 _ g2 =S6R2 0 C6R2 R20SR2 g3 R2 1 0 ft C8R2 0 SSR2 )3 R2 SSR2 0 C8R2 _ 1 2 $3 Figure 2. 10 Rudder 2 Fin Reference Frames Elevator 1 Figure 2.11 shows the sweepback and fin rotation for Elevator 1. The matrices for transformation from B to El can be derived as in Equation 2.8 and Equation 2.9. g2 =S6eEl C6El 0 E (2.8) g3 El 0 0 1 3 CEl 0 SSE1 CE S6l0 b SS6eEl C8El 0 (2.9) 3~~~ El ce S8E 0sE CSE 0 0 1 Elevator 2 Figure 2.12 shows the sweepback and fin rotation for Elevator 2. The matrices for transformation from B to E2 can be derived as in Equation 2.10 and Equation 2. 11. g2 =S6E2 C6E2 0 8, (2.10) g3 E2 0 0 1 8 I ft CSE2 0 SSE2 CE SE2 0 b }2 0 1 0 SOE2 C8E2 O Lo(2.11) 73 E2 S8E2 0 CSE2 1 8 Equations 2.2 to 2.11 will be used in later sections to transform the forces on fins and cavitator to the bodyfixed frame. 2.2.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 43,Z gR1z1 I~ II eE2 ~E2 22 Figure 2. 12 Elevator 2 Fin Reference Frames body can be found by knowing the rotation matrices. Thus, V = ubl t +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 2. 13 shows these parameters and their geometric interpretation. Angle of attack, a, can be defined as the angle between the projection of velocity V onto $2 3 plane and the $1 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 (fi, f2, 3) can be described by rotation of the B frame by an angle a. Thus the rotation Figure 2. 13 Angle of Attack (u) and Sideslip (P) matrix from Fbody to B can be written. by Ca 0 Sa b3 Sa 0 Ca f3 The angle of sideslip, P, is defined as the angle between the actual velocity V and the projection of VT onto babJ3 plane. Again, a frame Gbody can be defined by rotation of Fbody by an angle equal to P in negative f3 direction, thus giving a rotation matrix as in Equation 2. 13. 1l Cp Sp 0 f 2 = s cp o (2.13) Velocity of CG of the torpedo in the Gbody frame can be written as VAl, 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 2.14. by CpCa CuSp Sa g1 by CpSa SuSp Ca g3 body (.4 Using Equation 2.14 V can be rewritten as in Equation 2.15. V = Vit = VC~ubiV~pb+VCBub3(2.15) where V2 = V2 = 2 t 2 2. From Figure 2. 13, relations between the velocity compo nents and the angles of attack and sideslip can be derived. tan a = (2.16) v sin 0 = (2.17) Though the matrix GbodyB in Equation 2.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 GfinB 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 afin and Pfin is a two step process: 1. Obtain the velocity of center of pressure of fin. VCl'body = Veg +tE 0B X TEcgCP (2.18) where VCPbodv is velocity of CP of fin in B frame, Veg is the velocity of CG of the torpedo in E frame, E B is angular velocity of B in E, and regCP is position vector from CG to CPfin. Equation 2. 18 can be rewritten as in Equation 2. 19. fufin u by b2 b3 vfin v + p r (2.19) Wfin B c Xfin Yfin Zfin where regCP = Xfinby + Yfninb2 + Zfinb)3 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 2.3 to 2. 11. MR1 ),[R1 0 SSR CR1 0 S6~R1 VR1 vR1 0 R1 0 S6R CR1 0 C6R1 UR1 (.0 wR1 SSR1 0 C8R1 0 1 0 wR1 SuR2 C6R2 0 SSR2 COR2 0 S6eR2 UR2 PR20 0 SR2 CR2 PR (2.21) al :R2SSR2 0 C8R2 0 1 0 O vEl ),[0 1 0 S6El C8El 0 vEl2.2 wEl S8El 0 CSE1 0 0 1 wEl ( E2E2 2h 00 SSE2CE2 S6~E2 0 uE2 vE2 0 12 0 S6E2 C6E lE2 0 vE2 (.3 wEE2 2S8E2 0 CSE2 __ 0 0 1 wE2 B The left hand terms in Equations 2.20 to 2.23 give the velocity components at the CP for corresponding fins, in that fin frame. These can be used in Equations 2.16 and 2.17 to find the angle of attack and sideslip for a particular fin. 2.2.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. E B API 1t qi2 t r3 (2.24) 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 2.24 can be written as in Equation 2.25. E B = e3 t X2 1t~ (2.25) The rotation matrices from Equation 2.1 can be substituted into Equation 2.25 to obtain Equation 2.26. Eo B= (@ SOYG )bi+( C8SQ+8t CO)&2+( CCOC~8SA)&3 (2.26) Equations 2.24 and 2.26 can be equated to obtain Equation 2.27. p SO 0 1 Pr C8CO S~O0 Equation 2.27 can be rewritten as in Equation 2.28. p = 0 SOY q = YC8SQ +t OCO (2.28) r = YC8CO 8SO To apply Newton's Laws, acceleration of the CG is required. The values of pi, 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. d. d.. (v)l (v) +' oB xv (2.29) dt \,I dt \/B where, sub script I denotes Newtonian (inertial) frame, and B is the bodyfixed frame. lo3 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. u, 8, 82 E~ ~ CM = (2.30) = + ur pw (.1 wi +t pv uq Similarly, the rotational acceleration will be required in the frame E. This can be written analogous to Equation 2.30. (2.32) Pr 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. yx =( wv (2.33) Equation 2.1 can be substituted in Equation 2.33 to obtain Equation 2.34. x COCY C89P SO y = CYSSOS8 CSY SAS8SY +CYCO SCOC8 v (.4 zx CASOCY+SOSY CASOSY SOCY COC8 w 2.3 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 2.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 FP=P (2.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 2.3 1 in Equation 2.35, Newton's Laws for the torpedo can be rewritten as in Equation 2.36. m vt ur pw =F(.6 w +t pv uq ;I ut q v =F 2.3b Newton's laws can be extended to rotation. Equation 2.37 are the Newton's Laws for rotational motion. M c= II (2.37) = Icnnu+tE 0B x H where H (= IcMnE B) is the angular momentum, IcMn 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 2.32, the Newton's Law for rotation can be written as in Equation 2.38. li0 0 pi by b2 b3 0 I2 0 q + p r = M (2.38) 0 I3 7 P 293 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 extrap olated to calculate the numerical values for the forces. When the body of the vehicle other than the fins or the cavitator comes into contact with the cavity wall, the body experiences forces called as the planing forces during contact with the cavity wall. Details of these kind of hydrodynamic forces will be dealt with in Chapter 3. FHydrodynantic = FR1 tFR2 tFEl FE2 ec Fplaning (2.39) MHydrodynantic = MR1 tMR2 tME1 tME2 tMc Mplaning (2.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 byl 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. Immersion: This force results from the vehicle being submerged underwater. Since the center of buoyancy is assumed to be coincident with the center of gravity this force contributes only to the external force term. The total forces and moments are expressed in terms of these components. F = FHydrodynantic t Gray FProp tFlninersion (2.41) M = MHydrodynantic (2.42) 2.3.1 Forces on Cavitator Figure 2.14 shows the forces acting on the cavitator. Coefficient of lift (clc) and drag (cdc) for the cavitator are functions of angle of attack, occ, and halfangle, yl, of the cavi 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 cle and cde are determined using CFD and tabulated in [5]. These values have been extrapolated to calculate lift and drag. Le =2pif Secl,(ae,7 y) (2.43) De =~2plSecd (ae, y ) (2.44) where Sc is the crosssectional area of the cavitator base. Directions of lift (Le) and drag (Dc) are as shown in Figure 2.14(b). These can be transformed to the body axes using 2.2 and 2. 14 for the cavitator. b1, = C_B(sc) x Gcav_C(ue,,i) x gl82, (2.45 Thus the forces on cavitator, in body frame, can be written. Faxv = (2.46) Le (Rc,7 4) SGear C6c 0 S6c cpc~ac CcaSpc SW Dc(c = 1 0 Spc Ocp 0 0t S6c 0 C6c OcpSue SucSic Cac Le(Rc,7 4) where Fe is a 3xl 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. Me = rCPcar X Fec (2.47) where rCPcar is the position vector from that point to CP of cavitator. It is assumed that the CP lies on by when cavitator deflection is 0, and its coordinates with respect to body origin O, in this case, are (Xc, 0, 0). Thus from Figure 2.7, the coordinates of CP can be written for the case when the cavitator is deflected. rCPcar = bodv x c +PCGc 0 (2.48) The moments on the cavitator in bodyfixed can be obtained by substituting Equations 2.46 and 2.48 in Equation 2.47. Me = [(Xc t kPC~c~j by kPSrCcb3] x (2.49) Ccpuc , Sic C~cSuc CcaSic Oc SucSic Le/ De c (a) (b) Figure 2. 14 Cavitator: (a) Angle of Attack and Sideslip and (b) Hydrodynamic Forces C6c 0 S6c 0 10 S6c 0 C6c Suc Dc (Orc. y) 0 0( Cac Le (acG 4y~) L Lan S a ha Figure 2. 15 Fin Geometry 2.3.2 Forces on Fins Fin forces are also extrapolated from the CFD database [5], which gives the values of coefficients of lift (clf;,) and drag (cdf;,) for fins. These values are functions of angle of attack afin, fin sweepback 8fi;,, fin rotation 6f;,, fin immersion Ii,, and fin geometry. Figure 2. 15 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 2.15. 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 If,, can be defined as the ratio of the wetted length of the fin to its true length. If,, = (So/S) p;, (2.50) Ip;, determines the total hydrodynamic force on the fin. Fin Rotation (87;;,): As defined earlier, this is rotation about fin f2 axis. This deter mines the direction of the hydrodynamic force. Thus fin rotation is used as primary control for the torpedo. Fin Sweepback (8f;;;): 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 2. 15 and section 2.2.4. The database gives clfi, and cdfi, as a function of my;,,, 8p,, and Ip,,, thus lift and drag on the fins can be calculated by the normalizing factor. Lan = ~p2 pVS inc (2.51) 1 Din = 2pV S injcdpin (2.52) Where Sfin is the length of the fin as shown in Figure 2. 15 These forces have directions as shown in Figure 2.15. Before substituting in Equation 2.39, these forces have to be transformed to bodyfixed reference frame. This process involves following two rotations: 1. Rotate the frame Ffin (which has Lan and Dfin along its basis vectors) to align it with the finfixed frame using Equation 2.14 and 2. Rotate the above obtained finfixed frame to obtain the bodyfixed frame using Equa tions 2.3 to 2. 11. Thus the forces on the fins in bodyfixed frame axis can be obtained. *Rudder 1 ceR1 S6R1 0 0 0 1 S6R1 ceR1 0 _ cRIcaR1 caRISPR1 SpR1 c R1 cRISaR1 SaRISPR1 ceR2 S6R2 0 0 0 1 S6R2 C6R2 0 _ CPR2CtR2 CaR2SPR2 SBR2 c R2 CPR2SaR2 SaR2SPR2 C8R1 0 _SSR1 SaR1 0 CaR1 _ CsR2 0 SSR2 SaR2 0 CaR2 _ (2.53) (2.54) CSE1 0 SSE1 SaE1 0 CaE1 COEl S6El c 0 E10tcal SBE1 CE1SaE1 S 8 0 (2.55) FR1~ * Rudder 2 FR2,x * Elevator 1 (FE1,x 0 SSR1 1 0 0 C8R11 DR1L1 0 SSR2 1 0 0 C8R2 DR2 0 LR2 6El 0 El 0 c El PE SaEE1SE1 0 S8El 1 0 0 CSE1 DE *Elevator 2 FE2,x C6E2 S6E2 0 CSE2 0 S8E2 FE2L, S6eE2 8E2 FE2,z B 0 0 1 SSE2 0 CSE2 ( ) (2.56) SB2CE2 CESE2 0SE 0DE CE2SaE2 SaE2SPE2 CaE2_ LE2 Equations 2.53 to 2.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 2.47. Mfin =r nGCP x Ffin (2.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. r nG~root t, 1' ; 2 + Z$1 t3 (2.58) r otP ; 7 2(2.59) where (ft,}2,}) is finfixed coordinates for corresponding fin. Equations 2.58 and 2.59 can be added by using matrices given in Equations 2.3 to 2.11. Thus, the position vector from CG to CP of fins can be obtained. Rudder 1 XR1 ~ Yiot 8R1 S6R1 0 YR 1 oo'3tR10 0 ZR1 B Z ,ot B S6R1 ceR1 0 (2.60) C8R1 0 SSR1 ~: 0 1 0 Ayo SSR1 0 C8R1 0 R Rudder 2 XR2 X coPt C6R2 S6R2 0 YR2 = 7 ootR  ZR2 B Z~o~t B S6R2 ceR2 0 C8R2 0 SSR2 2XC ~ (.1 0 1 0 Ay( SSR2 0 C8R20 _ceEl S6El CSE1 0 0 1 SSE1 0 S6El cE1 CSE1 _ (2.62) YE2 = root +t S6E EX2X' CE2 S6E2 0 ZE2 BZrroot B0 0 1 (2.63) CSE2 0 S8E2 ~ E~ 0 10O SSE2 0 CSE20 Equations 2.60 to 2.63 give the position vector from CG to CP of the fins. These equa tions in conjunction with Equations 2.53 to 2.56, used in 2.57, gives the external moments on fins about the CG. Man, = by Xfin Ffin ,x Ffin~v Zfin Ffin ,z (2.64) 2.3.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 e^3 axis. Thus, the gravitational force can be written as in Equation 2.65. Fgray = mg^3 (2.65) * Elevator 1 SXEl ot El = 7 oot + X~O ZE1 Zroot B * Elevator 2 0 1 0 El Equation 2.65 can be reexpressed in body frame of reference using Equation 2. 1. COCY CGSY SO F gray = CY ~SOS8 CSY SASOSY +CCOC SCOC 0g CASOCY +SOSY CAS8SY SCY CIOC8 m E(2.66) = mg SOC8 SCOC 2.3.4 Equations of Motion Now that a mathematical formulation of forces on the torpedo has been achieved, these equations can be substituted into Equations 2.36 to 2.42 to obtain the dynamic equations of motion. Thus the force equations of motion can be summarized as in Equation 2.67. 2.3.4.1 Force equations ui +tqw vr prop FR1,x m vt ur pw = imrinFlnn+ 0 + FR1,y w + tpv uq 0FR1,z FR2,x FE1,x FE2 ,x SO FR2,y + FEyI +( FE~Sy +mg SOC8O + (2.67) FR 2 ,z FE 1, z FE 2 z CICO C6c 0 S6c Ccpuca CucSic Suc Dc(ac,T y) 0 1 0 Spc Ccp 0 0t S6c 0 C6c C~cSuc SucSic Cac Le(ac,71) 2C Some issues should be noted for Equation 2.67. *The planing forces have not yet been calculated in the equations of motion. The calculation of these forces will be explained in Chapter 3. In the case of fixed cavity model for the torpedo, these forces are neglected by assumption that the vehicle is centered in the cavity. COSY SG CASf SAS8SY +CYCO SCOC8 SOSY CAS8SY SOCY COC8 (2.70) *The propulsion force is constrained to be along negative by axis. 2.3.4.2 Moment equations li0 0 p1 by~ bg3 by bg b3 0 IO qr +p qr = Mplaning+t XR1 YR1 ZR1 0 0 Is r I p Igq I3r FR1,x FR1,y FR1,z by bg b3 by bg b3 XR2 YR2 ZR2 + XEl YEl ZE1 + (2.68) FR2, FRy F2,z FE1,x FElyI FE1,z by bg b3 by bg bi3 XE2 YE2 ZE2 +Xc +t AcCPC, 0 ACPSGc FE2x F2, FEzFc,x Fcy Fc,z Some issues should be noted for Equation 2.68. Some of the terms in Equation 2.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 moment due to planing forces have not yet been explained. Also, the moment due to planing forces is to be neglected in the case of fixed cavity model. 2.3.4.3 Orientation equations O O pocc e =~S 0 ~ C SO q~ (2.69) 2.3.4.4 Navigation equations x COCY y =I CYSSOS8 i CASOCY +t u v Equations 2.67 to 2.70 are a complete set of equations of motions for the supercavitating torpedo. CHAPTER 3 CAVITY AND PLANING DYNAMIC S If the body is immersed in a steady flow, it is reasonable to assume that the shape of the cavity is constant and axially symmetric. It may not, however, be so as cavities may pulsate due to reentrant jets located at the aft ends of a vehicle [4]. The cavity shape is generally specified by giving the length L and the maximum diameter d,, of the cavity, or the fineness ratio L/d,,. The diameter is specified as function of the arc length along the centerline of the cavity. Various models, analytical, empirical and numerical, have been developed over the years to determine the above mentioned parameters for a cavity. Waid [9] derived formulae for the cavity shapes of disks from measurements made on the cavities. May [4] lists various cavity models based on experimental data for geometries such as disks, cones, spheres, etc. Many models exist for supercavity behavior based on slenderbody theory, boundaryelement methods, and Reynoldsaveraged NavierStokes solvers. However, for timedomain simulation of vehicle dynamics, it is desirable to use simple models. This is particularly the case when we study the guidance, control and stability of the overall body. Two cavity models have been used in the simulations. They are the Munzerreichardt cavity model and the Logvinovich cavity model. Both cavity models assume the cross section of the cavity as being circular in shape. The radius of the cavity is specified as a function of the arc length along the cavity centerline. The modeling of the cavity, for simulation purposes, comprises of these cavity models and the description of the cavity centerline based on memory effects. 3.1 MunzerReichardt Model Reichardt showed that the fineness ratio of the cavity L/dinax is a function only of the cavitation number o. In other words, the ratio is independent of the cavitator shape and is only function of o. The formulae for the cavity dimensions according to the Munzer Reichardt [10] model are R = R,L, 1 Lazx 2 ,, [124 0. (3.1) R,, = RcClcdel.350.9 In the Equations 3.1 Lza is the length of the cavity, L is arc length along the cavity center line, Cd, is the cavitator drag coefficient, RC is the radius of the cavitator and R the cavity radius at that arc length on the cavity centerline. 3.2 Logvinovich Cavity Model 3.2.1 Logvinovich Theory of Independent Expansion Logvinovich has made the following fundamental observation: "Each crosssection of the cavity expands relative to the path of the bodycenter almost independently of the sub sequent or preceding motion of the body..." [1 l]. Logvinovich's cavity contour equations are another example of the kind which consider the cavity to be in the shape of an ellipsoid. The following formulas are taken from Logvi novich [11] for the cavity shape R = Rk RR 1 2k1 1tk RR2\ 1 t\ t R = k 1 1I 1 (3 .2) R, Rg i tkj tk In the Equations 3.2 Rk denotes the maximum cavity radius, R and Rt, denote cavity radius at a particular point in time on the cavity centerline, t time taken for cavity at that point to evolve from inception, tk is the time taken by cavitator to travel one body length of the vehicle, K is a correction factor and R is the time rate of change of the cavity radius. Explicit assumption is made of the cavitator being of disk shaped in deriving the above formulae. However, in our case the cavity due to a cone shaped cavitator is approximated with the one due to a disk shaped cavitator. The formulae in the above form make it possible to calculate the cavity profile also when t/tk > 1, but usually when t/tk > 1.5 the boundaries of the cavity are undeterminable. They start breaking up and foam begins to form. At t = 0 the contour expressed by Equation 3.2 is matched with that of the leading part of the empirical cavity 3x R = n (3.3) where Rn is the cavitator diameter. The maximum cavity radius Rk can be calculated using the formula 1=, +08 io(3 By selecting xl = 2Rn as our matching point we obtain R1 = 1.92Rn. A cavity "correction factor", K, is also required in the calculations. For the selection xl = 2Rn and K = 0.85, a good correlation can be obtained between the formulae in Equation 4.2 and the experi mental data. Under these conditions the cavity half length can be approximated using the formula in Equation 3.5 [l l]. L I (1.923.5 The cavity half length will be used to approximate the value of t/tk as L/Lk. The validity of this approximation can be seen from Equation 3.10. From the above formulae we can calculate the basic dimensions of the cavity. However, during simulation the coordinates of cavity centerline need to be located. The centerline is required in order to determine the final shape and position of the cavity with respect to the body. Effects like buoyancy and downwa~sh need to be considered for determining the final shape of cavity centerline. 3.2.2 Cavity Centerline The arc length of the actual path is assumed to be approximately equal to the distance from the cavitator to the planing point projected onto the body axis. The time lag in con sidering the memory effects shaping the cavity can be obtained by the expression z = L (3.6) Where, L is the distance along the cavity centerline. The above value is obtained through a 2tad order approximation = u~) +u~t) ut~r(3.7) The average velocity over the time period [tl t z] multiplied by the time period will give the distance covered, L, in time t. The distance L can be written as L= u t &TT&T+u=0(38 Now, solve the qluadratic equation for t, ui u2 4( ) (L) ui V'u2 2u I = (3.9r) Taking the limit of ii approaching zero and using the L'Hospital rule for calculating limits for indeterminate forms we have lim I=4 (3.10) iFromn one y... ..t of view, the variation in the cavity shape at a point in space is due to the memory effects associated with the i : :  of cavitator through that point at a i :; T us time. Thus to calculate the cavity i the cavitator position at each timestep is stored. The cavity centerline is formed by the i : .I II :: of the cavitator tip at current time and the previous (n1) time i i where 'n' is the number ofl 1...1= on the full length of the cavity. We can write yi = ye ( i 1At)(3.11) yizi z 1 i ) / where (xidyi,zi) are the coordinates of the ith point on the cavity centerline, (xc,,7,zc) are the coordinates of the cavitator tip and At is the time differential between points i and i 1. The buoyancy force acts on the cavity and the centerline is distorted over time. One way of representing the buoyancy force is by a constant acceleration acting against gravity. Now taking into account the effect of buoyancy on the cavity centerline, we have xi xc (t (i 1)At) zi zc (t (i ) At) b((i ) At) where b is the constant upward acceleration due to buoyancy on the cavity. In contrast Logvinovich [1 1] has a different form of equations for the buoyancy and downwash. They are given by the following expressions hbuoyancy (X) = 7 v2J RbXdx (3.13) F ,xdx hdownwash(X) = , (3.14) npf o R2(x) In these equations, X is the position on the body centerline in frame B, R is radius of cavity as function of position along body centerline, hbuovancy and hdownwash represent the vertical displacements due to buoyancy and downwash respectively at a point on the body centerline, Ok is the volume of the cavity enclosed by the crosssection at that point on the body centerline and V is the velocity of the centerofmass of the body. In the present work the cavity model is limited to the effects such as buoyancy, downwash and memory effects. Other complex phenomena associated with the shape of the cavity are ignored in the current investigation. 3.3 Planing Model The stability of the supercavitating vehicle does not depend simply on the hydrody namic coefficients as in the case of a fully wetted missile, but rather on the moments asso ciated with the nose and the aft. The contact of the aft of the vehicle with the cavity wall gives rise to the planing force. This force, depending on the strength of the tailmoment it generates, may send the body into tailslap. When sufficient restoring nose moments exist the body may cease oscillating but still have a part of the hull protruding outside the cavity. In this case part of the lift required by the vehicle to maintain it's flight is provided by the planing force. The paths an unpowered proj ectile may take for various such planing configurations has been given in May [4]. The lift needed to support the aftside of the body may be generated by either the fins, the planing force or a combination of the both. In the present work the Hassan's planing force model has been used. Figure 3.1 shows the vehicle undergoing tailslap. The planing force (Fp) and moment (Mp)are shown acting at transom of the planing section. Plaming Figure 3.1 Vehicle Undergoing TailSlap. 3.4 Planing Force Equations Wagner's planing theory [12] assumes planing and continuous immersion to be identi cal. The theory is mainly developed for 2dimensional bodies undergoing planing. Using potential flow, the immersion force on a flat wedge is calculated assuming that the force is due to fluid transport into a spray sheet. The force per unitlength of the cylinder is com puted and then integrated over the entire wetted area. The following are the equations that characterize this situation as given by Wagner [12].  a it w l Body Wall Figure 3.2 CrossSection of Cavity Assume h = (R rta2/ > 6 = 2tan h = hs A =R r (3.15) In these equations h is the depth of the penetration of the body outside cavity at transom. , R is the radius of the cavity at transom and r is the vehicle body radius at transom. Figure 3.2 defines these variables. The force is given by P = m*h;_ hl ~~t t2 (3.16) where dm = pnR2 1 + COS2 Sin2e = pniR2 (3.17) Letting dm* am* ah (3.18) dt ah at We can write Sh 1 +t 2h j~( 9 The term m*h~ is usually neglected. If we assume h(x) = ho xtan(u) I = (3.20) tan(u) lam* pnsR22 mf = h hhortan(,~dx = tan(u) 1 A h (3.21) It can be concluded that the planing force is given by ,=Y ( r +t ho 2n/(.2 Hassan [13]i, based on the Wagner's theory, has given the following set of planing force, moment and skinfriction equations. The term Fp denotes the planing force, Myp the planing moment and Ff the planing skinfriction force. 1 Rr rtho\ FP 2 ho+tR r) jr+2ho 1 hr +t ho Mlr = pw V2" (nR2) cosucostIh tRr t h (3.23) Dummy variables uc, us and S, are introduced into the skinfriction equation. These variables are defined by the following equations s r (3.24) Sw = 4r ((1+uc2)ATan(uc) M) + 2~(Rr)tn Os .5) ASin(u,) + 0.5us (1 us The final form of the skinfriction equation can be written as Ff. = 0.5pV'"cosolcosuSwCdt (3.25) Note that Hassan's equation for Fp can be derived from the Wagner's equation for Fp in Equation (6.8) by using the value h = Vsina. 3.5 Planing Kinematics 3.5.1 Calculation of Immersion The main assumption involved in the calculation of the immersion is that the cross section of the cavity is circular. Strictly speaking the cavity crosssection would undergo deformation due to gravity and other effects. It may become elliptical in shape [14]. Using this assumption we can, by simple geometric calculations, find the immersion depth, ho, at transom as follows ho = r R+Rt (3.26) In Equation 3.26, Rt is the distance between the centers of the circular crosssections of the cavity wall and the body. Because the crosssections are circular, a line passing through center of cavity and body also passes through the point of deepest immersion. 3.5.2 Method of Calculation of a This section describes the calculation of the angle a that appears in the planing theory. It is important not to confuse this a with the numerous places that a represents the angle of attack. Figure 3.3 shows the angle a as the angle of inclination between cavity centerline and body centerline. 3 Body Figure 3.3 Planing Section Define the following terms, y, : change in cavitator velocity over the interval [t z, T] expressed in dynamics frame B. Subscript 't' denotes transom or planing section. V : Velocity of cavitator at current time. Vc:Velocity of cavitator at the time: the: cavity at the i .1 ::1:: section formed. Rc : Location of cavity centerline in the dynamic reference frame. TF is the tangent to the cavity centerline at the planing point. It is also the normal to the crosssection. Note that calculating i1:; r i of deflections l::T:;1 that crosssections are elliptical; this will be ignored for the time being. The angle a that :i i :: in the planing theory is the angle between the body centerline and the cavity centerline. The theory assumes that cavity centerline and cavity surface are: parallel. This is not the case if the cavity radius is changing. The angle a is : :::i: . by finding the transformation matrix that rotates the: velocity vector Vc so that it is the unit normal i" = {1,010}"T in the frame B3. This transformation matrix rotates a vector in the dynamics frame BZ to the cavity normal frame. The: defining .;:: .0. :: is: o1 The rotation angles for T~ can be found using fx 1 cos(6) ... 9) defines 10tation sequence Introduce another rotation to account for the rotation of the cavity surface due to ex pansion of the cavity Rx, sin(ex) =   cl Letting cos(6x) 0 sin(6x) T3 = T3,2,l (0, 8x, 0) = 0 1 0 (3.28) sin(6x) 0 cos(6x) and T = T3 T2 we've got our transformation from dynamics coordinates to the frame in which the planing theory is defined. The angle a appearing in the Figure 3.3 is obtained from costu) = < Ti, l > > cos(u) = Tr y (3.29) This is the angle between the vehicle centerline and the cavity surface at the planing point Body v. Wal Figure 3.5 Velocity of Transom Relative to Frame B We need to include the effect of the velocity of the transom on the angle a. V,,= o x R, (3.30) In the above equation, V,, is the velocity of transom relative to the frame B. Rp is the position vector of the transom in the B given by { Lp,, 0, O} where LI, is the position of transom on the 1 axis. In order to include the effect of the velocity of transom, F),, on the angle a, F), is projected onto the vector R ~. Then the angle eR, as shown in Figure 3.3, is computed by the following equation sin6 = V(3.31) In the above equation V denotes the velocity of the vehicle. So, the angle alpha used in the planing force equation can be finally written down as a = u+6t (3.32) Cavity wall lfho Body Section Figure 3.6 Calculation of angle of plane inclination a with assumption that the cavity sur face is planar. A simple way to calculate a reasonably accurate value of a during computations is by the Equation (6.19). L tan(u) = (3.33) Here, we are assuming the cavity wall to be planar. Figure 3.6 shows the angle of plane inclination, a, with the planar assumption. 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 lin earization can be carried out similarly, as shown in Nelson [15]. The equations of motion, as in the case of a supercavitating torpedo, are represented by a set of firstorder differential equations. ri = f (x, u) (4.1) using f : 2" 2" as a nonlinear function of a timevarying vector x E "Z and u E 2". 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 2.67 and 2.70 can be written in simplified form using sums of total forces and moments acting on the body. m (u +tqw vr +t gS8) = X m (vtru pw gC8SA) = Y (4.2) m(wi!tpy qu gC8CO) = Z Ixt pt + r (I, ly) = L (4.3) ly@FFx )= M (4.4) Iz qI I)= N (4.5) 8O = )I 46 A 1 SO COr x COCY CGSY SO y = CYSSOS8 CSY SAS8SY +CYCO S~IOCv(47 z ~CASOCY +SOSY CASOSY SCY CCOCw E( wv 4 7 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, Y, 8, 4, x, y, z } (4.8) u = {61, 6R2, 6E1, 6E2, 8R1, 6R2, 8E 1, E2, 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,Po, go, ro, Yo, Oo, Go, xo,yo, zo } (4.9) = {uo, 0, 0, 0, 0, 0, 0, Go, 0, 0, 0, O } 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 +t 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. However, the following analysis is just an explanation of the method. Eventually, the linearization is carried out numerically. Thus, in the numerical method the planing forces are included, unless specified. The linearization procedure is resolved for the force equation in by direction. This equation relates the force, X, to the states. m (u tqw ru tgS8) =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. mgS~o = 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 +t AM) + (40 +t Af) t o +t A) (r0 +t A) (uo +t AM) dt (4.13) +tgS(Oo +t AG)] = Xo +t AK 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+ gAGC~o) = AK (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+ gAGC~o) = AK m (Av +t uoir g^0000) = AY (4.15) m (Au uo~q tgAGS~o) = AZ 4.1.3.2 Moment equations Ix Ap = AL ly~q = API (4.16) Izr = AN 4.1.3.3 Orientation equations C~o AG = Aq (4. 17) AQ = Ap+ T~o~r 4.1.3.4 Position equations Ax= SO,,an,,G+COI ~n u+S~o~w Ay = Av (4.18) Az = CO,,an,,G SOI ,Au+tC~o~w 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. K =X,Au+_tXvOv+tX,~ w +X,Ap+XAq+Xr, (4.19) + XeR1 6R1 +t X8R2 6R2 +t X8E1 OEl 1+ X8Ez 6E2 + XaR1 R1+XaR 6R2+tXaE16E 1+XaE 6E2+tXac i where the sub scripted X represents its partial derivative with respect to its sub script. aX x, = (4.20) au X=XO 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, SEl s E2, 6Els 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. B (FR1,x +t FR2,x +t FE1l,x +t FE2,x +t Fe,x +t Fprop ,x) (4.22) Expressions for each of the terms in Equation 4.22 have been derived in Chapter 2. For example, the expression for the force on cavitator is shown in Equation 4.23. Fc~x= Ci 0 Sc Sc Ci 0 (4.23) The velocity components in Equation 4.27 can be found using Equation 2.2. In Equation 4.23, ac, Pc, and thus Le and De are the only functions of u. Thus the partial derivatives with respect to ze can be obtained. acx= cCac 0~ r Sc, a 0~s, + 0~a ~ Le(ac,Y ) SC  Dc(ac,: S 0T N I 0 S6c Ocpuca CcaSpc Suc 1 0P Sic Oc 0 C6c OcpSuc SucSic Cca (4.24) It can be seen that a s, ,L and D are terms required to be calculated. These can be calculated from equations 2.16 and 2.17, which are restated in Equations 4.25 to Equation 4.27. tan(ac) = '" 21c (4.25) (4.26) (4.27) tan(Pc) = Ce V;2 =t + t + w; W C6c 0 0 1 u6 o w BI be obtained S6c 0 C6c (4.28) by $2 ~3 p q r (4.29) xc ye zc for the reference flight condition that is L IB Now the velocity components can stated in Equation 4.9. 2 ucdwe we due sec2 (ac)due= ucdwe we due due = u? +t w, CGcdwe SGcdue due = 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. (Vcdvc vcdYc) di = c Vc Y V dvc = at (xo, uo) (4.33) Now, using Equations 4.28 and 4.29, variation of velocity components of cavitator with respect to u can be obtained. auc = C6c = 0 = S~c (4.34) uc CGcuo The variation of ac can be obtained by differentiating Equation 4.25 at reference flight condition. Finally, combining Equations 4.32 to 4.34, the variations of ac and Pc with respect to u can be obtained at reference flight condition. ue C6Sc Sw S6,C Bu C~cSc S~C~c(4.35) uo uo apc 1 avc Bu uo 8"(4.36) _O 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. Le = ~plc2Sccle 2 (4.37) De = ~p~c2Scdc These forces can be differentiated by chain rule the derivative would be like in Equation 4.38. ta rpSc 2V~ccle + Vc (4.38) 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. ;uc +v +weu lr? (4.39)) u? + v + w Thus the derivative of the lift and drag forces can be obtained. aLe = pScYccle (4.40) aD, = pScYccdc (4.41) au Thus all the derivatives required to calculate righthand side of Equation 4.24 have been calculated. 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. The Mathematica code for linearization is described in Chapter Appendix C. 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. i = Ax+tBu y = Cx+tDu x E Zn,u E ~Z,ye Zm (4.42) Ae Znx Zn,Be Znx Zp C E Zmx Zn,De Zmx X F 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 =[n nA Aw Aq AG Av, Ap Ar Y Ax Ay Az (4.43) u = ic8El 6E2 6R1 6R2 8E1Bl BE R1 BR2 bfypro 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 {Y,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, AG correspond to longitudinal dynamics, which also means the dynamics in $1bi3 plane. The variables Av, Ap, Ar A correspond to lateral dynamics, which is the dynamics in $1b2 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 = AopAlong Acoupl~lt (4.44) 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 dynam ics. Acoupl and Acoup2 are coupling matrices. It is required that the coupling matrices become negligible 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. (4.49) Xsc m Blong . 0  From linearized equations, the four parts of the A matrix for the torpedo can be written as in Equation 4.45 to Equation 4.48. x, xw +x g~+x m 90 mt mw m gO tX o+z, z, uo +z gC~oS~o + ze~ Alon m m m m (4.45) 0 0 C~o 0 uo + gC~oC~o + Lr(l4 )qo Le Ix Ix C~oT.~o 'ocoose0,,1 ,,1~H vo +x ~P zv gS~oC~o $ m~ m Z)ro Mr(I4;oM S~o S ogo  C m Lp Ix Np lyx)90 1 Alatd Acoupl (4.46) (4.47) (4.48) ro +x Ix N, 0r  xz L, Ix Nw 0 gS~oSpo +~ Le No S~ogo +t C~oro z ~oro S4(o YTOo 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. XFprop,x 0x1 Tsc YFprop,x Blatd = (4.50) 4x10 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. Xlong = Ihu Aw Aq AG] xlatd = Av Ap Ar A' H = 16c El: 6E2 6R1 6R2 B8El 8E12 BR1 BR2 EFprop (4.51) Xilong Along Acoupl1 Xlong Blong 1x Ailatd Acu2Alt latd Blatd 8xl 8x8 8xl 81 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, 8 v, p, r, Control Sc, 6el, 6e2 6rl, 6r2 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. A deviation from these requirements will be observed in some of the controllers, mainly to improve performance and provide stability. There are various control methods, like linear quadratic regulator (LQR) synthesis, pu 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. pusynthesis 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. This chapter explains various problems associated with the control synthesis and the system model used for synthesis of the controller. 5.1 OpenLoop Performance for the Fixed Cavity Model Initially the equations of motion for the torpedo have been linearized for straight and level flight at a forward velocity of 75 msl x= {75, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, O} (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. H= ( cel,~e2,6rl, 6r2, Fprop (5.2) = {0.0067, 0.0106, 0.0106, 0, 0,4.0023e+t03} 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 2. 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, 6e, 8e,2, ir, in }. The matrices corresponding 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, Y}. 1 ~~~NL Equation~opd Integrator at For simt Clock Time Figure 5.1 ::::::i::. 1 Model for Open Loop Simulation state r Plotting Control ~ Trim 20 15 5 ' 0 01 0 03 000 20 40 6 8 00 times) oO 20 40 0 t9 4 3 5 P25 1 5 Fi:: 5.2, O  : rLoop i i..: for Torpedo: w?,p, q VVVVVVVVWVVWvywywywwyvi 4.5204 0.2616 0.0000 32.3010 406.0942 158.4675 0 12.0422 0.1813 1.1437  0 1.5417 15.7648 1.2077 0 1.3110 78.5888 3.5614 1.0000 9.8100 0 0 0 Along Blong Alatd = (5.3) (5.4) (5.5) 69.0608 69.0608 303.3736 303.3736 45.1531 45.1531 0 0 0.0002 54.2281 0.0025 1.0000 0 71.6011 0.3004 1.2528 0 1.0000 9.8100 0 0 0 0 0 14297.086 1.4129523 0 0 0 14297.086 1.4129523 0 0 366.60511 17276.994 54.561629 0 0 366.60511 17276.994 54.561629 0 0 Blatd (5.6) 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 +t 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 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. 1 Longitudinal Egienvalues for u=60:5:90 10 Lateral Egienvalues for u=60:5:90 0.55 .E E 0.5  1 10 0 20 10 ~~og 0 10 0 60 ~a (La 20 0 20 (a) Longitudinal (b) Lateral Figure 5.3 Variation of Eigenvalues with Change in Velocity 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. 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. Cavitator Deflection 15" I Sc I +t15" Cavitator Rate 25" s I Sc I +t25" s fins 60" I 6f I +60" Fin Rate 25" s I 8f I +25" s Thrust 0 < F,, < 30000N Table 5.2 Control Constraints 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 = ~. 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 5(t) = qachi(t) qcomm (t)  subj ect to ui i = Ax+tBu where, umin and umax are lower and upper bounds on control deflections. The quantities uin, and umax are lower and upper bounds on control deflection rates. The problem becomes a disturbance rej section 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 1800. 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. 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 Figure 5.4 Loop Gain 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: * +20% error in Cl of Cavitator. * +20% error in Cd of Cavitator. * +20% error in Cl of all the Fins. * +20% 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. 5.3.5 pu analysis: For the Hoo/p analysis [16] it is desired that the peak pu value for the closedloop system be close to 1. This ensures good robustness atleast equal to the values of uncertainties in the synthesis model. 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 [17]. 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)= (x'Te Qx +uTRu + 2x Nu)dt (6.1) 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. x = Ax +t 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, , x,,} (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.5. u = K (6.5) Equation 6.2 can be rewritten by substituting the new control law. i = Ax+tBu X1 r1 x2 (6.6) = Ax BK 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. i = Ax B kl, k2, , k;, (6.8) = Ax BKx+tBklrl = (A BK)x+tBkl rl 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. i( )= (A BK)x() )+ Bklrl (6.9) The error dynamics can be obtained by subtraction Equation 6.9 from Equation 6.8. i(t) i( ) = (A K) (x(t) x( )) (6. 10) e = (A BK)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 [17]. 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 = (ri xt) (6.12) The state space equation for the system with this modification can be written. i = Ax+tBu e"= ri xi (6.13) = r 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 i = {xl,x2, ,Xny e) T. So a new formulation of the state space equation can be written, AO B i= i+: u+ ri (6.14) > = Ai+tBu+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. 70 K= [kl,k2, , knkn 1] (6.16) Then, the new control law can be written as in Equation 6.17. u = Kxj =[k, k2 k Ilk~l] 11(6.17) = [kl, k2, , kn]x+t [kn 1]e" = Kx tkre which is represented in Figure 6.2. rl kltti ~y=x I K Figure 6.1 Controller for Tracking when Plant has an Integrator. rl + x + Bu X I K Figure 6.2 Controller for Tracking when Plant has no Integrator. 6.2 LQR Control for Fixed Cavity Model: 6.2.1 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. Long = diag( [0, 0, 0, 0, 10]) (6.18) Rlong = diag([ 5, 4]) The first four numbers in the Qlong 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 8. 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 S0.9681 (.9 K=I0.0000 0.0040 0.1016 0.0000 1.4195 1.1184 0.0000 0.0042 0.0995 0.0000 1.3981 1.1308 (.0 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 (.2 0.0071 K = 03 0.0005 0.1253 0.0132 0.0019 0.0000 (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.2.2 Nominal Closedloop Model 6.2.2.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. 10 ~~oo 84 5solP~~O 2 0 ~~0a70 (a) Longitudinal (b) Lateral Figure 6.3 Eigenvalues for the ClosedLoop System 6.2.2.2 Simulations The response of the vehicle to a longitudinal command is simulated and shown in Fig ures 6.4 to 6.6. 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 741 20 10 0 10 20 5 10 15 times) Figure 6.4 P~itch Command T~rackiing : q 1 20 0.5 ~ l1 10 0.5 ~ 1 20 1 5 times) 10 15 30 5 times) 10 15 Figure 6.5 Pitch C~ommand Tracking : C Sc only. This shows that the controller allows pure longitudinal motion to be uncoupled from the lateral m~otion. 1.15 1.5~ 10  00 5 times) 101 105 times) 105 Figure 6.6 Pitch Command Tracking : Sen Bel The response of the vehicle to a lateral, roll rate, command is shown in Figures 6.7 to 6.9. 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. 0 5 10 15 times) Figure 6.7 Roll Command Tracking: p O 5 10 15 05L5 10 15 times) times) Figure 6.8 Roll Command Tracking: Sc, Sc 0.4 2 0.3 2.5 0. 5 times) 101 05 times) 105 Figure 6.9 Roll Command Tracking: Se, Bel1 rl + y=x '3 K Figure 6.10 Breakpoints for Calculating the LoopGain for a Tracking Controller 6.2.2.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.10 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. 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 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) I ~Lateral Gain Margin(db) Phase Margin (deg) 1 22.964(at 0 rad/s) o 2 250.51 (at 0 rad/s) o 3 50.36 (at 0 rad/s) o 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 margmn. 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, 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 refl ect the robustness of the dominant dynamics. 6.2.3 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.2.3 and 6.3 show the percentage variation of the elements of A and B matrices for a 20% change in coefficients of lift cavitator. This comparison is done for cases with changes in other coefficients also. In all cases, 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 u w q Ovpr Y u 0.46 0.62 wl 5.52 6.86 1.58 q 1.05e5 34.8 13.58 v p Table 6.2 Percentage Variation in A Matrix due to 20% Variation in clc in a coefficient. For example, in Table 6.2.3 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. 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.2.3.1 Model Figure 6. 11 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 be seen that the longitudinal dynamics show some perturbation in the damping while the lateral system relatively unchanged. 6.2.3.2 Simulations 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 Sc 6el 6e2 r 6r2 w 20 q 20 v r Table 6.3 Percentage Variation in B Matrix due to 10% Variation in clc 15 6 10 5 4 10 6 11 00 80 6040020 10 0 60000 0 Real Real (a) Longitudinal (b) Lateral Figure 6. 11 Eigenvalues for the Perturbed ClosedLoop System: 20% Error in clfin 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 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 6.12 to 6.13 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 clfin. 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 C C 4 +20% error 4~  +20% error 20% error I 20% error 0. 2 a "2 ~ 1 I 4 44 6 CO10time(s) 203 0time(s)203 Figure 6. 12 Response for 20% Variation in clpn, P, 4 1 1.2  +20% error  +20% error S 20% error 1 20% error 0.5 9 90.8 o 0 0.6 0.50 1 20 30 0.40L 10 20 30 times) times) Figure 6.13 Response for 20% Variation in clpn: 6c, Bel 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. 6.2.3.3 Gain and phase margins Table 6.4 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 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. 81 Table 6.4 Gain and Phase Mtargin for Perturbed Closedloop System: : error in clfin CHAPTER 7 p/lHo SYNTHESIS CONTROL 7.1 Uncertainty The main parameters (besides the control) that dictate the forces on the the torpedo are the coefficient of lift, cl, and drag, cd, of the cavitator and fins. The values of these parameters are obtained using a CFD code in [5]. These values have also been obtained theoretically and experimentally in [4]. A comparison of the CFD and predicted results is shown in the figures 7. 1. It can be seen that there is a difference of at least 1 order of magnitude, in the value of cd, between the two datasets. Also, a supercavitating flow is a 2 phase flow with partial water and partial vapor. Thus the hydrodynamics for supercavitating vehicles are unique and need more investigation. It is clear that there is a need for inclusion of uncertainty for the values of cl and cd, which are main parameters that determine the hydrodynamics. A controller designed based on pu or Hoo theory is a robust controller, which deals with errors and uncertainties in the system, implicitly. Basic idea of pu synthesis is to reduce the gain from error or disturbance to the error in tracking. In the design of pu controller for torpedo, an uncertainty is assumed in the coefficients of lift, clc, and the coefficient of drag, cde, for the cavitator only. The following formulation for the synthesis model has been derived in detail for the longitudinal plant. The formulation for lateral plant is very similar, difference being only the states and controls. Also, the following formulation is for a pitch angle tracking controller. The formulation for tracking any other state, say pitch rate q, would be similar. The feedback to the controller will, in that case, be the tracking error mn q. Plot of Cd values from Database and Maytest Data : May Test Date SDataba~seImmersion=0J 1  Daab~aselImmrson0 3 DatabaseImerasron=0~ 5 DatabaseImmerston=0 7 DatabaSeImmer~con=0 9 ne i1 4 :: L c~z ; Y  Figure 7.1 Calculation of Uncertainty 7.2 Synthesis Model The state space formulation for longitudinal dynamics can be given as x = Ax +tBu y =Cx where x = {Au, Aw,A~q, AG} u = {0,8e} (7.1) (7.2) where, I4 represents an identity matrix of size 4 x 4. For simplicity, the B matrix is chosen to be 4 x 2, i.e., other controls are assumed to be fixed at their trim values. Thus the only controls are deflection of cavitator and the deflection of the elevators. The elevators are A = [aij]4x4 B = [b,7]4x2 assumed to have a symmetric deflection for a longitudinal tracking command. Thus only one elevator can be included in the synthesis model. Note that the symbol u is used to represent the control vector. u will also be be used to represent the forward velocity of the torpedo. As described earlier, parametric uncertainty is assumed in the coefficients of lift(clc, cle) and drag(cde, cde) of the cavitator and the elevators. Let W1 and W2 represent percentage errors in cle and cde respectively. Let W3 and W4 represent percentage errors in clc and cdc respectively. Thus a value of W1 = 0. 1 would mean a up to 10% error in value of cle and so on. Now, the true value of these coefficients can be written as cle = cle (l +t 6cieW1) cde = cde(1+8t cdeW2) (7.3) cle = cle (1 +t ScicW3) cdc = cdc(1+8t cdcW4) where, cle, cde, cle and cdc are the values of these coefficients from the database, i.e., their nominal values. 6ct and Scd are variables, whose numerical value determines the actual error in these coefficients. Thus (7.4) Elements ofA and B matrices will be functions of function of coefficients of lift and drag of cavitator, and some other terms which are of lesser importance for this analysis. In theory, the parametric form can be written as explicit function of cl and cd i.e. A = Ao tcleA t +cdeA2 +tclcA3 +tcdcA4 (7.5) B = Bo +t cleB1 +t cdeB2 +t clcB3 +t cdcB4 This formulation of A and B matrices can be tough to obtain for general flight condition, but it can be found numerically MATLAB, as described in the Appendix B. A symbolic derivation of these matrices is done using MATHEMATICA, but it is found that the terms are very long and the results runs through pages. The description of various codes has been given in Appendix C. The numerical values from both the codes match, thus leading to a verification of the numerical method. Now, using equation 7.3 in 7.5, we have A = Ao+ cle (1+ ScieW1 )A t +cde ( 1+8cde W2)A2 tC, +l( 1+8t ciW3)A3+ cdc (1+ tcdc W4)A4 B = Bo +t cle ( 1 +t Sci W1)B1 +t cde ( 1 +t cde W2)B2 +t clc ( 1 +t 8ci W3)B3 +t cdc ( 1 +t cdc W4)B4 (7.6) A = Ao +t6IciA t +8cdeA2 +8 cicA3 +8 cdcA4 (7.7) B = Bo +t6IciB1 +8 cdeB2 +8 cicB3 +8 cdcB4 where Ao = Ao +tcleA t +cdeA2 +tclcA3 +tcdcA4 Al = cleW1AI A2 = cde W2A2 (7.8) A3 = clcW3A3 A4 = cdcW4A4 and Bo = Bo +tcleB1 +tcdeB2 +tclcB3 +tcdcB4 B1 = cleW1BI B2 = cde W2B2 (7.9) B3 = clcW3B3 B4 = cdcW4B4 With this transformation the state space form shown in equation 7.1 can be written as x = (Ao +t6IciA t +8cdeA2 +8 cicA3 +8 cdcA4)X (7.10) + (Bo +t6cieB1 +8 cdeB2 +8 cicB3 +8 cdcB4) M In above equation, each of the terms corresponding to uncertainty, are like inputs to the system. For example, 6clAix is like an input, and is of size 4 x 1. This term can be replace by an input vector wl of same size. Also the term multiplying the uncertainty can be written as an output of the system. Thus lets make a further transformation as follows 11l 212 213 214 wit w13 W14 Aix (7.11) Scie 0 0 0 0 Scie 0 0 0 0 Scie 0 0 0 0 Scie similarly z2 = [Z2i] z3 = [Z3i] 4 4Zi T zg = [Zsi] Zg6 6Zi T z7 = [Z7i]" z8 8zi T [w2i]" [wati]" [wqi] T [wsi]" [wsi T [7i] T [W~i T SBiut A2x B2Ue = A3x = B3u A41 B421 = Diag[8ie]4x4 2 'Diag[6cde 4x4 3 : Diag[6cde 4x4 4 = Diag[8cic]4x4 5 = Diag[8cic]4x4 6 : Diag~[cT._, 1]Jx4 7 : Diag~[cT._, I]Jx4 8 (7.12) for i=1 to 4, and Diag[] represents a diagonal matrix with the diagonal elements given by those in [] and remaining elements as zeros. Using above transformation in equation 7.10, we get a new state space form as x = P~11x+Piu P = P21x+tP22u (7.13) where: U = w" w12 11'3 11'4 5/' M6 M7 *' 7 4x (7.14) I = :zl ~I z3 z4 5j 6 ;7 8g .716 and P?11[ = .x4 (7.15) pig = 4 4 4 4 4 4 4 4 ;13 (7.16) At 041x4 AZ 04x4 =A3 (7. 17) 04x4 A4 04Ix32 04x2 041x32 R1 04x32 x / = 04x32 04x2,(. 8 04x32 L3 04 x32 . x2 04x32 04x32 '. x2 I___ ___ ___ ____ ___I I_ x u Figur 7.2Linea Frationl Repesenatio Thus The ne sytmi 4iptad3 outu rsyse.Thsrpesnain oh systm i kno asLiner Factl Reprsenttion an a esonbabokdarma shown in figure 7.2.___ With thsfomlain,Z aW snhsscnrlecabefudsdsribe n[6.Fg ure 7.3 shows the sythsi model fr the system showsn i nation 71.T con o feedback signal uncertainty, a multiplicative uncertainty is added in the feedback signal for 8. The W, represents the multiplicative uncertain in 8, whose values is W,A,,. Thus, Oactual = S~oninal(1 +t WmA,,, (7.19) Various filters have been added to system to put constraints on the control surfaces and improve performance. These filters are like the Q and R matrices for and LQR control synthesis. These choice of these filters done both by trial and error and also some analysis of their frequency response. * W,: This is chosen so as to penalize the performance error. This is usually a low pass filter. Thus the filter has a high value at low frequencies and the value drops down at higher frequencies. The value is obtained by trial and error. WKc and WKe: These are the filter on control commands. These are usually high pass filter. The control surface motion allowed is the inverse of these filter. Thus a higher control motion will be allowed a low frequencies and low control motion at high frequencies. Isc 0 0 0 Is Sc4, O 0 0 0 Isc 0 0 0 I Scdc Auto Pilot Input Oo Figure 7.3 Synthesis Model for pu Controller \' :y: This signal is to add sensor noise in the system. In this synthesis, this is chosen as a constant value. A better approximation of this can be done: by knowt~ing some: sen sor character sticGs. Ac: This filter .. an actuator model. This represents the actuator dynamics in the frequency domain. The ::i.::i d and .: Hi.=:: e for the synthesis model are noise, ..7 (7.20) ek WK d= w e = z eg ek 