<%BANNER%>

Control strategies for supercavitating vehicles

University of Florida Institutional Repository

PAGE 1

CONTR OL STRA TEGIES FOR SUPERCA VIT A TING VEHICLES By ANUKUL GOEL A THESIS PRESENTED T O THE GRADU A TE SCHOOL OF THE UNIVERSITY OF FLORID A IN P AR TIAL FULFILLMENT OF THE REQ UIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE UNIVERSITY OF FLORID A 2002

PAGE 2

A CKNO WLEDGMENTS I w ould lik e to e xpress my sincere gratitude to my committee chairman, Dr Andre w K urdila, for his in v aluable guidance throughout the course of this project. I w ould also lik e to thank him for gi ving me this opportunity to w ork on such a f ascinating project. I w ould also lik e to thank my committee cochair Dr Richard C. Lind, for his in v aluable guidance and inspiration throughout the project. I w ould also lik e to sho w my sincere appreciation to Dr John Dzielski, Dr Nor man Fitz-Co y and Jammulamadaka Anand Kapardi for their v aluable contrib utions to this project. I w ould also lik e to e xpress my gratitude to all the members, past and present, of the Superca vitation Project. I w ould also lik e to thank the Of ce of Na v al Research for the support of the research grant for the project. On a personal note, I w ould lik e to thank all my friends and f amily members whose support helped me to aim to w ards my goals. ii

PAGE 3

T ABLE OF CONTENTS page A CKNO WLEDGMENTS . . . . . . . . . . . . ii LIST OF T ABLES . . . . . . . . . . . . . vi LIST OF FIGURES . . . . . . . . . . . . . ix ABSTRA CT . . . . . . . . . . . . . . x CHAPTER 1 INTR ODUCTION . . . . . . . . . . . . 1 1.1 Ca vitation . . . . . . . . . . . . . . . . . 1 1.2 T ypes of Superca vitating Projectiles . . . . . . . . . . 2 1.3 Related Research . . . . . . . . . . . . . . . 4 1.4 Ov ervie w of this Thesis . . . . . . . . . . . . . . 5 2 CONFIGURA TION OF VEHICLE . . . . . . . . . 7 2.1 Ca vitator . . . . . . . . . . . . . . . . . 7 2.2 Fins . . . . . . . . . . . . . . . . . . . 8 2.3 Maneuv ering . . . . . . . . . . . . . . . . 8 3 NONLINEAR EQ U A TIONS OF MO TION . . . . . . . 10 3.1 Kinematic Equations of Motion . . . . . . . . . . . . 10 3.1.1 Orientation of the T orpedo . . . . . . . . . . . 10 3.1.2 Orientation of the Ca vitator . . . . . . . . . . 13 3.1.3 Orientation of Fins . . . . . . . . . . . . . 14 3.1.4 Angle of Attack and Sideslip . . . . . . . . . . 18 3.1.5 Kinematic Equations . . . . . . . . . . . . 20 3.2 Dynamic Equations of Motion . . . . . . . . . . . . 23 3.2.1 F orces on Ca vitator . . . . . . . . . . . . . 25 3.2.2 F orces on Fins . . . . . . . . . . . . . . 27 3.2.3 Gra vitational F orces . . . . . . . . . . . . 31 3.2.4 Equations of Motion . . . . . . . . . . . . 31 4 LINEARIZA TION . . . . . . . . . . . . 34 4.1 Linearization . . . . . . . . . . . . . . . . 34 iii

PAGE 4

4.1.1 Need for Linearization . . . . . . . . . . . . 34 4.1.2 Generic F orm of Equations of Motion . . . . . . . . 35 4.1.3 Small Disturbance Theory . . . . . . . . . . . 35 4.1.4 Stability and Control Deri v ati v es . . . . . . . . . 37 4.2 State Space Representation . . . . . . . . . . . . . 42 5 CONTR OL DESIGN SETUP . . . . . . . . . . 46 5.1 Open-loop Performance . . . . . . . . . . . . . . 47 5.2 Closed-Loop Problem . . . . . . . . . . . . . . 50 5.3 Rob ustness of the Controller . . . . . . . . . . . . 52 5.3.1 Gain Mar gin . . . . . . . . . . . . . . 52 5.3.2 Phase Mar gin . . . . . . . . . . . . . . 52 5.3.3 Uncertainty In P arameters . . . . . . . . . . . 53 5.3.4 Controller Objecti v e . . . . . . . . . . . . 53 6 LQR CONTR OL . . . . . . . . . . . . 54 6.1 LQR Theory . . . . . . . . . . . . . . . . . 54 6.2 Control Synthesis . . . . . . . . . . . . . . . 58 6.3 Nominal Closed-loop Model . . . . . . . . . . . . 60 6.3.1 Model . . . . . . . . . . . . . . . . 60 6.3.2 Linear Simulations . . . . . . . . . . . . . 60 6.3.3 Gain and Phase Mar gins . . . . . . . . . . . 63 6.4 Perturbed Closed-loop Model . . . . . . . . . . . . 64 6.4.1 Model . . . . . . . . . . . . . . . . 66 6.4.2 Linear Simulations . . . . . . . . . . . . . 69 6.4.3 Gain and Phase Mar gins . . . . . . . . . . . 72 7 NONLINEAR SIMULA TIONS . . . . . . . . . 74 7.1 Nonlinear Simulations for Nominal System . . . . . . . . 74 7.2 Nonlinear Simulations for Perturbed System . . . . . . . . 79 8 CONCLUSION . . . . . . . . . . . . . 83 8.1 Summary . . . . . . . . . . . . . . . . . 83 8.2 Future W ork . . . . . . . . . . . . . . . . . 83 APPENDIX A REFERENCE FRAMES AND R O T A TION MA TRICES . . . . 84 B NUMERICAL TECHNIQ UES . . . . . . . . . . 86 B.1 Interpolation of F orce Data . . . . . . . . . . . . . 86 B.1.1 Extrapolation Scheme . . . . . . . . . . . . 86 i v

PAGE 5

B.1.2 Ca vitator . . . . . . . . . . . . . . . 87 B.1.3 Fins . . . . . . . . . . . . . . . . . 87 B.2 Numerical Linearization . . . . . . . . . . . . . 88 REFERENCES . . . . . . . . . . . . . 90 BIOGRAPHICAL SKETCH . . . . . . . . . . . 91 v

PAGE 6

LIST OF T ABLES T able page 5.1 Control P arameters . . . . . . . . . . . . . . . 46 5.2 Control Constraints . . . . . . . . . . . . . . . 51 6.1 Gain and Phase Mar gin with LQR Controller . . . . . . . . 64 6.2 Percentage V ariation in A Matrix due to 20% V ariation in cl c . . . . 67 6.3 Percentage V ariation in B Matrix due to 10% V ariation in cl c . . . . 67 6.4 Percentage V ariation in A Matrix due to 20% V ariation in cd c . . . . 67 6.5 Percentage V ariation in B Matrix due to 20% V ariation in cd c . . . . 68 6.6 Percentage V ariation in A Matrix due to 20% V ariation in cl f in . . . 68 6.7 Percentage V ariation in B Matrix due to 20% V ariation in cl f in . . . 68 6.8 Percentage V ariation in A Matrix due to 20% V ariation in cd f in . . . 69 6.9 Percentage V ariation in B Matrix due to 20% V ariation in cd f in . . . 70 6.10 Gain and Phase Mar gin for Perturbed Closed-loop System: 20% error in cl f in 73 B.1 Grid F or Experimental Ca vitator Data . . . . . . . . . . 88 B.2 Grid F or Experimental Fin Data . . . . . . . . . . . 88 vi

PAGE 7

LIST OF FIGURES Figure page 1.1 T ip V orte x Ca vitation . . . . . . . . . . . . . . 2 1.2 F ormation of Ca vity . . . . . . . . . . . . . . . 3 1.3 Superca vitating V ehicle . . . . . . . . . . . . . . 4 2.1 Superca vitating V ehicle . . . . . . . . . . . . . . 7 2.2 Ca vitator and Fins . . . . . . . . . . . . . . . 9 3.1 Body-x ed and Inertial Frames . . . . . . . . . . . . 11 3.2 Principle Planes of Symmetry for the T orpedo . . . . . . . . 12 3.3 Euler Angles of Rotation . . . . . . . . . . . . . 12 3.4 Ca vitator Reference Frame . . . . . . . . . . . . . 13 3.5 Rudder and Fin Reference Frames . . . . . . . . . . . 14 3.6 Rudder 1 Fin Reference Frames . . . . . . . . . . . 16 3.7 Rudder 2 Fin Reference Frames . . . . . . . . . . . 16 3.8 Ele v ator 1 Fin Reference Frames . . . . . . . . . . . 17 3.9 Ele v ator 2 Fin Reference Frames . . . . . . . . . . . 17 3.10 Angle of Attack ( a ) and Sideslip ( b ) . . . . . . . . . . 18 3.11 Ca vitator: (a) Angle of Attack and Sideslip and (b) Hydrodynamic F orces 25 3.12 Fin Geometry . . . . . . . . . . . . . . . . 28 5.1 Simulink Model for Open Loop Simulation . . . . . . . . 48 5.2 Open-Loop Response for T orpedo: wpq . . . . . . . . . 48 5.3 V ariation of Eigen v alues with Change in V elocity . . . . . . . 50 5.4 Loop Gain . . . . . . . . . . . . . . . . . 53 6.1 Controller for T racking when Plant has an Inte grator . . . . . . 57 6.2 Controller for T racking when Plant has no Inte grator . . . . . . 57 vii

PAGE 8

6.3 Eigen v alues for the Closed-loop System . . . . . . . . . 60 6.4 Pitch Command T racking for Linear System : qu . . . . . . . 61 6.5 Pitch Command T racking for Linear System : d c d c . . . . . . 61 6.6 Pitch Command T racking for Linear System : d e 1 d e 1 . . . . . . 62 6.7 Roll Command T racking for Linear System : pr . . . . . . . 62 6.8 Roll Command T racking for Linear System : d r 1 d r 1 . . . . . . 62 6.9 Breakpoints for Calculating the Loop-Gain for a T racking Controller . . 63 6.10 Gain and Phase Mar gin: Longitudinal Outer -loop . . . . . . . 64 6.11 Gain and Phase Mar gin: Longitudinal Inner -loop . . . . . . . 65 6.12 Gain and Phase Mar gin: Longitudinal All-loop . . . . . . . 65 6.13 Gain and Phase Mar gin: Lateral Outer -loop . . . . . . . . 66 6.14 Gain and Phase Mar gin: Lateral Inner -loop . . . . . . . . 66 6.15 Gain and Phase Mar gin: Lateral All-loop . . . . . . . . . 69 6.16 Eigen v alues for the Perturbed Closed-loop System: 20% Error in cl f in . 70 6.17 Pitch Command T racking for Perturbed Linear System : qu . . . . 71 6.18 Pitch Command T racking for Perturbed Linear System : d c d c . . . 71 6.19 Pitch Command T racking for Perturbed Linear System : d e 1 d e 1 . . . 71 6.20 Roll Command T racking for Perturbed Linear System : pr . . . . 72 6.21 Roll Command T racking for Perturbed Linear System : d r 1 d r 1 . . . 72 7.1 Complete Nonlinear Simulation with LQR Controller . . . . . . 75 7.2 Pitch Command T racking : pq . . . . . . . . . . . . 76 7.3 Pitch Command T racking : d c d c . . . . . . . . . . . 76 7.4 Pitch Command T racking : d e 1 d e 1 . . . . . . . . . . . 76 7.5 Pitch Command T racking : d r 1 xyzT rajectory . . . . . . . 77 7.6 Roll Command T racking: pq . . . . . . . . . . . . 77 7.7 Roll Command T racking: d c d c . . . . . . . . . . . . 78 7.8 Roll Command T racking: d e 1 d e 1 . . . . . . . . . . . 78 viii

PAGE 9

7.9 Roll Command T racking: d r 1 d r 1 . . . . . . . . . . . 78 7.10 Roll Command T racking:xyzT rajectory . . . . . . . . 79 7.11 Roll & Pitch Command T racking: pq . . . . . . . . . . 79 7.12 Roll & Pitch Command T racking: d c d c . . . . . . . . . 80 7.13 Roll & Pitch Command T racking: d e 1 d e 1 . . . . . . . . . 80 7.14 Roll & Pitch Command T racking: d r 1 d r 1 . . . . . . . . . 80 7.15 Roll & Pitch Command T racking:xyzT racking . . . . . . 81 7.16 Response for 20% V ariation in cl f in : uw . . . . . . . . . 81 7.17 Response for 20% V ariation in cl f in : pq . . . . . . . . . 82 7.18 Response for 20% V ariation in cl f in : d cd r 1 . . . . . . . . 82 7.19 Response for 20% V ariation in cl f in :xyzT rajectory . . . . . 82 A.1 Rotation of Frames . . . . . . . . . . . . . . . 84 B.1 Shape Function for One Dimensional Quadratic Scheme . . . . . 87 ix

PAGE 10

Abstract of Thesis Presented to the Graduate School of the Uni v ersity of Florida in P artial Fulllment of the Requirements for the De gree of Master of Science CONTR OL STRA TEGIES FOR SUPERCA VIT A TING VEHICLES By Anukul Goel December 2002 Chair: Andre w J. K urdila Cochair: Richard C. Lind Major Department: Mechanical and Aerospace Engineering Underw ater tra v el is greatly limited by the speed that can be attained by the v ehicles. Usually the maximum speed achie v ed by underw ater v ehicles is about 40 m/s. Superca vitation can be vie wed as a phenomenon that w ould help us to break the speed barrier in underw ater v ehicles. The idea is to mak e the v ehicle surrounded by w ater v apor while it is tra v eling underw ater Thus, the v ehicle actually tra v els in air and has v ery small skin friction drag. This allo ws it to attain high speeds underw ater Superca vitation is a phenomenon which is observ ed in w ater As the relati v e v elocity of w ater with respect to the v ehicle increases, the pressure decreases and subsequently it e v aporates to form w ater v apor Superca vitation has its dra wbacks. It is really hard to control and maneuv er a superca vitating v ehicle. The rst part of this w ork deals with modeling of a superca vitating torpedo. Nonlinear equations of motion are deri v ed in detail. The latter part of the w ork deals with nding a controller to maneuv er the torpedo. A controller is obtained by LQR synthesis. F or the synthesis, it is assumed that the ca vity is x ed and the torpedo is situated symmetrically in the ca vity The rob ustness analysis of the LQR controllers is carried out x

PAGE 11

by calculating the gain and phase mar gins. Simulations of the response for a perturbed system also ha v e been studied. xi

PAGE 12

CHAPTER 1 INTR ODUCTION Achie ving high speeds is an important issue for underw ater v ehicles. Ev en the common f astest underw ater v ehicles are restricted to tra v el at speeds around 40 ms1 The reason for this restriction is the drag induced by skin friction. When a body mo v es in a uid, a layer of the uid clings to the surf ace 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 w ater unlik e air is dominated by skin friction drag as compared to other sources such as pressure drag. Ov er the years, e xtensi v e research has been done to boost the speed of underw ater v ehicles. Ho we v er most of the studies were mainly focused on streamlining the bodies and impro ving their propulsion systems. Ev en though these ha v e pro v en to gi v e adv ancements in speed, there has not been a considerable reduction in skin friction drag. In the late 1970' s, the Russian Na vy w as able to engineer a torpedo, the Shkv al (squall) [ 1 ], that achie v ed a speed o v er 80 ms1 This phenomenal impro v ement in speed w as made possible by superca vitation. The idea w as to en v elop the torpedo in a gas so that it has little contact with w ater The Shkv al w as able to achie v e a tremendous reduction in skin friction drag and e xhibit high speed. 1.1 Ca vitation As the speed of an underw ater v ehicles increases, i.e., as the relati v e v elocity of w ater with respect to the v ehicle increases, the pressure decreases. The speed can be increased to the point the pressure goes belo w the v apor pressure of w ater and then b ubbles of w ater v apor are formed. This process is kno wn as ca vitation. Ca vitation is mostly observ ed at sharp corners of a body where the speed can reach high magnitudes. A classic e xample for ca vitation is at the tip of propellers, lik e the one sho wn in Figure 1.1 Since the propeller 1

PAGE 13

2 is rotating at high speed, the relati v e v elocity at the tips is high enough so that w ater at the tips v aporizes. Ca vitation has been e xtensi v ely researched, b ut remains a challenge for underw ater v ehicles. Figure 1.1 T ip V orte x Ca vitation [ 2 ] Ca vitation can be harmful in man y cases. The ca vitation re gion is usually turb ulent. When the ca vitation is not steady the pressure drop is momentary and v ery quickly the ca vitation b ubble encounters a re gion of high pressure that forces the b ubble to collapse lik e a tin y bomb This collapse causes high le v els of noise and also may cause considerable damage to the surf ace of the body Figure 1.2 sho ws the v arious stages of formation of ca vity It sho ws a b ullet-lik e body tra v eling in w ater at high speed. The ca vitation starts as v apor b ubbles near the re gion of high relati v e v elocity As the speed is further increased, the b ubbles mer ge to form a lar ge area of v apor On further increase in speed, the whole of the body is co v ered in v apor This stage is called the superca vitation. At this point, the condition is similar to tra v eling in air The skin friction drag is tremendously reduced, and thus high speed can be attained in this phase. 1.2 T ypes of Super ca vitating Pr ojectiles The rst v ersion of the Russian Shkv al w as a crude superca vitating v ehicle. It w as unguided and had a small range of about 5 miles. No w there are v arious adv anced forms of superca vitating bodies. One class of superca vitating bodies, called RAMICS (Rapid

PAGE 14

3 Figure 1.2 F ormation of Ca vity [ 3 ] Airborne Mine Clearance System), is a helicopter -borne weapon that destro ys surf ace and near -surf ace marine mines by ring superca vitating rounds at them. These are small b ulletlik e, at nosed projectiles that are designed to tra v el stably through air and w ater As the RAMICS can produce a ca vitation b ubble, the y can tra v el at high speed underw ater and can pierce the mines to destro y them. As the y are red from a distance in air the y need to be designed to tra v el in both phases. The RAMICS are better than con v entional b ullets, as con v entional b ullets are rapidly slo wed do wn by drag in w ater Another kind of a superca vitating projectile is a sub-surf ace gun system using Adaptable High-Speed Undersea Munitions (AHSUM). These are superca vitating “kinetic-kill” b ullets, red from guns tted under submer ged hull of submarines. These sonar directed b ullets w ould be used to intercept undersea missiles. The RAMICS and AHSUM are uncontrolled small range superca vitating projectiles. The ne xt le v el of superca vitating projectiles is lar ger torpedoes, with higher speeds and longer range. These v ehicles are much more comple x. The y require the design of a launch station. The y require detailed studies of hydrodynamics, acoustics, guidance and control, propulsion, etc. Ev en if the v ehicle is designed to be an uncontrolled torpedo, it is still a challenge to produce and maintain a ca vity around the v ehicle. The ca vity is usually produced by the nose of the v ehicle, which is specially shaped for the purpose. The nose is called a ca vitator The ca vitator may not be suf cient to produce the ca vity Thus air can be blo wn at the nose and v arious sections of the body so as to maintain and produce the

PAGE 15

4 ca vity Figure 1.3 sho ws a superca vitating torpedo tra v eling underw ater It can be seen that the whole of its body is en v eloped by a ca vity This is the kind of v ehicle that has been studied in this w ork. Cavity Fin Cavitator Figure 1.3 Superca vitating V ehicle [ 1 ] 1.3 Related Resear ch Research in the eld of superca vitation has been going on from the early 1900' s. But earlier research w as aimed at reduction of ca vitation so as to reduce noise and body damage. In the 90' s the focus shifted to e xploiting the ef fects of superca vitation. The w ork sho wn in May [ 4 ] has an e xtensi v e collection of e xperimental data for v ariation of forces on v arious superca vitating shapes. Coef cients of lift and drag are plotted with the v ariation of ca vitation number for shapes lik e disks, cones, ogi v es and wedges. The w ork done in this research mak es use of a CFD database pro vided in Fine [ 5 ]. This database has v alues for coef cients of lift and drag for conical ca vitators, which are functions of the half angle of the cone and the angle of attack. This database also has coef cients of lift and drag for wedges, which are functions of wetted surf ace of the wedge, angle of attack and sweepback angle (we discuss the denition of these geometric quantities such as the angle of attack and half angle later in this thesis). This information is useful to calculate the forces on ns of the torpedo. In late 90' s a lot of research has been done on the dynamics of superca vitating v ehicles. W ork sho wn in K ulkarni and Pratap [ 6 ] and Rand et al. [ 7 ] deals with studying dynamics of uncontrolled superca vitating projectiles. A dynamic model for RAMICS and AHSUM has been de v eloped. It is sho wn that these projectiles rotate inside the ca vity This rotation leads to impacts between the tail of the projectile and the ca vity w all. The frequenc y of the

PAGE 16

5 impact increases with time. These projectiles are v ery short range and ha v e a small time of ight, on the order of a fe w seconds. The w ork sho wn in Dzielski and K urdila [ 8 ] focuses on the formulation of a control problem for a superca vitating torpedo. A dynamical model for a n controlled torpedo has lik e wise been de v eloped. The model also includes a formulation for the ca vity It is observ ed that the weight of the body forces it to skip inside the w alls of the ca vity and the v ehicle is unstable. A control system is designed for the torpedo and results of closed-loop simulations ha v e been presented. 1.4 Ov er view of this Thesis This w ork aims at formulating the control design for a superca vitating torpedo. Equations of motion for the torpedo are deri v ed, and linear control methodologies are applied for pitch and roll control of the torpedo. The main purpose of this w ork is to analyze these controllers and ascertain their rob ustness to v arious errors and uncertainties. Chapter 2 of this thesis briey describes the conguration of the superca vitating torpedo for which the equations of motion ha v e been de v eloped. A detailed deri v ation of the equations of motion for the torpedo has been carried out in Chapter 3 V arious reference frames ha v e been used to obtain the kinematic equations of the torpedo. Dynamic equations are deri v ed using Ne wton' s La ws. V arious forces e xperienced by dif ferent re gions of the torpedo ha v e been e xplained. Chapter 4 describes linearization of the equations of motion using small disturbance theory Numerical methods are used for this purpose. It is observ ed that the linearization for a simple trim can be v ery complicated. Chapter 5 describes the control synthesis for the torpedo. Open-loop dynamics are sho wn. The closed-loop problem and v arious constraints on the torpedo ha v e been stated. Chapter 6 formulates a Linear Quadratic Re gulator (LQR) control design for the tor pedo. Controllers for pitch and roll rate control of the torpedo are obtained. The results for linear closed-loop system and a perturbed liner system ha v e been sho wn.

PAGE 17

6 The results for pitch and roll rate control for the nonlinear closed-loop system ha v e been presented in Chapter 7 The simulations for a perturbed system model also ha v e been presented.

PAGE 18

CHAPTER 2 CONFIGURA TION OF VEHICLE Although superca vitation can be a v ery helpful phenomenon, it presents signicant challenges in modeling and control of superca vitating v ehicles. As a signicant portion of the v ehicle is located in the ca vity the control, guidance and stability of the v ehicle has to be managed by v ery small re gions in front and aft of the v ehicle. Also generation of the ca vity can be a signicant problem. The major problems associated with the superca vitating v ehicles may be summarized as:generation and maintenance of ca vitybalancing the weight of the v ehiclecontrol and guidancestability Figure 2.1 is a gure of a superca vitating torpedo that is modeled in this w ork. The main parts of the torpedo are the ca vitator in the front and the four ns in the aft portion. The ca vitator is used to generate and maintain the ca vity The Ca vitator and the four ns together are also used for control and stability of the v ehicle. 2.1 Ca vitator The ca vitator is the element that generates a ca vity around the torpedo. Se v eral types of ca vitators, including cones, wedges and plates ha v e been in v estigated [ 4 ]. This project will consider a conical ca vitator as sho wn in Figure 2.1 The main parameter that denes Figure 2.1 Superca vitating V ehicle [ 8 ] 7

PAGE 19

8 a conical ca vitator is its half-angle. The coef cients of lift and drag, as functions of halfangle, for the ca vitator ha v e been generated using CFD and tab ulated in [ 5 ]. The ca vitator in this model has one de gree of freedom dened by a rotation angle about an axis perpendicular to its axis of symmetry The shape and location of the ca vity is a nonlinear function of this ca vitator deection angle and half angle of the cone. This shape determines the position where the ca vity intersects the body of the v ehicle. Thus, the amount of wetted area of the v ehicle is dependent on these angles, which in turn determines the ef cienc y of superca vitation achie v ed. As stated earlier a lar ge portion of the v ehicle is in the ca vity The lift produced by the ca vitator combined with the lift produced by the ns helps in balancing the weight of the body 2.2 Fins Although the ca vitator is capable of pro viding enough lift to sustain the body in w ater it does not usually pro vide suf cient forces and moments to stabilize and control the v ehicle. F or this purpose ns are required in the aft portion of the v ehicle. The ns help in counteracting the moments produced by the ca vitator and also pro viding some amount of lift to support the weight of the body There are four ns placed symmetrically along the girth of the v ehicle, near the tail. The ns also are the control surf aces, as the y can pro vide dif ferential forces, thus causing control moments. T w o of the ns sho wn in Figure 2.2 are horizontal (placed parallel to the axis of rotation of ca vitator), called ele v ators, are used to af fect the longitudinal dynamics for the v ehicle. The other tw o ns are called the rudders and are used to af fect the lateral dynamics to the v ehicle. The ns ha v e tw o de grees of freedom, a sweepback angle and an angle of rotation. These angles will be e xplained further in Chapter 3. 2.3 Maneuv ering The motion of the v ehicle is controlled by all v e control surf aces, the four ns and the ca vitator In a straight line motion the ca vitator and the ele v ators balance the weight of the

PAGE 20

9 Figure 2.2 Ca vitator and Fins v ehicle. A propulsion force at the tail pushes the v ehicle forw ard. The rudders usually ha v e a zero deection in such a case. The v ehicle depends on a bank-to-turn strate gy for making a turn, and cannot use traditional missile strate gies such as skid-to-turn. This dependence arises because the hull of the v ehicle is incapable of pro viding a lift force. The ns are the main lift generating surf aces. A dif ferential force from the ns can be used to generate a force to w ards the center of the turn.

PAGE 21

CHAPTER 3 NONLINEAR EQ U A TIONS OF MO TION The dynamics of the torpedo, whose conguration w as discussed in Chapter 2, can be mathematically formulated. A complete deri v ation of the equations of motion is gi v en in this chapter Section 3.1 deals with the setup of reference frames and deri v ation of the kinematic equations. The dynamic equations of motion are deri v ed in Section 3.2. 3.1 Kinematic Equations of Motion The denition of a suitable coordinate system and de grees of freedom is a precursor to an y formulation of equations of motion. The deri v ation of the equations of motion of multi-body systems is highly simplied by dening v arious reference frames and relations between them. Appendix A describes briey the concept of reference frames and rotation matrices. These concepts will be used e xtensi v ely in the deri v ation of equations of motion. The deri v ation of the equations of motion will be done in tw o parts. First, the kinematic equations will be deri v ed. These include the formulation of the position v ectors, v elocities and accelerations of v arious points on the torpedo. Ne xt, the dynamic equations will be deri v ed. Finally Ne wton' s La ws yield the dynamic equations of motion. 3.1.1 Orientation of the T or pedo Six de grees 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 gi v e the location of a point x ed on the torpedo. The other 3 DOFs gi v e the orientation of the torpedo in a x ed reference frame. The position of the torpedo in a reference frame can also be obtained by the inte gration of its v elocity in that reference frame. The torpedo is assumed to be mo ving in an earth-x ed reference frame E centered at an y con v eniently chosen point and described by the basis v ectorsˆ e 1ˆ e 2ˆ e 3. The earth-x ed 10

PAGE 22

11 ˆ e 1 ˆ e 3ˆ e2Oˆ b 1 ˆ b 2 ˆ b 3Figure 3.1 Body-x ed and Inertial Frames reference frame is sho wn in Figure 3.1 The v ector ˆ e 3 points in the do wnw ard direction, i.e., the direction of the gra vity The v ectors ˆ e 1 and ˆ e 2 are placed in the horizontal plane such that the basis v ectors form a right-handed coordinate system. As sho wn in the gure, ˆ e 1 points to the right and ˆ e 2 points outside the plane of the paper A body-x ed frame B is dened by the basis v ectorsˆ b 1ˆ b 2ˆ b 3so as to simplify the deri v ation of the equations of motion. The frame B is centered at O the center of gra vity of the torpedo, and mo v es with the torpedo. The reference frame B is sho wn in the Figure 3.1 It can be seen in Figure 3.2 that the torpedo has tw o planes of symmetry namely ˆ b 1 ˆ b 2 and ˆ b 1 ˆ b 3 The plane ˆ b 1 ˆ b 3 is called the longitudinal plane and plane ˆ b 1 ˆ b 2 the lateral plane. T ransformation from E to B can be achie v ed by 3 rotations. Man y such sets of rotations are possible. The rotation steps chosen here are the Euler angles of rotation, which are the ya w pitch and roll. As there are three rotations to be performed, tw o intermediate reference frames with basis v ectorsˆ x 1ˆ x 2ˆ x 3andˆ y 1ˆ y 2ˆ y 3will ha v e to be introduced to perform the transformation. The rotations, to be performed in order are listed belo w 1. Rotate the frame E about ˆ e 3 through a ya w angle, Y to obtain the frameˆ x 1ˆ x 2ˆ x 3. 2. Rotateˆ x 1ˆ x 2ˆ x 3about ˆ x 2 through a pitch angle, Q to obtain the frameˆ y 1ˆ y 2ˆ y 3. 3. Rotateˆ y 1ˆ y 2ˆ y 3about ˆ y 1 through a roll angle, F to obtain the frame B

PAGE 23

