OSCILLATIONS OF A BEAM IN FLUID FLOW: A MEANS FOR UNDERWATER PROPULSION By Sanders Wainwright Aspelund A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering – Doctor of Philosophy 2023 ABSTRACT Aquatic animals commonly oscillate their fins, tails, or other structures to propel and control themselves in water. These elements are not perfectly rigid, so the interplay between their stiffness and the fluid loading dictates their dynamics. We examine the propulsive qualities of a tail-like flexible beam over a range of flow speeds, with oscillations induced either by a known forcing function, feedback-based actuation, or a follower force. This is accomplished using the equations of fluid-immersed beams in combination with a set of tractable expressions for thrust and efficiency. We first compute these equations for a flexible beam actuated by a known forcing function over the external flow velocity and forcing frequency plane and show that the flexible propulsor has regions of both positive and negative thrust. We show the behavior of a sample underwater vehicle with fixed drag characteristics as an illustration of a realizable system. An alternate approach to generating the tail motion for underwater propulsion is to use a method where the oscillation of the flexible element is self-induced. We investigate a pinned-free beam in axial fluid flow, subjected to feedback-based actuation at the pinned end. The actuation may be a moment or a prescribed angle, and it is proportional to the state (curvature, slope, or displacement) of the beam at some point along its length. All equations and boundary condition terms are non-dimensionalized and the stability of the system is studied over a range of external flow velocities and sensing locations. For each combination of flow velocity and sensing location, the critical gain (positive or negative) for the onset of flutter is determined. This process, which is repeated for each combination of actuation and sensing modes, reveals that the closed-loop system exhibits a rich set of stability transitions, each associated with a traveling waveform in the flexible beam at the onset of flutter. With the intent of exploring the use of flexible fluttering beams for underwater propulsion, the efficiency of these waveforms is computed using slender-body theory. Additional insights into the efficiency of the waveforms are obtained through considerations of the smoothness of the traveling waveforms. Using a water tunnel at various flow speeds, we provide experimental validation of feedback-induced flutter of a beam for the specific case of moment actuation proportional to curvature. We show that the model, adapted for experimental considerations, results in flutter at very similar frequencies as the experimental results. We complete our investigation of flutter-based propulsion by examining how the onset of flutter is affected by the application of a follower force to a cantilevered beam in fluid flow. We consider the full range of fluid-mass to beam-mass ratios (from thick pipes to thin beams) given the possibility that the follower force could be generated by internal flow and a fluid jet. We follow the stability of the system as the flow velocity is increased to determine the onset of flutter. Over this full range of mass ratios and external flow velocities, we determine the critical follower force and critical frequency. This set of investigations, along with our investigations of forced oscillations and feedback-induced flutter, provide insights into how a flexible propulsor may be used as an alternative to the propeller. This thesis is dedicated to my parents who never gave up on me and provided me the strength to go on. Thank you. iv ACKNOWLEDGEMENTS I would like to begin by thanking my advisor, Dr. Ranjan Mukherjee for his unending patience and guidance through this process. Your continuous drive to make progress on this work gave me the strength to make it to the end and I would not be publishing this document if not for your support. I would also like to thank my Navy mentor, Dr. Aren Hellum for his tremendous advice and encouragement. Your enthusiasm and knowledge made it always a pleasure working with you. I also thank my remaining committee members, Dr. André Bénard and Dr. Firas Khasawneh for their support in getting me over the finish line. My utmost gratitude goes to my parents, Curtis and Hélène, for their inexhaustible support and encouragement. They got me through my darkest times and back on my feet and moving towards the final goal. Each step of the way, they were always there to celebrate my achievements and give me fuel to keep pushing. I would also like to thank my siblings, Keegan, Morgan, and Stanton, and also Arjun, for always being there to lend an ear and share tales of our lives with each other. To my wonderful girlfriend Rachel, I give my thanks for being the beacon on the horizon and giving me a reason to keep moving forward. And to Liam Faulkner, I thank for having given me the means with which to wrangle my mind and be able to make the progress needed to make it through. I cannot go without thanking Connor, Aakash, and Spencer, the best friends someone could ask for on this journey, and Krissy, who walked the path before me and provided a safe harbor for me on the way, both in sharing her home and always being there to talk to and empathize. My favorite memories of graduate school nearly always involve one of you! To my labmates, Anshul, Mahmoud, Sheryl, Robert, Shaede, Nilay, and Amer, and my undergraduates, Eddie, Nik, Jason, and Ben, I am grateful for all the moments we had, learning about each of you and helping me develop as a researcher, person, and mentor. A final shout-out goes to all of my fellow swing dancers and co-opers, for providing a space of camaraderie and fun. The support provided by the Office of Naval Research, ONR Grant N00014-19-1-2535, and by the National Science Foundation, NSF Grant 1703735, is gratefully acknowledged. v TABLE OF CONTENTS CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Forced Excitation of a Flexible Beam in Fluid Flow . . . . . . . . . . . . . . . . . 2 1.2 Feedback-Induced Flutter Instability of a Flexible Beam in Fluid Flow . . . . . . . 3 1.3 Effects of a Follower Force on the Inducement of Flutter in a Flexible Beam in Fluid Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Organization and Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 CHAPTER 2 FORCED EXCITATION OF A FLEXIBLE BEAM IN FLUID FLOW . . . 12 2.1 System Description and Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Dynamic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Method of Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Investigation of Propulsive Performance . . . . . . . . . . . . . . . . . . . . . . . 20 2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 CHAPTER 3 FEEDBACK-INDUCED FLUTTER INSTABILITY OF A FLEXIBLE BEAM IN FLUID FLOW . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 Method of Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3 Investigation of Flutter Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.4 Results of Feedback-Induced Instability . . . . . . . . . . . . . . . . . . . . . . . 37 3.5 Application to Underwater Propulsion . . . . . . . . . . . . . . . . . . . . . . . . 46 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 CHAPTER 4 EXPERIMENTAL VALIDATION OF FEEDBACK-INDUCED FLUTTER INSTABILITY OF A FLEXIBLE BEAM IN FLUID FLOW . . . . . . . . . 56 4.1 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.2 Experimental Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.4 Numerical simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.5 Conclusion: Comparison of Numerical and Experimental Results . . . . . . . . . . 86 CHAPTER 5 EFFECTS OF A FOLLOWER FORCE ON THE INDUCEMENT OF FLUTTER IN A FLEXIBLE BEAM IN FLUID FLOW . . . . . . . . . . . 88 5.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.2 Method of Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.3 Investigation of Flutter Instability . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 CHAPTER 6 CONCLUSION AND FUTURE DIRECTIONS . . . . . . . . . . . . . . . . 121 6.1 Forced Excitation of a Flexible Beam in Fluid Flow . . . . . . . . . . . . . . . . . 121 6.2 Feedback-Induced Flutter Instability of a Flexible Beam in Fluid Flow . . . . . . . 121 6.3 Experimental Validation of Feedback-Induced Flutter Instability of a Flexible Beam in Fluid Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 vi 6.4 Effects of a Follower Force on the Inducement of Flutter in a Flexible Beam in Fluid Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6.5 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 vii CHAPTER 1 INTRODUCTION Underwater robotic vehicles, also known as unmanned underwater vehicles (UUVs), are a rapidly growing field given their increasing utility due to autonomy, capabilities, and cost. These vehicles are used in a variety of applications, including military operations, environmental monitoring, exploration, and research [1]. In recent years, advances in robotics, propulsion, and materials science have enabled the development of highly capable UUVs with unique bio-inspired designs and new capabilities [2, 3]. UUVs can operate at extreme depths, in conditions that would be hazardous or impossible for humans to explore such as oil, gas, and mineral exploration [4–7]. UUVs are also being used to explore and map new and challenging environments such as underwater caves, which require highly maneuverable vehicles that can navigate through narrow passages [8]. Additionally, the military is one of the largest users of UUVs, with applications ranging from mine detection and disposal to submarine tracking and surveillance [9, 10]. UUVs can also be used for reconnaissance and intelligence gathering, and for delivering payloads to targeted locations [11]. Finally, UUVs are increasingly being used for environmental monitoring, including tracking pollution levels, monitoring fish populations, and studying the impacts of climate change on marine ecosystems [12–14]. One of the key challenges for underwater vehicles is propulsion. Water is denser and more viscous than air, making it difficult for traditional propulsion systems to achieve high speeds or maneuverability. To overcome this challenge, researchers have developed bio-inspired propulsion systems that mimic the movements of marine animals [15]. For example, some UUVs use flapping fins or undulating membranes to move through the water which offers improved maneuverability over traditional propeller-based systems [16, 17]. Efficiency is another important factor for UUVs, as they typically rely on batteries for power and may need to operate autonomously for extended periods of time, allowing for continuous data collection and analysis [18]. In this thesis, we explore how oscillations of a beam in fluid flow, whether due to a known sinusoidal input or as a result of flutter, can be harnessed to serve as a means of propelling an underwater vehicle. 1 Our research is different from standard propellers in that it explores the use of oscillating a flexible beam through forcing or flutter. Since we are using flexible beams, it behooves us to study the elastic response and stability of structures. There has been considerable work in the field of vibrations which focuses on the frequency dependent response of elastic structures to increase efficacy, efficiency, and in many cases, prevent damage from resonance [19, 20]. The study of stability focuses on the behavior of elastic structures when subject to differing system or environmental parameters, a field of study since as early as the late 1800s e.g. [21]. Stability can be lost drastically and permanently, such as the buckling of a column subject to too great of a compressive force, or through oscillation, such as with a garden hose whipping back and forth when unsecured [22]. Most research into oscillatory response behaviors has been for the purpose of minimizing or preventing them [23, 24], but here we seek take advantage of their amplifying power to generate thrust. 1.1 Forced Excitation of a Flexible Beam in Fluid Flow There is tremendous interest in bio-inspired underwater propulsion systems as alternatives to the standard propeller found on most undersea vehicles [25]. Nature has evolved a number of periodic propulsive techniques [26], e.g. thunniform, carangiform, anguilliform, batoid, that enable a large range of creatures to traverse the ocean [27, 28]. The periodic motions of fish allow them to attain tremendous speeds and agility during hunting and evasive maneuvers while also exhibiting impressive efficiencies over long distances. Given this impressive performance, a significant amount of analytical work has been committed to understanding the hydrodynamics of periodic motion, with the intention of adapting it to engineered systems [29–32]. Experimental work has also been done with foils and panels which are periodically pitching [33], heaving [34], or both pitching and heaving [35–37]; and using different shapes [38, 39], aspect ratios [40], and flexibilities [41, 42]. Because these engineered systems can be run repeatably and understood in detail, the analysis, measurement, and control [43] of these systems have provided scientists and engineers with a better understanding of periodic propulsion [28]. With a better understanding of periodic propulsion, a number of robotic implementations have 2 been built. These include the well-known RoboTuna [44], the recent designs described in [45–47], and many others reviewed in [48]. Many of these designs use fully-actuated, rigid-link systems for propulsion, but such fully-articulated systems have a drivetrain complexity that increases with the number of joints. For the purpose of simplification, flexible propulsors have been used in place of fully-actuated systems. Researchers have studied the effects of flexibility using distinct foils and panels [49,50], as well as continuous tails actuated with a dynamic moment [51] or with tail motion generated through instability [52]. Systems using a flexible tail for propulsion, as opposed to a fully-actuated rigid tail, are governed by equations modeling fluid-immersed beams [53–55]. In this work, we investigate the propulsion due to a flexible tail actuated by a dynamic moment at its base. This model can be studied over a large range of frequencies and external flow velocities, treating them as independent variables, to determine the tail motion. The tail motion of the forced system can be used to evaluate propulsive performance using Lighthill’s equations for thrust and Froude efficiency for slender bodies [29]. Together, these can be used to determine the resultant velocity, thrust, and efficiency of a robotic submersible driven at a specific tail frequency, based on its drag model and other system parameters. 1.2 Feedback-Induced Flutter Instability of a Flexible Beam in Fluid Flow An alternate method of producing oscillatory motion in a fluid-immersed beam for propulsion is through inducing flutter, rather than through the input of a known sinusoidal motion. The phenomenon known as flutter is a parametric self-excitation that occurs in many types of systems. Classic examples range from the casual, such as the flapping of flags in the wind [56], to the concerning, such as unwanted vibrations in the trans-Arabian oil pipeline [57], to the catastrophic, such as bridges collapsing (e.g. the infamous Tacoma Narrows Bridge [58]) and the loss of several aircraft [59, 60]. These examples, particularly the latter ones, have led to rigorous study of the phenomenon of flutter since at least the early 1950s [61–63] to improve these systems with respect to safety and efficiency. Work continues on these topics through today with more sophisticated methods and models. In general, elastic systems can lose stability through either flutter, characterized by oscillations of 3 increasing amplitude until a limit cycle is reached (or failure occurs), or divergence, characterized by near zero-frequency (often static) deflections away from the original point of equilibrium [64]. The mode of instability is a function of the type of loading (conservative vs non-conservative) [65, 66] as well as whether the system is damped or undamped [67]. Systems may also be studied from a discrete perspective [62, 68], with a limited number of natural frequencies and associated mode shapes, or from a continuous perspective [69], with infinite natural frequencies and associated mode shapes. For undamped systems (continuous or discrete), the natural frequencies lie on the axis of marginal stability and move towards or away from the origin as the system parameters change. If the complex conjugate pair of roots associated with the lowest frequency meet at the origin, then stability is lost through divergence [67], e.g. the buckling of a column subject to an axial compressive load. If two adjacent natural frequencies approach each other and meet, then the system loses stability through coupled-mode flutter. On the contrary, for undamped systems, the natural frequencies are not relegated to the axis of marginal stability and therefore stability can be lost through single-mode flutter in which a single complex conjugate pair crosses the axis of marginal stability onto the unstable half of the plane. Divergence can result from either conservative or non-conservative loading while flutter typ- ically1 results from non-conservative loading [72]. Elishakoff [73] presents a review of system stability subject to follower forces which are often used to generate non-conservative loading. This greater body of literature on the phenomenon of flutter of elastic systems includes discrete systems, but we focus on continuous systems for this work (which are surveyed by Langthjem and Sugiyam [74]), particularly with an emphasis on beams and pipes subjected to axial flow. These structures are also non-conservative systems and much of the early work surrounding their study is credited to Paidoussis [75] and Benjamin [68]. Investigation into this subject has been continued by several [76–80] to better understand the stability characteristics of more complex systems under varying conditions. 1In some recent works, it has been shown that flutter instability can occur in the presence of conservative loading when subject to non-holonomic constraints [70, 71]. 4 As alluded to in the examples given earlier in this section, most studies attempt to understand this behavior to avoid the instability. On the contrary, Hellum [52,81] investigated how to induce flutter instability for the purpose of underwater propulsion. Work by Abdullatif [82] further investigated the instigation of flutter through applying a dynamic moment to a beam with and without external flow. In this work, we extend this work by considering a significantly broader scope of methods to induce flutter through feedback. We also conduct experiments to validate our analytical and numerical findings. 1.3 Effects of a Follower Force on the Inducement of Flutter in a Flexible Beam in Fluid Flow As mentioned, another method of inducing flutter (or divergence) in beams and pipes is through the application of follower forces. Given the simplicity of analysis for the standard forms, but the impressive insights and lessons that can be learned from studying them, they have continued to be a very popular subject. Surveys of the field [83] were already being published in the early 1970s though the survey conducted by Langthjem and Sugiyama [74], provides a tremendous review of the field up through 2000. This survey focused on simple elastic structures subject to non-conservative follower loads, with a particular emphasis on canonical problems, including Beck’s column and Reut’s column, among others. Beck’s column is a simple, thin, flexible cantilevered beam that has a compressive force applied to its free end, with the force always aligned tangentially to the free end of the beam [84]. This differs from Reut’s column, for which the compressive force always acts in the direction of the central axis of the undeformed beam [61]. These seminal works have inspired several generations to continue studying the effects of non-conservative follower loads and how they can induced instability. Efforts ranged from studying the effects of a compressive follower force for asymmetric beams, showing how changing parameters could lead to flutter or buckling [85] to studies on how the frequency response of MEMS resonator could be tuned by varying the follower force. More recent work studies the non-linear limit cycle behavior of the 3D beam showing complex behavior [86]. The field has not been without criticism as Elishakoff [73] writes about, given the difficulty of 5 producing a tangential follower force without introducing significant additional effects such as fluid flows or changing tip masses. An early attempt by Done applied a follower force with a jet engine to show the onset of flutter instability, however, the study only lasted one second [87]. Sugiyama later attempted to show the same with a rocket engine, which burned for four seconds [88]. More recently Bigoni et al presented experimental evidence of flutter and divergence instabilities induced by dry friction, first with rigid links, and then with a continuous beam i.e. a Pflüger column (a Beck’s column with a concentrated tip mass) [89, 90]. Alongside this work on the stability of columns subjected to follower forces, both tangential and axial, there has been significant work on the stability of fluid conveying pipes and fluid-immersed beams, a field largely led by the work of Paidoussis and covered in his books on the topic [67, 91]. These systems follow similar equations of motion as elastic beams with follower forces as the effluent jet provides a compressive effect on the beam, an effect which can be amplified by the use of a convergent nozzle. However, the fluid imparts additional damping on the system which is why such lengths were gone to when trying to experimentally validate Beck, Reut, and Pflüger columns using jet and rocket engines [87, 88]. Different loading schemes were studied, such as pinned-pinned or clamped-pinned fluid-conveying pipes under axial load [92] or pinned-free fluid conveying pipe with an axial load but with the pinned end constrained by a rotational spring of varying stiffness [93]. Clamped-free fluid-conveying pipes with distributed tangential loads (akin to a very viscous, slow moving fluid being forced through a pipe) were studied [94] with Karimi- Nobandegani adding the condition of having the tube spin [95] which add additional inertial and coriolis forces. Ilgamov studied the effects on stability of the internal pressure from fitting a nozzle to the end of a vertical fluid-conveying pipe that also had a follower force seeking to validate the jet flow assumptions that had been made in previous works on the effects of the flow exiting the jet [96]. More recently, several studies have involved research into the behavior of nanorods and nanotubes subject to fluid flow conditions [97] with Bahaadini et al incorporating the nonlocal surface effects [98] and Zhou et al studying the dynamics of spinning nanotubes conveying magnetic nanofluid [99]. 6 Various methods have also been utilized to analyze the dynamics of fluid-conveying pipes with axial loads. Laithier studied the non-linear Timoshenko equations of a cantilevered and clamped- clamped fluid-conveying pipe through several derivations and compared them to the equations of motion derived from Newtonian methods, showing strong agreement [100]. Guran developed a fluid-conveying Reut’s column that accounts for shear stresses and extensibility, showing a loss of stability through both divergence and flutter [101]. Several studies were conducted as to the non-linear oscillations of fluid-conveying pipes, [102, 103] including vibration suppression in the post-critical regime [104]. A very thorough and elucidating survey of non-linear techniques can be found in [91]. 1.4 Organization and Contributions In this section we present the organization of this dissertation with an overview of each chapter and the research contributions presented therein. In Chapter 2, we delve into the propulsion capabilities of a flexible, tail-like beam that is driven by a dynamic moment over a broad frequency and flow speed range. We model an underwater vehicle that consists of a rigid body and a flexible beam connected by a revolute joint. It is here that we introduce the general structure of the flexible, tail-like beam in axial fluid flow that will be the focus of this thesis. The center of mass of rigid body of the model is limited to moving axially along a channel by a pinned joint about which it is free to rotate. The leading edge of the beam is attached by a revolute joint to the trailing edge of the rigid body while the trailing edge of the beam is free. To simplify the model and make it easier to solve, we establish a set of assumptions about the fluid-structure interactions and the physical properties of the rigid body and flexible tail-like beam. We then use the equation of motion for fluid-immersed beam oscillations and the dynamics of the rigid body acting as boundary conditions to mathematically define the system. Our method of solution involves finding the beam’s motion through sinusoidal angle actuation of the revolute joint. We present a series of tractable equations that allow us to calculate the thrust and efficiency of the system as it moves across the velocity-frequency plane. Lastly, we provide a demonstration of a sample underwater vehicle with fixed drag properties to give a practical example of how this 7 system might work in real-life situations. The contributions of this work presented in Chapter 2 are as follows. We determine the analytical model for a rigid body in flow and show that it is numerically tractable to solve over a large domain of model and state parameters. We show the utility of using slender body equations for thrust and power to analyze a large region of the frequency - flow velocity space. We show that the flexible propulsor exhibits areas of both positive and negative thrust. We discuss the effects of resonance on the system. By defining a drag model of the rigid body, we determine the equilibrium velocity of an example submersible propelled by a tail-like beam. We show that the efficiency exceeds 50% on this locus of dynamic equilibrium. In Chapter 3, we examine the use of different forms of feedback to induce flutter behavior in a beam in fluid flow. Our focus is on a model of a pinned-free beam that is actuated at its leading edge by a moment or angle proportional to the beam’s displacement, slope, or curvature at a specific point along its length. We describe the equation of motion and boundary conditions of the pinned joint which depend on the actuation and sensing method used. To simplify our analysis, we non-dimensionalize the equations of motion and boundary conditions. Next, we present a method of solution to determine the natural frequencies of the beam, taking into account parameters such as beam properties, fluid flow velocity, feedback gain, and the location of sensing. Initially, we examine the behavior of the system as the external flow velocity increases, with no feedback gain, for both the moment actuation and angle actuation cases. We then show how to determine the feedback gain at which the system becomes unstable. From the twelve possible combinations of actuation, sensing, and feedback gain, we select six cases to demonstrate the range of behaviors possible. We identify patterns in the critical stability and frequency surfaces and compare the behavior of the different cases. We then evaluate the potential of using flutter oscillations as a propulsor for an underwater vehicle, examining the thrust and efficiency characteristics and how they are influenced by the parameters. We also examine how the smoothness of the waveform affects efficiency, defining the Phase Smoothness Factor. Lastly, we study how the system behavior is affected certain system parameters and not others. 8 The contributions of this work presented in Chapter 3 are as follows. We model a pinned-free, flexible, tail-like beam immersed in fluid flow using fluid, subject to various leading-edge boundary conditions. We show how these non-conservative loading due to feedback causes the system to lose stability through flutter. We study various actuation-proportional-to-sensing combinations incorporating moment and angle actuation; displacement, slope, and curvature sensing; and both positive and negative sign of feedback gain. We determine the flow velocity at which the systems with zero feedback gain lose stability. We define the criteria for a system losing stability due to feedback from which we, in turn, define the critical stability surfaces to present the results. We present the results for six Cases, one for each combination of actuation and sensing over a large range of feedback locations and flow velocities which show a rich set of stability transitions. We determine the that the impact of flow velocity is much more mild than the significant impact of the location of sensing. We further determine the propulsive characteristics of the waveforms generated by the flutter instability through slender body equations. We then show that while the resulting waveform is composed of four individual traveling waves, we can define a Phase Smoothness Factor metric which can be used to reliably predict the efficiency based on the numerical waveform. Finally we determine that the propulsive characteristics of the beam do not depend on the combination of actuation and sensing by which flutter is produced; they depend only on the values of the dimensionless fluid velocity and critical frequency, which completely define the waveform. In Chapter 4, we aim to provide experimental validation for the feedback-induced flutter of a flexible beam in fluid flow, which was introduced in Chapter 3. To achieve this, we selected a suitable beam with the appropriate flexibility and density and installed strain gauges to measure curvature at a point along the centerline. We then constructed a setup for applying a moment to the beam’s leading edge and all necessary electronics for curvature measurement and torque control. We created a Simulink program to implement the feedback control scheme and ran it through the ControlDesk interface. We tested the system at various flow velocities and with both positive and negative feedback gain to determine the critical feedback gain and frequency at which flutter occurred. To validate our experimental findings, we adapted the moment-proportional-to- 9 curvature system from Chapter 3 to account for experimental considerations and solved it for points corresponding to those of the experiment. We then compare the experimental results and numerical simulations to test the validity of the model used. The contributions of this work presented in Chapter 4 are as follows. We design and construct a system for applying feedback-based actuation to the leading edge of a flexible, tail-like beam in fluid flow proportional to the curvature at a point along its centerline. This involves both the application, sensing, and amplification of strain gauges as well as the control and transmission of torque from a motor to the leading edge of the beam. We design a control scheme that enables the input of the feedback gain which we use to determine when the beam loses stability through flutter at different axial flow velocities. We show that the positive feedback gain behavior is significantly different to the negative feedback behavior for our system. We show that the behavior of the system is loosely correlated with the flow velocity. We modify the moment-proportional-to-curvature boundary condition in Chapter 3 to include the the inertia of the torque application hardware and the time delay that occurs due to the flexibility of the torque application hardware. We show that the numerical model of the experimental setup loses stability through flutter at frequencies very similar to those observed in the experiment. With this, we can state with confidence that the model used is a good approximation of the real system and that feedback-induced flutter of a flexible beam in axial fluid flow is achievable. In Chapter 5, we again examine the flutter characteristics of a beam in fluid flow subject to a follower force. This time, however, we focus on the impact of a follower force on the onset of flutter, given the possibility that the follower force could be generated by internal flow and a fluid jet emanating from a convergent nozzle, while retaining the same equations of motion. We start by defining the model of a cantilever-free beam with a tangential follower force in fluid flow, with the cantilevered end being upstream. Following the methods outlined in Chapter 3, we calculate the equations of motion and boundary conditions of the beam. We then use the same solution process to determine the natural frequencies of the beam, given parameters such as the beam material and geometric properties, fluid flow velocity, and follower force amplitude. We first study the behavior 10 of the system as the external flow velocity is increased with no follower force, over a range of beam-to-fluid mass ratios. We analyze these curves to identify the range of external flow velocities to examine further. We then show how to determine the follower force at which the system becomes unstable. We present surfaces over the mass ratio and external flow velocity plane, indicating the follower force required to induce flutter instability, the critical frequency at which the beam flutters, and the mode from which the beam lost stability. Additionally, we delve into the scaling of the system and observe some interesting convergence. We examine Ziegler’s paradox and its relevance to the problem, and compare the effects of a follower force on stability to the attachment of a nozzle to a fluid-conveying pipe using the same equations. The contributions of this work presented in Chapter 5 are as follows. We model a cantilever-free, flexible, tail-like beam in axial fluid flow with a tangential follower force applied at the free end. We define the equation of motion and boundary conditions of this system and show the method of determining the motion given a set of parameters. We explore the natural frequencies of the system without follower force over the full range of possible beam-to-fluid mass ratios as the fluid flow increases. These results present a rich set of behaviors which define the points from which our follower-force analysis begins. Our analysis of the critical follower force allows us to determine the frequency and mode at which the system loses stability through flutter given each combination of mass ratio and fluid flow velocity. We show that for any given flow velocity, increasing the mass ratio always requires a larger follower force to cause the system to lose stability. We show that by scaling the flow velocity axis by the square root of the mass ratio, we can get all of our results for the critical frequency surface and the critical mode surface to simplify tremendously. Our results for very low flow velocities confirm Ziegler’s paradox regarding the destabilizing nature of low damping on certain systems. Our analysis of incorporating a convergent nozzle to the downstream end of the system shows that the follower force can be feasibly simulated for higher flow velocities but not lower ones given the tremendous difficulty of actualizing a tangential follower force. Concluding remarks and directions for future work are given in Chapter 6. 11 CHAPTER 2 FORCED EXCITATION OF A FLEXIBLE BEAM IN FLUID FLOW In this chapter, we examine the propulsive qualities of a tail-like flexible beam actuated by a dynamic moment over a range of frequencies and flow speeds. We develop a model of an underwater vehicle composed of a rigid body connected to a flexible tail-like beam by a revolute joint. The rigid body is pinned about its center of mass and constrained to move along a channel. We lay out several assumptions used to simplify our model and assist with tractability of the solution. We use the equations of motion for fluid-immersed beams and the dynamics of the rigid body for the boundary conditions to defined the system mathematically. We present a method of solving for the motion of the beam under sinusoidal angle actuation of the revolute joint connecting the rigid body to the flexible beam. We use a set of tractable equations to solve for the thrust and efficiency of the system over the velocity-frequency plane and show that the flexible propulsor has regions of both positive and negative thrust. We also show the behavior of a sample underwater vehicle with fixed drag characteristics as an illustration of a realizable system. The layout of this chapter is organized as follows. The physical system is described in Section2.1 along with a set of assumptions that simplify the mathematical model. The mathematical model is presented in Section2.2 and the method of solution is discussed in Section2.3. For a specific set of system parameters, the thrust and efficiency characteristics of the underwater vehicle, associated with different driving frequencies, are discussed in Section2.4. 2.1 System Description and Assumptions This chapter was previously published as [105], and has been reformatted to meet the require- ments of this dissertation. Consider the underwater vehicle immersed in quiescent fluid, in Fig.2.1. It is comprised of a rigid body and a tail-like flexible beam. The center-of-mass of the rigid body is constrained to translate along a channel and the beam is connected to the rigid body by a revolute joint. The vehicle is propelled by the oscillatory motion of the flexible beam, which is generated by actively 12 controlling the revolute joint. We make several simplifying assumptions in modeling the vehicle. A1. The vehicle is in a state of dynamic equilibrium. The revolute joint of the vehicle translates with constant velocity 𝑈e along the negative 𝑋 axis in the quiescent fluid. We define the reference frame 𝑋𝑌 where the 𝑋 axis is aligned with the channel and the origin is located at the projection of the revolute joint onto the 𝑋 axis; this implies that the 𝑋𝑌 frame translates along the channel. The 𝑋𝑌 frame is therefore an inertial reference frame and 𝑈e denotes the external flow relative to this frame. The underwater vehicle is neutrally buoyant in a fluid of density 𝜌f and its motion is confined to the horizontal plane. A2. The 𝑥 1 𝑦 1 frame is fixed to the rigid body at its center-of-mass. The orientation of the rigid body relative to the 𝑋𝑌 frame is denoted by 𝛼, which is measured positive about the vertical axis. The 𝑥 2 𝑦 2 frame is located at the revolute joint where the 𝑥 2 axis is aligned with the slope of the beam. The angle between the 𝑥 1 and 𝑥 2 axes is 𝛿, which is measured positive about the vertical axis and is assumed to be small. A3. The rigid body is symmetric about the vertical plane that contains the 𝑥 1 axis. It has mass moment of inertia 𝐽 about its center-of-mass, which is located at a distance ℓ from the 𝑌 𝑦1 𝑥2 𝑦2 𝛿 𝑥1 ℓ 𝑈e 𝛼 𝑦(𝑥, 𝑡) 𝜃(𝑥, 𝑡) 𝑋 𝑥 𝐽 𝐿 Figure 2.1: An underwater vehicle comprised of a rigid body and a tail-like flexible beam. 13 revolute joint. For 𝛼 = 0, the projected area of the rigid body in the 𝑌 𝑍 plane is 𝐴𝑟 . The drag coefficient of the rigid body is 𝐶 𝐷 and the drag force acts in the positive 𝑋 direction at a distance 𝑙 𝑑 from the center-of-mass, measured along the 𝑥 1 axis. A4. The flexible beam satisfies Euler-Bernoulli assumptions. It has length 𝐿 and a rectangular cross-section with flexural rigidity 𝐸 𝐼, where 𝐸 and 𝐼 are the Young’s modulus and area moment of inertia, respectively. The fluid volume associated with the motion of the flexible beam yields an added mass per unit length 𝑚 e ; this fluid is assumed to have a uniform velocity profile [81], or an infinitely thin boundary layer. This assumption is made for tractability, see, e.g., [106]. A5. The displacement of a point on the beam at a distance 𝑥 from the revolute joint is denoted by 𝑦(𝑥, 𝑡) in the 𝑋𝑌 frame. The revolute joint has a displacement of 𝑦 0 , 𝑦(0, 𝑡). The slope of the beam is denoted by 𝜃(𝑥, 𝑡) , [𝜕𝑦(𝑥, 𝑡)/𝜕𝑥]. The slope of the beam at 𝑥 = 0, denoted by 𝜃 0 , 𝜃(0, 𝑡), is small. Since 𝜃 0 and 𝛿 are both small, 𝛼 , (𝜃 0 − 𝛿) is small. A6. The net drag on the underwater vehicle is entirely due to the drag on the rigid body, i.e., the flexible tail produces no drag. The drag force acts at a point on the longitudinal axis that lies behind the center-of-mass; this is consistent with hydrodynamically stable bodies. The net thrust produced by the vehicle is generated by the oscillating motion of the flexible tail. 2.2 Dynamic Model Rigid Body Dynamics The free-body diagrams of the rigid body and the flexible beam are shown in Fig.2.2. The reaction forces at the revolute joint are assumed to be F and S along the 𝑋 and 𝑌 axes. The drag force and the reaction force of the channel on the rigid body are denoted by D and R. The reaction moment about the 𝑍 axis is M. Summing the moments about the center-of-mass of the rigid body, we get: M − Sℓ cos 𝛼 + F (ℓ − ℓd ) sin 𝛼 = 𝐽 𝛼¥ (2.1) 14 Substituting cos 𝛼 ≈ 1 and (ℓ − ℓd ) sin 𝛼 ≈ 0, we get M − ℓS = 𝐽 𝛼¥ (2.2) One boundary condition at the revolute joint is geometric: 𝑦 0 = ℓ sin 𝛼 ≈ ℓ𝛼 = ℓ (𝜃 0 − 𝛿) (2.3) The other boundary condition arises from the shear force S and the moment M, which are dependent on the dynamics of the flexible beam. These will be described in the next sub-section. Fluid-Immersed Beam Dynamics The dynamics of the flexible beam in Fig.2.1 is identical to that of a beam in axial flow [52, 91]: 𝜕4 𝑦 2 2 2 𝐸𝐼 + 𝑚 𝑈 2 𝜕 𝑦 + 2𝑚 𝑈 𝜕 𝑦 + (𝑚 + 𝑚 ) 𝜕 𝑦 = 0 (2.4) e e e e e b 𝜕𝑥 4 𝜕𝑥 2 𝜕𝑥𝜕𝑡 𝜕𝑡 2 where 𝑚 e and 𝑚 b denote the mass per unit length of the surrounding fluid (added mass) and the beam, respectively. Including the wake of the rigid body or a boundary layer over the tail would 𝛿 𝑌 𝑋 M 𝛼 ℓ S F F M S 𝛼 D R ℓd Figure 2.2: Free body diagrams of the rigid body and flexible beam shown in Fig.2.1. 15 introduce a constant to the second term of (2.4) [106]. For a uniform velocity profile (A4), this constant is 1. The boundary conditions at the free end of the beam are: 𝜕2 𝑦 𝜕3 𝑦 𝐸𝐼 = 0, 𝐸 𝐼 =0 (2.5) 𝜕𝑥 2 𝑥=𝐿 𝜕𝑥 3 𝑥=𝐿 Using the moment equilibrium (2.2) and the pinned condition (2.3), the boundary conditions at the revolute joint can be expressed as:  3  𝜕2 𝑦 𝜕3 𝑦 𝜕 𝑦 𝑑2𝛿 𝐸𝐼 2 − ℓ𝐸 𝐼 3 =𝐽 − (2.6a) 𝜕𝑥 𝑥=0 𝜕𝑥 𝑥=0 𝜕𝑥𝜕𝑡 2 𝑥=0 𝑑𝑡 2   𝜕𝑦 𝑦 =ℓ −𝛿 (2.6b) 𝑥=0 𝜕𝑥 𝑥=0 We introduce the following change of variables: r s 𝑦 𝑥 𝑚e 𝐿2 𝐸𝐼 𝜈= , 𝑢= , 𝑢 e = 𝑈e , 𝜏=𝑡 𝐿 𝐿 𝐸𝐼 (𝑚 e + 𝑚 b )𝐿 4 to obtain the non-dimensional equation of motion of the beam 𝜕4𝜈 2 p 2 2 + 𝑢 2 𝜕 𝜈 + 2 𝛽𝑢 𝜕 𝜈 + 𝜕 𝜈 = 0 (2.7) e e 𝜕𝑢 4 𝜕𝑢 2 𝜕𝑢𝜕𝜏 𝜕𝜏 2 From (2.5), the non-dimensional boundary conditions at the free end of the beam are 𝜕2𝜈 𝜕3𝜈 = 0 =0 (2.8) 𝜕𝑢 2 𝑢=1 𝜕𝑢 3 𝑢=1 while (2.6a) and (2.6b) provide the boundary conditions at the revolute joint  3  𝜕2𝜈 𝜕3𝜈 𝜕 𝜈 𝑑2𝛿 −𝜅 3 =𝜆 − (2.9a) 𝜕𝑢 2 𝑢=0 𝜕𝑢 𝑢=0 𝜕𝑢𝜕𝜏 2 𝑢=0 𝑑𝜏 2   𝜕𝜈 𝜈 =𝜅 −𝛿 (2.9b) 𝑢=0 𝜕𝑢 𝑥=0 In (2.7) and (2.9), 𝛽, 𝜆, and 𝜅 are non-dimensional parameters: 𝑚e 𝐽 ℓ 𝛽= , 𝜆= , 𝜅= (2.10) 𝑚e + 𝑚b (𝑚 e + 𝑚 b )𝐿 3 𝐿 16 2.3 Method of Solution Exact Solution The system of differential equations (2.7,2.8,2.9) is solved using the procedure from [52,91]. Since the revolute joint is actively controlled, we assume  𝑚 + 𝑚  1/2 e b 𝛿(𝑡) = 𝛿0 e𝑖Ω𝑡 = 𝛿0 e𝑖𝜔𝜏 , 𝜔, Ω𝐿 2 (2.11) 𝐸𝐼 where 𝛿0 and Ω are the driving amplitude and frequency and 𝜔 is the non-dimensional excitation frequency. We assume the response of the flexible beam to have the form 𝜈(𝑢, 𝜏) = 𝜙(𝑢)e𝑖𝜔𝜏 (2.12) where 𝜙(𝑢) is a complex shape function with spatially-varying magnitude and phase. Substituting (2.12) into (2.7) results in p 𝜙′′′′ + 𝑢 2e 𝜙′′ + 2𝑖𝑢 e 𝛽𝜔 𝜙′ − 𝜔2 𝜙 = 0 (2.13) where ( )′ denotes the spatial derivative of ( ) with respect to 𝑢. Applying the boundary conditions in (2.8) and (2.9), we get 𝜙′′(1) = 0 𝜙′′′(1) = 0 (2.14) 𝜙′′(0) − 𝜅𝜙′′′(0) + 𝜔2𝜆𝜙′(0) = 𝜔2𝜆𝛿0 𝜙(0) − 𝜅𝜙′(0) = −𝜅𝛿0 where forcing terms appear on the right hand side. We assume the shape function to be of the form 𝜙(𝑢) = 𝐴1 e𝑧1 𝑢 + 𝐴2 e𝑧2 𝑢 + 𝐴3 e𝑧3 𝑢 + 𝐴4 e𝑧4 𝑢 (2.15) where 𝑧𝑖 = 𝑧𝑖 (𝜔) are the roots of the characteristic equation of (2.13). Given 𝑢 e , 𝛽, and 𝜔, substituting (2.15) into (2.14) yields 17       𝑧 21 e𝑧1 𝑧 22 e𝑧2 𝑧 23 e𝑧3 𝑧 24 e𝑧4   𝐴1   0            𝑧 1 e 1 𝑧 2 e 2 𝑧 3 e 3 𝑧 4 e 4   3 𝑧 3 𝑧 3 𝑧 3 𝑧 𝐴2   0     =  𝛿0 (2.16)   𝐴3   −𝜔2𝜆   𝜂1 𝜂2 𝜂3 𝜂4                𝜁1 𝜁2 𝜁3 𝜁4   𝐴4   −𝜅   where 𝜂𝑖 , 𝜁𝑖 , 𝑖 = 1, 2, 3, 4, are given by the relations 𝜂𝑖 = 𝑧𝑖2 − 𝜅𝑧𝑖3 + 𝜆𝜔2 𝑧𝑖 , 𝜁𝑖 = 1 − 𝜅𝑧𝑖 We can solve (2.16) for the 𝐴𝑖 to determine the shape function 𝜙(𝑢). The response of the flexible beam takes the form X 4 X 4 𝜈(𝑢, 𝜏) = 𝐴𝑛 e𝑧𝑛 𝑢 e𝑖𝜔𝜏 = 𝐴𝑛 eRe[𝑧𝑛 ]𝑢 e𝑖 ( Im[𝑧𝑛 ]𝑢+𝜔𝜏 ) . (2.17) 𝑛=1 𝑛=1 The forced response in (2.17) is the sum of four terms wherein each term is a product of two exponential terms: the first exponential term is bounded because 𝑢 is bounded and the second exponential term is periodic as its exponent is imaginary. Each term in (2.17) is a traveling wave. The forced response in (2.17) is also the steady-state response provided that the free response is transient and decays out with time. A Case Study We consider the rigid body in Fig.2.1 to be an ellipsoid with length along the 𝑥 1 , 𝑦 1 and 𝑧 1 axes equal to 0.4 m, 0.06 m, and 0.12 m, such that 𝐴𝑟 = 5.65 × 10−3 m2 , ℓ = 0.2 m. Its density is assumed to be the same as that of the surrounding fluid, which is water and equal to 𝜌f = 1000 kg/m3 ; this yields 𝐽 = 0.0123 kg m2 . The flexible beam has length 𝐿 = 0.45 m and a rectangular cross-section of width 0.001 m and height 0.1 m, such that 𝐼 = 8.333 × 10−12 m4 . The material of the beam is assumed to be Cirlex® , for which 𝜌b = 1420 kg/m3 , and 𝐸 = 2.7 GPa. For these dimensions, the linear density of the beam is 𝑚 b = 0.142 kg/m; following standard approximations in the literature [81, 107], the 18 0.0 𝑢 Re[𝜈(𝑢, 0)] Re[𝜈(𝑢, 𝜋/2𝜔)] Re[𝜈(𝑢, 𝜋/𝜔)] Re[𝜈(𝑢, 2𝜋/𝜔)] Re[𝜈(𝑢, 𝜋/4𝜔)] Re[𝜈(𝑢, 3𝜋/4𝜔)] Re[𝜈(𝑢, 5𝜋/4𝜔)] Re[𝜈(𝑢, 3𝜋/2𝜔)] Re[𝜈(𝑢, 7𝜋/4𝜔)] 1.0 Figure 2.3: A sequence of images showing the shape of the flexible beam Re[𝜈(𝑢, 𝜏)] during forced vibration for 𝜔 = 40.0 and 𝑢 e = 2.0 over one cycle of oscillation. The revolute joint is at 𝑢 = 0. linear density of the added mass of water is 𝑚 e = 7.854 kg/m. This corresponds to the mass of the cylinder of fluid surrounding the flapping beam. For the set of dimensional parameters provided above, the non-dimensional parameters in (2.10) are 𝛽 = 0.9822, 𝜆 = 0.0169, 𝜅 = 0.4444 A sequence of images showing the traveling wave nature of the flexible beam under forced vibration is shown in Fig.2.3 for the particular case of 𝜔 = 40.0 and 𝑢 e = 2.0. 19 2.4 Investigation of Propulsive Performance Thrust The method proposed by Lighthill for “slender fish" [29] can be used to estimate the thrust generated by the flexible beam in fluid flow. As in [31, 52], the more general case of a beam is used rather than the reduced case of a “tapered fish" that has neither mass nor area at its leading edge. The time-averaged net thrust produced by the tail-like flexible beam is   2  2  2  2  1  𝜕𝑦 𝜕𝑦 𝜕𝑦 𝜕𝑦  F̄ = 𝑚 e  − 𝑈e − − 𝑈e   (2.18) 2  𝜕𝑡 𝜕𝑥 𝑥=𝐿 𝜕𝑡 𝜕𝑥 𝑥=0    where the notation (·) is the cycle-average of (·). The thrust can be non-dimensionalized using the change of variables in Section 2.2 and computed over one cycle using 2 Z 2𝜋/𝜔 ( "   2  2# "  2  2# ) F̄𝐿 𝜔 𝜕𝜈 𝜕𝜈 𝜕𝜈 𝜕𝜈 F̄∗ = = − 𝑢e − − 𝑢e 𝑑𝜏 (2.19) 𝐸𝐼 4𝜋 0 𝜕𝜏 𝜕𝑢 𝜕𝜏 𝜕𝑢 𝑢=1 𝑢=0 The function 𝜈(𝑢, 𝜏) is solved using the procedure in Section 2.3 and has both real and imaginary parts. Since only the real part physically contributes to the thrust, Re[𝜈] is used in place of 𝜈 in (2.19) and has the form X4 Re[𝜈(𝑢, 𝜏)] = eRe[𝑧𝑛 ]𝑢 {Re [ 𝐴𝑛 ] cos (Im [𝑧 𝑛 ] 𝑢 + 𝜔𝜏) − Im [ 𝐴𝑛 ] sin (Im [𝑧 𝑛 ] 𝑢 + 𝜔𝜏)} 𝑛=1 (2.20) The amplitude of the sinusoidal excitation will determine the magnitude of the thrust; we chose 𝛿0 = 12 deg (≈ 0.2 rad). This does not universally guarantee small beam displacements due to resonance; to conform with Euler-Bernoulli assumptions, we restrict our investigation to the domain where the maximum elongation of the beam does not exceed 5%, i.e., Z1 q  sup 1 + Re[𝜈′] 2 𝑑𝑢 ≤ 1.05 (2.21) 𝜏∈[0,2𝜋/𝜔] 0 20 For the case study discussed in Section 2.3, Fig.2.4 shows a contour plot of the nondimensional time-averaged thrust produced by the flexible beam in the domain 𝜔 ∈ (0, 100] and 𝑢 e ∈ (0, 10]. The hatched regions show areas of negative thrust and solid grey regions show areas where the elongation exceeds 5% due to resonance. The general trend of Fig.2.4 shows higher thrust at higher frequency with regions of significantly greater thrust due to resonance. In the neighborhood of resonance, denoted by the solid grey areas in Fig.2.4, the tip of the tail moves significantly more which results in greater thrust - this follows from (2.19). As the frequency increases beyond a region of resonance, the tail movement, and thus the thrust, drop off quickly. It should be noted that 𝜔 and 𝑢 e are treated as independent variables in the results presented in Fig.2.4. In reality, 𝑢 e is a consequence of the dynamic equilibrium where the time-averaged thrust F̄ is equal to the time-averaged drag D̄ for a given 𝜔. The time-averaged drag is assumed to be 10 0.5 1.0 0.1 1.0 0.25 2.5 𝑢e 5 0.5 5.0 1.0 0.1 0.25 0 0 50 100 𝜔 Figure 2.4: Contour plot of the nondimensional time-averaged thrust F̄∗ in the 𝜔-𝑢 e plane; hatched areas signify regions with negative thrust and the three solid grey areas indicate regions around resonance where the elongation exceeds 5%. 21 1 0.8 0.6 0.4 0.2 0 0 5 10 15 20 25 0.6 0.58 0.56 0.54 0.52 0.5 0 5 10 15 20 25 Figure 2.5: (a) Locus of dynamic equilibrium for the underwater vehicle described in Section 2.3, (b) Froude efficiency along the locus of dynamic equilibrium. D̄ = (1/2)𝜌f 𝐶 𝐷 𝐴𝑟 𝑈e2 (2.22) with 𝐶 𝐷 = 0.1 [108]. In light of this, we now look at the locus of points where F̄ = D̄. A dimensional plot of the locus is shown in Fig.2.5(a). For a given Ω in this plot, if 𝑈e lies below (above) the locus, it will accelerate (decelerate) until its velocity reaches the locus. Consider the case of the vehicle accelerating from rest using Ω = 10 rad/s. Because its velocity initially lies below the locus, the vehicle accelerates to the speed on the locus which matches the frequency, 𝑈e = 0.41 m/s. By varying Ω, the vehicle can attain higher and lower 𝑈e for a given amplitude 𝛿0 . The general trend of Fig.2.5(a) shows higher 𝑈e at higher frequencies; 𝑈e increases rapidly as the system approaches resonance and decreases post-resonance. This is in accord with the trends observed in Fig.2.4. Each point on the locus can also be associated with a specific value of Froude efficiency [29]. The derivation of this efficiency follows. Efficiency To compute the propulsive efficiency, we use an expression derived by Lighthill for the time- averaged power P̄ required to generate the displacements 𝑦(𝑥, 𝑡) which produce the thrust: 22 "     # 𝜕𝑦 𝜕𝑦 𝜕𝑦 𝜕𝑦 𝜕𝑦 𝜕𝑦 P̄ = 𝑈e 𝑚 e + 𝑈e − + 𝑈e (2.23) 𝜕𝑡 𝜕𝑡 𝜕𝑥 𝑥=𝐿 𝜕𝑡 𝜕𝑡 𝜕𝑥 𝑥=0 This expression can be non-dimensionalized using the same change of variables as before to get Z 2𝜋/𝜔         P̄𝐿 3 𝜔 𝜕𝜈 𝜕𝜈 𝜕𝜈 𝜕𝜈 𝜕𝜈 𝜕𝜈 P̄∗ = = + 𝑢e − + 𝑢e 𝑑𝜏 (2.24) 𝐸𝐼 2𝜋 0 𝜕𝜏 𝜕𝜏 𝜕𝑢 𝑢=1 𝜕𝜏 𝜕𝜏 𝜕𝑢 𝑢=0 The Froude efficiency [29], defined as the ratio between the power of forward motion through the fluid and the power required to generate the thrust by the motion of fluid, is F̄𝑈e F̄∗ 𝑢 e 𝜂= = ∗ (2.25) P̄ P̄ The efficiency 𝜂 is defined when the vehicle is moving with a constant speed 𝑈e , so we compute the efficiency only for points on the locus of dynamic equilibrium. Figure 2.5 (b) shows the Froude efficiency as a function of Ω; the corresponding value of 𝑈e can be found from Fig.2.5 (a). For example, with a driving frequency of Ω = 10 rad/s, the underwater vehicle will have a steady state velocity of 𝑈e = 0.41 m/s with an efficiency of 𝜂 = 0.56. 2.5 Conclusion The dynamics of a flexible tail-like structure, connected to a rigid body by an actively controlled revolute joint, can be analyzed as a fluid-immersed beam in axial flow. The rigid body imposes boundary conditions at one end of the beam while the other end is free. Subject to simplifying assumptions, the dynamics of the flexible beam are analytically tractable and result in traveling waves. These traveling waves produce thrust that can propel the underwater vehicle by overcoming the drag of the rigid body. The efficiency of the thrust varies as a function of the flow velocity and oscillation frequency of the revolute joint. The locus of dynamic equilibrium points, where thrust and drag forces balance each other, was obtained for a sample vehicle; the efficiency values on the locus are found to exceed 50%. Since the analysis is based on Euler-Bernoulli beam theory, simulations were carried out using a small amplitude of the revolute joint. Nevertheless, the deflection of the flexible tail-like structure 23 becomes large near resonance and such regions were therefore excluded from our investigation. A more accurate model of the fluid-structure interaction is necessary to investigate the behavior of the dynamic system over the complete domain. In certain regions of the velocity-frequency plane, the flexible beam produces negative thrust, implying that it acts as a brake. It should be noted that negative thrust can potentially produce backward motion, but the current model needs to be expanded to account for negative flow velocity. Similar to the thrust, the power can also be negative. While this condition is not sustainable for a self-propelled underwater vehicle, it may be possible to exploit it for energy extraction if the underwater vehicle remains anchored in a flow such as a stream or river. 24 CHAPTER 3 FEEDBACK-INDUCED FLUTTER INSTABILITY OF A FLEXIBLE BEAM IN FLUID FLOW In this chapter we aim to study how different forms of feedback could be used to instigate flutter behavior in a beam in fluid flow. We begin by describing a model of a pinned-free beam which is actuated at its leading edge by an moment or angle condition that is proportional to the displacement, slope, or curvature of the beam at some point along its length. We formulate the equations of motion and boundary conditions which include three default boundary conditions as well as one that depends on which method of actuation and sensing is used. We provide a meethod of non-dimensionalizing the equations and boundary conditions to simplify their analysis. We work through the method of solution required to determine the natural frequencies of the beam given the set of parameters such as beam properties, location of sensing, fluid flow velocity, and feedback gain. We start by determining the behavior of the system as the external flow velocity is increased with no feedback gain which corresponds to a pinned-free beam for the moment actuation case and a cantilevered beam for the angle actuation case. From this, we show the method of determining the feedback gain at which the system loses stability. Of the twelve possible cases of actuation, sensing, and sign of feedback gain, we choose six cases that are illustrative of the behavior possible over the span of external flow velocities and locations of sensing. We describe several patterns that are observable in the critical stability and critical frequency surfaces and compare and contrast the behavior of the different cases. Following this, we explore the applicability of the flutter oscillations produced through feedback as a propulsor for an underwater vehicle. We explore the thrust and efficiency characteristics to see how they depend on the chosen parameters. We additionally explore how the smoothness of the waveform affects the efficiency, defining a Phase Smoothness Factor in the process. Finally, we study how the behavior of the system is determined by the external flow velocity and frequency of oscillation, independent of the other parameters of the system. The remainder of this chapter is organized as follows. The problem formulation is described in 25 Section 3.1 and the analytical method of solution is discussed in Section 3.2. The investigation of the flutter instability is defined in Section 3.3, and the numerical procedure is laid out in Section 3.3. Parameters for the system are given in Section 3.3 as well as the terminology for the results that are presented in 3.4. Section 3.5 lays out how this system could be applied to the propulsion of an underwater vehicle and also compares the dynamic and hydrodynamic characteristics at several points of the results. This chapter was previously published as [109], and has been reformatted to meet the require- ments of this dissertation. 3.1 Problem Formulation Consider the fluid-immersed flexible beam in Fig.3.1. The beam has length 𝐿, a rectangular cross-section with width 𝑤 and height ℎ, mass per unit length 𝑚 b , and Young’s modulus of elasticity 𝐸. The upstream end of the beam is connected to a fixed point by a revolute joint, which is actively controlled; the downstream end of the beam is free. The fluid is inviscid and flows with constant velocity 𝑈e . The equation of motion of the beam, ignoring gravitational, viscous, pressurization and tensile effects, is as follows [52, 91]: 𝜕 4 𝑦(𝑥, 𝑡) 2 2 2 𝐸𝐼 + 𝑚 𝑈 2 𝜕 𝑦(𝑥, 𝑡) + 2𝑚 𝑈 𝜕 𝑦(𝑥, 𝑡) + (𝑚 + 𝑚 ) 𝜕 𝑦(𝑥, 𝑡) = 0 (3.1) 𝑒 e 𝑒 e e b 𝜕𝑥 4 𝜕𝑥 2 𝜕𝑥𝜕𝑡 𝜕𝑡 2 𝑦 point of sensing 𝑤 𝑈e 𝑥 b 𝑥 𝐿 Figure 3.1: A flexible beam, connected at one end to a fixed point by a revolute joint and free at the other end, is immersed in a fluid flowing with constant velocity 𝑈e . 26 where 𝑦(𝑥, 𝑡) is the displacement of the beam, 𝐼 = (ℎ𝑤 3 /12) is the area moment of inertia of the beam, and 𝑚 𝑒 is the mass per unit length of the external fluid. The mass per unit length of the external fluid is approximated as the mass of water within the cylinder of unit length circumscribing the beam cross-section [107]. The displacement boundary condition of the pinned end of the beam and the shear and moment boundary conditions of the free end of the beam are given as 𝜕 2 𝑦(𝐿, 𝑡) 𝜕 3 𝑦(𝐿, 𝑡) 𝑦(0, 𝑡) = 𝐸 𝐼 = 𝐸 𝐼 =0 (3.2) 𝜕𝑥 2 𝜕𝑥 3 We assume that an actuator located at the revolute joint can either apply a bending moment or impose an angle condition on the beam at 𝑥 = 0. Furthermore, the bending moment or the imposed angle can be based on feedback: proportional to the curvature, slope, or displacement of the beam at some point along its length 𝑥 = b 𝑥 . Therefore, the final boundary condition at 𝑥 = 0 can take one of six forms depending on the two modes of actuation and three modes of feedback: 𝜕 2 𝑦(0, 𝑡) 𝜕 2 𝑦(b 𝑥 , 𝑡) 𝐸𝐼 = 𝐶m,𝑐 moment ∝ curvature (3.3a) 𝜕𝑥 2 𝜕𝑥 2 𝜕 2 𝑦(0, 𝑡) 𝜕𝑦(b 𝑥 , 𝑡) 𝐸𝐼 = 𝐶m,𝑠 moment ∝ slope (3.3b) 𝜕𝑥 2 𝜕𝑥 𝜕 2 𝑦(0, 𝑡) 𝐸𝐼 = 𝐶m,𝑑 𝑦(b 𝑥 , 𝑡) moment ∝ displacement (3.3c) 𝜕𝑥 2 𝜕𝑦(0, 𝑡) 𝜕 2 𝑦(b𝑥 , 𝑡) = 𝐶a,𝑐 angle ∝ curvature (3.3d) 𝜕𝑥 𝜕𝑥 2 𝜕𝑦(0, 𝑡) 𝜕𝑦(b 𝑥 , 𝑡) = 𝐶a,𝑠 angle ∝ slope (3.3e) 𝜕𝑥 𝜕𝑥 𝜕𝑦(0, 𝑡) = 𝐶a,𝑑 𝑦(b 𝑥 , 𝑡) angle ∝ displacement (3.3f) 𝜕𝑥 where 𝐶m,𝑐 , 𝐶m,𝑠 , 𝐶m,𝑑 , 𝐶a,𝑐 , 𝐶a,𝑠 and 𝐶a,𝑑 are feedback gains of appropriate dimensions. Remark 1 For the purpose of theoretical development, it is assumed that the curvature, slope, or displacement of the beam at an arbitrary point along its length can be measured directly using sensors or estimated from sensor measurements. 27 Remark 2 The actuator at the revolute joint can directly apply a moment proportional to the curvature, slope, or displacement at some point along the length of the beam. To impose an angle condition, whereby the revolute joint angle is proportional to the curvature, slope, or displacement, a feedback controller must be designed to drive the actuator to track the desired angle. We introduce the following change of variables: r s 𝑦 𝑥 b𝑥 𝑚e 𝐸𝐼 𝑣= , 𝑢= , 𝛾= , 𝑢 e = 𝑈e 𝐿 , 𝜏=𝑡 𝐿 𝐿 𝐿 𝐸𝐼 (𝑚 e + 𝑚 b )𝐿 4 to obtain the non-dimensional equation of motion of the beam 𝜕4𝑣 2 p 2 2 + 𝑢 2 𝜕 𝑣 +2 𝛽𝑢 𝜕 𝑣 + 𝜕 𝑣 = 0 (3.4) e e 𝜕𝑢 4 𝜕𝑢 2 𝜕𝑢𝜕𝜏 𝜕𝜏 2 where 𝛽 is the mass fraction: 𝑚e 𝛽= 𝑚e + 𝑚b From (3.2), the non-dimensional displacement boundary condition of the pinned end of the beam and the non-dimensional natural boundary conditions of the free end of the beam are given as 𝜕 2 𝑣(1, 𝜏) 𝜕 3 𝑣(1, 𝜏) 𝑣(0, 𝜏) = = = 0. (3.5) 𝜕𝑢 2 𝜕𝑢 3 From (3.3), the actuator-imposed boundary condition at the revolute joint takes one of the following six forms: 28 𝜕 2 𝑣(0, 𝜏) 𝜕 2 𝑣(𝛾, 𝜏) = 𝑐 m,𝑐 moment ∝ curvature (3.6a) 𝜕𝑢 2 𝜕𝑢 2 2 𝜕 𝑣(0, 𝜏) 𝜕𝑣(𝛾, 𝜏) = 𝑐 m,𝑠 moment ∝ slope (3.6b) 𝜕𝑢 2 𝜕𝑢 2 𝜕 𝑣(0, 𝜏) = 𝑐 m,𝑑 𝑣(𝛾, 𝜏) moment ∝ displacement (3.6c) 𝜕𝑢 2 𝜕𝑣(0, 𝜏) 𝜕 2 𝑣(𝛾, 𝜏) = 𝑐 a,𝑐 angle ∝ curvature (3.6d) 𝜕𝑢 𝜕𝑢 2 𝜕𝑣(0, 𝜏) 𝜕𝑣(𝛾, 𝜏) = 𝑐 a,𝑠 angle ∝ slope (3.6e) 𝜕𝑢 𝜕𝑢 𝜕𝑣(0, 𝜏) = 𝑐 a,𝑑 𝑣(𝛾, 𝜏) angle ∝ displacement (3.6f) 𝜕𝑢 where the non-dimensional feedback gains in (3.6) are related to their dimensional counterparts by the relations 𝐶m,𝑐 𝐿𝐶m,𝑠 𝐿 2𝐶m,𝑑 𝐶a,𝑐 𝑐 m,𝑐 = , 𝑐 m,𝑠 = , 𝑐 m,𝑑 = , 𝑐 a,𝑐 = , 𝑐 a,𝑠 = 𝐶a,𝑠 , 𝑐 a,𝑑 = 𝐿𝐶a,𝑑 𝐸𝐼 𝐸𝐼 𝐸𝐼 𝐿 3.2 Method of Solution To solve (3.4) for the boundary conditions in (3.5) and (3.6), we followed the procedure introduced in [91] and used in [52]. In particular, we assume the following separable form for 𝑣(𝑢, 𝜏): 𝑣(𝑢, 𝜏) = 𝑓 (𝑢)𝑒𝑖Ω𝜏 (3.7) where Ω is the non-dimensional frequency of oscillation. Substitution of (3.7) into (3.4) and (3.5) yields p 𝑓 ′′′′(𝑢) + 𝑢 2e 𝑓 ′′(𝑢) + 2𝑢 e 𝛽 𝑖Ω 𝑓 ′(𝑢) − Ω2 𝑓 (𝑢) = 0 (3.8) 𝑓 (0) = 𝑓 ′′(1) = 𝑓 ′′′(1) = 0 (3.9) 29 while substitution of (3.7) into (3.6) yields 𝑓 ′′(0) = 𝑐 m,𝑐 𝑓 ′′(𝛾) moment ∝ curvature (3.10a) 𝑓 ′′(0) = 𝑐 m,𝑠 𝑓 ′(𝛾) moment ∝ slope (3.10b) 𝑓 ′′(0) = 𝑐 m,𝑑 𝑓 (𝛾) moment ∝ displacement (3.10c) 𝑓 ′(0) = 𝑐 a,𝑐 𝑓 ′′(𝛾) angle ∝ curvature (3.10d) 𝑓 ′(0) = 𝑐 a,𝑠 𝑓 ′(𝛾) angle ∝ slope (3.10e) 𝑓 ′(0) = 𝑐 a,𝑑 𝑓 (𝛾) angle ∝ displacement (3.10f) Since (3.8) is an ordinary differential equation with constant coefficients, the solution of 𝑓 (𝑢) is assumed to be of the form 𝑓 (𝑢) = 𝐴𝑒 𝑧𝑢 ; this results in the characteristic equation p 𝑧 4 + 𝑢 2e 𝑧 2 + 2𝑢 e 𝛽 𝑖Ω𝑧 − Ω2 = 0 (3.11) For specific values of 𝑢 e and 𝛽, (3.11) provides four roots of 𝑧 𝑛 , 𝑛 = 1, 2, 3, 4, which are functions of Ω. The solution of 𝑓 (𝑢) takes the form 𝑓 (𝑢) = 𝐴1 𝑒 𝑧1 𝑢 + 𝐴2 𝑒 𝑧2 𝑢 + 𝐴3 𝑒 𝑧3 𝑢 + 𝐴4 𝑒 𝑧4 𝑢 (3.12) Substitution of the boundary conditions in (3.9) and (3.10) yields       1 1 1 1   𝐴1  0           𝑧 𝑒 1 𝑧 𝑒 2 𝑧 𝑒 3 𝑧 𝑒 4   𝐴2  0 2 𝑧 2 𝑧 2 𝑧 2 𝑧  1 2 3 4       =  (3.13) 𝑧3 𝑒 𝑧1 𝑧3 𝑒 𝑧2 𝑧3 𝑒 𝑧3 𝑧3 𝑧4       1 2 3 4 𝑒   𝐴3  0            𝛿1 𝛿2 𝛿3 𝛿4   𝐴4  0      | {z } Z where 𝛿𝑛 , 𝑛 = 1, 2, 3, 4, are defined as follows 30      𝑧 2𝑛 − 𝑐 m,𝑐 𝑧 2𝑛 𝑒 𝑧𝑛 𝛾 : moment ∝ curvature        𝑧 2𝑛 − 𝑐 m,𝑠 𝑧 𝑛 𝑒 𝑧𝑛 𝛾 : moment ∝ slope       𝑧 2𝑛 − 𝑐 m,𝑑 𝑒 𝑧𝑛 𝛾 : moment ∝ displacement 𝛿𝑛 ,     𝑧 𝑛 − 𝑐 a,𝑐 𝑧 2𝑛 𝑒 𝑧𝑛 𝛾 : angle ∝ curvature        𝑧 𝑛 − 𝑐 a,𝑠 𝑧 𝑛 𝑒 𝑧𝑛 𝛾 : angle ∝ slope      𝑧 𝑛 − 𝑐 a,𝑑 𝑒 𝑧𝑛 𝛾 : angle ∝ displacement  For each of the six modes of actuation and feedback combinations, non-trivial solutions of (3.13) can be obtained by solving the transcendental equation det (Z) = 0. For specific values of 𝑢 e , 𝛽, 𝛾, and the appropriate feedback gain 𝑐 ∈ {𝑐 m,𝑐 , 𝑐 m,𝑠 , 𝑐 m,𝑑 , 𝑐 a,𝑐 , 𝑐 a,𝑠 , 𝑐 a,𝑑 } (3.14) the transcendental equation can be solved numerically to get the complex frequencies Ω𝑖 , 𝑖 = 1, 2, · · ·, and the 𝑧 𝑛 terms, 𝑛 = 1, 2, 3, 4, for each Ω𝑖 . 3.3 Investigation of Flutter Instability Critical Stability While Section 3.2 provides the frequencies of oscillation, Ω𝑖 , 𝑖 = 1, 2, · · ·, for specific values of 𝑢 e , 𝛽, 𝛾, and 𝑐1, we seek to find the critical stability points where the system loses stability through flutter. For a particular Ω and corresponding 𝑧 𝑛 terms, 𝑛 = 1, 2, 3, 4, the solution of (3.4), (3.5), and (3.6) can be obtained by substituting (3.12) into (3.7): X 4 X 4 𝑣(𝑢, 𝜏) = 𝐴𝑛 𝑒 𝑧𝑛 𝑢 𝑒𝑖Ω𝜏 = 𝑒 −I𝑚[Ω]𝜏 𝐴𝑛 𝑒 R𝑒[𝑧𝑛 ]𝑢 𝑒𝑖 {I𝑚[𝑧𝑛 ]𝑢+R𝑒[Ω]𝜏 } (3.15) 𝑛=1 𝑛=1 where the coefficients 𝐴𝑛 , 𝑛 = 1, 2, 3, 4, can be obtained from the null space of Z in (3.13). It is clear from (3.15) that the stability of 𝑣(𝑢, 𝜏) is dependent on the exponential term outside the 1The discussion here is general and applies to all six modes of actuation and feedback, i.e., 𝑐 can be any element of the set in (3.14) 31 summation; if I𝑚[Ω] < 0, this term is unbounded as 𝑡 → ∞. Therefore, the point at which I𝑚[Ω] changes sign from positive to negative represents the onset of flutter instability. The first exponential term inside the summation is bounded because 𝑢 is bounded; the second exponential term yields periodic motion, because it has an imaginary exponent. It should be noted that (3.15) describes the solution for one non-dimensional frequency Ω. At the flutter instability point, one specific value of Ω, Ω = Ωc𝑟 , satisfies I𝑚[Ω] = 0 whereas all other Ω values satisfy I𝑚[Ω] > 0. The frequency Ωc𝑟 is real and is defined as the critical frequency. Since 𝑒 −I𝑚[Ω]𝜏 → 0 as 𝜏 → ∞ for all Ω 6= Ωc𝑟 , the complete solution at the flutter instability point takes the form X4 𝑣(𝑢, 𝜏) = 𝐴𝑛 𝑒 R𝑒[𝑧𝑛 ]𝑢 𝑒𝑖 {I𝑚[𝑧𝑛 ]𝑢+R𝑒[Ωc𝑟 ]𝜏 } (3.16) 𝑛=1 Since the imaginary exponent in (3.16) is a function of both 𝑢 and 𝜏, the above equation represents a traveling waveform. Remark 3 In the context of a fluid-immersed slender body, Lighthill [29] established that a trav- eling wave can generate positive thrust if the phase velocity of the wave is greater than the fluid velocity. Since the expression in (3.16) is comprised of four waveforms with different, spatially variable amplitudes and phase velocities, deriving a condition for positive thrust is not straight- forward. The propulsive characteristics of the waveform in (3.16) will be discussed in Section 3.5. Numerical Procedure We first determine the natural frequencies Ω𝑖 , 𝑖 = 1, 2, · · ·, for the unforced system, i.e., the system with 𝑢 e = 0 and 𝑐 = 0. The set of natural frequencies are determined separately for the two cases where the moment applied at the revolute joint is zero - pinned boundary condition; and the angle specified at the revolute joint is zero - cantilevered boundary condition. Unlike the cantilevered boundary condition, the pinned boundary condition will include the rigid-body mode; this requires 32 us to include Ω0 = 0 in the set of natural frequencies for the pinned case. We now introduce the following definition: Frequency Band: The range of frequencies in (0, Ω1 ] is defined as the first frequency band Π1 . The range of frequencies in (Ω 𝑗−1 , Ω 𝑗 ] is defined as the 𝑗-th frequency band Π 𝑗 , 𝑗 ≥ 2. To determine the critical stability points, we fix the value of 𝛽 and vary 𝑢 e and 𝛾 over some domain. For each point in this domain, we solve for the critical feedback gain 𝑐 = 𝑐 c𝑟 , which causes the system to lose stability through flutter. These points define a surface, which we refer to as the critical stability surface. Each point on the critical stability surface corresponds to a critical frequency Ωc𝑟 ; these points define a critical frequency surface. The critical stability and frequency surfaces are obtained as follows: For a specific value of 𝛾, we start with 𝑐 = 0 and 𝑢 e = 0.1. We use the first eleven natural frequencies of the beam, Ω 𝑘 , 𝑘 = 0, 1, 2, · · · , 10, for the pinned case and the first ten natural frequencies of the beam, Ω 𝑘 , 𝑘 = 1, 2, · · · , 10 for the cantilevered case, as the initial guesses to solve for the eigenfrequencies as the magnitude of 𝑐 is gradually increased. The process is continued until one of the Ω 𝑘 ’s satisfies the condition I𝑚[Ω 𝑘 ] = 0. This provides the value of 𝑐 c𝑟 and Ωc𝑟 for 𝑢 e = 0.1 and the specific value of 𝛾; the value of 𝑘 denotes the mode of flutter instability, which will be formally defined later. The process is repeated by gradually incrementing the value of 𝑢 e while keeping the value of 𝛾 fixed; the process is terminated when the value of 𝑐 c𝑟 is uniformly zero2. To obtain the critical stability and frequency surfaces, the overall process is repeated on a fine mesh grid for 𝛾. Simulation Environment We will investigate flutter instability for 𝛽 = 0.98223, 2This signifies that the external flow alone causes the beam to lose stability, much like a flag fluttering in the wind. 3This value of 𝛽 is chosen based on a dimensional example that we will consider later in Section 3.5. 33     [0.1, 0.9] : feedback based on curvature     𝛾 ∈ [0.1, 1.0] : feedback based on slope       [0.3, 1.0] : feedback based on displacement  and    [0.1, 9.0]  : actuator applies a bending moment 𝑢e ∈   [0.1, 19.0]  : actuator imposes an angle condition A different range of 𝛾 was chosen for each of the feedback modes. Since the beam has zero curvature at its free end, we restrict the upper bound to 0.9 for curvature feedback; since the beam has zero displacement at the revolute joint, we restrict the lower bound to 0.3 for displacement feedback. The other bounds were chosen such that the values of the critical feedback gain 𝑐 c𝑟 at the boundaries were not inordinately large compared to those within the bounds. The upper bounds on 𝑢 e were chosen based on the finding that the beam loses stability due to external flow alone at 𝑢 e = 𝑢 e,𝑐𝑟 = 8.99 when the actuator applies a bending moment equal to zero, i.e., pinned boundary condition; and 𝑢 e = 𝑢 e,𝑐𝑟 = 18.12 when the actuator imposes an angle equal to zero, i.e., cantilevered boundary condition. These critical velocities for the pinned and cantilevered boundary conditions are shown in the Argand diagrams in Fig.3.2. The Argand diagrams show the locus of the first few eigenfrequencies as 𝑢 e is increased from zero; each branch starts at a natural frequency of the system Ω 𝑘 . The procedure for computing the critical stability points, described in Section 3.3, can now be better explained with the help of the Argand digrams in Fig.3.2. A specific value of 𝑢 e = 𝑢 ∗e , corresponds to a specific point on each branch of the Argand diagram; note that these points correspond to 𝑐 = 0 and therefore the value of 𝛾 is immaterial. For a specific value of 𝛾 = 𝛾 ∗ , increasing the value of 𝑐 from zero results in eleven (ten) loci of the eigenfrequencies for the moment actuation (angle actuation) case that start on each of the branches of the appropriate Argand diagram at the points corresponding to 𝑢 e = 𝑢 ∗e and 𝑐 = 0. The critical value of 𝑐 = 𝑐 c𝑟 , 34 16 20 𝑘 =3 𝑘 =3 𝑘 =4 𝑘 =3 𝑘 =1 Im(Ω) Im(Ω) 𝑘 =1 𝑘 =2 8 10 𝑘 =2 𝑘 =2 𝑘 =2 𝑢 e,𝑐𝑟 = 8.99 𝑢e = 9 𝑢 e,𝑐𝑟 = 18.12 𝑘 =0 𝑘 =1 Ω2 Ω3 𝑘 =1 0 0 0 Ω0 Ω1 50 100 0 Ω1 Ω2 50 Ω3 100 Re(Ω) Re(Ω) (a) (b) Figure 3.2: Argand diagrams for the beam without feedback for (a) pinned boundary conditions for 𝑢 e = [0, 15] and (b) cantilevered boundary conditions for 𝑢 e = [0, 20]. For the pinned case, the locus originating at Ω0 moves up along the imaginary axis, then downward, and upward again; then, it meets the locus originating at Ω1 and breaks away into the complex plane resulting in instability for 𝑢 e,𝑐𝑟 = 8.99. corresponds to the lowest value of 𝑐 for which one of the eigenfrequencies satisfy I𝑚[Ω] = 0. We now introduce the following definition: Single Mode of Flutter Instability: The system loses stability through the 𝑘-th mode of flutter if the eigenfrequency satisfying I𝑚[Ω] = 0 originated on the 𝑘-th branch in the Argand diagram of Fig.3.2. The above definition implicitly assumes that the loci of the eigenfrequencies do not intersect each other prior to satisfying the condition I𝑚[Ω] = 0. To account for the possibility of intersection of loci, we introduce the following definition: Coupled Mode of Flutter Instability: The system loses stability through 𝑘 1 -𝑘 2 mode of flutter if the eigenfrequency satisfying I𝑚[Ω] = 0 can be traced back to the intersection of two loci that originated on the 𝑘 1 -th and 𝑘 2 -th branches of the Argand diagram of Fig.3.2. Based on the above definitions, in the absence of feedback, stability is lost through the 0-1 mode for the pinned boundary condition of Fig.3.2(a) and through the 1st mode for the cantilevered boundary condition of Fig.3.2(b). Remark 4 The coupled mode of flutter is the result of two loci intersecting on the imaginary axis of the Argand diagram. Alternatively, when two loci approach each other in the complex plane but 35 do not intersect before moving away, they exhibit the phenomenon of veering [110, 111]. It should be mentioned that both positive and negative values of the feedback gain 𝑐 can cause the system to lose stability. This means that for each of the three modes of feedback and two modes of actuation, there are two cases to be considered, namely 𝑐 > 0 and 𝑐 < 0. This results in twelve cases of which we present a subset of six illustrative cases which are categorized in Table 3.1. It should be mentioned that some of the cases not presented here exhibit divergence mode of instability over a large region of the 𝛾-𝑢 e domain; these include the two cases: moment ∝ slope with 𝑐 < 0, and moment ∝ displacement with 𝑐 < 0. Table 3.1: Six specific cases chosen for simulation. Case 1 2 3 4 5 6 Actuation Moment Angle Moment Angle Moment Angle Feedback Curvature Curvature Slope Slope Displacement Displacement Sign of 𝑐 𝑐<0 𝑐<0 𝑐>0 𝑐<0 𝑐>0 𝑐<0 We complete this Section by providing the first eleven (ten) natural frequencies of the beam for the pinned (cantilevered) boundary conditions in Table 3.2. These values will be useful when we present our results on the mode of flutter instability and the frequency band in which the system loses stability in the next few subsections. The Argand diagrams in Fig.3.2 indicate that, in the absence of feedback, external flow results in 0-1 mode of flutter instability in the first frequency band for the pinned boundary condition, and 1st mode of flutter instability in the fourth frequency band for the cantilevered boundary condition. Table 3.2: Natural frequencies of beam for pinned and cantilevered boundary conditions. Ω0 Ω1 Ω2 Ω3 Ω4 Ω5 Ω6 Ω7 Ω8 Ω9 Ω10 Pinned 0 15.4 50.0 104 178 272 386 519 672 844 1037 Cantilevered - 3.52 22.0 61.7 121 200 299 417 555 713 891 36 3.4 Results of Feedback-Induced Instability Case 1: Moment ∝ Curvature with Negative Feedback Gain For this case, the moment is proportional to the curvature with feedback gain 𝑐 < 0. For 𝑐 = 0, the moment is zero, which signifies the pinned boundary condition. Therefore, the values of the natural frequencies Ω 𝑘 , 𝑘 = 0, 1, 2, · · · , 10, are those of the pinned beam in Table 3.2. The critical frequency surface is shown in Fig.3.3(a). The different colors correspond to the frequencies in the color bar shown to the left and the lines demarcate the frequency bands Π 𝑗 , 𝑗 = 1, 2, · · · , 10. Figure 3.3(a) shows well-defined striations of constant frequency band for any given 𝛾 indicating that the critical frequency Ωc𝑟 is highly influenced by the value of 𝛾. It is evident that for a constant 𝑢 e , a change in the location of sensing (value of 𝛾) can result in discontinuous changes in Ωc𝑟 . In contrast, for a constant 𝛾, Ωc𝑟 changes gradually as 𝑢 e changes. In some instances, this results in a change in the frequency band of Ωc𝑟 , such as at 𝑢 e ≈ 6 and 𝛾 ≈ 0.8. The critical stability surface is shown in Fig.3.3(b). It exhibits crests and troughs corresponding to the striations in the critical frequency surface in Fig.3.3(a), with the crests corresponding to the boundaries of the striations. The 𝑐 c𝑟 values show an increasing trend as 𝑢 e increases and drop to 0.9 Ω10 3 Π10 Ω9 Π9 2 Ω8 0.7 Π8 | 𝑐 c𝑟 | Ω7 1 Π7 Ω6 Π6 𝛾 0 Ω5 9 Π5 Ω4 Π4 Ω3 0.3 6 Π3 Ω2 𝑢e Π2 3 Ω1 0.1 Π1 0.3 0 0.1 0.1 0.7 𝛾 0.1 3 𝑢e 6 9 0.9 (a) (b) Figure 3.3: Case 1: (a) Critical frequency surface (b) Critical stability surface. The colorbar pertains only to the critical frequency surface in (a) with the lines demarcating the frequency bands. To better illustrate the topography of the critical stability surface in (b), a suitable perspective view is provided with a color gradient. 37 0.9 Ω10 Π10 0.3 Ω9 Π9 0.2 Ω8 0.7 Π8 | 𝑐 c𝑟 | Ω7 0.1 Π7 Ω6 𝛾 Π6 0 Ω5 Π5 18 Ω4 Π4 Ω3 0.3 12 Π3 Π2 Ω2 𝑢e 6 Ω1 0.1 Π1 0.3 0 0.1 0.1 𝛾 0.7 0.1 6 12 18 0.9 𝑢e (a) (b) Figure 3.4: Case 2: (a) Critical frequency surface (b) Critical stability surface. The colorbar pertains only to the critical frequency surface in (a) with the lines demarcating the frequency bands. To better illustrate the topography of the critical stability surface in (b), a suitable perspective view is provided with a color gradient. zero for 𝑢 ≥ 𝑢 e,𝑐𝑟 ≈ 9, which corresponds to the point where the pinned beam becomes unstable solely due to the external flow - see Fig.3.2(a). The 𝑐 c𝑟 values are the lowest for low 𝑢 e and high 𝛾; they increase moderately with decreasing 𝛾 and more dramatically with increasing 𝑢 e . As 𝛾 approaches the pinned end of the beam from 0.3, the magnitude of 𝑐 c𝑟 increases dramatically, once near 𝛾 = 0.3 and again near 𝛾 = 0.1. As 𝛾 decreases from 0.3 to 0.1, it is noteworthy that the value of Ωc𝑟 in Fig.3.3(a) passes through increasingly higher frequency bands for the entire range of 𝑢 e . Case 2: Angle ∝ Curvature with Negative Feedback Gain For this case, the angle is proportional to the curvature with feedback gain 𝑐 < 0. For 𝑐 = 0, the angle is zero, which signifies the cantilevered boundary condition. Therefore, the values of the natural frequencies Ω 𝑘 , 𝑘 = 1, 2, · · · , 10, are those of the cantilevered beam in Table 3.2. The critical frequency surface is shown in Fig.3.4(a). Similar to Case 1, Ωc𝑟 shows a strong dependence on 𝛾. However, the striations are wider in general and do not extend for the full range of 𝑢 e . Compared to Case 1, Ωc𝑟 shows a greater dependence on 𝑢 e and the frequency band of Ωc𝑟 changes with 𝑢 e for all values of 𝛾 > 0.25. Indeed, for 𝛾 ≈ 0.8, Ωc𝑟 smoothly passes through three different frequency bands before exhibiting a sharp drop in frequency at 𝑢 e ≈ 16. Unlike Case 1, 38 30 𝑢 ∗e = 8.93 𝑢 ∗e = 8.87 20 𝛾 = 0.69 𝑢 ∗e = 10.6 𝑢 ∗e = 10.6 𝑘 =3 Im(Ω) Im(Ω) 𝑘 =8 10 𝑢 ∗e = 8.87 Ωc𝑟 = 29.4 10 𝛾 = 0.70 𝑘 =9 𝑐 c𝑟 = 0.21 𝑘 =1 𝛾 = 0.69 𝑢 ∗e = 8.93 Ω3 𝛾 = 0.70 𝑐 c𝑟 = 0.084 0 0 0 Ω1 60 520 700 Re(Ω) 𝑐 c𝑟 = 0.15 Re(Ω) (a) (b) Figure 3.5: Argand diagrams for the beam with feedback for Case 2: (a) a small change in 𝑢 e results in a change in the mode of flutter instability without significant change in Ωc𝑟 , and (b) a small change in 𝛾 results in a change in both the mode of flutter instability and Ωc𝑟 . for 𝛾 < 0.58, several of the striations are terminated “prematurely" as the system loses stability in low frequency bands for high values of 𝑢 e . The large swath of frequency band Π3 originating from 𝛾 ∈ [0.18, 0.24] at 𝑢 e = 0.1 is not associated with a single mode of flutter. For example, as can be seen from the Argand diagram in Fig.3.5(a), for 𝛾 = 0.2, the mode of flutter changes from 𝑘 = 3 for 𝑢 e = 8.87 to 𝑘 = 1 for 𝑢 e = 8.93 due to veering [110, 111]; neither Ωc𝑟 nor 𝑐 c𝑟 change significantly. The critical stability surface is shown in Fig.3.4(b); it shows crests and troughs along constant 𝛾 similar to the critical stability surface of Case 1 shown in Fig.3.3(b). Also, the value of 𝑐 c𝑟 becomes zero for 𝑢 e ≥ 𝑢 e,𝑐𝑟 ≈ 18.1, which corresponds to the point where the cantilevered beam becomes unstable solely due to the external flow - see Fig.3.2(b). Unlike Case 1, there are additional peaks of 𝑐 c𝑟 in the interior of the domain. The regions around these peaks are associated with changes in the mode and frequencies of flutter instability. For example, for 𝑢 e = 10.6, the mode of flutter changes from 𝑘 = 9 for 𝛾 = 0.69 with 𝑐 c𝑟 = 0.084 to 𝑘 = 8 for 𝛾 = 0.7 with 𝑐 c𝑟 = 0.15; the jump in the Ωc𝑟 and 𝑐 c𝑟 values can be seen from the Argand diagram in Fig.3.5(b). Remark 5 The magnitude of 𝑐 c𝑟 for Case 2 is O(0.1); this is an order of magnitude lower than that in Case 1. 39 1.0 Ω10 40 Π10 Ω9 Π9 Ω8 | 𝑐 c𝑟 | 20 Π8 0.7 Ω7 Π7 Ω6 Π6 𝛾 0 Ω5 9 Π5 Ω4 0.4 Π4 Ω3 6 Π3 Ω2 𝑢e Π2 3 0.1 Ω1 Π1 0.4 0 0.1 0.1 0.7 0.1 3 𝑢e 6 9 1.0 𝛾 (a) (b) Figure 3.6: Case 3: (a) Critical frequency surface (b) Critical stability surface. The colorbar pertains only to the critical frequency surface in (a) with the lines demarcating the frequency bands. To better illustrate the topography of the critical stability surface in (b), a suitable perspective view is provided with a color gradient. Case 3: Moment ∝ Slope with Positive Feedback Gain For this case, the moment is proportional to the slope with feedback gain 𝑐 > 0. The values of the natural frequencies Ω 𝑘 , 𝑘 = 0, 1, 2, · · · , 10, are those of the pinned beam in Table 3.2. The critical frequency surface is shown in Fig.3.6(a). As the first case to present slope-based feedback, Fig.3.6(a) shows significantly different behavior from the two previous cases. This is in large part due to the behavior of the eigenfrequency originating at Ω0 : this locus is restricted to the imaginary axis for Case 1 involving curvature feedback; for slope-based feedback, it breaks away from the imaginary axis into the complex plane and results in coupled-mode flutter. While Fig.3.6(a) shows some distinct striations of frequency bands along constant 𝛾, these are confined to low values of 𝛾. In this limited range of 𝛾, both Ωc𝑟 and 𝑐 c𝑟 decrease as 𝛾 increases. Contrary to the previous two cases, this trend continues as 𝛾 increases towards the free end of the beam (𝛾 = 1) resulting in a very large region of Π1 . Except for small regions where 𝛾 ≈ 0.4 and 𝛾 ≥ 0.95, stability is lost through flutter in low frequency bands for a large fraction of the 𝛾-𝑢 e domain. In the Π1 region, the system loses stability in the 1st mode for low values of 𝑢 e , in the 0-th mode for intermediate values of 𝑢 e , and in the 0-1 coupled mode for high values of 𝑢 e . Stability is 40 lost in the 0-th mode when a pair of loci originating at Ω0 intersect on the imaginary axis, break away, and cross the real axis at a non-zero frequency. On the other hand, stability is lost in the 0-1 coupled mode when the loci originating at Ω0 and Ω1 intersect on the imaginary axis, break away, and cross the real axis. This behavior is similar to what was observed in Fig.3.2(a), albeit for higher values of 𝑢 e in the absence of feedback. The critical stability surface is shown in Fig.3.6(b). Similar to the previous two cases, the value of 𝑐 c𝑟 is low for high values of 𝛾 and low values of 𝑢 e ; the value of 𝑐 c𝑟 increases with decrease in 𝛾 and increase in 𝑢 e . However, contrary to the previous two cases, the critical stability surface is smooth and contains region in the interior of the domain where 𝑐 c𝑟 ≈ 0. In particular, for 𝛾 ∈ [0.6, 1.0] and 𝑢 e ≈ 6.3, the system is marginally stable in the absence of feedback and loses stability with a negligible value of feedback gain 𝑐 c𝑟 ; the corresponding critical frequency Ωc𝑟 is also small. This region of low Ωc𝑟 and 𝑐 c𝑟 will be discussed further in Section 3.4 with the help of an Argand diagram. Remark 6 The magnitude of 𝑐 c𝑟 for Case 3 is O(10); this is an order of magnitude higher than that in Case 1 and two orders of magnitude higher than that in Case 2. Case 4: Angle ∝ Slope with Negative Feedback Gain For this case, the angle is proportional to the slope with feedback gain 𝑐 > 0. The values of the natural frequencies Ω 𝑘 , 𝑘 = 1, 2, · · · , 10, are those of the cantilevered beam in Table 3.2. The critical frequency and critical stability surfaces are shown in Figs.3.7(a) and (b). The striations over the critical frequency surface and crests and troughs over the critical stability surface are quite similar to those observed in Cases 1 and 2. Also similar to these cases, the value of 𝑐 c𝑟 jumps at 𝛾 ≈ 0.3 - see Fig.3.7(b). As with all cases discussed so far, Ωc𝑟 increases monotonically as 𝛾 decreases below 0.3 - see Fig.3.7(a). Similar to Case 2, there exists a range of 𝛾, 𝛾 ∈ [0.3, 1.0], for which Ωc𝑟 drops abruptly at a specific 𝑢 e , 𝑢 e < 𝑢 c𝑟 ; this is exhibited by a sudden change of the frequency bands from high to low in Fig.3.7(a). This sudden drop in Ωc𝑟 is accompanied by a 41 1.0 Ω10 Π10 2 Ω9 Π9 Ω8 Π8 0.7 | 𝑐 c𝑟 | 1 Ω7 Π7 Ω6 Π6 𝛾 0 Ω5 Π5 18 Ω4 0.4 Π4 Ω3 12 Π3 Ω2 𝑢e Π2 6 Ω1 0.1 Π1 0.4 0 0.1 0.1 0.7 0.1 6 12 18 1 𝛾 𝑢e (a) (b) Figure 3.7: Case 4: (a) Critical frequency surface (b) Critical stability surface. The colorbar pertains only to the critical frequency surface in (a) with the lines demarcating the frequency bands. To better illustrate the topography of the critical stability surface in (b), a suitable perspective view is provided with a color gradient. sharp transition in the slope of 𝑐 c𝑟 , from rapidly increasing with 𝑢 e to rapidly decreasing with 𝑢 e - see Fig.3.7(b). The small “island" of high-frequency band that appears at 𝑢 e ≈ 9 and 𝛾 ≈ 0.85 in Fig.3.7(a) is associated with the 𝑘 = 9 mode of flutter. At 𝑢 e ≈ 9, the 𝑘 = 1 locus is sufficiently far from the real axis in the Argand diagram in Fig.3.2(b) such that introduction of feedback causes the 𝑘 = 9 locus4 to reach the real axis prior to the 𝑘 = 1 locus because the 𝑘 = 9 locus is very sensitive to 𝑐 for 𝛾 values close to 0.85. When the 𝑘 = 1 locus starts sufficiently close to the real axis, introduction of feedback causes the 𝑘 = 1 mode to reach the real axis before the other modes - this explains the “sea" of low-frequency bands in a large fraction of the upper-right domain. The border of this low frequency region shows the effect of the interplay between 𝛾 and 𝑢 e on the stability characteristics of the system. In this low-frequency region, the magnitude of 𝑐 c𝑟 exhibits an undulatory behavior as 𝑢 e increases; this can be attributed to the oscillatory behavior of the 𝑘 = 1 locus in the Argand diagram of Fig.3.2(b). Remark 7 The magnitude of 𝑐 c𝑟 for Case 4 is O(1), the same magnitude as Case 1 but an order 4The 𝑘 = 9 locus is not shown in the Argand diagram in Fig.3.2(b) but shown in Fig.3.5(b) for 𝑢 e up to at least 10.6. 42 of magnitude higher than that shown in Case 2 and an order of magnitude lower than that shown in Case 3. Case 5: Moment ∝ Displacement with Positive Feedback Gain For this case, the moment is proportional to the displacement with feedback gain 𝑐 < 0. The values of the natural frequencies Ω 𝑘 , 𝑘 = 0, 1, 2, · · · , 10, are those of the pinned beam in Table 3.2. The critical frequency and critical stability surfaces are shown in Figs.3.8(a) and (b). Although based on displacement feedback, these plots have many similarities with those of Case 3 in Figs.3.6(a) and (b), which are based on slope feedback. In particular, both the critical frequency and stability surfaces show a nearly monotonic increase with decreasing 𝛾, with ridges on the critical stability surface demarcating the change in mode of flutter. Also, the frequency bands curve in the direction of increasing 𝛾 as 𝑢 e increases; this is distinctly different from the other cases where the narrow frequency bands or striations are largely independent of 𝑢 e . Similar to Case 3, there exists a large region of the 𝛾-𝑢 e domain where the system loses stability in the 1st mode for low values of 𝑢 e , in the 0-th mode for intermediate values of 𝑢 e , and in the 0-1 coupled mode for high values of 𝑢 e . For intermediate values of 𝑢 e where stability is lost in the 0-th mode, the behavior of the system is however distinctly different from that of Case 3. For slope-based feedback (Case 3), the locus originating at Ω0 curves towards the real axis immediately after breaking away from the imaginary axis; this results in low values of Ωc𝑟 and 𝑐 c𝑟 - see Fig.3.9 (a). For displacement-based feedback (this case), the locus moves away from the real axis and converges on it at a higher value of Ωc𝑟 ; the associated value of 𝑐 c𝑟 is also higher - see Fig.3.9 (b). Remark 8 The magnitude of 𝑐 c𝑟 for Case 5 is O(100). In comparison to the other two cases of moment actuation, the magnitude of 𝑐 c𝑟 is an order of magnitude higher than Case 3 and two orders of magnitude higher than Case 1. 43 0.9 Ω10 200 Π10 Ω9 Π9 Ω8 Π8 100 Ω7 Π7 Ω6 Π6 𝛾 0 Ω5 9 Π5 Ω4 Π4 6 Ω3 Π3 𝑢e Ω2 Π2 3 Ω1 0.3 Π1 0.3 0.6 0 0.1 𝛾 0.1 3 𝑢e 6 9 0.9 (a) (b) Figure 3.8: Case 5: (a) Critical frequency surface (b) Critical stability surface. The colorbar pertains only to the critical frequency surface in (a) with the lines demarcating the frequency bands. To better illustrate the topography of the critical stability surface in (b), a suitable perspective view is provided with a color gradient. Case 6: Angle ∝ Displacement with Negative Feedback Gain For this case, the angle is proportional to the displacement with feedback gain 𝑐 > 0. The values of the natural frequencies Ω 𝑘 , 𝑘 = 1, 2, · · · , 10, are those of the cantilevered beam in Table 3.2. The critical frequency and critical stability surfaces are shown in Figs.3.10(a) and (b); they resemble 15 𝑢 ∗e = 5 15 𝑢 ∗e = 5 Im(Ω) Im(Ω) 𝑘 =1 𝑘 =1 Ωc𝑟 = 23.20 𝑢 ∗e = 5 𝑢 ∗e = 5 𝑐 c𝑟 = 29.16 𝑘 =0 Ωc𝑟 = 3.94 𝑘 =0 𝑐 c𝑟 = 2.28 0 0 0 Ω0 10 Ω1 20 0 Ω0 10 Ω1 20 Re(Ω) Re(Ω) (a) (b) Figure 3.9: Argand diagrams for the beam with slope and displacement feedback. For 𝑢 ∗e = 5.0 and 𝛾 = 0.85, (a) moment actuation based on slope feedback (Case 3) results in Ωc𝑟 = 3.94 and 𝑐 c𝑟 = 2.28; (b) moment actuation based on displacement feedback (Case 5) results in Ωc𝑟 = 23.30 and 𝑐 c𝑟 = 29.16. 44 0.9 Ω10 Π10 20 Ω9 Π9 Ω8 Π8 | 𝑐 c𝑟 | 10 Ω7 Π7 Ω6 𝛾 Π6 0 Ω5 Π5 18 Ω4 Π4 Ω3 12 Π3 Π2 Ω2 𝑢e 6 0.3 Ω1 Π1 0 0.3 0.6 0.1 𝛾 0.1 6 12 18 0.9 𝑢e (a) (b) Figure 3.10: Case 6: (a) Critical frequency surface (b) Critical stability surface. The colorbar pertains only to the critical frequency surface in (a) with the lines demarcating the frequency bands. To better illustrate the topography of the critical stability surface in (b), a suitable perspective view is provided with a color gradient. those of Case 3 in Figs.3.6(a) and (b) and those of Case 5 in Figs.3.8(a) and (b) although the mode of actuation and sign of 𝑐 are different for both Cases 3 and 5 and the mode of sensing is different for Case 3. The tendency of the critical frequency to increase with increasing 𝑢 e and decreasing 𝛾 is particularly noticeable due to the “curving up" of the frequency bands in Fig.3.10(a); this behavior is present in all cases but is very distinct in Cases 3 and 5. Similar to Case 3, for low 𝑢 e , Fig3.10(a) shows small regions where stability is lost through flutter in high frequency bands. Otherwise, both the critical frequency and stability surfaces resemble Figs.3.8(a) and (b) pertaining to Case 5, increasing towards low 𝛾 and high 𝑢 e with ridges on the critical stability surface demarcating the changes in the mode of flutter. Remark 9 The magnitude of 𝑐 c𝑟 for Case 6 is O(10). A comparison of all six cases indicate that the magnitude of 𝑐 c𝑟 depends on the modes of actuation and sensing. A change in the mode of actuation from angle to moment increases 𝑐 c𝑟 by one order of magnitude on average. Similarly, changing the mode of sensing from curvature to angle as well as from angle to displacement increases 𝑐 c𝑟 by one order of magnitude on average - see Table 3.3. 45 Table 3.3: Order of magnitude for critical stability. Curvature Slope Displacement Sensing Sensing Sensing Moment Actuation 1.0 10 100 Angle Actuation 0.1 1.0 10 3.5 Application to Underwater Propulsion An Underwater Vehicle with a Flexible Propulsor All of the analysis presented in Sections 3.1 and 3.2 and the results presented in Section 3.4 can find a potential application in the propulsion of undersea vehicles. To motivate this, we consider a submersible comprised of a rigid body and a tail-like flexible beam, immersed in a quiescent fluid; the rigid body is connected to the flexible beam by an active revolute joint - see Fig.3.11. For the sake of simplicity, we assume that the drag of the submersible is entirely due to the rigid body and the thrust is produced entirely by the flexible tail. The submersible is assumed to move with constant velocity 𝑈e in a state of dynamic equilibrium, where thrust and drag forces are equal and opposite. The rigid body is assumed to have negligible rotational motion due to its inertia and consequently it translates with constant velocity 𝑈e . This is in conformity with the assumptions made in Section 3.1, namely, the pinned joint is at the origin of an inertial reference frame, and the beam is immersed in a fluid that moves with constant relative velocity 𝑈e . 𝑦 active joint flexible beam rigid body 𝑥 𝑈e 𝐿 Figure 3.11: A rigid body connected to a tail-like flexible beam by an active revolute joint. 46 Propulsive Characteristics Thrust, Power, and Efficiency For a “slender fish", Lighthill [29] estimated the thrust, power, and efficiency assuming that the fish has neither mass nor area at its leading edge. For the more general case, Hellum [52] adapted these expressions for computing the nondimensional thrust F and power P Z ("    # "    # ) Ω 2𝜋/Ω 𝜕𝑣 2 𝜕𝑣 2 𝜕𝑣 2 𝜕𝑣 2 F= − 𝑢e − − 𝑢e 𝑑𝜏 (3.17a) 4𝜋 0 𝜕𝜏 𝜕𝑢 𝜕𝜏 𝜕𝑢 𝑢=1 𝑢=0 Z 2𝜋/Ω         Ω 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑣 P= + 𝑢e − + 𝑢e 𝑑𝜏 (3.17b) 2𝜋 0 𝜕𝜏 𝜕𝜏 𝜕𝑢 𝑢=1 𝜕𝜏 𝜕𝜏 𝜕𝑢 𝑢=0 where Ω = Ωc𝑟 is the non-dimensional frequency of oscillation, and the thrust and power expres- sions are approximated by their averages over one period of oscillation. The function 𝑣(𝑢, 𝜏) is solved using the procedure outlined in Section 3.2 and has both real and imaginary parts. Since only the real part physically contributes to the thrust and power, Re[𝑣] is used in place of 𝑣 in (3.17) and has the form X4 Re[𝑣(𝑢, 𝜏)] = eRe[𝑧𝑛 ]𝑢 {Re [ 𝐴𝑛 ] cos (Im [𝑧 𝑛 ] 𝑢 + Ω𝜏) − Im [ 𝐴𝑛 ] sin (Im [𝑧 𝑛 ] 𝑢 + Ω𝜏)} 𝑛=1 (3.18) From the thrust and power values, the Froude efficiency [29] can be calculated as F𝑢 e 𝜂= (3.19) P Both the thrust and power expressions depend on the amplitude of the waveform of the flexible tail, which depends on the 𝐴𝑛 terms in (3.18). The 𝐴𝑛 terms are obtained from the null-space of Z in (3.13); therefore, the thrust and power will depend on the scaling of the null-space vector. This scaling is arbitrary because we are using a linear model of the system. Ideally, the amplitude of the waveform would be determined by a limit cycle analysis of the nonlinear model, which is outside 47 the scope of this work. Therefore, we focus on the waveform efficiency, which is not dependent on the amplitude of the waveform. Wave Speed and Phase Smoothness The efficiency of the tail-like propulsor, given by (3.19), assumes that the thrust generated is positive. This can be verified from the sign of F, computed using (3.17a) with an arbitrary amplitude of the waveform. For a tail oscillating with a waveform 𝑔(𝑢, 𝜏) = ℎ(𝑢) cos(Ω𝜏 − 𝑘𝑢) where 𝑘 is the non-dimensional wavenumber, Hellum [52] provided a simple condition for the tail to generate positive thrust. This condition, which was adapted from Lighthill [29], is given as Ω 𝑢 𝑒 𝛽−1/2 > 𝑢 𝑒 𝛽−1/2 ⇒ <1 (3.20) 𝑘 (Ω/𝑘) where (Ω/𝑘) is referred to as the non-dimensional phase velocity. When this condition is met, the efficiency can be alternately computed using the expression " # 1 𝑢 𝛽 −1/2 𝑒 𝜂∗ = 1 − 1− (3.21) 2 (Ω/𝑘) which has been adapted from [112] using non-dimensional variables. It was shown in [112] that 𝜂∗ is restricted to lie in the range [0.5, 1.0]. The motion of the flexible propulsor is comprised of four traveling waves - see (3.18); therefore the condition in (3.20) and the expression in (3.21) are inapplicable. We can however compute a value of the average non-dimensional wavenumber, which we denote by 𝑘; ¯ simulation results show that positive thrust is generated when (3.20) is satisfied with 𝑘 replaced by 𝑘. ¯ To compute 𝑘, ¯ we first recognize that 𝑣(𝑢, 𝜏) of (3.7) is a complex helix defined by the shape function 𝑓 (𝑢) which rotates with angular velocity Ω. At any given time, 𝑘¯ can be computed from the phase of the helix 𝜙(𝑢) as follows 48 𝑘¯ = 𝜙(0) − 𝜙(1), 𝜙(𝑢) = atan2(Im[ 𝑓 (𝑢)], Re[ 𝑓 (𝑢)]) (3.22) For the purpose of illustration, we plot 𝜙(𝑢) and 𝑘¯ for the specific operating point: 𝑢 𝑒 = 3.4 and 𝛾 = 0.90 for Case 6 - see Fig.3.12. It can be seen that 𝜙(𝑢) decreases as 𝑢 varies from 0 to 1, which signifies that the waveform travels from the hinged end to the free end of the beam. The value of 𝑘¯ is the negative of the slope of the straight line joining 𝜙(0) and 𝜙(1). The phase angle 𝜙 in Fig.3.12 is observed to exhibit undulations about the straight line, implying that the phase does not vary linearly. To characterize this variation, we plot the spatial derivative of the phase angle (𝑑𝜙/𝑑𝑢) in Fig.3.12 and define the phase smoothness factor (PSF): min|(𝑑𝜙/𝑑𝑢)| P𝑆𝐹 = (3.23) max|(𝑑𝜙/𝑑𝑢)| The value of PSF signifies the extent to which the traveling wave exhibit a stop-and-go motion as it moves from the hinged end (𝑢 = 0) to the free end of the flexible tail (𝑢 = 1). When 𝜙(𝑢) varies linearly, the value of PSF is equal to 1, which describes a wave traveling with constant velocity. In 0 -0.02 -1 -0.03 (𝑑𝜙/𝑑𝑢) 𝜙(𝑢) 𝜙(𝑢) -2 -0.04 -3 (𝑑𝜙/𝑑𝑢) slope:- 𝑘¯ -0.05 0.0 0.2 0.4 0.6 0.8 1.0 𝑢 Figure 3.12: Phase plot for Case 6 at 𝑢 𝑒 = 3.4 and 𝛾 = 0.9. 49 the next section, it will be shown that the Froude efficiency is highly correlated with the value of PSF, with higher efficiencies associated with PSF values closer to unity. Illustrative Examples of Traveling Waveforms Effect of Change in 𝑢 e and 𝛾 We illustrate the change in the traveling waveform and its propulsive characteristics due to changes in the external flow velocity 𝑢 e and the location of sensing 𝛾. We consider three operating points from Case 6, presented in Section 3.4: a nominal point (Point 1) and two other points, obtained by varying either 𝛾 (Point 2) or 𝑢 e (Point 3). The critical feedback gain, the critical frequency, the Froude efficiency, the average nondimensional wavenumber, the value of PSF, the wavespeed, and the value of efficiency computed using (3.21) at these three points are shown in Table 3.4. For the nominal point (Point 1: 𝑢 e = 3.2, 𝛾 = 0.30), the waveform is shown in Fig.3.13 (a). The waveform exhibits a “stop-and-go" motion with amplitudes varying considerably along the length of the tail, which corresponds to its low PSF value. The value of 𝜂∗ lies near the lower bound of the range [0.5, 1.0]; this follows from (3.21) since the wavespeed (Ω/ 𝑘) ¯ is significantly higher than the external flow velocity 𝑢 e 5. The value of 𝜂 computed using (3.19) is also near the lower bound of [0.5, 1.0]. When the location of sensing alone is changed (Point 2: 𝑢 e = 3.2, 𝛾 = 0.90), we notice a significant drop in the frequency with a less significant drop in the wavenumber, which results in a lower wavespeed - see Table 3.4. While still producing positive thrust as per (3.20), the drop in wavespeed results in a jump in the value of 𝜂∗ , which closely matches the computed value of 𝜂. 5For our simulations, we assumed 𝛽 = 0.9822 - see Section 3.3. This gives a value of 𝛽−1/2 = 1.0090 ≈ 1. Table 3.4: Propulsive characteristics at three operating points of Case 6. Point 𝑢e 𝛾 𝑐 c𝑟 Ω = Ωc𝑟 𝜂 𝑘¯ PSF ¯ (Ω/ 𝑘) 𝜂∗ 1 3.2 0.30 9.78 193.1 0.514 13.40 0.074 14.41 0.611 2 3.2 0.90 3.65 24.0 0.701 3.24 0.403 7.40 0.716 3 12.4 0.30 11.25 258.1 0.751 11.11 0.505 23.23 0.767 50 The waveform, shown in Fig.3.13 (b), indicates a fewer number of undulations compared to that of Fig.3.13 (a)in accordance with the the lower wavenumber. Additionally, it can be seen that the waveform travels much more smoothly along the length of the tail which is captured by its much higher value of PSF. When the external flow velocity alone is changed (Point 3: 𝑢 e = 12.4, 𝛾 = 0.30), we notice a moderate increase in the frequency but a slight decrease in the wavenumber, which results in a higher wavespeed - see Table 3.4. Though the wavespeed is increased, the system operates at a much higher external velocity which results in a jump in the value of 𝜂∗ , which, again, closely matches the computed value of 𝜂. The waveform, shown in Fig.3.13 (c), indicates a similar number of undulations compared to that of Fig.3.13 (a) due to the similar magnitude of the wavenumber. While there are many undulations, it can be seen that the waveform travels much more smoothly along the length of the tail, which is captured by its much higher value of PSF. Remark 10 The results in Table 3.4 indicate that there is a strong positive correlation between the value of PSF and the values of 𝜂 and 𝜂∗ . This trend has been observed for other points in the 𝛾-𝑢 e domain for Case 6, as well as other cases with different modes of sensing and actuation. Dependence of Propulsive Characteristics on 𝑢 e and Ω We select three operating points with identical values of 𝑢 e and similar values of Ω6 from three different cases, namely Case 1, Case 2, and Case 4. The propulsive characteristics of these three cases are found to be very similar - see Table 3.5. These characteristics, which include the efficiency, the average nondimensional wavenumber, the value of PSF, and the wavespeed, are particularly close for Cases 1 and 2, for which the frequencies differ by only 0.4%. The characteristics of Case 4 differ slightly more as its frequency is 4.5% lower than the other two cases. It is observed that for a constant 𝑢 𝑒 , a lower frequency is associated with higher values of 𝜂 and PSF; this agrees well with the trend illustrated by Points 1 and 2 in Section 3.5. The waveforms for the three cases 6For a fixed value of 𝑢 e , each critical frequency plot degenerates to a line where each point on the line corresponds to a different value of Ωc𝑟 . Through trial and error it is possible to find similar Ωc𝑟 values across multiple cases. 51 0 Ω𝜏 𝜋 2𝜋 0 𝑢 1 0 𝑢 1 0 𝑢 1 (a) (b) (c) Figure 3.13: Traveling waveforms over one complete cycle, shown at intervals of 𝜋/4, for the three operating points of Case 6 shown in Table 3.4: (a) Point 1, (b) Point 2, and (c) Point 3. are shown in Fig. 3.14. The waveforms for Cases 1 and 2 are nearly indistinguishable, while the waveform for Case 4 is slightly different from the other two, which can be attributed to the slightly different value of the frequency of oscillation. Remark 11 The data in Table 3.5 are a small set of results which indicate that the nature of a waveform and its propulsive characteristics depend solely on the values of 𝑢 e and Ω, and are Ω𝜏 = 0 0 𝑢 1 Figure 3.14: Traveling waveforms for the three operating points of Case 1 (dotted), Case 2 (dashed), and Case 4 (solid), shown in Table 3.5. The waveforms from Case 1 and 2 exhibit tremendous overlap while the waveform from Case 4 is slightly separated. 52 independent of the modes of actuation and sensing. Table 3.5: Propulsive characteristics at three operating points of Cases 1, 2, 4 with identical 𝑢 𝑒 and similar Ω values. Case 𝑢e 𝛾 𝑐 c𝑟 Ω = Ωc𝑟 𝜂 𝑘¯ PSF ¯ (Ω/ 𝑘) 𝜂∗ 1 6.9 0.30 1.46 152.8 0.720 9.51 0.318 16.073 0.715 2 6.9 0.11 0.11 152.2 0.723 9.47 0.318 16.074 0.715 4 6.9 0.92 1.24 145.7 0.746 9.06 0.324 16.079 0.715 The results observed in this section, that the waveforms are uniquely determined by the values of 𝑢 e and Ωc𝑟 and are independent of the modes of sensing and actuation, can be explained mathematically as follows. For convenience, we revisit (3.16), (3.11), and (3.13) from Sections 3.2 and 3.3: X4 𝑣(𝑢, 𝜏) = 𝐴𝑛 𝑒 R𝑒[𝑧𝑛 ]𝑢 𝑒𝑖 {I𝑚[𝑧𝑛 ]𝑢+R𝑒[Ωc𝑟 ]𝜏 } (3.16 revisited) 𝑛=1 p 𝑧 4 + 𝑢 2e 𝑧 2 + 2𝑢 e 𝛽 𝑖Ω𝑧 − Ω2 = 0 (3.11 revisited)       1 1 1 1   𝐴1  0           𝑧 2 𝑒 𝑧1 𝑧 2 𝑒 𝑧2 𝑧 2 𝑒 𝑧3 𝑧 2 𝑒 𝑧4   𝐴2  0  1 2 3 4       =  (3.13 revisited) 𝑧 3 𝑒 𝑧1 𝑧 3 𝑒 𝑧2 𝑧 3 𝑒 𝑧3 𝑧 3 𝑒 𝑧4   𝐴     1 2 3 4   3  0            𝛿1 𝛿2 𝛿3 𝛿4   𝐴4  0      | {z } Z The waveform in (3.16) is defined by three terms: Ωc𝑟 , 𝑧 𝑛 , and 𝐴𝑛 , of which Ωc𝑟 is the same for the different cases. From (3.11), we can see that for a fixed value of 𝛽, the solutions of 𝑧 𝑛 , 𝑛 = 1, 2, 3, 4, are uniquely defined by 𝑢 e and Ωc𝑟 and are independent of the modes of sensing and feedback. As a consequence, the first three rows of Z in (3.13) are solely functions of 𝑢 e and Ωc𝑟 . It can be shown that the first three rows of Z are linearly independent though Z is singular as that is how Ωc𝑟 was solved for. Therefore, the direction of the null space vector (𝐴1 , 𝐴2 , 𝐴3 , 𝐴4 )𝑇 is independent 53 of the last row of Z and thus independent of the modes of sensing and feedback. Therefore, each term of the waveform in (3.16) is solely a function of 𝑢 e and Ωc𝑟 . 3.6 Conclusion Beam flutter occurs due to non-conservative loading. Such loading is commonly generated by a follower force or through interaction of the beam with a fluid, flowing internally or externally. It is also possible to produce non-conservative loading by applying an actuation, proportional to some state of the beam, at one of the boundaries. This work provides a generalized description of flutter, generated using this method, in a pinned-free beam. The actuation can take the form of a moment or an angle prescribed at the pinned boundary. This actuation is proportional to some state of the beam (displacement, slope, or curvature) measured at any location along its length. The onset of flutter is not symmetric about zero gain, meaning that there are twelve potential flutter mechanisms identified here: two modes of actuation, three modes of sensing, and two signs of the critical gain for each combination. This opens up a wide range of physical mechanisms, beyond the standard follower force and fluid-flow mechanisms, that can be used to produce and study flutter, some of which we intend to realize in future work. For each combination of actuation and sensing, the critical gain was determined over a range of external flow velocities and sensing locations along the beam. For a majority of the twelve possible combinations, stability was lost through flutter; a representative sample of six cases (one for each combination of actuation and sensing) were investigated. These six cases illustrated a rich set of stability transitions that depend strongly on the location of sensing and mildly on the external flow velocity. It was observed that small changes in the location of sensing could result in very different modes of flutter with large jumps in the critical frequency thereby resulting in significantly different traveling waveforms. These traveling waveforms of the flexible beam could be exploited to develop a propulsion mechanism for underwater vehicles. Because the Euler–Bernoulli beam model is fourth-order in space, the solution naturally comprises four separate traveling waves. Based on the spatial derivative of the beam phase, we constructed a metric, “the phase smoothness factor”, which is 54 a measure of how closely the four traveling waves can be approximated by a single waveform. Waveforms which are smoother demonstrate higher propulsive efficiency. The same is observed in nature, where fish swim with optimized waveforms that are smooth and efficient. Interestingly, the propulsive characteristics of the beam do not depend on the combination of actuation and sensing by which flutter is produced; they depend only on the values of the dimensionless fluid velocity and critical frequency, which completely define the waveform. 55 CHAPTER 4 EXPERIMENTAL VALIDATION OF FEEDBACK-INDUCED FLUTTER INSTABILITY OF A FLEXIBLE BEAM IN FLUID FLOW In this chapter we aim to present experimental verification of the feedback-induced flutter of a flexible beam immersed in fluid flow as introduced in Chapter 3. In order to accomplish this, we procured a beam of suitable flexibility and density and affixed strain gauges to measure the curvature at a point along its centerline. We then built a setup for applying a moment to the leading edge of the beam as well as all of the electronics needed for the curvature measurement and torque control. We designed a Simulink program to implement the feedback control scheme which was run through the ControlDesk interface. We tested the system at several flow velocities with both positive and negative feedback gain to determine the critical feedback gain and the critical frequency at which flutter occurred. To validate the experimental findings, we adapted the moment-proportional-to- curvature system from Chapter 3 to account for experimental considerations. We then compare the results of the experiments and numerical simulations which show many similarities to give confidence to the model used. The layout of this chapter is organized as follows. The experimental hardware, setup, control scheme, and procedures are described in Sections 4.1 and 4.2. The experimental results are presented in Section 4.3 while the methods and results of the numerical model are in Section 4.4. The results of both are compared in Section 4.5. 4.1 Experimental setup Equipment The flexible beam used in this experiment had width 𝑤 = 0.9 mm, height ℎ = 100 mm, and length 𝐿 = 400 mm, and was made of Cirlex® with density 𝜌b = 1420 kg/m and Young’s modulus 𝐸 = 2.7 GPa. Mounted directly opposite to each other along the centerline of either face of the beam and 100 mm from the leading edge, were two strain gauges of 6 mm length, resistance ∆𝑅/𝑅0 𝑅 = 120 Ω, and 𝑘-factor 𝑘 = 2.35 𝜖 - see Fig.4.1. Each strain gauge was bonded to the 56 surface of the beam using Gorilla® superglue, soldered to lead wires, and all exposed electrical components were covered in several layers of clear nail polish to prevent infiltration of the fluid. Figure 4.1: CAD model of experimental setup on the aluminum plate. Shown are the flexible beam with attached strain gauge which is clamped by the rotating shaft. The shaft is actuated by the motor through the coupler with the encoder measuring the rotation. From each strain gauge, the pair of lead wires travelled along the face of the beam and up the rotating shaft to form half of Wheatstone Bridge - see Fig.4.2. This configuration is standard for measuring the curvature a beam while rejecting the normal strain and temperature fluctuations of the beam and strain gauges. The voltage across the Wheatstone bridge was filtered and amplified by an Omega® DMD4059 Bridge/Strain Gage Signal Conditioner powered by a 108 V DC signal1. The excitement voltage was 𝑉𝑒𝑥𝑐𝑖𝑡𝑒 = 1 V and the amplification was 𝑆 𝑎𝑚 𝑝 = 1000 mV / 30 mV. The amplified signal was passed to a dSPACE® DS1103 PPC Controller Board (hereafter dSPACE® board) which processed it as described in Section 4.1. The resultant torque command was sent to 1Powering the DMD4059 using the built-in 120 V AC power adapter caused an overpowering 60 Hz signal to dominate the strain gauge measurement. 57 Figure 4.2: Experimental setup: (a) Omega® DMD4059 Bridge/Strain Gage Signal Conditioner; (b) Advanced Motion Controls® 12A8 PWM servo drive; (c) dSPACE® DS1103 PPC Controller Board; (d) DC power supplies; (e) Plexiglass plates; (f) Wheatstone bridge; (g) Faulhaber® 3863H024 CR graphite commutation brushed DC motor; (h) Dynapar® E2320000541 quadrature encoder; (i) Clamping shaft; (j) Strain gauge (2x); (k) Flexible beam; (m) Grasshopper3 GS3-U3-28S4M camera. an Advanced Motion Controls® 12A8 PWM servo drive. This motor driver was powered by a 24 V DC power supply which amplified the signal by 𝑀𝑎𝑚 𝑝 = 1 A/V and applied a bi-directional PWM voltage to the Faulhaber® 3863H024 CR graphite commutation brushed DC motor with torque constant 𝜏𝑚 = 0.0398 N.m/A. This motor was connected by a custom-built shaft coupler to the 6.35 mm aluminum rod which had been split down the middle to clamp the leading edge of the beam. The rotation of the shaft was measured by a Dynapar® E2320000541 quadrature encoder (2000 ticks/rev) which was read by the dSPACE® board. The motor and encoder were mounted to the upper face of a 4.76 mm × 100 mm × 400 mm aluminum plate with the flexible beam held by the aluminum shaft below. 58 Testing Facility The plate was placed in the Turbulent Mixing and Unsteady Aerodynamics Laboratories (TMUAL) Large Water Tunnel Facility shown in Fig.4.3, marginally above the water level so as not to submerge the electronics but fully submerge the beam. The plate was secured to large suspended plexiglass plates that spanned the width of the tunnel, directly on the surface of the water. The water tunnel had dimensions width = 61 cm, height = 61 cm, and length = 244 cm, and could sustain a steady flow rate of 30 cm/s that has been shown to be very well conditioned. This flow velocity was controlled using a LabviewTM application. The beam was placed as centrally as possible along the width and length of the water tunnel to reduce the boundary layer effects of the tunnel walls. Below the water tunnel was a Grasshopper3TM GS3-U3-28S4M camera from Point Grey Research® which was used to record the videos and ensure consistent alignment of the beam. Figure 4.3: CAD model of the Turbulent Mixing and Unsteady Aerodynamics Laboratories (TMUAL) Large Water Tunnel Facility with a cross-section of 61 cm × 61 cm and a usable length of 244 cm. Located below the test section of the water tunnel is the camera used to record the Trials. 59 Control Scheme and Interface Figure 4.4: Simulink Model used to control the feedback on the dSPACE® board using the dSPACE® ControlDeskTM application interface. Shown are the encoder blocks, logic, and outputs (dashed line); the curvature-torque feedback loop including the dSPACE® blocks, the feedback gain, filters, and outputs (dot-dashed line); a cut-off function for stopping the system from damaging itself due to large deformation or torque input; and a PID program to return the system to the home orientation (dotted line). To implement the feedback control scheme, a SimulinkTM model, shown in Fig.4.4, was de- signed to be run on the dSPACE® board using the dSPACE® ControlDeskTM application interface - see Fig.4.5. The amplified curvature voltage 𝑉0 was read by the dSPACE® ADC block which was converted to the curvature using (4.1), a function of the beam width (𝑤), the strain gauge 𝑘-factor (𝑘), the Wheatstone bridge configuration (the factor of 4), and the Omega® Signal Conditioner (𝑉𝑒𝑥𝑐𝑖𝑡𝑒 and 𝑆 𝑎𝑚 𝑝 ) parameters – see Section 4.1. 4 Curvature = 𝑉0 (4.1) 𝑉𝑒𝑥𝑐𝑖𝑡𝑒 𝑘 𝑤 𝑆 𝑎𝑚 𝑝 This curvature reading went through a high-pass filter (4.2) with time constant 𝜏 = 1.5915, which corresponds to corner frequency 𝑓𝑐 = 2𝜋𝜏 1 = 0.1 Hz. This value was intentionally chosen to be quite slow so as to less affect the flutter instability dynamics. 𝜏𝑠 𝑉 𝑓 (𝑠) = Curvature(𝑠) (4.2) 1 + 𝜏𝑠 60 The curvature measurement was then multiplied by the desired feedback gain factor – see Section 4.2. Following this, to account for subsequent amplification and motor parameters, the signal was divided by the motor driver amplification 𝑀𝑎𝑚 𝑝 and the motor torque constant 𝜏𝑚 . This voltage was then transmitted to the Advanced Motion Controls® 12A8 PWM servo drive with the dSPACE® DAC block. The values transmitted by the Dynapar® encoder were read in by the dSPACE® encoder block as raw tick values. These were converted to radians and radians/sec given the 2000 ticks/rev of the encoder, the 16 : 20 gear ratio of the motor sprocket to the encoder sprocket, and the 1000 Hz frequency of the dSPACE® board real-time operating system. Together, the commanded torque and the maximum angular displacement were used as triggers to prescribe the maximum conditions that the shaft could experience before the system would automatically cut off the torque so as to prevent damage. Finally, a simple PID controller block was included to return the beam to its starting orientation between Trials. To interact with the feedback controller, the Simulink model was imported into the dSPACE® ControlDeskTM application interface - see Fig.4.5. This enabled the experimenter to view the curvature, encoder, and commanded torque values as well as to increment the feedback gain. For each Trial, the application would record the data at the 1000 Hz polling rate of the real-time operating system and then export that data to MATLAB® for further analysis. Additionally, the experimenter could set the home position and command the system to return to it between Trials. 61 Figure 4.5: ControlDesk application interface used for interacting with the dSPACE® board. (a) Controls for beginning and ending the data recording. (b) Commanded torque. (c) Raw curvature (red) and high-pass filtered curvature (green). (d) Encoder position (red) and velocity (green). (e) Controls for PID loop to return to home position. (f) Feedback gain interface. (g) Commanded torque (blue) and cutoff trigger statuses (red, green, yellow). 4.2 Experimental Methods Experimental Procedure For the initial setup, the beam was aligned with the direction of flow of the water tunnel using the viewport at the back of the water tunnel. The stress was removed from the strain gauge wire lead so that it would not act as a spring by applying a restoring force to the shaft or beam about the default position. The aluminum plate was then clamped down to the plexiglass plates of the water tunnel to prevent movement of the motor relative to the initial position. The ControlDeskTM program was then started at this point so as to set this orientation to the zero point. 62 The procedure for each Trial was as follows. The ControlDeskTM program was initialized and after waiting for the high-pass filter to remove the DC component of the strain readings, the measurement and video recording were begun. The gain was then slowly increased/decreased from 𝐶 = 0 in increments of ±0.01 until the onset of flutter instability was observed. The torque control was then cancelled, either by the operator or automatically by the program, and then the beam was commanded back to the original orientation. After this, the ControlDeskTM program was re-initialized and the process was repeated for the next Trial. Six Trials at different flow velocities and positive/negative gain were conducted - see Table 4.1. An additional Trial was conducted at a higher flow velocity without any feedback gain to ascertain that the system did not develop flutter when subject to fluid flow alone. Table 4.1: Trial scheme for experiments. Trial # Flow Velocity (cm/s) Gain Sign 1 0 + 2 0 - 3 10 + 4 10 - 5 20 + 6 20 - 7 30 N/A Numerical methods To analyze the data in MATLAB® , several tools were used and are briefly described here. A bidirectional (two-pass) first-order low-pass filter2 was applied to the high-pass-filtered curvature data to make the results smoother and allow the behavior to be more easily discerned. Additionally, for the Fast-Fourier Transform (FFT) results of the low frequency flutter, the weighted average of the three points nearest to the peak was taken as the flutter frequency to account for the very low resolution at such low frequencies. Finally, the built-in MATLAB® spectrogram function was utilized over the entirety of each Trial data, specifying the window to be ≈ 2 periods of the flutter 2Using the discrete low-pass filter equation 𝑋𝑛 = (1 − 𝜎)𝑋𝑛−1 + 𝜎𝑋𝑛𝑒𝑤 , with 𝜎 = 0.25. 63 frequency, but leaving the other arguments3 as default. The results of the spectrogram were used to determine the magnitude of the frequency of maximum amplitude in each window. It is important to note that this amplitude does not correspond directly to the amplitude of oscillation seen from the data directly unless there was only a single frequency present. 4.3 Results Trial #1 Results For this Trial, the flow velocity was 𝑈𝑒 = 0 cm/s and the feedback gain was incremented by +0.01. Increasing the feedback gain to 𝐶 = 0.15 over the first 54.511 s resulted in no discernible oscillatory behavior that can be seen from the high-pass-filtered curvature data in Fig.4.6. Derived from the spectrogram of the curvature data (described in Section 4.2), the amplitude of the highest-amplitude (peak) frequency shown in Fig.4.7 for this period is on the order of ≈ 10−6 m−1 . Upon increasing the feedback gain to 𝐶 = 0.16 at 𝑡 = 54.511 s, there begins to be a noticeable increase in the amplitude of the response which continues to grow as the feedback gain is further increased. For the period of 𝐶 = 0.22 from 𝑡 = 114.838 s to 𝑡 = 166.854 s, the system sustains a very slight but consistent oscillation as seen in Fig.4.6. The oscillation during this period has an amplitude of peak-frequency on the order of ≈ 10−3 m−1 with a frequency of 𝜔𝑐𝑟 = 0.3014 Hz as derived from the FFT procedure outlined in Section 4.2. This minimal but steady oscillatory behavior is evocative of that of a limit cycle, a nonlinear phenomenon in which system loses stability through flutter, but the magnitude of oscillation is held in check by the non-linear effects of the sensing and actuation system as well as the damping effects of the surrounding fluid. Upon raising the feedback gain to 𝐶 = 0.23 at 𝑡 = 166.854 s, the system gradually began oscillating with increasingly greater amplitude (see Fig.4.6) until the cut-off trigger was activated at 𝑡 = 193.547 and stopped the input torque. Over the course of these 26.7 s, the amplitude of the peak-frequency oscillation grew by two orders of magnitude (see Fig.4.7), clearly signifying that the system had lost stability through flutter. During this period, the frequency of oscillation 3nooverlap: samples of overlap between adjoining segments; and nfft: sampling points to calculate the discrete Fourier transform [https://www.mathworks.com/help/signal/ref/spectrogram.html] 64 0.5 0.4 C = 0.11 C = 0.16 C = 0.20 C = 0.23 0.3 0.2 Curvature (m-1 ) 0.1 0 -0.1 -0.2 -0.3 -0.4 0 20 40 60 80 100 120 140 160 180 200 Time (s) Figure 4.6: Curvature measurements for the full Trial of beam in 0 cm/s flow with positive feedback gain, that has been two-pass filtered - see Section 4.2. The solid vertical lines correspond to when the feedback gain was increased by 0.01. Prior to increasing the feedback gain to 𝐶 = 0.16 at 𝑡 = 54.511 s, there is no discernible oscillatory behavior, but increasing the feedback gain to 𝐶 = 0.16 causes the system to show a small but marked oscillatory behavior with a frequency of 𝜔𝑐𝑟 = 0.3014 Hz. Upon raising the feedback gain to 𝐶 = 0.23 at 𝑡 = 166.854 s, the system gradually began oscillating with increasingly greater amplitude at a frequency of 𝜔𝑐𝑟 = 0.426 Hz. The dashed vertical line denotes where the cut-off trigger was activated at 𝑡 = 193.547 s. 𝜔𝑐𝑟 = 0.426 Hz was greater than the previous, much smaller oscillations. 65 10-1 C = 0.11 C = 0.16 C = 0.20 C = 0.23 10-2 Magnitude of peak frequency 10-3 10-4 10-5 10-6 10-7 0 20 40 60 80 100 120 140 160 180 200 Time (s) Figure 4.7: Amplitude of the peak frequency for the full Trial of beam in 0 cm/s flow with positive feedback gain taken from a spectrogram with a window of 4.0 s (approximately 2 cycles), presented on a logarithmic scale. The solid vertical lines correspond to when the feedback gain was increased by 0.01 while the dashed vertical line denotes where the cut-off trigger was activated at 𝑡 = 193.547 s. Prior to increasing the feedback gain to 𝐶 = 0.16 at 𝑡 = 54.511 s, the magnitude of the peak frequency was on the order of ≈ 10−6 m−1 , but increasing the feedback gain to 𝐶 = 0.16 and beyond causes the magnitude of the peak frequency to increase by several orders of magnitude with increasing feedback gain. This behavior continues slowly until the feedback gain is set to 𝐶 = 0.23 at 𝑡 = 166.854 s at which point the magnitude increases much more rapidly until the cut-off trigger stops the Trial. Trial #2 Results For this Trial, the flow velocity was 𝑈𝑒 = 0 cm/s and the feedback gain was incremented by −0.01. Similarly to Trial 1, decreasing the feedback gain to 𝐶 = −0.25 over the first 74.649 s resulted in no discernible oscillatory behavior that can be seen from the high-pass-filtered curvature data in Fig.4.8. The amplitude of peak frequency for this period, shown in Fig.4.9, is on the order of ≈ 10−7 m−1 which is on the order of system noise. After decreasing the feedback gain to 66 0.25 0.2 C = --0.15 C = --0.20 C = --0.26 C = --0.29 C = --0.30 0.15 0.1 Curvature (m-1 ) 0.05 0 -0.05 -0.1 -0.15 -0.2 0 20 40 60 80 100 120 140 Time (s) Figure 4.8: Curvature measurements for the full Trial of beam in 0 cm/s flow with negative feedback gain, that has been two-pass filtered - see Section 4.2. The solid vertical lines correspond to when the feedback gain was decreased by 0.01. Prior to decreasing the feedback gain to 𝐶 = −0.26 at 𝑡 = 74.649 s, there is no discernible oscillatory behavior, but decreasing the feedback gain to 𝐶 = −0.26 causes the system to show a small but marked oscillatory behavior with a distinct frequency of 𝜔𝑐𝑟 = 4.786 Hz. Upon lowering the feedback gain to 𝐶 = −0.30 at 𝑡 = 104.737 s, the system gradually began oscillating with increasingly greater amplitude at a frequency of 𝜔𝑐𝑟 = 4.719 Hz. This limit cycle behavior remained steady until the Trial was manually ended at 𝑡 = 126.934 s, denoted by the dashed vertical line. 𝐶 = −0.26 at 𝑡 = 74.649 s, there is a noticeable jump of one order of magnitude in the amplitude of the response in Fig.4.9 (though indiscernible to one viewing the beam in person), which then continues to grow slowly as the feedback gain is further decreased to 𝐶 = −0.29. For the period of 𝐶 = −0.25 through 𝐶 = −0.29 from 𝑡 = 74.649 s to 𝑡 = 104.737 s, while the oscillation has an amplitude of peak-frequency on the order of ≈ 10−6 m−1 , it has a very distinct frequency of 𝜔𝑐𝑟 = 4.786 Hz as derived from the FFT procedure outlined in Section 4.2. This minuscule but steady oscillatory behavior is again akin to that of a limit cycle or of a system operating at marginal 67 10-2 C = -0.15 C = -0.20 C = -0.26 C = -0.29 C = -0.30 10-3 Magnitude of peak frequency 10-4 10-5 10-6 10-7 10-8 0 20 40 60 80 100 120 140 Time (s) Figure 4.9: Amplitude of the peak frequency for the full Trial of beam in 0 cm/s flow with negative feedback gain taken from a spectrogram with a window of 0.4 s (approximately 2 cycles), presented on a logarithmic scale. The solid vertical lines correspond to when the feedback gain was decreased by 0.01 while the dashed vertical line denotes where the Trial was manually ended at 𝑡 = 126.934 s. Prior to decreasing the feedback gain to 𝐶 = −0.26 at 𝑡 = 74.649 s, the magnitude of the peak frequency was consistently on the order of ≈ 10−7 m−1 . Decreasing the feedback gain to 𝐶 = −0.26 causes a noticeable jump of one order of magnitude in the amplitude of the response after which the magnitude slowly increases with decreasing feedback gain. When the feedback gain is set to 𝐶 = −0.30 at 𝑡 = 104.737 s, the magnitude increases much more rapidly to the order of ≈ 10−2 m−1 where upon it remains fairly constant until the Trial is ended manually at 𝑡 = 126.934 s. stability. Upon lowering the feedback gain to 𝐶 = −0.30 at 𝑡 = 104.737 s, the system gradually began oscillating with increasingly greater amplitude (see Fig.4.8), taking approximately 15.1 s (≈ 70 oscillations) to reach the maximum amplitude of ≈ 10−2 m−1 . It then held this amplitude of oscillation fairly steady in the manner of a limit cycle until the Trial was manually ended at 68 𝑡 = 126.934 s. For this period, the oscillation has a very distinct frequency of 𝜔𝑐𝑟 = 4.719 Hz which was insignificantly lower than that of the previous period of extremely small oscillation. The FFT also reports very noticeable higher harmonics of the critical frequency. Trial #3 Results 0.6 C = 0.16 C = 0.22 C = 0.25 C = 0.27 C = 0.29 0.5 0.4 Curvature (m-1 ) 0.3 0.2 0.1 0 -0.1 -0.2 0 10 20 30 40 50 60 70 80 90 100 Time (s) Figure 4.10: Curvature measurements for the full Trial of beam in 10 cm/s flow with positive feedback gain, that has been two-pass filtered - see Section 4.2. The solid vertical lines correspond to when the feedback gain was increased by 0.01. Prior to increasing the feedback gain to 𝐶 = 0.22 at 𝑡 = 32.327 s, there is no discernible oscillatory behavior, but increasing the feedback gain to 𝐶 = 0.22 causes the system to show a small but marked oscillatory behavior with a frequency of 𝜔𝑐𝑟 = 0.330 Hz. 𝐶 = 0.29 at 𝑡 = 97.793, the system very quickly lost stability through flutter before being stopped by the cut-off trigger at 𝑡 = 100.639 s (dashed vertical line). For this Trial, the flow velocity was 𝑈𝑒 = 10 cm/s and the feedback gain was incremented by +0.01. The feedback gain was quickly increased to 𝐶 = 0.15 as it was known that the flutter behavior would not manifest at such a low feedback gain, after which it was increased more slowly. 69 10-2 C = 0.16 C = 0.22 C = 0.25 C = 0.27 C = 0.29 10-3 Magnitude of peak frequency 10-4 10-5 10-6 10-7 0 10 20 30 40 50 60 70 80 90 100 Time (s) Figure 4.11: Amplitude of the peak frequency for the full Trial of beam in 10 cm/s flow with positive feedback gain taken from a spectrogram with a window of 4.0 s (approximately 2 cycles), presented on a logarithmic scale. The solid vertical lines correspond to when the feedback gain was increased by 0.01 while the dashed vertical line denotes where the cut-off trigger was activated at 𝑡 = 100.639 s. Prior to increasing the feedback gain to 𝐶 = 0.22 at 𝑡 = 32.327 s, the magnitude of the peak frequency was on the order of ≈ 10−6 m−1 , but increasing the feedback gain to 𝐶 = 0.22 and beyond causes the magnitude of the peak frequency to increase by several orders of magnitude with increasing feedback gain. This behavior continues slowly until the feedback gain is set to 𝐶 = 0.29 at 𝑡 = 97.793 s at which point the magnitude increases very rapidly until the cut-off trigger stops the Trial. Similar to the previous two Trials, increasing the feedback gain to 𝐶 = 0.21 over the first 32.327 s resulted in no discernible oscillatory behavior that can be seen from the high-pass-filtered curvature data in Fig.4.10. The amplitude of peak frequency shown in Fig.4.11 for this period is on the order of ≈ 10−6 m−1 . Upon increasing the feedback gain to 𝐶 = 0.22 at 𝑡 = 32.327 s, there begins to be a noticeable increase in the amplitude of the response which continues to grow as the feedback gain is further increased. Imperfections of the experimental setup cause the curvature data to 70 have non-perfect sinusoidal measurements. For the period of 𝐶 = 0.25 through 𝐶 = 0.27 from 𝑡 = 45.019 s to 𝑡 = 87.401 s, the system sustains a very slight but consistent oscillation as seen in Fig.4.10. The oscillation during this period has an amplitude of peak-frequency on the order of ≈ 10−4 m−1 with a frequency of 𝜔𝑐𝑟 = 0.330 Hz. Upon raising the feedback gain to 𝐶 = 0.29 at 𝑡 = 97.793 s, the system very quickly lost stability through flutter - see Fig.4.10. The increase in amplitude was very sharp and rapidly triggered the torque cut-off at 𝑡 = 100.639 s. Given the windowing of the spectrogram and the short time in which the sytem activated the cut-off trigger, the last point of the amplitude of peak frequency plot in Fig.4.11 is based in large part on the data from before the final increase in the feedback gain to 𝐶 = 0.29 at 𝑡 = 97.793 s. While there were not sufficient oscillations to apply the FFT method to the region of flutter, by measuring the half-period from the first peak to the first trough, the frequency can be approximated to be 𝜔𝑐𝑟 ≈ 0.458 Hz, again significantly higher that than of the oscillations of the region before. Trial #4 Results For this Trial, the flow velocity was 𝑈𝑒 = 10 cm/s and the feedback gain was incremented by −0.01. Overall, the behavior of the system is very similar to that of Trial 2 which had no flow and negative feedback gain. Decreasing the feedback gain to 𝐶 = −0.27 over the first 56.318 s resulted in no discernible oscillatory behavior that can be seen from the high-pass-filtered curvature data in Fig.4.12. The amplitude of the peak frequency shown in Fig.4.13 for this period is on the order of ≈ 10−7 m−1 . Upon decreasing the feedback gain to 𝐶 = −0.28 at 𝑡 = 56.318 s, there begins to be a noticeable increase in the amplitude of the response which continues to grow by about two orders of magnitude as the feedback gain is further decreased to 𝐶 = −0.30. For this period of 𝐶 = −0.28 through 𝐶 = −0.30 from 𝑡 = 56.318 s to 𝑡 = 96.823 s, while the oscillation grows to an amplitude of peak-frequency on the order of ≈ 10−5 m−1 , it has a very distinct frequency of 𝜔𝑐𝑟 = 4.814 Hz as derived from the FFT procedure outlined in Section 4.2. Upon lowering the feedback gain to 𝐶 = −0.31 at 𝑡 = 96.823 s, the system gradually began 71 0.25 0.2 C = --0.20 C = --0.25 C = --0.28 C = --0.31 0.15 0.1 Curvature (m-1 ) 0.05 0 -0.05 -0.1 -0.15 -0.2 0 20 40 60 80 100 120 Time (s) Figure 4.12: Curvature measurements for the full Trial of beam in 10 cm/s flow with negative feedback gain, that has been two-pass filtered - see Section 4.2. The solid vertical lines correspond to when the feedback gain was decreased by 0.01. Prior to decreasing the feedback gain to 𝐶 = −0.28 at 𝑡 = 56.318 s, there is no discernible oscillatory behavior, but decreasing the feedback gain to 𝐶 = −0.28 causes the system to show a small but marked oscillatory behavior with a distinct frequency of 𝜔𝑐𝑟 = 4.814 Hz. Upon lowering the feedback gain to 𝐶 = −0.31 at 𝑡 = 96.823 s, the system gradually began oscillating with increasingly greater amplitude at a frequency of 𝜔𝑐𝑟 = 4.741 Hz. This limit cycle behavior remained steady until the Trial was manually ended at 𝑡 = 118.549 s, denoted by the dashed vertical line. oscillating with increasingly greater amplitude (see Fig.4.12), taking approximately 5.5 s (approx- imately 25 oscillations) to reach the maximum amplitude on the order of ≈ 10−2 m−1 . It then holds this amplitude of oscillation fairly steady in the manner of a limit cycle until the Trial was manually ended at 𝑡 = 118.549 s. For this period, the oscillation has a very distinct frequency of 𝜔𝑐𝑟 = 4.741 Hz which was insignificantly lower than that of the previous period of extremely small oscillation. The FFT also reports very noticeable higher harmonics of the critical frequency. 72 10-2 C = --0.20 C = --0.25 C = --0.28 C = --0.31 10-3 Magnitude of peak frequency 10-4 10-5 10-6 10-7 10-8 0 20 40 60 80 100 120 Time (s) Figure 4.13: Amplitude of the peak frequency for the full Trial of beam in 10 cm/s flow with negative feedback gain taken from a spectrogram with a window of 0.4 s (approximately 2 cycles), presented on a logarithmic scale. The solid vertical lines correspond to when the feedback gain was decreased by 0.01 while the dashed vertical line denotes where the Trial was manually ended at 𝑡 = 118.549 s. Prior to decreasing the feedback gain to 𝐶 = −0.28 at 𝑡 = 56.318 s, the magnitude of the peak frequency was consistently on the order of ≈ 10−7 m−1 . Decreasing the feedback gain to 𝐶 = −0.28 and further causes a noticeable increases in the amplitude of the peak frequency with decreasing feedback gain. When the feedback gain is set to 𝐶 = −0.31 at 𝑡 = 96.823 s, the magnitude gradually increases to the order of ≈ 10−2 m−1 where upon it remains fairly constant until the Trial is ended manually at 𝑡 = 118.549 s. Trial #5 Results For this Trial, the flow velocity was 𝑈𝑒 = 20 cm/s and the feedback gain was incremented by +0.01. Increasing the feedback gain to 𝐶 = 0.22 over the first 60.887 s resulted in no discernible oscillatory behavior that can be seen from the high-pass-filtered curvature data in Fig.4.14. Derived from the spectrogram of the curvature data (described in Section 4.2), the amplitude of peak frequency shown in Fig.4.15 for this period is on the order of ≈ 10−6 m−1 . Upon increasing the feedback gain 73 0.4 0.3 C = 0.16 C = 0.23 C = 0.30 C = 0.37 0.2 0.1 Curvature (m-1 ) 0 -0.1 -0.2 -0.3 -0.4 -0.5 0 50 100 150 200 250 Time (s) Figure 4.14: Curvature measurements for the full Trial of beam in 20 cm/s flow with positive feedback gain, that has been two-pass filtered - see Section 4.2. The solid vertical lines correspond to when the feedback gain was increased by 0.01. Prior to increasing the feedback gain to 𝐶 = 0.23 at 𝑡 = 60.887 s, there is no discernible oscillatory behavior, but increasing the feedback gain to 𝐶 = 0.23 causes the system to show a small but marked oscillatory behavior with a frequency of 𝜔𝑐𝑟 = 0.385 Hz. Upon raising the feedback gain to 𝐶 = 0.37 at 𝑡 = 196.946 s, the system gradually began oscillating with increasingly greater amplitude at a frequency of 𝜔𝑐𝑟 = 0.426 Hz. This limit cycle behavior remained steady until the Trial was manually ended at 𝑡 = 196.936 s, denoted by the dashed vertical line. to 𝐶 = 0.23 at 𝑡 = 60.887 s, there begins to be a sporadic but noticeable increase of approximately one order of magnitude in the amplitude of the response which continues to grow slowly as the feedback gain is further increased. For the period of 𝐶 = 0.25 through 𝐶 = 0.36 from 𝑡 = 98.135 s to 𝑡 = 196.936 s, the system sustains a very slight but consistent oscillation as seen in Fig.4.14. The oscillation during this period has an amplitude of peak-frequency on the order of ≈ 10−5 m−1 with a frequency of 𝜔𝑐𝑟 ≈ 0.385 Hz though unlike the other FFT data sets, the peak frequency is not a clear single frequency. 74 100 C = 0.16 C = 0.23 C = 0.30 C = 0.37 Magnitude of peak frequency 10-2 10-4 10-6 10-8 0 50 100 150 200 250 Time (s) Figure 4.15: Amplitude of the peak frequency for the full Trial of beam in 20 cm/s flow with positive feedback gain taken from a spectrogram with a window of 4.0 s (approximately 2 cycles), presented on a logarithmic scale. The solid vertical lines correspond to when the feedback gain was increased by 0.01 while the dashed vertical line denotes where the cut-off trigger was activated at 𝑡 = 196.936 s. Prior to increasing the feedback gain to 𝐶 = 0.23 at 𝑡 = 60.887 s, the magnitude of the peak frequency was on the order of ≈ 10−6 m−1 , but increasing the feedback gain to 𝐶 = 0.23 and beyond causes the magnitude of the peak frequency to increase in a very inconsistent manner by an order of magnitude with increasing feedback gain. When the feedback gain is set to 𝐶 = 0.37 at 𝑡 = 196.946 s, the magnitude gradually increases to the order of ≈ 10−1 m−1 where upon it remains constant until the Trial is ended manually at 𝑡 = 233.120 s. Upon raising the feedback gain to 𝐶 = 0.37 at 𝑡 = 196.946 s, the system very quickly lost stability through flutter - see Fig.4.14. The increase in amplitude was of four orders of magnitude (≈ 10−5 m−1 to ≈ 10−1 m−1 ) and occurred in a little over 5 seconds (2.5 cycles of oscillation). Unlike the previous two Trials with positive feedback gain which triggered the safety cut-off due to their increasingly larger oscillations, the system for this Trial entered a large limit cycle oscillation which continued until the system was manually turned off at 𝑡 = 233.120 s. During this period, 75 the beam oscillated with a frequency of 0.577 Hz which is significantly higher than the period preceding the final increase in feedback gain. This frequency is also higher than that of the large flutter of the previous two Trials with positive feedback gain, Trials 1 and 3. Trial #6 Results 0.3 C = --0.21 C = --0.25 C = --0.29 C = --0.35 0.2 0.1 Curvature (m-1 ) 0 -0.1 -0.2 -0.3 0 20 40 60 80 100 120 140 Time (s) Figure 4.16: Curvature measurements for the full Trial of beam in 20 cm/s flow with negative feedback gain, that has been two-pass filtered - see Section 4.2. The solid vertical lines correspond to when the feedback gain was decreased by 0.01. Prior to decreasing the feedback gain to 𝐶 = −0.29 at 𝑡 = 55.017 s, there is no discernible oscillatory behavior, but decreasing the feedback gain to 𝐶 = −0.26 causes the system to show a small but marked oscillatory behavior with a distinct frequency of 𝜔𝑐𝑟 = 4.786 Hz. Upon lowering the feedback gain to 𝐶 = −0.35 at 𝑡 = 115.443 s, the system gradually began oscillating with increasingly greater amplitude at a frequency of 𝜔𝑐𝑟 = 4.719 Hz. This limit cycle behavior remained steady until the Trial was manually ended at 𝑡 = 126.934 s, denoted by the dashed vertical line. For this Trial, the flow velocity was 𝑈𝑒 = 20 cm/s and the feedback gain was incremented 76 10-2 C = --0.21 C = --0.25 C = --0.29 C = --0.35 10-3 Magnitude of peak frequency 10-4 10-5 10-6 10-7 10-8 0 20 40 60 80 100 120 140 Time (s) Figure 4.17: Amplitude of the peak frequency for the full Trial of beam in 20 cm/s flow with negative feedback gain taken from a spectrogram with a window of 0.4 s (approximately 2 cycles), presented on a logarithmic scale. The solid vertical lines correspond to when the feedback gain was decreased by 0.01 while the dashed vertical line denotes where the Trial was manually ended at 𝑡 = 126.934 s. Prior to decreasing the feedback gain to 𝐶 = −0.29 at 𝑡 = 55.017 s, the magnitude of the peak frequency was consistently on the order of ≈ 10−7 m−1 . Decreasing the feedback gain to 𝐶 = −0.29 causes a very slight increase in the magnitude of oscillation to ≈ 10−6 m−1 which was maintained with decreasing feedback gain. When the feedback gain is set to 𝐶 = −0.35 at 𝑡 = 115.443 s, the magnitude increases much more rapidly to the order of ≈ 10−2 m−1 where upon it remains fairly constant until the Trial is ended manually at 𝑡 = 136.961 s. by −0.01. Overall, the behavior of the system is similar to that of Trials 2 and 4 which had no or slower flow and negative feedback gain but there are distinct difference. As with the previous Trials, the initial decreasing of the feedback gain produced no discernible effect visible on the curvature data shown in Fig.4.16 and as further evinced by the magnitude of oscillation of the peak frequency (≈ 10−7 m−1 ) shown in Fig.4.17. At 𝑡 = 55.017 s the feedback gain was decreased to 𝐶 = −0.29 which caused a very slight increase in the magnitude of oscillation to ≈ 10−6 m−1 77 which was maintained until the feedback gain was decreased to 𝐶 = −0.35 at 𝑡 = 115.443 s. This slight increase prior to 𝐶 = −0.35 was significantly less apparent than that of the previous Trials. Fig.4.16 shows that, unlike all previous Trials, the magnitude of oscillation did not begin increasing immediately upon the final increment of feedback gain, but continued to maintain the same magnitude until 4.1 s later when it suddenly had a significantly larger oscillation. This magnitude of oscillation was maintained for 0.7 s after which it began increasing by several orders of magnitude as can be seen in both Fig.4.16 and Fig.4.17. The beam then fluttered at a relatively steady but slightly varying amplitude until the Trial was manually ended at 𝑡 = 136.961 s. During this period of the large limit cycle behavior, the system oscillated with a very distinct frequency of 4.932 Hz. Trial #7 Results The system in 30 cm/s flow and with no feedback did not exhibit any oscillatory behavior beyond what could be considered noise as can be seen in Fig.4.18. For the duration of this Trial, there was no distinct frequency of oscillation and the magnitude of curvature measurement was significantly smaller than even those of the small oscillatory periods of the previous Trials. This supports the claim that the flutter induced in the system was due entirely to the feedback rather than the fluid flow. 78 10-3 8 6 4 Curvature (m-1 ) 2 0 -2 -4 -6 0 2 4 6 8 10 12 14 Time (s) Figure 4.18: Curvature measurements for the full Trial of beam in 30 cm/s flow without feedback gain, that has been two-pass filtered - see Section 4.2. General Trends and Conclusions Table 4.2: Summary of results for all of the Trials. Values in parentheses denote low precision values not taken from a distinct peak of the FFT figures. Flow Gain Critical Critical Minor Oscillatory Trial # Velocity (cm/s) Sign Gain Frequency (Hz) Frequency (Hz) 1 0 + 0.23 0.426 0.301 2 0 - -0.30 4.719 4.810 3 10 + 0.29 (0.458) 0.331 4 10 - -0.31 4.741 4.778 5 20 + 0.37 0.577 0.390 6 20 - -0.35 4.932 (4.807) 7 30 N/A - - - There are a few interesting observations that can be taken from the results of the different Trials which are summarized in Table 4.2. The simplest is that the positive feedback gain Trials 1, 3, 79 and 5 are much more similar to each other than to the negative feedback gain Trials 2, 4, and 6, with respect to the frequency of flutter (≈ 0.5 Hz vs. ≈ 5 Hz) as well as the minor oscillation frequency before the larger flutter behavior (≈ 0.35 Hz vs. ≈ 5 Hz). This is despite the differences in flow velocity. Of interest is that the minor oscillation frequency before the large flutter behavior for the positive feedback gain Trials was significantly smaller than that of the frequency of flutter oscillations, while there was no distinct difference between the minor oscillation and larger flutter frequencies for the negative feedback gain Trials. For both cases, it is seen how the magnitude of the critical gain increased with flow velocity with both cases increasing at a greater than linear rate but with the positive feedback gain Trials increasing at a much greater rate than the negative feedback gain Trials. 4.4 Numerical simulation To verify that the theory used in Chapter 3 is consistent with the experimental results presented in Section 4.3, we use the parameters of the experimental beam in the model. In order to do this we first determine the non-dimensionalized values of the beam parameters, fluid flow, feedback gain, and flutter frequencies to input into the model. As the analysis will be done non-dimensionally, these values will also be used to re-dimensionalize the results attained. Additionally, we consider some extensions to the model that can help it better represent the experimental system. Using the process of extending the frequency loci from a specific flow velocity for a given curvature measurement location (see Section 3.3, particularly Figs.3.5 and 3.9), we shall determine the critical gain and critical frequency at which the model predicts the system loses stability. We shall then summarize the numerical results and compare them to those gathered from the experiment. Non-Dimensional Terms/Scaling Referring to Section 3.1, the non-dimensional values needed for the equation of motion (3.4) and the moment based on curvature boundary condition (3.6a) are 80 r s 𝑚e b𝑥 𝑚e 𝐸𝐼 𝐶m,𝑐 𝛽= , 𝛾= , 𝑢 e = 𝑈e 𝐿 , 𝜏=𝑡 , 𝑐 m,𝑐 = 𝑚e + 𝑚b 𝐿 𝐸𝐼 (𝑚 e + 𝑚 b )𝐿 4 𝐸𝐼 Using the values given in Section 4.1, we determine the following values or ratios: 𝑢e 𝜏 𝑐 m,𝑐 𝛽 = 0.984, 𝛾 = 0.25, = 8.21, = 0.317, = 57.2 𝑈e 𝑡 𝐶m,𝑐 Model Adjustments needed for Time Delay Referring again to Section 3.1, the model, assuming perfect and instantaneous curvature measure- ment and torque actuation, does not account for several aspects of the physical realization of the experiment. Particularly, the model does not account for the inertia of the shaft complex (motor, coupler, and shaft) which transmits the moment to the leading edge of the beam, nor the time delays inherent in acquiring, processing, and applying the moment to the beam. 𝑦 𝜕𝑦 𝜃 0 = 𝜕𝑥0 M 𝑥 𝐽 Figure 4.19: A flexible beam, connected at one end to a fixed point by a revolute joint and free at the other end, is immersed in a fluid flowing with constant velocity 𝑈e . In Fig.4.19 we show the actuated moment M which rotates the shaft complex (rotational inertia 𝐽) which in turn applies the moment to the leading edge of the beam. This actuated moment is proportional to the curvature measurement at 𝑥ˆ that has been delayed by a factor 𝑡0 . Summing the applied moment, the reaction moment of the beam (defined in Section 3.1), and rotational inertia of the shaft complex about the origin, we get the following equation: 81 𝜕 2 𝑦(0, 𝑡) 𝜕 2 𝑦(0, 𝑡) 𝜕 3 𝑦(0, 𝑡) 𝜕 2 𝑦(ˆ 𝑥 , 𝑡 − 𝑡0 ) M(𝐶m,𝑐 , 𝑥ˆ, 𝑡0 ) − 𝐸 𝐼 = 𝐽 𝜃¥0 → 𝐸𝐼 + 𝐽 = 𝐶 m,𝑐 𝜕𝑥 2 𝜕𝑥 2 𝜕𝑥𝜕𝑡 2 𝜕𝑥 2 (4.3) This equation is similar to that of (3.3a) but incorporates the shaft complex inertia and time delay. Following the same process as in Section 3.1, we non-dimensionalize (4.3) to get 𝜕 2 𝑣(0, 𝜏) 𝜕 3 𝑣(0, 𝜏) 𝜕 2 𝑣(𝛾, 𝜏 − 𝜏0 ) + 𝜆 = 𝑐 m,𝑐 (4.4) 𝜕𝑢 2 𝜕𝑢𝜕𝜏 2 𝜕𝑢 2 where 𝜆 is the non-dimensional inertia of the shaft complex and 𝜏0 is the non-dimensional time delay between the instantaneous curvature of the beam at 𝛾 and the proportional moment being applied to the leading edge of the beam. These terms are non-dimensionalized from the dimensional terms by s 𝐽 𝐸𝐼 𝜆= , 𝜏0 = 𝑡0 (𝑚 𝑒 + 𝑚 𝑏 )𝐿 2 (𝑚 e + 𝑚 b )𝐿 4 For 𝐽 = 1.349×10−5 we get 𝜆 = 1.128×10−5 while the ratio of the non-dimensional to dimensional time delay was 𝜏0 /𝑡0 = 0.3173, which corresponds to the same ratio as that in Section 4.4. While the curvature sensing to torque command process was shown to occur in the duration of a single control loop i.e. 1 ms, the flexibility of the shaft complex caused a noticeable lag between when the torque was commanded and actually applied to the leading edge of the beam. Evaluating an open-loop frequency sweep showed that the time delay could range from a minimum of 3 ms to as great as 50 ms depending on the amplitude of the oscillation which was itself a function of the frequency i.e. resonance. Through trial and error of the simulation model, the time delay for the positive feedback gain Trials was taken as 𝑡0 = 0 ms while the time delay for the negative feedback gain Trials was taken as 𝑡0 = 30 ms. This difference is possibly due to the higher frequency of oscillation of the negative feedback gain Trials causing the flexion of the shaft complex to have a relatively greater impact on the response behavior. 82 Continuing with the solution process laid out in Section 3.2, using the assumed form of the solution (3.7) i.e. 𝑣(𝑢, 𝜏) = 𝑓 (𝑢)𝑒𝑖Ω𝜏 , we get the following form for the boundary condition. This equation now explicitly includes the frequency unlike (3.10b) 𝑓 ′′(0) − 𝜆Ω2 𝑓 ′(0) = 𝑐 m,𝑐 𝑓 ′′(𝛾)𝑒 −𝑖Ω𝜏0 (4.5) Substitution of the assumed form of the solution of 𝑓 (𝑢) (see (3.12) in Section 3.2) into (4.5) gives the following equation       1 1 1 1   𝐴1  0           𝑧 2 𝑒 𝑧1 𝑧 2 𝑒 𝑧2 𝑧 2 𝑒 𝑧3 𝑧 2 𝑒 𝑧4   𝐴2  0  1 2 3 4       =  (4.6) 𝑧3 𝑒 𝑧1 𝑧3 𝑒 𝑧2 𝑧3 𝑒 𝑧3 𝑧3 𝑧4       1 2 3 4 𝑒   𝐴3  0            𝛿1 𝛿2 𝛿3 𝛿4   𝐴4  0      | {z } Z where the 𝑧 𝑛 terms are the the roots of the characteristic equation of motion (3.11) and the 𝛿𝑛 , 𝑛 = 1, 2, 3, 4, are defined as follows 𝛿𝑛 , 𝑧 2𝑛 − 𝑐 m,𝑐 𝑧 2𝑛 𝑒 𝑧𝑛 𝛾−𝑖Ω𝜏0 − 𝜆Ω2 𝑧 𝑛 Non-trivial solutions of (4.6), i.e. complex frequencies Ω𝑖 , 𝑖 = 1, 2, · · ·, can be obtained by numerically solving the transcendental equation det (Z) = 0 given specific parameter values and feedback gain. Numerical Results The Argand diagrams of the simulation results are shown in Fig.4.20 for the positive feedback gain and Fig.4.21 for the negative feedback gain. The units along both axes have been converted to dimensional Hz as described in Section 4.4. Each figure shows the loci of the first four modes for a series of configurations with the regions of frequency between the modes removed so as to better be able to discern the behavior of the system. One set of loci corresponds to increasing the 83 0.3 Flow velocity (m/s) Ue 0.3 0.25 Ue = 0.01 Ue = 0.10 0.2 Ue = 0.20 Im[ ] Hz 0.15 0.1 0.05 0 0 0.4 0.8 2 2.5 4.8 5 5.2 5.4 8.8 9 9.2 Mode 1 Mode 2 Mode 3 Mode 4 Re[ ] Hz Figure 4.20: Argand diagram of the first four modes for the simulation of the positive feedback gain Trials with time delay of 𝑡0 = 0.00 s. The x-axis has units of dimensional Hz with each of the four modes of interest shown in focus with the empty frequency ranges between them removed. We can see that the system loses stability in the first mode at each flow velocity: 𝑈𝑒 = 1 and 10 cm/s resulting in flutter at 0.27 Hz and 0.225 Hz, respectively, while 𝑈𝑒 = 20 cm/s results in divergence at the origin. The loci of the other modes do not approach the Re[Ω] axis nearly as quickly as the feedback gain is increased. flow rate from 𝑈𝑒 = 0 cm/s to 𝑈𝑒 = 30 cm/s with zero feedback gain, while the other three sets of loci correspond to increasing/decreasing the feedback gain while keeping the flow rate constant at the three values4 tested in the experiment: 𝑈𝑒 = 1, 10, and 20 cm/s. For these, the feedback gain was increased/decreased until one of the loci reached the Re[Ω] axis. Shared between both figures, the loci corresponding to increasing the flow rate from 𝑈𝑒 = 0 cm/s to 𝑈𝑒 = 30 cm/s with zero feedback gain, each continue to become more stable as 𝑈𝑒 increases. This corroborates the experimental observation in Trial 7 in Section 4.3 that the system did not lose stability due to flow alone. As can be seen in Fig.4.20, the simulation for the positive feedback gain system lost stability 4The simulated system corresponding to the experimental Trials with 𝑈𝑒 = 0 cm/s were simulated with 𝑈𝑒 = 1 cm/s so as to incorporate the damping term of the equation of motion (3.4). Otherwise, the system would be completely ignoring the damping effects of the surrounding fluid. 84 0.3 Flow velocity (m/s) U 0.3 e 0.25 Ue = 0.01 Ue = 0.10 0.2 Ue = 0.20 Im[ ] Hz 0.15 0.1 0.05 0 0.6 0.7 0.8 2.4 2.6 2.8 5.2 5.4 8.8 8.9 9 9.1 Mode 1 Mode 2 Mode 3 Mode 4 Re[ ] Hz Figure 4.21: Argand diagram of the first four modes for the simulation of the negative feedback gain Trials with time delay of 𝑡0 = 0.030 s. The x-axis has units of dimensional Hz with each of the four modes of interest shown in focus with the empty frequency ranges between them removed. We can see that the system loses stability in the third mode at each flow velocity: 𝑈𝑒 = 1, 10, and 20 cm/s crossing the Re[Ω] axis at frequencies of approximately 5.276, 5.315, and 5.343 Hz, respectively, signifying flutter. The lower modes also move towards the Re[Ω] axis while the fourth mode moves away. There is a mild but noticeable jump from the unforced loci given the inclusion of the rotor, coupler, and shaft inertia and time lag but this does not substantially affect the results shown. through the first mode for each flow velocity. The loci originating from the zero-feedback locus of the first mode for 𝑈𝑒 = 1 and 10 cm/s crossed the Re[Ω] axis at a very low frequency of approximately 0.27 Hz and 0.225 Hz, respectively, signifying low frequency flutter. The loci corresponding to 𝑈𝑒 = 20 cm/s reach the Re[Ω] axis at the origin signifying divergence rather than flutter. The loci of the higher modes with positive feedback gain also moved towards the Re[Ω] axis but not as quickly. The simulation for the negative feedback gain system shown in Fig.4.21 lost stability through the third mode for each flow velocity. The loci originating from the zero-feedback locus of the third mode for 𝑈𝑒 = 1, 10, and 20 cm/s crossed the Re[Ω] axis at frequencies of approximately 5.276, 5.315, and 5.343 Hz, respectively, signifying flutter. The loci of the lower modes with negative 85 feedback gain also moved towards the Re[Ω] axis but not as quickly while the loci of fourth mode moved away from the Re[Ω] axis. The initial abrupt movement from the zero-feedback loci was due to the introduction of the 𝑡0 = 0.030 s time delay introduced to have the results more closely match the experiment. 4.5 Conclusion: Comparison of Numerical and Experimental Results The results from the simulations are summarized in Table 4.3 alongside the pertinent results from the experimental results previously shown in Table 4.2 in Section 4.3. For both signs of the feedback, the stability was lost through flutter at approximately the correct frequency corresponding to the correct mode. Significantly, the loci of the other modes were either becoming more stable with increasing magnitude of feedback gain, or approaching the Re[Ω] axis considerably slower than the loci of the mode that lost stability. Additionally, as the flow velocity increased, a greater magnitude of feedback was required for both the experimental and simulated results, with respect to both the negative and positive feedback systems. In this regards, a point of strong disagreement was that the positive feedback gain jumped by considerable margins in the experimental results for the change in flow velocity but changed very slightly in the simulation results. The negative feedback gain showed the exact opposite comparison with the experimental results changing very Table 4.3: Summary of results for the numerical simulations of the system for all of the Trials alongside the re-presented experimental results from Section 4.3. Values in square brackets denote a simulation that lost stability through divergence. Values in parentheses denote an experimental frequency that was not derived from FFT analysis. Simulations for Trials with negative feedback had a 𝑡0 = 0.030 s time delay. Simulation Experimental Flow Gain Critical Critical Critical Critical Trial # Velocity (cm/s) Sign Gain Frequency (Hz) Gain Frequency (Hz) 1 0-1 + 0.0252 0.270 0.23 0.426 2 0-1 - -0.00061 5.276 -0.30 4.719 3 10 + 0.0262 0.225 0.29 (0.458) 4 10 - -0.00577 5.315 -0.31 4.741 5 20 + [0.0281] [0.0] 0.37 0.577 6 20 - -0.0110 5.343 -0.35 4.932 7 30 0 N/A N/A N/A N/A 86 slightly while the simulation results showed considerable increase with increase in flow velocity. The large disparity of the feedback gain values can be attributed in part to errors in the implementation of the scaling values described in Section 4.1. As each of these were only scaling terms, they can only account for an error by a fixed scale of the feedback gain, but would not have affected the behavior of the system. Additionally, the model used was of a linearized Euler-Bernouli beam in non-vortical flow which does not capture the large oscillation behaviour well, especially for the low frequency, large amplitude flutter response of the positive feedback gain Trials. Finally, the flexibility of the rotating shaft and time lag were not well understood, and while modelled simplistically, are another source of potentially disparity between the experimental and theoretical results. 87 CHAPTER 5 EFFECTS OF A FOLLOWER FORCE ON THE INDUCEMENT OF FLUTTER IN A FLEXIBLE BEAM IN FLUID FLOW In this chapter, we continue to explore the flutter characteristics of a beam in fluid flow. However, rather than feedback-induced flutter, we look into how a follower force affects the onset of flutter given varying beam geometry and fluid flow velocity. We define the model of a cantilever-free beam with a tangential follower force in fluid flow with the cantilevered end being upstream. This follower force is a standard method of inducing flutter used in the literature and could analogously be generated by internal pipe flow with fluid jet emanating from a convergent nozzle, resulting in the same equations of motion. We describe the system as a beam here for continuity. Following the methods in Chapter 3, we determine the equations of motion and boundary conditions of the beam. We work through the method of solution required to determine the natural frequencies of the beam given the set of parameters such as beam properties, fluid flow velocity, and follower force amplitude. We begin by determining the behavior of the system over a wide range of beam-to-fluid mass ratios as the external flow velocity is increased with no follower force. We study the behavior of these curves which are used for finding the range of external flow velocities to study further. From this, we show the method of determining the follower force at which the system loses stability. We present the surfaces over the mass ratio and external flow velocity plane which correspond to the follower force required to induce flutter instability, the critical frequency at which the beam flutters, and the mode from which the beam lost stability. We explore some insights on the scaling of the system and observe some interesting convergence. We delve into the phenomenon known as Ziegler’s paradox and explore how it relates to the problem at hand. Finally, using the identical equations between a fluid-conveying pipe and a fluid-immersed beam, we compare the effects of a follower force to attaching a nozzle to a fluid conveying pipe. The layout of this chapter is organized as follows. The problem formulation is described in Section 5.1 and the analytical method of solution is discussed in Section 5.2. The investigation of the flutter instability is defined in Section 5.3. Parameters for the system and an exploration of the 88 behavior of the system without follower force are given in Section 5.3 with an explanation of the method of solution of the critical follower force given in Section 5.3. The critical stability results are presented in Section 5.4 with further explorations detailed in Section 5.4 (Ziegler Paradox) and Section 5.4 (Effluent Jet with Nozzle Attachment). The Problem Formulation (Section 5.1) and Method of Solution (Section 5.2) of this chapter are adapted from those of Chapter 3 which was published as [109] and included here with permission. 5.1 Problem Formulation Consider the fluid-immersed flexible beam in Fig.5.1. The beam has length 𝐿, a rectangular cross-section with width 𝑤 and height ℎ, mass per unit length 𝑚 b , and Young’s modulus of elasticity 𝐸. The upstream end of the beam is fixed while the downstream end of the beam is free. At the free end of the beam, a follower force 𝑃 is applied and acts tangentially to the slope of the beam. The fluid is inviscid and flows with constant velocity 𝑈e . The equation of motion of the beam with follower force, ignoring gravitational, viscous, pressurization and tensile effects, is as follows [91]: 𝜕 4 𝑦(𝑥, 𝑡)   2 2 2 2 + 𝑃 𝜕 𝑦(𝑥, 𝑡) + 2𝑚 𝑈 𝜕 𝑦(𝑥, 𝑡) + (𝑚 + 𝑚 ) 𝜕 𝑦(𝑥, 𝑡) = 0 (5.1) 𝐸𝐼 + 𝑚 𝑒 𝑈e 𝑒 e e b 𝜕𝑥 4 𝜕𝑥 2 𝜕𝑥𝜕𝑡 𝜕𝑡 2 where 𝑦(𝑥, 𝑡) is the displacement of the beam, 𝐼 = (ℎ𝑤 3 /12) is the area moment of inertia of the beam, and 𝑚 𝑒 is the mass per unit length of the external fluid. The mass per unit length of the fluid is approximated as the mass of water within the cylinder of unit length circumscribing the 𝑦 𝑤 𝑈e 𝑃 𝑥 𝐿 Figure 5.1: A flexible beam, rigidly connected at one end and free at the other end, is immersed in a fluid flowing with constant velocity 𝑈e . At the free end of the beam, a follower force 𝑃 is applied and acts tangentially to the slope of the beam. 89 beam cross-section [107]. The boundary conditions of the fixed end of the beam and the shear and moment boundary conditions of the free end of the beam are given as 𝜕𝑦(0, 𝑡) 𝜕 2 𝑦(𝐿, 𝑡) 𝜕 3 𝑦(𝐿, 𝑡) 𝑦(0, 𝑡) = = 𝐸𝐼 = 𝐸 𝐼 =0 (5.2) 𝜕𝑥 𝜕𝑥 2 𝜕𝑥 3 We introduce the following change of variables: s r 𝑥 𝑦 𝐸𝐼 𝑚e 𝑃𝐿 𝑢= , 𝑣= , 𝜏=𝑡 , 𝑢 e = 𝑈e 𝐿 , P= 𝐿 𝐿 (𝑚 e + 𝑚 b )𝐿 4 𝐸𝐼 𝐸𝐼 to obtain the non-dimensional equation of motion of the beam 𝜕4𝑣  2  𝜕2𝑣 p 𝜕2𝑣 𝜕2𝑣 + 𝑢 e + P + 2 𝛽 𝑢 e + =0 (5.3) 𝜕𝑢 4 𝜕𝑢 2 𝜕𝑢𝜕𝜏 𝜕𝜏 2 where 𝛽 is the mass fraction of the fluid to the fluid and beam: 𝑚e 𝛽= 𝑚e + 𝑚b From (5.2), the non-dimensional geometric boundary conditions of the fixed end of the beam and the non-dimensional natural boundary conditions of the free end of the beam are given as 𝜕𝑣(0, 𝜏) 𝜕 2 𝑣(1, 𝜏) 𝜕 3 𝑣(1, 𝜏) 𝑣(0, 𝜏) = = = = 0. (5.4) 𝜕𝑢 𝜕𝑢 2 𝜕𝑢 3 5.2 Method of Solution To solve (5.3) for the boundary conditions in (5.4), we followed the procedure introduced in [91] and used in [52]. In particular, we assume the following separable form for 𝑣(𝑢, 𝜏): 𝑣(𝑢, 𝜏) = 𝑓 (𝑢)𝑒𝑖Ω𝜏 (5.5) where Ω is the non-dimensional frequency of oscillation while 𝑓 (𝑢) is the spatially varying shape function. Substitution of (5.5) into (5.3) and (5.4) yields 90   p 𝑓 ′′′′(𝑢) + 𝑢 2e + P 𝑓 ′′(𝑢) + 2𝑢 e 𝛽 𝑖Ω 𝑓 ′(𝑢) − Ω2 𝑓 (𝑢) = 0 (5.6) and 𝑓 (0) = 𝑓 ′(0) = 𝑓 ′′(1) = 𝑓 ′′′(1) = 0 (5.7) Since (5.6) is an ordinary differential equation with constant coefficients, the solution of 𝑓 (𝑢) is assumed to be of the form 𝑓 (𝑢) = 𝐴𝑒 𝑧𝑢 . This results in the characteristic equation for the system   p 𝑧 4 + 𝑢 2e + P 𝑧 2 + 2𝑢 e 𝛽 𝑖Ω𝑧 − Ω2 = 0 (5.8) For specific values of 𝑢 e , P, and 𝛽, (5.8) provides four roots of 𝑧 𝑛 , 𝑛 = 1, 2, 3, 4, which are functions of Ω. The quartic equation (5.8), can be solved analytically for each of these parameters and Ω. These four roots cause the solution of 𝑓 (𝑢) to take the form 𝑓 (𝑢) = 𝐴1 𝑒 𝑧1 𝑢 + 𝐴2 𝑒 𝑧2 𝑢 + 𝐴3 𝑒 𝑧3 𝑢 + 𝐴4 𝑒 𝑧4 𝑢 (5.9) Substitution of the boundary conditions in (5.7) into (5.9) yields       1 1 1 1   𝐴1  0        𝑧1 𝑧 𝑧 𝑧    0 4   𝐴2   2 3     =  (5.10) 𝑧 2 𝑒 𝑧1 𝑧 2 𝑒 𝑧2 𝑧 2 𝑒 𝑧3 𝑧 2 𝑒 𝑧4   𝐴  0  1 2 3 4   3         3 𝑧 3 𝑧 3 𝑧 3 𝑧     𝑧 1 𝑒 1 𝑧 2 𝑒 2 𝑧 3 𝑒 3 𝑧 4 𝑒 4   𝐴4  0      | {z } Z Non-trivial solutions of (5.10) can be obtained by solving the transcendental equation det (Z) = 0 (5.11) For specific values of 𝑢 e , P, and 𝛽, the transcendental equation can be solved numerically to get the complex frequencies Ω𝑖 , 𝑖 = 1, 2, · · ·, and from these, the 𝑧 𝑛 terms, 𝑛 = 1, 2, 3, 4, for each Ω𝑖 . 91 5.3 Investigation of Flutter Instability Critical Stability While Section 5.2 provides the frequencies of oscillation, Ω𝑖 , 𝑖 = 1, 2, · · ·, for specific values of 𝑢 e , P, and 𝛽, we seek to find the critical stability points where the system loses stability through flutter. For a particular Ω and corresponding 𝑧 𝑛 terms, 𝑛 = 1, 2, 3, 4, the solution of (5.3) and (5.4) can be obtained by substituting (5.9) into (5.5): X4 X 4 𝑣(𝑢, 𝜏) = 𝐴𝑛 𝑒 𝑧𝑛 𝑢 𝑒𝑖Ω𝜏 = 𝑒 −I𝑚[Ω]𝜏 𝐴𝑛 𝑒 R𝑒[𝑧𝑛 ]𝑢 𝑒𝑖 {I𝑚[𝑧𝑛 ]𝑢+R𝑒[Ω]𝜏 } (5.12) 𝑛=1 𝑛=1 where the coefficients 𝐴𝑛 , 𝑛 = 1, 2, 3, 4, can be obtained from the null space of Z in (5.10). It is clear from (5.12) that the stability of 𝑣(𝑢, 𝜏) is dependent on the exponential term outside the summation; if I𝑚[Ω] < 0, this term is unbounded as 𝑡 → ∞. Therefore, the point at which I𝑚[Ω] changes sign from positive to negative represents the onset of flutter instability. Additionally, the first exponential term inside the summation is bounded because 𝑢 is bounded while the second exponential term, because it has an imaginary exponent, yields periodic motion. It should be noted that (5.12) describes the solution for one non-dimensional frequency Ω representing a single solution to the transcendental equation of (5.11). At the flutter instability point, one specific value of Ω, Ω = Ωc𝑟 , satisfies I𝑚[Ω] = 0 whereas all other Ω values satisfy I𝑚[Ω] > 0. The frequency Ωc𝑟 is real and is defined as the critical frequency. Since 𝑒 −I𝑚[Ω]𝜏 → 0 as 𝜏 → ∞ for all Ω 6= Ωc𝑟 , the complete solution at the flutter instability point takes the form X4 𝑣(𝑢, 𝜏) = 𝐴𝑛 𝑒 R𝑒[𝑧𝑛 ]𝑢 𝑒𝑖 {I𝑚[𝑧𝑛 ]𝑢+R𝑒[Ωc𝑟 ]𝜏 } (5.13) 𝑛=1 As the imaginary exponent in (5.13) is a function of both 𝑢 and 𝜏, the above equation represents a traveling waveform; more specifically, a combination of four superimposed traveling waveforms. 92 Simulation Environment Our ultimate goal is to determine the critical follower force Pc𝑟 required to cause the system to lose stability through flutter over a domain of 𝛽 and 𝑢 e . For this work, we will investigate flutter instability for the following range of 𝛽1 𝛽 = {0.001, 0.005, 0.01} ∪ {0.025, 0.05, ..., 0.95, 0.975} ∪ {0.99, 0.995} (5.14) High values of 𝛽 → 1, correspond to systems with very large fluid mass to beam mass ratios such as tall ℎ and thin 𝑤 beams in a dense fluid flow or thin-walled pipes. Low values of 𝛽 → 0, correspond to systems with very small fluid mass to beam mass ratios such as thick, fluid-conveying pipes. As the equations of motion are identical between the two systems, the whole range is treated here for completeness. To determine the range of 𝑢 e we first determine the natural frequencies Ω𝑖 , 𝑖 = 1, 2, · · ·, for the cantilevered-unforced system in zero flow, i.e., the system with 𝑢 e = 0 and P = 0. We then select a value of 𝛽 and increase 𝑢 e from 𝑢 e = 0, solving for the natural frequencies of the Ω 𝑗 , 𝑗 = 1, 2, · · · of the unforced system in flow. This is continued until the system loses stability as described in Section 5.3 at 𝑢 e = 𝑢 c𝑟 (𝛽) for each 𝛽. From this process, a different range of fluid flow velocity 𝑢 e is chosen for each value of 𝛽 as determined by the value of 𝑢 c𝑟 which caused the system to lose stability through flutter. Here we define the value of 𝑘 as the mode through which the system lost stability through flutter. It is over this 𝛽-𝑢 e domain that we will the determine the critical follower force Pc𝑟 required to cause the system to lose stability through flutter. While the bounds of the 𝛽-𝑢 e domain have thus been defined, the behavior of the zero follower force P = 0 system as 𝑢 c𝑟 increases is important to explore. As will be described in Section 5.3, it is from intermediary points along each of these increasing 𝑢 e loci that the follower force P will be increased until the system loses stability through flutter. To begin, we consider the critical flow velocities 𝑢 c𝑟 , critical frequencies Ωc𝑟 , and mode, 𝑘, at which the system loses stability through flutter for the zero follower force P = 0 case as summarized in Fig.5.2. Of particular note are 1 𝛽 = 0 will be treated separately as a special case in Section 5.4. 93 18 ucr 70 16 cr 60 14 12 50 10 u e,cr 40 cr 8 30 6 20 4 10 2 k=2 k=3 k=2 k=1 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 5.2: Critical flow velocity 𝑢 𝑐𝑟 (solid line) and critical frequency Ω𝑐𝑟 (dotted line) corresponding to each 𝛽 for the zero follower force P = 0 system. The critical mode 𝑘 is denoted for each section. the periods where even though the mode of flutter 𝑘 changes as 𝛽 gradually increases, there are only relatively minor changes in 𝑢 c𝑟 and Ωc𝑟 . Conversely, there are the periods of relatively major changes in 𝑢 c𝑟 and Ωc𝑟 as 𝛽 gradually increases even though the mode of flutter 𝑘 remains the same. To understand the causes of this behavior, we follow the loci of increasing 𝑢 e for the first three2 natural frequencies corresponding to select 𝛽. Given the large number of 𝛽 values, we divide the set into seven subsets that exhibit distinctly similar behavior. From these an illustrative select few are chosen, so as to reduce the clutter in each figure as well as better show the behavior of the system as 𝑢 e increases from 0 for each 𝛽. The Argand diagrams for each subset of 𝛽 loci are shown 2It was determined that there was no behaviour of note above the third natural frequency. Higher frequencies were calculated throughout this process, but are here omitted. 94 14 12 10 ] 8  I 6 𝛽↑ 4 2 𝛽↑ 𝛽↑ 0 10 20 30 40 50 60 k=1 k=2 R  k=3 Figure 5.3: Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.001, 0.1, 0.2, 0.275}. Loci of increasing 𝛽 are denoted as well as the direction in which each loci moves as 𝑢 𝑒 increases. As previously shown in Fig.5.2, while the fluid flow velocity required for the system to lose stability, 𝑢 𝑐𝑟 , increases with 𝛽, the critical frequency Ω𝑐𝑟 decreases with 𝛽. in Figs.5.3 - 5.9. Representing the first subset of 𝛽, Fig.5.3 shows the Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.001, 0.1, 0.2, 0.275}. For this subset, as 𝑢 e increases, each loci, beginning from the natural frequencies of the cantilevered-unforced system in zero flow, i.e., the system with 𝑢 e = 0 and P = 0, initially becomes more stable by moving away from the Re[Ω] axis. However, the loci of the 𝑘 = 2 mode curve back towards the Re[Ω] axis and eventually cross it, denoting the point of the system becoming unstable. As seen in Fig.5.2, for this range of 𝛽, while the fluid flow velocity required for the system to lose stability, 𝑢 𝑐𝑟 , increases with 𝛽, the critical frequency Ω𝑐𝑟 with which it flutters decreases with 𝛽. Referring to Fig.5.2, we see that at 𝛽 = 0.3, the behavior of the system shows a sharp change 95 25 20 𝛽↑ 15    10 𝛽↑ 𝛽↑ 5 𝛽↑ 0 10 20 30 40 50 60 k=1 k=2 k=3 Figure 5.4: Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.3, 0.325, 0.35, 0.375, 0.4}. Loci of increasing 𝛽 are denoted as well as the direction in which each loci moves as 𝑢 𝑒 increases. All but the loci corresponding to 𝛽 = 0.4 are the solid lines. The loci of 𝛽 = 0.4 is included to show the veering effect which explains the change in mode from 𝑘 = 2 to 𝑘 = 3 seen in Fig.5.2. It is important to note, however, that this veering does not correspond to a significant change in either Ω𝑐𝑟 nor 𝑢 𝑐𝑟 . Additionally, this figure shows how the curvature of the locus corresponding to 𝑘 = 2 mode of 𝛽 = 0.3 just misses the Re[Ω] axis before crossing it at a significantly latter point with respect to both Ω𝑐𝑟 and 𝑢 𝑐𝑟 . This, as opposed to veering, is the phenomenon which causes the jumps in Ω𝑐𝑟 and 𝑢 𝑐𝑟 of Fig.5.2. 96 in critical flow velocity 𝑢 𝑐𝑟 and critical frequency Ω𝑐𝑟 while the mode of flutter remains constant at 𝑘 = 2. The cause of this behavior is shown in Fig.5.4 which is the Argand diagram for the cantilever beam with P = 0 for 𝛽 = {0.3, 0.325, 0.35, 0.375, 0.4}. This figure shows how the curvature of the locus corresponding to 𝑘 = 2 mode of 𝛽 = 0.3 just misses the Re[Ω] axis before crossing it at a significantly later point with respect to both Ω𝑐𝑟 and 𝑢 𝑐𝑟 . As 𝛽 continues to increase, both Ω𝑐𝑟 and 𝑢 𝑐𝑟 increase, as can also be seen in Fig.5.2. However, included in Fig.5.4 are the loci corresponding to 𝛽 = 0.4, which is when, while both Ω𝑐𝑟 and 𝑢 𝑐𝑟 continue to increase slightly, the mode of flutter changes from 𝑘 = 2 to 𝑘 = 3 due to the phenomenon of veering. Veering is when two loci approach each other in the complex plane but do not intersect before continuing in divergent directions [110, 111], at some points often dramatically swapping direction towards that previously taken by the other mode. Even with the dramatic change in the behavior of the loci, the values of Ω𝑐𝑟 and 𝑢 𝑐𝑟 in Fig.5.2 do not signify that anything in particular has happened. Looking at Fig.5.5, we have the Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.4, 0.45, 0.5, 0.525, 0.55}. For this range, the flow velocity required for the system to lose stability 𝑢 𝑐𝑟 continues to increase gradually with increasing 𝛽, but the critical frequency Ω𝑐𝑟 remains nearly unchanged. At 𝛽 = 0.55, we again see the veering phenomenon that causes a change in the mode of flutter, but this time from 𝑘 = 3 back to 𝑘 = 2, all while the values of Ω𝑐𝑟 and 𝑢 𝑐𝑟 in Fig.5.2 continue to hold steady or change as they had been. For Fig.5.6, we have the Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.55, 0.575, 0.6, 0.625}. Very similarly to Fig.5.5, the critical flow velocity 𝑢 𝑐𝑟 continues to increase gradually with increasing 𝛽 while the critical frequency Ω𝑐𝑟 stays nearly constant until the system again experiences veering at 𝛽 = 0.625, this time from 𝑘 = 2 to 𝑘 = 1. For Figs.5.7 - 5.9, the mode of flutter remains constant at 𝑘 = 1, but as Fig.5.2 shows, there are jumps in Ω𝑐𝑟 and 𝑢 𝑐𝑟 at 𝛽 = 0.7 and 𝛽 = 0.925. In Fig.5.7, we have the Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.625, 0.65, 0.675, 0.7}. We see the loci of the first mode exhibits several oscillations before contacting the Re[Ω] axis at nearly the same Ω𝑐𝑟 , but shallower and shallower angles before the loci corresponding to 𝛽 = 0.7 passes right above the Re[Ω] axis 97 25 20 15 𝛽↑    10 𝛽↑ 𝛽↑ 5 𝛽↑ 0 10 20 30 40 50 60 k=1 k=2   k=3 Figure 5.5: Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.4, 0.45, 0.5, 0.525, 0.55}. Loci of increasing 𝛽 are denoted as well as the direction in which each loci moves as 𝑢 𝑒 increases. All but the loci corresponding to 𝛽 = 0.55 are the solid lines. The loci of 𝛽 = 0.55 is included to show the veering effect as explained in Fig.5.4: The mode changes back to 𝑘 = 2 from 𝑘 = 3 though Ω𝑐𝑟 and 𝑢 𝑐𝑟 do not exhibit a jump in value. to cross at a much higher Ω𝑐𝑟 with a correspondingly higher 𝑢 𝑐𝑟 . This is the same behavior that marked the first jump of Ω𝑐𝑟 and 𝑢 𝑐𝑟 Fig.5.2 at 𝛽 = 0.3. In Fig.5.8, which shows the Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.7, 0.775, 0.825, 0.9, 0.925}, we again see the same behavior with the loci corresponding to 𝛽 = 0.925 skimming right above the Re[Ω] axis to cross at a much higher Ω𝑐𝑟 and higher 𝑢 𝑐𝑟 . Finally with Fig.5.9, we have the Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.925, 0.95, 0.975, 0.995}. These loci all exhibit similar behavior to each other but noticeably cross at a significantly higher Ω𝑐𝑟 with a higher 𝑢 𝑐𝑟 . 98 30 25 20  𝛽↑  15  𝛽↑ 10 𝛽↑ 𝛽↑ 5 𝛽↑ 0 10 20 30 40 50 60 k=1 k=2 e[  k=3 Figure 5.6: Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.55, 0.575, 0.6, 0.625}. Loci of increasing 𝛽 are denoted as well as the direction in which each loci moves as 𝑢 𝑒 increases. All but the loci corresponding to 𝛽 = 0.625 are the solid lines. In this set, the loci of 𝛽 = 0.625 is included to show the veering effect as explained in Figs.5.4 and 5.5. In this case however, the veering occurs between the loci of the 𝑘 = 1 and 𝑘 = 2 modes. 99 30 25 20   15 𝛽↑  𝛽↑ 10 𝛽↑ 5 0 10 20 30 40 50 60 k=1 k=2  ! " k=3 Figure 5.7: Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.625, 0.65, 0.675, 0.7}. Loci of increasing 𝛽 are denoted as well as the direction in which each loci moves as 𝑢 𝑒 increases. All but the loci corresponding to 𝛽 = 0.7 are the solid lines. The loci of 𝛽 = 0.7 is included to show how the loci of the 𝑘 = 1 mode transitions from crossing the Re[Ω] axis at around Re[Ω]= 27 to just missing the Re[Ω] axis before crossing at a much higher value of Ω𝑐𝑟 and 𝑢 𝑐𝑟 . This corresponds to the second jump of Ω𝑐𝑟 and 𝑢 𝑐𝑟 in Fig.5.2, similar to how the behavior shown in Fig.5.4 corresponds to the first set of jumps. 100 45 40 35 𝛽↑ 30 & 25 % $# 20 15 10 𝛽↑ 𝛽↑ 5 𝛽↑ 0 10 20 30 40 50 60 k=1 k=2 '() * k=3 Figure 5.8: Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.7, 0.775, 0.825, 0.9, 0.925}. Loci of increasing 𝛽 are denoted as well as the direction in which each loci moves as 𝑢 𝑒 increases. All but the loci corresponding to 𝛽 = 0.925 are the solid lines. The loci of 𝛽 = 0.925 is included to show how the loci of the 𝑘 = 1 mode transitions from crossing the Re[Ω] axis at around Re[Ω]= 45 to just missing the Re[Ω] axis before crossing at a much higher value of Ω𝑐𝑟 and 𝑢 𝑐𝑟 . This corresponds to the third jump of Ω𝑐𝑟 and 𝑢 𝑐𝑟 in Fig.5.2, similar to how the behavior shown in Figs.5.4 and 5.7 correspond to the first and second sets of jumps. 101 45 40 35 𝛽↑ 30 . 25 - ,+ 20 𝛽↑ 15 10 𝛽↑ 𝛽↑ 5 𝛽↑ 0 10 20 30 40 50 60 70 k=1 k=2 /01 2 k=3 Figure 5.9: Argand diagrams for the cantilever beam with P = 0 for 𝛽 = {0.925, 0.95, 0.975, 0.995}. Loci of increasing 𝛽 are denoted as well as the direction in which each loci moves as 𝑢 𝑒 increases. For such high values of 𝛽, the system loses stability at a very high Ω𝑐𝑟 though requiring a correspondingly high 𝑢 𝑐𝑟 to induce the loss of stability. 102 Numerical Procedure For each point in the 𝛽 − 𝑢 e domain3 defined in Section 5.3, we solve for the critical follower force P = Pc𝑟 , which causes the system to lose stability through flutter. These points define a surface, which we refer to as the critical follower force surface. Each point on the critical follower force surface corresponds to a critical frequency Ωc𝑟 ; these points define a critical frequency surface. The critical stability and frequency surfaces are obtained as follows: For a specific value of 𝛽, we start with P = 0 and 𝑢 e = 0.01. We use the first three natural frequencies of the beam, Ω 𝑗 , 𝑗 = 1, 2, 3, as the initial guesses to solve for the eigenfrequencies as the magnitude of P is gradually increased. The process is continued until one of the Ω 𝑗 terms satisfies the condition I𝑚[Ω 𝑘 ] = 0. This provides the value of Pc𝑟 and Ωc𝑟 for 𝑢 e = 0.01 and the specific value of 𝛽. Here, in keeping with Section 5.3, the value of 𝑘 denotes the mode of the unforced locus from which the forced locus which lost stability originated. The process is repeated by gradually incrementing the value of 𝑢 e while keeping the value of 𝛽 fixed; the process is terminated when the value of 𝑢 e = 𝑢 c𝑟 (𝛽)4. To obtain the critical stability and frequency surfaces, the overall process is repeated on a fine mesh grid for 𝛽. The procedure for computing the critical stability points, described in Section 5.3, can now be better explained with the help of the Argand diagrams in Figs.5.10 and 5.11. The thin lines represent the loci of increasing fluid flow 𝑢 e (𝛽 = 0.75) and no follower force as previously shown in Fig.5.8. For both figures, given the shared value of 𝛽 = 0.75 the zero follower force P = 0 loci are the same, and, as can be seen from Fig.5.8, eventually lose stability through the 𝑘 = 1 mode at Ω𝑐𝑟 = 26.08 and 𝑢 𝑐𝑟 = 8.78. At a specific point on each loci of the unforced Argand diagram corresponding to a specific value of 𝑢 e = 𝑢 ∗e , the follower force P is increased from P = 0. Increasing the value of P from zero results in a separate locus sprouting from each of the branches. These continue until the critical value of P = Pc𝑟 is reached when one of the eigenfrequencies satisfy Im[Ω] = 0, analogous to the P = 0 system. 3We have 0 ≤ 𝑢 e ≤ 𝑢 c𝑟 (𝛽) defining the range of 𝑢 e for each 𝛽 4This signifies that the fluid flow alone causes the beam to lose stability, much like a flag fluttering in the wind. 103 40 35 30 25 : 9 20 87 15 𝑢 ∗𝑒 = 4.0 𝑢 ∗𝑒 = 4.0 10 5 𝑢 ∗𝑒 = 4.0 0 0 10 20 30 40 50 60 70 k=1 k=2 345 6 k=3 Figure 5.10: Argand diagrams for the cantilever beam with 𝛽 = 0.75. The thin lines represent the loci of increasing fluid flow 𝑢 e and no follower force. The thicker lines represent the loci for the system with 𝑢 ∗e = 4.0 as the follower force P = 0 increases until the system loses stability when one of the loci crosses the Re[Ω] axis. Here, while the system with no follower force is stable at 𝑢 ∗e = 4.0 (though it eventually loses stability through the 𝑘 = 1 mode at Ω𝑐𝑟 = 44.58 and 𝑢 𝑐𝑟 = 13.14), the follower force causes the system to lose stability in the 𝑘 = 2 mode at Ω𝑐𝑟 = 13.80 and P𝑐𝑟 = 27.84. These Argand diagrams show how the behavior that results from increasing the follower force while holding the flow velocity constant differs greatly from the P = 0 loci. In Fig.5.10 the flow velocity is held at 𝑢 ∗e = 4.0, a point at which the system is initially stable. Increasing the follower force results in the new loci quickly branching from the P = 0 loci with the 𝑘 = 1 locus remaining on the Im[Ω] axis while moving away from the Re[Ω] axis, the 𝑘 = 2 locus heading quickly for the Re[Ω] axis, and the 𝑘 = 3 locus maintaining its Im[Ω] component while decreasing in Re[Ω]. Eventually, the system loses stability through flutter in the 𝑘 = 2 mode at Ω𝑐𝑟 = 13.80 104 and P𝑐𝑟 = 27.84, noticeably different than the 𝑘 = 1 mode through which the P = 0 system lost stability through flutter. In Fig.5.11 the flow velocity is held at 𝑢 ∗e = 7.0, from which the system is initially stable. Increasing the follower force results in the 𝑘 = 2 locus remaining nearly in place while performing a small loop. The loci of the 𝑘 = 1 and 𝑘 = 3 modes move nearly directly towards each other before veering apart with the 𝑘 = 1 locus moving away from the Re[Ω] axis while the 𝑘 = 1 locus moving towards the Re[Ω] axis. Eventually, the follower force causes the system to lose stability through flutter in the 𝑘 = 3 mode at Ω𝑐𝑟 = 26.44 and P𝑐𝑟 = 32.64. Of note, while the locus of increasing P of the 𝑘 = 1 mode begins closer to the Re[Ω] axis than that of the 𝑘 = 3 mode, the follower force imposes a stabilizing behavior on the locus of the 𝑘 = 1 mode, while that is not the case for the locus of the 𝑘 = 3 mode. The system does, however, still loses stability through flutter nonetheless. 105 40 35 30 25 B A 20 @? 𝑢 ∗𝑒 = 7.0 15 𝑢 ∗𝑒 = 7.0 10 5 𝑢 ∗𝑒 = 7.0 0 0 10 20 30 40 50 60 70 k=1 k=2 ;<= > k=3 Figure 5.11: Argand diagrams for the cantilever beam with 𝛽 = 0.75. The thin lines represent the loci of increasing fluid flow 𝑢 e and no follower force. The thicker lines represent the loci for the system with 𝑢 ∗e = 7.0 as the follower force P = 0 increases until the system loses stability when one of the loci crosses the Re[Ω] axis. Here, while the system with no follower force is stable at 𝑢 ∗e = 7.0 (though it eventually loses stability through the 𝑘 = 1 mode at Ω𝑐𝑟 = 44.58 and 𝑢 𝑐𝑟 = 13.14 - see Fig.5.8), the follower force causes the system to lose stability in the 𝑘 = 1 mode at Ω𝑐𝑟 = 26.44 and P𝑐𝑟 = 32.64. Additionally, we also see how increasing P can also cause the system to exhibit the phenomenon of veering as shown by the loci originating from the 𝑘 = 1 and 𝑘 = 3 modes. It is also interesting to note how, while the locus of increasing P of the 𝑘 = 1 mode begins closer to the Re[Ω] axis than that of the 𝑘 = 3 mode, the follower force imposes a stabilizing behavior on the locus of the 𝑘 = 1 mode, while that is not the case for the locus of the 𝑘 = 3 mode. The behavior of the locus of the 𝑘 = 2 mode, while not consequential is interesting nonetheless. 106 5.4 Results Critical Stability Surfaces Figure 5.12: Critical follower force P𝑐𝑟 surface. The 𝑢 𝑒 = 0.01 border is at a nearly constant P𝑐𝑟 = 17.557, from which P𝑐𝑟 increases or decreases for different values of 𝛽 as 𝑢 𝑒 increases. The entirety of the surface terminates at P𝑐𝑟 = 0 at 𝑢 𝑒 = 𝑢 𝑐𝑟 as this is where the system loses stability without any P as discussed in Section 5.3. For lower values of 𝛽, the surface slopes only downwards, while for higher values of 𝛽, the surface exhibits a series of ridges of greater P𝑐𝑟 emanating from the discontinuities also described in Section 5.3. To begin, we have the critical follower force P𝑐𝑟 surface shown in Fig.5.12. This figure shows how the value of P𝑐𝑟 changes for each 𝛽 as 𝑢 𝑒 increases. As will be discussed further in Section 5.4, at 𝑢 𝑒 = 0.01, for each 𝛽, the critical follower force is a nearly constant P𝑐𝑟 ≈ 17.557. As 𝑢 𝑒 increases, the entire surface eventually decreases to P𝑐𝑟 = 0 at 𝑢 𝑒 = 𝑢 𝑐𝑟 which corresponds to when the system loses stability through fluid flow alone. However, for higher values of 𝛽, there are a series of ridges of significantly higher P𝑐𝑟 seemingly emanating from each of the discontinuities. 107 50 45 40 35 30 Cc 𝛽↑ P 25 20 15 10 5 0 0 2 4 6 8 10 12 14 16 18 20 ue Figure 5.13: Curves of critical follower force P𝑐𝑟 with order of increasing 𝛽 denoted. Presented here is the same information as the surface of Fig.5.12 but as individual lines. Of particular note is that for a given 𝑢 𝑒 , the follower force required to cause the system to lose stability through flutter increases monotonically with increasing 𝛽. This view of the data really signifies the tremendous change in P𝑐𝑟 needed to cause the system to lose stability through flutter as 𝛽 and 𝑢 𝑒 vary. As 𝛽 is increased, the crests of these ridges curve backwards compared to 𝑢 𝑒 . The information provided in Fig.5.12 is also shown in Fig.5.13 as a set of distinct curves for each 𝛽 given in Section 5.3. The ridges are very pronounced from this view. For 𝛽 < 0.3, P𝑐𝑟 shows a monotonic decrease to zero as 𝑢 𝑒 increases. This is the region that corresponds to Fig.5.3, before the first discontinuity. For values of 𝛽 in 0.3 < 𝛽 < 0.475, the P𝑐𝑟 curves decrease in value before rising over a ridge of higher P𝑐𝑟 emanating from the first discontinuity, and then going down to P𝑐𝑟 = 0. For 𝛽 ≥ 0.475 the P𝑐𝑟 curves increase in value before cresting the ridge emanating from first discontinuity and then going down towards P𝑐𝑟 = 0. Of those, the P𝑐𝑟 curves 108 for 𝛽 ≤ 0.6750 continue falling to zero, while the P𝑐𝑟 curves for 𝛽 ≥ 0.7 (having made it past the second discontinuity), rise up another ridge of higher P𝑐𝑟 emanating from that second discontinuity. Finally, of the curves that make it past the second discontinuity, for values of 𝛽 ≤ 0.8750, the P𝑐𝑟 curves continue falling to P𝑐𝑟 = 0 while the P𝑐𝑟 curves for 𝛽 ≥ 0.9 (having made it past the third discontinuity), rise up a final ridge emanating from that third discontinuity. Of note is that for a given 𝑢 𝑒 , the follower force required to cause the system to lose stability through flutter increases monotonically with increasing 𝛽, i.e the P𝑐𝑟 curves never cross. Figure 5.14: Critical frequency Ω𝑐𝑟 surface. There are four large plateaus of nearly constant critical frequency Ω𝑐𝑟 corresponding to the plateaus of similar critical frequency Ω𝑐𝑟 of the unforced case in Fig.5.2. Rather than abrupt jumps in Ω𝑐𝑟 , these plateaus are separated by smooth transitions. Overall, this surface shows a significantly less varied dependence on 𝛽 and 𝑢 𝑒 than P𝑐𝑟 does as seen in Fig.5.12. Following this, Fig.5.14 shows the critical frequency Ω𝑐𝑟 surface at which the increasing follower force P𝑐𝑟 causes the system to lose stability for each 𝛽 and 𝑢 𝑒 . The first obvious 109 observation is that there are large plateaus of very similar critical frequency Ω𝑐𝑟 separated by smooth transitions which occur along slopes emanating from the discontinuities and corresponding to the ridges of P𝑐𝑟 seen in Figs.5.12 and 5.13. To reiterate, rather than having discontinuities, the follower force causes the system to exhibit smooth transitions between frequency plateaus. It is also noticeable how each of the plateaus of similar critical frequency Ω𝑐𝑟 very closely matches the frequency plateaus of Fig.5.2 which represents the unforced case with P = 0. Figure 5.15: The critical mode surface which shows from which unforced P = 0 loci, the forced loci which lost stability originated. Border corresponds to where the system loses stability without any Pval as discussed in Section 5.3. For the portion of lower 𝑢 𝑒 , the √ system lost stability through the 𝑘 = 2 mode which corresponds to where the damping term 𝑢 𝑒 𝛽 is relatively small. There is a curve stripe of 𝑘 = 3 mode in the middle of the space while the high 𝑢 𝑒 - high 𝛽 region is 𝑘 = 1 mode. Finally, we have the critical mode surface in Fig.5.15 which shows from which unforced P = 0 loci, the forced loci which lost stability originated. For the portion of lower 𝑢 𝑒 , the system lost √ stability through the 𝑘 = 2 mode which corresponds to where the damping term 𝑢 𝑒 𝛽 is relatively 110 small. There is a curve stripe of 𝑘 = 3 mode in the middle of the space while the high 𝑢 𝑒 - high 𝛽 region is 𝑘 = 1 mode. Of interest are the regions where the mode changes. As neither the critical follower force surface in Fig.5.12 nor the critical frequency surface in Fig.5.14 exhibit discontinuous jumps of any kind, the behavior of the critical mode surface indicates that each of these changes in critical modes is due to the phenomenon of veering as described in Section 5.3. √ Scaling by Damping Term 𝑢 𝑒 𝛽 While the discussion so far has been with respect to separate 𝛽 and 𝑢 𝑒 axes, interesting behavior √ becomes apparent when we scale our 𝑢 𝑒 axis by 𝛽 such that the x-axis is now exactly proportional to the coeffiecient of the damping term of (5.3). This leaves only 𝑢 2𝑒 + P as the other free parameter of the system. Figures 5.16(a) and (b) show how the system behavior collapses to much simpler √ surfaces when the x-axis is 𝑢 𝑒 𝛽. Fig.5.16(a) shows how the critical frequency Ω𝑐𝑟 at which √ √ the system loses stability is constant for any given 𝑢 𝑒 𝛽. For a given 𝑢 𝑒 𝛽, it is the sum of the forcing terms that uniquely define the state at which the system loses stability. Fig.5.16(b) shows how the critical mode 𝑘 𝑐𝑟 from which the system loses stability is also defined almost entirely by 𝛽 Ω𝑐𝑟 𝛽 𝑘 𝑐𝑟 (a) √ (b) √ 𝑢𝑒 𝛽 𝑢𝑒 𝛽 Figure √ 5.16: (a) Critical Frequency surface as presented in Fig.5.14, but with the x-axis scaled by 𝛽. If viewed from the XZ plane, there would only √ be a line visible. (b) Critical Mode surface as presented in Fig.5.15, but with the x-axis scaled by 𝛽. 111 √ 𝑢 𝑒 𝛽. There are, however, interesting deviations from this generalization: the 𝑘 = 2 mode stripe √ at 𝑢 𝑒 𝛽 ≈ 5 and a small region at the base of the 𝑘 = 3 mode stripe. This is in large part due to the veering phenomenon in which changing certain parameters slightly can cause the mode behavior to change, but the final state of the system at the point of instability to be largely unchanged. Ziegler’s Paradox As can be seen in Fig.5.13, the critical value of P which causes the system to lose stability through flutter for very small 𝑢 𝑒 , particularly 𝑢 𝑒 = 0.01, is P𝑐𝑟 ≈ 17.557, regardless of 𝛽 > 0. This, however, is significantly lower than the system without any damping at all as shown in Figs.5.17 and 5.18. In the 𝑢 𝑒 = 𝛽 = 0 (undamped) case5, the loci of the 𝑘 = 1 and 𝑘 = 2 modes move towards each other along the Re[Ω] axis until they meet at Re[Ω] = 11.0828 and separate off the axis as a Hamiltonian Hopf bifurcation [91]. One locus then moves into the Im[Ω] > 0 (stable) half-plane while the other locus moves into the Im[Ω] < 0 (unstable) half-plane, causing the system to lose stability through flutter. Fig.5.18 shows how Re[Ω 𝑘 ] changes with increasing P until the loci of the 𝑘 = 1 and 𝑘 = 2 modes meet at P𝑐𝑟 = 20.06 signifying the moment of flutter. Conversely, Figs.5.19 and 5.20 show how the system evolves with increasing P when 𝑢 𝑒 = 0.01 and 𝛽 = 0.001 are held constant (signifying very small damping). In this case, the loci of the 𝑘 = 1 and 𝑘 = 2 modes are very slightly above the Re[Ω] axis as they approach each other, in accordance with the damping having a stabilizing effect. However, as can be seen in Fig.5.19, the loci split before meeting. This behavior is shown in much greater detail in Fig.5.20. The loci of the damped system (solid) each have Im[Ω 𝑘 ] > 0 (stable) though they are very slowly diverging (in a Im[Ω 𝑘 ] sense) from each other as P increases. Eventually, the lowermost locus passes the Im[Ω] = 0 axis at P𝑐𝑟 = 17.56, which is a significantly lower value than the undamped case (denoted by the dashed lines). This destabilizing effect of damping is in accord with the Ziegler paradox as described by Ziegler in his seminal 1952 paper [62]. 5Both having the effect of making the third (damping/coriolis) term of (5.1) go to 0, while 𝑢 𝑒 = 0 entails no contribution to the second (forcing) term. 112 5 4 3 2 1 T S 0 QO HN HM HL HK HJ 0 10 20 30 40 50 60 70 DEF G Figure 5.17: Argand diagrams for the cantilever beam subject to increasing follower force P, without fluid flow or damping. Without damping, the roots of the system remain entirely real until the loci of the 𝑘 = 1 and 𝑘 = 2 modes meet at Re[Ω]= 11.0828 causing the system to lose stability through coupled-mode flutter by the mechanism of the Hamiltonian Hopf bifurcation [91]. With this in mind, it is interesting to take another look at Fig.5.13 and see how, as 𝑢 𝑒 increases, different values of 𝛽 have the system require increasing or decreasing amounts of P to destabilize the system. Of particularly peculiar note, is that for a significant amount of 𝑢 𝑒 , 𝛽 = 0.475 requires nearly constant P to lose stability. Curves of P𝑐𝑟 for 𝛽 < 0.475 decrease quadratically with increasing 𝑢 𝑒 while curves of P𝑐𝑟 for 𝛽 > 0.475 increase quadratically with increasing 𝑢 𝑒 . 113 70 60 50 40 Re[ ] 30 20 10 0 0 5 10 15 20 25 P Figure 5.18: Re[Ω] for the cantilever beam subject to increasing follower force P, without fluid flow or damping. This figure clearly shows how the real part of the natural frequencies of the 𝑘 = 1 and 𝑘 = 2 modes approach each other before meeting at Re[Ω] = 11.0828 when P = 20.06 and becoming complex conjugate pairs, signified by their both having the same real part. 114 0.6 0.4 0.2 f d 0 ba YZ\` YZ\_ YZ\^ 0 10 20 30 40 50 60 70 UVW X Figure 5.19: Argand diagrams for the cantilever beam subject to increasing follower force P, with fluid flow 𝑢 𝑒 = 0.01 and 𝛽 = 0.001. With this very small damping, the roots of the system travel just above the Re[Ω] axis until they very gradually veer apart before the point of meeting. The thin lines show how the system evolves for increasing 𝑢 𝑒 with P = 0 for comparison. Compared to the case of holding 𝑢 𝑒 constant while increasing P, this case experiences greater and greater damping. 115 10qr 5 4 3 2 1 p o 0 mn s 17.56 gl t -1.81073e-07 gk gj gi gh 15 16 17 18 19 20 21 22 P𝑐𝑟 Figure 5.20: Comparison of Im[Ω] for the cantilever beam subjected to increasing follower force P, both without fluid flow or damping (dashed) and also with 𝑢 𝑒 = 0.01 and 𝛽 = 0.001 (solid). This figure clearly shows how the imaginary part of the very lightly damped system crossed the Re[Ω] axis at P𝑐𝑟 = 17.56, well before the undamped system. 116 Comparison of Follower Force to Effluent Jet with Nozzle Attachment In the discussion so far, the value of 𝛽 has generally corresponded to the ratio of the added mass of the surrounding fluid to the mass of the beam. However, the system described by (5.1) was first approached using the same equation as a fluid conveying pipe whereby the fluid mass was not external, i.e. 𝑚 𝑒 , but internal, i.e. 𝑚𝑖 (with internal flow velocity 𝑈𝑖 analogous to external flow velocity 𝑈𝑒 ). As an extension of this (and for experimental convenience), Gregory and Païdoussis [55] analyzed what the effects would be of introducing a massless convergent nozzle to the downstream/effluent end of the fluid conveying pipe. As part of their work, they presented the variation of (5.1) with a convergent nozzle as 𝜕 4 𝑦(𝑥, 𝑡) 𝜕 2 𝑦(𝑥, 𝑡) 𝜕 2 𝑦(𝑥, 𝑡) 𝜕 2 𝑦(𝑥, 𝑡) 𝐸𝐼 + 𝑚 𝑈 𝑈 𝑖 i j + 2𝑚 𝑈 𝑖 i + (𝑚 i + 𝑚 b ) =0 (5.15) 𝜕𝑥 4 𝜕𝑥 2 𝜕𝑥𝜕𝑡 𝜕𝑡 2 where 𝐴𝑖 𝑈j = 𝑈i = 𝑈i 𝛼 𝑗 𝐴𝑗 and 𝛼 𝑗 = 𝐴𝑖 /𝐴 𝑗 is the ratio of the pipe cross-section to the cross-sectional area of the nozzle exit. Following the same process as before, this takes the nondimesnional form 𝜕4𝑣 𝜕2𝑣 p 𝜕2𝑣 𝜕2𝑣 + 𝑢 i 𝑢 j + 2 𝛽𝑖 𝑢 i + =0 (5.16) 𝜕𝑢 4 𝜕𝑢 2 𝜕𝑢𝜕𝜏 𝜕𝜏 2 where 𝑚i 𝛽𝑖 = 𝑚i + 𝑚b Using this, we can further analyze how the follower force 𝑃 would compare to different nozzle ratios. To do this, in view of the structure of (5.3), we further introduce P 𝛼𝑃 = 1 + 2 (5.17) 𝑢𝑒 117 Figure 5.21: Contour plots of pipe area to nozzle exit area ratio 𝛼 𝑗 (overlain on the full 3D surface) that would be equivalent to the effects of the follower force P which causes the system to lose stability to fluid flow 𝑢𝑖 . Values for 𝑢𝑖 < 2 are omitted because 𝛼 𝑗 becomes excessively large as 𝑢𝑖 → 0. The rightmost border of the surface is defined by 𝛼 𝑗 = 1 which corresponds to the system losing stability without any follower force, reverting to the original equation with only 𝑢 2𝑒 . which allows us to directly equate 𝛼𝑃 to 𝛼 𝑗 to see what nozzle ratio corresponds to which follower 2 force. These two terms, 𝛼𝑃 and 𝛼 𝑗 , are the coefficients multiplying the second term 𝑢 2e,𝑖 𝜕 2𝑣 of 𝜕𝑢 their respective equations. Of note is that neither of these two coefficients have an effect on the damping of the system which appears in the third term of each equation. Figure 5.21 shows the surface of 𝛼𝑃 = 𝛼 𝑗 (hereafter, just 𝛼 𝑗 ) as a function of 𝛽 and 𝑢 𝑒 = 𝑢𝑖 (hereafter, just 𝑢 𝑒 ). To begin, we omit the portion of the surface for 𝑢 𝑒 < 2 as the values of 𝛼 𝑗 become extremely large as 𝑢 𝑒 → 0 due to the division by 𝑢 2𝑒 . The physical meaning of this is that the nozzle would need to have a very fine opening relative to the nominal diameter of the pipe so as to generate a very high velocity of the effluent fluid. The dominance of 𝑢 2𝑒 continues over the whole 118 surface wherein the peaks of P as seen in Fig.5.12 are visible, but the division by 𝑢 2𝑒 results in a nearly monotonic decrease of 𝛼 𝑗 until the system loses stability through flutter without any follower force (𝛼 𝑗 = 1). However, it may be observed that there are regions where the value of 𝛼 𝑗 shows an increase with increasing 𝑢 𝑒 (signified by a contour line curving downward). These regions are always associated with the discontinuities of 𝑢 𝑐𝑟 (𝛽) wherein a decrease in 𝛽 post-discontinuity, would result in a continuously decreasing 𝑢 𝑐𝑟 . The regions of decreasing 𝛼 𝑗 correspond to places where the “relative" 𝛽 is smaller than the unforced (P = 0) system with equivalent 𝑢 𝑒 . On a related note, while it may seem obvious, it is significant that the 𝛼 𝑗 surface does not go below 𝛼 𝑗 = 1. This is in accord with the findings that a tensile force is always stabilizing in addition to analysis and experiments in which a diverging nozzle is attached to the end of the pipe which results in no axial effluence, and consequently, no loss of stability [113]. An alternative interpretation of the dips of the 𝛼 𝑗 contours of Fig.5.21, particularly the 𝛼 𝑗 = 1.1 contour and the implied 𝛼 𝑗 = 1 contour, is that a negative P follower force could re-stabilize the system through the discontinuities of 𝑢 𝑐𝑟 (𝛽). This implies that the 𝛼 𝑗 , and by inference, the P𝑐𝑟 , surfaces could be extended towards greater 𝑢 𝑒 by accepting P < 0 and 𝛼 𝑗 < 1. 5.5 Conclusion There are several methods of inducing flutter in non-conservative systems including due to a follower force or through interaction of the beam with a fluid, flowing internally or externally. In this work, we explore the interaction that develops from combining these two methods. We analyze a flexible, tail-like beam in axial fluid flow with a tangential follower force at the free end. The equation of motion and boundary conditions are defined, and the motion is determined using a set of parameters. Natural frequencies are explored without the follower force to establish a basis for our analysis. Critical follower forces are analyzed to determine the frequency and mode at which the system loses stability through flutter, given each combination of mass ratio and fluid flow velocity. Increasing the mass ratio monotonically requires a larger follower force to cause instability. Scaling the flow velocity axis by the square root of the mass ratio allows all results to collapse onto a single curve. Results confirm Ziegler’s paradox at very low flow velocities. 119 Finally, we investigate how incorporating a convergent nozzle to the downstream end of the system shows that the critical follower force can be simulated for higher flow velocities over a very wide range of mass ratios (though not feasibly for lower flow velocities). In conjunction with the work on fluid-conveying, fluid-immersed tails introduced by Hellum [52], this convergent nozzle would allow for the more controllable instigation of flutter. Utilizing both the effluent jet thrust as well as the oscillation of the fluttering tail could be used as a novel propulsive device for underwater vehicles. 120 CHAPTER 6 CONCLUSION AND FUTURE DIRECTIONS 6.1 Forced Excitation of a Flexible Beam in Fluid Flow In Chapter 2, we developed an analytical model for a rigid body in flow and demonstrated that it was possible to solve it numerically across a large domain of model and state parameters. Using slender body equations for thrust and power, we showed how they could be used to analyze a large region of the frequency-flow rate space. We observed that the flexible propulsor exhibited both positive and negative thrust and discussed how resonance affected the system. To find the equilibrium velocity of a submersible propelled by a tail-like beam, we defined a drag model of the rigid body. Finally, we found that the efficiency surpassed 50% on this dynamic equilibrium locus. 6.2 Feedback-Induced Flutter Instability of a Flexible Beam in Fluid Flow In Chapter 3, we modeled a pinned-free, flexible, tail-like beam immersed in fluid flow using fluid, subject to various leading edge boundary conditions. We showed how these non-conservative loading due to feedback caused the system to lose stability through flutter. We studied various actuation-proportional-to-sensing combinations incorporating moment and angle actuation; dis- placement, slope, and curvature sensing; and both positive and negative sign of feedback gain. We determined the flow velocity at which the systems with zero feedback gain lost stability. We defined the criteria for a system losing stability due to feedback from which we, in turn, defined the critical stability surfaces to present the results. We presented the results for six Cases, one for each combination of actuation and sensing over a large range of feedback locations and flow velocities which showed a rich set of stability transitions. We determined that the impact of flow velocity was much more mild than the significant impact of the location of sensing. We further determined the propulsive characteristics of the waveforms generated by the flutter instability through slender body equations. We then showed that while the resulting waveform was composed of four individual traveling waves, we could define a Phase Smoothness Factor metric which could be used to reliably predict the efficiency based on the numerical waveform. Finally, we determined that the propulsive 121 characteristics of the beam did not depend on the combination of actuation and sensing by which flutter was produced; they depended only on the values of the dimensionless fluid velocity and critical frequency, which completely defined the waveform. 6.3 Experimental Validation of Feedback-Induced Flutter Instability of a Flexible Beam in Fluid Flow In Chapter 4, we designed and constructed a system for applying feedback-based actuation to the leading edge of a flexible, tail-like beam in fluid flow proportional to the curvature at a point along its centerline. This involved both the application, sensing, and amplification of strain gauges as well as the control and transmission of torque from a motor to the leading edge of the beam. We designed a control scheme that enabled the input of the feedback gain which we used to determine when the beam lost stability through flutter at different axial flow velocities. We showed that the positive feedback gain behavior was significantly different from the negative feedback behavior for our system. We showed that the behavior of the system was loosely correlated with the flow velocity. We modified the moment-proportional-to-curvature boundary condition in Chapter 3 to include the inertia of the torque application hardware and the time delay that occurred due to the flexibility of the torque application hardware. We showed that the numerical model of the experimental setup lost stability through flutter at frequencies very similar to those observed in the experiment. With this, we could state with confidence that the model used was a good approximation of the real system and that feedback-induced flutter of a flexible beam in axial fluid flow was achievable. 6.4 Effects of a Follower Force on the Inducement of Flutter in a Flexible Beam in Fluid Flow In Chapter 5, we modeled a cantilever-free, flexible, tail-like beam in axial fluid flow with a tangential follower force applied at the free end. We defined the equation of motion and boundary conditions of this system and showed the method of determining the motion given a set of parameters. We explored the natural frequencies of the system without follower force over the full range of possible beam-to-fluid mass ratios as the fluid flow increased. These results presented a rich set of behaviors that defined the points from which our follower-force analysis began. Our 122 analysis of the critical follower force allowed us to determine the frequency and mode at which the system lost stability through flutter given each combination of mass ratio and fluid flow velocity. We showed that for any given flow velocity, increasing the mass ratio always required a larger follower force to cause the system to lose stability. We showed that by scaling the flow velocity axis by the square root of the mass ratio, we could get all of our results to collapse onto a single curve for each of the critical follower force surface, critical mode surface, and critical frequency surface. Our results for very low flow velocities confirmed Ziegler’s paradox regarding the destabilizing nature of low damping on certain systems. Our analysis of incorporating a convergent nozzle to the downstream end of the system showed that the follower force could be feasibly simulated for higher flow velocities (though not lower ones) given the tremendous difficulty of actualizing a tangential follower force. 6.5 Future Directions There are several avenues for future work regarding the forced excitation of flexible beams in fluid flow, particularly with respect to their performance for underwater propulsion. For example, each of the parameters of our model could be varied to study their impact on the propulsive performance of the system. The torque and power requirements at the revolute joint can be further investigated for feasibility in a real-world system. To accommodate torque and/or power limitations of a physical system, more complex trajectories through the velocity-frequency plane could be designed by changing the amplitude and frequency of the revolute joint. An underwater robotic vehicle utilizing this design would validate this model if implemented successfully. With respect to feedback-induced flutter of a flexible beam in fluid flow as a means of generating a traveling wave for propulsion, we envision a number of other directions along which this work can be extended. First, the observation that a waveform dominated by a single traveling wave appears to be more efficient than a mixed waveform is satisfying, but a proof has eluded us to date. It is also likely that analyzing a discrete analogy of the problem, in which flutter of an articulated system is produced through actuation of the base, can provide additional insights because of the reduction to finite dimension. Finally, we believe that the results can be extended to an interesting boundary 123 condition at the free end, that of a hydrodynamically “active” fin, which has mass, dimension, and experiences fluid forces. This conceptually mimics the situation of a fast-swimming thunniform fish, which has a large caudal fin at the end of its tail. While our experimental work validated the heart of feedback-induced flutter of a flexible beam in fluid flow, several lines of inquiry could be pursued in future research. The first would be to replicate the experiment with a more fine-grained incrementing of the feedback gain. As was seen several times in some of the Trials, the system went from very small oscillations to such large oscillations that the cut-off trigger was activated, while other Trials resulted in nice and steady limit cycles. These limit cycles resulted in beam motions that could have the desired potential as propulsive devices. Additionally, the energetics and hydrodynamic performance (i.e. thrust and efficiency) of the beams could be measured to ascertain their feasibility as a means of propelling a submersible. For the work regarding the effects of the follower force on the instigation of flutter, we investi- gated how incorporating a convergent nozzle to the downstream end of the system could be used to simulate the critical follower force for higher flow velocities over a very wide range of mass ratios. This work could be combined with the work on fluid-conveying, fluid-immersed tails introduced by Hellum [52], in which this convergent nozzle would allow for the more controllable instigation of flutter. Utilizing both the effluent jet thrust as well as the oscillation of the fluttering tail could be used as a novel propulsive device for underwater vehicles. With respect to this entire work, all aspects of it could be extended to incorporate non-linear dynamics and analysis, both for the flexible beams, and the fluid-structure interactions. This analysis could then be incorporated into control schema to better operate these flexible beams as propulsive devices of underwater vehicles. 124 BIBLIOGRAPHY [1] J. Yuh, G. Marani, D. R. Blidberg, Applications of marine robotic vehicles, Intelligent Service Robotics 4 (4) (2011) 221–231. [2] R. Bogue, Underwater robots: A review of technologies and applications, Industrial Robot: An International Journal 42 (3) (2015) 186–191. [3] A. Mendez, T. J. Leo, M. A. Herreros, Current state of technology of fuel cell power systems for autonomous underwater vehicles, Energies 7 (7) (2014) 4676–4693. [4] G. Foresti, Visual inspection of sea bottom structures by an autonomous underwater vehicle, IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics) 31 (5) (2001) 691–705. [5] P. Ridao, M. Carreras, D. Ribas, R. Garcia, Visual inspection of hydroelectric dams using an autonomous underwater vehicle, Journal of Field Robotics 27 (6) (2010) 759–778. [6] X. Xiang, B. Jouvencel, O. Parodi, Coordinated Formation Control of Multiple Autonomous Underwater Vehicles for Pipeline Inspection, International Journal of Advanced Robotic Systems 7 (1) (2010) 3. [7] N. Wakita, K. Hirokawa, T. Ichikawa, Y. Yamauchi, Development of Autonomous Underwa- ter Vehicle (AUV) for Exploring Deep Sea Marine Mineral Resources, Tech. rep., Mitsubishi Heavy Industries (Sep. 2010). [8] S.-F. Chen, J.-Z. Yu, Underwater cave search and entry using a robotic fish with embedded vision, in: Proceedings of the 33rd Chinese Control Conference, 2014, pp. 8335–8340. [9] D. D. Sternlicht, J. E. Fernandez, R. Holtzapple, D. P. Kucik, T. C. Montgomery, C. M. Loef- fler, Advanced sonar technologies for autonomous mine countermeasures, in: OCEANS’11 MTS/IEEE KONA, 2011, pp. 1–5. [10] R. W. Button, J. Kamp, T. B. Curtin, J. Dryden, A Survey of Missions for Unmanned Undersea Vehicles, Tech. rep., Defense Technical Information Center (Jan. 2009). [11] N. R. Council, D. o. E. a. P. Sciences, N. S. Board, C. o. A. V. i. S. o. N. Operations, Autonomous Vehicles in Support of Naval Operations, National Academies Press, 2005. [12] C. Stoker, D. Burch, B. Hine, J. Barry, Antarctic undersea exploration using a robotic submarine with a telepresence user interface, IEEE Expert 10 (6) (1995) 14–23. [13] M. Dunbabin, L. Marques, Robots for Environmental Monitoring: Significant Advancements and Applications, IEEE Robotics & Automation Magazine 19 (1) (2012) 24–39. [14] N. Agarwala, Monitoring the Ocean Environment Using Robotic Systems: Advancements, Trends, and Challenges, Marine Technology Society Journal 54 (5) (2020) 42–60. [15] D. S. Barrett, The design of a flexible hull undersea vehicle propelled by an oscillating foil, Ph.D. thesis, Massachusetts Institute of Technology (1994). 125 [16] P. Bandyopadhyay, Trends in biorobotic autonomous undersea vehicles, IEEE Journal of Oceanic Engineering 30 (1) (2005) 109–139. [17] J. L. Tangorra, S. N. Davidson, P. G. Madden, G. V. Lauder, I. W. Hunter, A Biorobotic Pectoral Fin for Autonomous Undersea Vehicles, in: 2006 International Conference of the IEEE Engineering in Medicine and Biology Society, 2006, pp. 2726–2729. [18] C. Eriksen, T. Osse, R. Light, T. Wen, T. Lehman, P. Sabin, J. Ballard, A. Chiodi, Seaglider: A long-range autonomous underwater vehicle for oceanographic research, IEEE Journal of Oceanic Engineering 26 (4) (2001) 424–436. [19] J. Bleck-Neuhaus, Mechanical resonance: 300 years from discovery to the full understanding of its importance (Nov. 2018). arXiv:1811.08353. [20] M. Buchanan, Going into resonance, Nature Physics 15 (3) (2019) 203–203. [21] L. Prandtl, Kipperscheinungen, Ph.D. thesis, Technische Hochschule Munich (1899). [22] R. C. Petrus, Dynamics of fluid-conveying Timoshenko pipes, Book, Texas A&M University (Aug. 2006). [23] D. Hill, A History of Engineering in Classical and Medieval Times, Routledge, London, 2013. [24] M. W. Kehoe, A historical overview of flight flutter testing, in: AGARD Structures and Materials Panel Meeting, Rotterdam, 1995. [25] G. V. Lauder, Swimming hydrodynamics: Ten questions and the technical approaches needed to resolve them, in: G. K. Taylor, M. S. Triantafyllou, C. Tropea (Eds.), Animal Locomotion, Springer, Berlin, Heidelberg, 2010, pp. 3–15. [26] T. Van Buren, D. Floryan, A. J. Smits, Bio-inspired underwater propulsors, in: Bio-Inspired Structures and Design, Cambridge University Press, 2020, pp. 113–137. [27] M. J. Lighthill, Hydromechanics of Aquatic Animal Propulsion, Annual Review of Fluid Mechanics 1 (1) (1969) 413–446. [28] G. V. Lauder, E. J. Anderson, J. Tangorra, P. G. A. Madden, Fish biorobotics: Kinematics and hydrodynamics of self-propulsion, Journal of Experimental Biology 210 (16) (2007) 2767–2780. [29] M. J. Lighthill, Note on the swimming of slender fish, Journal of Fluid Mechanics 9 (2) (1960) 305–317. [30] M. J. Lighthill, Large-amplitude elongated-body theory of fish locomotion, Proceedings of the Royal Society of London. Series B. Biological Sciences 179 (1055) (1971) 125–138. [31] T. Y.-T. Wu, Hydromechanics of swimming propulsion. Part 1. Swimming of a two- dimensional flexible plate at variable forward speeds in an inviscid fluid, Journal of Fluid Mechanics 46 (2) (1971) 337–355. 126 [32] K. W. Moored, D. B. Quinn, Inviscid Scaling Laws of a Self-Propelled Pitching Airfoil, AIAA Journal 57 (9) (2019) 3686–3700. [33] M. M. Koochesfahani, Vortical patterns in the wake of an oscillating airfoil, AIAA Journal 27 (9) (1989) 1200–1205. [34] G. C. Lewin, H. Haj-Hariri, Modelling thrust generation of a two-dimensional heaving airfoil in a viscous flow, Journal of Fluid Mechanics 492 (2003) 339–362. [35] F. S. Hover, Ø. Haugsdal, M. S. Triantafyllou, Effect of angle of attack profiles in flapping foil propulsion, Journal of Fluids and Structures 19 (1) (2004) 37–47. [36] D. A. Read, F. S. Hover, M. S. Triantafyllou, Forces on oscillating foils for propulsion and maneuvering, Journal of Fluids and Structures 17 (1) (2003) 163–183. [37] R. Zhu, J. Wang, H. Dong, D. Quinn, H. Bart-Smith, V. D. Santo, D. Wainwright, G. Lauder, Computational study of fish-shaped panel with simultaneously heaving and bending motion, in: AIAA Scitech 2019 Forum, American Institute of Aeronautics and Astronautics, 2019, p. 1655. [38] H. Park, Y.-J. Park, B. Lee, K.-J. Cho, H. Choi, Vortical structures around a flexible oscillating panel for maximum thrust in a quiescent fluid, Journal of Fluids and Structures 67 (2016) 241–260. [39] P. D. Yeh, Y. Li, A. Alexeev, Efficient swimming using flexible fins with tapered thickness, Physical Review Fluids 2 (10) (2017) 102101. [40] P. A. Dewey, B. M. Boschitsch, K. W. Moored, H. A. Stone, A. J. Smits, Scaling laws for the thrust production of flexible pitching panels, Journal of Fluid Mechanics 732 (2013) 29–46. [41] P. Riggs, A. Bowyer, J. Vincent, Advantages of a Biomimetic Stiffness Profile in Pitching Flexible Fin Propulsion, Journal of Bionic Engineering 7 (2) (2010) 113–119. [42] A. J. Richards, P. Oshkai, Effect of the stiffness, inertia and oscillation kinematics on the thrust generation and efficiency of an oscillating-foil propulsion system, Journal of Fluids and Structures 57 (2015) 357–374. [43] K. A. Morgansen, B. I. Triplett, D. J. Klein, Geometric Methods for Modeling and Control of Free-Swimming Fin-Actuated Underwater Vehicles, IEEE Transactions on Robotics 23 (6) (2007) 1184–1199. [44] M. S. Triantafyllou, G. S. Triantafyllou, An Efficient Swimming Machine, Scientific Amer- ican 272 (3) (1995) 64–70. arXiv:24980373. [45] B. P. Epps, P. Valdivia y Alvarado, K. Youcef-Toumi, A. H. Techet, Swimming performance of a biomimetic compliant fish-like robot, Experiments in Fluids 47 (6) (2009) 927. [46] G. Ozmen Koca, D. Korkmaz, C. Bal, Z. H. Akpolat, M. Ay, Implementations of the route planning scenarios for the autonomous robotic fish with the optimized propulsion mechanism, Measurement 93 (2016) 232–242. 127 [47] B. Lu, C. Zhou, J. Wang, Y. Fu, L. Cheng, M. Tan, Development and Stiffness Optimization for a Flexible-Tail Robotic Fish, IEEE Robotics and Automation Letters 7 (2) (2022) 834– 841. [48] R. Salazar, V. Fuentes, A. Abdelkefi, Classification of biological and bioinspired aquatic systems: A review, Ocean Engineering 148 (2018) 75–114. [49] J. Katz, D. Weihs, Hydrodynamic propulsion by large amplitude oscillation of an airfoil with chordwise flexibility, Journal of Fluid Mechanics 88 (3) (1978) 485–497. [50] A. P. Hoover, R. Cortez, E. D. Tytell, L. J. Fauci, Swimming performance, resonance and shape evolution in heaving flexible panels, Journal of Fluid Mechanics 847 (2018) 386–416. [51] P. Valdivia y Alvarado, K. Youcef-Toumi, Design of Machines With Compliant Bodies for Biomimetic Locomotion in Liquid Environments, Journal of Dynamic Systems, Measure- ment, and Control 128 (1) (2005) 3–13. [52] A. Hellum, R. Mukherjee, A. J. Hull, Flutter instability of a fluid-conveying fluid-immersed pipe affixed to a rigid body, Journal of Fluids and Structures 27 (7) (2011) 1086–1096. [53] M. P. Paidoussis, Dynamics of flexible slender cylinders in axial flow Part 2. Experiments, Journal of Fluid Mechanics 26 (4) (1966) 737–751. [54] M. J. Hannoyer, M. P. Paidoussis, Dynamics of Slender Tapered Beams With Internal or External Axial Flow—Part 1: Theory, Journal of Applied Mechanics 46 (1) (1979) 45–51. [55] R. W. Gregory, M. P. Paidoussis, W. R. Hawthorne, Unstable oscillation of tubular cantilevers conveying fluid II. Experiments, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 293 (1435) (1966) 528–542. [56] M. Argentina, L. Mahadevan, Fluid-flow-induced flutter of a flag, Proceedings of the National Academy of Sciences 102 (6) (2005) 1829–1834. [57] H. Ashley, G. Haviland, Bending Vibrations of a Pipe Line Containing Flowing Fluid, Journal of Applied Mechanics 17 (3) (2021) 229–232. [58] K. Y. Billah, R. H. Scanlan, Resonance, Tacoma Narrows bridge failure, and undergraduate physics textbooks, American Journal of Physics 59 (2) (1991) 118–124. [59] E. Schatzberg, Wings Of Wood, Wings Of Metal: Culture And Technical Choice In American Airplane Materials, 1914-1945, Books by Alumni (Jan. 1999). [60] M. F. Sturkey, Mayday: Accident Reports and Voice Transcripts from Airline Crash Inves- tigations, Heritage Press International, 2005. [61] VI. Reut, About the theory of elastic stability, Proceedings of the Odessa Institute of Civil and Communal Engineering 1 (1939) 126–138. [62] H. Ziegler, The stability criteria of elastomechanics, Archive of Applied Mechanics 20 (1) (1952) 49–56. 128 [63] V. V. Bolotin, Nonconservative Problems of the Theory of Elastic Stability, Macmillan, 1963. [64] R. L. Bisplinghoff, H. Ashley, R. L. Halfman, Aeroelasticity, Courier Corporation, 2013. [65] V. V. Bolotin, Dynamic Instabilities in Mechanics of Structures, Applied Mechanics Reviews 52 (1) (1999) R1–R9. [66] Y. Sugiyama, M. A. Langthjem, K. Katayama, Dynamic Stability of Columns under Noncon- servative Forces: Theory and Experiment, Vol. 255 of Solid Mechanics and Its Applications, Springer International Publishing, Cham, 2019. [67] M. P. Paidoussis, Fluid-Structure Interactions: Slender Structures and Axial Flow, Vol. 1, Academic Press, 2013. [68] T. B. Benjamin, G. K. Batchelor, Dynamics of a system of articulated pipes conveying fluid - I.Theory, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 261 (1307) (1962) 457–486. [69] M. P. Païdoussis, G. X. Li, Pipes Conveying Fluid: A Model Dynamical Problem, Journal of Fluids and Structures 7 (2) (1993) 137–204. [70] A. Cazzolli, F. Dal Corso, D. Bigoni, Non-holonomic constraints inducing flutter instability in structures under conservative loadings, Journal of the Mechanics and Physics of Solids 138 (2020) 103919. [71] A. Cazzolli, F. Dal Corso, D. Bigoni, Flutter Instability and Ziegler Destabilization Paradox for Elastic Rods Subject to Non-Holonomic Constraints, Journal of Applied Mechanics 88 (3) (Dec. 2020). [72] G. Herrmann, R. W. Bungay, On the Stability of Elastic Systems Subjected to Nonconser- vative Forces, Journal of Applied Mechanics 31 (3) (1964) 435–440. [73] I. Elishakoff, Controversy Associated With the So-Called “Follower Forces”: Critical Overview, Applied Mechanics Reviews 58 (2) (2005) 117–142. [74] M. A. Langthjem, Y. Sugiyama, Dynamic Stability of Columns Subjected to Follower Loads: A Survey, Journal of Sound and Vibration 238 (5) (2000) 809–851. [75] M. P. Paidoussis, Dynamics of flexible slender cylinders in axial flow Part 1. Theory, Journal of Fluid Mechanics 26 (4) (1966) 717–736. [76] X. Wang, F. Bloom, Dynamics of a Submerged and Inclined Concentric Pipe System with Internal and External Flows, Journal of Fluids and Structures 13 (4) (1999) 443–460. [77] E. D. Langre, M. P. Païdoussis, O. Doaré, Y. Modarres-Sadeghi, Flutter of long flexible cylinders in axial flow, Journal of Fluid Mechanics 571 (2007) 371–389. [78] Q. Qian, L. Wang, Q. Ni, Vibration and stability of vertical upward-fluid-conveying pipe immersed in rigid cylindrical channel, Acta Mechanica Solida Sinica 21 (5) (2008) 331–340. 129 [79] M. P. Paı¨doussis, T. P. Luu, S. Prabhakar, Dynamics of a long tubular cantilever conveying fluid downwards, which then flows upwards around the cantilever as a confined annular flow, Journal of Fluids and Structures 24 (1) (2008) 111–128. [80] M. Paak, M. P. Païdoussis, A. K. Misra, Nonlinear vibrations of cantilevered circular cylin- drical shells in contact with quiescent fluid, Journal of Fluids and Structures 49 (2014) 283–302. [81] A. Hellum, R. Mukherjee, A. Bénard, A. J. Hull, Modeling and simulation of the dynamics of a submersible propelled by a fluttering fluid-conveying tail, Journal of Fluids and Structures 36 (2013) 83–110. [82] M. Abdullatif, R. Mukherjee, A. Hellum, Critical Stability of a Hinged Beam With Dynamic Moment: With and Without External Flow, in: ASME International Design Engineering Technical Conferences, 2019, p. V008T10A066. [83] G. Herrmann, Dynamics and stability of mechanical systems with follower forces, Tech. Rep. NASA-CR-1782 (Nov. 1971). [84] M. Beck, Die Knicklast des einseitig eingespannten, tangential gedrückten Stabes, Zeitschrift für angewandte Mathematik und Physik ZAMP 3 (3) (1952) 225–228. [85] K.-H. Lin, S. Nemat-Nasser, G. Herrmann, Stability of a Bar Under Eccentric Follower Force, Journal of the Engineering Mechanics Division 93 (4) (1967) 105–115. [86] K. A. McHugh, E. H. Dowell, Nonlinear Response of an Inextensible, Cantilevered Beam Subjected to a Nonconservative Follower Force, Journal of Computational and Nonlinear Dynamics 14 (3) (Jan. 2019). [87] G. T. S. Done, Follower force instability of a pod-mounted jet engine, The Aeronautical Journal 76 (734) (1972) 103–107. [88] Y. Sugiyama, K. Katayama, S. Kinoi, Flutter of Cantilevered Column under Rocket Thrust, Journal of Aerospace Engineering 8 (1) (1995) 9–15. [89] D. Bigoni, G. Noselli, Experimental evidence of flutter and divergence instabilities induced by dry friction, Journal of the Mechanics and Physics of Solids 59 (10) (2011) 2208–2226. [90] D. Bigoni, O. N. Kirillov, D. Misseroni, G. Noselli, M. Tommasini, Flutter and divergence instability in the Pflüger column: Experimental evidence of the Ziegler destabilization paradox, Journal of the Mechanics and Physics of Solids 116 (2018) 99–116. [91] M. P. Paidoussis, Fluid-Structure Interactions, Volume 2: Slender Structures and Axial Flow, Elsevier, 2003. [92] R. H. Plaut, K. Huseyin, Instability of Fluid-Conveying Pipes Under Axial Load, Journal of Applied Mechanics 42 (4) (1975) 889–890. [93] A. Guran, R. H. Plaut, Stability boundaries for fluid-conveying pipes with flexible support under axial load, Archive of Applied Mechanics 64 (7) (1994) 417–422. 130 [94] G. U. O. Chang-qing, L. I. U. Hong-tao, W. Xiao-feng, Z. Chu-han, Vibration and Stability of Pipes Conveying Fluid with Distributed Follower Force, Engineering Mechanics 27 (4) (2010) 190–196. [95] A. Karimi-Nobandegani, S. A. Fazelzadeh, E. Ghavanloo, Effect of Uniformly Distributed Tangential Follower Force on the Stability of Rotating Cantilever Tube Conveying Fluid, Latin American Journal of Solids and Structures 13 (2016) 365–377. [96] M. A. Ilgamov, D. M. Tang, E. H. Dowell, Flutter and Forced Response of a Cantilevered Pipe: The Influence of Internal Pressure and Nozzle Discharge, Journal of Fluids and Structures 8 (2) (1994) 139–156. [97] Y. Xiang, C. M. Wang, S. Kitipornchai, Q. Wang, Dynamic Instability of Nanorods/Nanotubes Subjected to an End Follower Force, Journal of Engineering Mechanics 136 (8) (2010) 1054–1058. [98] R. Bahaadini, M. Hosseini, A. Jamalpoor, Nonlocal and surface effects on the flutter in- stability of cantilevered nanotubes conveying fluid subjected to follower forces, Physica B: Condensed Matter 509 (2017) 55–61. [99] Z.-X. Zhou, O. Koochakianfard, Dynamics of spinning functionally graded Rayleigh tubes subjected to axial and follower forces in varying environmental conditions, The European Physical Journal Plus 137 (1) (2021) 71. [100] B. E. Laithier, M. P. Païdoussis, The equations of motion of initially stressed Timoshenko tubular beams conveying fluid, Journal of Sound and Vibration 79 (2) (1981) 175–195. [101] A. Guran, T. M. Atanackovic, Fluid conveying pipe with shear and compressibility, European Journal of Mechanics - A/Solids 17 (1) (1998) 121–137. [102] M. P. Païdoussis, C. Semler, Non-linear dynamics of a fluid-conveying cantilevered pipe with a small mass attached at the free end, International Journal of Non-Linear Mechanics 33 (1) (1998) 15–32. [103] D. Meng, H.-Y. Guo, S.-P. Xu, Non-linear dynamic model of a fluid-conveying pipe under- going overall motions, Applied Mathematical Modelling 35 (2) (2011) 781–796. [104] Y. H. Lin, Y. K. Tsai, Non-Linear Active Vibration Control of a Cantilever Pipe Conveying Fluid, Journal of Sound and Vibration 202 (4) (1997) 477–490. [105] S. Aspelund, M. Abdullatif, R. Mukherjee, A. M. Hellum, Underwater Propulsion Using Forced Excitation of a Flexible Beam, ASME Journal of Vibration and Acoustics (2023). [106] J. Kutin, I. Bajsić, Fluid-dynamic loading of pipes conveying fluid with a laminar mean-flow velocity profile, Journal of Fluids and Structures 50 (2014) 171–183. [107] C. E. Brennen, A Review of Added Mass and Fluid Inertial Forces., Technical ADA110190, Defense Technical Information Center (1982). [108] F. M. White, Fluid Mechanics, McGraw Hill, 2011. 131 [109] S. Aspelund, R. Mukherjee, A. Hellum, Feedback-Induced Flutter Instability of a Flexible Beam in Fluid Flow, Journal of Sound and Vibration 547 (2023) 117488. [110] C. Pierre, Mode localization and eigenvalue loci veering phenomena in disordered structures, Journal of Sound and Vibration 126 (3) (1988) 485–502. [111] M. Abdullatif, R. Mukherjee, Effect of intermediate support on critical stability of a cantilever with non-conservative loading: Some new results, Journal of Sound and Vibration 485 (2020) 115564. [112] M. Sfakiotakis, D. Lane, J. Davies, Review of fish swimming modes for aquatic locomotion, IEEE Journal of Oceanic Engineering 24 (2) (1999) 237–252. [113] S. Rinaldi, M. P. Païdoussis, Dynamics of a cantilevered pipe discharging fluid, fitted with a stabilizing end-piece, Journal of Fluids and Structures 26 (3) (2010) 517–525. 132