Electromechanical retrofit of submarine control surface actuators

Copyright by Jonathan Robert LeSage 2010 Electromechanical Retrofit of Submarine Control Surface Actuators by Jonathan Robert LeSage, B.S.M.E. Thesis Presented to the Faculty of the Graduate School of The University of Texas at Austin in Partial Fulfillment of the Requirements for the Degree of Master of Science in Engineering The University of Texas at Austin May 2010 The Thesis Committee for Jonathan Robert LeSage Certifies that this is the approved version of the following thesis: Electromechanical Retrofit of Submarine Control Surface Actuators Approved by Supervising Committee: Raul G. Longoria Joseph J. Beaman To Stephanie for her continuous inspiration To my family for their ubiquitous support Acknowledgments The author would like to thank Raul Longoria for providing guidance and assistance during the course of this project. Additionally, the author would like to extend thanks to Bill Shutt for additional direction and the Applied Research Laboratories for the opportunity to work on this project. The author would also like to thank the Office of Naval Research and the Electric Ship Research and Design Consortium for providing the funding for this project. Jonathan Robert LeSage The University of Texas at Austin May 2010 v Electromechanical Retrofit of Submarine Control Surface Actuators Jonathan Robert LeSage, M.S.E. The University of Texas at Austin, 2010 Supervisor: Raul G. Longoria The current trend of the United States naval fleet has shifted towards the electrification of next generation vessels. In the case of large actuation systems such as the control surface actuators on submarines, the retrofitting of legacy hydraulic systems with modern electromechanical systems requires keen understanding of all actuation subsystems. This thesis explores two methods of characterization of the actuation and energy storage demands. Each technique relies on the utility of adopting a model basis, and specifically uses a bond graph approach in combination with two-port immittance relations. The conventional approach requires characterization and simulation of the entire submarine and control surface systems. This model basis, which incorporates the entire system, allows torque and power load estimations for known extreme maneuvers. Energy requirements for actuation in these scenarios provide insight vi into necessary energy storage. Two-port impedance synthesis techniques are also applied to the characterization of actuation subsystems in a more analytical process. These techniques explore how passive and active compensation may offer alternatives to established ways for representing actuation systems. A resulting simplified passive compensation function, which is physically realizable, is found for the control surface system for ocean disturbance rejection. Additionally, further contributions to the field of network synthesis as applied to physical systems include: further insight into realizable desired frequency response specifications, equivalent impedance for domain transformation, and model order reduction for simplified realizations. vii Contents Acknowledgments v Abstract vi List of Tables xi List of Figures xii Chapter 1 Introduction 1 1.1 Motivation for Actuation Retrofit . . . . . . . . . . . . . . . . . . . . 2 1.2 Overview of Synthesis Techniques . . . . . . . . . . . . . . . . . . . . 5 1.2.1 Impedance Based Approaches: Passive and Active Electrical Based Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 5 Extension of Impedance Based Methods: Multidisciplinary Synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.3 Active Element Synthesis . . . . . . . . . . . . . . . . . . . . 7 1.2.4 Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . 8 Chapter 2 Mathematical Model of Submarine Vehicular Dynamics 9 2.1 System Coordinates and State Variables . . . . . . . . . . . . . . . . 9 2.2 Vehicular Kinematics and Frame Transformation . . . . . . . . . . . 11 2.3 Nonlinear Submarine System Dynamics . . . . . . . . . . . . . . . . 12 2.4 Virtual Mass Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.5 Hydrodynamic Restoring Forces/Moments . . . . . . . . . . . . . . . 14 2.6 Potential Hydrodynamic Damping . . . . . . . . . . . . . . . . . . . 15 2.7 Control Inputs for Submarine Dynamics . . . . . . . . . . . . . . . . 17 viii 2.8 2.9 Full Vehicular Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.8.1 Parameter Calculation for Simulation . . . . . . . . . . . . . 18 Control Surface Moments/Forces . . . . . . . . . . . . . . . . . . . . 19 2.9.1 Moment Calculation from Nonlinear Theory of Airfoils . . . . 22 2.9.2 Hoerner Empirical Technique for Linearization . . . . . . . . 27 2.9.3 Summary of Submarine Dynamic Modeling . . . . . . . . . . 28 Chapter 3 Control Surface System Analyzation 29 3.1 Legacy Hydraulic System . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2 Control Surface Actuation System Modeling . . . . . . . . . . . . . . 33 3.3 Actuation Power Limitations and Power Curve . . . . . . . . . . . . 38 3.4 Application of System Model with Submarine Model . . . . . . . . . 39 3.4.1 Control Surface Load Estimation via Submarine Model . . . 40 3.4.2 Full Surface Sweep Maneuver Simulation . . . . . . . . . . . 40 3.4.3 Emergency Surfacing Maneuver Simulation . . . . . . . . . . 41 Energy Storage Requirements for Actuation . . . . . . . . . . . . . . 44 3.5 Chapter 4 Overview of Two-Port Synthesis 4.1 4.2 49 Essential Concepts for Synthesis Theory . . . . . . . . . . . . . . . . 49 4.1.1 Description of Two-port systems . . . . . . . . . . . . . . . . 50 4.1.2 Canonical Matrix Forms and the Corresponding Parameters . 51 4.1.3 Conversion between Canonical Forms . . . . . . . . . . . . . 53 4.1.4 Impedance Bond Graphs . . . . . . . . . . . . . . . . . . . . 55 4.1.5 Description of Synthesis . . . . . . . . . . . . . . . . . . . . . 56 Active Vehicle Suspension Example . . . . . . . . . . . . . . . . . . . 58 4.2.1 Bond Graph Model and Matrix Derivation . . . . . . . . . . 58 4.2.2 Specify Desired Frequency Response . . . . . . . . . . . . . . 62 4.2.3 Synthesized Impedance Function, Functional Expansion . . . 63 4.2.4 New Bond Graph Layout, Physical Realization . . . . . . . . 65 4.2.5 Simulation of Actuation System with Synthesized Compensator 69 4.2.6 Hydraulic Realizations via Secondary Synthesis . . . . . . . . 4.2.7 Simulation of Hydraulic Actuation System with Synthesized 4.2.8 71 Compensator . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 Summary of the Synthesis Methodology . . . . . . . . . . . . 79 ix Chapter 5 Synthesis of Control Surface Actuator Systems 80 5.1 Conventional and Synthesis Design Comparison . . . . . . . . . . . . 81 5.2 Bond Graph Model and Matrix Derivation with Unknown Impedance 81 5.3 Specify Desired Frequency Response . . . . . . . . . . . . . . . . . . 85 5.3.1 Natural System Response to Torque Disturbances . . . . . . 85 5.3.2 Ocean Disturbance Torque Spectrum . . . . . . . . . . . . . . 86 5.3.3 Notch Filter on Peak Ocean Torque Frequency . . . . . . . . 88 5.3.4 Band Reject on Ocean Torque Frequencies . . . . . . . . . . . 89 5.4 Synthesized Impedance Function . . . . . . . . . . . . . . . . . . . . 91 5.5 Synthesized Bond Graph and Physical Realization . . . . . . . . . . 94 5.5.1 Purely Active Realization of Synthesized System . . . . . . . 95 5.5.2 Separated Active and Passive Subsystem Realization . . . . . 101 5.6 Model Reduction of Synthesized Impedance . . . . . . . . . . . . . . 108 5.7 Synthesis Realization Comparison . . . . . . . . . . . . . . . . . . . 116 Chapter 6 Conclusions 120 6.1 Summary and Contributions . . . . . . . . . . . . . . . . . . . . . . . 120 6.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6.3 6.2.1 Extension of Synthesis Techniques . . . . . . . . . . . . . . . 122 6.2.2 Statistical Analysis of Power Variables . . . . . . . . . . . . . 122 Final Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 Appendix A Realizability Conditions for Passive Systems 125 A.1 One-port Passive System Conditions . . . . . . . . . . . . . . . . . . 125 A.2 Methods for determination of positive-realness . . . . . . . . . . . . 126 A.3 Extension to Two-port Systems . . . . . . . . . . . . . . . . . . . . . 127 Appendix B Balanced Realizations for Order Reduction 129 B.1 Conditions for a Balanced Realization . . . . . . . . . . . . . . . . . 129 B.2 Method to Determine Balanced Realizations . . . . . . . . . . . . . . 131 Bibliography 132 Vita 137 x List of Tables 1.1 Specified Performance Parameters . . . . . . . . . . . . . . . . . . . 4 2.1 Hydrodynamic and added mass coefficients . . . . . . . . . . . . . . 20 2.2 Control surface parameters to determine hydrodynamic forces . . . . 22 3.1 Control surface system parameter list . . . . . . . . . . . . . . . . . 34 xi List of Figures 1.1 Overall actuation system package with vital subdivisions. . . . . . . 3 2.1 Coordinate frame definitions of submarine vehicular dynamics. . . . 10 2.2 Depiction of the submarine system inputs and control surface hydrodynamic moment definition. . . . . . . . . . . . . . . . . . . . . . . . 21 2.3 Fin lift forces and moment arms . . . . . . . . . . . . . . . . . . . . 22 2.4 Center of pressure as a percentage of the length of the chord . . . . 23 2.5 Coefficient of lift for a NACA 0022 airfoil . . . . . . . . . . . . . . . 24 2.6 Coefficient of drag for a NACA 0022 airfoil . . . . . . . . . . . . . . 25 2.7 Hydrodynamic moments for a constant u = 15 m/s control surface for varying angles of attack. . . . . . . . . . . . . . . . . . . . . . . . 3.1 26 Legacy hydraulic system supply lines and primary accumulator energy storage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.2 Legacy hydraulic system and power transmission to control surfaces 31 3.3 Standard submarine spade actuation surface . . . . . . . . . . . . . . 32 3.4 Simplified layout of control surface system. . . . . . . . . . . . . . . 33 3.5 Bond graph form of the control surface system . . . . . . . . . . . . 34 3.6 Block diagram for control of control surface system . . . . . . . . . . 36 3.7 PID tuned frequency domain response comparison . . . . . . . . . . 37 3.8 Power limiting curve for actuation systems . . . . . . . . . . . . . . . 39 3.9 Control surface demands during full sweep simulation . . . . . . . . 41 3.10 Selected submarine states for full sweep maneuver. . . . . . . . . . . 42 3.11 Submarine states for rapid ascension maneuver . . . . . . . . . . . . 43 3.12 Control surface demands during emergency surfacing . . . . . . . . . 44 xii 3.13 Comparison of energy density to power density of energy storage devices 45 3.14 Potential centralized energy storage and distribution layout for submarine system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.15 Another potential energy storage and decentralized layout for submarine system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1 Examples of two-port system models . . . . . . . . . . . . . . . . . . 50 4.2 Finding z11 for a 2 port electrical network. . . . . . . . . . . . . . . . 51 4.3 Two-port system canonical forms . . . . . . . . . . . . . . . . . . . . 52 4.4 Relationships between canonical two-port systems . . . . . . . . . . 54 4.5 Basic elements of bond graph notation with causal definitions . . . . 55 4.6 Two-port transmission matrix forms for bond graph junctions . . . . 56 4.7 Analogous forms of system synthesis . . . . . . . . . . . . . . . . . . 57 4.8 Synthesis procedure for n-port systems to generate frequency domain compensators . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Active suspension system with unknown impedance function . . . . 60 4.10 Bond graph representation of the active suspension system . . . . . . 61 4.11 Experimentally obtained desired system frequency response (Vp /Vt ) . 62 4.9 4.12 Resulting frequency response for two-port matrix elements after system synthesis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.13 Synthesized bond graph elements separated into both active and passive components. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.14 Realization of the synthesized bond graph . . . . . . . . . . . . . . . 67 4.15 Flow source controller algorithm as determined by the synthesis expansion incorporated in the bond graph structure. . . . . . . . . . . 68 4.16 Velocity output at 10 rad/s comparison between suspension with only passive elements and a combination of active/passive. . . . . . . . . 69 4.17 Frequency response comparison between suspension with only passive elements and a combination of active/passive . . . . . . . . . . . . . 70 4.18 Force-speed of linear actuator at 35 rad/s input tire disturbance. . . 71 4.19 Force-speed of linear actuator at random tire disturbance. . . . . . . 72 4.20 Simple hydraulic realization of equivalent active vehicle suspension with a hydraulic piston coupled with a variable flow pump. . . . . . xiii 73 4.21 Bond graph of simple hydraulic realization with corresponding compensation control algorithm. . . . . . . . . . . . . . . . . . . . . . . . 73 4.22 Original simple hydraulic system with impedance function represented as an element and the new pure hydraulic form. . . . . . . . . . . . . 74 4.23 Pure hydraulic actuation realization in bond graph form with modified compensation controller algorithm. . . . . . . . . . . . . . . . . . 76 4.24 Pure hydraulic realization of the synthesized active vehicle suspension system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.25 Frequency response comparison between pure hydraulic realization and the desired frequency response. . . . . . . . . . . . . . . . . . . . 78 4.26 Pressure-flowrate curve for hydraulic actuator . . . . . . . . . . . . . 78 5.1 Flowchart comparison of the conventional design process in contrast with the synthesis procedure. . . . . . . . . . . . . . . . . . . . . . . 82 5.2 Control surface system layout with added unknown impedance element. 83 5.3 Bond graph representation of control surface model with specified unknown impedance function. . . . . . . . . . . . . . . . . . . . . . . 84 5.4 Unactuated frequency response of control surface system. . . . . . . 86 5.5 Magnitude of the s-Domain Representation of the Pierson-Moskowitz ocean spectrum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Proposed notch response alteration to impose on the natural control surface system response . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 89 Band reject response alteration to impose on the natural control surface system response . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 88 Comparison of the desired frequency response and the natural system response with a notch function. . . . . . . . . . . . . . . . . . . . . . 5.8 87 90 Comparison of the desired frequency response and the natural system response with a band reject function. . . . . . . . . . . . . . . . . . . 91 5.10 Comparison of the desired frequency response and the natural system response with a band reject function. . . . . . . . . . . . . . . . . . . 92 5.11 Frequency response of the synthesis derived actuation impedance function for the control surface system. . . . . . . . . . . . . . . . . . 93 5.12 Synthesized bond graph elements separated into both active and passive components. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv 94 5.13 Control surface system realization with ideal rotational actuator. . 95 5.14 Block diagram representation of control surface system and compensator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.15 Torque source controller algorithm as determined by the synthesis expansion incorporated in the bond graph structure. . . . . . . . . . 97 5.16 Control surface displacement and angular velocity with single actuation impedance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.17 Control surface actuator torque output and power flow velocity with single actuation impedance. . . . . . . . . . . . . . . . . . . . . . . . 99 5.18 Control surface actuator torque output and power flow velocity with single actuation impedance. . . . . . . . . . . . . . . . . . . . . . . . 100 5.19 Torque-Speed Characteristic Curve for single actuation impedance case compared to closed loop actuation. . . . . . . . . . . . . . . . . 100 5.20 Frequency response of actuation system with single actuation impedance actuator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.21 Control surface system realization with passive components separated. 102 5.22 Block diagram representation of control surface system and compensator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.23 Torque source controller algorithm with passive components separated.104 5.24 Control surface displacement and angular velocity with separated passive elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.25 Control surface actuator torque output and power flow with separated passive elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.26 Energy consumption by actuation system disturbance rejection with separated passive elements. . . . . . . . . . . . . . . . . . . . . . . . 106 5.27 Torque-Speed Characteristic Curve for separated passive elements case compared to closed loop actuation. . . . . . . . . . . . . . . . . 107 5.28 Frequency response of actuation system with separated passive elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.29 Comparison of original fourth order system to the reduced order second order approximation. . . . . . . . . . . . . . . . . . . . . . . . . 110 5.30 Reduced order impedance function in bond graph form yielding a purely passive system. . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.31 Reduced order mechanical realizations . . . . . . . . . . . . . . . . . 112 xv 5.32 Block diagram representation of control surface system and passive compensator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5.33 Unknown equivalent two port impedance function . . . . . . . . . . . 113 5.34 Passive compensation as realized in electrical form inside electromechanical actuator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.35 Electrical schematic of electrical compensation function. . . . . . . . 116 5.36 Control surface displacement and angular velocity with either reduced mechanical or electrical system. . . . . . . . . . . . . . . . . . . . . . 117 5.37 Control surface actuator torque output and power flow with either reduced mechanical or electrical system. . . . . . . . . . . . . . . . . 117 5.38 Energy consumption by actuation system disturbance rejection with either reduced mechanical or electrical system. . . . . . . . . . . . . 118 5.39 Reduced order torque-speed characteristic curve . . . . . . . . . . . . 118 5.40 Frequency response of actuation system with either reduced mechanical or electrical system. . . . . . . . . . . . . . . . . . . . . . . . . . 119 xvi Chapter 1 Introduction Current initiatives of the Office of Naval Research (ONR) have focused on the electrification of the naval fleet with a new generation of warships. With the increased interest in the next generation fleet, ONR established the Electric Ship Research and Development Consortium (ESRDC) to concentrate and centralize the various research on the fleet electrification project. The consortium includes five universities and the Naval Surface Warfare Center (NSWC). The University of Texas at Austin has focused on problems regarding thermal management, power systems and energy storage, and control system reconfiguration with fellow research institutions: the Applied Research Laboratories, the Center for Electromechanics, and the Microelectronics Research Center. This thesis focuses on the power and energy requirements for submarine control surface actuators during various emergency maneuvers. The model basis provides insight into potential demands on energy storage mechanisms in the event of power failure. In addition to a simulation basis for identification of potential load requirements, an impedance synthesis methodology is also utilized to procure impedance functions consisting of passive and/or active elements, which can be physically realized. A two-port synthesis methodology, more specifically, provides not only a model basis but also allows an analytical generation of potential system topologies due to the usage of the bond graph nomenclature. Comparisons of synthesized compensation functions are presented in conjunction with a tuned PID controller for disturbance rejection. 1 1.1 Motivation for Actuation Retrofit As a current priority, the Navy is looking into methods to reduce life cycle costs for future submarine deployment [1, 2]. One such aspect of cost reduction is the replacement of existing on board control surface hydraulic actuator designs, which have been in service since World War II, with electromechanical actuators. Control surface actuation for shipboard systems has relied heavily on hydraulic systems, taking advantage of legacy design knowledge. However, these hydraulic systems prove costly to construct and sustain. Initial studies suggest whole life costs (WLC) could be reduced by up to 50% [3] due to the conversion from hydraulic to electromechanical systems. A recent physical retrofitting of a hydraulic actuator on a Type 23 Frigate with an electromechanical system was performed [4]. Two linear actuators were used to control the aft fixed fin stabilizers on the vessel. To the knowledge of the author, this is the first instance of an electromechanical retrofit on a naval vessel in the research literature. The results of the study demonstrated not only the feasibility of the retrofitting process but also the viability of electromechanical systems in the naval setting. Further justification for continued deployment of these systems should show the ease of maintainability of electromechanical systems over hydraulics. Hydraulic systems entail additional piping and auxiliary power conversion systems, which add considerably to both construction and maintenance overheads. Hydraulic fluids are drained and the system is flushed in a process that takes upward of 48 hours in addition to any typical servicing which would be standard for any actuation system [4]. Moreover, modern electrical actuators employ intelligence maintenance monitoring with integrated sensors and the ability for prognostics on the fly. Additionally, retrofitting with electromechanical systems is motivated by electrification of shipboard systems, which promises improvements in actuation systems and resolution of problems associated with hydraulic systems. Not only do hydraulic systems add additional energy conversion paths in the submarine system, additional crew training beyond electrical knowledge is required for maintenance and emergency repairs. The study by Stafford [4] found that a electromechanical conversion required only two hours a month for planned maintenance requirements with other hourly reductions resulting from considerable time savings during dock- 2 Actuator System Power Conversion Prime Mover Control Power Electronics Power Source End Effect Motion Conversion Energy Storage Figure 1.1: Overall actuation system package with vital subdivisions. side maintenance [4]. In modern submarines, actuators serve as a synthesis of functions: motion conversion, power conversion, energy storage and system control as shown in Figure 1.1. Due to the hostile nature of the ocean, control surface actuators must function in strenuous circumstances and ensure functionality for surfacing in event of power failure. Control surface actuators must also be capable of meeting rigorous torques demands during diverse actuation maneuvers while maintaining high levels of reliability and impact tolerance. Other than situations of overheating and extended over loading conditions, electrical actuators for control surface replacement attained 98.5 percent availability [4]. Furthermore, electromechanical actuators, which are compact and lightweight, are ideal candidates to replace current hydraulic actuators. Due to the relative simplicity and compactness of electromechanical system, the retrofitted system should yield a reduction in total actuation system weight. The proposed electromechanical system tested by Stafford resulted in a 500 kg weight reduction over the existing hydraulic system [4]. The conversion of electrical voltage generated by the submarine’s power plant to hydraulic pressure results in energy loss. Since all energy conversion processes involve a loss in available energy, removal of the unnecessary conversion process results in energy conservation. One potential cause of concern for electromechanical actuation systems is the continuous power consumption present due to actuation position holding. In contrast, hydraulic systems provide continuous position holding with little power consumption. However, the experimental actuation fitting by Stafford 3 [4] resulted in the electromechanical actuator consuming 11.16 MJ in comparison to the standard hydraulic actuation system, which consumed 23.76 MJ producing a 53.03% reduction in energy consumption [4]. Since hydraulic actuation in submarines has an extensive history and many system designs and controls have been determined using hydraulics, a technique to determine actuator specifications and component topology is desired. By characterization of the existing hydraulic actuators in both time and frequency domains, electromechanical actuators with identical characteristics can utilize existing control algorithms with indistinguishable responses. Additionally, investigation of the actuation systems must consider submarine system-wide effects on power consumption and motion dynamics. Considerable concentration of mass towards the aft section of the ship results in stability issues for submarine motion. The Stafford report indicates strong evidence for the favorability of electromechanical system retrofitting. The report presents experimental confirmation of the possibility for a retrofit of existing hydraulic systems with positive benefits. However, many issues and concerns remain unanswered for a finalized retrofit solution such as system component topology, energy storage requirements/placement and power consumption during intense maneuvers. This thesis addresses each of these problems in turn and explores potential solutions for a final retrofit design in accordance to the specified parameters in Table 1.1. Table 1.1: Specified Performance Parameters [1] Parameter Specification Operating Torque Capacity 700, 000 ft-lbs (950 kN-m) 1, 000, 000 ft-lbs (1,360 kN-m) Peak Torque Capacity Operational Angular Displacement ±35 degrees Rate of Angular Displacement 5 ± 1 deg/s Position Accuracy ±0.25 degrees Energy Storage 1.4 MJ Since the retrofitting of the control surface actuators suggests system insertion, there is a need for methods to aid formulation of requirements for actuation subsystems that correspond to submarine operational requirements. Modern modeling and simulation tools enable a total systems approach that allows the possibility to propagate overall system response requirements into specified subsystems. How4 ever, early or concept design applications requires that physics-based models be developed with minimal design information. Detailed selection and sizing requires that requirements be formulated, and it can be helpful to have systematic ways to model and simulate systems for gaining insight during this process. The following section details the development of the impedance based synthesis techniques that are used in this thesis in conjunction with system modeling and parameterization techniques. 1.2 Overview of Synthesis Techniques Forms of the synthesis methodology have been present since the early twentieth century. Early circuit designers desired to produce circuit topologies given the desired performance requirements. Due to the ease of use, needing only a simple derivation and simple identification of elementary model elements, the transfer function became the ideal candidate for synthesizing results. Initially, only the ability to produce purely passive systems (inductors, capacitors and resistors) was known. However, the synthesis methodology has been expanded to include active systems (op amps, controlled voltage or current sources), as well as n-port systems. The synthesis methodology was eventually expanded, via the usage of the bond graph notation, to include not only electrical system but also mechanical systems (translational, rotational, hydraulic, etc.). Once extended to these new domains, many of the classical synthesis techniques were quickly adopted. 1.2.1 Impedance Based Approaches: Passive and Active Electrical Based Systems With the mathematical foundations of transfer functions, matrix algebra, and impedance modeling, the synthesis procedure resulted as an offshoot of a combination of mathematical techniques in both modeling and early control theory. Original work using the synthesis methodology was developed for the realization of electrical networks with specified performance specifications. Initial developments by Brune [5], Cauer [6], and Foster [7] laid the foundation to the classical synthesis procedures. Their work established many of the theorems for physical realization of derived impedance functions as well as establishing the necessary conditions for identification of passive 5 impedance elements. Additionally, their work established the procedure of decomposing resulting impedance function with functional expansion techniques such as partial fraction expansion. However, these initial techniques were limited in scope and only applicable to single input single output systems which contained purely passive elements. The synthesis methodology was extended in the 1950’s to active electrical systems with the realization of the forms of active impedance functions [8, 9, 10]. Active impedances were identified as elementary model components with negative impedances (such as negative value resistors, capacitors and inductors). These active impedance functions were found to be physically realizable via op-amps and/or controllable voltage/current sources. The extension to synthesize both active and passive functions is recapitulated in the works of Newcomb [8]. Additionally, at the same time, the synthesis methodology was extended from SISO systems to include n-port systems. The extension to include n-port systems, conducted by Bayard [10], Balabanian [11], Hazony [12], and Newcomb [8], laid out the framework for the field now known as classical synthesis. The potential for n-port models extended the capabilities of modeling to include any potential electrical system with both active and passive systems. Another result of the formation of classical synthesis techniques was the identification of canonical matrix model forms [13]. The canonical matrix forms commonly appear in physical systems and transformation between the canonical forms requires only simple algebraic manipulation. Classical synthesis of electrical networks provided the outline from which future researchers could extend to other physical domains. 1.2.2 Extension of Impedance Based Methods: Multidisciplinary Synthesis Since the existing methodology for electrical system synthesis relies primarily on impedance functions, an extension to multidisciplinary system is possible through bond graph formalism. Bond graph modeling, developed by Paynter [14] and detailed in [15], is a technique for visualization of energy flow throughout a system with the ability to indicate the causality of the energy flow. The fundamental elements of the bond graph represent generalized resistive, inductive and capacitive 6 elements and thus provide an ideal expansion to the classical synthesis techniques. Expanding on the early ideas of Brune and Cauer, bond graph formalism was used to synthesize a SISO mechanical system by Redfield and Krishnan [16, 17]. As an analogous mechanical system, they were able to produce a passive mechanical realization that provides a specified frequency response. This extension provided a modified methodology, which included bond graph representation of the system as well as the resulting synthesized impedance function. Due to the formalism inherent in bond graph notation, synthesized impedance functions, once functionally expanded, indicated internal connections between the components. 1.2.3 Active Element Synthesis With a basis for multidisciplinary passive synthesis, Connolly and Longoria further extended bond graph synthesis techniques to include active multidisciplinary elements [18, 19]. Again, much akin to the classical electrical network synthesis procedure, the active fundamental elements were modeled as negative impedance functions. Realization of these negative impedance functions produce a compensation algorithm from which a signal from the system is drawn to produce a certain input. In the active vehicle suspension synthesis in Connolly, the resultant active impedance function produced a damping torque in accordance to a measured velocity signal. Kim [20] made the extension to n-port networks with bond graph formalism. However, Kim [20] found that synthesized bond graph networks were difficult to realize for complex multiport networks. Multiport bond graphs could represent many of these complex systems, but again physical realization of the interconnected multiport bond graph elements was tenuous at best. However, with certain restrictions in play, the complexity of the synthesized results could be reduced. With the extension to bond graph methodology, the synthesis procedure is capable of producing a multidisciplinary impedance function which through bond graph formalism can be expanded and structuralized. The resulting synthesized topologies represent one of many potential results which satisfy the specified parameters. 7 1.2.4 Thesis Overview The overall goal of this thesis is to highlight the aspects of both conventional and synthesis techniques as applied to retrofitting systems. The subsequent chapters build the reader’s knowledge base of submarine dynamics, actuation systems and synthesis techniques for comprehension of the comparison of design techniques. Chapter 2 presents an overview of the derivation of the nonlinear submarine vehicular dynamics equations. Additionally, this section details the dynamics of the control surface stern planes for both the non-linear submarine model and the linearized actuation system model. After building a foundation on submarine dynamics, a conventional methodology for actuation sizing, based on model basis simulation, is presented in Chapter 3. In addition to an overview of the approach, several realistic submarine maneuver simulations are presented. Chapter 4 introduces the two-port methodology. Two-port systems are introduced and the synthesis of both active and passive elements with bond graph nomenclature is presented. Chapter 4 culminates with a full example of the synthesis methodology as applied to an active vehicle suspension system. Chapter 5 integrates the modeling of the control surface actuation system with the synthesis methodology. Synthesis is conducted for removal of flow-induced ocean disturbance torques. Several non-unique systems are analytically derived as a result. 8 Chapter 2 Mathematical Model of Submarine Vehicular Dynamics In this chapter, the time-invariant nonlinear submarine vehicular kinematics and dynamics are introduced. For the sake of the reader, the equations are partially derived to reveal the framework of the submarine motion and forces. Additionally, a complete presentation of all of the equations with definitions of every term can be found in Gertler [21] and Feldman [22] and a comprehensive derivation of all of the equations can be found in Watt [23, 24]. The notation of derivation follows that of Fossen [25] with some modifications. 2.1 System Coordinates and State Variables By establishing simple coordinate frame relationships, the overall complexity of a six degree of freedom (DOF) model can be reduced to the least complex non-linear form. In the case with a single moving body, two coordinate systems are used to depict the submarine kinematics, as shown in Figure 2.1. The earth-fixed frame, defined as the inertial frame, is fixed at some point in space, and the body-fixed frame is fixed to the submarine. The body-fixed coordinate system is selected to correspond with the principal axes of inertia to aid in the simplification of derivation. The state variables of the system are all measured in relation to these defined reference frames. With respect to the defined axis and the right hand rule, the Euler angles are (φ, θ, ψ) and the vehicular angular momen9 Earth-Fixed Reference Frame φ X θ Y ψ δr - rudder Body-Fixed Reference Frame Z δb - bowplanes δs - sternplane p - roll q - pitch x y v - sway u - surge r - yaw z w - heave Figure 2.1: Coordinate frame definitions of submarine vehicular dynamics. The rotations about the earth-fixed axis are defined as the Euler angles (φ, θ, ψ) and the rotations about the body-fixed axis are the vehicle angular momentums (p, q, r). 10 tums are (p, q, r). The position and orientation of the submarine are relative to the earth-fixed frame and the linear and angular velocities are relative to the body-fixed frame. Representing these expressions in terms of vectors, the system state variables are η and ν and the forces and moments are τ as illustrated by equation (2.1). η= ν= τ= 2.2 h h h x y z ϕ θ ψ u v w p q r X Y Z K M iT iT N (2.1) iT Vehicular Kinematics and Frame Transformation With the system state variable vectors defined, the transformation between the two frames is found by relating the angular orientations. In this derivation, the Euler angles are used due to the intrinsic relation to the notation of pitch, roll, and yaw. By applying the frame transformation matrix, the submarine states can be indicated in the global sense of the earth-fixed frame. The relation between state vectors is important for determining system response and eventual characterization from the earth-fixed frame. For a full derivation of the successive transformations that make up the resulting frame transformation matrix consult Watt [24]. The velocity vector of the earth-fixed frame and the velocities of the submarine fixed reference frame are related by equation (2.2). η̇ = J(η)ν (2.2) where J(η) is the rotation matrix expressed as equation (2.3). J(η) = " J1 (η2 ) 03×3 03×3 J2 (η2 )] # (2.3) The subordinate matrices J1 (η2 ) and J2 (η2 ) are given by equation (2.4) and equation (2.5) respectivly. 11  cos ψ cos θ J12 (η2 )  J1 (η2 ) =  sin ψ cos θ J22 (η2 ) J13 (η2 ) J23 (η2 ) − sin θ cos θ sin ϕ cos θ cos ϕ   1 sin ϕ tan θ cos ϕ tan θ   J2 (η2 ) =  0 cos ϕ − sin ϕ     (2.4) (2.5) 0 sin ϕ/cos θ cos ϕ/cos θ where the remaining components of matrix (2.4) are given by equation (2.6). J12 (η2 ) = − sin ψ cos ϕ + cos ψ sin θ sin ϕ J13 (η2 ) = sin ψ sin ϕ + cos ψ cos ϕ sin θ J22 (η2 ) = cos ψ cos ϕ + sin ϕ sin θ sin ψ (2.6) J23 (η2 ) = − cos ψ sin ϕ + sin θ sin ψ cos ϕ 2.3 Nonlinear Submarine System Dynamics To derive the dynamics of the system, Lagrange’s equation of motion is used. This allows the establishment of the general forces in the system, which opposes the fluid surrounding the vehicle. In the formulation of Lagrange’s equations, the independent generalized coordinates for each degree of freedom are required. In this case, the coordinates chosen are the state variable η(x y z φ θ ψ) which make the other state variable ν(u v w p q r) dependent variables through equation (2.2). Defining qk as each of the independent generalized coordinates, the equations of motion are found by equation (2.7).   ∂L d ∂L − = Qk dt ∂ q̇k ∂qk (2.7) where t is time, Qk are the fluid forces related to qk and L is the system Lagrangian. By expanding the terms of the Lagrange equation out, the equations in terms of the vectors η and ν can be found. For a full expansion of the Lagrange equations consult Watt [24]. Combining terms from the expansion of the Lagrange equations, the rigidbody dynamics of the submarine are denoted by equation (2.8). Mb ν̇ + Cb (ν)ν = τ 12 (2.8) where the vector ν represents the equation is defined with respect to the body. The rigid-body inertia matrix, Mb , is defined in equation (2.9). Cb represents the centripetal and Coriolis portions of the inertial matrix displayed in equation (2.10).  m 0 0 0 mzG   0 m 0 −mzG 0   0 0 m myG −mxG  Mb =   0 −mzG myG Ix −Ixy   mz 0 −mxG −Iyx Iy G  −myG mxG 0 −Izx −Izy −myG   mxG    0   −Ixz   −Iyz   Iz (2.9) where m is the mass of the submarine, I is the moment of inertia (or correspondingly, the product of inertia) with respect to the rotation, xg is the x-coordinate of the center of gravity, yg is the y-coordinate of the center of gravity, and zg is the zcoordinate of the center of gravity. Cb (ν) =  " 03×3 Cb1 (ν) T (ν) C (ν) −Cb1 b2 m(yG q + zG r) −m(xG q − w)  Cb1 (ν) =  −m(yG p + w) m(zG r + xG p) # −m(xG r + v) (2.10)   −m(yG r − u)  −m(zG p − v) −m(zG q + u) m(xG p + yG q)   0 −Iyz q + Iz r Iyz r − Iy q   Cb2 (ν) =  Iyz q − Iz r 0 −Ixz r + Ix p  −Iyz r + Iy q Ixz r − Ix p 0 (2.11) (2.12) Due to the components of vector ν in the Coriolis/centripetal term which couple with other terms in the ν vector, the dynamics of the submarine are nonlinear. Also note that the Cb matrix is in skew-symmetric form. 2.4 Virtual Mass Force The virtual mass (also known as the added mass) is the change in inertia of a system due to a body accelerating or decelerating through a fluid medium. This change in inertia arises due to the fact that the body and fluid cannot occupy the same space 13 simultaneously and fluid must be displaced. The virtual mass force assumption says that fluid accelerates with the body and thus adds inertia to the body. Since the virtual mass will contribute in a similar manner to the rigid-body dynamics, the equation for the force and moment contributions takes the form of equation (2.13). Mv ν̇ + Cv (ν)ν = τv (2.13) Much like the rigid-body dynamics which were derived above, the virtual mass contains both an inertial component and a Coriolis/centripetal component. From the derivation in Fossen [25], the virtual mass inertial matrix is given by equation (2.14).       Mv = −       Xu̇ Xv̇ Xẇ Xṗ Xq̇ Xṙ Yu̇ Yv̇ Yẇ Yṗ Yq̇ Zu̇ Zv̇ Zẇ Zṗ Zq̇ Ku̇ Kv̇ Kẇ Kṗ Kq̇  Yṙ   Zṙ    Kṙ   Mṙ   Nṙ Mu̇ Mv̇ Mẇ Mṗ Mq̇ Nu̇ Nv̇ Nẇ Nṗ Nq̇ (2.14) where the terms are the coefficients of virtual mass with respect to the corresponding components of the ν vector. Likewise in Fossen [25], the virtual Coriolis/centripetal matrix is equation (2.15).  0 0 0 0   0 0 0 Zẇ w   0 0 0 −Yv̇ v  Cv =   0 −Zẇ w Yv̇ v 0   Z w 0 −Xu̇ u Nṙ r  ẇ −Yv̇ v Xu̇ u 0 −Mq̇ q −Zẇ w 0 Xu̇ u −Nṙ r 0 Kṗ p Yv̇ v   −Xu̇ u    0   Mq̇ q   −Kṗ p   0 (2.15) where the terms are the coefficients of the Coriolis/centripetal contribution to the virtual force with respect to the corresponding components of the ν vector. 2.5 Hydrodynamic Restoring Forces/Moments Other forces on the submarine body are the effects of the weight and buoyancy of the vehicle. Since gravity always acts globally downward due to the earth-fixed 14 nature of gravity, the force due to gravity must be transformed to the center of mass in the vehicle frame. Similarly, the buoyancy must be transformed to the center of buoyancy in the submarine. The remaining force components of restoring moments can be expressed as a sum of the remaining forces/moments given by equation (2.16). τr = −g(η) (2.16) where g(η) are the restoring moments which are the gravity and buoyancy components. Note that g(η) is in terms of η. Since the buoyancy and the weight forces act due to the orientation of the submarine to the earth-fixed frame, a transformation between the earth-fixed frame and the body frame must be preformed. Using the transformation J1 (η2 ), the restoring moments/forces coefficients are given by equation (2.17).  (W − B) sin θ   −(W − B) cos θ sin ϕ   −(W − B) cos θ cos ϕ  g(η) =   −(yG W − yB B) cos θ cos ϕ + (zG W − zB B) cos θ sin ϕ   (zG W − zB B) sin θ + (xG W − xB B) cos θ cos ϕ  −(xG W − xB B) cos θ sin ϕ + (yG W − yB B) sin θ            (2.17) where (xg , yg , zg ) and (xb , yb , zb ) are the center of gravity and center of buoyancy respectively. 2.6 Potential Hydrodynamic Damping The remaining forces can be summarized as forms of hydrodynamic damping due to ocean wave disturbances, skin friction of the vehicle and damping due to vortex shedding due to form drag [25]. For this derivation, the submarine will be assumed to be significantly below the surface so that the interference from surface waves is minimal. The total potential damping force on the submarine body can be expressed as equation (2.18). τd = D(ν)ν 15 (2.18) where D is the matrix defining the potential damping from surface disturbances. The component of drag due to vortex shedding/form drag is a function of the frontal area shape of the submarine body. By integrating the pressure that is normal to the surface of the submarine, the force drag term is found and shown in equation (2.19). 1 f (U ) = − AρCd (Rn )U |U | 2 (2.19) where ρ is the density of the fluid medium, Cd is the coefficient of drag for the body, A is the frontal area, and U is the relative velocity of the viscous fluid over the body. Each degree of freedom of the body contains a second order force due to vortex shedding term. The coefficients of these terms are replaced with a single coefficient for simplicity. The potential damping coefficients due to the skin friction is found as in Watt [24] as the integral of the viscous stress tangent to the submarine body. As in the derivation of Watt [24], the terms higher than second order are negligible and can thus be dropped. Additionally, due to the planes of symmetry of the submarine body the damping matrices may be simplified to, D(ν) =   D11 (ν) =  " D11 (ν) D12 (ν) 03×3 D22 (ν) # (2.20) Xu + Xu|u| |u| 0 0 0 Yv + Yv|v| |v| 0 0 0 Zw + Zw|w| |w|   0 Xq|q| |q| 0   D12 (ν) =  Yp|p| |p| 0 Yr|r| |r|  0 Yq|q| |q| 0  Kp + Kp|p| |p| 0 Kr|r| |r|  D22 (ν) =  0 Mq + Mq|q| |q| 0 Yp|p| |p| 0  (2.21)   Nr + Nr|r| |r| (2.22)    (2.23) where the first terms in matrix are the coefficients of the skin friction factors and the second terms are the coefficients of the vortex shedding drag factors. 16 2.7 Control Inputs for Submarine Dynamics In addition to the passive dynamics of the submarine system, the model must incorporate system inputs, which allows active control of the system [25]. Due to the plethora of existing submarine designs, the control inputs to the system in this case are restricted to stern plane/bow plane control surfaces, aft rudders for yaw control and a propeller for surge velocity input. The control surfaces provide restoring forces into the system dynamics much akin to the hydrodynamic restoring forces, while the propeller force provides a simplified direct force in the surge plane. The force components of the control inputs can be expressed as a product of the input matrix and the input vector, τi = g(ν)ui (2.24) where g(ν) is the control parameter matrix which contains the parametric relationships of the input vector ui . The control parameter matrix g(ν) is given below,       g(ν) =       0 0 0 1 Yuuδr u2 0 0 0 Zuuδs u2 Zuuδb u2 0 0 0 0 Muuδs u2 Muuδb u2  0   0    0   0   0 Nuuδr u2 0 0 (2.25) where the control parameters are given by the size and shape of the control surfaces. The input vector, ui , is given as,  δr   δs ui =   δ  b Xprop       (2.26) where δr is the deflection of the rudder in radians, δs is the deflection of the stern plane in radians, δb is the deflection of the bow plane in radians and Xprop is the simplified thrust model which only applies force in the surge direction. 17 2.8 Full Vehicular Dynamics By compiling the forces and moments derived in the preceding analysis, the total forces and moments on the submarine body can be found. The total inertial matrix and Coriolis/centripetal matrix can be found by summing both the body dynamics and the added mass matrices, M = Mv + Mb (2.27) C = Cv + Cb (2.28) Combining these two summations with the total damping matrix D(ν), the total force/moment on the submarine body with respect to the body frame is found, M ν̇ + C(ν)ν + D(ν)ν + g(η) + g(ν)ui = τ (2.29) Nevertheless, the analysis of the motion of the submarine with respect to the earthfixed frame is of ultimate interest, so use the translation derived in equation (2.3). This transformation provides the vehicular dynamics in the earth-fixed frame, Mη (η)η̈ + Cη (η, ν)η̇ + Dη (η, ν)η̇ + gη (η) + g(ν)ui = τη (2.30) Mη (η) = J −T (η)Mtot J −1 (η) h i ˙ Cη (ν, η) = J −T (η) Ctot ν − Mtot J −1 (η)J(η) J −1 (η) (2.31) Dη (ν, η) = J −T (η)D(ν)J −1 (η) (2.33) gη (η) = J −T (η)g(η) (2.34) τη (η) = J −T (η)τ (2.35) where, 2.8.1 (2.32) Parameter Calculation for Simulation In the derivation of the submarine dynamic equations, a slew of hydrodynamic coefficients resulted, which model the motion of the rigid body in a particular direction. However, these coefficients exhibit considerable variation between various ocean faring crafts. For the estimation of the aforementioned coefficients, strip theory can be used to determine a slender body approximation of the hydrodynamic derivates. 18 The basis of strip theory is the separation of the vehicle body shape into many small strips [25]. The hydrodynamic coefficients are then calculated for each individual strip and the resulting coefficients are combined to form the total final coefficient. An example hydrodynamic derivative for a submerged slender vehicle is given by equation (2.36). −Yv̇ = Z L/2 A22 (y, z)dx (2.36) −L/2 where L is the length of the vehicle, and A22 is the strip area for a given y and z point on the vehicle cross section. However, due to the relative complexity of the submarine body, calculation of the hydrodynamic moments is ideal preformed by algorithms with simplified geometries. In Watt [23], a numerical strip technique is utilized to calculate scaling coefficients for each hydrodynamic derivative. These scaling factors are then combined with the density of the fluid medium and the length of the vessel to provide the total coefficient values. With the equations for the coefficients, each component is tabulated and indicated in Table 2.1. Since the focus of this thesis tends toward modern vessel retrofits, the submarine class of choice for this study was the Virginia class vessel [1]. A list of specifications and dimensions of the Virginia class submarine can be found in Table 2.2. For the values in Table 2.1, a length of 110 meters was used as indicated by Table 2.2. 2.9 Control Surface Moments/Forces In addition to the hydrodynamic forces that act on the submarine body, the moments and forces about the control surface rudders also play an important role in vehicular dynamics. Due to the structure of the control surface and the conditional assumption that the surface is completely immersed in a fluid media, the control surface experiences hydrodynamic lift identically as an airfoil. Likewise, the control surface will be subjected to a drag force and moments about the actuation shaft. Two techniques for the identifying the hydrodynamic forces on the control surface are outlined in this section. The first technique follows the theory of airfoils and utilizes the nonlinear coefficients of lift and drag to calculate the torque about the control surface [25, 26]. This fluid dynamics technique formulates an accurate method for identifying and modeling hydrodynamics. Alternatively, the Hoerner 19 Table 2.1: Hydrodynamic and added mass coefficients for submarine dynamics [23] Parameter Value Parameter Value Parameter Value xg 0.00 yg 0.00 zg 0.55 yb 0.00 zb -0.55 xb 0.00 Ixx 2.44 × 108 Iyy 8.13 × 109 Izz 8.13 × 109 Xu̇ Xẇ Xq̇ Kv̇ Kṗ Kṙ −3.69 × 105 −1.91 × 103 6.98 × 104 −2.16 × 107 −2.18 × 108 −1.67 × 108 Yv̇ Yṗ Yṙ Mu̇ Mẇ Mq̇ −1.32 × 106 −2.16 × 107 1.52 × 106 6.98 × 105 −1.90 × 107 −7.49 × 109 Zu̇ Zẇ Zq̇ Nv̇ Nṗ Nṙ −1.91 × 103 −1.09 × 107 −1.90 × 107 1.53 × 108 −1.67 × 108 −7.64 × 109 Xuu Xpp Xuq Yup Yr|r| Zq|q| Kur Mrr Mq|q| Nwp 1.00 × 104 −2.08 × 106 1.28 × 104 3.36 × 106 2.29 × 108 −2.29 × 108 −1.02 × 107 2.52 × 105 −1.06 × 101 0 −1.83 × 106 Yvv Xrr Xqq Yp|p| Zwp Kup Kr|r| Mwp Nup Nur 2.68 × 106 2.42 × 108 2.45 × 108 −1.92 × 105 −8.60 × 103 −3.13 × 107 −1.63 × 106 −2.44 × 105 −4.73 × 107 −4.77 × 108 Zww Xwp Xq|q| Yur Zuq Kp|p| Mpp Muq Np|p| Nr|r| 2.68 × 104 −1.65 × 105 3.70 × 104 7.43 × 106 −8.84 × 106 −1.63 × 106 2.11 × 107 4.49 × 108 −2.44 × 106 −1.06 × 101 0 Xδsδs Yδr −1.27 × 105 1.46 × 105 Xδrδr Nδr −1.27 × 105 −7.34 × 106 Zδs Mδs −1.58 × 105 −7.34 × 106 Xvr Ywr Ypq Xvp Kwp Kwr Mvr Nqr −Yv̇ Xẇ −Zq̇ −Yv̇ −Yṗ −Yṙ − Zq̇ Yṗ −Kṙ Xrp Yqr Zwq Zpp Kpq Kqr Mvp Nvq −Yṗ Xq̇ −Xẇ Yṗ Kṙ −Mq̇ + Nṙ −Yṙ −Xq̇ − Yṗ Xwq Ywp Zqq Zrp Kvq Mwq Mpr Npq Zẇ Zẇ Xẇ Yṙ Yṙ + Zq̇ Xq̇ Kṗ − Nṙ −Kṗ + Mq̇ empirical technique characterizes data from experimental studies on thin rigid bodies at small angles of attack [27]. This method of deriving hydrodynamic forces produces a form that is easily linearized for use in the impedance based synthesis techniques. Both techniques for identification of hydrodynamic forces requires a background of 20 δs Msp u - surge q - pitch w - heave Figure 2.2: Depiction of the submarine system inputs and control surface hydrodynamic moment definition. basic airfoil mechanics that are discussed in the following section [26]. Due to the surge momentum deviation by the cross-sectional area of the hydroplane and the resulting forces due to Newton’s second law, the control surface will generate a lift force. 1 FL = ρU 2 As CL 2 (2.37) where ρ is the fluid density, u is the surge velocity, A is the cross-sectional area and CL is the coefficient of lift. The coefficient of lift is dependent not only on the shape of the control surface but also on the angle of attack relative to the fluid medium. Likewise, the drag due to the fluid about the surface can be found, 1 Fd = ρU 2 As Cd 2 (2.38) where Cd is the coefficient of drag. The drag imposed by the control surfaces was also incorporated into the discussion of the derivation of the vehicular dynamics above. A full description of the drag induced by control surfaces can be found in Triantafyllou [28]. The drag and lift forces on the control surfaces result from the components of 21 Normal Force - N Lift force - L Control Rod δs Drag force - D x Fluid Velocity - Uo Cpx Figure 2.3: Fin lift forces and moment arms which produce hydrodynamic moments about the control surface. the three dimensional velocity vector. Each of the body-fixed velocities are coupled or have non-linear relationships, which ultimately result in the moment about the control surface. Formulation of the actuation torque requirements is based on the hydrodynamic torques induced by submarine motion and sub sea fluid flows. The parameters for the control surface were found for the Virginia class submarine vessel. Due to the lack of clear information on all parameters, reasonable estimates, established in Liu [29] for bearing damping, are used for simulation. The control surface parameters are summarized in Table 2.2. Table 2.2: Control surface parameters to determine hydrodynamic forces Parameter Value Description As 13.717 m2 Control surface planeform area c 2.00 m Chord length of control surface 0.3375 m Distance from leading edge to control rod x1 l 110 m Length of submarine 3 ρ 1010 kg/m Salt water fluid density 2.9.1 Moment Calculation from Nonlinear Theory of Airfoils The fluid dynamics technique approaches the fin hydrodynamics in a systematic manner similar to the airfoil dynamics techniques. In a similar manner, rather than attempting to analytically derive the gradients of fluid velocity and the fluid stress, empirical results of lift and drag coefficients are used. For this study, the 22 Center of Pressure - NACA 0022 Airfoil in Free Stream 45 Center of Pressure - CPc (%) 40 35 30 25 20 15 10 5 -40 -30 -20 -10 0 10 20 30 40 Angle of Attack (deg) Figure 2.4: Center of pressure as a percentage of the length of the chord from the leading edge of the control surface [26]. control surface profile was assumed to closely resemble that of a NACA 0022 airfoil [26]. Additionally, the control surfaces are modeled as have a spade rudder design, meaning that the entire surface deflects rather that a flap portion. One of the important components of the hydrodynamic forces is the center of pressure on the control surface. The location of the center of pressure determines the location of the hydrodynamic force lever arm. As the center of pressure shifts due to the change in angle of attack of the surface, the length of the force moment arm lengthens which in turn raises the torque imposed on the control surface actuator. The center of pressure as a percentage away from the leading edge of the control surface is illustrated in Figure 2.4. The center of pressure increases with respect to an increase in the angle of attack both positive and negative, as shown in Figure 2.4. Although the center of pressure in the x-direction does not change due to differing aspect ratios, the center of pressure in the y-direction depends on both angle of attack and aspect ratio. However, only the x-direction of the center of pressure is needed for hydrodynamic torques about the control rod. By looking up the center of pressure for a given angle 23 Coefficient of Lift - NACA 0022 Airfoil in Free Stream 1 0.8 Coefficient of Lift (CL) 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -40 -30 -20 -10 0 10 20 30 40 Angle of Attack (deg) Figure 2.5: Coefficient of lift for a NACA 0022 airfoil for differing angles of attack. [30] of attack, the lever arm from the center of pressure to the control rod is found by equation (2.39).   CPC (α) x̄ = c̄ − x1 100 (2.39) where c̄ is the chord of the surface and x1 is the distance from the leading edge to the control surface rod. The other angle of attack dependent coefficients determining the hydrodynamic forces are the coefficient of lift and drag. In a similar manner to the center of pressure, both coefficients rely on the airfoil profile and the angle of attack. Figures 2.5 and 2.6 display the empirically found coefficients of lift and drag for a NACA 0022 airfoil profile respectively [30]. The coefficient of lift typically increases with increased angle of attack until a peak lift angle is reached, as shown in Figure 2.5. Due to both the decrease in lift and increase in drag at angles greater than 30 degrees, the operation of submarine control surfaces is nominally restricted to plus or minus 30 degrees deflection. As seen in Figure 2.3, the summation of force vectors due to lift and drag produces a resulting normal force vector which is perpendicular to the control surface. This normal force 24 Coefficient of Drag - NACA 0022 Airfoil in Free Stream 0.7 Coefficient of Drag (CD) 0.6 0.5 0.4 0.3 0.2 0.1 0 -40 -30 -20 -10 0 10 20 30 40 Angle of Attack (deg) Figure 2.6: Coefficient of drag for a NACA 0022 airfoil for differing angles of attack. [30] vector, when crossed with the moment arm, results in the hydrodynamic moment. Rather than combining the force vectors, the coefficients can also be viewed as vectors, which can be combined. Combining the lift and drag coefficients yields equation (2.40). CN (α) = CL (α) cos α + CD (α) sin α (2.40) With the angle-dependent parameters defined, the force and moments on the surfaces can be investigated. Combining the equations of lift and drag and also accounting for the duality of the control surfaces, equation (2.41) results, FN,m (α) = ρU 2 As CN (α) (2.41) where As is the surface area of the control surface , ρ is the density of salt water and U is the coupled body-fixed velocity function which is dependent on both coupled velocities as well as control surface angles of deflection. Expanding the coupled body-fixed velocity function and accounting for the center of pressure moment arm, 25 Hydrodynamic torque - u = 15 m/s 2000 1500 Moment (kN-m) 1000 500 0 -500 -1000 -1500 -2000 -40 -30 -20 -10 0 10 20 30 40 Angle of Attack (deg) Figure 2.7: Hydrodynamic moments for a constant u = 15 m/s control surface for varying angles of attack. the hydrodynamic moment equation is found,   Mf (α) = ρAs CN (α)x̄(α) u2 α + uw + lp uq (2.42) Applying equation (2.42) to the NACA 0022 airfoil with parameters given by Table 2.2 with zero heave and pitch motion and a constant surge velocity of 15 m/s produces Figure 2.7. The hydrodynamic torques remain minute for small angles of attack and only contribute vigorously beyond a ±20 degree angle of attack. Note that the maximum specified torque output required of the electrical actuator is 1360 kN-m, so that angles of attack beyond 35 degrees are unobtainable at a surge velocity of 15 m/s. However since angles of attack of the surfaces are restricted to the ±30 degree range, all actuations lie within obtainable regions, outside of the inaccessible zone. 26 2.9.2 Hoerner Empirical Technique for Linearization Another technique, used to divulge the hydrodynamic moments imposed on the control surface by the surrounding fluid media, produces a simplified linear relationship between the angle of the control surface and hydrodynamic torque. The Hoerner realization of hydrodynamic forces results from empirical studies on thin rigid bodies at small angles of deflection [27]. This linearized form alleviates nonlinear system complications, which could result as a byproduct of the impedance synthesis techniques. However while this technique is ideal for the synthesis technique with a small angle assumption, large deflections are not properly reflected in the results. Moments are produced due to the offset of the actuation rod to the center of pressure of the lift and drag forces. By assuming a constant center of pressure, the moment arm for the hydrodynamic forces also remains constant. Finally, the coefficient of lift for small angles of attack can be represented by the Hoerner empirical fin lift coefficient (ch ) shown in equation (2.43) [28]. ch = As 1 + 1.8π πh2s (2.43) where As is the surface area of the control surface and hs is the height of the control surface. Since the deflection of the control surface is restricted to ±30 degrees, the Hoerner angle assumption should hold in this case. Altering the equation of lift and accounting for the duality of the control surfaces, equation (2.44) results, FL,m = ρU 2 As Cs (2.44) where U is the coupled body-fixed velocity function which is dependent on both coupled velocities as well as control surface angles of deflection. Expanding the coupled body-fixed velocity function and accounting for the center of pressure moment arm, the hydrodynamic moment equation is found,   Mf = ρAs Cs lp u2 δs + uw + lp uq (2.45) where lp is the moment arm to the center of pressure, δs is the control surface deflection, u is the surge velocity, w is the heave velocity and q is the pitch velocity. Coupled with the control surface dynamics, the hydrodynamic moments about the 27 control surface produce the required actuation torque for both steady state and transient actions. Since techniques presented in this thesis rely on linear impedance methods, the linearization of the hydrodynamic function proves necessary. The function contains both coupled and quadratic components, which should be removed for linearization. In the cases of interest, only motion in the surge direction is necessary while the other velocities are assumed to be zero. Additionally, the surge velocity is assumed to be constant so that both the quadratic and coupled components are eliminated. The linearized hydrodynamic moment is given by equation (2.46). Mf,l = ρAs Cs lp u2 δs (2.46) An interesting note on the linearized hydrodynamic moment is that the moment function behaves as a linear compliance. Combining terms to form the linearized moment compliance parameter, the function is simplified to a linear impedance element as shown by equation (2.47). Mf,l = ko θs ⇒ ko = ρAs Cs lp u2 2.9.3 (2.47) Summary of Submarine Dynamic Modeling The preceding section overviewed the derivation of the submarine vehicular dynamics as well as the selection of the submarine model parameters. With the submarine vehicle model derived, the following chapter extends the model basis to include not only the submarine and control surface hydrodynamics, but also the actuation system dynamics. Integration of these distinct system models provides a path for submarine maneuver simulation for characterization actuator and energy storage requirements. This conventional approach for requirement derivation relies heavily on a detailed model basis, hence the effort placed on model derivation. 28 Chapter 3 Control Surface System Analyzation 3.1 Legacy Hydraulic System Design of the hydraulic actuation systems located in modern submarines has remained relatively constant since the ships of the Second World War. Hydraulics remain the mainstay power source for most high force actuation demands such as the control surfaces and rudders. Additionally, the hydraulic systems operate the periscopes and masts, as well as many hull valves that require large forces to operate against the high pressures of the ocean. Rather than multiple power plants and compressors providing a supply of pressure to individual actuators, the hydraulic power plant is centrally located and provides a source of high pressure for all actuation demands. Rather than the original designs of the Second World War submarines which kept line pressures at approximately 1500 psi, modern submarines maintain a constant higher pressure around 3000 psi [31]. The higher constant supply pressure in conjunction with pressure accumulators, as seen in Figure 3.1 lessens the cycle demand on the hydraulic pumps, which significantly improved the silence of the vessels. Additionally due to the higher system pressures, the hydraulic ram actuators that drive the control surfaces are smaller [32]. The most prevalent layout of the actuation systems on modern submarines place the hydraulic components inside the pressure hull, see Figure 3.2 with the 29 Constant pressure supply line (~3000 psi) Return line Primary Accumulators Supply tank Pumps Figure 3.1: Legacy hydraulic system supply lines and primary accumulator energy storage [32] power transmission in the non-pressurized hull [1]. Hydraulic force action is produced by one or two hydraulic rams. These hydraulic rams position the control surface structures using displacement feedback, driven by electronically control hydraulic control valves. In the case of systematic power loss, manual control of hydraulic valves for positioning is also possible via a secondary manual control hydraulic loop. Manual override control draws energy from the hydraulic accumulation system, since the hydraulic pumps are driven by the electrical system [32]. The high pressure accumulators provide large capacitance for the pressure supply lines which allows the hydraulic pumps to remain off for considerable lengths of time. When the pressure supply in the hydraulic accumulators dwindles, the pumps cycle on and replenish the storage pressure. In the event of a power failure, the hydraulic systems operate independently of the central electrical power systems. With the energy storage provided by the hydraulic accumulators, the vessel can surface by the means of an emergency surfacing maneuver [32]. Several control surface and hull fin location designs exist for legacy systems. However, the most common design and the focus of this thesis is the standard layout of two aft rudders, two stern planes and two bow planes [1]. Additionally, the control surfaces have no standardized design. Permutations exist that include spade surfaces (all movable surface), skeg surfaces (fixed skeg with a single movable 30 Spade Control Surface Hydraulic Rams Tiller Arms Push Rods Pressure Hull Linear to Rotational Conversion Non-pressure Hull Figure 3.2: Legacy hydraulic system and power transmission to control surfaces. [1] rudder), or multiple flap designs. For this thesis, the control surface layouts are always assumed a spade surface design, as seen in Figure 3.3. While the hydraulic legacy systems have proven to be a reliable design for over 50 years, hydraulic systems present many sustainability and manufacturing issues. Due to the alignment critical components such as the linkages, pushrods, guide tubes, bearings, and hydraulic rams, the initial construction and continued maintenance cost remains steep [4]. The location of the hydraulic systems within the pressurized hull also restricts many variations on design. A potential outcome of a systematic change of the actuation systems would be the ability to locate the actuations systems outside of the pressurized hull. This could drastically reduce the complexity of the power transmission systems, which lowers life cycle costs and potentially decreases the audibility of the actuation. The hydraulic systems also require piping and ancillary systems that contribute significantly to construction and maintenance costs on top of the additional electrical infrastructure [4]. Separate operational systems also call for additional crew training and extra specialists for repairs and maintenance of hydraulics. While the hydraulic systems have preformed adequately for many years, there remains potential for design improvements on the existing actuation systems. 31 Figure 3.3: Standard submarine spade actuation surface for directional control [26] 32 Control Surface Submarine Pitch Motion - ωb Actuator Frame Control Surface Angular Velocity - ωcs Drive Shaft Stiffness - K B1 Linearized Hydrodynamic Compliance - Ko B2 Base (Submarine) Surface Rotational Inertia - Js Actuator Rotational Inertia - Jm Figure 3.4: Simplified layout of control surface system. 3.2 Control Surface Actuation System Modeling With a derived model of the submarine vehicular motion, the control surface input dynamics are now investigated. The two primary control force inputs to the submarine system are the control surfaces and the propulsion propeller. For this thesis, the input controls of interest are the control surfaces. Section 5.2 contains a complete derivation of the two-port control surface system equations. Before implementing the synthesis procedure to find potential design candidates, a simple closed-loop control surface control system serves as the baseline case for which all subsequent design iterations can be compared. The simple closed-loop case also represents a familiar design methodology for control systems especially when lacking expert knowledge of the underlying problem and/or special energy needs [26]. Simply closing the control loop does not address the issue of energy storage needs or realization of the physical components. Due to the relative simplicity of the control surface system model, the system equations relating the input torque from an actuation source to the control surface deflection are easily derived from Figure 3.4. This control surface system model assumes the actuation system will couple into rotational inertia, Jm , which is connected to the control surface by a drive shaft having stiffness, K. On each end of the drive shaft, bearing losses are accounted for by linear frictional damping coefficients B1 and B2 . The control surface is lumped as a single inertia with moment of inertia, Js . Additionally, the linearized ocean compliance, from equation (2.47) is 33 Sf : ωb 0 Se : Ta Ta R : b1 C : k/s 1 0 I : Jm s R : b2 C : ko /s 1 ωcs I : Js s Se : Td Figure 3.5: Bond graph form of the control surface system with ideal torque actuator. also incorporated into the two-port model as another applied torque on the control surface inertia [29]. The summation of the net acting forces in the control surface system produces the underlying dynamic equations of the control surface deflection as related to the input torque provided by the ideal rotational actuator. Since force components are assumed to be linear including the ocean hydrodynamic moments, the resulting equation is a linear ordinary differential equation. Figure 3.5 contains the topology of the energy flow in the control surface system in bond graph form. The parameters for this bond graph are given in Table 3.1. These values are calculated to model the motion of the control surface system in a Virginia class submarine. Values which were unavailable for simulation are estimated to correspond with known performance response of control surface systems [1]. Table 3.1: Control surface system parameter list Parameter Value Description Js 1.16 × 104 kg-m2 Control surface rotational inertia 1.67 × 103 N/s Internal bearing resistance B1 B2 2.26 × 104 N/s External bearing resistance and drag K 1.91 × 107 N-m/rad Control rod stiffness Jm 252 kg-m2 Actuator rotational inertia Ko 3.54 × 107 N-m/rad Linearized hydrodynamic coefficient Using the bond graph, the dynamic equations for the control surface system are derived. The simplified model provides equation (3.1), a differential equation 34 that models the control    ḣs 0     θ̇s   0     ḣ  =  −b /J 1 m  m   θ̇k 1/Jm surface movements and forces.    k −b2 /Js ko hs       0 0 1/Js    θs  +     −k 0 0   hm    0 −1/Js 0 θk 1 0 0     Td 0 0 0     ωb  (3.1)  0 0 1  Ta 0 0 0 Due to the linearity of the control surface differential equation in equation (3.1), the differential equation can be written as a transfer function relating the input torque to the output surface deflection. Equation (3.2) illustrates the control surface dynamics in transfer function form. s ωcs = 4 3 Ta n4 s + n3 s + n2 s2 + n1 s − ko (3.2) where the denominator coefficients are defined by, Js Jm k b1 Js b2 Jm n3 = + k k Jm k o b1 b2 + Js + Jm + n2 = k k b1 k o n 1 = b1 − +b k n4 = (3.3) Integrating equation (3.2), which alters the output to position rather than angular velocity, reveals that the current model is a fourth order system. Due to the high stiffness factor for the control rod, a second order high frequency harmonic would also be present in any simulations. For this derivation, the drive shaft stiffness K is allowed to approach infinity. This assumption simplifies the control surface model by eliminating high frequency harmonics from the drive shaft compliance and creates a straightforward model which can be more readily analyzed. By taking the limit as K approaches infinity, the new transfer function model of the system can be found, s ωcs = 2 Ta [Js + Jm ] s + [b1 + b2 ] s − ko (3.4) This simplified model accounts for only the second order motion of the control 35 Flow-induced Torque Disturbances Td Reference Angle + δr,cs - Control Surface Position Controller + Tact + Tnet Control Surface System ωcs 1 s Control Surface Angular Position δcs Figure 3.6: Block diagram layout for closed loop control of control surface system. surface inertia rather than dealing with high frequency harmonics. For this study on energy consumption and control, the elimination of harmonic effects is negligible. For the derivation of the closed loop control system, the linearized hydrodynamic Hoerner term provides a linear, time-invariant system. However, for the energy consumption simulations, the full non-linear hydrodynamic equations are used for a more accurate assessment of energy demand during actuation. With a system model, the next step in the conventional design process is to implement a controller scheme. Since all the poles of the derived transfer function have negative real values (in the left half plane), the control surface begins as a stable system and the system model is controllable by definition. Adding simple PID (Proportional, Integral, Derivative) control to the system to close the loop should therefore preserve the stability. Figure 3.6 shows an overview of the closed loop system in block diagram form. Closing the loop and implementing a PID controller scheme produces the transfer function shown in equation (3.5). Kd s2 + Kp s + Ki δcs = δref [Js + Jm ] s3 + [b1 + b2 + Kd ] s2 + [Kp − ko ] s + Ki (3.5) By adjusting the gains of the PID controller via Ziegler-Nichols tuning, the response of the control surface system is adjusted to meet the requirements as stated in Section 1.1. The gain parameters found are Kp = 2.99 × 106 , Ki = 7.59 × 105 , and Kd = 1.13 × 105 . The Bode plot shown in Figure 3.7 shows the frequency domain response of the control surface system. With the derived controller for the control surface system, the power requirements for the actuation system as designed in the conventional method can be estimated. By simulating extreme circumstance maneuvers of the control surface 36 Frequency Response - Open Loop Comparison 100 Magnitude (dB) 50 0 -50 -100 -150 -200 0 Plant Only PID Controller and Plant Phase (deg) -45 -90 -135 -180 10 -2 10 -1 10 0 10 1 10 2 10 3 Frequency (rad/sec) Figure 3.7: Frequency domain response comparison for PID tuned controller and open loop system. 37 system, the torque and angular velocity of the ideal actuation are used to size potential actuation candidates. A potential actuator must meet the demands of the most rigorous environments and maneuvers. 3.3 Actuation Power Limitations and Power Curve To further extend the control surface actuation system model, a power limiting curve is incorporated in the actuation model. All real actuators can only output a limited power, typically expressed by the manufacturer [33]. At high speeds the actuator is limited only by the total power of the system normally. However at low speeds and holding motions, the actuator is limited by the maximum torque output rather than the power. These limitations are especially true for electrical actuators where holding torques/forces remain one of the limiting factors. The power limitation curve is typically expressed as an inverse function of torque as a function of angular velocity as in equation (3.6). Ta (ωa ) = Pmax ωa (3.6) where Pmax is the constant maximum power of a given actuator and wa is the velocity of the actuator. During unimpeded holding, the above power limit function provides infinite torque at ωa = 0. To prevent this unrealizable actuator, a maximum torque is specified which lowers the potential power at lower velocities. Note that any torquespeed combination that lies beneath the power limit curve is physically realizable. By shifting the power limit until all simulated extreme maneuvers are realizable, the resulting power limit serves as an ideal estimate for a candidate actuator. For the case of the submarine control surface systems, the power limiting curve was determined by the specifications available in the ONR report summarized in Table 1.1. The specifications call for an operating torque capacity of 950 kNm with a preferred objective of 1,360 kN-m with an approximate tiller (moment) arm of 1.2 meters for linear actuators [1]. Additionally, the specifications identify a rate of angular displacement of 5 deg/s (0.087 rad/s). For over specification of the actuation device, the expected angular rate is moved to 10 deg/s (0.175 rad/s) [1]. At this angular velocity, the actuator can maintain the maximum torque load whereas at higher rates, the torque output diminishes. Combining these two terms 38 Actuator Power Limit 1400 Torque Output (kN-m) 1200 Maximum Power Curve Pmax = 237 kW 1000 800 600 400 200 0 5 10 15 20 25 30 35 40 Angular Velocity (deg/s) Figure 3.8: Power limiting curve for the control surface actuation system. provides the maximum power needed to drive the control surface system, 237 kW. Plotting the first quadrant torque-speed curve yields Figure 3.8. Figure 3.6 shows an overview of the closed loop system in block diagram form. The restrictions imposed by the power limitation incorporate into the model of the actuation system. The torque output of the actuator model is combined with the actuator velocity to determine the current actuator power demand. The actual system torque output is determined by a power check against the power curve in Figure 3.8. If the power exceeds the curve or the torque exceeds the torque limitation, the simulation stunts the final actuation power output. 3.4 Application of System Model with Submarine Model With the model and control system for the control surfaces, the entire submarine dynamics forms a model basis for the motions of the submarine in addition to the dynamics of the control surface systems. Rather than instantaneous deflections due to step inputs, the physics of the control surface inertia and feedback loops include the transients of the system. Additionally, the motion variables from the 39 submarine dynamics equations allow for proper calculation of the hydrodynamic moments encountered by the control surfaces. 3.4.1 Control Surface Load Estimation via Submarine Model With a model basis, simulation based design provides preliminarily insights into expected torque loads and potential energy storage demands. Of particular interest are the loads and energy required to perform sustained extreme maneuvers as well as an emergency ballast blow maneuver [1]. For a typical submarine mission, the actuation systems encounter the greatest hydrodynamic torque demands during sharp yaw or pitch banks [32]. Particular attention remains on the stern plane actuation system due to the location of the surfaces distally and external to the pressurized hull and due to the large moments required for actuation. The following two simulations explore the torque loads encountered by the stern plane actuation system in both rapid peak to trough deflection as well as the classic ballast blow maneuver. To demonstrate prediction of the dynamics of a submarine performing extreme pitch bank maneuvers, a full sweep and hold of the stern plane actuation system has been simulated. The deflection of the stern plane is swept from a peak of +30 degrees to the trough of -30 degrees. After a complete sweep of the range of the actuator, the deflection is held stable at -30 degrees. This maneuver is analogous to a rapid rising and falling of the submarine during operation without a change in ballast pressure. Additionally, the ONR utilizes this basic test for the quantification of actuation systems [1]. 3.4.2 Full Surface Sweep Maneuver Simulation With the combined submarine and actuation model as described in the previous sections, the full sweep and hold simulation can accurately estimate the torque loads encountered during this type of maneuver. In this particular simulation, the submarine travels at 15 m/s (30 knots) surge velocity representing near top speed for the Virginia class vessels [1]. The system thrusters attempt to maintain the constant velocity surge during the full sweep and hold and all other velocities begin with zero initial conditions. Additionally, a Pierson-Moskowitz spectrum is used to shape expected flow-induced disturbances [34]. Figure 3.9 shows the deflection of the control surface during the maneuver and displays the actuator torque requirements. 40 Control surface transient deflection Actuator Power Consumption (kW) 100 Deflection Reference 20 Power (kW) Deflection (deg) 40 0 −20 −40 0 20 40 60 Time (s) 80 100 0 −100 −200 120 0 Control surface moment during emergency surfacing Energy Consumption (kJ) Moment (kN−m) 500 0 0 20 40 60 Time (s) 80 100 40 60 Time (s) 80 100 120 Actuator Energy Consumption (kJ) 1000 −500 20 120 1500 1000 500 0 0 20 40 60 Time (s) 80 100 120 Figure 3.9: Control surface deflection and dynamic demands during full sweep simulation. As can also be seen from Figure 3.9 during the holding maneuver, a torque demand continues. From the persistent demand for torque actuation, the final actuation design must be capable of four quadrant operation, meaning both actuation and braking. In addition, Figure 3.10 conveys the important submarine model parameters of interest. 3.4.3 Emergency Surfacing Maneuver Simulation The other simulation of interest is the emergency surfacing maneuver [1]. Unlike the continuous sweep simulation above, the stern planes in the surfacing maneuver perform precision tracking to orientate the submarine to the proper global angle (θ) so that the vessel exits the water at the proper angle for structural integrity. Additionally, with the large change of buoyancy after blowing the ballast chamber, the stern plane actuators must ensure that the ascent is controlled rather than allowing the submarine system to destabilize. In the worst case, an uncontrolled ascension could exit the water perpendicularly causing bodily harm to the crew and severe damage to the structural systems. For the emergency surfacing maneuver simulation, the buoyancy doubles with respect to the initial value equivalent to the weight of the vehicle. For reference, the initial weight and buoyant force of the submarine are 7.06 × 105 N and after the 41 Submarine States during Full Control Surface Sweep Global Submarine Z Position Z Position (m) 1000 500 0 0 20 40 60 80 100 Time (s) Global Submarine Heading − Theta 120 100 0 −100 50 0 −50 0 20 40 60 80 Time (s) Heave Velocity − w 100 40 60 80 Time (s) Surge Velocity − u 100 120 0 20 40 60 80 100 Time (s) Pitch Angular Velocity − q 120 0 20 2 q (deg/s) w (m/s) 20 10 0 120 1 0 −1 0 20 u (m/s) Theta (degrees) X Position (m) Global Submarine X Position 0 20 40 60 Time (s) 80 100 0 −2 120 40 60 Time (s) 80 100 120 Figure 3.10: Selected submarine states for full sweep maneuver. ballast has filled with air the value of the buoyant force is 1.41 × 106 N. To deflect the submarine to the desired submarine global angle of 20 degrees, the stern planes begin initially deflected -30 degrees [23]. Another control loop for submarine angle control directs the control surfaces to drop to near zero deflection after reaching the target submarine angle. Once at the desired submarine angle and the proper buoyant force, the control surfaces act only reject disturbance and stabilize the path of the vessel. For further realism the Pierson-Moskowitz ocean spectrum is used to inject ocean disturbance torques to both the actuation system as well as the submarine dynamics. The important states of the submarine dynamics are given by Figure 3.11. Since the emergency rise maneuver only occurs in the vertical plane, the transverse plane parameters are unaffected by the maneuver and remain near zero. As can be seen by Figure 3.11, the submarine angle pitches to 20 degrees due to the control surface deflection. Additionally, the heave velocity rapidly raises due to the additional buoyant force. With the proper control surface actuation, the rapid rise remains fully stable. Figure 3.12 shows the deflection of the control surface during the emergency surfacing maneuver and displays the actuator torque 42 Global Submarine X Position Z Position (m) 400 200 0 0 10 20 30 Time (s) 40 50 60 0 10 30 Time (s) 40 50 q (deg/s) 0 10 20 30 Time (s) 40 10 50 30 Time (s) 40 50 60 0 10 20 30 Time (s) 40 50 60 50 60 Pitch Angular Velocity - q 0 -2 60 20 Surge Velocity - u 2 2 0 0 5 0 60 Heave Velocity - w 4 w (m/s) 20 0 10 20 0 200 Global Submarine Heading - Theta 40 Global Submarine Z Position 400 u (m/s) Theta (degrees) X Position (m) Submarine States during Emergency Surfacing 0 10 20 30 Time (s) 40 Figure 3.11: Selected submarine states for rapid ascension maneuver. requirements. As can be seen from Figure 3.12, the toque requirements never completely return to zero due to the additional hydrodynamic forces induced by motion. The actuation system encounters the largest torque requirements during the return actuation move after the submarine reaches the 20 degree angle. Since this maneuver is critical for maintaining stability, the second main actuation movement is the most difficult and important. Any energy storage must ensure enough energy to power both the first and second maneuvers, as not executing the second maneuver destabilizes the system. If the second actuation were to fail, the submarine would exit the ocean near vertically. With the known energy demands of these extreme actuation maneuvers, the energy demands yield sizing of potential energy storage devices. The following section discusses the amount of energy storage required for electromechanical actuation systems for the submarine application and potential energy storage mechanisms. 43 Control surface transient deflection Control surface moment during emergency surfacing 100 0 Moment (kN−m) Deflection (deg) 10 −10 −20 Deflection Reference −30 −40 0 10 20 30 Time (s) 40 50 50 0 −50 −100 60 0 Actuator Power Consumption (kW) Energy Consumption (kJ) Power (kW) 10 5 0 0 10 20 30 Time (s) 40 50 20 30 Time (s) 40 50 60 Actuator Energy Consumption (kJ) 15 −5 10 60 80 60 40 20 0 −20 0 10 20 30 Time (s) 40 50 60 Figure 3.12: Control surface deflection and dynamic demands during emergency surfacing. 3.5 Energy Storage Requirements for Actuation The previous section demonstrates the capabilities of simulation based design to provide insight into energy demands on actuation systems in extreme circumstances. With known energy requirements for these maneuvers, energy storage devices can be sized according to simulation results [4, 33]. However, note that although the simulation design technique relies on simulation-based techniques, the design process remains an ad hoc expert based process. Further design techniques such as the synthesis technique, which is investigated next, produces potential design based purely on analytical results. With the simulation results, potential energy storage devices and their capacity are presented. For storage of energy to power an electromechanical actuation system, many suitable technologies exist [35]. Figure 3.13 shows many of the potential energy storage technologies that could provide emergency power and regulate the voltages provided to the system electronics. Batteries provide energy storage by a conversion of the electrical potential to a chemical potential. On the other hand, the flywheel physically stores the potential energy in a rotating wheel requiring a conversion from electrical to rotational. 44 Ragone chart - Energy vs. Power delivery profile Enegy Density (Joules/kg) 10 6 Flywheels Batteries 10 5 Ultra Capacitors 10 4 Metal Oxide Capacitors 10 Carbon Capacitors 3 10 2 10 3 10 4 10 5 10 6 Power Density (Watts/kg) Figure 3.13: Comparison of energy density to power density of energy storage devices. [35] Capacitors store potentially energy in an electric field, directly storing charge on separated plates. Each of these technologies exhibit pros and cons, so final decision on the technologies requires a thorough comparison of energy and power density requirements, weight/balance issues, noise production in the submarine setting and the ease of implementation. The energy consumption predicted via simulation for the single actuator can support systematic selection and sizing of energy storage requirements. In the case of the full sweep simulation, a single actuator consumed 493 kJ and the emergency surfacing maneuver required 80 kJ. In the event of power failure, additional energy may also be needed in case actuations are incomplete. An askew control surface, for example, could potentially drive the entire submarine system into a state of instability. For a single actuator, the energy storage should be five times the encountered peak for an extension to the amount of time a maneuver could persist. The resulting predicted energy storage is 2465 kJ. After acquiring the estimated energy consumption, energy storage selection requires additional knowledge about the peak power draw. Comparing the peak power draw of 25 kW from the full sweep simulation to the estimated required energy storage of 2465 J, the energy storage device could take the form of batteries, ultra capacitors or flywheels in accordance with Figure 3.13. The energy storage 45 mechanism must contain a large energy density in comparison the relative power density of the system. Further investigation into the preferred energy storage for submarine applications must be conducted. The location of the energy storage device in the system topology also remains an ad hoc procedure [4, 35]. A potential design could have the energy storage in a centralized location much akin to the current hydraulic accumulator designs. However, centralized energy storage requires power transport and distribution that add losses to the system. Another possibility is to have localized energy storage for each actuation system. This setup provides quick response times but once the power in one area is exhausted, energy from other underutilized actuators proves difficult to relocate quickly to drained systems in the event of power failure. Also note that the voltage type of the ship power grid is also a crucial factor. Due to the strong push towards medium voltage DC power, the following design iteration utilizes medium voltage DC power systems [1]. Figure 3.14 illustrates the centralized energy storage design. In the centralized system, power buses distribute power amongst the various electrical tasks. For the ad hoc realization presented, power conversion systems lie in distributed positions for direct power supply to elements. Additionally in this case, a large central battery bank provides an energy storage mechanism that can recharge and discharge. Another potential design iteration, illustrated by Figure 3.15, provides individual energy storage mechanisms for each large actuation source. In this case, due to the decreased need for large energy density and the desire for quick response times, ultra capacitors supply power directly for energy demands. For additional systematic energy storage, a centralized battery bank provides additional power to the electrical grid in event of power failure. While the ad hoc design methodology produces feasible designs, the topology and energy storage mechanisms are not necessarily ideal for the submarine system. Completion of the design process requires expert and/or ad posteriori knowledge of energy systems as well as the submarine specific application. The system simulations provide insight into actuation energy requirements for physical submarine maneuvers which aid in energy storage sizing. In the following sections, the synthesis approach produces not only the required energy storage for a given actuation task, but also provides the topology of the energy storage mechanism. Rather than 46 Large Energy Storage Bank (Ba!eries) Power Source Power Conversion Control Surface Actuator Ballast Valve Actuator Mast and Periscope Actuator Control Surface Actuator Figure 3.14: Potential centralized energy storage and distribution layout for submarine system. Ultra Capacitor Power Electronics Control Energy Storage Bank Power Source Power Conversion Control Surface Actuator Ballast Valve Actuator Mast and Periscope Actuator Figure 3.15: Another potential energy storage and decentralized layout for submarine system. 47 ad hoc placement of components in a given system, a direct identification of both component location and size is found. 48 Chapter 4 Overview of Two-Port Synthesis This chapter introduces two-port systems and the application of synthesis techniques to these structures to obtain desired performance specifications. Two-port representation of systems has existed in the field of electrical engineering and only more recently have these techniques been applied to more general multi-domain problems. Additionally, utilizing two-port synthesis techniques in domains other than the electrical realm has only recently been applied more broadly. Synthesis techniques in the two-port domain add a layer of complexity to synthesizing results due to the possibility of extremely complex or non-physically realizable systems. The following chapter discusses in further detail the benefits of two-port synthesis and the techniques to deal with realizations of results. 4.1 Essential Concepts for Synthesis Theory The dynamics of complex engineering systems with multiple inputs and outputs can be modeled using a variety of techniques and formulations. However, the ability of many techniques for analysis of a variety of parameters and system definitions lacks both flexibility and robustness [16]. The most complete over-arching method of capturing system dynamics is a n-port system matrix representation [36, 37]. Using n-port system representation, the relationships between any number of inputs to any number of outputs can be analyzed individually which grants major flexibility to system analysis and control. 49 i2 i1 e1 v1 Two-port Network Two-port Network f1 v2 e2 f2 Effort = e = {voltage, force, torque, pressure} Flow = f = {current, velocity, angular velocity, flowrate} (a) (b) Figure 4.1: (a) Electrical two-port. (b) Generalized (multi-disciplinary) two-port with bonds conveying power, Pi = ei fi . 4.1.1 Description of Two-port systems More generally, the vast majority of complex systems can be represented as only a two-port system with only two inputs and outputs, as seen in Figure 4.1 [20]. In two-port systems effort and flow variables typically represent power variable pairs such as voltage-current, force velocity, etc. Consequently, any design problems including power transmission can be effectively modeled as an aggregate of two-port subsystems [13]. Another pleasing aspect of n-port systems (or the 2-port explicit case) is that simple matrices represent the models of a complex linear systems. This systematic approach improves the derivation speed of synthesis results and provides the ability to automate the synthesis methodology. For the simple two-port network in Figure 4.1, a potential grouping of parameters in matrix form is shown in equation (4.1). " e1 e2 # = " z11 z12 z21 z22 #" f1 f2 # (4.1) where the z parameters are impedance functions. Due to the system linearity, the matrix equation may also be represented V = ZI. The parameters of the n-port matrix models are interconnected by the linear relationships between the input and output parameters. Individually the parameters in the n-port matrix model represent the system in ideal conditions. For example in the electrical two-port in Figure 4.2, if i2 = 0 (second port is shorted) then z11 is 50 i2 = 0 + v1 i1 Two-port Network - v2 Figure 4.2: Finding z11 for a 2 port electrical network. found to be the ratio of v1 to i1 as shown by equation (4.2). z11 = v1 i2 = 0 i1 (4.2) These linear relationships hold for all systems including n-port multi-disciplinary systems. Rather than electrical networks as in the simple example above, the voltage measurement could be any effort variable (force, torque, pressure, voltage) and likewise the current measurement could be any flow variable (velocity, angular velocity, flow rate, current). 4.1.2 Canonical Matrix Forms and the Corresponding Parameters In addition to the impedance matrix form shown in equation (4.1), early developers of n-port matrix theory designated certain matrix forms as “canonical”. A detailed derivation of two-port matrix forms can be found in Huelsman [13] and the matrix forms are summarized in Figure 4.3. The synthesis methodology heavily relies on the canonical matrices and transformations between the various forms [38]. Since the primary matrix forms that are dealt with in the following sections are transmission matrices and immitance matrices, the following paragraphs outline these parameter forms in detail. Transmission matrices are heavily utilized in the derivation of models using bond graph impedances which is covered in this section. The transmission matrix allows conversion of impedance-based bond graphs (or networks) to symbolic matrix relations. The transmission matrix is written, " e1 f1 # = " n11 η12 y21 n22 51 #" e2 f2 # (4.3) e1 Transmission - M f1 e1 Impedance - Z e2 e1 Immittance - H f2 f1 Adpedance - G e2 f1 Admittance - Y f2 n11 η12 e2 y21 n22 f2 = η11 η12 η21 η22 f1 f2 = η11 n21 n12 f1 y11 e1 n21 n12 η22 y11 y12 e1 y21 y22 e2 = = = y22 e2 f2 [n ] dimensionless [η ] [e / f ] [ y ] [ f / e] Figure 4.3: Two-port system canonical forms with descriptions of transfer function relationships. 52 In the notation above, the n elements represent dimensionless quantities, the η elements represent impedance relationships (e/f ) and the y elements represent admittance (f /e). Each of the parameters are defined as, η11 = e1 f1 |e2 η11 = e1 f1 |e2 =0 n12 = e1 e2 |f1 n12 = e1 e2 |f1 =0 (4.4) =0 =0 As is discussed further in this section, the transmission matrix product for each 0 and 1 junction are used to determine the system transmission matrix. For system insertion, unknown impedance components can be added into the system model. Solving for the unknown impedance functions is the basis for the synthesis techniques [19]. Another important canonical matrix form is the immittance matrix which relates the effort of one input and the flow of the other to the converse flow and effort output. Since the derived models contain both a flow source and effort source, this form is required for synthesis. The immitance matrix has the form, " e1 f2 # = " η11 n12 n21 y22 #" # f1 e2 (4.5) where, η11 = e1 f1 |e2 η11 = e1 f1 |e2 =0 n12 = e1 e2 |f1 n12 = e1 e2 |f1 =0 (4.6) 4.1.3 =0 =0 Conversion between Canonical Forms Since the matrices remain in canonical form, transformation between any of the forms is a simple matter of algebraic manipulation. A table in Figure 4.4 summarizes transformations between the canonical forms of two-port matrices [13]. Using transformation relationships, the system matrix can be defined in multiple ways allowing many different transfer function relationships to be analyzed. The classic field of electrical systems primarily utilized two-port representation since the two-port representation accurately modeled the bulk of the common electrical systems. These system techniques have been extended to any physical domain for two-port representation. Due to the extension to multiple physical do53 A B e2 C D f2 = 1 A −∆ C 1 −D f1 = 1 B ∆ D 1 −C = 1 C A 1 = 1 D − ∆ e1 B 1 − A e2 e1 Transmission - M f1 e1 Impedance - Z e2 e1 Immittance - H f2 f1 Adpedance - G e2 f1 Admittance - Y f2 = f2 f1 e2 ∆ e1 −B f2 Figure 4.4: Relationships between canonical two-port systems. Conversions are all with respect to the transmission matrix, which can be easily derived from the system bond graph. 54 e f e f e f R R C e f C f ( s) =G e( s ) e f e( s ) =R f ( s) e f f ( s) = Cs e( s ) e f e( s ) 1 = f ( s ) Cs e f f ( s) 1 = e( s ) Is e( s ) = Is f ( s) I I n T r G Z e( s ) = n2Z ( s) f ( s) Z f ( s) = r 2 Z ( s) e( s ) Figure 4.5: Basic elements of bond graph notation with causal definitions. The casual definition implies the transfer function of an element. mains, the method of generating the two-port system matrices involves the use of bond graph topology. 4.1.4 Impedance Bond Graphs While synthesis and impedance methods are still applicable without bond graphs, the energy flow system ideally allows the designer to identify individual elements and the respective topology from the system as a whole [15]. Bond graph elements represent the basic linear or non-linear multi-disciplinary elements of realizable systems. For those unfamiliar with bond graphs, Karnopp [15], Redfield [16, 17, 39] and Connolly [18, 19] contain overviews of the subject. Since this synthesis approach uses impedances and the frequency domain, the reader needs only to be familiar with impedance-based bond graphs. Similar to impedance electrical networks, each bond graph element corresponds to a simple impedance function. Note, however, the causality of the bond graph alters how the impedance function should be interpreted. Figure 4.5 presents a summary of the potential bond graph elements with equivalent impedances based on the causality definitions of the bond. With the impedances of bond graph elements defined, the method of impedance model derivation can be derived using transmission matrix form. Bond graph junc- 55 e1 f1 1 e2 e1 f2 f1 = 1 Z e2 0 1 f2 1 0 e2 1/Z 1 f2 m e2 Z(s) e1 f1 0 e2 e1 f2 f1 = Z(s) e1 f1 e1 f1 T G e2 e1 f2 f1 e2 e1 f2 f1 = = 0 0 1/m f2 r e2 1/r 0 f2 0 Figure 4.6: Two-port transmission matrix forms for bond graph junctions tions also play a pivotal role in transmission matrix formulation. Zero junctions have common efforts, so a two port matrix correspondingly has a zero for the η12 parameter. One junctions have common flows, so a two port matrix likewise has a zero for the y21 parameter. The transmission matrix forms for common bond graph junctions are summarized in Figure 4.6 [38]. 4.1.5 Description of Synthesis Impedance-based bond graphs and n-port matrix models form the basis of the mathematical tools for the synthesis procedure. Synthesis methods can generate preliminary designs based on a desired system response between two variables of interest. Once an impedance-based bond graph has been derived for a given model, the de56 ω1 + T1 ωa - ω1 ω2 + - ω2 + T2 ωa Za(s) T1 ω1 0 T2 ω2 Ta ωa - Za(s) (a) (b) (c) Figure 4.7: (a) Electrical analog for ideal torque actuation. (b) Electrical analog showing active impedance replacing torque source. (c) Bond graph form of system in (b) signer inserts an unknown impedance function to represent system alteration or retrofitting. For example, an electrical analog model for the ideal torque source, Ta (t), is shown in Figure 4.7 (a). Retrofitting the torque source, the designer can remove the known ideal torque source and replace with an unknown impedance, Za (s) = Ta /ωa , shown in Figure 4.7 (b) as electrical analog or as in (c) using an equivalent bond graph form. The synthesis methodology determines the unknown impedance function in terms of specified system response requirements. As a result of the specified system response and system structure, this function may be purely passive, purely active, or a combination. The feasible realizations of the impedance functions result from the real component characteristics of the complex impedance. Positive real functions exhibit the most flexibility and are potentially realizable as purely passive system. Appendix A outlines the properties of positive real functions and the scenarios from which simple pure passive realizations result. Expansion of the resulting impedance function yields the individual impedance elements and the topology of the component interconnection. Standard methods of function expansion such as partial fractions and continuous fractions, amongst other techniques, represent common starting points for elements identification. For positive real systems under certain conditions discussed in Temes [40], systems composed of only RC, RI, and CI elements are realizable via Cauer or Foster synthesis. A positive real impedance function can always be realized as a RCIT system (R-element, 57 C-element, I-element, and transformers) using Brune synthesis shown in Temes [40] or Bott-Duffin synthesis to eliminate transformer coupled I-elements [41]. However, in the case of positive real function with pole removal at only s = 0 or s = ∞, a minimal physical realization is possible with only RCI elements [40]. Synthesis of multiport systems typically requires topology restrictions on the unknown impedance compensator. In the case of two-port systems, a two-port unknown impedance function could produce myriad answers many of which would be difficult to physically realize [20]. Restricting the form of the unknown impedance function to match the original connection topology (in shunt or parallel) reduces the number of answers produced by synthesis. A summary of the synthesis procedure is shown in Figure 4.8. Additionally, the content of section 4.2 demonstrates the synthesis procedure entirety. 4.2 Active Vehicle Suspension Example As an example of the two-port synthesis technique as utilized for realization of physical systems, the classical example of an active suspension system was chosen. The active suspension example suits the synthesis procedure since we are analytically searching for actuation impedance which will produce the desired response. Additionally, the breadth of solution methodology for this example in the literature is broad [19, 20, 42, 43, 44]. Verification of the system response via other experimental results or by intuition is also easy for the active suspension example. Analyzation of a four wheel active suspension system shows both medial and transverse symmetry so that the model reduces to a single active suspension system [42]. For this example, the model of the active suspension system will remain relatively simple with a non-complaint massless tire, as seen in Figure 4.9. Also, the mass of the quarter vehicle is lumped into a single block element with mass mp = 613kg. Finally, the connection element between the tire/road input velocity and the vehicle mass is unknown and represented by a black box. 4.2.1 Bond Graph Model and Matrix Derivation To begin the synthesis process, a model of the system including the unknown element must be found. As discussed previously, the models in this thesis will be constructed in two-port form and also represented with bond graphs. The bond graph model of 58 Incorporate known system into n-port bond graph model Assimilate unknown n-port impedance function Restrict form of unknown impedance accordingly Obtain desired transfer function for compensated response – dTF(s) Simulate system to find baseline response Find dTF(s) via curve fit of experimental data Form dTF(s) as modification on system baseline response Transform n-port matrix model to required canonical form Decompose Z(s) into basic impedances and set in bond graph form Equate n-port matrix to dTF(s) and solve Determine form of synthesized components (hydraulic, electrical, etc.) Incorporate all synthesized components in control algorithm Realize passive elements as separate entities Simulate results and repeat synthesis as necessary Figure 4.8: Synthesis procedure for n-port systems to generate frequency domain compensators 59 Fd vp vehicle body - mp Z(s) vt Figure 4.9: Active suspension system with unknown impedance function for specified response. Fd is applied disturbance force, vp is the body output velocity, vt is the tire input velocity and mp is the quarter vehicle mass. 60 I : mp . vp p Sf : vt Ft vt Ft vp 0 Fd vp 1 Se : Fd v z Ft Z(s) Figure 4.10: Bond graph representation of the active suspension system including the unknown two-port impedance function which has been restricted in form to a zero junction only. the quarter vehicle suspension is represented as in Figure 4.10. For this example, the form of the two-port unknown impedance function is restricted as a single impedance function off of a 1-junction. As discussed before, this restriction limits the quantity of designs but ensures a more simple realization procedure. Using the bond graph representation, the two-port matrices of each junction are formulated using transmission matrix form. Beginning with the right side of the bond graph model, the matrices are constructed for each element and parsed into series matrix multiplication as in equation (4.7). " 1 0 1/Z(s) 1 #" 1 mp s 0 1 # (4.7) From the series of element junctions, the impedance function that represents the entire two-port system is calculated. The total function is equation (4.8). " 1 1 Z(s) mp s mp s Z(s) 61 +1 # (4.8) Desired system (Vp /V t ) frequency response -10 Gain (dB) -15 -20 -25 -30 -35 1 10 10 2 Phase (degrees) 50 0 -50 -100 -150 1 10 Frequency (rad/s) 10 2 Figure 4.11: Experimentally obtained desired system frequency response (Vp /Vt ) 4.2.2 Specify Desired Frequency Response With the model of the active suspension system derived, the response of the system is specified. The relation of interest is between the input of the wheel/road velocity on the velocity of the body (vp /vt ) with no disturbance force applied to the vehicular body. The desired frequency response is chosen from experimental data from an actual active suspension system in [19], as show in Figure 4.11. The analytical representation of this desired frequency response represented in transfer function form is equation (4.11). This function relates the input of the wheel/road velocity to the velocity of the body. dT F (s) = a1 s + a0 s4 + b3 s3 + b2 s2 + b1 s + b0 (4.9) where a1 = 77, 254, a0 = 5.11 × 106 , b3 = 12.1, b2 = 56, 230, b1 = 628, 790, and b0 = 1.26 × 108 . 62 4.2.3 Synthesized Impedance Function, Functional Expansion With the specified response determined, the canonical form of the synthesis matrix must relate the same efforts and flows of the desired transfer function. The specified transfer function relates the flow input to the flow output of the two-port system. Looking over the canonical matrix forms, the H (immittance) matrix form contains the relationship f2 /f1 given that e2 = 0. Since the two-port matrices were derived via transmission matrices, a transformation between canonical forms is needed for the synthesis procedure. Using the canonical matrix relations discusses earlier, the H-form of the system equations is found as equation (4.10). " mp sZ(s) mp s+Z(s) Z(s) mp s+Z(s) Z(s) mp s+Z(s) −1 mp s+Z(s) # (4.10) Now that both the two-port model of the system and the desired transfer function are specified and ensured to be in proper form, the actuation impedance is synthesized. In the active suspension case, we have four equations with four unknowns to solve. The elements H11 , H12 , and H22 are equivalent to unspecified unknown frequency responses and M21 is equivalent to the desired frequency response, as shown in equation (4.11). H11 ⇒ H12 ⇒ H21 ⇒ H22 ⇒ mp sZ(s) mp s+Z(s) Z(s) mp s+Z(s) Z(s) mp s+Z(s) −1 mp s+Z(s) = U11 = U12 = s4 +b a1 s+a0 3 2 3 s +b2 s +b1 s+b0 (4.11) = U22 Solving the system of equations in equation (4.11) provides not only the resulting impedance function of the unknown impedance function Z(s), but also the resulting frequency responses of the other relations in the two-port system (U11 , U12 , U22 ). The resulting transfer functions for each element of the two-port matrix are shown in Figure 4.12. As seen in the graphs, the frequency response for H21 matches the specified response. Additionally, the resulting impedance function for Z(s) is found as in equation 4.12. This impedance takes the difference in velocity of the tire and the 63 Resultant system H 11 Resultant system H 12 80 -10 -15 Gain (dB) Gain (dB) 70 60 50 40 1 10 -20 -25 -30 10 Frequency (rad/s) -35 1 10 2 -75 -15 -80 -20 -25 -30 -35 1 10 10 2 Resultant system H 22 -10 Gain (dB) Gain (dB) Specified system H21 Frequency (rad/s) -85 -90 -95 Frequency (rad/s) 10 -100 1 10 2 Frequency (rad/s) 10 2 Figure 4.12: Resulting frequency response for two-port matrix elements after system synthesis. mass and produces a torque output to the system. Z(s) = nu2 s2 + nu1 s du4 s4 + du3 s3 + du2 s2 + du1 s + du0 (4.12) where nu2 = 4.75 × 108 , nu2 = 3.15 × 1010 , du4 = 10, du3 = 121, du2 = 5.62 × 105 , du1 = 5.52 × 106 , and du0 = 1.21 × 109 . Since equation (4.12) is not positivereal, physical realization of this system will contain active elements and possibly passive elements. The impedance function as a whole can be decomposed via a multitude of expansions techniques. For this example, partial fraction expansion of the impedance function resulted in equation (4.13). Z(s) = na1 s + na0 nb1 s + nb0 + s2 + da1 s + da0 s2 + db1 s + db0 (4.13) where na1 = 5.09 × 104 , na0 = −2.07 × 106 , da1 = 10.1, da0 = 2.24 × 103 , nb1 = −5.09 × 104 , nb0 = 4.98 × 107 , db1 = 1.96, and db0 = 5.40 × 104 . Since both terms of the resulting partial fraction expansion cannot be expressed via the basic bond graph elements, continuous fraction expansion further decomposes the original function to 64 equation (4.14). Z(s) = 1 k11 s + b11 + 1 m21 s+b21 + k (4.14) 1 31 +b 31 s where k11 = 3.91 × 106 , b11 = 854.4, m21 = −15.2, b21 = −1600, k31 = −6.2 × 104 , and b31 = 1596. Note, several of the derived parameters have negative values which implies active systems as expected. 4.2.4 New Bond Graph Layout, Physical Realization The expanded impedance function contains three bond junctions with both active and passive components. The bond graph representation of the expanded impedance function is Figure 4.13. The passive elements of the impedance function can be extracted into physically realizable components and the active components are lumped into a single flow source. With the negative impedance elements lumped into a flow source, the system is realized as Figure 4.14. As discussed earlier the flow or effort source elements are initially assumed to be ideal, meaning both lossless and energetic. Physical limitations of the actuator can be added and compensated for. The passive compliance and damping elements have identical parameters as derived from the expanded impedance function. For the flow source, the output flow (va ) is related to the effort signal (Ft ) via the transfer function of the active elements. Using the state equations derived for the active elements in equation (4.15), the compensator controller function is derived. b21 m21 p21 − k31 x31 k31 1 m21 p21 − b31 x31 ṗ21 = Ft − ẋ31 = (4.15) where m21 = −15.2, b21 = −1596, k31 = −6.2 × 104 , and b31 = 1595. Representing equation (4.15) in block diagram form results in a compensator algorithm which utilizes the parameters derived from the synthesis method. Figure 4.15 depicts the controller form integrated with the bond graph structure. From this form of bond graph, the state equations for the active suspension 65 Passive Elements C : k11 R : b11 I : mp . vp p Sf : vt Ft vt 0 Ft vp 1 Fd vp Se : Fd vz Ft -R : b21 1 -I : m21 -C : k31 0 R : b31 Active Elements Figure 4.13: Synthesized bond graph elements separated into both active and passive components. 66 Fd vp vehicle body - mp damper - b11 spring - k11 Flow source actuator vt Figure 4.14: Realization of the bond graph from Figure 4.13. The active elements are lumped together into a single ideal flow source and the passive elements are represented by simple linear elements. 67 C : k11 R : b11 I : mp . vp p Sf : vt Ft vt 0 va Ft Ft vp Fd vp 1 Se : Fd . p21 Ft - 1 s - Sf : va p21 b21 m21 va k31 x31 1 s . x31 1 m21 - k31 b31 1 m21 Figure 4.15: Flow source controller algorithm as determined by the synthesis expansion incorporated in the bond graph structure. 68 Active and passive chassis velocity comparison 1 vground 0.8 vactive vpas s ive 0.6 Velocity (m/s) 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Time (s) Figure 4.16: Velocity output at 10 rad/s comparison between suspension with only passive elements and a combination of active/passive. system are derived for performance analysis. The state equation is equation (4.16). " ṗ ẋk # = " 0 k11 − m1p 11 − kb11 #" p xk # + " −1 0 0 0 1 −1 #  Fd     vt  va (4.16) where k11 = 3.914 × 106 N-m, b11 = 854.4 N/s, mp = 613 kg, m21 = −15.2, b21 = −1596, k31 = −6.2x104 , and b31 = 1595. 4.2.5 Simulation of Actuation System with Synthesized Compensator For verification of the validity of the results of the synthesis procedure for the active suspension system, simulations are performed on the system state equations. The velocities of the vehicular body due to a sinusoidal input disturbance from the tire are compared with the flow source and without. In the case without the flow source, the only elements of the impedance function that remain are the passive spring and damper. Figure 4.16 depicts the comparison of the velocity cases. 69 Magnitude Gain Comparison -10 Active/Passive Passive Only Desired Gain (dB) -15 -20 -25 -30 -35 -40 1 10 Angular Velocity (rad/s) 10 2 Phase Comparison 20 Phase (degrees) 0 -20 -40 -60 -80 -100 -120 1 10 Angular Velocity (rad/s) 10 2 Figure 4.17: Frequency response comparison between suspension with only passive elements and a combination of active/passive The flow source has a clear effect of attenuation of the input to the system over the case with only the passive elements. A more clear visualization of the difference between the cases is seen through a swept frequency response simulation comparing the passive only system and the fully synthesized system. Figure 4.17 depicts the frequency response comparison on the magnitudes in comparison with the desired system response which was specified. With an ideal linear flow source actuator, the frequency response closely matches the specified response used for the synthesis method. The compensator shapes the output of the flow actuator to provide the ideal response when coupled with the passive elements that were synthesized. For actuation sizing, the power requirements of the synthesized actuator must be known. A simple method for characterization of the actuator is to plot a force-speed curve and find the bounds of the plot. The extreme values of the power conjugates, force and speed, provide the maximum power consumption of the actuator. Figure 4.18 shows the force-speed curve of the ideal actuation element for a constant sinusoidal input of 35 rad/s. However, a force-speed curve for a single input frequency does not fully 70 Force-speed Characteristic Curve of Actuator 3000 2000 Force Output (N) 1000 0 -1000 -2000 -3000 -4 -3 -2 -1 0 1 2 3 4 Velocity (m/s) Figure 4.18: Force-speed of linear actuator at 35 rad/s input tire disturbance. characterize the power consumption of the device at all expected inputs. A more practical simulation result that fully characterizes the actuator is the force-speed response to a random input, as in Figure 4.19. 4.2.6 Hydraulic Realizations via Secondary Synthesis The previous realization of the active suspension system assumed the availability of an ideal flow sources. This ideal device could be potentially realized as an electromechanical actuator or a hydraulic actuator assuming no transformation moduli. However in reality, there exists a scaling transformation between one physical domain to another. In the case of electrical to mechanical, there is a motor or actuator constant, and for hydraulic there is the area of the piston. Realization of the hydraulic actuation case is relatively straight forward. By inserting a transformation between the flow source and the remainder of the bond graph, the system effectively models a hydraulic system with a variable flow rate pump for control. Figure 4.20 shows the simple hydraulic realization of this system. However, it should be noted that by adding a transformation between the 71 Force-speed Characteristic Curve of Actuator 2000 1500 Force Output (N) 1000 500 0 -500 -1000 -1500 -2000 -1.5 -1 -0.5 0 0.5 1 1.5 Velocity (m/s) Figure 4.19: Force-speed of linear actuator at random tire disturbance. flow provided from the pump and the flow to the system, the controller signal must be adjusted accordingly to ensure two-port impedance equivalence. The new input signal to the variable control rate pump must have an additional scaling factor of the area of the hydraulic piston to cancel out the effect of the transformation modulus, as seen in Figure 4.21. Another potential hydraulic realization of the active vehicle suspension would have the capacitance and resistive elements in the hydraulic domain with the flow source. As in the simple hydraulic realization case, the two-port impedance of the modified bond graph must be equivalent to the original synthesized bond graph for identical performance. To find the equivalent impedance of the new modified bond graph via secondary synthesis, the original bond graph transmission matrix must be equivalent. Figure 4.22 shows the original bond graph with the flow source reverted to the original negative impedance function that was derived and the pure hydraulic form bond graph. Converting the controller algorithm back into an impedance function ensures that any necessary modifications to the controller are accounted for. For each bond graph, the two-port transmission matrices must be found for 72 Fd vp vehicle body - mp b11 k11 Hydraulic piston Variable flow pump vt Figure 4.20: Simple hydraulic realization of equivalent active vehicle suspension with a hydraulic piston coupled with a variable flow pump. C : k11 R : b11 I : mp . vp p Sf : vt Ft vt 0 va Ft Ft vp Fd vp 1 Se : Fd . p21 Ft - Ap : T p21 1 s - b21 m21 k31 Sf : Qin x31 1 s . x31 1 m21 - k31 b31 Qin va Ap 1 m21 Figure 4.21: Bond graph of simple hydraulic realization with corresponding compensation control algorithm. 73 I : mp C : k11 Sf : vt R : b11 Ap : Sf : 0 0 T Se : Fd 1 0 Za(s) Z12 Z21 Z22 I : mp 0 1 Ap : Z11 Sf : vt T Se : Fd Figure 4.22: Original simple hydraulic system with impedance function represented as an element and the new pure hydraulic form. the area of modification. Since no changes will occur to the mass of the vehicle due to the change of physical domain, the inertia element is not needed for the two-port representation. Equation (4.17) shows the matrix function for the original case with a place holder transformer element added and equation (4.18) shows the matrix function for the pure hydraulic case. " A 0 # 0  1 A 1 0 1 1 k11 +b11 +Za (s) s   = A 0 1 A  k11 +b11 +Za (s) s  1 A   (4.17) In the case for the pure hydraulic case, the equivalent two-port impedance is unknown and thusly represented as a matrix of unknowns. " Z11 Z12 Z21 Z22 #" A 0 0 1 A # = " AZ11 AZ21 Z12 A Z22 A # (4.18) For each two-port representation to have equivalent impedance matrices, the derived two-port matrices must be identically equal. Solving for the unknown Z matrices, 74 the equivalent structure is found in equation (4.19). " Z11 Z12 Z21 Z22 #  = 1 0 1 A2  k11 +b11 +Za (s) s  1  (4.19)  The resulting impedance function from equation (4.19) has the structure of a zero junction where the admittance is made up of the inverse sum of the elements off the junction. However due to the transformation modulus, there is an additional scaling factor on the two passive elements and the negative impedance function. To compensate for this scaling modulus, each term in the controller must be scaled by 1/A. Additionally, the passive component parameters must be scaled by the same 1/A factor. With these modified parameters and the structure from the two-port secondary synthesis, the new bond graph structure can generated. Since the equivalent transformation preserves the characteristics of the impedance function, the function remains non-positive-real requiring an active realization. Figure 4.23 shows the new bond graph layout of the purely hydraulic realization of the active vehicle suspension system with the altered controller algorithm. This straightforward method procures an active suspension system realization which is entirely hydraulic. From the bond graph, the realization of the physical structure is shown in Figure 4.24. The compliance element is now realized as a hydraulic accumulator and the resistive element represents a return pressure loss. From the pure hydraulic bond graph in Figure 4.23, the state equations for the hydraulic active suspension system are derived for performance analysis. The state equation is given by equation (4.20). " ṗ V̇a # = " 0 A − mpp Ap C11 − R111C11 #" p Va # + " −1 0 0 0 Ap −1 #  Fd     vt  va (4.20) where k11 = 3.914 × 106 /A2p N-m, b11 = 854.4/A2p N/s, mp = 613/A2p kg, m21 = −15.2/A2p , b21 = −1596/A2p , k31 = −6.2 × 104 /A2p , and b31 = 1595/A2p . 75 I : mp . vp p Sf : vt Ft vt 0 va Ft Ft vp Fd vp 1 Se : Fd . p21 Ft - Ap : T k C : 11 Ap b R : 11 Ap Pa . Va Pa p21 1 s - b21 m21 Qout Pa 0 Ap k31 x31 . x31 1 s Ap m21 - k31 b31 Qin Pa Sf : Qin Qin Ap m21 Figure 4.23: Pure hydraulic actuation realization in bond graph form with modified compensation controller algorithm. 76 Fd vp vehicle body - mp Hydraulic accumulator vt Hydraulic piston Variable flow pump Figure 4.24: Pure hydraulic realization of the synthesized active vehicle suspension system 4.2.7 Simulation of Hydraulic Actuation System with Synthesized Compensator For verification of the validity of the results of the secondary synthesis procedure for the hydraulic active suspension system, simulations are performed on the system state equations. The velocities of the vehicular body due to a sinusoidal input disturbance from the tire are compared with to the desired response. Figure 4.25 depicts the frequency response comparison on the magnitudes with the desired system response which was specified. The frequency response for both the pure hydraulic case and the linear actuation realization are analytically identical and the frequency response shows the correspondence. Again for actuation sizing, the power requirements of the synthesized hydraulic actuator must be known. Figure 4.26 shows the pressure-flowrate curve of variable flowrate pump (flow source) for a constant sinusoidal input of 35 rad/s the pressure-flowrate curve for a random input. 77 Magnitude Gain Comparison -10 Pure Hydraulic Desired Gain (dB) -15 -20 -25 -30 -35 1 10 Frequency (rad/s) 10 2 Phase Comparison Phase (degrees) 20 0 -20 -40 -60 -80 -100 -120 1 10 Frequency (rad/s) 10 2 Figure 4.25: Frequency response comparison between pure hydraulic realization and the desired frequency response. Pressure-Flowrate of Actuator - Sinusoid Input Pressure-Flowrate of Actuator - Random Input 400 600 300 400 200 Pressure (kPa) Pressure (kPa) 200 0 100 0 -100 -200 -200 -400 -300 -600 -0.02 -0.015 -0.01 -0.005 0 3 0.005 0.01 0.015 0.02 -400 -0.01 -0.008 -0.006 -0.004 -0.002 0 0.002 0.004 0.006 0.008 0.01 3 Flow rate (m/s) Flow rate (m/s) Figure 4.26: Pressure-flowrate of pure hydraulic actuator at 35 rad/s input tire disturbance and random tire disturbance. 78 4.2.8 Summary of the Synthesis Methodology The resulting active suspension design variation, as derived from the analytical synthesis procedure, illustrates the diversity of the resulting topologies. As an early stage analytical methodology, the synthesis technique could potentially be used to supplement design insight of experts in a field. Additionally, the procedure might aid in guiding a novice designer to identify multiple system topologies dissimilar to conventional designs. Obviously, designers must exercise prudence and ensure physical realizability of a system as well as monitor the complexity of any derived topology. The following chapter extends the synthesis procedure to the control surface actuation retrofit for another approach at system design. 79 Chapter 5 Synthesis of Control Surface Actuator Systems The prior coverage and examples of the synthesis technique illustrate the capabilities for additional analytical design iterations, which provide insight into potential design topologies. Expanding the application of the synthesis method from hypothetical examples to the real world applications provides further evidence for the justifiability of the technique. Additionally, the continued exploration of the synthesis technique provides insights into the limitations of the technique, as well as its effectiveness. Building on the previous displacement control system design that was covered in Section 3.2, the synthesis technique is applied to the control surface actuation system retrofit for additional insight into energy storage. As was introduced in this thesis, the trend towards electrification of the naval fleet has encouraged retrofitting and new system insertion in naval vessels [1]. A more traditional simulation based design techniques was outlined in Section 3.4. However, this technique relies heavily on a posteriori design experience and knowledge of existing actuation performance that is not always readily available. Additionally many design assumptions proved necessary with the submarine model based technique. The energy storage is assumed to be ideal and coupled directly into the actuation system as a lossless system. The synthesis technique, as demonstrated in the previous section, provides not only analytical component magnitudes but also potential system layouts. 80 5.1 Conventional and Synthesis Design Comparison As was discussed before, a crucial unknown constraint for the control surface actuation system is the required torques for both tracking control for deflections and disturbance rejection. Two different methods for exploring potential characterizations of actuation systems are explored and outlined in Figure 5.1. The conventional/direct method is to derive a closed-loop feedback model containing an ideal actuator model and to run case studies on required torques [25]. Conversely using synthesis, specifications for disturbance rejection can be used to produce a compensator design to reject ocean disturbance. Note that for command following, a feedback loop must also be derived for tracking control of the control surface on top of the compensator design. Using the conventional method with a closed-loop control model enables an analysis and characterization of the actuation torques when subjected to control surface deflections, ocean disturbance moments and submarine pitch motion. However, the controller design must be insured to be sufficiently robust as to control motion during both standard positioning maneuvers with ocean disturbances as well as emergency/ extreme submarine maneuvers. This direct approach allows the characterization of the torque-speed, power and energy requirements for the control surface actuation system. Coupling the conventional method with synthesis provides a methodical path to investigate the effect of specific compensation for ocean disturbance torques. Essentially, the torques required to reject disturbances are isolated and the impedance structure can be analyzed for potential separable passive structures. Active actuation torque and energy with a compensator can potentially be reduced due to passive components and concentrated disturbance rejection. The claims on the ability to synthesize a disturbance rejection compensator are investigated in the following section. 5.2 Bond Graph Model and Matrix Derivation with Unknown Impedance The system model is constructed with a two-port model topology with one input as the submarine pitch velocity and the other as ocean disturbance torque, as seen 81 Synthesis Implementaon Convenonal Implementaon Idenfy system Model system, actuaon model as unknown impedance funcon Model system, actuaon model as ideal source Synthesis provides compensator for adjusng system frequency response Implement controller algorithm Size actuaon using control simulaon Closed-loop controller implemented for reference tracking/ transient response Realize system as ideal actuaon source Size actuaon using control simulaon Realize system as actuaon source and/or synthesized passive elements Figure 5.1: Flowchart comparison of the conventional design process in contrast with the synthesis procedure. 82 Control Surface Actuator Frame Submarine Pitch Motion - ωb Control Surface Angular Velocity - ωcs Drive Shaft Stiffness - K Z(s) B1 Base (Submarine) Linearized Hydrodynamic Compliance - Ko B2 Actuator Rotational Inertia - Jm Surface Rotational Inertia - Js Figure 5.2: Control surface system layout with added unknown impedance element. in Figure 5.2. This control surface system model assumes the actuation system will couple into rotational inertia, Jm , which is connected to the control surface by a drive shaft having stiffness, K. On each end of the drive shaft, bearing losses are accounted for by linear frictional damping coefficients B1 and B2 . The control surface is lumped as a single inertia with moment of inertia, Js . The unknown active impedance, Za (s), is inserted between the submarine pitch input and the control surface model. For straightforward realization, the form of the two-port unknown impedance function is restricted as a single impedance function off of a 0-junction, represented by equation (5.1). As discussed before, this restriction limits the quantity of designs but ensures a more simple realization procedure. Physically the form restriction provides a system topology that is guaranteed to be realizable in a desirable form. # # " #" " 1 0 T2 T1 = ω2 1/Za (s) 1 ω1 (5.1) Another important consideration in the model of the control surface is the hydrodynamic torques on the surface due to the motion of the fluid about the surface body. As was shown in Section 2.9.2 on the derivation of the submarine vehicular dynamics, the hydrodynamic moment equation for the control surfaces can be linearized and approximated as compliance as shown in equation (5.2). This linearized ocean “compliance” is also incorporated into the two-port model as another applied torque on the control surface inertia. To = ko θs ⇒ ko = ρAs Cs lp u2 83 (5.2) C : k/s 1 0 Ta 0 Sf : ωb R : b1 Z(s) R : b2 C : ko /s 1 I : Jm s ωcs I : Js s Se : Td Figure 5.3: Bond graph representation of control surface model with specified unknown impedance function. Representation of the control surface system model in bond graph topology is shown in Figure 5.3. Using the bond graph representation, the two-port matrices of each junction are formulated using transmission matrix form. Beginning with the right side of the bond graph model, the matrices are constructed for each element and parsed into series matrix multiplication as in equation (5.3). " T1 ω1 # = " 1 0 Ya (s) 1 #" 1 B1 + J m s 0 1 #" 1 0 s/K 1 # (5.3) " 1 B2 + Js s − Ko /s 0 1 #" Td ωcs # For this derivation, the drive shaft stiffness K is allowed to approach infinity. This assumption simplifies the control surface model by eliminating high frequency harmonics from the drive shaft compliance. From the series of element junctions, the impedance function that represents the entire two-port system is calculated. The total function is represented by equation 5.4. " T1 ω1 # = " − Kso + (Js + Jm )s + B1 + B2   Ya (s) Ya (s) (Js + Jm )s − Kso + B1 + B2 + 1 1 #" Td ωcs # (5.4) For the synthesis procedure, the transfer function of interest must be iden- 84 tified. The specified transfer function relates the effort input (ocean disturbance torque) to the flow output (control surface deflection) of the two-port system. Looking over the canonical matrix forms, the H (immittance) matrix form contains the relationship wcs /Ta (f2 /e2 ) given that f1 = q = 0. Using the canonical matrix relations discusses earlier, the H-form, of the system equations, is found as equation 5.5. " T1 ωcs # = " Ã B̃ C̃ D̃ #" ωcs Td # (5.5) where the matrix parameters are described by equation (5.6). Ã = B̃ = C̃ = D̃ = 5.3 (Js +Jm )s2 +(B1 +B2 )s−Ko Ya (s)[(Js +Jm )s2 +(B1 +B2 )s−Ko ]+s s Ya (s)[(Js +Jm )s2 +(B1 +B2 )s−Ko ]+s s Ya (s)[(Js +Jm )s2 +(B1 +B2 )s−Ko ]+s s Ya (s) Ya (s)[(Js +Jm )s2 +(B1 +B2 )s−Ko ]+s (5.6) Specify Desired Frequency Response Due to the infinite specifications that could potentially be imposed as a desired response, a method of procuring a feasible desired response is needed. Rather than specifying a physically impossible response, the natural response of the system with respect to the input of interest poses a starting point. As long as the specified desired response does not radically deviate from the natural response, the synthesis procedure returns a sensible result. For example, a system that contains a natural first order low pass response cannot be transformed into a fifth order purely differential system without the addition of complex active systems. In this transformation, the inserted impedance function would result in a improper transfer function. 5.3.1 Natural System Response to Torque Disturbances The best way to generate realizable systems in the case of the control surface disturbance rejection compensator is to specify the desired response as an alteration of the baseline system response (i.e. the system with no actuation). With the bond graph model of the control surface system, the transfer function relating the input disturbance torque (Td ) to the output angular velocity of the control surface (ωcs ) with no actuation can easily be derived. Expanding D̃ in the immitance (H) matrix, 85 Natural system (w / T ) frequency response cs -100 d Gain (dB) -110 -120 -130 -140 -150 -2 10 10 -1 10 0 10 1 10 2 10 3 Phase (degrees) 100 50 0 -50 -100 -2 10 10 -1 0 10 10 Frequency (rad/s) 1 10 2 10 3 Figure 5.4: Unactuated frequency response of control surface system. the disturbance transfer function is found to be equation (5.7) which assumes the input submarine pitch remains constant zero. ωcs = Td (Js + Jm ) s − Ko s 1 + B1 + B2 + Z(s) (5.7) With Z(s) set to zero (natural response of system), alterations to the disturbance transfer function response spur from this baseline response. Figure 5.4 shows the control surface system natural response to ocean torque disturbances. 5.3.2 Ocean Disturbance Torque Spectrum With the natural system response to torque disturbances, the ocean disturbance torque spectrum should be analyzed to determine what alterations on the natural response are necessary for disturbance rejection. Many different models for ocean wave spectrum exist, but for the purpose of this derivation, the Pierson-Moskowitz ocean spectrum is used [34, 45]. Typically the Pierson-Moskowitz ocean spectrum 86 Pierson-Moskowitz s-Domain Sea Spectrum 40 20 Gain (dB) 0 -20 -40 -60 -80 -2 10 10 -1 10 0 10 1 10 2 Frequency (rad/s) Figure 5.5: Magnitude of the s-Domain Representation of the Pierson-Moskowitz ocean spectrum. is given by equation (5.8). S(ω) = 0.0081g −3.109/h2s ω4 e ω5 (5.8) where ω is the frequency in radians per second, g is the acceleration is gravity, and hs is the significant wave height in meters. However, since all models are derived in the s-domain, a s-domain representation of the Pierson-Moskowitz spectrum is required. A fourth order s-domain transfer function has been derived in [45] and is given by equation (5.9). G(s) = Kn s2 s4 + a 1 s3 + a 2 s2 + a 3 s + a 4 (5.9) where Kn = 1.341, a1 = 0.641, a2 = 0.648, a3 = 0.135, and a4 = 0.0641. The s-domain Pierson-Moskowitz ocean spectrum magnitude is shown in Figure 5.5. The ocean disturbance spectrum contains discernable torques from 0.18 rad/s to 1.25 rad/s with a peak at 0.4 rad/s. Referring back to Figure 5.4, the magnitude of the natural response of the control surface system is near the peak in this band of frequencies. The control surface system is thus susceptible to ocean torque disturbances. To minimize the effect of the ocean disturbance torques, the natural response is altered to attenuate ocean disturbances in the discernible torque range. 87 Notch frequency response Gain (dB) 0 -50 -100 -2 10 10 -1 10 0 10 1 Phase (degrees) 100 50 0 -50 -100 -2 10 10 -1 10 0 10 1 Frequency (rad/s) Figure 5.6: Proposed notch response alteration to impose on the natural control surface system response 5.3.3 Notch Filter on Peak Ocean Torque Frequency One potential alteration of form that can be imposed on the natural system response is a notch filter form. Rather than attempting to attenuate a large band of frequencies, the peak frequency of the ocean disturbance torque (0.4 rad/s) can be selected for the notch location. The form of the notch response is shown in equation (5.10). f s2 + 2ζωc s + ωc 2 = e (s + ωc )2 (5.10) where ζ = 10−5 which is the notch depth, and ωc = 0.4 rad/s is the notch frequency location. This function has a response as shown in Figure 5.6. By imposing the notch function to the natural system response, the desired system response function is determined. Combining these functions yields equation (5.11). dT F (s) = s2 + 2ζωc s + ωc 2 1 · (Js + Jm ) s − Kso + B1 + B2 (s + ωc )2 (5.11) Merging these functions produces a desired transfer function that will not 88 Desired system ( wcs /T d) frequency response Gain (dB) -100 -150 -200 -250 -2 10 10 -1 10 0 10 1 10 2 10 3 Phase (degrees) 150 Natural Response Desired Response 100 50 0 -50 -100 -2 10 10 -1 0 10 10 Frequency (rad/s) 1 10 2 10 3 Figure 5.7: Comparison of the desired frequency response and the natural system response with a notch function. radically alter the natural response of the system, but modifies the system to attenuate disturbance frequencies from 0.1 rad/s to 1.2 rad/s. The desired frequency response which is used for the synthesis procedure is shown in Figure 5.7. While the notch filter modification to the desired system response isolates primarily a single frequency, the ocean spectrum perturbs the system on a broader band of frequencies. Modifying the natural system response by a notch response would be ideal for synthesis of a system for 60 Hz noise removal or an equivalent design criteria. However, in the case for disturbance rejection, the notch filter does not provide ample spectrum attenuation. 5.3.4 Band Reject on Ocean Torque Frequencies Rather than dTF set as only a notched response or a band reject response, the desired function corresponds to the natural response combined with a filtered response. This combination eliminates overcompensation at higher frequencies which consumes energy voraciously. To filter out ocean disturbance torques, the natural response of the system is modified to have a band reject response. This modifying 89 Band-reject frequency response 0 Gain (dB) -5 -10 -15 -20 -25 -2 10 10 -1 10 0 10 1 10 2 10 3 Phase (degrees) 100 50 0 -50 -100 -2 10 10 -1 0 10 10 Frequency (rad/s) 1 10 2 10 3 Figure 5.8: Band reject response alteration to impose on the natural control surface system response function is shown in equation (5.12). f = e " s ατ1 + 1 s τ1 + 1 · #n s τ2 /α + 1 s τ2 + 1 (5.12) where τ1 = 0.1 is the lower limit of the band reject function, τ2 = 5 is the upper limit of the band reject function, α = 6 is the slope factor, n = 2 is the filter order. This function has a response as shown in Figure 5.8. By imposing the band rejection function to the natural system response, the desired system response function is determined. Combining these functions yields equations (5.13). 1 · dT F (s) = (Js + Jm ) s − Kso + B1 + B2 " s ατ1 + 1 s τ1 + 1 · #n s τ2 /α + 1 s τ2 + 1 (5.13) Merging these functions produces a desired transfer function that will not radically alter the natural response of the system, but modifies the system to attenuate disturbance frequencies from 0.1 rad/s to 5 rad/s. The desired frequency 90 Desired system ( wcs /Td ) frequency response -100 Gain (dB) -110 -120 -130 -140 -150 -2 10 10 -1 10 0 10 1 10 2 10 3 Phase (degrees) 100 Natural Response Desired Response 50 0 -50 -100 -2 10 10 -1 0 10 10 Frequency (rad/s) 1 10 2 10 3 Figure 5.9: Comparison of the desired frequency response and the natural system response with a band reject function. response which is used for the synthesis procedure is shown in Figure 5.9. Since the notch filter only marginally attenuates frequencies other than the selected natural frequency, the band reject form of filtering is ideal for this design. Using the resulting desired transfer function, the synthesis methodology is executed to produce an actuation impedance which provides the desired response. In this case, the desired response only attenuates the natural response, so a positive real synthesis is expected. 5.4 Synthesized Impedance Function Now that both the two-port model of the system and the desired transfer function are specified and ensured to be in proper form, the actuation impedance is synthesized. Due to the two-port model form, we have four equations with four unknowns to solve. The elements H11 , H12 , and H21 are equivalent to unspecified unknown frequency responses and H22 is equivalent to the desired frequency response, as shown in equation (5.14). 91 Resultant system H11 125 Gain (dB) Gain (dB) -10 115 110 105 -20 -30 -40 100 -3 10 10 -2 -1 10 10 Frequency (rad/s) 0 10 1 10 -50 -3 10 2 Resultant system H21 0 10 -2 10 -1 0 10 Frequency (rad/s) 10 1 10 2 10 3 Specified System H22 -110 -120 Gain (dB) -10 Gain (dB) Resultant system H12 0 120 -20 -30 -40 -130 -140 -150 -50 -3 10 10 -2 10 -1 0 10 Frequency (rad/s) 10 1 10 2 10 -160 -3 10 3 10 -2 10 -1 0 10 Frequency (rad/s) 10 1 10 2 10 3 Figure 5.10: Comparison of the desired frequency response and the natural system response with a band reject function. H11 ⇒ Ã = U11 H12 ⇒ B̃ = U12 (5.14) H21 ⇒ C̃ = U21 H22 ⇒ D̃ = 1 (Js +Jm )s− Kso +B1 +B2 +1/Ya ·  1 +1 ατ1 1 +1 τ1 · 1 +1 τ2 /α 1 +1 τ2 n Solving the system of equations in equation (5.14) provides not only the resulting impedance function of the unknown impedance function Z(s), but also the resulting frequency responses of the other relations in the two-port system (U11 , U12 , U21 ). The resulting transfer functions for each element of the two-port matrix are shown in Figure 5.10. As seen in the graphs, the frequency response for H22 matches the specified response. Additionally, the resulting impedance function for Z(s) is found as in equation (5.15). This impedance function measures the difference in velocity between the submarine base and the control surface and produces a torque output to the system. Za (s) = a 4 s4 + a 3 s3 + a 2 s2 + a 1 s + a 0 b 4 s 4 + b3 s 3 + b2 s 2 + b1 s + b0 (5.15) where a4 = 2.63 × 1018 , a3 = 3.42 × 1019 , a2 = 9.50 × 1019 , a1 = 4.55 × 1019 , a0 = 5.00×1018 , b4 = 1.55×1013 , b3 = 4.43×1013 , b2 = 4.72×1013 , b1 = 2.22×1013 , and b0 = 3.87 × 1012 . The actuation impedance function has a frequency response 92 Frequency Response of Actuator Impedance Function Gain (dB) 130 120 110 100 -2 10 10 -1 0 1 0 1 10 10 Frequency (rad/s) 10 2 10 3 Phase (degrees) 50 0 -50 -100 -2 10 10 -1 10 10 Frequency (rad/s) 10 2 10 3 Figure 5.11: Frequency response of the synthesis derived actuation impedance function for the control surface system. as seen in Figure 5.11. The actuation frequency response peaks at 0.6 rad/s which means that the torques provided by the actuation function are the greatest at that point. A positive real test reveals that equation (5.15) is a positive real function, which implies a purely passive realization of the function is permissible. However, further investigation indicates that positive real pole removal is required at values other than s = 0 or s = ∞. As a result the fourth order synthesized impedance requires either complex coupled Brune elements or active components for realization [5, 40]. Using partial fraction expansion, the actuation impedance function is expanded into separate terms so that individual impedance function components can be identified. Partial fraction expansion of equation (5.15) yields equation (5.16) which contains active elements resulting from the form of the positive real function. Za (s) = Zactive (s) + Zpassive (s) Zactive (s) = 1 s2 + ks + b 1 k22 21 21 + Zpassive (s) = b11 + 93 1 s2 + ks + b 1 k32 31 31 1 s + b1 k12 12 + 1 s + b1 k41 41 (5.16) I : Ja s Sf : ωb 0 R : b1 C : k/s 1 0 R : b2 C : ko /s 1 I : Js s Passive Elements R : b12 Z2(s) 1 Se : Td 0 C : k12 /s Z3(s) R : b11 -C : k41 0 -R : b41 Active Elements Figure 5.12: Synthesized bond graph elements separated into both active and passive components. where passive parameters are b11 = 1.70 × 105 N/s, k12 = 9.16 × 107 N-m/rad, b12 = 1.10 × 108 N/s, and where active parameters are k21 = 4.82 × 106 N-m/rad, b21 = 1.61 × 107 N/s, k22 = 5.79 × 106 N-m/rad, k31 = 1.04 × 107 N-m/rad, b31 = 2.49 × 107 N/s, k32 = 1.73 × 107 N-m/rad, k41 = −8.98 × 107 N-m/rad, and b41 = −1.50 × 108 N/s. 5.5 Synthesized Bond Graph and Physical Realization The components of the expanded actuation function have been separated into active and passive components. The active components are either negative impedances or non-realizable complex impedance elements which cannot be represented by simple impedance components. The passive actuation components have potential to be realized as simple impedance components such as rotational springs and dampers. The full bond graph representation of the derived impedance function within the framework of the control surface system is shown in Figure 5.12. Assuming ideal linear actuator realization, two realizations of the actuation 94 ωb Ideal Rotational Actuator Control Surface Inertia δcs ωcs Driven by controller and combined compensator Figure 5.13: Control surface system realization with ideal rotational actuator. function are investigated. A single actuator realization of the impedance with both the passive and active impedance provides a straightforward insight into a simple system. The other potential structure formulation is to separate the passive components from the actuation impedance into individually realized structures. Due to the analytical natural of the derived actuation impedance function, there are ample potential physical realizations of the synthesized function. A simple realization is to assume that all impedance elements are ideal, linear and purely rotational systems akin to the control surface system. In this case, any negative/nonsimple element impedances are realized with an ideal rotational actuator. As was shown in equation (5.16), the derived impedance function separates into both active and passive components. Physically separating the realization of these functions remains up to the discretion of the system designer. The repercussions of component separation should be investigated as part of the overall synthesis design procedure. Both designs will be investigated in the following section. 5.5.1 Purely Active Realization of Synthesized System If the passive and active components of the synthesized impedance function are not separated, then the derived impedance function is ultimately represented by a single modulated effort source, as shown in Figure 5.13. The modulation of the torque source is governed by the derived actuation impedance function. 95 Flow-induced Torque Disturbances Td Reference Angle + δr,cs - + Control Surface Position Controller + Tact + Tnet Control Surface System + Tactive + Active + ωcs 1 s Control Surface Angular Position δcs Tpassive ω11 Passive + Compensation ωb Submarine Pitch Motion Figure 5.14: Block diagram representation of control surface system and compensator. For the torque source, the output effort (Ta ) is related to the flow signal (ω11 ) via the transfer function of the active elements. Using the state equations derived for the active elements in equation 5.17, the compensator controller function is derived. θ̇12 = ω11 − θ̇21 = ω11 − θ̇31 = ω11 − θ̇41 = ω11 − k12 b12 θ12 k21 b21 θ21 k31 b31 θ31 k41 b41 θ41 − − k21 k22 θ21 k31 k32 θ31 (5.17) where passive parameters are b11 = 1.70 × 105 N/s, k12 = 9.16 × 107 N-m/rad, b12 = 1.10 × 108 N/s, and where active parameters are k21 = 4.82 × 106 N-m/rad, b21 = 1.61 × 107 N/s, k22 = 5.79 × 106 N-m/rad, k31 = 1.04 × 107 N-m/rad, b31 = 2.49 × 107 N/s, k32 = 1.73 × 107 N-m/rad, k41 = −8.98 × 107 N-m/rad, and b41 = −1.50 × 108 N/s. Using the above differential equations, the specified torque output of the ideal actuator as related to the difference angular velocity, ω11 (difference between the submarine pitch rate and the control surface rate) results from the summation of the impedances given by equation (5.18). Ta = b11 ω11 + k12 θ12 + k21 θ21 + k31 θ31 + k41 θ41 (5.18) Representing equation (5.18) in block diagram form results in a controller algorithm that utilizes the parameters derived from the synthesis method. Figure 5.15 depicts the controller form integrated with the bond graph structure. 96 I : Ja s 0 Sf : ωb Se : Ta R : b1 C : k/s 1 0 R : b2 C : ko /s I : Js s 1 Se : Td ω11 T11 b11 . θ12 Ta θ12 1 s - T12 k12 + + + + + 1 b12 . θ21 - θ21 1 s - k21 T21 1 b21 . θ31 - - 1 s 1 k22 θ31 1 s T31 k31 1 b31 1 s . θ41 - 1 k32 1 s θ41 k41 T41 1 b41 Figure 5.15: Torque source controller algorithm as determined by the synthesis expansion incorporated in the bond graph structure. 97 From this form of bond graph, the state equations for the control surface system are derived for performance analysis. Additionally, note that as before the control rod compliance is assumed to be near infinite stiffness. Due to the value of k approaching infinity, the two inertia elements are lumped together as a single inertia (Ja + Js ). The state equation is given by equation (5.19). " ḣs θ̇s # = " +b2 − Jbs1+J a 1 Js +Ja ko 0 #" hs θs # + " 1 0 1 0 0 0 #  Td     ωb  Ta+p (5.19) where Js = 23204 kg-m2 , b1 = 1.665 × 103 N/s, b2 = 2.260 × 105 N/s, Ja = 50.16 kg-m2 , ko = −8.819 × 104 N-m/rad. For verification of the validity of the results of the synthesis procedure for the control surface system, simulations are performed on the system state equations including the synthesized actuation functions. Since the compensation function rejects ocean disturbance torques but has no effect on the transient position control, the first simulations hold the control surface steady at zero degrees deflection and reject disturbance. To compare the efficacy of the actuation function, the synthesis results are contrasted to an identical closed loop system sans the actuation function. Rather than isolating and separating the passive elements from the derived actuation impedance function, the derived impedance is lumped as a single ideal actuation element. Again, the ocean torque disturbance is modeled using the s-domain Pierson-Moskowitz spectrum in equation (5.9). Figure 5.16 shows the control surface deflection and angular velocity for closed loop actuation systems with and without synthesized compensation. The energy consumed by each actuation case is given below in Figure 5.18. The total energy consumed accounts for both energy required by the active compensator as well as the command following PID actuation. In the actuation only case, only the PID actuation is utilized. The energy consumption of the synthesized compensation case is reduced in comparison to the closed loop actuation only case. For characterization of the actuation requirements for disturbance rejection, the torque speed curve the previous simulation is shown in Figure 5.19. The torque-speed characteristic curve reveals that the synthesized actuation function shapes the response by selectively utilizing more torque to filter out disturbance torques. 98 Control Surface Displacement Displacement (deg) 4 Actuation Only With Compensator Reference 2 0 -2 -4 0 20 40 60 80 100 120 140 160 180 200 160 180 200 Time (s) Control Surface Angular Velocity Angular Velocity (deg/s) 2 1 0 -1 -2 -3 0 20 40 60 80 100 120 140 Time (s) Figure 5.16: Control surface displacement and angular velocity with single actuation impedance. Active Actuation Torque Actuator Torque (kN-m) 30 Actuation Only With Compensation 20 10 0 -10 -20 0 20 40 60 80 100 Time (s) 120 140 160 180 200 140 160 180 200 Actuator Power Flow (kW) 0.2 Power (kW) 0.1 0 -0.1 -0.2 -0.3 -0.4 0 20 40 60 80 100 Time (s) 120 Figure 5.17: Control surface actuator torque output and power flow velocity with single actuation impedance. 99 Energy Consumption Comparison 10 Actuation Only With Compensation Energy (kJ) 8 6 4 2 0 0 20 40 60 80 100 120 140 160 180 200 Time (s) Figure 5.18: Control surface actuator torque output and power flow velocity with single actuation impedance. Torque-Speed Characteristic Curve of Actuator 25 Actuation Only With Comp 20 Torque Output (kN-m) 15 10 5 0 -5 -10 -15 -20 -2 -1.5 -1 -0.5 0 0.5 1 Angular Velocity (deg/s) 1.5 2 2.5 Figure 5.19: Torque-Speed Characteristic Curve for single actuation impedance case compared to closed loop actuation. 100 Control surface Magnitude Gain -110 Ideal Active Compensation Desired Response Gain (dB) -120 -130 -140 -150 -2 10 10 -1 10 0 10 1 10 2 10 3 Frequency (rad/s) Control surface Phase 100 Phase (degrees) 50 0 -50 -100 -150 -2 10 10 -1 10 0 10 1 10 2 10 3 Frequency (rad/s) Figure 5.20: Frequency response of actuation system with single actuation impedance actuator. As expected, the synthesized actuation impedance increases the “stiffness” of the control surface system. Since the synthesized elements are lumped into a single active impedance, the increased stiffness is due to the reaction of the actuator alone. This means that the actuator does not have to compensate for stiff compliances at lower frequencies. In the control surface disturbance rejection case, the actuation function used jointly is ideal. Similarly for verification of the system frequency response, sinusoidal disturbance torques at specific frequencies were used to test the system response. The following plot, Figure 5.20, shows the magnitude and phase characteristics of the control surface system with separated passive elements to ocean disturbance input torques. 5.5.2 Separated Active and Passive Subsystem Realization As in the active suspension example, separation of the passive elements reduces the control algorithm by eliminating passive elements from the function and physically 101 ωb Ideal Rotational Actuator Passive compensation as rotational spring and damping δcs ωcs Active compensation and controller only Control Surface Inertia Figure 5.21: Control surface system realization with passive components separated. Flow-induced Torque Disturbances Td Reference Angle + δr,cs - Control Surface Position Controller + Tact + + Tnet + + + Tactive Active ωcs Control Surface System 1 s Control Surface Angular Position δcs Tpassive Passive ω11 + ωb Submarine Pitch Motion Figure 5.22: Block diagram representation of control surface system and compensator. realizing them as passive components. For the sake of comparison, the energy consumption should be compared between both the pure active case and the separated passive case. By separating the active and passive components, a modulated effort source and a passive spring damper system represent the derived impedance functions, as shown in Figure 5.21. The modulation of the torque source is governed by the derived actuation impedance function. For the modified torque source, the transfer function of the active elements without the passive components relates the output effort (Ta ) via the flow signal 102 (ω11 ). Using the state equations derived from the active elements in equation (5.20), the compensator controller function is derived. θ̇21 = ω11 − θ̇31 = ω11 − θ̇41 = ω11 − k21 b21 θ21 k31 b31 θ31 k41 b41 θ41 − − k21 k22 θ21 k31 k32 θ31 (5.20) where the parameters remain identical to those given above for equation (5.17). Again, the specified torque output of the ideal actuator as related to the difference angular velocity, ω11 is the summation of the impedances given by equation (5.21). Ta = b11 ω11 + k21 θ21 + k31 θ31 + k41 θ41 (5.21) Representing the modified controller function equation (5.21) in block diagram form results in a modified controller algorithm, this utilizes the parameters derived from the synthesis method. Figure 5.23 depicts the modified controller form integrated with the bond graph structure. From this form of bond graph, the state equations for the control surface system are derived for performance analysis. The state equation is given by equation (5.22). +b2 − Jbs1+J a ko 0    1  θ̇s  =  Js +J a 1 − Js +J θ̇12 a 0 0 0 12 − kb12  ḣs    hs   1 0 1  Td         θ s  +  0 0 0   ωb  θ12 Ta 0 1 0 (5.22) where the state space parameters are the same as equation (5.19). The following simulations represent the case where the passive components are isolated and separated from the active elements. The ocean torque disturbance is modeled using the s-domain Pierson-Moskowitz spectrum in equation (5.9). Figure 5.24 shows the control surface deflection and angular velocity for closed loop actuation systems with and without synthesized compensation. The energy consumed by each actuation case is given in Figure 5.26. As can be seen by the results, separation of the passive elements does not always yield a desirable result. Energy consumption in the case where the elements have been separated is an order of magnitude greater than the closed loop only actuation. For 103 I : Ja s 0 Sf : ωb R : b1 C : k/s 1 0 R : b2 C : ko /s I : Js s 1 R : b12 1 T12 Se : Td 0 . θ12 Se : Ta C : k12 /s ω11 T11 + + + b11 . θ21 Ta - θ21 1 s - k21 + T21 1 b21 . θ31 - - 1 s 1 k22 θ31 1 s T31 k31 1 b31 1 s . θ41 - 1 k32 1 s θ41 k41 T41 1 b41 Figure 5.23: Torque source controller algorithm with passive components separated. 104 Control Surface Displacement 4 Actuation Only With Compensator Reference Displacement (deg) 3 2 1 0 -1 -2 -3 -4 0 20 40 60 80 100 120 140 160 180 200 160 180 200 Time (s) Control Surface Angular Velocity Angular Velocity (deg/s) 2 1 0 -1 -2 -3 0 20 40 60 80 100 120 140 Time (s) Figure 5.24: Control surface displacement and angular velocity with separated passive elements. 105 Active Actuation Torque 1000 Actuator Torque (kN-m) Actuation Only With Compensation 500 0 -500 -1000 0 20 40 60 80 100 Time (s) 120 140 160 180 200 140 160 180 200 Actuator Power Flow (kW) 6 5 Power (kW) 4 3 2 1 0 -1 0 20 40 60 80 100 Time (s) 120 Figure 5.25: Control surface actuator torque output and power flow with separated passive elements. Energy Consumption Comparison 140 Actuation Only With Compensation 120 Energy (kJ) 100 80 60 40 20 0 -20 0 20 40 60 80 100 120 140 160 180 200 Time (s) Figure 5.26: Energy consumption by actuation system disturbance rejection with separated passive elements. 106 Torque-Speed Characteristic Curve of Actuator 800 Actuation Only With Comp 600 Torque Output (kN-m) 400 200 0 -200 -400 -600 -800 -2 -1.5 -1 -0.5 0 0.5 Angular Velocity (deg/s) 1 1.5 2 2.5 Figure 5.27: Torque-Speed Characteristic Curve for separated passive elements case compared to closed loop actuation. characterization of the actuation requirements for disturbance rejection, the torque speed curve for ocean torque reject is shown in Figure 5.27. Again, the required torques are more than an order of magnitude greater in the separated elements case. The simulations illustrate that the separation of passive elements from the actuation function in not always ideal. Since the separated spring and damper combination have large coefficients, the passive elements radically increase the “stiffness” of the entire control surface system and create a parasitic effect. The active actuation function must compensate for the increased system stiffness by over-actuating to provide the desired system response. For verification of the system frequency response, sinusoidal disturbance torques at specific frequencies were used to test the system response. The following plot, Figure 5.28, shows the magnitude and phase characteristics of the control surface system with separated passive elements to ocean disturbance input torques. The simulated frequency response closely matches that of the desired response that was specified for the synthesis procedure. The compensation function helps to shape the system response over the range of specified frequencies. 107 Control surface Magnitude Gain -110 Ideal Active Compensation Desired Response Gain (dB) -120 -130 -140 -150 -2 10 10 -1 10 0 10 1 10 2 10 3 Frequency (rad/s) Control surface Phase 100 Phase (degrees) 50 0 -50 -100 -150 -2 10 10 -1 10 0 10 1 10 2 10 3 Frequency (rad/s) Figure 5.28: Frequency response of actuation system with separated passive elements. 5.6 Model Reduction of Synthesized Impedance The fourth order complexity of the synthesized impedance function complicates the realization procedure. Although the complex impedance function yields positive real characteristics, a purely passive realization of the function could not be completely derived due to function complexity. A purely passive realization of the fourth order system with complex roots requires coupled passive elements as shown in Temes [40]. Coupled passive elements are physically difficult to realize in many physical domains and should be avoided while using synthesis techniques. Another technique for realization of actuation impedances produces a reduced order model of the synthesized impedance function. Rather than attempting to realize the impedance function in the entirety, a balanced realization of the impedance model provides insight into the possibility of order reduction. Additionally, the state of positive-realness of a reduced order impedance function remains the same as the complex function. A reduced order impedance function contains many benefits. Not only is the positive-realness property of the function preserved, 108 but also the number of states is potentially reduced. The reduced order system also might eliminate the need for required coupled passive elements and provide a pure passive realization. The order of a given system can always be reduced via balanced realization with increasing error depending on the magnitude of the state contributions. The interested reader can find additional coverage of the balanced realization procedure in Appendix B. Conversion of the original resulting impedance function, shown again in equation (5.23), to a state space representation is necessary to produce the balanced realization. Za (s) = a 4 s4 + a 3 s3 + a 2 s2 + a 1 s + a 0 b 4 s 4 + b3 s 3 + b2 s 2 + b1 s + b0 (5.23) where a4 = 2.63 × 1018 , a3 = 3.42 × 1019 , a2 = 9.50 × 1019 , a1 = 4.55 × 1019 , a0 = 5.00×1018 , b4 = 1.55×1013 , b3 = 4.43×1013 , b2 = 4.72×1013 , b1 = 2.22×1013 , and b0 = 3.87×1012 . Conversion from transfer function form produces Za (s) in state space form as shown by equation (5.24).     Za (s) =     −2.867 −3.055 −1.433 −0.25 1 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1.725e6 5.628e6 2.698e6 2.873e5 1.702e5         (5.24) Computing the balanced realization of the impedance space state function produces both the balanced state space realization as well as the error G vector. The error G vector for Za(s) is given below, Gerror = h 1.5235 0.8543 0.0972 0.0026 iT The order of magnitude drop between the second and third entries in the error G vector implies that an accurate realization of a second order system exists for Za (s). The relative contributions from the balance third and fourth states are an order of magnitude less than the contributions of the first and second states. The reduced order Za (s) system is procured by elimination of the low contribution states shown 109 Frequency Response Comparison - Original Model and Reduced Order 130 Original System Reduced Order Magnitude (dB) 125 120 115 110 105 Phase (deg) 100 45 0 -45 -90 10 -2 10 -1 10 0 10 1 10 2 10 3 Frequency (rad/sec) Figure 5.29: Comparison of original fourth order system to the reduced order second order approximation. in equation (5.25).  −0.7043 −0.6297  Zred (s) =  0.6297 −1465 −0.0484 −287.7 −1465   287.7  1.702e5 (5.25) Conversion of the second order state space state representation to a SISO transfer function produces the reduced impedance function. Zred (s) = 1.702 × 105 s2 + 2.192 × 106 s + 6.498 × 105 s2 + 0.7527s + 0.4306 (5.26) The reduced actuation function has a frequency response as shown in Figure 5.29. For comparison, the frequency response plot also contains the original synthesized impedance function. Since the original impedance function is a positive real function, the reduced order system also is positive real such that a purely passive realization of the function must exist. Due to the order reduction, the realization of the impedance 110 R : b1 C : k/s 0 1 0 1 R : b11 R : b12 0 C : k12 /s I : m21 s 1 R : b21 I : Ja s Sf : ωb Passive Elements R : b2 C : ko /s 1 I : Js s Se : Td Figure 5.30: Reduced order impedance function in bond graph form yielding a purely passive system. function with desired response yields a less complicated result which can more readily be physically constructed. Continuous fraction expansion of the Zred (s) function produces equation (5.27). Zred (s) = b11 + s k12 + 1 b12 1 + 1 m21 s+b21 (5.27) where b11 = 1.70 × 105 , k12 = 2.06 × 106 , b12 = 4.36 × 106 , m21 = 6.917 × 106 , and b21 = 1.932 × 106 . Since all element coefficients from the resulting function expansion are positive and real values, the total function is purely passive. The full bond graph representation of the reduced impedance function within the framework of the control surface system is shown in Figure 5.30. One prospective realization of this passive system produces a purely mechanical linear system as shown in Figure 5.31 (a). The block m21 can be physically realized via an inerter element as proposed by Smith [46]. A more direct realization, due to the control surface system type, presents 111 b11 b12 k12 b11 b21 b12 m21 b21 k12 m21 (b) (a) Figure 5.31: (a) Linear mechanical realization of reduced impedance function. (b) One potential rotational realization of the reduced impedance. the passive elements in the rotational mechanical domain as seen in Figure 5.31 (b). The elements both in the linear and rotational domains must preserve the quantities resulting from the synthesis or the dynamic performance suffers. Since these quantities are sensitive to variability, a conversion to the electrical domain could prove more feasible for physical realization. Another potential passive realization of the control surface system shifts the linear mechanical components into the electrical domain. The physical domain transformation can be accomplished by moving the impedance function behind the rotational actuator source via a gyrator transformation. As in Section 4.2.6, the two-port impedance of the modified bond graph must be equivalent to the original synthesized bond graph for identical performance. To find the equivalent impedance of the new modified bond graph via secondary synthesis, the original bond graph transmission matrix must equate to the modified transmission matrix. Figure 5.33 shows the original bond graph with the flow source reverted to the original negative impedance function that was derived and the pure hydraulic form bond graph. Converting the compensator algorithm back into an impedance function ensures that any necessary modifications to the controller are accounted for. For each bond graph, the two-port transmission matrices must be found for modified portions of the bond graph. Since no changes occur beyond the synthesized 112 Flow-induced Torque Disturbances Td Reference Angle + δr,cs Control Surface Position Controller - + Tact + Tnet + ωcs Control Surface System + 1 s Control Surface Angular Position δcs Tpassive ω11 Passive + Submarine Pitch Motion ωb Figure 5.32: Block diagram representation of control surface system and passive compensator. Sf : ωb Kt : Se : 0 0 G Ta Z(s) R : b1 C : k/s 1 0 Z21 Z22 Kt : Z12 G 0 C : ko /s 1 I : Jm s Sf : ωb Z11 R : b2 I : Js s ωcs Se : Td R : b1 C : k/s 1 0 Ta I : Jm s R : b2 C : ko /s 1 ωcs I : Js s Se : Td Figure 5.33: Original pure mechanical domain system and the new partially electrical domain system with unknown equivalent two port impedance function. 113 impedance function, the remaining portions of the bond graph are not needed for the two-port representation. Equation (5.28) illustrates the matrix function for the original case with a placeholder gyrator element added, and equation (5.29) displays the matrix function for the electrical conversion. # " #" " 1 0 0 KT = 1 1 0 Zred (s) 1 K KT Zred (s) 1 KT T KT 0 # (5.28) In the case for the electrical conversion, the equivalent two-port impedance is unknown and thusly represented as a matrix of unknowns. " Z11 Z12 Z21 Z22 #" 0 KT 1 KT 0 # = " Z12 KT Z22 KT KT Z11 KT Z21 # (5.29) For each two-port representation to have equivalent impedance matrices, the derived two-port matrices must be identically equal. Solving for the unknown Z matrix, the equivalent structure is found in equation (5.30). " Z11 Z12 Z21 Z22 # = " 1 KT2 Zred (s) 0 1 # (5.30) The found impedance function from the above equation has the structure of a 1junction where the impedance is composed of the inverse of the elements off the junction. However due to the gyrator modulus, there is an additional scaling factor on all of the impedance elements. To compensate for the gyrator modulus, the coefficients of the impedance elements must be scaled accordingly. The resulting equivalent impedance function takes the form of equation (5.31). Zeq (s) = 1 b11 KT2 + (5.31) 1 1 m12 KT2 s+b12 KT2 + s k21 K 2 +s b21 K 2 T T / / where b11 = 1.70 × 105 , m12 = 2.06 × 106 , b12 = 4.36 × 106 , k21 = 6.917 × 106 , b21 = 1.932×106 and Kt is any positive real number. With these modified parameters and the structure from the two-port secondary synthesis, the new bond graph structure can generated. Figure 5.34 shows the new bond graph layout of the electrical conversion real114 I : Ja s Sf : ωb 0 : Kt G 1 2 R : KT b12 2 C : k21 KT R : b1 C : k/s 1 0 R : b2 C : ko /s 1 Sf : iact I : Js s Se : Td 2 0 R : b11/KT 1 I : KT m12 0 R : b21 KT 2 2 Figure 5.34: Passive compensation as realized in electrical form inside electromechanical actuator. ization of the control surface ocean disturbance rejection compensator. Realization of the gyrator and electrical domain portion of the bond graph yields Figure 5.35. The gyrator represents the ideal torque source actuation. Moving the compensation function from the mechanical domain to the electrical domain produces an internal voltage filter which indirectly filters the torque supplied to the control surface system. For verification of the validity of the results of the order reduction procedure for the control surface system, simulations are performed on the system state equations which include the synthesized reduced passive functions. To compare the efficacy of the reduced passive function, the synthesis results are contrasted to an 115 2 2 KT m12 KT b12 2 Tact KT b21 2 KT k21 2 Gyrator - KT b11/KT Figure 5.35: Electrical schematic of electrical compensation function. identical closed loop system sans the actuation function. Figure 5.36 shows the control surface deflection and angular velocity for closed loop actuation systems with the reduced order passive impedance. The energy consumed by each actuation case is given below in Figure 5.38. As can be seen by the results, a purely passive realization considerably reduces the energy consumption of the actuator for disturbance rejection. Energy consumption in the reduced order impedance function case is an order of magnitude less than any of the other actuation cases. For characterization of the actuation requirements for disturbance rejection, the torque speed curve for ocean torque reject is shown in Figure 5.39. Again, the required torques are more than an order of magnitude less in the pure passive elements reduced case. Similarly for verification of the system frequency response, sinusoidal disturbance torques at specific frequencies were used to test the system response. The following plot, Figure 5.40, shows the magnitude and phase characteristics of the control surface system with separated passive elements to ocean disturbance input torques. Since the deviation of the plot due to model reduction is minimal, the effects of the reduced impedance function are minute. 5.7 Synthesis Realization Comparison Due to the counteractive effect of separating the passive elements from the control algorithm, any permutations of the separated design will also consume excess energy. Attempting a hydraulic or electrical realization with separated passive elements 116 Control Surface Displacement Displacement (deg) 4 Actuation Only With Compensator Reference 2 0 -2 -4 0 20 40 60 80 100 120 140 160 180 200 160 180 200 Time (s) Control Surface Angular Velocity Angular Velocity (deg/s) 2 1 0 -1 -2 -3 0 20 40 60 80 100 120 140 Time (s) Figure 5.36: Control surface displacement and angular velocity with either reduced mechanical or electrical system. Actuator Torque (kN-m) Active Actuation Torque 20 Actuation Only With Compensation 10 0 -10 -20 0 50 100 Time (s) Actuator Power Flow (kW) 150 200 0 50 100 Time (s) 150 200 Power (kW) 0.4 0.2 0 -0.2 -0.4 Figure 5.37: Control surface actuator torque output and power flow with either reduced mechanical or electrical system. 117 Energy Consumption Comparison 10 Actuation Only With Compensation Energy (kJ) 8 6 4 2 0 0 20 40 60 80 100 120 140 160 180 200 Time (s) Figure 5.38: Energy consumption by actuation system disturbance rejection with either reduced mechanical or electrical system. Torque-Speed Characteristic Curve of Actuator 20 Actuation Only With Comp 15 Torque Output (kN-m) 10 5 0 -5 -10 -15 -2 -1.5 -1 -0.5 0 0.5 1 Angular Velocity (deg/s) 1.5 2 2.5 Figure 5.39: Torque-Speed Characteristic Curve for either reduced mechanical or electrical system compared to closed loop actuation. 118 Control surface Magnitude Gain -110 Ideal Active Compensation Desired Response Gain (dB) -120 -130 -140 -150 -2 10 10 -1 10 0 10 1 10 2 10 3 Frequency (rad/s) Control surface Phase 100 Phase (degrees) 50 0 -50 -100 -150 -2 10 10 -1 10 0 10 1 10 2 10 3 Frequency (rad/s) Figure 5.40: Frequency response of actuation system with either reduced mechanical or electrical system. yields identical results as the results of the separated compensation simulations above. As a result, the results of this synthesis procedure in this case should be restricted to the combined active and passive compensation controller or a reduced order purely passive compensator. Due to the resulting restriction for the synthesis results, no extrinsic energy storage is required for disturbance rejection in this case. However, the underlying compensation algorithm and pure passive case suggests intrinsic energy storage in the form of the hidden capacitive components. These underlying capacitive elements help to shape the response using the disturbance to aid in rejection. Since no energy storage is required for disturbance rejection, any added energy storage components supply energy to aid with transient control surface maneuvers and/or aid in the case of power failure. 119 Chapter 6 Conclusions 6.1 Summary and Contributions As part of the continued investigation for the feasibility of electrification of the naval fleet, this work began the investigation of retrofitting legacy hydraulic actuation systems in submarine vessels. In addition to required actuation moments, the potential need for energy storage was pursued utilizing two complementary techniques. Utilization of a detailed submarine model in conjunction with models of the control surface system formed the basis of the conventional approach. Simulation of detailed submarine maneuvers elucidated potential power requirements for actuation maneuvers. The synthesis design approach explored potential design realizations via specification of frequency response requirements on the control surface system. The two methodologies were also contrasted and the limitations and scope of each were detailed. As a preliminary for submarine control surface modeling, the non-linear time invariant equations of motion for a submerged vessel were derived. Additionally, the hydrodynamic moments encountered on the skeg control surfaces were discussed. While not necessary in the synthesis methodology, the submarine vehicular dynamics are critical for establishing dynamic submarine states using the conventional approach. With the submarine dynamics established, several extreme maneuver simulations provided insight into the vessel behavior as well as the torque demands on the actuator. Several critical realizations resulted from the simulations. In an 120 emergency surfacing maneuver, the actuation restoring the sternplane to maintain constant submarine pitch requires 78% more energy than the deflection maneuver. Furthermore, failure to perform this second actuation results in a near vertical vessel ascent. The full sweep maneuver proved to require the most energy per actuation and thusly should be utilized for energy storage size estimation. The two-port impedance bond graph synthesis methodology and required background knowledge were also introduced. Additionally, due to the developing nature of the active bond graph synthesis methods, several modifications to the methodology were introduced and discussed. Bond graph models for the control surface system were determined, allowing two-port impedance matrices to be derived. An actuation function was derived for ocean disturbance torque rejection, and functional expansion revealed various potential system topologies. Ultimately, order reduction was performed on the synthesized impedance function revealing a potential purely passive system realization. The primary original contributions to bond graph synthesis techniques of this research are briefly summarized as follows. 1. A technique for secondary synthesis was developed and demonstrated. Secondary synthesis provided a method for extracting equivalent two-port transmission matrices to produce equivalent structures across transformers or gyrators. This procedure proves important for conversion of systems across physical domains. (e.g. Linear mechanical system behind transformer to electrical system behind a gyrator.) 2. This research demonstrated the viability of order reduction on derived impedance functions for more physically realizable systems. As a result of model order and desired response order, a synthesized impedance function yields a higher order system. The complexity of the resulting system makes physical realization difficult and oftentimes causes positive real systems to result in active realizations. Order reduction of the impedance function produces simplified realizations and lower order passive real systems are often easily realizable. 121 6.2 Future Work Both the conventional method and the synthesis techniques presented in this thesis provide insight into energy demands of actuator systems under differing conditions. Of course, the developments suggested by this thesis move the field of synthesis and modeling forward, and the understanding of these techniques continues to expand. While both modeling and synthesis offer knowledge of system performance in the frequency domain, only the conventional modeling process generates direct insight for the transient response of an actuator subsystem. Bond graph impedance synthesis techniques are currently limited to frequency domain design. 6.2.1 Extension of Synthesis Techniques Extension of the bond graph synthesis methodology to incorporate transient requirements as well as frequency domain requirements would provide even more flexibility for the procedure. Designers could not only specify specific frequencies for attenuation, but also settling time and overshoot. For first or second order systems, simple conversions between transient and frequency domain requirements exist such as rise time, percent overshoot and settling time. However, for higher order systems simple relations are not nearly as lucid. Rather than modeling the system as an n-port impedance model, the state space bond graph model could be derived. Bond graphs in state space form provide the ability of system inversion [47]. An inverted system would allow the designer to specify a transient system output and the required inputs would result. This system inversion technique would allow desired transient responses rather than impedance functions as requirements. 6.2.2 Statistical Analysis of Power Variables In addition to extensions on the synthesis methodology, further work is needed in the characterization of energy storage demands. The current techniques outlined in this thesis reveal potential energy demands but offer no further insight on more specific types of energy storage. As was discussed in Section 3.5, power and energy draw requirements determine strong candidates for energy storage device type. A statistical technique for characterization of power draw and energy consumption 122 over a range of different simulation scenario could further spur energy storage device selection. Data collected over a range of potential simulation scenarios for both effort and flow from an actuator system can be collected and used to determine both transient power and energy consumed over time. Inspection of the power draw over many scenarios could potentially reveal separate zones of power consumption. A significant portion of the power demands makes up the standard power draw from the actuation system, peaks in the power consumption might identify zones of abnormally high power draw. If such a separated zone exists, power from the grid demanded by the actuator could be satisfied by lower standard consumption found while the remaining rare power peaks might be satisfied by high power density energy storage such as ultra capacitors or flywheels. 6.3 Final Remarks Both the conventional simulation based technique and the bond graph synthesis technique developed in this work were shown to be viable approaches for the early stage design of actuator subsystems. An approach for deriving candidate designs for actuator retrofit based on a conventional model basis and two-port synthesis methods have been described in this thesis. Currently, these techniques complement each other as preliminarily design techniques. The conventional detailed model basis reveals energy requirements at the expense of requiring an expert knowledge basis for accurate results. On the other hand, the synthesis technique requires only a basic system model requiring limited designer knowledge, making the method well suited for early design. However, synthesis is not without limitations as well. The technique only currently allows compensator synthesis rather than entire control system synthesis. The work of this thesis has determined how the synthesis methods benefits early stage design. Additionally, the model basis demonstrated how ship simulations with required maneuvers can directly generate the actuation system requirements. In addition, by comparing the different synthesis derived actuation models that augment the feedback-controlled torque, this work demonstrated that synthesis can lead to a notable reduction in energy consumption for the same torque requirement. This result suggests that synthesis methods may be helpful in guiding how passive and/or 123 active topologies for energy storage might be used within an actuation subsystem. Ultimately, this research has shown that utilization of the synthesis methods described can be used more universally to explore retrofitting existing technology. If a given actuation system exists, for example, deriving a model basis provides a direct means for synthesizing a generalized impedance function which can then be used to study realizations using other types of technologies. For example, dynamic performance specifications for a hydraulic system may be used to generate an equivalent electromechanical form. 124 Appendix A Realizability Conditions for Passive Systems Conditions exist which sufficiently indicate that a given impedance function is realizable via strictly passive elements [5]. However, note that satisfying these conditions does not imply that the passive realization of the system is succinctly derived with straightforward techniques. Complex manipulations and couplings of impedance functions are sometimes required to generate a purely passive realization, especially in the case of high order impedance functions. The following definitions provide the outlines of the composition of passively realizable systems in terms of all bond graph elements. A.1 One-port Passive System Conditions The location of the poles and/or zeros in a given polynomial function is an important concept when dealing with the complex s-domain. Many insights into the function response and type derived from the location of poles and zeros in the complex sdomain. Definition 1 A given polynomial function, A(s), in the s-domain is said to be a Hurwitz polynomial if, A(s) = a0 + a1 s + . . . + an sn 1. A(s) is a real polynomial, meaning that all values ai are real coefficients. 125 2. All of the zeros of A(s) are in the left half-plane. Any zeros that lie on the imaginary axis must be simple zeros. (Re(s) ≤ 0) The polynomial A(s) is said to be a strictly Hurwitz polynomial if all the zeros are in the left half plane (LHP) without any zeros which lie on the j-axis. (Re(s) < 0) Definition 2 A one-port impedance function, Z(s), is said to be passively realizable if the given impedance function found to be a positive real function [5]. 1. The impedance function, Z(s), must be a real rational function of s. Z(s) = a 0 + a 1 s + . . . + a n sn b0 + b1 s + . . . + bm s m where the coefficients ai and bi are all real values. 2. For all real values of ω, the following condition is satisfied. ReZ(jω) ≥ 0 3. All the poles of Z(s) are in the closed LHP of the s-plane. The third condition of positive realness on the pole locations can also be expressed in terms of Hurwitz polynomials. Depending on the polynomial type, the third condition is satisfied if the following conditions are met. 1. The denominator of Z(s) must be either a Hurwitz or strictly Hurwitz polynomial. 2. If the denominator polynomial is only Hurwitz, then the zeros on the jω-axis must be simple. Additionally, the poles as s → ∞ must also be simple. A.2 Methods for determination of positive-realness Simple tests have been found which provide insight into the positive real nature of a given function [40]. The first condition of positive realness is easily determined via inspection and verification of each polynomial coefficient being real valued. The second and third conditions are more difficult to ascertain and cannot be immediately found via inspection. 126 A technique which is a sufficient condition for proving Re(jω) ≥ 0 (condition 2 for positive real functions), is Strum’s theorem. Strum’s theorem provides a method through which the numbers of sign changes, meaning a simple real zero value) occur on the positive axis. Strum’s theorem presents the following condition. P (x) = an xn + an−1 xn−1 + . . . + a1 x + a0 ≥ 0 where x ≥ 0 and the function P(x) is found by division of the function of interest into even and odd polynomial values for both the numerator and denominator. The terms are then combined together as follows to produce P(x) where ω 2 = x. P (ω 2 ) = Neven (jω)Deven (jω) − Nodd (jω)Dodd (jω) The condition of Strum’s theorem requires that P(x) contains no simple real zeros on the positive interval, which satisfies Re(jω) ≥ 0. Sufficient conditions for the satisfaction of the third condition of positive realness can also be found via the manipulation of numerator and denominator of the function of interest. The following modification of the function Z(s) must have all zeros in the open LHP. Z mod (s) = N (s) + D(s) N (s) +1= D(s) D(s) The modified function above implies that N(s) + D(s) must be strictly Hurwitz in order to have all zeros in the open LHP. Therefore, the addition of the numerator and denominator provides a quick test of third condition of positive realness of a given function. A.3 Extension to Two-port Systems The prior discussion only dealt with the possibility of passive realization for a given one-port impedance transfer function. The following definition extends the previous definitions to encapsulate two-port functions [40]. As in the case of the one-port system, a two-port system must be a positive real function for a purely passive realization to exist. 127 Definition 3 A two-port impedance function system given by a two-port matrix Z(s) is positive real if the following conditions are satisfied. 1. Z(s) exists and is analytic on the interval of σ > 0, where σ is the real portion of s (s = σ + jω). 2. Z ∗ (s) = Z(s∗ ) on the interval σ > 0, where the superscript asterisk indicates the complex conjugate. 3. ZH (s) ≥ 0 on the interval σ > 0, where ZH (s) designates the Hermitian portion of the matrix. The above conditions for positive realness of two-port systems can also be expressed in a perhaps more readily understood form. The following conditions are equivalent to the initial three of the definition. 1. All elements of Z(s) are real rational functions. 2. Z(s) contains no poles in the range of σ > 0. 3. The poles of Z(s) on the jω-axis are simple values. 4. For the poles on the jω-axis, the residues of each element must be simple values. 5. ZH (jω) ≥ 0 for all real values of ω 128 Appendix B Balanced Realizations for Order Reduction In the case of several of the derived impedance functions found throughout this thesis, the order of the resulting function was the magnitude of both the order of the model function and the desired response function. For example, a second order model with and unknown impedance function that is run through the synthesis procedure with a desired response of a second order system procures a fourth order impedance function. Higher order impedance functions tend to be physically difficult to realize and even more difficult to physically manufacture. Additionally, high order functions that are found to be positive real often cannot easily be decomposed into purely passive elements without significant function manipulation. One potential method for impedance model simplification is order reduction of the function [48]. Reliable reduction of model order can be achieved by internally balancing the model system. This balanced realization contains the original number of states, but readjusts the state space values so that the relative contributions from each state are in descending order. B.1 Conditions for a Balanced Realization Consider a large order impedance transfer function Z(s) which must be either proper or strictly proper. By the definition of the transfer function, the function Z(s) can 129 be converted from transfer function form to the state space form. Z(s) = ẋ = Ax + Bu N (s) ⇒ D(s) y = Cx + Du A state coordinate transformation (x̂ = T x) can be performed on the state space representation of the function Z(s). This similarity transformation results in an equivalent system as shown. ˆ˙x = T AT −1 x̂ + T Bu y = CT −1 x̂ + Du In the case of a balanced realization, the unique transformation T when applied to both the observability gramians (Lo ) and controllability gramians (Lc ) results in modified gramians that are equivalent. The modified gramians are given as follows. L̂o = T −T Lo T −1 L̂c = T Lc T T In the case of an internally balanced realization when the modified gramians are equivalent, the resulting structure of the modified gramians is in diagonalized form. L̂o = L̂c = diag(g1 , g2 , . . . , gn ) where the diagonal terms are in descending order (g1 ≥ g2 ≥ ≥ gn ). Reduction of system order eliminates states from the balanced realization system according to the diagonalized gramians form. A proper model reduction identifies an order of magnitude drop in the diagonalized matrix as an indication of decreasing state contributions. The states of the balanced realization at this division are eliminated which reduces the model order. A reduction to any order is possible, however reduction errors become ever more present with increasing order reduction. kG(s) − Gapprox (s)k∞ = 2 n X λk k=m+1 The errors of the model reduction are bounded by the sum of the eliminated states from the diagonalized matrix gramian forms. 130 B.2 Method to Determine Balanced Realizations To determine the unique transformation T that yields the balanced realization, an algorithmic process that produces a decomposition of one of the gramians is found to produce equivalent controllability and observability gramians [48]. The following procedure for procuring the unique transformation is implementable in the Matlab environment with the “balreal” command. 1. The observability gramian must be factored into two equivalent matrices with one being orthogonal to the other. Lo = R T R where the matrix R cannot be singular. 2. Combining the R matrix found from the observability gramian with the controllability gramian produces an altered controllability matrix. Decomposing the new matrix into Eigenvalue and eigenvector matrices yields the following. RLc RT = U Λ2 U T where U are the eigenvectors which combine to form an orthogonal unitary matrix and Λ is a diagonal matrix comprised of the Eigenvalues. 3. The following procedure yields the unique transformation that produces the balanced realization. The matrix T is defined as follows. T = R−1 U Λ1/2 131 Bibliography [1] “ONR BAA Announcement Number 09-011,” Tech. Rep. Number 09-011, 2008. [2] H. Stevens, “Electric warship science and technology: A summary of onr 334 planning information,” Tech. Rep. MSD-50-TR-2002/23, 2002. [3] “Hydraulic v. electrical power for deepwell cargo pump systems,” tech. rep., Deltamarine for Hamworthy KSE. [4] B. Stafford and N. Osborne, “Technology development for steering and stabilizers,” Journal of Engineering for Maritime Environment, 2008. [5] O. Brune, “Synthesis of finite two-terminal network whose driving point impedance is a prescribed function of frequency,” Journal of Mathematics and Physics, vol. 10, pp. 191–236, 1931. [6] W. Cauer, Synthesis of Linear Communication Networks. New York: McGrawHill Book Company, 1958. [7] R. Foster, “Theorems on the driving point impedance of two-mesh circuits,” Bell Systems Technology Journal, vol. 3, p. 651, 1924. [8] R. W. Newcomb, Linear Multiport Synthesis. McGraw-Hill Book Company, Inc., New York, 1966. [9] A. Larky, “Negative-impedance converters,” Circuit Theory, IRE Transactions on, vol. 4, no. 3, pp. 124–131, 1957. [10] M. Bayard, “Synthesis of n-terminal pair network,” Proceedings of Brooklyn Polytech. Symposium on Modern Network Synthesis, vol. 1, pp. 66–83, 1952. 132 [11] N. Balabanian, Network Synthesis. New Jersey: Prentice-Hall, Inc., 1958. [12] D. Hazony, Elements of Network Synthesis. New York: Reinhold Publishing Corp., 1963. [13] L. P. Huelsman, Circuits, Matrices, and Linear Vector Spaces. McGraw-Hill Electronic Sciences Series, New York: McGraw-Hill Book Company, Inc., 1963. [14] H. M. Paynter, Analysis and Design of Engineering Systems. The MIT Press, 1960. [15] D. C. Karnopp, D. L. Margolis, and R. C. Rosenberg, System Dynamics: Modeling and Simulation of Mechatronic Systems. John Wiley & Sons, Inc., New York, 2001. [16] R. Redfield and B. Mooring, “Concept generation in dynamic systems using bond graphs,” ASME Winter Annual Meeting symposium, vol. 8, 1988. [17] R. C. Redfield, “Bond graphs as a tool in mechanical system conceptual design,” Automated Modeling, vol. 41, 1992. [18] T. J. Connolly, Synthesis of Multiple-Energy Active Elements for Mechanical Systems. PhD thesis, 2000. [19] T. J. Connolly and R. G. Longoria, “An approach for actuation specification and synthesis of dynamic systems,” Journal of Dynamic Systems, Measurement, and Control, vol. 131, no. 3, 2009. [20] S. Kim, A Bond Graph Approach to Analysis, Synthesis and Design of Dynamic Systems. PhD thesis, 2003. [21] M. Gertler and G. R. Hagen, “Standard equations of motion for submarine simulation,” tech. rep., Naval Ship Research and Development Center, June 1967 1967. [22] J. Feldman, “Dtnsrdc revised standard submarine equations of motion,” tech. rep., Naval Ship Research and Development Center, 1979. 133 [23] G. D. Watt, “Modelling and simulating unsteady six degrees-of-freedom submarine rising maneuvers,” tech. rep., Defence Research and Development Canada - Atlantic, February 2007 2007. [24] G. D. Watt, “Estimates for the added mass of a multi-component deeply submerged vehicle,” tech. rep., Defence Research and Development Canada - Atlantic, October 1988. [25] T. I. Fossen, Guidance and Control of Ocean Vehicles. Chichester: John Wiley & Sons., 1994. [26] A. F. Molland and S. R. Turnock, Marine Rudders and Control Surfaces. Oxford: Elsevier Ltd., 2007. [27] S. Hoerner, Fluid Dynamic Drag. Hoerner Fluid Dynamics, 1993. [28] M. S. Triantafyllou and F. S. Hover, “Maneuver and control of marine vehicles,” 2003. [29] S. Liu, L. Fang, and J.-l. Li, “Rudder/flap joint intelligent control for ship course system,” in IEEE International Conference on Mechatronics and Automation, (Harbin, China), pp. 1890 – 1895, IEEE, 2007. [30] I. H. Abbot, Theory of wing sections, including a summary of airfoil data. New York: Dover Publications, 1959. [31] USS Pampanito - Hydraulic System. 2009. [32] R. Burcher and L. Rydill, Concepts in Submarine Design. Ocean Technology, Cambridge: Cambridge University Press, 1994. [33] D. P. Rubertus, L. D. Hunter, and G. J. Cecere, “Electromechanical actuation technology for all-electric aircraft,” IEEE Transactions on Aerospace and Electronic Systems, vol. 20, no. 3, p. 243, 1983. [34] W. J. Pierson and L. Moskowitz, “A proposed spectral form for fully developed wind seas based on the similarity theory,” tech. rep., New York University, 1963. [35] J. Jensne, Fundamentals of energy storage. New York: Wiley, 1984. 134 [36] J. Wyatt, L. Chua, J. Gannett, I. Goknar, and D. Green, “Energy concepts in the state-space theory of nonlinear n-ports: Part i-passivity,” Circuits and Systems, IEEE Transactions on, vol. 28, no. 1, pp. 48–61, 1981. [37] J. Wyatt, L. Chua, J. Gannett, I. Goknar, and D. Green, “Energy concepts in the state-space theory of nonlinear n-ports: Part ii - losslessness,” Circuits and Systems, IEEE Transactions on, vol. 29, no. 7, pp. 417–430, 1982. [38] D. Karnopp and R. C. Rosenburg, Analysis and Simulation of Multiport Systems. Cambridge, Massachusetts: The MIT Press, 1968. [39] R. C. Redfield, “Configurational design of an internal velocity sensor using bond graphs,” Advances in Instrumentation, vol. 30, 1991. [40] G. C. Temes and J. W. LaPatra, Introduction to Circuit Synthesis and Design. New York: McGraw-Hill Book Company, 1977. [41] R. Bott and R. Duffin, “Impedance synthesis without use of transformers,” Journal of Applied Physics, vol. 20, p. 816, 1949. [42] D. Karnopp, “Power requirements for vehicle suspension systems,” Vehicle System Dynamics, vol. 21, no. 1, pp. 65–71, 1992. [43] D. Karnopp, “Active and semi-active vibration isolation,” Journal of Mechanical Design, Transactions of the ASME, vol. 117B, pp. 177 – 185, 1995. Vibration isolation;Semi active suspensions;Signal processors;. [44] M. C. Smith and G. W. Walker, “Performance limitations and constraints for active and passive suspensions: a mechanical multi-port approach,” Vehicle System Dynamics, vol. 33, no. 3, pp. 137–168, 2000. [45] C. K. Li and J. R. Leigh, “S-domain realisation of sea spectrum,” Electronics Letters, vol. 19, no. 5, 1983. [46] M. C.Smith, “Synthesis of mechanical networks: The inerter,” 2002. [47] R. F. Ngwompo, S. Scavarda, and D. Thomasset, “Inversion of linear timeinvariant siso systems modelled by bond graph,” Journal of the Franklin Institute, vol. 333, no. 2, pp. 157 – 174, 1996. 135 [48] B. Moore, “Principal component analysis in linear systems: Controllability, observability, and model reduction,” IEEE Transactions on Automatic Control, vol. 26, pp. 17–32, 1981. 136 Vita Jonathan Robert LeSage was born in Temple, Texas on December 4, 1985, the son of Gene and Jannette LeSage. Jonathan completed his undergraduate education at The University of Texas at Austin in 2008 with a Bachelors of Science in Mechanical Engineering. During his undergraduate studies, he participated as an undergraduate researcher working on projects on both design collaboration techniques and hybrid vehicle modeling and simulation. Additionally, Jonathan worked as a mechanical research intern for Schlumburger in both Houston, Texas and Clamart, France in 2007 and 2008 respectively. Returning to academics after his internship experience, Jonathan enrolled in the Dynamic Systems and Controls area of the Department of Mechanical Engineering at The University of Texas at Austin in the fall of 2008. Upon entering the program, he was hired as a graduate research assistant at the Applied Research Laboratories under the supervision of Bill Shutt and the academic guidance of Dr. Raul Longoria. Permanent Address: 207 Highland Ridge Terrace Johnson City, TN 37615 This thesis was typeset with LATEX 2ε 1 by the author. 1 A L TEX 2ε is an extension of LATEX. LATEX is a collection of macros for TEX. TEX is a trademark of the American Mathematical Society. The macros used in formatting this thesis were written by Dinesh Das, Department of Computer Sciences, The University of Texas at Austin, and extended by Bert Kay, James A. Bednar, and Ayman El-Khashab. 137