12 ˆ b 1 ˆ b 3 ˆ b 2 B Figure 3.2 Principle Planes of Symmetry for the T orpedo ˆ x 2 ˆ e 3ˆ x 3ˆ x1Y Y ˆ e 1ˆ e2ˆ x 3ˆ x 1Q Q ˆ y 3 ˆ y 1ˆ x 2ˆ y 2ˆ y 3 ˆ y 1ˆ b 1ˆ y 2F F ˆ b 3 ˆ b 2 Figure 3.3 Euler Angles of Rotation Figure 3.3 sho ws the abo v e rotations in order Based on the abo v e denition of rotations, the transformation matrix from E to B can be written as in equation 3.1 ˆ b 1 ˆ b 2 ˆ b 3 n r 1 0 0 0 C F S F 0S F C F r C Q 0S Q 0 1 0 S Q 0 C Q r C Y S Y 0S Y C Y 0 0 0 1 ˆ e 1 ˆ e 2 ˆ e 3 n r C Q C Y C Q S YS Q C Y S F S QC F S Y S F S Q S YC Y C F S F C Q C F S Q C YS F S Y C F S Q S YC Y S F C F C Q ˆ e 1 ˆ e 2 ˆ e 3 n E B ˆ e 1 ˆ e 2 ˆ e 3 n(3.1)

PAGE 24

13 3.1.2 Orientation of the Ca vitator As described earlier the ca vitator has only one de gree of freedom. It can rotate in the longitudinal plane about an axis parallel to the ˆ b 2 axis. The orientation of the ca vitator -x ed ax es with respect to the body x ed ax es is sho wn in Figure 3.4 B A ˆ b 1 ˆ b 3ˆ c1ˆ c3C d c ˆ b 3 ˆ b 1ˆ c 1P CP D C P A d c Figure 3.4 Ca vitator Reference Frame The deection of the ca vitator is gi v en by an angle, d c which is positi v e for a positi v e ca vitator rotation about the ˆ b 2 direction. The geometric center of rotation of ca vitator is denoted by P CP is the center of pressure for the ca vitator which is at a distance D C P from P along ˆ c 1 From Figure 3.4 the rotation matrix from B to ca vitator x ed frame C can be written as in Equation 3.2 ˆ c 1 ˆ c 2 ˆ c 3 n r C d c 0S d c 0 1 0 S d c 0 C d c ˆ b 1 ˆ b 2 ˆ b 3 n(3.2)

PAGE 25

14 3.1.3 Orientation of Fins Figure 3.5 sho ws the orientation of the n-x ed reference frames. F or con v enience, all the n frames are indicated by basis v ectorsˆ f 1ˆ f 2ˆ f 3. In te xt the y will be represented asˆ f 1ˆ f 2ˆ f 3f in where subscript f in corresponds to a particular n. Rudder 1 Rudder 2 Elevator 1 Elevator 2 FRONT VIEW TOP VIEW B ˆ f 1 ˆ f 1 ˆ b 3 ˆ b 1 ˆ b 2 ˆ f 2 ˆ f 3 ˆ f 3 ˆ f 2 B ˆ f 1 ˆ f 1 ˆ b 1 ˆ f 2 ˆ f 3 ˆ f 3 ˆ f 2 ˆ b 2 ˆ b 3 Figure 3.5 Rudder and Fin Reference Frames All the ns ha v e 2 DOFs, namely the sweepback angle ( q f in ) and the n rotation ( d f in ). These can be dened asSweepback angle ( q f in ): This parameter is the angle of rotation of a n about its ˆ f 3 axis. The direction of rotation for positi v e sweepback is such that the n is mo v ed to w ard the rear portion of the torpedo. Due to this denition, the positi v e sweepback is along the ne gati v e ˆ f 3 direction for all ns. Sweepback angle determines the amount of n that is en v eloped in the ca vity

PAGE 26

15Fin Rotation ( d f in ): This parameter is the angle of rotation of the n about its ˆ f 2 axis, in positi v e the ˆ f 2 direction. Fin rotation determines the magnitude and direction of the forces acting on the ns. The order of rotation in the abo v e case is important to obtain the correct equations. Sweepback has to be performed before n rotation. An intermediate reference frame G with basis v ectorsˆ g 1ˆ g 2ˆ g 3is introduced so as to simplify the deri v ation of rotation matrix from B to the n coordinates. Sweepback aligns the n-x ed frames with the intermediate frame G and then a n rotation puts the n frame in actual orientation with the ns. As the second rotation is identical in all cases, the transformation matrix from frame G to n frame F f in can be written as in Equation 3.3 ˆ f 1 ˆ f 2 ˆ f 3 nf in r C d f in 0S d f in 0 1 0 S d f in 0 C d f in ˆ g 1 ˆ g 2 ˆ g 3 nf in (3.3) The rotation matrix for sweepback and the transformation matrices from B to F f in frame for each of the ns can be deri v ed easily and are summarized belo w .Rudder 1 Figure 3.6 sho ws the sweepback and n rotation for rudder 1. The matrices for transformation from B to R 1 can be deri v ed as in Equation 3.4 and Equation 3.5 ˆ g 1 ˆ g 2 ˆ g 3 nR 1 r C q R 1 0 S q R 1S q R 1 0C q R 1 01 0 ˆ b 1 ˆ b 2 ˆ b 3 n(3.4) ˆ f 1 ˆ f 2 ˆ f 3 nR 1 r C d R 1 0S d R 1 0 1 0 S d R 1 0 C d R 1 r C q R 1 0 S q R 1S q R 1 0C q R 1 01 0 ˆ b 1 ˆ b 2 ˆ b 3 n(3.5)Rudder 2 Figure 3.7 sho ws the sweepback and n rotation for rudder 2. The matrices for transformation from B to R 2 can be deri v ed as in Equation 3.6 and Equation 3.7 ˆ g 1 ˆ g 2 ˆ g 3 nR 2 r C q R 2 0S q R 2S q R 2 0 C q R 2 0 1 0 ˆ b 1 ˆ b 2 ˆ b 3 n(3.6) ˆ f 1 ˆ f 2 ˆ f 3 nR 2 r C d R 2 0S d R 2 0 1 0 S d R 2 0 C d R 2 r C q R 2 0S q R 2S q R 2 0 C q R 2 01 0 ˆ b 1 ˆ b 2 ˆ b 3 n(3.7)

PAGE 27

16 ˆ g 3 ˆ g 1 ˆ g 2ˆ f 2 ˆ f 3 ˆ f 1 d R 1 q R 1 q R 1 q R 1 ˆ b 1 ˆ b 3ˆ b 2 ˆ g 1 ˆ g 2 ˆ g 3Figure 3.6 Rudder 1 Fin Reference Frames ˆ b 1 ˆ b 3 ˆ g 2ˆ b 2 ˆ g 1 q R 2 q R 2 ˆ g 1 ˆ g 2 ˆ f 1 ˆ f 2 d R 2 d R 2 ˆ g 3 q R 2 ˆ g 2 Figure 3.7 Rudder 2 Fin Reference FramesEle v ator 1 Figure 3.8 sho ws the sweepback and n rotation for Ele v ator 1 The matrices for transformation from B to E 1 can be deri v ed as in Equation 3.8 and Equation 3.9 ˆ g 1 ˆ g 2 ˆ g 3 nE 1 r C q E 1S q E 1 0S q E 1 C q E 1 0 0 01 ˆ b 1 ˆ b 2 ˆ b 3 n(3.8) ˆ f 1 ˆ f 2 ˆ f 3 nE 1 r C d E 1 0S d E 1 0 1 0 S d E 1 0 C d E 1 r C q E 1S q E 1 0S q E 1 C q E 1 0 0 01 ˆ b 1 ˆ b 2 ˆ b 3 n(3.9)Ele v ator 2 Figure 3.9 sho ws the sweepback and n rotation for Ele v ator 2 The matrices for transformation from B to E 2 can be deri v ed as in Equation 3.10 and

PAGE 28

17 ˆ g 3 ˆ g 1 ˆ g 2ˆ f 2 ˆ f 3 ˆ f 1 d E 1 q R 1 ˆ b 1 ˆ g 1 ˆ g 2 ˆ g 3 ˆ b 3 ˆ b 2 q E 1 q E 1 Figure 3.8 Ele v ator 1 Fin Reference Frames ˆ b 1 ˆ g 1 ˆ g 1 ˆ g 2 ˆ f 1 ˆ f 2 d R 2 ˆ g 3 ˆ g 2 ˆ b 2 q E 2 d E 2 q E 2 ˆ g 3ˆ b 3 q E 2 Figure 3.9 Ele v ator 2 Fin Reference Frames Equation 3.11 ˆ g 1 ˆ g 2 ˆ g 3 nE 2 r C q E 2 S q E 2 0S q E 2C q E 2 0 0 0 1 ˆ b 1 ˆ b 2 ˆ b 3 n(3.10) ˆ f 1 ˆ f 2 ˆ f 3 nE 2 r C d E 2 0S d E 2 0 1 0 S d E 2 0 C d E 2 r C q E 2 S q E 2 0S q E 2C q E 2 0 0 0 1 ˆ b 1 ˆ b 2 ˆ b 3 n(3.11) Equations 3.2 to 3.11 will be used in later sections to transform the forces on ns and ca vitator to the body-x ed frame.

PAGE 29

18 3.1.4 Angle of Attack and Sideslip The body-x ed reference frame has been dened, b ut the v elocity of v arious points on the body of the torpedo is yet to be dened. The torpedo is considered as a rigid body If the v elocity of a certain point on a rigid body is kno wn, the v elocity at an y other point on that body can be found by kno wing the rotation matrices. Thus, V u ˆ b 1v ˆ b 2w ˆ b 3 will be tak en as the v elocity of CG of the torpedo. The v elocity of other points on the torpedo can be dened subsequently T w o v ery useful parameters, angle of attack and angle of sideslip can be dened in conjunction with the orientation of the body axis with the v elocity of CG Figure 3.10 sho ws these parameters and their geometric interpretation. ˆ f 1 u wv ˆ b 1 ˆ b 2 ˆ b 3 a b V u ˆ b 1 v ˆ b 2 w ˆ b 3ˆ g1Figure 3.10 Angle of Attack ( a ) and Sideslip ( b ) Angle of attack, a can be dened as the angle between the projection of v elocity V onto ˆ b 2 ˆ b 3 plane and the ˆ b 1 axis. Angle of attack is positi v e when the nose of the torpedo points abo v e the v elocity v ector of the torpedo. As before, an intermediate frame F gi v en byˆ f 1ˆ f 2ˆ f 3can be described by rotation of the B frame by an angle a Thus the rotation matrix from F bod y to B can be written. ˆ b 1 ˆ b 2 ˆ b 3 n r C a 0S a 0 1 0 S a 0 C a ˆ f 1 ˆ f 2 ˆ f 3 nbod y (3.12) The Angle of sideslip, b is dened as the angle between the actual v elocity V and the projection of V onto ˆ b 2 ˆ b 3 plane. Again, a frame G bod y can be dened by rotation of F bod y

PAGE 30

19 by an angle equal to b in ne gati v e ˆ f 3 direction, thus gi ving a rotation matrix as in Equation 3.13 ˆ g 1 ˆ g 2 ˆ g 3 nbod y r C bS b 0 S b C b 0 0 0 1 ˆ f 1 ˆ f 2 ˆ f 3 nbod y (3.13) V elocity of CG of the torpedo in the G bod y frame can be written as V ˆ g 1 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 G bod y to B is important. It is gi v en by Equation 3.14 ˆ b 1 ˆ b 2 ˆ b 3 n r C b C a C a S bS aS b C b 0 C b S a S a S b C a ˆ g 1 ˆ g 2 ˆ g 3 nbod y (3.14) Using Equation 3.14 V can be re written as in Equation 3.15 V V ˆ g 1V C b C a u ˆ b 1V S b v ˆ b 2V C b S a w ˆ b 3 (3.15) where V 2V 2u 2v 2w 2 From Figure 3.10 relations between the v elocity components and the angles of attack and sideslip can be deri v ed. tan aw u (3.16) sin b v V (3.17) Though the matrix G bod y B in Equation 3.14 has been dened for the body-x ed reference frame and v elocity of CG of the torpedo, the equation is v alid for an y other rigid part of the body lik e the ns and ca vitator Thus, in case of a n, the v elocity V w ould correspond to a point (lik e the tip, center of pressure etc.) on that n, and G f in B matrix w ould correspond to the n-x ed reference frame. In this case the v elocity of center of pressure of the n will be used to dene the abo v e parameters. Thus, obtaining a f in and b f in is a tw o step process:

PAGE 31

20 1. Obtain the v elocity of center of pressure of n. V C P bod yV cgE w Br cg C P (3.18) where V C P bod y is v elocity of CP of n in B frame, V cg is the v elocity of CG of the torpedo in E frame, E w B is angular v elocity of B in E and r cg C P is position v ector from CG to CP f in Equation 3.18 can be re written as in Equation 3.19 u f in v f in w f in nB u v w ncg ˆ b 1 ˆ b 2 ˆ b 3 p q r X f in Y f in Z f in (3.19) where r cg C PX f in ˆ b 1Y f in ˆ b 2Z f in ˆ b 3 2. T ransform the v elocity (in E ) of CP of n from frame B to frame of corresponding n. This transformation is obtained by using rotation matrices deri v ed in Equations 3.3 to 3.11 u R 1 v R 1 w R 1 nR 1 r C d R 1 0S d R 1 0 1 0 S d R 1 0 C d R 1 r C q R 1 0 S q R 1S q R 1 0C q R 1 01 0 u R 1 v R 1 w R 1 nB (3.20) u R 2 v R 2 w R 2 nR 2 r C d R 2 0S d R 2 0 1 0 S d R 2 0 C d R 2 r C q R 2 0S q R 2S q R 2 0 C q R 2 01 0 u R 2 v R 2 w R 2 nB (3.21) u E 1 v E 1 w E 1 nE 1 r C d E 1 0S d E 1 0 1 0 S d E 1 0 C d E 1 r C q E 1S q E 1 0S q E 1 C q E 1 0 0 01 u E 1 v E 1 w E 1 nB (3.22) u E 2 v E 2 w E 2 nE 2 r C d E 2 0S d E 2 0 1 0 S d E 2 0 C d E 2 r C q E 2 S q E 2 0S q E 2C q E 2 0 0 0 1 u E 2 v E 2 w E 2 nB (3.23) The left hand terms in Equations 3.20 to 3.23 gi v e the v elocity components at the CP for corresponding ns, in that n frame. These can be used in Equations 3.16 and 3.17 to nd the angle of attack and sideslip for a particular n. 3.1.5 Kinematic Equations V elocity of the CG of the torpedo has been dened in the pre vious section. That v elocity denes the translational kinematics for the torpedo. Analogous to this quantity angular v elocity is required to dene the rotational kinematics. The angular v elocity of the hull has components p q and r in the frame B E w B Dp ˆ b 1q ˆ b 2r ˆ b 3 (3.24)

PAGE 32

21 As the transformation from E to B has already been dened in terms of rotational motions gi v e by Euler angles, the angular rates can also be obtained by dif ferentiation of Euler angles. Thus, another form of Equation 3.24 can be written as in Equation 3.25 E w B Y ˆ e 3 Q ˆ x 2 F ˆ b 1 (3.25) The rotation matrices from Equation 3.1 can be substituted into Equation 3.25 to obtain Equation 3.26 E w B FS Q Yˆ b 1 Y C Q S F Q C Fˆ b 2 Y C Q C F Q S Fˆ b 3 (3.26) Equations 3.24 and 3.26 can be equated to obtain Equation 3.27 p q r n r S Q 0 1 C Q S F C F 0 C Q C FS F 0 Y Q F n(3.27) Equation 3.27 can be re written as in Equation 3.28 p FS Q Y q Y C Q S F Q C F r Y C Q C F Q S F (3.28) T o apply Ne wton' s La ws, acceleration of the CG is required. The v alues of p q r will be the angular accelerations of torpedo in B The translational acceleration can be calculated by time dif ferentiation of V in Ne wtonian frame. A dif ferentiation formula can be used to nd the time deri v ati v e, in some frame, for a v ector dened in some other related frame. d d tv Id d tv BI w Bv (3.29) where, subscript I denotes Ne wtonian (inertial) frame, and B is the body-x ed frame. I w B is angular v elocity of the body (or body-x ed frame) in the I frame, v is the v elocity in I frame, of the point where acceleration is desired. Using the formula the acceleration of

PAGE 33

22 CM of torpedo in E can be obtained. E a C M u v w n ˆ b 1 ˆ b 2 ˆ b 3 p q r u v w (3.30) uqwvr vurpw wpvuq nˆ b (3.31) Similarly the rotational acceleration will be required in the frame E This can be written analogous to Equation 3.30 E a B p q r n ˆ b 1 ˆ b 2 ˆ b 3 p q r p q r p q r n(3.32) The position of torpedo can be found by inte grating the v elocity Letxyzrepresent the coordinates of CG in frame E Thus, the time deri v ati v e of these coordinates in E should represent the v elocity components of the torpedo in E frame. When these time deri v ati v es are transformed to body-x ed frame, the y w ould be equi v alent to the v elocity components in body-x ed frame. x y z nB u v w n(3.33) Equation 3.1 can be substituted in Equation 3.33 to obtain Equation 3.34

PAGE 34

23 x y z nE r C Q C Y C Q S YS Q C Y S F S QC F S Y S F S Q S YC Y C F S F C Q C F S Q C YS F S Y C F S Q S YS F C Y C F C Q u v w n(3.34) 3.2 Dynamic Equations of Motion No w that the accelerations of v arious parts of the torpedo are kno wn, Ne wton' s La ws can be used to deri v e the dynamic equations of motion. Ne wton' s la ws state that the rate of change of momentum is equal to the sum of e xternal force applied on the body Equation 3.35 can be obtained from Ne wton' s la ws by an assumption that the mass of the torpedo is constant. This assumption is v alid for a short period of time. The equations are F Pm a (3.35) where P is the linear momentum, m is mass of the body a is the acceleration and F is all the forces of the body Using Equation 3.31 in Equation 3.35 Ne wton' s La ws for the torpedo can be re written as in Equation 3.36 m uqwvr vurpw wpvuq nˆ bF (3.36) Ne wton' s la ws can be e xtended to rotation. Equation 3.37 are the Ne wton' s La ws for rotational motion. M HI C M aE w BH (3.37) where H (I C M E w B ) is the angular momentum, I C M 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 abo v e stated Ne wton' s la ws are only v alid when the quantities are calculated in an inertial frame of reference. Thus, the quantities ha v e been calculated in frame E Using Equation 3.32 the Ne wton' s La w for rotation can be written as in Equation 3.38

PAGE 35

24 I 1 0 0 0 I 2 0 0 0 I 3 p q r n ˆ b 1 ˆ b 2 ˆ b 3 p q r I 1 p I 2 q I 3 r M (3.38) T o deri v e the equations, the forces on each of the parts will be calculated indi vidually and then e xpressed in body-x ed reference frame, i.e., summation will be done in body reference frame. The rotation matrices deri v ed in pre vious sections will be used e xtensi v ely for this purpose. V arious types of forces are e xperienced by a mo ving torpedo in w ater These forces can be summarized as follo ws:Hydr odynamic F or ces : These are the forces e x erted by the uid on the torpedo. Thus the y e xist whene v er the uid is in contact with body F or superca vitating v ehicle, most of the body (hull) is inside the ca vity Only a portion of the ns and the ca vitator are in contact with the uid. Thus the lift and drag on ca vitator and ns are only hydrodynamic forces. The coef cients of lift and drag for the ns and ca vitator are functions of their angle of attack, immersion, sweepback angle, angle of rotation etc. and are tab ulated in a database [ 5 ]. This database will be interpolated and e xtrapolated to calculate the numerical v alues for the forces. Occasionally a part of the hull comes in contact with w ater and might e xperience some forces. These forces are kno wn as planing forces. It is assumed that the v ehicle is centered in the ca vity Thus the planing forces are ne glected in the summation of the hydrodynamics forces. F H yd r od ynamicF R 1F R 2F E 1F E 2F c (3.39) M H yd r od ynamicM R 1M R 2M E 1M E 2M c (3.40)Gra vitational F or ces : This is the gra vity forces on the body As the summation of moments will be tak en about the center of gra vity the gra vitational forces will not contrib ute to the summation on moments. The y will appear only in summation of forces.Pr opulsi v e : The torpedo has a propulsion system, which is similar to that of rock ets. The line of action of the propulsi v e force is assumed to be passing through center of gra vity and along ˆ b 1 axis. Thus this force will also contrib ute to the sum of forces, and not moments. The propulsi v e forces are pro vided by b urning the fuel, b ut for simplicity it will be assumed that the mass of the torpedo remains constant. The total forces and moments are e xpressed in terms of these components. FF H yd r od ynamicF Gr avF Pr o p (3.41) MM H yd r od ynamic (3.42)

PAGE 36

25 3.2.1 F or ces on Ca vitator Figure 3.11 sho ws the forces acting on the ca vitator Coef cient of lift ( cl c ) and drag ( cd c ) for the ca vitator are functions of angle of attack, a c and half-angle, g 1 2 of the ca vitator Half-angle, for a cone, is the angle made by axis of the cone with an y line passing through the v erte x and lying in the surf ace of the the cone. This parameter denes the main geometry of the conical ca vitator The v alues of cl c and cd c are determined using CFD and tab ulated in [ 5 ]. These v alues ha v e been e xtrapolated to calculate lift and drag. L c1 2 r V 2 c S c cl ca cg 1 2(3.43) D c1 2 r V 2 c S c cd ca cg 1 2(3.44) where S c is the cross-sectional area of the ca vitator base. Directions of lift ( L c ) and drag ( D c ) are as sho wn in Figure 3.11 (b). These can be transformed to the body ax es using 3.2 and 3.14 for the ca vitator ˆ b 1 ˆ b 2 ˆ b 3 nC Bd c G cav Ca cb c ˆ g 1 ˆ g 2 ˆ g 3 ncav (3.45) b c a c ˆ c 1 ˆ c 3 ˆ g 1 ˆ c 2 (a) ˆ c 1 ˆ g 2 ˆ g 3 M c D c L c ˆ g 1 (b) CP Figure 3.11 Ca vitator: (a) Angle of Attack and Sideslip and (b) Hydrodynamic F orces

PAGE 37

26 Thus the forces on ca vitator in body frame, can be written. F c F cx F cy F cz nB D ca cg 1 20L ca cg 1 2 nG cav r C d c 0 S d c 0 1 0S d c 0 C d c r C b c C a c C a c S b cS a cS b c C b c 0 C b c S a c S a c S b c C a c D ca cg 1 20L ca cg 1 2 n(3.46) where F c is a 3x1 matrix with each ro w corresponding to each direction in B basis. The forces are assumed to be acting at CP of the ca vitator Once the forces ha v e been calculated, the moment about an y point can be calculated. M cr C PcavF c (3.47) where r C Pcav is the position v ector from that point to CP of ca vitator It is assumed that the CP lies on ˆ b 1 when ca vitator deection is 0, and its coordinates with respect to body origin O in this case, areX c00. Thus from Figure 3.4 the coordinates of CP can be written for the case when the ca vitator is deected. r C Pcav X cD C P C d c 0D C P S d c nbod y (3.48) The moments on the ca vitator in body-x ed can be obtained by substituting Equations 3.46 and 3.48 in Equation 3.47

PAGE 38

27 M c X cD C P C d cˆ b 1D C P S d c ˆ b 3 r C d c 0 S d c 0 1 0S d c 0 C d c r C b c C a c C a c S b cS a cS b c C b c 0 C b c S a c S a c S b c C a c D ca cg 1 20L ca cg 1 2 n(3.49) 3.2.2 F or ces on Fins Fin forces are also e xtrapolated from the CFD database [ 5 ], which gi v es the v alues of coef cients of lift ( cl f in ) and drag ( cd f in ) for ns. These v alues are functions of angle of attack a f in n sweepback q f in n rotation d f in n immersion I f in and n geometry Figure 3.12 sho ws these parameters graphically and the y can be dened as follo ws:Fin Geometry: The geometry parameters for ns are L and S and wedge half angle ( h a ), as sho wn in Figure 3.12 These parameters are x ed according to the v alues gi v en by the database, so as to calculate the forces accurately .Fin Immersion: As the n is partially wetted by uid, the wetted length can be represented by parameter S 0 along n Y -axis. The immersion I f in can be dened as the ratio of the wetted length of the n to its true length. I f in S 0Sf in (3.50) I f in determines the total hydrodynamic force on the n.Fin Rotation ( d f in ): As dened earlier this is rotation about n ˆ f 2 axis. This deter mines the direction of the hydrodynamic force. Thus n rotation is used as primary control for the torpedo.Fin Sweepback ( q f in ): As dened earlier this is rotation about n ˆ f 3 axis. Sweepback determines the wetted surf ace of the n, thus the hydrodynamic force on the n.Angle of Attack: Angle of attack for the n is calculated as described in Figure 3.12 and section 3.1.4 The database gi v es cl f in and cd f in as a function of a f in q f in and I f in thus lift and drag on the ns can be calculated by the normalizing f actor L f in1 2 r V 2 S 2 f in cl f in (3.51) D f in1 2 r V 2 S 2 f in cd f in (3.52)

PAGE 39

28 wetted region D f in ˆ f f in 3 L f in ˆ f f in 1 ˆ f f in 2 V h a S L S 0 Figure 3.12 Fin Geometry Where S f in is the length of the n as sho wn in Figure 3.12 These forces ha v e directions as sho wn in Figure 3.12 Before substituting in Equation 3.39 these forces ha v e to be transformed to body-x ed reference frame. This process in v olv es follo wing tw o rotations: 1. Rotate the frame F f in (which has L f in and D f in along its basis v ectors) to align it with the n-x ed frame using Equation 3.14 and 2. Rotate the abo v e obtained n-x ed frame to obtain the body-x ed frame using Equations 3.3 to 3.11 Thus the forces on the ns in body-x ed frame axis can be obtained.Rudder 1 F R 1x F R 1y F R 1z nB r C q R 1S q R 1 0 0 01 S q R 1C q R 1 0 r C d R 1 0 S d R 1 0 1 0S d R 1 0 C d R 1 r C b R 1 C a R 1 C a R 1 S b R 1S a R 1S b R 1 C b R 1 0 C b R 1 S a R 1 S a R 1 S b R 1 C a R 1 D R 1 0L R 1 n(3.53)Rudder 2 F R 2x F R 2y F R 2z nB r C q R 2S q R 2 0 0 01S q R 2 C q R 2 0 r C d R 2 0 S d R 2 0 1 0S d R 2 0 C d R 2 r C b R 2 C a R 2 C a R 2 S b R 2S a R 2S b R 2 C b R 2 0 C b R 2 S a R 2 S a R 2 S b R 2 C a R 2 D R 2 0L R 2 n(3.54)

PAGE 40

29Ele v ator 1 F E 1x F E 1y F E 1z nB r C q E 1S q E 1 0S q E 1 C q E 1 0 0 01 r C d E 1 0 S d E 1 0 1 0S d E 1 0 C d E 1 r C b E 1 C a E 1 C a E 1 S b E 1S a E 1S b E 1 C b E 1 0 C b E 1 S a E 1 S a E 1 S b E 1 C a E 1 D E 1 0L E 1 n(3.55)Ele v ator 2 F E 2x F E 2y F E 2z nB r C q E 2S q E 2 0 S q E 2C q E 2 0 0 0 1 r C d E 2 0 S d E 2 0 1 0S d E 2 0 C d E 2 r C b E 2 C a E 2 C a E 2 S b E 2S a E 2S b E 2 C b E 2 0 C b E 2 S a E 2 S a E 2 S b E 2 C a E 2 D E 2 0L E 2 n(3.56) Equations 3.53 to 3.56 gi v e the components of hydrodynamic forces on ns in body-x ed frame. What remains is to nd the moment of these forces about CG of body The moments can be obtained in analogous to Equation 3.47 M f inr f in C GC PF f in (3.57) In this case, the root of ns is x ed to body and it can mo v e thus mo ving the CP of n relati v e to root. The position of CP from root is also kno w with respect to n coordinates. r f in C Gr oo tX f in r oo t ˆ b 1Y f in r oo t ˆ b 2Z f in r oo t ˆ b 3 (3.58) r f in r oo tC PD x f in C P ˆ f 1D y f in C P ˆ f 2 (3.59) whereˆ f 1ˆ f 2ˆ f 3is n-x ed coordinates for corresponding n. Equations 3.58 and 3.59 can be added by using matrices gi v en in Equations 3.3 to 3.11 Thus, the position v ector from CG to CP of ns can be obtained.

PAGE 41

30Rudder 1 X R 1 Y R 1 Z R 1 nB X r oo t R 1 Y r oo t R 1 Z r oo t R 1 nB r C q R 1S q R 1 0 0 01 S q R 1C q R 1 0 r C d R 1 0 S d R 1 0 1 0S d R 1 0 C d R 1 D x R 1 C P D y R 1 C P 0 nR 1 (3.60)Rudder 2 X R 2 Y R 2 Z R 2 nB X r oo t R 2 Y r oo t R 2 Z r oo t R 2 nB r C q R 2S q R 2 0 0 01S q R 2 C q R 2 0 r C d R 2 0 S d R 2 0 1 0S d R 2 0 C d R 2 D x R 2 C P D y R 2 C P 0 nR 2 (3.61)Ele v ator 1 X E 1 Y E 1 Z E 1 nB X r oo t E 1 Y r oo t E 1 Z r oo t E 1 nB r C q E 1S q E 1 0S q E 1 C q E 1 0 0 01 r C d E 1 0 S d E 1 0 1 0S d E 1 0 C d E 1 D x E 1 C P D y E 1 C P 0 nE 1 (3.62)Ele v ator 2 X E 2 Y E 2 Z E 2 nB X r oo t E 2 Y r oo t E 2 Z r oo t E 2 nB r C q E 2S q E 2 0 S q E 2C q E 2 0 0 0 1 r C d E 2 0 S d E 2 0 1 0S d E 2 0 C d E 2 D x E 2 C P D y E 2 C P 0 nE 2 (3.63) Equations 3.60 to 3.63 gi v e the position v ector from CG to CP of the ns. These equations in conjunction with Equations 3.53 to 3.56 used in 3.57 gi v es the e xternal moments on ns about the CG M f in ˆ b 1 ˆ b 2 ˆ b 3 X f in Y f in Z f in F f inx F f iny F f inz (3.64)

PAGE 42

