T H E DYNAMICS AND CONTROL OF A N UNDERWATER T O W E D VEHICLE By David A. Hopkin B. A. Sc. (Mechanical Engineering) Technical University of Nova Scotia A THESIS SUBM ITTED IN PARTIAL F U L F I L L M E N T O F T H E REQU IREMENTS F O R T H E D E G R E E O F M A S T E R O F A P P L I E D SC I E N C E in T H E F A C U L T Y O F G R A D U A T E STUDIES M E C H A N I C A L ENGINEER ING We accept this thesis as conforming to the required standard T H E UNIVERSITY O F BRITISH COLUMB IA April 1989 © David A . Hopkin, 1989 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of (Mec U.>. <Jp t^f-r«*.cy The University of British Columbia Vancouver, Canada Date DE-6 (2/88) 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 vehi-cle 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 con-sists 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 correc-tions for the airfoils, and the downwash effects of the front airfoils. The linear control strategy is based on linear, quadratic optimum control. Sim-ulation 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 Table of Contents Abstract ii List of Figures v i 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 Proposed Approach 6 2 Vehicle and Cable Model ing 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 W i n d Tunne 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 4 Trajectory Contro l 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 5 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 6 Simulation Results 68 6.1 Overview 68 6.2 State Feedback with Reference Input 69 6.2.1 Varying the L Q R Weightings 69 iv 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 7 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 List ing of M A T L A B Files 102 B Numeric Values Used in Simulations 116 C List ing of Cable Model ing Program 117 D Detailed Derivation of the Contro l Objectives 120 E W i n d Tunne l Test D a t a 122 F A C S L Program Listings 126 Bibl iography 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 106 52 4.19 Schematic for Process and Reference Input 53 vi 4.20 Schematic of Observer Design 55 4.21 Schematic of Composite System 56 5.22 Program Structure for SSCAN.CSL 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 an-gle 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 con-stant of 2 sec, body weighting of 60, downwash neglected 75 6.28 Linear(solid) and Non-linear(dashed) Results for Simulation 5, down-wash and non-linear body forces neglected 77 6.29 Linear(solid) and Non-linear(dashed) Results for Simulation 5, down-wash, non-linear body forces and cable angles neglected 78 6.30 Linear(solid) and Non-linear(dashed) Results for Simulation 7, time con-stant of 2 sec, body weighting of 1, with initial conditions 80 6.31 Linear(solid) and Non-linear(dashed) Results for Simulation 8, time con-stant 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 Ob-server Design with vertical position and front foil rotation measurements 84 6.34 Linear Process States(solid) and Observer States(dashed) for Second Ob-server Design, with vertical position and pitch angle measurements . . . 86 vii 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), Sec-ond Observer Design with vertical position and pitch angle measurements 89 6.37 Linear Process States(solid) and Observer States(dashed), Third Ob-server 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 Ax cross-sectional area at distance x from the towpoint ft2 Acf profile area subject to separated flow ft 2 ARi aspect ratio 6, span of airfoil ft. Cm normal force coefficient of body Cu airfoil lift coefficient Ccf crossflow drag coefficient of body Cd body drag coefficient CT normal cable drag coefficient Db body drag lb. D0 displaced volume of vehicle ft 3 d cable diameter ft. db diameter of body ft. Imth pitch and roll moments of inertia, including added moments of inertia Li hydrodynamic lift of airfoils lb. Lb hydrodynamic lift of body lb. lm distance from towpoint to center of mass ft. Ii distance from towpoint to the mounting point of the airfoils ft. M total mass of vehicle, including added masses slugs x Mo mass of the vehicle slugs Mb added mass of the vehicle (body and airfoils) slugs m o o lift coefficient slope for airfoil with infinite aspect ratio mf- 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 2 T cable tension lb. V,V0 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/ft3 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 Description of the Prob lem 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. All of this new devices, however, had at least one thing in common; the tasks they were performing were be-coming 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. Introduction 3 1.2 Current Technology and its Limitations 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 predomi-nant 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 sig-nals that are proportional to the measured ship motion, or more precisely the vertical Chapter 1. Introduction 4 Figure 1.1: Shipboard Motion Compensation Systems, reproduced from Berteaux[4] Chapter 1. Introduction 5 velocity of the ships' tow point. 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 com-pensation. 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. All 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. How-ever, 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 Fish-eries 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 A p p r o a c h 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 fol-lowing 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 de-scribed 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 hydro-dynamic 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 de-velopment. In addition, wind tunnel tests of a scaled model verified the values used for many of the hydrodynamic coefficients. Chapter 2 Vehicle and Cable Model ing 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 (AUV). The vehicle employs two pairs of control-lable 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 free-dom 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 fun-damental 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 modi-fications 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 14 2.2 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 ten-sion, 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 2d +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.5CrpdV2 where p is the fluid density, d is the cable diameter, V is the stream velocity, and Cr is a normal drag coefficient which is a function of the Reynolds Number and the stranding or roughness of the cable. For this work Cr is assumed to be 1.2, a typical value for a 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 17 2.3 Mathematical Model ing 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 P i t ch /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. L0 is the normal 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 + L2 cos j2 ~ T sin 9 + Lb cos a (2.3) T+ Hydrofoi 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 = CLiqiSi 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 19 ipi, and can be written as Cn = m,-0», where m,- is the slope of the lift coefficient curve 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 pV2 ft = ^ (2-6) where Vi = ^V02 + (y + (lm-li)ay (2.7) where lm 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(V2 + (y + (lm-li)aY) Li = m-i a + 8i- arctan f ^ - + ^ m (2.10) Chapter 2. Vehicle and Cable Modeling 20 B o d y Forces Lb and 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, Lp 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 Ax 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 Lp = 2.0 q0Ab(a-jb) (2.12) where 75 = arctan(y/Po) is the induced angle of attack for the body. The only con-tributing part of the body to Lp is the nose section, where dAb/dx ^ 0. Note that if the flow is assumed to be unseparated at the tail by the use of an aerodynamic tail section, J0L Axdx = 0 and no net lift is generated, and we have D'Alembert's paradox for invisid flow. The second component of Lb, LCJ, represents the normal force due to separated crossflow and is written as Lcf = Cefq0il>lAcf (2.13) where Ccf is the crossflow drag coefficient, and Acf is the profile area of the body subject to the separated crossflow. The profile area Acf is a function of body angle of 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 Lb = 2q0Ab (a - arctan + Ccfq0^Acf (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.5CdpAbVQ2 (2.15) 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. Cable 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 = — a (2.16) cos 9 where Tx is the net horizontal force. Since surge motion is assumed to be zero in this study, Tx is simply the sum of all the horizontal components of the forces shown in Figure 2.8, or Tx = 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,- + Lh sin a cos 9 A d d e d 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 M0 = O.8M0, 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. F i n a l F o r m of y Equat ion By substituting equations 2.14, 2.18, and 2.10, and the relationships for 7; and 7& into equation 2.3 we obtain rrii V A l • ( V r 0 2 - ^ • ( y - ^ ( / m - / l • ) A ) 2 ) , X « + « - — ( i + ^ ) ) x cos (arctan ^ + < ^ ^ j ) 0.6CdpAVf+ « + « , - a r c t a n U - + < i ^ L X + sin I arctan I — + —-'^ 1 ) V \Vo Vo )) (2qoA0 (a - arctan ( ^ r ) ) + Cc/SoV'b^c/^ sin a (2q0Ab (a - arctan + Cc/^oV'b^c/^ cos a F i n a l F o r m of a Equat ion tan0+ The simplest form of the ci equation is (2.19) M = Ima = Li(lm — li) cos 7J — Tlm sin 9 + L0lD cos a (2.20) Chapter 2. Vehicle and Cable Modeling 23 where (Zm — li), lb, and lm are the distances from the center of mass to the origin of the forces Li, Lb, and T respectively. The measures (lm — /,•) and lm 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(V02 + (y + (lm-li)a)2y a = — < m,(lm - U) x a + 6t -arctan ( ^ + cos I arctan I — + 0.5CdpAbVo2-r a + 6i - arctan U + V0 x y , (lm- U)a sin I arctan I ——h V \Vo Vo + > lm tan 0+ (2.21) (2q0Ab (a - arctan (^r)) + Cc/^oV'b^c/j sin a lb (2qoAb (a - arctan ^^-^j + Cc/go^Mc/^ c o s a where Im includes the added moment of inertia of the body and the airfoils. 2.3.2 Linear Coupled P i tch /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 x = y y a a (2.23) As shown in the previous section, these four variables completely describe the coupled pitch/heave motion. The input vector u is simply u = 62 (2.24) 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) | £ = 0 + d dx (Li cos 7,- — T sin 6 + Lb cos a) |s=0 (x — x) (2.25) Ima = (Li(lm - k) cos7,- - Tlm sin6 + Lbh cos a)\x=o + d -QZ (Li(lm - li) cos 7, — Tlm 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 where 0 1 0 0 0 0 Cyy Cyy Cya Cya B = Cyl 0 0 0 1 0 0 Cmy Cmy Cma Cma Cml Cm2 Cyy — Cyy — 20.0M -go (5,-m,- + 2.0 A) C y a — Cya — i f ( 5 , W i 9 0 - ^ + 2 - 0 q o A b , -Simiq0(lm - it) V0M Cyi — Cmy — Cmy Cma — Cma — Cmi = — Dblm 20.0/m -go V0Im 1 (Sirmilm - U) + 2.0Ablb) J- ySiTUiilm - h)qo - ^ -5,m,g 0 (Z m - Z,)2 + 2.0q0Ablb V0Im Si-miqQ(lm - U) (2.27) (2.28) (2.29) (2.30) (2.31) (2.32) (2.33) (2.34) (2.35) (2.36) (2.37) As a result of this linearization, the following nonlinear effects are lost: 1. The horizontal component of the cable tension Tx is now a constant, equal to 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 Cyy. 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 Lp. 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 0 1 0 0 -0.0256 -4.322 21.99 0 A = 0 0 0 1 -0.0305 -0.8963 4.51 -4.78 (2.38) for which the eigenvalues are [ -6.7 +0.18» - 0 . 1 8 » -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> ~ Cti4 + miln - m2ip2 = 0 (2.39) Chapter 2. Vehicle and Cable Modeling 27 -5 - 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 0 1 0 0 X = x + 0 mi — r a 2 u (2.40) with the state vector x and the input vector u as x = u Si S2 (2.41) 2.3.4 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 0 1 0 X = T X + 0 1 T . 0 1 u y = i 0 r 0 -T J (2.42) x 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 Tunnel 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 Ccf of the body lift, cannot. In ad-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 inter-actions 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 105. Therefore, tests are done at two values of 9ft, 1.0 x 105 and 1.53 x 105, to show the dependance of the measured data on 9ft. 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 Tunnel 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 .8 o JS 8 u | 13 5 5 4 3 2 1 0 -1 -2 -3 -4f -30 + - V = 12m/sec o - V = 18 m/sec x - Atraghji -20 20 -10 0 10 Body Rotation - Degrees Figure 3.11: Normal Force Coefficient for the Body 30 3.3 Test Results 3.3.1 B o d y N o r m a l Force Coefficient The normal force coefficient is defined as Cm = qoAb (3.43) 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 106 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 o a o U I 30 25 20 15 10 -5 --10 -15 + V =; 12rn/sec o - V = 18 m/sec -30 -20 20 30 -10 0 10 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 Pitching Moment 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 M CtAb = df, = body diameter (3.44) qoAbdb Again, the two sets of data represent two stream velocities. 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 Foil Rotation - Degrees -5 0 5 Foil Rotation - Degrees 10 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 Airfo 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 35 theoretical value for this slope, with an aspect ratio of 4, is 0.073 per degree. 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 Downwash 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 36 Degrees 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 37 Body Angle (Degrees) Measured Lift Coef-ficient Slope for the Rear Foils Theoretical Lift Co-efficient Slope for the Bear Foils 0 -0.0380 -0.0492 5 -0.0225 -0.0209 10 -0.0133 -0.0134 15 -0.0073 -0.0080 20 - -0.0058 25 - -0.0038 30 - -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 (L2)v = *T T L2 (3.45) V0lr 72 where T/Volr is the nondimensional vortex strength, £ 2 / 7 2 is the lift coefficient slope of the rear wing, and ix is the tail interference factor. The reference length lr is based on the rear wing dimensions and is chosen as lr = 27T&2. The vortex strength T can be expressed as R = J 0 " 1 ? ? 1 . , (3.46) A(yv - 0.5c4) where yv is the vertical vortex position as shown in Figure 3.15. Making these substi-tutions, and recalling AR+ = b 2/Si, equation 3.45 can be written in lift coefficient form as (m2)v = IT-—7^-7 n _ , v (3.47) IvAR^yy - 0.54) Implicit in the development of equation 3.47 are the following assumptions: Chapter 3. Wind Tunnel Tests 38 y ^ — r yv Figure 3.15: Assumed Paths of Wing Tip 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 Contro l 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 MIMO 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 MIMO 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 MIMO 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 MIMO 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 opti-mizes vehicle performance by directly weighting the importance of the errors in in-dividual 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 Contro l Objectives and 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 commer-cially 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 trans-ducers 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. An 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. An 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. An 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 This is the augmented system, where (4.48) x = Si S'2 y y a ct The first two states are r times the actual foil rotations, and result from the augmen-tation 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 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 Linear Contro l Design 4.3.1 Controllabil i ty 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. gain matrix G is chosen to minimize the quadratic performance integral to (4.49) Chapter 4. Trajectory Control 44 An 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 fol-lowing theorem is usually referred to as the controllability theorem, and applies to the more general time variant system. Controllabil i ty Theorem 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 e AtBB'e A,tdt (4.51) Jo 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 • • • Ak~*B ] (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) = fT eA'TC'CeArdT (4.53) Jo 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 JV = C CA OA*-1 (4.54) 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 B = 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 (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 105 1.71 x 104 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 observable1. However, since none of the above singular values is zero or close to zero (in the order of < 10 - 1 1 ) , 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 105 1.29 x 104 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 IEEE floating point 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 105, this means singular values in the order of 10 - 1 1 will result in a singular system Chapter 4. Trajectory Control 47 the singular values of the grammian become 8.17 x 105 1.29 x 104 1.26 7.22 x 10"1 6.22 x 10"3 8.14 x 10~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 105 1.32 x 104 1.49 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 r 1 0 0 0 0 0 0 0 1 0 0 0 and the singular values shown above. The second design uses vertical position and body rotation measurements, such that C = (4.58) C = 0 0 1 0 0 0 0 0 0 0 1 0 (4.59) and the singular values are 8.18 x 10s 1.29 x 104 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 0 1 0 (4.60) Chapter 4. Trajectory Control 48 u B o x C y Figure 4.16: Schematic of Regulator Feedback Law and the singular values are 8.18 x 105 1.29 x 104 1.27 1.00 7.43 x 10"3 7.14 x 10"4 4.3.2 State Feedback Design (4.61) A basic premise of linear state feedback design is that given a controllable, time-invariant system x = Ax + T5u y = 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 MIMO 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 ctn of the characteristic equation det(sl - A ) = sn + ctis"-1 + ... + <*„ and the desired eigenvalues of the system A i , . . . , A„. The eigenvalues can be used to determine the desired characteristic equation (s - - A 2) • • • (s - A„) = sn + ens"-1 +••• + «„ with which G can be expressed as G = [(QW)'] Ctl Ctr, (4.62) where Q is the controllability matrix and W = 1 a i ••• a„_ i 0 1 ••• Ctn-2 0 0 ••• 1 (4.63) If this method is considered for a MIMO system, the calculation of G is under-determined, with more gain values to solve for than eigenvalues. As previously men-tioned, 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 MIMO 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 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 (4.64) Voo = x'(t0)Mx(t0) where M satisfies the algebraic Riccati equation 0 = M A + A ' M - M B R - 1 B ' M + Q (4.65) Chapter 4. Trajectory Control 51 The resulting optimum gain matrix is expressed as G = R ^ B M (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 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 weight-ing 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 sta-ble 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 106. The gain matrix is determined for each weighting, and the resulting eigenvalues are plotted. Note that as R —> oo the eigenvalues of the closed loop system approach the open loop eigenvalues of x - (A - BG)x (4.67) y - Cx -0.5 -0.5 -2.41 -6.7 +0.18t -0.18t 4.3.3 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 106 Chapter 4. Trajectory Control 53 —e G o B . i o s X y 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 r a; r , where the error e is the difference between the state x and the reference input xr. The feedforward gain G r enables the controller to provide zero steady state error e with a nonzero input of — Grxr. Therefore, G 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 r is simply the inverse of the steady state of the process of 3.55. However, G r introduces a difficulty with the performance integral /o o [x'(r)Qx(r) + t i ' (r)R«(r)]rfr (4.68) If the input u is now not zero as a result of the contribution of G r « r , then the integral 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 - 1 B M (4.69) and GR = [C(A - BG)~LB]~1C(A - BG)~LA (4.70) With these gains defined, 3.55 can be re-written as x = A x + Bit where u = —Gx + ( G — G r ) x r or x = (A - B G ) x + B ( G - G r ) x r (4.71) y = Cx 4.3.4 Observer Design The controller designs of the previous two sections assume the state vector x is mea-sured and available for feedback. For this application, it is desirable to minimize the measurements to only the vertical position. The rate and angular measurements in-volve 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 y O y K u B o- x Figure 4.20: Schematic of Observer Design selected so as to minimize the error e = x — x. The result of this is A = A - K C B = B which allows 3.65 to be rewritten as x = (A - K C ) s + Bu + Ky (4.73) (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. An 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 ^ C H K K H i ~ B x. -, A «-Gr Xr r "Plant" B -o/>f A J J Figure 4.21: Schematic of Composite System must be expressed in the form X = AX + Bxr 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 r)cc r to yield x = Ax - BGx + B ( G - G r ) x r (4.75) In a similar manner, the expression in 4.72 is rewritten by making the substitutions for u = —Gx — ( G — G r ) x r and y = Cx to yield x = (A - B G - K C ) x + K C x + B ( G - G r ) x 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 57 A = B = B ( G — G r ) B ( G - G r ) (4.77) ( A - B G - K C ) K C - B G A 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 - Gr)xr (4.78) In a similar manner, since e = x — x, equation 4.76 becomes e = ( A - K C ) 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 pre-sented 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, nu-merical 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 er-rors 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 fol-lowing 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 SSCAN, 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 P R O G R A M I N I T I A L §. Statements to define initial conditions and constants E N D D Y N A M I C D E R I V A T I V E § Statements to describe the dynamic model E N D §. Statements executed every communication interval E N D E N D Figure 5.22: Program Structure for SSCAN.CSL 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 auto-sorting 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]. SSCAN.CSL 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 Tx of the cable tension is calculated. This value, along with 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 de-termined, 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 SSCAN.CSL through the A C S L processor before running the simulations. This procedure performs the auto-sorting of the state-ments, 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 com-plex elements. Fundamental math operators, elementary math functions, and more advanced scientific functions are all defined to operate on this fundamental matrix ob-ject. 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 Con-trol, System Identification, and Signal Processing toolboxes. The additional functions provided by the control toolbox are used extensively by the simulation program SIMU-L A T E . 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. Lqr_par Provides the option of preferentially weighting one of the states before per-forming 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, requir-ing only a time constant to be entered. Chapter 5. Numerical Simulation 65 Simulate Partlist S_state Heave Lqr. .par StJist Sim .opt Act.dyn Solve Observ New_st Pitch Newplot UJnput Refer Trans Attac Newplot P_vel P_pos P_omg P_ang Plotlist PJbi l l P_foil2 P.attacl P_attac2 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_input 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. Trans Performs the translation of the A C S L output file to defined M A T L A B variables. Newplot 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. Observ 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 , per-forms the linear uncoupled roll simulations. This program uses the script file ROLL_ST 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 per-form the pole placement for the system, and the simulation is performed with the Isim routine. Chapter 6 Simulation 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 ob-servable 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 with Reference Input 6.2.1 Vary ing 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 r . Note, the pitch angle steady state error for Chapter 6. Simulation Results sec. Pitch Angle .21 • 1 . 1 0 5 10 15 20 sec Attack Foil 1 sec .4I 1 L-i! , I 0 5 10 15 20 sec . 21 • 1 • 1 0 5 10 15 20 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 r is based on the linear model only. The 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 in-crease 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 sec Foil 1 Rotation l n Foil 2 Rotation ' ' 1 1 0 l 1 ' r sec sec Figure 6.25: Linear(solid) and Non-linear(dashed) Results for Simulation 2, Pitch anj Weighting is 60 Chapter 6. Simulation Results 73 6.2.2 Vary ing 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 Relative Importance of the Principle Non-l inear 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 calcula-tion 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 non-linear 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 . Attack Foil 2 sec sec Figure 6.26: Linear(solid) and Non-linear(dashed) Results for Simulation 3, Time Con-stant 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 con-stant 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 non-linear results of Figure 6.28. 6.3 State Feedback with Initial Conditions 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 shown in the following section. 6.4 Robustness of the State Feedback Design 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 80 sec sec Figure 6.30: Linear( solid) and Non-linear (dashed) Results for Simulation 7, time con-stant 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 con-stant 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 15 seconds Foil 1 Angle of Attack 5 10 15 seconds 20 20 10 15 seconds Foil 2 Angle of Attack 10 15 seconds 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 non-linear pitch angle of simulation 2. 6.5 Observer Design with 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 Vertical Velocity 5 10 15 seconds Body Rotation 5 10 15 seconds 20 20 5 10 15 seconds Body Rotation Rate 20 5 10 15 20 seconds Foil 1 Rotation Foil 2 Rotation 5 10 15 seconds Foil 1 Angle of Attack 5 10 15 seconds 20 20 5 10 15 seconds Foil 2 Angle of Attack 5 10 15 seconds 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 C = (6.81) The second observer design assumes that the vertical position and the pitch angle are available for measurement, such that 0 0 1 0 0 0 0 0 0 0 1 0 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 re-sponse 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 con-verged (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 Position Vertical Velocity 5 10 15 seconds Body Rotation 10 15 seconds 5 10 15 seconds Body Rotation Rate 5 10 15 seconds 20 20 Foil 1 Rotation Foil 2 Rotation 5 10 15 seconds Foil 1 Angle of Attack 5 10 15 seconds 20 20 5 10 15 seconds Foil 2 Angle of Attack 5 10 15 seconds 20 20 Figure 6.34: Linear Process States(solid) and Observer States(dashed) for Second Ob-server Design, with vertical position and pitch angle measurements Chapter 6. Simulation Results 87 Vertical Position Vertical Velocity 5 10 15 seconds Body Rotation 10 15 seconds T3 5 10 15 seconds Body Rotation Rate 20 Foil 1 Rotation Foil 2 Rotation 5 10 15 seconds Foil 1 Angle of Attack 20 5 10 15 20 seconds 5 10 15 20 seconds 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 0 0 1 0 0 0 0 0 0 1 0 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 non-linear 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 Vertical Velocity 5 10 15 seconds Body Rotation 20 -0.5 5 10 15 20 seconds Foil 1 Rotation Foil 2 Rotation 5 10 15 seconds Foil 1 Angle of Attack 20 5 10 15 seconds Foil 2 Angle of Attack 20 seconds seconds Figure 6.36: Non-linear Process States(solid) and Linear Process States(dashed), Sec-ond Observer Design with vertical position and pitch angle measurements Chapter 6. Simulation Results 90 Vertical Position 8 10 5 10 15 seconds Body Rotation -I L_ 10 15 seconds 20 20 Vertical Velocity 5 10 15 seconds Body Rotation Rate 5 10 15 seconds 20 20 Foil 1 Rotation Foil 2 Rotation 5 10 15 seconds 5 10 15 seconds Foil 2 Angle of Attack 5 10 15 seconds 20 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 Vertical Velocity 10 15 seconds 5 10 15 seconds Foil 1 Rotation Foil 2 Rotation 5 10 15 seconds Foil 1 Angle of Attack 5 10 15 seconds 20 20 10 15 seconds Foil 2 Angle of Attack 5 10 15 seconds 20 Figure 6.38: Non-linear Process States(solid) and Observer States(dashed), Third Ob-server 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 93 6.6 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 1 1 1 i 1 i i i i I 0 1 2 3 4 5 6 7 8 9 10 seconds Figure 6.40: Linear Uncoupled Roll Response, time constant of 2 sec, roll angle weight-ing of 8 Chapter 7 Conclusions 7.1 Contr ibut ion of 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 unrea-sonable, 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 strat-egy. 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 observ-ability 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 con-troller 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 2 , and makes extensive use of M A T L A B ' s programming capa-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 1Mitchell and Gauthier, Assoc., Inc. 2The Math Works, Inc. 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 Summary of the Simulation 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 actua-tor 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 feed-back 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 measur-ing 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 num-ber 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 100 7.3 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 List ing of M A T L A B Files '/,* Hop's Linearized Side Scan Sonar Simulation Program * ^ £ * * * * i | c g j c $ 3 l c * * * * * * 3 | c 3 f c * * * * * * $ ) | c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * fl 1- Calculate the basic state space representation using '/, the l i n e a r i z e d model '/, 2- Set LQR parameters '/, 3- Choose the various simulation options ( i e : root locus, '/, f u l l state simulation, include actuator dynamics, etc.) '/, 4- Calculate the input matrix and set i n i t i a l conditions */, 5- Perform simulation using f u l l state feedback '/, 6- Perform simulation using observer states '/, 7- Plot routine '/, 8- Load ACSL data (only must be done before plotting) */, 0- Quit •/********************* simulate.m************************** echo off */, This i s the main program for simulating the coupled pitch '/, and heave motion for the side scan sonar body. A simple '/, state feedback control i s presently used, employing LQR design techniques to determine the optimum state feedback '/, matrix. Actuator dynamics may be included. while 1 parts=[ 's_state' 'lqr.par' 'sim_opt' 'u_input' 'solve ' 'observ ' 'newplot' 'trans ' ] ; c l c help p a r t l i s t disp('Run each part i n sequence, followed by further ... choices'); qql=input('Enter the number of the part to execute:'); i f ((qql <= 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 '/, scr i p t f i l e that calculates the matrix values '/, of A,B,C and D i n the state space representation '/, of the coupled pitch 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; '/, Calc. the mass and the added mass of the body Ml=(W/32.2)+(pi*L*2.0*DIA*DIA/4.0); */, Calc. 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; Total 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)); '/, Calc. 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); '/, Calc. added moment of i n e r t i a of front f o i l s X mi.ff.cg = mi_ff + (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); '/, Calc. moment of i n e r t i a of rear 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 for added mom. of '/ i n e r t i a 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; '/, Originally a li n e a r approx. of the viscous force was used '/, This has now been replaced by a slender body normal force Appendix A. Listing of MATLAB Files '/, of 2 times alpha (alpha i s the induced angle of attack of '/, the body) This i s implimented with terms z6-z9 '/, k_x2=input('Enter slope of li n e a r approx. for x2 viscous '/. force='); '/. Z3=k_x2; Z3=0.0; Z4=0.4*DIA*(L-2.0*LM)"3; '/, A li n e a r approx. for the additional normal forces from */, omega i s s t i l l included, even though t h i s term has '/, been found to have l i t t l e effect on the solution. k_x4=input('Enter slope of the li n e a r approx. for x4 ... viscous force='); Z5=k_x4; '/, norm, force = nfc*Q*A Q=dynam. press., A=cross '/, sectional area, nfc=normal force coef., Munk's '/, approx. */. nf = (z7 * x3)-(z6 * x2) nfc=3.0; z6=nfc*V*area; z7=nfc*V*V*area; '/, pitching moment = nf * l s ls=dist. of nf from cm, '/, t y p i c a l l y at the nose/cylinder junction */, 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 ************************** '/, Script f i l e to calc. the coeff. i n the heave equation '/, of motion '/, to obtain a better approx. for cable angle, assume */, theta = y/20 instead of y/40 '/, Calc. c_yy due to vert, position c_yy=-BDR/20.0; '/, Calc. c_yv due to vert. v e l . c_yv=-((Zi*KCl/V)+(Z2*KC2/V)+Z3+z6); '/, Calc. c_ya due to pitch c_ya=((Z1*KC1)+(Z2*KC2)-(BDR*LM/20.0)+z7); '/, Calc. c_yq due to pitch rate c_yq=-((Z1*KC1*(LM-L1)/V)+(Z2*KC2*(LM-L2)/V)+Z4); '/, Form the 2x4 matrix from the heave equation heave_m=[0 MOO ;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 **************************** '/, sc r i p t f i l e that calculates the coeff. of the pitch '/, equation of motion '/, to obtain a better approx. f o r cable angle, assume '/, theta = y/20 instead of y/40 '/, Calc. the coeff. from v e r t i c a l position c_my=-BDR*LM/20.0; '/, Calc. c_mv due to vert. v e l . c_mv=-((Z1*KC1*(LM-L1)/V)+(Z2*KC2*(LM-L2)/V)+z8); Calc. c_ma due to pitch c_ma=(Zl*KCl*(LM-Ll))+(Z2*KC2*(LH-L2))-(BDR*LM*LM/20.0) ... +z9; '/, Calc. c_mq due to pitch rate c_mq=-((Z1*KC1*(LM-Li)~2/V)+(Z2*KC2*(LM-L2)"2/V)+Z5); '/, Calc. 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 pitch equation pitch_m=[0 0 0 IM ;c_my c_mv c_ma c_mq ]/IM; '/, Form the 2x2 input matrix from the pitch equation pitch_in=[0 0 ;c_mf c_mr ]/IM; '/,************************ lqr_par.m ************************ '/, Script f i l e to determine the LQR parameters '/, Input are: state 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 the state of Q to weight:'); st_wgt=input('Enter value of weighting:'); Q(st.lqr,st_lqr)=st_wgt; ^ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * l | l l | l * * * * * * * * ) | C ) | t ^ l l | l l | l 4 t * 4 c 4 C ! | c 4 C ) | C * * '/.* STATE LIST * '/.* 1- Front f o i l 2- Rear f o i l * '/,* 3- Heave 4- Heave rate * '/,* 5- Pitch angle 6- Pitch rate * 7 ! * * * $ $ $ $ $ $ * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * '/*********************** sim_opt.m ************************* '/, Script f i l e to determine the simulation options '/, Input are: Option f o r root locus plots or f u l l state Appendix A. Listing of MATLAB Files '/, simulation, option to solve Y as the input matrix '/, (must be selected to enable f o i l response p l o t s ) , '/, option to include actuator dynamics. '/, Option to include 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 off /, Script f i l e to include the actuator dynamics '/, F i r s t the model for the actuator i s made a_act=[-l/tc 0;0 -1/tc]; b_act=[l 0;0 1]; c_act=[l/tc 0;0 1/tc]; d_act=zeros(2); '/, Next t h i s model i s joined i n series to 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 labels as the o r i g i n a l '/, model y,********************** u_input.m ************************** '/, The following i s a scr 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 (sin) . '/, Given: a) the f i n a l time and in t e r v a l , '/, b) step amplitude and input state, '/, c) ramp amp., time of max. value and input state, '/, d) harm, amp., period, and input state. 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 for ramp input, 2 for 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. value:'); n2=t_rp/u_int; u_inc=a_rp/n2; u_in(l,l)=0.0; for i=l:n2. . , u_m(i+l, l)=u_in(i,l)+u_inc; end for i=n2+l:u_nn-l Appendix A. Listing of MATLAB Files u_in(i+l,l)=a_rp; end e l s e i f u_choice == 2, a_hm=input('enter the harm, amplitude:'); pd_hm=input('enter the period:'); 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 pick_ic=input('Inital conditions are zero. Okay?(y/n) ... : V s > ) ; i f p ick_ic == 'y', x_ic=[0 0 0 0 0 0]; else disp('Enter new i n i t i a l conditions matrix of the ... form'); x_ic=input('[betl bet2 y yd alpha alphad] : ' ) ; end y,************************* solve.m *********************** '/, Script f i l e which performs the solution part of the '/, analysis. '/, A LQR design i s used to f i n d the optimum state '/, feedback matrix. '/, Because of the cycling option, the c a l l to 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 routine, pack '/, perform linear-quadratic regulator design k=lqr(AA,BB,Q,R); '/, calc. new state 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'; '/, calc. angle of attack of f o i l s with attack.m attack newplot yt************************ new_st.m ************************ '/, Script f i l e to convert state space representation */, of the open loop system to a S.S. repr. with state '/, feedback gain k '/, i e : '/, X = A x + B u , Y = C x + D u '/, becomes '/, X = (A - B k) x + B k xref '/. or '/, X = Ac x + Be u '/, and Y w i l l be : '/, Y_option==0 Y = C x + D k xref Ac=AA-BB*k; Bc=BB*k; Appendix A. Listing of MATLAB Files 108 Cc=CC; Dc=DD*k; '/• Script f i l e to 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 rather '/, than 6x1 Bc=Bc(:,3); Cc=eye(6,6); Dc=zeros(6,l); ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ — ^ ^ ^ ^ 1 — VM ^ ^ ^ ^ ^ ^ ^ ^ ^ «V ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ I^f / - ^ T ^ *T* ^ ^ ^ ^ ^ *p 'T' T "I~ I »^ ^ ^ CLw JIW # III ^ T 1* ^ "X1 T* *P *p T* -I* T T ^ T ^ ^ A Script f i l e to calculate the angle of attack f o r each X f o i l 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************************ newplot.m ********************** X Plot script f i l e f o r selecting and plot t i n g either states X or inputs found from running cycle.m while 1 plots=[ 'p_pos ' 'p_vel ' 'p.ang ' 'p_omeg ' ' p . f o i l l ' 'p_foil2 ' 'p.attacl' 'p_attac2']; c l c help p l o t l i s t n=input('Enter number of plots per page:l,2,4 or 0 ... to q u i t ' ) ; i f ((n <= 0) I (n>8)) break end disp('Enter respective number of variable to p l o t ' ) ; nn=input('eg: 2 or [ 1 3 4 5] : ' ) ; nnn=input('Superimpose non-linear results? ... (y/n):','s'); i f n == 1, cl g plots=plots(nn,:); Appendix A. Listing of MATLAB 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 to make meta f i l e of current plot make_met end y,************************* p_pos.m ************************ '/, Script f i l e f o r pl o t t i n g the v e r t i c a l position x_cy(l,:) i f nnn == 'y', plot(t,x_cy(:,3),t,y_nl) t i t l e ( ' L and N-L Position'); 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. Velocity'); plot(t,x_cy(:,4)) t i t l e C V e r t i c a l v e l o c i t y ' ) ; end xlabel('sec'); y l a b e l ( ' f t / s e c ' ) ; y#***********************p.ang.m *************************** '/, Script f i l e to plot angular body rotation i f nnn =='y', plot(t,raddeg*x_cy(:,5),t,alph.nl) t i t l e ( ' L and N-L Body rotation'); 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 ************************* '/, Script f i l e to plot angular rotation rate 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)) title('Body Rotation Rate'); end xlabel('sec'); ylabel('degrees/sec'); yt************************* f o i l l . m ************************* '/, Script f i l e to plot f o i l 1 rotation and input to f o i l 1 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 Rotation'); 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,*************************** foil2.m *********************** '/, Script f i l e to plot 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 Rotation'); 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 ************************ '/, Script f i l e to plot 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 to plot angle of attack for 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'); 6 l S Q 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 plot ansl=input('Create a .met f i l e from current plot? y/n ... : \ ' s ' ) ; i f ansl == 'y', f_name=input('Input filename enclosed i n single ... quotes:','s'); eval(['meta',f_name]) ans2=input('Add to existing meta f i l e ? y/n :' , ' s ' ) ; i f ans2 == 'y', meta end end yt ********$*$********************************** V, * Hop's Plotting Routine * '/, 1- V e r t i c a l position 2- V e r t i c a l velocity '/, 3- Body rotation 4- Body rotation rate '/, 5- F o i l 1 response 6- F o i l 2 response '/, 7- Attack f o i l 1 8- Attack f o i l 2 •/,*********************** observ.m ************************** '/, Script f i l e which performs the solution part of the '/, analysis. A LQR design i s used to f i n d the optimum '/, state feedback matrix. The observer gains are then '/, c a l c , and a composite system matrix i s formed with '/, the new state vector being X=[xo x xr e f ] ' pack '/, perform linear-quadratic regulator design k=lqr(AA,BB,Q,R); '/, calc. new state space rep. with new_st.m new_st refer /, calc. observer gains using LQR design with weighting '/, of 70. The states 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 112 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 states i n the NL variables, which may '/, then be plotted against the actual, states by asking to '/, superimpose the non-linear results. However, loading '/, ACSL results after running t h i s simulation w i l l •/, overright the observer states. 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'; '/, calc. angle of attack of f o i l s with attack.m attack newplot y,*********************** trans.m ************************** '/, Script f i l e which w i l l convert a previously edited ACSL '/, data f i l e into separate matlab variables. '/, The edited f i l e must be i n the following form: •/.coll col2 col3 col4 col5 col6 col7 col8 col9 collO */, l i n e #, t ,y , yd , alph.alphd, b e t l , b e t l , atacl,atac2 '/, The matlab variables w i l l be the above names with the '/, extension _ n l (eg: y_nl ). disp('Has the ACSL data been previously converted to 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 to ... for020.mat.'); n=input('Okay? (y/n) : ','s'); i f n == 'n', disp('You may do so now...'); disp('Type !, your DOS copy command, [re t ] , ... [CTRL-Z] :') dispCeg: !copy [path]nlsim2.mat for020.mat ... [ret] [CTRL-Z] :') keyboard Appendix A. Listing of MATLAB Files end else !\matlab\translate end load 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); clear for020 '/,********************** uc_roll.m ************************ '/, Script f i l e to calc. the coeff. i n the uncoupled r o l l 7, equation of motion '/, Calculate the basic body parameters i n roll.in.m r o l l _ i n '/, Calc. 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; */, Calc. c _ l r cross derivative due to yaw '/, c _ l r w i l l be a linear approx. k_roll=SPANl~4*KCl*C0RDl; R_fit=0:0.002:.07; c _ r o l l = p o l y f i t ( R _ f i t , k _ r o l l * R _ f i t . ~ 2 , 1 ) ; c _ l r = c _ r o l l ( l , l ) ; */, Calc. c_lp 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; */, Calc. c _ l f l and c _ l f r from l e f t and right front f o i l s c.lf=KC1*Z1*(0.5*(SPAN1+DIA)); '/, Form the 2x2 matrix contributed by the r o l l equation aa=[ 0 1.1; 0 -c_lp ] / I _ l ; '/, Form the 2x1 input matrix from the r o l l equation bb=[0 ; c _ l f ] / I _ l ; 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 114 '/, 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'); '/,************************** roll_in.m ********************* '/, Script f i l e to calculate the basic body parameters '/, used by uc.roll.m to 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; '/, Calc. the mass and the added mass of the body Ml=(W/32.2)+(pi*L*2.0*DIA*DIA/4.0); '/, Calc. 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; Total 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)); '/, Calc. 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); '/, Calc. added moment of i n e r t i a of front f o i l s Appendix A . Listing of MATLAB Files '/, mi.ff.cg = mi.ff + (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 Calc. moment of i n e r t i a of rear 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 for added mom. of X i n e r t i a 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; '/, Originally a linear approx. of the viscous force was '/, used. This has now been replaced by a slender body '/, normal force of 2 times alpha (alpha i s the induced '/, angle of attack of the body) This i s implimented with */, terms z6-z9 '/, k_x2=input ('Enter slope of li 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 lin e a r approx. f o r the additional normal forces from X omega i s s t i l l included, even though t h i s term has been X found to have l i t t l e effect on the solution. k_x4=input(*Enter slope of the lin e a r approx. f o r x4 ... viscous force='); Z5=k_x4; X norm, force = nfc*Q*A Q=dynam. press., A=cross X sectional area, nfc=normal force coef., Munk's X approx. X nf = (z7 * x3)-(z6 * x2) nfc=3.0; z6=nfc*V*area; z7=nfc*V*V*area; X pitching moment = nf * Is ls=dist. of nf from cm, X t y p i c a l l y at the nose/cylinder junction X pm = (z8 * x3)-(z9 * x2) ls=2.0; z8=ls*z6; z9=ls*z7; d i s p C A l l done! !') ; A p p e n d i x B Numeric Values Used in Simulations Constant Description V = 5.1 ft/sec stream velocity m = 4.66 slug vehicle mass bx=b2 = 1.0 ft front and rear airfoil spans S1 = S2 = 0.5 ft front and rear airfoil areas h = 0.5 ft distance from towpoint to front foil mount l2 = 3.5 ft distance from towpoint to rear foil mount lm = 2.0 ft distance from towpoint to center of mass db = 0.75 ft body diameter Cd = 0.4 body drag coefficient Cr = 1.2 normal cable drag coefficient d = 0.0416 ft cable diameter W = 0.5 lb/ft cable weight per foot Polynomial Expressions for Body Normal Force and Moment Coefficients Cm = 20.26a3 - 0.9799a2 + 3.347a Cf4b -= -397.7a 3 + 6.311a2 + 131.7a 116 Appendix C List ing of Cable Model ing Program 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 T(1)=(LIFT**2+BDR**2)**0.5 TH(1)=ATAN(-LIFT/BDR) PRINT *,T(1),TH(1) DO J=l,102 c T(1)=TENSI0N AT VEHICLE, T(51)=TENSI0N AT HEAD OF c CABLE #1 c T(52)=TENSI0N AT DEPRESSOR, T(103)=TENSION AT c SURFACE c SIMILARLY FOR THETA c 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)) IF(J .EQ. 51) THEN c ... CALC. DEPRESSOR DRAG DDR=0.4*PI*(DD*V)**2/4.0 c 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) ENDDO ELSE ENDIF T(J+1)=T(J)+DT TH(J+1)=TH(J)-DTH c NOW WORK BACK DOWN THE CABLE TO CALC. THE (X,Y) c CQOR. X(1)=0.0 Y(1)=0.0 DO 1=1,102 IF(I .LT. 52) DS=DS2 IF(I .EQ. 52) DS=0.0 IF(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,'W23',W2,'DIA23',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 VEHICLE3',THC1) WRITECl,*) 'LIFT AND DRAG AT VEHICLE3'.LIFT.BDR WRITECl,*) 'TENSION AT DEPRESSOR3',TC52) TH(52)=THC52)*180.0/PI WRITECl,*) 'ANGLE AT DEPRESSOR3',THC52) WRITECl,*) 'TENSION AT END OF CABLE 1=',TC51) WRITECl,*) 'DEPRESSOR POSITION CX.Y)3',XC52),YC52) WRITECl,*) 'TOW VEHICLE POSITION CX.Y)3',X(103), c YC103) APPROX=ASINCCYC52)-YC103))/L1) APPROX=APPROX*180.O/PI WRITECl,*) 'APPROX. CABLE ANGLE AT VEHICLE3'.APPROX excur 3yC52)-yCl03) writeCi,*) 'Y excursion of ve h i c l e =',excur CLOSECD 10 FORMATCA1) 20 FORMATCF10.5) 40 END Appendix D Detailed Derivation of the Contro l Objectives Sonar Transducer Vo <E -4 FL ea oor 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 Vc = Irfp, where fp is the pulse 120 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 VA = V0 + ha where a is the body rotation rate. The maximum rotation rate is that for which VC = 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 mea-surement error for this situation is simply la — I, which is not to exceed 4Zj\ Using a small angle approximation, the geometry in Figure D.41 yields la = 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 m a x so as to introduce errors at both the leading and trailing edges of the target. Appendix E W i n d Tunnel Test D a t a Test Description Hemispherical JNose, front and rear foils re-moved, manometer ratio of: :10 Bodv lYont Rear Measurement s Angle Foil Foil b l Hi Velocity Angle Angle (mm Hg) -30 - - -84 t -143 Jt 110 -25 - - -60 -60 -102 | -99 110 -20 - - -42 -42 -70 -67 110 -15 - - -29 -30 -42 -43 110 -10 - - -19 -18 -25 -26 110 -5 - - -8 -8 -12 -11 110 5 - - 9| 10 13 12 110 10 - - 20 20 27 26 110 15 - - 31 31 42 43 110 20 - - 43 44 68 68 110 25 - - 64 63 102 100 110 30 - - 93 t 14: 5 t 110 Test Description Hemispherical JNose, front and rear foils removed as indicated, manometer ratio of 1:10 Bodv Front Rear Measurements Angle Foil Foil i i Velocity Angle Angle (mm Hg) U 10 - 109 128 118 0 5 70 J 7 1 86 I 85 118 0 -5 - -I h h 118 0 -10 - -110 -129 118 0 - 5 -13 J 1 -13 -1 11 0 140 0 - 10 i 17 140 0 -5 1 8 \ 18 -3 I."3 140 0 - -10 0 6 140 0 10 0 129 132 133 135 125 0 5 0 87 86 87 89 125 0 -5 0 -68 -68 -74 -73 125 0 -10 0 -125 -1 34 125 122 Appendix E. Wind Tunnel Test Data Hemispherical ?t.ip of "yin Test Description Nose, rear foils removed, manometer ra io Bodv Angle Rear Foil Angle Measurements 7TO" -5 0 0 5 5 5 5 5 10 10 10 10 10 ront Foil Angle ST 7T2T -71 -46 -93 -52 -20 32 88 120 -22 34 64 112 131 ^2 7Lo8" -89 -57 -117 -62 -23 40 108 146 -23 44 81 137 162 Velocity (mm Hg) 0 -5 -10 -15 -10 -5 0 5 -20 -15 -10 -5 -0 1W 110 110 110 140 110 110 110 108 108 108 110 110 110 Test Description" Hemispherical Nose, rear foils removed, manometer ?.tip of yin ratio Bodv Angle ont Foil Angle Rear Foil Angle "37 Measurements Velocity (mm Hg) "IDS" 108 108 110 110 108 108 108 108 110 108 108 108 108 108 15 15 15 15 15 20 20 20 20 20 25 25 25 25 25 2^ -20 -15 -10 -5 -30 -25 -20 -15 -10 -35 -30 -25 -20 -15 3^ 59 109 130 137 108 82 128 158 177 58 107 148 171 182 11 79 135 163 172 155 113 169 198 225 97 153 201 230 252 Appendix E. Wind Tunnel Test Data 124 Test Descri ption Hemispherical Nose, front foi Is removed, manometer r a t i o nj 1:1 n Bodv Front Rear Measurements Angle Foil Foil i i Velocity Angle Angle (mm Hg) 15 - -25 73 53 165 15 - -20 53 58 165 15 - -15 37 61 165 15 - -10 17 45 165 15 - -5 1 88 165 20 - -30 89 94 165 20 -25 73 92 165 20 - -20 55 98 165 20 - -15 41 103 165 20 - -10 37 96 165 25 - -35 111 134 165 25 - -30 98 140 165 25 - -25 79 141 165 25 - -20 70 143 165 25 - -15 61 145 165 Test Description" JNose, front foils removed, manometer Hemispherical i,tip of 1^ 0 rat o Bodv Angle ~5~ 5 5 5 5 10 10 10 10 10 ont Foil Angle Rear Measurements Foil i i i ' 2 Velocity Angle (mm Hg) -15 50 U 165 -10 21 20 165 -5 3 35 165 0 -21 38 165 5 -32 40 165 -20 61 25 165 -15 43 30 165 -10 17 45 165 -5 -7 63 165 0 -16 59 165 Test Description Hemispherical JNose, front and rear foils mounted, inmP^pr ra.tjn o f 1 ;1fl mane Body Angle LQ] Rear Foil Angle 0 0 0 0 -5 -5 -5 -5 -5 -10 -10 -10 -10 -10 Measurements S t ~0~ 0 0 0 0 5 5 5 5 5 10 10 10 10 10 :ont Foil Angle ST 7T25" -68 0 87 129 -137 -73 21 112 180 -70 -27 50 134 135 TT3T -74 0 87 133 -140 -68 34 132 208 -56 -5 78 170 168 Velocity (m m..Hs) -5 0 5 10 -15 -10 -5 0 5 -20 -15 -10 -5 0 125 125 125 125 125 170 184 137 137 167 168 172 137 137 125 Appendix E. Wind Tunnel Test Data 125 ~ Test Description Hemispherical Nose, front and rear foils mounted, ia,nornf»|er ra.t.jo of 1:1Q m n! Body Angle L Q l Re ront Foil Angle Lear Foil Angle ST Measurements ~72 -20 -15 -10 -5 -20 -15 -10 -15 -15 -15 -15 -20 -20 -20 Si Velocity (mm Hg) 15 15 15 15 15 20 20 20 ^37 48 97 183 166 163 180 220 95 146 244 218 240 288 300 T 6 T 167 173 174 137 164 172 174 Appendix F A C S L Program Listings ********************* Degr2.csl ******************* PROGRAM SSCAN COUPLED PITCH HEAVE MOTION INITIAL "—DEFINE ALL PRESET VARIABLES" LOGICAL MUNK , DWASH , LKUP constant s e t l = 1.0 , set2 = 1.0 BET1 = 0.0 , BET2 = 0.0 W = 150.0 , AMP = 0.0 , YIC = 0.0 , YDIC = 0.0 SPAN1 = 0.5 , C0RD1 = SPAN2 = 0.5 , C0RD2 = LD = 2.0 , LM - 2.0 LI = 0.5 , L2 = 3.5 , PI = 3.1415926 , MK = CINT =0.1 L = 4.0 , DIA V = 5.1 , Bl CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CINTERVAL CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT , DP = , TC = T l = 0 BET1 = 33 33 T2 = 0.0 LS = 2.0 3.0 0.0 1.0 0 O.C 0. 0. TSTOP ALPDIC Kl = 0 K4 = 0 K7 = 0 K10 = 60.0 = 0. 0 , 0 , 0 , 0.0 , Gl = 0.0 G6 = 0.0 G9 = 0.0 G12 = 0.0 R1=0.0 R6 = 0 LKUP = DEGRAD = RADDEG = Al > Q = 0 K2 = K5 = K8 = K l l G2 = G7 = G10 = G4 = , R2=0.0 0 , MUNK .TRUE. PI/180.0 180.O/PI (PI*DIA**2)/4.0 V**2 0.75 , TY = 0.0 0.0 , TZ = 0.0 Dl = 0.0 , BET2 = 0.0 ALPHIC =0.0 0.0 , K3 = 0.0 0 , K6 = 0.0 0 , K9 = 0.0 0.0 , K12 = 0.0 0 , G3 = 0.0 0 , G8 = 0.0 0.0 , G i l = 0.0 0.0 , G5 = 0.0 , R4=0.0 , R5 = 0.0 = .TRUE. , DWASH = .FALSE. 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 =.IMO+IMADD END $ "OF INITIAL" DYNAMIC DERIVATIVE "—CALC. 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) "—CALC. RAMP REF. INPUT" R3 = AMP*(RAMP(Tl)-RAMP(T2))/(T2-T1) "—CALC. 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 "—CALC. ACTUAL FOIL ROTATIONS USING ACTUATOR MODEL" BET1 = REALPL(TC,ERR1) BET2 = REALPL(TC.ERR2) "—CALC. EFFECTIVE ANGLE OF ATTACK OF EACH FOIL" ATTAC1 = ALPH+BET1-GAM1 ATTAC2 = ALPH+BET2-GAM2 "—CALC. 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) "—CALC. LIFT OF EACH FOIL" FL1 = 2.0*CL1*SPAN1*C0RD1*V1**2 FL2 = 2.0*CL2*SPAN2*CORD2*V2**2 : . :: v f LI FL1Y = FL1*C0S(GAM1) FL1X = FL1*SIN(GAM1) FL2Y = FL2*C0S(GAM2) FL2X = FL2*SIN(GAM2) "—CALC. 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 ... " CALC. X AND Y COMP. OF LIFT" 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) "—CALC. Y COMP. OF CABLE TENSION" TY = TX*TAN(TH) "—CALC. 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 "—CALC. MUNK'S PITCHING MOMENT" CM =-397.7*ATTAC3**3+6.311*ATTAC3**2+131.7*ATTAC3 "—NOTE CM IS IN IN-LB" MP = RSW(MUNK, LS*NF, CM*q*Al*DIA/12.0) "—CALC. Y DIR. DAMPING FORCE" DAMP = DP*1.2*L*DIA*YD*ABS(YD) "—CALC. PITCH DAMPING MOMENT" MOMENT =0.3*DIA*ALPHD*ABS(ALPHD)*((L-LM)**4+LM**4) "—CALC. 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" qi = 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_obs.csl ********************* PROGRAM SSCAN COUPLED PITCH HEAVE MOTION INITIAL "—DEFINE 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 Bl .0 .0 K2 = 0.0 , K5 = 0.0 , K8 = 0.0 , , K l l = 0.0 G2 = G7 = 0.0 0.0 LI = 0.5 , L2 = 3.5 , LS = 2.0 PI = 3.1415926 , MK = 3.0 CINT =0.1 L = 4.0 , DIA V = 5.1 TSTOP = 60 , Dl = 6.0 , BET2 = 0.0 ALPDIC = 0. Kl = 0.0 , K4 = 0.0 , K7 = 0.0 , K10 = 0.0 Gl = 0.0 , G6 = 0.0 , G9 = 0.0 , G10 = 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 , All= -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 0.75 , TY = 0.0 0.0 , TZ = 0.0 l 0ALPHIC =0.0 K3 = 0.0 K6 = 0.0 K9 = 0.0 , K12 = 0.0 G3 = 0.0 G8 = 0.0 G i l = 0.0 G5 = 0.0 CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT CONSTANT A46=0.0 , A51=0.0 , A52=0.0 A53=-.0542 , A54=0.0 , A55=-64.0525 A56=1.0 , A61=3.2789 , A62=-3.2789 A63=6.1095 , A64=-.7233 , A65=-47.7243 A66=-3.911 , KC13=29.6696 KC15=23.1549 , KC23=26.99 KC25=-28.7581 , KC33=64.2059 KC43=61.199 KC53=.0542 KC63=-6.1341 BGR13=-0.0200 KAPLHD =1.0 , KC11 = 7.878 , KC31 = 0.1927 KC51 = 0.1478 KC12 = -0.0005 KC32 = 0.1959 KC52 = -0.1943 KC35=0.0542 KC45=13.0878 KC55=64.0525 KC65=51.3639 BGR23=0.0093 LKUP KC21 KC41 KC61 KC22 KC42 KC62 DEGRAD RADDEG .TRUE -0.0005 3.4822 2.6915 7.877 , 3.2063 -3.1042 PI/180.0 180.O/PI Al = (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 = IMO+IMADD END $ "OF INITIAL" 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) "—CALC. MAG. OF FLOW FOR EACH AIRFOIL" VI = SQRT(q+YDFl**2) V2 = SQRT(Q+YDF2**2) "—CALC. RAMP REF. INPUT" , v v , , R3 = AMP*(RAMP(T1)-RAMP(T2))/(T2-T1) "—CALC. 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 "—CALC. ACTUAL FOIL ROTATIONS USING ACTUATOR MODEL1 BET1 = REALPL(TC.ERR1) BET2 = REALPL(TC.ERR2) "—CALC. EFFECTIVE ANGLE OF ATTACK OF EACH FOIL" ATTAC1 = ALPH+BET1-GAM1 ATTAC2 = ALPH+BET2-GAM2 "—CALC. 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) "—CALC. LIFT OF EACH FOIL" FL1 = 2.0*CL1*SPAN1*CORD1*V1**2 FL2 = 2.0*CL2*SPAN2*C0RD2*V2**2 "—CALC. X AND Y COMP, OF LIFT" FL1Y = FL1*C0S£GAM1) FL2Y = FL2*C0S(GAM2) FL1X = FL1*SIN(GAM1) FL2X = FL2*SIN(GAM2) "—CALC. 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 131 , 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 ... 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 - A T P T M n n r i r n 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)+ ... *LM)/20.0 .BLE TENSION" 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 tstop =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 A l l = A12 = A13 = A14 = A15 = A16 = A21 = A22 = = 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 S A23 = -0.4714 S A24 = -0.2325 S A25 = 5.7885 S A26 = 0.9004 S A31 = -0.1927 S A32 = -0.1959 S 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, Cal-ifornia, 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 Un-derwater 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. , Cam-bridge, 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 method-ology. In Proceedings of the 1986 American Control Conference, American Auto-matic 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. MIT Press, 1986. J .N. Nielson. Missile Aerodynamics. McGraw-Hill, 1960. B. Paul and A.I. Soler. Cable dynamics and optimum towing stategies for sub-mersibles. MTS Journal, 6(2), 1972. Kim 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 |
AggregatedSourceRepository | 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