Citation |

- Permanent Link:
- http://digital.auraria.edu/AA00001980/00001
## Material Information- Title:
- Finite element modeling for energy dissipation in a pressure flow
- Creator:
- Guo, Chwen-geng
- Place of Publication:
- Denver, CO
- Publisher:
- University of Colorado Denver
- Publication Date:
- 1988
- Language:
- English
- Physical Description:
- 126 leaves : illustrations ; 28 cm
## Subjects- Subjects / Keywords:
- Streamflow -- Experiments ( lcsh )
- Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Bibliography:
- Includes bibliographical references (leaves 102-104).
- Thesis:
- Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Civil Engineering
- Statement of Responsibility:
- by Chwen-geng Guo.
## Record Information- Source Institution:
- |University of Colorado Denver
- Holding Location:
- |Auraria Library
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Resource Identifier:
- 19783233 ( OCLC )
ocm19783233 - Classification:
- LD1190.E53 1988m .G86 ( lcc )
## Auraria Membership |

Downloads |

## This item has the following downloads: |

Full Text |

FINITE ELEMENT MODELING FOR ENERGY
DISSIPATION IN A PRESSURE FLOW by Chwen-geng Guo .S., National Cheng-Kung University, Taiwan, 1978 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Department of Civil Engineering 1988 This thesis for the Master of Science Degree by Chwen-geng Guo has been approved for the Department of civil Engineering by Date (Zb /fW Ronald C. Scherer Guo, Chwen-geng (M.S. Civil Engineering) Finite Element Modeling for Energy Dissipation in a Pressure Flow Thesis directed by Assistant Professor James C. Y. Guo. A finite element numerical model using quadratic isoparametric elements was developed to analyze a two- dimensional pressure flow and to predict the flow energy losses associated with the friction and cross sectional contraction. The unsteady flow governing equations were converted into stream function and vorticity diffusion equations. Initial conditions were determined by potential flow and boundary conditions included a uniform incoming flow and non-slip condition on the solid wall. Flow cases studied included the opening ratios ranging from 0.3 to 1 and Reynolds numbers varying from 50 to 200. In general, this numerical approach provides stable computations for flow to reach its steady state when >- Reynolds number is less than 200 and opening ratio is larger than 0.3. Although for high Reynolds number flows the exit velocity profile may slightly deviate from the exact velocity profile, the predicted friction loss and orifice loss still give good agreements with the theoretical solutions and laboratory data. ACKNOWLEDGMENTS This project was supported by the Denver Center for the Performing Arts and funded by the Helen G. Bmfils Foundation. I would like to thank Dr. Ronald C. Scherer, senior scientist, Denver Center for the Performing Arts, and Dr. James C. Y. Guo, assistant professor, University of Colorado at Denver, for their advice and encouragement on this project. Also, I wish to express my thanks to thesis committee, Dr. David W. Hubly, and Dr. John R. Mays, for introducing me to the Fluid Mechanics and Finite Element Method. Finally, I would like to give my special thanks to my parents, my wife and children, for their understanding throughout the entire project. CONTENTS CHAPTER I. INTRODUCTION ...................................... 1 Objective ....................................... 1 Scope of the study .............................. 2 II. LITERATURE REVIEW ................................. 3 Energy Losses in An incompressible Flow ...... 3 Numerical Modeling of An Incompressible Flow . 12 III. GOVERNING EQUATIONS OF FLOW MOTIONS .............. 16 Flow Motions and Descriptions .......... 16 Stream Function-Vorticity Equations of Flow Motion .................... 19 Nondimensionalization of the Governing Equations ................... 23 Boundary Conditions ............................ 25 IV. NUMERICAL MODELING BY FINITE ELEMENT SCHEME .... 29 Basic Concepts ................................ 29 Derivation oif Finite Element Equations of Flow Governing Equations ......... 35 Element Shape Function ........................ 40 Numerical Initial and Boundary Conditions .... 49 Computational Flow Chart ...................... 53 V. CASE STUDIES AND RESULTS ......................... 57 Case One : Study of Flow Patterns ............ 63 vi Case Two : Prediction of Energy Losses ....... 77 VI. CONCLUSIONS ...................................... 99 REFERENCE ................................ . ........ 102 APPENDICES A. THE DERIVATION OF ELEMENT EQUATION BY GAUSS-GREEN THEOREM ....................... 105 B. THE INTEGRATION OF SHAPE FUNCTIONS .............. 109 C. DEFINITION OF SYMBOLS ........................... 115 D. THE FORTRAN SOURCE CODE ......................... 117 TABLES Table 1. Computational Flow Chart .......................... 55 2. Summary of the Flow Parameter Used in the Case Studies ................................ 58 3. Computed Velocity Distribution for Opening Ratio = 0.7 and Reynolds number = 50 .. 73 4. Comparison of Peak Vorticity....................... 76 5. Summary of Computational Time Intervals and the Required Times for Flow to Reach Steady State .............................. 78 6. Computed Velocity Distribution in a Poiseuille Flow for Reynolds Number = 50 ........................ 81 7. Computed Velocity Distribution in a Poiseuille Flow for Reynolds Number = 100 ....................... 81 8. Comparison Between the Computed Exit Velocity Profile at X/B = 8.7 and The Analytical Profile for Poiseuille Flow............................. 84 9. Comparison of Friction Loss Between Computational and Analytical Predictions ........ 87 10. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile for the cases with a contraction .............. 94 11. Distribution of Flow Pressure with a Contraction .............................. 95 Comparison of Discharge Coefficients Between the Computed Results and Experimental Data .... 97 12. FIGURES Figure 1. Illustration of Sudden Expansion ................ 7 2. Illustration of the Flow Domain with Orifice Plate ........... 10 3. Rotational and Irrotational Flows between Two Infinite Plates .................... 18 4. Boundary Conditions of Poiseuille Flow ........... 27 5. Illustration of One-Dimensional Linear Finite Element .......................... 30 6. The Example of Eight-Node Quadratic Quadrilateral Elements ............... 34 7. The Eight-Node Quadratic Quadrilateral Isoparametric Element with Natural Coordinates ............................ 43 8. Finite Element Net Work for the Case with Opening Ratio of 1 ............. 59 9. Finite Element Net Work for the Case with Opening Ratio of 0.7 ........... 60 10. Finite Element Net Work for the Case with Opening Ratio of 0.5 ........... 61 11. Finite Element Net Work for the Case with Opening Ratio of 0.3 ........... 62 12. The Steady State Distribution of Stream Function and Contours of Vorticity for Opening Ratio = 0.7 and Reynolds Number = 50 .......................... 64 13. Illustration the Evolution of Flow Patterns Changing From Potential Flow to Its Steady State for Opening Ratio = 0.7 and Reynolds Number =50 ........................ 66 ix 14. The Variation of Velocity Profile along the Flow Direction for Opening Ratio = 0.7 and Reynolds number = 50 ..................... 72 15. The Distribution of the Wall Vorticity for Opening Ratio = 0.7 and Reynolds number = 50 ......................... 75 16. The Steady-State Velocity Profiles In A Poiseuille Flow for Reynolds Number = 50 ......................... 79 17. The Steady-State Velocity Profiles In A Poiseuille Flow for Reynolds Number = 100 ........................ 80 18. The Variation of Velocity Profile Along the Flow Direction for Opening Ratio = l and Reynolds number = 50 ..................... 82 19. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile for Poiseuille Flow with X/B = 8.7 ........... 83 20. Comparison of Friction Loss Between Computational and Analytical Predictions for Poiseuille Flow with X/B = 8.7 ........... 88 21. The Variation of Pressure Along the Flow Direction for Poiseuille Flow with X/B = 8.7 ............................... 89 22. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile for Opening Ratio = 0.7 .......................... 91 23. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile for Opening Ratio = 0.5 .......................... 92 24. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile for Opening Ratio = 0.3 .......................... 93 25. Computed Pressure Drops......................... 96 26. Comparison of Discharge Coefficients Between Computed Results and Experimental Data ....... 98 27. Coordinate Transformation .........*............ Ill CHAPTER I INTRODUCTION Objective Energy dissipation is an important subject to many professions. Examples include crude oil delivered in pipes water supply systems, air flow through the larynx and blood flow in arteries. Numerous studies have been conducted to develop analytical and empirical equations to estimate energy losses associated with the shear friction and eddy motions. Recently, energy loss in the larynx flow has attracted much attention because of the lack understanding of flows near the glottis. The Denver Center for the Performing Arts, in association with the University of Iowa has been engaged in a laboratory investigation of this phenomena. Since the energy dissipation in a flow is related to both flow and fluid properties, laboratory arrangements may not easily cover wide ranges of flow variables such as Reynolds number. As a result, numerical modeling is considered a possible helpful alternative. The purpose of this study is to develop a numerical model to analyze a two-dimensional pressure flow and to predict the flow energy losses associated with the wall 2 friction and cross sectional contraction. In this work, plate orifices are considered in order to establish the finite element technique before applying the technique to laryngeal flow problems. Scope of the Study After a review of numerical methods, the finite element method was selected to convert the nonlinear governing equations of unsteady flow into their numerical forms and to include the complicated boundary conditions in the numerical integrals. Solutions in terms of the stream function, vorticity, velocity and pressure drop at each time step were obtained from the numerical computations for unsteady flow. The steady state solutions were achieved by allowing the successive time-dependent solutions to converge. It is found that for low Reynolds number flows, comparisons between the predicted flow variables and analytical solutions are in good agreement. Although this developed numerical scheme is only for two-dimensional flow, it is believed that the approach is also suitable for a similar three-dimensional flow problem. CHAPTER TWO LITERATURE REVIEW Energy Losses in an Incompressible Flow The main objective of this study is to develop a numerical model for the prediction of flow patterns and energy losses in a two dimensional pressure flow. In general, there are two kinds of energy losses involved in a pressure flow. They are friction loss and minor loss. Friction in a flow is caused by the non-uniform velocity distribution and fluid viscosity. Friction results in energy dissipation that converts mechanical energy into thermal energy. This energy loss to the environment is irreversible. Friction loss is determined by the length of the conduits, flow Reynolds number, and wall roughness. Using dimensional analysis, the friction loss equation in a pressure flow may be derived as follows (Daugherty, Franzini, and Finnemore, 1985): where, hT= the head loss due to friction, L= the length of the conduit, R^ = hydraulic radius, V= cross sectional average flow velocity, g = acceleration of gravity, and C^= Darcy friction factor. 4 For a pipe flow, substituting Rh=D/4, in which D is the pipe diameter, into Eq. 2.1, the Darcy-Weisbach equation is obtained: h L f L V D 2g (Eq. 2.2) where D = the pipe diameter, and f = 4C^. To determine the pipe friction factor f the pipe flow experiments were conducted as early as the mid-1800s. For high Reynolds number flow, flow in rough pipe resulted in a higher energy loss than in smooth pipe. In 1933, J. Nikuradse experimented using artificially roughened pipes and obtained reliable information for friction factor (Gerhart and Gross, 1985). In 1944, the Moody Chart for Friction Factor f versus Reynolds number and relative roughness was published and, since then, it has been widely adopted by engineers. To expand the friction loss, Hagen and Poiseuille successfully derived analytical solutions in 1840 for laminar flows in a circular pipe and between two plates. They stated that for laminar flow in a circular pipe, the loss of head in friction is proportional to the first power of the cross sectional average velocity and is given by IL L h = 32------------= V (Eq. 2.3) L <5 D2 where ii = absolute viscosity, and <: = the density of fluid. By equating Eq. 2.2 and 2.3, we have 5 64 M> f = --------- expressed as (Eq. 2.4) Eq. 2.4 can be further 64 f = ------ (Eq. 2.5) Re in which Re= flow Reynolds number. From Eq. 2.5, it is found that in laminar flow the friction loss is only a function of Reynolds number and independent of the roughness of the pipe wall. On the Moody chart, the friction factor f for laminar flow is a straight line. This fact is not only true for a laminar pipe flow but also for a laminar flow in a noncircular duct. Jones (1976) expanded Eq. 2.4 to duct flow with the circular diameter D replaced by the hydraulic radius, D^, and conducted a series of experiments for various duct configurations. He suggested that for the laminar flow in a rectangular duct with width, a, and height, b, the friction factor may be described by k f = ------ (Eq. 2.6) Reh where k = 64 2 11 b 3 24 a (2 (Eq. 2.7) Reh - ___!i_ and Dh = hydraulic diameter. 6 For the laminar flow between two infinite plates, the value of variable a, in Eq. 2.7, is infinite. As a result, k is equal to 96 when is equal to 4B (Blevins 1984, and Gerhart and Gross i985). Losses due to the change of local pipe configurations are termed minor losses. The flow in a pipe system may need to pass through fittings, bends, or abrupt changes in the pipe configuration. These kinds of changes in local pipe geometries generate a great deal of disorderly motion and result in dissipation of flow momentum. Whenever the flow velocity is altered, eddy currents are induced and a loss of energy in addition to friction loss is created. Although minor losses are usually confined to a short length of flow path, the effects may not disappear for a considerable distance downstream. The minor loss may be expressed as v2 hT = K ----- (Eq. 2.8) L 2g where K = loss coefficient, which must be determined experimentally for each situation. The only loss coefficient which can be calculated by simple analytical methods is that the axisymmetric sudden - expansion as shown in Figure 1. This problem was first treated by Borda and Carnot. Using continuity, linear momentum, and energy equations and neglecting the shear stress along the wall, the head loss for sudden expansion Figure 1 Illustration of sudden expansion 8 can be expressed as hL = ( 1 - )2 ( ZL 2g (Eq. 2.9) where A1= the cross sectional area at position 1, A2= the cross sectional area at position 2, and V1= the average velocity at position 2. Comparing Eq. 2.9 with Eq. 2.8 yields A1 2 K = ( 1-----p- )z (Eq. 2.10) A2 Although we may be uneasy about the assumptions of the Borda-Carnot approach and feel that Eq. 2.10 is too simple to evaluate the accurate energy loss for sudden expansion, the laboratory experiments showed that the difference between observed and predicted values of K in Eq. 2.10 were within a few percent (Vennard 1969). Based on this fact, it can be assumed that minor losses may be considered for practical purposes to occur in a much more localized flow region. As to the minor losses other than sudden expansion, they all have to be determined by experiments. In 1885, Weisbach followed the approach used by Borda-Carnot to observe the energy losses of abrupt contraction and a fair - estimate of loss coefficients was obtained. In 1952, Gibson's studies on diffusers of conical form grasped the trend of the change of loss coefficient, which are defined as a function of cone angle and area ratio. In 1942, 9 Langhaar successfully calculated the loss of head at entrance by using boundary-layer theory. As to the losses in bends and elbows, there were studied by Beij and Hofmann in 1933. Hofmann gave loss coefficients for wide ranges of bend shape, roughness, and Reynolds number. Nigro (1978) and Pao(1961) indicated that for a flow with Reynolds number below 104, the discharge coefficient becomes increasingly sensitive to Reynolds number, first rising and then falling as Reynolds number decreases. The same phenomenon as indicated by Negro and Pao was observed by Scherer (1981) in his experiments on laryngeal flow for a convergent-divergent glottis and opening ratios d/B less than 0.05, where d = the glottal diameter, and B = the subglottal (tracheal) lateral diameter. The discharge coefficient for the pipe flow used by Nigro is Ud " 2 A p h [ i (-f-r ] (Eg. 2.11) where C = discharge coefficient, 11^= the mean flow velocity through the orifice throat, "a = the cross sectional area of orifice, A = the cross sectional area of the pipe, Ap = the pressure drop between the flange taps as shown in Fig 2, and = fluid density. Considering a unit width of a two-dimensional flow, Eq. 2.11 can further be simplified to evaluate the flow velocity between two infinite plates. Figure 2 Illustration of the flow domain with orifice plate 11 ud = 2 A p h [ 1 (-f-r ] (Eq. 2.12) where d/B = the opening ratio between two plates, d = half of the opening width at the orifice, B = half of the width between the two plates. As indicated by Nigro, for laminar flow through a thin square-edged orifice with Reynolds number between 100 and 200, the discharge coefficient for flange taps is approximately .70 to .80. The values of loss coefficients obtained by the pioneer researchers mentioned above were continually modified by researchers. There are presently many handbooks and technical papers published with experimental data (Idel'chek 1960, and ASHRAE Handbook of Fundamentals). Unfortunately, owing to different experiments conducted under different experimental conditions, data values for minor losses may be different for the same flow configuration. Therefore, the most representative one is difficult to determine. Furthermore, experiments for the evaluation of minor losses are common for turbulent flows, whereas the empirical data for laminar flows are less common. These are the main reasons why this numerical study for estimating the energy losses for a laminar flow is needed. In this study, the laminar flow patterns and energy losses associated with the plate orifice as shown in Figure 12 2, are numerically calculated. Before the selection of numerical methods, the following literature review was performed. Numerical Modeling of An Incompressible Flow There are three approaches by which we can solve fluid problems, namely 1. Experimental 2. Theoretical? and 3. Numerical The experimental approach has the capability of producing the most realistic answers for many flow problems. However, its cost is the highest one among these three methods. In theoretical approaches, assumptions should be made in order to make the problem tractable. Therefore, the problems that can be solved by theoretical approaches are restricted to simple boundary conditions and usually limited to linear problems. However, the shortcomings of experimental and theoretical methods do not exist in numerical approach. In the numerical approach, a limited number of assumptions are made and a digital computer is used to solve the governing eguations of a flow problem. There is no restriction to flow linearity and even complicated geometry can be treated. Furthermore, the computational costs are always much less 13 than that of laboratory equipment for experiments. It is why the numerical method has become more important when the computing facilities become more efficient and accessible. Most people attribute the first definite work of numerical approach to Richardson (1910) who introduced point iterative schemes to numerically solve Laplace's equation. His scheme used data available from the previous iteration to update each value of the unknown and was called the finite difference method (FDM). Richardson actually carried out calculations for the stress distributions in a masonry dam with great success. Up to the late 1930's, FDM was used for solving flow problems. Allen and Southwell (1955) applied relaxation scheme to solve the incompressible, viscous flow around a cylinder. Although they were not the first scientists to solve these kinds of problems numerically, their results did capture the flow phenomenon and had good accuracy. At the same time, John von Neumann developed his method for evaluating the stability of numerical methods for solving time-dependent problems. Douglas and Rachford (1950) developed an implicit method for solving parabolic and elliptic partial differential equations. From then on, the basic knowledge for solving fluid dynamics problems by FDM was complete. Before the finite element method (FEM) was fully developed, FDM was used to solve flow problems. The 14 advantage of FEM is that its rate of convergence is much faster than that of FDM. The finite element method was originally used by aircraft structural engineers in the 1950's to analyze a large system of structural elements in the aircraft. Usually, there are two ways to formulate element approximate equations: one is the variational method, which has a close relationship to the classical variational concepts of the Rayleigh-Ritz method (Reyleigh 1877;Ritz 1909); the other one is the weighted residual method, which is modeled after the well-known method of Galerkin (1915). In the weighted residual method, the error in the approximate satisfaction of the governing equations is not set to zero, but instead its integral with respect to the selected "weights" is required to vanish. The application of FEM to nonstructural problems such as fluid flows was initiated by Zienkiewicz and Cheung (1965). Pioneering research has employed weighted residual methods to formulate approximate element equations. Oden (1972) was among the first to derive the basic theoretical analog for the Navier-Stokes equations. Baker (1971, 1973) published his results for the analysis of incompressible flow using the Galerkin method. Olson (1974) used a pseudo- variational finite-element algorithm for solving the biharmonic equation for stream functions. Lynn(1974) utilized the least-squares methods to develop a finite- 15 element algorithm for laminar boundary-layer flow. In 1978, Backer successfully extended his studies from 2- dimensional laminar flow to 3-dimensional turbulent flow. From that time on, FEM was widely adopted for solving fluid dynamics problems. In this study, the numerical model for solving laminar, incompressible flow will be derived using FEM because of its flexibility of element shape and efficient computation. Galerkin's method is used to weight residuals in the approximate governing equations of flow. The procedure of formulation will be discussed in Chapter Four. CHAPTER THREE GOVERNING EQUATIONS OF FLOW MOTIONS Flow Motions and Descriptions Using viscosity as a criterion, fluid may be classified as real fluid and ideal fluid. For a real fluid, flow motion results in shear forces between fluid layers. But for an ideal flow, the shear force is assumed to be zero everywhere. Although this is an idealized situation, it does simplify the complexity of flow problems. The shear stress in a real flow is termed internal friction stress and results from the cohesion and momentum interchange between fluid layers (Daugherty, Franzini and Finnemore, 1985). The magnitude of shear stress is determined by fluid viscosity and the flow velocity gradient which is usually expressed by flow rotationality. For a viscous flow between two infinite flat plates, the relationship between viscosity and shear stress may be expressed as __ du r = (i ---- C Eg 3.1) dy where T = shear stress, y. = the absolute viscosity, and du/dy = the time rate of angular deformation or the velocity y gradient. From the above equation, it is noted that whenever 17 and wherever a flow velocity gradient exits, shear stress occurs correspondingly. Based on laboratory observations, the fluid in direct contact with a solid boundary has the same velocity as the boundary itself. Therefore, when the boundary is fixed, the no-slip condition is applied, i.e., the velocity of the fluid on the solid surface must be zero. As a result, the farther the fluid from the solid wall, the faster the flow velocity can be. This fact explains the reason why the existence of the wall results in flow velocity gradient and shear stress. As shown in Figure 3(a), the cross sectional velocity profile between two plates is parabolic for a laminar flow. This is a typical motion of a low Reynolds number viscous flow. Considering element 1 in Figure 3(a), the uneven velocity gradients in the x- and y-directions cause a rotation. This kind of flow motion is called rotational flow (Janna 1987). It means that the shape and orientation of any fluid element change from time to time and from one position to another as well. When the flow is assumed to be ideal, as shown in Figure 3(b), the flow velocity profile between two plates at any cross section is a straight line. As a result, the shape of a fluid element remains identical. It means that the fluid element does not rotate. This is so called irrotational flow (Janna 1987). Applying the concept of control volume to a two- dimensional, viscous flow moving on the xy-plane, the I 18 'ZA'//////S////////////V//////////;//////JS///////////;//////,///,, 4- ^ 7n7f/f/jj)//7Tn'nnnu/ui/ii////m/wwj/r/T777////M/wvM .< (a) Rotational flow ^ ////////// /// ////////////l/milLL W ri H. 4 (b) irrotational flow Figure 3 Rotational and irrotational flows between two infinite plates \ V \i 19 rotational velocity of a fluid element in the z axis may be expressed as 1 dv ou ez = ---- (------------) (Eq. 3.2) 2 dx 6y in which, z= rotational velocity in the z-direction, u = velocity component in the x-direction, and v = velocity component in the y-direction. The vorticity e used to express the rotational characteristic of flow, is defined as two times the rotational velocity, namely 6v 6u ez = 20z -------------- (Eq- 3.3) 6x 6y in which ez= vorticity in the z-direction. Although Eq's 3.2 and 3.3 are derived from mathematic concepts, they offer a feasible way to measure flow motion. The following section will further explain how to use vorticity to form the governing equations of the motion. Stream Function-Vorticity Equations of Flow Motion The basic governing equations for a two-dimensional incompressible ,viscous flow can be described by the principles of continuity and momentum. Conservation of momentum states: ou ou du dp d2u d2u 20 and, ou 6v 6v op 0 (--- + u --- + v ---) = - + p, ot ox oy 6y The continuity Equation states: (Eq. 3.4) 62v 62v ox2 6y2 (Eq. 3.5) 6u ov ---- + ----- = 0 (Eq. 3.6) 6x 6y where There are three different formulations used to solve Eq's 3.4 through 3.6 for flow problems. They are: 1. The stream function formulation combines all governing equations together and expresses all flow variables by stream function. As a result, a fourth order partial differential equation of stream function needs to be solved. Olson (1975) first applied this approach using finite element method to solve the flow around a circular cylinder. The results, including the separation phenomenon behind the cylinder, were found to be accurate. 2. The velocity-pressure formulation directly solves momentum and continuity equations with the pressure and velocity as unknowns. Kawahara et. al. (1974) used this approach to solve several steady state flow problems with Reynolds number up to 1200. 3. The stream function-vorticity formulation uses 21 stream function,
and pressure. The advantage of this approach is to reduce
approach to solve flows in a constricted channel with
u = ---- (Eq. 3.7)
0
UB ,
6 ( x0)2 6 ( y0)2
(Eq. 3.16)
e = 0 and
To estimate the vorticities generated from the wall,
Due to the occurrence of the maximum horizontal
The coefficients a^ and a2 can be determined by using the
L 1 L Z
where = the unknown nodal values, = shape functions.
D ----r + D ------ Gip + Q = 0 (Eq. 4.11)
, , 6 where
x oxz Y tyz
(Eq. 4.17)
(Eq. 4.21)
- e = ------ +----------------------- (Eq. 4.23)
partial differential equation. To overcome the nonlinearity
(Eq. 4.26)
Then, the element equation is reduced to
1 2 J 4 D b / o
where = the unknown nodal values, and = shape
mi+1 n
cose + d
sine ) ,
along line BC, 6(p/6n = 0,
(Eq. 4.52)
along line AD, tp = 1.
6B
oB |

Full Text |

PAGE 1 FINITE ELEMENT MODELING FOR ENERGY DISSIPATION IN A PRESSURE FLOW by Chwen-geng Guo B.S., National Cheng-Kung University, Taiwan, 1978 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Department of civil Engineering 1988 PAGE 2 This thesis for the Master of Science Degree by Chwen-geng Guo has been approved for the Department of civil Engineering by John R. Mays Ronald c. Scherer b. I PtPtft PAGE 3 Guo, Chwen-geng (M.S. civil Engineering) Finite Element Modeling for Energy Dissipation in a Pressure Flow Thesis directed by Assistant Professor James c. Y. Guo. A finite element numerical model using quadratic isoparametric elements was developed to analyze a twodimensional pressure flow and to predict the flow energy losses associated with the friction and cross sectional contraction. The unsteady flow governing equations were converted into stream function and vorticity diffusion equations. Initial conditions were determined by potential flow and boundary conditions included a uniform incoming flow and non-slip condition on the solid wall. Flow cases studied included the opening ratios ranging from 0.3 to 1 and Reynolds numbers varying from 50 to 200. In general, this numerical approach provides stable computations for flow to reach its steady state when Reynolds number is less than 200 and opening ratio is larger than 0.3. Although for high Reynolds number flows the exit velocity profile may slightly deviate from the exact velocity profile, the predicted friction loss and orifice loss still give good agreements with the theoretical solutions and laboratory data. PAGE 4 ACKNOWLEDGMENTS This project was supported by the Denver Center for the Performing Arts and funded by the Helen G. Bmfils Foundation. I would like to thank Dr. Ronald c. Scherer, senior scientist, Denver Center for the Performing Arts, and Dr. James c. Y. Guo, assistant professor, University of Colorado at Denver, for their advice and encouragement on this project. Also, I wish to express my thanks to thesis committee, Dr. David w. Hubly, and Dr. John R. Mays, for introducing me to the Fluid Mechanics and Finite Element Method. Finally, I would like to give my special thanks to my parents, my wife and children, for their understanding throughout the entire project. PAGE 5 CONTENTS CHAPTER I. INTRODUCTION . . . . . . . . . . . . . . . . . 1 Objective . . . . . . . . . . . . . . . . . . 1 Scope of the study . 2 II. LITERATURE REVIEW 3 Energy Losses in An Incompressible Flow . . . 3 Numerical Modeling of An Incompressible Flow 12 III. GOVERNING EQUATIONS OF FLOW MOTIONS . 16 Flow Motions and Descriptions . . . . . . . . 16 Stream Function-Vorticity Equations of Flow Motion . . 19 Nondimensionalization of the Governing Equations . 23 Boundary conditions . . . . . . . . . . . . . .25 IV. NUMERICAL MODELING BY FINITE ELEMENT SCHEME .. 29 Basic Concepts ....... 29 Derivation of Finite Element Equations of Flow Governing Equations ... 35 Element Shape Function ... 40 Numerical Initial and Boundary Conditions 49 Computational Flow Chart . 53 V. CASE STUDIES AND RESULTS 57 Case One : Study of Flow Patterns 63 PAGE 6 vi Case Two: Prediction of Energy Losses 77 VI. CONCLUSIONS 99 ........ 102 APPENDICES A. THE DERIVATION OF ELEMENT EQUATION BY GAUSS-GREEN THEOREM 105 B. THE INTEGRATION OF SHAPE FUNCTIONS 109 C. DEFINITION OF SYMBOLS 115 D. THE FORTRAN SOURCE CODE 117 PAGE 7 TABLES Table 1. Computational Flow Chart . 55 2. Summary of the Flow Parameter Used in the case studies .......................... 58 3. Computed Velocity Distribution for Opening Ratio= 0.7 and Reynolds number= 50 73 4. Comparison of Peak Vorticity 76 5. Summary of Computational Time Intervals and the Required Times for Flow to Reach Steady State . 78 6. Computed Velocity Distribution in a Poiseuille Flow for Reynolds Number= 50 81 7. Computed Velocity Distribution in a Poiseuille Flow for Reynolds Number= 100 . 81 8. Comparison Between the Computed Exit Velocity Profile at X/B = 8.7 and The Analytical Profile for Poiseuille Flow . .. 84 9. Comparison of Friction Loss Between Computational and Analytical Predictions 10. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile 87 for the cases with a contraction .. 94 11. Distribution of Flow Pressure with a Contraction . 95 12. Comparison of Discharge Coefficients Between the Computed Results and Experimental Data 97 PAGE 8 FIGURES Figure 1. Illustration of Sudden Expansion . 7 2. Illustration of the Flow Domain with Orifice Plate . . . . . 10 3. Rotational and Irrotational Flows between Two Infinite Plates 18 4. Boundary Conditions of Poiseuille Flow . 27 5. Illustration of one-Dimensional Linear Finite Element ........ 30 6. The Example of Eight-Node Quadratic Quadrilateral Elements ..... 34 7. The Eight-Node Quadratic Quadrilateral Isoparametric Element with Natural Coordinates ... 43 B. Finite Element Net Work for the case with Opening Ratio of 1 . 59 9. Finite Element Net Work for the Case with Opening Ratio of 0.7 ....... 60 10. Finite Element Net Work for the Case with Opening Ratio of 0.5 ...... 61 11. Finite Element Net Work for the case with Opening Ratio of 0.3 .... 62 12. The steady State Distribution of Stream Function and Contours of Vorticity for Opening Ratio = 0.7 and Reynolds Number = 50 ....... 64 13. Illustration the Evolution of Flow Patterns Changing From Potential Flow to Its Steady State for Opening Ratio = 0.7 and Reynolds Number= 50 ....... 66 PAGE 9 14. The Variation of Velocity Profile along the Flow Direction for Opening Ratio = 0.7 ix and Reynolds number= 50 ........ 72 15. The Distribution of the Wall Vorticity for Opening Ratio= 0.7 and number= 50 . ." ... 75 16. The steady-state Velocity Profiles In A Poiseuille Flow for Reynolds Number= 50 79 17. The Steady-state Velocity Profiles In A Poiseuille Flow for Reynolds Number= 100 80 18. The Variation of Velocity Profile Along the Flow Direction for Opening Ratio = 1 and Reynolds number= 50 ... 82 19. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile for Poiseuille Flow with X/B = 8.7 83 20. Comparison of Friction Loss Between Computational and Analytical Predictions for Poiseuille Flow with X/B = 8.7 88 21. The Variation of Pressure Along the Flow Direction for Poiseuille Flow with X/B = 8.7 89 22. Comparison Between The Computed Exit Velocity and.The Analytical Profile for Open1ng Rat1o = 0.7 .. 91 23. Comparison Between The Computed Exit Velocity Profile and The Analytical Profile for Opening Ratio = 0 5 9 2 24. Comparison Between The computed Exit Velocity Profile and The Analytical Profile for Opening Ratio= 0.3 93 '25. computed Pressure Drops .. 96 26. comparison of Discharge Coefficients Between Computed Results and Experimental Data 98 27. Coordinate Transformation 111 PAGE 10 CHAPTER I INTRODUCTION Objective Energy dissipation is an important subject to many professions. Examples include crude oil delivered in pipes, water supply systems, air flow through the larynx and blood flow in arteries. Numerous studies have been conducted to develop analytical and empirical equations to estimate energy losses associated with the shear friction and eddy motions. Recently, energy loss in the larynx flow has attracted much attention because of the lack understanding of flows near the glottis. The Denver Center for the Performing Arts, in association with the University of Iowa, has been engaged in a laboratory investigation of this phenomena. Since the energy dissipation in a flow is related to both flow and fluid properties, laboratory arrangements may not easily cover wide ranges of flow variables such as Reynolds number. As a result, numerical modeling is considered a possible helpful alternative. The purpose of this study is to develop a numerical model to analyze a two-dimensional pressure flow and to predict the flow energy losses associated with the wall PAGE 11 friction and cross sectional-contraction. In this work, plate orifices are considered in order to establish the finite element technique before applying the technique to laryngeal flow problems. Scope of the study 2 After a review of numerical methods, the finite element method was selected to convert the nonlinear governing equations of unsteady flow into their numerical forms and to include the complicated boundary conditions in the numerical integrals. Solutions in terms of the stream function, vorticity, velocity and pressure drop at each time step were obtained from the numerical computations for unsteady flow. The steady state solutions were achieved by allowing the successive time-dependent solutions to converge. It is found that for low Reynolds number flows, comparisons between the predicted flow variables and analytical solutions are in good agreement. Although this developed numerical scheme is only for two-dimensional fl.ow, it is believed that the approach is also suitable for a similar three-dimensional flow problem. PAGE 12 CHAPTER TWO LITERATURE REVIEW Energy Losses in an Incompressible Flow The main objective of this study is to develop a numerical model for the prediction of flow patterns and energy losses in a two dimensional pressure flow. In general, there are two kinds of energy losses involved in a pressure flow. They are friction loss and minor loss. Friction in a flow is caused by the non-uniform velocity distribution and fluid viscosity. Friction results in energy dissipation that converts mechanical energy into thermal energy. This energy loss to the environment is irreversible. Friction loss is determined by the length of the conduits, flow Reynolds number, and wall roughness. Using dimensional analysis, the friction loss equation in a pressure flow may be derived as follows (Daugherty, Franzini, and Finn.emore, 1985): ( Eq. 2.1) 2g where, hL= the head loss due to friction, L= the length of the conduit, Rh = hydraulic radius, V= cross sectional average flow velocity, g = acceleration of gravity, and Cf= PAGE 13 4 Darcy friction factor. For a pipe flow, substituting Rh=D/4, in which D is the pipe diameter, into Eq. 2.1, the Darcy-Weisbach equation is obtained: ( Eq. 2. 2) D 2g where D = the pipe diameter, and f = 4Cf. To determine the pipe friction factor f the pipe flow experiments were conducted as early as the mid-1800s. For high Reynolds number flow, flow in rough pipe resulted in a higher energy loss than in smooth pipe. In 1933, J. Nikuradse experimented using artificially roughened pipes and obtained reliable information for friction factor (Gerhart and Gross, 1985). In 1944, the Moody Chart for Friction Factor f versus Reynolds number and relative roughness was published and, since then, it has been widely adopted by engineers. To expand the friction loss, Hagen and Poiseuille successfully derived analytical solutions in 1840 for laminar flows in a circular pipe and between two plates. They stated that for laminar flow in a circular pipe, the loss of head in friction is proportional to the first power of the cross sectional average velocity and is given by 1.1. L ( Eq. 2. 3) where 1.1. = absolute viscosity, and = the density of fluid. By equating Eq. 2.2 and 2.3, we have PAGE 14 5 64 f = ----DV ( Eq. 2. 4) Since is flow Reynolds number, Eq. 2.4 can be further expressed as 64 f =---( Eq. 2. 5) Re .in which Re= flow Reynolds number. From Eq. 2.5, it is found that in laminar flow the friction loss is only a function of Reynolds number and independent of the roughness of the pipe wall. On the Moody chart, the friction factor f for laminar flow is a straight line. This fact is not only true for a laminar pipe flow but also for a laminar flow in a noncircular duct. Jones (1976) expanded Eq. 2.4 to duct flow with the circular diameter D replaced by the hydraulic radius, Dh' and conducted a series of experiments for various duct configurations. He suggested that for the laminar flow in a rectangular duct with width, a, and height, b, the friction factor may be described by k f =---( Eq. 2. 6) 64 where k = (Eq. 2.7) 2 + 11 _E._ (2 __ ) -3-24 a a and Dh = hydraulic diameter. PAGE 15 For the laminar flow between two infinite plates, the value of variable a, in Eq. 2.7, is infinite. As a result, k is equal to 96 when Dh is equal to 4B (Blevins 1984, and Gerhart and Gross i985). Losses due to the change of local pipe configurations are termed minor losses. The flow in a pipe 6 system may need to pass through fittings, bends, or abrupt changes in the pipe configuration. These kinds of changes in local pipe geometries generate a great deal of disorderly motion and result in dissipation of flow momentum. Whenever the flow velocity is altered, eddy currents are induced and a loss of energy in addition to friction loss is created. Although minor losses are usually confined to a short length of flow path, the effects may not disappear for a considerable distance downstream. The minor loss may be expressed as h = K L 2g (Eq. 2.8) where K = loss coefficient, which must be determined experimentally for each situation. The only loss coefficient which can be calculated by simple analytical methods is that the axisymmetric sudden expansion as shown in Figure 1. This problem was first treated by Borda and Carnot. Using continuity, linear momentum, and energy equations and neglecting the shear stress along the wall, the head loss for sudden PAGE 16 7 Figure 1 Illustration of sudden expansion '. PAGE 17 can be expressed as 2g ) (Eq. 2.9) where A 1= the cross sectional area at position 1, A 2= the cross sectional area at position 2, and v1= the average velocity at position 2. Eq. 2.9 with Eq. 2.8 yields 8 K A1 2 =(1---) A2 (Eq. 2.10) Although we may be uneasy about the assumptions of the Borda-carnot approach and feel that Eq. 2.10 is too simple to evaluate the accurate energy loss for sudden expansion, the laboratory experiments showed that the difference between observed and predicted values of K in Eq. 2.10 were within a few percent (Vennard 1969). Based on this fact, it can be assumed that minor losses may be considered for practical purposes to occur in a much more localized flow region. As to the minor losses other than sudden expansion, they all have to be determined by experiments. In 1885, Weisbach followed the approach used by Borda-Carnot to observe the energy losses of abrupt contraction arid a fair estimate of loss coefficients was obtained. In 1952, Gibson's studies on diffusers of conical form grasped the trend of the change of loss coefficient, which are defined as a function of cone angle and area ratio. In PAGE 18 Langhaar successfully calculated the loss of head at entrance by using boundary-layer theory. As to the losses in bends and elbows, there were studied by Beij and Hofmann in 1933. Hofmann gave loss coefficients for wide ranges of bend shape, roughness, and Reynolds number. 9 Nigro (1978) and Pao(1961) indicated that for a flow with Reynolds number below 104, the discharge coefficient becomes increasingly sensitive to Reynolds number, first rising and then falling as Reynolds number decreases. The same phenomenon as indicated by Negro and Pao was observed by Scherer (1981) in his experiments on laryngeal flow for a convergent-divergent glottis and opening ratios d/B less than 0.05, where d = the glottal diameter, and B = the subglottal (tracheal) lateral diameter. Nigro is The discharge coefficient for the pipe flow used by c -( [ 1 ] A 2 4 p ) (Eq. 2.11) where C = discharge coefficient, Ud= the mean flow velocity_ through the orifice throat, a = the cross sectional area of orifice, A = the cross sectional area of the pipe, p = the pressure drop between the flange taps as shown in Fig 2, and = fluid density. Considering a unit width of a two-dimensional flow, Eq. 2.11 can further be simplified to evaluate the flow velocity between two infinite plates. PAGE 19 'lo .. ... ..... .. I I __J i'l-I -Flange tap ...... -., ... _.. ..-: .( 'a,. ...... ,..,. ...... .. ,; -, _,. r.,.r .. ,., ... ,.-_.,. ........... ., _., .... -, ,!,; plate . __ .. :----:---Centerli!le ______ ... __ 2B 2d j_ If :: :-: .. .... -:. .. """11., ... _., .. ' ,, ,. ,. ,, tC' .... ........ ,. --. ;. ;. Solid wall Figure 2 Illustration of the flow domain with orifice plate Outflow ---..:,. _\ .. 1 ) / .. / ......... o PAGE 20 11 c ( [ 1 (_Q_)2 ] 2 .. p ) ( Eq. 2.12) B where d/B = the opening ratio between two plates, d = half of the opening width at the orifice, B = half of the width between the two plates. As indicated by Nigrq, for laminar flow through a thin square-edged orifice with Reynolds number between 100 and 200, the discharge coefficient for flange taps is approximately .70 to .ao. The values of loss coefficients obtained by the pioneer researchers mentioned above were continually modified by researchers. There are presently many handbooks and technical papers published with experimental data (Idel'chek 1960, and ASHRAE Handbook of Fundamentals). Unfortunately, owing to different experiments conducted under different experimental conditions, data values for minor losses may be different for the same flow configuration. Therefore, the most representative one is difficult to determine. Furthermore, experiments for the evaluation of minor losses are common for turbulent flows, whereas the empirical data for laminar flows are less common. These are the main reasons why this numerical study for estimating the energy losses for a laminar flow is needed. In this study, the laminar flow patterns and energy losses associated with the plate orifice as shown in Figure PAGE 21 2, are numerically calculated. Before the selection of numerical methods, the following literature review was performed. Numerical Modeling of An Incompressible Flow There are three approaches by which we can solve fluid problems, namely 1. Experimental 2. Theoretical; and 3. Numerical 12 The experimental approach has the capability of producing the most realistic answers for many flow problems. However, its cost is the highest one among these three methods. In theoretical approaches, assumptions should be made in order to make the problem tractable. Therefore, the problems that can be solved by theoretical approaches are restricted to simple boundary conditions and usually limited to linear problems. However, the shortcomings of experimental and theoretical methods do not exist in numerical approach. In the numerical approach, a limited number of assumptions are made and a digital computer is used to solve the governing equations of a flow problem. There is no restriction to flow linearity and even-complicated geometry can be treated. Furthermore, the computational costs are always much less PAGE 22 13 than that of laboratory equipment for experiments. It is why the numerical method has become more important when the computing facilities become more efficient and accessible. Most people attribute the first definite work of numerical approach to Richardson (1910) who introduced point iterative schemes to numerically solve Laplace's equation. His scheme used data available from the previous iteration to update each value of the unknown and was called the finite difference method (FDM). Richardson actually carried out calculations for the stress distributions in a masonry dam with great success. Up to the late 1930's, FDM was used for solving flow problems. Allen and Southwell (1955) applied relaxation scheme to solve the incompressible, viscous flow around a cylinder. Although they were not the first scientists to solve these kinds of problems numerically, their results did capture the flow phenomenon and had good accuracy. At the same time, John von Neumann developed his method for evaluating the stability of numerical methods for solving time-dependent problems. Douglas and Rachford (1950) developed an implicit method for solving parabolic and elliptic partial differential equations. From then on, the basic knowledge for solving fluid dynamics problems by FDM was complete. Before the finite element method (FEM) was fully developed, FDM was used to solve flow problems. The PAGE 23 advantage of FEM is that its rate of convergence is much faster than that of FDM. 14 The finite element method was originally used by aircraft structural engineers in the 1950's to analyze a large system of structural elements in the aircraft. Usually, there are two ways to formulate element approximate equations: one is the variational method, which ha.s a close relationship to the classical variational concepts of the Rayleigh-Ritz method 1909); the other one is the weighted residual method, which is modeled after the well-known method of Galerkin (1915). In weighted residual method, the error in the approximate satisfaction of the governing equations is not set to zero, but instead its integral with respect to the selected 11weights11 is required to vanish. The application of FEM to nonstructural problems such as fluid flows was initiated by Zienkiewicz and Cheung (1965). Pioneering research has employed weighted residual methods to formulate approximate element equations. Oden (1972) was among the first to derive the basic theoretical analog for the Navier-Stokes equations. Baker (1971, 1973) published his results for the analysis of incompressible flow using the Galerkin method. Olson (1974) used a pseudovariational finite-element algorithm for solving the biharmonic equation for stream functions. Lynn(1974) utilized the least-squares methods to develop a finite- PAGE 24 15 element algorithm for laminar boundary-layer flow. In 1978, Backer successfully extended his studies from 2-dimensional laminar flow to 3-dimensional turbulent flow. From that time on, FEM was widely adopted for solving fluid dynamics problems. In this study, the numerical model for solving laminar, incompressible flow will be derived using FEM because of its flexibility of element shape and efficient computation. Galerkin's method is used to weight residuals in the approximate governing equations of flow. The procedure of formulation will be discussed in Chapter Four. PAGE 25 CHAPTER THREE GOVERNING EQUATIONS OF FLOW MOTIONS Flow Motions and Descriptions Using viscosity as a criterion, fluid may be classified as real fluid and ideal fluid. For a-real fluid, flow motion results in shear forces between fluid layers. But for an ideal flow, the shear force is assumed to be zero everywhere. Although this is an idealized situation, it does simplify the complexity of flow problems. The shear stress in a real flow is termed internal friction stress and results from the cohesion and momentum interchange between fluid layers (Daugherty, Franzini and Finnemore, 1985). The magnitude of shear stress is determined by fluid viscosity and the flow velocity gradient which is usually expressed by flow rotationality. For a viscous flow between two infinite flat plates, the between viscosity and shear stress may be expressed as du r = (Eq. 3.1) dy shear stress, = the absolute viscosity, and dujdy = the time rate of angular deformation or the velocity I gradient. From the above equation, it is noted that whenever PAGE 26 17 and wherever a flow velocity gradient exits, shear stress occurs correspondingly. Based on laboratory observations, the fluid in direct contact with a solid boundary has the same velocity as the boundary itself. Therefore, when the boundary is fixed, the no-slip condition is applied, i.e., the velocity of the fluid on the solid surface must be zero. As a result, the farther the fluid from the solid wall, the faster the flow velocity can be. This fact explains the reason why the existence of the wall results in flow velocity gradient and shear stress. As shown in Figure 3(a), the cross sectional velocity profile between two plates is parabolic for a laminar flow. This is a typical motion of a low Reynolds number viscous flow. Considering element 1 in Figure 3(a), the uneven velocity gradients in the x-and y-directions cause a rotation. This kind of flow motion is called rotational flow (Janna 1987). It means that the shape and orientation of any fluid element change from time to time and from one position to another as well. When the flow is assumed to be ideal, as shown in Figure 3(b), the flow velocity profile between two plates at any cross section is a straight line. As a result, the shape of a fluid element remains identical. It means that the fluid element does not rotate. This is so called irrotational flow (Janna 1987). Applying the concept of control volume to a twodimensional, viscous flow moving on the xy-plane, the PAGE 27 (a) Rotational flow I .(b) flow .-Figure 3 Rotational and irrotational flows between two infinite plates 18 0 .............. i i L. S.. :{ -..,,,' PAGE 28 19 rotational velocity of a fluid element in the z axis may be expressed as 1 ov ou az = ( ---) (Eq. 3.2) 2 ox oy in which, e = rotational velocity in the z-direction, u z velocity component in the x-direction, and v = velocity component in the y-direction. The vorticity used to express the rotational characteristic of flow, is defined as two times the rotational velocity, namely ov ou = = 29 = z z ( Eq. 3. 3) ox oy in which z= vorticity in the z-direction. Although Eq's 3.2 and 3.3 are derived from mathematic concepts; they offer a feasible way to measure flow motion. The following section will further explain how to use vorticity to form the governing equations of the motion. Stream Function-Vorticity Eguations of Flow Motion The basic governing equations for a two-dimensional incompressible ,viscous flow can be described by the principles of continuity and momentum. Conservation of momentum states: ou ou ou op o 2 u o 2 u (-+ u + v -) = ---+ IJ. ( + ot ox oy ox ox2 oy2 ) PAGE 29 20 (Eq. 3.4) and, ou ov ov op o 2 v o 2 v (-+ u + v -) = --+ ( + ) ot ox oy oy ox2 oy2 (Eq. 3.5) The continuity Equation states: ou OV --+ = 0 ( Eq. 3. 6) ox oy where = the fluid density, and = the viscosity. There are three different formulations used to solve Eq's 3.4 through 3.6 for flow problems. They are: 1. The stream function formulation combines all governing equations together and expresses all flow variables by stream function. As a result, a fourth order partial differential equation of stream function needs to be solved. Olson (1975) first applied this approach using finite element method to solve the flow around a circular cylinder. The results, including the separation phenomenon behind the cylinder, were found to be accurate. 2. The velocity-pressure formulation directly solves momentum and continuity equations with the pressure and velocity as unknowns. Kawahara et. al. (1974) used this approach to solve several steady state flow problems with Reynolds number up to 1200. 3. The stream function-vorticity formulation uses PAGE 30 stream function, and vorticity, e, to replace velocity and pressure. The advantage of this approach is to reduce unknowns from three variables, u, v, and p, to two 21 variables, and e. Cheng (1972) successfully applied this approach to solve flows in a constricted channel with Reynolds numbers less than 100. Using the formulation, Eq's 3.4 and 3.5 are reduced to a diffusion equation of vorticity and Eq. 3.6 becomes a poisson equation of the stream function. In this study, the combination of and e is used to describe flow. The following section further presents the detailed derivations. From the principle of conservation of mass, the velocity components u and v are related to the stream function by the following equations, u = (Eq. oy v = -----(Eq. 3.8) ox Substituting Eq's 3.7 and 3.8 into Eq. 3.3 yields 2 = e (Eq. 3.9) in which + For an irrotational flow, Eq. 3.9 is further simplified to (Eq. 3.10) Taking the derivative ojoy of Eq. 3.1 and the PAGE 31 derivative -ojox of Eq. 3.2 and adding the resultant equations together, we obtain or < + 0 ) = + v ot ox oy 22 .,.2 = (Eq. 3.11) Dt Eq. 3.11 is a two-dimensional diffusion equation of vorticity. Eq's 3.9 and 3.11 represent two coupled second order equations governing incompressible viscous flows. The continuity equation is implicitly satisfied by the relationship between stream function and velocity. In order to predict energy loss, we need to calculate the pressure distribution through out the flow field. Taking the derivative of Eq. 3.4 with respect to x and the derivative of Eq. 3.5 with respect to y and then summing the resultant equations together, we have .,.2p = ou 2 OV 2 ou ov [ ( ) + ( ) + 2 ----] ox oy oy ox (Eq. 3.12) Substituting Eq. 3.6 into Eq. ;3.12 yields ou ov ou ov ----] (Eq. 3.13) OX oy oy ox Eq. 3.13 may further be expressed by stream function, namely PAGE 32 23 2 + -( ) ] (Eq. 3.14) Eq's 3.9, 3.11, and 3.14 are governing equations used in this study. It is noted that Eq's 3.9 and 3.14 take the same form as the Poisson equation provided that the terms on the right-handed side of each equation can be as constants. Nondimensionalization of the Governing Equations In fluid mechanics, the solution of flow problems may become more generalized when they are expressed in nondimensional forms. Without nondimensionalization, solutions of a numerical model may only be suitable for the specific conditions studied. In the process of nondimensionalization, certain flow dominating parameters like Reynolds number may be found. To perform nondimensionalization, some reference flow,quantities have to be chosen as characteristics of the flow system (Daily and Harleman, 1966). For a flow, the inflow velocity, u, may be chosen as the characteristic velocity and a half of the distance between two infinite plates, B, may be chosen as the characteristic length. A set of dimensionless flow quantities may be defined as follows: X y B B PAGE 33 24 0 u 0 v u = v = u u to tu 0 p = p = B u2 cp 0= cp 0 B = (Eq. 3.15) UB u Substitution of Eq. 3.15 into Eq. 3.9, the equation of stream function becomes 02( 0 ) 02( 0 ) cp cp + = 0 ( x0)2 0 ( y0)2 or 'f' 2 ( cpo ) = (Eq. 3.16) similarly, Eq 3.11 becomes + u + v = ( + ) o(yo) ua o(xo)2 o(yo)2 (Eq. 3.17) The term of in the above equation is the reciprocal of flow Reynolds number which delineates flow regimes, i.e., laminar or turbulent. UB Re = ---(Eq. 3.18) Reynolds number physically represents the ratio of inertia force of flow to viscous force. For a viscous flow between infinite plates, gravity does not affect the flow pattern. PAGE 34 25 It is also obvious that capillarity and surface tension are of no practical importance. Hence the significant forces affecting flow motions are inertia force and viscous force. Substituting Eq. 3.18 into Eq. 3.16 yields 1 (Eq. 3.19) Re Based on the similarity theorem, Eq's 3.17 and 3.19 shall depict the same flow dynamic behavior as long as the boundary conditions and Reynolds number remain the same. The dimensionless form of the pressure distribution equation is derived as follows: = -2 (-----------------) or o(u0 ) o(v0 ) o(u0l o(v0 ) = -2 ( ) (Eq. 3.20) o(yo) o(x0 ) o(x0 ) o(yo) Comparing Eq. 3.20 with Eq.3.14, it can be seen that the density of fluid is eliminated by nondimensionalization. In the fol.lowing sections, the flow problem will be solved using the nondimensional equations, and the superscript O's are neglected. Boundary Conditions When solving a partial differential equation, the boundary conditions play a determining role. In the study PAGE 35 26 of the flow between two infinite plates, as shown in Figure 2, flow is symmetric to its centerline, therefore, only a half of the flow field is used to define the computational domain. The inflow is assumed to have a uniform velocity profile with a constant horizontal velocity u. The boundary conditions associated with each governing equation of flow are shown in Figure 4 and summarized below : On the boundary AB : Flow enters with a uniform velocity profile, therefore, we have v = 0 u = 1 = uy = Y = 0 and the value of pressure is set to be a selected datum value such as 100. on the boundary BC Applying the non-slip condition to this boundary, we have = 0, and. v = u = 0 To estimate the vorticities generated from the wall, it is necessary to apply the Taylor expansion to Eq. 3.3. The detailed derivation will be discussed later. On the centerline AD : Based on the continuity and Eq. 3.15, the stream PAGE 36 p = Selected datum value .-. u=l centerline D : Y/B v=O cp=l' =0 ---u=l I v=O cp=y =0 c ...... .-"-_._.; ,.::-__:-:. ::. 6> ,: :-.. "'!-',:; .. ;, ._.._ :,,.-,.'--' 'J, . :.;._..._,....,,. -.--..--;t, I I : / __ ___. . -:;:--;::..-;.. r. Solid wall u=v=O cp=O f=fw Figure 4 Boundary conditions of Poiseuille flow X/B N -.J PAGE 37 function at the centerline must be = 1 Due to the occurrence of the maximum horizontal velocity, we can write 28 oujoy = o, (Eq. 3.21) Considering the symmetry of the flow about the centerline, we conclude that v = 0, and f = o, on this boundary. On the boundary DC This is the exit of the computational domain. Cheng (1972) suggested to introduce a velocity profile defined by a cubic equation of distance measured from the wall to calculate stream functions the line DC, to ensure that Poiseuille flow is generated. This suggestion forces the velocity profile on the line DC to agree with the analytic velocity profile, i.e., a parabolic curve. In this study, it is found that as long as the computational domain is long enough, the exit velocity profile agrees with the analytic solution quite well. Therefore, the condition that the vertical velocity vanishes on the boundary DC is adopted and used in the computation. PAGE 38 CHAPTER FOUR NUMERICAL MODELING BY FINITE ELEMENT.SCHEME Basic Concepts The finite element method is a numerical scheme and computational procedure for obtaining approximate solutions to many engineering problems (Segerlind, 1984). In general, the application of the method takes the following five steps: 1. Discretization and Selection of Element Configuration This step involves dividing the entire computational domain into a suitable number of "small" bodies called finite elements. The shapes of the elements may be triangular or rectangular, which is determined by the user. 2. Determination of the Approximatiqn Equation For each element, it is required to have an approximation equation for solutions in terms of the unknown nodal values. For example, a one-dimensional linear element AB, as shown in Figure 5, has an approximation equation expressed by a linear equation. The coefficients a1 and a2 can be determined by using the nodal conditions, namely cp = 1 cp = 2 at at X= X 1 X= X 2 (Eq. 4. 2) PAGE 39 cp = al. + a 2 x Tt2 tl 1 2 rxl L x2 Figure 5 Illustration of one-dimensional linear finite element 30 I X PAGE 40 Substituting Eq. X -x lp = L 1 4.2 into Eq. x-x + L 2 4.1 gives { Eq. 4. 3) Eq. 4.3 is a typical form for approximation equations used in the finite element method. The nodal values are 31 multiplied by linear functions of x, which are called shape functions or interpolation functions. To be more Eq. 4.3 may be rewritten as cp = N. l. l. ( i = 1 to n ) {Eq. 4.4) .. where = the unknown nodal values, Ni = shape functions. In the above linear case, the shape functions become X -x 2 N =--1 L and x-x 1 N =--2 L (Eq. 4.5) The subscript i in Eq. 4.4 is an index that indicates how many nodes are involved in an element. For a linear triangular element, the value of n equals 3. As to a quadratic quadrilateral element, the value of n equals a. 3. Derivation of Element Equations There are several ways to derive the element equations in the finite element method such as the variational method, and the weighted residual method. In this study, Galerkin's method is used, which is a kind of weighted residual method. The weighted residual method is used to minimize the residual resulting from.an approximate solution that is substituted into the differential equations governing a PAGE 41 32 field problem (Desai,l979). Usuc;llly, it requires that L I Wi(x)R(x) dx = o 0 ( Eq. 4. 6) where R(x) =the residual, Wi(x) =the ith nodal weighting function, and L = the length of integral. There are a number of schemes developed for the weighted residual method to determine the weighting functions, among which are collocation, subdomain, least squares and Galerkin's methods. When Galerkin's method is used, Eq. 4.6 is evaluated using the same functions for Wi(x) as were used in the approximate solution 1984). Thus, for any element, the weightihg functions for deriving element equation are exactly the same as its shape functions. When there are n nodes involved in an element, there are n shape fqnctions that can be used as weighting functions to form n equations as Eq. 4.6 for solving n unknown nodal values 'n For a two-dimensional field problem, Eq 4.6 may be rewritten in its matrix form I [N ]T R(x) dA = 0 A (i = 1 to n) (Eq. 4. 7) where [N]T = row vector containing the element shape functions. After integration, we have n simultaneous equations with n unknowns of ti. These simultaneous equations may be expressed as PAGE 42 33 = {F} (Eq. 4.8) where [K] = stiffness matrix or element property matrix, = the unknown nodal values, and {F} = force matrix or vector of nodal forces! which may be a zero vector. 4. Assembly of Element Equations to Obtain Global Equations and Introduction of Boundary Conditions The final aim is to obtain the global equations that approximately define the behavior of the entire system. As shown in Figure 6, when the eight-node quadratic quadrilateral elements are used ( for instance, the element one has nodes 1, 10, 15, 16, 17, 11, 3, and 2), the element equation can be written as where the superscripts represent the element number, and The sequence of node numbers takes the counterclockwise rule. Since the different elements have different nodes, we have to, based on the node numbers, incorporate all element equations into the global After the boundary conditions are introduced into the computation, the modified global equation similar to Eq. 4.8 may further be obtained. PAGE 43 4 3 2 1 (2) (6) 11 25 ... ( 1) (5) 10 15 24 Figure 6 The example of eight-node quadrilateral -elements l 34 .. PAGE 44 5. Computation and Solution Finding of the Global Equation The global equations are a set of linear or nonlinear simultaneous algebraic equations. They can be solved using the well-known Gaussian elimination method to obtain the value at each node. 35 In the following sections, a numerical model of twodimensional incompressible flow will be derived according to the steps stated above. Derivation of Finite Element Equations of Flow Governing Eguations To formulate element equations, a general form of the dimensional field equation can be used (Segerlind 1984). 02'P 02'P D --+ D Gcp + Q = 0 X OX2 y oy2 (Eq. 4.11) Comparing Eq. 4.11 with Eq. 3.9, we conclude that for the calculation of stream function, Dx = Dy = 1, G = o, and Q = e According to the Galerkin's method, the residual integral equation for any element e is given by ( ) I T o2m o2m {R e } = -[N] (Dx + D Gcp + Q)dA ax oy (Eq. 4 .12) where [N] is the row vector containing the element shape functions and is regarded as the weighting function in the Galerkin's method. To integrate the second-derivative terms in Eq. PAGE 45 36 4.12, the Green-Gauss theorem may be applied. After integration, Eq. 4.12 becomes where I = r [N]T( D ---cosa + D sine )dr x ox Y oy (Eq. 4.14) {k(e)} IA< ox o[N]T o[NJ o[N]T o[NJ = + Dy )dA OX ox oy oy + IAG[N]T[N] dA (Eq. 4.15) and (Eq. 4.16) As to two-dimensional time dependent problems, the general form of their governing equations may be expressed as o2n on Dx --=.-L + D --"' + Q = -"'-ox2 Y oy2 ot (Eq. 4.17) Eq. 4.17 can be rewritten as follows: + Q = ( Eq. 4.18) where = D --+ D --x ox2 Y oy2 The residual integral is PAGE 46 37 similarly, after integration, we have = 0 ( Eq. 4. 20) where {I(e)}, {k(e)}, and {f(e)} are the same as expressed in Eq. 4.13, and {t(e)} =the nodal values of The transpose {i(e)}T of {t(e)} is {t(e)}T = [ at n at ( Eq. 4. 21) ] (Eq. 4.22) The detailed derivation of Eq's 4.13 and 4.20 are presented in Appendix A. Next we apply the general forms to the flow problem. The governing equations for a two-dimensional, incompressible flow as derived in chapter two are: ., 2 = e: ( Eq. 3 .16) 1 .,2 De: e: = (Eq. 3 .19) Re Dt 2 ., p = 2 [ ou OV au ov -----] (Eq. 3.20) ox oy oy ox in which Eq's 3.16 and 3.19 are the two coupled governing equations of flow motion. Eq. 3.19 may be further expressed by stream function as follows: PAGE 47 38 1 +----(Eq. 4.23) Re ot oy ox ox oy Because both and e are unknowns, Eq. 4.23 is a nonlinear partial differential equation. To overcome the nonlinearity of Eq. 4.23, some researchers (Cheng, 1972) suggested that an unsteady flow problem can be regarded as linear in stream function and vorticity e at each time step. We may consider two solutions ( en ) at the nth time step and ( en+ 1 ) at a time increment t later. The steady-state solutions are achieved by allowing the successive timedependent solutions to converge. Therefore, the governing equations of flow motion become 2 =e n (Eq. 4. 24) and 1 ,2 oen OE'n E'n+1 = + Re ot oy ox OE'n ox oy (Eq. 4.25) It is noted that, at each time step, the and en in Eq. 4.25 are no longer unknowns. Thus Eq. 4.25 becomes solvable numerically. After the steady state is reached, the velocity components u and v at every node can be obtained by taking the first derivative of the stream function with respect to y and x .,respectiveiy. Then thevelocity g+adients in the right-hand side of Eq. can further be computed. As a result, Eq, 3.20 may be regarded as a Poisson equation. PAGE 48 39 If the shape functions used for calculating vorticities are the same as those used for stream functions, applying Eq. 4.13 to Eq's 4.24 and 3.20, and Eq. 4.20 to Eq. 4.25, the element equation can be obtained as where {I(e)} = {k(e)} = {k(e)} p (Eq. 4.26) ( Eq. 4. 27) ( Eq. 4. 28) --case + sine ) dr ox oy ( Eq. 4. 29) oEn+l OEn+l --'------'-case + sine ) dr ox oy (Eq. 4.30) op op ---case + sine )dr ox + / oy o[NJ oy ( Eq 4. 31) )dA ( Eq. 4. 32) (Eq. 4.33) ( Eq. 4. 34) i .. ( .. PAGE 49 40 [C(e)] = IA [N]T[N] dA (Eq. 4.35) {f(e)} = I ([N]T n) dA (Eq. 4.36) 'P .A {f(e)} I [N] T ( 6 'Pn+l On On )dA = oy ox ox oy A (Eq. 4.37) {f(e)} -2 I [N]T au ov au ov = (-------) dA p ox oy oy ox (Eq. 4.38) Eq's 4.26 through 4.38 are the element equations for solving a two-dimensional, incompressible flow problem. For an irrotational flow, the governing equation is y 2 'P = 0 (Eq. 3.10) Then, the element equation is reduced to ( Eq. 4. 26) where {I(e)} and [k(e)] can be referred to Eq1s 4.29 and 4.32 respectively. In the sections, the element shape functions will further be discussed and determined. Element Shape Function To apply the finite element method to a flow problem, we need to decide what.kind of elements and element equations should be used. The elements may be triangular, rectangular, and other shapes. The element equations may be PAGE 50 !'i1 linear, quadratic, or cubic. There are two requirements that should be satisfied in the determination of element shape function. The first requirement is the continuity between adjacent elements. In Eq. 4.6, R(x) represents the governing equation. suppose that R(x) contains derivatives up to the r-th order. To avoid the integral becoming infinite, the compatibility requirement (Huebner, 1975) requires that its (r-1)-th order derivative must be continuous. This ensures that the piecewise continuity is maintained at element interfaces. In this study, the linear element shape function may be chosen. ;Because the governing equations of flow motion contain the second-derivative terms, a linear element with continuity of first order can sufficiently satisfy this requirement. The second requirement, criterion of completeness, results from the approximation theorem ( Zienkiewicz, 1977). finite element method has an implicit assumption that convergence occurs as the size of each element decreases. To ensure that this assumption is satisfied well, at least a constant value can be obtained in a local element when the size of any element tends to zero. Therefore, if a derivative of order r exists in the governing equation, it is necessary for the element shape function to be at least of the order r so that, in the limi.t, such a constant value can be acquired. Based on this criterion, a quadratic PAGE 51 42 element shape function was chosen. It is noted that the criterion of completeness is not as necessary condition as the requirement of compatibility, for determining the order of element shape function. But, in this study, since a time-dependent equation of vorticity (Eq. 3.19) needs to be dealt with, it is better to use a quadratic element shape funqtion to ensure an efficient convergence. In addition, considering the application of this numerical model to a flow problem with curved boundary, the quadratic quadrilateral isoparametric elements were selected and Qsed in this study. Isoparametric elements are useful in modeling a field problem with the curved boundary and in grading a mesh from coarse to fine ( Cook, 1981). They make it possible to use nonrectangular quadrilateral elements for solving a field problem with the curved boundary. The interpolation equation for the eight-node quadratic quadrilateral isoparametric element shown in Figure 7 is = a1 + a2 a + a3 B + a4 aB + a5 a 2 + a6s2 + a 7 a 2 B + a8aA2 (Eq. 4. 39) with nodal conditions : a = -1 and B = -1 at node 1, a= 0 and B = -1 at node 2, a = 1 and B = -1 at node 3, a= 1 and B = 0 at node 4, a= 1 and B = 1 at node 5, PAGE 52 a J Figure 7 The eight-node quadratic quadrilateral r isoparametric element with natural coordinates 43 PAGE 53 44 Q = 0 and B = 1 at node 6, Q = -1 and B = 1 at node 7, and Q = -1 and B = 0 at node a. (Eq. 4.40) in which a and B are natural coordinates. The element shape functions can be obtained by applying their basic characteristic to Eq. 4.39. The basic characteristic of the shape functions is that each shape function has a value of one at its own node and zero at the .other node. Applying this characteristic to Eq. 4.39 gives a standard form as follows: cp = N.t. 1 1 ( i = 1 to a ) where t. = the unknown nodal values, 1 functions, in which N1 = N2 = N3 = \(1+a)(1-B)(a-B-1) N4 = N5 = \{1+a)(1+B)(a+B-1) N6 = N 7 = \(1-a)(1+B)(-a+B-1) N 8 = (Eq. 4.4) and N. 1 = shape (Eq. 4.41) Applying Eq 4.41 to Eq's 4.29 through 4.38, the integration can be performed. There are two types of integrations: 1. The integrations occur on the boundary, such as PAGE 54 45 Eq's 4.29, 4.30, and 4,31. This kind of derivative boundary condition will be discussed in section 4.4. 2. The integrations occur within the elements. They can further be subdivided into two different types: the integrations of the shape functions such as Eq's 4.35, and 4.36, and the containing partial derivative terms such as Eq's 4.32, 4.33, and 4.34, that involve the first derivatives of the shape functions, and Eq's 4.37 and 4.38, that include the first derivatives of stream functions, vorticities, and velocities. For the purpose of integration, the change of variables is required. For a two-dimensional integral, the change of variables equation is I I A 1 1 f{x,y)dA =I I f[x{a,B),y{a,B)] I det[J] I dadS -1 -1 (Eq. 4. 42 )" in which f(x,y) = the function expressed by global coordinates, f[x{a,B),y{a,B)] =the function expressed by natural coordinates, and [J] = Jacobian matrix, which is defined as [J] = o(x,y) o(a,B) = ox oa ox 68 {Eq. 4.43) and det[J] = the determinant of the Jacobian matrix. Applying Eq. 4.43 to Eq. 4.35 yields PAGE 55 46 (Eq. 4.44) Similarly, Eq. 4.36 becomes 1 {f(:)} = JJ_1([N(a,B)]T[N(a,B)] {e }) ldet[JJI dadS ( Eq. 4. 45) Eq's 4.44 and 4.45 can be evaluated by Newton-Cotes or Gauss-Legendre quadrature method. For those integrals containing the partial derivatives of the shape functions, we need to change the local derivatives in terms of a and B to global derivatives. The relationships between local derivatives and global derivatives are oN. oNi 1 (a,B) (a,B) ox -1 oa = [J(a,B)] (Eq. 4.46) oN. oN. 1 (a,B) l (a,B) oy os where [J] =Jacobian matrix, and =inverse of [J]. However, in Eq. 4.46 an explicit form of [J]-1 cannot be obtained (Zienkiewicz,1977). Numerical schemes may be used to evaluate element integrals. In this study, GaussLegendre quadrature is adopted. The Gauss-Legendre quadrature selects the sampling points between the integral intervals to achieve the best accuracy. For double it states that PAGE 56 47 R = J l Jl m -l 1f(a,B) dadB (Eq. 4. 47) where m is the number of the selected sampling points with respect to B, which can be determined by equating (2m-1) to the highest power of B existing in the function f(a,B). If the resulting m is not an integer, the next larger integer is used. n is the number of the selected sampling points with respect to a, which can be determined the same way as that used for determining m. wi is the weighting factor for the sampling points with respect to B, and Wj is the weighting factor for the sampling points with respect to a. The detailed derivation of Eq's 4.42 and 4.46 and the further explanations of Gauss-Legendre quadrature are presented in Appendix B. How many sampling points should be chosen when using Gauss-Legendre quadrature to maintain the ? In this study, referred to Eq. 4.44, the highest power of variables. appearing in [N]T[N] are a 4 and s4 Therefore, the value of 3 is used for both m and n. It means that the order of quadrature is 3 and there are 3x3 po1nts used for computation. Cook (1974) suggested for plane elements, either linear or quadratic, the 2x2 points are usually the best, except for a markedly nonrectangular or elongated PAGE 57 48 quadratic element, where 3x3 points may give better accuracy. In this study, the numerical model developed may further be applied to the analysis of laryngeal flow. Considering the irregularity of larynx geometries, the nonrectangular and elongated elements were expected. Therefore, 3x3 points are considered the most reasonable choice. Applying Eq. 4.47 to Eq. 4.44, it gives [C(e)] 3 3 )]T[N(a = I: I: {W W [N(a ,B ,B )]} 1 det[JJ 1 i=l j=l i j j i j i (Eq. 4.48) The local coordinates of each sampling point and its corresponding weighting factor in the apove equation are presented in Appendix B. After the coordinates of each sampling point were determined, det[J] at each point can be obtained using Eq. B.3. Then Eq. 4.48 can further be evaluated numerically. The same procedure is applied to Eq. 4.45. As to Eq 4.32, which contains the first derivatives of the shape functions, with the aid of Eq's 4.46 and.4.47, it becomes {k (e)} 3 3 o[N(aj,Bi)J o[N(a.,B.)J = I: I: w. w. ( J ]. cp i=l j=l ]. J ox ox o[N(a.,B.)] o[N(aj,Bi)] + J ]. ) 1 det[JJ I -oy oy (Eq. 4.49) Similar approach to the forementioned can be applied to Eq's PAGE 58 49 4.33 and 4.34. The term of in Eq 4.37, may be evaluated by applying Eq. 4.46 to the center-point of each element. For any element, equals a constant. By applying Gauss-Legendre quadrature, Eq. 4.37 becomes {f(e)} 3 3 T = :E :E (w.w. [N(aj,Bi)] } 1 det[JJ I Cc i=l j=l 1 J (Eq. 4.50) 0 0 where Cc = ( n n ) evaluated at oy OX ox oy the center of each element. Eq. 4.38 can be solved by the same way as used for ( ou ov ou ov Eq. 4.37. The term of ox oy ox ), 1n Eq. 4.37 are also evaluated at the center of each element. Numerical Initial and Boundary Conditions Since Eq. 3.19 is a time-dependent equation, the initial conditions have to be determined. In this study, the initial flow field was considered to be a potential flow (ideal flow) with a uniform inflow of velocity u at the entrance. As shown in Figure 4, the initial vorticity is zero everywhere in the computational domain. To determine the distribution of stream function, the governing equations for the irrotational flow were used. There are two kinds of boundary conditions involved in this study. They are, Dirichlet boundary conditions, .. PAGE 59 50 that are the values of specified on the boundaries, and Neumann boundary conditions, that are the values of derivatives and specified on the boundaries. Recall that in the derivation of element equations, the boundary conditions have already included as a part of element equations (Eq's 4. 29, 4. 30, 4. 31) and take the following form ( D ---cose + D Y sine ) x ox oy (Eq. 4.51) where e = the angle to the outward normal to the boundaries of each element. If Dx = Dy = D, then Eq. 4.51 becomes D = 0 on where is the derivative normal to the boundary. To evaluate Eq. 4.29, in which theNeumann boundary conditions of the stream functions were involved, the derivative boundary conditions of stream functions in Figure 4 are listed below: along line AB, = o, and e = 1ao0 along line BC, = o, along line AD, = o, and e = 90, along line DC, = o, and e = oo, Substituting Eq. 4.52 into Eq. 4.29 yields {I (e) } = 0 (Eq. 4. 52) (Eq. 4. 53) Eq. 4.53 implies that the Neumann boundary conditions for PAGE 60 51 stream functions are self-consistent in the element equations. If boundaries AB and DC in Figure 4 are vertical straight lines, the contribution of the derivative boundary conditions of the stream functions to the element equations automatically vanishes. Therefore, for the stream function, only Dirichlet boundary conditions need to be considered. As a result, Eq. 4.29 can be removed from the element equations. The Dirichlet boundary conditions of the stream functions as shown inFigure 4 are listed below: Along line BC, = o, and along line AD, = l How to impose Dirichlet boundary conditions by modifying the global equations can be found in Reference 7, 9, and 29. To consider the numerical boundary conditions of vorticity, a linear combination of E and its normal derivative is required to be specified throughout the entire computational domain (Backer, 1983). On the fixed wall, BC, the value of vorticity is a priori unknown. To evaluate.the wall vorticity, the no-slip boundary condition alone is insufficient. From Eq. 4.24, at a point (x0,y0 ) on the wall, the vorticity may be by (Eq. 4. 54) where B is the coordinate normal to the wall. To evaluate PAGE 61 52 joB2 we can apply Taylor series expansion to a point (x1,y1 ) which is located a small distance from the wall along the B direction. ( Eq. 4. 55) since the no-slip condition indicates that OB ( xo Yo ) = o we have 2 en(xo,Yo) = 8 2 (Eq. 4 . 56) Wall vorticity is then expressed in terms of the stream functions on the wall and the node at a small distance from the wall. Because B is the coordinate normal to the wall, the construction of the grids for computing wall vorticity do need special attention. J As to the derivative boundary condition of vorticity, the values of oen+l/ox along boundaries AB, BC, CD and oen+l/oy along boundary BC are undeterminable. Eq. 4.30 may be considered as a negligible term which is further ignored from the global equations . In this study, it is shown that neglecting the derivative boundary conditions of vorticity can still give acceptable results which agree to analytic solution well. The same situation may be applied PAGE 62 53 to the derivative boundary conditions of pressure. By adding case times (Eq.3.4) to sine times (Eq.3.5) and using the appropriate boundary conditions, we have ( Eq. 4. 57) Substituting Eq. 4.57 into Eq. 4.31, it is found that the integrations of the second derivative of velocity along the element boundaries are difficult to evaluate. Like Eq. 4.30, Eq. 4.31 is considered as a numerically negligible term and ignored from the global equation. Computational Flow Chart To develop a digital computer program for analyzing a two-dimensional, incompressible, viscous flow between two infinite plates, the initial conditions are determined using the governing equations for ideal flow and then Eq's 4.26 and 4.27 are used to calculate stream functions and vorticities for successive time steps. Neglecting the derivative boundary conditions, Eq's 4.2 and 4.27 are reduced to [k(e)] cp {cp(e)} = n+1 {f(e)} cp ( Eq. 4. 58) [k(e)]{f:(e)} + [c(e)]{;(e)} = {f(e)} n+1 n (Eq. 4.59). in which d fn+1 n {f:n } = { dt } = { .&t } (Eq. 4.60) PAGE 63 54 Substituting Eq. 4.60 into Eq. 4.59, the recurrence relation is (Eq. 4.61) At each time step, Eq. 4.58 is first used to calculate Substituting new values of to Eq. 4.58, the corresponding wall vorticities can be obtained. Using Eq. 4.61, the vorticity can be updated for the next computational time step. The criterion for determining the steady state was defined to be whem the largest value of OE/Ot in the entire computational domain was less than 0.04 (Backer, 1983). After the steady state is reached, using Eq. 4.28, the corresponding pressure distribution is further computed. Comparing pressure drops between downstream and upstream flows, we can determine the energy dissipation. A flow chart illustrated the computational procedures is presented in Table 1. The source code and FORTRAN program can be found in Appendix D. PAGE 64 FLOW CHART A I Define Flow Field 1 Discretize and Select the Element configuration Specify the Approximation Equation The Isoparametric Quadratic Quadrilateral Elements are Used Set Up Initial condition .,2 cp= 0 The Boundary Conditions for Initial State Are Shown in Figure 4. (Continued) Table 1 computational Flow Chart. 55 PAGE 65 56 Flow Chart B t = 0 .,.2 cp = 0 and = 0 l J t = t + &t J 'l l Determination of wall vorticity .. 2 ... f = -(cp cp ] , .... w 2 w w+l B l 1 .,.2 D Re = Dt l .... 2 cp = l Check steadiness no 0/0t < 0.04 I yes Calculation of Pressure "2P 2 [ ou ov ou OV ] = CSX6Y-6Y6X I End I PAGE 66 CHAPTER FIVE CASE STUDIES AND RESULTS There are four different contraction geometries used in this study. They can be distinguished by the opening ratio, r, which is defined as r = d B (Eq. 5.1) where d = a.half width of the opening, and B =a half width of the distance between two plates. When r = 1, flow does not have any contraction. Results of these cases present predictions of friction losses only. When r < 1, flow is subject to minor losses due to contraction. The opening ratios used in this study were 0.3; 0.5, 0.7, and 1.0. Table 2 summarizes the total number of nodal elements, and Reynolds numbers used in the studies of these four cases. The corresponding element network for each case is shown in Figures 8 through 11. Since the abrupt changes of stream function and vorticity are expected to occur near the orifice, those cases with a smaller opening should use a larger number of smalier elements around the orifice. Thus, the number of nodal points and elements around the orifice used in the case with an opening ratio of 0.3 has as nearly 1.5 times nodal points and elements as used in the case with an opening ratio of 0.7. The total length of the computational domain also PAGE 67 Table 2 Summary of the Flow Parameters used in the case Studies. Number of Number Opening Reynolds nodal of Length Separation ratio number points element L/B eddy ==--=================--=====--==================== 1.0 50 413 120 8.7 No 1.0. 75 413 120 8.7 No 1.0 100 413 120 8.7 No 1.0 150 413 120 8.7 No 1.0 200 413 120 8.7 No 0.7 50 235 64 6.5 Yes o. 7 100 235 64 6.5 Yes .7 200 235 64 6.5 Yes 0.5 50 309 86 8.7 Yes 0.5 100. 309 86 8.7 Yes 0.5 20" 0 309 86 8.7 Yes 0.3 5 0 305 84 9.3 Yes 0.3 100 305 84 9.3 Yes 0.3 150 305 84 9 3 Yes ====--== ... U1 Q) PAGE 68 FINITE (0.0, 1.0) y/9 (0.0 ,0.0)" ELEMENT GRID FOR PAGE 69 FINITE ELEMENT GRID NET'w'DRK FOR R=0.7 -. PAGE 70 ._ .. .l. :: .. . FINITE ELEMENT GRID NETV/DRK FOR R=0.5 (CENTERLINE OF (0.0, 1.0) y/B 1111 DNJ( (0.0 ,0.0) -ORIFICE .. (SOLID VIALL) x/B. Figure 10 Finite element grid net work for the case with opening ratio of o.s \ \ (8.7, 1.0) (8.7, 0.0)" --_j ..r.;,'-'! ) l 0'1 .... ..... \:; ,., ',. ...., ___ ,.. PAGE 71 FINITE ELEMENT GRID NET\JORK FOR 'R=0.3 co.o, 1.0) CCENTERLINE OF FLO"'w'). .y/B II II ll if (0.0 ,0.0) .Orifice (SOLID 'w'ALL) Figure 11 Finite element grid net work for the case with opening_ratio of o.J x/B. (9.3, 1.0) <9.3, o.o> 0\ N /-,-. i ... PAGE 72 63 needs special attention. After the flow passing the orifice, it takes a longer distance for the flow with more contraction to return to its parabolic velocity profile than that with a smaller contraction. When the computational domain downstream of the orifice is not long enough, the numerical procedure may become unstable and results are not reliable. In this study, it was found that the range of computational domain, in terms of L/B varied from 6.5 to 9.3. Predictions using these computational domains give good agreements to the laboratory and analytical results. Predictions and comparisons for each case study are presented as follow: Case One : Study of Flow Patterns As may be expected, eddies are generated behind the orifice plate. Table 2 indicates that as long as the opening ratio is less than unity, eddies always occur. For high Reynolds number flows, eddies immediately downstream of the orifice plate may further shed out. Figure 12 presents steady state flow with an opening ratio of 0.7 and Reynolds number of 50. The computational time interval, At, was 0.01. It took 640 time steps, i.e. t 0 = 6.4, for flow to reach the steady state. The upper half of Figure 12 shows the distribution of stream function and the lower half presents the contours of vorticity. It is noted that a reverse flow occurs immediately downstream of the PAGE 73 L------:----====()-------------J --cp=. -20--------------_J o.S Figure 12 'The steady state distribution of stream function and contours of vorticity for Opening Ratio = 0.7 and Reynolds Number = 50 .... .. / 0\ PAGE 74 orifice plate and the highest vorticity is produced near the orifice. On this figure, the standing eddies are clearly shown. The general flow patterns predicted in this study are very similar to that reported by White (1974), Chung (1978), and, Kwon and Pletcher (1981), although the numerical approaches were different. 65 Evolutions of flow patterns with respect to time changing from potential to steady state flow are shown in Figure 13. It can be seen that the eddies are growing in size as time goes on. Using the steady state criterion mentioned in Chapter four, flow becomes steady when t 0 = 6.4. However, it is found that when t 0 = 3.2, the variations of streamline patterns have already been numerically negligible. This may imply that the criterion for steady state suggested by Baker (1983) seems to be conservative. Due to the inadequate accuracy of software for plotting streamlines and the coarseness of elements used, some flow patterns in Figure 13 are not shown to be separating from the edge of the orifice. Figure 14 displays the plot of velocity profiles along the flow direction at X/B = o., 1.1, 2.2, 2.5, 2.9, 4.7, and 6.5. Near the entrance, the velocity profile shows that the flow has not been fully developed. Behind the. orifice, negative velocity occurs. This means that eddy motions in the area of wake produce a reverse flow. Table 3 PAGE 75 .:..: .. (. .. t o.o 1 Y/B -:::::::: 61 I I I I I I I I I I I I I 0 1 2 3 Lt. 5 X/B i 0 ... t 0.2 1 Y/B i :::::=::::: 6$. I I I I I ' i fJ 1 2 3 Lt 5 X/B to 0.4 1 } Y/B t I I I I I I I I I I 0 1 't 5 X/B to 0.6 .. 1 Y/B (continued> Figure 13 Illustration the evolution of flow patterns changing from potential flow to its steady state for Opening Ratio 0.7 and Reynolds Number= 50 ) --66 I I 6 I G ' ... I 6 PAGE 76 67 ______ _.::,_: J ,0 . . : 1. 2 . : : _. 3 .. 't . 5 X/B 6 .. .. --.--.-. -. I I : 'l-c=================================::::::' to 1 o -.. .. ,. Y/B 2 X/8 ::-.. 0 t l. 2 l 1 3 X/B .. tL-----:-___ ---::===:---------.0 I I I I I I I I 1 .. 1 II I I I 1 . 2 I . 3 't 5 X/B ... ---. PAGE 77 I. 0 t 1.6 . . ... I( 1 2 3 5 I I I G t0 1.a l'-c==========:======================== y/B 0 t 2.0 Y/B 1 2 3 Y/B 2 3" ___ 5 ._X/B G PAGE 78 69 0 t 2.4 ..:;;::'. y/B f:---:::pr l \-2._ -.... .. ,,, 0 1 2 3 tt 5 X/B G 0 t II 2.6 I' I I -1., I e 1 2 3 tt 5 x;a G t0 11 2.a r?... 1 . Y/B C--.--:------:---======-------. .!: .. 0 1 2 3 lt. 5 X/B I I I I : . G 0 t II 3o0 .. 1-c==================================== Y/B I 1 = I I I I I I I .. ---... -.. ---. .:. -----. -----------. 5 X/B. G PAGE 79 70 0 t 3.2 Y/B :r::: -:-I I I I I If I I I I I 0 1 2 3 ... '+ to 3 .4 I I I I I I I 5 X;B 6 . 1--c:=============================== Y/B 6 t::----?f&:_ -, 1''' 'I 0 1 2 3 '+ 0 t 3.6 I I I I I I 5 X/B I G Y/B .-r 1 2 3 Lf 5 6 0 .t 3.8 1-c=================================== Y/B 6 I 0 a 1 2 3 PAGE 80 D 1 2 3 o. t 4.8 I I I .lt I ' 5 xis I I 6 71' y/B. '. 1 2 3 0 : t .. 5.2 : I j I 4 I I I I I 5 6 .X/B Y/B i---------:-: 0 I 0 1 2 3 4 0 t .4 I I I I 5 X/B I I I G Y/B I I I I I I i I I I I I I 0 1 2 3 4 5 G X/B ."' 0 ,' j I : .. PAGE 81 Y/B _f//// 4u=l 2 Figure 14 The variation of velocity profile along the flow direction for Opening Ratio = 0.7 and Reynolds number = 50 X/B ....,J N PAGE 82 73 Table 3 Computed Velocity Distribution For opening ratio= 0.7, and Reynolds number= 50 ============================================================ o.oo X/B 1.10 2.50 2.90 4.70 6.50 2.20 ============================================================ Distance Velocity Distribution At orifice yjB ujUo ujUo ============================================================ 0.00 1.00 0.00 0.00 0.00 0.00 o.oo 0.05 1.00 0.18 -0.09 -o.o8 0.08 0.14 0.10 1.00 0.36 -0.17 -0.13 0.15 0.27 0.20 1.00 0.65 0.03 0.05 0.36 0.54 0.30 1.00 0.92 0.33 0.28 0.58 0.77 o.oo 0.45 1.00 1.11 1.03 0.92 0.96 1.06 1.27 0.60 1.00 1.27 1.67 1.61 1.37 1.29 1.50 0.80 1.00 1.30 1.74 1.80 1.58 1.45 1.57 1.00 1.00 1.33 1.80 1.96 1.77 1.57 1.64 ============================================================ ----------------------------------------------------------------------------Exit Profile Analytical Profile ====================================== Distance Y/B Velocity Profile ujUc U/UC ====================================== o.oo 0.05 0.10 0.20 0.30 0.45 0.60 0.80 1.00 o.oo 0.09 0.18 0.34 0.49 0.68 0.82 0.92 1.00 0.00 0.10 0.19 0.36 0.51 0.70 0.84 0.96 1.00 Note: X= horizontal distance. y= distance from solid wall. B= half width between two plates. Uo= entrance velocity. Uc= centerline velocity PAGE 83 74 presents velocity profiles at different locations and a comparison between the exit velocity profile and the analytic parabolic velocity profile. It is noticed that although the existence of the orifice does disturb velocity profile, flow still returns to Poiseuille profile far downstream. In this study, the velocity profile for Poiseuille flow can be expressed as (Gerhart and Gross, 1895) U/Uc = 1 -( 1-y2 ) where u= velocity in x-direction, Uc= u at centerline, y= distance from the solid wall. The wall vorticity is also worthy of studying because the wall shear stress is directly proportional to its value. Figure 15 depicts the distribution of the wall vorticity for the flow with an opening ratio of 0.7 and Reynolds number of 50. It shows that the peak wall vorticity occurs at the contraction section and the wall shear stress changes its direction corresponding to the flow direction near the wall. Table 4 presents the values of peak vorticity for various cases under different Reynolds numbers. It is noted that for the with an opening ratio = 1, the peak vorticity, highest shear stress, occurs at the entrance because the flow experiences the steepest velocity gradient at the beginning of the entrance length and its value is proportional to the Reynolds numbers. On the contrary, for those cases with an opening ratio less PAGE 84 22 20 18 16 14 .a 12 ] 1:! 10 0 8 6 4 2 0 -2 0 2 4 6 X/B Figure 15. The distribution of the wall vorticity for opening Ratio= 0.7 and Reynolds number 50 Ul PAGE 85 Table 4 Comparison of Peak Vorticity. ============================================================= The Values Opening. Reynolds of Peak occuring ratio number Vorticity Position ============================================================= 1.0 50 16.62 Near Entrance 1.0 75 17.04 Near Entrance 1.0 100 17.36 Near Entrance 1.0 150 17.84 Near Entrance 1.0 200 18.19 Near Entrance 0.7 50 20.78 Near Orifice 0.7 100 23.49 Near Orifice 0.7 200 26.43 Near Orifice 0.5 50 34.15 Near Orifice 0.5 100 .38. 98 Near Orifice 0.5 200 42.40 Near Orifice 0.3 50 60.89 Near Orifice 0.3 100 66.80 --. Near_Orifice ---o-:3-150 73.77 Near Orifice ============================================================= -.J O'o PAGE 86 77 than one, the highest shear stress occurs at the orifice and their values are a function of both Reynolds number and the opening ratio. With the same Reynolds number, the value of peak vorticity with an opening ratio of 0.3 is as nearly three times as that with an opening ratio of 0.7. Computational time intervals and the required times for flow to reach the steady state are presented in Table 5. It is noted that the higher the Reynolds number is, the more time steps are required. Case Two : Prediction of Energy Losses Energy dissipation in a pressure flow includes friction loss and minor loss. In this portion, flows are studied with and without an orifice. When the opening ratio is unity, the energy loss is essentially a friction Figures 16 and 17 present the steady-state velocity profiles at various locations with Reynolds number equal to 50, and 100 respectively. The corresponding data can.be found in Tables 6 and 7. Flow starts with a uniform velocity profile at the entrance anq gradually changes to a fully developed flow through the entrance length. Figure 18 depicts the variations of velocity profiles along the flow direction with Reynolds number of 50. Figure 19 and Table 8 show comparison between the exit velocity profiles under Reynolds. numbers equal to 50, 75, 100, 150, and 200 with the analytic velocity profile of a fully developed flow. It is noted that the exit velocity profiles with PAGE 87 Table 5 summary of Computational Time Intervals and the Required Times for Flow to Reach steady State. ========================================================== Time Interval Opening. Reynolds Used for ratio number Computation Time Steps for Reaching Steady State ========================================================== 1.0 50 0.0100 960 1.0 75 0.0100 1200 1.0 100 0.0100 1200 1.0 150 0.0100 1900 1.0 200 0.0050 2400 0.7 50 0.0200 320 0.7 100 0.0100 780 0.7 200 0.0050 2015 0.5 50 0.0100 633 0.5 100 0.0100 754 0.5 200 0.0050 1992 0.3 50 0.0100 567 0.3 100 0.0100 7 .37 0.3 150 0.0025 3376 '.J (%) PAGE 88 0 :J :J ,... u 0 cu > 1.5 1.4 -, _m .:t 2.0 0 0 v ).9 ----+ 1.3 1.2 1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 A 6.3 X s. 7 0.2 0.4 Distance y /B 0.6 0.8 Figure 16 The steady-state velocity profiles ln a Poiseuilie flow for Reynolds Number 50 1.0 -.J \D PAGE 89 2 1.9 1.8 1.7 ..6 1.5 1.4 1.3 0 1.2 :::> 1.1 :::J p 1 'ij 0.9 0 Ill 0.8 > 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0.00 X/B o o + 0 3.9 .1::.. 6. 3 v a. 1 I I--0.20 0.40 Distance y/B 0.60 0.80 Figure 17 The steady-state velocity profiles in a Poiseuille flow for Reynolds Number.= 100 .. 1.00 CXI 0 PAGE 90 81 Table 6 Computed Velocity Distribution in a Poiseuille Flow Reynolds Number= so. ============================================================ X/B o.oo 2.00 3.90 6.30 8.70 Exit Analytical Profile Profile ============================================================ Distance Velocity distribution Velocity Profile yjB U/UO u;uc U/UC ==================.========================================== o.oo 1.00 o.oo o.oo o.oo o.oo 0.00 0.00 0.05 1.00 0.18 0.16 0.15 0.15 0.15 0.98 0.10 1.00 0.35 0.32 0.30 0.30 0.30 0.19 0.15 1.00 0.49 0.45 0.43 0.42 0.42 0.28 0.20 1.00 0.64 0.59 0.57 0.56 0.56 0.36 0.40 1.00 1.07 1.01 0.99 0.98 0.98 0.64 0.60-1.00 1.27 1.27 1.27 1.27 1.27 0.84 0.80 1.00 1.33 1.39 1.43 1.44 1.44"0.96 1.00 1.00 1. 33 1.43 1.48 1.49 1.49 1.00 ============================================================ Table 7 Computed Velocity Distribution in a Poiseuille Flow Reynolds Number= 100 ============================================================ X/B Exit Analytical 0.00 2.00 3.90 6.30 8.70 Profile Profile ============================================================ Distance Velocity Distribution Velocity Profile Y/B U/UO U/UC ujUc ============================================================ 0.00 1.00 0.00 o.oo 0.00 0.00 o.oo o.oo 0.05 1.00 0. 21-0.18 0.17 0.16 0.11 0.01 0.10 1.00 0.40 0.36 0.33 0.32 0.23 0.19 0.15 1.00 0.56 0.50 0.47 0.46 0.32 0.28 0.20 1.00 0.74 0.66 0.61 0.60 0.42 0.36 0.40 1.00 1.13 1.08 1.03 1.02 0.72 0.64 0.60 1.00 1.24 1.27 1.27 1.27 0.89 0.84 0.80 1.00 1. 24 1.32 1.37 1.39 0.98 0.96 1.00 1.00 1.24 1.32 1.39 1.42 1.00 1.00 ============================================================ Note: Uo = entrance veloctiy. Uc= centerline velocity. X horizontal distance. y = distance from the plate. B = half width between two piates. PAGE 91 _y)s 1 g u=1 ... 0 .. I I I I I I _1 I w f "1'' h' U.I.IU.I.Ii1.111.1/.l.l 'I 11//// .1/r.I.IIN F /N F/1'11" f.l /.'/.II' .I /NF'f/'-I'H".IUF 2 4-6 Figure 18 The variation of velocity profile along the flow direction for Opening Ratio = 1 and Reynolds number = 50 I I J I r r .. ,, I """"" 8 .. Q) liJ PAGE 92 2 1.9 1.8 1.7 1.6 1.5 1.4 1.3 u 1.2 :::> 1.1 ::J .0 1 u 0.9 0 Ill 0.8 > 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 .. 0.00 Re 0 50 + 75 <> 100 )(. 150 /:). 200 ."'\1 Analytical 0.20 0.40 0.60. 0.80 Distance y/B Figure 19 -Comparison between the computed exit velocity profile and the analytical profile for Poiseuille flow with X/B = 8.7 .(--_ 1.00 co .w :... PAGE 93 Table 8 Comparison between the Computed Exit Velocity Profile at X/B = 8.7 and the Analytical Velocity Profile for Poiseuille Flow ============================================================ Reynolds Number 50 75 100 150 200 Analytical Velocity Profile ============================================================ Distance Exit Velocity at X/B=8.7 Y/B u;uc U/UC ============================================================ o.oo 0.00 o.oo 0.00 0.00 o.oo o.oo 0.05 0.10 0.11 0.11 0.13 0.14 0.01 0.10 0.20 0.21 0.23 0.24 0.27 0.19 0.15 0.28 0.29 0.32 0.35 0.38 0.28 0.20 0.38 0.38 0.42 0.45 0.51 0.36 0.40 0.66 0.67 0.72 0.76 0.82 0.64 0.60 0.85 0.86 0.89 0.92 0.96 0.84 0.80 0.97 0.96 0.98 0.98 0.99 0.96 1.00 1.00 1.00 1.00 1.00 1.00 1.00 ============================================================ Note: Uo = entrance veloctiy. Uc= centerline velocity. X = horizontal distance. y = distance from the plate. B = half width between two plates. Q) PAGE 94 85 Reynolds number less than 100 give a better agreement to the exact solution. The discrepancy between the predicted exit velocity profiles and analytic solution for high Reynolds number flows may be due to the length of computational flow domain. When it is hot long enough, flow may not become fully developed at the exit. Theoretically, an infinite distance is required for a flow with uniform entrance velocity to become fully developed, but it has been proved both by theory and observation that for a pipe flow, the maximum velocity at the centerline of the pipe may reach 99 percent of its ultimate value in the distance determined as follows: L = 0.058 RePo (Eq. 5.2) where Rep is the Reynolds number for pipe flow and defined as UD IJ. in which U =average flow velocity, and D =the diameterof the pipe (Daugherty, Franzini and Finnemore). Applying Eq 2.7 to Eq. 5.2 and using Reynolds numper determined by hydraulic radius we have L = .348 Re B (Eq. 5.3) For Re = 100, the required length to reach 99% of fully developed flow is 35 B. This is impractical to adopt this computational domain because of the elongation of PAGE 95 computational time. After many tests, the length of 8.7 B is used in this study based on the judgement of the agreement between the predicted exit velocity profile and the exact velocity profile for flows with Reynolds number less than 100. 86 As to the prediction of friction loss, computational domain of L/B=8.7 can still give fair accurate predictions even for the flows with higher Reynolds number. Numerical predictions of friction loss expressed by flow pressure drop are presented in Table 9 and Figure 20. The analytical pressure drop in Table 9 was computed by Eq 5.3 that can be derived by applying Eq. 2.7 to Eq. 2.1 and using hydraulic radius for Reynolds number. 3 L 4 p = ( Eq. 5. 3) Re B where L/B = 8.7 in this study. It can be seen that the computed friction loss is in a good agreement to the analytical solutions. Figure 21 presents the variations of flow pressure along the flow direction with various Reynolds numbers. It shows that independent of Reynolds number, kinetic effects dominate over the first 30 percent of the entrance length. Thereafter, viscous effects gradually increase in importance and become the sole influence at the end of the duct, where the pressure gradient is a constant. To predict the energy losses due to the contraction, PAGE 96 Table 9 Comparison of Friction Loss between Computational and Analytical Predictions ================================================ Reynold Number Computed Analytical Pressure Drop Pressure Drop 50 75 100 150 200 0.50 0.36 0.25 0.18 0.14 0.52 0.35 0.26 0.17 0.13 ================================================ Q) '-.J PAGE 97 0.55 0.5 0.45 0.4 Ill '1 _3 0.35 c 0 :;::; 0 0.3 .:: LL. 0.25 0.2 0.15 0.1 50 70 90 110 130 150 Reynolds Number [] Computed + Analytical .Figure 20 Comparison of friction loss between comput"ational and analytical predictions for Poiseuille flow with X/B = s.7 170 190 (X) (X) PAGE 98 p uo """2 QJ Ill Ill QJ 0.8 ..... ----' ------' ,, ------------------., ,____ ----,_. ...... .......... ----------------------------!: .. 0. 6 0 2 '+' 6 Distance X/B Figure 21. The variation of pressure along the flow directiort for Poiseuille flow with X/B = 8 Re 200 150 100 75 50 .> (X) \0 PAGE 99 90 opening ratios of 0.3, 0.5 and 0.7 are used. The exit velocity profiles for different cases are shown in Figures 22, 23, and 24 and Table 10. It seems that the computational domain was not long enough for the flow with r = 0.3. Table 11 presents the computed pressure drop. Figure 25 plots the computed pressure drops for the flows with the opening ratios of 1, 0.7, 0.5, and 0.3. Using Eq. 2.11, the orifice pressure drop can be converted to discharge coefficient which can be compared with the experimental data observed by Nigro (1978) for pipe flow with a circular orifice. It is noted that the orifice pressure drop measured by Nigro is the pressure difference between two flange taps that were located at one inch on both sides of circular orifice. In this the nearest nodal points at one inch on both sides of the plate orifice were chosen to calculate orifice pressure drops. Comparison between the predicted pressure drop and experimental data is shown in Table 12 and Figure 26. Although energy dissipation of a plate orifice may be somewhat different from a circular orifice, these two sets of pressure drops capture a similar trend and even numerical deviations are fairly small. PAGE 100 0.9 0.8 0.7 0 0.6 :::> :J >. 0 5 u 0 Ill 0.4 > 0.3 0.2 0.1 0 0 0.2 0.4 0.6 0.8 Distance y/B 0 Re=50 + Re=100 <> Re=200 1::. Analytical Figure 2"2 comparison between the computed exit. velocity profile and the analytical Profile for \D -Opening Ratio = 0.7 -. . PAGE 101 1 0.9 0.8 0.7 0 0.6 ::::> "-.. ::J >. 0.5 u 0 (I) 0.4 > 0.3 0.2 0.1 0 0 0 Re=SO I 0.2 0.4 0.6 0.8 Distance y /8 + Re=100 <> Re=200 h. Analytical Figure 23 between the computed exit velocity profile and the analytical profile opening Ratio = 0.5 1 U) N PAGE 102 0.9 0.8 0.7 0 0.6 :J ......... :J >. 0.5 u 0 Gl 0.4 > 0.3 0.2 0.1 0 ,. 0 0 Re=50 0.2 0.4 Distance y /B + Re=100 <> Re=150 0.6 0.8 A Analytical \D w PAGE 103 94 Table. 10 Comparison between the Computed Velocity Profiles and the Analytical Profile for cases with Contraction. =========================================================== Exit Velocity Profile ujUc Analytical d/B = 0.7 d/B = 0.5 ujUc =========================================================== Distance Reynolds number yjB 50 100 200 50 100 200 =========================================================== o.oo o.oo 0.00 o.oo 0.00 o.oo o.oo o.oo 0.05 0.09 0.08 o.o8 0.09 0.09 0.12 0.10 0.10 0.17 0.17 0.16 0.19 0.19 0.23 0.19 0.20 0.34 0.33 0.32 0.37 0.37 0.40 0.36 0.30 0.49 0.48 0.48 0.52 0.52 0.51 0.51 0.45 0.68 0.66 0.66 0.69 0.69 0.69 0.70 0.60 0.82 0.84 0.84 0.84 0.84 0.84 0.84 0.80 0.92 0.93 0.94 0.96 0.96 0.95 0.96 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 =========================================================== ================================================== Exit Velocity Profile u;uc d/B = 0.3 Analytical Profile ujUc ================================================== Distance Reynolds Number yjB 50 100 150 ================================================== o.oo o.oo 0.00 0.00 0.00 0.05 0.12 0.11 0.08 o.1o 0.10 0.28 0.24 0.22 0.19 0.25 0.52 0.44 0.44 0.44 0.40 0.70 0.57 0.60 0.64 0.70 0.87 0.85 0.84 0.91 1.00 1.00 1.00 1.00 1.00 ================================================== Note: Uo = entrance veloctiy. Uc= centerline velocity. X = horizontal distance. y = distance from the B = half width between two plates. PAGE 104 Table 11 Pressure Drop of the Flow with a contraction. ============================================== b/B Ratio Reynolds Number 0 0.3 0.5 0.7 ============================================== 50 75 100 150 200 0.50 0.36 0.25 0.18 0.14 10.31 9.75 10.03 --2.94 0.98 2.74 0.84 2.89 0.81 ============================================== Note: --= not computed d/B = opening ratio between two plates computed in this study \0 ur PAGE 105 p uo a. e Q "QI Ill, Ill Gl 10 5 d/8 = O.J d/8 = o.s )( )( d/8 1:1: 0.7 1 d/8 = 1 ---Reynolds Number Figure Computed pressure drops U) 0\ PAGE 106 Table 12 Comparison of Discharge Coefficients Between the Computed Results and Experimental Data. ============================================================ d/B Ratio d/D Ratio Reynolds Number 0.3 0.5 0.7 0.3 0.5 0.7 ============================================================ 50 100 150 200 0.70 0.72 0.71 0.71 0.74 0.72 0.73 0.78 0.80 0.70 0.71 0.70 0.69 0.71 0.73 0.72 0.71 0.74 0.79 0.81 0.82 ============================================================ Note: --= not computed d/B = opening ratio between two plates computed in this study. d/D = opening ratio in a circular pipe observed in lab. \0 -.,J PAGE 107 0.9 ..., ; 0.8 ..... u ..... Ql 0 u CD 00 as .c: o. 7 ..... Q X A + .. d/B 0.3 o.s 0.7 experimental data d/8 .. 0.7 d/8 .. 0.5 d/B 0.3 10 20 40" 60 100 200 400 600 Reynolds Number Figure 26 Comparison of.discharge coefficients between computed results and experimental \0 co PAGE 108 CHAPTER SIX CONCLUSIONS A finite element numerical model using quadratic isoparametric elements was developed to analyze a Poiseuille flow with and without a contraction. In this study, the unsteady flow governing were converted into stream function and vorticity transport equations. A general numerical form of the Poisson equation was developed to integrate finite element equations. Initial conditions were determined by potential flow and boundary conditions included a uniform incoming flow and non-slip condition on the plate. The range of computational dimensionless time increment was from 0.0025 to 0.01. Flow cases studied include the opening ratios ranging from 0.3 to 1 and Reynolds number varying from 50 to 200. In general, this numerical approach provides stable computations for flow to reach its steady when Reynolds number is less than 200. It was found that the computational domain of L/B =. 8.7 gives, in general, reasonably good predictions compared with analytical solutions. Although for high Reynolds number flows the exit velocity profile may slightly deviate from the exact velocity profile, the predicted friction loss and PAGE 109 orifice loss still give good agreements to theoretical solutions and laboratory data. From the comparisons and discussions above, this numerical model is capable of computing flow patterns and predicting energy losses for a two dimensional flow with Reynolds number less than 200 and orifice opening ratio 100 larger than 0.3. This range may cover most blood flows in the artery and oil movements in lubrication devices. To extend its application to high Reynolds number flows, some improvements are recommended as follows: 1. Efficiency of Computation: In this study, for the case with the opening ratio of 0.7 and Reynolds number of 50, it took 4 hours of computer time to reach the steady state. Requirement of such a longcomputation time may be due to the use of quadratic isoparametric elements. If the linear triangular elements were used, it might reduce the computer time to one half of that used in this study. In addition, the values of wall vorticities play an importance role in the determination of numerical convergence. The. value of the wall .vorticity, as shown in Eq. 4.54, is proportional to the reciprocal to B 2 At each time step, the wall vorticity has to be updated first. If B is too large, the value of the wall vorticity becomes too small. It decelerates the rate of convergence; as a result, it requires more computer time to reach the steady state. Conversely, if B is too small, the numerical procedure may PAGE 110 101 become unstable due to the abrupt increment in the values of the wall vorticity. Therefore, how to select the size of wall elements and determine the value of 8 are subjects for the future studies. 2. Expansion of Its Usefulness: It has been found that when Reynolds number is greater than 200 or the opening ratio is less than 0.3, this numerical model becomes unstable and the steady state can never be reached. As indicated in the beginning of this study, the development of this model was aimed at the application of the finite element scheme to the numerical modeling of larynx flows. How to extend the application of this numerical model from lower Reynolds number flows to turbulent flows may present another challenge to the future study. PAGE 111 REFERENCE 1. Anderson, D. A., Tannehill, J. c., and Pletcher, R. H., Computational Fluid Mechanics and Heat Transfer, McGrawHill, N.Y., 19a4. : Chaps 1, and 6. 2. Baker, A. J., Finite Element Computational Fluid Mechanics, Hemisphere Publishing, N.Y., 19a3. : Chaps 1, 2, 3, 4, and 5. 3 .. Baker, A. J. ''Finite Element Solution Algorithm for Incompressible Fluid Dynamics" in Finite Elements in Fluids 2. Ed. Gallagher, Oden, Taylor, and Zienkiewicz. John Wiley, London, 1975 : 67:a2. 4. Bathe, K. J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, N.J, 19a2. : Chap 5. 5. Blevins, R. D., Applied Fluid Dynamics Handbook, Nostrand Reinhold, N.Y., i9a4. : Chaps 5, 6, and 7. 6. Cheng, R. T., "Numerical Solution of Navier-stokes Equations by Finite Element Method" The Phy. of Fluids 12, Dec. 1972. : 209a:2105. 7. Chung, T. J., Finite Element Analysis in Fluid Dynamics, McGraw-Hill, 197a. : Chaps 1, 2, 3, 4, and 5. a. Connor, J. J. and Brebbia, c. A., Finite Element Techniques for Fluid Flow, Newnes-Butterworths, LOndon, 1976. : Chap 9. 9. Cook, R. D., concepts And Applications of Finite Elenient. Analysis, John Wiley and Sons, N.Y., 19a1. : Chap 5. 10. Daily, J.W. and D.R.F., Flyid Dynamics, Addison-Wesley, Reading, Mass., 1966 : 137:143. 11. Daugherty, R. L., Franzini, J. B., and Finnemore, Fluid Mechanics with Engineering Applications, McGraw-Hill, N.Y., 19a5. : Chaps 1, 7, a, and 12. Desai, c. s., Elementary Finite Element Method, Prentice-Hall, Englewood N.J., 1979. :Chap 1. 13. Fox, R. w. and McDonald, A. T., Intro. to Fluid Mechanics, John Wiley and Sons, N.Y., 19a5 :Chaps 6, and a. PAGE 112 14. Gerald, c. F., Applied Numerical Analysis, AddisonWesley, Reading, Mass., 1970 : Chaps 7, and 9. 103 15. Gerhart, P. M. and Gross, R. J., Fundamentals of Fluid Mechanics, Addison-Wesley, Reading, Mass., 1985. : Chap 7. 16. Huebner, K. H., The Finite Element Method for Engineers, John Wiley and Sons, N.Y., 1975. : Chap 9. 17. Idel'chek, I. E., Handbook of Hydraulics ResistanceCoefficients of Local Resistance and of Friction, u.s National Technical Information Service, 1960. 18. Janna, W. s., Intra. to Fluid Mechanics, PWS Publishers, Boston, Mass., 1987 : 554:557. 19. Jones, o. c., "An Improvement in the Calculation of Turbulent Friction in Rectangular Ducts" ASME J. of Fluids Eng., June 1976. 20. Kawahara, M., Yoshimura, N., and Nakagawa, K., "Analysis of steady Incompressible Viscous Flow" in Finite Element Methods in Flow Problems. Ed. Oden, Zienkiewicz, Gallagher, and Taylor. The University of Alabama Press, Huntsville, Ala, 1974 : 107:120. 21. Kersten, R. D., Engineering Differential Systems, McGraw-Hill, N.Y., 1969 : Chaps 3, and 7. 22. Kreyszig, E., Advanced Engineering Mathematics, John Wiley and Sons, N.Y., 1983 : 412:447. 23. Kwon, o. K. and Pletchcher, R. H. "Prediction of the Incompressible Separated Boundary Layers Including Viscous-Inviscid Interaction" J. Fluids Eng. 101,1979.: 466 :472. 24. Longwell, P. A., Mechanics of Fluid Flow, McGraw-Hill, N.Y., 1980 : Chap. 25. Nigro, F. E. B., et al., 11A numerical study of the Laminar Viscous Incompressible Flow Through a Pipe Orifice" J. of Fluid Eng. 100, 1978. : 467:472. 26. Olson, M.D. "Variational Finite Element Methoc;is for TwoDimensional and Axisymmetric Navier-Stokes Equations" in Finite Elements in Fluids 1. Ed. Gallagher, Oden, Taylor, and Zienkiewicz. John Wiley, London, 1975 : 57:72. PAGE 113 27. Pao, R. H. F. Fluid Mechanics, John-Wiley and Sons, N.Y., 1961. 28. Pin Tong, "On the Solution of the Navier-Stokes Equations in Two-Dimensional and Axial-Symmetric Problems" in Finite Element Methods in Flow Problems. Ed. Oden, Zienkiewicz, Gallagher, and Taylor. The University of Alabama Press, Huntsville, Ala, 1974 57:66. 104 29. Segerlind, L. J., Applied Finite Element Analysis, John Wiley and Sons, Canada, 1984. 30. Scherer, R., Laryngeal Fluid Mechanics:Steady Flow Considerations Using Static Models, University of Iowa, diss., 1981. 31. Streeter, v. L., Fluid Dynamics, McGraw-Hill, N.Y.,1948 : Chap 1, 2, 10, 11, and 12. 32. Vennard, J. K., "One-Dimensional Flow" in Handbook of Fluid Dynamics. Ed. Streeter. McGraw-Hill, N.Y., 1969 : Chap 3. 33. White, F. M., Viscous Fluid Flow, McGraw-Hill, N.Y., 1974 : Chap 4. 34. Zienkiewics, o. c., The Finite Element Method, McGrawHill, London, 1977 : Chaps 1,3,7, and 8. 35. Zienkiewics, o. c. and Cheung, Y. K., "Finite Element in the Solution of Field Problems", The Engineer, 1965 : 507:510. PAGE 114 APPENDIX A THE DERIVATION OF ELEMENT EQUATION BY GAUSS-GREEN THEOREM The Gauss-Green divergence theorem in a xy-plane states (Kreyszig, 1983): Let R be a closed bounded region in the xy-plane whose boundary r consists of finitely many smooth curves. Let U(x,y) be a vector function which is continuous and has continuous first partial derivatives in some domain containing R. Then, we can write I I R di v u dxdy = I r u n dr (Eq. A.l) where div = o +-o-6X oy and n = the outer unit normal vector. If we write u and n in terms of unit components, namely ( Eq. A. 2) n = cose i + sine j (Eq. A 3) Substituting Eq's .A.2 and A.3 into Eq. A.l gives + (Eq. A. 4) The integral equation we have to deal with is ( ) I T 02'" (R e } = -[N] (Ox _,.._2 + Dy Gp + Q)dA ox oy PAGE 115 Using the product rule for differentiation in the above equation gives 106 { R (e) } = -J [ 0 _!:____ ([ N] T ocp ). + 0 ( [ N] T ocp ) ] dA A x ox ox Y oy oy o[NJT ocp o[NJT ocp + JA( 0 -+ 0 +G[N]Tcp + Q[N]T)dA x ox ox Y oy oy ( Eq. A. 5) Applying Eq. A.4 to the first integral of the right-hand side of Eq. A.5, it becomes = ocp ocp Jr( o [N]T -cose + oy [NJT -sine )dr x ox oy (Eq. A. 6) Eq. A.6 can be transformed to its final form using the following relationship. (Eq. A. 7) Substitution and rearranging yield J T ocp ocp {R(e)} =-r.( ox [NJ -cose + oy [N]T -.-sine )dr ox oy JAc ox o[N]T o[NJ o[N]T o[NJ )dA ] {t(e)) + [ + 0 ox ox y oy oy + ( JAG[N]T[N] dA ) { t (e)} -JAQ[N]T dA (Eq. A. 8) PAGE 116 107 As to a two-dimensional time-dependent problem, the residual integral is {R(e)} = -J [N]T ( cp(x,y) + Qgf )dA (Eq. 4.19) D + o2tn where cp(x,y) = D Gcp x ox2 Y oy2 Comparing Eq. 4.19 with Eq. 4.13, the only additional term is I ( [N]T )dA ot ( Eq. A. 9) Assume that ocpjot varies quadratically within any element. ocpjot can be expressed as = N ;. (e) ot i 'i!i (Eq. A.10). where t = nodal values of ocpjot. Substituting Eq. A.10 into Eq. A.9 gives (Eq. A.11) Adding Eq. A.ll to Eq. A.S, the final form of the residual integral of a time-dependent field problem is {R(e)} ocp ocp = -J ( D [N]T -case + Dy [N]T -sine )dr r x ox oy + r JAc o[N]T o[NJ o[N]T o[Nl {t(e)} D + D )dA ] X ox ox y oy oy PAGE 117 + ( JAG[N]T[N] dA ) {t(e)} JAQ[N]T dA + ( J [N]T[N] dA ) {;(e)} 108 A.12) PAGE 118 APPENDIX B THE INTEGRATION OF SHAPE FUNCTIONS Since the shape functions are expressed in terms of natural coordinates, to perform integration, the change of variables is required. In case of a double integral J I f(x,y)dxdy R the new variables of integration a,B can be introduced by setting x = x(a,B), and y = y(a,B) where the functions x(a,B) and y(a,B) are continuous and have continuous first derivatives in some region R* in the aB-plane so that each point (a,B) in R* corresponds to a point [x(a,B), y(a,B)] in Rand conversely. Kreyszig (1983) indicated that I I f(x,y)dxdy =I I f[x(a,B),y(a,B)] R R* o(x,y) I dadB o(a,B) (Eq. B.l) where o(x,y)jo(a,B) is the Jacobian matrix, which can further be expressed as [Jl = o(x,y) o(a,B) = ox OQ ox 613 (Eq. B. 2) To prove Eq. B.l, the two coordinate systems as PAGE 119 110 shown in Fig 27 were considered (Chung; 1978). The directions of the cartesian coordinates x-y and the arbitrary nonorthogonal (possibly curvilinear ) isoparametric coordinates a-B are expressed by qnit vectors i1 i2 and the tangent vectors t1 t2 respectively. The relationship between i and t are t = 2 ox oy --oa ox oB The differential area (shaded) is ox _l!y oa oa o dadB or dxdy i3 = I det[J] I dadB i3 Thus, Eq. B.l can be obtained. ox ..21.. oB oB 0 (Eq. B.J) For an isoparametric quadratic quadrilateral element, the transformation equations are given by y = N. (a,B) Y. l. l. where N1 are shape functions, x.,Y. are the global l. l. coordinates of the nodes. The Jacobian matrix can be written as the matrix product as PAGE 120 i \ 1 1 1 i__ y X Figure 27 Coordinate transformation PAGE 121 112 x1 y2 oN1 oN2 oN3 oN4 oN5 oN6 oN7 oN8 x2 y2 oa oa oa oa oa oa oa oa XJ YJ [J] = oN1 oN2 oN3 oN4 oN5 oN6 oN7 oN8 x4 y4 oB oB oB os oB as a a as xs y5 x6 y6 x7 y7 xa Ya (Eq. B.J) in which oN1 joa = \(1-B)(2a+B), oN2 ;oa = and so on. In Eq. B.J, the terms of oNifoa and oNi/68 can be evaluated. The global coordinates of nodes, i.e., X. and l. Yi, are already known. Therefore, the Jacobian matrix for a element can be obtained, which is a 2x2 matrix. Substituting the determinant of the Jacobian matrix into the integrate equations, i.e., Eq's 4.35 and 4.36, the integral without the partial derivative terms can be performed. For those integral containing the partial derivative terms, we need to change the local derivative in terms of a and B to global derivatives in terms of x and y. Applying the chain rule for differentiation toa shape function Ni(a,B)t it gives oNi(a,B) oNi(a,B) 6x oNi(a,B) oy = --+ oa ox oa oy oa oNi(a,B) oNi(a,B) ox oNi(a,B) oy oB = ox ""68 + oy 08 or in a matrix form '. ; PAGE 122 113 i oN. ox oy oN. l. (a,B) l. (a, B) oa 6a ox = oN. ox oy oN. l. (a, B) l. (a,B) oB 68 68 oy (Eq. B.4) The coefficient matrix is the Japobian matrix [J]. To find the global derivatives we invert [J] and write oN. oN. l. (a,B) l. (a,B) 6x -1 oa = [J(a,B)] (Eq. B.S) oN. oN. l. (a,B) l. (a,B) oy oB However, an explicit form of [J]-1 pannot be obtained (Zienkiewicz 1977). The numerical methods have to be used. In this study, the Gauss method is recommended because it has been proved to be in finite element method. Gauss-Legendre Quadrature To approximate the integral 1 R = J f(x) dx -1 (Eq. B. 6) we can select f(x) at the midpoint, i.e., f(O) and multipty by the length of the interval, i.e., 2. Thus we find R 2f(O). This result is an exact solution if the function y is a straight line. Generalization of Eq. B.6 yields 1 R = J f ( x) dx :E ( w i f (xi ) ) -1 ( Eq. B. 7) where = number of points selected. Thus, to approximate PAGE 123 114 R, we evaluate f(xi) at each location xi, multiply f(xi) by an appropriate weighting value Wi' and sum the resultants up. The Gauss-Legendre quadrature locates the sampling points to achieve the greatest accuracy. Sampling points are located symmetrically with respect to the center of the interval. Symmetrically paired points have the same weighting value Wi. In general, a polynomial of degree 2n-1 is integrated exactly by n sampling points. further R = Applied to two-dimensional be expressed as 1 1 1 J_1J-1f(a,B) dadB I_1 problems, Eq. B.7 may n [,:E(Wj f (aj,B))] dB ]=1 m n :E :E [W.W.f(a.,B.)] i=1 j=1 1 J J 1 (Eq. B.8) The locations of sampling points and weighting coefficients can be found in Reference 4, a, 9, and 29. PAGE 124 Symbol A a a B b c g hL [J] k K L N. l. p Re APPENDIX C DEFINITION OF SYMBOLS Definition Cross-sectional area The width of a rectangular duct The cross sectional area of orifice Half of the width between two plates The height of a rectangular duct Discharge coefficient Darcy friction factor Pipe diameter Hydraulic diameter Half of the opening width at the orifice Darcy friction factor for pipe flow Acceleration of gravity Head loss Jacobian matrix Friction coefficient ( k = f Re ) Loss coefficient Length of conduit Shape function Pressure Pressure difference Flow Reynolds number PAGE 125 r t u u v v X y w. 1 Q B r r 1 e a Reynolds number for the flow in a rectangular duct Hydraulic radius Opening ratio Time 116 Inflow velocity for a Poiseuille flow Mean flow velocity through the orifice Velocity in x-direction Cross sectional average flow velocity Velocity in y-direction Cartesian coordinate, usually parallel to main flow direction Cartesian coordinate, usually perpendicular to main flow Weighting function Natural coordinate Natural coordinate Element boundary Shear stress Absolute viscosity Density Nodal value The angle to the outward normal Rotatoional velocity Stream function Vorticity PAGE 126 APPENDIX D THE FORTRAN SOURCE CODE { PAGE 127 c program TDICF common/elmatx/esmC8,8) ,efCB> ,phiCB) ,nsCS) ,wsmC8,8) ,p(200) common/matl/ge,qe,ivor,vk,wvC600) common/pdxy/vnCB),pnxC8l,pnyC8>,xxC8),yy(8),xd,yd,det common/crd/xc(600),ycC600),velx(600),velyC600) common/tle/title(20) common/av/aC40000),jqf,jqsm,np,nbw dimension nel(200,8),ick(600),w(8),q(200),sv(600) dimension iwc(l50,2),1pc(20),u(8),v(8) dimension ibc(l50,2),bc(l50),cc(40000) data in/60/,io/61/ openCin,f1le='FLOW.INP',status='old',form='formatted', $ access='sequential') open( io,file= 'FLOW. OUT' ,status= 'unknown' ,form=' formatted', $ access='sequential') c INPUT OF THE TITLE CARD AND PARAMBI'ERS c c writeC",l) 1 format(/lx,'Plaese enter Reynolds number :') readC",">re vk=l./re wr1teC",2) 2 format(/lx,'T1me interval for every computational step?') readC","Jdeltat readC1n,5)t1tle 5 format(20a4) read(in,")np,ne,nbc,nwc,npc read(in,")Cxc(i),i=l,np) read(in,")(yc(i),1=l,np) writeCio,8)title,np,ne,nbc,nwc,npc,re 8 format(////10x,20a4//lOx,5hNP =,16/lOx,ShNE =,16/ $ lOx,ShNBC =,i6/10x,ShNWC =,i6 $ /lOx,ShNPC =,16//lOx,ShRE =,16) do 9 kk::l,ne read(in,")n,(ne1Cn,i),1=1,8) 9 continue read(in,") (1bc(i,l),ibcC1,2),bc(1),1=l,nbc) read(in,"> CiwcCi,l),iwcCi,2),1=l,nwc) read(in,") (ipc(i),i=l,npc) c CALCULATION OF BANDWIDTH c nbw::Q do 20 kk,;.l,ne do 22 1=1,8 22 ns(i)=nel(kk,i) lk=7 do 21 1=l,lk ij=i+l do nj=ij,8 nb=iabs(ns(i)-ns(j)) ifCnb.le.nbw) qoto 21 inbw=kk nhw=nh 21 continue 20 continue nhw=nhw+l writeCio,24) nhw,inbw 24 format(//10x,l2hBANDWIDTH IS,I4,11H IN ELEMENT,I4) 118 PAGE 128 c c INITIALIZATION OF THE COLUMN VECTOR A() c -jgf=np jgsm=jgf+np j1=jend-jgf if(jend.gt.40000)then 25 format(l0x,30hDIMENSION OF A VECTOR EXCEEDED/lOX, $ 20HEXECUTION TERMINATED) goto 100 end if c -, c GENERATION OF THE SYSTEM OF EQUATIONS c c do 26 i=l,np 26 wvCi)=O. istep=O ivor=O istop=O 30 continue do 27 i=l,jend a(i)=O. cc(i)=O. 27 continue do 32 kk=l,ne do 31 1=1,8 ns(i)=nel(kk,i) 31 j=ns(i) ge=O. qe=O. ifCivor.eq.l)then qe=qCkk) ge=l. end if iftivor.eq.2)then ge=O. qe=pCkk) end if CALL ELSTMX C KK) ----c DIRECT STIFFNESS PROCEDURE c c:: do 33 i=1,B ii=ns(i) aCjgf+ii)=aCjgf+ii)+ef(i) do 34 j=l,B jj=ns(j)+l-11 if(jj.1e.O) goto 34 cc(ji)=cc(ji)+wsm(i,j) 34 continue 33 continue 32 continue c:: MODIFICATION AND SOLUTION OF THE SYSTEM OF EQUATIONS c:: if(ivor.ne.l)then if(ivor.eq.2)goto 75 119 PAGE 129 do 40 i=l,nbc if(ibc(i,2).eq.O)then ib=ibc(i,l) bv=bc(i) end if CALL MODIFYCib,bv) 40 continue CALL DCMPBD CALL SLVBD do 56 i=l,np 56 sv(i)=a(i) 98 continue do 76 i=l,np velx(i)=O. 76 vely(i)=O. do 70 kk=l,ne do 71 1=1 ,8 ns(i)=nel(kk,i) j=ns(i) ,. 11 phi( i> =a( j) ELGRAD ( KK) 70 continue 81 do 80 kk=l,ne do 81 1=1,8 ns(i)=nel(kk,i) j=ns(i) xx(i)=xc(j) yy(i)=yc(j) u(i)=velx(j) v(i)=vely(j) continue gdxu=O. gdxv=O. gdyu=O. gdyv=O. call pderv(O.,O.) do 83 j=l,8 83 gdyv=gdyv+pny(j)*v(j) continue 80 continue ivor=2 goto 30 75 continue do 54 i=l,npc ib=ipcCi) bv=lOO. CALL MODIFY{ib,bv) 54 continue CALL DCMPBD CALL SLVBD .. '.'. 44 format(//lOx,'step =',f8.4/) write(io,45) 45 formatCllx,'NO.',lx,'Streamline',4x,'Vorticity' ,5x'Pressure', $ 9x, 'Vel-x' ,6x, 'Vel-y' I) 120 PAGE 130 write(io,43) (i,sv(i),wv(i),a(i),velx(i),velyCil,i=l,np) 43 format(l0x,13,2x,fl0.5,2x,fl0.5,2x,fl2.2,4x,f9.5,2x,f9.5) end if if(istep.qt.3000)qoto 100 if(istop.eq.l)goto 100 do 47 i=l,nwc if(iwc(i,2).ne.Olthen j=iwc(1,1) k=iwc(1,2) & end if 47 continue do 51 kk=l,ne do 52 1=1,8 ns(i)=nelCkk,i) j=ns(i) xxC1)=xc(j) yy(i)=yc:(j) w(i)=wv(j) 52 phi(1)=sv(j) qdxphi=O. qdyphi=O. qdxw=O. qdyw=O. PDERV(O.,O.) do 53 j=l,8 53 continue 51 continue ivor=l istep=istep+l goto 30 end if do 60 i=l,np sum=o. ilimi t=nbw+i-1 if(i11mit.gt.np)11imit=np do 61 j=l,111mit if(i.le.j)then k=j-1 else k=i-j+1 ij=iij+jgsm ifCiij.qt.nnc)goto 61 end if 61 continue aCjgf+i)=aCjgf+i)+sum/deltat 60 continue do 6.2 1=1 ,nne 121 PAGE 131 c j=i+jgsm a(j)=a(j)+cc(j)/de1tat 62 continue do. 63 i=1,nwc if(iwc(i,2).eq.O)then ib=iwc(i,1) hv=o. else ib=iwcC i, 1) bv=wv(ib) endif CALL MODIFY(ib,bv) 63 continue CALL DCMPBD CALL SLVBD do 90 i=1,np ick(i)=l def=aCi)-wv(i) deft=def/de1tat if(deft.lt.0.04) ick(i)=O 90 continue do 99 i=l,np if(ick(i).eq.1)then do 64 i=1,np 64 wv(i)=a(i) ivor=O goto 30 end if 99 continue istop=l goto 98 100 continue stop end SUBROUTINE ELSTMXCKK) c EVALUATES THE ELEMENT STIFFNESS MATRIX AND ELEMENT FORCE VECTOR c common/e1matx/esm(8,8),ef(8),phi(8),ns(8),wsm(8,8),p(200) common/matl/ge,qe,ivor,vk,wv(600) common/pdxy/vn(B),pnx(B),pny(B),xx(B),yyCB),xd,yd,det common/ipts/vxC9r,vyC9) ,wc(9) common/crd/xc(600),yc(600),velx(600),ve1y(600) dimension w(B),ff(8,8) do 1 1=1,8 ef(i)=O.O do 1 j=l,8 ff(i,j)=O. wsm(i,j)=O.O 1 esm(i,j)=O.O do 10 1=1,8 j=ns(i) w( i) =wv( j) xx(i)=xc(j) 10 yy(i)=yc(j) 122 PAGE 132 CALL INGPTS do 3 ii=l,9 call pderv(vx(ii),vy(ii)) do 2 i=l,B do 2 j=l,B if ( 1 vor. eq. 0) ff ( i, j) =ff ( i, j) +vn( 1) j) ( i1) 2 continue 3 continue if(ivor.eq.l)then. do 4 i=l,B do 4 j=l,B 4 continue end if if(ivor.eq.O)then do 5 i=l,B do 5 j=l,B 5 continue end if return end SUBROUTINE ELGRADCKK) comrnon/elrnatx/esrn(8,B),ef(B),phi(8),ns(B),wsm(B,8),p(200) common/pdxy/vn(8),pnx(8),pny(8),xx(B),yy(8),xd,yd,det cornmon/ipts/vx(9),vy(9),wc(9) cornmon/crd/xc(600),yc(600),velx(600),vely(600) dimension xq(9),yq(9),gdx(9),gdy(9) xq/-l.,O.,l.,l.,l.,0.,-1.,-1.,0./ data yq/-l.,-l.,-l.,O.,l.,l.,l.,O.,O./ data io/611 do 1 i=l,B j=ns(i) xx(i)=xc(j) yy(i)=yc(j) 1 continue do 10 1=1,9 gdx(i)=O. gdy(i)=O. CALL PDERV PAGE 133 c velx(j)=Cvelx(j)+qdx(1))/2. endif ifCvely(j).eq.O.)then vely(j)=qdy(i) else vely(j)=(vely(j)+qdy(1))/2. end if 15 continue return end SUBROUTINE INGPTS common/ipts/vx(9),vy(9),wc(9) dimension a(3),b(3) data data b/5.,8.,5./ c GENERATION OF THE NINE INTERGATION POINTS c c n=O do 1 1=1,3 do 1 j=l,3 n=n+l vxCn)=aC1) vyCn)=aCj) 1 return end SUBROUTINE PDERV(Xl,X2) C EVALUATES THE JOCOBIAN TRANSFORMATION MATRIX C AND USES THE INVERSE JOCOBIAN MATRIX TO OBTAINTHE DERIVATIVES C OF THE SHAPE FUNCTIONS WITH TO X AND Y. c common/pdxy/vnn(8),pnxC8),pny(8),x(8),y(8),xd,yd,det common/derv/vn(8),pnsCB),pne(8) real jocb(2,2) do 1 1=1,2 do 1 j=l,2 1 jocb(i,j)=O.O xd=O.O "yd=O.O CALL QDSHFNCXl,X2) do 2 1=1,8 2 .. jocbC 2, 2 > =joc:bc 2, 2 > +pnec 1 > b=jocb(l ,1) jocb(l,l)=jocb(2,2)/a jocb(l,2)=-jocb(l,2)/a jocb(2,1)=-jocb(2,1)/a jocbC 2, 2) =b/a det=abs(a) 124 PAGE 134 c do 3 1=1,8 vnn(i)=vn(i) pnx(i)=jocb(l,l)*pns(i)+jocb(1,2)*pne(i) pny(i)=jocbC2,1)*pns(i)+jocbC2,2)*pne(i) return end SUBROUTINE QDSHFNCSI,ETA) 125 c CALCULATES THE VALUE OF THE SHAPE FUNCTIONS AND THEIR DERIVATIVES c common/derv/vn(B),pns(B),pne(B) dimension sq(B),eqCB) data sq/-1.,0.,1.,1.,1.,0.,-1.,-1./ data eq/-1.,-1.,-1.,0.,1.,1.,1.,0./ do 5 1=1,7,2 vn(i)=.25*C1.+si*sq(i))*C1.+eta*eq(i))*(si*sq(i)+etalleq(i)-1.) pns(i)=.25*C1.+eta,i:eqCi))ll(2.,r,:si+eta*eq(i)*sqCi)) 5 pne( i) =. 25* C l.+si*sq( i)) *C 2. *eta+si*sq( i) lleq( i) )_ do 6 1=2,6,4 vn(i)=.5*(1.-s1**2>*C1.+eta*eq(i)) pns(i)=-si*C1.+eta*eq(i)) 6 pneCi)=.5*eqCi)*C1.-si**2) do 7 1=4,8,4 vnC i) =. 5* ( 1. +si*sq( i)) *< 1. -eta**2) pns(i)=.5*sq(i)*(1.-eta**2) 7 pne(i)=-eta*C1.+si*sq(i)) return end .. SUBROUTINE MODIFY(ib,bv) common/av/a(40000),jqf,jqsm,np,nbw data inl60/,io/61/ k=ib-1 do. 211 j=2,nbw m=ib+j-1 if(m.qt.np) qoto 210 -aCjgf+m)=aCjqf+m)-a(ij)*bv aCij)=O.O 210 if(k.le.O) qoto 211 kj=jgsm+(j-l)*np+k-Cj-l)*Cj-2)/2 a(jgf+k)=aCjqf+k)-a(kj)*bv a(kj)=O.O k=k-1 211 continue a(jqf+ib)=a(jqsm+ib)*bv 221 continue return end -SUBROUTINE DCMPBD common/av/a(40000),jqf,jgsm,np,nbw io=61 npl=np-1 do 226 i=1,npl mj=i+nbw-1 if(mj.qt.np) mj=np PAGE 135 nj=i+1 mk=nbw if((np-i+1).1t.nbw) mk=np-i+1 nd=O do 225 j=nj,mj mk=mk-1 nd=nd+1 n1=nd+1 do 225 k=1,mk nk=nd+k ii=jgsm+i 225 a(jk)=a(jk)-a(in1)"a(ink)/a(ii) 226 continue return end SUBROUTINE SLVBD common/av/a(40000),jgf,jqsm,np,nbw npl=np-1 do 250 i=1,np1 mj=i+nbw-1 if(mj.gt.np) mj=np nj=i+1 1=1 do 250 j=nj,mj 1=1+1 250 continue a(npl=a(jqf+np)/a(jgsm+np) do 252 k=1,np1 i=np-k mj=nbw if((i+nbw-1).gt.np) mj=np-1+1 sum=O.O do 251 j=2,mj n=i+j-1 251 252 a(i)=(a(jgf+i)-sum)/a(jgsm+i) return end 126 |