31 3.2.3 Gra vitational F or ces F or simplicity mass ( m ) of the torpedo is assumed to be constant o v er time. This is not the case in reality b ut the approximation is reasonable for considering short time maneuv ers. Acceleration due to gra vity g is also assumed to be constant as torpedo mo v es in space. The direction of gra vity is gi v en by ˆ e 3 axis. Thus, the gra vitational force can be written as in Equation 3.65 F gr avmg ˆ e 3 (3.65) Equation 3.65 can be re-e xpressed in body frame of reference using Equation 3.1 F gr av r C Q C Y C Q S YS Q C Y S F S QC F S Y S F S Q S YC Q C F S F C Q C F S Q C YS F S Y C F S Q S YS F C Y C F C Q 0 0 mg nEmg S Q S F C Q C F C Q nB (3.66) 3.2.4 Equations of Motion No w that a mathematical formulation of forces on the torpedo has been achie v ed, these equations can be substituted into Equations 3.36 to 3.42 to obtain the dynamic equations of motion. Thus the force equations of motion can be summarized as in Equation 3.67

PAGE 43

32 3.2.4.1 F or ce Equations m uqwvr vurpw wpvuq nBF immer sion F pr o p 0 0 nB F R 1x F R 1y F R 1z nB F R 2x F R 2y F R 2z nB F E 1x F E 1y F E 1z nB F E 2x F E 2y F E 2z nBmg S Q S F C Q C F C Q nB r C d c 0 S d c 0 1 0S d c 0 C d c r C b c C a c C a c S b cS a cS b c C b c 0 C b c S a c S a c S b c C a c D ca cg 1 20L ca cg 1 2 nC (3.67) Some issues should be noted for Equation 3.67 .The planing forces ha v e not been included in the equations of motion. These forces are ne glected by assumption that the v ehicle is centered in the ca vity .The propulsion force is constrained to be along ne gati v e ˆ b 1 axis. 3.2.4.2 Moment Equationsr I 1 0 0 0 I 2 0 0 0 I 3 p q r n ˆ b 1 ˆ b 2 ˆ b 3 p q r I 1 p I 2 q I 3 r ˆ b 1 ˆ b 2 ˆ b 3 X R 1 Y R 1 Z R 1 F R 1x F R 1y F R 1z ˆ b 1 ˆ b 2 ˆ b 3 X R 2 Y R 2 Z R 2 F R 2x F R 2y F R 2z ˆ b 1 ˆ b 2 ˆ b 3 X E 1 Y E 1 Z E 1 F E 1x F E 1y F E 1z ˆ b 1 ˆ b 2 ˆ b 3 X E 2 Y E 2 Z E 2 F E 2x F E 2y F E 2z ˆ b 1 ˆ b 2 ˆ b 3 X cD C P C d c 0D C P S d c F cx F cy F cz (3.68) Some issues should be noted for Equation 3.68

PAGE 44

33Some of the terms in Equation 3.68 are sho wn as a determinant. The y need to be e xpanded and written as components in body-x ed frame so as to equate the lefthand and right-hand terms.Moments due to gra vitation do not appear because the moments are tak en about CG .Again, the moments due to planing forces ha v e been ne glected. 3.2.4.3 Orientation Equations Y Q F n r 0 S F C Q C F C Q 0 C FS F 1 S F S Q C Q C F S Q C Q p q r n(3.69) 3.2.4.4 P osition Equations x y z nE r C Q C Y C Q S YS Q C Y S F S QC F S Y S F S Q S YC Y C F S F C Q C F S Q C YS F S Y C F S Q S YS F C Y C F C Q u v w n(3.70) Equations 3.67 to 3.70 are a complete set of equations of motions for the superca vitating torpedo.

PAGE 45

CHAPTER 4 LINEARIZA TION 4.1 Linearization 4.1.1 Need f or Linearization The equations of motion for the torpedo are identical to airplane equations of motion b ut the forces terms on the right-hand side of these equations are dif ferent. Thus, the linearization can be carried out similarly as sho wn in Nelson [ 9 ]. The equations of motion, as in the case of a superca vitating torpedo, are represented by a set of rst-order dif ferential equations. xfxu(4.1) using f : n n as a nonlinear function of a time-v arying v ector x n and u n F or control design, the system dynamics are observ ed at some trim conditions by gi ving perturbations to states of the system at that trim. The dynamics associated with these perturbations are obtained by linearization. An adv antage 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. Y et, there are fe w disadv antages for this processLinearized equations can predict the system performance only in a small range of operations. The linearized equations change as the operating point of system changes, thus making it dif cult for simulating true beha vior of system.Information relating to nonlinearities lik e hysteresis, backlash, coulomb friction, discontinuities etc. may be lost by linearizing the system.A controller that is good for linearized system might ha v e v ery poor performance for the nonlinear equations. An iterati v e process may be needed to nd a controller that is good for nonlinear equations. 34

PAGE 46

35 4.1.2 Generic F orm of Equations of Motion The equations of motion in Equations 3.67 and 3.70 can be written in simplied form using sums of total forces and moments acting on the body m uqwvrgS Q X m vr upwg C Q S F Y m wpvqug C Q C F Z (4.2) I x pqrI zI y L (4.3) I y qr pI xI z M (4.4) I z rpqI yI x N (4.5) Y Q F n r 0 S F C Q C F C Q 0 C FS F 1 S F S Q C Q C F S Q C Q p q r n(4.6) x y z nE r C Q C Y C Q S YS Q C Y S F S QC F S Y S F S Q S YC Y C F S F C Q C F S Q C YS F S Y C F S Q S YS F C Y C F C Q u v w n(4.7) These equations of motions are coupled by the state v ector x and are dependent on the control v ector u x uvwpqrYQFxyzu q R 1q R 2q E 1q E 2d R 1d R 2d E 1d E 2d cF pr o p(4.8) 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 ight condition), i.e., the equations thus deri v ed will be v alid for the torpedo operating at and near that v alue of v ector x The operating point is chosen to correspond to the trim condition in Equation 4.9

PAGE 47

36 x 0 u 0v 0w 0p 0q 0r 0Y 0Q 0F 0x 0y 0z 0 u 0000000Q 00000(4.9) This corresponds to straight and le v el ight with constant v elocity As the torpedo may be tra v eling at other ight conditions, the linearization at those conditions w ould be carried out numerically which will be e xplained in later sections. A v alue of u 0 is found numerically that satises the equations of motion for a gi v en v alue of x 0 Then a disturbance of D x is gi v en to the equations of motion from the reference condition thus changing the ight conditions to x 0D x Se v eral assumptions are made to carry out the linearization:The disturbances from reference ight condition are small. Thus the terms in v olving higher order of disturbances D will be ne glected. Furthermore rst order terms in v olving D will be approximated as in Equation 4.10 S inD DC osD 1 (4.10)The propulsi v e forces and mass are assumed to be constant.Planing and immersion forces are ne glected for this analysis. The linearization procedure is presolv ed for the force equation in ˆ b 1 direction. This equation relates the force, X to the states. m uqwr ugS Q X (4.11) Using the ight condition from Equation 4.9 in Equation 4.11 gi v es the v alue of force at the reference trim condition. mgS Q 0X 0 (4.12) The perturbation equation, i.e., the equation with ight condition disturbed by D x can be obtained by substituting the perturbed ight condition into Equation 4.11 md d tu 0D u q 0D q w 0D w r 0D r u 0D u gSQ 0DQ X 0D X (4.13) Equations 4.12 and 4.13 can be combined to gi v e the linearized form of Equation 4.11 for straight and le v el ight condition.

PAGE 48

37 mD ug DQ C Q 0 D X (4.14) Proceeding in a similar w ay all other equations of motion can be linearized. The linearized equations for straight le v el ight are sho wn in Equation 4.15 to Equation 4.18 4.1.3.1 F or ce Equations mD ug DQ C Q 0 D X mD vu 0 D rg DF C Q 0 D Y mD uu 0 D qg DQ S Q 0 D Z (4.15) 4.1.3.2 Moment Equations I x D pD L I y D qD M I z D rD N (4.16) 4.1.3.3 Orientation Equations D YD r C Q 0 D QD q D FD pT Q 0 D r (4.17) 4.1.3.4 P osition Equations D x S Q 0 u 0 DQC Q 0 D uS Q 0 D w D yD v D z C Q 0 u 0 DQS Q 0 D uC Q 0 D w (4.18) 4.1.4 Stability and Contr ol Deri v ati v es The v ariations in total force and moment are often dif cult to compute.These v ariations in forces can be calculated by chain rule for deri v ati v es. As stated in Equation 4.8 these are functions of state v ariables x and control v ariables u Thus for e xample D X can be written by chain rule as in Equation 4.19

PAGE 49

38 D XX u D uX v D vX w D wX p D pX q D qX r D rX Y DYX Q DQX F DFX pr o p D F Pr o pX q R 1 q R 1X q R 2 q R 2X q E 1 q E 1X q E 2 q E 2X d R 1 d R 1X d R 2 d R 2X d E 1 d E 1X d E 2 d E 2X d c d c (4.19) where the subscripted X represents its partial deri v ati v e with respect to its subscript. X u X u xx 0 (4.20) Each of these subscripted v ariables that ha v e a subscript of state v ariable are kno wn as stability deri v ati v es and the ones with a control v ariable as subscript are kno wn as a control deri v ati v e. There can be as man y stability deri v ati v es as there are forces and state and control v ariables. Man y of these are ne gligible, depending on the reference ight condition. These dependencies are kno wn usually by e xperimentation or numerical calculations. F or e xample, the force, X is observ ed to depend mainly on a subset of the state and control v ariable. Thus only the stability deri v ati v es that correspond to these v ariables ha v e to be retained in Equation 4.19 when straight and le v el ight is considered. Xf unc tuwd E 1d E 2q E 1q E 2d cF pr o p(4.21) The ne xt problem is calculating numerical v alues of these deri v ati v es. In most cases it is easy to calculate these numerically or using a symbolic manipulation softw are. F or some reference points, it is possible to do the calculation manually The calculation of X u will be done manually for straight and le v el ight. X u uF R 1xF R 2xF E 1xF E 2xF cxF pr o px(4.22) Expressions for each of the terms in Equation 4.22 ha v e been deri v ed in Chapter 3 F or e xample, the e xpression for the force on ca vitator is sho wn in Equation 4.23 F cx C d c 0 S d c r C b c C a c C a c S b cS a cS b c C b c 0 C b c S a c S a c S b c C a c D ca cg 1 20L ca cg 1 2 nC (4.23)

PAGE 50

39 In Equation 4.23 a c b c and thus L c and D c are the only functions of u Thus the partial deri v ati v es with respect to u can be obtained. u F cx C d c 0 S d c r S b c C a c b c uC b c S a c a c uS a c S b c a c uC a c C b c b c uC a c a c uC b c b c uS b c b c u 0S b c S a c b c uC b c C a c a c u S a c C b c b c uC a c S b c a c uS a c a c u D ca cg 1 20L ca cg 1 2 nC r C d c 0 S d c 0 1 0S d c 0 C d c r C b c C a c C a c S b cS a cS b c C b c 0 C b c S a c S a c S b c C a c u D ca cg 1 20 u L ca cg 1 2 nC (4.24) It can be seen that a c u b c u L c u and D c u are terms required to be calculated. These can be calculated from equations 3.16 and 3.17 which are restated in Equations 4.25 to Equation 4.27 tana c w c u c (4.25) tanb c v c V c (4.26) V 2 cu 2 cv 2 cw 2 c (4.27) The v elocity components in Equation 4.27 can be found using Equation 3.2

PAGE 51

40 u c v c w c nC r C d c 0S d c 0 1 0 S d c 0 C d c u c v c w c nB (4.28) u c v c w c nB u v w nB ˆ b 1 ˆ b 2 ˆ b 3 p q r x c y c z c (4.29) No w the v elocity components can be obtained for the reference ight condition that is stated in Equation 4.9 u c v c w c nC C d c u 0 0 S d c u 0 n(4.30) The v ariation of a c can be obtained by dif ferentiating Equation 4.25 at reference ight condition. sec 2a cd a cu c d w cw c d u c u 2 cd a cu c d w cw c d u c u 2 cw 2 c (4.31) d a cC d c d w c u 0S d c d u c u 0 atx 0u 0(4.32) Similarly the v ariation of b c can be obtained by dif ferentiating Equation 4.26 at reference ight condition. d b c V c d v cv c dV c V c V 2 cv 2 c d v c u 0 atx 0u 0(4.33) No w using Equations 4.28 and 4.29 v ariation of v elocity components of ca vitator with respect to u can be obtained. u c uC d c v c u0 w c uS d c (4.34)

PAGE 52

41 Finally combining Equations 4.32 to 4.34 the v ariations of a c and b c with respect to u can be obtained at reference ight condition. a c uC d c u 0 w c uS d c u 0 u c uC d c S d c u 0S d c C d c u 00 (4.35) b c u 1 u 0 v c u0 (4.36) Clearly tw o of the deri v ati v es that are required to calculate stability deri v ati v es ha v e been obtained. It w as pre viously stated that lift and drag are calculated using the coef cient of lift and drag. L c1 2 r V 2 c S c cl c D c1 2 r V 2 c S c cd c (4.37) These forces can be dif ferentiated by chain rule the deri v ati v e w ould be lik e in Equation 4.38 L c u1 2 r S c2 V c cl c V c uV 2 c cl c a c a c u(4.38) The Equation 4.38 requires tw o deri v ati v es. One of the deri v ati v es is calculated in Equation 4.35 The other deri v ati v e can be calculated using Equation 4.27 V c u u u 2 cv 2 cw 2 c u c u c uv c v c uw c w c u u 2 cv 2 cw 2 c1 (4.39)

PAGE 53

