A MODAL FINITE ELEMENT METHOD FOR THE NAVIER-STOKES EQUATIONS by ZLATKO SAVOR B . S c . ( E n g . ) , U n i v e r s i t y o f Zagreb, 1969 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE i n THE FACULTY OF GRADUATE STUDIES Department o f C i v i l Engineer ing We accept t h i s t h e s i s as conforming to the r equ i red standard THE UNIVERSITY OF BRITISH COLUMBIA September 1977 (aT) Z l a t k o Savor , 1977 In p r e s e n t i n g t h i s t h e s i s in p a r t i a l f u l f i l m e n t o f the requirements fo r an advanced degree at the U n i v e r s i t y of B r i t i s h Columbia, I agree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e fo r reference and study. I f u r t h e r agree tha t permiss ion fo r e x t e n s i v e copying of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the Head of my Department or by h i s r e p r e s e n t a t i v e s . It i s understood that copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l ga in s h a l l not be a l lowed without my w r i t t e n p e r m i s s i o n . V Department of C i v i l Engineering The U n i v e r s i t y of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date September .30, 1977 i i ABSTRACT A modal f i n i t e element method i s presented f o r the steady s t a t e and t r a n s i e n t analyses o f the plane f low of incompress ib le Newtonian f l u i d . The governing r e s t r i c t e d f u n c t i o n a l i s d i s c r e t i z e d w i t h a high p r e c i s i o n t r i a n g u l a r stream f u n c t i o n f i n i t e element. Eigenvalue a n a l y s i s i s c a r r i e d out on the r e s u l t i n g d i s c r e t i z e d problem, under the assumption tha t the non l inea r convec t ive term i s equal to z e r o . A f t e r t r u n c a t i n g at va r ious l e v e l s o f approximat ion to o b t a i n a reduced number o f modes, the t r ans fo rmat ion to the new vec to r space, spanned by these modes i s performed. Advantage i s taken o f the ..symmetric and the an t i symmetr ic p r o p e r t i e s o f the modes i n order to s i m p l i f y the c a l c u l a t i o n s . The Lagrange m u l t i p l i e r s technique i s employed to {incorporate the nonhomo-geneous boundary c o n d i t i o n s . The steady s t a t e a n a l y s i s i s c a r r i e d out by u t i l i z i n g the Newton-Raphson i t e r a t i v e procedure. The a l g o r i t h m f o r t r a n s i e n t a n a l y s i s i s based upon backward f i n i t e d i f f e r e n c e s i n t ime . Numerical r e s u l t s are presented f o r the f u l l y developed plane P o i s e u i l l e f l o w , the f low i n a square c a v i t y , and the f low over a c i r c u l a r c y l i n d e r problems. These r e s u l t s c f o r the steady s t a t e are compared w i t h the r e s u l t s obta ined by d i r e c t f i n i t e element approach on the same g r i d s and the r e s u l t s obta ined by f i n i t e d i f f e r e n c e s technique . I t i s concluded tha t the number o f modes, which are to be r e t a i n e d i n the a n a l y s i s i n order to ach ieve reasonable r e s u l t s , inc reases w i t h the ref inement o f the f i n i t e element g r i d . Furthermore, the cho i ce of modes to be used depends on the problem. F i n a l l y i t i s e s t a b l i s h e d , tha t t h i s new modal method i i i y i e l d s good r e s u l t s i n the range of moderate Reynolds numbers w i t h about 50% or l e s s o f the modes of the problem. T h i s , i n t u r n , means tha t the time i n t e g r a t i o n s can be performed on a g r e a t l y reduced number o f equations and hence p o t e n t i a l savings i n computer t ime are s i g n i f i c a n t . i v TABLE OF CONTENTS Page ABSTRACT i i TABLE OF CONTENTS • i v LIST OF TABLES v i LIST OF FIGURES v i i NOTATION i x ACKNOWLEDGEMENTS x CHAPTER 1. INTRODUCTION . 1 2. BASIC THEORY 3 2.1 P a r t i a l D i f f e r e n t i a l Equation 3 2.2 Restri c t e d V a r i a t i o n a l P r i n c i p l e 4 3; FINITE ELEMENT FORMULATION . . . . . . . . 8 3.1 Functional on Element Level . . . 8 3.2 Global Functional . 11 4- MODAL APPROACH 14 4.1 Transformation by Condensation . 16 4.2 Transformation i n Conjunction with Lagrange M u l t i p l i e r s Technique 20 4.3 F i n a l Remarks . 22 5; SOLUTION OF NONLINEAR EQUATIONS IN TIME DEPENDENT AND STEADY STATE ANALYSES 24 5.1 Time Dependent Approach . . 25 5.2 Steady State Approach 29 5.3 F i n a l Remarks 31 V CHAPTER 6. PRESSURE FIELD 33 7. NUMERICAL IMPLEMENTATION 36 8. EXAMPLE APPLICATIONS 40 8.1 F u l l y Developed Plane P o i s e u i l l e Flow . 40 8.2 C i r c u l a t o r y Flow i n a Square Cavity 42 8.3 Flow Around a C i r c u l a r Cylinder . 48 9. CONCLUSIONS 59 BIBLIOGRAPHY 61 v i LIST OF TABLES Table(s) Page I Number of Symmetric and Antisymmetric Modes and Lowest and Highest Eigenvalues f o r Three F i n i t e Element Grids ' Used f o r Square Cavity Flow Problem Simulation . . . . 64 II Error i n Stream Function at MidnodefS for Various Number of Modes and Various Reynolds Numbers for SQCA 12-13 Grid 65 III Comparison of Numerical Results for Square Cavity Flow Problem 66 TV Comparison of Numerical Results f o r Steady Flow Around a C i r c u l a r Cylinder 73 v i i LIST OF FIGURES Figure(s) Page 1 Local ( C,n ) Coordinates of the Triangular Element . . . 74 2 The 18 DOF Truncated Quintic Element 74 3 P o i s e u i l l e Flow Problem Configuration 75 4 Mode Shapes & Eigenvalues f or P o i s e u i l l e Flow . . . . . . 76 5 Streamlines & V e l o c i t i e s at Midnode 3 for P o i s e u i l l e Flow 77 6 C i r c u l a t o r y Flow Induced by a Moving L i d Over a Square Cavity 78 7 F i n i t e Element Grids f o r Square Cavity Flow 79 8-1 Mode Shapes & Eigenvalues for SQCA 12-13 Grid 80 8-2 F i r s t 30 Mode Shapes & Eigenvalues for SQCA 36-29 Grid . . 81 9 Streamlines for Square Cavity Flow f o r Various Reynolds Numbers R 83 10 E q u i - V o r t i c i t y Lines f o r Square Cavity Flow for Various Reynolds Numbers R 91 11 Equi-Pressure Lines f o r Square Cavity Flow for Various Reynolds Numbers R . . 99 12-1 F i n i t e Element Grids for a Flow Around a C i r c u l a r Cylinder CYLFL 84-58 Grid 107 12- 2 F i n i t e Element Grids f o r a Flow Around a C i r c u l a r Cylinder CYLFL 92-63 Grid 108 13- 1 F i r s t 15 Antisymmetric Mode Shapes & Eigenvalues f o r CYLFL 84-58 Grid 109 13-2 F i r s t 20 Antisymmetric and F i r s t 10 Symmetric Mode Shapes & Eigenvalues for CYLFL 92-63 Grid 110 v i i i Figure(s) Page 14 Streamlines f o r Steady Flow Around a C i r c u l a r Cylinder for Various Reynolds Numbers R 112 15 E q u i - V o r t i c i t y Lines f o r Steady Flow Around a C i r c u l a r Cylinder for Various Reynolds Numbers R . . . 121 16 Streamlines for Transient Flow Around a C i r c u l a r Cylinder at Reynolds No. R = 140 for Various Time Instants . . . 130 ix NOTATION The s p e c i f i c use and meaning of symbols are defined i n the text where they are introduced. I n d i c i a l notation based upon the summation convention i s adopted throughout the text. A l t e r n a t i v e l y , the more common matrix notation i s also used when i t r e s u l t s i n equations written more concisely, and for f i n a l expressions. Vector quantities are indicated by a lower bar, matrices by two lower bars, and t r i p l y subscripted arrays by three lower bars. When the arrays have to be written out i n f u l l , vectors are denoted by { } brackets, matrices by [ ] brackets, and trans-posed vectors by < > brackets. A comma followed by an index, appearing as a subscript, designates a p a r t i a l d e r i v a t i v e with respect to a s p a t i a l d e r i v a t i v e i n the d i r e c t i o n of that index. The Greek symbol e implies 'belongs to' unless s p e c i f i e d otherwise. X ACKNOWLEDGEMENTS The author wishes to express h i s gratitude to h i s adviser, Dr. M.D. Olson, for h i s invaluable advice and guidance during the research and preparation of t h i s t h e s i s . He also expresses thanks to Dr. S.Y. Tuann for h i s valuable suggestions during the research work. The f i n a n c i a l support of the Research Council of SRH, Yugoslavia, and the Department of C i v i l Engineering at UBC i s g r a t e f u l l y acknowledged. F i n a l l y , many thanks to my wife Iva for a l l her help and encourage-ment . CHAPTER 1 Introduction The numerical analysis of the plane flow of incompressible Newtonian f l u i d constitutes an important area i n computational mechanics and engin-eering p r a c t i c e . T r a d i t i o n a l l y , the most popular and thoroughly studied methods for tr e a t i n g problems of t h i s type have been the various techniques based on f i n i t e d i f f e r e n c e d i s c r e t i z a t i o n . The great progress i n a p p l i c a t i o n of f i n i t e element procedures has been made only i n the l a s t few years. [1] . The purpose of t h i s thesis i s to e s t a b l i s h whether i t i s f e a s i b l e to reduce the number of d i s c r e t e v a r i a b l e s appearing i n f i n i t e element approach by performing the transformation to modal coordinates. This i s e s p e c i a l l y important i n applications to transient problems where the cost of integrating the d i s c r e t e equations i n time soon becomes p r o h i b i t i v e l y expensive. Stream function alone approach i s used for f i n i t e element modelling. The r e s t r i c t e d functional governing the problem i s d i s c r e t i z e d with a high p r e c i s i o n f i n i t e element of C 1 c l a s s . Eigenvalue analysis i s c a r r i e d out on the r e s u l t i n g d i s c r e t i z e d problem under the assumption that the nonlinear convective term i s equal to zero. A f t e r truncating at various l e v e l s of approximation to obtain a reduced number of modes, transformation to the new basis defined by these modes i s performed. It i s shown that i t i s computationally more e f f i c i e n t to employ transformation i n conjunction with Lagrange m u l t i p l i e r s technique than transformation by condensation, even though the system s i z e i s expanded and the system matrix i s nonpositive-d e f i n i t e . Steady state analysis i s c a r r i e d out by u t i l i z i n g the Newton-Raphson i t e r a t i v e procedure. Algorithms for transient analysis are based upon.backward f i n i t e differences i n time and d i f f e r only i n how the nonlinear 2 term i s treated. The pressure d i s t r i b u t i o n i s obtained by solving the Poisson equation a f t e r a l l modal stream function variables i n the f i n i t e element coordinates have been found. Numerical r e s u l t s are presented and discussed for the f u l l y developed P o i s e u i l l e flow, the flow i n the square cavity and the flow over a c i r c u l a r cylinder problems. The r e s u l t s f o r the steady state are compared with the r e s u l t s obtained by d i r e c t f i n i t e element approach on the same grids and the r e s u l t s obtained by f i n i t e d i f f e r e n c e procedures. 3 CHAPTER 2 Basic Theory In t h i s chapter the p a r t i a l d i f f e r e n t i a l equation governing the problems of plane flow lafi incompressible Newtonian f l u i d and the r e s t r i c t e d v a r i a t i o n a l p r i n c i p l e derived from i t are presented to re c a p i t u l a t e the t h e o r e t i c a l basis on which modal approach s h a l l l a t e r be b u i l t . 2.1 P a r t i a l D i f f e r e n t i a l Equation The d e s c r i p t i v e set governing the plane flow of incompressible Newtonian f l u i d consists of: (a) the s p a t i a l equation of motion, derived from the dynamics of flow considerations using the p r i n c i p l e of conserva-t i o n of l i n e a r momentum, and c a l l e d Navier-Stokes equation: 3u. + u u = -p,i, + — (u + u ) , i xefi, t>0 (1) pt 3 !>3 1 K e !»3 3* 1 3 (b) the continuity equation derived from the kinematics of f l u i d flow v i a the conservation of mass: V • u = u. . = 0 xefi, t>0. (2) x , i -The appropriate boundary and i n i t i a l conditions also have to be prescribed. The equations are w r i t t e n i n an Eulerian frame of reference f i x e d i n space. The normalizing system used i n above equations i s based on the convective time scale L/U, where L i s a c h a r a c t e r i s t i c length and U i s a c h a r a c t e r i s t i c v e l o c i t y of the problem. Pressure i s normalized with respect to the reference pressure pU 2. Re i s the dimensionless Reynolds number defined as Re = ^p-. The other parameters of the problem are the density of the f l u i d p, which i s assumed constant and independent of temperature, the absolute v i s c o s i t y of the f l u i d y, and the kinematic v i s c o s i t y defined as v = —. u . ( i = 1,2) are the components of v e l o c i t y i n x and y d i r e c t i o n s P i r e s p e c t i v e l y , and 0, i s the domain under consideration with the boundary V : Viscous stresses for an i s o t r o p i c Newtonian f l u i d are defined a x. . = y(u. . + u. .) (3) i j i , J l o i By introducing the stream function concept the continuity equation becomes exactly s a t i s f i e d and the pressure i s completely e l i -minated, leaving the stream function as the only dependent v a r i a b l e . The stream function i s defined as: f = u f = -v, (4) y x Then by c r o s s - d i f f e r e n t i a t i n g ,Eq. -:(1). J f or each of the two d i r e c t i o n s x and y to eliminate the pressure and by expressing the v e l o c i t i e s i n terms of the stream function v i a Eq. (4), the following s i n g l e - v a r i a b l e fourth order p a r t i a l d i f f e r e n t i a l equation i s obtained: ^ - - + * h 'V <« As 2.2 R e s t r i c t e d V a r i a t i o n a l P r i n c i p l e As proven by Finlayson [2], no exact v a r i a t i o n a l p r i n c i p l e e x i s t s for t h i s equation due to the non-selfadjoint convective terms. Therefore Olson [3] resorted to r e s t r i c t e d v a r i a t i o n a l p r i n c i p l e . Yamada et a l . [4] point out that the paper [3] was " i t s e l f a break-through i n the a p p l i c a t i o n of the f i n i t e element method to f i e l d of f l u i d dynamics". The d e r i v a t i o n of the p r i n c i p l e i s repeated here, extended for the time dependent term, as t h i s p r i n c i p l e y i e l d s the f u n c t i o n a l , which i n i t s d i s c r e t i z e d form serves as the basis for the modal approach. If Eq. (5) i s m u l t i p l i e d by Sty and integrated over the domain Q the following varied i n t e g r a l i s obtained: 61 = Jf[±- vSjj + * VIM) - + V 2 ( i ) - -f- (vH)]6^dQ. (6) 2 ^ K e x y y x dt Integrating Eq. (6) by parts twice for the viscous term and once for the convective and the time dependent terms and employing the Green-Gauss theorem y i e l d s the following v a r i a t i o n a l statement: +Vip0V6i!>]'dP. + f[^~ d ^ 2 ^ - 72iJ)0nx7i|>0 (7) where the notation r e f e r s to v a r i a t i o n i n 1J1 only, while IJJ° i s kept f i x e d . A f t e r the v a r i a t i o n i s set equal to IJJ° and the governing d i f f e r e n t i a l equation (5) recovered as Euler equation. If the v a r i a t i o n a l operation 6 i s p u l l e d i n front of the i n t e g r a l expression, and the boundary i n t e g r a l s neglected for now, the following f u n c t i o n a l i s obtained: = // [-5iT ( V 2* 2) + ( A 2 f ) i|/ (8) 0 Z K y x - (i))°V2f) + Vif)°Vi|)]dfi. 6 This f u n c t i o n a l y i e l d s the boundary conditions: 1 H]pn) either &\b = 0 or - — V 2 T J J + V 2 i M + — = 0 (9) Re n s 3t and e i t h e r S\\i = 0 or z, = V 2 I J J = 0 (10) where n i s the d i r e c t i o n of the outward normal on the boundary, San. u i s the d i r e c t i o n tangential to the boundary, and t, represents the v o r t i c i t y . The l e f t hand side conditions are the l o g i c a l ' r i g i d ' boundary conddta.6n>ss corresponding to s p e c i f i e d stream function and i t s normal d e r i v a t i v e . The 'natural' boundary conditions on the r i g h t hand side are not the ones that we would l i k e to have, namely those of constant pressure and zero shear s t r e s s . These boundary conditions can be e a s i l y implemented though, by simply adding the appropriate boundary i n t e g r a l s to the fu n c t i o n a l Eq. (8). As we use tr i a n g u l a r f i n i t e elements to be covered i n more d e t a i l l a t e r , the boundary i n t e g r a l s are added i n such a way, that the natural boundary conditions can be approximated only on one edge of the t r i a n g l e . This e f f e c t i v e l y means that only i f a p a r t i c u l a r f i n i t e element i s to be used as boundary element with two of i t s v e r t i c e s l y i n g on the part of the boundary, where any one or both 'natural' boundary conditions are prescribed, these boundary i n t e g r a l s are added. In order to sim p l i f y the equations we switch now to l o c a l coordinates £ ,n as shown i n F i g . (1), and obtain the modified governing f u n c t i o n a l as: i e 0 M ° ) = / / ( v 2 r0 - (4>°vV) ^ n + O ^ v V ) 4>5 (ID + v^v*] d ? d n + r [iL. * j * - i r r + ^ ° n ) * ] n = 0 d 5 . - D 7 The f u n c t i o n a l y i e l d s the ' r i g i d ' boundary conditions; 6^ = 0 and <5^ = 0 n (12) and the 'natural' boundary conditions: nnn ) + \i> 4v - = 0 (13) and T = Re v%n = 0. Both ' r i g i d ' boundary conditions can be exactly s a t i s f i e d for a s t r a i g h t boundary, since stream function i s uniquely determined by and at the two v e r t i c e s l y i n g on the boundary, and itsnnormal d e r i v a t i v e , a cubic, i s uniquely determined by and ij; at these v e r t i c e s . For curved boundaries, however, s p e c i a l provisions are necessary, and the desired boundary conditions are r e a l i z e d exactly only at d i s c r e t e nodes. These provisions i n the element area integrations are included i n t h i s work, but they w i l l not be reported, as they are well documented i n Tuann and Olson [5,6]. 8 CHAPTER 3 F i n i t e Element Formulation By performing the aforementioned i n t e g r a t i o n by parts the order of the function space i s lowered. This r e l a x a t i o n of the continuity require-ments on the stream function i s referred to as a 'weak' formulation of the problem. I t permits us to seek the s o l u t i o n i n the Sobolev space W^ , which contains a l l functions, whose second derivatives are square integrable, i . e . f i n i t e . In other words using c l a s s i c a l f i n i t e element vocabulary we need a f i n i t e element of c l a s s . Such an element s a t i s f i e s the 'compati-b i l i t y ' conditions, as i t provides interelement continuity of \j> and i t s normal d e r i v a t i v e and the governing fu n c t i o n a l contains derivatives of at most second order. In t h i s chapter the de r i v a t i o n of the fun c t i o n a l f o r an element dEclass i s presented. The procedure of f i n d i n g the extremum of the global functional r e s u l t i n g i n the d i s c r e t i z e d equations of- motion i s indicated. 3.1 Functional on Element Level Here the tria n g u l a r element derived by Cowper et a l . [7] i s adopted. This element has since been recognized as one of the most accurate f i n i t e elements a v a i l a b l e f o r plate bending problems, Gallagher [8]. I t should be noted, though, that t h i s element requires a su b s t a n t i a l formu-l a t i v e e f f o r t . Olson was the f i r s t to modify i t for a p p l i c a t i o n to the steady state two-dimensional viscous flow problems. Only the relevant steps i n the element d e r i v a t i o n areiaincluded and for more d e t a i l s i t i s referred to Olson [1,3]. The element and i t s nodal variables are shown i n Fig. (2). The f i e l d v a r i a b l e i s interpolated within the element by a truncated q u i n t i c polynomial: 2 0 m. n. i|> = E a.?" V 1 ( 1 4 ) 1 = 1 1 such that tangential v e l o c i t y component along the edge 1 - 2 i s a cubic i n £. Eighteen r e l a t i o n s expressing nodal variables ^ i n terms of polynomial c o e f f i c i e n t s A can be established by d i f f e r e n t i a t i n g ip with respect to the E,,r) coordinates and by s u b s t i t u t i n g the l o c a l coordinates of the v e r t i c e s . Two a d d i t i o n a l r e l a t i o n s are obtained by constraining the tangential v e l o c i t y to vary as a cubic along the remaining two edges, 2 - 3 and 3 - 1 . This i s written i n matrix form as : < ^ , 0 , 0 T > = T A ( 1 5 ) Inverting the above equation the vector of polynomial c o e f f i c i e n t A can be expressed i n terms of the vector of nodal variables jp£ as: A = | " 1 < ^ I , 0 , 0 > T = T 2 ^ ( 1 6 ) where i s 2 0 x 1 8 matrix consisting of the f i r s t 1 8 columns of T The introduction of the transformation from l o c a l coordinates £,n, to global coordinates x,y i n the form: Jfe-L = R * ( 1 7 > where ip i s the vector of nodal variables i n the global coordinate system into Eq. ( 1 6 ) y i e l d s the following r e l a t i o n between jp and A: A = T 2 R I ( J = S jfc. ( 1 8 ) 10 The element matrices are then obtained by su b s t i t u t i n g Eq. (14) into the governing fu n c t i o n a l Eq. (11) and carrying out the integrations. It should be noted, that the terms i n the functional are f i r s t arranged i n such a way, that the r e s u l t i n g arrays are symmetric i n t h e i r f i r s t two indices. The fun c t i o n a l on the element l e v e l can then be written i n d i s -c r e t i z e d form as: A M T ) = ii- j ^ v j + Q i J k ^ ? k + <3 V5 (19) where K?. = S , S .k (20) k3 rk S3 rs i s the element l i n e a r d i s s i p a t i v e matrix Q?., = S . S .S art . (21) ^13 k r.k S3 t f ^ r s t i s the element nonlinear convective matrix and Mf. = S , S .1 (22) K3 rk S3 rs i s the element consistent mass matrix, and ij; represents a time de r i v a t i v e . In the above equations i , j , k take values from 1 to 18, and r , s , t from 1 to 20, resp e c t i v e l y . The other arrays appearing i n these equations are defined as: k = m m 0Cm~l7(10-(m A-1);LDX -F(m + m- = 4-t n )+ n ) rs r s r r • 3 v s r V r s r ' r s + n n (n - 1) (n - 1) x F(m + m , n + n - 4Y r s r s r s r s + [m n (m - 1) (n - 1) + m n (m - 1) (n - 1) ] r s r s s r s r x F(m +m - 2 , n + n - 2) + [m n (m - 1 ) r s r s r s r + m n ( m - l ) ] x G ( m + m -•2, n + n - 1 ) s r s r s r s (23) 11 q r s t = \ ( n r m t " m r I 1 t ) m s ( m S " 1 } + ( n S m t " m S n t ) m r ( m r " 1 } x F(m +m + m - 3 , n + n + n ) + ( n m - m n ) n. ( n - 1 ) r s t ' r s t 2 v r t r t s s + ( n ni - m n ) n ( n - l ) x F ( m +m + m - 1, n + n + n - 3 ) s t s t r r v r s t r s t - m m (m + m + m - 2 ) x G ( m +m + m - 3 , n + n + n ) 2 r s r s t r s t r s t - i n n (m + m ) x G(m + m + m - l , n + n + n - 2 ) (24) 2 r s r s r s t r s t £ = m m x F ( m +m - 2 , n + n ) + n n x F ( m + m , n + n rs r s r s r s r s r s r s - 2) (25) where F and G represent the exact i n t e g r a l formulas obtained by carrying out the integrations of the general term £ m n n over the area of the t r i a n g l e and along the n = 0 edge, respectively. w \ m + ± r / / i , \m+ln min! M n F(m,n) = c [(a) - (-b) J ( m + n + 2)! ^ ' G(m,n) = ^ [ ( a ) m + 1 - (-b) 1^ 1] i f n = 0 (27) G(m,n) = 0 i f n + 0. 3.2 Global Functional The global governing f u n c t i o n a l f o r the enti r e problem i n d i s c r e t i z e d form i s obtained by following the usual f i n i t e element assem-blage process. That i s , the element matrices are appropriately summed into the glo b a l matrices, taking into account the symmetry and handedness properties. The homogeneous boundary conditions are introduced during the process as w e l l , by simply s t r i k i n g out rows and columns of the global 12 arrays corresponding to the variables on the boundary with s p e c i f i e d zero value. The problem obtained i n that manner, that i s , the problem of the order equal to the number of the nonzero nodal variables w i l l be referred to as the 'net' problem. It should be noted, that so f a r , the nonhomogeneous boundary conditions have not been taken into account. The aforementioned f u n c t i o n a l can be written as: K M H = h V A + W W k + V 3*k i>j»k 1, . . . , r (28) where r i s the''net' problem s i z e , and iji and represent the g l o b a l vectors of nonzero nodal v a r i a b l e s . The global d i s s i p a t i v e matrix and the g l o b a l consistent mass matrix M^... are each stored columnwise i n an one-dimensional array of s i z e l b . x r. The nonlinear convective matrix Q „ k i s also stored i n f u l l , because i t was found to be numerically more e f f i c i e n t , than to r e c a l c u l a t e i t element by element. The matrix cannot be accommodated i n the core, because of i t s s i z e , so i t i s stored s l i c e by s l i c e columnwise i n the order of the index k on the high speed disk. Each one-dimensional s l i c e i s of s i z e lb x (2 lb - 1), or lb x r, depending on whether the 'net' s i z e of the problem i s bigger or smaller than the bandwidth 2 l b - 1. The half-bandwidth of the problem i s defined as l b . The procedure of seeking the extremum of the g l o b a l f u n c t i o n a l Eq. (28) with respect to while keeping f i x e d , y i e l d s : u1 = (i V J + v * ; ^ 0 + v? ^ • ° (29) where the i n d i c e s , appearing as subscripts, i , j , and k take on, successively, values from 1 to r , as before. The r e s t r i c t e d v a r i a t i o n a l p r i n c i p l e now permits to replace 13 jf)° by and as Sip^ i s a r b i t r a r y , the set of f i r s t order nonlinear d i f f e r e n -t i a l equations i n di s c r e t e form f or the r gl o b a l variables ip i s obtained: K . ip . + Q. .. ip . ip . + M . ip . = 0 (30) where the indices take on successively the same values, as i n the previous two equations. The f u n c t i o n a l Eq. (28), or a l t e r n a t i v e l y the system of d i f f e r e n -t i a l equations Eq. (30),, w i l l serve as the basis f o r employing mbdauVvapproach. -14 CHAPTER 4 Modal Approach The idea to employ modal approach originates i n the f i e l d of s t r u c t u r a l dynamics where t h i s approach has been applied s u c c e s s f u l l y f o r years to l i n e a r and mildly nonlinear problems and recently extended to cover strongly nonlinear problems. In the category of l i n e a r s t r u c t u r a l dynamics problems the advantages of modal approach are quite obvious. Not only do the equations of motion uncouple, under the assumption that damping matrix can be repre-sented as a l i n e a r combination of s t i f f n e s s and mass matrices, but also, due to the f a c t that most of the frequency content of the loading i s con-tained i n the lowest modes f o r many types of p r a c t i c a l loadings, only a f r a c t i o n of the t o t a l number of uncoupled equations i n generalized coordinates need be considered i n order to obtain a reasonable approximation to the actual response of the system [9]. We define Emidifclhy: nonlinear systems as those systems for which nonlinear deformation mechanisms do not cause major changes i n the d e f l e c t i o n patterns. For these systems the equations of motion are no longer uncoupled due to off-diagonal terms appearing i n the generalized s t i f f n e s s matrix. The response may s t i l l be evaluated, though, by d i r e c t i n t e g r a t i o n of a l i m i t e d set of equations of motion i n generalized coordinates. In the category of strongly nonlinear problems modal analysis has been t r i e d , to my knowledge, only by N i c k e l [10]. He casts the dynamic equations i n the incremental form. Af t e r f i n d i n g the i n i t i a l modes and frequencies he proceeds to compute the subsequent modal spectrum for non-l i n e a r states, employing an extremely fast and e f f i c i e n t eigenproblem solver that involves matrix m u l t i p l i c a t i o n s only and uses the most recently computed 15 spectrum as i n i t i a l estimate. So he obtains the equations of motion i n the uncoupled form even when strong n o n l i n e a r i t i e s are present. The advantage of h i s procedure gets more pronounced as the bandwidth of the problem i n the o r i g i n a l f i n i t e element coordinates increases. I t should be noted though, that again the assumption i s made, that the lowest natural frequencies and associated modes dominate the incremental motion. This assumption, which i s e s s e n t i a l l y equivalent to the previous statement, that the frequency content of the loading i s contained the lowest modes, i s empirically v e r i f i e d f o r great many s t r u c t u r a l dynamics problems, notably those of earthquake e x c i t a t i o n . Had the same been true for the problems of f l u i d dynamics, we would have t r i e d to perform the modal decomposition for nonlinear states as w e l l , employing only the lowest modes. Unfortunately, p r i o r to t h i s work nothing was known about the a p p l i c a t i o n of modal analysis to f l u i d s . So the f i r s t fundamental question that we have attempted to answer was how many modes were to be included i n order to obtain a reasonable s o l u t i o n . That i s why we i n t e n t i o n a l l y r e s t r i c t e d the class of the problems that we have been t r y i n g to solve to mildly nonlinear problems, that i s , to the range of low to moderate Reynolds numbers. Thus we needed to perform the l i n e a r eigenvalue analysis only, but f o r a l l the modes. The eigenvalue analysis f o r the subsequent nonlinear states, based upon the tangent matrix obtained as a l i n e a r combination of the l i n e a r viscous matrix and the non-l i n e a r convective matrix, evaluated with stream function values at a previous time step, for a l l the modes, would have been p r o h i b i t i v e l y ex-pensive, of course. The second fundamental question that we have t r i e d to answer was, what was the maximum Reynolds number governing the nonlinear behaviour, for which the l i n e a r modes could s t i l l p redict the accurate enough s o l u t i o n . 16 In t h i s chapter two procedures f o r employing modal a n a l y s i s , transformation by condensation and transformation i n conjunction with Lagrange m u l t i p l i e r s technique w i l l be presented. Both procedures s h a l l i n turn be applied, the f i r s t one to Eq. (30), the second one to the f u n c t i o n a l Eq. (28). I t w i l l be assumed, as already stated previously, that the homo-geneous boundary conditions have been taken account of, but not the non-homogeneous ones. As a matter of f a c t , as must be obvious from Eq. (30), which i s homogeneous, the only 'loading' comes p r e c i s e l y from these non-homogeneous boundary conditions, that i s , from the prescribed nonzero variables on the boundary. 4.1 Transformation by Condensation In the following d e r i v a t i o n the indices appearing as subscripts i , j , and k take on successively values from 1 to s, where s i s the t o t a l number of unconstrained variables equal to the s i z e of the problem, which w i l l be referred to as the 'net net' problem, the indices m and n values from 1 to p, where p Is the number 6f constrained variables, on the kinematic and mixed boundaries, and the indices a,3, and y values from 1 to w, where w i s the number of the eigenvectors used. Throughout the de r i v a t i o n i t w i l l be assumed that the number of the eigenvectors used w i s equal to the s i z e of the 'net net' problem, so that the equations of motions i n the generalized coordinates are the exact equivalent of the same equations i n f i n i t e element coordinates. Later i n actual applications we s h a l l obviously attempt to reduce the number of eigenvectors used. The equations of motion w i l l be written for the unconstrained variables only as follows: 17 V j + fe V j + Q u k * l * J + V m + fe U njk n r j i m k r r m nmk n m (31) As there are p prescribed variables il> and ip the m u l t i p l i c a t i o n s m n involving these variables may be performed. Taking advantage^of the symmetry of the array Q i n the f i r s t two subscripts i and j , we can write: Q ..ip + Q. = Q ty. + Q V . = C..iJ).. (32) nik n j i m k I m njk n j njk n T j j k T j By transposing a l l the known terms to the r i g h t hand side, the following set of equations i s obtained: 11 .iji. + K , + C., M -H- .;^ Q. ip.\p. = F (33) V T J Re k j r j j k r j ^ i j k y x r j k where F^ can be viewed as the 'loading' term defined by: F. = -M. i|» - | - K. ij) - Q , t|> i|> . (34) k l m m Re km m nmkTnTm We now propose to transform Eq. (33) using the following transformation on the f i n i t e element stream function vector: * J ( t ) - V y a C t ) ( 3 5 ) where E. i s a square matrix and y (t) i s a time dependent vector of 3a a the order s. The matrix E i s established by solving f o r : V J + V J = ° * ( 3 6 ) The s o l u t i o n can be postulated to be of the form: 18 * =+ e'^ (37) where t i s the time var iab le and X i s a constant defined to represent the time decay of the vector <j>. By subs t i tu t ing Eq. (37) in to Eq. (36) the generalized eigen-value problem i s obtained, from which the eigenvectors cj) and the associated eigenvalues X are to be determined: ^ ^ - v v 0 - (38) I t i s evident now that E i s to be defined as a matrix whose columns are jus t the eigenvectors <j>: E = < i 1 , * 2 , i 3 . • • • ,£s*.- (39) By introducing Eq. (35) into Eq. (33), we obtain: 1 M. .E . y +~— K . . E . y + C . . E . y + Q E E. y y = F . kj j a a Re kj j a a j k j a a xjk i f j a a p k (40) By mul t ip ly ing the above equation by E the equations of motion kj i n the eigenvector basis are obtained i n the form: M . . E . E. y + I - K . . E . E, y + C.. E. 'E, y KJ j a ky a Re KJ j a ky a j k j a ky a + Q. .. E E . E y y . = F. E. . ^xj k xB j a ky a 3 k ky (41) The above equation can be s impl i f i ed taking into account the orthonormal properties of the eigenvectors with respect to the mass matrix M^j . These properties can be expressed as fo l lows: 19 M, .E, E. kj ky jet M, .E, E. KJ ky j a K. .E, E. KJ ky jet K, .E. E. kj ky j a The equation of motion for the o-th generalized coordinate y can then be w r i t t e n as: o y + k~ x y + C.. E..E. y + Q. E..E. E. y y = F. E. (43) ^o Re o-'o j k j X j y a i j k ig j a ko^cr 3 k ko or i n the matrix form f o r the whole system: where J. i s the i d e n t i t y matrix, y i s the vector of f i r s t order time derivatives of generalized coordinates, X i s the diagonal matrix of T eigenvalues, y i s the vector of generalized coordinates, C = E CE i s the T square matrix, Q = E QEE i s the condensed nonlinear convective matrix, T and F = E F i s the load vector. I t should be noted that the procedure l a i d out above i s completely analogous to the treatment of an a r b i t r a r y support e x c i t a t i o n i n s t r u c t u r a l dynamics. This transformation by condensation procedure seems quite appeal-ing. I t uses the c l a s s i c a l eigenvectors which s a t i s f y the homogeneous boundary conditions only. Unfortunately we have to deal with the nonlinear convective matrix, as w e l l . Although i t i s numerically quite easy to code an integer vector to keep track of the prescribed values of the stream func-t i o n on the boundary, even the process of forming the nonlinear matrix i n the f i n i t e element coordinates i s an extremely complicated one, and breaking i f y = ct = 0 i f Y ^ ct (42) cr i f Y = ctY = 0 i f Y ^ ct. 20 t h i s matrix into four d i f f e r e n t arrays, as indicated i n Eq. (31), to carry out the m u l t i p l i c a t i o n s of Eqs. (32) and (34), seemed quite hopeless. I t i s because of these numerical d i f f i c u l t i e s that the procedure had to be abandoned and we concentrated instead on the transformation i n conjunction with Lagrange m u l t i p l i e r s technique. 4.2 Transformation i n Conjunction with Lagrange M u l t i p l i e r s Technique In t h i s d e r i v a t i o n the indices i , j , and k take on successively values from 1 to r , where r i s the s i z e of the 'net' problem, the index m values from 1 to p, where p i s the number of constrained v a r i a b l e s on the kinematic and mixed boundaries, and the indices a,$, and y values from 1 to w, where w i s the number of eigenvectors used to approximate the so l u t i o n . Throughout the de r i v a t i o n w w i l l be taken equal to the s i z e of the 'net' problem r i n order to obtain the exact equivalent of the formulation of the problem i n the f i n i t e element coordinates. The p l i n e a r constraints imposed upon the stream function can be expressed i n Che form: Gmk*k " Tm " ° C 4 5 ) where G ^ i s the rectangular matrix of constraints of s i z e p x r with a l l the entries equal to zero except f o r the diagonal entries corresponding to the prescribed values of the stream function on the boundary, ip, i s the vector of stream function variables of s i z e r, and T i s the vector of m these prescribed values of the stream function on the boundary of s i z e p. Then Eq. (45), m u l t i p l i e d by the vector of Lagrange m u l t i p l i e r s h of s i z e p, i s added to the global f u n c t i o n a l Eq. (28) and the augmented fu n c t i o n a l obtained i n the form: 21 + G , h - T h = 0. mk k m mm (46) The v a r i a t i o n of the augmented f u n c t i o n a l defined by the above equation, f i r s t with respect to ip while keeping ip0 and h f i x e d and then with respect to h while keeping i|> and f i x e d , y i e l d s : 3ijj£ V ) r j Re k j r j H i j k r i r j mk m (47) 81 , h = G Ji " T = 0 9h mi ] m m J J where i n the f i r s t subset of the above set of equations has been replaced by as permitted by the r e s t r i c t e d v a r i a t i o n a l p r i n c i p l e , and i n the second subset the dummy subscript k has been replaced by j . By employing the coordinate transformation indicated i n Eq. (35), with the matrix E established v i a Eqs. (36) to (39), and by introducing t h i s transformation into Eq. (47) we obtain: M E. y +^-K..E. y + Q . E..E. y y . + G.h = 0 Tcj j a a Re kj j a a i j k l g j a a B mk m (48) G .E. y - T = 0. mj j a a m By m u l t i p l y i n g the f i r s t subset of the above set of equations by E^ and by taking advantage of the orthonormality properties of the eigenvectors with respect to the mass matrix M^ _. as expressed i n Eq. (42) the following system of equations i s obtained: y + ^ ~ ^ y + Q. ., E.„E. E, y y„ + G . E. h = 0 Y Re Y Y i j k l g j a ky a B mk ky m (49) G E. y - T = 0 mj j a y a m^ 22 or i n the matrix form: Iy + Xy + Qyy + | T h = 0 (50) Ey - T = 0 where I i s the i d e n t i t y matrix, A i s the diagonal matrix of eigenvalues, T Q = E QEE i s the transformed nonlinear convective matrix, E i s the matrix -T of eigenvectors, and E EG i s the matrix consisting of the entries i n eigenvectors at the constrained degrees of freedom. The matrices I and X are each of s i z e w x w, the matrix Q i s of s i z e w x w x w, and the matrix E i s of s i z e p x w. 4.3 F i n a l Remarks It should be emphasized again, that while i n the f i r s t method, the transformation by condensation, the general eigenvalue analysis i s run a f t e r the nonhomogeneous boundary conditions have been taken account of, i n the second method employing Lagrange m u l t i p l i e r s technique the generalized eigenvalue analysis i s run p r i o r to taking account of the given values of stream function on the boundary. Also, the f i r s t method contracts the s i z e of the system of equations to be solved to s = r - p, while the second method expandsthe s i z e to q = r +.p, where r i s the 'net' problem s i z e , and p i s the number of constraints. This does not e f f e c t our choice of the second method to a great extent, however, because i n the f l u i d problems the number of constraints i s r e l a t i v e l y small compared to the s i z e r and i t can be further reduced by our numerical procedure through the use of 'master slave' option. If the number of eigenvectors i s equal to the number of f i n i t e 23 element degrees of freedom, mathematically the same space i s spanned by the eigenvectors as by the nodal point f i n i t e element stream function variables and consequently the same s o l u t i o n must be obtained by both analyses. But, employing a l l the eigenvectors of the l i n e a r eigenvalue problem, would a c t u a l l y be a step backward, because although the global consistent mass matrix and the l i n e a r global d i s s i p a t i v e matrix become uncoupled, the equations i n generalized coordinates are s t i l l coupled through the g l o b a l nonlinear convective matrix, which becomes f u l l , whereas i t was banded i n f i n i t e element coordinates. So, the modal approach can only be more e f f i c i e n t , i f a reasonable approximation of the s o l u t i o n of f i n i t e element equations of motion Eq. (30) can be obtained by using a s i g n i f i c a n t l y reduced number of eigenvectors. We note, that so far we have only been concerned with the exact and approximate solutions of these discretetequations. Whether a good approximation to the s o l u t i o n of the a c t u a l continuum problem w i l l be obtained, depends on the f i n i t e elements employed, the f i n i t e element meshes, and the boundary conditions. 24 CHAPTER 5 Solution of Nonlinear Equations i n Time Dependent and Steady State Analyses I t i s appropriate to note here, that i n the c l a s s i c a l computational f l u i d dynamics, based on f i n i t e d i f f e r e n c e s , the studiesoof even steady state problems are mostly based on the time dependent equations, because, f i r s t l y , t h i s time dependent approach does not postulate the existence of a steady state s o l u t i o n , secondly, the procedure i s more f l e x i b l e i n the sense, that the transient s o l u t i o n can be achieved, i f so desired, and t h i r d l y and most important, the unsteady equations i n f i n i t e differences are easier to handle and more stable than t h e i r steady counterparts. The steady state s o l u t i o n i s obtained, i f i t e x i s t s , as the asymptotic l i m i t of the time i n t e g r a t i o n . In the f i n i t e element d i s c r e t i z a t i o n of f l u i d dynamics problems, the contrary seems to be true, that i s , i t i s computationally more e f f i c i e n t to seek the steady state s o l u t i o n only, i f the transient solutions are of no i n t e r e s t . As we are also confident, that the steady state solutions for the problems, that we intend-to solve, do e x i s t , and indeed, excellent r e s u l t s have been obtained employing the same boundary conditions, the same 18 d.o.f. t r i a n g u l a r f i n i t e elements, and the same f i n i t e element meshes by Olson [3], and Tuann and Olson [5,6], we s h a l l implement the time dependent, but also the steady state approach. The time dependent approach w i l l be used only f o r the numerical study of the flow around a c i r c u l a r c y l i n d e r , where the transient solutions are desired, as o s c i l -l a t o r y behaviour i s to be expected at the c r i t i c a l value of the Reynolds numb er. As we propose to solve the governing set of nonlinear equations i n generalized coordinates Eq. (50) simultaneously, rather than to attempt 25 a p a r t i t i o n e d s o l u t i o n , we s h a l l f i r s t cast Eq. (50) i n the form of one nonlinear time dependent matrix equation as follows: ~ t -2 y 0 0 h _ _ + Re= E E T 0 y Q y y + h -T — = 0. (51) Then the time dependent approach s h a l l be applied d i r e c t l y to the above equation, while the steady state approach s h a l l be applied to the steady equivalent of the same equation obtained by simply s e t t i n g the time dependent term equal to zero. 5.1 Time Dependent Approach We s h a l l assume l i n e a r time dependence of the vector of generalized coordinates y over the time i n t e r v a l t , which can be written as: T-t v (T) = v + — - (y + A t — y ) I K l J i t At Vit+At ltJ (52) D i f f e r e n t i a t i n g the above equation with respect to x, we obtain: v (T) = — (y + A t — y ) i ^ ; At ^t+At y-tJ (53) or evaluating at time T = t + At y (t + At) = — ( Z t + f o A t ~ y t ) (54) We note, that the above equation represents a backward f i n i t e d i f f e r e n c e scheme, which i s unconditionally stable for l i n e a r problems. This very useful feature cannot be ascertained f or nonlinear problems, however, and consequently the c r i t i c a l time step for those problems must be determined by numerical experiments. 26 I f we introduce Eq. (54) into Eq. (51), evaluate at time t+At, and add the appropriate terms, we obtain the fol lowing set of now algebraic nonlinear equations: 1 1 -T TT I + ^ E At= Re= = ^t+At i^t+At^t+At It f-Y At- t T -t+At (55) Three time integrat ion algorithms s h a l l i n turn be applied to the above equation. These algorithms d i f f e r only i n how the nonlinear term i s treated. A l l b'fsthem are quite crude, but r e l a t i v e l y cheap, and we consider them suf f i c i en t for th i s work. In the f i r s t algorithm we simply evaluate the nonlinear term for the previous time step and move i t to the r ight hand side to obtain: 1 1 -T f-rl- + —-X E At= Re= = ^t+At ^t+At At^t i^t^t T -t+At (56) The system matrix of the above set of l i nea r equations i s square, r e a l , symmetric, and nonsingular, but i t has an associated quadratic form which i s i n d e f i n i t e , so that Cholesky decomposition cannot be used. As th is matrix remains constant i t needs be inverted only once and for the con-secutive steps only the backsubsti tut ion has to be performed. In the second procedure the nonlinear term i s evaluated as 9. Y tY t +^ t» so that the fol lowing system of l inea r equations for the generalized coordinates y and h i s obtained: 2.7 E 2 t+At T -t+At (57) I t should be noted, that the nonlinear term i s treated i n the above algorithm i n the same way as an Picard i t e r a t i o n or the method of successive s u b s t i -t u t i o n . Here the i t e r a t i o n s s h a l l not be performed, because the cost would be quite high, and consequently the dynamic equilibrium i s again not exactly s a t i s f i e d . The algorithm can be very e a s i l y modified, though, to incorpor-ate the i t e r a t i v e procedure. The system matrix, which has now become un-symmetric, has to be updated at each time step, but that allows the use of much larger time step than i n the f i r s t algorithm. As a matter of f a c t , our experience has been that t h i s second algorithm i s numerically more e f f i c i e n t . [11], the nonlinear term i s treated by transposing i t to the ri g h t hand side as a d d i t i o n a l 'pseudoload' vector, s i m i l a r l y as i n the f i r s t algorithm. But then the whole r i g h t hand side including the nonlinear term i s expanded into a f i r s t order Taylor ser i e s about the previous time step. Denoting the r i g h t hand side vector as F, defined as: In the t h i r d algorithm, following a proposal by S t r i c k l i n et. a l . (58) we can wri t e : = F + At 9t - f (59) 28 By using a f i r s t order backward diffe r e n c e expression to approximate the time d e r i v a t i v e i n the above equation, i t becomes: F = T2F - F -t+1 - t £ t - l (60) We note, that the use of Eq. (59) corresponds to a l i n e a r extrapolation of the 'loads' at the two previous time steps. By introducing Eq. (60) into Eq. (55), i t can be written as: 1 1 -T At= Re= = -t+At ^t+At = F -t+At (61) The above procedure cannot be started d i r e c t l y , so we s t a r t i t by solving for the f i r s t two time steps with the f i r s t algorithm. In a l l three algorithms a test i s included on whether the steady state has been achieved. In the f i r s t algorithm a l l the entries i n the two consecutive s o l u t i o n vectors are successively scanned, and i f none of the absolute differences i n the two entries corresponding to the same coordinate exceeds a preassigned value, we assume that the steady state has been achieved. In the second algorithm the steady state has been obtained i f the test on the two consecutive determinants of the system matrix of Eq. (57) | D t + A t - D ' l < e (62) i s satisfied,where e i s a preassigned small value. In the t h i r d procedure the entries i n the two consecutive 'load' vectors are compared and the maximum absolute difference i n the two entries 29 corresponding to the same coordinate i s found. Then i f t h i s maximum difference i s smaller than a prescribed small value, we assume that the steady state s o l u t i o n has been achieved. 5.2 Steady State Approach We propose to solve the steady equivalent of Eq. (51) by Newton-Raphson method, which seeks to exactly s a t i s f y the equations of equilibrium by i t e r a t i n g u n t i l a s p e c i f i e d l e v e l of accuracy i s attained. We denote an approximate vector of generalized coordinates as y and an approximate vector of Lagrange m u l t i p l i e r s by h. Then the vector of residuals F(y,h) may be written as: F(y,h) = 1 -T Re= = y h + gZZJ -T (63) A Taylor ser i e s expansion of the vector of residuals around the p o s i t i o n (y,h) y i e l d s the following expression for the vector of residuals at an adjacent state (y + Ay, h + Ah) 3 F 3 F ? 2 F(y + Ay, h + Ah) = F(y,h) + |-Ay + -^Ah *[0[£AyVAK) . ]. (64) d-In the above equation the vector notation i s used for s i m p l i c i t y . P a r t i a l d e r i vatives may be written more rigorously as: 3F 3F V~Ay = — A y . and -rr-Ah 8y J- y^ Jj 9h -F. T - ^ A h h, k k (65) where the index i takes on successively values from 1 to w+p, the index j values from 1 to w, and the index k values from 1 to p, where w i s the number of eigenvectors used, and p i s the number of constraints i n the problem. 30 The conventional Newton-Raphson procedure retains only the terms up to the f i r s t order p a r t i a l d e r i v a t i v e s i n the Taylor expansion Eq. (64). We assume, i n addition, that the vector of residuals corresponding to the state (y + Ay, h + Ah) i s zero. These assumptions allow to rewrite Eq. (64) as: 9F dF (66) The p a r t i a l derivatives appearing i n the above equation may be obtained by d i f f e r e n t i a t i n g the steady equivalent of Eq. (51), as follows: 9F 9F • § 7 A y + 9 h ^ " Ay Ah (67) By introducing Eqs. (63) and (67) into Eq. (66) we obtain the following system of l i n e a r equations: 1 R e * + 2 2 Z = Ay n+1 "1 1 , ^T —X E Re= 1 0 _ • J Mn Qyy -T T - E y , (68) where n i s the i t e r a t i o n number. The above equation i s solved to determine the (n+1) increments i n generalized coordinates and Lagrange m u l t i p l i e r s . These increments are then used to determine an improved vector of generalized coordinates y n+^, and an improved vector of Lagrange m u l t i p l i e r s l} n +^> where: Zn+l = Y-n + A y-n+l h , T = h + Ah 1 . -n+1 -n -n+1 Equations (68) and (69) comprise the set of recurrence r e l a t i o n s needed i n the Newton-Raphson procedure. Beginning with an i n i t i a l estimate of the vectors y and h, the equations (68) and (69) are successively applied to y i e l d better and better approximation. The i n i t i a l vector of generalized coordinates i s obtained by transforming to the eigenvector basis the i n i t i a l guess on the stream function vector, which consists of a l l zero entries except f o r the con-strained ones, when the problem i s started up with Re = 1. The Jacobian of the system, that i s , the determinant of the system matrix, i s used as a test on the convergence of the procedure. At each i t e r a t i o n the Jacobian i s recorded and compared with the Jacobian of the previous i t e r a t i o n . The process i s stopped when the test | j n + 1 - J n | a _ L < e (70) i s s a t i s f i e d , where e i s the preassigned accuracy c r i t e r i o n . 5.3 F i n a l Remarks In both transient and steady state analyses, a f t e r the generalized coordinates have been computed, the corresponding nodal stream function variables are obtained v i a Eq. (35). Then the stream function subvector, con s i s t i n g of the stream function and i t s f i r s t and second d e r i v a t i v e s , can be r e a d i l y obtained at any point i n the domain. It i s done, element 31 (69) 3 2 by element, by computing the polynomial c o e f f i c i e n t s for the element under consideration using Eq. (18), then c a l c u l a t i n g the stream function subvector i n the l o c a l coordinate system v i a Eq. (14), and f i n a l l y transforming t h i s subvector back to the global f i n i t e element coordinates. The interpolated values of the stream function are needed, i n p a r t i c u l a r , to p l o t the streamlines. A rectangular g r i d i s used for contour p l o t t i n g . 33 C H A P T E R 6 P r e s s u r e F i e l d I n o u r a p p r o a c h , b a s e d o n t h e s t r e a m f u n c t i o n o n l y , t h e d i m e n s i o n l e s s p r e s s u r e f i e l d c a n b e c a l c u l a t e d w h e n a l l t h e n o d a l v a r i a b l e s i n t h e g l o b a l f i n i t e e l e m e n t c o o r d i n a t e s h a v e b e e n f o u n d . T h e p r e s s u r e i s g o v e r n e d b y t h e P o i s s o n e q u a t i o n o b t a i n e d b y t a k i n g t h e d i v e r g e n c e o f t h e m o m e n t u m e q u a t i o n s . B y e m p l o y i n g t h e c o n -t i n u i t y e q u a t i o n E q . (2) t h e s i m p l i f i e d e x p r e s s i o n f o r p r e s s u r e i s o b t a i n e d i n t h e f o r m : V 2 p = - 2 ( u v - u v ) = 2 0 ip - ) 2 = f . (71) r y x x y x x y y x y T h e e q u i v a l e n t v a r i a t i o n a l p r i n c i p l e f o r t h e a b o v e e q u a t i o n c a n b e s t a t e d a s : n = / / ( - j ( V p ) 2 + f p ) d f i - j> |2- pdr (72) a r w h e r e i s c a l c u l a t e d f r o m t h e m o m e n t u m e q u a t i o n E q . (1) w r i t t e n f o r dn t h e n d i r e c t i o n : | £ = - ^ - 0 + $ ) + V> - i») 1J1 ). (73) 9n R e s n n s s s n s s s s n T h e f u n c t i o n a l o n t h e e l e m e n t l e v e l i n t h e l o c a l c o o r d i n a t e s , ( F i g . 1), i s d e f i n e d b y t h e f o l l o w i n g e x p E g s s i p n : n 6 = / / (y(P? + P 2 ) + fp)d5dn + | E - p | d £ . (74) 34 We adopt the same f i n i t e element mesh f o r pressure as we have used f o r stream function, hence the source function f can be expressed, v i a Eq. (14), as: 2 0 2 0 f = 2 E E a a m ( m - l ) n (n - 1) - m m n n ] r s r r s s r s r s r s (75) Cm +m - 2 n +n - 2 _ r s r s The pressure gradient along the boundary edge 1-2 i s given by: \[ 0 • t- h <•{«+ w + <v« - hh^- (76) I t should be noted, that as the fu n c t i o n a l Eq. (73) contains derivatives of p up to the f i r s t order, i t would have been s u f f i c i e n t to use any element of c l a s s . Here i n order to conform with the so l u t i o n the same.truncated q u i n t i c polynomial i s used f o r i n t e r p o l a t i o n of pressure within the element, that i s : 2 0 m. n. p = Z h.K \ \ (77) i = l 1 By repeating the steps equivalent to those, indicated i n Eqs. (15), (16), (17), and (18), s u b s t i t u t i n g then Eqs. (75), (76), and (77) into the fun c t i o n a l Eq. (74), carrying out the integrations and trans-forming to global f i n i t e element coordinates, the d i s c r e t i z e d f u n c t i o n a l on the element l e v e l i s obtained i n the form: n e = i g TKP + (G + F)P (78) 35 where K i s the element matrix, G i s the load vector of 'body' load terms, and F i s the load vector of 'surface' load terms. By forming the global problem i n the usual manner the following system of equations i s obtained: KP + F = 0 (79) where K i s symmetric and p o s i t i v e d e f i n i t e matrix, and P i s the s o l u t i o n vector. We note, that the f u n c t i o n a l Eq. (72) y i e l d s the Poisson equation Eq. (71) as Euler equation, and the boundary conditions of ei t h e r I2- = or <5p = 0. (80) an dn Since the pressure gradient i s known everywhere, v i a Eq. (73), a l l the boundary conditions are of the Neumann type. That makes the so l u t i o n of Eq. (79) nonunique, however, because the system matrix K i s singular. To avoid t h i s d i f f i c u l t y we have to impose a D i r i c h l e t boundary condition at an a r b i t r a r y node. For more d e t a i l s , concerning the c a l c u l a t i o n of pressure f i e l d i t i s re f e r r e d to Tuann and Olson [5,6]. 3 6 CHAPTER 7 Numerical Implementation A l l our programs are written i n Fortran i n double p r e c i s i o n arithmetic. We employ the dynamic storage option, so that we can control the s i z e of a l l the arrays at execution time, rather than at compilation time. Thus we are able to run problems of various sizes without wasting any of the v i r t u a l memory and we can accomodate i n core a much bigger transformed nonlinear matrix Q. As we have already mentioned previously we need to compute a l l the modes of the general l i n e a r eigenvalue problem Eq. (38) i n an e f f i c i e n t manner. For t h i s purpose we use the very f a s t d i r e c t eigenvalue solver contained i n the program DRSGAL. The eigenvalue problem i s solved as follows. F i r s t l y , as the global consistent mass matrix i s p o s i t i v e d e f i n i t e , i t s inverse can be found by LU decomposition. The global l i n e a r d i s s i p a t i v e matrix K i s premultiplied by t h i s inverse to transform the general eigenvalue problem to the simpler form K lp = A lp (81) The symmetric matrix K, of the order N, i s reduced to a symmetric t r i d i a g o n a l matrix, a f t e r N-2 orthogonal s i m i l a r i t y transformations, using the Householder method. The eigenvalues and the eigenvectors of the t r i d i a g o n a l matrix are found by QL transformations, and transformed back to the eigenvectors of K. -5 3 The CPU time, i n seconds, for the eigenvalue a n a l y s i s , i s about 1.-8 x 10 x N , where N i s the s i z e of 'net' problem. Then we solve the l i n e a r steady equivalent of Eq. (51) i n order to determine which modes give a reasonable representation of the solu t i o n of the l i n e a r 37 problem with the convective terms taken equal to zero. This i s conveniently done by a d i r e c t solver which i s based on Gaussian elimination. The CPU —6 3 time f o r t h i s procedure, i n seconds, i s approximately 5.2 x 10 (NT) , where NT i s equal to the number of the eigenvectors used plus the number of constraints of the problem. We note, that i t i s p a r t i c u l a r l y important to perform t h i s procedure i n the square cavity flow problem, where, as i t turns out, some higher modes must be included i n order to approximate even the so l u t i o n of the l i n e a r problem. The flow around a c i r c u l a r c y l i n d e r , on the other hand, due to the l e s s stringent boundary conditions, behaves much l i k e the problems encountered i n s t r u c t u r a l dynamics, that i s , the higher modes need not be included. A f t e r having established which modes we are going to use as the new basis, we proceed to set up the complete nonlinear problem i n the f i n i t e element coordinates. ThewCPU time f o r s e t t i n g up the matrices i s about 2.6 seconds per each new group of elements, where a group consists of elements, which have the same dimensions and the same o r i e n t a t i o n i n space. We note, that for the square cavity flow and the flow around a c i r c u l a r cylinder problems we choose the o r i g i n of the global coordinate system, so that there e x i s t s one axis of symmetry. Consequently, the modes are eit h e r symmetric or antisymmetric, which r e s u l t s i n some s p e c i a l properties of the transformed nonlinear convective matrix Q. = Q.., E.,E. E, . Namely, i f a l l three of lmn l j k i i jm kn J the modes 1, m, and n are symmetric, or i f any two of them are antisymmetric and the t h i r d i s symmetric, the corresponding entry i n the transformed matrix Q i s zero, otherwise i t i s nonzero. Consequently, only the m u l t i p l i c a t i o n s i n the transformation, which r e s u l t i n the nonzero e n t r i e s , are performed increasing the e f f i c i e n c y of the transformation procedure. The CPU time, i n seconds, for the transformation, can be estimated by the formula 38 7.0 x 10~ 6 x N x [NE x LB3 x (LB3-1) + NE 2 x LB3 + NE 2 x (NE+1) / 2] , where N i s the 'net' s i z e of the problem, NE i s the number of the modes used i n the transformation, and LB3 = 2 x LB - 1 i s the bandwidth of the problem i n the f i n i t e element coordinates. We note, that the dynamic storage option increased the number of modes, that can be used i n the anal y s i s , to 62. This l i m i t of maximum 62 modes i s imposed by the capacity of the v i r t u a l memory of the IBM 370/168 machine at UBC, that i s , the maximum order of the transformed nonlinear matrix Q that can be kept i n core i f double p r e c i s i o n arithmetic i s used, even taking account of the symmetry i n the f i r s t two in d i c e s , i s 62. Otherwise, a u x i l i a r y storage locations have to be employed. Haying stored the global consistent mass matrix-* the eigenvalues and the corresponding modes, and the transformed nonlinear matrix Q on tape, we have a l l the arrays indicated i n Eq. (50) and can run both steady state and time dependent analyses. In both cases, due to the Lagrange m u l t i p l i e r s technique, the system matrix i s nonpositive d e f i n i t e , although i t i s non-singular. Accordingly, the Gaussian elimination with p a r t i a l p i v o t i n g , and forward and backward s u b s t i t u t i o n i s used to solve the r e s u l t i n g set of l i n e a r i z e d algebraic equations. The CPU time f o r one i t e r a t i o n of the —6 3 Newton-Raphson process, i n seconds, i s about 4.0 x 10 x (NT) . For the f i r s t time dependent algorithm, with the nonlinear term transposed to the ri g h t hand side, and evaluated at the preceding time step, the CPU time per —6 3 time step, i n seconds, i s approximately 2.4 x 10 x (NT) . The CPU time for the second time dependent algorithm, where the system matrix i s updated —6 3 at every time step, i n seconds, i s about 3.4 x 10 x (NT) . As before, NT denotes the number of eigenvectors plus the number of constraints. 2 The scalar v o r t i c i t y C, = V ^ i s obtained by the i n t e r p o l a t i o n program, 39 which i s used f o r p l o t t i n g . A f t e r the stream function subvector has been obtained at a l l points of a s p e c i f i e d regular i n t e r p o l a t i o n mesh, the v o r t i c i t y at each of those points i s simply equal to the sum of the second derivatives of the stream function il; and ib . The mesh i s , i n general, xx yy nonuniform with v a r i a b l e spacings i n x- and y- d i r e c t i o n s . The whole domain can be covered by the mesh, or j u s t some regions of i n t e r e s t . The stream-l i n e s , the pressure f i e l d , and the v o r t i c i t i e s are then plotted using the standard contour subroutines. We note, that while the stream function and the v e l o c i t i e s are continuous, the v o r t i c i t y i s only piecewise continuous, which accounts f o r not so good v o r t i c i t y p l o t s , e s p e c i a l l y for crude f i n i t e element g r i d s . The pressure f i e l d representation i s excellent, however, because the method f o r c a l c u l a t i n g pressures has a tendency of 'smearing out'. 40 CHAPTER 8 Example Applications The foregoing modal f i n i t e element method has been used to solve several example flow problems. Our c r i t e r i o n i n choosing example problems has been a v a i l a b i l i t y of known r e s u l t s , so that ready checks would be provided. Numerical r e s u l t s f o r f u l l y developed plane P o i s e u i l l e flow, c i r c u l a t o r y flow i n a square cavity and the flow around a c i r c u l a r c y l i n d e r problems, obtained by modal f i n i t e element method, are presented i n t h i s chapter. The re s u l t s f o r P o i s e u i l l e flow problem are compared with the exact closed form s o l u t i o n . The steady state r e s u l t s f o r c i r c u l a t o r y flow i n a square cavity and the flow around a c i r c u l a r cylinder problems are compared to the r e s u l t s obtained by d i r e c t f i n i t e element approach using the same f i n i t e element grid and to the r e s u l t s obtained by various f i n i t e d i f f e r e n c e techniques, considered exact herein, as no closed form solutions e x i s t for these problems. 8.1 F u l l y Developed Plane P o i s e u i l l e Flow The f i r s t example chosen was that of a f u l l y developed flow between p a r a l l e l w alls. The exact sol u t i o n for t h i s problem shows the flow laminar and d i s t r i b u t e d p a r a b o l i c a l l y between the walls, with a corresponding l i n e a r pressure f i e l d decreasing downstream. The f i n i t e element g r i d used for t h i s problem i s shown i n F i g . (3). The following boundary conditions were used, a l l of them kinematic i n the sense of the r e s t r i c t e d v a r i a t i o n a l p r i n c i p l e : 2 3 (a) on the stream function ; ip = ip(y) = 3y - 2y on the upstream section, ip = 1 on the upper w a l l , and ip = 0 on the lower w a l l ; (b) on i t s normal d e r i v a t i v e ; i p x = - v = 0 on the upstream and the downstream sections, and ip = u = 0 on therupper and the lower w a l l . On the downstream section the 41 natural boundary conditions of zero pressure gradient across the flow P y = 0 were to be approximated. The 'net' s i z e of the problem was 12 with 4 constraints ; i p = 1 at nodes 2 and 5, i p = 6 at node 1, and i p = -6 at node 2. We note, yy yy that because the f i n i t e element used contains a complete cubic for v e l o c i t y i n t e r p o l a t i o n , i t was capable of exactly representing the parabolic v e l o c i t y p r o f i l e and the downstream natural boundary conditions of constant pressure. The 12 mode shapes from the l i n e a r eigenvalue analysis are shown i n F i g . (4), along with the corresponding eigenvalues f o r Re = 1. The curves represent equal steps i n stream function values. These modes were then used to represent the nonlinear equations indicated i n Eq. (50). Both steady state and transient analyses were then performed using 5,6,7,8,9,10,11, and 12 modes, re s p e c t i v e l y . We note, that the lowest number 5 a c t u a l l y gave only one free mode, because of the 4 constraints imposed on the problem. The c a l c u l a t i o n s were started with a l l free nodal variables equal to zero i n both the analyses. The Newton-Raphson steady state procedure converged i n 3 or les s i t e r a t i o n s i n a l l cases. This simple problem was also used to check the e f f i c i e n c y of the three time i n t e g r a t i o n algorithms. In the f i r s t and the t h i r d algorithms we used the time step of PAt = 2/X , where X was the eigenvalue associated max max with the highest mode kept i n the a n a l y s i s . In the second algorithm, due to i t s greater numerical s t a b i l i t y , we could use a bigger time step, so we a r b i t r a r i l y chose At-= 0.15 seconds. The f i r s t algorithm, with the nonlinear term transposed to the r i g h t hand side and evaluated at the previous time step, converged to the steady state i n maximum 67 increments, when a l l 12 modes were used. We assumed that the steady state was achieved, when the maximum absolute diff e r e n c e i n 42 the two e n t r i e s corresponding to the same coordinate did not exceed e = 10 In the second algorithm the steady state s o l u t i o n was attained i n 7 increments, with e of Eq. (62) prescribed as 10 ^. When the exact s o l u t i o n was used as the i n i t i a l guess the t h i r d algorithm was able to reproduce t h i s r e s u l t i n 3 steps. When a l l the free nodal variables were zeroed, however, the algorithm did not converge. On the basis of these c a l c u l a t i o n s we decided tosemploy only the f i r s t two algorithms i n t h i s work, although the t h i r d one might have converged, as w e l l , had the i n i t i a l conditions been closer to the true s o l u t i o n . The r e s u l t i n g streamlines and the predicted maximum v e l o c i t i e s i n the d i r e c t i o n of the flow at the midnode 3 f o r 5,6,7,9, and 12 modes, res p e c t i v e l y , are shown i n F i g . (5). We see, that although the 5 mode r e s u l t i s rather poor, the 6 mode one i s already quite acceptable, and compares very well with the e x a c t r r e s u l t . As a matter of f a c t , as long as at l e ast f i r s t 9 modes are used the exact r e s u l t i s obtained and only the generalized coordinates associated with modes 1,6, and 9 are f i n i t e , so that j u s t these modes contribute to the s o l u t i o n . The exact s o l u t i o n i s namely antisymmetric i n stream function with respect to the axis z of F i g . (3), so that none of the symmetric modes can a f f e c t i t . Taking that into account we reproduced the exact s o l u t i o n with j u s t 5 modes, but we had to include modes 1,6, and 9, whichppossessathe antisymmetric properties, while the remaining 2 modes were a r b i t r a r y . The generalized coordinate associated with mode 1 i s always predominant. A l l r e s u l t s are independent of the Reynolds number. 8.2 C i r c u l a t o r y Flow i n a Square Cavity As the second problem we chose the flow within a square cavity which i s bounded by three fixed walls and an upper l i d moving with constant 43 v e l o c i t y i n i t s own plane, shown i n Fig.(6). Mathematical s i n g u l a r i t i e s are present i n the problem i n the regions of the upper v e r t i c e s , 2 and 4, where a f i x e d wall meets a moving l i d , due to the flow having to move without s l i p at the speed of the moving l i d and yet be zero at the f i x e d w a l l . To avoid t h i s d i f f i c u l t y two extra nodes were introduced close to the upper v e r t i c e s , thus allowing for t r a n s i t i o n from zero normal v e l o c i t y at the f i x e d walls to the prescribed tangential v e l o c i t y of the moving l i d at these new nodes. Two a d d i t i o n a l nodes were needed i n the neighbourhood of the upper v e r t i c e s on the fixed walls, as w e l l , because of e x i s t i n g asymmetric pressure s i n g u l a r i t i e s . A l l boundary conditions imposed on the problem were ' r i g i d ' i n the sense of the r e s t r i c t e d v a r i a t i o n a l p r i n c i p l e . Accordingly, the only free nodal v a r i a b l e along the f i x e d walls 1-2 and 3-4 was ip , along the bottom f i x e d w a l l 1-3 ' P y y * while on the upper l i d if^ and T|I were the only free nodal v a r i a b l e s , so that the three fixed walls and the moving l i d a c t u a l l y were zero streamlines. The s i z e of the cavity was assumed to be unity, and the given v e l o c i t y of the moving l i d was also taken to be unity i n the d i r e c t i o n to the l e f t . Thus, the Reynolds number was n a t u r a l l y defined as Re = ^. The three f i n i t e element grids used are shown i n F i g . (7) together with the number of elements (NE), the number of nodes (NO), the s i z e of the 'net' problem (NN), and the h a l f bandwidth i n the f i n i t e element coordinates (LB). The number of constraints for a l l three grids was equal to 1, because by using the 'master-slave' option we forced a l l the tangential v e l o c i t y degrees of freedom ip^ on the upper l i d to be the same, thereby leaving only free v a r i a b l e s . A l l 15 mode shapes with associated eigenvalues for SQCA 12-13 g r i d 44 and the f i r s t 30 modes for SQCA 36-29 g r i d are shown i n F i g . (8), arranged i n the order of the ascending eigenvalues. We note, that the shape of the highest modes of the SQCA 12-13 g r i d i s already influenced by the f i n i t e element layout. The most noteworthy feature, v a l i d f o r a l l g r i d s , i s the strk i n g increase i n complexity of the modes with only one order of magnitude increase i n the corresponding eigenvalues. The lowest and the highest eigenvalue f o r a l l three g r i d s , together with the number of symmetric and antisymmetric modes and the t o t a l number of modes are given i n Table ( I ) . The i n t e r v a l bounded by the lowest and the highest eigenvalue of the f i n e r grids also includes a l l the eigenvalues of the cruder g r i d s , which serves as a check on the eigenvalue sol v e r . The lowest eigenvalue only changes s l i g h t l y from g r i d to g r i d , which seems to suggest that i t i s predicted accurately enough even f o r the very crude SQCA 12-13 g r i d . The corresponding modes are also almost i d e n t i c a l . The modes are subdivided into the symmetric and the antisymmetric ones with respect to the v e r t i c a l y axis. Only the symmetric modes are d i r e c t l y loaded by the v e l o c i t y y constraint on the moving l i d , while the antisymmetric modes are excited only through nonlinear coupling. It i s then obvious that the antisymmetric modes do not contribute to the l i n e a r -4 sol u t i o n f o r Re = 10 , while without these modes no nonlinear behaviour can be simulated. When the modal approach was t r i e d on t h i s problem, we found that reasonable r e s u l t s could only be obtained i f some of the higher modes were included i n the an a l y s i s . This was true even f o r the l i n e a r problem with Re = 0, so we included a l l the symmetric modes, whose generalized coordinates were greater than some prescribed small value e for Re = 0 decomposition. The antisymmetric modes used corresponded to the lowest eigenvalues. 45 On the SQCA 12-13 g r i d we used 15 and 11 modes, on the SQCA 36-29 gr i d 36 and 42 modes, and on the SQCA 26-53 g r i d 60 modes, r e s p e c t i v e l y . Only the Newton-Raphson steady state procedure was run. The c a l c u l a t i o n s -4 were started with the Reynolds number Re = 10 using the n u l l s o l u t i o n , with a l l free v a r i a b l e s equal to zero, as anxini,tia ,l guess i n the f i n i t e element coordinates. The transformation to the eigenvector basis using the orthonormal properties of the eigenvectors had obviously to be performed before the i t e r a t i o n s were started to y i e l d an^nitti-aU.guess i n the generalized coordinates. A f t e r the s o l u t i o n i n the eigenvector basis had been obtained the transformation back to the f i n i t e element coordinates was performed v i a Eq. (35) to y i e l d the stream function s o l u t i o n . This s o l u t i o n was then used as an i n i t i a l guess for the next higher Reynolds number and the process -4 repeated. These steps were ca r r i e d out at RS = 10 , 1, 10, 20, 40, 100, 200 and 400 for a l l three g r i d s . For SQCA 76-53 g r i d the rangee o>f Reynolds numbers was extended to Re = 3000 with the a d d i t i o n a l steps c a r r i e d out at Re = 600, 1000, 1400, 1700, 2000, 2200, 2400, 2800 and 3000. The accuracy —6 test on the Jacobian e of Eq. (70) was set equal to 10 . Regardless of the gr i d and the number of modes used, that i s the problem s i z e , convergent solutions for the whole range of Reynolds numbers under consideration were obtained i n 3 to 6 i'iterations. We note that t h i s held true even f or Reynolds numbers up to 3000, so i t i s concluded that numerical s t a b i l i t y of modal method i s quite high. When the Poisson equation, Eq. (79), was solved to obtain the pressure f i e l d , the D i r i c h l e t boundary condition p = 0 was imposed at the node located i n the middle of the bottom wa l l . This e f f e c t i v e l y means that the pressure,distribution i s referenced to the pressure at the middle of the bottom w a l l . A l l the other pressure nodal variables were l e f t free, 46 and accordingly, for a l l t r iangular f i n i t e elements with two of the ver t ices located on.any part of the boundary, the boundary in tegra l of Eq, (74) had to be included i n the element funct ional expression to conform with the v a r i a t i o n a l p r i n c i p l e . Error i n stream function at midnode 5 of the SQCA 12-13 g r id for various combinations of modes i s l i s t e d i n Table ( I I ) . We note that jus t one antisymmetric mode i s enough to produce resul ts with maximum 1% error for Reynolds numbers up to about 40, when compared to the d i rec t f i n i t e element so lu t ion on the same g r i d . Complete resu l t s for Reynolds numbers. !Re = 0, 10, 20, 40, 100, 200 and 400 are plot ted i n F igs . (9) - (11). The t o t a l number of modes and the number of symmetric and antisymmetric modes, together with predicted coordinates of the vortex centre and predicted values of stream function v o r t i c i t y and pressure at the vortex centre are l i s t e d i n Table ( I I I ) . These resul ts are compared with the d i rec t f i n i t e element resul ts on the same grids and with the Burg:graf's f i n i t e differences resul ts [12]. Some small d i s -crepancies between the d i rec t f i n i t e element results reported here and the ones of [5] are due to the different in te rpo la t ion grids used. Namely, we used crude 20x20 in te rpola t ion gr id i n this work, while [5] employed more accurate i r r egu la r in te rpo la t ion g r i d . Only the general shape of streamlines from modal method and d i rec t f i n i t e element approach can be compared with Burggraf's r e s u l t s , because contour l eve l s do not match. The v o r t i c i t i e s and the pressures can be compared d i r e c t l y , however. The SQCA 12-13 gr id i s too crude to produce accurate resul ts even for low Reynolds numbers and i f the d i r ec t f i n i t e element approach i s used. This i s espec ia l ly true for the v o r t i c i t y predic t ions , but a l so , for the stream function and the pressure predic t ions , which are i n general too high. 47 The 11 modes r e s u l t s with j u s t 1 antisymmetric mode do agree quite w e l l with the 15 modes r e s u l t s and the d i r e c t f i n i t e element approach r e s u l t s on the same g r i d , which are of course i d e n t i c a l , up to Re of about 40. For higher Re the agreement r a p i d l y d e t e r i o r a t e s . In the 15 modes r e s u l t s the change of p o s i t i o n of the vortex centre with increasing Re i s followed to some extent,while 11 modes r e s u l t s show the vortex centre stationary f o r the whole range of Re under consideration. In the 60 modes r e s u l t s on the SQCA 76-53 g r i d the lower r i g h t secondary vortex, which appears i n the exact s o l u t i o n , does show at Re = 40, then increases gradually at Re = 100 and 200, but disappears at Re = 400. For a l l these Reynolds numbers the seconary vortex i s s h i f t e d a l i t t l e to the l e f t . The stronger lower l e f t vortex of the exact s o l u t i o n does not show at a l l . The v o r t i c i t y contours are somewhat wavy i n appearance, e s p e c i a l l y i n the region close the bottom w a l l , but the general trend seems to be preserved. The pressure contours agree xjuite w e l l with the exact r e s u l t s except for the region close to the bottom w a l l . The position of the vortex centre i s predicted very w e l l for the whole range of Reynolds numbers under consideration. The same i s true for the stream function and the pressure values at' the vortex centre. The v o r t i c i t y value at the vortex centre agrees well with the exact r e s u l t f or Re up to about 20, while for higher Re the predicted v o r t i c i t y values are much too high. B Best modal r e s u l t s are achieved with 36 modes and 42 modes on the SQCA 36-29 g r i d . Both the 36 modes and the 42 modes r e s u l t s appear to reproduce the streamlines.predicted by the d i r e c t f i n i t e element approach on the same g r i d very w e l l for the whole range of Reynolds numbers considered, e s p e c i a l l y near the upper l i d where the gradients are high. They also compare reasonably well with Burggr.af's r e s u l t s , considered exact herein, reproducing the lower r i g h t 48 s e c o n d a r y v o r t e x , b u t n o t t h e s t r o n g e r l o w e r l e f t s e c o n d a r y v o r t e x . The v o r t i c i t y p r e d i c t i o n s a l s o compare s a t i s f a c t o r i l y w i t h t h e d i r e c t f i n i t e e l ement a p p r o a c h r e s u l t s and even B u r g g f f f ' s r e s u l t s , c o n s i d e r e d e x a c t a g a i n , f o r t h e w h o l e r ange o f R e y n o l d s numbers unde r c o n s i d e r a t i o n . The 42 modes r e s u l t s a r e s l i g h t l y b e t t e r t h a n t h e 36 modes ones as c o u l d have been e x p e c t e d f o r b o t h t h e s t r e a m f u n c t i o n and t h e v o r t i c i t y . The p r e s s u r e r e s u l t s a r e i n e x c e l l e n t agreement w i t h t h e r e s u l t s o b t a i n e d f o r t h e f u l l SQCA 36-29 g r i d and even w i t h t h e f u l l SQCA 76-53 g r i d r e s u l t s , c o n s i d e r e d e x a c t h e r e i n , up t o Re o f about 100 . F o r h i g h e r Re t h e agreement i s s t i l l f a i r , e s p e c i a l l y be tween moda l r e s u l t s and t h e d i r e c t f i n i t e e l emen t a p p r o a c h on t h e same g r i d w i t h t h e 42 modes r e s u l t s a g a i n h a v i n g a s l i g h t edge o v e r 36 modes r e s u l t s . F rom a l l t h e s e r e s u l t s we c o n c l u d e t h a t good agreement w i t h t h e fexact' s o l u t i o n o f t h e d i s c r e t i z e d p r o b l e m , o b t a i n e d by d i r e c t f i n i t e e l ement a p p r o a c h , f o r t h i s r a n ge o f R e y n o l d s number s , c an be a c h i e v e d by e m p l o y i n g abou t 50% o f t h e modes, t h a t i s , t h e number o f modes, w h i c h a r e r e q u i r e d , goes up w i t h t h e r e f i n e m e n t o f t h e f i n i t e e l ement mesh. Tha t e x p l a i n s why 36 modes and 42 modes on t h e c r u d e r SQCA 36-29 g r i d y i e l d b e t t e r r e s u l t s t h a n 60 modes on t h e f i n e r SQCA 76-53 g r i d , w h i c h c o n t a i n s much more d e g r e e s o f f r e e d o m . 8.3 F l o w A round a C i r c u l a r C y l i n d e r The f i n a l example c o n s i d e r e d was t h e c l a s s i c a l p r o b l e m o f t h e f l o w a round a c i r c u l a r c y l i n d e r . P h y s i c a l l y t h e p r o b l e m i n v o l v e s an i n f i n i t e l y l o n g c y l i n d e r immersed i n a f l u i d medium o f i n f i n i t e e x t e n t . We a d o p t e d f i n i t e d o m a i n s , h oweve r , as i t i s u s u a l l y done i n c o m p u t a t i o n a l f l u i d dynamic s i n v o l v i n g e x t e r i o r f l o w . 49 The f i r s t domain was used for numerical simulation of the steady flow, while the second domain, obtained by extending the f i r s t domain for si x u nits downstream, was used for the transient flow simulation. Both domains are shown i n F i g . (12), together with the adopted f i n i t e element g r i d s , the number of elements (NE), the number of nodes (NO) and the c e n t r a l angle (AQ). The centre of the cylinder i s located at the o r i g i n of the (x,y) plane for both domains. The boundary i s divided into the inflow section I\, the outflow section T Q , the c y l i n d e r wall and the top and the bottom sections p a r a l l e l to the incident uniform stream flow of un i t v e l o c i t y along the p o s i t i v e x-axis. The boundary conditions at the inflow section r were a l l ' r i g i d ' , ft obtained by spec i f y i n g the uniform stream flow defined by u=l and v=0. The stream function nodal v a r i a b l e s , as the d i r e c t r e s u l t of these ' r i g i d ' boundary conditions, were given as ip=y, ip =0, il> =1, \b =0, and ip =0, while was l e f t J x ' r y ' Txy yy xx fre e . On the outflow section T q the 'rigid'boundary condition v=o was used to pair up with the 'natural' boundary condition of constant pressure Py=0 i n the steady case. These conditions allowed the u=- v e l o c i t y to develop on t h i s section. The v= -ty^® ' r i g i d ' boundary condition implied that i»^=Q, which i n turn indicated that ip^ was to be unknown. The 'natural' boundary condition jp^ =0 required that ib and lp were to be unknown. Accordingly, on thi s section the boundary conditions i n terms of nodal stream function v a r i a b l e s were given as \b =0 and \b =0 with lp, \b , ip and \b l e f t f r e e . The ' r i g i d ' x xy y xx yy boundary conditions of n o - s l i p i>n=0 and <P.g=0 were d i r e c t l y imposed at each of the nodes located- on the c y l i n d e r w a l l T ^ . The <Pg=0 condition implied that ib =0 and ft =0 on r , and also that ip was a constant, which was set equal to rns r s s w r zero. Accordingly, at any node on of the s i x nodal v a r i a b l e s i n the l o c a l (s,n) system, which was used i n place of the global (x,y) system, only was retained while a l l others were zeroed and eliminated. This had the e f f e c t 50 of s a t i s f y i n g the n o - s l i p conditions exactly only at d i s c r e t e nodes on the cylinder w a l l , unlike the s t r a i g h t boundaries, where the prescribed boundary conditions can be r e a l i z e d pointwise. On.the top and the bottom sections the ' r i g i d ' boundary condition of u=l was set to pair up with the 'natural' boundary condition of constant pressure p =0 for the steady state a n a l y s i s . These conditions allowed the v- component to develop along T^. The u=iji =1 ' r i g i d ' boundary condition implied that ^ X y = 0 which required i p x to be unknown. The 'natural' boundary condition pv=0 indicated that i l l and iii had to be y xx unknown. Consequently, the boundary conditions i n terms of nodal stream function v a r i a b l e s were given as \b =1 and ip =0 with i l> , i i i , iji and i i i l e f t ° y xy x xx yy fr e e . A l t e r n a t i v e l y , f o r the steady flow the assumption can be made that the flow i s symmetric with respect to the x-axis and so only flow i n the upper ha l f domain y > 0 need be considered. The symmetryoof the flow implies that u i s an even function of y and v i s an odd function of y, so that both v and v are odd and vanish on the x-axis. It follows that both i i i and £ are y x odd and therefore also vanish on y=0. So the ' r i g i d ' boundary condition ip=0 i s imposed on the x-axis, which has become the part of the boundary, to pair up with the 'natural' boundary condition of zero shear stress or ?=0. The 'natural' boundary condition requires that iii =0 and 1J1 =0, the l a t t e r of xx yy which i n turn implies that \\i i s unknown. Hence i i i - i s also unknown, and the boundary conditions i n terms of nodal stream function v a r i a b l e s are given as il>=0, <JJ =0, i l l =0 and i i i =0, with i i i and i i i l e f t f r e e . We note that the x xx yy y xy 'natural' boundary condition £=0 holds only at the nodes.on x-axis. The transient analysis requires the f u l l domain to be represented, since the symmetry has vanished. The same boundary conditions were imposed on the inflow boundary r and the cylinder w a l l as i n the steady case. 51 The top and the bottom sections T were assumed to be zero f r i c t i o n smooth r u walls and so the ' r i g i d ' boundary conditions ip=y and v=ip^=l were imposed there. The nodal stream function v a r i a b l e s on T , as the d i r e c t r e s u l t u of these ' r i g i d ' boundary conditions, were given as ^=y, =0, ^ =1, tp =0 and ii; =0, while \li was l e f t f r e e . On the outflow s e c t i o n T we needed Txy r y y o the flow to be as unconstrained as possi b l e , so that the ' r i g i d ' boundary condition v=0 was d e f i n i t e l y out of question. The best that we could do was to leave a l l nodal stream function v a r i a b l e s on T f r e e , thus i n l i e u of the governing r e s t r i c t e d v a r i a t i o n a l p r i n c i p l e e f f e c t i v e l y s pecifying the 'natural' boundary conditions of zero shear stress T=0 and constant pressure Py =0. We note, that even these less stringent boundary conditions are erroneous. That i s why we extended the computational domain for s i x units downstream thus hoping to reduce the influence of the outflow boundary conditions on the region of s p e c i a l i n t e r e s t behind the obstacle, where unsteady behaviour was expected at a c e r t a i n value of Reynolds number. When the eigenvalue problem was solved the i|> along the boundaries I\ and were constrained to be equal by using the 'master-slave' option. For the GYLFL 92-63 g r i d the stream function ip was also constrained to be equal along the boundary using the same option. The number of unknowns for the CYLFL 84-58 g r i d was 223 and for the CYLFL 92-63 g r i d 226. The f i r s t 15 mode shapes, antisymmetric i n stream function, with the associated eigen-values, are shown i n F i g . (13-1) for the CYLFL 84-58 g r i d . The f i r s t 20 mode shapes, antisymmetric i n stream function with the associated eigenvalues, and the f i r s t 10 mode shapes, symmetric i n stream function, with the associated eigenvalues, are shown i n F i g . (13-2). The contour l e v e l s represent equal steps i n stream function and are plotted for the upper ha l f domain y > 0 to save space. The net f o r p l o t t i n g was 26 x 48 for the CYLFL 84-58 g r i d and 52 26 x 54 f o r t h e CYLFL 92-63 g r i d , w i t h x - and y - s p a c i n g s v a r y i n g f r o m 0.4 t o 1 .0 , i n c r e a s i n g o u t w a r d l y f r o m t h e c e n t r e o f t h e c y l i n d e r , f o r b o t h g r i d s . We o b s e r v e t h a t t h e e i g e n v a l u e s i n c r e a s e q u i t e r a p i d l y w i t h mode number compared t o t h e s q u a r e c a v i t y f l o w p r o b l e m . When moda l a p p r o a c h was t r i e d on t h e CYLFL 84 -58 g r i d f o r t h e s t e a d y f l o w we o b s e r v e d t h a t no h i g h e r modes we re needed i n o r d e r t o a t t a i n r e a s o n a b l e r e s u l t s . We a l s o o b s e r v e d t h a t t h e g e n e r a l i z e d c o o r d i n a t e s a s s o c i a t e d w i t h t h e s y m m e t r i c modes were a l l p r a c t i c a l l y z e r o f o r t h e w h o l e r ange o f R e y n o l d s numbers c o n s i d e r e d , c o n f i r m i n g t h a t , a s e x p e c t e d , t h e s y m m e t r i c modes c o u l d no t c o n t r i b u t e t o t h e s o l u t i o n w h i c h was a n t i s y m m e t r i c i n s t r e a m f u n c t i o n . H e n c e , i n a l l s ub sequen t c a l c u l a t i o n s r e p o r t e d h e r e i n , we u sed t h e uppe r h a l f doma in y > 0 o f F i g . ( 13 -1 ) o n l y , a p p l y i n g on t h e bounda ry d e f i n e d by t h e x - a x i s t h e bounda r y c o n d i t i o n s as c o v e r e d p r e v i o u s l y . The d i s c r e t e p r o b l e m , t h u s o b t a i n e d , had 104 d e g r e e s o f f r eedom i n t h e f i n i t e e l e m e n t c o o r d i n a t e s , w h i l e t h e b a n d w i d t h was 48 . The number o f e l e m e n t s was 42 and t h e number o f nodes 34. 62 modes were u sed i n t h e m o d a l r e p r e s e n t a t i o n , a l l o f them e v i d e n t l y a n t i s y m m e t r i c i n s t r e a m f u n c t i o n . The number o f c o n s t r a i n t s i n t h e p r o b l e m was 3 , ip=3.0 and #=20.0 a t t h e two nodes on t h e I\ bounda ry o b t a i n e d f r o m t h e bounda ry c o n d i t i o n ifj=y t h e r e , and i|> =1 a t a l l nodes l o c a t e d on t h e T. and r b o u n d a r i e s . The c a l c u l a t i o n s were s t a r t e d w i t h Re = 1 u s i n g x u t h e n u l l s o l u t i o n w i t h a l l v a r i a b l e s , e x c e p t f o r t h e c o n s t r a i n e d o n e s , z e r o e d a s an i n i t i a l g ue s s i n t h e f i n i t e e l ement c o o r d i n a t e s . T h i s i n i t i a l g ue s s v e c t o r was t h e n t r a n s f o r m e d t o t h e e i g e n v e c t o r b a s i s t o y i e l d an i n i t i a l gues s v e c t o r i n t h e g e n e r a l i z e d c o o r d i n a t e s . A f t e r t h e s o l u t i o n i n t e rms of g e n e r a l i z e d c o o r d i n a t e s had been o b t a i n e d by Newton-Raphson i t e r a t i v e p r o c e d u r e , i t was t r a n s f o r m e d b a c k t o t h e f i n i t e e l ement c o o r d i n a t e s v i a E g . (35) t o y i e l d t h e s t r e a m f u n c t i o n s o l u t i o n . T h i s s o l u t i o n was t h e n u sed 53 as an i n i t i a l guess f o r next higher Re, and the process repeated. These steps were c a r r i e d out at Re = 1,5,7,10,20,40,70 and 100. Maximum 5 i t e r a t i o n s were needed to a t t a i n an accuracy of e=10 ^ i n the whole range of Reynolds numbers considered. No o s c i l l a t o r y solutions or i n s t a b i l i t i e s were encountered. The pressure f i e l d was computed only for the steady flow. The upper hal f domain y > 0 only was considered. When the Poisson equation (79), was solved for pressure, a l l the boundary conditions were of the Neumann type, as already mentioned previously i n Chapter (6). Tuann and Olson [6] found, however, that the predicted pressure on the f a r boundaries P., T and r r x u o was very nearly zero, as i t should have been according to boundary conditions imposed on the stream function IJJ. As a l l our boundary conditions on the stream function lb are the same as those of [6], for the steady flow, we forced p to be i d e n t i c a l l y zero on these boundaries. The symmetry condition Py=0 was also enforced a l l along the x-axis. The dimensionless t o t a l drag c o e f f i c i e n t C^ was obtained as the sum of the f r i c t i o n drag c o e f f i c i e n t C^ calculated as C,. = % k s i n 6d6, and the pressure drag c o e f f i c i e n t f Re o nn ' n=0 0^ , given by the expression C^ = -/ p cost,0d6. These l i n e i n t e g r a l s were only approximated, however, because the integrations were not performed along the c y l i n d e r , but along the polygonal segments of the f i n i t e element g r i d approximating the c y l i n d e r . The number of unknowns i n the d i s c r e t e problem was 155, the bandwidth 50, while the number of elements and the number of nodes were 42 and 34, r e s p e c t i v e l y , that i s the same as for the stream function c a l c u l a t i o n s . Complete stream function and v o r t i c i t y r e s u l t s f o r steady flow f o r Reynolds numbers Re = 1,5,7,10,20,40,70 and 100 are plotted i n F i g s . (14) and (15), r e s p e c t i v e l y . Only the region of i n t e r e s t i s shown extending from 3 units upstream to 9 units downstream i n the x d i r e c t i o n measured from the centre of the c y l i n d e r , and for 3 units i n the d i r e c t i o n of the p o s i t i v e 54 y-axis measured from the x-axis. The i n t e r p o l a t i o n net was 25 x 66 with the x-spacings^varying from 0.1 to 0.4, and the y—spacings from 0.1 to 0.2, both increasing outwardly from the centre of the c y l i n d e r . The net functions, ib and t,, at the i n t e r i o r points of the cylinder were set to zero. The contour l e v e l s were s p e c i f i e d to match those reported by Tuann and Olson [6], so that d i r e c t comparisoncwould be p o s s i b l e . These modal r e s u l t s are compared with the d i r e c t f i n i t e element r e s u l t s on the same g r i d [6] and the f i n i t e d i f f e r e n c e r e s u l t s by Dennis and Chang [13] or Takami and K e l l e r [14,15]. The agreement between modal stream function r e s u l t s and stream function r e s u l t s by the d i r e c t f i n i t e element method on the same g r i d i s excellent up to Re of about 40, although the f i r s t appearance of a negative valued stream function becomes v i s i b l e at Re = 10, whereas i t i s observable at Re = 7 i n the d i r e c t f i n i t e element r e s u l t s . At Re = 70 the zero streamline does not extend f a r enough downstream and t h i s trend gets even more pronounced at Re = 100 r e s u l t i n g i n the wake being too short. The agreement between modal and f i n i t e d i f f e r e n c e s r e s u l t s i s excellent up to Re of about 20. For higher Re the length of the wakes and the p o s i t i o n s of the vortex centre predicted by modal method do not match those predicted by f i n i t e d i f f e r e n c e s . This discrepancy increases with Re, so that while modal r e s u l t s for Re = 40 are s t i l l acceptable, the predicted wakes at Re = 70 and 100 are much too short. Modal v o r t i c i t y predictions agree quite well with the d i r e c t f i n i t e element approach*,predictions again up to Re of about 40. The agreement deteriorates with increasing Re, though, and while for higher Re i n t h i s range modal values i n the regions of most i n t e r e s t behind and immediately above the cylinder are very close to f i n i t e element values, i n the region i n front and farther above the cylinder modal method predicts values which 55 do not seem to exist i n d i r e c t f i n i t e element r e s u l t s . At Re = 70 and 100 the predicted values i n these regions get even worse, hut, while the values i n the regions of most i n t e r e s t behind and close above the c y l i n d e r are not very close to those predicted by d i r e c t f i n i t e element approach, the general trend of e q u i - v o r t i c i t y l i n e s i s s t i l l reproduced. E q u i - v o r t i c i t y l i n e s predicted by modal method agree with the f i n i t e d ifferences r e s u l t s up to Re = 20 i n the regions of most i n t e r e s t . For higher Re modal r e s u l t s are somewhat wavy and some e q u i - v o r t i c i t y l i n e s have kinks showing the influence, of the g r i d . Again the discrepancies, as expected, increase with increasing Re. Equi-pressure l i n e s were not plotted because no comparisons were a v a i l a b l e . The f r i c t i o n drag c o e f f i c i e n t C^ ,, the pressure drag c o e f f i c i e n t 0^, the t o t a l drag c o e f f i c i e n t and the pressure values at the leading edge and the t r a i l i n g edge of the cylinder are l i s t e d i n Table (IV), however, and compared to the d i r e c t f i n i t e element approach r e s u l t s [6] and the f i n i t e d i fferences r e s u l t s [13,14], A l l these r e s u l t s are c o n s i s t e n t l y lower than the d i r e c t f i n i t e element r e s u l t s on the same g r i d for the whole range of Reynolds numbers under consideration. The d i f f e r e n c e i s drag c o e f f i c i e n t s , C^, 0^ and C^, increases with increasing Reynolds numbers with the maximum diff e r e n c e of about 30% at Re = 100. The d i f f e r e n c e i n the pressure value at the leading edge of the cylinder i s c o n s i s t e n t l y much higher than the d i f f e r e n c e i n the pressure value at the t r a i l i n g edge. Modal r e s u l t s are not c o n s i s t e n t l y higher than the f i n i t e differences r e s u l t s , considered exact herein, as the d i r e c t ' f i n i t e element'results are. Amazingly enough, though, a l l modal r e s u l t s f o r Re up to 40 are closer to the f i n i t e d i f f e r e n c e s r e s u l t s than the d i r e c t f i n i t e element approach r e s u l t s on the same g r i d . This can only be explained by reasoning that some approximation errors 56 introduced i n the pressure c a l c u l a t i o n s by f i n i t e element d i s c r e t i z a t i o n were cancelled by further approximation errors introduced by modal a n a l y s i s . From these r e s u l t s we conclude that good agreement with the d i r e c t f i n i t e element approach r e s u l t s f o r the steady case can be achieved by modal method for Reynolds numbers up to about 20 by employing about 50% of the modes. The agreement between modal r e s u l t s and the 'exact' r e s u l t s obtained by f i n i t e d i fferences i s good f o r Reynolds numbers only up to about 40, as i t i s also a function of f i n i t e element d i s c r e t i z a t i o n . For the transient analysis on the CYLFL 92-63 g r i d 62 modes were employed. F i r s t 52 modes, antisymmetric i n stream function, and f i r s t 10 modes symmetric i n stream function. The transformation of the nonlinear convective Q matrix to the eigenvector basis, spanned by these modes, would have been very expensive for the f u l l g r i d . We succeeded, however, to reduce the CPU time needed for the transformation:. 3.5 times by making use of the s p e c i a l properties of the transformed nonlinear convective § matrix, covered i n Chapter (7). I t was done as follows. F i r s t l y , the § matrix was formed for the upper h a l f of the CYLFL 92-63 g r i d . A l l boundary conditions imposed on the P., T , T and E boundaries were the same as for the f u l l g r i d , while l o u w a l l nodal v a r i a b l e s at the nodes located on the x-axis, which became a boundary, were l e f t f r e e . The number of degrees of freedom for t h i s g r i d was 139 and the bandwidth was 48, while for the f u l l g r i d they would have been 226 and 78, r e s p e c t i v e l y . Then the modes of the f u l l g r i d were truncated, so that only entries corresponding to the degrees of freedom of the half g r i d were retained and the transformation performed using these truncated modes. By multi p l y i n g a l l entries of the r e s u l t i n g array by a factor of 2.0, the transformed nonlinear convective matrix Q for the f u l l g r i d was f i n a l l y obtained. The mass matrix, needed i n the time algorithms for transformation from f i n i t e element coordinates 57 to generalized coordinates and v i c e versa, was computed for the f u l l g r i d . The constraints on the stream function, obtained from the boundary condition \\)=y imposed on the r and boundaries, were ip= 3.0 and ip= -3.0 at the two nodes located on the r boundary, IJJ= 20.0 at a l l nodes located on the top T u boundary and ip= -20.0 at a l l nodes located on the bottom boundary. In addition, the constraint ib =1 was imposed at a l l nodes situated on both the T and r boundaries, so that the t o t a l number of constraints i n the 1 u problem was 5. -The time integrations were started with an i n i t i a l s o l u t i o n obtained by perturbing the steady state r e s u l t s for a p a r t i c u l a r Reynolds number under consideration. The constrained v a r i a b l e s on the boundary were kept f i x e d throughout the time integrations which i n e f f e c t amounted to s p e c i f y i n g time independent boundary conditions. The c a l c u l a t i o n s at Re=20 were performed with two time int e g r a t i o n algorithms, the f i r s t one defined by Eq. (56) and the second one by Eq. (57). From these t r i a l c a l c u l a t i o n s we found, that while both algorithms yielded the same r e s u l t s the second algorithm was more e f f i c i e n t , because a much larger time step could be used without endangering the numerical s t a b i l i t y . So, i n a l l subsequent c a l c u l a t i o n s we e x c l u s i v e l y used t h i s algorithm. When the time analysis was performed at Re=20, 40 and 70 the steady state was reached i n a l l three cases. We would have expected that to happen at Re=20 and 40, but not at Re=70, where the flow should have become unsteady with o s c i l l a t i o n s of the downstream part of the wake. We reasoned, that t h i s was the r e s u l t of f i n i t e element d i s c r e t i z a t i o n errors and further truncation errors introduced by modal analysis, which employed only 10 modes, symmetric i n stream function. These errors introduced an ' a r t i f i c i a l v i s c o s i t y ' e f f e c t which lowered the ' e f f e c t i v e ' Reynolds number. This ' a r t i f i c i a l 58 v i s c o s i t y ' was, we thought, the main cause of numerical o v e r - s t a b i l i t y of our d i s c r e t e problem, Eq. (50), as compared to the actual hydrodynamic s t a b i l i t y inherent i n the governing p a r t i a l d i f f e r e n t i a l equation, Eq. (5). So we decided to increase the Reynolds number, thus also hoping to increase the ' e f f e c t i v e ' Reynolds number, and to perform the time integrations at Re=140, although judging from the steady r e s u l t s reported i n [6], the f i n i t e element g r i d was too crude at such a high Re. The integrations were performed with a time step At=0.15 seconds, while the test on whether the steady state had been achieved, Eq. (62), was set to be e=10 - 6. The i n i t i a l s o l u t i o n was a r b i t r a r i l y s p e c i f i e d by perturbing the steady state stream function values at the nodes located i n the upper ha l f domain y>0 upstream i n the neighbour-hood of the cylinder by 10% and at the mirror images of these nodes located i n the lower ha l f domain y<0 by 15%. The stream function r e s u l t s at time T=3, 6, 9, 12, and 16.5 seconds are plotted i n F i g . (16). These r e s u l t s seems to indi c a t e expected o s c i l l a t o r y behaviour, but f i n a l l y the steady state s o l u t i o n , antisymmetric i n stream function, was obtained, anyway. We speculate that t h i s happened p r i m a r i l y because of the wrong downstream boundary conditions whose influence was r e f l e c t e d to the computational domain f o r c i n g the flow back to steady state. Unfortunately, we could not come up with any better downstream boundary conditions and t h i s i s l e f t for some future study. We also note that the extension of the domain downstream, indicated i n F i g . (12-2) might not have been enough to reduce t h i s boundary conditions e f f e c t on.the region of i n t e r e s t . To summarize, i t seems that the i n a b i l i t y of our procedure to predict the expected unsteady flow was caused to some extent by truncation errors introduced by f i n i t e element d i s c r e t i z a t i o n and a d d i t i o n a l modal truncation errors, but p r i m a r i l y by the boundary conditions s p e c i f i e d on the downstream boundary Tn. 59 CHAPTER 9 Conclusions A modal f i n i t e element method for the steady state and the transient analyses of the plane flow of incompressible Newtonian f l u i d has been presented. The governing r e s t r i c t e d f u n c t i o n a l was d i s c r e t i z e d with a high p r e c i s i o n t r i a n g u l a r stream function f i n i t e element of C.^ class _p3jiS.«L Eigenvalue analysis was c a r r i e d out on the l i n e a r part of the problem obtained by deleting the nonlinear convective term. I t was found that the Lagrange m u l t i p l i e r s technique was computationally more e f f i c i e n t for incorporating the nonhomogeneous boundary conditions than the condensation procedure. In the l a t t e r procedure numerical d i f f i c u l t i e s were encountered when dealing with the nonlinear convective matrix. The matrix equations to be solved, when the Lagrange m u l t i p l i e r s technique i s applied, are nonsingular but i n d e f i n i t e . We found that t h i s posed no computational d i f f i c u l t i e s as the Gaussian elimination with p a r t i a l p i v o t i n g and forward and backward s u b s t i t u t i o n could conveniently be used to solve such equations. The number of modes was r e s t r i c t e d to 62 because of the computer core capacity and t h i s l i m i t a t i o n has to be further explored. We found that the computer time for the transformation of the nonlinear convective matrix to modal coordinates could be s i g n i f i c a n t l y reduced by taking advantage of the symmetric and antisymmetric properties of the modes. The transformation procedure was s t i l l quite expensive. Therefore, i t i s concluded that there w i l l be some p r a c t i c a l l i m i t on the s i z e of the problem that can be solved by modal method. 60 When modal method was applied to several flow problems we found that the number of modes, which were to be retained i n the analysis i n order to achieve reasonable r e s u l t s , increased with the refinement of the f i n i t e element g r i d . Furthermore, the choice of modes depended on the problem. In the square c a v i t y flow problem some higher modes had to be included i n order to approximate even the l i n e a r part, while for the flow around a c i r c u l a r cylinder no higher modes needed to be included. In the steady state analysis the convergent so l u t i o n was obtained i h 6 or l e s s i t e r a t i o n s , for the whole range of Reynolds numbers considered, regardless of the gr i d and the number of modes used. That i s , numerical i n s t a b i l i t i e s frequently encountered i n the f i n i t e d ifference method at higher Reynolds numbers were never experienced. It i s concluded, that t h i s new modal f i n i t e element method i n general y i e l d s good r e s u l t s i n the range of moderate Reynolds numbers with about 50% or less of the t o t a l modes. When the time dependent analysis was applied to the flow around a c i r c u l a r c y l i n d e r i t was concluded that the i n a b i l i t y to predict unsteady behaviour, expected at higher Reynolds numbers, was pri m a r i l y caused by the outflow boundary conditions. This lends hope that once these boundary conditions have been corrected or at le a s t t h e i r influence reduced by extending the computational domain further downstream, i t w i l l be possible to perform the time integrations on a greatly reduced number of equations by employing modal an a l y s i s . Hence, s i g n i f i c a n t savings i n computer costs can be achieved. F i n a l l y , as good r e s u l t s have been obtained i n t h i s thesis for moderate Reynolds numbers employing a greatly reduced number of l i n e a r modes, we speculate, that the extension of modal method to higher Reynolds numbers i s quite f e a s i b l e . It could be achieved, as suggested by Nickel;/[10] , by introducing modal decompositions for the subsequent nonlinear states based upon the tangent"matrix. 61 BIBLIOGRAPHY 1. Olson, M.D., "Comparison of Various F i n i t e Element Solution Methods for the Navier-Stokes Equations", Proceedings of International Conference on F i n i t e Elements i n Water Resources, Princeton U n i v e r s i t y , Princeton, U.S.A., Ju l y 12-16, 1976. 2. Finlayson, B.A., The Method of Weighted Residuals and V a r i a t i o n a l P r i n c i p l e s , Academic Press, New York, 1972. 3. Olson, M.D., "A V a r i a t i o n a l F i n i t e Element Method for Two-Dimensional Steady Viscous Flows", Proceedings of McGill-Engineering I n s t i t u t e of Canada Conference on F i n i t e Element Methods i n C i v i l Engineering, M c G i l l U n i v e r s i t y , Montreal, June 1-2, 1972, pp. 518-618. Also i n F i n i t e Elements i n Fl u i d s - Volume 1, (Ed. Gallagher, R.H., Oden, J.T., Taylor, C , and Zienkiewicz, O.C.), John Wiley & Sons, New York, 1975. 4. Yamada, Y.,L't!o, K. , Yokouchi, Y., Tamano, T. and Ohtsubo, T., " F i n i t e Element Analysis of Steady F l u i d and Metal Flow", F i n i t e Elements i n F l u i d s - Volume 1, (Ed. Gallagher, R.H., Oden, J.T., Taylor, C , and Zienkiewicz, O.C.), John Wiley & Sons, New York, 1975. 5. Tuann, S.Y., and Olson, M.D., "A Study of Various F i n i t e Element Methods for the Navier-Stokes Equations", S t r u c t u r a l Research Series, Report No. 14, Department of C i v i l Engineering, U n i v e r s i t y of B r i t i s h Columbia, Vancouver, Canada, May 1976. 6. Tuann, S.Y., and Olson, M.D., "Numerical Studies of the Flow Around a C i r c u l a r Cylinder by a F i n i t e Element Method", Structural Research 62 Series, Report No. 16, Department of C i v i l Engineering, U n i v e r s i t y of B r i t i s h Columbia, Vancouver, Canada, December 1976. 7. Cowper, G.R., Kosko, E., Lindberg, G.M., and Olson, M.D., "A High P r e c i s i o n Triangular Plate Bending Element','" Aeronautical Report LR-514, National Research Council of Canada, Ottawa, 1968. 8. Gallagher, R.H., F i n i t e Element Analysis, P r e n t i c e - H a l l , Englewood C l i f f s , New Jersey, 1975. 9. Clough, R.W., and Penzien, J . , Dynamics of Structures, McGraw-Hill, 1975. 10. N i c k e l , R.E., "Nonlinear Dynamics by Mode Superposition", Computer Methods i n Applied Mechanics and Engineering, 7 (1976), pp. 107-129, North-Holland. 11. S t r i c k l i n , J.A., Martinez, J.E., T i l l e r s o n , J.R., Hong, J.H., and Ha i s l e r , W.E., "Nonlinear Dynamic Analysis of Shells of Revolution by the Matrix Displacement Method", American I n s t i t u t e of Aeronautics and Astronautics Journal, 9 (4), pp. 629-636, A p r i l 1971. 12. Burggglij O.R., " A n a l y t i c a l and Numerical Studies of the Structure of Steady Separated Flows1,', Journal of F l u i d Mechanics, V o l . 24, Part 1, pp. 113-151, 1966. 13. Dennis, S.C.R., and Chang, G.Z., "Numerical Solutions for Steady Flow East a C i r c u l a r Cylinder at Reynolds Numbers Up to 100", Journal of F l u i d Mechanics, Vol. 42, Part 3, pp. 471-489, 1970. 63-14. Takami, H., and K e l l e r , H.B., "Steady Two-Dimensional Viscous Flow of an Incompressible F l u i d Past a C i r c u l a r Cylinder", Phys. F l u i d s Supplement I I , Vol. 12, pp. 11-51, 1969. 15. K e l l e r , H.B., and Takami, H., "Numerical Studies of Steady Viscous Flow Around Cylinders", Numerical Solutions of Nonlinear D i f f e r e n t i a l Equations, (Ed. Greenspan D.), John Wiley & Sons, pp. 115-140, 1966. Grid Number of Modes Eigenvalues Symmetric Antisymmetric T o t a l Symm. Antisymm. Tot a l Lowest Highest Lowest Highest Lowest Highest SQCA 12-13 10 5 15 53.3062 12799.6 108.352 9516.34 53.3062 12799.6 SQCA 36-29 40 31 71 51.6295 45973.8 92.4626 44085.0 51.6295 45973.8 SQCA 76-53 94 81 175 50.7516 49398.6 92.1604 48683.4 50.7516 49398.6 Number of Symmetric and Antisymmetric Modes and Lowest and Highest Eigenvalues for Three F i n i t e Element Grids Used for Square Cavity Flow Problem Simulation Number of Modes Reynolds Number, R 1 10 20 40 100 200 400 Symm. Antisymm. Total Error i n % at Midnodes5 10 5 15 0. P- 0. 0. 0. 0. 0. 9 5 14 0.76 0.78 0.85 1.08 1.63 1.96 1.94 8 5 13 1.40 1.42 1.47 1.65 2.06 2.26 1.31 8 4 12 1.40 1.42 1.49 1.75 2.36 2.75 9.18 10 1 11 0. 0. 0.27 1.02 11.33 17.22 42.73 TABLE II Error i n Stream Function at Midnode 5 for Various Number of Modes and Various Reynolds Numbers for SQCA 12-13 Grid a-Method Grid Number of Modes X v.c. y v . c . v.c. ^v. c. ^v. c. Grid for P l o t t i n g Symm. Antisymm , Total Modal SQCA 12-13 10 5 15 0.0 0.25 0.1355 5.5225 0.0 20 x 20 10 1 11 0.0 0.25 0.1355 5.5225 0.0 SQCA 36-29 27 15 42 0.0 0.25 0.0989 2.8059 -0.0892 26 10 36 0.0 0.25 0.0984 1.8501 -0.0835 SQCA 76-53 39 21 60 0.0 0.25 0.1040 2.7424 -0.0003 DMeet F i n i t e Element SQCA 12-13 - - - 0.0 0.25 0.1355 5.225 0.0 SQCA 36-29 - - - 0.0 0.25 0.0986 3.033 -0.0862 Approacl SQCA 76-53 - - - 0.0 0.25 0.0996 2.9775 -0.0928 [12] 50 x 50 - - - 0.0 0.27 0.100 3.20 0.0 TABLE III R = 0 Comparison of Numerical Results for Square Cavity Flow Probl Method Grid Number of Modes X V. c. ^v. c. v.c. Cv.c. P Grid for P l o t t i n g Symm. Antisymm, Total v. c. 10 5 15 0.0 0.25 0.1355 5.5196 -2.3096 SQCA 12-13 10 1 11 0.0 0.25 0.1355 5.5207 -2.3092 Modal • 27 15 42 0.0 0.25 0.0988 2.8050 -0.8911 SQCA 36-29 26 10 36 0.0 0.25 0.0984 1.8494 -0.8347 SQCA 76-53 39 21 60 0.0 0.25 0.1040 2.7276 -0.9828 20 x 20 SQCA 12-13 - — _ 0.0 0.25 0.1355 5.5196 92.3096 Dire c t F i n i t e Element SQCA 36-29 - - - 0.0 0.25 0.0985 3.0090 -0.8607 Approach SQCA 76-53 - - - 0.0 0.25 0.0995 2.9657 -0.9270 • R = 10 TABLE III (cont.) Comparison of Numerical Results f o r Square Cavity Flow Problem Method Grid Number of Modes X v.c. y v . c . v.c. C v.c. ^v. c. Grid for P l o t t i n g Symm. Antisymm Total Modal SQCA 12-13 10 5 15 0.0 0.25 0.1355 5.5107 -4.6126 20 x 20 10 1 11 0.0 0.25 0.1354 5.5157 -4.6096 SQCA 36-29 27 15 42 -0.05 0.25 0.0987 2.2439 -1.2252 26 10 36 -0.05 0.25 0.0986 1.9595 -1.1424 SQCA 76-53 39 21 60 -0.05 0.25 0.1046 2.8044 -1.5328 Dire c t F i n i t e Element Approach SQCA 12-13 - - - 0.0 0.25 0.1355 5.5107 -4.6126 SQCA 36-29 - - - -0.05 0.25 0.0985 2.2887 -1.1057 SQCA 76-53 - - - -0.05 0.25 0.0996 3.1068 -1.2662 R = 20 TABLE III (cont.) Comparison of Numerical Results f o r Square Cavity Flow Problem Method Grid Number of Modes " v . c . X v.c. " v ' v . c . ^v.c. v.c, v.c. 7" v.c. ? v.c. n ^v. c. Grid for P l o t t i n g Symm. Antisymm Total Modal SQCA 12-13 10 5 15 0.0 0.25 0.1354 5.4729 -9.1711 20 x 20 10 1 11 0.0 0.25 0.1350 5.4965 -9.1518 SQCA 36-29 27 15 42 -0.05 0.25 0.0995 2.2381 -3.1196 26 10 36 -0.05 0.25 0.0995 1.9080 -2.9461 SQCA 76-53 39 21 60 -0.10 0.25 0.1062 6.9153 -3.3526 Dire c t F i n i t e Element Approach SQCA 12-13 - - - 0.0 0.25 0.1354 5.4729 -9.1711 SQCA 36-29 - - - -0.05 0.25 0.0994 2.2729 -3.0044 SQCA 76-53 - - - -0.05 0.25 0.1003 3.1277 -3.2471 R = 40 TABLE III (cont.) Comparison of Numerical Results f o r Square Cavity Flow Probl em Method Grid Number of Modes X v.c. y v . c . v.c. ? v.c. '"'v.c. Grid for P l o t t i n g Symm. Antisymm Total Modal SQCA 12-13-10 5 15 0.0 0.20 0.1359 5.1545 -21.0284 20 x 20 10 1 11 0.0 0.25 0.1331 5.4002 -22.022 SQCA 36-29 27 15 42 -0.10 0.20 0.1045 2.3892 -9.2827 26 10 36 -0.10 0.20 0.1061 2.4479 -9.4048 SQCA 76-53 39 21 60 -0.10 0.25 0.1127 7.1780 -12.3943 Dire c t F i n i t e Element Approach SQCA 12-13 - - - 0.0 0.20 0.1359 5.1545 -21.0284 SQCA 36-29 - - - -0.10 0.20 0.1054 3.0346 -9.8903 SQCA 76-53 - - - -0.10 0.25 0.1037 2.8612 -9.7013 [12] 50 x 50 - - - -0.13 0.24 0.101 3.14 -18.1 TABLE .III (cont.) R = 100 Comparison of Numerical Results f o r Square Cavity Flow Problem Method Grid Number of Modes X v.c. y v . c . v.c. G v.c. P v . c . Grid for P l o t t i n g Symm. Antisymm , Total Modal SQCA 12-13 10 5 15 0.0 0.20 0.1428 4.9600 -42.3894 20 x 20 10 1 11 0.0 0.25 0.1307 5.2643 -41.5979 SQCA 36-29 27 15 42 -0.05 0.15 0.1152 2.7952 -25.2413 26 10 36 -0.10 0.15 0.1186 1.9509 -23.4828 SQCA 76-53 39 21 60 -0.15 0.15 0.1202 4.7300 -26.6176 Dire c t F i n i t e Element Approach SQCA 12-13 - - - 0.0 0.20 0.1428 4.9600 -42.3894 SQCA 36-29 - - - -0.10 0.15 0.1182 2.9340 -25.6447 R = 200 TABLE III (cont.) Comparison of Numerical Results f o r Square Cavity Flow Problem Method. Grid Number of Modes X v.c. y v . c . W v.c. v.c. P v . c . Grid for P l o t t i n g Symm. Antisymm Total Modal SQCA 12-13 10 5 15 0.0 0.15 0.1703 7.3450 -117.7910 20 x 20 10 1 11 0.0 0.25 0.1790 5.1650 -78.950 SQCA 36-29 27 15 42 -0.05 0.10 0.1246 1.9322 -57.0767 26 10 36 -0.05 0.05 0.1279 2.1609 -51.2457 SQCA 76-53 39 21 60 -0.05 0.05 0.1315 1.2139 -54.7102 D i r e c t F i n i t e Element Approach SQCA 12-13 - - - 0.0 0.15 0.1703 7.3450 -117.7910 SQCA 36-29 - - - -0.10 0.10 0.1319 4.6136 -58.6027 [5] SQCA 76-53 - - • - -0.056 0.083 0.1213 2.5099 -49.8779 i r r e g u l a r [12] *QC40 x 40 - - - -0.06 0.12 0.102 2.15 -71.7 TABLE III (cont.) R = 400 Comparison of Numerical Results for Square Cavity Flow Problem 73 Method & Grid Reynolds Number, R 1 5 7 10 20 40 70 100 F r i c t i o n Drag C o e f f i c i e n t , C f Modal CYLFL m~S% 7.525 2.301 1.816 1.408 0.848 0.506 0.325 0.239 [6] 42 7.573 2.386 1.918 1.529 0.997 0.653 0.447 0.347 [13] - 1.917 1.553 1.246 0.812 0.524 0.360 0.282 Pressure Drag C o e f f i c i e n t , C P Modal CYLFL 42-34 7.557 2.538 2.089 1.712 1.199 0.896 0.737 0.654 [C>6] 42 7.837 2.704 2.263 1.906 1.443 1.149 0.965 0.874 [13] - 2.199 1.868 1.600 1.233 0.998 0.852 0.774 Drag C o e f f i c i e n t , Modal CYLFL 42-34 15.082 4.839 3.905 3.120 2.047 1.402 1.062 0.893 [6] 42 15.410 5.091 4.181 3.435 2.440 1.802 1.412 1.221 [13]2 - 4.116 3.421 2.846 2.045 1.522 1.212 1.056 [14] 10.109 - • 3.303 2.800 2.013 1.536 - -Pressure at Leading Edge, P(if) Modal CYLFL '42-34 5.602 2.004 1.666 1.375 0.965 0.763 0.708 0.690 [6] 42 5.829 2.228 1.919 1.678 1.418 1.351 1.315 1.282 [13] - 1.872 1.660 1.489 1.269 1.144 1.085 1.060 [14] 3.905 - 1.637 1.474 1.261 1.141 — — Pressure at T r a i l i n g Edge, -P(0) Modal CYLFL %$rl% 3.876 1.222 1.013 0.845 0.638 0.545 0.502 0.479 [6] 42 3.845 1.242 1.050 0.896 0.698 0.580 0.488 0.436 [13] - 1.044 0.8700 0.742 0.589 0.509 0.439 0.393 [14] 2.719 - 0.783 0.670 0.537 0.512 - -TABLE IV Comparison of Numerical Results for Flow Around .a C i r c u l a r Cylinder 74 75 i r / = 3 y 2 - 2 y 3 ^ z (parabolic velocity A prof i le ) - -^ I {unknown velocity profile ) FIGURE 3 POISEUILLE FLOW PROBLEM CONFIGURATION 76 E V - 9 . 8 7 4 7 1 5 . 7 2 7 4 0 . 9 5 3 6 2 . 4 2 5 9 8 . 0 1 0 1 1 4 . 7 0 1 2 6 . 7 4 1 2 6 . 9 2 1 7 2 . 9 9 2 0 4 . 4 3 2 4 7 . 2 2 F I G U R E 4 M O D E S H A P E S & E I G E N V A L U E S F O R P O I S E U I L L E F L O W 5 M O D E S V 3 = - 1 . 5 7 8 0 8 6 M O D E S V 3 = 1 . 5 7 5 2 4 7 M O D E S 1 2 M O D E S V 3= 1.51168 V 3 = 1 . 5 0 0 0 0 F I G U R E 5 S T R E A M L I N E S & V E L O C I T I E S A T M I D N O D E 3 F O R P O I S E U I L L E F L O W 78 FIGURE 6 CIRCULATORY FLOW INDUCED BY A MOVING LID OVER A SQUARE CAVITY 79 (a) S Q C A 13-12 NE = 1 2 , NO = 13 NNS=15, N N P = 7 7 L B S =13, L B P = 5 4 ( b ) S Q C A 3 6 - 2 9 N E = 3 6 , N O = 2 9 NNS = 71 NNP=173 L B S = 32 L B P = 7 8 ( c ) S Q C A 7 6 - 5 3 N E = 7 6 , N O = 5 3 NNS = 175 L B S =44 NNP = 317 L B P = 7 8 1 \0A. 0 . 4 NE= NO. O F E L E M E N T S , NO = NO. O F NODES NN = NET NO. O F UNKNOWNS L B = H A L F BANDWIDTH S = STREAM FUNCTION P = PRESSURE FIGURE 7 FINITE E L E M E N T GRIDS FOR SQUARE CAVITY FLOW 80 n r 1 1 is fit. (cm) . E V = 5 3 . 3 0 6 9 9 . 8 0 4 1 0 8 . . 3 5 2 1 5 3 . 9 4 5 i l l (ID 1 7 0 . 2 5 6 1 8 5 . 5 2 4 3 7 2 . 2 4 1 3 9 1 , 7 7 8 4 0 8 . 2 8 4 7 9 6 . 3 8 9 1 5 7 4 . 9 2 9 0 2 4 . 0 2 i ^ ^ ^ ^ ^ ^ ( S T 9 5 1 6 . 3 4 1 0 0 1 5 . . 2 1 2 7 9 9 . 6 F I G U R E , 8 - 1 M O D E S H A P E S & E I G E N V A L U E S F O R S Q C A 1 2 - 1 3 G R I D 81 EV= 51.629 8 1 . 9 8 0 9 2 . 4 6 2 ail 0 1 5 3 . 5 6 9 ,--vv-, 1 6 5 . 7 3 5 3 . 1 2 4 2 4 2 . 6 9 7 2 5 5 . 3 8 5 2 5 7 2 1 1 m m 1 2 9 . 1 6 4 fit 1 9 3 . 1 0 8 IP m 2 8 2 . 7 7 8 2 9 2 . 5 9 9 3 4 4 . 6 7 3 3 4 8 . 5 7 9 3 5 6 . 8 6 2 ~ F I G U R E 8 - 2 F I R S T 3 0 M O D E S H A P E S & E I G E N V A L U E S F O R S Q C A 3 6 - 2 9 G R I D 82 3 7 5 . 3 0 6 4 2 4 . 8 1 8 4 2 6 . 4 8 5 Star 4 4 1 . 9 7 5 479.444 4 8 2 . 6 3 9 5 1 8 . 4 9 8 5 5 3 . 0 7 9 SSI i i H 592.345 5 9 3 . 5 5 2 5 9 5 . 2 . 7 1 6 0 3 . 9 4 4 6 5 9 3 2 3 6 7 9 . 2 8 4 F I G U R E . 8 - 2 ( C O N T ) F I R S T 3 0 M O D E S H A P E S & E I G E N V A L U E S F O R S Q C A 3 6 - 2 9 G R I D 83 Figure 9 Streamlines for Square Cavity Flow for Various Reynolds Numbers R Contours represent equal steps i n stream function, unless s p e c i f i e d otherwise. The vortex centre i s marked by a cross. The numbers below the fig u r e s r e f e r to the number of modes used i n the c a l c u l a t i o n s , while FE denotes the d i r e c t f i n i t e element approach. 84 F I G U R E 9 -1 R = 0 85 FIGURE 9 - 2 R = 10 86 FIGURE 9 - 3 R=20 87 88 FIGURE 9 - 5 R = 100 89 91 Figure 10 E q u i - V o r t i c i t y Lines for Square Cavity Flow for Various Reynolds Numbers R The vortex centre i s marked by a cross. The numbers below the figures r e f e r to the number of modes used i n the c a l c u l a t i o n s , while FE denotes the d i r e c t f i n i t e element approach. Contour l e v e l s are l a b e l l e d only on the c e n t r a l f i g u r e . For R = 0, 10,20,40 and 100, i|> = -1.0,0.0,1.0,3.0,5.0 contours are plotted, and for R = 200 and 400, i> = -1.0,0.0,1.0,2.0,2.2, and 3.0 contours. ( c ) 6 0 F E R E F 1 2 F I G U R E 1 0 - 1 R = 0 93 9 6 97 (a) 11 15 FE (b) 6 0 FIGURE 10-6 R = 200 98 (a) 11 15 FE (b) 36 42 FE ( O 6 0 REF 5 REF 12 FIGURE 10 -7 R = 4 0 0 99 F i g u r e 11 E q u i - P r e s s u r e L i n e s f o r Square C a v i t y F l o w f o r V a r i o u s R e y n o l d s Numbers R The v o r t e x c e n t r e i s marked by a c r o s s . The numbers b e l o w t h e f i g u r e s r e f e r t o t h e number o f modes u sed i n t h e c a l c u l a t i o n s , w h i l e FE d e n o t e s t h e d i r e c t f i n i t e e l emen t a p p r o a c h . C o n t o u r l e v e l s a r e l a b e l l e d o n l y on t h e c e n t r a l f i g u r e . F o r R = 0 , 1 0 , 2 0 , and 4 0 , p = - 2 0 . 0 , - 1 0 . 0 , - 5 . 0 , - 1 . 0 , 0 . 0 , 1 . 0 , 5 . 0 , 1 0 . 0 , and 2 0 . 0 , c o n t o u r s a r e p l o t t e d , f o r R = 100 p = - 1 5 . 0 , - 7 . 5 , - 0 . 5 , 0 .0 and 15 .0 c o n t o u r s , f o r R = 200 p = - 2 0 . 0 , - 1 0 . 0 , - 5 . 0 , - 1 . 0 , 0 . 0 , 3 0 . 0 , and 60 .0 c o n t o u r s and f o r R = 400 p = - 6 0 . 0 , - 3 0 . 0 , 0 . 0 , 3 0 . 0 and 60 .0 c o n t o u r s . 103 (a) 11 15 FE ( b ) 60 F E FIGURE 11-4 R = 40 (a) 36 42 FE ( c ) 60 FE REF 12 FIGURE 11 - 5 R=100 105 (a) 11 15 FE (b) 36 42 FE \ ( c ) 6 0 FIGURE 1 1 - 6 R=200 FIGURE 12-1 FINITE ELEMENT GRIDS FOR A FLOW AROUND A CIRCULAR CYLINDER C Y L F L 8 4 - 5 8 GRID FIGURE 12-2 FINITE ELEMENT GRIDS FOR A FLOW AROUND A CIRCULAR CYLINDER C Y L F L 92 -63 GRID ' 109 EV= 0.020 0.248 0.517 0.826 1.193 FIGURE 13-1 0.070 0.378 0.986 1.424 0.191 0.428 0.669 1.085 1.503 FIRST 15 ANT ISYMMETRIC MODE SHAPES & EIGENVALUES FOR CYLFL 8 4 - 5 8 GRID 110 EV= 0.016 0.186 0.259 0.397 0.510 0.614 0 .774 0.795 1.020 1.180 1.487 1.503 1.730 1.849 2.051 FIGURE 13-2 FIRST 20 ANTISYMMETRIC AND FIRST 10 SYMMETRIC MODE SHAPES & EIGENVALUES FOR C Y L F L 9 2 - 6 3 GRID I l l EV= 2.450 2.844 2.955 3.481 3.607 0 . 0 8 6 0.137 0.186 0 . 3 8 9 0.417 0.476 0.597 0.706 0 . 8 9 0 0.962 FIGURE 13-2 ( CONT. ) FIRST 20 ANTISYMMETRIC AND FIRST 10 SYMMETRIC MODE SHAPES & E IGENVALUES FOR C Y L F L 9 2 - 6 3 GRID 112 Figure 14 Steamlines for Steady Flow Around a C i r c u l a r Cylinder f o r Various Reynolds Numbers R Values of the dimensionless stream function ib are shown for each streamline on the bottom f i g u r e . Values of i|> for closed streamlines are given below for a s p e c i f i e d value of the Reynolds number (4) R - 10 : *c = -0.0002 (5) R = 20 : *c =' -0.0080, -0.0058 (6) R = 40 : V = -0.0328, -0.0246, -0.0164, -0 (7) R = 70 : *c -0.07, -0.06,-0.035, -0.023 (8) R = 100: -0.1, -0.08, -0.05, -0.035 The number below the top f i g u r e denotes the number of modes used i n the c a l c u l a t i o n s , while Ref. 6 r e f e r s to the d i r e c t f i n i t e element approach using the same g r i d . 113 REF 14 FIGURE 14-1 R=1 114 REE 13 FIGURE 14-2 R=5 115 REF 13 FIGURE 14-3 R = 7 116 REF 13 FIGURE 14-4 R = 1 0 117 REF 13 FIGURE 14 - 5 R=20 62 REF 13 FIGURE 14-6 R=40 119 REF 13 FIGURE 1 4 - 7 R = 70 120 REF 13 FIGURE 14-8 R=100 121 Figure 15 Equi-Vorticity Lines for Steady Flow Around a Circular Cylinder for Various Reynolds Numbers R Values of the negative dimensionless vort ic i ty £ are shown for each equi-vorticity l ine . The number below the top figure denotes the number of modes used in the calculations, while Ref. 6 refers to the direct f in i te element approach using the same gr id . 122 123 L, Jo.s REF 6 FIGURE 1 5 - 2 R = 5 124 REF 14 FIGURE 1 5 - 3 R=7 1 2 5 REF 15 FIGURE 1 5 - 4 R = 10 126 REF 14 FIGURE 15-5 R = 20 127 FIGURE 15-6 R= 40 128 REF 13 FIGURE 1 5 - 7 R = 70 129 REF 13 FIGURE 15-8 R = 100 130 Figure 16 Streamlines f o r Transient Flow Around a C i r c u l a r Cylinder at Reynolds Nonce'-. 1T;=140 for Various Time Instants The following streamlines are plotted ii = -0.631, -0.4115, -0.129, -0.03, -0.0175, -0/0115 0.0, 0.0115, 0.0175, 0.03, 0.129, 0.4115, 3nd30.631 131 T = 9 SECONDS FIGURE 16 STREAMLINES FOR TRANSIENT FLOW AROUND A CIRCULAR CYLINDER AT REYNOLDS NO. R = 140 FOR VARIOUS TIME INSTANTS 132 T = 16.5 SECONDS FIGURE FLOW NO. R 16 ( CONT. ) STREAMLINES FOR TRANSIENT AROUND A CIRCULAR CYLINDER AT REYNOLDS = 140 FOR VARIOUS TIME INSTANTS
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Modal finite element method for the navier-stokes equations
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Modal finite element method for the navier-stokes equations Savor, Zlatko 1977
pdf
Page Metadata
Item Metadata
Title | Modal finite element method for the navier-stokes equations |
Creator |
Savor, Zlatko |
Date Issued | 1977 |
Description | A modal finite element method is presented for the steady state and transient analyses of the plane flow of incompressible Newtonian fluid. The governing restricted functional is discretized with a high precision triangular stream function finite element. Eigenvalue analysis is carried out on the resulting discretized problem, under the assumption that the nonlinear convective term is equal to zero. After truncating at various levels of approximation to obtain a reduced number of modes, the transformation to the new vector space, spanned by these modes is performed. Advantage is taken of the ..symmetric and the antisymmetric properties of the modes in order to simplify the calculations. The Lagrange multipliers technique is employed to {incorporate the nonhomo-geneous boundary conditions. The steady state analysis is carried out by utilizing the Newton-Raphson iterative procedure. The algorithm for transient analysis is based upon backward finite differences in time. Numerical results are presented for the fully developed plane Poiseuille flow, the flow in a square cavity, and the flow over a circular cylinder problems. These resultscfor the steady state are compared with the results obtained by direct finite element approach on the same grids and the results obtained by finite differences technique. It is concluded that the number of modes, which are to be retained in the analysis in order to achieve reasonable results, increases with the refinement of the finite element grid. Furthermore, the choice of modes to be used depends on the problem. Finally it is established, that this new modal method yields good results in the range of moderate Reynolds numbers with about 50% or less of the modes of the problem. This, in turn, means that the time integrations can be performed on a greatly reduced number of equations and hence potential savings in computer time are significant. |
Subject |
Navier-Stokes equations |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-02-19 |
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. |
IsShownAt | 10.14288/1.0062684 |
URI | http://hdl.handle.net/2429/20518 |
Degree |
Master of Science - MSc |
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_1977_A7 S39.pdf [ 6.28MB ]
- Metadata
- JSON: 831-1.0062684.json
- JSON-LD: 831-1.0062684-ld.json
- RDF/XML (Pretty): 831-1.0062684-rdf.xml
- RDF/JSON: 831-1.0062684-rdf.json
- Turtle: 831-1.0062684-turtle.txt
- N-Triples: 831-1.0062684-rdf-ntriples.txt
- Original Record: 831-1.0062684-source.json
- Full Text
- 831-1.0062684-fulltext.txt
- Citation
- 831-1.0062684.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062684/manifest