N o n l i n e a r A n a l y s i s o f R i g i d B o d y - V i s c o u s F l o w I n t e r a c t i o n by Pareshkumar Gordhandas Pattani B.Tech., Indian Institute of Technology, Bombay, 1981 M.Eng., The University of Brit ish Columbia, 1983 A Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies Department of Civil Engineering We accept this thesis as conforming to the required standard The University of British Columbia Apr i l 1986 © Pareshkumar Gordhandas Pattani 1986 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. I t i s understood that copying or p u b l i c a t i o n of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. 2 ? The University of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Abstract i i A b s t r a c t This thesis decribes the work on extending the finite element method to cover interaction between viscous flow and a moving body. The problem configuration of interest is that of a two-dimensional incompressible flow over a solid body which is elastically supported or alternatively undergoing a specified harmonic oscillation. The problem addressed in this thesis is that of an arbitrarily shaped body undergoing a simple harmonic motion in an otherwise undisturbed fluid. The finite element modelling is based on a velocity-pressure primitive variable representation of the Navier-Stokes equations using curved isoparametric elements with quadratic interpolation for velocities and bilinear for pressure. The problem configuration is represented by a fixed finite element grid but the body moves past the grid. The nonlinear boundary conditions on the moving body are obtained by expanding the relevant body boundary terms to first order in the body amplitude ratio to approximate the velocities at the finite element grid points. The method of averaging is used to analyse the resulting periodic motion of the fluid. The stability of the periodic solutions is studied by introducing small perturbations and applying Floquet theory. Numerical results are obtained for three different body shapes, namely, (1) a square body oscillating parallel to one of its sides, (2) an oscillating circular body and (3) a symmetric Joukowski profile oscillating parallel to the line of symmetry. The latter case is considered to investigate the flow pattern around an asymmetrical body. In all cases, results are obtained for steady streaming, instantaneous velocity vectors in the fluid domain, added mass, added damping, added force and stability of the flow. A comparision is made between the numerical and published experimental steady streamlines. Very good agreement is obtained for the basic nonlinear phenomenon of steady streaming. i i i C o n t e n t s Abstract i i Figures v Tables Acknowledgement x~" Chapter 1 Introduction and Literature Review 1 1.1 General Remarks and Previous Studies 1 1.1.1 Introduction .1 1.1.2 Fluid Mechanics 3 1.1.3 Fluid-Structure Coupling .4 1.1.4 Analysis in Time Domain 5 1.1.5 Stability of Steady State Solutions 5 1.2 Scope of Present Investigation 7 Chapter 2 Fundamentals of Fluid-Structure Interac-tion 0 2.1 General Description of Governing Equations .9 2.1.1 General Remarks 9 2.1.2 Conservation Equations and Boundary Con-ditions 10 2.1.3 Comparison of Two Associated Flows 12 2.1.4 Non-dimensional Form of Conservation Equa-tions 16 2.2 Finite Element Formulation and Matrix Equations 19 2.2.1 General Remarks 19 2.2.2 Restricted Variational Principle 19 2.2.3 Finite Element Formulation 20 2.3 Moving Body Boundary Conditions 23 2.4 Approximate Method for Equations having Periodic Solutions 27 2.4.1 General Remarks 27 2.4.2 Modified Method of Averaging 27 Chapter 3 Stability of Steady State Solutions 31 3.1 Introduction 31 3.2 General Description of the Governing Equations 32 3.3 Floquet Theory 33 3.4 Numerical Determination of the Transition Matrix 35 3.4.1 General Remarks 35 3.4.2 A Single Pass Integration Scheme 36 Chapter 4 Numerical Investigations 38 Contents iv 4.1 I n t r o d u c t i o n . 38 4.2 C o m p u t a t i o n o f t h e B o d y B o u n d a r y V e l o c i t y 40 4.3 D e t e r m i n a t i o n o f S t reaml ines f o r S teady P a r t o f t h e S o l u t i o n 4 1 4 .3 .1 G o v e r n i n g E q u a t i o n s 4 1 4.3.2 F i n i t e E l e m e n t D i s c r e t i z a t i o n a n d M a t r i x E q u a t i o n s 4 1 4.4 D e t e r m i n a t i o n o f A d d e d M a s s , A d d e d D a m p i n g a n d A d d e d Force 42 4 .4 .1 G e n e r a l R e m a r k s 42 4.4.2 D e t e r m i n a t i o n o f F l u i d Forces 43 4.4.3 C o m p u t a t i o n o f A d d e d M a s s , A d d e d D a m p -i n g a n d A d d e d Force 46 4.5 Resu l ts f o r O s c i l l a t i n g Square B o d y 4 7 4 .5 .1 G e n e r a l R e m a r k s 4 7 4.5.2 F l o w Resu l t s 51 4.5.3 A d d e d Mass a n d A d d e d D a m p i n g 6 1 4.5.4 S t a b i l i t y A n a l y s i s 66 4.6 Resu l ts f o r O s c i l l a t i n g C i r c u l a r B o d y . 7 1 4 .6 .1 G e n e r a l R e m a r k s 71 4.6.2 F l o w Resu l t s 73 4.6.3 A d d e d Mass a n d A d d e d D a m p i n g 87 4.6.4 S t a b i l i t y A n a l y s i s 9 1 4.7 Resu l t s f o r O s c i l l a t i n g S y m m e t r i c J o u k o w s k i P ro f i l e 96 4.7.1 G e n e r a l R e m a r k s 96 4.7.2 F l o w Resu l t s . 98 4.7.3 A d d e d M a s s , A d d e d D a m p i n g a n d A d d e d Force 107 4.7.4 S t a b i l i t y A n a l y s i s 108 C h a p t e r 5 C o n c l u s i o n s 1 1 4 5.1 C o n c l u d i n g R e m a r k s 114 5.2 Suggest ions f o r F u r t h e r D e v e l o p m e n t 115 R e f e r e n c e s 1 1 6 A p p e n d i x 1 M a t r i c e s o f S e c t i o n 2 . 3 1 1 9 A p p e n d i x 2 N e w t o n - R a p h s o n P r o c e d u r e 1 2 1 A 2 . 1 N e w t o n - R a p h s o n I t e r a t i o n P r o c e d u r e 121 A 2 . 2 D e t a i l s o f M a t r i c e s o f Sec t ion 2.4.2 . . 122 A p p e n d i x 3 M a t r i x o f S e c t i o n 3 . 2 1 2 6 V F i g u r e s Figure 2.1 Fluid Domain 10 Figure 2.2 Comparison of Two Associated Flows 13 Figure 2.3 Problem Configuration 17 Figure 2.4 Isoparametric element with quadratic interpola-tion for Velocity and Bilinear Interpolation for Pressure 21 Figure 2.5 Motion of Body 24 Figure 2.6 Body-Fluid Interface Element 26 Figure 4.1 Stress Components 43 Figure 4.2 Cross Section of Body 45 Figure 4.3 Finite Element Grids for Square Body. 48 Figure 4.4 Steady Streamlines for Oscillating Square Body Dn for = 1.99, -r = 30 52 Figure 4.5 Steady Streamlines for Oscillating Square Body for R„ = 26, ^ = 30 53 o Figure 4.6 Steady Streamlines for Oscillating Square Body for R„ = 120, ^ = 30 54 o Figure 4.7 Velocity Vectors for Oscillating Square Body. Rw = 1.99, p = 0.1, = 30 57 o Figure 4.8 Velocity Vectors for Oscillating Square Body. i? w = 26, P = 0.1, ^ = 30 58 Figure 4.9 Velocity Vectors for Oscillating Square Body. Ru = 120, p = 0.1, ^ = 30 59 o Figure 4.10 Added Mass, m 0 for Oscillating Square Body for P = 0 63 Figure 4.11 Added Damping, ca for Oscillating Square Body for/? = 0 65 Figure 4.12 Finite Element Grids for Circular Body. 72 Figure 4.13 Steady Streamlines for Oscillating Circular Body for Ru = 3.834, ^ = 7.7 75 Figures v i Figure 4.14 Steady Streamlines for Oscillating Circular Body for R„ = 2 1 . 3 4 , ^ = 15.5 76 o Figure 4.15 Steady Streamlines for Oscillating Circular Body for = 2 7 8 . 2 , ^ = 30 77 o Figure 4.16 Steady Streamlines for Oscillating Circular Body for j3 = 0.1, ^ = 30 80 Figure 4.17 Steady Streamlines for Oscillating Circular Body for R„ = 2 1 . 3 4 , ^ = 30 81 o Figure 4.18 Velocity Vectors for Oscillating Circular Body. = 3.834, P = 0.1, ^ = 30 84 o Figure 4.19 Velocity Vectors for Oscillating Circular Body. i ? w = 278.2, p = 0.043, ^ = 30 85 Figure 4.20 Velocity Vectors for Oscillating Circular Body. Rw = 324, p = 0.1, ^ = 30 86 b Figure 4.21 Added Mass, ma for Oscillating Circular Body for P = 0 88 Figure 4.22 Added Damping, ca for Oscillating Circular Body for P = 0 89 Figure 4.23 Added Mass, ma for Oscillating Circular Body 92 Figure 4.24 Added Damping, ca for Oscillating Circular Body 93 Figure 4.25 Finite Element Grids for Joukowski Profile 96 Figure 4.26 Steady Streamlines for Oscillating Joukowski Profile for R^ = 55.1, ^ = 18.3 99 o Figure 4.27 Steady Streamlines for Oscillating Joukowski Profile for = 206, Q = 18.3 100 6 Figure 4.28 Velocity Vectors for Oscillating Joukowski Profile. R„ = 55.1, P = 0.01, ^ = 18.3 104 o Figure 4.29 Velocity Vectors for Oscillating Joukowski Profile. R„ = 55.1, P = 0.1, ^ = 18.3 105 b Figure 4.30 Velocity Vectors for Oscillating Joukowski Profile. R„ = 206, P = 0.1, ^ = 18.3 106 b Figure 4.31 Added Mass, ma for Oscillating Joukowski Profile for p = 0 109 Figures Figure 4.32 Added Damping, ca for Oscillating Joukowski Profile for /3 = 0 .110 v i i i T a b l e s Table 4 .1 Maximum Steady Stream Function for Oscillating Square Body 51 Table 4.2 Maximum Steady Component of Velocity for Os-cillating Square Body 56 Table 4.3 Added Mass, ma and Added Damping, ca for Os-cillating Square Body for Ru large and /? = 0. 62 Table 4 .4 Added Mass, ma for Oscillating Square Body for P = 0 62 Table 4.5 Added Damping, ca for Oscillating Square Body for p = 0 64 Table 4.6 Added Mass, m0 and Added Damping, ca for Os-cillating Square Body 66 Table 4 .7 Symmetric Mode Stability Analysis for Oscillat-ing Square Body. = 1.99, = 0.1 6 7 Table 4.8 Symmetric Mode Stability Analysis for Oscillat-ing Square Body. Ru = 1.99, /? = 0.1 68 Table 4 .9 Symmetric Mode Stability Analysis for Oscillat-ing Square Body. R„ = 120, p = 0 .1 69 Table 4 .10 Antisymmetric Mode Stability Analysis for Oscil-lating Square Body. R„ = 120, p = 0.1 70 Table 4 .11 Maximum Steady Stream Function for Oscillating Circular Body. 74 Table 4.12 Maximum Steady Component of Velocity for Os-cillating Circular Body. 82 Table 4 .13 Added Mass, ma Added Damping, ca for Oscil-lating Circular Body for P = 0 87 Table 4.14 Added Mass, ma and Added Damping, ca for Os-cillating Circular Body for j R w large and P = 0 90 Table 4.15 Added Mass, ma for Oscillating Circular Body 90 Table 4 .16 Added Damping, ca for Oscillating Circular Body 9 1 Table 4 .17 Symmetric Mode Stability Analysis for Oscillat-ing Circular Body. Ru = 324 , P = 0 . 1 . 94 Table 4.18 Antisymmetric Mode Stability Analysis for Oscil-lating Circular Body. R^ = 324 , p = 0.1 95 Tables Table 4.19 Table 4.20 Table 4.21 Table 4.22 Table 4.23 Table 4.24 Table 4.25 ix Maximum and Minimum Steady Stream Function for Oscillating Joukowski Profile. 98 Maximum Steady Component of Velocity for Os-cillating Joukowski Profile : 103 Added Mass, ma and Added Damping, ca for Os-cillating Joukowski Profile for i2w large and P = 0 107 Added Mass, ma and Added Damping, ca for Os-cillating Joukowski Profile for P = 0 108 Added Mass, m0, Added Damping, c0 and Added Force, fa for Oscillating Joukowski Profile . . . . I l l Symmetric Mode Stability Analysis for Oscillat-ing Joukowski Profile. = 206, P = 0.0251 I l l Antisymmetric Mode Stability Analysis for Oscil-lating Square Body. Rw = 206, p = 0.0251 112 X A c k n o w l e d g e m e n t The author wishes to express his immense gratitude to his research supervisor Professor M. D. Olson for his advice and encouragement throughout the course of research work and in the preparation of this thesis. The author also wishes to thank Professors M. de St. Q. Isaacson and I. S. Gartshore for numerous helpful discussions. Financial support in the form of a research assistantship from the Natural Sciences and Engineering Research Council of Canada and a university graduate fellowship from The University of British Columbia is gratefully acknowledged. The author finally wishes to thank his wife, Kalyani, for her constant en-couragement and patience throughout the preparation of this thesis. x i To My Mother, Hansa My Father, Gordhandas My Wife, Kalyani C H A P T E R 1 I N T R O D U C T I O N A N D L I T E R A T U R E R E V I E W 1 . 1 G E N E R A L R E M A R K S A N D P R E V I O U S S T U D I E S 1.1.1 I n t r o d u c t i o n . In recent times, significant attention has been devoted to the development of fluid-structure interaction analysis procedures. However, the complexity of this class of problems is such that very few analytical solutions exist and therefore the necessity for numerical methods of analysis has arisen. Of these methods, the finite element method allows the investigator to tackle the problem most elegantly in all its aspects, that is, complicated geometry and nonlinearity in fluids as well as structural representation. Basically, there are two main classes of problems encountered in fluid-str-ucture interaction. They are either (1) structures surrounded by fluids, or (2) fluids con-tained within structures. Examples of (1) are wave forces on offshore structures, flow induced vibrations in aerospace structures and pressurized water reactor cores. Examples of (2) are the seismic response of ground supported, cylindrical liquid storage tanks, vibra-tions of biomechanical structures such as artificial heart and aerospace applications such as rocket fuel tanks. The finite element methods have been applied to fluid-structure interaction problems by several investigators. A large volume of research in this area is based on 1 INTRODUCTION AND LITERATURE REVIEW / 1.1 2 inviscid fluid both compressible or incompressible with emphasis on aerodynamic and offshore structure applications [1-5]. The work of applying finite element methods to viscous fluid-structure interaction problems have been carried out by Olson et al. [6,7] and Liu [8]. Liu [8] treats the problem of fluid-structure interaction in a very general way and therefore it does not lend itself to economical utilization for important subclasses of problems. One such subclass, found in offshore technology, is wave flow past a section of a structural member. There are strong similarities between the wave flow past a section of a struc-tural member and a two-dimensional harmonic flow past a similar section. The latter reference flow, which is conveniently represented analytically and obtained under labora-tory conditions, is generally used to generate a variety of pertinent results. For most applications of practical interest the forces acting on a structural section in two-dimensional flow are given by the so called Morison's equation. This equa-tion is empirical in nature but most extensively used in view of the fact that even in a steady, two-dimensional, separated flow past a smooth cylinder one does not have either a theoretical or a numerical solution which explains all the characteristics of the flow as functions of Reynolds number. The solution of harmonic flow past a section of a structural member is closely related to the solution of a harmonically oscillating section of a structural member in an otherwise still fluid [9]. This in turn is related to the solution of a rigid structural member which is elastically supported and subject to an external harmonic forcing function in an otherwise still fluid. Solution of such a system for an incompressible real (viscous) fluid would involve the solution of the complete Navier-Stokes equations with appropriate no-slip boundary conditions. The presence of nonlinear convective acceleration terms in the Navier-Stokes equations makes the solution difficult. Stokes [10] obtained the solution for such a system for a circular section by neglecting the nonlinear convective acceleration INTRODUCTION AND LITERATURE REVIEW / 1.1 3 terms. This can be done for the flow for which viscous forces are dominant, that is, when the velocity is very small, or, speaking more generally, when Reynolds number is very small. For very large Reynolds numbers the two-dimensional flow of a real (viscous) fluid past a body began to be understood only in the present century following Prandtl's discovery of the boundary layer in 1904. Analytic solutions have been obtained for bound-ary layer equations for the flow around an oscillating circular cylinder [11]. These solutions give the same surface skin friction as obtained by Stokes using a linearized form of Navier-Stokes equations; but the pressure distribution is assumed to be, and remains, that of the inviscid theory. Hence when the primary interest is in the force and therefore the pressure, on the body, the boundary layer theory cannot be used. The ability to solve this type of fluid-structure interaction problem by finite element methods would depend on the ability to solve the corresponding uncoupled fluid problem by finite elements, and also the ability to interface fluid and structural subdo-mains. In engineering applications it is of interest to determine critical conditions for the occurence of instability for infinitesimal disturbances. The ability to determine stability of the flow in this type of fluid-structure interaction problem would involve the ability to determine the stability of time dependent phenomenon. It is the author's goal in this thesis to obtain the solution for the flow around an elastically supported rigid body and investigate the stability of this solution. 1.1.2 F l u i d M e c h a n i c s . In the area of fluid mechanics, much work has been done in recent years to develop numerical solutions for two-dimensional flow around bodies using finite element and finite difference techniques. Results of acceptable accuracy have been obtained for flows and forces on fixed bodies for a finite range of Reynolds number. INTRODUCTION AND LITERATURE REVIEW / 1.1 4 T h e r e are essent ia l ly t h r e e f i n i t e e lement approaches t o dea l w i t h i n c o m -press ib le v iscous flow, based o n (1 ) t h e p r i m i t i v e va r iab les u , v , p ; (2 ) t h e s t r e a m f u n c t i o n a n d v o r t i c i t y (3 ) t h e s t r e a m f u n c t i o n V a lone . Essent ia l l y , i n each case t h e d e p e n d e n t va r iab les are a p p r o x i m a t e d p iecewise across t h e d o m a i n b y f i n i t e e lement i n t e r p o l a t i o n s t o o b t a i n t h r o u g h v a r i a t i o n a l o r w e i g h t e d res idua l f o r m u l a t i o n a n o n l i n e a r s y s t e m o f m a t r i x e q u a t i o n s . O l s o n a n d T u a n n [12] have presented an excel lent c o m p a r i s i o n be tween v a r i o u s m e t h o d s . T u a n n a n d O l s o n [13] have also o b -t a i n e d finite e lement s o l u t i o n s f o r s teady flow a r o u n d c i r c u l a r c y l i n d e r s i n t h e R e y n o l d s n u m b e r range f r o m 1 t o 100. T h e same t h r e e approaches m e n t i o n e d i n c o n t e x t w i t h f i n i t e e lement m e t h o d s have been t r a d i t i o n a l l y e m p l o y e d i n finite d i f ference m e t h o d s . Fo r a d e t a i l e d ana lys is o f t h e t h r e e f o r m u l a t i o n s a n d t h e finite d i f fe rence techn iques , t h e reader is r e f e r r e d t o Roache 's b o o k Computational Fluid Dynamics [14] . M o r e r e c e n t l y exce l lent resu l t s have been o b t a i n e d b y D a v i s a n d M o o r e [15] f o r u n s t e a d y flow a r o u n d squares f o r R e y n o l d s n u m b e r be tween 100 a n d 2800 u s i n g finite d i f ference m e t h o d s . O f t h e n u m e r i c a l m e t h o d s i n t r o d u c e d i n t o mechan ics w i t h i n t h e last severa l decades t h e finite d i f fe rence m e t h o d has become p r o m i n e n t i n t h i s field. T h e finite e lement m e t h o d is a r e l a t i v e l y recent a r r i v a l b u t has a l ready become v e r y p o p u l a r a n d seems l i ke l y t o d isp lace t h e finite d i f ference m e t h o d i n m a n y p r o b l e m areas because o f t w o i nhe ren t advan tages . F i r s t l y , a c o m p l e x b o u n d a r y g e o m e t r y is m o r e r e a d i l y a c c o m m o d a t e d i n a finite e lement f o r m u l a t i o n , a n d secondly , a finite e lement p r o g r a m is m o r e r e a d i l y a d a p t a b l e t o subsequent p r o b l e m s o f d i f fe ren t c o n f i g u r a t i o n a n d spec i f i ca t i on . 1 .1 .3 F l u i d - S t r u c t u r e C o u p l i n g . F i n i t e e lement m e t h o d s c a n be used t o analyse b o t h fluid a n d s t r u c t u r a l p a r t s o f a c o u p l e d fluid-structure p r o b l e m . T h e equa t ions o f m o t i o n i n t h e fluid a n d so l id INTRODUCTION AND LITERATURE REVIEW / 1.1 5 are all expressed in terms of nodal variables and are of the same form. The equations of motion governing the coupled system may thus be solved simultaneously and a strong cou-pling between fluid and solid responses is achieved. However, suitable conditions must be prescribed along fluid-solid interfaces to match the velocities, tangential and normal to the solid without gross distortion to the finite element grid. To achieve this Hughes et al. [16] have proposed a mixed Langragian-Eulerian description in which each degree-of-freedom may be assigned to move at a fraction of the fluid particle velocity. For fluid-structure problems of the type being dealt here it has been indicated by Olson and Pattani [7] that it is possible to keep the finite element grid fixed, but allow the body to move past the grid and expand the relevant boundary terms by the Taylor series to approximate the velocities at the finite element grid points. 1.1.4 Ana l y s i s i n Time Domain. The techniques available for determining the steady-state behaviour of forced oscillations of nonlinear systems can be divided into two different groups. The first group includes the methods of averaging and multiple scales. With this group one determines the equations describing the amplitudes and phases. The second group includes the Lindst£dt-Poincare1 technique, the method of harmonic balance and the Galerkin procedure. With this group one determines directly the steady-state solutions. For the details of these techniques the reader is referred to Nayfeh and Mook's book Nonlinear Oscillations [17]. The quadratic nonlinear terms in the Navier-Stokes equations are responsible for the phenomenon of steady streaming in which the solution for velocities have a steady part as well as an oscillating part. The method of averaging needs to be modified, as will be outlined later on in this thesis, to obtain the steady streaming part. 1.1.5 S t a b i l i t y of Steady State Solutions. The transition from steady state to transient laminar flow with increasing Reynolds number is considered to be an instability and represents bifurcation of the solu-tion. In fact the complete transition to turbulent flow is thought to be a series of succesive INTRODUCTION AND LITERATURE REVIEW / 1.1 6 b i f u r c a t i o n s . T h e s t a b i l i t y o f a l a m i n a r f low is o f p r i m e i m p o r t a n c e i n m a n y eng ineer ing a p p l i c a t i o n s . T o u n d e r s t a n d t h e phys ics o f t h e p r o b l e m i t is necessary t o k n o w t h e c r i t i c a l R e y n o l d s n u m b e r a t w h i c h t h e l a m i n a r f low loses i t s s t a b i l i t y , a n d t h e n a t u r e o f i n s t a b i l i t y . T h e n u m e r i c a l s t e a d y s t a t e s o l u t i o n m e t h o d s f o r v iscous incompress ib le flow p r o b l e m s gove rned b y t h e Nav ie r -S tokes e q u a t i o n s is f a i r l y w e l l d e v e l o p e d . However , t h e assoc ia ted s t a b i l i t y o f these flows has seen m u c h less e x p l o r a t i o n . L i n [18] has w r i t t e n a m o n o g r a p h d e a l i n g w i t h some genera l concep ts o f s t a b i l i t y . M o r e r e c e n t l y Joseph [19] has dea l t w i t h t h e s t a b i l i t y p r o b l e m i n t h e l i g h t o f b i f u r c a t i o n t h e o r y . O n e o f t h e s i m p l e r s t a b i l i t y p r o b l e m s , governed b y t h e O r r - S o m m e r f e l d equa-t i o n w h i c h co r responds t o a l i n e a r p e r t u r b a t i o n o f t h e Nav ie r -S tokes equa t ions a n d reduc-t i o n t o a o n e - d i m e n s i o n a l e igenvalue p r o b l e m f o r Po iseu i l le flow, has been ex tens ive ly i n v e s t i g a t e d . T h e p r o b l e m has been t r e a t e d u s i n g an a s y m p t o t i c m e t h o d b y L i n [18] a n d u s i n g finite d i f ferences b y T h o m a s [20] . G r o s c h a n d Sa lwen [21] have used a n e igenfunc-t i o n e x p a n s i o n a n d O r z a g [22] has used C h e b y s h e v p o l y n o m i a l s . M o r e recen t l y , S a r a p h et a l . [23] have used t h e finite e lement a p p r o a c h f o r s t a b i l i t y s tud ies u s i n g l i near i zed g o v e r n i n g equa t i ons i n s t e a d o f r e s o r t i n g t o t h e u s u a l O r r - S o m m e r f e l d e q u a t i o n . T h e O r r - S o m e r f e l d e q u a t i o n has also been a p p l i e d t o t h e b o u n d a r y layer s t a b i l i t y p r o b l e m . G r o s c h a n d Sa lwen [24] e x t e n d e d t h e i r p rev ious w o r k a n d m o r e r e c e n t l y V a n S t i j n a n d V a n de V o o r e n [25] have dea l t w i t h t h e p r o b l e m . Fasel [26] has also s t u d i e d t h e b o u n d a r y layer p r o b l e m b u t b y a n e n t i r e l y d i f fe ren t a p p r o a c h . H e s t a r t s w i t h a finite d i f fe rence r e p r e s e n t a t i o n o f t h e u n s t e a d y Nav ie r -S tokes equa t i ons a n d i n t r o d u c e s a p e r i o d i c d i s t u r b a n c e a n d t h e n examines i t s b e h a v i o u r b y i n t e g r a t i n g i n t i m e . Recen t l y , O l s o n a n d B a l d w i n [27] have used t h e finite e lement a p p r o a c h t o dea l w i t h a genera l l i nea r i zed p e r t u r b a t i o n o f t w o - d i m e n s i o n a l Nav ie r -S tokes equa t ions f o r s t e a d y s o l u t i o n s . T h e genera l l i near i zed p e r t u r b a t i o n o f a finite e lement rep resen ta t i on of t h e t w o - d i m e n s i o n a l u n s t e a d y Nav ie r -S tokes equa t i ons has n o t been p r e v i o u s l y dea l t w i t h . INTRODUCTION AND LITERATURE REVIEW / 1.2 7 1 . 2 S C O P E O F P R E S E N T I N V E S T I G A T I O N This thesis describes the development work on extending the finite element method to cover interaction between viscous flow and a moving body and to investigate the stability of the resulting flow. The problem configuration of interest is that of two-dimensional incompressible flow over a solid body which is elastically supported or alter-natively undergoing a specified harmonic oscillation. In Chapter 2, the fundamentals of fluid-structure interaction and the finite element formulation are discussed. The chapter begins by stating the governing differ-ential equations for the fluid and choosing the non-dimensional parameters to effectively characterize the present flow situation. The flow situation in which the fluid remote from the stationary body has a uniform acceleration will be compared to the related situation in which the identical body has an equal acceleration in the opposite direction, while the fluid is otherwise stationary. It will be shown that it is not necessary to solve both these cases independently as the results from one can be used to obtain pertinent results for the other. The restricted variational principle governing the flow situation will be discussed and finite element matrix equations will be obtained for a suitable choice of an element. The discretized form of the boundary conditions are then obtained. The method of slowly varying amplitudes is applied to obtain the approximate periodic solutions of the finite ele-ment ordinary differential equations. A Newton-Raphson iteration scheme will be outlined for the solution of the resulting algebraic equations. In Chapter 3, the stability of the steady state solutions obtained in Chapter 2, is discussed. The chapter begins by obtaining the governing differential equations for small superimposed disturbances to the steady state solutions. Floquet theory is then applied to characterise these small disturbances. An efficient numerical integration scheme is presented to obtain the transition matrix characterizing the solution of the governing differential equations. INTRODUCTION AND LITERATURE REVIEW / 1.2 8 The numerical investigations are described in Chapter 4. A finite element scheme is outlined to obtain streamlines for the steady part of the velocity field solutions. A numerical integration scheme to compute added mass and added damping is presented next. Numerical results are then obtained for three different body cylinders, respectively, (1) a square cylinder oscillating parallel to one of its sides, (2) an oscillating circular cylinder and (3) a symmetric Joukowski profile cylinder oscillating parallel to the line of symmetry. The latter case is considered in order to investigate the flow pattern around an asymmetrical body. In all cases results are obtained for the steady streaming, the added mass and added damping, and for the stability of the periodic flow. Conclusions and suggestions for further developments are presented in Chap-ter 5. C H A P T E R 2 F U N D A M E N T A L S O F F L U I D - S T R U C T U R E I N T E R A C T I O N 2 . 1 G E N E R A L D E S C R I P T I O N O F G O V E R N I N G E Q U A T I O N S 2 .1 .1 G e n e r a l R e m a r k s In this section the governing equations for the coupled viscous fluid-solid body system are presented. For the fluid domain the velocity-pressure primitive variable v form of the Navier-Stokes equations for two-dimensional incompressible viscous flow is used. The equation of motion of an elastically supported rigid body with a single transla-tion degree of freedom is also presented. This elastically supported rigid body is subject to an externally specified forcing function which alternatively may be considered to be undergoing a specified motion. The most common applications would involve predicting the forces on and hence the response of a structural section represented by a spring mass system in a two-dimensional flow. The structure itself might be initially stationary, and subjected to a specified time dependent flow, or it can be subjected to an external force in an otherwise still fluid. Isaacson[9] has indicated that these problems are interrelated and the solution of one leads to the solution of the other without much additional computational effort. There are strong similarities between the wave flow past a section of a struc-tural member and a two-dimensional harmonic flow past a similar section (see Sarpkaya 9 FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.1 10 and Isaacson[33]). The solution of the latter reference flow is closely related to the solu-tion of a harmonically oscillating section of a structural member in an otherwise still fluid. Therefore, the problem configuration considered here is that of a elastically supported rigid body of arbitrary shape undergoing specified harmonic oscillations in an otherwise still fluid. The equations of motion are non-dimensionalised to effectively characterize this flow situation. 2 .1 .2 C o n s e r v a t i o n E q u a t i o n s a n d B o u n d a r y C o n d i t i o n s flow are the Navier-Stokes equations, the equations of continuity, and the boundary and initial conditions. Solution of the flow problem will be sought within a plane domain fi bounded by a contour T which, in general, is composed of two distinct parts denoted by r « and T „ respectively (see Figure 2.1). The equations of motion in Cartesian coordinates (z,y) are given by (see Schlichting[ll]) F I G U R E 2 .1 Fluid Domain. The basic equations governing two-dimensional, viscous incompressible fluid Du _ 1 (dax Dt ~ p\dx Dv _ 1 (drn Dt ~ p\dx (2.1) FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.1 11 where u, v are the fluid velocities in the z,y directions, respectively; ax,<Ty,Txy are the normal stresses in z, y direction and the shear stress respectively; p is the fluid density; ^ is the material derivative given by Dt ~ dt + U d x + v 6 y The constitutive relation for stresses in terms of velocities is given by <rx = -p + 2 / i — ox Txy - / X \dx + dy) where p is the fluid pressure and ft is the fluid absolute viscosity. Introducing equations (2.2) into the fundamental equations of motion (2.1), we obtain Du _ _ldp / dhi cPu d2v \ D t ' p d x + U \ dx2 +dy2 +dxdy) , , n » . A » > + j * + i £ + * l { ( i , ! 0 e n ' < > o ( 2 - s ) Dt pdy \dx2 dy2 dxdy) where u = ^ is the fluid kinematic viscosity. These very well known differential equations are usually referred to as the Navier-Stokes equations. The equation of continuity for two-dimensional incompressible flow assumes the following form: dv „ ,„ ,x fa**-0 <2-4> Boundary conditions are given on the velocity components over the part boundary r„, referred to as the kinematic boundary, and traction over the remaining part Tf, referred to as the mechanical boundary. Therefore, the boundary conditions are u = u, v = C (i, y ) € r w > < > 0 r (d d M r a 1 - (*.sOer„oo FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.1 12 where (u, v) are the specified velocities on r „ , (X.Y) are the specified tractions on Tf, and(n 1 , n 2 ) are the direction cosines of the outward pointing normal to the boundary. The equation of motion for an elastically supported rigid body of arbitrary shape is given by: ms + ks = f{t) (2.6) where m is the mass of the body, k the elastic spring constant, s the dispacement motion and the time dependent loading consists of two components f(t) = Ff(t)+Fe(t) (2.7) where Ff(t) is the fluid force and Fe(t) is the external force on the body. 2 . 1 . 3 C o m p a r i s o n o f T w o A s s o c i a t e d F l o w s We consider an associated pair of two-dimensional flows as defined in Fig-ure 2.2 (see also Isaacson[9]). In Case I the fluid remote from a stationary body has a uniform time dependent velocity U{t) and in Case II the identical body has an equal velocity U(t) in the opposite direction, while the fluid remote from the moving body is stationary. The cartesian axes are so defined that these velocities lie in the x direction. The flow quantities in the co-ordinate system fixed to the body in Case I are denoted by subscript a, that fixed relative to the body in Case II by subscript b and fixed relative to the undisturbed fluid in Case I f by subscript c. The equations of motion and continuity, given by (2.3) and (2.4) for any fluid particle in Case I are: Dua 1 dpa "a __*_£ Dt ~ pdxa Dt V dx\ dyl dxadya) = _ I ^ + „ ( f £ + 2 f ^ + / ^ ) ( 2 .8 ) pdya \dxl dyl dxadyaJ dua dva dxa dya where the material derivative j j j is given by D_ _ _3_ a Dt~ d t + U a d x a + V a d y a FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / Zl Case I Fluid remote from body with uniform time dependent velocity U(t). Inertial co-ordinate system (xai ya) is fixed with respect to the body. 7b Case II Stationary fluid remote from the body moving with velocity —U(t). Noniner-tial co-ordinate system (ij,, yj,) fixed wi-th respect to the body is moving with velocity — U(t) relative to inertial co-ordinate system (xc,yc) fixed with re-spect to the undisturbed fluid. F I G U R E 2.2 Comparison of Two Associated Flows. FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.1 14 and the boundary conditions are: « a = 0 va = 0 on the body boundary ua = U(t) va=0 in the far field For a fluid particle in Case II, the equations of motion and the continuity equation in the inertial co-ordinate system [xc, yc) are: Dt pdx, \ dx\ * dyl + dxcdyj pdye +V\dxl + 'dy* +dxedyj ( > ^ £ + ^ = 0 dxe dye Dt where the material derivative j j j is given by Z> _ d_ _d_ _d_ Dt ~ dt + Ucdxc +Vcdyc' Now the velocities of any point in terms of co-ordinate systems (xb, yb) and (xCJ yc) are related by: ue = ub-U »c = Vb And the pressure at any point is independent of the frame of reference: (2.10a) Pc = Pb (2.106) The axes are related by: xc = xb— I U(t) dr /o (2.10c) ye = yb The absolute acceleration of a fluid particle, in the noninertial co-ordinate system (xbt yb) moving with velocity — U(t) relative to the inertial co-ordinate system (xc,yc), is given by (see Batchelor[28]) Due = Dub dU{t) d t <>•»> Dt Dt FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.1 15 Hence, upon substituting (2.10) and (2.11) into (2.7), we obtain dt • Dt pdxb V dx\ dy\ ^ dxbdyb) Dv, ~Dt I P9yb \dx2b dyl + dxbdyb) dub dvb _ Q or, on rearrangement Dt p \dxb P d t ) ^ \ dxy dy\ + dxbdyj Dvb_ ldPb (d\ d\ d*ub \ ( -Dt--pdyb + l'{dxl +2dyl+dxbdyb) ( W 2 ) dub dvb _ Q where the material derivative is given by D__d_ _d_ d Dt~ dt+Vbdxb + Vbdya and the appropriate boundary conditions are: ub = 0 ffe = 0 on the body boundary ub = U(t) vb = 0 in the far field The differential equations and the boundary conditions are the same in both Case I (equations (2.8)) and Case II (equations (2.12)). Therefore, assuming a unique solution for all t, we have va = n dPa^dPb (2-13) dya dyb dpa _ d£b _ dU_ dxa dxb ^ dt The consequence of equations (2.13) is that the velocities, and hence the flow pattern are identical in both Case J and Case II. The only difference between the two cases as related to pressure is the additional uniform pressure gradient — p^ in Case / and acting in the direction of acceleration. FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.1 16 I t is poss ib le t o cons ider t h e wave flow pas t a s t r u c t u r a l m e m b e r b y cons id-e r i n g a t w o - d i m e n s i o n a l h a r m o n i c flow pas t a s i m i l a r sec t ion ( S a r p k a y a a n d Isaacson [33] ) . T h i s s i t u a t i o n c a n be o b t a i n e d b y l e t t i n g U(t) be a h a r m o n i c f u n c t i o n i n Case / . I i a p r o b l e m is ana lysed i n t h e i n e r t i a l c o o r d i n a t e s y s t e m ( x c > yc) as p resc r ibed i n Case II f o r a h a r m o n i c v e l o c i t y — U(t)t t h e n i t is poss ib le t o o b t a i n t h e flow q u a n t i t i e s i n t h e n o n -i n e r t i a l c o o r d i n a t e s y s t e m (x^J t o ) b y u s i n g equa t i ons (2 .10 ) . F u r t h e r m o r e , t o o b t a i n t h e flow q u a n t i t i e s i n Case / f o r a h a r m o n i c U(t) t h e t r a n f o r m a t i o n s g i ven b y equa t i ons (2.13) need t o b e p e r f o r m e d . T h e r e f o r e , i t is poss ib le t o ana lyse t h e case o f a sec t ion o f s t ruc -t u r a l m e m b e r u n d e r g o i n g a h a r m o n i c m o t i o n i n an o the rw ise s t i l l fluid a n d o b t a i n t h e flow q u a n t i t i e s f o r t h e case o f a t w o - d i m e n s i o n a l h a r m o n i c flow pas t a s i m i l a r sec t ion b y t h e a p p r o p r i a t e t r a n s f o r m a t i o n s . 2 . 1 . 4 N o n - d i m e n s i o n a l F o r m o f C o n s e r v a t i o n E q u a t i o n s T h e p r o b l e m u n d e r c o n s i d e r a t i o n is t h a t o f a b o d y o f a r b i t r a r y shape w i t h a cha rac te r i s t i c l e n g t h b i n t h e d i r e c t i o n o f m o t i o n , u n d e r g o i n g a s i m p l e h a r m o n i c d isp lace-m e n t m o t i o n i n a n o t h e r w i s e s t i l l fluid ( F i g u r e 2 .3 ) . T h e s i m p l e h a r m o n i c d i sp lacement m o t i o n t h e b o d y is u n d e r g o i n g is g i ven b y s = Sq s in ut w h e r e s 0 is t h e d i sp lacemen t a m p l i t u d e a n d u t h e impressed f requency . T h e n t h e char -ac te r i s t i c v e l o c i t y i n t h e flow is t h a t o f t h e b o d y , n a m e l y Uq = us0. I n t r o d u c i n g t h e n o n - d i m e n s i o n a l va r iab les u - v z P U = V = P— 7T i n t o Nav ie r -S tokes e q u a t i o n s (2.3) a n d c o n t i n u i t y e q u a t i o n ( 2 . 4 ) , we o b t a i n du ( du du\ _ _ 1 _ ( dhx_ cPu d 2v _ dp\ dt + ^ \Udx + Vdy) ~ \ dx* + dy 2 + dxdy dx) dv a ( dv dv\ 1 (d 2v d 2v d 2u dp\ fo . A , du dv dx dy FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.1 FIGURE 2.3 Problem Configuration. FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.1 18 where all quantities are their respective non-dimensional values and the overtilda has been dropped for convenience. It should be noted here that the pressure has been non-dimension alised with respect to the characteristic shear stress as opposed to the character-istic dynamic pressure since it is anticipated that the flow regime would be a slow, shear dominated one. There are two natural Reynolds numbers in the problem, namely the usual Reynolds number and the frequency Reynolds number given respectively by Re=*± and i C = — • V V Note that P K b is the body amplitude ratio. Introducing non-dimensional variables u = u V = V P = X y -X = b y = b u = V = 0 u 0 x = X fiun/b Y = fMU0/b u Y fiu0/b into boundary conditions (2.5), we obtain u = u, v = v (x,y) € r „ , i > 0 {x,y)eT„t>0 P + 2 ^ ] n i + [ ( ^ + f i ) ] r i 2 = = X [(f^ + l ; ) ] ^ + h + 2 ^ ] ^ = F (2.15) where all quanties are their respective non-dimensional values and the overtilda has again been dropped for convenience. Introducing non-dimensional variables u = u V = V P~ P Uq Uq fiuo/b j _ T*y Hu0/b °y - {MU0/b x y fiu0/b X y b y = b into the constitutive equations (2.2), we obtain du °* = -*+2irx dv <Ty = -p + 2n7, (2-16) dy Tx» ~ \dx + dy) FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.2 19 where again all quantities are their respective non-dimensional values. Introducing non-dimensional variables 3 - m s = - t = ut m = b pAbl 7 k "rT _ -^e ~ Ft k = tt = 1 f = —L_ where Ab is the cross-sectional area of the body cylinder and / the length of the cylinder, into equation (2.6), we obtain ms + ks = Jb(Fe + Ff) (2.17) where again all quantities are their respective non-dimensional values, and r = J^o_ = l_ ( V \ 6 uHAb K \AbJ' 2.2 F I N I T E E L E M E N T F O R M U L A T I O N A N D M A T R I X E Q U A T I O N S 2.2.1 General Remarks. In this section the Navier-Stokes equations presented in section 2.1 are dis-cretized by the finite element methodology. The starting point in the finite element methodology is usually the functional equivalent form of the differential equations. Olson and Tuann[12] have indicated that it is convenient to recast into a restricted variational principle, which is exactly equivalent to the Galerkin method. This process is outlined here. The finite element discretization of the restricted variational form of the Navier-Stokes equations is carried out using curved isoparametric elements with quadratic interpolation for velocities and bilinear interpolation for pressure. 2.2.2 Restricted V a r i a t i o n a l P r i n c i p l e . FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.2 20 It is now well accepted (Finlayson[29]) that no variational principle in the classical sense can be constructed for the full Navier-Stokes problem represented by equa-tions (2.14) and (2.15). Nevertheless, it is possible to construct a restricted variational principle equivalent to the field equations (2.14) and (2.15) in the form (2.18) - £ ( s + S ) } " - J U < * - + * > « By restricted we mean there are two types of velocity, (u, v) and (u°, y°), where (u°, v°) is associated with inertial terms and (u, v) with the remaining terms. When the first variation is taken, the (u°, v°) is held constant, and at the end of the process the two velocities are equated, thereby restoring the original field equations (2.14) and (2.15). This is identical to the Galerkin method but allows one to think in the usual terms of seeking a stationary point to some functional. 2 . 2 . 3 F i n i t e E l e m e n t F o r m u l a t i o n . Olson and Tuann[12] have shown that the finite element interpolation for pressure should be no higher than that for the dilation | | + |^ in order to avoid spurious singularities for the (u, v, p) integrated formulation. Therefore, the finite element interpo-lation for pressure p should be at least one degree less than that for the velocity components (u, v). Also, the form of it in equation (2.18) reveals that only C° continuity is required of the finite element interpolations for each of the (u, v,p). To satisfy both these conditions the finite element representation of the Navier-Stokes equations (2.14) is achieved here us-ing curved isoparametric elements with quadratic interpolation for velocities and bilinear interpolation for the pressure. FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.2 21 FIGURE 2 .4 Isoparametric Element with Quadratic Interpolation for Velocity and Bilinear Interpolation for Pressure. (a) Element in (s.t) space, (b) Mapping of the element into (x, y) space. In isoparametric elements the shape function used are the same both for interpolating the variables u,v,p and for transformation from 8,t natural co-ordinates to x,y co-ordinates (see Figure 2.4). The velocities and pressure are represented by u = JV,u,(0 t = l,2,...,8 u = NiVi{t) s = l,2,...,8 (2.19) p = MiPi(t) • = !,...,4 FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.2 where iV,-,M,- are the shape functions given by 22 *1 = - 1 ( 1 - . ) ( ! - + • + « ) isr2 = -1(1+ s ) ( l -0 ( l - s + 0 N3 = - l ( l + 5 ) ( l _ 0 ( l - s - 0 JV4 = - l ( l - S ) ( l + 0(l + s - 0 Nb = \(l-s7)(l-t) tf6 = \{l + s){l-t2) N7 = l ( l - s 2 ) ( l + 0 JV8 = 1 ( 1 - , ) ( ! - « » ) . M1 = 1(1-,)(!-<) M 2 = l ( l + 3 ) ( l -0 M 3 = l ( i + *)(i + 0 M 4 = 1 ( 1 - « ) ( ! + «) and u,-,v,-,p,- are the time dependent nodal variables as shown in Figure 2.4(b) Substituting equations (2.19) into equation (2.18), carrying out the integra-tion over one element and following the restricted variational principle, yields FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.3 23 where the arrays are M? = f j NiNJdA = MJf *>jt k — 1,2,..., 8 dy dx „* ft JA „„ f f 3JV,- JA » = 1 | 2 , . . . , 8 2.3 M O V I N G B O D Y B O U N D A R Y C O N D I T I O N S The body in an otherwise still fluid is describing a motion given by (see Figure 2.5) s(t) = s0 sin ut (2.21) or in terms of velocities u(s(<)) = u(t) = u0cosu;£ (2.22) v(s(t)) = v{t) = 0 where u0 = Sqoj. It is assumed that the finite element grid remains fixed at the mean position of the body but the body moves past the grid. To obtain velocities at the mean position of the body for any t > 0, we expand the relevant boundary terms by the Taylor series: «(«) = «(«)+ «(?••:) +••• )a\ ( 2 ' 2 3 ) •W—(o) + . ( ^ ) +••• FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.3 24 y • FIGURE 2.5 Motion of Body. where the subscript 0 indicates that the derivatives are evaluated at s = 0, that is, at the mean position of the body. Truncating the series at order s and substituting from equations (2.21) and (2.22) into equations (2.23) and rearranging, we obtain (2.24) To non-dimensionalize equations (2.24), we introduce the variables t = into equations (2.24), to obtain (2.25) FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.3 25 where all quantities are their respective non-dimensional values and the overtilda has been dropped for convenience. The zero position in equation (2.25) indicates that the boundary conditions need to be applied to the finite element grid edges representing the boundary of the body in its mean position. Note that the above procedure can be applied to any description of s(t) so long as s(t) is small compared to the body dimension b. It is seen from equations (2.14) that /? is the obvious perturbation param-eter in that it governs the nonlinear convection terms. Therefore, representation of the boundary condition to the order ft in the equations (2.25) is consistent with the order of nonlinearity in the Navier-Stokes equations. To discretize the boundary condition in equations (2.25), we substitute equa-tions (2.19) into them for nodal variables at the edge of the finite element grid interfacing with the mean position of the body, to obtain u(0) = cos t - B ( ^ u , J sin t ON \ (2' 26) v(0) = - / 9 ^ « / , j sin* For one element on the boundary as shown in Figure 2.6, the equations (2.26) for the nodal variables on the edge interfacing with the mean position of the body reduce to u,- = cos t — pC^Uj sin t - /?C i k sin t cos t (2.27) t;,- = — pC^Vj-s'mt where is summed over the velocity degrees of freedom other than those on the edge interfacing with the mean position of the body, k is summed over the velocity degrees of freedom on the edge interfacing with the mean position of the body, and • m, C 8 = I " X T ) (2-28) dN-where the subscript i indicates that is evaluated at the location of u,- or whichever is appropriate. When the finite element equations (2.20) are assembled there are (say) net u and v nodal variables numbering n each and net p nodal variables numbering m. Out FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.3 26 Body Boundary FIGURE 2.6 Body-Fluid Interface Element. of these, the u and v nodal variables each numbering r are located on the boundary and are known from equations (2.27). These variables are segregated by suitable matrix manipulation and transformed to the right hand side. This results in the matrix equations of the form [M + p sin iP]d + [K + p cos tQ + P sin <R]d = F sin t + G cos t + p{H cos 21 + L sin t cos t + J sin 21} where M is the moss matrix given by M = and K is the stiffness matrix given by Mg 0 0 0 MJf 0 0 0 0 Kg -Pl -P*. -py. 0 (2.29) FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.4 27 and d = {uj u2 . . . u„_ r vx v2 ... «n_r P l p2 ... Pm}T is the nodal vector of unknowns. Details of all other matrices are given in Appendix 1. 2.4 A P P R O X I M A T E M E T H O D F O R E Q U A T I O N S H A V I N G P E R I O D I C S O L U T I O N S 2.4.1 General Remarks. The techniques available for determining the steady-state behaviour of the forced oscillation of a nonlinear system can be divided into two basic groups. The first group includes the methods of averaging and multiple scales. With this group, the equa-tions describing the amplitudes and phases are first determined. These equations are then transformed into an autonomous system, where the steady-state solutions correspond to the singular points of this autonomous system. The second group includes the Lindstedt-Poincare" technique, the method of harmonic balance and the GaleYkin procedure. With this group one determines directly the steady-state solutions. For the details of these tech-niques the reader is referred to the excellent book Nonlinear Oscillations by Nayfeh and Mook[17]. The quadratic nonlinear terms in the Navier-Stokes equations are responsible for the phenomenon of steady streaming in which the velocities have a steady part as well as a time dependent part. The method of averaging needs to be modified to obtain the steady streaming part of the solution from equations (2.29). The modified method of averaging is illustrated by obtaining the steady state solution of the first order ordinary differential equations (2.29) governing the fluid-structure interaction problem under consideration. 2.4.2 Modified Method of Averaging. Let us assume that the solution of equations (2.29) for small /? takes the form d = A+B(t)cos* + C(r)sin* (2.30) FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.4 28 where B(£)>C(£) are assumed to be slowly varying functions of non-dimensional time t. The first term A represents the steady streaming part of the solution which naturally arises for the system with quadratic non-linearities as is encountered here. From equation (2.30), we obtain d = -B{t) sin t + C(t) cos t + B(i) cos t + C(t) sin t (2.31) To obtain an autonomous system of equations governing the amplitudes A , B(tf)> and C(£)> equations (2.30) and (2.31) are first substituted into equations (2.29). Then to obtain the three sets of equations for the three sets of unknowns A , B(<) and C(i), the resulting equations are (1) averaged over the period 2ir, (2) multiplied by sin t and averaged over the period 2k and (3) multiplied by cos t and averaged over the period 2jt. By hypothesis, for small /?, B(i) and C(<) vary much more slowly with t than t itself. Therfore B(r), C(t), B(i) and C(t) are considered to be constant while performing the FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.4 29 averaging. Through this process the following equations are obtained "o 0 fp] 0 M 0 \ B + 0 0 M K | [P-Q] §R PQ K M PR - M K Is) * 5 * [ajaI + B)BU2 + C J C I / 2 ) + * 5 * lAiAl + B)BU2 + cviCU2) ? { H + J } = < (2.32) where b = { b ; 5; b>.} t, c = { c j c j c ; } T , are average values (over one period) of A,B(t) and C(t). Here the superscripts u, v,p indicate that the amplitudes are associated with u,v,p degrees of freedom respectively. The steady state solution of the problem corresponds to the singular points of the autonomous system of equations (2.32), that is, when B = C = 0. This results in a set of nonlinear algebraic equations for A, B , C. These equations are solved using the Newton-Raphson iteration procedure as outlined in Appendix 2. The equations obtained for the increments to the solution vector, following FUNDAMENTALS OF FLUID-STRUCTURE INTERACTION / 2.4 30 the procedure of Appendix 2, are of the form [J]{Ax} = {-/} (2.33) where [J] is a coefficient matrix, {Ax} is an increment to the solution vector and {—/} is a load vector. The details of these matrices corresponding to the nonlinear algebraic equations under consideration are also given in Appendix 2. The modified Newton-Raphson procedure was adopted to obtain the steady-state solution of equations (2.32). In this procedure the coefficient matrix [J] in equa-tion (2.33) is not computed after every iteration but is computed only after a specified number of iterations. C H A P T E R 3 S T A B I L I T Y O F S T E A D Y S T A T E S O L U T I O N S 3.1 I N T R O D U C T I O N In the present study a steady state solution was obtained to the fluid-structure interaction problem by assuming the solution to be of periodic nature with slowly varying amplitudes. Therefore, it was thought necessary to conduct a stability analysis of this solution to ascertain whether it is physically realizable or not. The general description of the governing equations for stability analysis are obtained by carrying out a perturbation analysis of a finite element representation of the two-dimensional unsteady Navier-Stokes equations. This is done by introducing a small perturbation to the steady state solution, and this results in a parametrically excited sys-tem. The nature of the small perturbations is examined by applying Floquet theory. This in turn requires a transition matrix to be determined, the eigenvalues of which characterise the nature of the small perturbations. The main limitation for numerical implementation of Floquet theory has been the computational effort necessary for evaluating the transition matrix for large systems. This evaluation requires n passes for the calculation of n lin-early independent solutions over the period T of the system. To overcome this deficiency Friedmann, Hammond and Woo[30] have developed a numerical scheme that yields the 31 STABILITY OF STEADY STATE SOLUTIONS / 3.2 32 transition matrix in one pass rather than n passes. A similar numerical scheme has been developed here. 3.2 G E N E R A L D E S C R I P T I O N O F T H E G O V E R N I N G E Q U A T I O N S In Chapter 2, an approximate steady state solution of equations (2.20) was obtained using the method of averaging. This steady state solution from equations (2.32) may be written as d = {uy Vj p y }r = A + B cos f. + C sin t (3.1) Once the steady state solution is calculated its stability is usually investigated by super-imposing a small perturbation on d, that is, by letting d = d + e (3.2) where the subscripts u y , vj, pj indicate the perturbations to the Uy , «y, py degrees of freedom, respectively. Substituting equation (3.2) into equations (2.20), using the fact that d is the steady state solution of equations (2.20) and linearizing equations in e, we obtain Me + Z{t)s = 0 (3.3) where Z(t) is a time dependent periodic matrix, that is, Z(t) = Z(t + 2rr). The details of Z(t) are given in Appendix 3. The boundary conditions to be applied are eUj.t eVj. = 0 on the outer boundary r „ as well as at the mean position of the moving body where velocities are specified. STABILITY OF STEADY STATE SOLUTIONS / 3.3 33 3.3 F L O Q U E T T H E O R Y Floquet theory is used for characterizing the functional behaviour of a system of linear ordinary differential equations with periodic coefficients as in equations (3.3) (see Nayfeh and Mook[17]). For the system in equations (3.3), it is possible to define a fundamental set of solutions where N is the order of matrices in equations (3.3) after applying the boundary conditions. This fundamental set of solutions can be expressed in the form of an N x N matrix [s] as W = [^ i e2 • • • eN] (3.4) It is clear that [e] satisfies the matrix equation M[£] + Z(t)[e] = 0 (3.5) Since Z(t) = Z ( * + 2tt), M[i(t + 2ic)] + Z{t + 2ir)[e(t + 2ir)] (3.6) = M[e(* + 2tt)] + Z{t)[e{t + 2ir)} = 0 from which it follows that [e(t + 2ir)] is also a fundamental set of solutions. Hence it is related to [s[t)\ by [i(t + 2*)] = •[*(«)] (3-7) where 9 is a nonsingular constant N x N matrix called the transition matrix. Introducing a transformation [7(01 = » W 0 1 (3-8) where 9 is a nonsingular constant N x N matrix, equations (3.7) can be expressed as [7(<.+ 2ir)l = l - ^ h f W l = T[7(*)l (3-9) STABILITY OF STEADY STATE SOLUTIONS / 3.3 34 The matrix can be chosen so that the matrix T has a Jordan canonical form. This form depends on the eigenvalues of 9 which are the N roots of the characteristic equation * - A [ f ] = 0 (3.10) When the roots of equation (3.10) are distinct, then the Jordan canonical form of T is a diagonal matrix given by A x 0 0 ... 0 0 A 2 0 ... 0 T = 0 0 A 3 . . . 0 0 0 0 . . . A j v . Consequently, writing equations (3.9) in component form, it follows that (3.11) (3.12) 7<(< + 2ir) = At-7,-(0 for i = 1,2,..., AT It follows from equations (3.12) that 7,-(* + n2ir) = A?7.-(0 (3.13) where n is an integer. As t — • oo (ie., n — » oo), it follows from equations (3.13) that From equations (3.13) it also follows that if A,- = 1, 7,- is periodic with period 2ir while if A,- = —1, 7,- is periodic with period 4tt. Therefore, stability of the steady state solution requires that + A}4 < 1 (3-15) where A# . and A/, are the real and imaginary parts of A,-, that is, A,- = Aj j . + iXj.. The boundedness criteria of equation (3.14) applies even when roots A,- of equation (3.10) are not distinct. For further information about this and a Jordan canonical form of T when roots A,- of equation (3.10) are not distinct, the reader is referred to the book by Nayfeh and Mook[17]. STABILITY OF STEADY STATE SOLUTIONS / 3.4 35 3.4 N U M E R I C A L D E T E R M I N A T I O N O F T H E T R A N S I T I O N M A T R I X 3.4.1 General Remarks. Of the various methods available for studying the stability of a system with periodic coefficients, Floquet theory is the most general one. Until very recently the major difficulty has been the computational effort required in evaluating the transition matrix. With the advent of modern computers and the development of efficient numerical treatment schemes, this difficulty has been overcome, to a certain extent. In the past, direct numerical integration has been carried out to evaluate the transition matrix (eg., Friedmann and Silverthorn[31]). This is done by using the initial conditions [e(0)j = [T\. Then the transition matrix 9 = [s(27r)] according to equa-tions (3.7). The main shortcoming of this approach has been that for matrix 9 of order N, N integration passes would be required with each integration pass yielding one column of Friedmann, Hammond and Woo[30] have suggested an efficient numerical integration scheme which yields the transition matrix * in a single pass without resorting to the application of initial conditions. Their integration scheme is based upon the fourth order Runge-Kutta scheme with Gill coefficients. In the present study a similar scheme is developed, but based on the trape-zoidal rule integration scheme. An implicit integration scheme was neccesary here since the matrix M in equations (3.5) is a singular matrix. The trapezoidal rule is chosen since it is the only member of the trapezoidal family of methods which is second order accurate and there is no time step restriction imposed by stability requirements (see Hughes[32]). STABILITY OF STEADY STATE SOLUTIONS / 3.4 36 3.4.2 A Single Pass Integration Scheme. U s i n g t h e t r a p e z o i d a l r u l e t o s o l v e e q u a t i o n s ( 3 . 3 ) c o n s i s t s o f t h e f o l l o w i n g s e t o f e q u a t i o n s : M v „ + 1 + Z B + 1 e n + 1 = 0 = e„ + A r v n + 0 . 5 ( 3 . 1 6 ) •n+o .5 = 0 . 5 v B + 0 . 5 v B + 1 w h e r e en a n d v n a r e a p p r o x i m a t i o n s t o e(tn) a n d £(tn) r e s p e c t i v e l y ; Z B = Z ( * B ) ; At i s t h e t i m e s t e p a s s u m e d c o n s t a n t i n t h e p r e s e n t c i r c u m s t a n c e s . C o m b i n i n g e q u a t i o n s ( 3 . 1 6 ) o n e o b t a i n s [ M + 0 . 5 A t Z B + 1 ] e B + 1 = [ M - 0 . 5 A i Z n ] e B ( 3 . 1 7 ) w h i c h c a n a l s o b e w r i t t e n i n t h e f o r m e B + 1 = D B e n ( 3 . 1 8 ) w h e r e t h e d y n a m i c m a t r i x D „ = D ( t B ) = [ M + 0.5AtZ^]-1^ - 0 . 5 A * Z „ ] U s i n g e q u a t i o n s ( 3 . 1 8 ) t h e f o l l o w i n g e x p r e s s i o n s c a n b e i m m e d i a t e l y w r i t t e n o u t s2 = D l £ l = D i D 0 e 0 ( 3 . 1 9 ) sn = D B _ i D n _ 2 • • • D 0 £ o F o r c o n s t a n t s t e p s i z e At = 2n/n a n d t0 = 0 a n d tn = 2ir, t h e l a s t o f e q u a t i o n s ( 3 . 1 9 ) r e d u c e s t o e(27r) = D ( 2 t t - A * ) D ( 2 * r - 2 A < ) • • • D ( 0 ) e ( 0 ) = f[D(2ir-ibA0e(0) ( 3 - 2 0 ) k=i U s i n g e q u a t i o n s ( 3 . 7 ) f o r a n y o n e f u n d a m e n t a l se t o f s o l u t i o n s a n d f o r t = 0 g i v e s e(2n) = * e ( 0 ) ( 3 . 2 1 ) STABILITY OF STEADY STATE SOLUTIONS / 3.4 37 Comparing equations (3.20) and equations (3.21) the transition matrix $ is given by n * = n D ( 2 , r ~ * A t ) (3-22) k=i where it is understood that the order in the product is important and the k th term is supposed to be before the k + 1 th term. After numerically evaluating the transition matrix using equations (3.22), the eigenvalues are evaluated using equations (3.10) and checked for the stability requirement given by equation (3.15) for t = 1,2,..., N. C H A P T E R 4 N U M E R I C A L I N V E S T I G A T I O N S 4 . 1 I N T R O D U C T I O N For a specified rigid body shape and specified values of frequency Reynolds number R^ and the body amplitude ratio B, the flow quantities such as the velocity and pressure can be obtained by determining the average values of the amplitudes A, B, C from equations (2.32) for a finite element discretization of the fluid domain. The steady state solution of the problem corresponds to the singular points of the autonomous system of equations (2.32), that is, when B = C = 0. This results in a set of nonlinear algebraic equations for A, B, C. These equations are solved using the Newton-Raphson iteration procedure as outlined in Appendix 2. A computer program in FORTRAN is developed to carry out the entire solution procedure. The flow solution is then used to determine the velocities at the mean position of the body using equations (2.27). This results in velocities with superharmonics as outlined in Section 4.2 later in the chapter. The overall view of the flow pattern at any instant of time can be obtained by plotting the velocity vectors in the fluid domain. Since the velocity solution has a steady part as well as a time dependent part, it is possible to obtain the stream function \P for the steady part of the velocity and, thus, obtain the flow pattern for the basic 38 NUMERICAL INVESTIGATIONS / 4.1 39 nonlinear phenomenon of steady streaming. This flow pattern is then compared with published experimental results. The finite element procedure for obtaining this flow pattern is outlined later in this chapter. A postprocessing package is developed in FORTRAN to compute the stream function for the steady part of the velocity using the finite element technique and to plot the stream function contours and the velocity vectors at any instant of time. Forces on structures are of central importance in engineering applications. The procedure to determine the resulting forces on the rigid body cylinder due to the fluid flow is described in some detail in this chapter. These forces are characterised by the concept of added mass, added damping and added force. A FORTRAN computer program is developed to evaluate these quantities. There are two main concerns regarding the numerical solution obtained for the present nonlinear fluid-structure interaction problem, namely, (1) the fluid flow is assumed to be periodic since the motion of the body is periodic and (2) the stable flow situation would be physically realizable but a numerically converged solution may not be stable, and therefore, physically unrealizable. Thus, a method for determining the stability of the numerically converged solution is developed as described in Chapter 3. A computer program in FORTRAN is developed to implement this procedure. All the computer programs are implemented on a 48-megabyte Amdahl 5840 system with double accelerator at the University of British Columbia. Double precision arithmetic is used throughout to reduce the effect of round off errors. The solution of the linear algebraic equations and the matrix inversions are performed using a sparse matrix solving package called SPARSPAK (University of Waterloo). The numerical results are obtained for three different body shapes, namely, (1) a square body oscillating parallel to one of its sides, (2) an oscillating circular body and (3) a symmetric Joukowski profile oscillating parallel to the line of symmetry. The latter case is considered to investigate the flow pattern around an asymmetrical body. In all NUMERICAL INVESTIGATIONS / 4.2 40 cases results are obtained for the basic nonlinear phenomenon of steady streaming, the in-stantaneous velocity vectors in the fluid domain, added mass, added damping, added force and stabilty of the flow. The numerical results are compared with published experimental observations wherever possible. 4.2 C O M P U T A T I O N O F T H E B O D Y B O U N D A R Y V E L O C I T Y The velocities at the mean position of the body are computed using equa-tions (2.27) after the flow solution has been obtained. From equations (2.30) the velocity at any node in the fluid domain other than at the mean position of the body consists of three components, namely, u y = A] + B) cost + C] sin t (4.1) Vj =A^ + B) cos t + C) sin t For one element on the body boundary as shown in Figure 2.6, equations (4.1) are substituted into equations (2.27) to obtain the nodal velocities at the mean position of the body. This results in the following equations: u{ = - i sCgCJ + cost - pCiiAvjs\nt + cos2t + X-p (-CyBJ - C a ) sin 2t (4.2) = " ^ C S C J - PC-^A) sint + jj8C sCJ cos2t - \pC^B) sin 2t where ; is summed over the velocity degrees of freedom other than those on the edge interfacing with the mean position of the body and k is summed over the velocity degrees of freedom on the edge interfacing with the mean position of the body. Note that the velocities in equations (4.2) have superharmonics and at t = 0, when the body is at its mean position, the nodal velocities correspond to the body velocity as they should. NUMERICAL INVESTIGATIONS / 4.3 41 4.3 D E T E R M I N A T I O N O F S T R E A M L I N E S F O R S T E A D Y P A R T O F T H E S O L U T I O N 4 . 3 . 1 G o v e r n i n g E q u a t i o n s . The A component of the solution in the equations (2.30) represents the steady flow part of the solution which has been referred to as the streaming flow in classical fluid mechanics (Schlichting [11]). It is possible to obtain the stream function \P for this steady flow simply by solving the Poisson's equation du dv , v * = -< = ay-Tx <4-3> where u, v are the steady flow components of the velocities and f is the steady flow component of the vorticity. The equation (4.3) is represented by the functional where fl is the domain of the problem and g(S) = is the tangential velocity specified on the natural part of the boundary Ct. 4 . 3 . 2 F i n i t e E l e m e n t D i s c r e t i z a t i o n a n d M a t r i x E q u a t i o n s In this work the finite element interpolation chosen for the stream function \& is the same as that used for the velocities in Subsection 2.2.3. This is done for ease of computation and ease of setting up the matrix equations. The stream function is represented by * = JVi*,- t = 1,2,3,..., 8 (4.5) where N{ are the shape functions given by Nx = -7(1 - a)(l - 0(1 + • + 0 JV2 = --(1 + «)(1 - 0(1 -» + t) 4 4 A 3 = ~ ( l + « ) ( l - « ) ( l - « - < ) A4 = - i ( l - « ) ( l + *)(! + • - ' ) * 5 = ^ ( i - s 2 ) ( i - 0 J\r6 = i ( i + S )( i-* 2 ) ^ 7 = ^ ( i - O ( i + 0 J\r8 = i ( i - S ) ( i - * 2 ) NUMERICAL INVESTIGATIONS / 4.4 42 and are the nodal variables. Equations (4.5) are substituted into the functional in equation (4.4). The functional is then minimised with respect to the nodal values for a single element. For g(S) = 0, the corresponding element equations are given by di where the arrays are = [%,] {*>} " {Pi) = {0} (4-6) t — 1,2,...,8 dA and fle denotes that the corresponding integration is taken over a single element only. All such element equations are assembled to produce a set of global equations pertaining to the fluid domain. These equations may be written in the matrix form as Z * = P (4.7) The solution of equations (4.7) substituted back into equation (4.5) gives the finite element representation of the stream function ^ in the fluid domain. 4.4 D E T E R M I N A T I O N O F A D D E D M A S S , A D D E D D A M P I N G A N D A D D E D F O R C E 4.4.1 General Remarks. A procedure to determine the fluid forces acting on the oscillating body is outlined here. For a dynamic system these forces are generally characterised by the concept of added mass, added damping and added force. The details of the concept and evaluation of the added mass are described in an excellent book by Sarpkaya and Isaacson [33]. The basic definition of the added mass is NUMERICAL INVESTIGATIONS / 4.4 43 7 FIGURE 4.1 Stress Components. given as the quotient of the additional force required to produce accelerations throughout the Quid divided by the acceleration of the body. The fluid forces on a rigid body performing harmonic oscillations will in gen-eral consist of three components, namely, (1) the component in phase with the acceleration of the body, (2) the component in phase with the velocity of the body and (3) a constant component. These components are each associated with the added mass, added damping and added force, respectively. The procedure for numerical determination of this quantities is outlined in this section. The nondimensional values of the added mass, the added damping and the added force are evaluated in the following sections for various body shapes and for a range of nondimensional parameters, the frequency Reynolds number and the body amplitude ratio B. 4.4.2 Determination of Fluid Forces. Knowing the stress components o~x1 <ryt at any point P in the two di-mensional fluid domain, the stress on any plane BC through this point and inclined at an angle 9 to the x-axis (see Figure 4.1) can be calculated from the equations of statics and are given by o — ox sin2 9 + <ry cos2 9 — 2tx1) sin 9 cos 9 (4.8) r = r,y(sin2 9 - cos2 9) + (<r, - ay) sin 9 cos 9 NUMERICAL INVESTIGATIONS / 4.4 44 where a and r are the normal and shear components of stresses, respectively, on the plane BC. In matrix form equations (4.8) can be written as f a i r sin20 cos20 -2sin0cos0 \ t J — [sin 0 cos 0 — sin 0 cos 0 sin2 8 — cos2 6 •xy } (4.9) The velocities and pressure are represented by the finite elements by u = AT(«t- t = l , 2 , ...,8 v = NiV{ t = l , 2 , . . . , 8 (4.10) p = MiPi » ' = 1 , . . . , 4 where iV,-,M,- are the shape functions, given in Subsection 2.2.3 and u,-,v,-,p,- are the nodal variables. For an isoparametric element the coordinates of a point in the element are also given by x = JV,x,- i = 1,2,..., 8 (4.11) ,- , y = iV,y,- t = l , 2 , . . . , 8 where x,-,yt- are the nodal coordinates. From equations (4.10), we obtain du dNi dx Ui dNi dv dNi dx ~ du dx ~ dv dx dNi dy~ 9y~ dy Vi (4.12) Using equations (4.11) and the chain rule, the derivatives of the shape functions are given by dNi_dNidx dNi dy da dx ds dy ds _ dNi^dN^ dNidNj dx ds Xj + dy ds Vj dNi_dNidx dNi dy dt dx dt+ dy dt dK dNj dNi dNj — * -X • A - i-tl . dx dt 3 dy dt V] (4.13) or in matrix form I dN,-ON,- ' d'Ni dNi I ~dTx) -dTvj (4.14) NUMERICAL INVESTIGATIONS / 4.4 7 ' A 45 F I G U R E 4.2 and alternatively written as Cross Section of Body. i f H t ) (4.15) where the matrix J is the Jacobian of the transformation and is given by the the corre-sponding matrix in equations (4.14). From equations (4.15), the following are obtained ( H M D (4.16) The equations (4.10) are substituted into equations (2.16) to obtain 0 -Mi 2 ^ ~Mt oy ox IS) (4.17) Therefore, to obtain the normal and shear components of the stresses acting on plane BC at point P in the fluid domain, first the derivatives of the shape functions are evaluted using equations (4.16), then the stress components ax, ayt rxy are evaluated using equations (4.17) and lastly equations (4.9) are used. This whole process can be implemented in a computer program using matrix manipulations after the solution for the velocities and pressure in the fluid domain have been obtained. For a body of general cross section (see Figure 4.2) and axial length /, the resultant fluid force in the x-direction acting on it is given by Ff = -I j{osm9+ rcos9)d$ (4.18) NUMERICAL INVESTIGATIONS / 4.4 46 where $ is the coordinate along the boundary and 0 is the angle which S makes with the positive x-axis. Introducing non-dimensional variables b fiu0/b - Ff - r f iiu0l fwo/b into equation (4.18), we obtain Ff = -j{asmO + r cos 0)dS (4.19) where all quantities are their respective non-dimensional values and the overtilda has been dropped for convenience. In general the resultant fluid force Ff on a body oscillating with velocity u0coswi in the fluid domain will consist of three components, namely, (1) the component in phase with the acceleration of the body, (2) the component in phase with the velocity of the body and (3) the component which is constant. This is so because the solution for the velocities and pressure in the fluid domain consists of a steady as well as oscillating components. Therefore, the fluid force Ff obtained from equations (4.19) after integration would be of the following form Ff = P + Qcost + Raint (4.20) where t is the non-dimensional time and P. Q, R are evaluated using equation (4.19). It should be noted here that for a body with at least one axis of symmetry, the resultant fluid force Ff in the direction perpendicular to the axis of symmetry will have P = 0. 4.4.3 Computation of A d d e d Mass, Added Damping and Added Force. The equation of motion for the body is given by equation (2.17). The non-dimensional displacement of the body, s = Bs'mt. and the fluid force Ff from equa-tion (4.20) are substituted into equation (2.17) to obtain Fe = (^^-R) Bint-Q cost+ ^smt-P (4.21) V Jb J Jb NUMERICAL INVESTIGATIONS / 4.5 47 The equation of motion can also be written in the form {m + ma)s + cas + ks + fa = JbFe (4.22) where ma, ca and fa are the non-dimensional values of the added mass, added damping and added force, respectively. The non-dimensional displacement of the body, s = fi sin t, is substituted to obtain F e = f = ^ _ ! ^ ) s i n ( + ^ C 0 S , + M s i n ( + £ ( 4 . 2 3 ) \ Jb Jb J Jb Jb Jb Comparing equation (4.21) and equation (4.23), we obtain RJb m * = fi Q Jb (4.24) fa = -PJb Therefore, to evaluate the non-dimensional values of the added mass, added damping and added force, we first carry out the numerical integration as indicated in the Subsection 4.4.2 to obtain P, Q and R and then use equations (4.24). The added force is zero for a body which has at least one axis of symmetry and the resultant fluid force Fj is evaluated perpendicular to the axis of symmetry, since P = 0 for such a body. 4.5 RESULTS FOR O S C I L L A T I N G SQUARE B O D Y 4.5.1 General Remarks. Tatsuno[34] has presented excellent experimental results for an oscillating square body, and these results are used for the verification of the present method. The results of Tatsuno are for the basic nonlinear phenomenon of steady streaming which are compared with the numerically obtained stream function for the steady component of the flow. NUMERICAL INVESTIGATIONS / 4.5 48 The numerical'results are obtained for the square body having sides of unity with the centroid of its area located at the origin of the x-y coordinate system and sides NUMERICAL INVESTIGATIONS / 4.5 49 parallel to the axes. The body is performing harmonic oscillations parallel to the x-axis. Using symmetry, only one half of the domain is modelled. Two finite element grids shown in Figure 4.3 are developed. For the flow situation, the boundary conditions are specified by the velocities u = v = 0 at the outer boundary and the velocities specified according to equations (2.27) at the body boundary. The symmetry condition is obtained by letting the velocity v = 0 and the shear stress rxy = 0 on the symmetry line. For a finite element scheme only kinematic boundary conditions need to be specified and the homogenous natural boundary conditions come out as a part of the solution in the limit of grid refinement. Therefore, on the symmetry line only v = 0 is enforced. The pressure needs to be specified only at one location in the domain to obtain the pressure datum. In the present case the pressure distribution is symmetric, therefore the pressure is specified to be zero at the two symmetric locations, namely, (1) and (2) shown in Figure 4.3. For the stability analysis, the velocity perturbation boundary conditions are given by sn = ev = 0 at the outer boundary and on the body boundary since the velocities are specified at these locations. Since only one half of the domain is modelled, there are two cases to be investigated, namely, the symmetric and antisymmetric modes of stability. The symmetric case is obtained by forcing the velocity perturbation ev = 0 on the symmetry line and the antisymmetric case is obtained by forcing the velocity perturbation en = 0 on the symmetry line. The pressure perturtion ep — 0 at the same two locations (1) and (2) shown in Figure 4.3. The finite element grid shown in Figure 4.3(a) is developed to perform pre-liminary calculations and to debug the program. The grid has 16 elements and 65 nodes. There are 121 net degrees of freedom, 98 of which are for the velocities and 23 for the pressure. Out of the 121 net degrees of freedom, 18 velocity degrees of freedon lie on the body boundary and are specified from the moving body boundary condition. For the flow situation, there are three types of coefficients A, B, C and therefore this results in 309 NUMERICAL INVESTIGATIONS / 4.5 50 variables in total. For the stability analysis, the total number of variables is 103. The re-sults from this coarse grid seemed to exhibit the correct trends for both the flow situation and the stability analysis, but the accuracy was rather poor. The finite element grid shown in Figure 4.3(b) is developed to represent Tatsuno's experimental setup and has ^ = 30, where DQ is the diameter of the outer fixed cylinder and 6 is the characteristic length of the body. The grid has 40 elements and 147 nodes. There are 294 net degrees of freedom, 242 of which are for the velocities and 52 for the pressure. Out of the 294 net degrees of freedom, 34 velocity degrees of freedon lie on the mean position of the body boundary and are specified from the moving body boundary condition. For the flow situation, there are three types of coefficients A, B, C and this results in 780 variables in total. For the stability analysis, the total number of variables is 260. The flow results are presented in two ways, namely, the instantaneous veloc-ity vectors in the flow field and the streamline contours for the steady part of the velocity. The streamline contours obtained from the steady component A of the solution will here-after be referred to as steady streamlines. The numerical steady streamlines are compared with the corresponding experimental streamlines. The added mass and added damping quantities are tabulated for all the flow situations considered. The added force is zero for the present case since the fluid force Ff is evaluated in the direction perpendicular to the axis of symmetry of the body. The results for the linear case, that is for the amplitude parameter (3 = 0, are plotted for the whole range of frequency Reynolds number Rw. The symmetric mode stability analysis for one particular flow situation is performed for the whole range of time step intervals using the grid in Figure 4.3(a). This is done to check the numerical scheme, to observe the trends with the decreasing time step interval and to obtain guidelines for a suitable time step interval for the more refined grid. The results for this particular flow situation are also obtained using the grid in NUMERICAL INVESTIGATIONS / 4.5 51 K fi Location of Maximum Stream Function Maximum Steady Stream Function X y 1.99 0.06 1.2 1.05 0.02247 1.99 0.10 1.2564 1.0636 0.02012 26.00 0.06 1.0636 0.8218 0.03694 26.00 0.10 1.16 0.8218 0.02971 120.00 0.08 -2.7 2.5 0.05744 120.00 0.10 -2.3684 2.4164 0.02996 T A B L E 4.1 Maximum Steady Stream Function for Oscillating Square Body. Figure 4.3(b). One more flow situation is investigated using the grid in Figure 4.3(b), where the stability analysis is performed both in symmetric and antisymmetric modes. 4.5.2 Flow Results. The numerical and experimental streamlines are shown in Figures 4.4, 4.5 and 4.6 for values of of 1.99, 26, 120, respectively. In each case, the numerical results are shown for two values of the body amplitude parameter 6, whereas the experimental result is for only one of these values. The numerical results are obtained by plotting 31 equally spaced contours between the maximum and minimum value of the steady stream function in the plotted domain. For each case the location and the value of the maximum stream function is described in Table 4.1. Since the body in the present case is symmetric about the y-axis the value of the minimum steady stream function is the negative of the maximum value and is located at the symmetric point about the y-axis. Tatsuno[34] obtained his experimental pictures by a stroboscopic technique which effectively eliminated the periodic part of the flow. Tatsuno showed that there are three distinct flow regimes with different flow patterns and these figures show one example in each regime. The flows clearly become more complex as Rw increases. The example shown in Figure 4.4 has been classified as Type I by Tatsuno. In this type of flow there is one main vortex system in each quadrant. In this case, the streaming flow does not separate at the edge of the body. The flow pattern of this vortex NUMERICAL INVESTIGATIONS / 4.5 52 NUMERICAL INVESTIGATIONS / 4.5 53 NUMERICAL INVESTIGATIONS / 4.5 54 NUMERICAL INVESTIGATIONS / 4.5 55 system is symmetric with respect to both the axis of oscillation and the axis perpendicular to it. Overall the numerical and experimental results agree very well. However, there are some small differences. In Figure 4.4, the experimental result which is for fi = 0.1 compares better with the numerical result for fi = 0.06. The numerical result for fi = 0.1 shows the start of another pair of vortices above (and below by symmetry) the body which are not evident in the experimental result. The example shown in Figure 4.5 has been classified as Type II by Tatsuno. In this example the flow separates from the square body at the edge and the twin vortices are generated enclosed between the primary vortex systems and the square body sides parallel to the oscillation direction. The twin vortices due to the flow separation have been observed to grow with increasing value of Ru or the amplitude parameter fi and ultimately the flow pattern shifts to that of Type HI. Overall the numerical and experimental results agree very well for this case, too. However, there are some small differences. In Figure 4.5, the experimental result which is for fi = 0.06 compares better with the numerical result for /? = 0.1. The numerical result for fi = 0.06 does not show the appearance of the secondary twin vortices. As the amplitude parameter fi is increased to 0.1 the secondary twin vortices appear but the location of the primary vortices remains essentially the same and this case compares better with the experimental result. The example shown in Figure 4.6 has been classified as Type i n by Tatsuno. There are two vortex systems close to the body in each quadrant and flow separates at the edge of the square body. Unlike the flow pattern of Type II, this example does not enclose the second set of twin vortices between the primary vortices and the sides parallel to the direction of oscillation. It has been observed that the two vortex systems in each quadrant become small with the increasing value of Ru. Overall the numerical and experimental results agree very well in this case too, though, small discrepancies do appear. In Figure 4.6, the experimental result which is for fi = 0.08 compares better with the numerical result for fi = 0.1. The numerical result at fi = 0.06 does not clearly show the appearance of the second set of twin vortices, but both the vortex systems are clearly visible at fi = 0.1. NUMERICAL INVESTIGATIONS / 4.5 56 fi Location of Maximum Steady Velocity umax x y «o 1.99 0.06 0.625 0.625 0.081 1.99 0.10 0.5 0.25 0.282 26.00 0.06 0.625 0.3125 0.177 26.00 0.10 0.5 0.25 0.565 120.00 0.08 0.5 0.25 0.583 120.00 0.10 0.5 0.25 0.555 TABLE 4.2 Maximum Steady Component of Velocity for Oscillating Square Body. The photographs for the experimental results in Figures 4.4, 4.5 and 4.6 appear in Reference[34] and the original photographs have been obtained directly from Tatsuno. The steady flow components illustrated in these figures are called secondary flows, because they are expected to be small compared to the main periodic flow. For example, the forced single degree of freedom system with quadratic nonlinearities (same as the Navier-Stokes equation) described in the book by Nayfeh and Mook[17] show the steady part of the response to be of the order of the perturbation parameter (fi in the Navier-Stokes equations) compared to the primary harmonic response. The actual maximum of the resultant of the steady component of the velocities and its location in the flow domain are shown in Table 4.2. Since the body is symmetric about the y-axis these maximum velocities occur at two locations symmetric about the y-axis. These maximum velocities are obtained only from a check of nodal values, so they may not be the global maxima. The steady state solution obtained from equations (2.32) can be written as d = {«_,• vj pj}T = A + Bcosi + Csini (4.25) Some of the results obtained from equations (4.25) are shown in Figures 4.7, 4.8 and 4.9 for of 1.99, 26 and 120, respectively. For all cases, the body amplitude parameter fi = 0.1. These figures show velocity vectors over a uniform mesh in the vicinity of the body for t = 0, JT/4, TT/2, 3TT/4, and ir, that is equal time steps over one-half cycle of body NUMERICAL INVESTIGATIONS / 4.5 57 t = 0 scale = 0.26 •* * * ** -* > t = | aco/e = 0.36 ' '• ::::::::::: * = § sca/e = 0.66 * = 3* 5 c a / e = o.30 t = ir sca/e = 0.26 FIGURE 4.7 Velocity Vectors for Oscillating Square Body, = 1.99, /? = 0.1, ^ = 30. NUMERICAL INVESTIGATIONS / 4.5 58 \ 1 \ UJ -> 5S = ? * = 0 scale = 0.23 t = | sca/e = 0.28 i = | sca/e = 0.36 -i = 2f aco/e = 0.24 > / > > * * . . t = ir scale = 0.23 FIGURE 4.8 Veiocity Vectors for Oscillating Square Body. i? w — ^ = 30. 26, 8 = 0.1, NUMERICAL INVESTIGATIONS / 4.5 59 t = 0 scale = 0.24 t fc f * = f scale = 0.30 .' f - r t > ; •5 : • * = | aca/e = 0.49 t = sca/e = 0.26 i i ' r : : r n ; t = ir scale = 0.24 FIGURE 4 . 9 Velocity Vectors for Oscillating Square Body. R„ = 120, B = 0.1, ^ = 30. NUMERICAL INVESTIGATIONS / 4.5 60 motion. The size of each vector indicates the velocity magnitude. Note that the scale in each illustration is different. The magnitude of the velocities is obtained by multiplying the length of the arrows in centimetre by 1.5 and dividing by the scale. It is clear from these figures that there is one main vortex in the upper half of the fluid flow and although it does move slightly from side to side it is located just above (and below by symmetry) the oscillating body. As the frequency Reynolds number RU increases, this vortex moves closer to the body and decreases in size, thereby increasing the velocity gradients in that region. Consequently, the boundary layer thickness is also decreased. It should be noted that the velocities at the mean position of the body shown in Figures 4.7, 4.8 and 4.9 do not include the superharmonics of equations (4.2) since they were not computed at the time these figures were plotted. But this will not affect the velocities in the fluid domain except at the nodes located at the mean position of the body and that too to a minor extent. Picard iteration* was used for solving the nonlinear algebraic equations in the foregoing calculations. I t seems to work reasonably, well. Al l but one case, namely, for = 120 and fi = 0.08, were successful using the linear solution as the initial guess. The number of iterations required gradually increased from 4 to 10 as R„ and fi were increased. The aforementioned case for R^ = 120 and fi = 0.08 did not converge when a linear solution was used to start it. However, a converged result is obtained by using steps from RU = 26 and fi = 0.06 to i£w = 90 and fi = 0.08 and subsequently to the desired result. Overall the present method of representing the interaction between a solid body and viscous flow seems to work well and good agreement is obtained between the numerical and experimental results. However, the comparison between the experimental t This is one of the methods of successive linearization. Here the gness of an init ia l solution is made, the nonlinear terms computed and moved over to right hand side, and then the resulting linear equations are solved to get a new solution. This process is repeated until the converged solution is obtained. NUMERICAL INVESTIGATIONS / 4.5 61 and the numerical steady streamlines show small differences. It is generally observed that the secondary vortices seen in the experimental flow do not occur at exactly the same value of the body amplitude parameter /? for the numerical case, but occurs at a value close to it. Furthermore, the flow becomes more complex with increasing R^. Al l of this indicates the need for further grid refinement near the body before the secondary vortices appear at the same /3 as the experiment and before R^ can be increased further. The present computing facilities at the University of British Columbia do not allow for much further grid refinement due to the enormous memory requirements. 4.5.3 Added Mass and Added Damping. The results for the added mass m 0 and added damping ca are obtained using both the finite element grids in Figure 4.3. The grids in Figure 4.3(a) and Figure 4.3(b) will henceforth be referred to as Grid 4.3(a) and Grid 4.4(b), respectively. First, the results are obtained for the linear case, that is, when the body amplitude parameter /? = 0. Due to the nonavailability of any analytical results to make a direct comparison, finite element results for ma obtained by Irani[35] are made use of for this purpose. Irani used the stream function formulation of the Navier-Stokes equations and the 18 degrees of freedom quintic element for the finite element interpolation of the fluid domain. The comparison is also made with the limiting case of inviscid flow. The added mass and added damping values for a bluff body oscillating at a high frequency with a small amplitude are mainly given by inviscid flow considerations and are very little affected by viscosity. Hence the values of ma and ca for large values of iZ w and /? = 0 should approach the values predicted by inviscid flow theory. The values of the added mass and added damping for large values of Ru and p = 0 for both grids are presented in Table 4.3. For R^ = 10000, Grid 4.3(a) gives ma = 1.2655 and Grid 4.3(b), 1.2371 compared to 1.19 for inviscid flow[33]. Theoretically, in the limit of very high frequencies ma should approach the value predicted for inviscid flow. The limiting value ma = 1.24 predicted using Grid 4.3(b) and ma = 1.27 predicted NUMERICAL INVESTIGATIONS / 4.5 62 Rw Added Mass, ma Added Damping, ca Grid 4.3(a) Grid 4.3(b) Grid 4.3(a) Grid 4.3(b) 500 5000 10000 1.2749 1.2656 1.2655 1.3163 1.2379 1.2371 0.0614 0.0064 0.0032 0.2567 0.0302 0.0151 T A B L E 4.3 Added Mass, ma and Added Damping, ca for Oscillating Square Body for Ru large and (3 = 0. Added Mass, m 0 Grid 4.3(a) Grid 4.3(b) Irani Results l 6.2587 7.5857 8.05 3.384 4.0893 4.5531 6.05 21.34 2.0595 2.4472 3.3 43 1.7026 1.9926 2.4 100 1.4579 1.6506 1.645 144 1.3886 1.5531 1.383 177 1.3577 1.5047 278.2 1.3103 1.4091 324 1.2998 1.3809 T A B L E 4.4 Added Mass, m 0 for Oscillating Square Body for /3 = 0. using Grid 4.3(a) are within 4% and 6.5%, respectively, of the inviscid value of 1.19. The value for the added damping, ca appproaches zero for both the grids as it should from inviscid flow considerations. This indicates that the accuracy and representation of the flow phenomena for the linear case is very good, and results improve with grid refinement as they should. For moderate values of R^, the added mass for the linear case obtained using Grids 4.3 and those obtained by Irani[35] are presented in Table 4.4 and graphically in Figure 4.10. The values of ma obtained using Grid 4.3(a) are consistently lower than that obtained using Grid 4.3(b) except in the limiting case of very large R^, where the result for Grid 4.3(b) is closer to the inviscid value and lower than for Grid 4.3(a). The values of ma obtained using Grid 4.3(b) are close to the results obtained by Irani except for the limiting case of very large Rw, where Irani's result appproaches unity rather than 1.19 for NUMERICAL INVESTIGATIONS / 4.5 O 63 © CO o o to W CO a S o © © I LEGEND v = Grid 4.3(a) • = Grid 4.3(b) o = I r a n i Results — # = 9 — i r — i 1 ' I r " — I ' 1 • i 0.0 2.0 4 .0 6.0 8.0 10.0 12.0 — | . 1 1 1 r 14.0 16 .0 18 .0 2 0 . 0 F I G U R E 4.10 Added Mass, m. for Oscillating Square Body for B = 0. NUMERICAL INVESTIGATIONS / 4.5 64 1 3.384 21.34 43 100 144 177 278.2 324 Added Damping, ca Grid 4.3(a) 11.2686 4.8642 1.5425 0.9663 0.5143 0.4154 0.3541 0.2425 0.2116 Grid 4.3(b) 11.2261 4.9243 1.6652 1.1254 0.6825 0.5528 0.4928 0.3813 0.3471 T A B L E 4.5 Added Damping, ca for Oscillating Square Body for 3 = 0. the inviscid case. For moderate values of R^, the added damping for the linear case obtained using Grids 4.3 are presented in Table 4.5 and graphically in Figure 4.11. The values of ca obtained using Grid 4.3(a) are consistently lower than that obtained using Grid 4.3(b) except for = 1. For the limiting case of large Ru, the results from both grids approach zero, same as for the inviscid case. The results for the added mass and added damping presented in Table 4.6 are for finite values of body amplitude parameter B. For this case there are no published results to compare with. The results in Table 4.6 show a general trend of decreasing ma and ca with increasing R^, same as for the linear case. Also, it is observed that the values of ma and c„ show only a small change with a small change in B and values are close to the values for the linear case for the B values considered here. The results presented in Table 4.6 are for the flow situations described earlier in Subsection 4.5.2 except for the case of R^ = 566.44. The numerical solution for i? w = 566.44 and B = 0.058 did not converge using the Picard iteration even though several different initial guesses were tried. However, it did converge after 12 iterations when the full Newton-Raphson solver was used with the solution for Ru = 566.4 and B = 0.0376 as an initial guess solution. NUMERICAL INVESTIGATIONS / 4.5 O N - I F I G U R E 4 .11 Added Damping, c« for Oscillating Square Body for p = 0. NUMERICAL INVESTIGATIONS / 4.5 66 Grid 4.3(b) P Added Mass, Added Damping, ca 1.99 0.06 5.3975 7.1963 1.99 0.1 5.6237 6.7802 26.01 0.06 2.1886 1.7092 26.01 0.1 2.3074 1.2300 120 0.08 1.6697 0.2790 120 0.1 1.6410 0.4115 566.44 0.0376 1.3078 0.2346 566.44 0.058 1.3140 0.2468 T A B L E 4.6 Added Mass, ma and Added Damping, ca for Oscillating Square Body. Overall, the results for the added mass and added damping for the linear case of P = 0 compare very well with Irani's finite element results and approach the correct limiting values of the inviscid flow. The results also exhibit the correct trends for finite values of body amplitude parameter p. Al l these indicate that the procedure for computing the added mass and added damping is correct and the finite element representation of the fluid domain and the Taylor series expansion of the nonlinear boundary conditions adequately represent the basic characteristics of the physical phenomena. 4.5.4 S t a b i l i t y Analysis. For the stability analysis the same two grids are used. Grid 4.3(a) has 23 pressure degrees of freedom and 80 velocity degrees of freedom for both the symmetric and antisymmetric modes of stability analysis. Grid 4.3(b) has 52 pressure degrees of freedom and 208 velocity degrees of freedom for both the symmetric and antisymmetric modes. Preliminary calculations are performed using Grid 4.3(a). The symmetric mode stability analysis is carried out for the oscillating square body for i 2 w = 1.99 and p = 0.1 and the results are presented in Table 4.7 for the top 52 eigenvalues in descending order of magnitudes. The flow situation considered here has also been observed experimentally (see Subsection 4.5.2) and, therefore, it is expected to be stable. The results show that there are 46 eigenvalues (twice the number of pressure degrees of freedom) with magnitude equal NUMERICAL INVESTIGATIONS / 4.5 67 Grid 4.3(a) Eigenvalue No. At=*f 1*1 1*1 At — — |A| At — — |A| |A| . . so |A| 1 1 1 1 1 1 1 2 1 1 1 1 1 1 3 1 1 1 1 1 1 46 1 1 1 1 1 1 47 0.84714 0.51039 0.22444 0.021145 0.021158 0.021167 48 0.84527 0.51492 0.22002 0.011728 0.011739 0.011747 49 0.83926 0.49601 0.20630 0.011375 0,011386 0.011394 50 0.80885 0.42786 0.14784 0.001679 0.001684 0.001687 51 0.80592 0.42170 0.14309 0.000219 0.000220 0.000221 52 0.80072 0.40445 0.13495 0.000031 0.000031 0.000031 TABLE 4.7 Symmetric Mode Stability Analysis for Oscillating Square Body. R„ = 1.99, B = 0.1. to 1 and the rest are all less than 1. This, according to the criteria given by equation (3.15), indicates that the numerical steady state solution is stable for the symmetric perturbation. As the non-dimensional time step interval At is reduced from At = ^ to A* = 1^ , it is observed that the magnitude of the eigenvalues (other then |A| = 1) decrease and converge to a constant value at At = ^ onwards. The converged results show that the magnitudes of the eigenvalues other than |A| = 1 are much less than 1, qualitatively indicating that this is a strongly stable situation. The results presented in Table 4.8 are for the same flow situation as in the previous case. However, the calculations are performed using Grid 4.3(b) and results are presented for the top 115 eigenvalues in descending order of magnitude. The time step for convergence is expected to be much smaller for Grid 4.3(b) as compared to the coarser Grid 4.3(a) (see Friedmann et al [30]). The calculations are performed for At = ^ and At = only, because of the large CPU time requirements. The results show that there are 104 eigenvalues (twice the number of pressure degrees of freedom) with magnitude equal to 1 and the rest are all less than 1. This, according to the criteria given by equation (3.15) NUMERICAL INVESTIGATIONS / 4. 5 68 Eigenvalue No. Grid 4.3(b) |A| At - 2K . . 4 0 |A| 1 1 1 2 1 1 3 1 1 104 1 1 105 0.87084 0.58450 106 0.86875 0.55128 107 0.80963 0.45520 108 0.77856 0.36658 109 0.77583 0.28828 110 0.77286 0.24335 111 0.76806 0.19939 112 0.75991 0.17017 113 0.74788 0.11284 114 0.73748 0.10934 115 0.73032 0.10521 T A B L E 4.8 Symmetric Mode Stability Analysis for Oscillating Square Body. Ru = 1.99, p = 0.1. indicates that the numerical steady state solution is stable for the symmetric perturbation. As the non-dimensional time step interval is reduced from At = |£ to At = ^ it is observed that the magnitude of eigenvalues (other than |A| = 1) decreases, similar to the previous case. The results for the top 115 eigenvalues in descending order of magnitude are presented in Table 4.9 and Table 4.10 for the symmetric and antisymmetric perturbation modes, respectively. The flow situation considered is = 120 and P = 0.1. This flow situation is very close to the flow situation i ? w = 120 and P = 0.076, where the numerically converged results could not be obtained. The results show that there are 104 eigenvalues (twice the number of pressure degrees of freedom) with the magnitude equal to 1 and rest are all less than 1. This, according to the criteria given by equation (3.15) indicates that the NUMERICAL INVESTIGATIONS / 4.5 69 Grid 4.3(b) Eigenvalue No. At — — |A| 1*1 1 1 1 2 1 1 3 1 1 104 1 1 105 0.99080 0.99080 106 0.99015 0.99015 107 0.98694 0.98694 108 0.98320 0.98320 109 0.97966 0.97966 110 0.97641 0.97641 111 0.97332 0.97332 112 0.97104 0.97104 113 0.96457 0.96457 114 0.95692 0.95692 115 0.95261 0.95261 TABLE 4.9 Symmetric Mode Stability Analysis for Oscillating Square Body. Ru = 120, 0 = 0 .1 . numerical steady state solution is stable for both symmetric mode and the antisymmetric modes. The results also show that the eigenvalues have essentially converged at At = and further reduction of the time step interval to At = ^ brings about no change in them. Many of the eigenvalues in the present case have magnitudes close to 1 which qualitatively indicates a weakly stable situation. This may be due to the fact that the flow under consideration is very close to the flow for which a numerically converged steady state solution was not obtained, and therefore, this may represent a bifurcation point for transition to different flow patterns. For all the cases considered here, the number of eigenvalues with magnitude 1 are twice the number of pressure degrees of freedom. No explanation for this can be offered at this time, but it may be due to the fact that for incompressible flow the mass NUMERICAL INVESTIGATIONS / 4 15 70 Grid 4.3(b) Eigenvalue No. |A| |A| 1 1 1 2 1 1 3 1 1 104 1 1 105 0.99605 0.99605 106 0.99007 0.99007 107 0.99005 0.99005 108 0.98724 0.98724 109 0.98592 0.98592 110 0.97984 0.97984 111 0.97664 0.97664 112 0.97421 0.97420 113 0.97182 0.97182 114 0.97014 0.97014 115 0.96493 0.96493 T A B L E 4.10 Antisymmetric Mode Stability Analysis for Oscillating Square Body. R„ = 120, p = 0.1. matrix M is rank deficient by the number of pressure degrees of freedom because the continuity equation does not involve the pressure term explicitly. Overall, the implicit numerical integration scheme using the trapezoidal rule to obtain the transition matrix seems to work very well. The convergence to the transition matrix and thereby to the eigenvalues is obtained by increasing the number of time steps (that is, decreasing the time step interval) and the scheme is unconditionally stable with respect to the size of the time step. The computer time required for obtaining the converged results is very high because the number of matrix inversions required are equal to the number of time steps. Therefore, any increase in the number of time steps would involve a substantial increase in the computer time required. NUMERICAL INVESTIGATIONS / 4.6 71 4.6 RESULTS FOR O S C I L L A T I N G C I R C U L A R B O D Y 4.6.1 General Remarks. Tatsuno[36,37] has presented excellent experimental results for the oscillating circular body which are used for the verification of the present method. The results of Tatsuno are for the basic nonlinear phenomenon of steady streaming and are compared with the numerically obtained stream function for the steady component of the flow. The numerical results are obtained for a circular body of characteristic length 6 of unity (that is, diameter = 1 for the present case) with its centroid located at the origin of the x-y coordinate system. The body is performing oscillations parallel to the x-axis. Using symmetry, only one-half of the domain is modelled. Three finite element grids shown in Figure 4.12 are developed. The grids in Figures 4.12(a) and (b) are developed to represent Tatsuno's experimental setup. For the flow situation, the boundary conditions are specified by the velocities u = v = 0 at the outer boundary and the velocities specified according to equations (2.27) at the body boundary. The symmetry condition is again obtained by setting the velocity v = 0 on the symmetry line. The pressure needs to be specified only at one location in the domain to obtain the pressure datum. In the present case the pressure distribution is symmetric, therefore the pressure is specified to be zero at the two symmetric locations, namely, (1) and (2) shown in Figure 4.12. For the stability analysis , the boundary conditions are given by the velocity perturbations e„ = ev = 0 at the outer boundary and the body boundary since the velocities are specified at these locations. Since only one half of the domain is modelled, there are two cases to be investigated namely, the symmetric and antisymmetric modes of stability. The symmetric case is obtained by setting the velocity perturbation ev = 0 on the symmetry line and the antisymmetric case is obtained by the velocity perturbation NUMERICAL INVESTIGATIONS / 4.6 FIGURE 4 . 1 2 Finite Element Grids for Circular Body. NUMERICAL INVESTIGATIONS / 4.6 73 eu = 0 on the symmetry line. The pressure perturbation ep = 0 at the same two locations (1) and (2) shown in Figure 4.12. The three finite element grids shown in Figure 4.12 all have a different ^ ratio, where DQ is the diameter of the outer fixed cylinder and 6 is the characteristic length of the body (diameter in the present case of the circular body). The grids in Figure 4.12(a), Figure 4.12(b) and Figure 4.12(c) have ^ = 7.7, ^ = 15.5 and ^ = 30, respectively. Al l these grids have 48 elements and 173 nodes. There are 351 net degrees of freedom, 290 of which are for the velocities and 61 for the pressure. Out of the 351 net degrees of freedom, 34 velocity degrees of freedon lie on the body boundary and are specified from the moving body boundary condition. For the flow situation, there are three types of coefficients A , B, C and therefore this results in 951 variables in total. For the stability analysis, the total number of variables is 317. The flow results are presented in two ways, namely, the instantaneous veloc-ity vectors in the flow field and the streamline contours for the steady part of the velocity. The streamline contours obtained from the steady component A of the solution will here-after be referred to as steady streamlines. The numerical steady streamlines are compared with the corresponding experimental streamlines. The added mass and added damping quantities are presented for a range of p and i? w values including the linear case of /? = 0. The added force is zero for the present case since the fluid force Ff is evaluated in the direction perpendicular to the axis of symmetry of the body. The stability analysis for the symmetric and antisymmetric modes is per-formed only for one flow situation. More extensive investigation has already been carried out for the square body and the same characteristic trends are expected for the circular body. Furthermore, the stability computations are very expensive. 4.6.2 Flow Results. NUMERICAL INVESTIGATIONS / 4.6 74 Location of Maximum Maximum Steady K P A ) Stream Function Stream Function b X y -3.834 0.1 7.7 -0.9 0.9 0.0044 3.834 0.13 7.7 -0.9 0.9 0.0059 3.834 0.161 7.7 0.45 0.3 0.0024 21.34 0.0925 15.5 -0.45 0.525 0.0414 278.2 0.038 30 -1.3 1.3 0.0129 278.2 0.043 30 -1.3 1.25 0.0119 1 0.1 30 1.2 1.2 0.0278 3.383 0.1 30 1.2 1.2 0.0440 21.34 0.1 30 1.65 1.65 0.0798 43 0.1 30 1.95 2.025 0.0969 100 0.1 30 2.35 2.575 0.0878 177 0.1 30 2.35 2.65 0.0285 278.2 0.1 30 -2.9 2.45 0.0712 324 0.1 30 -2.9 2.325 0.1158 21.34 0.03 30 0.65 0.625 0.0088 21.34 0.06 30 0.75 0.425 0.0356 21.34 0.1 30 1.65 1.65 0.0798 21.34 0.13 30 -2.7 2.4 0.0426 TABLE 4.11 Maximum Steady Stream Function for Oscillating Circular Body. The numerical and the experimental steady streamlines are shown in Fig-ures 4.13, 4.14 and 4.15 for values of R„ of 3.834, 21.34 and 278.2, respectively. The numerical steady streamlines shown in Figure 4.16 are for 0 = 0.1 and in Figure 4.17 are for Ru = 21.34. The numerical results are obtained by plotting 31 equally spaced contours between the maximum and the minimum value of the steady stream function in the plotted domain. For each case the value of the location of the maximum stream function and its value are presented in Table 4.11. Since the body in the present case is symmetric about the y-axis the value of minimum steady stream function is the negative of the maximum value and located at the corresponding symmetric point about the y-axis. Tatsuno has experimentally observed that when R^ is less than about 64 the whole flow field in each quadrant is occupied by only one circulatory streaming vortex (the NUMERICAL INVESTIGATIONS / 4.6 75 NUMERICAL INVESTIGATIONS / 4.6 76 0 = 0.0925 0 = 0.0925 F I G U R E 4.14 Steady Streamlines for Oscillating Circular Body for R = ^ = 15.5. 21.34, NUMERICAL INVESTIGATIONS / 4.6 77 NUMERICAL INVESTIGATIONS / 4.6 78 inner circulation). Also, at a critical value of the body amplitude parameter fi the inner circulation shrinks down sharply and this critical value depends on R„. In Figure 4.13 the experimental results are shown for = 3.834 and fi = 0.161 whereas the numerical results are shown for the same value of Rw and three values of fi, namely, 0.1, 0.13 and 0.161. The experimental result is for Ru less than 64 and clearly shows the inner circulation occupying the whole flow field in each quadrant. This indicates that the critical value of fi for inner circulation to shrink is greater than 0.161. The numerical results for fi = 0.1 and 0.13 show the inner circulation occupying the whole flow field in each quadrant, similarly to the experiment. The flow pattern of the inner circulation and the location of its core is very close to the experiment. But the numerical result for fi = 0.161 shows that the inner circulation has sharply shrunk. This indicates that the critical value of fi is between 0.13 and 0.161. Overall, the numerical results exhibit the same phenomena as observed experimentally. In Figure 4.14 the experimental and numerical are both shown for Ru = 21.34 and fi = 0.0925. For both cases the core of the inner circulation is closer to the cylinder than in the previous case. Nevertheless, the circulation expands over the whole flow field. The core of the inner circulation appears to be closer to the cylinder for the experimental case as compared to the numerical case. Overall, the agreement between the experimental and the numerical results is very good. For both cases, the inner circulation occupies the whole flow field in each quadrant. This indicates two things, namely, (1) the value of R^ for transition to the case where the secondary vortex system appears is greater than J2W = 21.34 and (2) at R^ = 21.34 the critical value of fi at which the inner circulation shrinks is greater than fi = 0.0925. In Figure 4.15 the experimental results are shown for Ru = 278.2 and fi = 0.046 whereas the numerical results are shown for the same value of Rw and two values of fi, namely, 0.038 and 0.043. The numerical and experimental results are for Ru greater than 64 and both clearly show the appearance of the secondary vortices close to the body. NUMERICAL INVESTIGATIONS / 4.6 79 There is very good agreement between the experimental and numerical results regarding the flow pattern and the location of both the vortices in each quadrant. However, the numerical results show the secondary vortex more tightly packed to the circular body than the experimental one. The numerical result for fi = 0.046 did not converge even when the solution for fi = 0.043 was used as an initial guess solution. This may represent a bifurcation point close to fi = 0.046 for transition to different flow patterns. This may also be indicated by the fact that the stream function for fi = 0.043 exhibits a sharper curvature close to the horizontal stagnation points as compared to the case of fi = 0.038. The photographs for the experimental results in Figures 4.13 and 4.14 were obtained from Reference[36] whereas the photograph in Figure 4.15 appears in Refer-ence^?] and the original photograph was obtained directly from Tatsuno. The numerical steady streamlines shown in Figure 4.16 are for fi of 0.1 and Ru of 1, 3.383, 21.34, 43, 100, 177, 278.2 and 324. The stream function results show that the core of the primary vortex system moves away from the circular body as Ru is increased from 1 to 324. This is also indicated in Table 4.11, which shows that the location of the maximum value of the steady stream function moves away from the body for the present case as Rw is increased. The numerical streamlines also show that the inner circulation occupies the whole flow field for R„ of 1, 3.383, 21.34, 43 and 100. The results for Ru of 177, 278.2 and 324 show the appearance of the secondary vortices close to the body. This indicates that the value of Ru for transition to the case where the secondary vortices appear is between R^ = 100 and R^ = 177. The fact that the numerical streamlines for Ru = 177 and fi = 0.03 (not shown here) also show the appearance of secondary vortices indicate that at Rw = 177 the secondary vortices appear irrespective of the value of fi. The numerical steady streamlines shown in Figure 4.17 are for Ru of 21.34 and fi of 0.03, 0.06, 0.1 and 0.13. The stream function results show that the inner circu-lation occupies the whole flow field for fi of 0.03, 0.06 and 0.1. For fi = 0.13 the inner circulation has shrunk close to the body. This indicates that the critical value of fi for inner circulation to shrink is between 0.1 and 0.13 for R„ = 21.34. NUMERICAL INVESTIGATIONS / 4.6 80 NUMERICAL INVESTIGATIONS / 4.6 81 NUMERICAL INVESTIGATIONS / 4.6 82 K P D0 b Location of Maximum Steady Velocity X y - «o 1 0.1 30 0.5 0 0.2859 3.834 0.1 7.7 0.0975 0.4904 0.0119 3.834 0.1 30 0.5 0 0.4942 3.834 0.13 7.7 0.0975 0.4904 0.0163 3.834 0.161 7.7 0.4904 0.0975 0.0582 3.834 0.161 30 0.3536 0.3536 0.2486 21.34 0.03 30 0.3536 0.3536 0.1010 21.34 0.06 30 0.3536 0.3536 0.7288 21.34 0.0925 15.5 0.1913 0.4619 0.4794 21.34 0.0925 30 0.5 0 1.0424 21.34 0.1 30 0.5 0 0.8141 21.34 0.13 30 0.3536 0.3536 0.6353 43 0.1 30 0.1913 0.4619 0.8819 100 0.1 30 0.1913 0.4619 0.8281 177 0.1 30 0.1913 0.4619 0.6835 278.2 0.038 30 0.3536 0.3536 0.2317 278.2 0.043 30 0.4619 0.1913 0.2587 278.2 0.1 30 0.3536 0.3536 0.5581 324 0.1 30 0.4619 0.1913 0.5678 TABLE 4.12 Maximum Steady Component of Velocity for Oscillating Circular Body. The steady flow components illustrated in these figures are expected to be small compared to the main periodic flow. The actual maximum of the resultant of the steady component of the velocities, its location in the flow domain and ^ ratio for which the numerical computation is performed are shown in Table 4.12. Since the body is sym-metric about the y-axis these maximum velocities occur at two locations symmetric about the y-axis. These maximum velocities are obtained only from a check of nodal values, and hence may not be the global maximum. In all cases the maximum velocity occurred at the body mean position. It should be noted that for one anomalous case, that is, when = 21.34, /? = 0.0925 and ^ = 30, the maximum steady velocity is greater than the characteristic body velocity u 0 . Even this case does not violate the body boundary condition because NUMERICAL INVESTIGATIONS / 4.6 83 of the superharmonics obtained in the computation of the velocities at the body mean position given by equations (4.2). In Table 4.12, the effect of the change of ^ is observed for three flow sit-uations, namely, (1) Ru = 3.834 and p = 0.1, (2) Ru = 3.834 and p = 0.161 and (3) R„ = 21.34 and p = 0.1. In each case there is a change in the location and magnitude of the maximum steady velocity with a change in The flow pattern for the steady streaming looks essentially the same for the flow situations (1) and (3). However some differences were observed for the flow situation (2). This may be due to a shift in the critical value of p for the inner circulation to shrink. Further investigations are necessary to examine the effect of The steady state solution obtained from equations (2.32) can be written as shown in equations (4.25). Some of the results from equations (4.25) are shown in Figures 4.18, 4.19 and 4.20 for = 3.834, 278.2 and 324, respectively. The results in Figures 4.18 and 4.20 are shown for the body amplitude parameter P = 0.1 and the one in Figure 4.19 for p = 0.043. These figures show velocity vectors over a uniform mesh in the vicinity of the body for t = 0, JT/4, ir/2 and 3TT/4, that is, equal time steps over one-half cycle of the body motion. The size of each vector indicates the velocity magnitude. Note that the scale in each illustration is different and the magnitude of the velocities is obtained by dividing the length of the arrows in centimetres, by the scale. It is clear from these figures that there is one main vortex in the upper half of the fluid flow and although it does move slightly from side to side it is located just above (and below by symmetry) the oscillating body. As the frequency Reynolds number Ry, increases, this vortex moves closer to the body and decreases in size, thereby increasing the velocity gradients in that region. Consequently, the boundary layer thickness is also decreased. A full nonlinear Newton-Raphson solver with the Jacobian [J] updated ev-ery alternate iteration is used for solving the nonlinear algebraic equations. Overall, the NUMERICAL INVESTIGATIONS / 4.6 84 \\\\\\\\\m : : : : : ; 111111 ; | ! | | | * %. \ % > ' M i ' P I ^-J-ji '• ::• = 5 = 22 = 33: : : : : : : 3 5 = 5 3 = 33k \\\[\\\\\\\\\\\\\\ mNmjil f * * > 5 5 i Ijij HP w ttj £JJ: :::::::: : * = 0 sca/e = 0.2292 < = f sco/c = 0.1771 t = f sca/e = 0.1775 < = *f Sca/e = 0.3123 F I G U R E 4.18 Velocity Vectors for Oscillating Circular Body. = 3.834, 0 = 0.1, ^ = 3 0 . NUMERICAL INVESTIGATIONS / 4.6 85 t = 0 scale = 0.2990 t = \ scale = 0.3085 t = § scale = 0.5672 t = scale = 0.2755 F I G U R E 4.19 Velocity Vectors for Oscillating Circular Body. R„ = 278.2, P = 0.043, ^ = 30. NUMERICAL INVESTIGATIONS / 4.6 86 t = 0 scale = 0.2760 FIGURE 4.20 Velocity Vectors for Oscillating Circular Body. R„ = 324, £ = 01. T = 3 0-NUMERICAL INVESTIGATIONS / 4.6 87 K Finite Element Results Exact Results (Stuart[38]) Added Mass, ma Added Damping, ca Added Mass, ma Added Damping, ca 1 7.2981 12.1030 6.6569 13.6569 3.383 4.2444 5.1392 4.0756 5.4403 21.34 2.2017 1.5967 2.2244 1.5993 43 1.8329 1.0587 1.8627 1.0487 100 1.5093 0.6730 1.5657 0.6457 . 177 1.3505 0.4869 1.4252 0.4704 278.2 1.2592 0.3711 1.3392 0.3679 324 1.2343 0.3377 1.3143 0.3390 T A B L E 4.13 Added Mass, ma and Added Damping, ca for Oscillating Circular Body for 0 = 0. nonlinear solver works very well. The converged results are obtained within five iterations for all cases presented here except one. For the case of i ? w = 177, 0 = 0.13 and ^ = 30 the converged solution is obtained after eleven iterations when the case of = 3.383 and 0 = 0.03 is used as an initial guess solution. In some cases the converged solution is not obtained at all. For example, for R^ = 278.2, 0 = 0.0463 and ^ = 30 the converged solu-tion is not obtained even when the case of R^ = 278.2 and 0 = 0.043 is used as an initial guess solution. This numerical instability may represent a bifurcation point to different flow patterns. 4.6.3 A d d e d Mass and A d d e d Damping . The results for the added mass, ma and added damping, ca are obtained using the finite element grid in Figure 4.12(c) with ^ = 30. First, the results are obtained for the linear case, that is, when the body amplitude parameter 0 = 0. For this case, Stuart[38] has presented analytical results for a circular body. Therefore, a direct comparison is possible for this case and the results are tabulated in Table 4.13 and presented graphically in Figures 4.21 and 4.22. The agreement between the exact analytical results and the finite element results is very good. The finite element results are within 6% of the exact results for all cases except one. For the case of Ru = 1, the finite element results for the added mass and added damping differ from the NUMERICAL INVESTIGATIONS / 4.6 NUMERICAL INVESTIGATIONS / 4.6 NUMERICAL INVESTIGATIONS / 4.6 90 K Added Mass, ma Added Damping, ca 500 1.1768 0.2564 5000 1.0631 0.0393 10000 1.0607 0.0198 T A B L E 4.14 Added Mass, ma and Added Damping, ca for Oscillating Circular Body for i2w large and 0 = 0. Added Mass, m„ 0 = 0 0 = 0.03 13 = 0.1 l 7.2981 7.2373 7.0753 3.383 4.2444 4.2003 4.0845 21.34 2.2017 2.1731 2.0903 43 1.8329 1.8091 1.7383 100 1.5093 1.4933 1.4634 177 1.3505 1.3417 1.3445 278.2 1.2592 1.2563 1.2759 324 1.2343 1.2332 1.2522 T A B L E 4.15 Added Mass, ma for Oscillating Circular Body. exact results by 9.6% and 11.4%, respectively. This may be due to the fact that the exact analytical results exhibit a singularity at = 0. That is, the value of the added mass and added damping go to infinity as Ru approaches zero. The comparison is also made with the limiting case of inviscid flow. The added mass and added damping values for a bluff body oscillating at a high frequency with a small amplitude is not affected much by viscosity. Hence, the values of m a and c„ for large values of R^ and 0 = 0 should approach the values predicted by inviscid flow. The results obtained for this case are presented in Table 4.14. For R„ = 10000, the value of the added mass obtained using the finite element method is 1.0607 as compared to 1 for the inviscid flow [33]. Therefore, the value of ma obtained using the finite element method differs only by about 6% from the limiting case of inviscid flow. The value of the added damping approaches zero as it should from the inviscid flow considerations. For the amplitude parameter 0 = 0 and for the finite values of 0 of 0.03 NUMERICAL INM5STIGATI0NS / 4.6 91 Added Damping, ca 0 = 0 0 = 0.03 0 = 0.1 1 12.1030 12.1369 12.2025 3.383 5.1392 5.1713 5.2158 21.34 1.5967 1.6290 1.6232 43 1.0587 1.0919 1.0601 100 0.6730 0.7027 0.6414 177 0.4869 0.5100 0.4356 278.2 0.3711 0.3865 0.3087 324 0.3377 0.3498 0.2825 T A B L E 4.16 Added Damping, ca for Oscillating Circular Body. and 0.1, the results for the added mass and added damping are tabulated in Tables 4.15 and 4.16 and presented graphically in Figures 4.23 and 4.24, respectively. For the case of finite values of the body amplitude parameter 0, no analytical or experimental results are available for direct comparison. The result for 0 = 0.03 shows that the added mass is lower and the added damping is higher than for the case of 0 = 0. The same is true for the more complex flow situation of j3 = 0.1 except at high values of Ru where the situation is reversed, that is, the added mass terms are higher and the added damping terms lower than the case of 3 = 0. Overall, the results for the added mass and added damping for the linear case of 0 = 0 compare very well with the exact analytical results of Stuart and approach the correct limiting values of the inviscid case. The results for small but finite values of the body amplitude parameter 0 show very little departure from the linear case and show the correct trends from the viscous flow considerations. Al l this indicates that the procedure for computing the added mass and added damping is correct and the finite element repre-sentation of the fluid domain and the Taylor series expansion for the nonlinear boundary conditions are adequate to capture the basic characteristics of the physical phenomena. 4.6.4 Stability Analysis. The results for the stability analysis are obtained by using the finite element NUMERICAL INVESTIGATIONS / 4.6 NUMERICAL INVESTIGATIONS / 4.6 93 NUMERICAL INVESTIGATIONS / 4.6 94 Eigenvalue No. \M \M 1 I I 2 l I 3 l I 122 l I 123 0.99597 0.99597 124 0.99594 0.99594 125 0.99372 0.99372 126 0.99147 0.99147 127 0.98686 0.98686 128 0.98472 0.98472 129 0.98353 0.98353 130 0.97974 0.97941 131 0.97941 0.97941 132 0.97645 0.97645 133 0.97352 0.97352 134 0.97244 0.97244 135 0.96951 0.96951 T A B L E 4.17 Symmetric Mode Stability Analysis for Oscillating Circular Body. Ru = 324, p = 0.1. grid in Figure 4.12(c) with ^ = 30. This grid has 61 pressure degrees of freedom and 256 velocity degrees of freedom for both the symmetric and antisymmetric modes of stability analysis. The stability analysis is performed for the oscillating circular body for R^ = 324 and P = 0.1. The results for the top 135 eigenvalues are presented in descending order of their magnitudes in Tables 4.17 and 4.18 for the symmetric and antisymmetic modes, respectively. The flow situation considered here is very close to the flow situations with the same R^ of 324 and P of 0.06 and 0.13, where the numerical results did not converge. The results show that there are 122 eigenvalues (twice the number of pressure degrees of freedom) with the magnitude equal to 1 and the rest are all less than 1. This, according to the criteria given by equation (3.15) indicates that the steady state solution NUMERICAL INVESTIGATIONS / 4.6 95 Eigenvalue No. |A| |A| 1 1 1 2 1 1 3 1 1 122 1 1 123 0.99846 0.99846 124 0.99504 0.99504 125 0.99490 0.99490 126 0.99123 0.99123 127 0.99006 0.99006 128 0.98816 0.98816 129 0.98676 0.98676 130 0.98421 0.98421 131 0.98370 0.98371 132 0.98168 0.98168 133 0.98105 0.98105 134 0.98012 0.98012 135 0.97781 0.97781 T A B L E 4.18 Antisymmetric Mode Stability Analysis for Oscillating Circular Body. R_, = 324, 0 = 0.1. is stable for both symmetric and antisymmetric perturbations. The results also show that the eigenvalues have essentially converged at At = JQ and further reduction of the time step interval to At = ^ brings about no change in them. Many of the eigenvalues in the present case have magnitudes close to 1 which qualitatively indicates a weakly stable situation. This may be due to the fact that the flow under consideration is very close to the flow for which a numerically converged steady state solution was not obtained, and therefore, may represent a bifurcation point for transition to different flow patterns. I t should be noted that in the present case also, the number of eigenvalues with magnitude 1 are twice the number of pressure degrees of freedom. NUMERICAL INVESTIGATIONS / 4.7 96 ID J O U K O W S K I P R O F I L E 4.7.1 General Remarks. Tatsuno[37] has presented excellent experimental results for the oscillating symmetric Joukowski profile which are used for the verification of the present method. Tatsuno again presents the results for the basic nonlinear phenomenon of steady streaming. This case is considered to investigate the flow pattern around an asymmetrical body. The numerical results are obtained for a symmetric Joukowski profile with a characteristic length b of unity. The body is performing oscillations parallel to the x-axis. Using symmetry, only one half of the domain is modelled. The finite element grid shown in Figure 4.25 is developed for this purpose. To compare with Tatsuno's experimental setup, the ratio of thickness to chord length for the Joukowski profile is taken to be 0.6 and the grid has ratio ^ = 18.3, where DQ is the diameter of the outer fixed cylinder. NUMERICAL INVESTIGATIONS / 4.7 97 For the flow situation, the boundary conditions are specified by the velocities u = v = 0 at the outer boundary and the velocities specified according to equations (2.27) at the body boundary. The symmetry condition is obtained by setting the velocity v = 0 on the symmetry line. The pressure needs to be specified only at one location in the domain to obtain the pressure datum. In the present case the pressure is specified to be zero at location (1) shown in Figure 4.25. For the stability analysis, the boundary conditions are given by the velocity perturbations eu = sv = 0 at the outer boundary and on the body boundary since the velocities are specified at these locations. Since only one half of the domain is modelled, two cases are to be investigated namely, the symmetric and antisymmetric perturbations. The symmetric case is obtained by setting the velocity perturbation sv = 0 on the symmetry line and the antisymmetric case is obtained by setting the velocity perturbation e„ = 0 on the symmetry line. The pressure perturbation ep = 0 at the same location (1) shown in Figure 4.25. The finite element grid shown in Figure 4.25 has 48 elements and 173 nodes. There are 352 net degrees of freedom, 290 of which are for the velocities and 62 for the pressure. Out of the 352 net degrees of freedom, 34 velocity degrees of freedon lie on the body boundary and are specified from the moving body boundary condition. For the flow situation, there are three types of coefficients A , B, C, and therefore, this results in 954 variables in total. For the stability analysis, the total number of variables is 318. The flow results are presented in two ways, namely, the instantaneous veloc-ity vectors in the flow field and the streamline contours for the steady part of the velocity. The streamline contours obtained from the steady component A of the solution will here-after be referred to as steady streamlines. The numerical steady streamlines are compared with the corresponding experimental streamlines. The added mass, added damping and added force quantities are presented for a range of and 3 values including the linear case of 0 = 0. The added force is zero NUMERICAL INVESTIGATIONS / 4.7 98 K P Location of Minimum Stream Function Minimum Steady Stream Function Location of Maximum Stream Function Maximum Steady Stream Function X y X y 55.1 0.01 -0.7 0.35 -0.0015 0.05 0.425 0.0016 55.1 0.064 -1.15 1.05 -0.0341 -0.05 0.325 0.0138 55.1 0.1 0.7 2.275 -0.0469 -0.25 0.375 0.0047 206 0.0091 -0.6 0.275 -0.0013 -0.2 1.575 0.0018 206 0.0131 -0.6 0.275 -0.0182 -0.05 1.575 0.0029 206 0.0188 -0.06 0.275 -0.002686 0.05 1.575 0.005125 206 0.0251 -0.6 2.275 -0.0037 0.25 1.575 0.0096 206 0.0325 -0.6 0.275 -0.0049 0.4 1.575 0.0234 206 0.0622 1.05 1.125 -0.0337 -0.05 0.3 0.0112 206 0.1 1.5 2.125 -0.1205 -0.25 0.35 0.0049 TABLE 4.19 Maximum and Minimum Steady Stream Function for Oscillating Joukowski Profile. for the linear case of P = 0 but is not necessarily zero for the finite values of (3 because of the asymmetry of the Joukowski profile. The stability analysis for the symmetric and antisymmetric modes is again performed only for one flow situation. 4.7.2 Flow Results. The numerical and experimental steady streamlines are shown in Figures 4.26 and 4.27 for values of R^ of 55.1 and 206, respectively. The numerical results are obtained by plotting 31 equally spaced contours between the maximum and minimum value of the steady stream function in the plotted domain. For each case the location and value of both the maximum and minimum values of the stream function are presented in Table 4.19. An example of the flow at = 55.1, a relatively small value, is given in Figure 4.26. In this case the experimental result is for the body amplitude parameter p of 0.06 whereas the numerical results are for three values of (3, namely, (1) 0.1, (2) 0.06 and (3) 0.01. The flow field around the body is almost totally occupied by the inner circulatory streaming for the experimental case. The numerical case that compares very well with the FIGURE 4.26 Steady Streamlines for Oscillating Joukowski Profile for Rw Q = 18.3. = 55.1, NUMERICAL INVESTIGATIONS / 4.7 100 Ik) FIGURE 4.27 Steady Streamlines for Oscillating Joukowski Profile for J2 W = 206, - 18.3 (a) fi = 0.0091, (b) fi = 0.0131, (c) fi = 0.0188, (d) fi = 0.0251, (e) fi = 0.0325, (f) fi = 0.0622, (g) fi = 0.1, (h) fi = 0.0091, (i) fi = 0.0131, 0) fi = 0.0251, fk) fi = 0.0457. NUMERICAL INVESTIGATIONS / 4.7 101 experiment is for P = 0.01. For this case, too, the flow field around the body is almost occupied by the inner circulatory streaming. The pair of vortices above the body (and below by symmetry) are located almost symmetrically about the y-axis with the vortex on the cuspidal edge of the profile occupying slightly a larger portion of the domain. As P is increased , the asymmetry of the location of the vortices becomes more pronounced. At 0 = 0.06, same as the experiment, the numerical steady streamlines show that the vortex on the obtuse edge occupies a larger portion of the plotted domain and the vortex on the cuspidal edge has shrunk. At P = 0.1, the vortex which was on the obtuse edge for P = 0.06, has now shifted its location towards the cuspidal edge and completely occupies the second vortex, which itself has shrunk even more. The numerical result for P = 0.01 compares very well with the experiment for P = 0.06. Overall, the numerical results show the characteristics of the experimental flow. The example in Figure 4.27 shows the steady streamlines for the numerical case at seven different values of P and for the experimental case at four different values of P, all for a moderate value of of 206. The numerical results are for P of 0.0091, 0.0131, 0.0188, 0.0251, 0.0325, 0.0622, and 0.1 and the experimental results are for p of 0.0091, 0.0131, 0.0251 and 0.0457. The experimental results show a more complicated flow pattern for small values of p. At P = 0.0091, it is observed that there are three pairs of vortices very close to the body. Two of them are located adjacent to the body while the third one is located further away in the flow field on the cuspidal edge. As the value of p is increased, it is observed that the central vortex shrinks and the vortex on the cuspidal edge becomes more dominant, and at P = 0.0457 the central vortex has almost disappeared. In all the experimental results shown, the vortex on the obtuse edge remains essentially at the same location. The numerical result at P = 0.0091 shows three main vortices, two of which are located adjacent to the body while the third one is located further away in the flow field. As the value of p is increased the vortex adjacent to the body and close to the cuspidal edge shrinks and the vortex further away in the flow field becomes more dominant and NUMERICAL INVESTIGATIONS / 4.7 102 moves towards the cuspidal edge. At /? = 0.0251 and /? = 0.0325 the vortex adjacent to the body and close to the cuspidal edge has disappeared and the vortex further away in the flow field has moved to the cuspidal edge. For all values of /? from 0.0091 to 0.0325, the vortex adjacent to the body and close to the obtuse edge remains at essentially the same location. The numerical results for = 206 and ft = 0.0457 did not converge even when the converged results for = 206 and two separate values of /? of 0.0325 and 0.0622 were used as initial guess solutions. This may indicate a change in flow pattern which is indeed seen at /? = 0.0622. In this case again three vortices appear with two vortices located close to the obtuse edge and tightly packed adjacent to the body while the third one spreads more into the flow field and is located at the cuspidal edge. As (3 is increased to 0.1, the vortices located close to the body shrink and the vortex on the cuspidal edge becomes more dominant and spreads over the whole plotted domain. Unfortunately, the experimental results at higher values of /? are not available to observe the existence of a critical value of (3 for transition to a different flow pattern. Overall, the numerical results exhibit the trends observed experimentally. The photographs for the experimental results in Figures 4.25 and 4.26 appear in Reference[37] and original photographs were obtained directly from Tatsuno. The steady flow components illustrated in these figures are again expected to be small compared to the main periodic flow. The actual maxima of the resultant velocities and their locations in the flow domain are shown in Table 4.20. These maximum velocities were obtained only from a check of nodal values, so they may not be the global maxima. In all cases, except for R„ = 55.1 and j3 = 0.1, the maximum velocity occurred at the body mean position. For the aforementioned case, the maximum steady velocity occurs at a position very close to the body mean position. The steady state solution obtained from equations (2.32) can be written as shown in equations (4.25). Some of the results from equations (4.25) are shown in NUMERICAL INVESTIGATIONS / 4.7 103 K P Location of Maximum Steady Velocity umax X y «0 55.1 0.01 -0.05105 0.2313 0.0407 55.1 0.064 -0.0567 0.2466 0.6910 55.1 0.1 -0.7 0 0.5134 206 0.0091 -0.5105 0.2313 0.0582 206 0.0131 -0.5105 0.2313 0.0839 206 0.0188 -0.5105 0.2313 0.1206 206 0.0251 -0.5105 0.2313 0.1610 206 0.0325 0.1134 0.1431 0.3073 206 0.0622 0.1134 0.1431 0.8377 206 0.1 -0.5105 0.2313 0.5473 TABLE 4.20 Maximum Steady Component of Velocity for Oscillating Joukowski Profile. Figures 4.28,4.29 and 4.30. The results for R„ = 55.1 are shown in Figures 4.28 and 4.29 for 0 of 0.01 and 0.1, respectively, and the results for Ru = 206 and 0 = 0.1 are shown in Figure 4.30. These figures show velocity vectors over a uniform mesh in the vicinity of the body for t = 0, JT/4, TT/2, 3TT/4 T T, 5TT/4, 3TT/2 and 7JT/4, that is, equal time steps over one full cycle of body motion. Unlike the case of the square and the circular body where it is necessary to show the results only for one-half cycle of body motion because of the symmetry, the results for the Joukowski profile have to be shown for one full cycle. The size of each vector indicates the velocity magnitude. Note that the scale in each illustration is different. The magnitude of the velocities is obtained by multiplying the length of the arrows in centimetres by 1.5 and dividing by the scale. It is clear from these figures that there is one main vortex in the upper half of the fluid flow and although it does move slightly from side to side it is located just above (and below by symmetry) the oscillating body. For Ru = 55.1 and 0 of 0.01 and 0.1 the location of the main vortex is essentially the same. As i2 w is increased to 206, this vortex system moves closer to the body and decreases in size, thereby increasing the velocity gradients in that region. Consequently, the boundary layer thickness is also decreased. N U M E R I C A L INVESTIGATIONS / 4.7 104 1 M 1 t = 0 scale = 0.3131 t = \ scale = 0.4289 t = | scale = 0.7698 < = ^ scafe = 0.4493 * = ir scale = 0.3151 t = ^ sca/e = 0.4298 < = ^ scale = 0.7718 < = If scale = 0.4454 FIGURE 4.28 Velocity Vectors for Oscillating Joukowski Profile. Ru 0.01, ^ = I 8 - 3 -= 55.1, B = NUMERICAL INVESTIGATIONS / 4.7 105 t = 0 scale = 0.2100 t = \ scale = 0.1785 t = ir scale = 0.3065 i b i i ' * = ^ scale = 0.3368 . . . : : | •'I'.. < = § sca/e = 0.2402 < = ^ scale = 0.3258 * = % sco/c = 0.7597 < = 2f. scale = 0.2841 FIGURE 4.29 Velocity Vectors for Oscillating Joukowski Profile, = 55.1, fi 0.1, ^ = 18.3. NUMERICAL INVESTIGATIONS / 4.7 106 ii t = 0 scale = 0.2753 33-t = f sca/c = 0.2851 t = | sca/e = 0.2458 < = ^ scale = 0.3161 «t = ir scale = 0.3194 t = ^ sco/e = 0.3079 t = ^ scale = 0.5794 < = ?f scale = 0.2222 FIGURE 4.30 Velocity Vectors for Oscillating Joukowski Profile. = 206, 0 = 0.1, ^ = 18.3. NUMERICAL INVESTIGATIONS / 4.7 107 K Added Mass, m a Added Damping, ca 500 0.8993 0.3971 5000 0.6951 0.0793 10000 0.6830 0.0426 1 x 106 0.6775 0.0033 T A B L E 4.21 Added Mass, ma and Added Damping, ca for Oscillating Joukowski Profile for large and 0 = 0. 4.7.3 A d d e d Mass, Added Damping and Added Force. The results for the added mass (m a ) , added damping (c a) and added force (/„) are obtained using the finite element grid in Figure 4.25 with ^ = 18.3. Unlike the previous cases of the oscillating square and circular bodies, the oscillating symmetric Joukowski profile gives the added force due to its asymmetry. Since the added force is caused by the asymmetry in the steady component of the velocity, the linear case of 0 = 0 does not give the added force for any body profile. First, the results are obtained for the linear case, that is, when the body amplitude parameter 0 = 0. For this simple case too, no analytical results are available to make a direct comparison. The values ma and ca for large values of Ru and 0 = 0 are presented in Table 4.21. The limiting value obtained for ma is 0.6775 and ca is 0.0033. These limiting values should approach the values predicted by the inviscid flow. The results for symmetric Joukowski profile from the inviscid flow considerations are not available but the added mass can be approximately taken to be equal to pit ( | ) , where t is the thickness of the Joukowski profile, using the results for ellipse or plate from Reference[33]. This gives the non-dimensional value of the added mass, m B to be 0.72 which is within 6% of the limiting value 0.6775 obtained numerically. The value of ca approaches zero as it should from inviscid flow considerations. The values of the added mass and added damping for moderate values of Ru and for the linear case of 0 = 0 are presented in Table 4.22 and graphically in Figures 4.31 and 4.32. NUMERICAL INITSTIGATIONS / 4.7 108 Ru Added Mass, ma Added Damping, ca 1 9.9114 20.0877 3.383 5.5558 8.3242 21.34 2.3333 2.5329 43 1.7894 1.5888 100 1.3987 0.9696 177 1.1789 0.7252 278.2 1.0331 0.5656 324 0.9925 0.5172 T A B L E 4.22 Added Mass, ma and Added Damping, ca for Oscillating Joukowski Profile for p = 0. The results for the added mass, added damping and added force presented in Table 4.23 are for the finite values of body amplitude parameter 0. For this case too there are no published results to compare with. The results in Table 4.23 show a general trend of decreasing values of ma and ca with increasing Rw same as for the linear case. It is also observed that ma and ca show only a small change with a small change in 0 and the values are close to the linear case for the values of P considered here. Generally, the added force increases with increasing except for one case. For the case of Rw = 206 and P = 0.0325, the added force becomes negative. I t is also observed from Table 4.23 that for a constant P (for example, p = 0.06 and P = 0.1) the value of fa decreases with the increasing i2 w . The results presented in Table 4.23 are for the flow situations considered in Subsection 4.7.2. Overall, the results for the linear case of p = 0 show the trends observed for the square and circular body and approach a limiting value for ma and ca when R^ becomes very large. The values for the finite values of p exhibit the correct trends. 4.7.4 Stabil i ty Analysis. The results for the stability analysis are obtained by using the finite element grid in Figure 4.25 with = 18.3. This grid has 62 pressure degrees of freedom and 256 velocity degrees of freedom for both the symmetric and antisymmetric modes of stability analysis. NUMERICAL INVESTIGATIONS / 4.7 109 FIGURE 4.31 Added Mass, ma for Oscillating Joukowski Profile for 0 = 0. NUMERICAL INVESTIGATIONS / 4.7 110 NUMERICAL INVESTIGATIONS / 4.7 Ul K fi Added Mass, ma Added Damping, ca Added Force, fa 55.1 0.01 1.6546 1.3613 0.6264 x 10 - 5 55.1 0.064 1.5565 1.4289 0.9697 x 10 - 3 55.1 0.1 1.6851 1.3470 0.3173 x 10~ 2 206 0.0091 1.1242 0.6715 0.2953 x 1 0 - 5 206 0.0131 1.1221 0.6742 0.5914 x 10 - 5 206 0.0188 1.1173 0.6806 0.1117 x 10" 4 206 0.0251 1.1082 0.6942 0.1609 x 10~ 4 206 0.0325 1.0826 0.7445 -0.3959 x 10~ 5 206 0.0622 1.1409 0.6040 0.2276 x 10~ 3 206 0.1 1.1868 0.6102 0.1047 x 10" 2 T A B L E 4.23 Added Mass, ma, Added Damping, c 0 and Added Force, fa for Oscil-lating Joukowski Profile. Eigenvalue No. A * — — 1*1 At - 2s. " L — 40 |A| 1 1 1 2 1 1 3 1 1 124 1 1 125 0.98638 0.98638 126 0.98478 0.98478 127 0.97876 0.97876 128 0.97069 0.97069 129 0.96840 0.96840 130 0.96319 0.96319 131 0.96020 0.96020 132 0.95102 0.95102 133 0.94683 0.94683 134 0.93428 0.93428 135 0.93314 0.93314 T A B L E 4.24 Symmetric Mode Stability Analysis for Oscillating Joukowski Profile. Rw = 206, fi = 0.0251. The stability analysis is performed for the oscillating Joukowski profile for = 206 and fi = 0.0251. The results for the top 135 eigenvalues in descending order NUMERICAL INVESTIGATIONS / 4.7 112 Eigenvalue No. |A| |A| 1 1 1 2 1 1 3 1 1 124 1 1 125 0.99460 0.99460 126 0.98569 0.98569 127 0.98478 0.98478 128 0.98154 0.98154 129 0.97874 0.97874 130 0.97059 0.97059 131 0.96721 0.96721 132 0.96319 0.96319 133 0.96045 0.96045 134 0.95995 0.95995 135 0.95067 0.95067 T A B L E 4.24 Antisymmetric Mode Stability Analysis for Oscillating Joukowski Pro-file. Rw = 206, B - 0.0251. of magnitude are presented in Tables 4.24 and 4.25 for the symmetric and antisymmetic perturbation modes, respectively. The flow situation considered here is very close to the flow situation with the same Rw of 206 and B of 0.0457, where the numerically converged results were not obtained. The results show that there are 124 eigenvalues (twice the number of pressure degrees of freedom) with magnitude equal to 1 and the rest are all less than 1. This, according to the criteria given by equation (3.15) indicates that the steady state solution is stable for both symmetric and antisymmetric modes. The results also show that the eigenvalues have essentially converged at A* = and further reduction of the time step interval to At = | | brings about no change in them. Many of the eigenvalues in the present case have magnitudes close to 1 which qualitatively indicates a weakly stable situation. This may be due to the fact that the flow under consideration is very close to the flow for which the numerical results did not converge, and therefore, may represent a NUMERICAL INVESTIGATIONS / 4.7 113 bifurcation point for transition to different flow patterns. The number of eigenvalues with magnitude 1 is again twice the number of pressure degrees of freedom. CHAPTER 5 CONCLUSIONS 5.1 C O N C L U D I N G R E M A R K S The present method of representing the interaction between a solid body and viscous flow seems to work very well. In particular, the Taylor series expansion of the oscillating body boundary conditions with respect to the fixed finite element grid is apparently adequate to capture the basic nonlinear phenomenon of secondary flow. The full nonlinear Newton-Raphson procedure for solving the nonlinear al-gebraic equations is successful for the range of frequency Reynolds number considered in the present study. The overall agreement between the numerical and the experimental results for the basic nonlinear phenomenon of steady streaming is very good. The results for the instantaneous velocity vectors in the fluid domain show one main vortex in the upper half ( and lower half by symmetry ) of the transient flow. As the frequency Reynolds number Ru increases, this vortex moves closer to the body and decreases in size, thereby increasing the velocity gradients in that region. This indicates that the boundary layer thickness decreases with the increasing R^, as it should. The results for the added mass and added damping for the linear case of 0 = 0 compare very well with the other finite element results and exact analytical results. 114 CONCLUSIONS / 5.2 115 The values of the added mass and added damping for the linear case of /? = 0 approach the correct limiting values for the inviscid flow when Ru is large. The values of the added mass, added damping and added force are also obtained for finite values of the body amplitude parameter /?. But no results are available for direct comparision for this case. Overall, the implicit numerical integration scheme using the trapezoidal rule to obtain the transition matrix for stability analysis seems to work very well. The con-vergence to the transition matrix and thereby to the eigenvalues is obtained by increasing the number of time steps (that is, decreasing the time step interval). The scheme is un-conditionally stable with respect to the size of the time step. The computer time required for obtaining the converged results is very high because the number of matrix inversions required are equal to the number of time steps. Therefore, any increase in the number of time steps would correspond to a substantial increase in the computer time required. 5.2 S U G G E S T I O N S F O R F U R T H E R D E V E L O P M E N T Some specific recommendations for further studies based on the work in this thesis are: (1) Extend the analysis presented here to the case of a body oscillating in fluid with currents. This will provide a method for analysing the interaction between waves, currents and structures. (2) Extend the analysis presented here to the case of a more realistic multi-degree of freedom model for the structure. (3) Incorporate turbulence modelling for the fluid to obtain results for large values of Reynolds number suitable for engineering applications. 116 References [1] Brebbia, C. A. and Walker, S. (1979) Dynamic Analysis of Offshore Struc-tures, Newnes-Butterworths, London. [2] Zienkiewicz, O. C , Bettess, P., and Kelly, D. W. (1978) The Finite Element Method for Determining Fluid Loadings on Rigid Structures. Two- and Three-Dimensional Formulations. In Numerical Methods in Offshore Engi-neering, eds. O. C. Zienkiewicz, P. Lewis, and K. G. Stagg, Wiley, Chichester, England, pp. 141-183. [3] Armen, Harry (ed.), Stiansen S. (ed.) (1980) Computational Methods for Offshore Structures, ASME, AMDV 3 7 . [4] Chan, S. T. K. and Brashears, M. R. (1976) Analysis of Transonic Flow Over Lift ing and Oscillating Airfoils, 2nd Symp. on Finite Element Meths. in Flow Probs., Santa Margherita Ligure, Italy. [5] Shen, S. F. (1978) Transonic Aerodynamic Computation with Finite Ele-ment Method, in Finite Elements in Fluids, voi.3, eds. Gallagher, R. H., Zienkiewicz, O. C , Oden, J. T., Cecchi, M. M., Taylor, C , Wiley, London, pp 183-204. [6] Olson, M. D.and Irani, M. (1983) Finite Element Analysis of Viscous Flow-Solid Body Interaction, 2nd Int. Conf. on Num. Meths. in Laminar and Turbulent Flow, Seattle, U.S.A. [7] Olson, M. D.and Pattani, P. G. (1984) Nonlinear Analysis of Rigid Body-Viscous Flow Interaction, 5th Symp. on Finite Element Meths. in Flow Probs., Austin, Texas, also in Finite Elements in Fluids, vol.6, eds. Gal-lagher, R. H., Carey, G. F., Oden J. T., Zienkiewicz, O. C , Wiley, (1985). [8] Liu, W. K. (1980) Development of Finite Element Procedures for Fluid-Structure Interaction, Ph.D. Thesis, California Inst, of Tech., Pasadena. [9] Isaacson, M. de St. Q. (1974) The Forces on Circular Cylinders in Waves, Ph.D. Thesis, University of Cambridge, England. [10] Stokes, G. G. (1851) On the Effect of the Internal Friction of Fluids on the Motion of the Pendulums, Trans. Camb. Phil. Soc. 9 , Pt. II, pp 8-106. See also Math, and Phys. Pap. 3 , pp 1-141, Cambridge Univ. Press. 117 [11] Schlichting, H. (1968) Boundary Layer Theory, McGraw-Hill, 6th Ed. [12] Olson, M. D. and Tuann, S. Y. (1978) Primitive Variables versus Stream Function Finite Element Solutions of the Navier-Stokes-Equations, in Finite Element in Fluids, vol. 3, eds. Gallagher, R. H., Zienkiewicz, O. C., Cecchi, M. M., and Taylor, C., Wiley, London. [13] Tuann, S. Y.and Olson, M. D. (1978) Numerical Studies of the Flow around a Circular Cylinder by a Finite Element Method, Computers and Fluids, Vol.6, pp 219-240. [14] Roache, P. J. (1972) Computational Fluid Dynamics, Hermosa Publishers, Albuquerque, New Mexico. [15] Davis, R. W. and Moore, E. F. (1981) The Numerical Solutions of Flow Around Squares, Proc. 2nd Int. Conf. on Num. Meth. in Lam. and Turb. Flow, Venice, Italy, pp 279-290. [16] Hughes, T. J. R., Liu, W. K., and Zimmermann, T. K. (1978) Langragian-Eulerian Finite Element Formulation for Incompressible Viscous Flows, U.S.-Japan Conf. on Interdisciplinary Finite Element Analysis, Cornell Univer-sity, Ithaca N.Y. [17] Nayfeh, A. H. and Mook, D. T. (1979) JVoniinear Oscillations, Wiley, New-York. [18] Lin, C. C. (1955) The Theory of Hydrodynamic Stability, Cambridge Uni-versity Press, London. [19] Joseph, D. D. (1953) Stability of Fluid Motions I, Springer-Verlag, New-York. [20] Thomas, L. H. (1953) The Stability of Plane Poiseulle Flow, Phys. Rev., Vol.91, pp 780-783. [21] Grosch, C. E. and Salwen, H. (1968) The Stability of Steady and Time-Dependant Plane Poiseulle Flow, J. Fluid Mech., Vol.34, Part 1, pp 177-205. [22] Orszag, S. A. (1971) Accurate Solutions of Orr-Sommerfeld Stability Equa-tion, J. Fluid Mech., Vol.50, Part 4 , pp 689-703. [23] Saraph, V. R., Vasudev Rao, B., and Panikar, J. T. (1979) Stability of Parallel Flows by Finite Element Method, Int. J. JVum. Meths. Eng., Vol.14, pp 1257-1270. [24] Grosch, C. E. and Salwen, H. (1978) The Continous Spectrum of the Orr-Sommerfeld Equation Part 1: The Spectrum and Eigenfunction, J. Fluid Mech., Vol. 87, Part 1, pp 33-54. References 118 [25] Van Stijn, T H . L. and Van de Vooren, A. I. (1980) An Accurate Method of Solving the Orr-Sommerfeld Equation, J. Fluid Mech, Vol. 114, Part 1. [26] Fasel, H. (1976) Investigation of the Stability of Boundary Layers by a Finite Difference Model of the Navier-Stokes Equations, J. Fluid Mech, Vol. 78, Part 2, pp 355-383. [27] Olson, M. D. and Baldwin, J. F. (1981) Eigenvalue Perturbation of Navier-Stokes Equations by Finite Elements, Proc. 2nd Int. Conf. Num. Meth. Lam. Turb. Flows, Venice, Italy [28] Batchelor, G. K. (1974) An Introduction to Fluid Dynamics, Cambridge University Press, London, pp 139-140. [29] Finlayson, B. A. (1972) The Method of Weighted Residuab and Variational Principles, Academic Press, New York. [30] Friedmann, P., Hammond, C. E. and Woo, T. H. (1977) Efficient Numerical Treatment of Periodic Systems with Application to Stability Problems, Int. J. Num. Meths. Eng., Vol. 11 , pp 1117-1136. [31] Friedmann, P. and Silverthorn, L. J. (1974) Aeroelastic Stability of Periodic Systems with Application to Rotor Blade Flutter, AIAA J. 12, pp 1559-1565. [32] Hughes, T. J. R. (1983) Analysis of Transient Algorithms with Particular Reference to Stability Behavior, in Computational Methods for Transient Analysis, Vol. 1, eds. Belytschko, T. and Hughes, T. J. R., North-Holland. [33] Sarpkaya, T. and Isaacson, M. (1981) Mechanics of Wave Forces on Offshore Strctures, Van Nostrand Reinhold Company, New York. [34] Tatsuno, M. (1974) Circulatory Streaming in the Vicinity of an Oscillating Square Cylinder, J. Physical Soc. Japan, Vol. 36, No. 4, pp 1185-1191. [35] Irani, M. B. (1982) Finite Element Analysis of Viscous Flow and Rigid Body Interaction, Master's Thesis, Department of Civil Engineering, Univ. of British Columbia, Canada. [36] Tatsuno, M. (1973) Circulatory Streaming Around an Oscillating Circular Cylinder at Low Renolds Numbers, J. Physical Soc. Japan, Vol. 35, JVo. 3, pp 915-920. [37] Tatsuno, M. (1980) Secondary Flow Around Oscillating Asymmetric Cylin-ders, Reports of Research Institute of Applied Mechanics, Kyushu University, Vol. XXVm, No. 89, pp 35-47. [38] Stuart, J . T. (1963) Periodic Boundary Layers, Fluid Motion Memoirs-Laminar Boundary Layers, ed. Rosenhead, L., Oxford University Press, pp 390. APPENDIX 1 MATRICES OF SECTION 2.3 The net global u and v degrees of freedom for finite element discretization of viscous fluid is each numbering n including r degrees of freedom at the interface between the finite element grid and the mean position of the moving body. There are m net global p degrees of freedom. Also the u and v degrees of freedom at the edge of an element which is at the interface between the viscous fluid and the mean position of the body are each numbering q. Therefore, we can write the matrices in section 2.3 as given below: P = -i = n - r + l 0 0 0 0 E M^ch- 0 l=n-r+l 0 0 Q = P + y [6x. +6* .) y 6y. o »=n—r+1 m=n—r+1 0 0 m=n—r+1 0 tmj R = -E Kych- E K^ch- o l-n-r+l l=n-r+l E K^ch- E K^Ch- 0 / = n - r + l / = n - r + l - E P£Ch- - E Pf&j o / = n - r + l l=n-r+l 119 MATRICES OF SECTION 2.3 / 1.0 120 F = j = n - r + l 0 Af"" tj E A?/ j=n—r+1 E j = n — r + 1 - E j ' =n—r+1 P ? . H = < E Mt? E Qy - E f = n - r + l j = l y = n - r + l f c = n - r + l 0 E E « 3 " E c , y J = n - r + l j ' = l L=< E K^tch-l-n-r+l j = l - E E Clj / = n - r + l j = l E Kftcij J = - I l=n-r+l ^ j = l o APPENDIX 2 NEWTON-RAPHSON PROCEDURE A2.1 N E W T O N - R A P H S O N I T E R A T I O N P R O C E D U R E The standard Newton-Raphson iteration procedure is adopted here. The scheme is outlined for a set of r simultaneous nonlinear algebraic equations for r unknowns xx, x2,.. •, xr. These equations can be written in the form fi{xi,x2)...,xr) = 0 / 2 ( i i , z 2 , . . . , i r ) = 0 (A2-1) fr( xl>x2i' --.^r) = 0 Expanding the jth equation around x 1 0 , x20,..., x r 0 using the Taylor series we obtain fj{xU X2, • • • . xr) « //(*10i *20, • • • . * ro) + ( | ^ ) where the derivatives are evaluated at (XIQ,X2Q,. .. ,xRQ) a n d A i j = xi — x 1 0 A x 2 = x 2 — X20 A l f X^ *^r0* 121 NEWTON-RAPHSON PROCEDURE / A2.2 122 Similarly, carrying out the expansion for all equations in (A2.1) and rear-ranging, we obtain ( i & ) o ( ^ ) o " ( i £ ) 0 ( & ) o •• • • • • ( i £ ) 0 ( ^ ) o A x 2 . A x r t Equation (A2.3) can be written in the form ~flixlOiX20i-'-i xro) -f2ix10,x20i--->xro) — fr(x10>x20t' •• >xro) (A2.3) [ J ] {Ax} = { - / } (A2.4) where the matrices denote the corresponding terms in equation (A2.3). The algorithm for obtaining the solution vector { x } is outlined below: (1) Guess the solution vector {X}Q. (2) Compute [J]. (3) Compute { - / } . (4) Obtain { A x } from [ J ] {Ax } = { - / } . (5) Obtain {x}x = { x } 0 + { A x } . (6) If {x}i is close to { x } 0 within the specified tolerance, go to (9). (7) Set { x } 0 = { x } x . (8) If [J] needs to be modified in this iteration go to (2), otherwise go to (3). (9) Set {x } = { x } i and stop. A 2 . 2 D E T A I L S O F M A T R I C E S O F S E C T I O N 2.4.2 The steady state solution of the fluid-structure interaction problem un-der consideration corresponds to the singular points of the autonomous system of equa-tions (2.32), that is, when B = C = 0. This results in a set of nonlinear algebraic equations for A , B and C which are solved using the Newton-Raphson iteration procedure described NEWTON-RAPHSON PROCEDURE / A2.2 123 in section A2.1. The equations obtained for the increment to the solution vector are of the form of equation (A2.4) or equation (2.33) where [J} = A B C p e r QUI {Ax} = I AA AB AC {-/} = K PQ BR ' [ P - Q ] E - M fR M K P{ ' ( A ) A l + + C^Cl/l) + 6!.k (AJAJ + B°B°k/2 + C ; C « / 2 ) *5* [ A ? A l + B 1 B U 2 + c l C V 2 ) + 6!.k (A«AI + B;BI/2 + c ; q / 2 ) 6?jk i^Bl + AlBJ) + Sfjk (A«BI + AlBfj % (AV *5* 0 The details of submatrices of [J] are given below: °imjAm 4 - 8X A* + °imjAm A = K + B 61 + SX- A* tjm m 'ijmA™ cx *u °imjAm + 6y AVM + &ijmAm 0 NEWTON-RAPHSON PROCEDURE / A2.2 fc=n-r+l 6?,-tkj B = fi 2 + E Sf, k=n-r+l 0 'ijk 6? B L tmj rn Jk=n-r+ l t E **/ * = n - r + l 0 0 0 1 tmj m + S i j m B m 6X. B L tjm m0 tjm »i + * ? j m B V m 0 0 2 2 5? - C l tmj m + SfmjCm + t'jmCm 6* C L tjm m0 6?. C l tjm m tiimjCm + 6? im; m 'tjm rn 0 0 E % * = n - r + l E k-n-r+i 0 0 tmj rn + S i m j B m + S i j m B m S i j m B m f • E *, fc=n-r+l 5?. fl« tjm rn S i m j B m + 6? B L ' tmj m + B L 1 tjm m 0. 0 0 0 NEWTON-RAPHSON PROCEDURE / A2.2 125 r6? C i tmj tn 9 = fiR + fi tjm tn 0 6?. C* tjm tn ex fiu + 6? CL ' tmj m + 6?- Cl 0 0 0 I = A where m is summed from 1 to (n — r). APPENDIX 3 MATRIX OF SECTION 3.2 The matrix Z(t) in equations (3.3) is given by Z{t) = K + 0 6LjVm + 6 L j u m + Sijm um tjm Tr* 6f. um tjm *n Simj um + 6imj Vm + 6ijm Vm 0 126
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Nonlinear analysis of rigid body-viscous flow interaction
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Nonlinear analysis of rigid body-viscous flow interaction Pareshkumar Gordhandas, Pattani 1986
pdf
Page Metadata
Item Metadata
Title | Nonlinear analysis of rigid body-viscous flow interaction |
Creator |
Pareshkumar Gordhandas, Pattani |
Publisher | University of British Columbia |
Date Issued | 1986 |
Description | This thesis decribes the work on extending the finite element method to cover interaction between viscous flow and a moving body. The problem configuration of interest is that of a two-dimensional incompressible flow over a solid body which is elastically supported or alternatively undergoing a specified harmonic oscillation. The problem addressed in this thesis is that of an arbitrarily shaped body undergoing a simple harmonic motion in an otherwise undisturbed fluid. The finite element modelling is based on a velocity-pressure primitive variable representation of the Navier-Stokes equations using curved isoparametric elements with quadratic interpolation for velocities and bilinear for pressure. The problem configuration is represented by a fixed finite element grid but the body moves past the grid. The nonlinear boundary conditions on the moving body are obtained by expanding the relevant body boundary terms to first order in the body amplitude ratio to approximate the velocities at the finite element grid points. The method of averaging is used to analyse the resulting periodic motion of the fluid. The stability of the periodic solutions is studied by introducing small perturbations and applying Floquet theory. Numerical results are obtained for three different body shapes, namely, (1) a square body oscillating parallel to one of its sides, (2) an oscillating circular body and (3) a symmetric Joukowski profile oscillating parallel to the line of symmetry. The latter case is considered to investigate the flow pattern around an asymmetrical body. In all cases, results are obtained for steady streaming, instantaneous velocity vectors in the fluid domain, added mass, added damping, added force and stability of the flow. A comparision is made between the numerical and published experimental steady streamlines. Very good agreement is obtained for the basic nonlinear phenomenon of steady streaming. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-07 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0062904 |
URI | http://hdl.handle.net/2429/27181 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1986_A1 P37_7.pdf [ 11.94MB ]
- Metadata
- JSON: 831-1.0062904.json
- JSON-LD: 831-1.0062904-ld.json
- RDF/XML (Pretty): 831-1.0062904-rdf.xml
- RDF/JSON: 831-1.0062904-rdf.json
- Turtle: 831-1.0062904-turtle.txt
- N-Triples: 831-1.0062904-rdf-ntriples.txt
- Original Record: 831-1.0062904-source.json
- Full Text
- 831-1.0062904-fulltext.txt
- Citation
- 831-1.0062904.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062904/manifest