42 Thus the deri v ati v e of the lift and drag forces can be obtained. L c ur S c V c cl c (4.40) D c ur S c V c cd c (4.41) Thus all the deri v ati v es required to calculate right-hand side of Equation 4.24 ha v e been calculated. All the terms on right-hand side of the Equation 4.20 can be calculated in a similar manner It is tedious to calculate the deri v ati v es analytically in such a w ay The comple xity increases for other ight conditions. F or practical purposes these deri v ati v es are calculated using symbolic manipulation softw ares lik e `Mathematica' or by using numerical methods. The numerical methods used to calculate the deri v ati v es ha v e been described in Appendix B 4.2 State Space Repr esentation Equations 4.15 to 4.18 are a complete set of linearized equations of motion for the torpedo. The y can be represented in a more con v enient form kno wn as the State Space F orm. The state space equations are a set of rst-order dif ferential equations. xAxBu yC xDu x nu py m A n nB n p C m nD m p (4.42) Equation 4.42 is a generalized form of state space representation for an y system. Each of the terms in the equations has a particular importance for describing the dynamics of the system.State V ariable x : The state v ariables for a system are a set of v ariables, when kno wn at time t 0 and along with input u are suf cient to determine the state of the system at an y time tt 0 All the states of the system need not be measurable.

PAGE 54

43Input V ariable u : This is the control surf ace deections.Output V ariable y : The output v ariables are the measured parameters. These may or may not be same as the state v ariables. The output v ariables are usually considered to be measurable b ut sometimes the y are estimated. The matrices ABC and D may either be constant or time-v arying functions. In the case of the superca vitating torpedo, the state v ector is of size 12 (n) and the control v ector is of size 10 (p). x D u D w D q DQ D v D p D r DF DY D x D y D zu d c d E 1 d E 2 d R 1 d R 2 q E 1 q E 2 q R 1 q R 2 D F pr o p(4.43) Some of these controls may not be needed for some maneuv ers. From the linearized equations it can be observ ed that the state v ariables are not coupled by the statesYxyz. These four states can be remo v ed from the analysis for no w The system becomes a 8 state system. These states can be further di vided into longitudinal and lateral-directional dynamics. The v ariables D uD wD qDQ correspond to longitudinal dynamics, which also means the dynamics in ˆ b 1 ˆ b 3 plane. The v ariables D vD pD rDF correspond to lateral dynamics, which is the dynamics in ˆ b 1 ˆ b 2 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 v ariables. The plant matrix A can be di vided into four parts. A r A l ong A cou p 1 A cou p 2 A l a t d (4.44) When A is di vided as in equation 4.44 where each element is a 44 matrix, A l ong w ould correspond to longitudinal dynamics and A l a t d w ould correspond to lateral dynamics. A cou p 1 and A cou p 2 are coupling matrices. It is required that the coupling matrices become ne gligible for the equations to be decoupled. If these parts are not ne gligible, the equations cannot be decoupled, and a 8 state model will be required to be considered. From linearized

PAGE 55

44 equations the four parts of the A matrix for the torpedo can be written as in Equation 4.45 to Equation 4.48 A l ong r X u mq 0X w mw 0X q mg C Q 0X Q m q 0Z u m Z w m u 0Z q mg C F 0 S Q 0Z Q m M u I y M w I y M q I y M Q I y 0 0 C F 0 0 (4.45) A l a t d r Y v m Y p mu 0Y r mg C Q 0 C F 0Y F m L v I x L p I x L r I zI yq 0 I x L F I x N v I z N p I yI xq 0 I z N r I z N F I z 0 1 C F 0 T Q 0 q 0 C F 0 S Q 0r 0 S F 0 S q 0 C Q 0 (4.46) A cou p 1 r r 0X v m X p m v 0X r m X F m p 0Z v mv 0Z p m Z r mgS F 0 C Q 0Z F m M v I y M p I xI zr 0 I y M r I xI zp 0 I y M F I y 0 0S F 0S F 0 q 0C F 0 r 0 (4.47) A cou p 2 r r 0Y u m p 0Y w m Y q mgS Q 0 S F 0Y Q m L u I x L w I x L q I zI yr 0 I x L Q I x N u I z N w I z N q I yI xp 0 I z N Q I z 0 0 S F 0 T Q 0 S F 0 q 0C F 0 r 0 (4.48) Similarly B is a 810 matrix, whose elements are just the control deri v ati v es according to their locations in the matrix. The rst 4 ro ws correspond to longitudinal dynamics and the last 4 correspond to lateral dynamics. B l ong r X d c mX F pr o px m . . . . 00 410 (4.49)

PAGE 56

45 B l a t d r Y d c mY F pr o px m . . . . 00 410 (4.50) No w the complete state space representation for the torpedo can be written as in Equation 4.51 which gi v es tw o sets of equations. The rst set is the longitudinal equations and the second set is the lateral-directional equations. x l ong D u D w D q DQT x l a t d D v D p D r DFT u d c d E 1 d E 2 d R 1 d R 2 q E 1 q E 2 q R 1 q R 2 D F pr o p x l ong x l a t d n81 r A l ong A cou p 1 A cou p 2 A l a t d 88 x l ong x l a t d n81 r B l ong B l a t d 810 u 101 (4.51)

PAGE 57

CHAPTER 5 CONTR OL DESIGN SETUP This chapter deals with the control design for the torpedo described in pre vious chapters. V arious parameters associated with the control are restated in T able 5.1 T able 5.1 Control P arameters Longitudinal Lateral State u w q Q v p r F Y Control d c d e 1 d e 2 d r 1 d r 2 It should be noted that Y has been included in the states though it w as observ ed in the state matrices that all other v ariables are independent of Y It will be seen later that the inclusion of Y in the feedback states helps in impro v ement of performance. Also, for longitudinal control, tw o ele v ators and the ca vitator are required. Similarly for lateral direction, the rudders should be enough for control. There are v arious control methods, lik e linear quadratic re gulator (LQR) synthesis, synthesis etc., which can be used to design a controller Each of these methods has adv antages and disadv antages. LQR method gi v es a constant gain controller which is based on minimization of a quadratic performance inde x and considers the problem of rob ustness only in terms of gain and phase mar gins. -synthesis deals with rob ustness with respect to a wide v ariety of uncertainties to minimize an innity-norm matrix b ut the resulting controller can be of high order Re gardless of comple xity and rob ustness, each design method presents dif culties. The LQR method w as chosen for controller synthesis for the torpedo. This method w as chosen because 1. It is easy to implement in a nonlinear simulation. 2. The resulting rob ustness for the LQR controller w as acceptable. 46

PAGE 58

47 3. It is straight-forw ard to v ary some design parameters to achie v e performance goals. This chapter e xplains v arious problems associated with the control synthesis and the system model used for synthesis of the controller 5.1 Open-loop P erf ormance Initially the equations of motion for the torpedo ha v e been linearized for straight and le v el ight at a forw ard v elocity of 75 ms1 x 7500000000000(5.1) It is found that the ca vitator and tw o ele v ators are suf cient to maintain the abo v e trim. It is also assumed that the v alue of propulsi v e force required to maintain this trim is x ed at the required v alue. u d cd e 1d e 2d r 1d r 2F pr o p 0006700106 001060040 023 e03(5.2) where the deections are gi v en in radians, and F pr o p is in Ne wtons. As these parameters are obtained numerically it may not be a perfect trim. The system may ha v e some nonzero accelerations, and consequently may tend to de viate from straight and le v el ight. T o check this, the open-loop dynamics are simulated at this condition without an y feedback. Figure 5.1 sho ws the Simulink model used for open-loop simulation. The control surf ace deections are x ed at their trim v alues for this simulation. The closed-loop system is obtained using the equations of motion that were deri v ed in Chapter 3 The state deri v ati v es are then inte grated to obtain the state at the ne xt instant. Figure 5.2 sho ws the open-loop response for the torpedo at this trim condition. It can be seen that the open-loop system is unstable. The simulation is carried out at the trim that is sho wn in Equation 5.1 i.e., all the v alues sho wn in these gures ha v e to be zero. The system is seen to ha v e oscillation about non-zero states. The state and control matrices obtained for this trim condition are sho wn in Equations 5.3 to 5.6 There are 5 control v ariables,d cd e 1d e 2d r 1d r 2. The matrices corresponding

PAGE 59

48 State Feedback In1 state For Plotting simt Time MATLAB Function NL Equation Torpedo 1 s Integrator Out Control at Trim Clock Figure 5.1 Simulink Model for Open Loop Simulation 0 20 40 60 80 100 0 5 10 15 20 time(s)w (ms-1) 0 20 40 60 80 100 -0.04 -0.03 -0.02 -0.01 0 0.01 time(s)p (rad s-1) 0 20 40 60 80 100 0 0.5 1 1.5 2 2.5 3 3.5 4 time(s)q (rad s-1) Figure 5.2 Open-Loop Response for T orpedo: wpq

PAGE 60

49 to the lateral dynamics are of dimension 5 because the state Y is included in the lateral dynamics. Thus the lateral states are no wvprFY. A l ong r 45204 15417 131109810002616157648 785888 0 00000 1207735614 0 0 0 10000 0 (5.3) B l ong r 323010 690608690608 0 040609423033736 3033736 0 0 1584675451531 451531 0 0 0 0 0 0 0 (5.4) A l a t d r 12042200002716011 98100 001813542281 03004 0 0 114370002512528 0 0 0 10000 0 0 0 0 0 10000 0 0 (5.5) B l a t d r 0 0 036660511 36660511 014297086142970861727699417276994 01412952314129523 5456162954561629 0 0 0 0 0 0 0 0 0 0 (5.6) The longitudinal eigen v alues are 211414 4513718262 00178and the lateral eigen v alues are0 54228900002 6647272683 6647272683. The eigen v alues clearly sho w that the system is unstable. It can also be seen that the longitudinal dynamics ha v e no oscillatory modes. Figure 5.3 sho ws the v ariation of eigen v alues

PAGE 61

50 -30 -20 -10 0 10 -1 -0.5 0 0.5 1 Longitudinal Egienvalues for u=60:5:90 Real ( l long )Imag ( llong) (a) Longitudinal -80 -60 -40 -20 0 20 -10 -5 0 5 10 Lateral Egienvalues for u=60:5:90 Real ( l latd )Imag ( llatd) (b) Lateral Figure 5.3 V ariation of Eigen v alues with Change in V elocity for the torpedo for dif ferent v elocities. State v alues are x ed e xcept for forw ard v elocity which is v aried from 60 ms1 to 90 ms1 The gures sho w that the v ariation is small and most of the eigen v alues stay in ne gati v e half of comple x plane. 5.2 Closed-Loop Pr oblem As stated earlier the control problem can be subdi vided into v arious problems. Each can be solv ed to get a nal controller The ultimate goal of the controller design is to achie v e a desired trajectory or reach a particular point with minimization of some perfor mance criteria. The achie v ement of this goal requires addressing maneuv ering, trimming, guidance and na vigation. This thesis will consider the basic problem of maneuv ering. 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:T rack a step command in pitch or roll rate of size up to 30 d e gs .Maintain an o v ershoot less than 15%.Ha v e a rise time of less than 0.6 sec.Ha v e a steady state error of less than 5%. Besides meeting the abo v e mentioned performance criteria, the controller is also required to stabilize the closed-loop system.

PAGE 62

51 T able 5.2 Control Constraints Ca vitator Deection 15d c 15 Ca vitator Rate 25 s d c 25 s ns 60d f 60 Fin Rate 25 s d f 25 s Thrust 0F pr o p30000 N V arious bounds are placed on the control surf ace deections and rates. These bounds are listed in T able 5.2 The bounds are included in the nonlinear simulations and it is required that there is no saturation of deection or the rates at least for the commands ha ving the rate 30 d e gs An actuator model is included in the system so as to constrain the rates of the control surf ace motion. This actuator is realized as a lo w pass lter A c80 s80 Addition of this lter synthesizes a controller that has slo wer control deections. Let q commtbe a function of time, dening the desired pitch rate prole. The aim of the controller is to nd a control la w utthat yields an achie v ed pitch rate, q ac hit, to minimize the optimization problem stated in Equation 5.7 nd utthat minimizes zt q ac hit q commt subject to u minuu max u min u u max xAxBu (5.7) where, u min and u max are lo wer and upper bounds on control deections. The quantities u min and u max are lo wer and upper bounds on control deection rates. The problem becomes a disturbance rejection problem, when the commanded v ariable is 0 for all time. This is an optimization problem,where the state space equation is an equality constraint and the control surf ace bounds are inequality constraints.

PAGE 63

52 5.3 Rob ustness of the Contr oller A control system that is designed to accommodate the uncertainties in a mathematical model is called a rob ust control system. Usually the response of a model does not accurately match the response of the true system. A rob ust control system should gi v e the desired performance not only during the simulations of the model, b ut for the true system also. V arious 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 mar gins are generally used to predict the rob ustness of a control system. Physically these are just the f actors by which the feedback gain can be increased and yet ha v e a stable real system. A formal denition of these can be gi v en by using a frequenc y analysis for a feedback system. 5.3.1 Gain Mar gin Figure 5.4 sho ws a typical feedback system in v olving a plant, P and a controller K. The gain for the system in dotted re gion is kno wn 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 mar gin is the lar gest f actor by which this gain can be increased and still ha v e a stable system. Physically it means if the response of the torpedo for a gi v en ele v ator input is higher by a f actor of the gain mar gin, the torpedo is still stable. The gain mar gin is usually e xpressed in decibel (db) units, and can be easily obtained from the Bode plots for the system. The gain mar gin is a measurement of the magnitude on the Bode plot, at the point where the phase is 180 o 5.3.2 Phase Mar gin Gain is a v alid rob ustness criteria when the system has real eigen v alues. But usually the eigen v alues ha v e imaginary components and thus the phase is also a concern. Phase mar gin is the measure of the maximum possible phase lag before the system becomes unstable. From the Bode plot, phase mar gin is the phase when the magnitude of the gain is zero.

PAGE 64

53 + K P Figure 5.4 Loop Gain 5.3.3 Uncertainty In P arameters Another f actor that can determine the rob ustness of a controller is its response to errors in kno wn parameters. As stated earlier the coef cients of lift and drag are calculated from a CFD database. The accurac y of the model depends on accurac y of this data. Rob ustness of a controller can be assessed by introducing errors in the data and checking ho w it performs. The follo wing v ariations ha v e been introduced in the system to check for performance of the system with intrinsic uncertainties: 20% error in C l of Ca vitator 20% error in C d of Ca vitator 20% error in C l of all the Fins. 20% error in C d of all the Fins. 5.3.4 Contr oller Objecti v e In terms of rob ustness, the controller is required to meet v arious performance objecti v es. These objecti v e can be summarized as:The closed-loop system should ha v e a gain mar gin of at least 6 dB.The closed-loop system should ha v e a phase mar gins of at least 45 d e g

PAGE 65

CHAPTER 6 LQR CONTR OL 6.1 LQR Theory The tracking problem, lik e the one gi v en in Equation 5.7 can be solv ed by using a combination of feedback and feedforw ard control [ 10 ]. The Linear Quadratic Re gulator (LQR) problem is to nd an optimal feedback matrix K such that the state-feedback la w u K x minimizes the linear quadratic cost function sho wn in Equation 6.1 Ju 0x T Qxu T Ru2 x T N ud t (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 v ariable x represents the state v ector y is the output v ector and u is the input v ector xAxBu (6.2) yx (6.3) The LQR controller is realized by a constant gain matrix K such that the feedback u K x mak es x go to zero. By a modication to this la w the LQR method can also be used for tracking. The state v ector x is of size n x x 1x 2 x nT (6.4) Let the tracking problem be for the state x 1 to track a step command r 1 The idea is to mak ex 1r 1go to zero using a LQR controller The ne w control la w can be chosen as in Equation 6.17 54

PAGE 66

55 u K x 1r 1 x 2 . x n (6.5) Equation 6.2 can be re written by substituting the ne w control la w xAxBuAxBK x 1r 1 x 2 . x n (6.6) F or simplicity assume that there is only one control, u (this is dif ferent from v elocity u ). The controller K is of size n1 and it can be e xpanded in its elements. K k 1k 2 k n(6.7) Equation 6.6 can be re written by substituting the K in its e xpanded form. xAxBk 1k 2 k n x 1r 1 x 2 . x n AxBK xBk 1 r 1 ABKxBk 1 r 1 (6.8) It should be noted that the command r 1 is a step command. The steady-state dynamics of the system can be obtained from Equation 6.8 x ABKx Bk 1 r 1 (6.9)

PAGE 67

56 The error dynamics can be obtained by subtraction Equation 6.9 from Equation 6.8 xt x ABKxt x(6.10) e ABKe (6.11) where e xt x So, the tracking problem is cast as a re gulator problem. The ne w state v ector is the steady-state error e which is made zero using the re gulator Figure 6.1 sho ws the block diagram for this closed-loop system. It is required that the closed-loop system has an inte grator some where so as to mak e the steady-state error go to 0 [ 10 ]. That is, e has to go to zero rather than e so as to achie v e good tracking. Figure 6.2 sho ws the ne w conguration of a system that has no inte grator and thus an inte grator has to been included during design. Thus, the inte gral of the actual error has to be made to go to zero so as to achie v e a good tracking. e r 1x 1(6.12) The state space equation for the system with this modication can be written. xAxBu er 1x 1r 1C x (6.13) where x 1C x 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 no w has n1 states with state v ector x x 1x 2 x n eT So a ne w formulation of the state space equation can be written, x r A 0C 0 x r B 0 u r 0 1 r 1 (6.14) x A x Bu I r (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 deri v ed. The LQR controller K will be a constant matrix of size n1 as the system no w is of size n1.

PAGE 68

57 K k 1k 2 k nk n1(6.16) Then, the ne w control la w can be written as in Equation 6.17 u K x k 1k 2 k nk n1 r x e k 1k 2 k nx k n1 e K xk I e (6.17) which is represented in Figure 6.2 + k 1 K r 1 xAxB h yx x Figure 6.1 Controller for T racking when Plant has an Inte grator + + k I K xAxBu yxr 1 x C Figure 6.2 Controller for T racking when Plant has no Inte grator

PAGE 69

58 6.2 Contr ol Synthesis The torpedo system does not ha v e an inte grator 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 le v el ight at 75ms1 The linearized dynamics are rst separated into the longitudinal and lateral dynamics as gi v en in T able 5.1 The controls used in longitudinal direction are the ca vitator and 2 ele v ators. The controls in lateral direction are the 2 rudders. It is observ ed that for straight and le v el ight the longitudinal and lateral dynamics are practically decoupled. Once the state and control matrices ha v e been obtained, the main v ariables 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 v alues respecti v ely A higher v alue in Q matrix w ould cause a better tracking. A lar ger R w ould constrain the control surf ace deection. An optimum combination of the matrices is obtained iterati v ely so as to get good tracking with achie v able control deections. The matrices for longitudinal pitch rate tracking are gi v en in Equation 6.18 Q l ongd ia g 000010 R l ongd ia g 54 (6.18) The rst four numbers in the Q l ong correspond to weightings on the four longitudinal states. The y are chosen to be 0. W e do not w ant to restrict the states from changing. This is especially important for weightings on q and Q A weighting on these v ariables w ould 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 v alue of weighting w ould gi v e a better tracking, b ut it is seen that it w ould require v ery high control rates. The rst number in the weighting matrix R l ong 5, corresponds to ca vitator weighting and the number 4 is for

PAGE 70

59 ele v ator weighting. Ele v ator weighting is chosen to be smaller so as to encourage the controller to use more ele v ator than the ca vitator This gi v es a more stable performance. The control matrices obtained for the longitudinal dynamics are gi v en in Equations 6.19 and 6.20 k I r 1118209681 (6.19) K r 00000 00040 0101600000 1419511184 0000000042009950000013981 11308 (6.20) Similar process is in v olv ed in the design of the lateral controller Initially the lateral controller is designed with only four state feedback, and Y is ne glected in the feedback. In this case it is found that the torpedo has high side w ash and de viates considerably from the original path, e v en when the pitch angle is 0. T o a v oid this, Y is included in the feedback states. It is also observ ed that a penalty on the ya w motion causes the controller to command a v ery f ast control surf ace motion. Also, a continuous correction of control surf ace deection is required to pre v ent the ya w motion entirely Thus an optimum combination of the weighting matrices is obtained that w ould pre v ent a v ery high ya w motion b ut w ould still ha v e slo w control surf ace motion. Q l a t dd ia g 00000 1 R l a t dd ia g 10001000 (6.21) The rst 5 numbers correspond to 5 states and the last number is weighting for the error The R l a t d is of dimension 2 as only the rudders are included in the synthesis. The weighting on the rudders is high as it is observ ed that the roll rate is v ery sensiti v e to the rudder deection. The control matrices obtained for the lateral dynamics are gi v en in the Equations 6.22 and 6.23

PAGE 71

60 k I r 0007100071 (6.22) K103 r 000050125300132 0001900000 000050125400132 0002600000 (6.23) The feedback matrix K for lateral dynamics is of size 25, which is sho wn in Equation 6.23 6.3 Nominal Closed-loop Model 6.3.1 Model Figure 6.3 sho ws the eigen v alues for the closed-loop longitudinal and lateral systems. It can be seen that both systems are stable as all the eigen v alues are in the left half of the comple x plane. Also, each of the dynamics has one eigen v alue at the origin, which is introduced due the inte grator in the system. -1000 -800 -600 -400 -200 0 200 -15 -10 -5 0 5 10 15 Real ( l long )Imag ( llong) (a) Longitudinal -1000 -800 -600 -400 -200 0 200 -8 -6 -4 -2 0 2 4 6 8 Real ( l latd )Imag ( llatd) (b) Lateral Figure 6.3 Eigen v alues for the Closed-loop System 6.3.2 Linear Simulations The response of the closed-loop linear system has been sho wn in Figures 6.4 to 6.8 The simulations for lateral and longitudinal systems ha v e been carried out separately as the linear system is decoupled into lateral and longitudinal.

PAGE 72

61 Figure 6.4 sho ws the tracking obtained for a 15 d e gs pitch rate command. The command is achie v ed in 0.17s with an o v ershoot of 3.95% and with no steady-state error Figures 6.5 and 6.6 sho w the control surf ace deections and rates required to achie v e the commands. Though there are some quick deections, the rates are still under the constraints. 0 5 10 15 -20 -10 0 10 20 Pitch Rate(deg s-1)Time(s) command Achieved 0 5 10 15 72.5 73 73.5 74 74.5 75 75.5 time(s)u (m s-1) Figure 6.4 Pitch Command T racking for Linear System : qu 0 5 10 15 -1.5 -1 -0.5 0 0.5 1 time(s)dc (deg) 0 5 10 15 -40 -30 -20 -10 0 10 20 time(s)rate dc(deg s-1) Figure 6.5 Pitch Command T racking for Linear System : d c d c Figure 6.7 sho ws the roll rate tracking obtained for a 15 d e gs roll rate command. The command is achie v ed in 0.53s with no o v ershoot and a 0 steady-state error The v ariation of the ya w rate is also sho wn in the gure and it can be seen that the ya w motion is coupled with the roll. At the end of the roll doublet, the torpedo has a non-zero ya w angle thus changing the direction of motion. The control surf ace deections required for the roll rate command are sho wn in Figure 6.8 The rudder deection is small for reasons e xplained in the ne xt chapter

PAGE 73

62 0 5 10 15 -1 -0.5 0 0.5 1 1.5 time(s)de1 (deg) 0 5 10 15 -20 -10 0 10 20 30 time(s)rate de1 (deg s-1) Figure 6.6 Pitch Command T racking for Linear System : d e 1 d e 1 0 5 10 15 -20 -10 0 10 20 Roll Rate(deg/s)Time(s) Commanded Achieved 0 5 10 15 -2 0 2 4 6 8 time(s)r (deg s-1) Figure 6.7 Roll Command T racking for Linear System : pr 0 5 10 15 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 time(s)dr1 (deg) 0 5 10 15 -0.5 0 0.5 time(s)rate dr1(deg s-1) Figure 6.8 Roll Command T racking for Linear System : d r 1 d r 1

PAGE 74

63 6.3.3 Gain and Phase Mar gins The LQR tracking system sho wn in Figure 6.2 is ob viously more comple x than the system sho wn in Figure 5.4 Thus, the loop gain can be dened in man y w ays in this case. The block diagram can be brok en at dif ferent points so as to simplify it to the form sho wn in Figure 5.4 Figure 6.2 is redra wn in Figure 6.9 which sho ws the possible breakpoints for this system. F or understanding, the output of plant P is di vided into tw o parts, one is the achie v ed v alue of the commanded v ariable ( r a ) and the other is remaining states of the plant P ( x ). The break points are numbered 1 to 3. The system can be brok en at each of these points to gi v e a loop gain. These gains will be named outer -loop, inner -loop and all-loop gains respecti v ely + 1 3 2 + K x r a k 1r 1 yx xAxBu Figure 6.9 Breakpoints for Calculating the Loop-Gain for a T racking Controller Gain and Phase mar gins for each of the abo v e possible break points ha v e been calculated for both the longitudinal and lateral controllers. T able 6.1 lists the gain and phase mar gins for the torpedo with LQR controller that w as obtained in pre vious sections. All mar gins are quite high and meet the desired conditions of 6dB for gain and 45 d e g for phase mar gin. Figures 6.10 to 6.15 sho w the corresponding bode plots for the data gi v en in T able 6.1 Also, the lateral controller is unable to stabilize the unstable spiral mode. Thus the closed-loop system is inherently unstable due to this pole and w ould consequently ha v e ne gati v e gain mar gin. Numerous simulations sho w that the af fect of spiral mode is ne gligible

PAGE 75

64 T able 6.1 Gain and Phase Mar gin with LQR Controller Longitudinal Gain Mar gin(db) Phase Mar gin (de g) 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 57.606(at 20.845 rad/s) Lateral Gain Mar gin(db) Phase Mar gin (de g) 1 22.964(at 0 rad/s) 2 250.51 (at 0 rad/s) 3 50.36 (at 0 rad/s) i.e., the time to double for the instability is considerably lar ger than the maneuv ering time of the torpedo. So, the closed-loop system model is reduced by remo ving the spiral mode from the model. The gain and phase mar gins in T able 6.1 are for this reduced-order system and reect the rob ustness of the dominant dynamics. Bode Diagram (rad/sec)Phase (deg) Magnitude (dB) -100 -50 0 Gm = 21.056 dB (at 47.498 rad/sec), Pm = 64.846 deg (at 9.0625 rad/sec) 10 0 10 1 10 2 10 3 -270 -225 -180 -135 -90 Figure 6.10 Gain and Phase Mar gin: Longitudinal Outer -loop 6.4 P erturbed Closed-loop Model A perturbed system model is formed by adding an error to the v alues of coef cients of lift and drag for the ns and ca vitator Ne w v alues of trim deection are obtained for the perturbed model and thus a ne w set of A and B matrices is obtained. T ables 6.4 to 6.9

PAGE 76

65 Bode Diagram (rad/sec)Phase (deg) Magnitude (dB) -50 0 50 Gm = 327.87 dB (at 0 rad/sec), Pm = 77.118 deg (at 25.925 rad/sec) 10 0 10 1 10 2 10 3 -180 -135 -90 -45 0 45 90 Figure 6.11 Gain and Phase Mar gin: Longitudinal Inner -loop Bode Diagram (rad/sec)Phase (deg) Magnitude (dB) -100 0 100 Gm = Inf, Pm = 57.606 deg (at 20.845 rad/sec) 10 0 10 2 -180 -135 -90 Figure 6.12 Gain and Phase Mar gin: Longitudinal All-loop sho w the percentage v ariation of the elements of A and B matrices for a 20% change in coef cients of lift and drag of ca vitator and ns. Fe w elements in the state and control matrices change. In most cases, the change in elements of A and B matrices is a linear function of the change in a coef cient. F or e xample, in T able 6.4 there are 8 terms that sho w a v ariation due to a 20% v ariation in coef cient of lift of the ca vitator The term A(3,1) sho ws a lar ge v ariation b ut its numerical v alue is ne gligible. The term A(3,2) sho ws a 34% v ariation b ut this term is also small compared to other terms. Remaining terms in the matrix sho w v ery small v ariation. Some terms in the B matrix sho w a 20% v ariation.

PAGE 77

66 Bode Diagram (rad/sec)Phase (deg) Magnitude (dB) -50 -40 -30 -20 Gm = 22.964 dB (at 0 rad/sec), Pm = Inf 10 0 10 1 10 2 10 3 90 135 180 Figure 6.13 Gain and Phase Mar gin: Lateral Outer -loop Bode Diagram (rad/sec)Phase (deg) Magnitude (dB) -50 -40 -30 -20 Gm = 250.51 dB (at 0 rad/sec), Pm = Inf 10 0 10 1 10 2 10 3 -90 -45 0 45 90 Figure 6.14 Gain and Phase Mar gin: Lateral Inner -loop Thus some terms in controllability matrix change considerably This w ould mean that for an error in these coef cients, the response w ould sho w some dif ference in control surf ace deection. As it is observ ed that the closed-loop system has good gain and phase mar gins, this ef fect on B matrix should not be of much concern. 6.4.1 Model Figure 6.16 sho ws the eigen v alues for the perturbed closed-loop longitudinal and lateral systems. An error of -20% is included in the v alue of coef cient of lift for the ns. It can

PAGE 78

67 T able 6.2 Percentage V ariation in A Matrix due to 20% V ariation in cl c u w q Q v p r F Y u 0.46 0.62 w 5.52 6.86 1.58 q 1.05e5 34.8 13.58 Q v p r F Y T able 6.3 Percentage V ariation in B Matrix due to 10% V ariation in cl c d c d e 1 d e 2 d r 1 d r 2 u w 20 q 20 Q v p r F Y T able 6.4 Percentage V ariation in A Matrix due to 20% V ariation in cd c u w q Q v p r F Y u 10 -6 7.6 w 1.54 0.36 q 658 7.8 3 Q v 2 20 0.4 p 2 -15.6 r -10 0.8 8.4 F Y

PAGE 79

68 T able 6.5 Percentage V ariation in B Matrix due to 20% V ariation in cd c d c d e 1 d e 2 d r 1 d r 2 u 20 w q 0.12 Q v p r F Y T able 6.6 Percentage V ariation in A Matrix due to 20% V ariation in cl f in u w q Q v p r F Y u 1.22 0.64 w 14.4 10.0 -0.9 q -1e5 -60 3 Q v 16.2 -1.2 p 19 36.0 r 25.4 0.94 10.2 F Y T able 6.7 Percentage V ariation in B Matrix due to 20% V ariation in cl f in d c d e 1 d e 2 d r 1 d r 2 u w 20 20 q 20 20 Q v 20 20 p 20 20 20 20 r 20 20 F Y

PAGE 80

69 Bode Diagram (rad/sec)Magnitude (dB) -80 -70 -60 -50 Gm = 50.36 dB (at 0 rad/sec), Pm = Inf 10 0 10 1 10 2 10 3 90 135 180 Phase (deg) Figure 6.15 Gain and Phase Mar gin: Lateral All-loop T able 6.8 Percentage V ariation in A Matrix due to 20% V ariation in cd f in u w q Q v p r F Y u 9.4 24.0 12.4 w 1.34 -0.12 q -68 -2.6 0.4 Q v 20 1.76 -2.3 -0.13 p -2.3 1.12 -0.62 r 2.74 18.26 1.12 F Y be seen that the longitudinal dynamics sho w some perturbation in the damping while the lateral system relati v ely unchanged. 6.4.2 Linear Simulations Figures 6.17 to 6.19 sho w the response of the perturbed system for a 15 d e gs pitch rate doublet command. The perturbation to the system is a 20% error in the v alue of the coef cient of lift of the ns. It can be seen that the performance criteria are met e v en in the case of the perturbed system. It can also be observ ed that the o v ershoot is increased for the perturbed system. The performance is achie v ed at the cost of small perturbations in

PAGE 81

70 T able 6.9 Percentage V ariation in B Matrix due to 20% V ariation in cd f in d c d e 1 d e 2 d r 1 d r 2 u 20 20 w q 1.36e-2 Q v p r 20 20 12.6e-3 2.6e-3 F Y -1000 -800 -600 -400 -200 0 -15 -10 -5 0 5 10 15 RealImag ( llong) (a) Longitudinal -1000 -800 -600 -400 -200 0 200 -8 -6 -4 -2 0 2 4 6 8 RealImag ( llatd) (b) Lateral Figure 6.16 Eigen v alues for the Perturbed Closed-loop System: 20% Error in cl f in other states of the system. As the control ef fecti v eness of the control surf aces is changed, the amount of control surf ace deection is also changed by a constant f actor Figures 6.20 to 6.21 sho w the response of the perturbed system for a 15 d e gs roll rate doublet command. The perturbation to the system is a 20% error in the v alue of the coef cient of lift of the ns. It can be seen that the performance criteria are met e v en in case of perturbed system. In this case also, it can be observ ed that there is a perturbation in other states.

PAGE 82

71 0 5 10 15 -20 -15 -10 -5 0 5 10 15 20 Pitch Rate(deg s-1)Time(s) No Uncertainty +20% in cl f -20% in cl f 0 5 10 15 72 72.5 73 73.5 74 74.5 75 75.5 time(s)u (m s-1) No Uncertainty +20% in clf-20% in clf Figure 6.17 Pitch Command T racking for Perturbed Linear System : qu 0 5 10 15 -1.5 -1 -0.5 0 0.5 1 time(s)dc (deg) No Uncertainty +20% in cl f -20% in cl f 0 5 10 15 -40 -30 -20 -10 0 10 20 time(s)rate dc(deg s-1) No Uncertainty +20% in clf-20% in clf Figure 6.18 Pitch Command T racking for Perturbed Linear System : d c d c 0 5 10 15 -1 -0.5 0 0.5 1 1.5 time(s)de1 (deg) No Uncertainty +20% in cl f -20% in cl f 0 5 10 15 -20 -10 0 10 20 30 time(s)rate de1 (deg s-1) No Uncertainty +20% in clf-20% in clf Figure 6.19 Pitch Command T racking for Perturbed Linear System : d e 1 d e 1

PAGE 83

72 0 5 10 15 -20 -15 -10 -5 0 5 10 15 20 Roll Rate(deg s-1)Time(s) No Uncertainty +20% in cl f -20% in cl f 0 5 10 15 -1 0 1 2 3 4 5 6 7 time(s)r (deg s-1) No Uncertainty +20% in clf-20% in clf Figure 6.20 Roll Command T racking for Perturbed Linear System : pr 0 5 10 15 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 time(s)dr1 (deg) No Uncertainty +20% in cl f -20% in cl f 0 5 10 15 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0.2 0.25 time(s)rate dr1(deg s-1) No Uncertainty +20% in clf-20% in clf Figure 6.21 Roll Command T racking for Perturbed Linear System : d r 1 d r 1 6.4.3 Gain and Phase Mar gins T able 6.10 lists the gain and phase mar gins for the perturbed closed-loop system. The perturbed system also has good gain and phase mar gins. Comparing the v alues with T able 6.1 it can be seen that there are small changes in the v alues e xcept for the lateral all-loop. The last v alue is increased to sho wing an impro v ement for the perturbed system. From the analysis of the perturbed closed-loop system it can be said that the linear model is rob ust to v arious uncertainties in the system.

PAGE 84

73 T able 6.10 Gain and Phase Mar gin for Perturbed Closed-loop System: 20% error in cl f in Longitudinal Gain Mar gin(db) Phase Mar gin (de g) 1 21.193(at 49.599 rad/s) 68.981(at 8.7966 rad/s) 2 320.33(at 0 rad/s) 77.605(at 25.925 rad/s) 3 60.305(at 22.552 rad/s) Lateral Gain Mar gin(db) Phase Mar gin (de g) 1 24.391(at 0 rad/s) 2 278.23 (at 0 rad/s) 3 (at 0 rad/s)

PAGE 85

CHAPTER 7 NONLINEAR SIMULA TIONS The controller for longitudinal and lateral dynamics ha v e been obtained separately That is, the longitudinal controller is to achie v e a required pitch rate and the lateral controller achie v es a gi v en roll rate. Once the controllers ha v e been found for linear systems, the y are emplo yed with the nonlinear torpedo model and the performance is check ed using numerical simulation. Figure 7.1 sho ws the complete nonlinear simulation model for the torpedo with the LQR controller The nonlinearity in the model is pro vided by both the aerodynamics and control surf ace constraints. These constraints, gi v en in T able 5.2 are not directly included in the linear model. So it is important to nd their ef fects on the nonlinear simulation. It should be noted that there is a constraint on thrust, b ut thrust is assumed to be constant with time. The thrust constraint is used to nd a trim v alue for v arious trajectories. 7.1 Nonlinear Simulations f or Nominal System The simulations ha v e been carried out for v arious commands with the nominal system. The `Nominal system' is the nonlinear torpedo model assuming that it is completely accurate. No uncertainty is included in the model. The simulations sho w good tracking response while meeting all the performance criteria. The response of the v ehicle to a longitudinal command is simulated and sho wn in Figures 7.2 to 7.5 These gures sho w the response for a pitch rate doublet of 15 d e gs The rise time for the pitch rate command of 15 d e gs is 0.18s and there is an o v ershoot of 11.53%. The steady-state error is .8%. The controller is able to command pitch rates as high as 30 d e gs It is observ ed that the v ehicle motion is conned to longitudinal plane 74

PAGE 86

75 Low Pass Filter 4 Control Rates 3 Controls 2 Position 1 States In1 state pos state rates K*u size1 K*u size command_roll roll command In1 long latd p q re-arrange state K*u r2d6 K*u r2d1 command_pitch pitch command 1 s int1 K*u d2r1 K*u d2r In1 In2 Out1 Out2 change units K*u change size of output K*u change sige of output simt To Workspace MATLAB Function NL Equations Torpedo K*u Klongfwd K*u Klongbck K*u Klatdfwd K*u Klatdbck 1 s Int_state 1 s Int_pos 1 s Int du/dt Derivative in Out Control Limiters 0.3438 Clock In1 Controls Actuator state longitudinal 14 9 9 9 9 5 5 5 5 5 5 5 5 4{4} 4{4} 6 6 2 2 2 5 5 5 5 9 12 3 2 9 3 2{2} 2{2} 3 5 5 5 5 5 5 5 5 5 Figure 7.1: Complete Nonlinear Simulation with LQR Controller

PAGE 87

76 0 5 10 15 -5 0 5 x 10 -5 time(s)p (deg s-1) 0 5 10 15 -20 -10 0 10 20 time(s)q (deg s-1) Figure 7.2 Pitch Command T racking : pq 0 5 10 15 -1 -0.5 0 0.5 1 time(s)dc (deg) 0 5 10 15 -30 -20 -10 0 10 20 time(s)rate dc (deg s-1) Figure 7.3 Pitch Command T racking : d c d c only This sho ws that the controller allo ws pure longitudinal motion to be uncoupled from the lateral motion. 0 5 10 15 0 0.5 1 1.5 2 time(s)de1 (deg) 0 5 10 15 -15 -10 -5 0 5 10 15 20 time(s)rate de1 (deg s-1) Figure 7.4 Pitch Command T racking : d e 1 d e 1

PAGE 88

77 0 5 10 15 -5 0 5 10 15 20 x 10 -7 time(s)dr1 (deg) 0 200 400 600 800 -50 0 50 -300 -200 -100 0 100 x(m) y (m)z (m)starting point Figure 7.5 Pitch Command T racking : d r 1 xyzT rajectory The response of the v ehicle to a lateral, roll rate, command is sho wn in Figures 7.6 to 7.10 A roll rate command of 15 d e gs is achie v ed in .52s with an o v ershoot of 0% and a steady-state error of 0.09%. The controller is able to command a roll rate motion of as high as 50 d e gs before a saturation of control surf ace rate is reached. It is observ ed 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 deection for a roll rate command is small. This is e xpected as the terms corresponding to roll rate from rudder are an order of 3 times lar ger than the terms corresponding to pitch rate from ele v ators. It is assumed that the control surf ace deection is achie v able. 0 5 10 15 -20 -10 0 10 20 time(s)p (deg s-1) 0 5 10 15 -2.5 -2 -1.5 -1 -0.5 0 0.5 time(s)q (deg s-1) Figure 7.6 Roll Command T racking: pq Figures 7.11 to 7.15 sho w the response of the torpedo for a combined roll and pitch rate command, similar to a windup turn. In this case, the torpedo is gi v en a 5 d e gs roll

PAGE 89

78 0 5 10 15 0 0.2 0.4 0.6 0.8 1 time(s)dc (deg) 0 5 10 15 -0.5 0 0.5 1 1.5 2 2.5 3 time(s)rate dc (deg s-1) Figure 7.7 Roll Command T racking: d c d c 0 5 10 15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 time(s)de1 (deg) 0 5 10 15 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 time(s)rate de1 (deg s-1) Figure 7.8 Roll Command T racking: d e 1 d e 1 0 5 10 15 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 time(s)dr1 (deg) 0 5 10 15 -0.5 0 0.5 time(s)rate dr1 (deg s-1) Figure 7.9 Roll Command T racking: d r 1 d r 1

PAGE 90

79 0 500 1000 1500 0 500 1000 0 500 1000 x(m) y(m)z(m) Starting Point Figure 7.10 Roll Command T racking:xyzT rajectory rate command from 2 to 12 seconds, 5 d e gs pitch rate command from 12 to 22 seconds, -5 d e gs pitch rate command from 22 to 32 seconds and then -5 d e gs roll rate command from 32 to 42 seconds. As the v ehicle motion is a little dif ferent from the actual trim, it can be seen the v ehicle has considerable side w ash. Despite this side w ash, the controllers gi v e a good tracking performance. The rise times for roll and pitch commands are .5s and 0.22s respecti v ely The o v ershoot for the same are 0% and 0% respecti v ely 0 10 20 30 40 50 -6 -4 -2 0 2 4 6 time(s)p (deg s-1) 0 10 20 30 40 50 -6 -4 -2 0 2 4 6 time(s)q (deg s-1) Figure 7.11 Roll & Pitch Command T racking: pq 7.2 Nonlinear Simulations f or P erturbed System The performance of the controllers is studied using the simulation with a perturbed system model. An error is assumed in the v alues of v arious coef cients and a correction f actor is added. Response of the closed-loop nonlinear system is not much af fected by the v ariations in coef cients of lift and drag. It is observ ed that the controller commands the

PAGE 91

80 0 10 20 30 40 50 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 time(s)dc (deg) 0 10 20 30 40 50 -8 -6 -4 -2 0 2 4 time(s)rate dc (deg s-1) Figure 7.12 Roll & Pitch Command T racking: d c d c 0 10 20 30 40 50 0.2 0.4 0.6 0.8 1 time(s)de1 (deg) 0 10 20 30 40 50 -4 -2 0 2 4 6 time(s)rate de1 (deg s-1) Figure 7.13 Roll & Pitch Command T racking: d e 1 d e 1 0 10 20 30 40 50 -0.01 -0.005 0 0.005 0.01 time(s)dr1 (deg) 0 10 20 30 40 50 -0.04 -0.02 0 0.02 0.04 time(s)rate dr1 (deg s-1) Figure 7.14 Roll & Pitch Command T racking: d r 1 d r 1

PAGE 92

81 0 500 1000 1500 -1000 0 1000 2000 -1000 0 1000 2000 x(m) y(m)z(m) Starting Point Figure 7.15 Roll & Pitch Command T racking:xyzT racking 0 5 10 15 20 74.9 75 75.1 75.2 75.3 75.4 75.5 75.6 75.7 timeu(m/s) No Uncertainty 20% in cl f -20% in cl f 0 5 10 15 20 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 timew(m/s) No Uncertainty 20% in cl f -20% in cl f Figure 7.16 Response for 20% V ariation in cl f in : uw system to a ne w trim state which is also a straight and le v el ight, with change in speed and control deections. After that, the system follo ws a pitch or roll command as well as before. Figures 7.16 to 7.19 sho w the response for one such case. In this case a roll doublet is commanded to the system, and there is an error of20% in the v alue of cl f in It can be clearly seen that the v ehicle has gone to another trim state and then it follo ws the command equally well. There is almost no change in the trajectory of the v ehicle. The control surf ace deections are similar with a constant of fset. Such response has also been check ed for other cases. The af fect of error is similar in all cases. By abo v e analysis it can be said that the LQR controller designed for the torpedo for pitch and roll rate tracking commands is f airly stable and can be e xpected to achie v e good performance for the real torpedo.

PAGE 93

82 0 5 10 15 20 -20 -15 -10 -5 0 5 10 15 20 timep(deg/s) No Uncertainty 20% in cl f -20% in cl f 0 5 10 15 20 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 timeq(deg/s) No Uncertainty 20% in cl f -20% in cl f Figure 7.17 Response for 20% V ariation in cl f in : pq 0 5 10 15 20 -0.5 0 0.5 1 timecavitator (deg) No Uncertainty 20% in cl f -20% in cl f 0 5 10 15 20 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 timerudder (deg) No Uncertainty 20% in cl f -20% in cl f Figure 7.18 Response for 20% V ariation in cl f in : d cd r 1 0 500 1000 1500 0 500 1000 0 500 1000 x(m) y (m)z (m) No Uncertainty 20% in cl f -20% in cl f Figure 7.19 Response for 20% V ariation in cl f in :xyzT rajectory

PAGE 94

