T H E D Y N A M I C S A N D C O N T R O L OF A N U N D E R W A T E R T O W E D VEHICLE By David A . Hopkin B. A . Sc. (Mechanical Engineering) Technical University of Nova Scotia A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E D E G R E E OF MASTER OF APPLIED SCIENCE in T H E F A C U L T Y O F G R A D U A T E STUDIES MECHANICAL ENGINEERING We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A April 1989 © David A . Hopkin, 1989 In presenting this degree at the thesis University in partial of freely available for reference copying of department publication by his or and her purposes may representatives. <Jp The University of British Columbia Vancouver, Canada DE-6 (2/88) be It shall not permission. Date requirements for an advanced Library shall make it study. I further agree that permission for extensive of this thesis for financial gain Department of (Mec U.>. the British Columbia, I agree that the this thesis for scholarly or fulfilment of t^f-r«*.cy granted is be by understood the that head of my copying or allowed without my written Abstract A linear, quadratic optimum control strategy is applied to a non-linear dynamic model of an underwater towed vehicle in an effort to minimize the pitching motion of the vehicle during the tracking of a reference input. The dynamic model is a two-dimensional, coupled pitch/heave model. The cable interactions with the vehicle are simplified to steady state, and formed into a multi-variable look-up table used in the non-linear model. The normal force resulting from the body of the vehicle is non-linear and consists of two components, an invisid slender-body theory component, and a separated crossflow component. In addition, the dynamic model includes the non-linear effect of the front airfoil's downwash acting upon the rear airfoils. Aerodynamic testing of a scaled vehicle provides the expressions for the non-linear normal body force and moment. These tests also verify the finite aspect ratio corrections for the airfoils, and the downwash effects of the front airfoils. The linear control strategy is based on linear, quadratic optimum control. Simulation results show that proper selection of the state and input weighting matrices result in minimizing the pitch angle of the vehicle to within the control objectives. In addition, simulations of various observer designs show how the tracking and attitude control varies with the selected measurements. ii T a b l e of Contents Abstract ii List o f F i g u r e s vi Acknowledgments ix Nomenclature x 1 Introduction 1 1.1 Description of the Problem 1 1.2 Current Technology and its Limitations 3 1.2.1 Shipboard Compensation 3 1.2.2 Vehicle Compensation 5 1.3 2 Proposed Approach 6 Vehicle a n d C a b l e M o d e l i n g 12 2.1 Related Publications 12 2.2 Cable Model 14 2.3 Mathematical Modeling of the Vehicle 17 2.3.1 Nonlinear Coupled Pitch/Heave Model 17 2.3.2 Linear Coupled Pitch/Heave Model 23 2.3.3 Linear Uncoupled Roll Model 26 2.3.4 Actuator Dynamics 27 iii 3 4 5 6 W i n d T u n n e l Tests 29 3.1 Overview of Aerodynamic Testing 29 3.2 Wind Tunnel Setup 30 3.3 Test Results 32 3.3.1 Body Normal Force Coefficient 32 3.3.2 Body Pitching Moment Coefficient 33 3.3.3 Airfoil Lift Coefficient 34 3.3.4 Downwash Effects 35 Trajectory Control 40 4.1 Related Publications 40 4.2 Control Objectives and Approach 41 4.3 Linear Control Design 43 4.3.1 Controllability and Observability 43 4.3.2 State Feedback Design 48 4.3.3 Reference Model 51 4.3.4 Observer Design 54 Numerical Simulation 58 5.1 Overview 58 5.2 A C S L Environment 60 5.3 P C - M A T L A B Environment 63 S i m u l a t i o n Results 68 6.1 Overview 68 6.2 State Feedback with Reference Input 69 6.2.1 69 Varying the L Q R Weightings iv 7 6.2.2 Varying Actuator Time Constant 73 6.2.3 Relative Importance of the Principle Non-linear Effects 73 6.3 State Feedback with Initial Conditions 76 6.4 Robustness of the State Feedback Design 79 6.5 Observer Design with Reference Input 83 6.6 Roll Simulation 93 Conclusions 95 7.1 Contribution of the Thesis 95 7.2 Summary of the Simulation Results 98 7.3 Needed Further Research 100 Appendices 102 A L i s t i n g of M A T L A B Files 102 B N u m e r i c Values U s e d i n Simulations 116 C L i s t i n g of C a b l e M o d e l i n g P r o g r a m 117 D D e t a i l e d D e r i v a t i o n of the C o n t r o l Objectives 120 E W i n d T u n n e l Test D a t a 122 F A C S L P r o g r a m Listings 126 Bibliography 134 v List of Figures 1.1 Shipboard Motion Compensation Systems, reproduced from Berteaux[4] 4 1.2 Proposed Towing Geometry . . 7 1.3 Towed Body Motion from Chapman 8 1.4 Cable Profile Used in the Hydrodynamic Modelling 10 1.5 Vehicle Geometry 10 2.6 Steady State Cable Model 15 2.7 Cable Angle at Vehicle 16 2.8 Forces Acting on the Vehicle 18 2.9 Root Locus for Open Loop Response 27 3.10 Wind Tunnel Test Setup 30 3.11 Normal Force Coefficient for the Body 32 3.12 Moment Coefficient for the Body 33 3.13 Lift Coefficients for Front and Rear Airfoils 34 3.14 Lift Coefficient Curves for the Front Foil/Body configuration(solid) and the Front/Rear Foil/Body conflguration(dashed) 36 3.15 Assumed Paths of Wing Tip Shed Vortices 38 4.16 Schematic of Regulator Feedback Law 48 4.17 Regulator Root Locus, Q(5,5) varied from 3 to 32769 52 4.18 Regulator Root Locus, R varies from 1 to 1.04 x 10 52 4.19 Schematic for Process and Reference Input vi 6 53 4.20 Schematic of Observer Design 55 4.21 Schematic of Composite System 56 5.22 Program Structure for S S C A N . C S L 60 5.23 Flow Chart for P C - M A T L A B program S I M U L A T E . M 65 6.24 Linear(solid) and Non-linear(dashed) Results for Simulation 1 70 6.25 Linear(solid) and Non-linear(dashed) Results for Simulation 2, Pitch angle Weighting is 60 72 6.26 Linear(solid) and Non-linear(dashed) Results for Simulation 3, Time Constant Reduced to 0.5 sec 74 6.27 Linear(solid) and Non-linear(dashed) Results for Simulation 4, time constant of 2 sec, body weighting of 60, downwash neglected 75 6.28 Linear(solid) and Non-linear(dashed) Results for Simulation 5, downwash and non-linear body forces neglected 77 6.29 Linear(solid) and Non-linear(dashed) Results for Simulation 5, downwash, non-linear body forces and cable angles neglected 78 6.30 Linear(solid) and Non-linear(dashed) Results for Simulation 7, time constant of 2 sec, body weighting of 1, with initial conditions 80 6.31 Linear(solid) and Non-linear(dashed) Results for Simulation 8, time constant of 2 sec, body weighting of 60, with initial conditions 81 6.32 Non-linear results for Simulation 9 at 2.5 ft/sec(solid) and Simulation 10 at 7.5 ft/sec(dashed), time constant of 2 sec, body weighting of 60 . 82 6.33 Linear Process States(solid) and Observer States(dashed), First Observer Design with vertical position and front foil rotation measurements 84 6.34 Linear Process States(solid) and Observer States(dashed) for Second Observer Design, with vertical position and pitch angle measurements . . . vii 86 6.35 Non-linear Process States(solid) and Observer States(dashed), Second Observer Design with vertical position and pitch angle measurements . 87 6.36 Non-linear Process States(solid) and Linear Process States(dashed), Second Observer Design with vertical position and pitch angle measurements 89 6.37 Linear Process States(solid) and Observer States(dashed), Third Observer Design with both foil rotations, vertical position, and pitch angle measurements 90 6.38 Non-linear Process States(solid) and Observer States(dashed), Third Observer Design with both foil rotations, vertical position, and pitch angle measurements 91 6.39 Linear Process States(solid) and Non-linear Process States(dashed), Third Observer Design with four measurements 92 6.40 Linear Uncoupled Roll Response, time constant of 2 sec, roll angle weighting of 8 94 D.41 Side Scan Sonar Beam Geometry 120 D.42 Target Measurement Error Resulting From Body Rotations 121 viii Acknowledgments I would like to extend my sincerest appreciation to my thesis supervisors, Dr. Ian Gartshore and Dr. Michael Davies, for their guidance and support thoughout this work. It has been a great pleasure working with them. I would also like to thank the staff of the drafting department and the model shop at Defence Research Establishment Pacific, for their fast and proficient work in providing the aerodynamic test model. Also, a very special thanks to my wife and children, for their love and support. ix Nomenclature Af, cross-sectional area of cylindrical body ft 2 A cross-sectional area at distance x from the towpoint ft 2 Af profile area subject to separated flow ft 2 ARi aspect ratio 6, span of airfoil Cm normal force coefficient of body Cu airfoil lift coefficient Cf crossflow drag coefficient of body Cd body drag coefficient C normal cable drag coefficient x c c T ft. D body drag lb. D displaced volume of vehicle ft d cable diameter ft. db diameter of body ft. Imth pitch and roll moments of inertia, b 0 3 including added moments of inertia Li hydrodynamic lift of airfoils lb. Lb hydrodynamic lift of body lb. l distance from towpoint to center of mass ft. m Ii M distance from towpoint to the mounting point of the airfoils ft. total mass of vehicle, including added masses slugs x Mo mass of the vehicle slugs M added mass of the vehicle (body and airfoils) slugs moo lift coefficient slope for airfoil with infinite aspect ratio m- lift coefficient slope for each airfoil P(d) tangential hydrodynamic loading function Q{B) normal hydrodynamic loading function Qo dynamic pressure lb/ft 2 Qi airfoil dynamic pressure lb/ft 2 R normal cable drag force per unit length lb/ft Si airfoil area ft T cable tension lb. stream velocity ft/sec Vi velocity at each airfoil ft/sec W cable weight per unit lenght lb/ft y vertical vehicle position ft. p fluid density slug/ft b f V,V 0 2 3 cable friction factor a vehicle pitch angle radians Si airfoil rotation angle radians 7i induced angle of attack of the airfoils radians 76 induced angle of attack of the body radians e cable angle radians i>i angle of attack of the airfoils radians $b angle of attack of the body radians <t> roll angle radians Chapter 1 Introduction 1.1 D e s c r i p t i o n of the P r o b l e m The concept of towing various devices from a floating vessel is certainly not new, and likely dates back several centuries. However, over the last century, the type of device and its sophistication has changed enormously. Many of these devices, such as un- derwater sonar, sprang from military applications. Others were developed as a result of increased interest in physical and chemical oceanography. A l l of this new devices, however, had at least one thing in common; the tasks they were performing were becoming more and more complex. This increasing demand on performance meant it was no longer sufficient to simply tow the devices behind the ship. For many applications, it was desirable to completely control the tow path and attitude of the towed body. This new requirement has precipitated many studies in the last two decades, all of which specifically or partially address a similar, fundamental problem: how to tow an underwater vehicle from a surface ship in a random sea and maintain complete control over the vehicle's trajectory and attitude. In order to understand this problem further, it is necessary to identify the major obstacles to achieving complete control over a towed vehicle's trajectory and attitude. They can be divided into two general categories, environmental and human related. The principal environmental problems are of two forms. First is the induced motion at the towed vehicle due to surface ship motion. This is a result of the surface ship 1 Chapter 1. Introduction 2 motion propagating along the tow cable to the underwater vehicle, and is a function of the tow cable profile and the vehicle body design. The induced motion principally affects the attitude of the vehicle, although in high sea states vehicle performance may also be compromised due to the random trajectory. The second environmental problem is the danger of vehicle collision with the sea bed due to changes in the sea floor contour. This problem is very prevalent in geophysical applications where the vehicle, often towed on 1000 meters or more of cable, is required to be 50 meters or less from the sea floor. The human related problems are of many forms, and depend on the complete towing setup. However, two sources of problems can be identified, ship helmsmanship and the competence of the oceanographic winch operator. Fluctuation in the tow ship speed will result in fluctuation of the underwater vehicle depth. For certain oceanographic experiments this may be very undesirable, and in cases where the vehicle is only meters from the sea bed, there may be catastrophic results. Unexperienced winch operators, or even a lapse of concentration by a competent operator, could have similar undesirable results. The remainder of this chapter reviews the current technology used in attempting to solve these problems. The limitations of these techniques are identified, and a proposed approach to satisfy better the requirements of the current, sophisticated underwater towed vehicle is outlined. Chapter 1. 1.2 Introduction 3 C u r r e n t Technology a n d its L i m i t a t i o n s 1.2.1 Shipboard Compensation Shipboard compensation is one of two principal methods used for achieving control over a towed vehicle's trajectory and attitude. In the past this has been the predominant method used to provide ship motion compensation, thereby reducing the induced motion at the towed vehicle. However, in certain cases, as discused later, this method can also be used to achieve a certain degree of bottom following. A recent publication by Berteaux [5] contains a very good review of the principal techniques which have been used in trying to achieve shipboard motion compensation. Figure 1.1, reproduced from [5], illustrates three of these techniques. The ram tensioner is a passive device, and works on the principal of maintaining a constant cable tension using a precharged hydraulic cylinder. It typically has a high frequency response due to the low inertia of the hydraulic cylinder. The boom bobber, which works in a similar manner to the ram tensioner, tends to have a lower frequency response due to the added inertia of the boom. However, the boom bobber tends to fatigue the cable less by employing fewer sheaves over which the cable passes. These passive compensation methods have the advantage of low power requirement during operation. However, both techniques have the common disadvantages of occupying a large deck space and, more importantly, of providing only a finite distance over which they can compensate. The controllable winch is an active device. It tends to have a moderate to high frequency response, and achieves its control signals from two basic methods. The first method employes the winch as a constant tensioning device, where the control signal is proportional to the difference between the measured cable tension and a reference cable tension. A second, and more common method, is to generate winch control signals that are proportional to the measured ship motion, or more precisely the vertical Chapter 1. Introduction Figure 1.1: Shipboard Motion Compensation Systems, reproduced from Berteaux[4] 4 Chapter 1. Introduction velocity of the ships' tow point. 5 Saunders [27] outlines this latter method in more detail, and provides some experimental results. The controllable winch generally takes up less deck space, and is limited only by cable length when providing motion compensation. However, active devices tend to require substantial amounts of power for proper operation. Shipboard compensation techniques provide several advantages over the vehicle compensation methods described below in 1.2.2. A l l three of the described techniques provide good motion compensation in their operating range. The control equipment is located entirely onboard the ship, so ensuring a lower replacement cost of the towed vehicle. Onboard equipment also simplifies maintenance and troubleshooting. However, this method of compensation does not provide good bottom following capabilities. Only for a near vertical cable profile is it possible to achieve partial bottom following, and even for this case the nonlinear dynamics of the tow cable limit the method. 1.2.2 Vehicle Compensation Vehicle compensation is the second method used to maintain control over the towed vehicle's trajectory and attitude. This method employs a somewhat more sophisticated towed vehicle, capable of generating a variable lift or depressing force which effectively flies the vehicle up and down. There are several systems of this type commercially available in Canada. The best documented of these is the Batfish, which was developed for the Bedford Institute of Oceanography by Desserault [13] at the Technical University of Nova Scotia. The vehicle, which uses a single set of airfoils and a clever hydraulic system, is capable of flying a preprogrammed flight trajectory. The vehicle is predominantly used for C T D profiling. Chapter 1. Introduction 6 Another of these vehicles, the Brutiv, was developed by the Department of Fisheries to fly over scallop and oyster beds and maintain photographic records of their development or depletion. To achieve bottom following, the vehicle employs a simple feedback signal proportional to the error of the measured height off the bottom and a preset height. These two vehicles illustrate this method's advantages over shipboard compensation. They provide reasonably good bottom following and trajectory following capabilities. The cable profiles for towing these vehicles are usually such that the ship motion is partially damped out. However, neither of these, nor any other vehicle today, is capable of providing attitude control. 1.3 Proposed Approach The previous two sections outline the advantages and disadvantages of using shipboard and towed vehicle compensation. Shipboard compensation is better suited for ship motion isolation applications, while vehicle compensation provides good trajectory following capabilities. Since it is desirable to provide both ship motion isolation and trajectory following, a hybrid of these two compensation methods is proposed. The towing arrangement for this hybrid system is shown in Figure 1.2. It consists of two cable parts,a negatively buoyant cable joining to a depressor weight, and a neutrally buoyant cable joining to an actively controlled vehicle. By assuming the mass of the depressor is an order of magnitude larger than the mass of the neutrally buoyant cable, the resulting sharp discontinuity of cable angles at the depressor is assumed to partially decouple the two parts. In this way, the depressor remains strongly coupled to the vertical ship motion, but the vehicle motion can be considered uncoupled. * The first part of the cable system, if assumed uncoupled from the second part, is a Chapter 1. Introduction 7 Figure 1.2: Proposed Towing Geometry configuration that has been examined and simulated extensively in the past decade, as in [23,11,8,26,2,18,9,12,10]. Figure 1.3, reproduced from Chapman [10], shows how a towed body might be expected to respond to a sinusoidal surface ship motion comprising of both a surge and a heave component. The upper part of the figure shows the general profile of the cable, with the applied ± 2 m of tangential motion at the top of the cable. The lower part of the figure is an exploded view of the body trajectory at the lower end of the cable. The motion predicted is elliptical, with the major axis motion close to 4 meters and the minor axis motion less than 0.5 meters. The second part of the cable system, and interactions of the neutrally buoyant cable with a dynamically controlled vehicle, have not previously been examined. This Figure 1.3: Towed Body Motion from Chapman Chapter 1. Introduction 9 interaction is important in developing a hydrodynamic model of the vehicle, and is modelled here using the following simplifying assumptions: 1. The induced heave and surge at the depressor of the two cable system does not differ from that of a single cable system described in Chapman [10], 2. The magnitude of the induced heave is small compared to the length of the neutrally buoyant tow cable, 3. As a result of hydrodynamic damping the heave motion of the depressor is damped out part way along the neutrally buoyant cable (Bowden cable hypothesis described in Chapman [10]), 4. To facilitate a straightforward linear control design, the induced surge at the vehicle due to the induced surge at the depressor is neglected. These assumptions leave the vehicle travelling at a constant velocity, subject to the cable profiles shown in Figure 1.4. The following work develops a non-linear model of this simplified cable/vehicle system. A steady state model of the cable is developed and combined with a non-linear hydrodynamic model of the towed vehicle. This non-linear model is linearized about an equilibrium towing point. The resulting linear model is then used to develop a linear controller to provide bottom following and attitude control for the vehicle. This linear controller is then used in simulations of both the linear and non-linear hydrodynamic models. The geometry of the towed vehicle used in the development of the non-linear hydrodynamic model is shown in Figure 1.5. This geometry is such that an existing vehicle developed at the Defence Research Establishment Pacific could be easily modified to produce a working model for possible future experiments. The vehicle is also very Chapter 1. Introduction 10 Figure 1.5: Vehicle Geometry Chapter 1. Introduction 11 symmetric, eliminating the effects of coupled inertia terms in the non-linear model development. In addition, wind tunnel tests of a scaled model verified the values used for many of the hydrodynamic coefficients. Chapter 2 Vehicle a n d C a b l e M o d e l i n g 2.1 Related Publications Publications relating to the work of this chapter can be divided into two groups, those relating to the development of a mathematical model of an autonomous underwater vehicle, and those relating to the modelling of towed cable systems. A paper by Shupe and McGeer[28] was recently presented at a Military Robotic Applications Workshop. Part of this paper deals with the mathematical modelling of an autonomous underwater vehicle ( A U V ) . The vehicle employs two pairs of controllable hydrofoils for maintaining depth and pitch/roll stability. A controllable rudder provides directional stability. The mathematical model is linear, and consists of two uncoupled parts, a two degree of freedom longitudinal model and a three degree of freedom lateral/directional model. The effect of surge is neglected. The hydrodynamic lift coefficients for the airfoils are based on steady state airfoil theory. These calculations include a finite aspect ratio correction and a body upwash correction. Although details are not provided, a body moment coefficient is also included in the model. The thesis by Boncal[6] also deals with the modelling and simulation of a selected A U V . The model developed is a complete six degree of freedom non-linear model, which is linearized about a straight line trajectory. The vehicle uses two pairs of controllable hydrofoils for roll/pitch/heave stability and a controllable rudder for sway/yaw stability. The hydrodynamic coefficients and equations for this particular vehicle were developed 12 Chapter 2. Vehicle and Cable Modeling 13 previous to Boncals' work, and are therefore only briefly outlined. However, Appendix A of [6] does include a complete listing of the non-linear equations of motion. The papers on cable modelling can be divided into dynamic and steady state models. The dynamic models, such as presented by Chapman[10] and Burroughs and Benz[7], tend to use discretization to solve the set of nonlinear partial differential equations that describe this continuum problem. However, the method of generating the resulting set of ordinary differential equations varies. Chapman[10] uses Newtonian dynamics to generate a set of recursive equations, while Burroughs and Benz[7] use Lagrangian dynamics to generate a set of equations of the form 1(0)8 = F(0), where 6 is the vector of generalized coordinates, F is vector of the forcing functions, and I is the inertial matrix. Of fundamental importance to all dynamic studies are the hydrodynamic loading coefficients applied to the cable. These coefficients are the result of extensive steady state modelling and experimentation. In this area, Steady-State Theory of Towing Cables by Eames[15] is a well recognized paper. The paper reviews the accepted fundamental assumptions made in steady state analysis, presents the governing equations, and provides a very comprehensive review of hydrodynamic loading functions for both faired and nonfaired tow cables. A paper by Anderson[3] provides a few further modifications to the loading functions proposed by Eames[15]. Calkins[8] proceeds one step further than Eames[15] and Anderson[3], and takes both experimental and theoretical relationships for the loading functions and presents a complete steady state analysis of a high speed marine towed system. Chapter 2. Vehicle and Cable Modeling 2.2 14 Cable M o d e l As mentioned in section 1.3, the interaction between the neutrally buoyant cable and the vehicle is important for developing the dynamic model of the vehicle. The dynamic model requires a relationship for cable angle and tension at the vehicle as a function of other known variables. To obtain this relationship in a form which can be included in both the linear and nonlinear model, it is necessary to make the additional following assumptions: 1. Inertial effects of the cable are neglected, 2. The cable is completely flexible, and supports no internal forces other than tension, 3. The induced hydrodynamic forces due to the cable motion are small compared to the hydrodynamic forces due to the constant stream velocity, and may be neglected. These assumptions simplify the cable/vehicle interaction to a steady state problem. The steady state cable model is developed in a similar manner to Eames[15], with a cable element shown in Figure 2.6. By summing the forces normal and tangential to the cable element ds and using small angle approximations, the following governing equations are derived: dT = P(6)ds Td9 = -Q(8)ds (2.1) where P(0) and Q(0) are the tangential and normal components of the hydrodynamic force acting on the cable element ds. As proposed in [15], these hydrodynamic loading functions can be written as: Q(0) = R{(1-p)sin d 2 +psm6) Chapter 2. Vehicle and Cable Modeling 15 Figure 2.6: Steady State Cable Model P(6) = Rpcos0 (2.2) where R is the drag of a cable element normal to the flow, and p is referred to as the friction ratio, which for a bare cable is assumed to be 0.05. The normal drag force R is expressed as: R = 0.5C pdV 2 r where p is the fluid density, d is the cable diameter, V is the stream velocity, and C is a r normal drag coefficient which is a function of the Reynolds Number and the stranding or roughness of the cable. For this work C is assumed to be 1.2, a typical value for a r stranded wire rope cable with no fairing. The equations in 2.1 are recursive and, given the boundary condition at the vehicle, can be solved along the cable to the depressor. The depressor is then chosen as the origin and, using the array of incremental cable angles dd, the profile of the cable and the position of the vehicle are determined. A Fortran program 2 C A B L E . F O R , listed Chapter 2. Vehicle and Cable Modeling 16 Vehicle Position - Ft Figure 2.7: Cable Angle at Vehicle in Appendix C , solves these equations for input cable parameters and geometry, and plots the resulting cable profile. The boundary condition at the vehicle end of the cable consist of lift and drag components. By determining the steady state cable profiles for various values of lift and drag, it is possible to establish a relationship for the cable angle at the vehicle as a function of the vertical position y and the body drag. This relationship is presented in Figure 2.7. Figure 2.7 is translated to a lookup table and used in the nonlinear model to provide the necessary relationship between cable angle at the vehicle and the known variables, vehicle position and drag. In the linearized model, the drag of the vehicle is assumed to be a constant. The linear model, therefore, simply uses a straightline approximation to one curve of the Figure 2.7 that corresponds to this constant drag value. Chapter 2. Vehicle and Cable Modeling 2.3 17 M a t h e m a t i c a l M o d e l i n g of the Vehicle The following mathematical models are developed for the vehicle shown in Figure 1.5. This vehicle has two pairs of symmetric airfoils located equidistant fore and aft of the center of the vehicle. The body is cylindrical in shape, with a hemispherical nose and a blunt tail section. The mass is assumed to be uniformly distributed along the length of the vehicle. With this geometry, the vehicle is assumed to have three point symmetry about its center of mass, which is located at the midpoint of the vehicle. This symmetry allows the inertial cross product terms to be set to zero, simplifying the mathematical model to an inertially uncoupled system. One additional assumption is made to further simplify the cable/vehicle interaction. It is assumed that the maximum excursion of the vehicle from its equilibrium towing position is small compared to the length of the neutrally buoyant tow cable. This allows the induced surge due to this motion to be neglected. For models of this type, the longitudinal (pitch/heave) motions are uncoupled from the lateral (yaw/sway/roll) motions. Therefore, the following three independent models are developed: a nonlinear coupled pitch/heave model, a linear coupled pitch/heave model, and a linear coupled yaw/sway/roll model which is reduced to an uncoupled roll model. The coordinate system for the coupled pitch/heave motion is shown in Figure 2.8. 2.3.1 Nonlinear Coupled Pitch/Heave M o d e l As previously mentioned, the body symmetry reduces the modeling problem to that of modeling an inertially uncoupled system. Newtonian dynamics are used to separately solve for y and <ii, with Figure 2.8 showing the forces acting upon the vehicle. In the Figure, L\ and Li are the hydrodynamic lifts from each pair of foils. L is the normal 0 Chapter 2. Vehicle and Cable Modeling 18 x Figure 2.8: Forces Acting on the Vehicle force generated by the body of the vehicle, and Dy, is the associated drag of the body. T is the magnitude of the cable tension at the vehicle, acting at an angle 9. Summing forces in the y direction, the simplest form of the y equation is written as: 53 F — My = Li cos 7x + L cos j ~ T sin 9 + L cos a 2 2 b (2.3) T+ H y d r o f o i l Lift Forces In equation 2.3, the hydrodynamic lift forces L\ and are derived from linear, steady- state airfoil theory. For this case, Li = C qiSi Li t = l,2 (2.4) where Cu is the lift coefficient of the airfoils, g,- is the dynamic pressure for the airfoil, and Si is the airfoil area. The lift coefficient Cn is a function of the angle of attack Chapter 2. Vehicle and Cable Modeling ipi, and can be written as Cn = 19 m,-0», where m,- is the slope of the lift and ipi is the angle of attack. For a finite aspect ratio airfoil, m, = ——- (2.5) where the aspect ratio ARi is defined as bf/Si, where 6, and Si are the span and the area of foil respectively. The correction factor accounts for the reduction of lift as ARi is reduced from its usual theoretical value of infinity to a smaller, finite value. Kuethe/Chow[19] present a more complete description of this effect and its derivation. For the numeric simulation, ARi is based on 6,- equal to the total span of the pair of foils less the body diameter, and is verified with wind tunnel tests. The dynamic pressure for the airfoil can be written as pV 2 ft = ^ (2-6) where Vi = ^V where l m 2 0 + (y + (l -l )ay m (2.7) i is the distance from the tow point to the center of mass, and I,- is the distance from the tow point to the mounting point of the airfoil. Vo is the stream velocity. The angle of attack ^ is expressed as i> = a + 6i - 7i (2.8) { The induced angle of attack 7; is comprised of two components, so that r .: = arctan ( ± + (2.9) With equations 2.9, 2.8, 2.7 and 2.6, equation 2.4 becomes 'pSi(V + (y + 2 Li = m-i a + 8i- arctan f ^ - + ^ m (l -li)aY) m (2.10) Chapter 2. Vehicle and Cable Modeling 20 B o d y Forces Lb a n d Db Lb is the force generated by the body. This force acts normal to the body, and is modeled using a method similar to that presented by Nielson[25] and Atraghji[4], and is usually referred to as Allen and Perkin's crossflow theory. Lb is assumed to have two components. The first component, L p represents the calculated lift using an invisid slender-body theory, and results from the integration of < -»> 2 where qo is the dynamic pressure, ipb is the angle of attack of the body, and A x is the cross-sectional area at a point x from the front of the vehicle. Integration of 2.11 over the length of the vehicle and substitution of if)b = & ~ 76 yields L = 2.0 q Ab(a-jb) p (2.12) 0 where 75 = arctan(y/Po) is the induced angle of attack for the body. The only contributing part of the body to L is the nose section, where dAb/dx ^ 0. Note that if p the flow is assumed to be unseparated at the tail by the use of an aerodynamic tail section, J A dx = 0 and no net lift is generated, and we have D'Alembert's paradox L 0 x for invisid flow. The second component of Lb, L , CJ represents the normal force due to separated crossflow and is written as Lf = c C q il>lA f ef 0 (2.13) c where C f is the crossflow drag coefficient, and A f is the profile area of the body c c subject to the separated crossflow. The profile area A f is a function of body angle of c attack ifrb and the body and nose geometry. For this reason, experiments are required to accurately represent this component of Lb. Section 3 outlines the wind tunnel testing performed for this particular vehicle. Chapter 2. Vehicle and Cable Modeling 21 Combining equations 2.12 and 2.13 yields L = 2q A b 0 b (a - arctan + C q ^A cf 0 cf (2.14) The body force Db is assumed to be a constant and results from the separated flow at the blunt tail, such that Db = 0.5C pA V d b (2.15) 2 Q where Ab is the cross-sectional area of the cylindrical body, and Cd is the drag coefficient. Any additional induced body drag has been neglected. C a b l e Forces In equation 2.3, T represents the cable tension at the vehicle acting at an angle 9. From simple trigonometric relations, T can be written as T = — cos 9 (2.16) a where T is the net horizontal force. Since surge motion is assumed to be zero in this x study, T is simply the sum of all the horizontal components of the forces shown in x Figure 2.8, or T x = Db cos a + Li sin 7; + Lb sin a (2-17) The cable angle 9 is determined from the lookup table outlined in section 2.2. Thus, T _ Db cos a + Li sin 7,- + L sin a cos 9 h Added Mass In equation 2.3, M includes the added mass of the body and the two pairs of airfoils. To determine the added mass of the body, Mb, Newman's[24] added mass coefficients for spheroids are used. This method yields a value of Mb = 0.8poDo, where DQ is the Chapter 2. Vehicle and Cable Modeling 22 volume of the body. Since the vehicle is neutrally buoyant, this means M = O.8M0, 0 where Mo is the vehicle mass. This value of 0.8 is slightly less than the commonly used value of 1.0 for cylinder since, for a body of this length to diameter ratio, the end effects are significant in determining Mb. The added mass of the airfoils is assumed to be equal to the mass of fluid contained in a circular cylinder of volume 7r6j 5,/4. Final F o r m of y Equation By substituting equations 2.14, 2.18, and 2.10, and the relationships for 7; and 7& into equation 2.3 we obtain VA •(V l rrii r 2 0 -^•(y- (/ -/ •)A) ) 2 ^ m « + « - — ( i cos (arctan ^ + , X l + ^ ) ) x < ^ ^ j ) 0.6C pAVf+ (2.19) d « + «,-arctan U - sin I arctan I — + V (2qoA \Vo + < i ^ L X —-'^ 1 ) Vo + tan0+ )) (a - arctan ( ^ r ) ) + Cc/SoV'b^c/^ sin a 0 (2q A 0 (a - arctan b + Cc/^oV'b^c/^ cos a Final F o r m of a Equation The simplest form of the ci equation is M = I a = Li(l m m — li) cos 7J — Tl sin 9 + L l cos a m 0 D (2.20) Chapter 2. Vehicle and Cable Modeling where (Z — li), lb, and l are the distances from the center of mass to the origin of m m 23 the forces Li, Lb, and T respectively. The measures (l — /,•) and l m m are constant for a given body geometry. However, lb is a function of the body angle of attack if>b, and requires wind tunnel tests to be accurately represented. Since the form of 2.20 is very similar to 2.3, the final form of the a equation is simply fpAi(V + (y + 2 0 m,(l - U) (l -li)a) y 2 m x a + 6 -arctan ( ^ + m t cos I arctan I — + V 0 0.5C pA V -r 2 d b o (2.21) a =— < a + 6i - arctan U + x y , (l sin I arctan I ——h m V (2q A 0 b \Vo where I m 2.3.2 + > l tan 0+ m Vo (a - arctan (^r)) + Cc/^oV'b^c/j sin a l (2qoAb (a - arctan b U)a ^^-^j + Cc/go^Mc/^ c o s a includes the added moment of inertia of the body and the airfoils. Linear Coupled Pitch/Heave M o d e l The linear coupled pitch/heave model is the second of three models developed, and is the fundamental building block for the control work of Chapter 3.0. For this reason it is desirable to obtain the linear model in the following form x = Ax + Bit (2.22) Chapter 2. Vehicle and Cable Modeling 24 Here, x is referred to as the state vector and u the input vector, while A and B are suitable matrices. The state vector for this linear model is y y x= (2.23) a a As shown in the previous section, these four variables completely describe the coupled pitch/heave motion. The input vector u is simply (2.24) u = 62 the rotation angles of the airfoils. With these vectors defined, A must be a 4 x 4 matrix and B must be a 4 x 2 matrix. In order to define these two matrices, the nonlinear model of section 2.3.1 is linearized. Linearization of the model is achieved by retaining only the first order terms of a Taylor series expansion about an equilibrium state x for equations 2.3 and 2.20. This yields My = (Li cos 7; — T sin 6 + Lb cos ct) | d (Li cos 7,dx Ia m = (Li(l m — T sin 6 + Lb cos a) | s=0 - k) cos7,- - Tl m d -QZ (Li(l m £ = 0 + (x — x) (2.25) sin6 + L h cos a)\ =o + b - li) cos 7, — Tl m x sin 6 + Lbh cos a) \x=o(x - x) (2.26) Performing the above partial differentiation and evaluating the terms at x = 0 yields Chapter 2. Vehicle and Cable Modeling 25 the following form of A and B 1 0 0 0 0 Cyy Cyy Cya Cya 0 0 0 0 Cyl B = 1 0 Cmy Cmy CmaCma (2.27) 0 Cl C2 m m where Cyy —20.0M -go Cyy — C y a — i f ( (2.28) 5 , W i 9 - 0 -Simiq (l 0 Cya — (2.29) (5,-m,- + 2.0 A) ^ + 2 - 0 q o A (2.30) , b - it) m (2.31) VM 0 Cyi — (2.32) — Dblm 20.0/ -go (Sirmilm - U) + 2.0A l ) V Im 1 J- ySiTUiilm - h)qo ^ + 2.0q A lb Cmy — (2.33) m Cmy (2.34) b b 0 Ca — m Ca — m 0 - 5 , m , g ( Z - Z,) 0 (2.35) b 2 m (2.36) V Im 0 Ci m = Si-miq (lm - U) (2.37) Q As a result of this linearization, the following nonlinear effects are lost: 1. The horizontal component of the cable tension T is now a constant, equal to x the body drag Db. As mentioned in section 2.2, this means the cable angle is proportional to the y excursion, as shown in equation 2.28 for C . yy 2. The foil velocity is now a constant. The induced velocity terms involve the square of the states, and for a linear control design, must be neglected. Chapter 2. Vehicle and Cable Modeling 26 3. The crossflow component of the body force is now neglected, leaving the body force equal to the potential flow contribution L . p The importance of the loss of the above effects is discussed further in Chapter 5.0. With the above linear model, it is now possible to examine the stability of the open loop system. Appendix B details the numeric values selected for the model coefficients denned by equations 2.28 through to 2.37. For these values A = 0 1 0 0 -0.0256 -4.322 21.99 0 0 0 0 1 -0.0305 -0.8963 4.51 -4.78 (2.38) for which the eigenvalues are [ -6.7 +0.18» -0.18» -2.41 ] These eigenvalues show the open loop system is slightly oscillatory at a stream velocity of 1.5 m/sec. Varying the stream velocity from 0.15 to 2.0 m/sec produces the root locus plot in Figure 2.9. This shows the pair of oscillatory poles move outwards along the imaginary axis, with the other two poles moving outwards along the negative real axis. 2.3.3 Linear Uncoupled R o l l M o d e l The linear uncoupled roll model is the last of the three mathematical models. This model is developed to only provide insight into the time constants and magnitudes of airfoil defections associated with roll correction and stability. For this reason, a linear model is adequate. The governing equation is simply Ii4> ~ C i4 + miln - m ip = 0 t 2 2 (2.39) Chapter 2. Vehicle and Cable Modeling -5 - 27 4 - 3 - 2 - 1 0 1 Real Part Figure 2.9: Root Locus for Open Loop Response where C^ represents the resistance to roll due to the airfoils, and tpi is the angle of attack of the foils. In this particular case, ipi = the foil rotation. This equation, when cast into the conventional linear, time-invariant form appears as X = 0 0 1 x + 0 0 u mi — r a (2.40) 2 with the state vector x and the input vector u as x = u Si S 2.3.4 (2.41) 2 Actuator Dynamics The dynamics of the actuator are modeled as a first order system with a given time constant r . While the selection of r is at present completely arbitrary, it does provide the flexibility to simulate the vehicles performance for various actuator time constants. This flexibility is very important, since the value of r is often dictated by criteria not Chapter 2. Vehicle and Cable Modeling 28 related to the vehicle performance, such as electrical power or physical space limitations. By selecting large time constants it may be possible to maintain selected performance criteria while minimizing sudden motions of the vehicle. Expressed in matrix form, the actuator model is 1 X = 0 1 0 T 0 X + 0 T . i y= u 1 (2.42) 0 r 0 1 x T J where the state vector x is 1/r times the actual foil rotation, and the input vector u is the input signal to the actuator. For the linear simulations, the linear vehicle model is augmented with this actuator model to yield a sixth order system. Chapter 3 W i n d T u n n e l Tests 3.1 Overview of Aerodynamic Testing The accuracy of the three mathematical models presented in the previous section is partly dependent upon the relationships developed for the hydrodynamic coefficients. While many of the coefficients can be calculated using basic theoretical equations, some coefficients, such as the crossflow component C f of the body lift, cannot. In adc dition, even the basic theoretical equations become less representative in the presence of complex body/airfoil interactions. Many of these body/airfoil and airfoil/airfoil interactions are particular to the vehicle geometry, and can only be determined accurately with proper aerodynamic or hydrodynamic testing. Aerodynamic testing is selected over hydrodynamic testing because of the existing and readily available aerodynamic test facility. The tests are performed in the large boundary layer wind tunnel in the Mechanical Engineering Department's Aerodynamic Laboratory. The test section of the tunnel is about 1.5 x 2.0 meters, with a maximum wind speed of approximately 25 m/sec. This large test section means a 62% scaled model results in 0.53% blockage at zero angle of attack and 2.6% blockage at a 30 degree angle of attack. For these low values of blockage, blockage corrections can be neglected. The maximum speed of 25 m/sec provides a Reynolds Number of 3£ = 2.0 X 10 . s This Reynolds Number, which is based on the body diameter, is less than the desired 29 Chapter 3. Wind Tunnel Tests 30 Figure 3.10: Wind Tunnel Test Setup value of 9ft = 3.0 x 10 . Therefore, tests are done at two values of 9ft, 1.0 x 10 and 5 5 1.53 x 10 , to show the dependance of the measured data on 9ft. 5 The scope of the testing is limited to measuring static aerodynamic coefficients. The normal force coefficient of the body and the lift coefficient of the airfoils are determined, and data is taken to examine the airfoil downwash relationship. No provisions are made to examine any rate coefficients. 3.2 W i n d T u n n e l Setup The setup for the wind tunnel testing is shown in Figure 3.10. The body of the vehicle is made from 5 inch, Schedule 80 P V C pipe, with a machined P V C spherical nose piece. The airfoils are N A C A 0015 sections[l], cut from P V C plate and mounted at their quarter cord location. The airfoils are adjustable in 5 degree increments. A machined aluminum insert piece is lightly press fit into the tail of the vehicle. The Chapter 3. Wind Tunnel Tests 31 entire assembly is then lightly press fit onto the tapered end of the sting balance, resulting in the mounted configuration shown in Figure 3.10. The balance is a specially designed apparatus built by Aerolab for the Mechanical Engineering Department. The sting type of mounting arrangement allows the body to pitch and yaw and still maintain the body essentially in the middle of the tunnel. The balance provides five signals which can be used to obtain the pitch, yaw and drag forces, as well as pitch and yaw moments. A typical test session is as follows: 1. The balance is assembled and the model is mounted. If the airfoils are mounted they are set to their initial angles. The electronics are then zeroed, and the wind tunnel is turned on and brought up to speed. 2. Using the balance, the body is then rotated in 5 degree increments, up and down through a set of values which always keep the airfoils in their linear operating range. This results in two sets of data points for each angle combination. 3. The airfoil angle is then adjusted, and the process is repeated. Periodically, the velocity is measured to verify it is remaining constant. 4. The velocity is then changed, and the above steps 2, and 3 are repeated. A total of over 750 data points are recorded for this particular vehicle. The following foil/body combinations are set up and tested in the tunnel: 1. the body by itself with no airfoils, 2. the body with only the front foils, 3. the body with only the rear foils, 4. the body with both sets of airfoils. Chapter 3. Wind Tunnel Tests 32 5 4 3 .8 o JS 8 u | 13 5 2 1 0 -1 + - V = 12m/sec o - V = 18 m/sec x - Atraghji -2 -3 -4f -30 -20 -10 0 10 20 30 Body Rotation - Degrees Figure 3.11: Normal Force Coefficient for the Body 3.3 3.3.1 Test Results B o d y N o r m a l Force Coefficient The normal force coefficient is defined as Cm = qoA (3.43) b Figure 3.11 shows the empirical normal force coefficient Cm for the body with no foils attached. The straight solid line represents the invisid slender body theory. The deviation from this theory is the separated crossflow component. The '+' and the 'o' represent data from two velocities, with each marked point representing the average of at least two data points. The data points marked with an 'x' are from Atraghji[4], and show good agreement. The data from Atraghji is for a Reynolds Number of 1.4 x 10 6 based on maximum body diameter, and total body length to diameter ratio of 17 (Atraghji test configuration 1001). A second order least square polynomial fit of this Chapter 3. Wind Tunnel Tests 33 30 25 20 15 o a 10 I o U + V =; 12rn/sec o - V = 18 m/sec -5 -10 -15 -30 -20 -10 0 10 20 30 Body Rotation - Degrees Figure 3.12: Moment Coefficient for the Body data is used by the non-linear model. The polynomial coefficients are presented in Appendix B . 3.3.2 B o d y P i t c h i n g M o m e n t Coefficient For this work, the pitching moment coefficient is assumed to be about the center of mass . Figure 3.12 shows the measured moment coefficient CMI resulting from the body normal force, where the moment coefficient is defined as CtAb = M qoAbdb df, = body diameter Again, the two sets of data represent two stream velocities. (3.44) The tendency for Cj^b to approach zero for higher angles of attack is a result of the crossflow component dominating over the invisid component. At higher angles of attack, the point on the body for which the flow starts to separate moves from the tail of the body forward Chapter 3. Wind Tunnel Tests 34 •G B o o U -5 Foil Rotation - Degrees 0 5 10 Foil Rotation - Degrees Figure 3.13: Lift Coefficients for Front and Rear Airfoils towards the nose. As more of the body produces separated flow, the net crossflow force moves towards the middle of the body (see [25] pg. 90), which in this case is the center of mass. Therefore, at the higher angles of attack the moment coefficient tends towards zero. The non-linear model employs a least square polynomial fit to the data shown in Figure 3.12. These polynomial expressions are included in Appendix B . 3.3.3 A i r f o i l Lift Coefficient Figure 3.13 shows the measured lift coefficients for the front and rear airfoils. Each set of data was recorded in the absence of the other airfoil, but in the presence of the body. The solid lines represent linear, least square fits to the data, adjusted to pass through the origin. The slope of the lift coefficient curve for the front foil is 0.087 per degree, and the slope of the rear foil lift coefficient is 0.075 per degree. The calculated Chapter 3. Wind Tunnel Tests theoretical value for this slope, with an aspect ratio of 4, is 0.073 per degree. 35 The slightly larger slope of the front foil lift coefficient curve is expected to result from the front foil operating in the flow field of the nose section. This flow field creates a velocity greater than the stream velocity over part of the front foils, resulting in the increased lift. However, the rear foils appear to have a more uniform velocity of Vo over their span. 3.3.4 D o w n w a s h Effects In order to evaluate the downwash effects of the front foils on the rear foils, the body is tested with both sets of airfoils mounted. During the testing, the front foil and the body are rotated through a sequence of angles such that the angle of attack of the front foil does not exceed 10 degrees. The rear foil is rotated to maintain an angle of attack of 0 degrees. By comparing the lift coefficient of this configuration with the lift coefficient of the front foil only configuration, it is possible to determine the lift coefficient of the rear foils in the presence of the downwash of the front foils. Figure 3.14 shows the lift coefficients of these two configurations for various body angles. The solid lines in each plot represent the lift coefficient curve for the front foils only configuration. The dashed lines represent the results for the front and rear foils combination. A linear, least squares fit for these two sets of curves gives the lift coefficient slope for the front foils only, and a lift coefficient slope for the front and rear foil combination. The difference between these two slopes gives the lift coefficient slope of the rear foil in the presence of the downwash of the front foil for various body angles. These slopes are shown in Table 3.1. Note that as the body angle increases, the effect of the downwash diminishes. Also, since the effect of the downwash is a function of both the angle of attack of the front foils and the body, this effect can only be included in the non-linear model. Chapter 3. Wind Tunnel Tests Degrees 36 Degrees Figure 3.14: Lift Coefficient Curves for the Front Foil/Body configuration(solid) and the Front/Rear Foil/Body configuration(dashed) Chapter 3. Wind Tunnel Tests Body Angle (Degrees) 0 5 10 15 20 25 30 37 Measured Lift Coefficient Slope for the Rear Foils -0.0380 -0.0225 -0.0133 -0.0073 Theoretical Lift Coefficient Slope for the Bear Foils -0.0492 -0.0209 -0.0134 -0.0080 -0.0058 -0.0038 -0.0026 - Table 3.1: Measured and Theoretical Lift Coefficient Slopes for Rear Foils resulting from the Downwash of the Front Foils Table 3.1 also shows the theoretical lift coefficient slopes of the rear foil in the presence of the front foil downwash. These calculated slopes are based on the wing-tail interference material presented by Nielson[25]. In this work, Nielson shows how the lift on a tail section in the presence of a fixed vortex can be expressed as (L )v 2 = *T L T 2 (3.45) V l 72 0 r where T/Vol is the nondimensional vortex strength, £ 2 / 7 2 is the lift coefficient slope r of the rear wing, and ix is the tail interference factor. The reference length l is based r on the rear wing dimensions and is chosen as l = 27T&2. The vortex strength T can be r expressed as = J " ?? ., 0 1 (3.46) 1 A(y - 0.5c4) R v where yv is the vertical vortex position as shown in Figure 3.15. Making these substitutions, and recalling AR+ = b /Si, equation 3.45 can be written in lift coefficient form 2 as (m )v 2 = IT-—7^-7 IvAR^yy _ , v - 0.54) n Implicit in the development of equation 3.47 are the following assumptions: (3.47) Chapter 3. Wind Tunnel Tests 38 y ^ — r yv Figure 3.15: Assumed Paths of Wing T i p Shed Vortices 1. A single vortex model is sufficiently accurate to predict the rear wing loading, 2. The single vortex is shed at the wing tip of each front foil, as shown in Figure 3.15, 3. The vortex strength is time independent, and moves in the stream direction, as shown in Figure 3.15, 4. The effects of body shed vortices are neglected. From the assumed vortex geometry, it is possible to determine values of ix (see [25] pg. 193) for particular body geometries, and using equation 3.47 determine the various theoretical lift coefficient slopes shown in Table 3.1. The measured results agree very Chapter 3. Wind Tunnel Tests 39 well with the theory. The large error at a body angle of 0 degrees results from the shed vortex passing very close to the rear airfoil wing tip. A singularity exists at the wing tip, and it is therefore impossible to select a value for ix accurately. For the non-linear model, a lookup table is used to describe the downwash lift coefficient as a function of the front foil and body angles of attack. Chapter 4 Trajectory Control 4.1 Related Publications As a result of the increased complexity of underwater vehicles and of the tasks that they perform, modern control strategies are starting to be examined for use in these vehicles. Because most of these vehicles have several control variables, they are classified as multi-input multi-output, or M I M O systems. Such systems employ one of a variety of control strategies, depending on the particular requirements of the vehicle and the desired complexity of the system. One of the simplest strategies for controlling M I M O systems, known as decoupling, assigns each control variable to its own task. This strategy is presented in a paper by Shupe and McGeer[28]. By separating the tasks of the fore and aft pairs of control surfaces to heave and pitch control respectively, the state feedback design is reduced to a simple pole placement problem. The resulting controller is very straightforward to implement. However, this technique does not take into account the possibility of interaction between the various loops. A M I M O system allows the designer not only to select the closed-loop pole locations as in Shupe[28], but to also optimize the performance of the system in some manner. The particular type of optimization depends upon vehicle design requirements. Two recent theses by Harris[17] and Dreher[14], and a publication by Martin[20], present control strategies based on the Linear Quadratic Gaussian (LQG) with Loop Transfer 40 Chapter 4. Trajectory Control 41 Recovery (LTR) design methodology. This technique shapes the vehicle response in the frequency domain, based on low frequency design criteria and a particular crossover frequency. The crossover frequency is chosen to minimize the effects of high frequency modeling errors, sensor noise, and, in some cases, surface wave excitation. A third control strategy for M I M O systems is linear quadratic optimal control. A masters thesis by Boncal[6] presents the use of this control strategy to study model based maneuvering controls for autonomous underwater vehicles. This method optimizes vehicle performance by directly weighting the importance of the errors in individual states and inputs of the model. This strategy is particularly well suited for tracking applications, as Boncal shows by using a reference model in conjunction with the optimal control to achieve tracking of multiple inputs. 4.2 C o n t r o l Objectives a n d A p p r o a c h In order to outline the specific control objectives, a brief overview of the use of this particular vehicle is necessary. The vehicle and towing configuration, while general in nature, is specially suited for side scan sonar applications. Side scan sonar is commercially used for providing geo-physical surveys of the sea floor. The sonar transducers used for this work tend to have very narrow, fan shaped beam patterns. The transducers are located on the vehicle so that the beams look down and to each side of the vehicle, and at any point in time survey a narrow strip of the sea floor oriented normal to the attitude of the vehicle. A constant attitude of the vehicle produces consecutive, parallel strip surveys which, when combined, produce a coherent survey image. In addition, the height of the vehicle above the sea floor affects the quality of the survey record. A n optimal height is in the order of 10 meters. Appendix D provides a detailed derivation of the maximum pitch angle and pitch Chapter 4. Trajectory Control 42 rate allowed for side scan sonar operation. The maximum pitch angle ensures that the size of detected objects is within a chosen tolerance, and the maximum pitch rate ensures there are no gaps in the survey record. Based on this criteria, the control objectives can be summarized as follows: 1. A trajectory-following capability is required to maintain the vehicle at a constant height above the sea floor. A n estimated maximum climb angle of 11 degrees requires a maximum average vertical velocity of 1 ft /sec. 2. During bottom following, the induced pitching motion should be minimized. A n estimated maximum pitch angle is 6 degrees for conventional sonar, 2 degrees for specialized applications. 3. During bottom following, the pitch rate must be minimized. Estimated maximum pitch rates are 20 degrees/sec for conventional sonar, 4 degrees/sec for specialized applications. The linearized model presented in section 2.3.2 is of the form x = Ax + y = Cx -f D M (4.48) This is the augmented system, where x= Si S' y y a ct 2 The first two states are r times the actual foil rotations, and result from the augmentation of the actuator dynamics. A is a 6 x 6 matrix, and B is a 6 x 2 matrix. Since y is the measured output, C varies according to the measured states, and D is zero. This forms a classic multi-input multi-output system. As previously mentioned, this Chapter 4. Trajectory Control 43 allows the closed loop dynamics of the system to be selected, and at the same time optimization of the performance of the system with respect to a chosen criterion. Linear, quadratic optimum control is chosen to provide this system optimization. In using an optimum control strategy, a feedback law of the form u = Gx is used. The gain matrix G is chosen to minimize the quadratic performance integral (4.49) to where Q and R are symmetric positive definite weighting matrices of the state and the input respectively. By selecting appropriate values of Q and R, the desired control objectives are achieved while maintaining the inputs below saturation levels. The composite control system is developed in four stages. Section 3.3.1 examines the controllability and observability of the open loop system. Section 3.3.2 presents the derivation of the optimum control strategy that determines the gain matrix G . A reference input is then introduced in section 3.3.3 to provide the tracking capability. Finally, in section 3.3.4, an observer design is presented. The observer design is needed for situations in which not all of the states are available to be used for the feedback law u = Gx. The reference input is also included in the observer design. 4.3 4.3.1 Linear Control Design Controllability and Observability The concepts of controllability and observability were first introduced by Kalman in the late 1950's. Controllability determines the ability of the input(s) to directly or indirectly effect all of the states. Observability determines the ability to estimate all of the states from examination of the output y. Friedland[16] provides formal definitions of these concepts. Chapter 4. Trajectory Control 44 A n uncontrollable system is therefore one in which at least one state can not be affected by the input. If this uncontrollable subsystem is stable, the system is referred to as stabilizable. However, if the subsystem is unstable, so is the system, and a control system is of no use. It is therefore necessary to examine the controllability of a system before pursuing the development of a control strategy. The condition of observability is important for the development of an observer. As mentioned, the purpose of the observer design is to provide an estimate of all the states from the output y. By definition, if the system is not observable, then the design of an observer is not possible. There are several means for determining the controllability of a system. The following theorem is usually referred to as the controllability theorem, and applies to the more general time variant system. C o n t r o l l a b i l i t y T h e o r e m 1 A system is controllable if and only if the matrix P(T, t) = £ $(T, A)B(A)B'(A)$'(T, t)d\ (4.50) is nonsingular for some T > t, where $(T, t) is the state-transition matrix of the system. P(T, t) is the controllability Grammian, and for time-invariant systems can be simplified to P(T) = [ T Jo e BB'e dt At (4.51) A,t A second, and somewhat simpler, method commonly used for the time-invariant case is to show that the rank of the controllability matrix Q = [ B AB • • • A ~*B ] k (4.52) is equal to k, the order of the system. If the rank of Q is less than k, the system is not controllable. Chapter 4. Trajectory Control 45 To determine the observability of a system there exists a similar theorem to that of the controllability theorem. For the time-invariant case, this observability theorem can be simplified to M(T) = f e ' C'Ce dT Jo T A T (4.53) Ar where M(T) is referred to as the observability Grammian. In this case, if the singular values of M ( T ) are non-zero the system is observable, and the closeness of the singular values to zero indicate the systems degree of observability. A second method also exists for determining the observability of time-invariant cases. A system is observable if the rank of the observability matrix C JV = CA (4.54) OA*- 1 is equal to k, the order of the system. If the rank of N is less than k, the system is not observable. To determine the controllability and observability of the linear model presented in section 2.3.2, it is necessary to substitute the numeric values listed in Appendix B for the model coefficients. For these values A = -0.5 0 0 0 0 0 1 0 0 -0.5 0 0 0 0 0 1 0 0 0 1 0 0 0 0 6.06 6.06 -0.026 -5.51 28.0 0 0 0 0 0 0 0 0 1 0 0 5.42 -5.42 -0.03 -0.89 4.51 -6.44 0 0 B = (4.55) Chapter 4. Trajectory Control 46 This is now the augmented system, with the actuator dynamics included in the A matrix. The singular values of the controllability grammian are 6.96 x 10 5 1.71 x 10 4 2.72 1.07 6.3 x 10" 2 5.49 x 10" 3 and the rank of the controllability matrix is 6. Note, the examination of the singular values only provides a qualitative understanding of the controllability and observability of the system. It is difficult to establish a singular value below which the system is not controllable or observable . However, since none of the above singular values is zero 1 or close to zero (in the order of < 1 0 -11 ) , and the rank of the controllability matrix is equal to the order of the system, the system is considered controllable. For the observability, C is initially chosen to represent the measurement of only the vertical position, so that C = [0 0 1 0 0 0 ] (4.56) In this case, the singular values of the observability grammian are 8.17 x 10 5 1.29 x 10 4 1.26 6.3 x 10" 3 8.32 x 10~ 5 2.48 x 10" 11 and the rank of the observability matrix is 5. Since one singular value is very close to zero, and the rank of the observability matrix is less than the order of the system, the system is not observable. By including a measurement of one of the foil deflections, such that 1 0 0 0 0 0 0 0 1 0 0 0 (4.57) Numerically, on computers using IEEEfloatingpoint arithmetic, the relative accuracy of numbers is about 16 significant decimal digits. As a general rule of thumb, the number of digits lost to numerical roundoff during gaussian elimination is equal to the exponent of the ratio of the largest singular value to the smallest. With the largest singular value always about 10 , this means singular values in the order of 10 will result in a singular system 5 -11 Chapter 4. Trajectory Control 47 the singular values of the grammian become 8.17 x 10 1.29 x 10 5 1.26 4 7.22 x 10" 6.22 x 10" 1 8.14 x 10~ 3 5 and the rank of the observability matrix is 6. This value of C yields an observable system. Note, that the additional measurement has moved the smaller singular values away from zero. In the limit of measuring all the states, the singular values are 8.37 x 10 5 1.32 x 10 1.49 4 1.02 1.44 x 10" 1 5.76 x 10~ 2 and have been moved as far from zero as possible. This process is often referred to as increasing the systems degree of observability. For the numerical simulations, three observer designs are examined. The first design is with C = r 1 0 0 0 0 0 0 1 0 0 0 0 (4.58) and the singular values shown above. The second design uses vertical position and body rotation measurements, such that C = 0 0 1 0 0 0 0 0 0 0 1 0 (4.59) and the singular values are 8.18 x 10 s 1.29 x 10 4 1.27 7.78 x 10" 2 7.21 x 10" 3 5.8 x 10~ 4 The final observer design adds the foil measurements to the second design, such that 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 (4.60) Chapter 4. Trajectory Control u 48 o B x C y Figure 4.16: Schematic of Regulator Feedback Law and the singular values are 8.18 x 10 4.3.2 5 1.29 x 10 4 1.27 1.00 7.43 x 10" 3 7.14 x 10" 4 State Feedback D e s i g n A basic premise of linear state feedback design is that given a controllable, timeinvariant system x = Ax + T5u y = (4.61) Cx it is possible to choose an input of the form u = —Gx that will place the resulting closed loop poles at any desired location. The resulting regulator design is schematically shown in Figure 4.16. While this is a multi-input multi-output system (MIMO), it is possible to gain considerable insight into the advantages of M I M O systems by briefly reviewing the methods of pole placement for single-input single-output systems (SISO). For a SISO system, the Bass-Gura formula[16] can be used for determining G. The method Chapter 4. Trajectory Control 49 is based on using the polynomial coefficients ct of the characteristic equation n det(sl - A ) = s + ctis"- + ... + <*„ n 1 and the desired eigenvalues of the system A i , . . . , A„. The eigenvalues can be used to determine the desired characteristic equation (s - - A ) • • • (s - A„) = s + ens"- +••• + «„ n 1 2 with which G can be expressed as Ctl G = [(QW)'] (4.62) Ct , r where Q is the controllability matrix and W = 1 ai ••• a„_i 0 1 • • • Ct -2 0 0 ••• n (4.63) 1 If this method is considered for a M I M O system, the calculation of G is underdetermined, with more gain values to solve for than eigenvalues. As previously mentioned, a simple way to deal with this problem is to set some of the gains to zero, as done in the paper by Shupe[28]. This results in a simple pole placement problem, and is straightforward to implement. However, if this situation is viewed from a different perspective, the additional freedom should be taken advantage of in the design process. Consider a gain matrix G for a n input M I M O system. From the Bass-Gura formula, it has been shown that only k gain values of the n x k values provided are needed to place the poles at any desired location. The remaining (n — !)& gains can be considered Chapter 4. Trajectory Control 50 extra, and if they can be solved for can be used to provide system optimization of some type. Linear, quadratic optimum control provides the necessary extra relationships needed to be able to solve for these gains, and in doing so, not only places the closed loop poles at the desired location, but also provides system optimization. Again consider the linear, time-invariant system of 3.55. Linear quadratic optimum control provides a strategy for finding a state feedback matrix G that minimizes the performance integral (4.64) where Q and R are symmetric positive definite state and control weighting matrices respectively. The matrix Q weights the importance of the deviation of each state from the origin. Therefore, a relatively large weighting of one state will produce a controller that will preferentially keep that state close to the origin. This capability is particularly well suited to this application, since two of the control objectives are to minimize the pitch angle and pitch rate of the vehicle. For this application, the control weighting matrix R also plays an important role. By heavily weighting the states with Q , the resulting controller may produce very large control signals. R represents the cost of the control effort. By selecting larger weightings of R , the resulting control signals can be maintained at levels that will not cause the airfoils to stall and possibly create an unstable situation. To obtain an expression for the optimum gain matrix G , consider the performance integral again. For this case it can be shown that Voo = x'(t )Mx(t ) 0 0 where M satisfies the algebraic Riccati equation 0 = M A+ A'M - MBR - 1 B'M + Q (4.65) Chapter 4. Trajectory Control 51 The resulting optimum gain matrix is expressed as G = R^BM (4.66) Friedland[16] presents a complete derivation of the gain matrix G for this infinite time solution. With G defined, the closed loop regulator system can be expressed as x - (A - BG)x (4.67) y - Cx The response of this regulator is a function of the eigenvalues of ( A — B G ) , which are determined by the selection of G . To gain some insight into how changes in the weighting matrix Q affect the calculation of G , and in turn the eigenvalues of the system, a root locus plot is presented in Figure 4.17. The figure shows how the eigenvalues move as the pitch angle weighting of Q is varied from 3 to 32769. The remaining values of Q and the values of R are held constant. Figure 4.17 shows that the system remains stable at the larger state weightings, and that the response time of the system decreases. In contrast to weighting the states, the inputs can also be weighted using R to produce the root locus plot in Figure 4.18. In this case, both terms of R are varied from 1 to 1.04 x 10 . The gain matrix is determined for each weighting, and the resulting 6 eigenvalues are plotted. Note that as R —> oo the eigenvalues of the closed loop system approach the open loop eigenvalues of -0.5 4.3.3 - 0 . 5 -2.41 - 6 . 7 +0.18t -0.18t Reference M o d e l In the previous section, an expression is presented for G based on a optimum regulator design. This satisfies two of the three control objectives, namely the minimization of Chapter 4. Trajectory Control an M Figure 4.17: Regulator Root Locus, Q(5,5) varied from 3 to 32769 Figure 4.18: Regulator Root Locus, R varies from 1 to 1.04 x 10 6 Chapter 4. Trajectory Control —e G 53 o B o .i y X s Figure 4.19: Schematic for Process and Reference Input pitch angle and pitch rate. However, the third control objective has not been addressed. The requirement of tracking a reference signal in order to maintain a constant height from the sea floor is not satisfied with a regulator design. To satisfy this tracking objective, it is necessary to introduce a reference model. The reference model is schematically shown in Figure 4.19. The control law is now u = —Ge — G a ; , where the error e is the difference between the state x and r r the reference input x . r The feedforward gain G r enables the controller to provide zero steady state error e with a nonzero input of — G x . Therefore, G » r r r r represents the control signals required to maintain the vehicle at some steady state equilibrium position other than the origin. Clearly these signals cannot be zero. It can also be shown that G is simply the inverse of the steady state of the process of 3.55. r However, G introduces a difficulty with the performance integral r oo / [x'(r)Qx(r) + t i ' ( r ) R « ( r ) ] r f r (4.68) If the input u is now not zero as a result of the contribution of G « , then the integral r r Chapter 4. Trajectory Control 54 goes to infinity as T goes to infinity. Friedland[16] presents a complete derivation of the solution to this problem, assuming the reference signals are constant. The results show that for optimum control with a reference input G = R- BM (4.69) 1 and G R = [C(A - BG)~ B]~ C(A L 1 - BG)~ A (4.70) L With these gains defined, 3.55 can be re-written as x = A x + Bit where u = —Gx + ( G — G ) x r r or (4.71) x = (A - BG)x + B ( G - G )x r r y = Cx 4.3.4 Observer Design The controller designs of the previous two sections assume the state vector x is measured and available for feedback. For this application, it is desirable to minimize the measurements to only the vertical position. The rate and angular measurements involve more expensive equipment for both sensing and data transmission. If the system is shown to be observable, as in 3.3.1, then an observer can be designed that will provide an estimate of the state vector based on only the measured output y. The classic observer design is due to Luenberger in the early 1960's, and is shown schematically in Figure 4.20. Luenbergers' method obtains an estimate of the state x from the system x = A x + Bu + Ky (4.72) which is excited by the measurement y and the input u. The matrices A , B and K are Chapter 4. Trajectory Control 55 u B y O K o- x y Figure 4.20: Schematic of Observer Design selected so as to minimize the error e = x — x. The result of this is A = A - KC B = B (4.73) which allows 3.65 to be rewritten as x = ( A - K C ) s + Bu + Ky (4.74) If the state estimate x is now substituted for the state vector x in the previous optimum state feedback design, the composite system shown in Figure 4.21 is formed. A n underlying strength of this formulation of the composite system results from the separation principle. This principal states that the gains G and K may be determined separately, without consideration of each other, and without destroying the optimality of the composite system. Note, this may no longer be true when a non-linear process is involved. From the standpoint of the numerical simulation, the composite system Chapter 4. Trajectory Control 56 TJompensaFor r B r ^ C H K K Hx. -, "Plant" -o/>f B i ~ A «- G A J r J X r Figure 4.21: Schematic of Composite System must be expressed in the form X = AX + Bx r where X is a composite state vector formed as X - x x The expression in 4.71 is rewritten, replacing the input with u = —Gx — ( G — G )cc r r to yield x = Ax - BGx + B ( G - G ) x r (4.75) r In a similar manner, the expression in 4.72 is rewritten by making the substitutions for u = —Gx — ( G — G ) x r r and y = C x to yield x = (A - B G - KC)x + K C x + B ( G - G )x r r (4.76) Combining expressions 4.75 and 4.76 into matrix form yields the final expression X = AX + B x r Chapter 4. Trajectory Control A = 57 ( A - B G - K C ) KC -BG A B(G — G ) r B = (4.77) B(G - G ) r It is interesting to note that the eigenvalues for this composite system are simply the combined eigenvalues of the closed loop process ( A — B G ) and the observer ( A — K C ) . This becomes clear by denning the observer error e — x — x. With this, equation 4.75 becomes x = ( A - B G ) « + B G e + B ( G - G )x r r (4.78) In a similar manner, since e = x — x, equation 4.76 becomes e = (A-KC)e (4.79) While these two equations above completely describe the composite system, they also clearly represent two separate dynamic systems connected in series. Therefore, the response of the composite process will be determined by the eigenvalues of ( A — B G ) and ( A - K C ) . Chapter 5 Numerical Simulation 5.1 Overview With the mathematical models developed in Chapter 2 and the control strategies presented in Chapter 3, it is now possible to consider a numerical simulation of the system. However, it is first necessary to consider the limitations of this, and indeed any, numerical simulation. The intent of this simulation is not to attempt to predict the exact response of the towed vehicle during operational towing. There have been numerous assumptions made during the derivation of both the linear and nonlinear models. Some of these assumptions relate to tow cable geometry and response, some relate to modelling errors of the vehicle's hydrodynamic coefficients. In addition, the exact nature of the required tracking input is unknown, and measurement noise and system disturbances have been neglected. However, even with all of these unmodelled effects, a numerical simulation will show the qualitative, and to some degree the quantitative, nature of the system response. If a simulation predicts reasonable system behaviour, then a proper implementation of the system is likely to behave in the same, reasonable manner. It is, therefore, the intent of these simulations to show that the modeled system has a stable, well behaved response, and with the proper selection of the gain matrices, will satisfy the required control objectives. The simulations will also show how the selection of the measurement states for the observer design affect the performance of the vehicle. 58 Chapter 5. Numerical Simulation 59 In order to perform this numerical analysis, two simulation environments are used. The first of these is P C - M A T L A B . P C - M A T L A B is a scientific and engineering program developed by MathWorks Inc. to run on IBM and other MS-DOS compatible personal computers. P C - M A T L A B provides all of the necessary control design features, and provides an excellent environment for the linear simulations. While P C - M A T L A B is capable of nonlinear simulations, the features are not as extensive as provided by A C S L . A C S L , or Advanced Continuous Simulation Language, was developed by Mitchell and Gauthier, Associates. The language was designed to simulate continuous systems described by time dependant, non-linear differential equations. The program is cur- rently running on the Mechanical Engineering Department's V A X 750. This increased processing speed and the comprehensive programming language make A C S L the better environment for the non-linear simulations. The general procedure for a particular simulation run can be divided into the following four stages: 1. The linear state feedback design is performed in the P C - M A T L A B environment using the program S I M U L A T E , described in the following section. Either a state feedback or an observer design is performed. 2. The simulation of the linear model is performed using S I M U L A T E . 3. The gain matrices K and G are transferred to the A C S L program S S C A N , and the simulation of the non-linear system is performed. 4. The results of the non-linear simulation are transferred back to the P C - M A T L A B environment and directly compared to the linear simulation results. Comparing the non-linear and linear results indicates how well the non-linear model behaves with a linear control strategy. Chapter 5. Numerical Simulation 60 PROGRAM INITIAL §. Statements to define initial conditions and constants END DYNAMIC DERIVATIVE § Statements to describe the dynamic model END §. Statements executed every communication interval END END Figure 5.22: Program Structure for S S C A N . C S L 5.2 A C S L Environment A C S L is the continuous simulation language used for the non-linear simulations. The language provides a programming environment similar to that of F O R T R A N , yet has several differences. The language provides for free form input, and employs an autosorting routine to ensure that variables are not used before they are calculated. There is also a very comprehensive collection of special functions defined, and the option of using one of five predefined integration routines. A complete description of these functions and integration routines is presented in the A C S L Users Manual[21]. S S C A N . C S L is the A C S L program developed to provide the non-linear simulations (see Appendix F for program listing). The program structure is shown in Figure 5.22, and is typical of most A C S L programs. In the first section, INITIAL, the following is performed: • All of the constants are denned. This includes the stream velocity, body geometry, feedback gains, initial conditions, integration parameters, and logicals that control program execution. Chapter 5. Numerical Simulation 61 • Any variables that remain constant during the integration are determined here. This includes the mass and moment of inertia of the vehicle, the added mass and moment of inertia, and the finite aspect ratio correction. Three logicals are used in this program, munk, Ikup, and dwash. If munk is set to be true, the program uses only the slender body theory component of the body lift. If munk is set to be false, the least square polynomial fits of the experimental results in section 2.4.3 are used to determine the body lift and moment. The logical Ikup is used to determine how the cable angle theta is calculated. If Ikup is set true, a lookup table based on the steady state cable modeling is used. If Ikup is set false, the non-linear model will use the same linear approximation used by the linear model. If dwash is set to be true, the program calculates the downwash affects of the front foil on the rear foil. This calculation uses a lookup table based on the experimental results presented in section 2.4.3. If dwash is set to be false, downwash affects are neglected. In the Derivative section the dynamic model is described. The following outlines the main calculations required to define the non-linear model: • The induced angle of attack 7; is calculated for each pair of foils. • The velocity magnitude Vi is calculated. • The reference input is calculated. This is selected to be a ramped input to a constant value. • The foil rotations are calculated based on the gain values K and G determined from S I M U L A T E . A first order model of the actuator is modeled using the A C S L function REALPL. Chapter 5. Numerical Simulation 62 • The angle of attack ipi is determined. • The lift coefficient slope is calculated, based on the logical setting dwash. • The net lift L,- of each foil is determined and separated into an x and y component. • The x component T of the cable tension is calculated. This value, along with x the value of y, is used to enter into the lookup table for the cable angle 8. The vertical component of the cable tension is then calculated. • The body normal force and moment are determined, based on the logical munk. • The net vertical force and moment acting at the body center of gravity are determined, and the vertical and angular accelerations are calculated. • Using the default fourth order Runge-Kutta integration routine, the vertical and angular velocities and positions are calculated. The observer model performs the same calculations as above, with the addition of the following: • Using the gains G and K, and the linear process matrices A, B, and C the derivatives of the observer states are determined. • Using the default integration routine, the observer states are calculated. These linear observer states replace the non-linear process states in the calculation of the foil rotations. The last part of the program, which is outside of the Derivative section, contains the statements which control the termination of the program. In this case, if the time exceeds the input duration of the simulation, the program will terminate. Chapter 5. Numerical Simulation 63 It is necessary to run the A C S L program S S C A N . C S L through the A C S L processor before running the simulations. This procedure performs the auto-sorting of the statements, and produces a standard Fortran program which is then compiled and linked to any additional graphics libraries. The resulting executable program inputs a data file which is used to perform the following: • Any predefined constants or logicals may be assigned new values. • The output variables are chosen. • Any output plots or file printouts may be specified. The output plots are displayed on Tektronix 4010 compatible displays, while the output files are transferred to the P C - M A T L A B environment for direct comparison with the linear simulation results. 5.3 P C - M A T L A B Environment P C - M A T L A B is used to develop the linear process model and, using this model, design the linear controller described in Chapter 3. To facilitate this model development, P C M A T L A B creates a scientific, interactive programing environment. The fundamental object in this environment is an undimensioned rectangular matrix with possibly complex elements. Fundamental math operators, elementary math functions, and more advanced scientific functions are all defined to operate on this fundamental matrix object. This provides for a very powerful programing environment. For a complete listing of the P C - M A T L A B functions, see the P C - M A T L A B Users Guide[22]. In addition to the basic scientific operators of the kernal P C - M A T L A B program, there are additional toolboxes provided by MathWorks to enhance the capabilities of Chapter 5. Numerical Simulation 64 P C - M A T L A B in particular engineering disciplines. Currently available are the Control, System Identification, and Signal Processing toolboxes. The additional functions provided by the control toolbox are used extensively by the simulation program SIMULATE.M. S I M U L A T E . M is the main program developed in M A T L A B . This program, and the associated subprograms are referred to as script files. When these program names are entered in the M A T L A B environment, the files are simply executed one line at a time. Any variables defined during the execution of the program are stored in memory, and available for use by any of the subprograms. This provides for very fast and flexible program development, at the expense of slower execution speeds and a maximum size and number of defined variables. Figure 5.23 shows a flow chart for S I M U L A T E . M . The following briefly describes the tasks of the principle subprograms. Simulate Provides for the selection of the underlying 8 principle subprograms. These must initially be executed in order to correctly define all the variables. S_state This script file defines the body parameters, calculates the added mass and moment of inertia, and determines the aspect ratio corrections. Then, using the files heave.m and pitch.m, the linear dynamic model is assembled. L q r _ p a r Provides the option of preferentially weighting one of the states before performing the linear quadratic regulator (LQR) design. The state and the weighting value are entered. Sim_opt Offers the option of augmenting the fourth order system created by S jstate.m with actuator dynamics. The actuator is modeled as a first order system, requiring only a time constant to be entered. Chapter 5. Numerical Simulation 65 Simulate Partlist S_state Lqr. .par Sim .opt Solve Observ UJnput Heave StJist Act.dyn New_st Refer Trans Newplot Attac Newplot Pitch P_vel P_pos PJbill P_foil2 P_omg P_ang P.attacl P_attac2 Plotlist Makemet Figure 5.23: Flow Chart for P C - M A T L A B program S I M U L A T E . M Chapter 5. Numerical Simulation 66 U _ i n p u t This script file is used to generate the reference input matrix. The user selects either a ramp to a constant value input, or a sine function input. Non-zero initial conditions may also be set. T r a n s Performs the translation of the A C S L output file to defined M A T L A B variables. N e w p l o t This is the kernal plotting routine. The number of plots and the variables to be plotted are first selected. It is also possible to have a set of non-linear results superimposed on the linear results. The appropriate P_*.m script files are called to plot the selected results. After plotting, a M A T L A B meta file may be created from the plot on the screen. This device independent meta file is then used by a post-processing graphics routine to create a device dependant printable file. Solve This is one of the main computational script files. First, the L Q R design is performed using M A T L A B ' s Iqr function. New_st.m takes the regulator gain and solves for the resulting closed loop matrices, while Refer.m adds the reference input. The final system is then simulated using the function him. The angles of attack of the airfoils are calculated in Attack.m, and Newplot.m is called for the plotting of the results. O b s e r v This script file is similar to Solve.m, with the addition of performing the observer design. The function Iqr is used to determine the observer gain matrix C. For this calculation, state weightings are much larger than the weightings for the regulator design, ensuring the response of the observer is faster than the system. The simulated system is the twelfth order composite system presented in section 3.3.4. This file also loads the observer states into the non-linear variables, allowing them to be plotted by requesting the superposition of the non-linear results. Chapter 5. Numerical Simulation 67 In addition to the main program S I M U L A T E , a shorter routine, U C - R O L L , performs the linear uncoupled roll simulations. This program uses the script file R O L L _ S T to calculate the vehicle parameters. The linear model is then assembled and augmented with the selected actuator dynamics. M A T L A B ' s L Q R design routine is used to perform the pole placement for the system, and the simulation is performed with the Isim routine. Chapter 6 S i m u l a t i o n Results 6.1 Overview The simulation results are presented in five sections. The results of the first section are for the full state feedback design. This assumes that all of the states are available for measurement. For all of the simulations except roll, the reference input is a ramp from 0 to 10 feet in 10 seconds, which then remains constant for a further 10 seconds. First, the linear quadratic regulator (LQR) weightings for the state feedback design are varied. These results indicate the order of magnitude of the L Q R weightings and the resulting state gains necessary to achieve the desired control objectives outlined in section 3.2. Next, the actuator time constant is varied. These results indicate how the vehicles performance varies for fast and slow actuator time constants, and is a useful criteria in selecting the vehicle actuators. The final set of results for the full state feedback model present the relative importance of the three principle non-linear affects; airfoil downwash, body normal force and moment, and cable angle. Three consecutive simulation runs are presented. The first run neglects the downwash affect, the second neglects downwash and non-linear body forces, and the third run neglects all three affects. The second section of results present the vehicle's ability to recover from a given set of initial conditions and return to the desired trajectory. The results are presented for two L Q R weightings, showing the ability to recover as a function of the chosen state 68 Chapter 6. Simulation Results 69 weightings. Section three presents two simulation results that indicate the robustness of the full state feedback design to velocity fluctuations. While this is not an exhaustive examination of the vehicle robustness, the velocity fluctuation is considered to be the largest and most common parameter fluctuation. The results presented in Section four show how the selection of the measurement states affect the vehicle performance. While both simulations were shown to be observable systems in section 3.3.1, the first simulation showed poor vehicle response. The final section presents a simulation of the linear roll model. As previously mentioned, this model is only intended to provide insight into the magnitude of foil rotations required to correct for roll rotations. 6.2 State Feedback w i t h Reference Input 6.2.1 V a r y i n g the L Q R Weightings As described in section 3.3.2, the optimum control strategy allows the individual states and inputs to be preferentially weighted. Figure 6.24 shows the results of the first simulation run in which all of the states and inputs have a weighting of one. The actuator time constant for this and the following run is arbitrarily chosen as 2 seconds. In the figure, the solid line represents the linear model results, and the dashed line represents the non-linear A C S L simulation results. The results show that the non-linear model behaves well with the linear controller, and that even for this first selection of weightings the control objectives for the general side scan application can be meet. In addition, the following is shown: 1. The steady state response for the linear model is shown to be zero. This is a result of the feedforward gain G . Note, the pitch angle steady state error for r Chapter 6. Simulation Results sec. Pitch Angle .4I 0 .21 0 • 5 1 10 sec . 15 1 20 .1 2 0 1 5 • 5 L-i! 10 sec 1 10 sec , I 15 20 • 1 15 20 Attack Foil 1 sec sec Figure 6.24: Linear (solid) and Non-linear( dashed) Results for Simulation 1 Chapter 6. Simulation Results 71 the non-linear model is not zero, since G is based on the linear model only. The r non-linear vertical position steady state error is close to zero, and is not a tracking problem. 2. The angles of attack of the foils remain well within their linear range, and are possibly even too small. Since this particular maneuver represents close to the maximum desired performance, the foil areas could be reduced somewhat to increase the angles of attack. Note, the difference in the linear and non-linear angles of attack for foil 2 is a result of the downwash effect. In the non-linear model, foil 2 is operating in the downwash of foil 1. Therefore, to achieve the required lift it must be at a greater angle of attack than is predicted in the linear model. In simulation 2, the L Q R weighting for the pitch angle is increased to 60. The results in Figure 6.25 show that the pitch angle can be decreased to the lowest value specified in the control objectives. In order to achieve this improved performance, the angles of attack have increased. However, the attack angles are still well within their linear range, and again suggest the initial foil size could be decreased. The effect of downwash is more noticeable in simulation 2. Foil 1 is generating greater lift than in the first simulation, and produces a stronger downwash effect on foil 2. This is evident from the increased separation of the linear and non-linear curves for foil 2 angle of attack in Figure 6.25. The pitch angle weighting of 60 is chosen arbitrarily to yield pitch angles less than the control objectives. While larger weightings will produce what might appear to be even better results, this may not be the case. Neither of these simulation models consider process or measurement noise. The high gain system resulting from very large weightings may have an undesirable response when noise is added to the system. Chapter 6. Simulation Results sec Foil 1 Rotation ' sec sec ' 1 Foil 2 Rotation l n 1 0 l 1 ' r sec Figure 6.25: Linear(solid) and Non-linear(dashed) Results for Simulation 2, Pitch anj Weighting is 60 Chapter 6. Simulation Results 6.2.2 73 Varying Actuator T i m e Constant The time constant for the first two simulations is arbitrarily chosen as 2 seconds. To see the effect of varying this time constant, simulation 3 uses a time constant of 0.5 seconds. The pitch angle weighting is 60. The results in Figure 6.26 show that this reduced time constant further minimizes the pitch angle. However, as the pitch angle is continually reduced from the results of simulation 1 to simulation 3, there is more of a difference between the linear and non-linear pitch angles. While reducing the time constant improves the simulated vehicle response, there is a point at which the validity of both simulations break down. Neither simulation model considers the affects of unsteady aerodynamics. In its simplest approximation, unsteady affects will limit the rate at which lift can be generated by the airfoils. Actuator time constants or state weightings which require the airfoils to create lift at a rate exceeding this limit will produce misleading results. 6.2.3 R e l a t i v e I m p o r t a n c e of the P r i n c i p l e N o n - l i n e a r Effects This section considers the importance of including the downwash, non-linear body force, and non-linear cable angle calculation in the non-linear process model. By sequentially omitting each non-linear effect, it is possible to observe the nature in which each one influences the response. This is useful in determining how uncertainties in the calculation of these non-linear effects is likely to influence the system response. In addition, if further optimization of the response is required, these results will show which nonlinear effects need to be accounted for in a more comprehensive control strategy. For the following three simulations, the actuator time constant is set at 2.0 seconds and the pitch angle weighting is 60. In simulation 4, the downwash affects have been removed from the non-linear simulation. These results are shown in Figure 6.27, and indicate Chapter 6. Simulation Results 74 sec Attack FoU 1 sec . Attack Foil 2 sec Figure 6.26: Linear(solid) and Non-linear(dashed) Results for Simulation 3, Time Constant Reduced to 0.5 sec Chapter 6. Simulation Results 75 sec Figure 6.27: Linear(solid) and Non-linear(dashed) Results for Simulation 4, time constant of 2 sec, body weighting of 60, downwash neglected Chapter 6. Simulation Results 76 that the downwash influences the results in the following way: 1. The apparent increase in the angle of attack for foil 2 in Figure 6.25 is a result of the downwash. In Figure 6.27, the linear and non-linear curves have collapsed close to each other. 2. The differences in attack angles for the foils introduces a substantial difference between the linear and non-linear pitch angle. Neglecting downwash reduces this difference, as shown in Figure 6.27. In simulation 5, both the downwash and the non-linear body forces have been neglected from the non-linear model. These results are shown in Figure 6.28, and show that the non-linear component of the body force affects the vehicle motion greatest when the vertical velocity is the greatest. Therefore, the non-linear component of the body force has little affect on the steady state response, as evident from comparison of Figures 6.27 and 6.28. The last of these three simulations, simulation 6, neglects downwash, non-linear body forces, and non-linear cable angles. These results are shown in Figure 6.29, and show that the non-linear cable affect is greatest as the vehicle moves farther from the linearization point. For this particular reference input, this is reflected in the steady state error, which for the non-linear results shown in Figure 6.29 is less than the nonlinear results of Figure 6.28. 6.3 State Feedback w i t h Initial C o n d i t i o n s Simulations 7 and 8 show how well the vehicle responds to non-zero initial conditions. The reference input is still applied, with an actuator time constant of 2 seconds. For both simulations, the vertical position is given a 1 foot initial condition, and the initial Chapter 6. Simulation Results 77 Figure 6.28: Linear(solid) and Non-linear(dashed) Results for Simulation 5, downwash and non-linear body forces neglected Chapter 6. Simulation Results 78 Figure 6.29: Linear (solid) and Non-linear(dashed) Results for Simulation 5, downwash, non-linear body forces and cable angles neglected Chapter 6. Simulation Results 79 condition for body rotation is 10 degrees. The results for simulation 6 are shown in Figure 6.30, and are for a pitch angle weighting of 1. These results indicate that the vehicle recovers well from the initial conditions, and quickly resumes the desired trajectory. Except for the first few seconds of recovery, the maximum pitch angle does not exceed the maximum values observed for the zero initial condition simulation 1 shown in Figure 6.24. For simulation 8, the pitch angle weighting is increased to 60. The results are shown in Figure 6.31, and indicate that for the full state feedback design, increasing the state weightings does not have a detrimental affect on the results. For the observer design, this is not always the 6.4 shown in the following section. Robustness o f the State Feedback D e s i g n As mentioned, this section is not an exhaustive examination of the controller robustness. However, it is very useful to examine the performance range resulting from the largest fluctuation of any one parameter. For underwater vehicle simulations, the parameter with the greatest fluctuation is usually the stream velocity. Simulation 9 is run with a stream velocity of 2.5 ft/sec, and simulation 10 is run at 7.5 ft/sec. This velocity range represents the most likely minimum and maximum velocities that might occur for a mean towing velocity of 5 ft/sec. For both simulation runs, the actuator time constant is 2 seconds, and the pitch angle weighting is 60. The linear controller for these two simulations is designed at a stream velocity of 5 ft/sec. The results for simulation 9 are shown as the solid lines in Figure 6.32. The overall sluggish response of the vehicle is a result of the foils having to rotate to a greater angle of attack in order to generate the required amount of lift. While the resulting maximum pitch angle is still within the loosest control objectives, it may not be adequate for the Chapter 6. Simulation Results sec 80 sec Figure 6.30: Linear( solid) and Non-linear (dashed) Results for Simulation 7, time constant of 2 sec, body weighting of 1, with initial conditions Chapter 6. Simulation Results 81 Figure 6.31: Linear(solid) and Non-linear(dashed) Results for Simulation 8, time constant of 2 sec, body weighting of 60, with initial conditions Chapter 6. Simulation Results 82 Vertical Position Vertical Velocity Foil 1 Rotation Foil 2 Rotation 5 10 seconds 15 20 Foil 1 Angle of Attack 5 10 seconds 15 10 seconds 15 Foil 2 Angle of Attack 20 10 seconds 15 Figure 6.32: Non-linear results for Simulation 9 at 2.5 ft/sec(solid) and Simulation 10 at 7.5 ft/sec(dashed), time constant of 2 sec, body weighting of 60 Chapter 6. Simulation Results 83 higher performance applications. The dashed lines in Figure 6.32 are for simulation 10. The higher velocity of 7.5 ft/sec has the affect of slightly improving the vehicle performance. The maximum pitch angle of simulation 10 is slightly less than the nonlinear pitch angle of simulation 2. 6.5 O b s e r v e r D e s i g n w i t h Reference Input This section examines the performance of three observer designs. For all three designs, the actuator time constant is again 2 seconds, and the pitch angle state weighting is 60. For the simulations, the vertical position and the pitch angle observer states are given initial conditions of 1 ft. and 10 degrees, respectively. The reference input is unchanged from that used in the state feedback designs. The first design assumes that the front foil rotation and the vertical position are available for measurement, such that C = 1 0 0 0 0 0 0 0 1 0 0 0 (6.80) Practically, this provides the minimum number of measurement states by measuring two of the easiest states. While this value of C is shown to yield an observable system, the observability Grammian does have one singular value close to zero. Figure 6.33 shows how this singular value produces a poorly observable system, where the poor state estimates give a controller that needs to overcompensate. The results are for the linear simulation, where the solid lines are the process states and the dashed lines are the observer states. The system has a stable response, with satisfactory tracking. However, the pitch angle and the pitch rate exceed the control objectives. In addition, the foils are required to operate outside their linear range. Therefore, even though this is an observable system, the simulation results show it is an unsatisfactory design. Chapter 6. Simulation Results 84 Vertical Position 5 10 seconds 15 Vertical Velocity 20 Body Rotation 5 10 seconds 5 10 seconds 15 20 10 seconds 20 5 15 15 10 seconds 15 20 Foil 2 Rotation 20 Foil 1 Angle of Attack 5 15 Body Rotation Rate Foil 1 Rotation 5 10 seconds 5 10 seconds 15 Foil 2 Angle of Attack 20 5 10 seconds 15 Figure 6.33: Linear Process States(solid) and Observer States(dashed), First Observer Design with vertical position and front foil rotation measurements Chapter 6. Simulation Results 85 The second observer design assumes that the vertical position and the pitch angle are available for measurement, such that C = 0 0 1 0 0 0 0 0 0 0 1 0 (6.81) The simulation results in Figure 6.34 are for the linear model, and show the linear process states (solid lines) and observer states (dashed lines). The overall vehicle response is better than the first observer design, with the pitch angle and pitch rate remaining within the control objectives. However, the maximum pitch rate occurring in the first few seconds is greater than that occurring in the full state feedback design (solid lines of Figure 6.25). This is clearly a result of the observer state having not yet converged to the process state. When the process and observer states have converged (after about 5 seconds), the process states for the observer design match closely to the process states of the full state feedback design. Figure 6.35 shows the results from the non-linear A C S L model, where solid lines are the non-linear process states and the dashed lines are the linear observer states. The results are also well behaved, and follow the trends of the linear simulation. The maximum pitch angle and pitch rates have increased slightly, but still remain within the control objectives. The sharp peaks in the airfoil angles of attack are also due to the observer states initially not matching the non-linear process states. This affect is slightly more exaggerated than in the linear simulation, where the maximum angle of attack for foil 2 is only about 1.75 degrees. It is interesting to note that the linear observer states do not all converge to the non-linear process states, as expected. This results in the non-linear process states of the observer design slightly differing from the non-linear process states of the full state feedback design in Figure 6.25. However, this does not appear to strongly effect the tracking of the system. Figure 6.36 shows a comparison between the states of the non-linear simulation (dashed lines) and the states of the linear simulation (solid lines) Chapter 6. Simulation Results 86 Vertical Velocity Vertical Position 5 10 seconds 15 5 Body Rotation 10 seconds 10 seconds 15 15 5 10 seconds 15 20 10 seconds 15 20 Foil 2 Rotation 20 Foil 1 Angle of Attack 5 15 Body Rotation Rate Foil 1 Rotation 5 10 seconds 5 10 seconds 15 20 Foil 2 Angle of Attack 20 5 10 seconds 15 20 Figure 6.34: Linear Process States(solid) and Observer States(dashed) for Second Observer Design, with vertical position and pitch angle measurements Chapter 6. Simulation Results 87 Vertical Position 5 10 seconds Vertical Velocity 15 5 Body Rotation 10 seconds 15 20 Body Rotation Rate T3 10 seconds 15 Foil 1 Rotation 5 10 seconds 15 Foil 2 Rotation 20 Foil 1 Angle of Attack 5 10 seconds 15 20 5 10 seconds 15 20 Figure 6.35: Non-linear Process States(solid) and Observer States(dashed), Second Observer Design with vertical position and pitch angle measurements Chapter 6. Simulation Results 88 for this observer design. The difference in the tracking error is negligible, with only a slight increase in the maximum pitch angle. The final observer design assumes that the front and rear foil rotations, the vertical position, and the pitch angle are all available for measurement, such that C = 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 (6.82) The results in Figure 6.37 are for the linear simulation, and show the linear process states (solid lines) and the linear observer states (dashed lines). The results indicate that with four measurements, the observer states converge to the process states faster than in the previous two designs. This has the affect of smoothing the response of the system and reducing the large oscillations that occur in the first few seconds of the previously described observers. Figure 6.38 shows the non-linear process states (solid lines) and the observer states (dashed lines) from the non-linear simulation. The nonlinear process states have a smoother response than the non-linear states of the second observer design. Adding two more measurements forces the observer foil rotations to converge to the non-linear state values. However, three of the remaining observer states are now farther from their non-linear process state than in the previous observer design. Overall, four measurements appear to improve the observer design. Figure 6.39 shows the linear (solid lines) and non-linear (dashed lines) states for this design. These results compare closely to the results of the full state feedback results in Figure 6.25. However, for an actual implementation, the additional cost of measuring these two extra states must be weighed against this predicted improved performance. Chapter 6. Simulation Results 89 Vertical Position 5 10 seconds 15 Vertical Velocity 20 -0.5 5 10 seconds 15 20 15 20 Body Rotation Foil 1 Rotation 5 10 seconds 15 Foil 2 Rotation 20 5 10 seconds Foil 1 Angle of Attack Foil 2 Angle of Attack seconds seconds Figure 6.36: Non-linear Process States(solid) and Linear Process States(dashed), Second Observer Design with vertical position and pitch angle measurements Chapter 6. Simulation Results 90 Vertical Position Vertical Velocity 8 5 10 seconds 15 20 Body Rotation 10 -I 10 seconds L_ 15 5 10 seconds 15 20 Body Rotation Rate 20 5 Foil 1 Rotation 10 seconds 15 20 Foil 2 Rotation 5 10 seconds 15 20 Foil 2 Angle of Attack 5 10 seconds 15 5 10 seconds 15 20 Figure 6.37: Linear Process States(solid) and Observer States(dashed), Third Observer Design with both foil rotations, vertical position, and pitch angle measurements Chapter 6. Simulation Results 91 Vertical Position 10 seconds Vertical Velocity 15 5 Foil 1 Rotation 5 10 seconds 15 10 seconds 15 15 Foil 2 Rotation 20 Foil 1 Angle of Attack 5 10 seconds 10 seconds 15 Foil 2 Angle of Attack 20 5 10 seconds 15 20 Figure 6.38: Non-linear Process States(solid) and Observer States(dashed), Third Observer Design with both foil rotations, vertical position, and pitch angle measurements Chapter 6. Simulation Results 92 Figure 6.39: Linear Process States(solid) and Non-linear Process States(dashed), Third Observer Design with four measurements Chapter 6. Simulation Results 6.6 93 R o l l Simulation The linear uncoupled roll model described in section 2.3.3 is simulated in the P C M A T L A B environment with the script file U C . R O L L . M . The basic vehicle geometry is the same as in the previous simulations, with an actuator time constant of 2 seconds. While the model is a single input system, M A T L A B ' s linear quadratic regulator design routine is still used to calculate the feedback gain matrix. This routine is used here simply for convenience. The results in Figure 6.40 are for a roll angle state weighting of 8. These results show that the roll response is well behaved, and with the proper selection of the weighting matrix, the controller can maintain the foils in their linear operating range. There is no roll control objective that relates to side scan sonar operation. However, the sum of the foil angles of attack during the tracking of the reference input and the foil rotations during roll correction must remain within the linear operating range of the foils. Chapter 6. Simulation Results 94 Roll Response _151 0 1 1 1 1 2 3 i 4 1 5 i 6 i i 7 8 i 9 I 10 seconds Figure 6.40: Linear Uncoupled Roll Response, time constant of 2 sec, roll angle weighting of 8 Chapter 7 Conclusions 7.1 C o n t r i b u t i o n o f the Thesis The contributions of this thesis can be divided into five areas of effort. The first of these is in the understanding and application of underwater tow cable modeling. The complete dynamic modeling of an underwater tow cable system is a complex problem, and has occupied the efforts of research groups for many years. It is therefore unreasonable, and perhaps even unnecessary, to try to include a complete dynamic model into a linear control strategy. For this reason, a simplified, steady state cable model is developed and included in the non-linear model of the process. In developing this steady state model, a more representative linear model is also incorporated into the linear control design. To implement this steady state model, a simple lookup table is developed. This lookup table is generated from the results of many runs of a steady state cable profile modeling program. This program determines a steady state cable profile given a set of boundary conditions, which in this case is the lift and drag applied at the free end of the cable. From the profile, it is possible to determine the vertical position and the cable angle at the end of the cable, and in turn relate this to the horizontal drag. This information is then compiled into a lookup table, and used by the non-linear model to calculate the vehicle cable angle during the simulation. The inclusion of this type of cable model with a dynamic towed vehicle model appears to be unique, and has not 95 Chapter 7. Conclusions 96 been presented in any of the reviewed literature. The second area of contribution is in the development of the two dynamic models of the vehicle. The non-linear model is developed using Newtonian Dynamics, and assumes the system is inertially uncoupled. The forces acting on the vehicle are limited to the hydrodynamic lift of the airfoils, a non-linear hydrodynamic body normal force, a body drag, and a cable tension. The calculation of the airfoil lift includes the effects of finite aspect ratios, induced angles of attack from body motions, and the downwash effects of the front airfoils. The body normal force is comprised of linear viscous component and a non-linear separated flow component. The cable tension is calculated from the net horizontal force on the vehicle, and the cable angle supplied from the previously mentioned lookup table. The body drag is simple pressure drag. The linear model is developed by retaining only the first order terms of a Taylor series expansion for the complete non-linear model. The linearization point is the equilibrium towing position, where the vertical position and body rotation are both zero. As a result of this linearization, the cable angle calculation reduces to a linear function of the vertical position, the body normal force is simplified to a linear viscous lift, and the airfoil downwash effects are lost completely. This linear model is used in the development of the various control strategies. The third area of contribution is in the wind tunnel test results. These tests are performed to provide an expression for the non-linear body normal force and body moment. Expressions for this particular body geometry do not appear to be available from current literature. These tests, therefore, have been very valuable in providing new aerodynamic test data for a slender body operating at high angles of attack in subsonic speeds. In addition, the tests are used to verify the finite aspect ratio correction for the airfoil lift, and the effect of the rear foils operating in the downwash of the front airfoils. To facilitate these experiments, a scaled model of the vehicle is constructed. Chapter 7. Conclusions 97 Also, before testing, the calibration of the newly installed test apparatus is verified. The fourth area of contribution is the development of the trajectory control strategy. The linear dynamic model is shown first to be a controllable, and with the proper selection of measurement states, an observable, system. This examination of observability is also unique to this work. The reviewed literature assumes all of the states are available for measurement. This may not always be the case, and as shown here, the measurement states may have a large affect on the performance of the system. Linear, quadratic optimum control is used to develop a regulator design. The input to this optimum controller is then modified to include a reference input. This optimum controller is then combined with the observer design to form the composite system used for the simulations. The final area of contribution is in the development of the linear and non-linear simulation capabilities. The advanced continuous simulation language A C S L 1 is used for the non-linear simulations. The comprehensive programming environment of A C S L provides a multi-dimensional lookup table routine. This routine is used to implement the cable angle calculation, and the front foil downwash effects. The use of logical switches provides the option of neglecting any of the three principle non-linear effects, and the wide selection of integration routines simplifies the necessary programming. Simulation of the linear model is performed with the scientific and engineering program P C - M A T L A B , and makes extensive use of M A T L A B ' s programming capa2 bilities. A complete menu driven system is developed to allow the user to assemble a linear model of the vehicle, select a linear quadratic state weighting, an actuator time constant, a desired reference input, and then perform either a full state feedback design or a linear observer design with pre-specified measurement states. The user then has Mitchell and Gauthier, Assoc., Inc. The Math Works, Inc. 1 2 Chapter 7. Conclusions 98 the option of plotting any of the output states, observer states and airfoil angles of attack. In addition, non-linear simulation results may be loaded and superimposed with the linear results. The linear model is developed from a Taylor series expansion of the complete non-linear model, where only the first order terms are retained. 7.2 S u m m a r y o f the S i m u l a t i o n Results The simulation results address five topics of interest. The first topic is the performance of a full state feedback controller. The results of the first two simulations, Figures 6.24 and 6.25, show that the control objectives can be met, and that increasing the state weighting further minimizes the maximum body rotation and rotation rate. The actuator time constant is reduced from 2 seconds to 0.5 seconds in simulation 3, Figure 6.26. This has the effect of further improving the response of the vehicle. However, reducing the time constant does push the simulations to the edge of their validity. Actuator time constants that require the foils to generate lift at too large a rate (see Section 5.2.2) will produce misleading results. The results from simulations 4, 5, and 6, examine the relative importance of the three principle non-linear effects. The results show that the non-linear cable effect is most prominent when the vehicle is far from the linearization point. For these simulations, this introduces a steady state error. The non-linear body normal force has no effect at steady state, and is greatest during the transient motion when the induced angle of attack for the body is greatest. The downwash effects are evident continuously during the simulation, but are particularly noticeable at steady state. This is where the front foil is generating substantial lift to maintain the vertical position, and since the body angle is close to zero, the rear foil is directly immersed in the front foil's downwash. Second, the simulations address the effect of providing the vehicle a non-zero initial Chapter 7. Conclusions 99 condition, and observing the ability of the vehicle to return to the desired trajectory. The results for simulations 7 and 8 in Figures 6.30 and 6.31, show that the introduction of an initial condition has a dramatic effect on the maximum rotation rate. This rate has increased from a maximum of 4 degrees/second for the full state feedback of simulation 1, to just under 20 degree/second. However, the peak rotation rate only lasts for 2 seconds, after which time the response closely follows the full state feedback results. Third, two simulations are performed to determine the robustness of the state feedback design to large fluctuations in the stream velocity. The results in Figure 6.32 show that even for a 50% reduction in the stream velocity, the vehicle maintains good tracking, with the body rotation remaining under 8 degrees. Also, since the foils remain well within their linear range, the initially selected foil sizes are to large. Fourth, several alternative observer designs are examined. It is shown that measuring just the vertical position and the front foil rotation yields a poor performing system, even though the system is observable. Replacing the front foil measurement with the more difficult pitch angle measurement improves the results substantially. If the number of measurements is increased to four by adding the foil rotation measurements, the results are only marginally improved. The extra two measurements produce smoother changing angles of attack of the foils. This is an improvement from the standpoint of unmodeled unsteady effects (see section 5.5). However, the maximum body rotation is slightly greater for four measurements than for two measurements. The final simulation shows the nature of the roll response of the vehicle, using one pair of foils for the control. The vehicle recovers from a 10 degree roll angle in about 5 seconds, while requiring a maximum of 2 degrees in foil rotation. Chapter 7. Conclusions 7.3 100 Needed Further Research This work shows that a proper implementation of the control strategies presented in Chapter 4 will result in a system with a well behaved response and good tracking ability. It also shows how proper selection of the state weighting will , under ideal operating conditions, maintain the control objectives. However, ideal operating condition do not often prevail, and further work is needed to extend the control strategies and vehicle models to account for these uncertainties. The control strategies should be extended to include the introduction of process and measurement noise. This would consist of determining exactly what is the best way of modeling these disturbances, and then changing the linear quadratic regulator design used in determining the observer gains, to a linear quadratic estimator design. A second useful extension to the control strategies would be to investigate the use of gain scheduling to compensate for stream velocity fluctuations. This might consist of using the lookup table capability of A C S L to determine a set of regulator gains, depending on the magnitude of the stream velocity. In this way, a closer to optimal system performance could be maintained over the large speed range. In addition, before any implementation can be considered, the effects of actuator saturation must be examined. The actual lift generated by an airfoil drops off quickly for angles of attack in excess of about 10 degrees. If the control strategy does not account for this, the controller will detect the effect of this reduced lift and simply rotate the foils even farther past their linear limit, which will in turn produce even less lift. This will produce an unstable situation. Setting a limit on the angle of attack is the most straightforward method to ensure a stable response. However, simulations would be required to determine the effect of actuator saturation (ie: when the maximum 10 degree angle of attack is reached). In addition, this places increased emphasis on the Chapter 7. Conclusions 101 ability of the non-linear process to calculate a true angle of attack. The non-linear dynamic model should also be extended to include three dimensional coupling effects. This would require research into the nature and magnitude of the coupling effects, and developing a similar control strategy for the control of sway. To effectively control sway it may be necessary to add two pairs of vertically mounted foils, and also attempt to satisfy the very stringent restrictions for yaw rotation and yaw rate. The addition of vertical fins may also precipitate further wind tunnel tests to establish the effect of airfoils operating in the wake of the body. Initially, any extensions to the model into three dimensions could simply use the existing steady state cable model. However, at some point in time it may be necessary to extend this cable model to a complete dynamic model of the tow cable system. The potential improvement as a result of this extension will need to be carefully examined, since this is considered to be a major undertaking. Appendix A L i s t i n g of M A T L A B Files '/,* Hop's L i n e a r i z e d Side Scan Sonar Simulation Program * ^ £ * * * * i | c g j c $ 3 l c * * * * * * 3 | c 3 f c * * * * * * $ ) | c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * fl 1'/, '/, 2'/, 3'/, '/, 4*/, 5'/, 6'/, 7'/, 8*/, 0- C a l c u l a t e the b a s i c s t a t e space r e p r e s e n t a t i o n u s i n g the l i n e a r i z e d model Set LQR parameters Choose the various s i m u l a t i o n options ( i e : r o o t l o c u s , f u l l s t a t e s i m u l a t i o n , i n c l u d e actuator dynamics, etc.) C a l c u l a t e the input matrix and set i n i t i a l c o n d i t i o n s Perform s i m u l a t i o n u s i n g f u l l s t a t e feedback Perform s i m u l a t i o n using observer s t a t e s Plot routine Load ACSL data (only must be done before p l o t t i n g ) Quit •/********************* simulate.m************************** echo o f f */, T h i s i s the main program f o r s i m u l a t i n g the coupled p i t c h '/, and heave motion f o r the s i d e scan sonar body. A simple '/, s t a t e feedback c o n t r o l i s p r e s e n t l y used, employing LQR design techniques to determine the optimum s t a t e feedback '/, matrix. Actuator dynamics may be included. while 1 parts=[ ' s _ s t a t e ' 'lqr.par' 'sim_opt' 'u_input' 'solve ' 'observ ' 'newplot' 'trans ' ] ; clc help p a r t l i s t disp('Run each p a r t i n sequence, followed by f u r t h e r ... choices'); qql=input('Enter the number of the p a r t to execute:'); i f ( ( q q l <= 0) | (qql>8)) break end parts=parts(qql,:); eval(parts); end y***ii(*:ic**>ii*****************s_state .m *********************** 102 Appendix A. '/, '/, '/, Listing of MATLAB Files s c r i p t f i l e t h a t c a l c u l a t e s the matrix values of A,B,C and D i n the s t a t e space r e p r e s e n t a t i o n of the coupled p i t c h and heave motion, where: '/. '/. X = A x + B u Y = C x + D u raddeg=57.29578; W=150.0; SPAH1-1.0; CORD1=0.5; SPAN2-1.0; CORD2=0.5; spant=l.0; cordt=0.5; areat=spant*cordt; LD=2.0; LM=2.0; Ll=0.5; L2=3.5; L=4.0; DIA=0.75; area=0.25*pi*DIA*DIA; V=5.1; q=V~2; '/, C a l c . the mass and the added mass of the body Ml=(W/32.2)+(pi*L*2.0*DIA*DIA/4.0); */, C a l c . the added mass of the f o i l s mff=2.0*(pi*CORDl~2)*SPANl/4.0; mrf=2.0*(pi*C0RD2~2)*SPAN2/4.0; T o t a l mass =Ml+mff+mrf; BDR=0.4*pi*DIA*DIA*V*V/4•0; AR1=SPAN1/C0RD1; AR2=SPAN2/C0RD2; art=2.0*spant/cordt; KCl=2.0*pi/(1.0+(2.0/ARl)); KC2=2.0*pi/(1.0+(2.0/AR2)); kct=2.0*pi/(1.0+(2.0/art)); '/, C a l c . moment of i n e r t i a f o r the body IM0=W*(DIA*DIA/4.0+L*L/3.0)/(4.0*32.2); '/, C a l c . added moment of i n e r t i a of f r o n t f o i l s X m i . f f . c g = m i _ f f + (mff * d_ff~2) mi_ff=(0.5*mff*C0RD1~2)/4.0; d_ff=LM-Ll-(CORDl/4.0); mi.ff_cg=mi_ff+(mff*d_ff~2); '/, C a l c . moment of i n e r t i a of r e a r f o i l s mi.rf=(0.5*mrf*C0RD2~2)/4.0; d_rf=L2-LM+(CORD2/4.0); mi.rf_cg=mi_rf+(mrf*d_rf"2); '/, Correct moment of i n e r t i a of body f o r added mom. of '/ inertia IM=IM0*1.8; '/, Add moment of i n e r t i a of f o i l s IM=IM+(2.0*mi_ff_cg)+(2.0*mi_rf_cg); Zl=2.0*SPAN1*CORD1*V*V; Z2=2.0*SPAN2*CORD2*V*V; '/, O r i g i n a l l y a l i n e a r approx. of the viscous f o r c e was used '/, This has now been r e p l a c e d by a slender body normal f o r c e Appendix A. Listing of MATLAB Files '/, of 2 times alpha (alpha i s the induced angle of attack of '/, the body) T h i s i s implimented with terms z6-z9 '/, k_x2=input('Enter slope of l i n e a r approx. f o r x2 viscous '/. force='); '/. Z3=k_x2; Z3=0.0; Z4=0.4*DIA*(L-2.0*LM)"3; '/, A l i n e a r approx. f o r the a d d i t i o n a l normal f o r c e s from */, omega i s s t i l l included, even though t h i s term has '/, been found t o have l i t t l e e f f e c t on the s o l u t i o n . k_x4=input('Enter slope of the l i n e a r approx. f o r x4 ... viscous force='); Z5=k_x4; '/, norm, f o r c e = nfc*Q*A Q=dynam. p r e s s . , A=cross '/, s e c t i o n a l area, nfc=normal f o r c e coef., Munk's '/, approx. */. nf = (z7 * x3)-(z6 * x2) nfc=3.0; z6=nfc*V*area; z7=nfc*V*V*area; '/, p i t c h i n g moment = nf * l s l s = d i s t . of nf from cm, '/, t y p i c a l l y at the n o s e / c y l i n d e r j u n c t i o n */, pm = (z8 * x3)-(z9 * x2) ls=2.0; z8=ls*z6; z9=ls*z7; heave pitch '/, Form the composite A and B matrices A=[heave_m;pitch_m]; B=[heave_in;pitch_in]; C=eye(4); D=zeros(4,2); d i s p C A l l done! ! ') ; y^********************** heave.m ************************** '/, S c r i p t f i l e t o c a l c . the c o e f f . i n the heave equation '/, of motion '/, t o o b t a i n a b e t t e r approx. f o r cable angle, assume */, t h e t a = y/20 i n s t e a d of y/40 '/, C a l c . c_yy due t o v e r t , p o s i t i o n c_yy=-BDR/20.0; '/, C a l c . c_yv due t o v e r t . v e l . c_yv=-((Zi*KCl/V)+(Z2*KC2/V)+Z3+z6); '/, C a l c . c_ya due t o p i t c h c_ya=((Z1*KC1)+(Z2*KC2)-(BDR*LM/20.0)+z7); '/, C a l c . c_yq due t o p i t c h r a t e c_yq=-((Z1*KC1*(LM-L1)/V)+(Z2*KC2*(LM-L2)/V)+Z4); '/, Form the 2x4 matrix from the heave equation heave_m=[0 M O O ;c_yy c_yv c_ya c_yq ]/M; '/, Form the input matrix c_yf=Zl*KCl; c_yr=Z2*KC2; Appendix A. Listing of MATLAB Files heave_in=[0 0 ;c_yf c_yr]/M; y , * * * * * * * * * * * * * * * * * * * * * * p i t c h . i i i **************************** '/, s c r i p t f i l e t h a t c a l c u l a t e s the c o e f f . of the p i t c h '/, equation of motion '/, t o o b t a i n a b e t t e r approx. f o r cable angle, assume '/, t h e t a = y/20 i n s t e a d of y/40 '/, C a l c . the c o e f f . from v e r t i c a l p o s i t i o n c_my=-BDR*LM/20.0; '/, C a l c . c_mv due t o v e r t . v e l . c_mv=-((Z1*KC1*(LM-L1)/V)+(Z2*KC2*(LM-L2)/V)+z8); Calc. c_ma due t o p i t c h c_ma=(Zl*KCl*(LM-Ll))+(Z2*KC2*(LH-L2))-(BDR*LM*LM/20.0) +z9; ... '/, C a l c . c_mq due t o p i t c h r a t e c_mq=-((Z1*KC1*(LM-Li)~2/V)+(Z2*KC2*(LM-L2)"2/V)+Z5); '/, C a l c . c_mfl, c_mfr, and c_mr c_mf=Zl*KCl*(LM-Ll); c_mr=Z2*KC2*(LM-L2); '/, Form the 2x4 matrix from the p i t c h equation pitch_m=[0 0 0 IM ;c_my c_mv c_ma c_mq ]/IM; '/, Form the 2x2 input matrix from the p i t c h equation pitch_in=[0 0 ;c_mf c_mr ]/IM; '/,************************ lqr_par.m ************************ '/, S c r i p t f i l e t o determine the LQR parameters '/, Input a r e : s t a t e and value of weighting [h,hh]=size(A); [s,ss]=size(B); Q=eye(h+ss); R=eye(ss); help s t . l i s t st_lqr=input('Input t h e s t a t e of Q t o weight:'); st_wgt=input('Enter value of weighting:'); Q(st.lqr,st_lqr)=st_wgt; ^***********************************l|ll|l********)|C)|t^ll|ll|l4t*4c4C!|c4C)|C** '/.* '/.* '/,* '/,* 1- Front f o i l 3- Heave 5- P i t c h angle STATE LIST 2- Rear f o i l 4- Heave r a t e 6- P i t c h r a t e * * * * 7!***$$$$$$************************************************** '/*********************** sim_opt.m ************************* '/, S c r i p t f i l e t o determine the s i m u l a t i o n options '/, Input are: Option f o r root locus p l o t s o r f u l l s t a t e Appendix A. Listing of MATLAB '/, '/, '/, Files s i m u l a t i o n , option t o solve Y as the input matrix (must be s e l e c t e d t o enable f o i l response p l o t s ) , option t o include actuator dynamics. '/, Option t o i n c l u d e f o i l actuator dynamics act_opt=inputC Include actuator dynamics? y/n : ' , ' s ' ) ; i f act_opt == 'y', tc=input('Input actuator time constant:'); act_dyn end y,****************** act_dyn.m ***************************** echo o f f /, S c r i p t f i l e t o i n c l u d e the actuator dynamics '/, F i r s t the model f o r the actuator i s made a _ a c t = [ - l / t c 0;0 - 1 / t c ] ; b_act=[l 0;0 1]; c _ a c t = [ l / t c 0;0 1/tc]; d_act=zeros(2); '/, '/, Next t h i s model i s j o i n e d i n s e r i e s t o the dynamic open loop model [AA,BB,CC,DD] */, '/, =series(a_act,b_act,c_act,d_act,A,B,C,D); The f i n a l system uses the same l a b e l s as the o r i g i n a l model y,********************** u_input.m ************************** '/, The f o l l o w i n g i s a s c r i p t f i l e which generates '/, two types of input f i l e s : a ramp input and a */, harmonic input ( s i n ) . '/, Given: a) the f i n a l time and i n t e r v a l , '/, b) step amplitude and input s t a t e , '/, c) ramp amp., time of max. value and input s t a t e , '/, d) harm, amp., p e r i o d , and input s t a t e . u_choice=l; u_ft=input('enter the f i n a l time:'); u_int=input('enter the time i n t e r v a l : ' ) ; t=0:u_int:u_ft; u_nn=l+(u_ft/u_int); u_in=zeros(u_nn,1); while u_choice, u_choice=input('Enter 1 f o r ramp input, 2 f o r s i n input, ... 0 to quit : ' ) ; i f u.choice == 1, a_rp=input('enter the ramp amplitude:'); t_rp=inputCenter time of max. v a l u e : ' ) ; n2=t_rp/u_int; u_inc=a_rp/n2; u_in(l,l)=0.0; f o r i=l:n2. ., u_m(i+l, l)=u_in(i,l)+u_inc; end f o r i=n2+l:u_nn-l Appendix A. Listing of MATLAB end Files u_in(i+l,l)=a_rp; e l s e i f u_choice == 2, a_hm=input('enter t h e harm, amplitude:'); pd_hm=input('enter t h e p e r i o d : ' ) ; ul_in=ones(u_nn,1); ul_in=a_hm*ul_in.*sin(2.0*pi*t/pd_hm)'; u_in(:,l)=ul_in; end end p i c k _ i c = i n p u t ( ' I n i t a l c o n d i t i o n s are zero. Okay?(y/n) ... :Vs>); i f p i c k _ i c == 'y', x_ic=[0 0 0 0 0 0 ] ; else d i s p ( ' E n t e r new i n i t i a l c o n d i t i o n s matrix of the ... form'); x _ i c = i n p u t ( ' [ b e t l bet2 y yd alpha alphad] : ' ) ; end y,************************* solve.m *********************** '/, S c r i p t f i l e which performs the s o l u t i o n p a r t of t h e '/, analysis. '/, A LQR design i s used t o f i n d the optimum s t a t e '/, feedback matrix. '/, Because of t h e c y c l i n g option, t h e c a l l t o the '/, p l o t t i n g routine '/, i s made d i r e c t l y from t h i s r o u t i n e , pack '/, perform l i n e a r - q u a d r a t i c r e g u l a t o r design k=lqr(AA,BB,Q,R); '/, c a l c . new s t a t e space rep. with new_st.m new_st refer x_cy=lsim(Ac,Bc,Cc,Dc,u_in,t,x_ic); y_cy=-k*x_cy'+(k(:,3)-G_ref(:,3))*u_in'; '/, c a l c . angle of attack of f o i l s with attack.m attack newplot y ************************ new_st.m ************************ '/, S c r i p t f i l e t o convert s t a t e space r e p r e s e n t a t i o n */, of the open loop system t o a S.S. r e p r . with s t a t e '/, feedback g a i n k '/, ie: '/, X = Ax + B u , Y = Cx + Du '/, becomes '/, X = (A - B k) x + B k x r e f '/. or '/, X = Ac x + Be u '/, and Y w i l l be : '/, Y_option==0 Y = C x + D k xref t Ac=AA-BB*k; Bc=BB*k; Appendix A. Listing of MATLAB Files 108 Cc=CC; Dc=DD*k; '/• '/, S c r i p t f i l e t o introduce a feedforward path f o r the reference input Ar=zeros(6,6); E=AA-Ar• Cr=[0 0*1 0 0 0;0 0 0 0 1 0 ] ; B_ref=inv(Cr*inv(Ac)*BB)*Cr*inv(Ac); G_ref=B_ref*E; Bc=BB*(k-G_ref); '/, Modify Be so the r e f . input need only be l x l r a t h e r '/, than 6x1 Bc=Bc(:,3); Cc=eye(6,6); Dc=zeros(6,l); ^^^^ /- ^ T ^^^^^^^^^^^^^^^ ^ *T* ^ ^ ^ ^ ^ *p 'T' T "I~ I »^ — ^ ^ ^ ^ 1— ^ ^ CLw VM JIW # III ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ T 1* «V ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^If ^ "X 1 T* *P *p T* -* I T T ^T^^ A S c r i p t f i l e t o c a l c u l a t e the angle of attack f o r each X foil ydfl=(x_cy(:,4)+((LM-Ll)*x_cy(:,6)))/V; ydf2=(x.cy(:,4)+((LM-L2)*x_cy(:,6)))/V; gaml=atan(ydf1); gam2=atan(ydf2); attacl=x_cy(:,5)+(x_cy(:,l)/tc)-gaml; attac2=x_cy(:,5)+(x_cy(:,2)/tc)-gam2; X************************ X X newplot.m ********************** P l o t s c r i p t f i l e f o r s e l e c t i n g and p l o t t i n g e i t h e r s t a t e s or inputs found from running cycle.m while 1 p l o t s = [ 'p_pos ' 'p_vel ' 'p.ang ' 'p_omeg ' 'p.foill ' 'p_foil2 ' 'p.attacl' 'p_attac2']; clc help p l o t l i s t n=input('Enter number of p l o t s p e r page:l,2,4 or 0 ... to q u i t ' ) ; i f ( ( n <= 0) I (n>8)) break end d i s p ( ' E n t e r r e s p e c t i v e number of v a r i a b l e t o p l o t ' ) ; nn=input('eg: 2 or [ 1 3 4 5] : ' ) ; nnn=input('Superimpose n o n - l i n e a r r e s u l t s ? ... (y/n):','s'); i f n == 1, clg plots=plots(nn,:); Appendix A. Listing of MATLAB '/, end Files eval(plots); pause e l s e i f n == 2, clg plotl=plots(nn(l,1),:); plot2=plots(nn(l,2),:); subplot(121).eval(plotl); subplot(122),eval(plot2); pause else clg plotl=plots(nn(1,1),:); plot2=plots(nn(l,2),:); plot3=plots(nn(l,3),:) plot4=plots(nn(l,4),:) subplot(221),eval(plotl); subplot(222),eval(plot2); subplot(223),eval(plot3); subplot(224),eval(plot4); pause end C a l l make_met t o make meta f i l e of current p l o t make_met y,************************* p_pos.m ************************ '/, S c r i p t f i l e f o r p l o t t i n g the v e r t i c a l p o s i t i o n x _ c y ( l , : ) i f nnn == 'y', plot(t,x_cy(:,3),t,y_nl) t i t l e ( ' L and N-L P o s i t i o n ' ) ; plot(t,x_cy(:,3)); t i t l e ( ' V e r t i c a l Position'); end xlabeK'sec.') ; ylabel('feet'); y,*************************** p.vel.m ********************** '/, Script f i l e to plot v e r t i c a l velocity i f nnn == 'y', plot(t,x_cy(:,4),t,yd.nl) t i t l e ( ' L and N-L Vert. V e l o c i t y ' ) ; plot(t,x_cy(:,4)) t i t l e C V e r t i c a l velocity'); end xlabel('sec'); ylabel('ft/sec'); y ***********************p.ang.m *************************** '/, S c r i p t f i l e t o p l o t angular body r o t a t i o n # i f nnn =='y', plot(t,raddeg*x_cy(:,5),t,alph.nl) t i t l e ( ' L and N-L Body r o t a t i o n ' ) ; else Appendix A. Listing of MATLAB Files plot(t,raddeg*x_cy(:,5)) title('Body Rotation'); end xlabel('sec'); ylabel('degrees'); y ************************p_omeg *m ************************* '/, S c r i p t f i l e t o p l o t angular r o t a t i o n r a t e # i f nnn == 'y', plot(t,raddeg*x_cy(:,6),t,alphd.nl) t i t l e ( ' L and N-L Body Rot. Rate'); plot(t,raddeg*x_cy(:,6)) t i t l e ( ' B o d y Rotation Rate'); end xlabel('sec'); ylabel('degrees/sec'); y ************************* f o i l l . m ************************* '/, S c r i p t f i l e t o p l o t f o i l 1 r o t a t i o n and input t o f o i l 1 t i f nnn == 'y', plot(t,raddeg*x_cy(:,i)/tc,t,betl.nl) t i t l e ( ' L and N-L F o i l 1 R o t a t i o n ' ) ; else plot(t,raddeg*x_cy(:,l)/tc,t,raddeg*y_cy(1,:)') t i t l e ( ' F o i l 1 Response'); end xlabel('sec'); ylabel('degrees'); y,*************************** f o i l 2 . m *********************** '/, S c r i p t f i l e t o p l o t f o i l 2 response i f nnn == 'y', plot(t,raddeg*x_cy(:,2)/tc,t,bet2_nl) t i t l e ( ' L and N-L F o i l 2 R o t a t i o n ' ) ; else plot(t,raddeg*x_cy(:,2)/tc,t,raddeg*y_cy(2,:)') t i t l e ( ' F o i l 2 Response'); end xlabel('sec'); ylabel('degrees'); X************************* attacl.m ************************ '/, S c r i p t f i l e t o p l o t angle of attack f o r f o i l 1 i f nnn == 'y', plot(t,raddeg*attacl,t,atacl.nl) t i t l e C L and N-L Attack F o i l 1'); else plot(t,raddeg*attacl) t i t l e ( ' A t t a c k F o i l 1'); end xlabel('sec'); ylabel('degrees'); Appendix A. Listing of MATLAB Files '/Script f i l e t o p l o t angle of attack f o r f o i l 2 i f nnn == 'y', plot(t,raddeg*attac2,t,atac2_nl) t i t l e d and N-L Attack F o i l 2'); 6l SQ plot(t,raddeg*attac2) t i t l e ( ' A t t a c k F o i l 2'); end xlabeK'sec'); ylabel('degrees'); y ********************** make.met ************************** '/.Script f i l e f o r producing a .met f i l e from the current p l o t # ansl=input('Create a .met f i l e from current p l o t ? y/n ... :\'s'); i f ansl == 'y', f_name=input('Input filename enclosed i n s i n g l e ... quotes:','s'); eval(['meta',f_name]) end ans2=input('Add t o e x i s t i n g meta f i l e ? y/n : ' , ' s ' ) ; i f ans2 == 'y', meta end y V, ********$*$********************************** * Hop's P l o t t i n g Routine * '/, '/, '/, '/, 1357- t Vertical position Body r o t a t i o n F o i l 1 response Attack f o i l 1 2468- Vertical velocity Body r o t a t i o n r a t e F o i l 2 response Attack f o i l 2 •/,*********************** observ.m ************************** '/, S c r i p t f i l e which performs the s o l u t i o n p a r t of t h e '/, a n a l y s i s . A LQR design i s used t o f i n d the optimum '/, s t a t e feedback matrix. The observer gains are then '/, c a l c , and a composite system matrix i s formed with '/, the new s t a t e v e c t o r being X=[xo x x r e f ] ' pack '/, perform l i n e a r - q u a d r a t i c r e g u l a t o r design k=lqr(AA,BB,Q,R); '/, c a l c . new s t a t e space rep. with new_st.m new_st refer /, c a l c . observer gains u s i n g LQR design with weighting '/, of 70. The s t a t e s y and alpha are used i n the observer */, design co=[0 0 1 0 0 0 ;0 0 0 0 1 0 ] ; ao=AA; Appendix A. Listing of MATLAB Files qo=1000*eye(6,6); ro=eye(2,2); lo=lqr(ao',co',qo.ro); CSl=(Ac-lo'*co); CS2=lo'*co; BSl=BB*(k(:,3)-G_ref(:,3)); CS3=-BB*k; CS4=AA; CS=[CS1 CS2 ;CS3 CS4 ] ; Ac=CS; Bc=[BSl;BSl] ; Cc=eye(12); Dc=zeros(12,1); i f p i c k _ i c == 'y' x_ic=zeros(l,12); end x_cy=lsim(Ac,Bc,Cc,Dc,u_in,t,x_ic); */, Place the observer s t a t e s i n the NL v a r i a b l e s , which may '/, then be p l o t t e d against the actual, s t a t e s by asking t o '/, superimpose the n o n - l i n e a r r e s u l t s . However, l o a d i n g '/, ACSL r e s u l t s a f t e r running t h i s s i m u l a t i o n w i l l •/, o v e r r i g h t the observer s t a t e s . betl_nl=raddeg*x_cy(:,1)/tc; bet2_nl=raddeg*x_cy(:,2)/tc; y_nl=x_cy(:,3); yd_nl=x_cy(:,4); alph_nl=raddeg*x_cy(:,5); alphd_nl=raddeg*x_cy(:,6); x_cy(:,[l:6]) = []; y_cy=-k*x_cy'+(k(:,3)-G_ref(:,3))*u_in'; '/, c a l c . angle of attack of f o i l s with attack.m attack newplot y,*********************** trans.m ************************** '/, S c r i p t f i l e which w i l l convert a p r e v i o u s l y e d i t e d ACSL '/, data f i l e i n t o separate matlab v a r i a b l e s . '/, The e d i t e d f i l e must be i n the f o l l o w i n g form: •/.coll c o l 2 c o l 3 c o l 4 c o l 5 c o l 6 c o l 7 col8 c o l 9 c o l l O */, l i n e #, t ,y , yd , alph.alphd, b e t l , b e t l , a t a c l , a t a c 2 '/, The matlab v a r i a b l e s w i l l be the above names with the '/, extension _ n l (eg: y _ n l ) . disp('Has the ACSL data been p r e v i o u s l y converted t o a ... .mat f i l e ' ) ; nn=input('(y/n) : ' , ' s ' ) ; i f nn == 'y', disp('The .mat f i l e must have been copied t o ... for020.mat.'); n=input('Okay? (y/n) : ' , ' s ' ) ; i f n == 'n', disp('You may do so now...'); disp('Type !, your DOS copy command, [ r e t ] , ... [CTRL-Z] :') d i s p C e g : !copy [path]nlsim2.mat for020.mat ... [ret] [CTRL-Z] :') keyboard 112 Appendix A. Listing of MATLAB else Files end !\matlab\translate end l o a d for020 y_nl2=for020(:,3); yd_nl2=for020(:,4); alph_nl2=for020(:,5); alphd_nl2=for020(:,6); betl_nl2=for020(:,7); bet2_nl2=for020(:,8); atacl_nl2=for020(:,9); atac2_nl2=for020(:,10); c l e a r for020 '/,********************** u c _ r o l l . m ************************ '/, S c r i p t f i l e t o c a l c . the c o e f f . i n the uncoupled r o l l 7, equation of motion '/, C a l c u l a t e the b a s i c body parameters i n r o l l . i n . m r o l l _ i n '/, C a l c . I _ l the r o l l moment of i n e r t i a '/, Body alone I_ll=W*DIA~2/(8.0*32.2); */, Then add f o i l s I_lf2=(pi*mff*SPANl~2*C0RDl/4.0)+(0.5*mff*SPAN1~2); I_lr2=(pi*mrf*SPAN2~2*C0RD2/4.0)+(0.5*mrf*SPAN2~2); I.l=I_ll+I_lf2+I_lr2; */, C a l c . c _ l r cross d e r i v a t i v e due t o yaw '/, c _ l r w i l l be a l i n e a r approx. k_roll=SPANl~4*KCl*C0RDl; R_fit=0:0.002:.07; c_roll=polyfit(R_fit,k_roll*R_fit.~2,1); c_lr=c_roll(l,l); */, C a l c . c _ l p from the wings and t a i l k_wing=4.0*SPANl"3*KCl*q*CORDl/3.0; k_tail=2.0*spant~3*kct*q*cordt/3.0; c_lp=k_wing+k_tail; */, C a l c . c _ l f l and c _ l f r from l e f t c.lf=KC1*Z1*(0.5*(SPAN1+DIA)); and r i g h t f r o n t foils '/, Form the 2x2 matrix c o n t r i b u t e d by the r o l l equation aa=[ 0 1.1; 0 - c _ l p ] / I _ l ; '/, Form the 2x1 input matrix from the r o l l bb=[0 ; c _ l f ] / I _ l ; equation cc=eye(2); dd=[0;0]; '/, Add actuator time cons. tc=input('Enter actuator time constant: ' ) ; a_act=-l/tc; b_act=1.0; c_act=l/tc; d_act=0.; [ac,bc,cc,dc]=series(a_act,b_act,c_act,d_act,aa,bb,cc,dd); Appendix A. Listing of MATLAB Files '/, Perform LQR design q=[l 0 0;0 8 0;0 0 1]; r=1.0 k=lqr(ac,bc,q,r); ac=ac-bc*k; bc=[0;0;0] cc=eye(3); dc=bc; u_input xx=lsim(ac,bc,cc,dc,u_in,t,x_ic); yy=-k*xx'; plot(t,raddeg*xx(:,2),t,yy*raddeg/2.,t,raddeg*xx(:,1)/ ... (tc*2.0)) t i t l e ( ' R o l l Response'); xlabel('seconds'); ylabel('Degrees'); '/,************************** r o l l _ i n . m ********************* '/, S c r i p t f i l e t o c a l c u l a t e the b a s i c body parameters '/, used by u c . r o l l . m t o model the uncoupled r o l l motion raddeg=57.29578; W=150.0; SPAN1=1.0; CORD1=0.5; SPAN2=1.0; C0RD2=0.5; spant=l.0; cordt=0.5; areat=spant*cordt; LD=2.0; LH=2.0; Ll=0.5; L2=3.5; L=4.0; DIA=0.75; area=0.25*pi*DIA*DIA; V=5.1; q=V~2; '/, C a l c . the mass and the added mass of the body Ml=(W/32.2)+(pi*L*2.0*DIA*DIA/4.0); '/, C a l c . the added mass of the f o i l s mff=2.0*(pi*C0RDl~2)*SPANl/4.0; mrf=2. 0* (pi*C0RD2~2) *SPAN2/4.0; T o t a l mass =Mi+mff+mrf; BDR=0.4*pi*DIA*DIA*V*V/4.0; AR1=SPAN1/C0RD1; AR2=SPAN2/C0RD2; art=2.0*spant/cordt; KCl=2.0*pi/(1.0+(2.0/ARl)); KC2=2.0*pi/(1.0+(2.0/AR2)); kct=2.0*pi/(1.0+(2.0/art)); '/, C a l c . moment of i n e r t i a f o r the body IM0=W*(DIA*DIA/4.0+L*L/3.0)/(4.0*32.2); '/, C a l c . added moment of i n e r t i a of f r o n t f o i l s 114 Appendix A . Listing of MATLAB Files '/, m i . f f . c g = m i . f f + (mff * d_ff~2) mi_ff=(0.5*mff*C0RD1"2)/4.0; d_ff=LM-Ll-(C0RDl/4.0); mi.ff _cg=mi_ff+(mff*d_ff2); X C a l c . moment of i n e r t i a of r e a r f o i l s mi.rf=(0.5*mrf*C0RD2~2)/4.0; d_rf=L2-LM+(CORD2/4.0); mi.rf_cg=mi_rf+(mrf*d_rf*2); '/, Correct moment of i n e r t i a of body f o r added mom. of X inertia IM=IM0*1.8; '/, Add moment of i n e r t i a of f o i l s IM=IM+(2.0*mi_ff_cg)+(2.0*mi_rf_cg); Z1=2.0*SPAN1*C0RD1*V*V; Z2=2.0*SPAN2*CORD2*V*V; '/, O r i g i n a l l y a l i n e a r approx. of the viscous f o r c e was '/, used. T h i s has now been r e p l a c e d by a slender body '/, normal f o r c e of 2 times alpha (alpha i s the induced '/, angle of attack of the body) T h i s i s implimented with */, terms z6-z9 '/, k_x2=input ('Enter slope of l i n e a r approx. f o r x2 viscous X force='); X Z3=k_x2; Z3=0.0; Z4=0.4*DIA*(L-2.0*LM)"3; X A l i n e a r approx. f o r the a d d i t i o n a l normal f o r c e s from X omega i s s t i l l i n c l u d e d , even though t h i s term has been X found t o have l i t t l e e f f e c t on the s o l u t i o n . k_x4=input(*Enter slope of the l i n e a r approx. f o r x4 ... viscous force='); Z5=k_x4; X norm, f o r c e = nfc*Q*A Q=dynam. press., A=cross X s e c t i o n a l area, nfc=normal f o r c e coef., Munk's X approx. X nf = (z7 * x3)-(z6 * x2) nfc=3.0; z6=nfc*V*area; z7=nfc*V*V*area; X p i t c h i n g moment = nf * I s l s = d i s t . of nf from cm, X t y p i c a l l y at the n o s e / c y l i n d e r j u n c t i o n X pm = (z8 * x3)-(z9 * x2) ls=2.0; z8=ls*z6; z9=ls*z7; d i s p C A l l done! !') ; Appendix B N u m e r i c Values U s e d i n Simulations Constant V = 5.1 ft/sec m = 4.66 slug b =b = 1.0 ft S = S = 0.5 ft h = 0.5 ft l = 3.5 ft l = 2.0 ft d = 0.75 ft C = 0.4 C = 1.2 x 2 1 2 2 m b d r d = 0.0416 ft W = 0.5 lb/ft Description stream velocity vehicle mass front and rear airfoil spans front and rear airfoil areas distance from towpoint to front foil mount distance from towpoint to rear foil mount distance from towpoint to center of mass body diameter body drag coefficient normal cable drag coefficient cable diameter cable weight per foot Polynomial Expressions for Body Normal Force and Moment Coefficients Cm = 20.26a - 0.9799a + 3.347a Cf4b - = -397.7a + 6.311a + 131.7a 3 2 3 2 116 Appendix C L i s t i n g of C a b l e M o d e l i n g P r o g r a m 2CABLE.F0R CALCULATES THE STATIC EQUILIBRIUM PROFILE OF A TWO CABLE TOWING SYSTEM. CABLE 1 IS NEUTRALLY BOUYANT AND CABLE TWO IS ARMOUR. A DEPRESSOR IS LOCATED AT THE JUNCTION OF THE TWO. V=VELOCITY,T=CABLE TENSION,TH=CABLE ANGLE,DT & DTH=INC. D=CABLE DIA.,W=CABLE WEIGHT, BD=BODY DIA.,WD=DEPR. WEIGHT X,Y=CABLE COOR.,DD=DEPR. DIA.,Ll=BUOYANT CABLE LENGTH L2=ARM0UR CABLE LENGTH,K=CL/ALPHA RATIO,ALPHA=FOIL ANGLE BDR=B0DY DRAG,NU=FRICTION RATIO USED TO CALC. LOADING FUNC. F=N0RMAL LOADING FUNC, G=TANGENTIAL LOADING FUNC. REAL T(103),DT,TH(103),DTH,V,DIA,W,BD REAL WD,X(103),Y(103) REAL DD,L1,L2,K,NU,CR,ALPHA,BDR,PI,DS1,DS2 REAL DS,D1,D2,WC REAL LIFT,CORD,SPAN,F,G,DDR,RX,RY,DIA1,DIA2,W1,W2 REAL NU1.NU2,APPROX,excur PI=3.1415926 0PEN(UNIT=1, FILE='2CABLE.DAT',STATUS='OLD') READ GENERAL DATA READ(1,*T V,DD,BD,WD,CR,K READ FOIL DATA READ(1,*) CORD,SPAN READ CABLE 1 DATA READ(1,*) Ll.Wl.DIAl.NUl READ CABLE 2 DATA READ(1,*J L2,W2,DIA2,NU2 CLOSE(1) PRINT *, 'V=',V,'DD=',DD,'BD=',BD,'WD=',WD,'CR=', c CR,'K=',K PRINT *, 'CORD=>,CORD,'SPAN='.SPAN PRINT *, 'Ll=',L1,'W1=',W1,'DIA1=',DIAi,'NUl=>,NU1 PRINT *, 'L2=',L2,'W2=',W2,'DIA2=',DIA2,'NU2=',NU2 PRINT *,'ARE THESE VALUES CORRECT?(Y/N)' READ(*,10) REPLY IF(REPLY .EQ. 'N') GO TO 40 PRINT *,'ENTER A VALUE OF ALPHA' READ *, ALPHA CALC. INC. LENGTHS DS1 AND DS2 DSl=Ll/50. DS2=L2/50. CALC. TOTAL LIFT FROM BOTH FOILS LIFT=2*K*ALPHA*C0RD*SPAN*V**2 CALC. TOW VEHICLE DRAG-SPHERE IN UNIFORM FLOW CD=0.4 BDR=0.4*PI*(BD*Vy**2/4.0 CALC. INITIAL TENSION AND CABLE ANGLE 117 Appendix C. Listing of Cable Modeling Program c c c c c c c c c c T(1)=(LIFT**2+BDR**2)**0.5 TH(1)=ATAN(-LIFT/BDR) PRINT *,T(1),TH(1) DO J=l,102 T(1)=TENSI0N AT VEHICLE, T(51)=TENSI0N AT HEAD OF CABLE #1 T(52)=TENSI0N AT DEPRESSOR, T(103)=TENSION AT SURFACE SIMILARLY FOR THETA CHOOSE CABLE PARAMETERS IFCJ .GT. 51) THEN DS=DS2 W=W2 NU=NU2 DIA=DIA2 ELSE DS=DS1 W=W1 NU=NU1 DIA=DIA1 ENDIF F=(1-NU)*SIN(TH(J))**2+ABS(NU*SIN(TH(J))) G=NU*COS(TH(J)) D1=CR*DIA*DS*F*V**2 D2=CR*DIA*DS*G*V**2 DT=D2+W*SIN(TH(J)) DTH=(D1-W*C0S(TH(J)))/(2.0*(T(J)+DT)) I F ( J .EQ. 51) THEN ... CALC. DEPRESSOR DRAG DDR=0.4*PI*(DD*V)**2/4.0 CALC. TENSION AND ANGLE AT DEPRESSOR RX=(DDR+TCJ)*COS(TH(J)))**2 RY=(WD+T(J)*SIN(TH(J)))**2 T(J+1)=(RX+RY)**.5 TH(J+1)=ATAN(RY/RX) ELSE T(J+1)=T(J)+DT TH(J+1)=TH(J)-DTH ENDIF ENDDO NOW WORK BACK DOWN THE CABLE TO CALC. THE (X,Y) CQOR. X(1)=0.0 Y(1)=0.0 DO 1=1,102 I F ( I .LT. 52) DS=DS2 I F ( I .EQ. 52) DS=0.0 I F ( I .GT. 52) DS=DS1 X(I+i)=X(I)+DS*C0S(TH(103-I)) Y(I+i)=Y(I)-DS*SIN(TH(103-I)) ENDDO OPEN(UNIT= i,FILE='XXXI.DAT',STATUS='OLD») OPEN(UNIT=2,FILE='YYY1.DAT',STATUS='OLD') OPEN(UNIT=3,FILE='THETA.DAT»,STATUS='OLD') DO 1=1,103 WRITECl,20) X(I) WRITE(2,20) Y(I) WRITE(3,20) TH(I) ENDDO CLOSE(l) CL0SE(2) CL0SE(3) OPEN(UNIT=1,FILE='2CABLE.OUT',STATUS=»OLD') WRITECl,*) 'V=',V,'DD=',DD,'BD=',BD,'WD=',WD,'CR=' c CR,'K=',K Appendix C. Listing of Cable Modeling Program WRITE(1,*) 'CORD=' ,CORD,'SPAN=',SPAN WRITE(1,*) 'Ll=',L1,'Wl=',W1,'DIA1=',DIA1,»NU1=»,NU1 WRITECl,*) 'L2=',L2,'W2 ',W2,'DIA2 ',DIA2,'NU2=',NU2 WRITECl,*) 'ANGLE OF ATTACK OF F0ILS='.ALPHA WRITECl,*) 'TENSION AT SURFACE=',T(103) TH C103)=TH C103)*180.O/PI WRITECl,*) 'ANGLE AT SURFACED,TH(103) WRITECl,*) 'TENSION AT VEHICLE-',T(1) THCD=THCD*180.0/PI WRITECl,*) 'ANGLE AT VEHICLE ',THC1) WRITECl,*) 'LIFT AND DRAG AT VEHICLE '.LIFT.BDR WRITECl,*) 'TENSION AT DEPRESSOR ',TC52) TH(52)=THC52)*180.0/PI WRITECl,*) 'ANGLE AT DEPRESSOR ',THC52) WRITECl,*) 'TENSION AT END OF CABLE 1=',TC51) WRITECl,*) 'DEPRESSOR POSITION CX.Y) ',XC52),YC52) WRITECl,*) 'TOW VEHICLE POSITION CX.Y) ',X(103), c YC103) APPROX=ASINCCYC52)-YC103))/L1) APPROX=APPROX*180.O/PI WRITECl,*) 'APPROX. CABLE ANGLE AT VEHICLE '.APPROX excur yC52)-yCl03) 3 3 3 3 3 3 3 3 3 3 w r i t e C i , * ) 'Y e x c u r s i o n of v e h i c l e =',excur 10 20 40 CLOSECD FORMATCA1) FORMATCF10.5) END Appendix D D e t a i l e d D e r i v a t i o n o f the C o n t r o l Objectives Sonar Transducer Vo <E ea -4oor FL Near Field Resolution Figure D.41: Side Scan Sonar Beam Geometry Figure D.41 shows the location of the sonar transducer on the towed vehicle, and the resulting sonar beam geometry. The darkened area of the sea floor is the width of the strip surveyed by one pulse of the transducers. The length of this strip varies, but is typically 50 meters, or 25 to either side of the vehicle. For this work it is assumed that the vehicle height from the sea floor h is always within the so called near field of the sonar. Here, the resolution is limited to the sonar transducer length IT. Using this near field resolution, it is possible to obtain simple expressions for the maximum allowed body rotation angle and rotation rate. The maximum rotation rate is based on ensuring that there is minimal overlap between two consecutive pulses. To do this, consider a point A located at the midpoint of each pulse, as in Figure D.41. The propagation rate of this point must not exceed V = Irfp, where f is the pulse c 120 p Appendix D. Detailed Derivation of the Control Objectives 121 la h - 4"— Target Length Figure D.42: Target Measurement Error Resulting From Body Rotations rate, which is assumed to be 10 per second. The actual propagation rate of point A is V A = V + ha 0 where a is the body rotation rate. The maximum rotation rate is that for which V = VA' Therefore, if h = 10, Vo = 1.55, and IT = 0.5, the maximum rotation rate is approximately 20 degrees/second. For this work, the maximum rotation angle is based on how accurately the size of observed sea floor objects can be determined. As a result of the sonar beam geometry, there is already an inherent minimum accuracy of IT- It is desirable to have the added error resulting from the body rotation angle not exceed 4/x- Therefore, to determine ctmax, consider only the error resulting from a, as depicted in Figure D.42. The measurement error for this situation is simply l — I, which is not to exceed 4Zj\ Using a small angle approximation, the geometry in Figure D.41 yields C a l a = l + 2ha so that Error = 2ha For h = 10, the maximum body rotation is approximately 6 degrees. This value for ctmax is clearly for the worst possible situation, where the body is rotated at a so as to introduce errors at both the leading and trailing edges of the target. m a x Appendix E W i n d T u n n e l Test D a t a Test Description Hemispherical JNose, front and rear foils removed, manometer ratio of: :10 Bodv lYont Rear Measurement s bl Angle Foil Foil Velocity Hi Angle Angle (mm Hg) -84 -30 -143 Jt 110 t -25 -60 -60 -102 | -99 110 -20 -15 -10 -5 5 10 15 20 25 30 - - -42 - -19 -8 9| 20 31 43 64 - - - -29 93 -42 -30 -18 -8 10 20 31 44 63 t -70 -42 -25 -12 13 27 42 68 102 -67 -43 -26 -11 12 26 43 68 100 t 14: 5 110 110 110 110 110 110 110 110 110 110 Test Description Hemispherical JNose, front and rear foils removed as indicated, manometer ratio of 1:10 Bodv Front Rear Measurements Angle Foil Foil ii Velocity Angle Angle (mm Hg) U 0 0 0 0 0 0 0 0 0 0 0 10 - -5 -10 - - 5 10 5 -5 -10 5 10 -5 -10 0 0 0 0 109 70 7 1 128 86 -110 -129 -1 11 J -13 i17 J1 -13 -I h 1 8 \ 018 129 132 87 86 -68 -68 -125 122 Ih 85 0 -3 I." 6 133 135 87 89 -74 -73 -134 3 118 118 118 118 140 140 140 140 125 125 125 125 Appendix E. Wind Tunnel Test Data Test Description Hemispherical Nose, rear foils removed, manometer ratio ?t.ip of "yin Measurements ront Rear Bodv Foil Velocity Foil ST ^2 Angle Angle Angle (mm Hg) 7TO" 7T2T 7Lo8" 1W 0 -5 -71 -89 110 0 -5 -46 -57 110 0 -10 -93 -117 110 5 -15 -52 -62 140 5 -10 -20 -23 110 5 32 110 -5 40 5 88 110 0 108 5 120 108 5 146 -22 108 -20 10 -23 34 108 -15 10 44 64 110 -10 10 81 112 110 -5 10 137 131 110 -0 10 162 Test Description" Hemispherical Nose, rear foils removed, manometer ?.tip of yin ratio ont Rear Measurements Bodv Foil "37 Velocity Foil Angle Angle Angle (mm Hg) ^2 3^ 11 15 "IDS" -20 15 59 79 108 -15 15 108 109 135 -10 15 110 130 163 -5 15 110 137 172 -30 20 108 108 155 -25 20 108 82 113 -20 20 108 128 169 -15 20 108 158 198 -10 20 110 177 225 -35 25 108 97 58 -30 25 108 107 153 -25 25 108 148 201 -20 25 108 171 230 -15 25 108 182 252 Appendix E. Wind Tunnel Test Data 124 Test Descri ption Hemispherical Nose, front foi Is removed, manometer r a t i o nj Bodv Angle 15 15 15 15 15 20 20 20 20 20 25 25 25 25 25 1:1 n Front Foil Angle - - Rear Foil Angle -25 -20 -15 -10 -5 -30 -25 -20 -15 -10 -35 -30 -25 -20 -15 ii 73 53 37 17 1 89 73 55 41 37 111 98 79 70 61 Measurements Velocity (mm Hg) 53 58 61 45 88 94 92 98 103 96 134 140 141 143 145 165 165 165 165 165 165 165 165 165 165 165 165 165 165 165 Test Description" Hemispherical JNose, front foils removed, manometer i,tip of 1^0 ratio ont Rear Measurements Bodv Foil Angle Foil ii i' Velocity Angle Angle (mm Hg) 2 ~5~ 5 5 5 5 10 10 10 10 10 -15 -10 -5 0 5 -20 -15 -10 -5 0 50 21 3 -21 -32 61 43 17 -7 -16 U 20 35 38 40 25 30 45 63 59 165 165 165 165 165 165 165 165 165 165 Test Description Hemispherical JNose, front and rear foils mounted, i n m P ^ p r ra.tjnLQ] o f 1 ;1fl mane :ont Rear Measurements Body Foil Angle St ST Velocity Foil Angle Angle ( m m .. H s) ~0~ 0 0 0 0 5 5 5 5 5 10 10 10 10 10 -5 0 5 10 -15 -10 -5 0 5 -20 -15 -10 -5 0 0 0 0 0 -5 -5 -5 -5 -5 -10 -10 -10 -10 -10 7T25" TT3T -68 -74 0 0 87 87 133 129 -137 -140 -68 -73 34 21 112 132 180 208 -56 -70 -27 -5 50 78 134 170 135 168 125 125 125 125 125 170 184 137 137 167 168 172 137 137 125 Appendix E. Wind Tunnel Test Data ~ Test Description Hemispherical Nose, front and rear foils mounted, ia,nornf»|er ra.t.joL of Q l 1:1Q man! Re ront Measurements Body Lear Foil Angle Velocity ST Si Foil Angle Angle (mm Hg) 15 ~72 ^37 T6T 15 -20 -15 48 167 95 15 -15 -15 97 173 146 15 -10 -15 183 244 174 15 -5 -15 166 218 137 20 -20 -20 163 240 164 20 -15 -20 180 288 172 20 -10 -20 220 300 174 125 Appendix F A C S L P r o g r a m Listings ********************* Degr2.csl ******************* PROGRAM SSCAN COUPLED PITCH HEAVE MOTION INITIAL " — D E F I N E ALL PRESET VARIABLES" LOGICAL MUNK , DWASH , LKUP constant s e t l = 1.0 , set2 = 1.0 , DP = 0.0 CONSTANT BET1 = 0.0 , BET2 = 0.0 , TC = 1.0 CONSTANT W = 150.0 , AMP = 0.0 , T l = 0 0 CONSTANT YIC = 0.0 , YDIC = 0.0 BET1 = O.C CONSTANT SPAN1 = 0.5 , C0RD1 = 0.33 CONSTANT SPAN2 = 0.5 , C0RD2 = 0.33 CONSTANT LD = 2.0 , LM - 2.0 T2 = 0.0 CONSTANT LI = 0.5 , L2 = 3.5 , LS = 2.0 CONSTANT PI = 3.1415926 , MK = 3.0 CINTERVAL CINT =0.1 CONSTANT L = 4.0 , DIA 0.75 , TY = 0.0 CONSTANT V = 5.1 , B l 0.0 , TZ = 0.0 CONSTANT TSTOP 60.0 Dl = 0.0 , BET2 = 0.0 CONSTANT ALPDIC = 0.0 ALPHIC =0.0 CONSTANT K l = 0 0 , K2 = 0.0 , K3 = 0.0 CONSTANT K4 = 0 0 , K5 = 0 , K6 = 0.0 CONSTANT K7 = 0 0 , K8 = 0 , K9 = 0.0 CONSTANT K10 = 0.0 , K l l 0.0 , K12 = 0.0 CONSTANT Gl = 0.0 G2 = 0 , G3 = 0.0 CONSTANT G6 = 0.0 G7 = 0 , G8 = 0.0 CONSTANT G9 = 0.0 G10 = 0.0 , G i l = 0.0 CONSTANT G12 = 0.0 G4 = 0.0 , G5 = 0.0 CONSTANT R1=0.0 , R2=0.0 CONSTANT R6 = 0 0 , MUNK , R4=0.0 , R5 = 0.0 = .TRUE. , DWASH = .FALSE. CONSTANT LKUP = .TRUE. DEGRAD = PI/180.0 RADDEG = 180.O/PI A l > (PI*DIA**2)/4.0 Q = V**2 LL = L/DIA — C A L C . MASS, INCLUDING ADDED MASS" M = (W/32.2)+((PI*L*2.0*DIA**2)/4.0) — C A L C . BODY DRAG BDR" BDR = 0.4*Al*q — C A L C . ASPECT RATIO CORRECTION FOR FOIL 1" ARI = (2.0*SPAN1)**2/(2.0*SPAN1*C0RD1) KC1 = 2.0*PI/(1.0+(2.0/ARl)) — C A L C . ASP. RATIO CORREC. FOR FOIL 2" AR2 = (2.0*SPAN2)**2/(2.0*SPAN2*CORD2) KC2 = 2.0*PI/(1.0+(2.0/AR2)) — C A L C . MOMENT QF INERTIA" IMO = W*(DIA**2?4.0+L**2/3.0)/(4*32.2) — C A L C . ADDED MOMENT OF INERTIA" IMADD = 0.8*IM0 126 Appendix F. ACSL Program Listings IM END $ "OF INITIAL" =.IMO+IMADD DYNAMIC DERIVATIVE " — C A L C . INDUCED ANGLE OF FLOW FOR EACH FOIL" YDF1 = YD+(LM-L1)*ALPHD YDF2 = YD+ (LM-L2)*ALPHD GAM1 = ATAN(YDE1/V) GAM2 = ATANCYDF2/V) " CALC. MAG. OF FLOW FOR EACH AIRFOIL" VI = SQRT(Q+YDF1**2) V2 = SQRT(Q+YDF2**2) " — C A L C . RAMP REF. INPUT" R3 = AMP*(RAMP(Tl)-RAMP(T2))/(T2-T1) " — C A L C . FOIL ROTATIONS FROM GAINS fc" KU1 =(K1*BET1*TC)+(K2*BET2*TC)+(K3*Y)+ . (K4*YD)+(K5*ALPH)+(K6*ALPHD) KU2 =(K7*BET1*TC)+(K8*BET2*TC)+(K9*Y)+ . (K10*YD)+(K11*ALPH)+(K12*ALPHD) KR1 =(K1*R1)+(K2*R2)+(K3*R3)+(K4*R4)+ ., (K5*R5)+(K6*R6) KR2 =(K7*R1)+(K8*R2)+(K9*R3)+(K10*R4)+ . (K11*R5)+(K12*R6) GR1 = (G1*R1)+(G2*R2)+(G3*R3)+(G4*R4)+ . (G5*R5)+(G6*R6) GR2 = (G7*R1)+(G8*R2)+(G9*R3)+(G10*R4)+ (G11*R5)+(G12*R6) ERR1 = KR1-KU1-GR1 ERR2 = KR2-KU2-GR2 " — C A L C . ACTUAL FOIL ROTATIONS USING ACTUATOR MODEL" BET1 = REALPL(TC,ERR1) BET2 = REALPL(TC.ERR2) " — C A L C . EFFECTIVE ANGLE OF ATTACK OF EACH FOIL" ATTAC1 = ALPH+BET1-GAM1 ATTAC2 = ALPH+BET2-GAM2 " — C A L C . COEFF. OF LIFT FOR FOILS" TABLE SLWO , 1, 6 / 0, 0.0873, 0.1745, 0.2618 ... , 0.3491, 0.4363, 5.111, 4.973, 4.165 ... , 3.575, 2.979, 3.781/ TABLE SLW , 1, 6 / 0, 0.0873, 0.1745, 0.2618 ... , 0.3491, 0.4363, 2.933, 3.684, 3.405 ... , 3.151, 2.649, 3.561/ SL0PE1=RSW(DWASH, SLW0(ATTAC3), KC1) SL0PE2=RSW(DWASH, SLW0(ATTAC3), KC2) DWRF=RSW(DWASH, (SLW0(ATTAC3) ... -SLW(ATTAC3))*ATTAC1, 0.0) CL1 = setl*SLOPEl*ATTACl CL2 = set2*(SL0PE2*ATTAC2-DWRF) " — C A L C . LIFT OF EACH FOIL" FL1 = 2.0*CL1*SPAN1*C0RD1*V1**2 FL2 = 2.0*CL2*SPAN2*CORD2*V2**2 " CALC. X AND vY COMP. OF LI LIFT" : . :: f FL1Y = FL1*C0S(GAM1) FL1X = FL1*SIN(GAM1) FL2Y = FL2*C0S(GAM2) FL2X = FL2*SIN(GAM2) " — C A L C . X COMP. OF CABLE TENSION" TX = BDR+FL1X+FL2X —CALC. CABLE ANGLE" TABLE ANGLE, 2, 6, 8 ... / 4.0 , 5.23 , 6.61 , 8.17 , 9.88 , 11.76 , 0 , 2 , 4 , 6 , 8 , 10 , 12 , 14 ... , 0.0 , 0.0 , 0.0 . 0.0 . 0.0 . 0.0 ... Appendix F. ACSL Program Listings , 3.83 , 3.74 , 3.65 , 3.52 , 3.39 , 3.31 ... , 10.12 ,9.0 ,8.25 , 7.74 ,7.39 ,7.10 , 18.78 ,15.88 ,14.1 ,12.88 ,12.11 ,11.5 ... , 28.88 ,24.17 ,21.1 ,19.0 ,17.5 ,16.4 , 39.37 ,33.25 ,29.0 ,25.8 ,23.5 ,22.0 , 48.93 ,42.15 ,37.0 ,33.0 ,30.0 ,27.8 , 55.96 ,49.5 ,44.2 ,39.8 ,36.6 ,33.9/ TH = RSW(LKUP,ANGLE(TX,Y)*DEGRAD,(Y+ALPH ... *LM)/20.0) " — C A L C . Y COMP. OF CABLE TENSION" TY = TX*TAN(TH) " — C A L C . MUNK'S NORMAL FORCE" GAM3 = ATANCYD/V) ATTAC3 = ALPH-GAM3 CN = RSW(MUNK,MK*ATTAC3, (20.26*ATTAC3**3)- ... (0.9799*ATTAC3**2)+(3.347*ATTAC3)) NF = CN*Q*A1 " — C A L C . MUNK'S PITCHING MOMENT" CM =-397.7*ATTAC3**3+6.311*ATTAC3**2+131.7*ATTAC3 " — N O T E CM IS IN IN-LB" MP = RSW(MUNK, LS*NF, CM*q*Al*DIA/12.0) " — C A L C . Y DIR. DAMPING FORCE" DAMP = DP*1.2*L*DIA*YD*ABS(YD) " — C A L C . PITCH DAMPING MOMENT" MOMENT =0.3*DIA*ALPHD*ABS(ALPHD)*((L-LM)**4+LM**4) " — C A L C . INDUCED NET FORCE FROM ABOVE" " EQUATION OF MOTION FOR Y" YDD = (FL1Y+FL2Y-TY-DAMP+NF)/M YD = INTEG(YDD.YDIC) Y = INTEG(YD,YIC) " EQUATION OF MOTION FOR ALPH" q i = FL1Y*(LM-L1) q2 = FL2Y*(LM-L2) Q3 = LM*TY ALPHDD = (qi+Q2-Q3-M0MENT+MP)/IM ALPHD = INTEG(ALPHDD,ALPDIC) ALPH = INTEG(ALPHD,ALPHIC) "—CONVERT TO DEGREES FOR PLOTTING" PALPH = RADDEG*ALPH PTAC1 = RADDEG*ATTAC1 PTAC2 = RADDEG*ATTAC2 PBET1 = RADDEG*BET1 PBET2 = RADDEG*BET2 PALPHD = RADDEG*ALPHD END $ "OF DERIVATIVE" II SPECIFY TERMINATION CONDITIONS" TERMT(T.GE.TSTOP; END $ "OF DYNAMIC" END $ "OF PROGRAM" I k * * * * * * * * * * * * * * * * * * * L _ o b s . c s l ********************* PROGRAM SSCAN COUPLED PITCH HEAVE MOTION INITIAL " — D E F I N E ALL PRESET VARIABLES" LOGICAL LKUP , DWASH , MUNK CONSTANT MUNK = .FALSE. , DWASH = .TRUE, constant s e t l = 1.0 , set2 = 1.0 , DP = 0.0 CONSTANT BET1 = 0.0 , BET2 = 0.0 , TC = 1.0 CONSTANT W = 150.0 , AMP = 0.0 , T l = 0.0 CONSTANT YIC = 0.0 , YDIC = 0.0 , BET1 =0.0 CONSTANT SPAN1 = 0.5 , C0RD1 =0.33 CONSTANT SPAN2 = 0.5 , C0RD2 =0.33 CONSTANT LD = 2.0 , LM = 2.0 , T2 = 0.0 ix F. ACSL Program Listings CONSTANT CONSTANT CINTERVAL CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT LI = 0.5 , L2 = 3.5 , LS = 2.0 PI = 3.1415926 , MK = 3.0 CINT =0.1 L = 4.0 , DIA 0.75 , TY = 0.0 Bl 0.0 , TZ = 0.0 V = 5.1 .0 , D D ll = 6.0 0.0 , BET2 TSTOP = 60.0 BET2 = 0.0 ALPHIC =0.0 ALPDIC = 0..0 K l = 0.0 , K2 = 0.0 , K3 = 0.0 K4 = 0.0 , K5 = 0.0 , K6 = 0.0 K7 = 0.0 , K8 = 0.0 , K9 = 0.0 K10 = 0.0 , K l l = 0.0 , K12 = 0.0 G3 = 0.0 Gl = 0.0 , G2 = 0.0 G8 = 0.0 G6 = 0.0 , G7 = 0.0 G i l = 0.0 G9 = 0.0 , G10 = 0.0 G5 = 0.0 G12 = 0.0 , G4 = 0.0 R1=0.0 , R2=0.0 , R4=0.0 , R5 = 0.0 R6 = 0.0 , 0B1IC=0.0 , OB2IC=0.0 OB3IC=0.0 , 0B4IC=0.0 , OB5IC=0.0 0B6IC=0.0 , A l l = -3.4818 , A12=1.1072 A13=-30.6417 , A14=-.1199 , A15=-31.85 A16=-1.6919 , A21=1.1072 , A22=-2.7272 A23=-27.0896 , A24=-.2244 , A25=31.7882 A26=.8441 , A31=0.0 , A32=0.0 A33=-64.2059 , A34=1.0 , A35=-.0542 A36=0.0 , A41=4.1846 , A42=4.1846 A43=-61.2225 , A44=-3.9744 , A45=7.1344 CONSTANT A46=0.0 , A51=0.0 , A52=0.0 CONSTANT A53=-.0542 , A54=0.0 , A55=-64.0525 CONSTANT A56=1.0 , A61=3.2789 , A62=-3.2789 CONSTANT A63=6.1095 , A64=-.7233 , A65=-47.7243 CONSTANT A66=-3.911 , KC13=29.6696 CONSTANT KC15=23.1549 , KC23=26.99 CONSTANT KC25=-28.7581 , KC33=64.2059 KC43=61.199 CONSTANT KC35=0.0542 KC53=.0542 CONSTANT KC45=13.0878 KC63=-6.1341 CONSTANT KC55=64.0525 BGR13=-0.0200 CONSTANT KC65=51.3639 KAPLHD =1.0 CONSTANT BGR23=0.0093 .TRUE , KC11 = 7.878 CONSTANT LKUP -0.0005 , KC31 = 0.1927 CONSTANT KC21 KC51 = 0.1478 3.4822 CONSTANT KC41 KC12 = -0.0005 2.6915 CONSTANT KC61 7.877 , KC32 = 0.1959 CONSTANT KC22 KC52 = -0.1943 3.2063 CONSTANT KC42 KC62 -3.1042 CONSTANT DEGRAD PI/180.0 RADDEG 180.O/PI A l = (PI*DIA**2)/4.0 Q = v**2 LL = L/DIA — C A L C . MASS, INCLUDING ADDED MASS" M = (W/32.2)+((PI*L*2.0*DIA**2)/4.0) — C A L C . BODY DRAG BDR" BDR = 0.4*Ai*q — C A L C . ASPECT RATIO CORRECTION FOR FOIL 1" ARI = (SPAN1**2)/(SPAN1*C0RD1) KC1 = 2.0*PI/(1.0+(2.0/ARl)) — C A L C . ASP. RATIO CORREC. FOR FOIL 2" AR2 = (SPAN2**2)7(SPAN2*C0RD2) KC2 = 2.0*PI/(1.0+(2.0/AR2)) CALC. MOMENT OF INERTIA" , IMO = W*(DIA**2/4.0+L**2/3.0)/(4*32.2) — C A L C . ADDED MOMENT OF INERTIA" IMADD = 0.8*IM0 Appendix F. ACSL Program Listings IM END $ "OF INITIAL" = IMO+IMADD DYNAMIC DERIVATIVE " CALC. INDUCED ANGLE OF FLOW FOR EACH FOIL" YDF1 = YD+(LM-L1)*ALPHD YDF2 = YD+(LM-L2)*ALPHD GAM1 = ATAN(YDF1/V) GAM2 = ATAN(YDF2/V) " — C A L C . MAG. OF FLOW FOR EACH AIRFOIL" VI = SQRT(q+YDFl**2) V2 = SQRT(Q+YDF2**2) " — C A L C . RAMP REF. INPUT" , ,, R3 = AMP*(RAMP(T1)-RAMP(T2))/(T2-T1) " — C A L C . FOIL ROTATIONS FROM GAINS K" KU1 =(K1*0B1)+(K2*0B2)+(K3*0B3)+ ... (K4*OB4)+(K5*0B5)+(K6*0B6) KU2 =(K7*0B1)+(K8*0B2)+(K9*0B3)+ ... (K10*0B4)+(K11*0B5)+(K12*0B6) KR1 =(K1*R1)+(K2*R2)+(K3*R3)+(K4*R4)+ ... (K5*R5)+(K6*R6) KR2 =(K7*Ri)+(K8*R2)+(K9*R3)+(K10*R4)+ ... (K11*R5)+(K12*R6) GR1 = (G1*R1)+(G2*R2)+(G3*R3)+(G4*R4)+ ... (G5*R5)+(G6*R6) GR2 = (G7*R1)+(G8*R2)+(G9*R3)+(G10*R4)+ .. (G11*R5)+(G12*R6) ERR1 = KR1-KU1-GR1 ERR2 = KR2-KU2-GR2 " — C A L C . ACTUAL FOIL ROTATIONS USING ACTUATOR MODEL BET1 = REALPL(TC.ERR1) BET2 = REALPL(TC.ERR2) " — C A L C . EFFECTIVE ANGLE OF ATTACK OF EACH FOIL" ATTAC1 = ALPH+BET1-GAM1 ATTAC2 = ALPH+BET2-GAM2 " — C A L C . COEFF. OF LIFT FOR FOILS" TABLE SLWO , 1, 6 / 0, 0.0873, 0.1745, 0.2618 ... , 0.3491, 0.4363, 5.111, 4.973, 4.165 ... , 3.151, 2.649, 3.561/ TABLE SLW , 1, 6 / 0, 0.0873, 0.1745, 0.2618 ... , 0.3491, 0.4363, 2.933, 3.684, 3.405 ... , 3.151, 2.649, 3.561/ SL0PE1=RSW(DWASH, SLW0(ATTAC3), KC1) SL0PE2=RSW(DWASH, SLW0(ATTAC3), KC2) DWRF=RSW(DWASH, (SLW0(ATTAC3) ... -SLW(ATTAC3))*ATTAC1, 0.0) CL1 = setl*SLOPEl*ATTACl CL2 = set2*(SL0PE2*ATTAC2-DWRF) v v 1 " — C A L C . LIFT OF EACH FOIL" FL1 = 2.0*CL1*SPAN1*CORD1*V1**2 FL2 = 2.0*CL2*SPAN2*C0RD2*V2**2 " — C A L C . X AND Y COMP, OF LIFT" FL1Y = FL1*C0S£GAM1) FL2Y = FL2*C0S(GAM2) FL1X = FL1*SIN(GAM1) FL2X = FL2*SIN(GAM2) " — C A L C . X COMP. OF CABLE TENSION" TX = BDR+FL1X+FL2X " — C A L C . CABLE ANGLE" TABLE ANGLE, 2, 6, 8 ... / 4.0 , 5.23 , 6.61 , 8.17 , 9.88 , 11.76 . , 0 , 2 , 4 , 6 , 8 , 10 , 12 , 14 . . . , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ... Appendix F. ACSL Program Listings " , 3.83 , 3.74 , 3.65 , 3.52 , 3.39 , 3.31 ... , 10.12 ,9.0 ,8.25 , 7.74 ,7.39 ,7.10 , 18.78 ,15.88 ,14.1 ,12.88 ,12.11 ,11.5 ... , 28.88 ,24.17 ,21.1 ,19.0 ,17.5 ,16.4 , 39.37 ,33.25 ,29.0 ,25.8 ,23.5 ,22.0 , 48.93 ,42.15 ,37.0 ,33.0 ,30.0 ,27.8 , 55.96 ,49.5 ,44.2 ,39.8 ,36.6 ,33.9/ CHOOSE TH BASED ON LOGICAL LHUP" TH = RSWCLKUP, ANGLE(TX,Y)*DEGRAD,(Y+ALPH ... *LM)/20.0 .BLE TENSION" CN = RSW(MUNK,MK*ATTAC3, (20.26*ATTAC3**3)- ... (0.9799*ATTAC3**2)+(3.347*ATTAC3)) NF = CN*Q*A1 "—CALC. MUNK'S PITCHING MOMENT" CM =-397.7*ATTAC3**3+6.311*ATTAC3**2+131.7*ATTAC "—NOTE CM IS IN IN-LB" MP = RSW(MUNK, LS*NF, CM*q*Al*DIA/12.0) "—CALC. PITCH DAMPING MOMENT" MOMENT =0.3*DIA*ALPHD*ABS(ALPHD)*((L-LM)**4+LM**4) II " I-ATP TMnnrirn MTTT trnnri? lrnnM Anrnnr" EQUATION OF MOTION FOR Y" YDD = (FL1Y+FL2Y-TY+NF)/M YD = INTEG(YDD.YDIC) Y = INTEG(YD.YIC) "—EQUATION OF MOTION FOR ALPH" Qi = FL1Y*(LM-L1) Q2 = FL2Y*(LM-L2) Q3 = LM*TY ALPHDD o (Q1+Q2-Q3-M0MENT+MP)/IM ALPHD = INTEG(ALPHDD,ALPDIC) ALPH = INTEG(ALPHD,ALPHIC) " CONVERT TO DEGREES FOR PLOTTING" PALPH = RADDEG*ALPH PTAC1 = RADDEG*ATTAC1 PTAC2 = RADDEG*ATTAC2 PBET1 = RADDEG*BET1 PBET2 = RADDEG*BET2 PALPHD = RADDEG*ALPHD P0B5 = RADDEG*0B5 P0B5D = RADDEG*0B5D P0B6 = RADDEG*0B6 "—CALC. OBSERVER USING INPUT A AND B MATRICES" 0B1D = (A11*0B1)+(A12*0B2)+(A13*0B3)+ ... (A14*0B4)+(A15*0B5)+(A16*0B6)+ ... (KC11*BET1*TC)+(KC12*BET2*TC)+ ... (KC13*Y)+(KC15*ALPH)- ... (BGR13*R3)+(K3*R3) 0B2D = (A21*0Bi)+(A22*0B2)+(A23*0B3)+ ... (A24*0B4)+(A25*0B5)+(A26*0B6)+ ... (KC21*BET1*TC)+(KC22*BET2*TC)+ ... (KC23*Y)+(KC25*ALPH)- ... (BGR23*R3)+(K9*R3) 0B3D = (A31*0B1)+(A32*0B2)+(A33*0B3)+ ... (A34*0B4)+(A35*0B5)+(A36*0B6)+ ... (KC31*BET1*TC)+(KC32*BET2*TC)+ ... (KC33*Y)+(KC35*ALPH) 0B4D = (A41*0B1)+(A42*0B2)+(A43*0B3)+ ... (A44*0B4)+(A45*0B5)+(A46*0B6)+ ... 131 Appendix F. ACSL Program Listings (KC41*BET1*TC)+(KC42*BET2*TC) (KC43*Y)+(KC45*ALPH) 0B5D = (A51*0B1)+(A52*0B2)+(A53*0B3)+ (A54*0B4)+(A55*0B5)+(A56*0B6) (KC51*BET1*TC)+(KC52*BET2*TC) (KC53*Y)+(KC55*ALPH) 0B6D = (A61*0B1)+(A62*0B2)+(A63*0B3)+ (A64*0B4)+(A65*0B5)+(A66*0B6) (KC61*BET1*TC)+(KC62*BET2*TC) (KC63*Y)+(KC65*ALPH) OBI = INTEG(OBID.OBIIC) 0B2 = INTEG(0B2D,0B2IC) 0B3 = INTEG(0B3D,0B3IC) 0B4 = INTEG(0B4D,0B4IC) 0B5 = INTEG(0B5D,0B5IC) 0B6 = INTEG(0B6D,0B6IC) END $ "OF DERIVATIVE" " SPECIFY TERMINATION CONDITIONS" TERMT(T.GE.TSTOP) END $ "OF DYNAMIC" END $ "OF PROGRAM" ******************** L_obs.dat ********************* SET TITLE="NON-LINEAR PITCH AND HEAVE MOTION" S TCWPRN =72,DIS=9 $" FORCE 3 COLUMN OUTPUT WIDTH" S LKUP = .TRUE. $".TRUE.=LOOKUP, .FALSE.=LINEAR" S DWASH = .TRUE. S MUNK = .FALSE. s t s t o p =20.0 S C0RD1 = 0.5 S SPANi =1.0 S C0RD2 = 0.5 S SPAN2 =1.0 S LI = 0.5 $"FOIL 1 POS. FROM T.P." S L2 = 3.5 $"FOIL 2 POSITION FROM TOW POINT" s dp = 0.0 $"viscous cross flow 0=no, l=yes" S TC = 2.0 $"TIME CONS. FOR FOIL ACTUATOR" S 0B3IC =1.0 S 0B5IC = 0.1745 s s e t l =1.0 S MK = 3.0 S KI = 3.3547 $"q(5,5)=60" S K2 = -1.4061 S K3 = 0.9440 S K4 = 0.0824 S K5 = 10.1240 S K6 = 1.3461 S K7 = -1.4061 S K8 = 2.9425 S K9 = 0.2756 S K10 = 0.2325 S K l l = -5.5942 S K12 = -0.9004 S Gl = 3.3547 S G2 = -1.4061 S G3 = -0.0200 S G4 = 0.0824 S G5 = 7.8516 S G6 = 1.3461 S G7 = -1.4061 S G8 = 2.9425 S G9 = 0.0093 S G10 = 0.2325 S G i l = 2.6904 S G12 = -0.9004 Appendix F. ACSL Program Listings KC13 KC23 KC33 KC43 KC53 KC63 KC15 KC25 KC35 KC45 KC55 KC65 All = A12 = A13 = A14 = A15 = A16 = A21 = A22 = S S S S S S S = 0.1927 = 0.1959 = 8.9518 = 5.5722 = 0.9666 = 0.0155 = 0.1478 = -0.1943 = 0.9666 = 16.9162 = 8.5757 = 2.2682 -11.7327 1.4066 -1.1367 -0.0824 -10.2717 -1.3461 1.4066 -11.3195 A23 = -0.4714 A24 = -0.2325 A25 = 5.7885 A26 = 0.9004 A31 = -0.1927 A32 = -0.1959 A33 = -8.9518 S A34 = 1.00 S A35 = -0.9666 S A36 = 0.0 S A41 = 2.5854 S A42 = 2.8613 S A43 = -5.5978 S A44 = -5.5118 S A45 = 11.1426 S A46 =0.0 S A51 = -0.1478 S A52 = 0.1943 S A53 = -0.9666 S A54 = 0.0 S A55 = -8.5757 S A56 =1.0 S A61 = 2.7261 S A62 = -2.3134 S A63 = -0.0460 S A64 = -0.8963 S A65 = 2.2420 S A66 = -6.4400 S BGR13 = -0.0200 S BGR23 = 0.0093 S T l = 0.0 $"TIME OF RAMP BEGINNING" S T2 = 10.0 $"TIME OF MAX. VALUE" S AMP = 10.0 $"MAX. AMP. OF RAMP" S CALPLT=.TRUE. S PRNPLT=.FALSE. PREPAR T,Y,YD,YDD,MOMENT,PALPHD,PALPH PREPAR PBET1,PBET2,PTAC1,PTAC2,R3 PREPAR NF,MP,OB i,0B2,0B3,0B4,0B5,0B6 PREPAR P0B5.P0B6.P0B5D START PLOT R3,Y,Palph $"PLOT PTAC1,PTAC2" PLOT YD,PALPHD PLOT P0B5.P0B6.P0B5D $"PLOT PBET1 PBET2" PRINT "NCIPRN"=i,T,0B3,0B4,P0B5,P0B6,OBI,0B2,PTAC1,PTAC2 STOP Bibliography I.H. Abbott and A . E . VonDoenhoff. Theory of Wing Sections. Dover, 1959. J . M . Abel. Cable interactions in a depth controlled submersible. J. Hydronautics, 6(2), 1972. G . F . Anderson. Tow cable loading functions. AIAA Journal, 5(2), 1967. E . Atraghji. A method for estimating the loading distribution on long slender bodies of revolution at high angles of attack in incompressible flow. In AGARD Conference Proceedings No. 204, 1976. H . O. Berteaux. Improvement of Intermediate Oceanographic Winches. Technical Report, Woods Hole Oceanographic Institution, 1985. R . J . Boncal. A Study of Model Based Maneuvering Controls for Autonomous Underwater Vehicles. Master's thesis, Naval Postgraduate School, Monterey, California, 1987. J.D. Burroughs and R . C . Benz. Computer support in the design and testing of undersea towed systems. Marine Technology, 1974. D . E . Calkins. Hydrodynamic analysis if a high-speed marine towed system. J. Hydronautics, 13(1), 1979. T . C . Cannon and J . Genin. Three-dimensional dynamical behaviour of a flexible towed cable. Aeronautical Quarterly, 1972. D . A . Chapman. Effects of ship motion on a neutrally-stable towed fish. Ocean Engineering, 9(3), 1982. Y . Choo and M . J . Casarella. Configuration of a towline attached to a vehicle moving in a circular path. J. Hydronautics, 6(1), 1972. T . H . Dawson. Dynamics of shallow towing cables. Ocean Engineering, 4, 1977. Jean-Guy Dessureault. The Lateral and Directional Stability of Batfish, an Underwater Towed Body with Wings. Master's thesis, Technical University of Nova Scotia, Halifax, Nova Scotia, 1976. [14] L . J . Dreher. Robust Rate Control System Designs for a Submersible. Master's thesis, M.I.T., Cambridge, Massachusetts, 1984. 134 Bibliography 135 [15] M . C . Eames. Steady-State Theory of Towing Cables. Technical Report, Royal Institution of Naval Architects, 1967. B. Friedland. Control System Design. McGraw-Hill, 1986. K . J . Harris. Automatic Control of a Submersible. Master's thesis, M.I.T., Cambridge, Massachusetts, 1984. N . E . Jeffrey. Influence of design features on underwater towed system stability. J. Hydronautics, 2(4), 1968. A . M . Kuethe and C . Chow. Foundations of Aerodynamics. John Wiley and Sons, third edition edition, 1976. R . J . Martin. Multivariable control of a submersible using lqg/ltr design methodology. In Proceedings of the 1986 American Control Conference, American Automatic Control Council, 1986. Mitchell and Gauthier, Assoc., Inc. Advanced Continuous Simulation Language (ACSL)- Users Guide. Mithcell and Gauthier, Associates, 1981. C . Moler and J . Little. PC-MATLAB Users Guide. The Mathworks, Inc., 1987. J.D. Mudie and W . D . Ivers. Simulation studies of the response of a deeply towed vehicle to various towing ship maneuvers. Ocean Engineering, 3, 1975. J . N . Newman. Marine Hydrodynamics. M I T Press, 1986. J . N . Nielson. Missile Aerodynamics. McGraw-Hill, 1960. B. Paul and A.I. Soler. Cable dynamics and optimum towing stategies for submersibles. MTS Journal, 6(2), 1972. K i m Saunders. Design and initial testing of a motion compensating winch system. In Oceans 85, 1985. L . Shupe and T . McGeer. A fundamental mathematical model of the longitudinal and lateral/directional dynamics of the dolphin type unmanned semi-submersible. In DREP/RRMC Military Robotic Applications Workshop, 1988.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The dynamics and control of an underwater towed vehicle
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The dynamics and control of an underwater towed vehicle Hopkin, David A. 1989
pdf
Page Metadata
Item Metadata
Title | The dynamics and control of an underwater towed vehicle |
Creator |
Hopkin, David A. |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | A linear, quadratic optimum control strategy is applied to a non-linear dynamic model of an underwater towed vehicle in an effort to minimize the pitching motion of the vehicle during the tracking of a reference input. The dynamic model is a two-dimensional, coupled pitch/heave model. The cable interactions with the vehicle are simplified to steady state, and formed into a multi-variable look-up table used in the non-linear model. The normal force resulting from the body of the vehicle is non-linear and consists of two components, an invisid slender-body theory component, and a separated crossflow component. In addition, the dynamic model includes the non-linear effect of the front airfoil's downwash acting upon the rear airfoils. Aerodynamic testing of a scaled vehicle provides the expressions for the non-linear normal body force and moment. These tests also verify the finite aspect ratio corrections for the airfoils, and the downwash effects of the front airfoils. The linear control strategy is based on linear, quadratic optimum control. Simulation results show that proper selection of the state and input weighting matrices result in minimizing the pitch angle of the vehicle to within the control objectives. In addition, simulations of various observer designs show how the tracking and attitude control varies with the selected measurements. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080681 |
URI | http://hdl.handle.net/2429/27889 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1989_A7 H67.pdf [ 7.18MB ]
- Metadata
- JSON: 831-1.0080681.json
- JSON-LD: 831-1.0080681-ld.json
- RDF/XML (Pretty): 831-1.0080681-rdf.xml
- RDF/JSON: 831-1.0080681-rdf.json
- Turtle: 831-1.0080681-turtle.txt
- N-Triples: 831-1.0080681-rdf-ntriples.txt
- Original Record: 831-1.0080681-source.json
- Full Text
- 831-1.0080681-fulltext.txt
- Citation
- 831-1.0080681.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080681/manifest