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

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 FitzCo 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 Openloop Performance . . . . . . . . . . . . . . 47 5.2 ClosedLoop 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 Closedloop 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 Closedloop 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 Closedloop 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 Bodyx 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 OpenLoop 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 Closedloop 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 LoopGain 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 Allloop . . . . . . . 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 Allloop . . . . . . . . . 69 6.16 Eigen v alues for the Perturbed Closedloop 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 ulletlik 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 subsurf ace gun system using Adaptable HighSpeed Undersea Munitions (AHSUM). These are superca vitating Â“kinetickillÂ” 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 closedloop 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. Openloop dynamics are sho wn. The closedloop 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 closedloop 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 closedloop 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 halfangle. 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 banktoturn strate gy for making a turn, and cannot use traditional missile strate gies such as skidtoturn. 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 multibody 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 earthx 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 earthx ed 10
PAGE 22
11 Âˆ e 1 Âˆ e 3Âˆ e2OÂˆ b 1 Âˆ b 2 Âˆ b 3Figure 3.1 Bodyx 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 righthanded 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 bodyx 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 nx 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 nx 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 bodyx ed frame.
PAGE 29
18 3.1.4 Angle of Attack and Sideslip The bodyx 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 bodyx 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 nx 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 bodyx ed frame. I w B is angular v elocity of the body (or bodyx 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 bodyx ed frame, the y w ould be equi v alent to the v elocity components in bodyx 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 bodyx 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 halfangle, g 1 2 of the ca vitator Halfangle, 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 crosssectional 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 bodyx 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 bodyx 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 nx ed frame using Equation 3.14 and 2. Rotate the abo v e obtained nx ed frame to obtain the bodyx ed frame using Equations 3.3 to 3.11 Thus the forces on the ns in bodyx 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 bodyx 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 nx 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 ree 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 bodyx ed frame so as to equate the lefthand and righthand 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 righthand 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 rstorder dif ferential equations. xfxu(4.1) using f : n n as a nonlinear function of a timev 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 righthand side of Equation 4.24 ha v e been calculated. All the terms on righthand side of the Equation 4.20 can be calculated in a similar manner It is tedious to calculate the 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 rstorder 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 timev 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 lateraldirectional 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 lateraldirectional 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 innitynorm 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 straightforw 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 Openloop 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 openloop dynamics are simulated at this condition without an y feedback. Figure 5.1 sho ws the Simulink model used for openloop simulation. The control surf ace deections are x ed at their trim v alues for this simulation. The closedloop 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 openloop response for the torpedo at this trim condition. It can be seen that the openloop system is unstable. The simulation is carried out at the trim that is 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 nonzero 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 (ms1) 0 20 40 60 80 100 0.04 0.03 0.02 0.01 0 0.01 time(s)p (rad s1) 0 20 40 60 80 100 0 0.5 1 1.5 2 2.5 3 3.5 4 time(s)q (rad s1) Figure 5.2 OpenLoop 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 ClosedLoop 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 closedloop 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 closedloop system should ha v e a gain mar gin of at least 6 dB.The closedloop 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 statefeedback 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 steadystate 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 steadystate error e which is made zero using the re gulator Figure 6.1 sho ws the block diagram for this closedloop system. It is required that the closedloop system has an inte grator some where so as to mak e the steadystate 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 Closedloop Model 6.3.1 Model Figure 6.3 sho ws the eigen v alues for the closedloop 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 Closedloop System 6.3.2 Linear Simulations The response of the closedloop 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 steadystate 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 s1)Time(s) command Achieved 0 5 10 15 72.5 73 73.5 74 74.5 75 75.5 time(s)u (m s1) 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 s1) 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 steadystate 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 nonzero 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 s1) 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 s1) 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 s1) 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 allloop gains respecti v ely + 1 3 2 + K x r a k 1r 1 yx xAxBu Figure 6.9 Breakpoints for Calculating the LoopGain 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 closedloop 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 closedloop 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 reducedorder 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 Closedloop 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 Allloop 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 closedloop 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 closedloop 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 Allloop 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.36e2 Q v p r 20 20 12.6e3 2.6e3 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 Closedloop 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 s1)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 s1) No Uncertainty +20% in clf20% 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 s1) No Uncertainty +20% in clf20% 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 s1) No Uncertainty +20% in clf20% 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 s1)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 s1) No Uncertainty +20% in clf20% 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 s1) No Uncertainty +20% in clf20% 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 closedloop 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 allloop. The last v alue is increased to sho wing an impro v ement for the perturbed system. From the analysis of the perturbed closedloop 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 Closedloop 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 steadystate 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 rearrange 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 s1) 0 5 10 15 20 10 0 10 20 time(s)q (deg s1) 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 s1) 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 s1) 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 steadystate 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 s1) 0 5 10 15 2.5 2 1.5 1 0.5 0 0.5 time(s)q (deg s1) 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 s1) 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 s1) 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 s1) 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 s1) 0 10 20 30 40 50 6 4 2 0 2 4 6 time(s)q (deg s1) 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 closedloop 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 s1) 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 s1) 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 s1) 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 openloop 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 closedloop 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 openloop 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 closedloop 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 xaxis. 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 i1 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 vitationCa 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 CavityRunning Behavior of Missiles Arlington, V A, SEAHA C TR 752, Na v al Sea Systems Command, 1975. [5] N. Fine, Six De gr eeofF 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 parttime basis. His research interests include controls and dynamics and optimization. 91