CHAPTER 8 CONCLUSION 8.1 Summary A dynamical model for a superca vitating v ehicle has been obtained. The v ehicle is found to be open-loop unstable, and a controller for stabilizing the pitch and roll rate motion has been obtained. The LQR controller sho ws good tracking performance for the v ehicle and all the control objecti v es are met. The controller is also found to be rob ust to errors in ca vity prediction and v elocity changes. This rob ustness is further demonstrated by the f act that the closed-loop system has high gain and phase mar gins. 8.2 Futur e W ork The dynamical analysis of the v ehicle has been deri v ed with an assumption that the ca vity shape is x ed. The open-loop ca vity dynamics need to be modeled and included in the synthesis. Rob ust control methodologies lik e -synthesis can be applied to obtain a more rob ust controller A rob ust control design could include the uncertainties in the model during the control synthesis. The LQR controllers obtained are typically kno wn as the `inner -loop' controllers. An outer -loop controller is also needed for guidance and na vigation. The idea is that the outer loop controller can be modeled for tracking the trajectory in space, based on the closed-loop dynamics of the inner -loop model. 83

PAGE 95

APPENDIX A REFERENCE FRAMES AND R O T A TION MA TRICES x 2 y 2 x 3 y 3 x 1 x 2 q Figure A.1 Rotation of Frames Figure A.1 sho ws tw o frames Xx 1 x 2 x 3and Yy 1 y 2 y 3. Y is rotated from X by an angle q about x-axis. Thus the basis v ectors of frame Y can be written in terms of basis v ectors of X frame. y 2x 2 cosq x 3 sinqy 3 x 2 sinq x 3 cosq(A.1) This relation can also be e xpressed in terms of matrices. y 1 y 2 y 3 n r 1 0 0 0 cosqsinq0sinqcosq x 1 x 2 x 3 n(A.2) This w as a case of simple rotation. The matrix abo v e in square brack ets is kno wn as the rotation matrix from X to Y and is represented as X Y The rotation matrix can be generalized for a case when the tw o reference frames are arbitrarily oriented. 84

PAGE 96

85 X Y r y 1x 1 y 1x 2 y 1x 3 y 2x 1 y 2x 2 y 2x 3 y 3x 1 y 3x 2 y 3x 3 (A.3) where ( ) means the dot product of the tw o v ectors. Thus, YX YX (A.4)

PAGE 97

APPENDIX B NUMERICAL TECHNIQ UES B.1 Inter polation of F or ce Data This section describes the numerical technique used to obtain the v alues of coef cients of lift and drag for ca vitator and ns. B.1.1 Extrapolation Scheme F or a better result, a quadratic interpolation/e xtrapolation scheme is used. Thus 3 points w ould be required to obtain an interpolated or e xtrapolated data v alue. Figure B.1 sho ws the shape functions used for one dimensional interpolation. Say pointsx i1x ix i1are used to nd the v alue of a function f at point x The v alue of fxw ould be gi v en by a parameter a and the three shape function N 1, N 2 and N 3. N 11 2 x i1x ix i1 x ix i1a x i1x i1 x ix i1a 2 N 2 x i1x i12 x i1x i x ix i1a x i1x i12 x i1x i x ix i1a 2 N 3 x i1x i x ix i1a x i1x i1 x ix i1a 2 (B.1) where, the v alue of a shape function can be obtained by nding the v alue of a axx i1 x i1x i1 (B.2) a 01f or x x i1x i1a0 f or xx i1 a1 f or xx i1 Thus a 01for interpolation and it is greater than 1 or less than 0 for e xtrapolation. 86

PAGE 98

87 0 0.2 0.4 0.6 0.8 1 -0.5 0 0.5 1 AlphaShape FunctionN1 N2 N3 X i-1 X X i i+1 Figure B.1 Shape Function for One Dimensional Quadratic Scheme fx N 1fx i1 N 2fx i N 3fx i1(B.3) This method can be e xtended for 2D and 3D as in case for ca vitator and ns respecti v ely B.1.2 Ca vitator The coef cients of lift ( cl c ) and drag ( cd c ) for the ca vitator are functions of half angle ( h a ) of ca vitator cone and angle of attack for ca vitator ( a c ). The CFD data [ 5 ] is a v ailable for combination of points gi v en in T able B.1 Equation B.3 can be e xtended for 2D ca vitator fa ch a 3 i1 3 j1 N1iN2jfa ci h aj (B.4) a ciV alue of a c at i t h node h ajV alue of h a at j t h node N1ii t h Shape function for a c N2jj t h Shape function for h a B.1.3 Fins The coef cients of lift ( cl f in ) and drag ( cd f in ) for the ns are functions of angle of attack ( a f ) for n, immersion ( S f ) and sweepback angle ( q f ). The CFD data is a v ailable for combination of points gi v en in T able B.2 Equation B.3 can be e xtended for 3D n.

PAGE 99

88 T able B.1 Grid F or Experimental Ca vitator Data Half Angle ( h a de g) 15, 30, 45, 60, 75, 90 Angle of Attack ( a c de g) 0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15 T able B.2 Grid F or Experimental Fin Data Immersion ( S f ) 0.1,0.3,0.5,0.7,0.9 Sweepback ( q f de g) 0,15,30,45,60,70 Angle of Attack ( a f de g) 0,1,2,3,4,5,6,7,8,9,10,12,15 Data not A v ailable for 1. S0 & S01 & q0 2. S01 & S03 & q30 3. S03 & S05 & q45 4. S05 & S07 & q60 5. S07 & S09 & q70 6. S09 & S1 & q0 7. a0 8. a15 fS fq fa f 3 i1 3 j1 3 k1 N1iN2jN3kfS fi q fj a fk (B.5) SiV alue of S f at i t h node q fjV alue of q f at j t h node a fkV alue of a f at k t h node N1ii t h Shape function for S f N2jj t h Shape function for q f N3jk t h Shape function for a f B.2 Numerical Linearization Numerical linearization can be done by the `linmod' command in the Matlab Simulink toolbox. This can also be done by noting that, the terms in the A and B matrices are the deri v ati v es of state rates with respect to states and controls. F or e xample, suppose x 0 and u 0

PAGE 100

89 represent the state and control v alues at trim. It should be noted that x 0 is a 91 (e xcluding the positionsxyz) v ector and u 0 is 51 (ca vitator and four ns). x 0 u 0 w 0 q 0 Q 0 v 0 p 0 r 0 F 0 Y 0T u 0 d c 0 d e 1 0 d e 2 0 d r 1 0 d r 2 0T (B.6) The equations of motion are of the form as in equation B.7 xfxu(B.7) where the function f is a v ector function ha ving 9 outputs and there are 9 states. The code for nonlinear equation of motion, tak es x and u as inputs and gi v e the v alue of x as output. Let e dene a v ery small change. No w the element Aijcan be calculated as in equation B.8 Aij fx 0ej u 0ifx 0u 0i e 1ij9 (B.8) where, ejmeans a matrix of size x 0 with all zeros e xcept j t h element, which is equal to e and f i represents the i t h element of v ector f An element Bijalso can be obtained in a similar w ay Bij fx 0u 0ej ifx 0u 0i e 1i91j5 (B.9) where, ejmeans a matrix of size u 0 with all zeros e xcept j t h element, which is equal to e

PAGE 101

REFERENCES [1] S. Ashle y W arpdrive Underwater 2002, http://www .diodon349.com/K urskMemorial/W arpdri v e underw ater .htm accessed: March 2002. [2] M. Billet, Cavitation Applied Research Laboratory at the Pennsylv ania State Univ ersity 2000, http://www .arl.psu.edu/areas/ca vitation/ca vitation.h tml accessed: March 2002. [3] D. R. Stinebring, M. L. Billet, J. W Lindau and R. F K unz, “De v eloped Ca vitation-Ca vity Dynamics, ” VKI Special Course on Superca vitating Flo ws, http://www .arl.psu.edu/areas/compmech/publications.h tml February 2001 [4] A. May W ater Entry and Cavity-Running Behavior of Missiles Arlington, V A, SEAHA C TR 75-2, Na v al Sea Systems Command, 1975. [5] N. Fine, Six De gr ee-of-F r eedom F in F or ces for the ONR Super cavitating T est Bed V ehicle Anteon Corporation, 2000, http://www .anteon.com accessed: September 2000. [6] S. S. K ulkarni and R. Pratap, “Studies on Dynamics of a Superca vitating Projectile, ” Applied Mathematical Modeling v ol. 24, pp. 113–129, 2000. [7] R. Rand, R. Pratap, D. Ramani, J. Cipolla and I. Kirschner “Impact Dynamics of a Superca vitating Underw ater Projectile, ” in Pr oceedings of the Thir d International Symposium on P erformance Enhancement for Marine Applications Ne wport, RI, pp. 215–223, 1997, http://tam.cornell.edu/randpdf/Rand pub .html accessed: March 2002. [8] J. Dzielski and A. J. K urdila, “ A Benchmark Control Problem for Superca vitating V ehicles and an Initial In v estigation of Solutions, ” accepted for publication in J ournal of V ibr ation and Contr ol [9] R. C. Nelson, Flight Stability and A utomatic Contr ol Boston, MA, McGra w Hill, 1997. [10] K. Ogata, Modern Contr ol Engineering Upper Saddle Ri v er NJ, Prentice Hall, 2002. 90

PAGE 102

BIOGRAPHICAL SKETCH Anukul Goel w as born in Luckno w India, on March 3 rd 1978, and raised in Hyderabad, India. Anukul attended the Indian Institute of T echnology located in Mumbai, India, where he recei v ed a Bachelor of T echnology de gree in aerospace engineering in 2000. Since 2000, Anukul has attended the Colle ge of Engineering at the Uni v ersity of Florida, Gainesville, to pursue his M. S. in aerospace engineering. During this time he w ork ed as a teaching assistant and a research assistant in the Mechanical and Aerospace Engineering Department on a part-time basis. His research interests include controls and dynamics and optimization. 91


Permanent Link: http://ufdc.ufl.edu/UFE0000537/00001

Material Information

Title: Control strategies for supercavitating vehicles
Physical Description: Mixed Material
Creator: Goel, Anukul ( Author, Primary )
Publication Date: 2002
Copyright Date: 2002

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0000537:00001

Permanent Link: http://ufdc.ufl.edu/UFE0000537/00001

Material Information

Title: Control strategies for supercavitating vehicles
Physical Description: Mixed Material
Creator: Goel, Anukul ( Author, Primary )
Publication Date: 2002
Copyright Date: 2002

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0000537:00001


This item has the following downloads:


Full Text









CONTROL STRATEGIES FOR SUPERCAVITATING VEHICLES


By

ANUKUL GOEL





















A THESIS PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE


UNIVERSITY OF FLORIDA


2002













ACKNOWLEDGMENT S


I would like to express my sincere gratitude to my committee chairman, Dr. Andrew

Kurdila, for his invaluable guidance throughout the course of this project. I would also like

to thank him for giving me this opportunity to work on such a fascinating project.

I would also like to thank my committee cochair, Dr. Richard C. Lind, for his invaluable

guidance and inspiration throughout the project.

I would also like to show my sincere appreciation to Dr. John Dzielski, Dr. Nor-

man Fitz-Coy and Jammulamadaka Anand Kapardi for their valuable contributions to this

project. I would also like to express my gratitude to all the members, past and present, of

the Supercavitation Project.

I would also like to thank the Office of Naval Research for the support of the research

grant for the project.

On a personal note, I would like to thank all my friends and family members whose

support helped me to aim towards my goals.













TABLE OF CONTENTS


page


ACKNOWLEDGMENTS .

LIST OF TABLES .


. . .. ii

. . . . v i


LIST OF FIGURES ............

ABSTRACT . . . .

CHAPTER


1 INTRODUCTION .


. . . . 1


1.1 Cavitation .... .. . . ...
1.2 Types of Supercavitating Projectiles .
1.3 Related Research .............
1.4 Overview of this Thesis . .....

2 CONFIGURATION OF VEHICLE .

2.1 Cavitator ........
2 .2 F in s . . . . .
2.3 Maneuvering ..............


3 NONLINEAR EQUATIONS OF MOTION .....

3.1 Kinematic Equations of Motion ................
3.1.1 Orientation of the Torpedo ...............
3.1.2 Orientation of the Cavitator .............
3.1.3 Orientation of Fins ...................
3.1.4 Angle of Attack and Sideslip .......
3.1.5 Kinematic Equations ..................
3.2 Dynamic Equations of Motion ................
3.2.1 Forces on Cavitator .. ................
3.2.2 Forces on Fins . . . . . .
3.2.3 Gravitational Forces .. ...............
3.2.4 Equations of M otion .. ...............

4 LINEARIZATION .........


4.1 Linearization .......


2
4
5

7
. . 1
. . 2


. . 34









4.1.1 Need for Linearization ................. 34
4.1.2 Generic Form of Equations of Motion . . .. . 35
4.1.3 Small Disturbance Theory . . . . . 35
4.1.4 Stability and Control Derivatives . . .. . 37
4.2 State Space Representation . . . . . . 42

5 CONTROL DESIGN SETUP . . . . . 46

5.1 Open-loop Performance . . . . . . . 47
5.2 Closed-Loop Problem . . . . . . . 50
5.3 Robustness of the Controller . . . . . . 52
5.3.1 G ain M argin . . . . . . . 52
5.3.2 Phase M argin . . . . . . . 52
5.3.3 Uncertainty In Parameters . . . . . 53
5.3.4 Controller Objective . . . . . . 53

6 LQR CONTROL . . . . . . 54

6.1 LQR Theory.......... ... . . ............ 54
6.2 Control Synthesis . . . . . . . 58
6.3 Nominal Closed-loop Model . . . . . . 60
6.3.1 M odel . . . . . . . . 60
6.3.2 Linear Simulations . . . . . . 60
6.3.3 Gain and Phase Margins . . . . . .... 63
6.4 Perturbed Closed-loop Model . . . . . . 64
6.4.1 M odel . . . . . . . . 66
6.4.2 Linear Simulations . . . . . . 69
6.4.3 Gain and Phase Margins . . . . . .... 72

7 NONLINEAR SIMULATIONS . . . . 74

7.1 Nonlinear Simulations for Nominal System . . . . 74
7.2 Nonlinear Simulations for Perturbed System . . . . 79

8 CONCLUSION . . . . . . 83

8.1 Sum m ary . . . . . . . . 83
8.2 Future Work.......... ... . . ............ 83

APPENDIX

A REFERENCE FRAMES AND ROTATION MATRICES . 84

B NUMERICAL TECHNIQUES. . . . . 86

B.1 Interpolation of Force Data . . . . . . 86
B.I.1 Extrapolation Scheme .. . . . . . 86









B .1.2 C avitator . . . . . . . 87
B .1.3 F in s . . . . . . . . 87
B.2 Numerical Linearization .................. ........ 88

REFERENCES . . . . . . 90

BIOGRAPHICAL SKETCH . . . . . 91














LIST OF TABLES


Table

5.1

5.2

6.1

6.2

6.3

6.4

6.5

6.6

6.7

6.8

6.9

6.10

B.1

B.2


page

. . 46

. . 5 1

. . 64

. . 67

. . 67

. . 67

. . 68

. . 68

. . 68

. . 69

. . 70


Gain and Phase Margin for Perturbed Closed-loop System: 20% error in clfin

Grid For Experimental Cavitator Data .. .................

Grid For Experimental Fin Data .. ...................


Control Param eters . . .. . . .

Control Constraints . . . . . .

Gain and Phase Margin with LQR Controller . ......

Percentage Variation in A Matrix due to 20% Variation in cc .

Percentage Variation in B Matrix due to 10% Variation in c .

Percentage Variation in A Matrix due to 20% Variation in cdc .

Percentage Variation in B Matrix due to 20% Variation in cdc .

Percentage Variation in A Matrix due to 20% Variation in clfin .

Percentage Variation in B Matrix due to 20% Variation in clfin .

Percentage Variation in A Matrix due to 20% Variation in cdfin .

Percentage Variation in B Matrix due to 20% Variation in cdfin


1













LIST OF FIGURES

Figure page

1.1 Tip Vortex Cavitation . . . . . . . 2

1.2 Formation of Cavity. . . ................... 3

1.3 Supercavitating Vehicle . . . . . . . 4

2.1 Supercavitating Vehicle . . . . . . . 7

2.2 Cavitator and Fins . . . . . . . 9

3.1 Body-fixed and Inertial Frames .. . . . . . 11

3.2 Principle Planes of Symmetry for the Torpedo . . .. . 12

3.3 Euler Angles of Rotation ................... . . 12

3.4 Cavitator Reference Frame . . . . . . 13

3.5 Rudder and Fin Reference Frames . . . . . 14

3.6 Rudder 1 Fin Reference Frames . . . . . 16

3.7 Rudder 2 Fin Reference Frames . . . . . 16

3.8 Elevator 1 Fin Reference Frames . . . . . .... 17

3.9 Elevator 2 Fin Reference Frames . . . . . 17

3.10 Angle of Attack (a) and Sideslip () . . . . . 18

3.11 Cavitator: (a) Angle of Attack and Sideslip and (b) Hydrodynamic Forces 25

3.12 Fin G eom etry . . . . . . . . 28

5.1 Simulink Model for Open Loop Simulation . . . . 48

5.2 Open-Loop Response for Torpedo: w,p,q . . . . 48

5.3 Variation of Eigenvalues with Change in Velocity . . . 50

5.4 Loop G ain . . . . . . . . 53

6.1 Controller for Tracking when Plant has an Integrator. . . . 57

6.2 Controller for Tracking when Plant has no Integrator. . . . 57









6.3 Eigenvalues for the Closed-loop System . . .. . 60

6.4 Pitch Command Tracking for Linear System : q,u . . . 61

6.5 Pitch Command Tracking for Linear System : c, c . .. 61

6.6 Pitch Command Tracking for Linear System : ei, e . . .. 62

6.7 Roll Command Tracking for Linear System : p,r . . . 62

6.8 Roll Command Tracking for Linear System : 6r1, 6 . . 62

6.9 Breakpoints for Calculating the Loop-Gain for a Tracking Controller . 63

6.10 Gain and Phase Margin: Longitudinal Outer-loop . . . 64

6.11 Gain and Phase Margin: Longitudinal Inner-loop . . . 65

6.12 Gain and Phase Margin: Longitudinal All-loop ..... . . 65

6.13 Gain and Phase Margin: Lateral Outer-loop.. .. ....... .. 66

6.14 Gain and Phase Margin: Lateral Inner-loop . . . . 66

6.15 Gain and Phase Margin: Lateral All-loop . . . . 69

6.16 Eigenvalues for the Perturbed Closed-loop System: 20% Error in clfin 70

6.17 Pitch Command Tracking for Perturbed Linear System : q, u . . 71

6.18 Pitch Command Tracking for Perturbed Linear System : 8c, . 71

6.19 Pitch Command Tracking for Perturbed Linear System : 6ei, e . 71

6.20 Roll Command Tracking for Perturbed Linear System : p,r . . 72

6.21 Roll Command Tracking for Perturbed Linear System : 1, . 72

7.1 Complete Nonlinear Simulation with LQR Controller . . .. 75

7.2 Pitch Command Tracking : p, q .. . . . . . 76

7.3 Pitch Command Tracking: . . . . . 76

7.4 Pitch Command Tracking: e,,pe . . . . . 76

7.5 Pitch Command Tracking : 81, {x,y,z} Trajectory . . 77

7.6 Roll Command Tracking: p, q .. . . . . . 77

7.7 Roll Command Tracking: . . . . . 78

7.8 Roll Command Tracking: 8ei,e . . . . . 78









7.9 Roll Command Tracking: ,1,,1 . . . . . 78

7.10 Roll Command Tracking: {x,y,z} Trajectory . . . . 79

7.11 Roll & Pitch Command Tracking: p, q . . . . . 79

7.12 Roll & Pitch Command Tracking: .. .. .. ..... 80

7.13 Roll & Pitch Command Tracking: 1, . . . . 80

7.14 Roll & Pitch Command Tracking: ,1, 1 . . . . 80

7.15 Roll & Pitch Command Tracking: {x,y,z} Tracking .. . ...... 81

7.16 Response for 20% Variation inclfi: u,w . . .. . 81

7.17 Response for 20% Variation in clfi: p, q . . .. . 82

7.18 Response for 20% Variation in clfi,: ,1 . . . . 82

7.19 Response for 20% Variation in clfi,: {x,y,z} Trajectory . . 82

A.1 Rotation of Frames.. . . ................... 84

B.1 Shape Function for One Dimensional Quadratic Scheme . . 87















Abstract of Thesis Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Master of Science

CONTROL STRATEGIES FOR SUPERCAVITATING VEHICLES

By

Anukul Goel

December 2002


Chair: Andrew J. Kurdila
Cochair: Richard C. Lind
Major Department: Mechanical and Aerospace Engineering


Underwater travel is greatly limited by the speed that can be attained by the vehicles.

Usually, the maximum speed achieved by underwater vehicles is about 40 m/s. Supercav-

itation can be viewed as a phenomenon that would help us to break the speed barrier in

underwater vehicles. The idea is to make the vehicle surrounded by water vapor while it

is traveling underwater. Thus, the vehicle actually travels in air and has very small skin

friction drag. This allows it to attain high speeds underwater. Supercavitation is a phe-

nomenon which is observed in water. As the relative velocity of water with respect to

the vehicle increases, the pressure decreases and subsequently it evaporates to form water

vapor. Supercavitation has its drawbacks. It is really hard to control and maneuver a su-

percavitating vehicle. The first part of this work deals with modeling of a supercavitating

torpedo. Nonlinear equations of motion are derived in detail. The latter part of the work

deals with finding a controller to maneuver the torpedo. A controller is obtained by LQR

synthesis. For the synthesis, it is assumed that the cavity is fixed and the torpedo is situated

symmetrically in the cavity. The robustness analysis of the LQR controllers is carried out









by calculating the gain and phase margins. Simulations of the response for a perturbed

system also have been studied.















CHAPTER 1
INTRODUCTION


Achieving high speeds is an important issue for underwater vehicles. Even the common

fastest underwater vehicles are restricted to travel at speeds around 40 ms-1. 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 ms-1. This phenomenal improvement in speed was made

possible by supercavitation. The idea was to envelop the torpedo in a gas so that it has little

contact with water. The Shkval was able to achieve a tremendous reduction in skin friction

drag and exhibit high speed.

1.1 Cavitation

As the speed of an underwater vehicles increases, i.e., as the relative velocity of water

with respect to the vehicle increases, the pressure decreases. The speed can be increased

to the point the pressure goes below the vapor pressure of water and then bubbles of water

vapor are formed. This process is known as cavitation. Cavitation is mostly observed at

sharp covers of a body where the speed can reach high magnitudes. A classic example for

cavitation is at the tip of propellers, like the one shown in Figure 1.1. Since the propeller






2

is rotating at high speed, the relative velocity at the tips is high enough so that water at

the tips vaporizes. Cavitation has been extensively researched, but remains a challenge for

underwater vehicles.














Figure 1.1 Tip Vortex Cavitation [2]

Cavitation can be harmful in many cases. The cavitation region is usually turbulent.

When the cavitation is not steady, the pressure drop is momentary and very quickly the

cavitation bubble encounters a region of high pressure that forces the bubble to collapse

like a tiny bomb. This collapse causes high levels of noise and also may cause considerable

damage to the surface of the body.

Figure 1.2 shows the various stages of formation of cavity. It shows a bullet-like body

traveling in water at high speed. The cavitation starts as vapor bubbles near the region of

high relative velocity. As the speed is further increased, the bubbles merge to form a large

area of vapor. On further increase in speed, the whole of the body is covered in vapor. This

stage is called the supercavitation. At this point, the condition is similar to traveling in air.

The skin friction drag is tremendously reduced, and thus high speed can be attained in this

phase.

1.2 Types of Supercavitating Projectiles

The first version of the Russian Shkval was a crude supercavitating vehicle. It was

unguided and had a small range of about 5 miles. Now, there are various advanced forms

of supercavitating bodies. One class of supercavitating bodies, called RAMICS (Rapid













PFigure 1.2 Formall tion of Cavity [3]










Airborne Mine Clearance System), is a helicopter-borne weapon that destroys surface and

near-surface marine mines by firing supercavitating rounds at them. These are small bullet-

like, flat nosed projectiles that are designed to travel stably through air and water. As the

RAMICS can produce a cavitation bubble, they can travel at high speed underwater and

can pierce the mines to destroy them. As they are fired from a distance in air, they need to

be designed to travel in both phases. The RAMICS are better than conventional bullets, as
conventional bullets are rapidly slowed down by drag in water.













Another kind of a supercavitating projectile is a sub-surface gun system using Adapt-

able High-Speed Undersea Munitions (AHSUM). These are supercavitating "kinetic-kill"

bullets, fired from guns fitted under submerged hull of submarines. These sonar directed

bullets would be used to intercept undersea missiles.

The RAMICS and AHSUM are uncontrolled small range supercavitating projectiles.

The next level of supercavitating projectiles is larger torpedoes, with higher speeds and

longer range. These vehicles are much more complex. They require the design of a launch

station. They require detailed studies of hydrodynamics, acoustics, guidance and control,

propulsion, etc. Even if the vehicle is designed to be an uncontrolled torpedo, it is still

a challenge to produce and maintain a cavity around the vehicle. The cavity is usually

produced by the nose of the vehicle, which is specially shaped for the purpose. The nose

is called a cavitator. The cavitator may not be sufficient to produce the cavity. Thus air can

be blown at the nose and various sections of the body so as to maintain and produce the
IJmillcd t.-aUl


--- r^Try-



Figure 1.2 Formation of Cavity [3]


Airborne Mine Clearance System), is a helicopter-bome weapon that destroys surface and

near-surface marine mines by firing supercavitating rounds at them. These are small bullet-

like, flat nosed projectiles that are designed to travel stably through air and water. As the

RAMICS can produce a cavitation bubble, they can travel at high speed underwater and

can pierce the mines to destroy them. As they are fired from a distance in air, they need to

be designed to travel in both phases. The RAMICS are better than conventional bullets, as

conventional bullets are rapidly slowed down by drag in water.

Another kind of a supercavitating projectile is a sub-surface gun system using Adapt-

able High-Speed Undersea Munitions (AHSUM). These are supercavitating "kinetic-kill"

bullets, fired from guns fitted under submerged hull of submarines. These sonar directed

bullets would be used to intercept undersea missiles.

The RAMICS and AHSUM are uncontrolled small range supercavitating projectiles.

The next level of supercavitating projectiles is larger torpedoes, with higher speeds and

longer range. These vehicles are much more complex. They require the design of a launch

station. They require detailed studies of hydrodynamics, acoustics, guidance and control,

propulsion, etc. Even if the vehicle is designed to be an uncontrolled torpedo, it is still

a challenge to produce and maintain a cavity around the vehicle. The cavity is usually

produced by the nose of the vehicle, which is specially shaped for the purpose. The nose

is called a cavitator. The cavitator may not be sufficient to produce the cavity. Thus air can

be blown at the nose and various sections of the body so as to maintain and produce the









cavity. Figure 1.3 shows a supercavitating torpedo traveling underwater. It can be seen that

the whole of its body is enveloped by a cavity. This is the kind of vehicle that has been

studied in this work.







Figure 1.3 Supercavitating Vehicle [1]

1.3 Related Research

Research in the field of supercavitation has been going on from the early 1900's. But

earlier research was aimed at reduction of cavitation so as to reduce noise and body damage.

In the 90's the focus shifted to exploiting the effects of supercavitation.

The work shown in May [4] has an extensive collection of experimental data for vari-

ation of forces on various supercavitating shapes. Coefficients of lift and drag are plotted

with the variation of cavitation number for shapes like disks, cones, ogives and wedges.

The work done in this research makes use of a CFD database provided in Fine [5]. This

database has values for coefficients of lift and drag for conical cavitators, which are func-

tions of the half angle of the cone and the angle of attack. This database also has coefficients

of lift and drag for wedges, which are functions of wetted surface of the wedge, angle of

attack and sweepback angle (we discuss the definition of these geometric quantities such as

the angle of attack and half angle later in this thesis). This information is useful to calculate

the forces on fins of the torpedo.

In late 90's a lot of research has been done on the dynamics of supercavitating vehicles.

Work shown in Kulkami and Pratap [6] and Rand et al. [7] deals with studying dynamics

of uncontrolled supercavitating projectiles. A dynamic model for RAMICS and AHSUM

has been developed. It is shown that these projectiles rotate inside the cavity. This rotation

leads to impacts between the tail of the projectile and the cavity wall. The frequency of the






5

impact increases with time. These projectiles are very short range and have a small time of

flight, on the order of a few seconds.

The work shown in Dzielski and Kurdila [8] focuses on the formulation of a control

problem for a supercavitating torpedo. A dynamical model for a fin controlled torpedo

has likewise been developed. The model also includes a formulation for the cavity. It is

observed that the weight of the body forces it to skip inside the walls of the cavity and the

vehicle is unstable. A control system is designed for the torpedo and results of closed-loop

simulations have been presented.

1.4 Overview of this Thesis

This work aims at formulating the control design for a supercavitating torpedo. Equa-

tions of motion for the torpedo are derived, and linear control methodologies are applied

for pitch and roll control of the torpedo. The main purpose of this work is to analyze these

controllers and ascertain their robustness to various errors and uncertainties.

Chapter 2 of this thesis briefly describes the configuration of the supercavitating torpedo

for which the equations of motion have been developed.

A detailed derivation of the equations of motion for the torpedo has been carried out in

Chapter 3. Various reference frames have been used to obtain the kinematic equations of the

torpedo. Dynamic equations are derived using Newton's Laws. Various forces experienced

by different regions of the torpedo have been explained.

Chapter 4 describes linearization of the equations of motion using small disturbance

theory. Numerical methods are used for this purpose. It is observed that the linearization

for a simple trim can be very complicated.

Chapter 5 describes the control synthesis for the torpedo. Open-loop dynamics are

shown. The closed-loop 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 closed-loop system and a perturbed liner system have been shown.






6

The results for pitch and roll rate control for the nonlinear closed-loop system have

been presented in Chapter 7. The simulations for a perturbed system model also have been

presented.















CHAPTER 2
CONFIGURATION OF VEHICLE


Although supercavitation can be a very helpful phenomenon, it presents significant

challenges in modeling and control of supercavitating vehicles. As a significant portion of

the vehicle is located in the cavity, the control, guidance and stability of the vehicle has to

be managed by very small regions in front and aft of the vehicle. Also generation of the

cavity can be a significant problem. The major problems associated with the supercavitat-

ing vehicles may be summarized as:

* generation and maintenance of cavity
* balancing the weight of the vehicle
* control and guidance
* stability

Figure 2.1 is a figure of a supercavitating torpedo that is modeled in this work. The main

parts of the torpedo are the cavitator in the front and the four fins in the aft portion. The cav-

itator is used to generate and maintain the cavity. The Cavitator and the four fins together

are also used for control and stability of the vehicle.

2.1 Cavitator

The cavitator is the element that generates a cavity around the torpedo. Several types

of cavitators, including cones, wedges and plates have been investigated [4]. This project

will consider a conical cavitator as shown in Figure 2.1. The main parameter that defines





Figure 2.1 Supercavitating Vehicle [8]

Figure 2.1 Supercavitating Vehicle [8]









a conical cavitator is its half-angle. The coefficients of lift and drag, as functions of half-

angle, for the cavitator have been generated using CFD and tabulated in [5].

The cavitator in this model has one degree of freedom defined by a rotation angle about

an axis perpendicular to its axis of symmetry. The shape and location of the cavity is

a nonlinear function of this cavitator deflection angle and half angle of the cone. This

shape determines the position where the cavity intersects the body of the vehicle. Thus, the

amount of wetted area of the vehicle is dependent on these angles, which in turn determines

the efficiency of supercavitation achieved.

As stated earlier, a large portion of the vehicle is in the cavity. The lift produced by the

cavitator combined with the lift produced by the fins helps in balancing the weight of the

body.

2.2 Fins

Although the cavitator is capable of providing enough lift to sustain the body in water,

it does not usually provide sufficient forces and moments to stabilize and control the ve-

