UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Finite element analysis of periodic viscous flow in a constricted pipe Singh, Rajesh Kumar 1989

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1989_A7 S55.pdf [ 6.34MB ]
Metadata
JSON: 831-1.0062609.json
JSON-LD: 831-1.0062609-ld.json
RDF/XML (Pretty): 831-1.0062609-rdf.xml
RDF/JSON: 831-1.0062609-rdf.json
Turtle: 831-1.0062609-turtle.txt
N-Triples: 831-1.0062609-rdf-ntriples.txt
Original Record: 831-1.0062609-source.json
Full Text
831-1.0062609-fulltext.txt
Citation
831-1.0062609.ris

Full Text

FINITE E L E M E N T ANALYSIS OF PERIODIC VISCOUS F L O W IN A CONSTRICTED PIPE by Rajesh K u m a r Singh B . Tech., Indian Institute of Technology, Kanpur , 1987 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF M A S T E R O F A P P L I E D S C I E N C E in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF CIVIL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA June, 1989 © Rajesh K u m a r Singh, 1989 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of Br i t i sh Columbia , I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my depart-ment or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of C i v i l Engineering The University of Br i t i sh Columbia Vancouver, B . C . , C A N A D A V 6 T 1W5 June, 1989 Abstract The finite element method is extended to analyze axisymmetric periodic viscous flow i n a constricted pipe. This method is used to solve the pulsatile fluid flow through an artery which may be stenosed due to impingement of extravascular masses or due to intravascular atherosclerotic plaques developed at the wal l of the artery. The artery geometry is approximated as an axisymmetric channel. B lood is assumed to be a Newtonian fluid, i.e., it follows a linear relationship between the rate of shear strain and the shear stress. The numerical model uses the finite element method based on velocity-pressure pr imit ive variable representation of the Navier-Stokes equations. Eight-noded iso-parametric elements wi th quadratic interpolation for velocities and bilinear for pres-sure are used. A truncated Fourier series is used to approximate the unsteadiness of flow and a modified method of averaging is used to obtain a periodic solution. Two non-dimensional parameters are used to characterize the flow: the frequency Reynolds number i ^ , and the steady Reynolds number Re. The pulsatile flow is modeled as a linear combination of steady, cosinusoidal and sinusoidal components. The mag-nitude of each component is determined by minimizing the error in approximation through Galerkin's procedure. Results are obtained for various values of and RE for laminar flow. Results are also presented for l imit ing cases and, wherever possi-ble, numerical results thus obtained are compared wi th analytical and experimental results published in the literature. i i Contents Abstract ii Tables vi Figures vii Notation xii Acknowledgement xvi 1 Introduction and Literature Review 1 1.1 General Remarks 1 1.2 Stenosis and Pulsatile Flow 3 1.3 Analysis in Time Domain 6 1.4 Scope of Present Investigation 7 2 Derivation of Governing Equations 8 2.1 Conservation Equations and Boundary Conditions 8 2.2 Non-dimensional Form of Governing Equations 12 2.3 Stream Function and Poisson's Equation 14 3 Finite Element Formulation 16 3.1 Discretized Form of Governing Equations 16 3.2 Boundary Conditions 21 3.3 Steady State Periodic Solution 22 in 3.4 Characterization of Arter ia l F l u i d F low 25 3.5 Streamlines and Streaklines 30 4 Numerical Results 33 4.1 Steady Flow Through a Pipe 34 4.1.1 Theoretical Results 34 4.1.2 F in i te Element Results 35 4.2 Oscil lat ing F low Through a Pipe 36 4.2.1 Theoretical Results 36 4.2.2 F in i te Element Results 41 4.3 Developing Flow in the Inlet Length of a Pipe 49 4.3.1 Theoretical Results 49 4.3.2 F in i te Element Results 51 4.4 F low Through an Axisymmetr ic Orifice 52 4.5 F low Through a Stenosed Pipe 65 4.5.1 Geometry of the Stenosis 65 4.5.2 Steady Flow Through Stenosis 67 4.5.3 Pulsati le Flow Through Stenosis 84 5 Conclusions 113 5.1 Concluding Remarks 113 5.2 Suggestions for Further Developments 115 Bibliography 116 Appendixes 121 A Newton-Raphson Procedure 121 B Details of Tangent Stiffness Matrix and Nonlinear Load Vector 123 iv C Determination of the Stream Function Values D Determination of Stress Components Tables 4.1 Comparison of Finite Element Solution with Exact Solution 37 4.2 Comparison of exact and finite element solution for oscillating flow for i% = 1.0 42 4.3 Comparison of exact and finite element solution for oscillating flow for 'Re = 10.0 42 4.4 Comparison of exact and finite element solution for oscillating flow for Re = 20.0 43 4.5 Instantaneous rate of flow through a pipe under cosinusoidal pressure difference 45 4.6 Parameters for the finite element grids for stenotic flow 70 4.7 Maximum stream function values for steady flow through the stenosis 84 vi Figures 2.1 Fluid Domain 9 3.1 Isoparametric Element (a) Element in (£, n) space, (b) Element in (x, y) space 18 3.2 Typical arterial flow waveforms (Young, 1979) 27 3.3 Function H(t) 28 4.1 Parallel Flow with Parabolic velocity distribution 34 4.2 One Element Representation of the Fully Developed Flow 36 4.3 Flow in a pipe under oscillating pressure 40 4.4 Finite element grid for oscillating flow through a pipe 41 4.5 Variation of cosinusoidal and sinusoidal components for unsteady par-allel flow for Ru, = 1, 10 and 20 44 4.6 Streaklines for R^ = 1 for t = 0, —, TT and — 46 4.7 Streaklines for R^ = 10 for t = 0, —, it and — 47 7T 37T 4.8 Streaklines for Rw = 20 for t = 0, - , TT and — 48 4.9 Flow in the inlet section of a Pipe 50 4.10 Finite element grid for developing flow 53 4.11 Streamline contours for developing flow for Re =0, 1 and 10 54 4.12 Streamline contours for developing flow for Re =50, 100 and 150 . . . 55 4.13 Streamline contours for developing flow for Re =200, 250 and 300 . . 56 4.14 Streamline contours for developing flow for Re =350 and 400 57 4.15 ^ versus Re 58 R vn 4.16 Velocity distribution along the centerline for different Re 59 4.17 Flow through an axisymmetric orifice 60 4.18 Finite element grid for flow through an axisymmetric orifice 61 4.19 Streamline contours for orifice flow for Re =0, 1 and 5 . . . 62 4.20 Streamline contours for orifice flow for Re =10, 15 and 20 63 4.21 Comparison of streamline contours for orifice flow for Re = 10 . . . . 64 4.22 Geometry of Stenosis in an axisymmetric pipe 66 4.23 Grid 1 used for analysis of stenotic flow 68 4.24 Grid 2 used for analysis of stenotic flow 69 4.25 Streamline contours for steady flow through stenosis for grid 1 for Re = 0, 1 and 10 72 4.26 Streamline contours for steady flow through stenosis for grid 1 for Re = 50, 100 and 150 73 4.27 Streamline contours for steady flow through stenosis for grid 1 for Re = 200, 250 and 300 74 4.28 Streamline contours for steady flow through stenosis for grid 2 for Re = 0, 10 and 50 75 4.29 Streamline contours for steady flow through stenosis for grid 2 for Re = 100, 150 and 200 76 4.30 Streamline contours for steady flow through stenosis for grid 2 for Re = 300, 400 and 500 77 4.31 Streamline contours for steady flow through stenosis for grid 2 for Re = 600, 700 and 800 . 78 4.32 Streamline contours for steady flow through stenosis for grid 2 for Re = 900, 1000 and 1250 79 4.33 Streamline contours for steady flow through stenosis for grid 2 for Re = 1500, 1750 and 2000 ' . . . '80 viii 4.34 Shear stress along the wall of the pipe for grid 1 for Re = 0 (linear), 50, 100, 200 and 300 81 4.35 Shear stress along the wall of the pipe for grid 2 for Re = 0 (linear), 50, 100, 200 and 300 82 4.36 Shear stress along the wall of the pipe for grid 2 for Re = 400, 500, 1000, 1500 and 2000 83 4.37 Streaklines for pulsatile flow through stenosis for Rw = 5, Re = 0 7T 37T (linear) for t = 0, —, it and — . 86 Li L» 4.38 Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 1 for . 7T 37T t = 0, —, 7r and — 87 Li Li 4.39 Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 5 for „ 7T , 37T * = 0, - , 7T and — 88 Lt Lt 4.40 Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 10 for 7T 3TT < = 0, | , 7T and y 89 4.41 Streaklines for pulsatile flow through stenosis for R^, = 5, Re = 20 for „ 7T , 3ff < = 0, - , 7T and y 90 4.42 Streaklines for pulsatile flow through stenosis for R^, = 5, Re = 50 for „ 7T , 37T t = 0, - , 7T and — 9 1 4.43 Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 100 for 7T 37T —, 7r and — 2' 2 4.44 Streaklines for pulsatile flow through stenosis for Rw = 5, Re = 150 for ^ 7T 37T f = 0, - , 7T and — . 9 3 4.45 Streaklines for pulsatile flow through stenosis for R^ = 10, Re = 0 (linear) for t = 0, —, 7r and — 94 Li Li 4.46 Streaklines for pulsatile flow through stenosis for R^ = 10, Re = 1 for . 7T 37T t = 0, —, 7r and — 95 f = 0, - ,  92 IX 4.47 Streaklines for pulsatile flow through stenosis for Rw = 10, Re = 5 for t = 0, —, 7r and — 96 Li Lt 4.48 Streaklines for pulsatile flow through stenosis for R^ = 10, Re = 10 for t = 0, —, 7r and — 97 Li Lt 4.49 Streaklines for pulsatile flow through stenosis for R^ = 10, Re = 20 for t = 0, - , 7T and — 98 Lt 4.50 Streaklines for pulsatile flow through stenosis for Rw — 10, Re = 50 for „ i" , 3 7 T * = 0, - , 7T and —- 99 Li Lt 4.51 Streaklines for pulsatile flow through stenosis for R^ = 10, Re = 100 7T 37T for t = 0, - , 7r and — . . . . 100 Li Lt 4.52 Streaklines for pulsatile flow through stenosis for Rw = 10, Re = 150 for t = 0, —, 7r and — . . . 101 Li Li 4.53 Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 0 7T 3TT (linear) for £ = 0, —, 7r and — 102 Lt Lt * 4.54 Streaklines for pulsatile flow through stenosis for Rw = 20, Re = 1 for „ ir , 3ir t = 0, - , 7T and — . 103 4.55 Streaklines for pulsatile flow through stenosis for Rw = 20, Re = 5 for * = 0, - , x and — 104 4.56 Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 10 for 7T , 3TT —, 7r and — 2' 2 4.57 Streaklines for pulsatile flow through stenosis for Rw = 20, Re = 20 for t = 0, - , ir and — 106 4.58 Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 50 for 7T 3?T t = 0, - , ir and — 107 4.59 Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 100 „ , 3 7 T for t = 0, - , 7T and —- 108 t = 0, —, 7r and — 105 x 4.60 Streaklines for pulsatile flow through stenosis for i ? w = 20, Re = 150 for t = 0, - , 7T and — 109 Lt Lt 4.61 Instantaneous shear stress along the wall of the pipe for grid 1 for Rw = 5 and Re = 0 (linear), 10, 50, 100 and 150 . . . . . . . 110 4.62 Instantaneous shear stress along the wall of the pipe for grid 1 for R^ = 10 and Re = 0 (linear), 10, 50, 100 and 150 I l l 4.63 Instantaneous shear stress along the wall of the pipe for grid 1 for Rw = 20 and Re = 0 (linear), 10, 50, 100 and 150 112 D . l Use of transformation (x = x(x', y'); y = y(x', y')) to determine stresses along a curved wall 130 xi Notation Lower case Greek symbols a "periodicity factor", real part of a complex number ft imaginary part of a complex number 5 radius at the minimum section e non-dimensional parameter for developing flow £ vorticity of the flow 77 non-dimensional radial coordinate 9 angle of rotation K grid efficiency factor fi viscosity of the fluid v kinematic viscosity £ non-dimensional axial coordinate p density of the fluid c^ azimuthal stress in three dimensions cTr radial stress in three dimensions crx axial stress in two dimensions o~y radial stress in two dimensions o~z axial stress in three dimensions T0 Z shear stress on azimuthal face in radial direction TR<P shear stress on radial face in azimuthal direction TTZ shear stress on radial face in axial direction rxy shear stress on axial face in radial direction <t> azimuthal coordinate in three dimensions xii yj stream function u frequency u>~1 characteristic time Upper case Greek symbols T fluid boundary Ts traction boundary T u kinematic boundary T e finite element boundary finite element traction boundary T„ finite element kinematic boundary A increment (prefix) Q fluid domain fle finite element domain Lower case Roman symbols a pre-stenosis length, pre-orifice length a(y) steady part of pulsatile flow b post-stenosis length, post-orifice length b(y) magnitude of cosinusoidal part of pulsatile flow d nodal vector of unknowns f consistent load vector / length of stenosis, length of orifice IE inlet length n a axial cosine of unit outward normal to the boundary n 2 radial cosine of unit outward normal to the boundary p pressure &Padd additional pressure drop xiii r radial coordinate in three dimensions r(x) instantaneous radius t time u axial component of velocity UQ characteristic velocity v radial component of velocity vr radial component of velocity in three dimensions azimuthal component of velocity in three dimensions vz axial component of velocity in three dimensions x axial coordinate in two dimensions y radial coordinate in two dimensions z axial coordinate in three dimensions Upper case Roman symbols A steady streaming solution Ao(y) steady part of pulsatile approximation BQ(y) magnitude of cosinusoidal part of pulsatile approximation B(<) magnitude of cosinusoidal component Co(y) magnitude of sinusoidal part of pulsatile approximation C(t) magnitude of sinusoidal component C ° continuity of zeroth order derivatives C 1 continuity of first order derivatives Fr body force in radial direction body force in azimuthal direction Fz body force in axial direction J Jacobian matrix Jo Bessel's function of the first kind of order zero xiv K stiffness matrix M mass matrix Mj j—th. bilinear isoparametric shape function Ni i—th quadratic isoparametric shape function Q(t) instantaneous rate of flow R radius of the pipe Rw frequency Reynolds number R0 characteristic length Re Reynolds number 72,/, error in approximating Poisson's equation 7£i error in approximating x—momentum equation 7^ 2 error in approximating y—momentum equation 72.3 error in approximating continuity equation Tin error in approximating x—component of traction boundary 72-r2 error in approximating y—component of traction boundary T tangent stiffness matrix U axial component of specified kinematic boundary Uo (maximum) entrance velocity V radial component of specified kinematic boundary X axial component of specified traction boundary Y radial component of specified traction boundary YQ Bessel's function of the second kind of order zero ZQ development length D ' . , , . . d d d — total material derivative = — + u——V v— Dt at ox oy ' differentiation non-dimensional value differentiation with respect to time xv Acknowledgement My special thanks to my supervisor Dr. M.D. Olson, for his guidance and encour-agement throughout the course of my research work and in the preparation of this thesis. I would also like to thank Mr. Tom Nicol of the UBC Computing Center for many valuable suggestions in using the SPARSPAK package. Financial support in the form of a Research Assistantship from the Natural Sciences and Engineering Research Council of Canada and a Graduate Fellowship from the University of British Columbia is gratefully acknowledged. xvi To My Parents Indu and Madan Singh xvii C H A P T E R 1 Introduction and Literature Review 1.1 General Remarks F l u i d mechanics poses one of the most difficult problems for applied mathematics: the solution of the Navier-Stokes equations and simplifications thereof for flow conditions dictated by practical applications of the laws of fluid mechanics and technological de-velopments, rather than by academic questions. It was only after the advent of digital electronic computers, that the solutions for such problems could be developed with success. In the beginning of fluid mechanics research, l imit ing situations—offering permissible simplifications of the governing equations to solvable ones—were primar-i ly the goal of research. Thus came the Prandtl ' s lift line theory, the Prandtl-Meyer flow, or the Blasius series solution. Even though these problems were relatively quite simple, a l l of them already involved numerical work and a substantial amount of t ime i n which those solutions were derived. Prandtl 's lifting line theory yielded an integro-differential equation for the circulation, the Prandtl-Meyer flow a transcen-dental equation i n terms of M a c h number to be solved for the Prandtl-Meyer angle and the Blasius series solution lead to a nonlinear two point boundary value problem. T h e n came the desk calculator which enabled the computation of more complex flows such as supersonic nozzle flow by the method of characteristics. The first electronic computers immediately enabled the solution of the in i t ia l value two point boundary value problem for two-dimensional boundary layers and of the potential equation for flows wi th rather arbitrary boundary shapes, a substantial step forward. Krause [11] 1 Chapter 1: Introduction and Literature Review 2 has compiled a short history of early works in this field. W i t h the advent of high speed digital computers, researchers started developing numerical techniques for the solving hitherto unsolvable differential equations. The basic idea behind most of these techniques is to try to 'simulate' the physical phe-nomena in some approximate way and to improve the approximation in successive attempts. Computat ional fluid dynamics has always been at the forefront of the de-velopment of these numerical techniques [12]. Several reasons have contributed to this leadership among many disciplines some of which are: • the inherently nonlinear behavior of fluids (e.g., convection, turbulence) which must be accounted for, • the mixed hyperbolic-elliptic character of the governing part ial differential equa-tions, and • the need for engineers designing new airplanes to obtain much more accurate performance estimates and therefore more accurate results for the simulation than, say, their colleagues designing bridges. Another reason that follows from above is the size of typical problems. A problem in structural mechanics is considered 'large' if it exceeds 5 x 10 3 nodes, whereas large in fluid problems means more than 5 x 10 5 grid points [12]. For about 25 years the computational fluid dynamics grew around F in i te Dif-ference Methods, as these were relatively simple to understand and code, easy to vectorize and the structured grids typically associated wi th them described appro-priately the simple geometric complexity of the fields that were solved. However, as the computers became bigger and faster, attempts were made to simulate more and more complex flow domains and it soon became clear that structured grids were not flexible enough to describe these domains. It was at this point in time that unstruc-tured grids and Fini te Element Methods—a natural way of discretizing operators on Chapter 1: Introduction and Literature Review 3 them—entered the scene. Since then the finite element methods have become a de jure standard techniques, although some researchers s t i l l prefer the finite difference methods. Since the in i t ia l development efforts for the finite element methods for fluid me-chanics, much attention has been devoted to the treatment of the primit ive variable (velocities and pressure) formulations. The dominant approach used in the finite el-ement methods is the mixed order interpolation scheme which attempts to eliminate the spurious pressure modes [31], i.e., the tendency to produce 'checkerboard' pres-sure distributions [24]. This scheme is analogous to the 'staggered gr id ' used in the finite difference methods. However, mixed order interpolation is not total ly effective at el iminating the spurious pressure modes for simple elements, i.e., bilinear triangle or bilinear quadrilaterals where constant element pressures are used [24]. De Bre-maecker [5] has suggested use of penalty functions to avoid this problem; Rice and Schnipke [24] have presented an equal order velocity pressure interpolation that does not show the spurious pressure modes at a l l ; Zinser and B e n i m [37] have suggested a segregated velocity pressure formulation of Navier-Stokes equations which results in an approximate pressure equation and allows uncoupling of pressure and velocity solutions. Numerous other suggestions have been made. For the present investiga-t ion, we have used a mixed order interpolation wi th quadratic bilinear element where pressure is not constant over the element; along the element boundaries the variation is linear whereas in the domain the variation is linear plus one quadratic cross term. This eliminates the above mentioned 'checkerboard' pressure distributions. 1.2 Stenosis and Pulsatile Flow The world "stenosis" is a generic medical term which means a narrowing of any body passage, tube, or orifice. Thus, an arterial stenosis refers to a narrowed or con-stricted segment of an artery and may be caused by impingement of extravascular Chapter 1: Introduction and Literature Review 4 masses or due to intravascular atherosclerotic plaques which develop at the walls of the artery. In the arterial systems of humans, it is quite common to find narrowings, or stenoses, some of which are, at least approximately, axisymmetric or 'collar-like.' Pa r t i a l occlusion of blood vessels due to stenosis is one of the most frequently oc-curring abnormalities in the cardiovascular system of humans; it is a well established fact that about 75% of al l deaths are caused due to circulatory diseases [2]. Of these, atherosclerosis is the most frequent. Considerable studies related to the circulatory flow in blood vessels were given major attention towards the beginning of the present century. Quite a few analytical and experimental investigations related to blood flow w i t h different perspectives have already been carried out. Interest in studies of this particular domain of biomechanics has increased wi th the discovery that many car-diovascular diseases are associated wi th the flow conditions in the blood vessels; one major type of flow disorder results from stenosis. Diagnosis of arterial diseases is usually carried out by the method of X-ray angiog-raphy which involve considerable risk of morbidity [3]. So any attempts to develop non-invasive techniques for detecting arterial disease like stenosis are very worthwhile. As a result, numerical modeling of arterial blood flow has long been of interest, but development of realistic models has occurred only in recent times [22]. Causes of the slow progress include nonlinearities in the model equations, unsteady flow in the vas-cular system, and elastic properties of arteries. Most of the studies have been carried out assuming that blood is a Newtonian fluid. However it has now been accepted that blood behaves more like a non-Newtonian fluid under certain conditions, in particular at low shear rate. Srivastava [28] has presented a study of the effects of stenosis on the flow of b lood when blood is represented by a couple stress fluid, a special type of non-Newtonian fluid which takes particle size into account. This , however, is beyond the scope of the current investigation. Young [35] has presented a survey of some of the early works done in this field. Chapter 1: Introduction and Literature Review 5 Over the years, Daly [4], Deshpande [6,7], Yoganathan et al. [34], Phi lpot et al. [21], Porenta, Young and Rogge [22], Forrester and Young [8,9] and others have shed considerable light on the complex problem of blood flow through an artery in the presence of stenoses. Daly studied the effects of pulsatile flow through stenosed canine femoral arteries for constrictions i n the range 0-61% and presented a comparison of in-vivo and in-vitro results. Local flow reversal was observed in the wake of the stenosis during systole and during diastolic flow reversal. He also concluded that the rate of development of the local reverse flow is very sensitive to the height of the stenosis, a result supported by Forrester as well. Deshpande obtained numerical solutions for steady flow through axisymmetric, contoured constrictions in a rigid tube uti l izing full Navier-Stokes equations i n cyl in-drical coordinates. Laminar boundary conditions, normally applicable only at x = ± o o , were applied at finite values of x by using a transformation of the type r\ = tanh(Kx) , where is K is a constant chosen to group the grid points i n an efficient manner. The stenosis geometry was modeled as a full cosine wave. Yoganathan et al. obtained steady flow velocity measurements in a pulmonary artery model with varying degrees of pulmonic stenosis. Phi lpot et al. used a 7-mW He-Ne laser as the light source for flow visualization in a Pulmonary artery model. They photographed the setup over the entire systolic period, thus observing the development and dissi-pation of the flow. Porenta et al. developed a one dimensional analysis which included taper, branches and obstructions and used Kantovich-like numerical technique and yielded stable nu-merical solution after 3 cycles, independent of in i t i a l values. The arterial cross-section was of the form a • exp(—ax), where a is a small positive number. Low frequency behavior was approximated wi th seven cosine and seven sine components of the gen-eralized Fourier series. Propagation of a pressure wave was evident from his results. Forrester modeled the stenosis as a full cosine wave and used the Karman-Pohlhausen Chapter 1: Introduction and Literature Review 6 method as described by Schlichting [26] to model the flow using a five parameter fourth order polynomial and obtained an analytical solution. He has also reported some experimental results which supported his solution scheme. In particular, his experiments clearly reveal the development of separation downstream of the stenosis and the reattachement points as well as variation of incipient separation point and crit ical Reynolds number wi th height of stenosis. 1.3 Analysis in Time Domain Pulsati le flow in the arteries is highly variable wi th time and to understand al l the effects of stenosis, one has to solve the full Navier-Stokes equations for any time t. In general, this can be done by obtaining the finite element solution for the steady part and then carrying out simultaneous numerical integration, choosing a sufficiently small step so as not to allow divergence of the results. This , however, is extremely costly in terms of C P U time and storage since, at each step, we must solve for the effects of the quadratic, nonlinear, convection terms. McLach lan [14] and Nayfeh and Mook [17] have pointed out that the solution of the differential equations wi th such nonlinearities has a steady part as well as an oscillating part. The techniques available for determining the steady state behavior and the slowly varying periodic behavior of forced nonlinear systems can be divided into two groups: the first group includes the method of averaging and multiple scales. W i t h this group one determines the equa-tions describing phases and amplitudes which are to be solved simultaneously. The second group includes the Lindstedt-Poincare technique and the method of harmonic balance. Details of these methods are presented in [14,17]. Pat tan i [19,20] has modified the method of averaging and shown that it can be used successfully for this k ind of problem. Using the Floquet theory, he has also shown the solution to be stable. This method is used here. Chapter 1: Introduction and Literature Review 7 1.4 Scope of Present Investigation This thesis describes the development work on extending the finite element method to cover the effects of stenoses on pulsatile flow in arteries. The problem configuration of interest is that of an axisymmetric, or 'collar-like' artery. The pulsatile flow is approximated as a combination of steady, sinusoidal and cosinusoidal velocities. In Chapter 2, the fundamentals of fluid flow i n cyl indrical coordinates are pre-sented. The chapter begins by stating the governing differential equations for the fluid, derived in [26] based-on the principle of conservation of mass and momen-t u m , and choosing non-dimensional parameters to effectively characterize the flow through a pipe. Kinemat ic and traction boundary conditions are also presented and non-dimensionalized. In the last section, the Poisson's equation for stream function is derived. Chapter 3 is about discretization of the various differential equations and the boundary conditions using a family of locally defined shape functions. Residuals are formed and Galerkin's procedure is used to minimize those residuals. The modified method of averaging is used to solve for the steady streaming part and slowly varying periodic parts. A Newton-Raphson iteration scheme is outlined for the solution of the resulting set of nonlinear, algebraic equations. Numerical results are described in Chapter 4. Besides the pulsatile flow through stenoses, several other l imit ing problems are also solved to check the code. These results compare favorably wi th published analytical and experimental results in the literature. Streaklines and streamlines are plotted for different values of the non-dimensional parameters Re and Rw. In the last section, conclusions and suggestions for further investigations are presented. C H A P T E R 2 Derivation of Governing Equations In this chapter, the governing equations for axisymmetric fluid flow are presented. The velocity-pressure primitive variable form of the Navier-Stokes equations for two dimensional incompressible viscous flow is used. The equations of motion are non-dimensionalized to effectively characterize the flow situation. Final ly , an approximate representation for pulsatile flow is presented. A l l the analysis is based on incompress-ib i l i ty assumptions which hold well for liquids. 2.1 Conservation Equations and Boundary Conditions The basic equations governing viscous, incompressible fluid flow are the conservation equations: conservation of mass and conservation of momentum. These equations along w i t h the specified boundary and ini t ia l conditions, in general, need to be solved over a domain f i bounded by a contour T which is composed of two distinct parts— kinematic boundary and traction boundary—denoted as Tu and T s , respectively, as shown in figure (2.1) If r, <p and z denote the radial , azimuthal and axial coordi-nates, respectively, of a three-dimensional system of cylindrical coordinates and vT, v^, vz denote the velocity components in the respective directions, then the governing equations for incompressible fluid flow are given by Schlichting [26] as 8 er 2: Derivation of Governing Equations r u Figure 2.1: Fluid Domain (dvr dvr Vfidvr V&2 dvr\ 0[lT  + v'-dT-  + 7-d+-T +  v*-d7)  = „ dp (d2vT 1 dvr vr 1 d2vr dr dr2 r or 2 dv* a y r 2 ' r 2 d<f>2 r2 d<f> dz2 ^ dt  T dr r d<j> r  2 dz J * r5<A V dr2 + r dr r2 r2 dcf?2 r2 d<f> dz2 (dvz dvz [-d7  + V'-dr-V4, dvz r d<j> dvz + — — + "z o dz Chapter 2: Derivation of Governing Equations 10 The stress components assume the form crr = -p + 2n — or 10 f ldv<f> ,  VA <r, = - P + 2u [--Qj + j ) „2 = -p + 2n— Tr4> = V dr \ r J r d(p /<9ur <9u2\ (2.5) where, p is the fluid density, v = — is the kinematic viscosity, p is the fluid pressure, Fr, Fj, and Fz are the body forces in the radial, azimuthal, and axial directions, respectively, crr, CT^ and crz are the stress in the radial, azimuthal and axial directions, respectively, and Tr<j>, Tj,z and rrz the corresponding shear stresses. Equations for axisymmetric flow are obtained from equations (2.1-2.5) by letting 1. U0 = 0, 2. § = 0,and 3. x, y refer to axial and radial directions, respectively. Furthermore, for the present investigation the body forces are assumed to be zero. This results in the following two momentum equations for the x and y directions, D d d d . , , where — = — + u——\- v— is the total material derivative. Dt dt ox dy Du Dt Dv Dt pdx \dx2 d2u d2u 1 du^ dy2 ydy. I dp + v pdy \dx2 d2v d2v I dv vN + ^ + " ~ o (i,y)eO, t>0 (2.6) dy2 p dy y2 Chapter 2: Derivation of Governing Equations 11 Equat ion (2.4) reduces to du dv v ox oy y (z,y)6f2, t>0 which is the equation of continuity for axisymmetric flow. . .,. d (du dv v\ , d (du dv v' A d d i n g I — h -5—I— and ^ ~ + 7 T + ~ ox \ox oy yJ oy \ox oy yt sides of the 1st and 2nd of equations (2.6) we obtain (2.7) to the right hand Idp p dx Du _ ~Dt Dv _ Idp Ut  = ~o"dy  + V ^ 2u d 2u d2u 1 (du <9t>N dx 2 dy 2 dxdy y \dy dx. d 2v d 2v d 2u 2 + 2 — + + -'dv dx 2 dy 2 dxdy y \dy y (x,y)eft, t>0. (2.8) Equations (2.8) are referred to as the Navier-Stokes equations and are preferred over equations (2.6) since the latter lead to symmetric matrices when integrated over a finite element. Of the total six stresses given by equations (2.5), only the following three are non-tr ivial . These are also referred to as constitutive relations. • xy n du rt dv _ (dv du ^ \dx dy, (2.9) Boundary conditions are given on the velocity components over Tu, referred to as the kinematic boundary and on traction over Ts, referred to as the mechanical boundary or traction boundary. These conditions are u = U , u = V (x,y)er„, t > 0 o~xnx + Txyn2 = X (x,y)eT3, t>0 Txynx + o-yn2 = Y (2.10) Chapter 2: Derivation of Governing Equations 12 Substituting equations (2.9) into equations (2.10), we obtain the boundary conditions in terms of pressure and velocity components. „ du (du dv ^ \dy dx ni + u = u, v = y (x,y)eru, o o 'du dv^ ^ {dy  + dx, -p + 2fi dv dy n 2 = X n2 = Y (ar,y)€r„ O O where, U , V are the specified velocities on Tu, X , Y are the specified tractions on F3, and n 1 ; T i 2 are the direction cosines of the unit outward normal to the boundary. (2.11) 2.2 Non-dimensional Form of Governing Equations The Navier-Stokes equations, the continuity equation and the boundary conditions can be non-dimensionalized by choosing suitable characteristic values for the coor-dinates, the velocities and the time-scale variables and defining non-dimensional pa-rameters. Choosing the non-dimensional variables as x _ y x = — y = — t = ut Ro Ro u v „ p u = — v = — p = UQ UQ (2.12) fj.u0/Ro where, RQ, U0 and u>~ 1 are the characteristic values of the coordinates, the velocities and time-scale, respectively, and introducing equations (2.12) into equations (2.8) and (2.7), we obtain the non-dimensional form of the Navier-Stokes equations and the continuity equation for axisymmetric flow. Chapter 2: Derivation of Governing Equations 13 „ du n ( du du\ dp R„ — + Re [U— + V— = - 7 p + dt \ dx dy J dx „ dv _ / dv dv\ dp Rw— + Re [u— + v— ) = + dt \ dx dy J dy ^d 2u d 2u d 2v 1 fdu dv^ dx 2 dy 2 dxdy y \dy dx/ d 2v ^d 2v d 2u 2 (dv v\ dx 2 dy 2 dxdy y \dy y I du dv v T + T + -  = 0 dx oy y (2.13) A l l quantities are their respective non-dimensional values and the tilde has been UJR 2 omitted for convenience. The two non-dimensional parameters, R^ = — - , the v UQRQ frequency Reynolds number and RE = , the steady Reynolds number, com-v pletely characterize the flow. It should be noted that the pressure has been non-dimensionalized wi th respect to the characteristic shear stress rather than the dynamic pressure as it is anticipated that the flow domain would be a slow, shear dominated one. Introducing equations (2.12) and the following additional non-dimensional vari-ables <7 = < T x <j = < T y f = T x y UUQ/RQ  Y UUQ/RQ  XY UUQ/RQ UUQ/RQ UUQ/RQ UQ UQ (2.14) into the constitutive equations (2.9) and boundary conditions (2.11), we obtain ndu a x = - p + 2d-x ndv <Ty = - p + 2d7l (2-15) dy _ dv du Txy = [fa  + dy', Chapter 2: Derivation of Governing Equations 14 and du (du dv' \dy dxt nx + du dv \dy dx) u = U , v = V (x,y)eTu, £ > 0 n 2 = X dv (x,y)eTs, t>0 (2.16) n 2 = Y where, again, a l l the quantities are their respective non-dimensional values and the tilde has been omitted for convenience. 2.3 Stream Function and Poisson's Equation Equations (2.13) must be solved for the boundary conditions (2.16) i n order to find the velocities and the pressure at any point in the domain f2. The first difficulty en-countered by any solution method arises naturally from the nonlinear non-self-adjoint convective terms i n the Navier-Stokes equations. A further significant difficulty aris-ing in al l numerical formulations of the governing equations in primitive variables (u ,u,p) representation is due to the interaction of the pressure and dilatation terms. This difficulty arises from the fact that the continuity equation merely represents a constraint on the divergence of the velocities rather than a full th i rd equation cou-pling the pressure wi th the velocities. As a result, it is difficult to develop solution methods which are capable of determining the primitive variable (w,u, p), simulta-neously, wi th the same degree of accuracy—the pressure cannot be obtained equally accurately as the velocities—when a single restricted variational functional is con-structed as the basis of formulation. A n y effort to do so leads to spurious pressure modes, i.e., 'checkerboard' pressure distributions. Furthermore, an exact variational formulation is non-existent due to the presence of the non-self-adjoint, nonlinear, con-vective terms. Tuann and Olson [31,32] have shown that the only appearance of p in a restricted variational functional is coupled wi th the dilation term which is one order less than the velocities. Chapter 2: Derivation of Governing Equations 15 One way to avoid this difficulty is to modify the formulation of primitive equa-tions, so that either the system consists of the Navier-Stokes equations and an induced equation which involves pressure in a relatively dominant role, or else, the diver-gence constraint is satisfied explicit ly and exactly, dropping the primitive variables approach. Alternatively, one fourth order equation can be obtained for the stream function. Relative merits of these formulations in avoiding the difficulties arising out of incomplete coupling between (u,i>) and p are discussed i n [31]. Another way to satisfy the continuity equation exactly is to define the velocities i n terms of yj, the stream function, as ldyj u = - — y oy _ld_V> y dx (2.17) Hence, dx dy y dx \y dy J dy \ y dx J y \ y dx) 1 d 2yb t 1 dyj I d 2yj 1 dtp y dxdy y2 dx y dxdy y2 dx = 0, i.e., the continuity equation is satisfied point-by-point. The vorticity £ is given as du dv , Substituting equations (2.17) into (2.18), we obtain i fd2i> a2_0_ia_v>\ y\dx2 dy2 y dy) ~ du_dv  (2- 19) dy dx This is called Poisson's equation and is.used to compute the stream function, once the pr imit ive variables have been determined. (Right hand side is known.) C H A P T E R 3 Finite Element Formulation In this chapter, the Navier-Stokes equations and the continuity equation (2.13) along w i t h the boundary conditions (2.16) are discretized using the finite element method and a set of locally defined polynomial shape functions to represent the primit ive vari-able u, v and p. Galerkin's procedure is used to minimize the error in approximation by making the base vectors orthogonal to the residuals. Discretization of boundary conditions is also considered. A steady state periodic solution is obtained using the modified method of averaging and a Newton-Raphson iteration procedure is outlined for solving the resulting set of nonlinear algebraic equations. A truncated Fourier series is used to approximate the unsteady solution as a slowly varying periodic so-lut ion. The Galerkin's procedure is also used to find a representation of pulsatile flow. Brief details of computation of stream functions and plott ing of streamlines and streaklines are also presented 3.1 Discretized Form of Governing Equations Olson and Tuann [31] have shown that the finite element interpolation for the pres-du dv sure should be no higher than that for the dilation term ——I- TT " in order to avoid dy dy spurious singularities for the (u,v,p) integrated formulation. Therefore, the finite element interpolation function for the pressure p should be at least one degree less than that for the velocity components (u,v). The highest order of derivative in the equations (2.13) is two for the velocity components (u,v) and one for the pressure 16 Chapter 3: Finite Element Formulation 17 p; as a result C 1 continuity is required of the interpolation functions for (u,v). This difficulty can, however, be avoided by integrating the first two of equations (2.13) by parts to reduce the highest order of derivative to one for the velocity components (u, v) and to zero for the pressure p. Now, only C° continuity is required of the inter-polat ion functions for equations (2.13). Furthermore, for the present investigation, isoparametric elements are used, i.e., the same shape functions are used for inter-polating the variable (u,v,p) and for transformation from (£,77) natural coordinates to {x,y) coordinates. Taking these conditions into consideration, curved isopara-metric elements shown in figure (3.1), wi th quadratic interpolation for the velocities and bilinear interpolation for the pressure are used to carry out the finite element discretization of the governing equations. The velocities u, v represented by u = NiUi{t) v = NiVi(t) 8 where Ni are the quadratic isoparametric shape functions given as N7 = N5 = N, = 4< 1 -0 ( l - ' 7 ) ( l + * + '7) - 1 ( 1 + 0 ( 1 + »7)(1 N6 N8 N, N2 - 0(1 + + ( - V) | ( i + 0 ( i -v2) ^ ( l - O ( W ) . (3.2) Chapter 3: Finite Element Formulation Figure 3.1: Isoparametric Element (a) Element in (£,77) space, (b) Element (x,y) space. Chapter 3: Finite Element Formulation 19 The pressure p is represented by p = MjPj(t) ;' = ! , . . . , 4 where Mj are the bilinear isoparametric shape functions given as Ma = 1(1-0(1-17) M2 = 1(1 + 0 ( 1 - 1 / ) M 3 = i ( l + 0 ( l + » 7 ) M4 = 1 (1-0 (1+ .7). (3.3) Substituting the interpolation functions (3.2) and (3.3) into the governing equa-tions (2.13) and the boundary conditions (2.16), we obtain Tlx = Rw—g^-Uj + Re " dNk ^MdNk dy dMj dN, 'd 2Ni d 2Nj d 2Ni 1 (dNj R  dN'v +R Tiw—QJ-Vj + lte ' dNk ^MdNk dMj d*Nj 2 d 2 N j v  9 2 N j u - - ^• dx 2 1 dy2 1 dxdy  j y \ dy  j y  j (3.4) 72, = dNj dNj Nj _ t t. + — V j + 7 2 n = ndNj j P j + ~dTUj ni + 'dNj dNj  y + I T * n 2 — X KV2 = dNj dy •Uj + -dx~V> n x + n2-Y where 7Z\, 722 and 72.3 a r e the errors in approximating the x-momentum equation, the ^-momentum equation and the continuity equation, respectively; 72-ri a n d 72p2 are the errors in approximating the x-component and y-component of the traction boundary condition, respectively. These are called residuals and are identically zero Chapter 3: Finite Element Formulation 20 only when the interpolation functions are capable of representing the exact solution. To minimize the errors, we make the residuals orthogonal to their base vectors—i.e., the shape functions JVt- and Mi—over the finite element. This leads to II K1Nidto+ [ KnNidT = 0 JJ w J r» // K2Ni<m+ I KriNidT = 0 JJ u e J r» JJ n^MidSl = 0, (3.5) where f i e is the domain and Tes is the traction part of the boundary of the element under consideration. Substituting equations (3.4) into (3.5) and integrating the ap-propriate terms by parts, we obtain 0 0" Uj Kij Uuv Kij -Pij ' Uj " 0 0 < Kji Kij y -Pij < Vj 0 0 0. .Pj. .-Pji 0 .Pj . '  sfjkuJuk + 8ijkvjuk ' 0 ' + Re < SijkUjVk + tijkVjVk • = < 0 -o . 0 . (3.6) Chapter 3: Finite Element Formulation 21 where dot denotes differentiation w i t h respect to time and the arrays are //,. 'dNidNj i ^dN^d^ 2_ i,j,k = 1,...,8 Luu 'mm, + ^ dn n e \ ax c7z dy oy J = 11 Qs \ ox ox oy oy yz 1 V.uv °ijk JJ n e dx dy = If wf-r JJ ue dx IIQENTN3 dy an * - IL ( ^•Mi + -NiMA dQ dx y 7 i = l , . . . , 8 j = l , . . . , 4 . (3.7) Vector f on the right hand side of equation (3.6) is the consistent load vector and is given as f = < ( f Xiv.dr-) J r; / Yiv4-dr J r« 3.2 Boundary Conditions For a l l the nodes that are on the kinematic boundary Tu, the nodal values are forced to the specified values of the corresponding variables. Thus, in the present investigation, the kinematic boundary conditions are satisfied exactly—node-by-node, that is. As is evident from equations (3.5), the traction boundary condition, i n general, is satisfied only in the "mean"— TZn and TZr2 are not identically zero. The specified stresses X and Y are integrated over the traction boundary, Tes, to give the consistent load vector f. Chapter 3: Finite Element Formulation 22 The computer program used for the present investigation is set up for homogeneous traction boundary conditions, X = 0 and Y = 0 only. Consequently, the consistent load vector f = 0 and the traction boundary conditions are approximately satisfied. Equat ion (3.6) reduces to [M]{d} + [K]{d} + i2e< ' sfjkuJuk + Sljkviuk ' 6fjkujVk + Sijkvjvk \ = 0 0 (3.8) where, M is the mass matr ix given by 0 0" M = 0 0 0 0 0 K is the stiffness matr ix given by Kij •L.UV Kij -Pij ' K = Uuv Kji Uw Kij -Pi 0 and is the nodal vector of unknowns. d = u. [Pj ) 3.3 Steady State Periodic Solution Due to the presence of the quadratic nonlinear terms in the Navier-Stokes equations, the velocities as well as the pressure have a steady component and a time-dependent Chapter 3: Finite Element Formulation 23 component. One way to obtain the complete solution is to find a steady state finite element solution and then use simultaneous time-scale integration of the governing equations. Thi s , however, is very costly and requires large storage and C P U time. Quite often the purpose is well served wi th an approximate slowly varying periodic solution. Nayfeh and M o o k [17] have classified the available techniques for deter-mining the steady state behavior of the forced oscillation of a nonlinear system into two groups: the first group includes the method of averaging and multiple scales and the second group includes such methods as Poincare technique and the method of harmonic balance. Pa t tan i [19] has modified the method of averaging, as given by Nayfeh and Mook [17], to obtain the steady streaming part, as well as the slowly varying pe-riodic part, of the solution from equations (3.8). The starting point in this method is to assume a form of the solution as d = A + B ( * ) c o s * + C ( * ) s i n i (3.9) where B ( i ) , C(t) are assumed to be slowly varying functions of non-dimensional time t. The first term A represents the steady streaming part of the solution which nat-ural ly arises for systems of equations wi th quadratic nonlinearities as is encountered here. Alternatively, equation (3.9) can be interpreted as the first three terms of the Fourier expansion of d. Differentiating equation (3.9), we obtain d = - B ( t ) s i n i + C ( t ) c o s i + B ( £ ) c o s i + C ( f ) s i n i (3.10) To obtain an autonomous system of equations governing the amplitudes A , B (£ ) and C ( i ) , equations (3.9) and (3.10) are first substituted into the equations (3.8). Then to obtain the three sets of equations for the three sets of unknowns A , B ( i ) and C(t), the resulting equations are 1. averaged over the period 27T, Chapter 3: Finite Element Formulation 24 2. mult ipl ied by cos t and averaged over the period 2w, and 3. mult ipl ied by s i n i and averaged over the period 27T. B y hypothesis, for small values of Re and R^, B ( i ) and C(t) vary much more slowly wi th time than the functions sin/j, C O S T ; and time t itself. Therefore, B(t), C ( t ) , B(t) and C(t) are considered to be constant while performing the averaging. This process yields "0 0 0 " ' A ' " K 0 0 " ' A ' 0 M 0 < B > + 0 K M < B • .0 0 M . . c . . 0 - M K. .c. (- 6fjk(A»At + BVBJI/2 + C y C J / 2 ) + 8 yjk(A»Al + + C « Q / 2 ) ) « % , ( A ? A £ + B ? B £ / 2 + C)Cl/2) + Sljk{K)A.l + B]B%/2 + C J Q / 2 ) 0 -^•*(AJBJ| + AJfBJ) + «? i f c (AJBH + A £ B ? ) . + i? e < ^ * ( A J B j ; + AJBJ) + Sfik(A)B% + AlB]) 0 % ( A ? C £ + AJCJ ) + - ^ ( A J C X + AJfCV) ^ ( A j q ; + AJCJ) + ^ f c ( A y c j s + A j q ) 0 ' 0 < 0 . 0 (3.11) Chapter 3: Finite Element Formulation 25 where A = < r A y Ay B = <! ( B y l By B p i J ( c y ) cy y. are the average values over one period (27r) of A , B(rj) and C ( i ) . The superscripts u, v and p indicate that the amplitudes are associated wi th the u, v and p degrees of freedom respectively. The steady state solution corresponds to the singular points of the autonomous system of equations when B = C = 0. This results in a set of nonlinear algebraic equations for A , B and C which are solved using the modified Newton-Raphson itera-t ion procedure as outlined i n Appendix A . The equations obtained for the increments to the solution vector are of the form [T] {Ax} = {-f} (3.12) where, [T] is the tangent stiffness matr ix , { A x } is the incremental solution vector and {—f} is the unbalanced load vector. Details of equations (3.12) are given in Appendix B . The above matrices are calculated for each finite element and are assembled into the corresponding global matrices. This is carried out by programs A X I L I N and A X I N O L I N on the S U N 4/260 computer system and the resulting set of sparse, linear, algebraic equations are solved using the University of Waterloo solution package S P A R S P A K . 3.4 Characterization of Arterial Fluid Flow In the application of the principles and concepts of fluid mechanics to the stenosis problem, several general characteristics play a major role. These include the Theo-logical properties of the fluid, the nature of flow (e.g., steady or unsteady, laminar Chapter 3: Finite Element Formulation 26 or turbulent) , the geometry of stenosis and the distribution of flow across the cross-section after it has sufficiently "stabil ized" over a long development length. B lood consists of formed elements, red blood cells, white blood cells, platelets etc. suspended in plasma. Al though plasma is a Newtonian fluid, Young [35] has reported that blood exhibits non-Newtonian behavior at low shear rates. However, at higher shear rates commonly found in larger arteries, blood is generally assumed to behave as a Newto-nian fluid. The scope of the current investigation is to analyze the effects of stenosis on the fluid flow, without particular reference to the actual artery stenosis caused due to the impingement of extravascular masses, or due to the intravascular atherosclerotic plaques; therefore no effort is made to further justify the Newtonian fluid assumption. Rather, an attempt is made to represent the actual velocity distr ibution across the artery cross-section. Ar ter i a l blood flow is highly pulsatile w i t h the instantaneous flow rate varying over a wide range during a flow cycle. Young [35] has reported representative waveforms for the left circumflex coronary artery and the femoral artery which are reproduced here i n figure (3.2). These are often approximated by the following equation. Uo(y,t) = a(y) + (b(y)cost)'H(t) where, 1 0 < t < 7T, H(t) = I 0 7T < t < OCK, and H(t + a) = H(t). To get r id of the discontinuity arising from H(t) in figure (3.3), we approximate Chapter 3: Finite Element Formulation 27 Figure 3.2: Typical arterial flow waveforms (Young, 1979) Chapter 3: Finite Element Formulation 28 Hit) 1 _ 0 IT OCT (a + l)7T 2a7T ' ( 2 a + l ) 7 T Figure 3.3: Function H(t) UQ(y, t) as Uo{y,i) ~ A0(y) + B0{y)cost + C0(y)sint ' My)' ' l " T 'i4o(y) ' (1 cos t s'mt) < B0(y) > = < cos t > < B0(y) . c0(y). . s int , . C0(y) . A0(y), B0(y) and Co(y) are determined by minimizing the error in approximation by Galerkin's procedure. The resulting residual ' 7lu = AQ(y) + B0(y) cos t + C0(y) sin t - [a(y) + (b(y) cost )- H(t)} is made orthogonal to the base vectors 1, cost and s int . Thi s gives ' 1 ' ' 1 " T 'A0' ' 1 ' ' 0 ' < cos t > < cos t • dt< Bo • - / < ./ 0 cos i - (a + (6cosi ) • H(t)) dt = < 0 • k s int , > s int , . Co . , s int , . 0 . Chapter 3: Finite Element Formulation 29 or, /ai 0 " 1 cost sin t COS t cos2t sin t cos t d t < B0 > smt sin t cos t sm2t . Co . /ai 0 or, Q7T s m a i OJ7T sin 2a7r sin car —— H 2 4 sin a?r — COS OJ7T 1 COS t s int — cos OCK s i n 2 a x air sin 2a7r _ cost fir - d t + b / < J 0 cosH , sin t cos 11 \ d t Bo > = a < Q7T sin air — COS OJ7T 2 2 4 Equat ion (3.13) can be solved for different values of a. In particular, . 37T tor a = —, 2 ' 2 1 0 ir K 2 ' (3.13) ' Mv)' '1 - ' 0.433897 ' < B0(y) • = a(y) < 0 . + b(y) < 0.0446932 . C 0 (y ) . . 0 . k 0.675718 t for a = 2ir, (Mv Bo(y Co(y and for a = Sir, > = a(y) < • + b(y) < 1 ir 0 1 ' 2 ' ' Mv)' ' 1 " '0.18091221 ' < B0(y) » = a(y) < 0 . + b(y) « 0 . C 0 (y ) , . 0 . . 0.29494257 , B o t h a(y) and b(y) represent the fully developed solutions for the steady state laminar flow and cosinusoidal laminar flow, respectively. Schlichting [26] has given Chapter 3: Finite Element Formulation 30 y the exact solution for both these parameters. a(y) is a quadratic function of — and b(y) is a function of > — w h e r e R is the radius and J 0 is Bessel's function of the first k ind of order zero. For ease of programming, though, both a(y) and b(y) are assumed to vary parabolically across the cross-section. To compensate for this approximation, enough entrance length is provided to let the flow develop naturally. 3.5 Streamlines and Streaklines A line i n the fluid whose tangent is everywhere parallel to the fluid velocity u(x,y,t) instantaneously, is a streamline. The family of streamlines at any time t are solutions of the equation 1 dx I dy y u(x,t) y v(x,t) where u and v are the component of u(x, y, t) velocity parallel to the axes x and y, respectively. The path of a material element of fluid does not, in general, coincide wi th a streamline. The pathline and streamline are the same only when the fluid motion is steady, i.e., only the A component is nonzero. A streakline is that line on which lie those fluid elements that, at some earlier instant, passed through a certain point of space. Thus, when a marking material is dropped in the flow, the visible path produced thus is a streakline. The family of streamlines and streaklines are found by solving the Poisson's equa-t ion (2.19). After the governing equations for the primitive variables (u,v,p) have been solved, the vorticity du dv dy dx can be computed at any point. In order to be consistent, the finite element in-terpolation functions chosen for the stream function ?/> to discretize the Poisson's equation (2.19) are the same as those used for the velocities. This also simplifies the computation and setting up of the matr ix equation. The stream function is Chapter 3: Finite Element Formulation 31 represented by tb = Ntfi i = l , 2 , . . . , 8 (3.14) where Ni are the shape functions given by equations (3.2) and yji are the nodal values. Substituting equation (3.14) and _ du dv dy dx dNj dNj dy  3 dx into equation (2.19), we obtain 1 TU = -d 2Nj d 2Nj _ IdNj dx 2 dy 2 y dy dNj dy dN< where TZ^, is the residual due to error in approximation. To minimize the error in approximation, we make 7Z^ orthogonal to the base vectors, i.e., the shape functions Jj' Tl^NidQ. = 0 or, JJ n ey 1 \d 2Nj d 2Nj 1 dNj dx 2 dy 2 y dy dNi dNj  dy  i dvj Integrating the appropriate terms by parts and simplifying, we obtain NW = (Pi) where the arrays are NidCl = 0. (3.15) * - II 'ldWdN± 1 dNj dNj y dx dx y dy dy dn 'dNj n e I dy  U^ dx  VJ dN4 NidCl + I -gNdT J r%y ij = 1 , 2 , . . . , 8. Cle is the domain of a single element, Tes is the traction part of the boundary of the element and 9 = d^_ dn dN3 dx dNj n 2 Chapter 3: Finite Element Formulation 32 is the tangential derivative of tp specified on the traction part of the boundary r*. The above process is repeated for each element and the matrices are assembled into the corresponding global matrices. The resulting set of sparse, linear, algebraic equations are solved using the University of Waterloo solution package S P A R S P A K . Substitut-ing back the solution i n equation (3.14), we obtain the finite element representation of the stream function ^ in the fluid domain. C H A P T E R 4 Numerical Results In this chapter numerical results are presented for different flow configurations. The results are obtained for various values of Re and in the viscous flow regime and, wherever possible, these results are compared wi th analytical , experimental or numer-ical results published in the literature. From the basic flow results, stream function values are computed and streamline and streakline plots are obtained. Streamlines are the stream function contours of the steady part of the velocity field and streaklines are the stream function contours of the total velocity field at one instant i n time. The steady state periodic solution of the problem corresponds to the singular points of the autonomous system of equations (3.11), that is, when B = C = 0. This results in a set of nonlinear algebraic equations in A , B, C. These equations are solved using the modified Newton-Raphson procedure as outlined in Appendix B . A l l the computer programs are implemented on the A m d a h l 5840 and on a depart-mental Sun 4/260 computer at the University of Br i t i sh Columbia . Double precision (Real*8) arithmetic is used throughout to reduce the effect of round-off errors. The resulting set of sparse, linear, algebraic equations are factorized and solved using the University of Waterloo solution package called S P A R S P A K . A 100 x 100 plot grid is used, in conjunction wi th the U B C surface visualization program *SURFACE, to obtain contour plots of streamlines and streaklines. Stream function values at the plot grid nodes are obtained by the procedure outlined in Appendix C . 33 Chapter 4: Numerical Results 34 4.1 Steady Flow Through a Pipe 4.1.1 Theoretical Results The steady flow through a straight pipe of circular cross-section is the simplest case wi th rotational symmetry for which an exact solution can be obtained. Let the x-axis be selected along the axis of pipe, figure (4.1), and let y denote the radial coordinate measured from the axis outwards. For fully developed flow the velocity component i n the radial direction, denoted by v, is zero; the velocity component parallel to the axis, denoted by u, depends on y alone and the pressure is constant at every cross-section. The governing equations (2.13) reduce to d 2u 1 du dp dy 2 y dy dx 2R Figure 4.1: Paral le l Flow wi th Parabolic velocity distr ibution The second and th i rd equations of (4.1) give p = p(x) and u = u(y). Rearranging the first equation, we obtain d 2u 1 du _ dp dy 2 y dy dx' ^ Chapter 4: Numerical Results 35 The left hand side of equation (4.2) is a function of y alone and the right hand side is a function of x alone; both these conditions can be true if and only if each side is equal to a constant. Hence, d 2u 1 du dy 2 y dy dp — = constant dx or, ld_ ydy du 'dy dp — — = constant dy or, du dy Spy2 , r faJ + C l or, For u to be bounded along the axis of the pipe, C\ = 0 and for u = 0 at y = R, R 2 C2 = — — • This gives u (y) dp' where \/Q~J  IS ^ n e pressure gradient which can be deduced from the given boundary conditions. 4.1.2 Finite Element Results We note that the exact solution obtained above is quadratic in y for the velocity component u and linear in x for the pressure p. This means that the quadratic isoparametric interpolation chosen for the current investigation should be capable of reproducing the results 'exactly'—up to the accuracy of the arithmetic used, that is. For this purpose, a single element, as shown in figure (4.2), is used to represent the fully developed flow. W i t h node 1 as the origin, we choose the radius R = 1 and length / = 4. The following kinematic boundary conditions were applied: Chapter 4: Numerical Results 36 • along x = 0: u(0,y) = u0(y) = (1 — y 2) and i> = 0; • along x = 4: i> = 0 and p = 0; • along y = 0: u = 0; and V • along y = 1: u = 0 and u = 0. 2/ u = ( l - y 2 ) u = 0 u = 0,v = 0 7^  4.0 t; = 0 p = 0 v = 0 Figure 4.2: One Element Representation of the Ful ly Developed Flow du The traction boundary condition along the centerline, y = 0, namely — = 0, was oy not enforced—it was left for the program to seek a satisfactory approximation to this dp condition. Also by letting u(0,y) = (1 — y 2 ) , impl ic i t ly we let — = —4. This wi th dx the condition that p = 0 at x = 4 gives p = (16 — 4x). The finite element results are given i n table (4.1). They match the exact solution up to the precision of the R e a l * 8 arithmetic used. No iterations were required to converge to the true solution. 4.2 Oscillating Flow Through a Pipe 4.2.1 Theoretical Results The case of the flow of a fluid through a pipe under the influence of a periodic pressure difference affords another example where an exact solution can easily be computed. This type of flow occurs, e.g., under the influence of a reciprocating piston, and its Chapter 4: Numerical Results 37 Table 4.1: Comparison of F in i te Element Solution w i t h Exact Solution Node Coordinates Variable F E solution Exact solution 1 (0.0,0.0) P 16.0000000001713 16.0 2 (4.0,0.0) u 1.0000000000046 1.0 4 (0.0,1.0) P 15.9999999999584 16.0 5 (2.0,0.0) u 1.0000000000055 1.0 6 (4.0,0.5) u 0.7500000000023 0.75 theory is presented in Schlichting [26]. It is assumed that the pipe is very long and circular i n cross-section. Under the assumption that the pipe is very long, the flow is taken to be independent of x . This means that the axial velocity component u is also independent of x , i.e., u = u(y,t). This simplifies the continuity equation to dv v 7T + - = 0 dy y or, or, ydy yv = constant = g(t) for some arbitrary function g{i). That is impossible unless v = 0. dp The y—momentum equation reduces to — = 0, i.e., p = p(x, t); the x—momentum equation reduces to j^9u _ dp ^ d 2u ^1 du dt dx dy 2 y dy (4.3) which is 'exact' as it implies no additional simplifications. The boundary conditions are u(y,t) = 0 for y = R at the wal l . Assume that the pressure gradient caused by the motion of the piston is harmonic and is given by dp - - = P o c o s i , Chapter 4: Numerical Results 38 where p0 denotes a constant. For this problem it is convenient to use complex notation and to put dp - nattr ibuting physical significance only to the real part (or, only to the imaginary part, if the pressure difference is sinusoidal). Note that the minus sign indicates that pressure drops i n the positive x direction. Assuming that the velocity function has the form u(y,t) = f(y)e lt, we obtain, % -e2 dy 2 Substituting equations (4.4) into the equation (4.3), and canceling the common term e lt from both left hand side and right hand side, we obtain the following differential equation for the function f(y): f"(y) + l-f\y) + (-iRMy) = -PO. To transform the differential equation above, we use f'(y) = v = yy(-iRJ)-dfiv) dn dn dy df(rj) f"(y) = (-iR.) dn d2f(y) drj2 ' Thi s gives, „2d2f(v) , df(n) 2 2 2f P o \ which has the particular solution fp(n) = ^—*^T^ a n d homogeneous solution A(7?) = AJ0(r)) + BYo(n), where A and B are, in general, functions of time only, and Chapter 4: Numerical Results 39 JQ and YQ are the Bessel's function of the first and second k i n d , respectively, of order zero. A d d i n g the two solutions and switching back to the y—variable gives For the solution to be bounded, B = 0 since YQ [yy/—iR^) is unbounded at y = 0. A is determined using u = 0 at y = R at the wall . Hence, = » [ / ( y ) - e " ] = 3? #u, \ Jo (Ry/^iRZ) The above can be further simplified by splitting the bracketed quantity into sine dp and cosine terms. Hence for cosinusoidal pressure difference — — = po cos t, ox u(y, t) = ^-(0 cos t + a sin t) Rw dp and for sinusoidal pressure difference — — = po s i n i , dx u(y, t) = a cos t + f3 sin i) Ru where a + i j Jo{y^iRZ)\ \ Jo (RV 1^) J ' A general discussion of the nature of the solution for arbitrary values of R^ is im-possible owing to the presence of Bessel's function wi th a complex argument, but for the two l imit ing cases of very low Ru and very large Rw, the solution above can be simplified. For very small values of R^, the above can be simplified by retaining only the first two terms of the series expansion Similarly using the asymptotic expansion of the Bessel function • J 0 ( z ) « , / A e - r l / 2 Chapter 4: Numerical Results 40 Figure 4.3: Flow in a pipe under oscillating pressure Chapter 4: Numerical Results 41 we can obtain an expression for very large values of R^. In general, the velocity distribution across a cross-section of the pipe under the influence of cosinusoidal pressure difference is as shown in figure (4.3). 4.2.2 Finite Element Results We choose a pipe of radius R = 1 and length / = 4 represented by a 5 x 8 grid as shown i n figure (4.4). We compute the fully developed velocity distr ibution using the series solution presented in the last section and specify these values at the entrance x = 0. The program seeks its own approximation to the flow downstream and stablizes before the exit x = 4. The velocity distributions thus obtained at the exit are compared wi th the exact solution. u = ( a cos t -+- (3 sin t) v = 0 J, u = 0, v = 0 1.0 v = 0 p = 0 4.0 •*-x Figure 4.4: F ini te element grid for oscillating flow through a pipe The exact solution presented in the last section is independent of the Reynolds number Re and hence we expect the numerical solution to converge without any iterations; this indeed is the case. Furthermore, we choose p0 = 1.0. Thus u(y, i) = - j - ( acos t + f3sint) Rw where a and /? are functions of y and are computed using the series solution. The program is run for three different frequency Reynolds numbers Rw = 1, 10 and 20. Chapter 4: Numerical Results 42 The results are compared in tables (4.2), (4.3) and (4.4). Table 4.2: Comparison of exact and finite element solution for oscillating flow for Re = 1.0 R = 1.0, po = 1.0 Rw = 1.0 y Exact Solution F in i te Element Solution a a P 0.0 0.2419938822 0.0454856408 0.2419887660 0.0455076871 0.125 0.2382643889 0.0445439938 0.2382631459 0.0445506155 0.25 0.2270649264 0.0417627709 0.2270629093 0.0417728360 0.375 0.2083631157 0.0372732688 0.2083631013 0.0372747261 0.5 0.1821068921 0.0312947819 0.1821058316 0.0313012409 0.625 0.1482273619 0.0241352748 0.1482278532 0.0241334887 0.75 0.1066428141 0.0161922692 0.1066424134 0.0161955928 0.875 0.0572639012 0.0079538898 0.0572645050 0.0079488352 1.0 0.0 0.0 0.0 0.0 Table 4.3: Comparison of exact and finite element solution for oscillating flow for Re = 10.0 R = 1.0, po = 1.0, Ru, = 10.0 Exact Solution Fini te Element Solution y a P a 0 0.0 0.0456409109 0.1109031291 0.0455297253 0.1108901990 0.125 0.0460493859 0.1091161974 0.0460123330 0.1091272587 0.25 0.0470648468 0.1037100366 0.0470122009 0.1037222211 0.375 0.0480518359 0.0945759605 0.0480434236 0.0945908185 0.5 0.0479364625 0.0816249082 0.0479068297 0.0816443909 0.625 0.0451974541 0.0649275362 0.0452096686 0.0649315983 0.75 0.0378803750 0.0449123607 0.0378715381 0.0449288345 0.875 0.0236620859 0.0226194555 0.0236821472 0.0225941990 1.0 0.0 0.0 0.0 0.0 Chapter 4: Numerical Results 43 Table 4.4: Comparison of exact and finite element solution for oscillating flow for Re = 20.0 R = 1.0, po = 1-0, Rw = 20.0 y Exact Solution Fini te Element Solution a a 0.0 0.0042213053 0.0601667061 0.0041474265 0.0600827801 0.125 0.0050090036 0.0598214600 0.0049651334 0.0598033790 0.25 0.0072867953 0.0586030834 0.0072308490 0.0585803199 0.375 0.0107517756 0.0559870567 0.0107261714 0.0560003541 0.5 0.0147491753 0.0511898489 0.0147068992 0.0512030557 0.625 0.0180745128 0.0433473344 0.0180846562 0.0433671069 0.75 0.0187627664 0.0318412776 0.0187522740 0.0318704518 0.875 0.0139564340 0.0168281658 0.0139112653 0.0167978695 1.0 0.0 0.0 0.0 0.0 The solution becomes increasingly wavy in the y—direction as Rw increases, and for this reason grid points lie closer to each other in this direction than along the pipe axis. The aspect ratio is 6.4 which is much less than 20 and hence no nu-merical difficulties arising out of the large aspect ratio were encountered. A t low frequency Reynolds number, the cosinusoidal component is dominant—one full order of magnitude higher than the sinusoidal part. However, as shown i n figure (4.5), wi th increasing frequency Reynolds number, the sinusoidal component increases whereas the cosinusoidal component decreases. At Rw = 10, both components are of the same order of magnitude, and at 7^, = 20 the sinusoidal component becomes the dominant one. F rom the velocity field results, stream function values at grid points were calcu-lated. These values were then extrapolated to a 100 x 100 grid using the defined interpolation shape functions and streakline plots are shown i n figures (4.6), (4.7), Chapter 4: Numerical Results 44 • ^  0.00 Q Cosinusoidal component — ~ Sinusoidal component 0.05 Veloci ty profi le 7 ^ = 10 0.15 0.00 O Cosinusoidal component . ~ —> Sinusoidal component 0.03 0.05 0.08 Veloci ty profi le Ru = 20 0.10 Figure 4.5: Variation of cosinusoidal and sinusoidal components for unsteady parallel flow for K = 1, 10 and 20 Chapter 4: Numerical Results 45 7T 37T (4.8) for t = 0, —, 7T and — . Twenty contours spaced equally between the instan-Lt Li taneous min imum and the instantaneous max imum values of the stream function are plotted. Note that a l l the contours are almost completely horizontal indicating that the finite element representation of fully developed flow under cosinusoidal pressure is very good. The instantaneous rate of flow, denoted by Q(t), is found by integrating the it—velocity across the cross-section of the pipe as shown below and is p r = 1 since u(y , t) = Q(t) = / u(y,t)-{2iry)dy J o = 2 ^ ( 1 , t) y dy Table 4.5: Instantaneous rate of flow through a pipe under cosinusoidal pressure difference R = 1.0, po = 1.0 Frequency Reynolds Number time Ru = 1.0 Rw = 10.0 i L 2 0 . 0 0 0.38126 0.10945 0.04144 7T 2 0.06357 0.16707 0.10706 7T -0.38126 -0.10945 -0.04144 3TT 2 -0.06357 -0.16707 -0.10706 (Negative instantaneous rate of flow means that the flow is to the left) Chapter 4: Numerical Results 46 Chapter 4: Numerical Results 47 Chapter 4: Numerical Results 48 Chapter 4: Numerical Results 49 4.3 Developing Flow in the Inlet Length of a Pipe 4.3.1 Theoretical Results As a th ird example of axisymmetric flow, we consider the case of flow in the inlet length of a straight circular pipe. A t a large distance from the inlet the velocity distr ibution across the cross-section becomes parabolic. We assume that the velocity i n the inlet section is uniformly distributed over the circular cross-section of radius R and that its magnitude is U o . Owing to viscous friction, a boundary layer wi l l be formed on the wal l and the depth of the boundary layer w i l l increase in the downstream direction. A t the beginning, i.e., at small distances from the inlet section, the boundary layer w i l l grow in the same way as it would along a flat plate at zero incidence [26]. The resulting velocity profile w i l l consist of an axisymmetric boundary layer near the wal l joined in the center by an area of constant velocity. Since the volume of the flow must be the same for every cross-section, the decrease i n the rate of flow near the wal l which is due to friction must be compensated by a corresponding increase near the axis. Thus the boundary layer is formed under the influence of an accelerated external flow, as distinct from the case of a flat plate. A t larger distances from the inlet section the central region of constant velocity asymptotically shrinks to zero, and finally the velocity profile is transformed into the parabolic distribution of Poiseuille flow. This process can be analyzed mathematically in one of two ways. Firs t , the in-tegration can be performed in the downstream direction so that the boundary layer growth is calculated for an accelerated external stream. Secondly, it is possible to analyze the progressive deviation of the profile from its asymptotic parabolic distri-bution, i.e., integration can proceed in the upstream direction. Having obtained both solutions, say in the form of series expansions, we can retain a sufficient number of terms i n either of them, joining the two solutions at a section where both are st i l l Chapter 4: Numerical Results 50 applicable. In this way the solution for the whole inlet is obtained. Figure 4.9: Flow in the inlet section of a Pipe As shown i n figure (4.9), this problem has a singularity at x = 0, y = R, namely u(Q-,R) = U 0 u(0 + , i r : ) = 0. We start by satisfying the continuity equation: fR(2Try)-u(x,y)dy = U 0 • ( r f ) J o or, fuydy = \VoR2. JO z Next we introduce a non-dimensional parameter ll_ _ x_ 6~\l R e ' R and expand U(x) = u (x ,0) in a series form and balance the momentum equation. Further details are given in the book by Schlichting [26]. He has reported the following values for the inlet length and the additional pressure drop Apadd-Chapter 4: Numerical Results 51 lE '= 0.04(2R)-Re &Padd = -0.3005i? e Using equations (3.14) and (3.15), we can then compute stream function values. A t the entrance, y dx and I ^ = „ = l • y dy Hence, y2 y2 Choosing ib = 0 along the centerline of the pipe gives C = 0 and yj'entrance = —• W e also know that at the exit, the velocity is fully developed. Hence v = 0 and u = A{1 — y2) where A is a constant. Aga in integrating the velocity fields and choosing tp = 0 along the centerline, yjexn = ~^y2(^ ~ y2)- Since the volume of flow must remain constant at any section, tbentrance(y = 1) = i>exit(y = 1)- This gives A = 2, i.e., the centerline velocity is doubled. 4.3.2 Finite Element Results As noted i n subsection 4.3.1, this problem has a singularity at x = 0, y = R. For this reason, the finite element grid used to solve this problem needs to be considerably finer near this point as compared to the grid downstream. A grid of 20 x 10 elements shown i n figure (4.10) is used wi th the length and width of successive elements chosen according to an arithmetic progression. Another patch of 4 x 10 elements precedes the pipe entrance representing the undisturbed uniformly distributed flow of magnitude U 0 = 1.0. The finite element program is run for Re = 0, 1, 10, 50, 100, 150, 200, 250, 300, 350 and 400 with a tolerance level of 1.0 x 10~ 5 . Results converge in 4 Chapter 4: Numerical Results 52 to 8 iterations wi th more iterations required as the Reynolds number is increased progressively. The distance along the centerline of the pipe between the entrance and a point ZQ is called the development length where ZQ is chosen such that u(oo,0) - u ( Z o , 0 ) < 0.001. u(oo,0) Twenty streamline contours equally spaced between the max imum and the min imum values of the stream function are plotted in figures (4.11), (4.12), (4.13) and (4.14). Aga in the y—dimension has been enlarged four times to clearly reveal the contours. Twenty streamline contours equally spaced between the max imum and the min imum values of the stream function are plotted i n figures (4.11), (4.12), (4.13) and (4.14). Aga in the y—dimension has been enlarged four times to clearly reveal the contours. The relationship between ZQ and RE is shown i n figure (4.15). Change in u velocity along the centerline as the flow develops is presented in figure (4.16) for different values of R.. 4.4 Flow Through an Axisymmetric Orifice A steady flow problem of interest to both engineers and mathematicians is that of a viscous, incompressible fluid through an orifice. Greenspan [10] used a finite difference method wi th a simple smoothing which yielded diagonally dominant systems of linear algebraic equations. Olson [18] used a pseudo-variational finite element theory with a high precision 18 degrees-of-freedom triangular element to get good results for Reynolds numbers up to Re = 20. Figure (4.17) shows a step change in diameter of the pipe over the range a < x < (a + 2/). For the purposes of comparison, we chose R = 1.0, / = 0.5 , 6 = 0.5. a = 5.25 and b = 15.55 such that at the exit, the flow was fully developed. The finite element grid used for a section bounded by the wal l and the centerline is shown in figure (4.18). The rigid boundary conditions were imposed as follows: Chapter 4: Numerical Results Finite Element Grid for Developing Flow R = 1.0 length of pipe=15.0 Number of nodes = 789 Number of elements = 280 Number of equations for steady flow = 1665 x—axis scale: y—axis scale: 1 cm = 1.104 units 1 cm = 0.276 units Figure 4.10: Finite element grid for developing flow Chapter 4: Numerical Results RE= 0 RE= 1 RE= 10 Figure 4.11: Streamline contours for developing flow for Re =0, 1 and 10 Chapter 4: Numerical Results RE- 50 Figure 4.12: Streamline contours for developing flow for Re =50, 100 and 150 Chapter 4: Numerical Results RE- 200 56 RE= 300 Figure 4.13: Streamline contours for developing flow for Re =200, 250 and 300 Chapter 4: Numerical Results RE= 350 Figure 4.14: Streamline contours for developing flow for Re =350 and 400 Chapter 4: Numerical Results 5 0 -i 0 - j — i — \ — i — | — i — i — i — i — i — i — i — | — i — i — i — | — i — i — i — | 0 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 Reynold number (Re) Figure 4.15: — versus Re Chapter 4: Numerical Results 2 . 0 o 1.6 H ^ > 1 .4 CD Co 1-2 1 .0 / / ' / / / I I I I / i — r 1 r R e =0 (l inear R e =10 R e - 5 0 R e =100 R e - 1 5 0 i 1 1 1 r ~> 1 ' 1 1 1 3 5 7 9 11 13 15 Distance from entrance Figure 4.16: Velocity distribution along the centerline for different Re Chapter 4: Numerical Results 60 Figure 4.17: Flow through an axisymmetric orifice • u = 4(1 — y2), v = 0 on the upstream edge (Poiseuille conditions), • u = 0, v = 0 all along the wall, • v = 0 along the centerline, and • v == 0, p = 0 on the downstream edge. du The natural boundary condition — = 0 was left for the program to approximate oy along the centerline. Since at the entrance v = 0, 4>{Q,y) = f\{\-y2)-y dy J o = y2(2-y2) and thus V(0,1) = 1.0. The program was run for Re = 0 (linear), 1, 5, 10, 15 and 20 and the stream function solution was obtained from the primitive variable solution. Streamline con-tour plots based on a 100 x 100 rectangular grid using the interpolation functions defined in Chapter 3, are shown in figures (4.19), (4.20) and (4.21). The first ten Chapter 4: Numerical Results 61 Finite Element Grid for Orifice Flow a = 4.75 b = 17.05 / = 0.5 R= 1.0 8 = 0.5 Number of nodes = 289 Number of elements = 80 Number of equations for steady flow = 479 x-axis scale: 1 cm = 1.540 units y-axis scale: 1 cm = 0.385 units Figure 4.18: Finite element grid for flow through an axisymmetric orifice Chapter 4: Numerical Results 62 R E = 0 R E = 1 \ 1 =| R E = 5 \ r~ = Figure 4.19: Streamline contours for orifice flow for Re =0, 1 and 5 Chapter 4: Numerical Results 63 R E = 1 0 — • — —Cv R E = 1 5 5fc R E = 2 0 Figure 4.20: Streamline contours for orifice flow for i2e =10, 15 and 20 Chapter 4: Numerical Results 64 (a) Greenspan [10], 1973 (c) present study Figure 4.21: Comparison of streamline contours for orifice flow for Re = 10 Chapter 4: Numerical Results 65 contours have an increment of 0.1, the value of the 11th contour is ip = 1.001 and the subsequent contours have an increment of 0.005. Aga in the y—dimension has been enlarged four times in the first two figures. In figure (4.21), the finite differ-ence results of Greenspan [10], and the pseudo-variational finite element results of Olson [18] are shown along wi th the results of the present study. It is seen that the contours obtained using the 8-noded isoparametric elements compare well with those of Greenspan's finite difference results and Olson's high-precision finite element results using the same contour values thus indicating that the present axisymmetric program is accurately representing the phenomenon. The small bubble upstream of the step observed by Greenspan is missing here, probably because of grid coarseness. 4.5 Flow Through a Stenosed Pipe In this section numerical results are presented for the parametric study performed for the case of flow through a stenosed pipe. A number of researchers have been working on this problem for the past few years and a brief summary of the past work is given i n Chapter 1. The geometry of the stenosis, various geometrical parameters, the nature of the constriction and the choice of finite element element grids (coarse and fine) are presented i n the first subsection. Results for the steady flow, and for periodic flow for Rw = 5, Rw = 10 and Rw = 20 are presented in the subsequent subsections. Shear stress along the wal l of the pipe is determined by the method outlined in Appendix D . 4.5.1 Geometry of the Stenosis Figure (4.22) represents a model of the stenosis. The various geometric parameters are as follows: a is the length of the unconstricted region before the stenosis; b is the length of the unconstricted region after the stenosis; 21 is the length of the stenosis; Chapter 4: Numerical Results 66 u = U 0 ( y , t ) v = Q 1 — -*-x Figure 4.22: Geometry of Stenosis i n an axisymmetric pipe r(x) is the effective radius of the pipe at any cross-section x; R is the radius at any unconstricted cross-section; and S is the maximum radial constriction so that the radius of the pipe at its m i n i m u m is (R — 6). Furthermore, the stenosis is symmetric about x = (a + I). B y suitably varying the a, b, I, 8, and choosing appropriate representation of r(x) we can represent almost any k ind of stenosis. Chakravarty [2,3], Forrester [8,9], Deshpande [6] and many others have chosen to model the constriction wi th a half cosine wave, i.e., 7r(x — (a + /)) r(x) = < R 1 + cos a<x< (a + 21); otherwise. Thi s model , however, has a drawback v iz . the beginning and end of the stenosis are dr not smooth: that is, ^— has a jump at x = a and x = (a + 21). To get around this difficulty, some researchers use a full cosine wave to model the stenosis. Another common technique is to assume that the variation of radius between x = a and x = (a + 2/) is circular, i.e., for some values for or, /3 and 7, ( x - a ) 2 + ( j / - / ? ) 2 = 7 2 . Chapter 4: Numerical Results 67 To determine a , /? and 7, we use y = R at x = a and x = (a + 2/), and y = (R — 8) at x = (a + /). This gives r(x) = < (*2-/2) (tS2 + /2) 1 2 - [x - (a + Of } a < x < ( a + 2Z); 1 2 2(5 otherwise. This model also suffers from the same drawback, although to a lesser extent. We have adopted this model for the present study. The choice of 8 is governed by the area reduction desired. Normal ly one restricts 8 to a m a x i m u m of 0.5JR. We chose 8 = 0.707i?; this results i n a 50% area reduction at the max imum constriction. Choosing the other parameters is not so easy; b, for example, should be long enough to allow the flow to develop again after separation. As a result, b is highly dependent on the Reynolds number. Likewise, a should be long enough to compensate for any separation occurring near the beginning of stenosis— especially at higher frequency Reynolds number, Rw. The half-length of the stenosis, /, affects the rate of change of diameter. Thus a very small / means that the radius decreases, and subsequently increases, at a very fast rate and the resulting flow might be similar to the flow over an orifice. The grids shown in figure (4.23) and figure (4.24) were used for this problem with the parameters as given i n the table (4.6) below. The number of elements across the cross-section remains constant; between a < x < (a + 2/), the width of each element T\ X) is reduced by a factor of ———. The number of elements along the y—direction and in R , the pre-stenosis zone are chosen the same for both the grids. Over the length of the stenosis and in its wake, the grid points are closely spaced. 4.5.2 Steady Flow Through Stenosis For the case of steady flow through a stenosed pipe, both grids were used. A t the entrance section, x = 0, the flow is specified as fully developed, i.e., u(0,y) = (1 — Chapter 4: Numerical Results 68 Coarse Grid for Pulsatile Flow (Grid 1) a = 5.0 b = 5.0 / = 3.0 R = 1.0 8 = 0.707 Number of nodes = 381 Number of elements = 112 Number of equations for steady flow = 772 Number of equations for unsteady flow = 2316 x—axis scale: 1 cm = 1.052 units y—axis scale: 1 cm = 0.263 units Figure 4.23: Grid 1 used for analysis of stenotic flow Chapter 4: Numerical Results 69 Fine Grid for Pulsatile Flow (Grid 2) Number of nodes = 719 Number of elements = 216 Number of equations for steady flow = 1487 x-axis scale: 1 cm = 1.564 units y-axis scale: 1 cm = 0.391 units Figure 4.24: Grid 2 used for analysis of stenotic flow Chapter 4: Numerical Results 70 Table 4.6: Parameters for the finite element grids for stenotic flow Parameter Grid 1 Grid 2 Designation short/coarse long/fine No. of nodes 381 719 No. of elements in x—direction 14 27 No. of elements in y—direction 8 . 8 Unconstricted radius R 1.0 1.0 Pre-stenosis length a 5.0 5.0 Post-stenosis length b 5.0 10.0 Half-stenosis length / 3.0 3.0 M a x i m u m constriction S 0.707 0.707 y 2) and v = 0. Results converged for Reynolds numbers up to Re = 500 for the coarse grid (grid 1). For Reynolds numbers higher than 500, grid 1 did not lead to converged solutions due to either, shorter post-stenosis length or due to inabil i ty of the finite elements to represent the rapid changes in primitive variables because of the coarseness of the grid; perhaps a combination of both the above. G r i d 2, however, yielded converged solutions up to very high Reynolds number and in relatively less number of iterations. We present these results for Reynolds numbers up to Re = 2000. We note that the pre- and post-stenosis lengths for grid 1 are the same and hence the linear problem, i.e., the flow through the stenosis at zero Reynolds number, is symmetric about the line x = (a +1). As such we expect the solution to be symmetric about this line; this indeed is the case. Furthermore, because of this symmetry, the y—component of fluid velocity along the grid line passing through x = (a + l) should be zero. Since the area of cross-section is effectively halved at this point, u must be doubled. This , too, is satisfied by the finite element results. A t higher Reynolds numbers, beginning at Re = 300, a very characteristic region of separation is observed in the wake of the stenosis. F rom the streamline contours for the coarse grid, we can see the onset of the disturbance which subsequently leads Chapter 4: Numerical Results 71 to the separation. Streamline contours for grid 2 very clearly reveal this phenomenon. A t Reynolds numbers higher than Re = 2000, the contours plots are a l i t t le zig-zagged which suggests that the results thus obtained are grid dependent. Thi s numerical in-accuracy results from the fact that in the current formulation the continuity equation is not satisfied exactly, whereas the streamline formulation assumes it to be satisfied exactly. This can be avoided by analyzing the flow wi th a very fine grid. Streamline contours for selected values of Re are given in figures (4.25-4.27) for the coarse grid and i n figures (4.28-4.33) for the finer grid. A 100 x 100 rectangular grid in conjunction wi th the scheme outlined in Appendix C is used to plot the contours. Twenty streamline contours equally spaced between the m i n i m u m and the max imum values are plotted for grid 1. For the finer grid, contours 1-10 have an increment of 0.02, contours 11-14 have an increment of 0.01, contours 15-16 have an increment of 0.005 and subsequent contours have an increment of 0.0005. M a x i m u m stream function values calculated at the rectangular grid points used for plott ing for different values of Re are summarized in table (4.7). Because of the no slip condition along the wal l , the fluid experiences a high a shear stress there. The shear acting on the wal l is equal and opposite of the shear acting on the fluid at y = r(x). Computat ion of the shear stress from the primitive variables solution is similar to the computation of vorticity £ and further details are given i n Appendix D. A t a fully developed, non-stenosed, cross-section the non-dimensional shear stress equals 2. Over the length of the stenosis, wi th the increase in velocity gradient, there is a corresponding increase in the shear stress. In the wake of the stenosis, as separation occurs, the velocity gradients, and consequently the shear stress changes sign. This is shown in figure (4.34) for the coarse grid and in figures (4.35) and (4.36) for the fine grid. Furthermore, up to Re = 300, the results obtained using the coarse grid are very close to those obtained using the fine grid. This suggests that grid 1 is sufficiently accurate up to that Reynolds number. Chapter 4: Numerical Results RE- 0 72 Figure 4.25: Streamline contours for steady flow through stenosis for grid 1 for Re 1 and 10 Chapter 4: Numerical Results 73 R E = 5 0 R E = 1 0 0 R E = 1 5 0 Figure 4.26: Streamline contours for steady flow through stenosis for grid 1 for Re = 50, 100 and 150 Chapter 4: Numerical Results RE- 2 0 0 74 Figure 4.27: Streamline contours for steady flow through stenosis for grid 1 for Re = 200, 250 and 300 Chapter 4: Numerical Results 75 R E = 0 R E = 1 0 R E = 5 0 Figure 4.28: Streamline contours for steady flow through stenosis for grid 2 for Rt = 0, 10 and 50 Chapter 4: Numerical Results 76 R E = 1 0 0 R E = 2 0 0 Figure 4.29: Streamline contours for steady flow through stenosis for grid 2 for Re = 100, 150 and 200 Chapter 4: Numerical Results 77 RE= 300 RE-= 5 0 0 Figure 4.30: Streamline contours for steady flow through stenosis for grid 2 for Re = 300, 400 and 500 Chapter 4: Numerical Results RE= 600 78 RE= 700 RE= 800 Figure 4.31: Streamline contours for steady flow through stenosis for grid 2 for Re = 600, 700 and 800 Chapter 4: Numerical Results 79 RE= 900 RE= 1250 Figure 4.32: Streamline contours for steady flow through stenosis for grid 2 for Re = 900, 1000 and 1250 Chapter 4: Numerical Results 80 RE = 1500 RE= 1750 Figure 4.33: Streamline contours for steady flow through stenosis for grid 2 for Re = 1500, 1750 and 2000 Chapter 4: Numerical Results 81 Grid 1 — 2~I i i i i i i i i i i i i i i i i 0 4 8 1 2 1 6 Distance from entrance Figure 4.34: Shear stress along the wall of the pipe for grid 1 for Re = 0 (linear), 50, .100, 200 and 300 Chapter 4: Numerical Results 82 - 2 -Grid 2 5 -/ \ / i </ V fc> 8 - '»' » 1' '.' A? / v: g>6- / \\! / \\ o \\ 4 -co i \ \ CO £ / M k 1 i * \ «i) ^ ? f / "' \ •So ^ '-' \.\ » \ \ i >'. V,\ CO 0 1 \ V, S o -V \ \ ^ -CO - l i n e a r - R e =50 - R e =100 - - R e =200 - R e - 3 0 0 T — i — r 0 T — i — i — | — r 12 — i — i — i — i — | 16 20 4 8 Distance from entrance Figure 4.35: Shear stress along the wall of the pipe for grid 2 for i?e = 0 (linear), 50, 100, 200 and 300 Chapter 4: Numerical Results 83 Grid 2 8 " CO - 4 -CO " —S~\—i—i—i—i—i—i—i—i—i—i—i—i—i i i — i — i — i — i — i 0 4 8 12 16 20 Distance from entrance Figure 4.36: Shear stress along the wall of the pipe for grid 2 for Re = 400, 500, 1000, 1500 and 2000 Chapter 4: Numerical Results 84 Table 4.7: M a x i m u m stream function values for steady flow through the stenosis G r i d 2 a = 5.0 6 = 10.0 / = 3.0 R = 1.0 8 = 0.707 Reynolds Number Re Maximum Stream Function Value Location X y 0 (linear) 0.2500 5.04 0.99 50 0.2500 5.04 0.99 100 • 0.2500 5.04 0.99 200 0.2500 12.18 0.99 300 0.2503 11.13 0.93 400 0.2508 11.13 0.90 500 0.2513 11.13 0.89 600 0.2517 11.13 0.88 700 0.2522 12.18 0.89 800 0.2527 12.18 0.88 900 0.2531 12.18 0.88 1000 0.2535 12.18 0.87 1250 0.2543 12.18 0.86 1500 0.2549 12.18 0.86 1750 0.2554 12.18 0.85 2000 0.2559 12.39 0.85 4.5.3 Pulsatile Flow Through Stenosis As pointed out i n section 3.4, the pulsatile nature of arterial fluid flow is approximated by the first three components of the Fourier series. Thus for any non-zero frequency Reynolds number R^, the problem size increases approximately nine-fold over that for the corresponding steady problem. For this reason, the finer grid (grid 2) could not be used for the results discussed in this sub-section; al l the solutions were obtained using grid 1 wi th a = 5.0, b = 5.0, / = 3.0, R = 1.0 and 8 = 0.707. A parametric study was conducted for frequency Reynolds numbers Ru = 5, 10 and 20, and steady Reynolds number Re = 0 (linear), 1, 5, 10, 20, 50, 100 and 150, thus resulting in a total of 24 Chapter 4: Numerical Results 85 cases. For Reynolds numbers higher than 150, grid 1 was not able to represent the flow satisfactorily and after a few iterations, the solution went divergent. Stream function values are computed from the primitive variable solutions and a 100 x 100 rectangular grid in conjunction wi th the scheme outline i n Appendix C is used to plot the contours. Twenty streakline contours equally spaced between the m i n i m u m and the max imum values are plotted at t = 0, —, ir and — . Computation L* L* of the shear stress from the primitive variables solution is similar to the computation of vorticity £; further datails are given in Appendix D . A long the curved portion of the wal l , axial stress <rx, radial stress o~y and the shear stress rxy are computed. To find the shear stress acting on the wall , r ^ , Mohr ' s circle construction is used. For any fixed Reynolds number, when the steady and the periodic components of the pulsatile flow are additive, the "true" Reynolds number is higher and consequently the 'deviation' from the fully developed flow is more. Stream function contours shown in figures (4.37-4.60) clearly reveal the effects of the stenosis. The instantaneous shear values along the wal l of the pipe are shown i n figures (4.61-4.63) for t = 0, —, -K and 3TT . . — . A t the beginning of the cycle, t = 0, the steady and the periodic components Lt are additive and hence the shear stress and the stream function values along the wal l are max imum at this moment. A t quarter cycle, t = —, only the steady and L\ the sinusoidal components are effective and hence the maximum shear stress and the m a x i m u m value of stream function along the wal l drop; at half cycle, t = TT, the cosinusoidal component is substractive and hence the maximum shear stress and the m a x i m u m value of stream function along the wal l is min imum at this instance. During the second half of the cycle, this process is reversed and the max imum value is again obtained at t = 2ir. Due to the coarseness of the grid, however, the instantaneous shear stress values along the wal l at various moments during a complete cycle showed a very small amount of oscillation near the min imum value. These oscillations, too, damp out at half cycle t = ir, when the cosinusoidal component is substractive. Chapter 4: Numerical Results 86 RE= 0 RW= 5 t = o RE=0 RW=5 RE= 0 RW= 5 t = 7T RE=0 RW= 5 Figure 4.37: Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 0 (linear) , . 7T 37T for t = 0, —, 7r and — Chapter 4: Numerical Results 87 RE = 1 RW= 5 t = o RE= 1 RW= 5 RE= 1 RW= 5 Figure 4.38: Streaklines for pulsatile flow through stenosis for Ru, = 5, Re = 1 for „ 7T , 37T t = 0, - , 7T and y Chapter 4: Numerical Results 88 RE= 5 RW= 5 t = o RE= 5 RW= 5 7T t = 2 RE= 5 RW= 5 t = 7T RE= 5 RW= 5 Figure 4.39: Streaklines for pulsatile flow through stenosis for R\j = 5, Re = 5 for „ TC , 37T f = 0, -, 7T and — Chapter 4: Numerical Results 89 RE = 10 RW= 5 t = 0 RE= 10 RW = 5 RE= 10 RW= 5 t = X RE= 10 RW= 5 Figure 4.40: Streaklines for pulsatile flow through stenosis for R^ t = 0, - , x and — = 5, Re = 10 for Chapter 4: Numerical Results 90 RE= 20 RW= 5 •RE=20 RW= 5 Figure 4.41: Streaklines for pulsatile flow through stenosis for R^, « = 0, - , 7r and — = 5, Re = 20 for Chapter 4: Numerical Results 91 RE= 50 RW= 5 t = o RE= 50 RW= 5 RE= 50 RW= 5 t = 7T RE= 50 RW= 5 Figure 4.42: Streaklines for pulsatile flow through stenosis for Ru, ' 2 ' T T = 5, Re = 50 for Chapter 4: Numerical Results RE- 100 RW= 5 t = o RE= 100 RW= 5 RE= 100 RW= 5 t = TT RE- 100 RW = 5 Figure 4.43: Streaklines for pulsatile flow through stenosis for Ru, 7T - 37T t = 0, - , and y = 5, Re = 100 for Chapter 4: Numerical Results 93 RE= 150 RW= 5 t = 0 RE= 150 RW= 5 7T <=2 RE= 150 RW= 5 t = 7T RE= 150 RW= 5 Figure 4.44: Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 150 for „ 7r 3TT t = 0, - , TT and y Chapter 4: Numerical Results 94 RE= 0 RW= 10 t = o RE= 0 RW= 10 RE = 0 RW= 10 t = 7T RE= 0 RW= 10 Figure 4.45: Streaklines for pulsatile flow through stenosis for R\, = 10, Re = 0 (linear) for t = 0, ^ , TT and — Chapter 4: Numerical Results 95 RE= 1 RW= 10 Figure 4.46: Streaklines for pulsatile flow through stenosis for R^, = < = 0, - , 7T and — 10, Re = 1 for Chapter 4: Numerical Results 96 RE= 5 RW= 10 • i = 0 RE- 5 RW= 10 RE= 5 RW= 10 ' Figure 4.47: Streaklines for pulsatile flow through stenosis for R\, = 10, Re = 5 for ' 2' * T Chapter 4: Numerical Results 97 RE= 10 RW= 10 t = o RE= 10 RW= 10 RE= 10 RW= 10 Figure 4.48: Streaklines for pulsatile flow through stenosis for = t = 0, - , 7T and — 10, Re = 10 for Chapter 4: Numerical Results 98 RE= 20 RW= 10 t = o RE= 20 RW= 10 RE= 20 RW= 10 RE= 20 RW= 10 Figure 4.49: Streaklines for pulsatile flow through stenosis for R\j = 10, Re = 20 for t = 0, —, 7r and — 2 2 Chapter 4: Numerical Results 99 RE = 50 . RW= 10 t = o RE= 50 RW= 10 RE= 50 RW= 10 Chapter 4: Numerical Results 100 RE= 100 RW= 10 t = o RE= 100 RW= 10 RE = 100 RW= 10 t = 7T RE= 100 RW= 10 Figure 4.51: Streaklines for pulsatile flow through stenosis for R^ — 10, Re = 100 for „ 7T , 37T t = 0, - , 7T and y Chapter 4: Numerical Results 101 RE = 150 RW= 10 t = o RE= 150 RW= 10 RE= 150 RW= 10 RE= 150 RW= 10 Figure 4.52: Streaklines for pulsatile flow through stenosis for R^, = 10, Re = 150 for Chapter 4: Numerical Results 102 RE= 0 RW= 20 t = o REr 0 RW= 20 RE= 0 RW= 20 t = IT RE= 0 RW= 20 Figure 4.53: Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 0 (linear) for t = 0, ^ , ir and — Chapter 4: Numerical Results RE- 1 RW= 20 103 Chapter 4: Numerical Results 104 RE= 5 RW= 20 RE= 5 RW= 20 RE= 5 RW= 20 t = V RE- 5 RW= 20 Figure 4.55: Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 5 for Chapter 4: Numerical Results 105 RE = 10 RW= 20 t = o RE= 10 RWr 20 Figure 4.56: Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 10 for t = 0, - , 7T and — Chapter 4: Numerical Results RE- 20 RW= 20 106 Figure 4.57: Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 20 for t = 0, - , JT and — Chapter 4: Numerical Results RE= 50 RW= 20 107 t = 0 RE=50 RW=20 RE= 50 RW= 20 t = TX RE= 50 RW= 20 Chapter 4: Numerical Results 108 RE= 1 0 0 RW= 2 0 t = o RE= 1 0 0 RW= 2 0 RE- 1 0 0 RW= 2 0 t = IT RE= 1 0 0 RW- 2 0 Figure 4.59: Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 100 for t = 0, —, ir and — Chapter 4: Numerical Results RE= 150 RW= 20 109 t = 0 RE= 150 RW= 20 RE= 150 RW= 20 t = 7T RE= 150 RW= 20 Figure 4.60: Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 150 for t = 0, —, 7r and — Chapter 4: Numerical Results 110 < = 0 5 -- 8 - | — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i — i 0 4 . 8 12 16 Distance from entrance -8 I i i i , , i 0 * 8 12 16 Distance from entrance 0) 40' 5 24-Co CO ca 5-CO S» a CD CO 16-8-0--8' t = 7T — linear — R.= 10 — R.-50 — R.-100 - - R.-150 - 1 — 12 16 0 4 ' ' ' 8 Distance from entrance 8 12 16 Distance from entrance Figure 4.61: Instantaneous shear stress along the wall of the pipe for grid 1 for R\j = 5 and Re = 0 (linear), 10, 50, 100 and 150 Chapter 4: Numerical Results 111 t = 0 5 •• -8-|—i—i—i—r—i—'—i—i—1—1—i—i—'—i—<—i 0 4 6 12 16 Distance from entrance 7T to —8 [ I I I — i i i i i—i—i i i i—i—i—i 0 4- 8 12 16 Distance from entrance lineor R.= 10 - R.=50 R.= 100 R.= 150 Distance from entrance t = >7T lineor R.-10 R.=50 R.= 100 6 * 8 12 16 Distance from entrance Figure 4.62: Instantaneous shear stress along the wall of the pipe for grid 1 for Ru, = 10 and Re = 0 (linear), 10, 50, 100 and 150 Chapter 4: Numerical Results 112 •8-T—' r—T 1 1 r—I 1 1 1—I 1 r—T 1—| —8"| 1 1 I I 1 1 1—1 1 1 1 1 1 1 I | 0 4 8 12 16 0 4 8 12 16 Distance from entrance Distance from entrance Figure 4.63: Instantaneous shear stress along the wall of the pipe for grid 1 for R^, = 20 and Re = 0 (linear), 10, 50, 100 and 150 C H A P T E R 5 Conclusions 5.1 Concluding Remarks The present method of representing the complete solution as a sum of steady, cosi-nusoidal and sinusoidal components seems to work very well . The modified method of averaging used is also observed to be quite accurate. In particular, the circular arc approximation of stenoses works well. The overall agreement between the present study and other published results in the literature is good. The modified Newton-Raphson procedure for solving the nonlinear algebraic equations is very successful and converges fast for the range of frequency Reynolds numbers and Reynolds numbers considered i n this study. The results obtained for the fully developed flow under the influence of a cos-inusoidal pressure difference agree well wi th the numerical and graphical values re-ported by Schlicting. The streakline plots are almost completely horizontal indicating that even though the solution becomes increasingly wavy wi th increasing frequency Reynolds number, the finite element approximation is very good. For the devel-oping flow, the program consistently exhibited good stability, fast convergence and good accuracy, even for higher Reynolds numbers. The effect of Reynolds number on the development length was in good agreement wi th the theoretical values reported in [26]. The additional pressure drop, Apadd obtained from the present study showed some discrepancies probably due to the inadequacies of the bilinear pressure shape functions and coarseness of the grid. 113 Chapter 5: Conclusions 114 The pr imit ive variable solution as well as the stream function contour plots for axisymmetric flow through the orifice converged fast for Reynolds number up to Re = 30 even though at the min imum section area of cross section of the pipe was reduced by 75% and the u—velocity was quadrupled. The stream function contours for Re = 10, plotted to 1:1 scale, compares well wi th the pseudo-variational finite element results of Olson [18] and the finite difference results of Greenspan [10]. It may be noted that Greenspan's method converged in about 40 to 90 'outer' iterations compared to 4 to 8 iterations required by the present method. Greenspan's method also required a very fine finite difference mesh. Steady results for flow through a stenosed pipe were obtained for two different grids and were used to check for the consistency of the results. A l l the major characteristics of the flow such as separation and re-attachment points, variation of shear stress along the wal l , including the max imum and the min imum values, max imum stream function values, etc., show a monotonic improvement wi th the finer grid. The separation and the re-attachment points as obtained from the stream function results and from the shear stress results are i n total agreement wi th each other. For unsteady flow through a stenosed pipe only the coarse grid was used and satisfactory results were obtained for R^ = 5, 10 and 20, and for Reynolds numbers up to Re = 150. For any fixed Reynolds number, when the steady and the periodic components of the pulsatile flow are additive, the "true" Reynolds number is higher and consequently the 'deviation' from the fully developed flow is more. Stream func-tion contours clearly reveal the effects of the stenosis and instantaneous shear values along the wal l of the pipe indicate the init iat ion of the separation phenomenon; for higher Reynolds numbers, the flow w i l l separate and the region of back flow w i l l have negative shear stress. Due to the coarseness of the grid, however, the instantaneous shear stress values along the wal l at various moments during a complete cycle showed a very small amount of oscillation near the min imum value. Chapter 5: Conclusions 115 5.2 Suggestions for Further Developments Some specific recommendations for further studies based on the work in this thesis are: • perform the analysis wi th the transformations V = f(x,y) to map the curved boundary onto a straight line £ = £o, and n = tanh(i^x) to map the infinite domain to the finite domain. • perform the analysis for multiple stenoses in a series. • conduct a grid refinement study to determine its effects on the numerical results. • perform the analysis w i t h high precision elements using the stream function and the vorticity as the primitive variables and compare and contrast it wi th the present study. Bibliography [1] Ca i lk , D . A . , and Naghdi , P . M . , Axisymmetric Motion of a Viscous Fluid In-side a Slender Surface of Revolution, Transactions of A S M E : Journal of Appl ied Mechanics, V o l . 54, 1987, pp. 190-196. [2] Chakravarty, S., Effects of Stenosis on the Flow Behavior of Blood in Artery, International Journal of Engineering Science, V o l . 25, 1987, pp. 1003-1016. [3] Chakravarty, S., and Mis ra , J .C . , Flow in Arteries in the Presence of Stenosis, Journal of Biomechanics, V o l . 19, 1986, pp. 907-918. [4] Daly, Bart J . , A Numerical Study of Pulsatile Flow Through Stenosed Canine Femoral Arteries, Journal of Biomechanics, V o l . 9, 1976, pp. 465-475. [5] De Bremaecker, J.C1., Penalty Solution of the Navier-Stokes Equations, Com-puters & Fluids , V o l . 15, 1987, pp. 275-280. [6] Deshpande, M . D . , Giddens, D .P . , and M a b o n , R . F . , Steady Laminar Flow Through Modeled Vascular Stenosis, Journal of Biomechanics, V o l . 9, 1976, pp. 165-174. [7] Deshpande, M . D . , Steady Laminar and Turbulent Flow Through Vascular Steno-sis Models, P h . D . Thesis, Georgia Institute of Technology, At lanta , G A 30332, U S A , 1977. [8] Forrester, J . H . , and Young, D . F . , Flow Through a Converging-Diverging Tube and its Implications in Occlusive Vascular Disease—I: Theoretical Development, Journal of Biomechanics, V o l . 3, 1970, pp. 297-305. 116 Bibliography 117 [9] Forrester, J . H . , and Young, D . F . , Flow Through a Converging-Diverging Tube and its Implications in Occlusive Vascular Disease—77; Theoretical and Exper-imental Results and Their Implications, Journal of Biomechanics, V o l . 3, 1970, pp. 307-316. [10] Greenspan, D . , Numerical Studies of Viscous, Incompressible Flow Through an Orifice for Arbitrary Reynolds Number, International Journal for Numerical Methods i n Engineering, V o l . 6, 1973, pp. 489-496. [11] Krause, E . , Computational Fluid Dynamics: Its Present Status and Future Di-rection, Computers & Fluids , V o l . 13, 1985, pp. 239-269. [12] Lohner, R. , Finite Elements in CFD: What Lies Ahead, International Journal for Numerical Methods i n Engineering, V o l . 24, 1987, pp. 1741-1756. [13] MacDona ld , D . A . , Pulsatile Flow in Catheterized Artery, Journal of Biomechan-ics, V o l . 19, 1986, pp. 239-249. [14] McLacha lan , N . W . , Ordinary Non-linear Differential Equations in Engineering and Physical Sciences, Oxford University Press, 1950. [15] Moorty , S., A Parametric Study of Rigid Body-Viscous Flow Interaction, M . A . S c . Thesis, The University of Br i t i sh Columbia , Vancouver, Canada, 1987. [16] Moorty , S., A Numerical Study of Low Reynolds Number Fluid-Structure Inter-action, Journal of Fluids and Structures, V o l . 3, 1989, pp. 37-60. [17] Nayfeh, A . H . , and Mook , D . T . , Nonlinear Oscillations, Wi ley International, New York , 1979. [18] Olson, M . D . , Variat ional-Finite Element Methods for Two-Dimensional and A x -isymmetric Navier-Stokes Equations, Finite Elements in Fluids, V o l . 1, pp. 57-72, editors Gallagher, R . H . , et al . , John Wiley, London, 1975. Bibliography 118 [19] Pat tani , P . G . , Nonlinear Analysis of Rigid Body-Viscous Flow Interaction, P h . D . Thesis, The University of Br i t i sh Columbia , Vancouver, Canada, 1986. [20] Pat tani , P . G . , and Olson, M . D . , Periodic Solution of Rigid Body-Viscous Flow Interaction, International Journal for Numerical Methods i n Fluids , V o l . 7, 1987, pp. 653-695. [21] Phi lpot , E . , Yoganathan, A . P . , Sung, H . - W . , Woo, Y . - R . , Franch, R . H . , Sahn, D . J . , and Valdez-Cruz, L . , In-Vitro Pulsatile Flow Visualization Studies in a Pulmonary Artery Model, Transactions of A S M E : Journal of Biomechanical E n -gineering, V o l . 107, 1985, pp. 368-375. [22] Porenta, G . , Young D . F . , and Rogge, T . R . , A Finite-Element Model of Blood Flow in Arteries Including Taper, Branches and Obstructions, Transactions of A S M E : Journal of Biomechanical Engineering, V o l . 108, 1986, pp. 161-167. [23] Pozrikidis , C , A Study of Peristaltic Flow, Journal of F l u i d Mechanics, V o l . 180, 1987, pp. 515-527. [24] Rice , J . G . , and Schnipke, R . J . , An Equal Order Velocity Pressure Formulation That Does Not Exhibit Spurious Pressure Modes, Computer Methods in Appl ied Mechanics and Engineering, V o l . 58, 1986, pp. 135-149. [25] Robichaud, M . P . , and Tanguy, P . A . , Finite Element Solution of Three-Dim-ensional Incompressible Fluid Problems by a Preconditioned Conjugate Residual Method, International Journal for Numerical Methods in Engineering, Vo l . 24, 1987, pp. 447-457. [26] Schlichting, H . , Boundary Layer Theory, M c G r a w H i l l , 6th edition, 1968. [27] Scott, P.S., and M i r z a , F . A . , A Finite Element Analysis of Laminar Flows Through Planar and Axisymmetric Abrupt Expansions, Computers & Fluids , Bibliography 119 Vol. 14, 1986, pp. 423-432. [28] Srivastava, L.M., Flow of Couple Stress Fluid Through Stenotic Blood Vessels, Journal of Biomechanics, Vol. 18, 1985, pp. 479-485. [29] Talukder, N., Karayannacos, P.E., Nerem, R.M., and Vasco, J.S., An Experimen-tal Study of the Fluid Dynamics of Multiple Noncritical Stenoses, Transactions of the ASME: Journal of Biomechanical Engineering, Vol. 99, 1977, pp. 74-82. [30] Temam, R., Navier-Stokes Equations, North-Holland Publishing Company, 1st revised edition, 1979. [31] Tuann, S.-Y., and Olson, M.D., Numerical Studies of the Flow around a Circular Cylinder by a Finite Element Method, Computers in Fluids, Vol. 6, 1978, pp. 219-240. [32] Tuann, S.-Y., and Olson, M.D., A Study of Various Finite Element Solution Methods for the Navier-Stokes Equations, Structural Research Series, Report No. 14, Department of Civil Engineering, The University of British Columbia, Vancouver, B.C., CANADA, May 1976. [33] Utku, M., and Carey G.F., Stream Function Solution of Navier-Stokes Problems with Inter-Element Penalties, International Journal for Numerical Methods in Fluids, Vol. 7, 1987, pp. 191-193. [34] Yoganathan, A.P., Ball, J., Woo, Y.-R., Philpot, E.F., Sung, H.-W., Franch, R.H., and Sahn, D.J., Steady Flow Velocity Measurements in a Pulmonary Artery Model with Varying Degrees of Pulmonic Stenosis, Journal of Biomechanics, Vol. 19, 1986, pp. 129-146. [35] Young, D.F., Fluid Mechanics of Arterial Stenoses, Journal of Biomechanical Engineering, Vol. 101, 1979, pp. 157-175. Bibliography 120 [36] Zienkiewicz, O.C., The Finite Element Method in Engineering Science, McGraw Hill, London, 1971. [37] Zinser, W., and Benim, A.C., A Segregated Formulation of Navier-Stokes Equa-tions with Finite Elements, Computer Methods in Applied Mechanics and Engi-neering, Vol. 57, 1986, pp. 223-237. A P P E N D I X A Newton-Raphson Procedure The modified Newton-Raphson iteration procedure adopted in this study is outlined here. The scheme is outlined for a set of n simultaneous nonlinear algebraic equations i n n unknowns xx, x 2 , . . . , xn. These equations can be written in the form / i ( x i , x 2 , . . . , x n ) = 0 f2(xx,x2, . . . , x n ) = 0 ( A . l ) fn[xX,X2, . . . ,Xn) = 0. Expanding the i t h equation around Xio, x 2 0 , • - • , xnQ using the Taylor series, we obtain fi(xx,x2,...,xn) « fi{xx0,x20,... ,xn0) + A x i where the derivatives are evaluated at ( x i 0 , x 2 o, • • • , xn0) and A x i = Xi — Xio A x 2 = x 2 — X 2 Q (A.2) (A.3) A x n — Xji XJIQ. Similarly, carrying out the expansion for a l l the equations in ( A . l ) and rearranging, we obtain 121 Appendix A: Newton-Raphson Procedure 122 rfdfA fdfA _ fdfA \dxj0 \dx2)Q \dxn)Q fdfA (dfA (dfA VSa J 0 \dx2)0 \ d x n ) 0 fdfA (dfA _ (dfA \dxj0 \dx2)Q "* \dxn)Q Equations (A.4) can be written in the form [T]{Ax} = {—/} where the matrices denote the corresponding terms in equation (A.4) . The algorithm for obtaining the solution vector {x} is outlined below: 1. Guess the solution vector {x}0. 2. Compute [T]. 3. Compute {—/}. 4. Obta in {Ax} from [T]{Ax} = {-/}. 5. Obta in {x}x = {x}0 + {Ax}. 6. If {x}i is close to {x}0 w i th in the specified tolerance go to 8, otherwise set {x}0 = {x}x. 7. If [T] needs to be modified in this iteration go to 2, otherwise go to 3. 8. Set {x} = {x}i and stop. ' A x i ' ' —fi(xio,x2o,.. A x 2 -f2(xw,X20, • • • , Xno) < > = = < > (A.4) . A x n , • —fn(x10, #20) • • • > ^ n o ) > A P P E N D I X B Details of Tangent Stiffness Matrix and Nonlinear Load Vector The left hand side of equation (3.12) can be re-written as K 0 0 " 'A B C " [T] = 0 K M + V A 0 . 0 - M K . Z 0 A. where K, and M are the stiffness and mass matrices. Sub-matrices A to £, using index notation, are given as Au 4- fix Au 4- FH Av f)y Au 0' A = Re SfmjAm + &\mjAvm + 6 i j m A n 0 B = Re 2 ' £x Do i fix D U I cy r>v ijm m ' imj m. ' imi m 8X- Bv Bu 8X Bu + 8y Bv +6?- Bv 0 Re 2 • sx rvu. I cx fiu i sy /~iv 8ijm C771 &imjCm + 8fmjCm + 8\jmCvm 0 123 Appendix B: Details of Tangent Stiffness Matrix and Nonlinear Load Vector V = RK Cx DU i SX JDU I CI) Dt) cx TDV ijm m Uijm-Dm ^imj B m + fiimjBm + 0~ViimBVm 0 £ = RE The incremental solution vector is {Ax} = < r A A A B A C and the load vector is {-f} = K 0 0 0 0 K M - M K (A) B C \ 6?jk{A]At + B?B£/2 + C?C2/2) + 5? i fc(A?Aj£ + BVB£/2 + CvC^/2) ) 8fjk(A]Al + BJBJI/2 + C?C£/2) + ^ ( A J A J + B?B£/2 + C*Cj;/2) 0 ^•*(AJBJ + AJfBJ) + S*ik(A?Bl + AJJBy) -Re ( Sfjk(A]Bt + AlBJ) + Sfjk(A]Bl + AJBJ) 0 6fjk(A]Ct + AtC]) + 8yjk{A)Cl + A j q ) ^• fc(AJCj; + AZC?) + ^ f c (AVCJ + Aj£CJ) 0 A P P E N D I X C Determination of the Stream Function Values Consider an eight-noded isoparametric element as shown i n the figure (3.1). The shape functions iV,- and M,- for the element are presented in Chapter 3. The value of the stream function $ and the coordinates (x, y) at a point can be obtained in terms of the nodal values \1/,- and nodal coordinates (x,-,y,-). * = x = NiXi * = 1, 2,..., 8 V = Nan The following procedure is adopted to obtain the value of accurately at each node of the plott ing grid which lies wi th in the finite element. 1. Assume a quadratic interpolation for £ and n of the form £ = Ai + A2x + A3y + A4x2 + A5xy + A6y2 + ATx 2y + A8y 2x ( C . l ) 7] = B1 + B2x + B3y + B4x2 + B5xy + B6y2 + B7x2y + B8y2x The (x, y) and (£, n) values at each of the eight nodes on the element are known and hence eight equations each in A 1 , A 2 , . . . , A 8 and Bx, B2,..., B% are ob-tained and solved for. Thus for any point (x0,y0) i n the finite element, using equation ( C . l ) , corresponding ( ^ o ? 7 ? © ) coordinates can be calculated. 2. Determine the maximum and minimum value of x and y coordinates of the finite element. The approximate area occupied by this element is a rectangle of 125 Appendix C: Determination of the Stream Function Values 126 length xmax - xmin along the x - a x i s and height ymax - ymin along y - a x i s ; this approximation works reasonably well for an element which is not too skewed. 3. Determine (£, 7?) value of each point of the plotting grid wi th coordinates (x p , yp) w i t h i n the finite element as outlined below: (a) Take the in i t i a l guess value of (£P,77P) to be that obtained from equa-t ion ( C . l ) , i.e., iv = M + A2xp + A3yp + A j x 2 , + A5xpyp + A6y2p + A7x\yp + Asy2pxp np = B1 + B2xp + B3yp + B4xl + B5xpyp + B&y2v + B7x2pyp + B8y2xp (b) Calculate the shape functions Ni for i = 1 , 2 , . . . , 8 for these values of £ p , np using the equations given i n Chapter 3. (c) Calculate the Jacobian J. (d) Determine the (a/, y') coordinates of the point using the calculated value of shape functions. dNj dNj J = ds dNj ds .7 = 1,2,..., 8 (e) Determine the error in the (x, y) coordinates. Ax = xp — x' Ay = yp - y' Determine the error in the (£,77) coordinates. Appendix C: Determination of the Stream Function Vaiues 127 (g) The new value of (£, n) is given as ? = + ^ V' = VP + A??. W i t h these new values of (£, 77), repeat steps (b) to (f) unt i l desired convergence is achieved. 4. Once convergence has occurred in step 3, the stream function value can be obtained using \I> = N^i. A P P E N D I X D Determination of Stress Components For axisymmetric flow, the three non-trivial components of the stress tensor in their non-dimensional form are given as ndu 0-y = -P + 2 'xy dv dy _ (dv du s \dx dy. (D.l) Consider an eight-noded isoparametric element as shown in the figure (3.1). Using the shape functions AT,- and M t - presented in Chapter 3, equations ( D . l ) are discretized as -PjMj + 2Ui -PjMj + 2vi dNj dx dNi dNi + Ui dy dNi' * = 1,2,...,8; ; = 1,2,...,4. (D.2) dx ' dy Stress components crx, ery and rxy at some point (x, y) G Q can be obtained in terms of the nodal values of the primitive variables (u,v,p). The following procedure is adopted to obtain the values of o~x, o~y and rxy accurately at any given point (x0,y0). 1. Determine the element wi th in which point (x0, yQ) lies x . 1 If the point lies on a boundary, stresses are obtained for each element that borders at the given point and a weighted average is obtained; in case of a stress discontinuity all the values are reported. 128 Appendix D: Determination of Stress Components 129 2. Assume a quadratic interpolation for £ and rj of the form £ = At + A2x + Azy + A4x 2 + A5xy + A6y 2 + A7x 2y + A8y 2x (D.3) rj = B1 + B2x + B3y + B4x2 + Bhxy + B6y2 + B7x 2y + B&y 2x The (x, y) and (£, r?) values at each of the eight nodes on the element are known and hence eight equations each in A\, A2,..., As and B\, B2,..., B$ are obtained and solved for. 3. Determine (£, n) value of each point of the plott ing grid wi th coordinates (x 0 , y0) within the finite element as outlined below: (a) Take the in i t i a l guess value of (£ 0 , n0) to be that obtained from equa-tion (D.3), i.e., to = Ax + A2x0 + A3y0 + A4x 20 + A5x0y0 + A6y 20 + A7x 20y0 + A8y 2x0 r]0 = Bi + B2x0 + B3y0 + B4x\ + B5x0y0 + B6yl + B7x2Qy0 + B 8 y 2 x 0 dNi d N (b) Calculate the shape function derivatives -—- and -r-^- for i = 1 ,2 , . . . ,8 o£ or) for these values of £ 0 , 770 using the equations given in Chapter 3. (c) Calculate the Jacobian J. rdNi dN n ds  X% dtXi dN dN 1 ds  Vi dt  V i l (d) Determine the (x', y') coordinates of the point using the calculated value of shape functions. x' = NiXi i = l ,2, . . . ,8 y' = Niyi (e) Determine the error in the (x,y) coordinates. Ax = x0 — x' Ay = y0- y' Appendix D: Determination of Stress Components 130 (f) Determine the error in the (£,?/) coordinates. J An } ^  J Ay J (g) The new value of (£, n) is given as £' = & + A£ n' = n0 + An. W i t h these new values of (£, 77), repeat steps (b) to (f) u n t i l desired convergence is achieved. dNi dNi . 4. Once convergence has occurred m step 3, compute —— and —— using ox oy dN dNi N dx dN K dy ' dt dNi dn These derivatives can now be substituted into equations (D.2) to give stress results. Figure D . l : Use of transformation (x = x(x',y');y = y(x',y')) to determine stresses along a curved wal l Appendix D: Determination of Stress Components 131 To obtain stresses in any direction other than the global axes, we use the Mohr's circle construction. As shown in the figure ( D . l ) , for equil ibrium crx, o~y and rxy acting on unit surface area must balance crxi and Txiyi acting on y/2 surface area. Thus [crx + o~y) (crx - (jy) 2 2 I yvx — Vy . crxi = 1 — — cos 20 + rxy s i n 29 ((Tx + 0-y) (CTX — (Jy) <?y< = — — - — — cos 29 + rxy s i n 29 T x , y l = _ ( ^ - ( 7 i / ) s i n 2 g + T ^ c o s 2 9 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062609/manifest

Comment

Related Items