hicle. For this purpose fins are required in the aft portion of the vehicle. The fins help in

counteracting the moments produced by the cavitator and also providing some amount of

lift to support the weight of the body. There are four fins placed symmetrically along the

girth of the vehicle, near the tail. The fins also are the control surfaces, as they can provide

differential forces, thus causing control moments. Two of the fins shown in Figure 2.2 are

horizontal (placed parallel to the axis of rotation of cavitator), called elevators, are used to

affect the longitudinal dynamics for the vehicle. The other two fins are called the rudders

and are used to affect the lateral dynamics to the vehicle.

The fins have two degrees of freedom, a sweepback angle and an angle of rotation.

These angles will be explained further in Chapter 3.

2.3 Maneuvering

The motion of the vehicle is controlled by all five control surfaces, the four fins and the

cavitator. In a straight line motion the cavitator and the elevators balance the weight of the















Figure 2.2 Cavitator and Fins


vehicle. A propulsion force at the tail pushes the vehicle forward. The rudders usually have

a zero deflection in such a case.

The vehicle depends on a bank-to-turn strategy for making a turn, and cannot use tradi-

tional missile strategies such as skid-to-turn. This dependence arises because the hull of the

vehicle is incapable of providing a lift force. The fins are the main lift generating surfaces.

A differential force from the fins can be used to generate a force towards the center of the

turn.















CHAPTER 3
NONLINEAR EQUATIONS OF MOTION


The dynamics of the torpedo, whose configuration was discussed in Chapter 2, can be

mathematically formulated. A complete derivation of the equations of motion is given in

this chapter. Section 3.1 deals with the setup of reference frames and derivation of the

kinematic equations. The dynamic equations of motion are derived in Section 3.2.

3.1 Kinematic Equations of Motion

The definition of a suitable coordinate system and degrees of freedom is a precursor

to any formulation of equations of motion. The derivation of the equations of motion of

multi-body systems is highly simplified by defining various reference frames and relations

between them. Appendix A describes briefly the concept of reference frames and rotation

matrices. These concepts will be used extensively in the derivation of equations of motion.

The derivation of the equations of motion will be done in two parts. First, the kinematic

equations will be derived. These include the formulation of the position vectors, velocities

and accelerations of various points on the torpedo. Next, the dynamic equations will be

derived. Finally, Newton's Laws yield the dynamic equations of motion.

3.1.1 Orientation of the Torpedo

Six degrees of freedom (DOF) are required to describe the position and orientation of

the torpedo when it is considered a rigid body. Of these, three DOFs give the location of a

point fixed on the torpedo. The other 3 DOFs give the orientation of the torpedo in a fixed

reference frame. The position of the torpedo in a reference frame can also be obtained by

the integration of its velocity in that reference frame.

The torpedo is assumed to be moving in an earth-fixed reference frame E, centered at

any conveniently chosen point and described by the basis vectors (el e^, 3). The earth-fixed






















3 e


Figure 3.1 Body-fixed and Inertial Frames

reference frame is shown in Figure 3.1. The vector 63 points in the downward direction,

i.e., the direction of the gravity. The vectors e\ and j2 are placed in the horizontal plane

such that the basis vectors form a right-handed coordinate system. As shown in the figure,

ei points to the right and ^2 points outside the plane of the paper. A body-fixed frame B is

defined by the basis vectors (bi, b2, b3) so as to simplify the derivation of the equations of

motion. The frame B is centered at 0, the center of gravity of the torpedo, and moves with

the torpedo. The reference frame B is shown in the Figure 3.1. It can be seen in Figure

3.2 that the torpedo has two planes of symmetry namely b1b2 and b1b3. The plane b1b3 is

called the longitudinal plane and plane b1b2, the lateral plane.

Transformation from E to B can be achieved by 3 rotations. Many such sets of rotations

are possible. The rotation steps chosen here are the Euler angles of rotation, which are the

yaw, pitch and roll. As there are three rotations to be performed, two intermediate reference

frames with basis vectors (21,42,43) and (C1i,2,j3) will have to be introduced to perform

the transformation. The rotations, to be performed in order, are listed below.

1. Rotate the frame E about 63 through a yaw angle, ?, to obtain the frame (21,42,23).
2. Rotate (h1,22,h3) about x2 through a pitch angle, 0, to obtain the frame (\1,Y2,Y3).
3. Rotate (i1,J2,j3) about Uf through a roll angle, (, to obtain the frame B.

























Figure 3.2 Principle Planes of Symmetry for the Torpedo


3.3
Figure 3.3 Euler Angles of Rotation


Figure 3.3 shows the above rotations in order. Based on the above definition of rotations,

the transformation matrix from E to B can be written as in equation 3.1.


1 0 0 CO 0 -so CY SwT

0 C' SO 0 1 0 -ST CY

0 -SO Cc SO 0 CO 0 0

COCY COST -SO

CYTSSO COST SISOST + CYCc SICO

CISOCY + SO ST CSSOST CS C(CCO


0 fe

0 j2

1 83

j ]
e2


e3


bi

b2

b3


(3.1)










3.1.2 Orientation of the Cavitator

As described earlier, the cavitator has only one degree of freedom. It can rotate in the

longitudinal plane about an axis parallel to the b2 axis. The orientation of the cavitator-fixed

axes with respect to the body fixed axes is shown in Figure 3.4.



A.) bi


Ap e
ACP
b3 ^ \








A
b3
Figure 3.4 Cavitator Reference Frame

The deflection of the cavitator is given by an angle, 8c, which is positive for a positive

cavitator rotation about the b2 direction. The geometric center of rotation of cavitator is

denoted by P. CP is the center of pressure for the cavitator, which is at a distance Acp from

P, along c1.

From Figure 3.4, the rotation matrix from B to cavitator fixed frame C, can be written

as in Equation 3.2.
ei C8c 0 -SSc bh

c2 = 0 1 0 b2 (3.2)

C3 SSc 0 CSc h3
k ^ k >










3.1.3 Orientation of Fins

Figure 3.5 shows the orientation of the fin-fixed reference frames. For convenience, all

the fin frames are indicated by basis vectors (11,12,~3). In text they will be represented as

(f1,f2,f3) fn, where subscript fin corresponds to a particular fin.


Rudder (1-


-fRu



Rudder 22


FRONT VIEW
f2


TOP VIEW


Elevator l






E-----------At



Elevator 2


Figure 3.5 Rudder and Fin Reference Frames


All the fins have 2 DOFs, namely the sweepback angle (Ofi,) and the fin rotation (fi,).

These can be defined as

Sweepback angle (Ofi,): This parameter is the angle of rotation of a fin about its f3
axis. The direction of rotation for positive sweepback is such that the fin is moved
toward the rear portion of the torpedo. Due to this definition, the positive sweepback
is along the negative 3 direction for all fins. Sweepback angle determines the amount
of fin that is enveloped in the cavity.


a, f








* Fin Rotation (8,i,): This parameter is the angle of rotation of the fin about its f2 axis,
in positive the f2 direction. Fin rotation determines the magnitude and direction of
the forces acting on the fins.
The order of rotation in the above case is important to obtain the correct equations.
Sweepback has to be performed before fin rotation. An intermediate reference frame G
with basis vectors (gI, 2,g 3) is introduced so as to simplify the derivation of rotation ma-
trix from B to the fin coordinates. Sweepback aligns the fin-fixed frames with the interme-
diate frame G and then a fin rotation puts the fin frame in actual orientation with the fins.
As the second rotation is identical in all cases, the transformation matrix from frame G to
fin frame Ffi, can be written as in Equation 3.3.

l C8fin 0 -S8fin gl
f2 = 0 1 0 I 2 (3.3)
f3 fin S8fn 0 Cafn1 g3 fin

The rotation matrix for sweepback and the transformation matrices from B to Ffin frame
for each of the fins can be derived easily, and are summarized below.
* Rudder 1 Figure 3.6 shows the sweepback and fin rotation for rudder 1. The matrices
for transformation from B to R1 can be derived as in Equation 3.4 and Equation 3.5.

Si [ -COR1 0 SOR1 I1
S2 = -SORI 0 -COR i 2 (3.4)
g3 R1 0 -1 0 J3
fC-R1 0 -SR1 -C 0 SOR f
i 2 0 1 0 -SOR1 0 -COR1 2 (3.5)
S3 SSR1 0 CR1 0 -1 0 3

Rudder 2 Figure 3.7 shows the sweepback and fin rotation for rudder 2. The matrices
for transformation from B to R2 can be derived as in Equation 3.6 and Equation 3.7.

g -COR2 0 -SOR2 ~j
S2 = -SOR2 0 COR2 2 (3.6)
g3 1 0 b3
SCs8 0 -S8R2 -COR2 0 -S
2 0 0 1 0 -SSR2 0 COR2 2 (3.7)
1i3 R SSR2 0 CR2 0 -1 0 ( 3












OR1


gz,f2

6R1






f
4


g3 f3


Figure 3.6 Rudder 1 Fin Reference Frames


Figure 3.7 Rudder 2 Fin Reference Frames

* Elevator 1 Figure 3.8 shows the sweepback and fin rotation for Elevator 1. The
matrices for transformation from B to El can be derived as in Equation 3.8 and
Equation 3.9.


0 b\
0 b2


-COEI -SOE1
-SOE1 COE1
0 0


0 b2
-1 3 3


* Elevator 2 Figure 3.9 shows the sweepback and fin rotation for Elevator 2. The
matrices for transformation from B to E2 can be derived as in Equation 3.10 and


61 3
g2

g3 El


-COE1
-SOE1
0


CSEl
0
S8EI


-SOE1
COEI
0

-S8E1
0
C E1


(3.8)


(3.9)









g2,f2







fl
\


g3 f3


Figure 3.8 Elevator 1 Fin Reference Frames

b,2

E2 3 31
E
>E2





0E2




Figure 3.9 Elevator 2 Fin Reference Frames
; ________^^ yI








Figure 3.9 Elevator 2 Fin Reference Frames


Equation 3.11.


SOE2 0 r
-COE2 0 b2
0 1 3
-SSE2 -COE2
0 -SOE2
CSE2 0


SOE2
-COE2
0


10 b2
1 I 3


Equations 3.2 to 3.11 will be used in later sections to transform the forces on fins and


cavitator to the body-fixed frame.


g2
g83 E2

f E2
I j3 P 2


S-COE2
-SOE2
0
C8E2 0
0 1
S8E2 O


(3.10)


(3.11)









3.1.4 Angle of Attack and Sideslip

The body-fixed reference frame has been defined, but the velocity of various points on

the body of the torpedo is yet to be defined. The torpedo is considered as a rigid body. If the

velocity of a certain point on a rigid body is known, the velocity at any other point on that

body can be found by knowing the rotation matrices. Thus, V = ub1 + vb2 + wb3 will be

taken as the velocity of CG of the torpedo. The velocity of other points on the torpedo can

be defined subsequently. Two very useful parameters, angle of attack and angle of sideslip

can be defined in conjunction with the orientation of the body axis with the velocity of CG.

Figure 3.10 shows these parameters and their geometric interpretation.




V ubi +( v)(b ,)+wb,
-V





b3
Figure 3.10 Angle of Attack (a) and Sideslip (3)

Angle of attack, a, can be defined as the angle between the projection of velocity V

onto b2b3 plane and the b1 axis. Angle of attack is positive when the nose of the torpedo

points above the velocity vector of the torpedo. As before, an intermediate frame F given

by (, f12,13) can be described by rotation of the B frame by an angle a. Thus the rotation
matrix from Fbody to B can be written.

by Ca 0 -Sa fi

b2 = 0 1 0 f2 (3.12)

b3 Sa 0 Ca f3
J body

The Angle of sideslip, 3, is defined as the angle between the actual velocity V and the

projection of V onto b2b3 plane. Again, a frame Gbody can be defined by rotation of Fbody









by an angle equal to 3 in negative f3 direction, thus giving a rotation matrix as in Equation

3.13.


Sg2 S= S C3 0 /2 (3.13)

3 body 1 3 bo

Velocity of CG of the torpedo in the Gbody frame can be written as Vgl, where V is

magnitude of V. It will be seen that drag and lift on the torpedo can be obtained in this

frame. Thus a transformation from Gbody to B is important. It is given by Equation 3.14.

b1 CPCaC C cSp -S( ]i

b2 = -S3 C3 0 g2 (3.14)

b3 C3Sc SS3 Ca J body3

Using Equation 3.14 V can be rewritten as in Equation 3.15.

V = Vgl
(3.15)
= VCp3Ca b VS b2 VCS b3(3.15)
H V W

where V2 = V2 = 2 +2 +w2. From Figure 3.10, relations between the velocity compo-

nents and the angles of attack and sideslip can be derived.

w
tan = (3.16)
u
-v
sin 3 = (3.17)

Though the matrix Gbody-B in Equation 3.14 has been defined for the body-fixed ref-

erence frame and velocity of CG of the torpedo, the equation is valid for any other rigid

part of the body like the fins and cavitator. Thus, in case of a fin, the velocity V would

correspond to a point (like the tip, center of pressure etc.) on that fin, and Gfin_B matrix

would correspond to the fin-fixed reference frame. In this case the velocity of center of

pressure of the fin will be used to define the above parameters. Thus, obtaining afi, and

3fin is a two step process:








1. Obtain the velocity of center of pressure of fin.

VCPbody = Vcg +E oB X rcgCP (3.18)

where Vcpbody is velocity of CP of fin in B frame, Vcg is the velocity of CG of the
torpedo in E frame, E o is angular velocity of B in E, and rcgcp is position vector
from CG to CPf,. Equation 3.18 can be rewritten as in Equation 3.19.

{ fin bl b2 b3
Vfpin, v + p q r (3.19)
Wfin B w g Xfin Yin Zin

where rcgCp = Xfinbl + Yinb + Zfinb3.

2. Transform the velocity (in E) of CP of fin from frame B to frame of corresponding
fin. This transformation is obtained by using rotation matrices derived in Equations
3.3 to 3.11.


UR1
VR1
WR1
UR2
VR2
il12


}:


R2


C6R1
0
S6R1
C8R2
0
S8R2


-SSR1
0

-S8R2
0
C6R2


UEI CsE1 0 -SE1I
VE 1 0 1 0
WE1 l E SE1 0 C8E1
UE2 C6E2 0 -S8E2
VE2 0 1 0
1, E2 S8E2 0 C8E2

The left hand terms in Equations 3.20


-COR1
-SOR1
0
-COR2
-SOR2
0
-COEI
-SOEI
0
-COE2
-SOE2
0


0 SORI UR1
0 -CORI VR1
-1 0 WR1 B
0 -SOR2 U R2
0 COR2 VR2
--1 0 12 RB
-SE1 o0 UEI
COE1 0 VE1
0 -l1 WE1 B
SOE2 0 UE2
-COE2 0 VE2
0 1 lI B


to 3.23 give the velocity components at the CP


for corresponding fins, in that fin frame. These can be used in Equations 3.16 and 3.17 to

find the angle of attack and sideslip for a particular fin.

3.1.5 Kinematic Equations

Velocity of the CG of the torpedo has been defined in the previous section. That velocity

defines the translational kinematics for the torpedo. Analogous to this quantity, angular

velocity is required to define the rotational kinematics. The angular velocity of the hull has

components p, q and r in the frame B.


EoB Apb + qb2 + rB3


(3.24)


(3.20)


(3.21)


(3.22)


(3.23)








As the transformation from E to B has already been defined in terms of rotational motions

give by Euler angles, the angular rates can also be obtained by differentiation of Euler

angles. Thus, another form of Equation 3.24 can be written as in Equation 3.25.

E-B = Ye3 + Ox2 + [1 (3.25)

The rotation matrices from Equation 3.1 can be substituted into Equation 3.25 to obtain

Equation 3.26.

E oB = (6 SO )bi + (YC SO + OCD)&b2 + (COCO OS) b3 (3.26)

Equations 3.24 and 3.26 can be equated to obtain Equation 3.27.


p -So 0 1 Y
q = CSO COS 0 (3.27)

r COCO -SO 0 6

Equation 3.27 can be rewritten as in Equation 3.28.

p = -SOY

q = TCOSS + OCO (3.28)

r = TCOCO OSO

To apply Newton's Laws, acceleration of the CG is required. The values of p, q, r

will be the angular accelerations of torpedo in B. The translational acceleration can be

calculated by time differentiation of V in Newtonian frame. A differentiation formula can

be used to find the time derivative, in some frame, for a vector defined in some other related

frame.
S(v) = (v) +I o xv (3.29)
dt I dt B

where, subscript I denotes Newtonian (inertial) frame, and B is the body-fixed frame. I'B

is angular velocity of the body (or body-fixed frame) in the I frame, v is the velocity in

I frame, of the point where acceleration is desired. Using the formula the acceleration of









CM of torpedo in E can be obtained.

f b Z2 3
EaCM + p q r (3.30)

w u v w

u +qw vr

= + ur pw (3.31)

w + pv- uq

Similarly, the rotational acceleration will be required in the frame E. This can be written

analogous to Equation 3.30.

b, bi2 b3
EOI= + p q r

r p q r
(3.32)
P
= q

r

The position of torpedo can be found by integrating the velocity. Let (x,y, z) represent

the coordinates of CG in frame E. Thus, the time derivative of these coordinates in E should

represent the velocity components of the torpedo in E frame. When these time derivatives

are transformed to body-fixed frame, they would be equivalent to the velocity components

in body-fixed frame.
x u

y = v (3.33)

Ez w

Equation 3.1 can be substituted in Equation 3.33 to obtain Equation 3.34.









x COCY COsT -so u



z C)SOCY +SST CS OSY- SOCY CCOC w

3.2 Dynamic Equations of Motion

Now that the accelerations of various parts of the torpedo are known, Newton's Laws

can be used to derive the dynamic equations of motion. Newton's laws state that the rate of

change of momentum is equal to the sum of external force applied on the body. Equation

3.35 can be obtained from Newton's laws by an assumption that the mass of the torpedo is

constant. This assumption is valid for a short period of time. The equations are

YF=P
(3.35)
= ma

where P is the linear momentum, m is mass of the body, a is the acceleration and F is all the

forces of the body. Using Equation 3.31 in Equation 3.35, Newton's Laws for the torpedo

can be rewritten as in Equation 3.36.


ui + qw -vr
m v+ur -pw =F (3.36)

w + pv uq


Newton's laws can be extended to rotation. Equation 3.37 are the Newton's Laws for

rotational motion.


(3.37)
=Ic + E (B x H

where H (= IcMEoB) is the angular momentum, IcM is moment of inertia matrix of the

body, a is the angular acceleration and M is all the moments on the body. It should be

noted that above stated Newton's laws are only valid when the quantities are calculated in

an inertial frame of reference. Thus, the quantities have been calculated in frame E. Using

Equation 3.32, the Newton's Law for rotation can be written as in Equation 3.38.









1I 0 0 1p bI b2 b3

0 12 0 q + p q r =AM (3.38)

0 0 13 r Ilp I2q I3r

To derive the equations, the forces on each of the parts will be calculated individually,

and then expressed in body-fixed reference frame, i.e., summation will be done in body

reference frame. The rotation matrices derived in previous sections will be used extensively

for this purpose.

Various types of forces are experienced by a moving torpedo in water. These forces can

be summarized as follows:

* Hydrodynamic Forces: These are the forces exerted by the fluid on the torpedo.
Thus they exist whenever the fluid is in contact with body. For supercavitating ve-
hicle, most of the body (hull) is inside the cavity. Only a portion of the fins and the
cavitator are in contact with the fluid. Thus the lift and drag on cavitator and fins are
only hydrodynamic forces. The coefficients of lift and drag for the fins and cavitator
are functions of their angle of attack, immersion, sweepback angle, angle of rotation
etc. and are tabulated in a database [5]. This database will be interpolated and ex-
trapolated to calculate the numerical values for the forces. Occasionally, a part of the
hull comes in contact with water and might experience some forces. These forces are
known as planing forces. It is assumed that the vehicle is centered in the cavity. Thus
the planing forces are neglected in the summation of the hydrodynamics forces.

FHydrodynamic = FR1 + FR2 + FE1 + FE2 + Fc (3.39)
MHydrodynamic = MR1 +MR2 +ME1 +ME2 +Mc (3.40)

Gravitational Forces: This is the gravity forces on the body. As the summation of
moments will be taken about the center of gravity, the gravitational forces will not
contribute to the summation on moments. They will appear only in summation of
forces.

Propulsive: The torpedo has a propulsion system, which is similar to that of rockets.
The line of action of the propulsive force is assumed to be passing through center of
gravity and along b1 axis. Thus this force will also contribute to the sum of forces,
and not moments. The propulsive forces are provided by burning the fuel, but for
simplicity it will be assumed that the mass of the torpedo remains constant.

The total forces and moments are expressed in terms of these components.

F = FHydrodynamic + FGrav + FProp (3.41)


M = MHydrodynamic


(3.42)









3.2.1 Forces on Cavitator

Figure 3.11 shows the forces acting on the cavitator. Coefficient of lift (clc) and drag

(cdc) for the cavitator are functions of angle of attack, ac, and half-angle, yi, of the cavi-
2
tator. Half-angle, for a cone, is the angle made by axis of the cone with any line passing

through the vertex and lying in the surface of the the cone. This parameter defines the main

geometry of the conical cavitator. The values of cl/ and cdc are determined using CFD and

tabulated in [5]. These values have been extrapolated to calculate lift and drag.


Lc = pV Sccl (c, ) (3.43)
2 2
Dc = 2pV2Secdc(oc,7y) (3.44)

where Sc is the cross-sectional area of the cavitator base. Directions of lift (Lc) and drag

(Dc) are as shown in Figure 3.11(b). These can be transformed to the body axes using 3.2

and 3.14 for the cavitator.




b2 = CB(c) x Gcav-C(ac, 3c) x &2 (3.45)




SC^2 M- n


Lc Dc
g3

(a) (b)

Figure 3.11 Cavitator: (a) Angle of Attack and Sideslip and (b) Hydrodynamic Forces









Thus the forces on cavitator, in body frame, can be written.


F JR

Fc = Fcy
Fc, B


-D c(ac, Yi)
= 0 (3.46)

-Lc(ac,,7 )
20 cav

C8c 0 SSc Cp cCac CacXSc -Sac -Dc (ac 7)

= 0 1 0 -Sc CPc 0 0

-sac 0 C6c CcSac ScSc Cuc -Lc((c, 7)

where Fc is a 3x1 matrix with each row corresponding to each direction in B basis. The

forces are assumed to be acting at CP of the cavitator. Once the forces have been calculated,

the moment about any point can be calculated.

Mc = rcpcav X Fc (3.47)

where rcPcav is the position vector from that point to CP of cavitator. It is assumed that the

CP lies on b1 when cavitator deflection is 0, and its coordinates with respect to body origin

O, in this case, are (Xc, 0, 0). Thus from Figure 3.4, the coordinates of CP can be written

for the case when the cavitator is deflected.


Xc +ACPC8c

rcpcav = 0 (3.48)

-AcPSc body
body

The moments on the cavitator in body-fixed can be obtained by substituting Equations 3.46

and 3.48 in Equation 3.47.








C c 0 SSc

M = [(Xc+ AcpCc) bi- AcpS8c3] x 0 1 0

-SSc 0 CSc
(3.49)
CC(,ca, Ccas -Sa, -D(ac, i)

-SIp cpc 0 0
CPcSoc( SacSPIc Cac -Lc (oc,7i)

3.2.2 Forces on Fins

Fin forces are also extrapolated from the CFD database [5], which gives the values of

coefficients of lift (clfin) and drag (cdfin) for fins. These values are functions of angle

of attack Cafi,, fin sweepback Ofi,, fin rotation 8fi,, fin immersion Ifi, and fin geometry.

Figure 3.12 shows these parameters graphically, and they can be defined as follows:

* Fin Geometry: The geometry parameters for fins are L and S, and wedge half angle
(ha), as shown in Figure 3.12. These parameters are fixed according to the values
given by the database, so as to calculate the forces accurately.

Fin Immersion: As the fin is partially wetted by fluid, the wetted length can be rep-
resented by parameter So along fin Y-axis. The immersion Ifi, can be defined as the
ratio of the wetted length of the fin to its true length.
fin = (So/S)fin (3.50)
Ifin determines the total hydrodynamic force on the fin.
Fin Rotation (8fi,): As defined earlier, this is rotation about fin j2 axis. This deter-
mines the direction of the hydrodynamic force. Thus fin rotation is used as primary
control for the torpedo.
Fin Sweepback (Ofi,,): As defined earlier, this is rotation about fin f3 axis. Sweepback
determines the wetted surface of the fin, thus the hydrodynamic force on the fin.

Angle of Attack: Angle of attack for the fin is calculated as described in Figure 3.12
and section 3.1.4.

The database gives clfin and cdfin as a function of (fi,, Ofin and Ifin, thus lift and drag on

the fins can be calculated by the normalizing factor.

1 2 2
Lfin = 2 V SinCfin (3.51)

Dfin = p2Sfincdfin (3.52)








tions 3.3 to 3.11.










Figure 3.12 Fin Geometry

Where Sfi is the length of the fin as shown in Figure 3.12 These forces have directions
as shown in Figure 3.12. Before substituting in Equation 3.39, these forces have to be
transformed to body-fixed reference frame. This process involves following two rotations:

1. Rotate the frame Ffin (which has Lfi, and Dfin along its basis vectors) to align it with
the fin-fixed frame using Equation 3.14 and
2. Rotate the above obtained fin-fixed frame to obtain the body-fixed frame using Equa-
tions 3.3 to 3.11.
Thus the forces on the fins in body-fixed frame axis can be obtained.

* Rudder 1
FRl,x -COR1 -SOR1 0 CSR1 0 S8R1
FRIy 0 0 -1 0 1 0
FRI,z B SOR1 -COR1 0 -SSR1 0 C6SR
S CPRICR1 CaRISPR1 -SaR1 -DR1I
_-SR1 CfR1 0 0
CPR1SaR1 SaR1SPR1 CaR1 -LR1

Rudder 2
FR2,x -COR2 -SOR2 0 F c8R2 0 S R2
{FR2,- Y 0 0 --1 0 1 0
FR2,z B -SOR2 COR2 0 -SSR2 0 C8R2
SCR2CaR2 Ca(R2SR2 -S(R2 -DR2
-SCR2 CSR22 0 0-
CIR2SaR2 SaR2SPR2 CUR2 -LR2 J








* Elevator 1
FEI,x -COEI -SOE1 0 CE1 0 S8Ei
FE1, = -SOE1 COE1 0 0 1 0
FEI,z B 0 0 -1 -SSEI 0 CSE
CE1CEcoE1 CoESIE1 -SE(1 -DE 1
-S3E1 CPEl 0 0
C3E1SaE1 SaE1SE1 CaE1 -LE1

Elevator 2
FE2,x1 -COE2 -SOE2 0 CE2 0 S8E2
FE2,y = SE2 -COE2 0 0 1 0
FE2z 0 0 1 -S8E2 0 CSE2
(3.56)
SCE2CaE2 C(E2SIE2 -S(E2 -DE2
-SPE2 CPE2 0 0
CPE2S(E2 S(E2SPE2 C(E2 -LE2

Equations 3.53 to 3.56 give the components of hydrodynamic forces on fins in body-fixed
frame. What remains is to find the moment of these forces about CG of body. The moments
can be obtained in analogous to Equation 3.47.

MIfin = rf_Cp X Ffi (3.57)

In this case, the root of fins is fixed to body, and it can move thus moving the CP of fin
relative to root. The position of CP from root is also know with respect to fin coordinates.

rC-root =Xotb I+Yrb2+Zr t3 (3.58)

oiot-cp = Aj1 +A 2 (3.59)

where (fi,2,/3) is fin-fixed coordinates for corresponding fin. Equations 3.58 and 3.59
can be added by using matrices given in Equations 3.3 to 3.11. Thus, the position vector
from CG to CP of fins can be obtained.








* Rudder 1
XR1 r1ot -COR1 -SOR1 0
YR1 + 0 0 -1
ZR1 Z ot R1 -COR1 0
SC8R1 0 SSR1 ~ Ay(3
0 1 0 Ayj
-S8R1 0 C8R1 0 R1

Rudder 2

SYR2 2 Yt 0 0 -1
ZR2 B Z -SR2 CR2 0
OL.2 (3.61)
C8R2 0 SSR2 (
0 1 0 Ayp
-SSR2 0 CsR2 0 R2

Elevator 1

YE I y( oot -S E1I COEI 0
ZE 1 B Z oEo B O 0 -
0 E (3.62)
C8E1 0 S8E A
0 1 0 AyCp
-SSE1 0 CSE1 0 El

Elevator 2
XE2 Xr 20t -C[ E2 -SOE2 0
YE2 Y t + SOE2 -COE2 0
ZE2 J R I ot I 0 1
B C8E2 0 8E2 (3.63)

C-8E2 0 CSE2 0 E

Equations 3.60 to 3.63 give the position vector from CG to CP of the fins. These equa-
tions in conjunction with Equations 3.53 to 3.56, used in 3.57, gives the external moments
on fins about the CG.
b1 b2 b3
Mfi = Xfi. Yfin Zfin (3.64)
Ffin,x Ffiny Ffin,z











3.2.3 Gravitational Forces

For simplicity, mass (m) of the torpedo is assumed to be constant over time. This

is not the case in reality but the approximation is reasonable for considering short time

maneuvers. Acceleration due to gravity, g, is also assumed to be constant as torpedo moves

in space. The direction of gravity is given by 63 axis. Thus, the gravitational force can be

written as in Equation 3.65.
Frav = mgd3 (3.65)

Equation 3.65 can be re-expressed in body frame of reference using Equation 3.1.


COCY COST -SO 0

Fgrav = CYS4SO COSY SISOST + COC SICO 0

CISOCY +S ST CISOST SICY CODCO mg
E (3.66)
-So
= mg SOCO

C(CO
B

3.2.4 Equations of Motion

Now that a mathematical formulation of forces on the torpedo has been achieved, these

equations can be substituted into Equations 3.36 to 3.42 to obtain the dynamic equations of

motion. Thus the force equations of motion can be summarized as in Equation 3.67.








3.2.4.1 Force Equations

i + c qw-vr -Fprop FRx
m r -+ u pWr = Fimmersion + 0 + FR y +

w+ pv uq 0 FR1,z

FR2,x FE1,x FE2,x -SO
FR2,y + FEly + FE2y +mg SOCO + (3.67)
FR2,z FE1,z FE2,z COI C

CSc 0 S81 CpcCac CacSpc -Sc1 -Dc(aey)

0 1 0 -SPc CPc 0 0

-Sc 0 C C CP ccSc S4cSPc Cac -Lc(ac,7y)

Some issues should be noted for Equation 3.67.

The planing forces have not been included in the equations of motion. These forces
are neglected by assumption that the vehicle is centered in the cavity.

The propulsion force is constrained to be along negative b\ axis.

3.2.4.2 Moment Equations

Il 0 0 p b b 2 b3 bI b2 b3

0 12 0 Y + p q r = XR YR1 ZR +

0 0 13 Ilp Iq I3r FRI,x FRIy FR1,z

b1 b2 b3 b1 b2 b3

XR2 YR2 ZR2 + XE1 YE1 ZE1 + (3.68)

FR2,x FR2, FR2,z FEI,x FEly FEI,z

b1 b2 b3 b1 b2 b3

XE2 YE2 ZE2 + Xc + ApC8c 0 -AcpSSc

FE2,x FE2y FE2,z Fc,x Fcy Fc,z

Some issues should be noted for Equation 3.68.









* Some of the terms in Equation 3.68 are shown as a determinant. They need to be
expanded and written as components in body-fixed frame so as to equate the left-
hand and right-hand terms.

Moments due to gravitation do not appear because the moments are taken about CG.

Again, the moments due to planing forces have been neglected.

3.2.4.3 Orientation Equations

j } 0 C- P



S1 SSO CO r

3.2.4.4 Position Equations

x COCY COsT -so u

y = CS~YSO-C SCY SISOSY +CYCI SCO v (3.70)

z CASOCY +SASY C(SDSY SOCY COCO w

Equations 3.67 to 3.70 are a complete set of equations of motions for the supercavitating

torpedo.















CHAPTER 4
LINEARIZATION


4.1 Linearization

4.1.1 Need for Linearization

The equations of motion for the torpedo are identical to airplane equations of motion

but the forces terms on the right-hand side of these equations are different. Thus, the

linearization can be carried out similarly, as shown in Nelson [9]. The equations of motion,

as in the case of a supercavitating torpedo, are represented by a set of first-order differential

equations.

x f(x, u) (4.1)

using f : 9)" 9" as a nonlinear function of a time-varying vector x E 9" and u E 9".

For control design, the system dynamics are observed at some trim conditions by giving

perturbations to states of the system at that trim. The dynamics associated with these

perturbations are obtained by linearization.

An advantage by linearization is that most of the control methodology is based on linear

equations of motion. A controller is designed initially for the linear system and then tested

for the actual nonlinear system. Yet, there are few disadvantages for this process

* Linearized equations can predict the system performance only in a small range of op-
erations. The linearized equations change as the operating point of system changes,
thus making it difficult for simulating true behavior of system.

Information relating to nonlinearities like hysteresis, backlash, coulomb friction, dis-
continuities etc. may be lost by linearizing the system.

A controller that is good for linearized system might have very poor performance for
the nonlinear equations. An iterative process may be needed to find a controller that
is good for nonlinear equations.









4.1.2 Generic Form of Equations of Motion

The equations of motion in Equations 3.67 and 3.70 can be written in simplified form

using sums of total forces and moments acting on the body.

m(u + qw vr +gSO) = X

m(v + ru -pw gCOSc) = Y (4.2)

m(w + pv qu gCOCO) = Z

IP +qr(I-Iy) =L (4.3)

Iyq+rp( ) =M (4.4)

Ir + pq(Iy I) =N (4.5)

0 C0

0 = 0 CO -SO q (4.6)

{ 1 SO CO r


x COCY COS -So u
y = CTSISOS-C SY SISOSY+CTCO SICO v (4.7)

z COSOCY +SOSY CSOSY SOCY COC w
E

These equations of motions are coupled by the state vector, x, and are dependent on the

control vector, u.

x = {u, v, w,, p, q,r, T, 0, (, x,y, z}
(4.8)
U = {OR1, R2, OE1, OE2, 8R1, 2 R2, 8E1, 8E2, 6c,Fprop

4.1.3 Small Disturbance Theory

The small disturbance theory will be used for linearization of equations of motion.

According to the theory the linearization will be carried about an operating point (reference

flight condition), i.e., the equations thus derived will be valid for the torpedo operating at

and near that value of vector x. The operating point is chosen to correspond to the trim

condition in Equation 4.9.









xo = {uo, vo,wo,, qo, o,Yo, o, o,xo,yo,zo}
(4.9)
= {uo,0, 0 0, 0, 0, 0, 0, 0, 0,}

This corresponds to straight and level flight with constant velocity. As the torpedo may be

traveling at other flight conditions, the linearization at those conditions would be carried out

numerically, which will be explained in later sections. A value of uo is found numerically,

that satisfies the equations of motion for a given value of xo. Then a disturbance of Ax

is given to the equations of motion from the reference condition thus changing the flight

conditions to xo + Ax. Several assumptions are made to carry out the linearization:

* The disturbances from reference flight condition are small. Thus the terms involv-
ing higher order of disturbances A will be neglected. Furthermore first order terms
involving A will be approximated as in Equation 4.10.

Sin(A) = (A)
(4.10)
Cos(A) = 1

The propulsive forces and mass are assumed to be constant.

Planing and immersion forces are neglected for this analysis.

The linearization procedure is resolved for the force equation in bf direction. This

equation relates the force, X, to the states.

m(u + qw -ru+gSO) = X (4.11)

Using the flight condition from Equation 4.9 in Equation 4.11 gives the value of force at

the reference trim condition.
mgSOo = Xo (4.12)

The perturbation equation, i.e., the equation with flight condition disturbed by Ax can be

obtained by substituting the perturbed flight condition into Equation 4.11.

m[ (uo + Au) + (qo + Aq) (wo + Aw) (ro + Ar) (uo + Au)
(4.13)
+gS(o+ A)]= Xo + AX

Equations 4.12 and 4.13 can be combined to give the linearized form of Equation 4.11 for

straight and level flight condition.









m(Au+gAOCOo) = AX (4.14)

Proceeding in a similar way all other equations of motion can be linearized. The lin-

earized equations for straight level flight are shown in Equation 4.15 to Equation 4.18.

4.1.3.1 Force Equations
m(Au +gAOCOo) = X

m(Av + uoAr gAcCOo) = AY (4.15)

m(Au uoAq+gAOSOo) = AZ

4.1.3.2 Moment Equations
IAp = AL

IyAq = AM (4.16)

IAr = AN

4.1.3.3 Orientation Equations


COo

A = Aq (4.17)

A = Ap +TOoAr

4.1.3.4 Position Equations

Ax = -SOn,, AO + COoAu + SOoAw

Ay =Av (4.18)

Az = -CO,,n, ,A SO ,A/ + COoAw

4.1.4 Stability and Control Derivatives

The variations in total force and moment are often difficult to compute.These variations

in forces can be calculated by chain rule for derivatives. As stated in Equation 4.8, these are

functions of state variables x and control variables u. Thus for example AX can be written

by chain rule as in Equation 4.19.









AX =X,,A +XAv +X,, Ai +XpAp +XqAq + XAr

+ XAT + XAO + X~A( + XpropAFrop
(4.19)
+ XO eR1 R+ +XOR 2 + XOEO E1 + xO E2E

+XaR6R1 +XaR28R2 +Xa6E1 +X26E22 +X&6cc

where the subscripted X represents its partial derivative with respect to its subscript.

X, = (4.20)
Iux=x0

Each of these subscripted variables that have a subscript of state variable are known as

stability derivatives and the ones with a control variable as subscript are known as a control

derivative. There can be as many stability derivatives as there are forces and state and

control variables. Many of these are negligible, depending on the reference flight condition.

These dependencies are known usually by experimentation or numerical calculations. For

example, the force, X, is observed to depend mainly on a subset of the state and control

variable. Thus only the stability derivatives that correspond to these variables have to be

retained in Equation 4.19, when straight and level flight is considered.

X = funct(u,w, 6E1,6E2, OE1, E2, 6c,Fprop) (4.21)

The next problem is calculating numerical values of these derivatives. In most cases it

is easy to calculate these numerically or using a symbolic manipulation software. For some

reference points, it is possible to do the calculation manually. The calculation ofX,, will be

done manually for straight and level flight.

X, = (FR,x + FR2,x + FE,x + FE2,x + F,x + Fprop,x) (4.22)

Expressions for each of the terms in Equation 4.22 have been derived in Chapter 3. For

example, the expression for the force on cavitator is shown in Equation 4.23.

cp>cCac CacSc -Sac -Dc (ac, yi)
Fc,x = C8c 0 SSc -SIc CIc 0 0 (4.23)

cp3cSac ScacSpc Cac -Lc (Cc, 7)








In Equation 4.23, ac, Pc, and thus Lc and Dc are the only functions of u. Thus the partial

derivatives with respect to u can be obtained.


-Fc,
au


-ScLcS13j~7- CCkCI3 ~~Ps cc~


-SacSPL^+C(aXCC -C(a, U
-au -u a
-S3 aC 0
S(XCCaPCL I+C(- aa -S-a^-_

0 SSc CPcC c CScSPc -Sac

1 0 -Sic C3c 0

0 Cc CcSocs SacSPc Cac


It can be seen that e-a, aP-, L and Dc are terms required to be calculated. These can be

calculated from equations 3.16 and 3.17, which are restated in Equations 4.25 to Equation

4.27.


tan(ac) = -
Uc

ye
tan(3c) = c+

V; /2+g+w;


(4.25)

(4.26)

(4.27)


The velocity components in Equation 4.27 can be found using Equation 3.2.


CIc 0 SSc

-SICcs LaP- -+ CP3S o


-SP3Sa^CL +CPCCQL-

-Dc (ac,) C6c

0 + 0

-Lc(c c,7) c -S8c

Dc( 0c, l

0

-Lkc( ocY)7C
TH -1 ^r


(4.24)










uc C8c 0

Vc = 0 1

We S8c 0

Uc U

Vc = V +

Wc W
)IB )B


-SSc uc

0 vc

C8c w B

b1 b2 b3

p q r

Xc Yc Zc


Now the velocity components can be obtained for the reference flight condition that is


stated in Equation 4.9.


uc C 0cUO

vc g= 0

WJ SScuo
( ); c-


(4.30)


The variation of ac can be obtained by differentiating Equation 4.25 at reference flight

condition.


sec2 (c)dac


ucdwc wcduc
U2
c


ucdwc wcduc
dac =

CScdw SScduc
d U =
Uo Uo


(4.31)



(4.32)


at (xo, uo)


Similarly the variation of Pc can be obtained by differentiating Equation 4.26 at reference

flight condition.

S (VYdvc vcdVc)
Vc V (4.33)
dvc
=- at (xo,uo)


Now, using Equations 4.28 and 4.29, variation of velocity components of cavitator with

respect to u can be obtained.


CSc


au


(4.34)


(4.28)


(4.29)









Finally, combining Equations 4.32 to 4.34, the variations of ac and 3P with respect to u can

be obtained at reference flight condition.

(xc Cs8c dw S6Sc uc
au uo au uo au
C CSSSc S8cC8c (4.35)
uo uo
=0

apc 1 v,
au Uo au (4.36)
=0

Clearly, two of the derivatives that are required to calculate stability derivatives have been

obtained. It was previously stated that lift and drag are calculated using the coefficient of

lift and drag.

1
L, = pV2S cclc
2 (4.37)
1
Dc = -pV2Sccd
2

These forces can be differentiated by chain rule the derivative would be like in Equation

4.38.
aLC 1 Vc a\cc c1
u- = pSc 2Vclu + Vc u (4.38)
au 2 [ du duc ou c

The Equation 4.38 requires two derivatives. One of the derivatives is calculated in Equation

4.35. The other derivative can be calculated using Equation 4.27.


-Vc a \U2+V2+W21
au au [ c ] w
uc L v, + VC a +WC3
uc + au a+ (4.39)
uV V1









Thus the derivative of the lift and drag forces can be obtained.


3L
c pSScVcl, (4.40)
au
aDc
u= pScVccdc (4.41)
du

Thus all the derivatives required to calculate right-hand side of Equation 4.24 have been cal-

culated. All the terms on right-hand side of the Equation 4.20 can be calculated in a similar

manner. It is tedious to calculate the derivatives analytically in such a way. The complexity

increases for other flight conditions. For practical purposes these derivatives are calculated

using symbolic manipulation software like 'Mathematica' or by using numerical methods.

The numerical methods used to calculate the derivatives have been described in Appendix

B.

4.2 State Space Representation

Equations 4.15 to 4.18 are a complete set of linearized equations of motion for the

torpedo. They can be represented in a more convenient form known as the State Space

Form. The state space equations are a set of first-order differential equations.

x = Ax +Bu

y = Cx +Du

x E 9n,u E Wy E 9m (4.42)

A E n x 9n,B E n x p

C E Z x 9~,D 9Z x W

Equation 4.42 is a generalized form of state space representation for any system. Each of

the terms in the equations has a particular importance for describing the dynamics of the

system.

* State Variable x: The state variables for a system are a set of variables, when known
at time to and along with input u, are sufficient to determine the state of the system at
any time t > to. All the states of the system need not be measurable.








* Input Variable u: This is the control surface deflections.

* Output Variable y: The output variables are the measured parameters. These may
or may not be same as the state variables. The output variables are usually considered
to be measurable but sometimes they are estimated.
The matrices A, B, C and D may either be constant or time-varying functions.

In the case of the supercavitating torpedo, the state vector is of size 12 (n) and the

control vector is of size 10 (p).


x= [A Aw Aq Ag Av Ap Ar A AT Ax Ay Az]
SJ (4.43)
u = c 6E1 8E2 6R1 8R2 OEl 0E2 OR1 0R2 A-prop


Some of these controls may not be needed for some maneuvers. From the linearized equa-

tions it can be observed that the state variables are not coupled by the states {(,x,y,z}.

These four states can be removed from the analysis for now. The system becomes a 8 state

system. These states can be further divided into longitudinal and lateral-directional dynam-

ics. The variables Au,Aw, Aq,AO correspond to longitudinal dynamics, which also means

the dynamics in b1b3 plane. The variables Av, Ap, Ar, Ac correspond to lateral dynamics,

which is the dynamics in flb2 plane. Sometimes the lateral and longitudinal equations can

be decoupled. Thus if the torpedo is making a pull climb/descent to a certain depth, usually

its dynamics can be represented by longitudinal state variables. The plant matrix A can be

divided into four parts.


A F= Along Acoupl (4.44)
A A 1 (4.44)
Acoup2 Alatd


When A is divided as in equation 4.44, where each element is a 4 x 4 matrix, Along would

correspond to longitudinal dynamics and Alatd would correspond to lateral dynamics. Acoupl

and Acoup2 are coupling matrices. It is required that the coupling matrices become negligi-

ble for the equations to be decoupled. If these parts are not negligible, the equations cannot

be decoupled, and a 8 state model will be required to be considered. From linearized







44


equations, the four parts of the A matrix for the torpedo can be written as in Equation 4.45

to Equation 4.48.


x. -qo + x-
xu Zx

qo + -
Mu Mw


0 0
o, oP




Lv Lp

Nv Np-(Iy-Ix)qo
T 1
0 1


Along








Alatd








Acoupl


w -wo+

uo + L
T/ m


MZq
Cy
Cco


Lr




C


xp
m

-vo +
Mp-(Ix--Iz)r
iy
0


Po+
Lw
Ix
Nw
0--


-gCGo+x-

-gCQoSo + -

ME
ly
0


-uo + -gCOoCoo + r.
-(Iz-Iy)0q Le
Ix Ix
Nr N

kI@0 TTo qoC(oSO, -,, ,,,,

vo +x, x
o m m

ZI -gS0oCeo +

0 Mr-(Ix-Iz)po Me
Iy Iy
-Sco -Scoqo Ccoro

I -gSOoSQco +
Lq-(Iz-Iy)ro Lo

Nq-(Iy-Ix)PO No
IS IS
SIoTO0 Scoqo +C4oro


Similarly B is a 8 x 10 matrix, whose elements are just the control derivatives according to

their locations in the matrix. The first 4 rows correspond to longitudinal dynamics and the

last 4 correspond to lateral dynamics.


X8c XFpropX
m m

Blong .

0 *.. 0


(4.49)


4x10


ro+

Po +
Mv

0



Lu
-ro +

Ix
N,

0


(4.45)








(4.46)








(4.47)








(4.48)











YSc YFprop,
m m
Blatd= : (4.50)

0 *** 0
S4x10

Now the complete state space representation for the torpedo can be written as in Equation

4.51 which gives two sets of equations. The first set is the longitudinal equations and the

second set is the lateral-directional equations.



T1
Xlong = Au Aw Aq AO


xlatd = Av Ap Ar AI]

U= 8c 8E1 6E2 8R1 8R2 OEl 0E2 OR1 OR2 A o (4.51)


SXlong Along Acoupl Xlong + Blong X10

Xlatd Acoup2 Alatd Xlatd Blatd
8x1 L 8x8 8x1 L 8x10















CHAPTER 5
CONTROL DESIGN SETUP


This chapter deals with the control design for the torpedo described in previous chap-

ters. Various parameters associated with the control are restated in Table 5.1.

Table 5.1 Control Parameters


Longitudinal Lateral
State u, w, q, 0 v, p, r, 0, Y
Control 6c, ei, 8e2 8rl, 8r2

It should be noted that Y has been included in the states though it was observed in

the state matrices that all other variables are independent of Y. It will be seen later that

the inclusion of Y in the feedback states helps in improvement of performance. Also,

for longitudinal control, two elevators and the cavitator are required. Similarly for lateral

direction, the rudders should be enough for control.

There are various control methods, like linear quadratic regulator (LQR) synthesis, u-

synthesis etc., which can be used to design a controller. Each of these methods has advan-

tages and disadvantages. LQR method gives a constant gain controller which is based on

minimization of a quadratic performance index and considers the problem of robustness

only in terms of gain and phase margins. p-synthesis deals with robustness with respect to

a wide variety of uncertainties to minimize an infinity-norm matrix but the resulting con-

troller can be of high order. Regardless of complexity and robustness, each design method

presents difficulties. The LQR method was chosen for controller synthesis for the torpedo.

This method was chosen because

1. It is easy to implement in a nonlinear simulation.
2. The resulting robustness for the LQR controller was acceptable.









3. It is straight-forward to vary some design parameters to achieve performance goals.

This chapter explains various problems associated with the control synthesis and the system

model used for synthesis of the controller.

5.1 Open-loop Performance

Initially the equations of motion for the torpedo have been linearized for straight and

level flight at a forward velocity of 75 ms-1

x= {75, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 (5.1)

It is found that the cavitator and two elevators are sufficient to maintain the above trim.

It is also assumed that the value of propulsive force required to maintain this trim is fixed

at the required value.

U = {Sc, 5el, 5e2, 8rl, 5r2,Fprop}
(5.2)
= {0.0067, 0.0106, -0.0106,0,0, 4.0023e + 03}

where the deflections are given in radians, and Fprop is in Newtons. As these parameters

are obtained numerically, it may not be a perfect trim. The system may have some non-

zero accelerations, and consequently may tend to deviate from straight and level flight. To

check this, the open-loop dynamics are simulated at this condition without any feedback.

Figure 5.1 shows the Simulink model used for open-loop simulation. The control surface

deflections are fixed at their trim values for this simulation. The closed-loop system is

obtained using the equations of motion that were derived in Chapter 3. The state derivatives

are then integrated to obtain the state at the next instant.

Figure 5.2 shows the open-loop response for the torpedo at this trim condition. It can

be seen that the open-loop 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 non-zero states.

The state and control matrices obtained for this trim condition are shown in Equations

5.3 to 5.6. There are 5 control variables, {Sc, 8el, 8e2, 8rl, r2}. The matrices corresponding
























S NL Equation Integrator L
J Torpedo
at Foi

-- simt

Clock Time
Figure 5.1 Simulink Model for Open Loop Simulation

Figure 5.1 Simulink Model for Open Loop Simulation


vIvvvv~vvwovvvv'fvvv


I0 20 40time(s40
4-
35
3
1--2 5

S2
0-1 5
1 -


80 100 -0 040
80 100 0


20 40 60
times)


Figure 5.2 Open-Loop Response for Torpedo: w,p, q


Control
Trim


state
r Plotting


20


15


5"


5


80 100









to the lateral dynamics are of dimension 5 because the state Y is included in the lateral

dynamics. Thus the lateral states are now {v,p,r,1, Y}.


-4.5204 1.5417 1.3110 -9.8100

-0.2616 -15.7648 78.5888 0
Along = (5.3)
0.0000 1.2077 -3.5614 0

0 0 1.0000 0



-32.3010 69.0608 -69.0608 0 0

-406.0942 -303.3736 303.3736 0 0
Blong= (5.4)
158.4675 -45.1531 45.1531 0 0

0 0 0 0 0


-12.0422 -0.0002 -71.6011 9.8100 0

-0.1813 -54.2281 0.3004 0 0

Alatd= 1.1437 -0.0025 -1.2528 0 0 (5.5)

0 1.0000 0 0 0

0 0 1.0000 0 0


0 0 0 -366.60511 366.60511

0 -14297.086 -14297.086 -17276.994 -17276.994

Blatd= 0 -1.4129523 -1.4129523 54.561629 -54.561629 (5.6)

0 0 0 0 0

0 0 0 0 0

The longitudinal eigenvalues are {-21.1414,-4.5137,1.8262,-0.0178} and the lat-

eral eigenvalues are {0, -54.2289,0.0002,-6.6472 + 7.2683,-6.6472 7.2683}. The

eigenvalues clearly show that the system is unstable. It can also be seen that the longi-

tudinal dynamics have no oscillatory modes. Figure 5.3 shows the variation of eigenvalues










Longitudinal Egienvalues for u=60:5:90 Lateral Egienvalues for u=60:5:90
1 1

0.5- 5







-0 -20 Rea-10 long) 0 10 0 -60 -20 0 20
Reral (oneal (at)

(a) Longitudinal (b) Lateral

Figure 5.3 Variation of Eigenvalues with Change in Velocity


for the torpedo for different velocities. State values are fixed except for forward velocity,

which is varied from 60 ms-1 to 90 ms-1. The figures show that the variation is small and

most of the eigenvalues stay in negative half of complex plane.

5.2 Closed-Loop 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 closed-loop system.









Table 5.2 Control Constraints


Cavitator Deflection -15 < 8 c < +15
Cavitator Rate -25/s <, 8 < +25/s
fins -60 < 8f < +600
Fin Rate -25/s < 8f < +250/s
Thrust 0 < Fprop < 30000N


Various bounds are placed on the control surface deflections and rates. These bounds

are listed in Table 5.2. The bounds are included in the nonlinear simulations and it is

required that there is no saturation of deflection or the rates at least for the commands

having the rate 30 deg/s.

An actuator model is included in the system so as to constrain the rates of the control

surface motion. This actuator is realized as a low pass filter, Ac = 80 Addition of this

filter synthesizes a controller that has slower control deflections.

Let qcomm(t) be a function of time, defining the desired pitch rate profile. The aim of

the controller is to find a control law u(t) that yields an achieved pitch rate, qachi(t), to

minimize the optimization problem stated in Equation 5.7.

find u(t)

that minimizes (t) = |qachi(t) qcomm(t)

subject to Umin < U < Umax (5.7)

Umin < ll <- lmax

x = Ax + Bu

where, umin and Umax are lower and upper bounds on control deflections. The quantities

Umin and Umax are lower and upper bounds on control deflection rates.
The problem becomes a disturbance rejection problem, when the commanded variable

is 0 for all time. This is an optimization problem,where the state space equation is an

equality constraint and the control surface bounds are inequality constraints.









5.3 Robustness of the Controller

A control system that is designed to accommodate the uncertainties in a mathematical

model is called a robust control system. Usually the response of a model does not accurately

match the response of the true system. A robust control system should give the desired

performance not only during the simulations of the model, but for the true system also.

Various parameters can be introduced in the model to simulate uncertainties. Random

noise can be added to output signal to simulate measurement errors, the gains in signals

can be changed to model uncertainty in response etc. Gain and phase margins are generally

used to predict the robustness of a control system. Physically, these are just the factors

by which the feedback gain can be increased and yet have a stable real system. A formal

definition of these can be given by using a frequency analysis for a feedback system.

5.3.1 Gain Margin

Figure 5.4 shows a typical feedback system involving a plant, P, and a controller, K.

The gain for the system in dotted region is known as the loop gain. The loop gain is

important because it determines stability. Errors in the predicted loop gain could cause

errors in predicted stability. The gain margin is the largest factor by which this gain can be

increased and still have a stable system. Physically, it means if the response of the torpedo

for a given elevator input is higher by a factor of the gain margin, the torpedo is still stable.

The gain margin is usually expressed in decibel (db) units, and can be easily obtained from

the Bode plots for the system. The gain margin is a measurement of the magnitude on the

Bode plot, at the point where the phase is 180.

5.3.2 Phase Margin

Gain is a valid robustness criteria when the system has real eigenvalues. But usually

the eigenvalues have imaginary components and thus the phase is also a concern. Phase

margin is the measure of the maximum possible phase lag before the system becomes

unstable. From the Bode plot, phase margin is the phase when the magnitude of the gain is

zero.




















Figure 5.4 Loop Gain

5.3.3 Uncertainty In Parameters

Another factor that can determine the robustness of a controller is its response to errors

in known parameters. As stated earlier, the coefficients of lift and drag are calculated from a

CFD database. The accuracy of the model depends on accuracy of this data. Robustness of

a controller can be assessed by introducing errors in the data and checking how it performs.

The following variations have been introduced in the system to check for performance of

the system with intrinsic uncertainties:

* 120% error in C1 of Cavitator.
* 120% error in Cd of Cavitator.
* 120% error in C1 of all the Fins.
* 120% error in Cd of all the Fins.

5.3.4 Controller Objective

In terms of robustness, the controller is required to meet various performance objec-

tives. These objective can be summarized as:

* The closed-loop system should have a gain margin of at least 6 dB.
* The closed-loop system should have a phase margins of at least 45 deg.















CHAPTER 6
LQR CONTROL

6.1 LQR Theory

The tracking problem, like the one given in Equation 5.7, can be solved by using a

combination of feedback and feedforward control [10]. The Linear Quadratic Regulator

(LQR) problem is to find an optimal feedback matrix K such that the state-feedback law

u = -Kx minimizes the linear quadratic cost function shown in Equation 6.1.


J(u)= (xTQx + uTRu + 2xTNu)dt (6.1)
0

The basic idea of LQR control is to bring the state of the system close to zero. A linear

system can be represented in the state space form as in Equations 6.2 and 6.3. The matrices

A and B are the state and control matrices. The variable x represents the state vector, y is

the output vector and u is the input vector.

S= Ax+Bu (6.2)

y=x (6.3)

The LQR controller is realized by a constant gain matrix K, such that the feedback

u = -Kx makes x go to zero. By a modification to this law, the LQR method can also be

used for tracking. The state vector x is of size n.

x= {X1,X2,"--, XnT (6.4)

Let the tracking problem be for the state xl to track a step command rl. The idea is to

make (xl rl) go to zero using a LQR controller. The new control law can be chosen as in

Equation 6.17.










x1 rl

X2
u= -K (6.5)


S x

Equation 6.2 can be rewritten by substituting the new control law.

x = Ax + Bu

x1 r1

X2 (6.6)
= Ax BK


\ x,


For simplicity, assume that there is only one control, u (this is different from velocity u).

The controller K is of size n x 1 and it can be expanded in its elements.

K =[kl,k2,... ,k,] (6.7)

Equation 6.6 can be rewritten by substituting the K in its expanded form.


xi -ri

X2
k = Ax B[kl, k2, ,kn]

(6.8)
\ xn

=Ax BKx +Bklrl

=(A BK)x +Bkirl

It should be noted that the command rl is a step command. The steady-state dynamics of

the system can be obtained from Equation 6.8.


(o) = (A -BK)x(-) +Bkirl


(6.9)









The error dynamics can be obtained by subtraction Equation 6.9 from Equation 6.8.

k(t) 2(-) = (A BK) (x(t) x(-)) (6.10)

e = (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 steady-state error e, which is made zero using the regulator. Figure 6.1

shows the block diagram for this closed-loop system. It is required that the closed-loop

system has an integrator somewhere so as to make the steady-state error go to 0 [10]. That

is, e has to go to zero rather than e so as to achieve good tracking. Figure 6.2 shows the new

configuration of a system that has no integrator and thus an integrator has to been included

during design. Thus, the integral of the actual error has to be made to go to zero so as to

achieve a good tracking.

e= (ri-xl) (6.12)

The state space equation for the system with this modification can be written.

S=Ax+Bu

e= ri -xi (6.13)

= rl Cx

where xl = Cx. It can be seen that the error equation is similar to state equation. Thus e

can be considered as another state, .i.e, the system now has n + 1 states with state vector

x = {xl,x2, ,Xn,, T. So a new formulation of the state space equation can be written,

A 0 B 0
X= x+ u+ ri (6.14)
-C 0 0 1

x=Ax+Bu+Ir (6.15)

which is similar to Equation 6.8. The error dynamics of this system represent the form of

state space equations, for which a LQR controller can be derived. The LQR controller K,

will be a constant matrix of size n + 1 as the system now is of size n + 1.






57

K= [kl,k2,.. ,kn lk,+I] (6.16)

Then, the new control law can be written as in Equation 6.17.

u = -K

x
= -[k, k2, kn I kn+1
S(6.17)

= -[ki,k2,--- ,kn]x + [-kn+l ]

=-Kx +ki

which is represented in Figure 6.2.


r +0 x =Ax+Brk X
y=x




K


Figure 6.1 Controller for Tracking when Plant has an Integrator.



k1 =Ax+Bu
+0 f+ y=x





K



C -


Figure 6.2 Controller for Tracking when Plant has no Integrator.









6.2 Control Synthesis

The torpedo system does not have an integrator in the system. A tracking controller can

be obtained from LQR method by the process described in Section 6.1. The controller is

obtained for a trim state of straight and level flight at 75ms -1. The linearized dynamics are

first separated into the longitudinal and lateral dynamics as given in Table 5.1. The controls

used in longitudinal direction are the cavitator and 2 elevators. The controls in lateral

direction are the 2 rudders. It is observed that for straight and level flight the longitudinal

and lateral dynamics are practically decoupled.

Once the state and control matrices have been obtained, the main variables that the LQR

controller depends on are the weighting matrices Q, R and N. In this case the cross coupling

matrix N is chosen to be 0. The matrices Q and R penalize the cost function for higher state

and control values respectively. A higher value in Q matrix would cause a better track-

ing. A larger R would constrain the control surface deflection. An optimum combination

of the matrices is obtained iteratively, so as to get good tracking with achievable control

deflections.

The matrices for longitudinal pitch rate tracking are given in Equation 6.18.

Qlong = diag([0, 0, 0, 0, 10])
(6.18)
Rlong = diag([5,4])

The first four numbers in the Qiong correspond to weightings on the four longitudinal states.

They are chosen to be 0. We do not want to restrict the states from changing. This is

especially important for weightings on q and 0. A weighting on these variables would

restrict them from changing. The last number, 10, is weighting on the error. This is chosen

to be high to penalize the tracking error. A higher value of weighting would give a better

tracking, but it is seen that it would require very high control rates. The first number in

the weighting matrix Rlong, 5, corresponds to cavitator weighting and the number 4 is for









elevator weighting. Elevator weighting is chosen to be smaller so as to encourage the

controller to use more elevator than the cavitator. This gives a more stable performance.

The control matrices obtained for the longitudinal dynamics are given in Equations 6.19

and 6.20.

1.1182
k1 = (6.19)
-0.9681

-0.0000 0.0040 0.1016 -0.0000 1.4195 -1.1184
K = (6.20)
0.0000 -0.0042 -0.0995 -0.0000 -1.3981 1.1308

Similar process is involved in the design of the lateral controller. Initially the lateral

controller is designed with only four state feedback, and Y is neglected in the feedback. In

this case it is found that the torpedo has high sidewash and deviates considerably from the

original path, even when the pitch angle is 0. To avoid this, Y is included in the feedback

states. It is also observed that a penalty on the yaw motion causes the controller to com-

mand a very fast control surface motion. Also, a continuous correction of control surface

deflection is required to prevent the yaw motion entirely. Thus an optimum combination of

the weighting matrices is obtained that would prevent a very high yaw motion but would

still have slow control surface motion.

Qlatd = diag([0, 0, 0, 0, 0,.1])
(6.21)
Rlatd = diag([1000, 1000])

The first 5 numbers correspond to 5 states and the last number is weighting for the error.

The Rlatd is of dimension 2 as only the rudders are included in the synthesis. The weighting

on the rudders is high as it is observed that the roll rate is very sensitive to the rudder

deflection. The control matrices obtained for the lateral dynamics are given in the Equations

6.22 and 6.23.










-0.0071
k, = (6.22)
-0.0071


=0.0005 -0.1253 -0.0132 0.0019 -0.0000 (6.23)
K=10.3 (6.23)
0.0005 -0.1254 -0.0132 0.0026 -0.0000


The feedback matrix K for lateral dynamics is of size 2 x 5, which is shown in Equation

6.23.

6.3 Nominal Closed-loop Model


6.3.1 Model

Figure 6.3 shows the eigenvalues for the closed-loop longitudinal and lateral systems.

It can be seen that both systems are stable as all the eigenvalues are in the left half of

the complex plane. Also, each of the dynamics has one eigenvalue at the origin, which is

introduced due the integrator in the system.

15 8
6
10
4
5 2


-5 -E -2
-5
-4
-10 -6
-6

-10o0 -800 -600 4P(lon o00 0 200 -1o00 -800 -600 4 -d00 0 200


(a) Longitudinal (b) Lateral

Figure 6.3 Eigenvalues for the Closed-loop System


6.3.2 Linear Simulations

The response of the closed-loop linear system has been shown in Figures 6.4 to 6.8.

The simulations for lateral and longitudinal systems have been carried out separately as the

linear system is decoupled into lateral and longitudinal.









Figure 6.4 shows the tracking obtained for a 15 deg/s pitch rate command. The com-

mand is achieved in 0.17s with an overshoot of 3.95% and with no steady-state error. Fig-

ures 6.5 and 6.6 show the control surface deflections and rates required to achieve the com-

mands. Though there are some quick deflections, the rates are still under the constraints.


5 Time(s) 10 15 0 5 times) 10

Figure 6.4 Pitch Command Tracking for Linear System : q, u


5 times) 10 15 0 5 times) 10

Figure 6.5 Pitch Command Tracking for Linear System : 8c, 8c


Figure 6.7 shows the roll rate tracking obtained for a 15 deg/s roll rate command. The

command is achieved in 0.53s with no overshoot and a 0 steady-state error. The variation

of the yaw rate is also shown in the figure and it can be seen that the yaw motion is coupled

with the roll. At the end of the roll doublet, the torpedo has a non-zero yaw angle thus

changing the direction of motion. The control surface deflections required for the roll rate

command are shown in Figure 6.8. The rudder deflection is small for reasons explained in

the next chapter.






























































-0.01

-0.02

-0.03
0


20
I/)
ao) 1
10



-10
-1" 0


S5 times) 10 15 0 5 times) 10

Figure 6.6 Pitch Command Tracking for Linear System : 8el, 8el


Commanded
- Achieved


5 10 15 -) 5 10
Time(s) times)

Figure 6.7 Roll Command Tracking for Linear System : p,r


U)
0)
0)

t0)
g!
e


-n .


5 times) 10 15 0 5 times) 10

Figure 6.8 Roll Command Tracking for Linear System : 6r1, 6r


rf-


-r--L r









6.3.3 Gain and Phase Margins

The LQR tracking system shown in Figure 6.2 is obviously more complex than the

system shown in Figure 5.4. Thus, the loop gain can be defined in many ways in this case.

The block diagram can be broken at different points so as to simplify it to the form shown

in Figure 5.4. Figure 6.2 is redrawn in Figure 6.9 which shows the possible breakpoints for

this system. For understanding, the output of plant P is divided into two parts, one is the

achieved value of the commanded variable (ra) and the other is remaining states of the plant

P (x). The break points are numbered 1 to 3. The system can be broken at each of these

points to give a loop gain. These gains will be named outer-loop, inner-loop and all-loop

gains respectively.

rl + 1 x -Ax + Bu ra
-- K -+v






2



Figure 6.9 Breakpoints for Calculating the Loop-Gain for a Tracking Controller

Gain and Phase margins for each of the above possible break points have been calcu-

lated for both the longitudinal and lateral controllers. Table 6.1 lists the gain and phase

margins for the torpedo with LQR controller that was obtained in previous sections. All

margins are quite high and meet the desired conditions of 6dB for gain and 45 deg for phase

margin. Figures 6.10 to 6.15 show the corresponding bode plots for the data given in Table

6.1.

Also, the lateral controller is unable to stabilize the unstable spiral mode. Thus the

closed-loop system is inherently unstable due to this pole and would consequently have

negative gain margin. Numerous simulations show that the affect of spiral mode is negligible,









Table 6.1 Gain and Phase Margin with LQR Controller


Longitudinal
Gain Margin(db) Phase Margin (deg)
1 21.056(at 47.498 rad/s) 64.846(at 9.0625 rad/s)
2 327.87(at 0 rad/s) 77.118(at 25.925 rad/s)
3 o 57.606(at 20.845 rad/s)
Lateral
Gain Margin(db) Phase Margin (deg)
1 22.964(at 0 rad/s)
2 250.51 (at 0 rad/s)
3 50.36 (at 0 rad/s)


i.e., the time to double for the instability is considerably larger than the maneuvering time

of the torpedo. So, the closed-loop system model is reduced by removing the spiral mode

from the model. The gain and phase margins in Table 6.1 are for this reduced-order system

and reflect the robustness of the dominant dynamics.

Bode Diagram
Gm = 21.056 dB (at 47.498 rad/sec), Pm = 64.846 deg (at 9.0625 rad/sec)
0-

S-50



1 -
-100
-90
a-135
P-180
--225
-27 01 /s-c102 3
10 10 (rad/sec) 10 10

Figure 6.10 Gain and Phase Margin: Longitudinal Outer-loop


6.4 Perturbed Closed-loop Model

A perturbed system model is formed by adding an error to the values of coefficients

of lift and drag for the fins and cavitator. New values of trim deflection are obtained for

the perturbed model and thus a new set of A and B matrices is obtained. Tables 6.4 to 6.9









Bode Diagram
Gm = 327.87 dB (at 0 rad/sec), Pm = 77.118 deg (at 25.925 rad/sec)
50

0

-50

45
0-
-45
-90_
-135-
-1804 1---3
100 10 (rad/sec) 102 10

Figure 6.11 Gain and Phase Margin: Longitudinal Inner-loop

Bode Diagram
Gm = Inf, Pm = 57.606 deg (at 20.845 rad/sec)
100


0


100----
-90


1-135


-180
100 (rad/sec) 102

Figure 6.12 Gain and Phase Margin: Longitudinal All-loop


show the percentage variation of the elements of A and B matrices for a 20% change in

coefficients of lift and drag of cavitator and fins. Few elements in the state and control

matrices change. In most cases, the change in elements of A and B matrices is a linear

function of the change in a coefficient. For example, in Table 6.4 there are 8 terms that

show a variation due to a 20% variation in coefficient of lift of the cavitator. The term

A(3,1) shows a large variation but its numerical value is negligible. The term A(3,2) shows

a 34% variation but this term is also small compared to other terms. Remaining terms in

the matrix show very small variation. Some terms in the B matrix show a 20% variation.










Bode Diagram
Gm = 22.964 dB (at 0 rad/sec), Pm = Inf
-20

-30

8-40

-50
180


S135


90
100 10 (rad/sec) 102 103

Figure 6.13 Gain and Phase Margin: Lateral Outer-loop

Bode Diagram
Gm = 250.51 dB (at 0 rad/sec), Pm = Inf
-20

--30

1-40

-50
90
145
0-
Q-45 -
-90 M o
10 101 (rad/sec) 102 103

Figure 6.14 Gain and Phase Margin: Lateral Inner-loop


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 closed-loop system has good gain and phase margins,

this effect on B matrix should not be of much concern.

6.4.1 Model

Figure 6.16 shows the eigenvalues for the perturbed closed-loop longitudinal and lateral

systems. An error of -20% is included in the value of coefficient of lift for the fins. It can











Table 6.2 Percentage Variation in A Matrix due to 20% Variation in clc


u w q O v p r c
u 0.46 0.62
w 5.52 6.86 1.58
q 1.05e5 34.8 13.58


v
p
r
~( ~~
Ty


Table 6.3 Percentage Variation in B Matrix due to 10% Variation in cc


8c 8el 8e2 861 8r2
U
u
w 20
q 20
0
v
p
r
(D-
Ty


Table 6.4 Percentage Variation in A Matrix due to 20% Variation in cdc


u w q O v p r (Y
u 10 -6 7.6
w 1.54 0.36
q 658 7.8 3


v 2 20 0.4
p 2 -15.6
r -10 0.8 8.4
~( ~
Ty










Table 6.5 Percentage Variation in B Matrix due to 20% Variation in cdc


8c 8el 8e2 8r1 6r2
u 20
w
q 0.12
0
v
p
r
~(~
Ty


Table 6.6 Percentage Variation in A Matrix due to 20% Variation in clfi,


u w q O v p r (
u 1.22 0.64
w 14.4 10.0 -0.9
q -le5 -60 3


v 16.2 -1.2
p 19 36.0
r 25.4 0.94 10.2






Table 6.7 Percentage Variation in B Matrix due to 20% Variation in clfi,


8c 8el 6e2 6r1 6r2
u
w 20 20
q 20 20
O
v 20 20
p 20 20 20 20
r 20 20
~(~
Ty









Bode Diagram
Gm = 50.36 dB (at 0 rad/sec), Pm = Inf
-50

-60

C -70

-80
180


S135



10 10 (rad/sec) 10 10

Figure 6.15 Gain and Phase Margin: Lateral All-loop


Table 6.8 Percentage Variation in A Matrix due to 20% Variation in cdfi,


u w q O v p r ( Y
u 9.4 24.0 12.4
w 1.34 -0.12
q -68 -2.6 0.4

v 20 1.76 -2.3 -0.13
p -2.3 1.12 -0.62
r 2.74 18.26 1.12






be seen that the longitudinal dynamics show some perturbation in the damping while the

lateral system relatively unchanged.

6.4.2 Linear Simulations

Figures 6.17 to 6.19 show the response of the perturbed system for a 15 deg/s pitch

rate doublet command. The perturbation to the system is a 20% error in the value of the

coefficient of lift of the fins. It can be seen that the performance criteria are met even in

the case of the perturbed system. It can also be observed that the overshoot is increased

for the perturbed system. The performance is achieved at the cost of small perturbations in











Table 6.9 Percentage Variation in B Matrix due to 20% Variation in cdfi,


8c 8el 8e2 8rl 8r2
u 20 20
w
q 1.36e-2


v

p
r 20 20 12.6e-3 2.6e-3

_(____ ______


15 8
6
10
4
5-
2-


E -2
-5
-4
-10
-6

-1100 -800 -600 -400 -200 0 -1400 -800 -600 -400 -200 0 200
Real Real

(a) Longitudinal (b) Lateral

Figure 6.16 Eigenvalues for the Perturbed Closed-loop System: 20% Error in clfi,


other states of the system. As the control effectiveness of the control surfaces is changed,

the amount of control surface deflection is also changed by a constant factor.

Figures 6.20 to 6.21 show the response of the perturbed system for a 15 deg/s roll

rate doublet command. The perturbation to the system is a 20% error in the value of the

coefficient of lift of the fins. It can be seen that the performance criteria are met even in

case of perturbed system. In this case also, it can be observed that there is a perturbation in

other states.
















- No Uncertainty
+20% in cl
20% in clf


/ -- No Uncertainty
V + +20% in cl
...... -20% in cl
5 times) 10
times)


Figure 6.17 Pitch Command Tracking for Perturbed Linear System : q, u


o0 05


-40
5 10 15 0 5 10
times) times)


Figure 6.18 Pitch Command Tracking for Perturbed Linear System : 8c, 8c


5 times) 10
times)


20

U)
) 10

-0
1 0

-10


S-20
15 0


- No Uncertainty
+20% in cl
...... -20% in clf


5time(s) 10


Figure 6.19 Pitch Command Tracking for Perturbed Linear System: 6ei, 5e1


I 74
E
"-73.5
:3


5 Time(s) 10
Times)


72.5

S 720
15 0


1 5


1


S05


-" 0


-05











20
15
10
o 5
0
" o
-5
<-10


- No Uncertainty
_ +20% in clf
-20% in clf


Figure 6.20 Roll Command Tracking for Perturbed Linear System : p,r


-- No Uncertainty
-002 ___ +20% In clf -0.1
... -20% In clf
-003 -0.15
0 5 10 15 0 5 10 15
times) times)


Figure 6.21 Roll Command Tracking for Perturbed Linear System : 8rl, 6r



6.4.3 Gain and Phase Margins

Table 6.10 lists the gain and phase margins for the perturbed closed-loop 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 all-loop.

The last value is increased to o showing an improvement for the perturbed system.

From the analysis of the perturbed closed-loop system it can be said that the linear

model is robust to various uncertainties in the system.






73


Table 6.10 Gain and Phase Margin for Perturbed Closed-loop System: 20% error in clfin


Longitudinal
Gain Margin(db) Phase Margin (deg)
1 21.193(at 49.599 rad/s) 68.981(at 8.7966 rad/s)
2 320.33(at 0 rad/s) 77.605(at 25.925 rad/s)
3 Io 60.305(at 22.552 rad/s)
Lateral
Gain Margin(db) Phase Margin (deg)
1 24.391(at 0 rad/s) o
2 278.23 (at 0 rad/s) C
3 (at 0 rad/s) Co















CHAPTER 7
NONLINEAR SIMULATIONS


The controller for longitudinal and lateral dynamics have been obtained separately. That

is, the longitudinal controller is to achieve a required pitch rate and the lateral controller

achieves a given roll rate. Once the controllers have been found for linear systems, they are

employed with the nonlinear torpedo model and the performance is checked using numeri-

cal simulation. Figure 7.1 shows the complete nonlinear simulation model for the torpedo

with the LQR controller.

The nonlinearity in the model is provided by both the aerodynamics and control surface

constraints. These constraints, given in Table 5.2, are not directly included in the linear

model. So it is important to find their effects on the nonlinear simulation. It should be

noted that there is a constraint on thrust, but thrust is assumed to be constant with time.

The thrust constraint is used to find a trim value for various trajectories.

7.1 Nonlinear Simulations for Nominal System

The simulations have been carried out for various commands with the nominal sys-

tem. The 'Nominal system' is the nonlinear torpedo model assuming that it is completely

accurate. No uncertainty is included in the model. The simulations show good tracking

response while meeting all the performance criteria.

The response of the vehicle to a longitudinal command is simulated and shown in Fig-

ures 7.2 to 7.5. These figures show the response for a pitch rate doublet of 15 deg/s. The

rise time for the pitch rate command of 15 deg/s is 0.18s and there is an overshoot of

11.53%. The steady-state error is .8%. The controller is able to command pitch rates as

high as 30 deg/s. It is observed that the vehicle motion is confined to longitudinal plane


















































Derivative r2dl


Figure 7.1: Complete Nonlinear Simulation with LQR Controller











x10- 20


10


0 0-






0 5 times) 10 15 0 5 times) 10 15


Figure 7.2 Pitch Command Tracking : p, q

1 20


0.5 10I

0) 0
00


-0.5 -
-20

-10 -30
0 5 times) 10 15 0 5 times) 10 15


Figure 7.3 Pitch Command Tracking : c, 8c


only. This shows that the controller allows pure longitudinal motion to be uncoupled from

the lateral motion.


2 20

15
1.5 10

5
Cn
z 0
0.5 16 -5
-10

0 -15
0 5 times) 10 15 0 5 times) 10 15


Figure 7.4 Pitch Command Tracking : e,1, e1










2 10-7
-300
15
-200
1-20 starting point
0100-




800
600
5 10 15 0 0 4
times) 10o 20 0 x(m)

Figure 7.5 Pitch Command Tracking : 8ir, {x,y,z} Trajectory


The response of the vehicle to a lateral, roll rate, command is shown in Figures 7.6 to

7.10. A roll rate command of 15 deg/s is achieved in .52s with an overshoot of 0% and

a steady-state error of 0.09%. The controller is able to command a roll rate motion of as

high as 50 deg/s before a saturation of control surface rate is reached. It is observed that

there is some longitudinal motion in this case. This longitudinal motion has been reduced

by inclusion of Y in the feedback states to the controller. It can be seen that the rudder

deflection for a roll rate command is small. This is expected as the terms corresponding to

roll rate from rudder are an order of 3 times larger than the terms corresponding to pitch

rate from elevators. It is assumed that the control surface deflection is achievable.

20 0.5
120

00
SI -0.5

0
i-1.5
-10
-2

-20 5 times) 10 15 2 5 times) 10 15



Figure 7.6 Roll Command Tracking: p, q


Figures 7.11 to 7.15 show the response of the torpedo for a combined roll and pitch

rate command, similar to a windup turn. In this case, the torpedo is given a 5 deg/s roll



























5 times) 10 15 5 times) 10


Figure 7.7 Roll Command Tracking: c, 8c


-2

V -2.5

1 -31
5 times) 10 15 0 5 times) 10 15


Figure 7.8 Roll Command Tracking: 8el, 8el


-0.5
5 time(s)10 15 times)


Figure 7.9 Roll Command Tracking: 8,1, 6,1


-0
to
0.4

0.3

0.2


-0.01

-0.02

-0.03
-*













1000


500-
Starting Point
0
1000 1500
500 1000
500
y(m) 0 0 x(m)


Figure 7.10 Roll Command Tracking: {x,y,z} Trajectory


rate command from 2 to 12 seconds, 5 deg/s pitch rate command from 12 to 22 seconds,

-5 deg/s pitch rate command from 22 to 32 seconds and then -5 deg/s roll rate command

from 32 to 42 seconds. As the vehicle motion is a little different from the actual trim, it can

be seen the vehicle has considerable sidewash. Despite this sidewash, the controllers give

a good tracking performance. The rise times for roll and pitch commands are .5s and 0.22s

respectively. The overshoot for the same are 0% and 0% respectively.



4 I 4


(n i(n
0 0,

_-_2 '--2

-4 -4

-6o 10 20 30 40 50 0 10 20 30 40 50
times) times)

Figure 7.11 Roll & Pitch Command Tracking: p, q


7.2 Nonlinear Simulations for Perturbed System

The performance of the controllers is studied using the simulation with a perturbed

system model. An error is assumed in the values of various coefficients and a correction

factor is added. Response of the closed-loop nonlinear system is not much affected by the

variations in coefficients of lift and drag. It is observed that the controller commands the














2i
0 --- -- ------


10 20 30 40 50
times)


Figure 7.12 Roll & Pitch Command Tracking: 8c, 6c




6


4
(n
2-

0

-2


40 10 20 30 40 50
times)


Figure 7.13 Roll & Pitch Command Tracking: e,1, e,


-0.01 -0.04
0 10 20 tme(30 40 50 0
times)


10 20time(s 0 40 50
times)


Figure 7.14 Roll & Pitch Command Tracking: 8,r, 6












2000-

1000-
Starting Point
0-

-1000>
2000
1000 1500
1000
0 500
y(m) -1000 0 x(m)


Figure 7.15 Roll & Pitch Command Tracking: {x,y, z} Tracking

757 02

757-6 ____------------------ 01 ------- ----------
t0-

755
756 / 014 .
753 E
3-02
75 2


.75. __ 20% In cl. -0. 4 1 __ 20% In cl
72 in 04 20% in ci
-20% in cl-20% in cl
74 -051
074 1 0 5 10 15 20
time time


Figure 7.16 Response for 20% Variation in clfi,: u,w



system to a new trim state which is also a straight and level flight, with change in speed

and control deflections. After that, the system follows a pitch or roll command as well

as before. Figures 7.16 to 7.19 show the response for one such case. In this case a roll

doublet is commanded to the system, and there is an error of 20% in the value of clfi,.

It can be clearly seen that the vehicle has gone to another trim state and then it follows the

command equally well. There is almost no change in the trajectory of the vehicle. The

control surface deflections are similar with a constant offset. Such response has also been

checked for other cases. The affect of error is similar in all cases.

By above analysis it can be said that the LQR controller designed for the torpedo for

pitch and roll rate tracking commands is fairly stable and can be expected to achieve good

performance for the real torpedo.















- No Uncertainty
20% in cl
.....-20% in cl


-0) 5


--1.5

-2

-2.5

0 -3
20 0


Figure 7.17 Response for 20% Variation in clfin: p, q


15 20


Figure 7.18 Response for 20% Variation in clfin: 8c, 861






No Uncertainty
20% inclf
1000 -20% in clf


E 500
N

0
1000


500
y(m)


1000
500
n x(m)


Figure 7.19 Response for 20% Variation in clfin: {x,y,z} Trajectory


No Uncertainty
20% in clf
-20% in clf

15


5 10
time















CHAPTER 8
CONCLUSION


8.1 Summary


A dynamical model for a supercavitating vehicle has been obtained. The vehicle is

found to be open-loop unstable, and a controller for stabilizing the pitch and roll rate motion

has been obtained. The LQR controller shows good tracking performance for the vehicle

and all the control objectives are met. The controller is also found to be robust to errors in

cavity prediction and velocity changes. This robustness is further demonstrated by the fact

that the closed-loop system has high gain and phase margins.



8.2 Future Work


The dynamical analysis of the vehicle has been derived with an assumption that the

cavity shape is fixed. The open-loop cavity dynamics need to be modeled and included in

the synthesis.

Robust control methodologies like p-synthesis can be applied to obtain a more robust

controller. A robust control design could include the uncertainties in the model during the

control synthesis.

The LQR controllers obtained are typically known as the 'inner-loop' controllers. An

outer-loop controller is also needed for guidance and navigation. The idea is that the outer-

loop controller can be modeled for tracking the trajectory in space, based on the closed-loop

dynamics of the inner-loop model.















APPENDIX A
REFERENCE FRAMES AND ROTATION MATRICES


X3
Y3




YX2



X1 ,X2

Figure A. 1 Rotation of Frames


Figure A.1 shows two frames X{xix2X3} and Y{yly2y3}. Y is rotated from X by an

angle 0 about x-axis. Thus the basis vectors of frame Y can be written in terms of basis

vectors of X frame.

y2 = X2Cos(O) +X3sin(0)
(A.1)
y3 = -x2sin(O) +X3cos(O)

This relation can also be expressed in terms of matrices.

y1 1 0 0 x1

Y2 = 0 cos(O) sin(O) x2 (A.2)

y3 0 -sin(O) cos(O) x3

This was a case of simple rotation. The matrix above in square brackets is known as

the rotation matrix from X to Y and is represented as X_Y The rotation matrix can be

generalized for a case when the two reference frames are arbitrarily oriented.






85

(yl,xi) (y, x2) (y1, x3)
XJ = (y2,X1) (y2,X2) (y2,X3) (A.3)

(Y3,X1) (Y3,X2) (Y3,X3)

where (,) means the dot product of the two vectors. Thus,
Y = XY*X (A.4)














APPENDIX B
NUMERICAL TECHNIQUES


B.1 Interpolation of Force Data

This section describes the numerical technique used to obtain the values of coefficients

of lift and drag for cavitator and fins.

B.1.1 Extrapolation Scheme

For a better result, a quadratic interpolation/extrapolation scheme is used. Thus 3 points

would be required to obtain an interpolated or extrapolated data value. Figure B.1 shows

the shape functions used for one dimensional interpolation. Say, points {xi-1, i,xxi+1} are

used to find the value of a function f at point x. The value of f(x) would be given by a

parameter a and the three shape function N1, N2 and N3.

N = 1 + a + or
(xi- x --ii) (Xi- Xi-1)

(xi+ xi-1 )2 (xi+ 1 1)2 a2 (B. 1)
(xi- --xi) (xi -xi+) (i- -xi)_ (xi -,Yi+)
N3 (xi-1 -xi) (x+i-1 -xi+) j

(xi xi+I) (xi xi+I)

where, the value of a shape function can be obtained by finding the value of a.

x xi-1 (B.2)
a= -- (B.2)
xi+1 xi-


aE[0,1] for x [xi-1,xi+1]

a < for x < xil

a > 1 for x > xi+l

Thus a E [0, 1] for interpolation and it is greater than 1 or less than 0 for extrapolation.











N2

N\ / \ 3
0
t 0.5
(D I-
O / \





-0.5 I +
0 0.2 0.4 0.6 0.8 1
Alpha


Figure B. Shape Function for One Dimensional Quadratic Scheme


f(x) = N1 *f(xi-) + N2 f(xi) +N3 *f(xi+l) (B.3)

This method can be extended for 2D and 3D as in case for cavitator and fins respectively.

B.1.2 Cavitator

The coefficients of lift (clc) and drag (cdc) for the cavitator are functions of half angle

(ha) of cavitator cone and angle of attack for cavitator (Xc). The CFD data [5] is avail-

able for combination of points given in Table B.1. Equation B.3 can be extended for 2D

cavitator. 3 3
f(xc, ha) = I N(1, i)N(2,j)f(oc (i),ha(j)) (B.4)
i=1j=1

ac (i) Value of ac at ith node

ha (j) Value of ha at jth node

N(1, i) ith Shape function for ac

N(2,j) jth Shape function for ha

B.1.3 Fins

The coefficients of lift (clfin) and drag (cdfin) for the fins are functions of angle of

attack (af) for fin, immersion (Sf) and sweepback angle (Of). The CFD data is available

for combination of points given in Table B.2. Equation B.3 can be extended for 3D fin.









Table B.1 Grid For Experimental Cavitator Data


Half Angle (ha deg) I {15, 30, 45, 60, 75, 90}
Angle of Attack (ac deg) {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15}


Table B.2 Grid For Experimental Fin Data


Immersion (Sf) {0.1,0.3,0.5,0.7,0.9}
Sweepback (Of deg) {0,15,30,45,60,70}
Angle of Attack (af deg) {0,1,2,3,4,5,6,7,8,9,10,12,15}
1. S>0&S<0.1&0>0
1. S > 0 & S < 0. 1 & 0 > 0
2. S>0.1 & S< 0.3 & > 30
3. S>0.3 &S< 0.5 & > 45
4. S> 0.5 &S< 0.7 & 0 > 60
Data not Available for <0.9
5. S > 0.7 & S < 0.9 & 0 > 70
6. S>0.9&S< 1 &0>0
7. a<0
8. a> 15


f(Sf,, O ), = I I IN(1, i)N(2, j)N(3,k)f(Sf(i), O(j),ayf(k))
i=lj=1 k=l


S(i)

Of(jW

af (k)

N(1, i)

N(2,j)

N(3,j)


(B.5)


Value of Sf at ith node

Value of Of at jth node

Value of af at kth node

ith Shape function for Sf

jth Shape function for f/
kth Shape function for af


B.2 Numerical Linearization

Numerical linearization can be done by the 'linmod' command in the Matlab Simulink

toolbox. This can also be done by noting that, the terms in the A and B matrices are the

derivatives of state rates with respect to states and controls. For example, suppose xo and uo








represent the state and control values at trim. It should be noted that xo is a 9 x 1 (excluding

the positions {x,y, z}) vector and uo is 5 x 1 (cavitator and four fins).


xo = uo wo qo Oo vo Po ro (o To

uo = {c0 e,10 e,20 ,10 ,r20 (B.6)


The equations of motion are of the form as in equation B.7.

x = f(x,u) (B.7)

where the function f is a vector function having 9 outputs and there are 9 states. The code

for nonlinear equation of motion, takes x and u as inputs and give the value of x as output.

Let e define a very small change. Now the element A (i, j) can be calculated as in equation

B.8.
S f (xo + (j), uo)i f(xo, U)i
A(i,j) = 1 < i,j < 9 (B.8)

where, 7(j) means a matrix of size xo with all zeros except jth element, which is equal to

e, and f. represents the ith element of vector f. An element B(i, j) also can be obtained in

a similar way.

f (xo, uo +e (j)) f/(xo, uo)i
B(i,j) = 1 f

where, 7(j) means a matrix of size uo with all zeros except jth element, which is equal to

F.