FINITE E L E M E N T ANALYSIS OF PERIODIC VISCOUS F L O W IN A C O N S T R I C T E D P I P E by Rajesh K u m a r Singh B . Tech., Indian Institute of Technology, K a n p u r , 1987 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E D E G R E E OF MASTER OF APPLIED SCIENCE in T H E FACULTY OF G R A D U A T E STUDIES DEPARTMENT OF CIVIL ENGINEERING We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA June, 1989 © Rajesh K u m a r Singh, 1989 In presenting this thesis i n p a r t i a l fulfillment of the requirements for an advanced degree at the University of B r i t i s h C o l u m b i a , I agree that the L i b r a r y shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes m a y be granted b y the head of my department or b y 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 B r i t i s h C o l u m b i a Vancouver, B . C . , C A N A D A V 6 T 1 W 5 June, 1989 Abstract T h e finite element method is extended to analyze axisymmetric periodic viscous flow i n a constricted pipe. T h i s 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 w a l l of the artery. T h e artery geometry is approximated as an axisymmetric channel. B l o o d is assumed to be a Newtonian fluid, i.e., it follows a linear relationship between the rate of shear strain and the shear stress. T h e numerical model uses the finite element method based on velocity-pressure p r i m i t i v e variable representation of the Navier-Stokes equations. Eight-noded iso- parametric elements w i t h quadratic interpolation for velocities and bilinear for pressure 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. T w o non-dimensional parameters are used to characterize the flow: the frequency Reynolds number i ^ , and the steady Reynolds number Re. T h e pulsatile flow is modeled as a linear combination of steady, cosinusoidal and sinusoidal components. T h e mag- nitude of each component is determined by m i n i m i z i n g the error i n approximation through Galerkin's procedure. Results are obtained for various values of and RE for laminar flow. Results are also presented for l i m i t i n g cases and, wherever possible, numerical results thus obtained are compared w i t h analytical and experimental results published i n the literature. ii Contents Abstract ii Tables vi Figures vii Notation xii Acknowledgement xvi 1 2 3 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 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 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 4 3.4 Characterization of A r t e r i a l F l u i d F l o w 25 3.5 Streamlines and Streaklines 30 Numerical Results 33 4.1 Steady F l o w T h r o u g h a P i p e 34 4.1.1 Theoretical Results 34 4.1.2 F i n i t e Element Results 35 4.2 4.3 5 Oscillating F l o w T h r o u g h a P i p e 36 4.2.1 Theoretical Results 36 4.2.2 F i n i t e Element Results 41 Developing F l o w i n the Inlet Length of a P i p e 49 4.3.1 Theoretical Results 49 4.3.2 F i n i t e Element Results 51 4.4 F l o w T h r o u g h an A x i s y m m e t r i c Orifice 52 4.5 F l o w T h r o u g h a Stenosed P i p e 65 4.5.1 Geometry of the Stenosis 65 4.5.2 Steady Flow Through Stenosis 67 4.5.3 Pulsatile F l o w Through Stenosis 84 Conclusions 113 5.1 C o n c l u d i n g 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 4.2 Comparison of exact and finite element solution for oscillatingflowfor 42 i% = 1.0 4.3 Comparison of exact and finite element solution for oscillating flow for 'Re = 10.0 4.4 42 Comparison of exact and finite element solution for oscillating flow for Re = 20.0 4.5 37 43 Instantaneous rate of flow through a pipe under cosinusoidal pressure difference 45 4.6 Parameters for the finite element grids for stenotic 4.7 Maximum stream function values for steady flow through the stenosis vi flow 70 84 Figures 2.1 Fluid Domain 3.1 Isoparametric Element 9 (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 parallel 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 4.15 ^ versus Re . . 57 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 Li L» (linear) for t = 0, —, it and — . 86 4.38 Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 1 for . 7T 37T t = 0, —, 7r and — Li 87 Li 4.39 Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 5 for „ 7T * = 0, - , Lt , 7T and 37T — 88 Lt 4.40 Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 10 for 3TT 7T < = 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 4.42 Streaklines for pulsatile flow through stenosis for R^, = 5, Re = 50 for „ 7T t = 0, - , 90 , 37T 7T and — 9 1 4.43 Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 100 for 7T 37T 7T 37T f = 0, —, - , 7r and and — 2' 2 4.44 Streaklines for pulsatile flow through stenosis for Rw = 5, Re = 150 for ^ f = 0, - , 7T and — 92 .93 4.45 Streaklines for pulsatile flow through stenosis for R^ = 10, Re = 0 (linear) for t = 0, —, 7r and — Li 94 Li 4.46 Streaklines for pulsatile flow through stenosis for R^ = 10, Re = 1 for . 7T 37T t = 0, —, 7r and — 95 IX 4.47 Streaklines for pulsatile flow through stenosis for Rw = 10, Re = 5 for t = 0, —, 7r and — Li 96 Lt 4.48 Streaklines for pulsatile flow through stenosis for R^ = 10, Re = 10 for t = 0, —, 7r and — Li 97 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" 37T * = 0, - , 7T and —Li 99 Lt 4.51 Streaklines for pulsatile flow through stenosis for R^ = 10, Re = 100 7T 37T for t = 0, - , 7r and — . . . . 100 4.52 Streaklines for pulsatile flow through stenosis for Rw = 10, Re = 150 for t = 0, —, 7r and — . . . 101 Li Lt Li Li 4.53 Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 0 3TT 7T (linear) for £ = 0, —, 7r and — Lt 102 Lt * 4.54 Streaklines for pulsatile flow through stenosis for Rw = 20, Re = 1 for „ ir t = 0, - , , 3ir 7T and — . 103 4.55 Streaklines for pulsatile flow through stenosis for Rw = 20, Re = 5 for *t = 0, - , x and — = 0, —, 7r a n d — 104 105 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 — 4.59 Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 100 „ , 37T for t = 0, - , 7T and —- 107 108 x 4.60 Streaklines for pulsatile flow through stenosis for i ? w = 20, Re = 150 for t = 0, - , 7T and — Lt 109 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 Ill 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 D.l 112 Use of transformation (x = x(x', y'); y = y(x', y')) to determine stresses along a curved wall 130 xi Notation Lower case G r e e k 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 U p p e r case G r e e k symbols T T s fluid boundary traction boundary Tu kinematic boundary Te finite element boundary finite element traction boundary T„ finite element kinematic boundary A increment Q fluid domain fl (prefix) finite element domain e Lower case R o m a n symbols a a(y) pre-stenosis length, pre-orifice length steady part of pulsatile flow b b(y) post-stenosis length, post-orifice length 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 na axial cosine of unit outward normal to the boundary n2 radial cosine of unit outward normal to the boundary p &Padd pressure additional pressure drop xiii r r(x) radial coordinate in three dimensions instantaneous radius t time u axial component of velocity UQ v v r characteristic velocity radial component of velocity radial component of velocity in three dimensions azimuthal component of velocity in three dimensions v z 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 B (y) magnitude of cosinusoidal part of pulsatile approximation Q B(<) Co(y) C(t) magnitude of cosinusoidal component magnitude of sinusoidal part of pulsatile approximation magnitude of sinusoidal component C° continuity of zeroth order derivatives C1 continuity of first order derivatives F r body force in radial direction body force in azimuthal direction F z J Jo body force in axial direction Jacobian matrix 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 at ox oy total material derivative = — + u——V v— — Dt ' 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 encouragement 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 CHAPTER 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 developments, rather t h a n by academic questions. It was only after the advent of digital electronic computers, that the solutions for such problems could be developed w i t h success. In the beginning of fluid mechanics research, l i m i t i n g situations—offering permissible simplifications of the governing equations to solvable ones—were primarily the goal of research. Thus came the P r a n d t l ' s lift line theory, the P r a n d t l - M e y e r flow, or the Blasius series solution. E v e n though these problems were relatively quite simple, a l l of t h e m already involved numerical work and a substantial amount of t i m e i n which those solutions were derived. Prandtl's lifting line theory yielded an integro-differential equation for the circulation, the P r a n d t l - M e y e r flow a transcendental equation i n terms of M a c h number to be solved for the P r a n d t l - M e y e r 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. T h e first electronic computers immediately enabled the solution of the i n i t i a l value two point boundary value problem for two-dimensional boundary layers and of the potential equation for flows w i t h 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 i n 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 t r y to 'simulate' the physical phenomena i n some approximate way and to improve the approximation i n successive attempts. C o m p u t a t i o n a l fluid dynamics has always been at the forefront of the development 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 m i x e d hyperbolic-elliptic character of the governing p a r t i a l differential equations, 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. A n o t h e r reason that follows from above is the size of typical problems. A problem i n structural mechanics is considered 'large' if it exceeds 5 x 10 3 nodes, whereas large i n fluid problems means more than 5 x 10 5 grid points [12]. For about 25 years the computational fluid dynamics grew around F i n i t e Difference Methods, as these were relatively simple to understand and code, easy to vectorize and the structured grids typically associated w i t h t h e m described appropriately 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 i n time that unstructured grids and F i n i t e Element M e t h o d s — 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 still prefer the finite difference methods. Since the i n i t i a l development efforts for the finite element methods for fluid mechanics, m u c h attention has been devoted to the treatment of the p r i m i t i v e variable (velocities and pressure) formulations. T h e dominant approach used i n the finite element methods is the m i x e d order interpolation scheme which attempts to eliminate the spurious pressure modes [31], i.e., the tendency to produce 'checkerboard' pressure distributions [24]. T h i s scheme is analogous to the 'staggered g r i d ' used i n the finite difference methods. However, m i x e d order interpolation is not totally effective at eliminating 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; R i c e 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 i n an approximate pressure equation and allows uncoupling of pressure and velocity solutions. Numerous other suggestions have been made. For the present investiga- t i o n , we have used a m i x e d order interpolation w i t h quadratic bilinear element where pressure is not constant over the element; along the element boundaries the variation is linear whereas i n the domain the variation is linear plus one quadratic cross term. T h i s eliminates the above mentioned 'checkerboard' pressure distributions. 1.2 Stenosis and Pulsatile Flow T h e world "stenosis" is a generic medical t e r m which means a narrowing of any body passage, tube, or orifice. T h u s , 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.' P a r t i a l occlusion of blood vessels due to stenosis is one of the most frequently occurring abnormalities i n the cardiovascular system of humans; it is a well established fact that about 75% of all deaths are caused due to circulatory diseases [2]. O f these, atherosclerosis is the most frequent. Considerable studies related to the circulatory flow i n blood vessels were given major attention towards the beginning of the present century. Quite a few analytical and experimental investigations related to b l o o d flow w i t h different perspectives have already been carried out. Interest i n studies of this particular domain of biomechanics has increased w i t h the discovery that many cardiovascular diseases are associated w i t h the flow conditions i n the blood vessels; one major type of flow disorder results from stenosis. Diagnosis of arterial diseases is usually carried out by the m e t h o d of X - r a y angiography which involve considerable risk of m o r b i d i t y [3]. So any attempts to develop non-invasive techniques for detecting arterial disease like stenosis are very worthwhile. A s a result, numerical modeling of arterial b l o o d flow has long been of interest, but development of realistic models has occurred only i n recent times [22]. Causes of the slow progress include nonlinearities i n the model equations, unsteady flow i n the vascular system, and elastic properties of arteries. Most of the studies have been carried out assuming that b l o o d is a Newtonian fluid. However it has now been accepted that b l o o d behaves more like a non-Newtonian fluid under certain conditions, i n particular at low shear rate. Srivastava [28] has presented a study of the effects of stenosis on the flow of b l o o d when b l o o d is represented by a couple stress fluid, a special type of non-Newtonian fluid which takes particle size into account. T h i s , however, is beyond the scope of the current investigation. Y o u n g [35] has presented a survey of some of the early works done i n this field. Chapter 1: Introduction and Literature Review 5 Over the years, D a l y [4], Deshpande [6,7], Yoganathan et al. [34], P h i l p o t et al. [21], Porenta, Y o u n g and Rogge [22], Forrester and Y o u n g [8,9] and others have shed considerable light on the complex problem of blood flow through an artery i n the presence of stenoses. D a l y 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. L o c a l flow reversal was observed i n 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 i n a rigid tube utilizing full Navier-Stokes equations i n cylindrical coordinates. L a m i n a r boundary conditions, n o r m a l l y applicable only at x = ± o o , were applied at finite values of x by using a transformation of the type r\ = t a n h ( K x ) , where is K is a constant chosen to group the grid points i n an efficient manner. T h e stenosis geometry was modeled as a full cosine wave. Yoganathan et al. obtained steady flow velocity measurements i n a pulmonary artery model w i t h varying degrees of pulmonic stenosis. P h i l p o t et al. used a 7 - m W He-Ne laser as the light source for flow visualization i n a P u l m o n a r y artery model. T h e y photographed the setup over the entire systolic period, thus observing the development and dissipation 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 numerical solution after 3 cycles, independent of i n i t i a l values. T h e arterial cross-section was of the form a • exp(— ax), where a is a small positive number. L o w frequency behavior was approximated w i t h seven cosine and seven sine components of the generalized 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 p o l y n o m i a l 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 critical Reynolds number w i t h height of stenosis. 1.3 Analysis in Time Domain Pulsatile flow i n the arteries is highly variable w i t h time and to understand all 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. T h i s , however, is extremely costly i n terms of C P U time and storage since, at each step, we must solve for the effects of the quadratic, nonlinear, convection terms. M c L a c h l a n [14] and Nayfeh and M o o k [17] have pointed out that the solution of the differential equations w i t h such nonlinearities has a steady part as well as an oscillating part. T h e 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 equations 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 i n [14,17]. P a t t a n i [19,20] has modified the method of averaging and shown that it can be used successfully for this k i n d of problem. Using the Floquet theory, he has also shown the solution to be stable. T h i s method is used here. Chapter 1.4 1: Introduction and Literature Review 7 Scope of Present Investigation T h i s thesis describes the development work on extending the finite element method to cover the effects of stenoses on pulsatile flow i n arteries. T h e problem configuration of interest is that of an axisymmetric, or 'collar-like' artery. T h e pulsatile flow is approximated as a combination of steady, sinusoidal and cosinusoidal velocities. In Chapter 2, the fundamentals of fluid flow i n cylindrical coordinates are presented. T h e chapter begins b y stating the governing differential equations for the fluid, derived i n [26] based-on the principle of conservation of mass and moment u m , and choosing non-dimensional parameters to effectively characterize the flow through a pipe. K i n e m a t i c 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. T h e 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. N u m e r i c a l results are described i n Chapter 4. Besides the pulsatile flow through stenoses, several other l i m i t i n g problems are also solved to check the code. These results compare favorably w i t h published analytical and experimental results i n 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. CHAPTER 2 Derivation of Governing Equations In this chapter, the governing equations for axisymmetric fluid flow are presented. T h e velocity-pressure primitive variable form of the Navier-Stokes equations for two dimensional incompressible viscous flow is used. T h e equations of m o t i o n are nondimensionalized to effectively characterize the flow situation. Finally, an approximate representation for pulsatile flow is presented. A l l the analysis is based on incompressi b i l i t y assumptions which hold well for liquids. 2.1 Conservation Equations and Boundary Conditions T h e 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 initial conditions, i n 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 T and T s , respectively, as u shown i n figure (2.1) If r, <p and z denote the radial, azimuthal and axial coordinates, respectively, of a three-dimensional system of cylindrical coordinates and v , T v^, v denote the velocity components i n the respective directions, then the governing z equations for incompressible fluid flow are given b y Schlichting [26] as 8 er 2: Derivation r of Governing Equations u Figure 2.1: Fluid Domain (dvr [lT dvr 0 + '-dT- v „ ^ dt T * + (d 2vT dr dr r dvz '-dr- = vr 1 dr 2 r r2 ' r2 d<j> r dz J 2 V4, + *-d7) v 1 dvr V dr V dvr\ 7-d+-T r5<A + V& 2 + dp (dvz [-d7 Vfidvr — or 2 r2 r dr + dvz dvz — r d<j> + "z o dz r2 2 dv* d 2vr d<f> 2 dcf? 2 r2 r2 d<f> d<f> ay dz 2 dz 2 Chapter 2: Derivation of Governing Equations 10 The stress components assume the form cr = -p + 2n r <r, = - P — or Tr4> = V dr \ r J r d(p 10 f ldv<f> , VA + 2u [--Qj + j ) „2 = -p + (2.5) /<9ur 2n— <9u2\ where, p is the fluid density, v = — is the kinematic viscosity, p is thefluidpressure, 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 axisymmetricfloware 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 Dt dt ox dy . where — = — + u——\- v— Du d 2u 1 du^ dy ydy. pdx \dx Dv I dp dv pdy 2 I dv 2 + \dx 2 dv 2 + v , d 2u Dt Dt , is the total material derivative. 2 ^ dy 2 + " ~ p dy vN o y2 (i,y)eO, t>0 (2.6) Chapter 2: Derivation of Governing Equations 11 E q u a t i o n (2.4) reduces to du dv v ox oy y (2.7) t>0 (z,y)6f2, w h i c h is the equation of continuity for axisymmetric flow. . .,. Adding d (du dv v\ I — h -5—I— ox \ox oy yJ , and d (du ^~ oy \ox dv v' + 7T + ~ oy yt to the right hand sides of the 1st and 2nd of equations (2.6) we obtain Du _ ~Dt ^ 2u Idp dx 2 p dx Dv _ Ut = Idp ~o"dy d 2v + dx V 2 d 2u d2u dy 2 dxdy d 2v d 2u + 2 — + dy 2 dxdy 1 (du <9t>N y \dy + 2 'dv y \dy dx. (x,y)eft, t>0. y (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. O f the total six stresses given by equations (2.5), only the following three are n o n - t r i v i a l . These are also referred to as constitutive n rt _ • xy (dv ^ \dx relations. du dv (2.9) du dy, B o u n d a r y conditions are given on the velocity components over Tu, as the kinematic boundary and on traction over Ts, boundary or traction boundary. o~ n + T n = X Txynx = Y x xy referred to as the mechanical These conditions are u = U, u = V x referred to 2 + o-yn2 (x,y)er„, t > 0 (2.10) (x,y)eT3, t>0 Chapter 2: Derivation of Governing Equations 12 Substituting equations (2.9) into equations (2.10), we obtain the boundary conditions i n terms of pressure and velocity components. = u, u „ du 'du ^ {dy (du dv ^ \dy dx ni + dv^ + dx, -p + 2fi dv dy v = y n2 = X n2 = Y (x,y)eru, o o (2.11) (ar,y)€r„ O O where, U , V are the specified velocities on T , u X , Y are the specified tractions on F , and 3 n 1 ; T i 2 are the direction cosines of the unit outward n o r m a l to the boundary. 2.2 Non-dimensional Form of Governing Equations T h e Navier-Stokes equations, the continuity equation and the boundary conditions can be non-dimensionalized by choosing suitable characteristic values for the coordinates, the velocities and the time-scale variables and defining non-dimensional parameters. Choosing the non-dimensional variables as x _ y x = — y = — t = ut Ro Ro u v „ p u = — v = — p = fj.u0/Ro UQ UQ (2.12) 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 „ du n ( du R„ — + Re [U— dt \ dx of Governing du\ + V— dy J Equations dp = -7p „ dv _ / dv dv\ Rw— + Re [u— + v— ) = dt \ dx dy J ^d 2u + dx dp dy 13 dx 2 d 2v + dx 2 du dv v T + T + dx oy y- = d 2u d 2v 1 fdu dv^ dy dxdy y \dy dx/ 2 ^d 2v d 2u 2 (dv v\ dy 2 dxdy y \dy yI 0 (2.13) A l l quantities are their respective non-dimensional values and the tilde has been o m i t t e d for convenience. T h e two non-dimensional parameters, R^ = frequency Reynolds number and RE pletely characterize the flow. UQRQ = UJR 2 —-, v the , the steady Reynolds number, com- v It should be noted that the pressure has been non- dimensionalized w i t h 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 variables <7 = < T x <j UUQ/RQ UUQ/RQ f < T y = UUQ/RQ Y = XY UUQ/RQ UQ T x y UUQ/RQ (2.14) UQ into the constitutive equations (2.9) and boundary conditions (2.11), we obtain n = - a x p du d-x dv + 2 n <Ty = _ Txy = - p d7dy l + 2 dv [fa du + dy', ( 2- 15) Chapter 2: Derivation of Governing Equations 14 and u = U, du du \dy (du dv' \dy dxt nx + dv dx) dv n2 n2 v = V = = (x,y)eTu, £>0 (x,y)eTs, t>0 X (2.16) 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 i n the domain f2. T h e first difficulty encountered b y any solution method arises naturally from the nonlinear non-self-adjoint convective terms i n the Navier-Stokes equations. A further significant difficulty arising i n all numerical formulations of the governing equations i n p r i m i t i v e variables ( u , u , p ) representation is due to the interaction of the pressure and dilatation terms. T h i s difficulty arises from the fact that the continuity equation merely represents a constraint on the divergence of the velocities rather than a full t h i r d equation coupling the pressure w i t h the velocities. A s a result, it is difficult to develop solution methods which are capable of determining the primitive variable (w,u, p), simultaneously, w i t h 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, convective terms. T u a n n and Olson [31,32] have shown that the only appearance of p i n a restricted variational functional is coupled w i t h the dilation t e r m 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 equations, so that either the system consists of the Navier-Stokes equations and an induced equation which involves pressure i n a relatively dominant role, or else, the divergence constraint is satisfied explicitly and exactly, dropping the primitive variables approach. Alternatively, one fourth order equation can be obtained for the stream function. Relative merits of these formulations i n 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 u ldyj - — y oy _ld_V> y dx = (2.17) Hence, dx dy y dx \y dy J 1 d yb 2 y dxdy = t dy \ y dx J y \ 1 dyj I d yj 1 dtp y 2 dx y dxdy y 2 dx 2 y dx) 0, i.e., the continuity equation is satisfied point-by-point. T h e vorticity £ is given as du dv , Substituting equations (2.17) into (2.18), we obtain i fd 2i> a _0_ia_v>\ y\dx 2 dy 2 2 y dy) ~ du_dv dy - (2 19) dx T h i s is called Poisson's equation and is.used to compute the stream function, once the p r i m i t i v e variables have been determined. (Right hand side is known.) CHAPTER 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 p r i m i t i v e variable u, v and p. Galerkin's procedure is used to minimize the error i n approximation by m a k i n g 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 solution. T h e Galerkin's procedure is also used to find a representation of pulsatile flow. Brief details of computation of stream functions and p l o t t i n g of streamlines and streaklines are also presented 3.1 Discretized Form of Governing Equations Olson and T u a n n [31] have shown that the finite element interpolation for the presdu dv sure should be no higher than that for the dilation t e r m ——I- TT" i n 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). T h e highest order of derivative i n the equations (2.13) is two for the velocity components (u,v) 16 and one for the pressure Chapter 3: Finite Element p; as a result C Formulation 17 continuity is required of the interpolation functions for (u,v). 1 This difficulty c a n , 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. N o w , only C ° continuity is required of the interp o l a t i o n functions for equations (2.13). Furthermore, for the present investigation, isoparametric elements are used, i.e., the same shape functions are used for interp o l a t i n g the variable (u,v,p) to {x,y) coordinates. and for transformation from (£,77) n a t u r a l coordinates T a k i n g these conditions into consideration, curved isopara- metric elements shown i n figure (3.1), w i t h quadratic interpolation for the velocities and bilinear interpolation for the pressure are used to carry out the finite element discretization of the governing equations. T h e velocities u, v represented b y where Ni are the quadratic N, = u = NiUi{t) v = NiVi(t) 8 isoparametric shape functions given as 4 < - 0 ( l - ' 7 ) ( l + * + '7) N2 - 1 ( 1 + 0 ( 1 + »7)(1 N, 1 - 0(1 + N5 = N6 |(i + 0(i N7 = N8 ^(l-O(W). + ( - V) -v 2) (3.2) Chapter 3: Finite Element Formulation Figure 3.1: Isoparametric Element (x,y) space. (a) Element i n (£,77) space, (b) Element Chapter 3: Finite Element Formulation 19 T h e pressure p is represented by p = MjPj(t) ;' = ! , . . . , 4 where Mj are the bilinear isoparametric shape functions given as Ma M 3 = 1(1-0(1-17) M2 = i(l+0(l+»7) M4 = 1(1 + = 0(1-1/) (3.3) 1 ( 1 - 0 ( 1 + .7). Substituting the interpolation functions (3.2) a n d (3.3) into the governing equations (2.13) a n d the boundary conditions (2.16), we obtain Tlx = + Rw—g^-Uj R " 'd 2Ni R 'v +R Ti —QJ-Vj + lt d*Nj dx 2 72, = 72n = K2 = V ' dy 2 d 2 N j dy v dN, ^MdNk 9 2 N j 2 V dNj ~dT •Uj + 1 (dNj dNk tt dNj d 2Ni dMj e 1 n dMj dy dNj _ . +dNj — j P j + ^MdNk d 2Nj dN w dNk e Uj -dx~ V> 1 j dxdy u j y \ dy - j ^• y j (3.4) +Nj ni + nx + 'dNj dNj + IT* y n2 — X n2-Y where 7Z\, 72 2 and 72.3 a r e the errors i n approximating the x - m o m e n t u m equation, the ^-momentum equation and the continuity equation, respectively; 72-ri a n d 72p2 are the errors i n 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 m i n i m i z e the errors, we make the residuals orthogonal to their base vectors—i.e., the shape functions JVt- and Mi—over the finite element. T h i s leads to II K1Nidto+ JJ w [ KnNidT J r» = 0 // K2Ni<m+ JJ u e I KriNidT J r» = 0 n^MidSl = JJ where f i e is the domain and T e s (3.5) 0, is the traction part of the boundary of the element under consideration. Substituting equations (3.4) into (3.5) and integrating the appropriate terms b y parts, we obtain 0 0" 0 0 0 0 Uj < 0. ' fjk J k s + Re < u u ijk Uj Vk S -Pij ij -Pij K ji K K .Pj. Uuv ij ij K 0 .-Pji + + o 8ijkvjuk tijkVjVk 0 ' ' • = < 0 .0. ' Uj " y < Vj .Pj . (3.6) Chapter 3: Finite Element Formulation 21 where dot denotes differentiation w i t h respect t o time and the arrays are Luu 'mm, //,. n \ ax e 'dNidNj = 11 Q \ s V.uv JJ °ijk n e ^ + c7z i dy ^dN^d^ ox ox dx dn oy J oy oy 2_ y z 1 i,j,k = 1,...,8 dy = If wf-r JJ u e dx an i = l,...,8 IIQE NTN3 dy * - IL ( ^•Mi dx j = l,...,4. + -NiMA dQ y 7 (3.7) Vector f on the right hand side of equation (3.6) is the consistent load vector and is given as ( f Xiv.dr-) J r; f =< / Yiv -dr 4 J r« 3.2 Boundary Conditions For a l l the nodes that are on the kinematic boundary T , the nodal values are forced to u the specified values of the corresponding variables. Thus, i n the present investigation, the kinematic boundary conditions are satisfied exactly—node-by-node, that is. A s is evident from equations (3.5), the traction boundary condition, i n general, is satisfied only i n the " m e a n " — TZn and TZr2 are not identically zero. T h e specified stresses X and Y are integrated over the traction boundary, T , to give the consistent load e s vector f. Chapter 3: Finite Element Formulation 22 T h e 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 a n d the traction boundary conditions are approximately satisfied. E q u a t i o n (3.6) reduces to ' fjk J k s [M]{d} + [K]{d} + i2e< u + Sljk i k u v u 6fjk jVk + Si vjv u jk k ' \ = 0 (3.8) 0 where, M is the mass m a t r i x given b y 0 M = 0 0" 0 0 0 0 K is the stiffness m a t r i x given b y K K = ij Uuv K ji •L.UV ij -Pij ' Uw -Pi K K ij 0 and u. d = [Pj ) is the n o d a l vector of unknowns. 3.3 Steady State Periodic Solution Due to the presence of the quadratic nonlinear terms i n 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. T h i s , however, is very costly and requires large storage and C P U time. Quite often the purpose is well served w i t h an approximate slowly varying periodic solution. Nayfeh and M o o k [17] have classified the available techniques for deter- m i n i n g 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. P a t t a n i [19] has modified the method of averaging, as given by Nayfeh and M o o k [17], to obtain the steady streaming part, as well as the slowly varying periodic part, of the solution from equations (3.8). T h e starting point i n this method is to assume a form of the solution as d = A + B(*)cos* + C(*)sini (3.9) where B ( i ) , C(t) are assumed to be slowly varying functions of non-dimensional time t. T h e first t e r m A represents the steady streaming part of the solution which naturally arises for systems of equations w i t h 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)cosi + B(£)cosi + C(f)sini (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). T h e n to obtain the three sets of equations for the three sets of unknowns A , B ( i ) and the resulting equations are 1. averaged over the period 27T, C(t), Chapter 3: Finite Element Formulation 24 2. multiplied by cos t and averaged over the period 2w, and 3. multiplied 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 w i t h time than the functions sin/j, COST; and time t itself. Therefore, B(t), C ( t ) , B(t) and C(t) are considered to be constant while performing the averaging. T h i s process yields "0 0 0 M 0 .0 0 M. 0 " 'A ' < B > + .c. (- 6fjk(A»At "K 0 0 K . 0 -M 0 " M 'A ' < B • K. .c. + B V B J I / 2 + C y C J / 2 ) + 8 yjk(A»Al « % , ( A ? A £ + B ? B £ / 2 + C)Cl/2) + Sljk{K)A.l + + B]B%/2 + C«Q/2) ) + CJQ/2) 0 -^•*(AJBJ| + AJfBJ) + « ? i f c ( A J B H + A £ B ? ) . +i?e < ^ * ( A J B j ; + A J B J ) + Sf (A)B% ik + AlB]) 0 % ( A ? C £ + A J C J ) + - ^ ( A J C X + AJfCV) ^ ( A j q ; + AJCJ) + ^fc(Aycjs + A j q ) 0 '0 <0 .0 (3.11) Chapter 3: Finite Element where Formulation 25 (By l r Ay A = < Ay (cy) B = <!By cy p y. B J A , B(rj) and C ( i ) . T h e superscripts are the average values over one period (27r) iof u, v and p indicate that the amplitudes are associated w i t h the u, v and p degrees of freedom respectively. T h e steady state solution corresponds to the singular points of the autonomous system of equations when B = C = 0. T h i s results i n a set of nonlinear algebraic equations for A , B and C which are solved using the modified Newton-Raphson iteration procedure as outlined i n A p p e n d i x A . T h e equations obtained for the increments to the solution vector are of the form [ T ] { A x } = {-f} (3.12) where, [T] is the tangent stiffness m a t r i x , { A x } is the incremental solution vector and {—f} is the unbalanced load vector. Details of equations (3.12) are given i n A p p e n d i x B . T h e above matrices are calculated for each finite element and are assembled into the corresponding global matrices. T h i s 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 Theological 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 crosssection after it has sufficiently "stabilized" over a long development length. B l o o d consists of formed elements, red blood cells, white blood cells, platelets etc. suspended i n plasma. A l t h o u g h 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 i n larger arteries, blood is generally assumed to behave as a Newtonian fluid. T h e 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 distribution across the artery cross-section. A r t e r i a l b l o o d flow is highly pulsatile w i t h the instantaneous flow rate varying over a wide range during a flow cycle. Y o u n g [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, H(t) and H(t + a) = = 1 0 < t < 0 7T < 7T, I t < OCK, H(t). To get r i d of the discontinuity arising from H(t) i n figure (3.3), we approximate Chapter 3: Finite Element Formulation Figure 3.2: Typical arterial flow waveforms (Young, 1979) 27 Chapter 3: Finite Element Formulation 28 Hit) 1 _ IT 0 OCT ( a + l)7T 2a7T ' (2a+l)7T Figure 3.3: F u n c t i o n H(t) UQ(y, t) as Uo{y,i) ~ A0(y) + B0{y)cost + C0(y)sint ' (1 cos t s'mt) < B0(y) . A0(y), B0(y) ' My)' > = < l cos t . sint , c0(y). T " > 'i4o(y) ' < B0(y) . C0(y) . and Co(y) are determined by m i n i m i z i n g the error i n approximation by Galerkin's procedure. T h e resulting residual ' 7lu = AQ(y) + B0(y) cos t + C0(y) sin t - [a(y) + (b(y) c o s t ) - H(t)} is made orthogonal to the base vectors 1, cost and s i n t . T h i s gives ' 1 ' < cos t k sint , ' > < 1 T " cos t > sint , • 'A0' dt< Bo . Co . •- / ./ 0 ' < 1 ' '0 ' cos i - (a + ( 6 c o s i ) • H(t)) dt = < 0 • , sint , .0. Chapter 3: Finite Element or, " COS 0 sin t cos 2t sin t cos t sin t cos t 1 sm 2t t smt ai / 29 cost 1 ai / Formulation 0 or, sin car smai OJ7T —— H 2 t - dt + b / < J 0 sin 2a7r , sin t cos 11 Bo > = a < 4 air _ 2 \ dt cosH 2 1 Q7T — cos OCK sin2ax sin a?r 2 — COS OJ7T . Co . cost fir COS sint Q7T d t < B0 > sin 2a7r 0 sin air ir — COS OJ7T 4 K 2 ' (3.13) E q u a t i o n (3.13) can be solved for different values of a. In particular, . 37T tor a = —, 2 ' '1 ' Mv)' < B0(y) • = ' 0.433897 ' - a(y) < 0 . + b(y) < 0.0446932 .0. . C0(y) . k for a = 2ir, 0.675718 t 1 (Mv ir Bo(y > = a(y) < • + b(y) < 0 1 Co(y ' 2 ' and for a = Sir, ' Mv)' < B0(y) » . C0(y) , B o t h a(y) ' 1 " = a(y) < 0 . + b(y) « .0. '0.18091221 ' 0 . 0.29494257 , 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 b o t h 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 i n d of order zero. For ease of programming, though, b o t h a(y) and b(y) are assumed to vary parabolically across the cross-section. T o 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. T h e family of streamlines at any time t are solutions of the equation 1 y dx I dy 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. T h e path of a material element of fluid does not, i n general, coincide w i t h a streamline. T h e pathline and streamline are the same only when the fluid m o t i o n 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. T h u s , when a m a r k i n g material is dropped i n the flow, the visible p a t h produced thus is a streakline. T h e family of streamlines and streaklines are found by solving the Poisson's equation (2.19). After the governing equations for the p r i m i t i v e variables (u,v,p) have been solved, the vorticity can be computed at any point. du dv dy dx I n order to be consistent, the finite element i n - terpolation functions chosen for the stream function ?/> to discretize the Poisson's equation (2.19) are the same as those used for the velocities. T h i s also simplifies the computation and setting up of the m a t r i x equation. T h e stream function is Chapter 3: Finite Element Formulation 31 represented b y 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 dNj dx dy dNj dx 3 into equation (2.19), we obtain 1 TU = d 2Nj d 2Nj dx dy - 2 _ IdNj dN< dNj y dy 2 dy where TZ^, is the residual due to error i n approximation. To m i n i m i z e the error i n approximation, we make 7Z^ orthogonal to the base vectors, i.e., the shape functions Jj' Tl^NidQ. = 0 or, JJ 1 \d 2Nj d 2Nj 1 dNj dNj dx 2 dy 2 y dy dy ny e dNi i dvj NidCl = 0. Integrating the appropriate terms by parts and simplifying, we obtain NW (3.15) (Pi) = where the arrays are * - II 'ldWdN± y dx dx 1 dNj dNj y dy dy 'dNj dN4 n I dy U^ dx e e element and 9 = ij e s is the traction part of the boundary of the d^_ dn dN3 dx = 1 , 2 , . . . , 8. NidCl + I -gNdT J r%y VJ Cl is the domain of a single element, T dn dNj n2 Chapter 3: Finite Element Formulation is the tangential derivative of tp specified on the traction part of the boundary r*. 32 The above process is repeated for each element and the matrices are assembled into the corresponding global matrices. T h e resulting set of sparse, linear, algebraic equations are solved using the University of Waterloo solution package S P A R S P A K . Substituting back the solution i n equation (3.14), we obtain the finite element representation of the stream function ^ i n the fluid domain. CHAPTER 4 Numerical Results In this chapter numerical results are presented for different flow configurations. T h e results are obtained for various values of Re and i n the viscous flow regime and, wherever possible, these results are compared w i t h analytical, experimental or numerical results published i n the literature. F r o m 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. T h e 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. T h i s results i n a set of nonlinear algebraic equations i n A , B, C. These equations are solved using the modified Newton-Raphson procedure as outlined i n A p p e n d i x B . A l l the computer programs are implemented on the A m d a h l 5840 and on a departmental Sun 4/260 computer at the University of B r i t i s h C o l u m b i a . Double precision (Real*8) arithmetic is used throughout to reduce the effect of round-off errors. T h e 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, i n conjunction w i t h 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 i n Appendix C. 33 Chapter 4.1 4.1.1 4: Numerical Results 34 Steady Flow Through a Pipe Theoretical Results T h e steady flow through a straight pipe of circular cross-section is the simplest case w i t h rotational symmetry for which an exact solution can be obtained. L e t the x-axis be selected along the axis of pipe, figure (4.1), a n d let y denote the radial coordinate measured from the axis outwards. For fully developed flow the velocity component i n the radial direction, denoted b y v, is zero; the velocity component parallel to the axis, denoted b y u, depends on y alone a n d the pressure is constant at every cross-section. T h e governing equations (2.13) reduce to d 2u 1 du dp dy y dy dx 2 2R Figure 4.1: P a r a l l e l F l o w w i t h P a r a b o l i c velocity distribution T h e second and t h i r d equations of (4.1) give p = p(x) a n d u = u(y). Rearranging the first equation, we obtain d 2u 1 du _ dp dy 2 y dy dx' ^ Chapter 4: Numerical Results 35 T h e left h a n d side of equation (4.2) is a function of y alone and the right hand side is a function of x alone; b o t h these conditions can be true if and only if each side is equal to a constant. Hence, d 2u 1 du dy y dy 2 dp — = constant dx or, ld_ du ydy 'dy dp — — = constant dy or, du dy Spy 2 faJ , + C r l or, C2 For u to be bounded along the axis of the pipe, C\ = 0 and for u = 0 at y = R, R2 = — — • T h i s 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 W e note that the exact solution obtained above is quadratic i n y for the velocity component u and linear i n x for the pressure p. T h i s 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 i n 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. T h e 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 = u = 0,v = 0 7^ (l-y2) u = 0 t; = 0 p = 0 4.0 v = 0 Figure 4.2: One Element Representation of the F u l l y Developed F l o w du T h e 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. A l s o by letting u(0,y) = (1 — y ) , i m p l i c i t l y we let — = —4. T h i s w i t h dx 2 the condition that p = 0 at x = 4 gives p = (16 — 4x). T h e finite element results are given i n table (4.1). T h e y match the exact solution up to the precision of the R e a l * 8 arithmetic used. N o iterations were required to converge to the true solution. 4.2 4.2.1 Oscillating Flow Through a Pipe Theoretical Results T h e 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. T h i s 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 i n i t e Element Solution w i t h E x a c t Solution Node 1 2 4 5 6 Coordinates (0.0,0.0) (4.0,0.0) (0.0,1.0) (2.0,0.0) (4.0,0.5) Variable P u P u u F E solution 16.0000000001713 1.0000000000046 15.9999999999584 1.0000000000055 0.7500000000023 E x a c t solution 16.0 1.0 16.0 1.0 0.75 theory is presented i n 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 . T h i s means that the axial velocity component u is also independent of x , i.e., u = u(y,t). T h i s simplifies the continuity equation to dv v 7T + - = 0 dy y or, ydy or, yv = constant = g(t) for some arbitrary function g{i). T h e y—momentum equation reduces to T h a t is impossible unless v = 0. dp equation reduces to — = 0, i.e., p = p(x, t); the x—momentum j^9u _ dt dp ^ d 2u ^1 dx dy 2 du y dy (4.3) which is 'exact' as it implies no additional simplifications. T h e boundary conditions are u(y,t) = 0 for y = R at the w a l l . Assume that the pressure gradient caused by the motion of the piston is harmonic and is given by dp --= P o cosi, Chapter 4: Numerical Results 38 where p0 denotes a constant. For this problem it is convenient to use complex notation and to put dp - n a t t r i b u t i n g 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. A s s u m i n g that the velocity function has the form u(y,t) = f(y)e lt, we obtain, % e 2 dy 2 Substituting equations (4.4) into the equation (4.3), and canceling the common t e r m 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 v = f'(y) = yy(-iRJ)dfiv) dn dn dy df(rj) dn f"(y) = (-iR.) d 2f(y) drj 2 ' T h i s gives, „2d 2f(v) , df(n) which has the particular solution fp(n) A(7?) = AJ0(r)) + BYo(n), 2 f 2 = ^—*^T^ 2 a n d P o \ homogeneous solution 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 For the solution to be bounded, B = 0 since YQ [yy/—iR^) gives 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) T h e above can be further simplified b y splitting the bracketed quantity into sine dp and cosine terms. Hence for cosinusoidal pressure difference — — = po cos t, ox u(y, t) = ^-(0 Rw cos t + a sin t) 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 Jo{y^iRZ)\ j \ Jo (RV 1^) J' A general discussion of the nature of the solution for arbitrary values of R^ is i m possible owing to the presence of Bessel's function w i t h a complex argument, but for the two l i m i t i n g 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 • J0(z)«,/Ae-rl/2 Chapter 4: Numerical Results Figure 4.3: Flow in a pipe under oscillating pressure 40 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 i n figure (4.3). 4.2.2 Finite Element Results W e choose a pipe of radius R = 1 and length / = 4 represented by a 5 x 8 grid as shown i n figure (4.4). W e compute the fully developed velocity distribution using the series solution presented i n the last section and specify these values at the entrance x = 0. T h e program seeks its own approximation to the flow downstream and stablizes before the exit x = 4. T h e velocity distributions thus obtained at the exit are compared w i t h the exact solution. J, u = 0, v = 0 u = ( a cos t -+- (3 sin t) 1.0 v = 0 v = 0 p = 0 •*-x 4.0 Figure 4.4: F i n i t e element grid for oscillating flow through a pipe T h e exact solution presented i n 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. T h u s u(y, i) = - j - ( a c o s t + Rw f3sint) where a and /? are functions of y and are computed using the series solution. T h e program is r u n for three different frequency Reynolds numbers Rw = 1, 10 and 20. Chapter 4: Numerical Results 42 T h e results are compared i n 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 E x a c t Solution y 0.0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0 a 0.2419938822 0.2382643889 0.2270649264 0.2083631157 0.1821068921 0.1482273619 0.1066428141 0.0572639012 0.0 0.0454856408 0.0445439938 0.0417627709 0.0372732688 0.0312947819 0.0241352748 0.0161922692 0.0079538898 0.0 F i n i t e Element Solution a P 0.2419887660 0.0455076871 0.2382631459 0.0445506155 0.2270629093 0.0417728360 0.2083631013 0.0372747261 0.1821058316 0.0313012409 0.1482278532 0.0241334887 0.1066424134 0.0161955928 0.0572645050 0.0079488352 0.0 0.0 Table 4.3: Comparison of exact and finite element solution for oscillating flow for Re = 10.0 R y 0.0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0 = 1.0, po = 1.0, E x a c t Solution a P 0.0456409109 0.1109031291 0.0460493859 0.1091161974 0.0470648468 0.1037100366 0.0480518359 0.0945759605 0.0479364625 0.0816249082 0.0451974541 0.0649275362 0.0378803750 0.0449123607 0.0236620859 0.0226194555 0.0 0.0 Ru, = 10.0 F i n i t e Element Solution a 0 0.0455297253 0.1108901990 0.0460123330 0.1091272587 0.0470122009 0.1037222211 0.0480434236 0.0945908185 0.0479068297 0.0816443909 0.0452096686 0.0649315983 0.0378715381 0.0449288345 0.0236821472 0.0225941990 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 E x a c t Solution y 0.0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1.0 a 0.0042213053 0.0050090036 0.0072867953 0.0107517756 0.0147491753 0.0180745128 0.0187627664 0.0139564340 0.0 0.0601667061 0.0598214600 0.0586030834 0.0559870567 0.0511898489 0.0433473344 0.0318412776 0.0168281658 0.0 F i n i t e Element Solution a 0.0041474265 0.0049651334 0.0072308490 0.0107261714 0.0147068992 0.0180846562 0.0187522740 0.0139112653 0.0 0.0600827801 0.0598033790 0.0585803199 0.0560003541 0.0512030557 0.0433671069 0.0318704518 0.0167978695 0.0 T h e solution becomes increasingly wavy i n the y—direction as Rw increases, and for this reason grid points lie closer to each other i n this direction than along the pipe axis. T h e aspect ratio is 6.4 which is much less t h a n 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), w i t h increasing frequency Reynolds number, the sinusoidal component increases whereas the cosinusoidal component decreases. A t Rw = 10, b o t h components are of the same order of magnitude, and at 7^, = 20 the sinusoidal component becomes the dominant one. F r o m the velocity field results, stream function values at grid points were calculated. 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 Cosinusoidal c o m p o n e n t — ~ Sinusoidal c o m p o n e n t 7 ^ = 10 •^ Q 0.00 0.15 0.05 Velocity profile Cosinusoidal c o m p o n e n t . ~ —> Sinusoidal c o m p o n e n t Ru = 20 O 0.00 0.03 Velocity 0.05 0.08 0.10 profile Figure 4.5: Variation of cosinusoidal and sinusoidal components for unsteady parallel flow for K = 1, 10 and 20 Chapter 4: Numerical Results 7T 37T Lt Li 45 (4.8) for t = 0, —, 7T and — . Twenty contours spaced equally between the instantaneous m i n i m u m and the instantaneous m a x i m u m 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. T h e instantaneous rate of flow, denoted by Q(t), is found b y integrating the it—velocity across the cross-section of the pipe as shown below and is p Q(t) = = since u ( y , t) = Table 4.5: difference y / r=1 J o u(y,t)-{2iry)dy 2 ^ ( 1 , t) dy Instantaneous rate of flow through a pipe under cosinusoidal pressure R = 1.0, po = 1.0 time 0 7T 2 7T 3TT 2 Frequency Reynolds Number iL20.0 Ru = 1.0 Rw = 10.0 0.38126 0.10945 0.04144 0.06357 0.16707 0.10706 -0.38126 -0.10945 -0.04144 -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: 4.3 49 Numerical Results Developing Flow in the Inlet Length of a Pipe 4.3.1 Theoretical Results A s a t h i r d example of axisymmetric flow, we consider the case of flow i n the inlet length of a straight circular pipe. A t a large distance f r o m the inlet the velocity distribution across the cross-section becomes parabolic. W e 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 . O w i n g to viscous friction, a boundary layer w i l l be formed on the w a l l and the depth of the boundary layer w i l l increase i n the downstream direction. A t the beginning, i.e., at small distances from the inlet section, the boundary layer w i l l grow i n the same way as it would along a flat plate at zero incidence [26]. T h e resulting velocity profile w i l l consist of an axisymmetric boundary layer near the w a l l joined i n 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 wall 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. T h i s process can be analyzed mathematically i n one of two ways. F i r s t , the i n tegration can be performed i n 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 distribution, i.e., integration can proceed i n the upstream direction. H a v i n g obtained both solutions, say i n the form of series expansions, we can retain a sufficient number of terms i n either of t h e m , joining the two solutions at a section where b o t h are still Chapter 4: Numerical 50 Results applicable. In this way the solution for the whole inlet is obtained. Figure 4.9: F l o w i n the inlet section of a P i p e A s shown i n figure (4.9), this problem has a singularity at x = 0, y = R, namely u(Q-,R) = u(0+,ir:) = U0 0. W e start b y satisfying the continuity equation: f (2Try)-u(x,y)dy = R J o U0 • (rf) or, fuydy JO = \VoR . 2 z N e x t we introduce a non-dimensional parameter ll_ 6 and expand U(x) ~\l _ x_ R ' e R = u ( x , 0 ) i n a series form and balance the m o m e n t u m equation. Further details are given i n 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 &Padd '= 0.04(2R)-Re = -0.3005i?e U s i n g 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, y 2 y Choosing ib = 0 along the centerline of the pipe gives C = 0 and W e also know that at the exit, the velocity is fully developed. u = A{1 — y 2) where A is a constant. choosing tp = —• Hence v = 0 and A g a i n integrating the velocity fields and 0 along the centerline, yjexn must remain constant at any section, yj'entrance = 2 = ~^y (^ tb trance(y en 2 ~ y )2 = 1) = Since the volume of flow i> it(y ex = 1)- T h i s gives A = 2, i.e., the centerline velocity is doubled. 4.3.2 Finite Element Results A s 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 g r i d of 20 x 10 elements shown i n figure (4.10) is used w i t h the length and w i d t h of successive elements chosen according to an arithmetic progression. A n o t h e r patch of 4 x 10 elements precedes the pipe entrance representing the undisturbed uniformly distributed flow of magnitude U 0 = 1.0. T h e finite element program is r u n for Re = 0, 1, 10, 50, 100, 150, 200, 250, 300, 350 a n d 400 w i t h a tolerance level of 1.0 x 1 0 ~ 5 . Results converge i n 4 Chapter 4: Numerical Results 52 to 8 iterations w i t h more iterations required as the Reynolds number is increased progressively. T h e 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 ) u(oo,0) < 0.001. T w e n t y streamline contours equally spaced between the m a x i m u m and the m i n i m u m values of the stream function are plotted i n figures (4.11), (4.12), (4.13) and (4.14). A g a i n the y—dimension has been enlarged four times to clearly reveal the contours. T w e n t y streamline contours equally spaced between the m a x i m u m and the m i n i m u m values of the stream function are plotted i n figures (4.11), (4.12), (4.13) and (4.14). A g a i n the y—dimension has been enlarged four times to clearly reveal the contours. T h e relationship between ZQ and RE is shown i n figure (4.15). Change i n u velocity along the centerline as the flow develops is presented i n figure (4.16) for different values of R.. 4.4 Flow Through an Axisymmetric Orifice A steady flow problem of interest to b o t h engineers and mathematicians is that of a viscous, incompressible fluid through an orifice. Greenspan [10] used a finite difference method w i t h a simple smoothing which yielded diagonally dominant systems of linear algebraic equations. Olson [18] used a pseudo-variational finite element theory w i t h 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 i n 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. T h e finite element grid used for a section bounded by the wall and the centerline is shown i n figure (4.18). T h e 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 = Number of elements = Number of equations for steady flow = x—axis scale: y—axis scale: 789 280 1665 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 R =50, 100 and 150 e Chapter 4: Numerical 56 Results RE- 200 RE= 300 Figure 4.13: Streamline contours for developing flow for R =200, 250 and 300 e Chapter 4: Numerical Results RE= 350 Figure 4.14: Streamline contours for developing flow for Re =350 and 400 Chapter 50 0 4: Numerical Results -i 0 1 0 0 2 0 0 3 0 0 4 0 0 - j — i — \ — i — | — i — i — i — i — i — i — i — | — i — i — i — | — i — i — i — | Reynold number Figure 4.15: — versus Re (R ) e 5 0 0 Chapter 4: Numerical Results 2.0 o 1.6 R R R R R H ' ^ > 1 .4 / / I I I / / / e e e e e = 0 (linear =10 -50 =100 -150 / I CD Co 1-2 1 .0 i—r 1 1 1 3 Distance r i 5 1 7 from 1 1 9 r ~> 11 1 ' 13 entrance Figure 4.16: Velocity distribution along the centerline for different Re 1 15 Chapter 4: Numerical Results 60 Figure 4.17: Flow through an axisymmetric orifice • u = 4(1 — y ), v = 0 on the upstream edge (Poiseuille conditions), 2 • 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\{\-y )-y 2 J o dy y 2(2-y ) 2 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 contour 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 R= b = 17.05 1.0 / = 0.5 8 = 0.5 Number of nodes Number of elements Number of equations for steady flow x-axis scale: y-axis scale: = = = 289 80 479 1 cm = 1.540 units 1 cm = 0.385 units Figure 4.18: Finite element grid for flow through an axisymmetric orifice Chapter 4: Numerical Results 62 RE= 0 RE= 1 1 \ =| \ r~ RE= 5 = Figure 4.19: Streamline contours for orifice flow for Re =0, 1 and 5 Chapter 4: Numerical Results 63 RE= — • — 10 —Cv RE= 15 5fc RE= 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. been enlarged four times i n the first two figures. A g a i n the y—dimension has 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 w i t h the results of the present study. It is seen that the contours obtained using the 8-noded isoparametric elements compare well w i t h 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. T h e small bubble upstream of the step observed b y 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. T h e 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 subsections. = 5, Rw = 10 and Rw = 20 are presented i n the subsequent Shear stress along the w a l l of the pipe is determined by the method outlined i n A p p e n d i x D . 4.5.1 Geometry of the Stenosis F i g u r e (4.22) represents a model of the stenosis. T h e 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 = U0(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 m a x i m u m 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 i n d of stenosis. Chakravarty [2,3], Forrester [8,9], Deshpande [6] and many others have chosen to m o d e l the constriction w i t h a half cosine wave, i.e., 1 + cos r(x) = < 7r(x — (a + /)) a<x< R (a + 21); otherwise. T h i s m o d e l , however, has a drawback v i z . the beginning and end of the stenosis are dr not smooth: that is, ^— has a j u m p at x = a and x = (a + 21). T o get around this difficulty, some researchers use a full cosine wave to model the stenosis. A n o t h e r 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 , /? a n d 7 , we use y = R at x = a a n d x = (a + 2/), and y = (R — 8) at x = (a + /). T h i s gives 1 (* -/ ) 2 2 r(x) = < (tS + / ) 2 2 12 2(5 2 - [x - (a + O f } a < x < ( a + 2Z); otherwise. T h i s m o d e l also suffers from the same drawback, although to a lesser extent. W e have adopted this m o d e l for the present study. T h e choice of 8 is governed b y the area reduction desired. N o r m a l l y one restricts 8 to a m a x i m u m of 0.5JR. W e chose 8 = 0.707i?; this results i n a 50% area reduction at the m a x i m u m 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. A s 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. T h e half-length of the stenosis, /, affects the rate of change of diameter. Thus a very small / means that the radius decreases, a n d subsequently increases, at a very fast rate a n d the resulting flow might be similar to the flow over an orifice. T h e grids shown i n figure (4.23) and figure (4.24) were used for this p r o b l e m w i t h the parameters as given i n the table (4.6) below. T h e number of elements across the cross-section remains constant; between a < x < (a + 2/), the w i d t h of each element T\ X) is reduced b y a factor of ———. T h e number of elements along the y—direction R , and i n the pre-stenosis zone are chosen the same for both the grids. Over the length of the stenosis a n d i n 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, b o t h 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 R = 1.0 b = 5.0 / = 3.0 8 = 0.707 Number of nodes Number of elements Number of equations for steady flow Number of equations for unsteady flow x—axis scale: y—axis scale: = = = = 381 112 772 2316 1 cm = 1.052 units 1 cm = 0.263 units Figure 4.23: Grid 1 used for analysis of stenotic flow Chapter 4: Numerical 69 Results Fine Grid for Pulsatile Flow (Grid 2) Number of nodes Number of elements Number of equations for steady flow x-axis scale: y-axis scale: = = = 719 216 1487 1 cm = 1.564 units 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 Designation N o . of nodes N o . of elements i n x—direction No. of elements i n y—direction Unconstricted radius R Pre-stenosis length a Post-stenosis length b Half-stenosis length / M a x i m u m constriction S y 2) and v = 0. Grid 1 Grid 2 short/coarse 381 14 8 . 1.0 5.0 5.0 3.0 0.707 long/fine 719 27 8 1.0 5.0 10.0 3.0 0.707 Results converged for Reynolds numbers up to Re = 500 for the coarse g r i d (grid 1). For Reynolds numbers higher than 500, grid 1 d i d not lead to converged solutions due to either, shorter post-stenosis length or due to inability of the finite elements to represent the rapid changes i n primitive variables because of the coarseness of the grid; perhaps a combination of b o t h the above. G r i d 2, however, yielded converged solutions up to very high Reynolds number and i n relatively less number of iterations. W e present these results for Reynolds numbers up to Re = 2000. W e 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). A s 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. T h i s , 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 i n the wake of the stenosis. F r o m the streamline contours for the coarse g r i d , 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 little zig-zagged which suggests that the results thus obtained are grid dependent. T h i s numerical i n accuracy results from the fact that i n the current formulation the continuity equation is not satisfied exactly, whereas the streamline formulation assumes it to be satisfied exactly. T h i s can be avoided by analyzing the flow w i t h a very fine grid. Streamline contours for selected values of Re are given i n 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 i n conjunction w i t h the scheme outlined i n A p p e n d i x C is used to plot the contours. T w e n t y streamline contours equally spaced between the m i n i m u m and the m a x i m u m 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 plotting for different values of Re are summarized i n table (4.7). Because of the no slip condition along the wall, the fluid experiences a high a shear stress there. T h e shear acting on the w a l l is equal and opposite of the shear acting on the fluid at y = r(x). C o m p u t a t i o n of the shear stress from the primitive variables solution is similar to the computation of vorticity £ and further details are given i n A p p e n d i x D. A t a fully developed, non-stenosed, cross-section the nondimensional shear stress equals 2. Over the length of the stenosis, w i t h the increase i n velocity gradient, there is a corresponding increase i n the shear stress. In the wake of the stenosis, as separation occurs, the velocity gradients, and consequently the shear stress changes sign. T h i s is shown i n figure (4.34) for the coarse grid and i n 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. T h i s suggests that grid 1 is sufficiently accurate up to that Reynolds number. Chapter 4: Numerical 72 Results RE- 0 Figure 4.25: Streamline contours for steady flow through stenosis for grid 1 for Re 1 and 10 Chapter 4: Numerical Results 73 RE= 50 RE= 100 RE= 150 Figure 4.26: Streamline contours for steady flow through stenosis for grid 1 for R = 50, 100 and 150 e Chapter 4: Numerical 74 Results RE- 200 Figure 4.27: Streamline contours for steady flow through stenosis for grid 1 for R = 200, 250 and 300 e Chapter 4: Numerical 75 Results RE= 0 RE= 10 RE= 50 Figure 4.28: Streamline contours for steady flow through stenosis for grid 2 for Rt = 0, 10 and 50 Chapter 4: Numerical 76 Results RE= 100 RE = 2 0 0 Figure 4.29: Streamline contours for steady flow through stenosis for grid 2 for R = 100, 150 and 200 e Chapter 4: Numerical Results 77 RE= 300 RE-= 500 Figure 4.30: Streamline contours for steady flow through stenosis for grid 2 for R = 300, 400 and 500 e Chapter 4: Numerical Results 78 RE= 600 RE= 700 RE= 800 Figure 4.31: Streamline contours for steady flow through stenosis for grid 2 for R = 600, 700 and 800 e Chapter 4: Numerical Results 79 RE= 900 RE= 1250 Figure 4.32: Streamline contours for steady flow through stenosis for grid 2 for R = 900, 1000 and 1250 e 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 — 2~I 0 i i i i i 4 Distance i i i 1 i 8 from i i i 1 2 i i i i 1 6 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 Grid 5 - / / \ i </ V fc> '»' » 1' 8- 2 '.'/Av:? / / g>6- \\! \\ - linear -- R =50 R =100 R =200 Re-300 e e e o \\ 4co CO «i) ^ •So CO S i \ \ £/ k1 ? ^ M i* \ f/ '-' "' \ \.\ »\\ i >'. V,\ 0\ 1 V \\ V, o- ^ CO -2- T — i — r 0 4 Distance T—i—i—|—r 8 12 from —i—i—i—i—| 16 20 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 8 2 " CO - 4 - CO " —S~\—i—i—i—i—i—i—i—i—i—i—i—i—i 0 4 Distance 8 12 from i i—i—i—i—i—i 16 20 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 a = 5.0 6 = 10.0 Reynolds Number Re 0 (linear) 50 100 • 200 300 400 500 600 700 800 900 1000 1250 1500 1750 2000 4.5.3 Grid 2 / = 3.0 R = 1.0 8 = 0.707 Maximum Stream Function Value 0.2500 0.2500 0.2500 0.2500 0.2503 0.2508 0.2513 0.2517 0.2522 0.2527 0.2531 0.2535 0.2543 0.2549 0.2554 0.2559 Location X y 5.04 5.04 5.04 0.99 0.99 0.99 0.99 0.93 0.90 0.89 0.88 0.89 0.88 0.88 0.87 0.86 0.86 0.85 0.85 12.18 11.13 11.13 11.13 11.13 12.18 12.18 12.18 12.18 12.18 12.18 12.18 12.39 Pulsatile Flow Through Stenosis A s pointed out i n section 3.4, the pulsatile nature of arterial fluid flow is approximated b y the first three components of the Fourier series. T h u s 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 i n this sub-section; all the solutions were obtained using grid 1 w i t h 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 i n 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 p r i m i t i v e variable solutions and a 100 x 100 rectangular grid i n conjunction w i t h the scheme outline i n A p p e n d i x C is used to plot the contours. T w e n t y streakline contours equally spaced between the m i n i m u m and the m a x i m u m values are plotted at t = 0, —, ir and — . C o m p u t a t i o n L* L* of the shear stress from the primitive variables solution is similar to the computation of vorticity £; further datails are given i n A p p e n d i x D . A l o n g the curved portion of the wall, axial stress <rx, radial stress o~y and the shear stress rxy are computed. To find the shear stress acting on the wall, r ^ , M o h r ' 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 i n figures (4.37-4.60) clearly reveal the effects of the stenosis. T h e instantaneous shear values along the w a l l of the pipe are shown i n figures (4.61-4.63) for t = 0, —, -K and 3TT —. Lt . . A t the beginning of the cycle, t = 0, the steady and the periodic components are additive and hence the shear stress and the stream function values along the w a l l are m a x i m u m at this moment. A t quarter cycle, t = —, only the steady and L\ the sinusoidal components are effective and hence the m a x i m u m shear stress and the m a x i m u m value of stream function along the w a l l drop; at half cycle, t = TT, the cosinusoidal component is substractive and hence the m a x i m u m shear stress and the m a x i m u m value of stream function along the w a l l is m i n i m u m at this instance. D u r i n g the second half of the cycle, this process is reversed and the m a x i m u m value is again obtained at t = 2ir. Due to the coarseness of the grid, however, the instantaneous shear stress values along the wall at various moments during a complete cycle showed a very small amount of oscillation near the m i n i m u m value. These oscillations, too, damp out at half cycle t = ir, when the cosinusoidal component is substractive. Chapter 4: Numerical 86 Results RE= 0 RW= 5 t = RE=0 RW=5 RE= 0 RW= 5 o t = 7T RE=0 RW= 5 Figure 4.37: Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 0 (linear) , for . t 7T 37T = 0, —, 7r and — Chapter 4: Numerical 87 Results RE = 1 RW= 5 t RE= 1 RW= 5 RE= 1 RW= 5 =o Figure 4.38: Streaklines for pulsatile flow through stenosis for Ru, = 5, R = 1 for e „ t 7T = 0, - , , 37T 7T and y Chapter 4: Numerical Results RE= 5 88 RW= 5 =o t RE= 5 RW= 5 7T t RE= 5 = RW= 5 t = RE= 5 2 7T RW= 5 Figure 4.39: Streaklines for pulsatile flow through stenosis for R\j = 5, Re = 5 for „ , 37T TC f = 0, -, 7T and — Chapter 4: Numerical 89 Results RE = 10 RW= 5 t RE= 10 RW = 5 RE= 10 RW= 5 = 0 t = X RE= 10 RW= 5 Figure 4.40: Streaklines for pulsatile flow through stenosis for R^ = 5, R = 10 for e t = 0, - , x and — Chapter 4: Numerical Results 90 RE= 20 RW= 5 •RE=20 RW= 5 Figure 4.41: Streaklines for pulsatile flow through stenosis for R^, = 5, R = 20 for e « = 0, - , 7r and — Chapter 4: Numerical Results RE= 50 91 RW= 5 RE= 50 RW= 5 RE= 50 RW= 5 RE= 50 t =o t = 7T RW= 5 Figure 4.42: Streaklines for pulsatile flow through stenosis for Ru,= 5, R = 50 for e ' 2' T T Chapter 4: Numerical Results RE- 100 RW= 5 t RE= 100 RW= 5 RE= 100 RW= 5 =o t = TT RE- 100 RW = 5 Figure 4.43: Streaklines for pulsatile flow through stenosis for Ru,= 5, R = 100 for e 7T t = 0, - , - 37T and y Chapter 4: Numerical 93 Results RE= 150 RW= 5 t = 0 RE= 150 RW= 5 7T <=2 RE= 150 RW= 5 t RE= 150 = 7T RW= 5 Figure 4.44: Streaklines for pulsatile flow through stenosis for R^ = 5, Re = 150 for „ 7r t = 0, - , TT 3TT and y Chapter 4: Numerical Results RE= 0 94 RW= 10 t RE= 0 RE = 0 =o RW= 10 RW= 10 t = RE= 0 7T 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 RE= 1 95 RW= 10 Figure 4.46: Streaklines for pulsatile flow through stenosis for R^, = 10, R = 1 for e < = 0, - , 7T and — Chapter 4: Numerical Results RE= 5 96 RW= 10 •i =0 RE- 5 RE= 5 RW= 10 RW= 10 ' Figure 4.47: Streaklines for pulsatile flow through stenosis for R\, = 10, Re = 5 for ' 2' * T Chapter 4: Numerical 97 Results RE= 10 RW= 10 t RE= 10 RW= 10 RE= 10 RW= 10 Figure 4.48: Streaklines for pulsatile flow through stenosis for t = 0, - , 7T and — =o = 10, R = 10 for e Chapter 4: Numerical Results RE= 20 98 RW= 10 t = RE= 20 RW= 10 RE= 20 RW= 10 RE= 20 RW= 10 o Figure 4.49: Streaklines for pulsatile flow through stenosis for R\j = 10, Re = 20 for t = 0, —, 2 7r and — 2 Chapter 4: Numerical Results RE = 50 . 99 RW= 10 t RE= 50 RW= 10 RE= 50 RW= 10 =o Chapter 4: Numerical Results RE= 100 100 RW= 10 t = RE= 100 RW= 10 RE = 100 RW= 10 t RE= 100 o = 7T RW= 10 Figure 4.51: Streaklines for pulsatile flow through stenosis for R^ — 10, Re = 100 for „ 7T t = 0, - , , 37T 7T and y Chapter 4: Numerical Results RE = 150 101 RW= 10 t RE= 150 RW= 10 RE= 150 RW= 10 RE= 150 RW= 10 =o Figure 4.52: Streaklines for pulsatile flow through stenosis for R^, = 10, Re = 150 for Chapter 4: Numerical Results RE= 0 102 RW= 20 t = REr 0 RW= 20 RE= 0 RW= 20 o 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 103 Results RE- 1 RW= 20 Chapter 4: Numerical 104 Results 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, R = 5 for e Chapter 4: Numerical 105 Results RE = 10 RW= 20 t RE= 10 RWr =o 20 Figure 4.56: Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 10 for t = 0, - , 7T and — Chapter 4: Numerical 106 Results RE- 20 RW= 20 Figure 4.57: Streaklines for pulsatile flow through stenosis for R^ = 20, Re = 20 for t = 0, - , JT and — Chapter 4: Numerical 107 Results RE= 50 RW= 20 t = RE=50 RW=20 RE= 50 RW= 20 0 t = TX RE= 50 RW= 20 Chapter 4: Numerical Results RE= 100 108 RW= 2 0 t RE= 100 RW= 2 0 RE- 100 RW= 2 0 =o t = IT RE= 100 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 109 RW= 20 t =0 RE= 150 RW= 20 RE= 150 RW= 20 t RE= 150 = 7T 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 t = entrance -8 I 0 i * 8 Distance from Distance from i 12 i , , i 16 entrance 7T 0) 40' 5 Co CO ca — — — — -- 24- linear R.= 10 R.-50 R.-100 R.-150 16- 5CO S» a CD CO 8- 0- -8' 0 4 Distance ' ' ' 8 from - 1 — 12 16 entrance 8 12 16 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 7T 0 t = 5 •• -8-|—i—i—i—r—i—'—i—i— — —i—i—'—i—<—i 0 4 6 12 16 1 Distance from 1 to —8 [ I I I — i i i i i—i—i i i i—i—i—i 0 48 12 16 Distance entrance t from = >7T lineor R.-10 R.=50 R.= 100 lineor R.= 10 - R.=50 R.= 100 R.= 150 6 Distance from entrance entrance * Distance 8 from 12 16 entrance Figure 4.62: Instantaneous shear stress along the wall of the pipe for grid 1 for Ru, = 10 and R = 0 (linear), 10, 50, 100 and 150 e Chapter 4: Numerical Results •8-T—' r—T 1 1 r—I 1 1 1—I 1 r—T 1—| 0 4 8 12 16 Distance from entrance 112 —8"| 1 1 I I 1 1 1—1 1 1 1 1 1 1 I | 0 4 8 12 16 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 CHAPTER 5 Conclusions 5.1 Concluding Remarks T h e present method of representing the complete solution as a sum of steady, cosinusoidal and sinusoidal components seems to work very well. T h e modified method of averaging used is also observed to be quite accurate. In particular, the circular arc approximation of stenoses works well. T h e overall agreement between the present study and other published results i n the literature is good. T h e 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. T h e results obtained for the fully developed flow under the influence of a cosinusoidal pressure difference agree well w i t h the numerical and graphical values reported by Schlicting. T h e streakline plots are almost completely horizontal indicating that even though the solution becomes increasingly wavy w i t h 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. T h e effect of Reynolds number on the development length was i n good agreement w i t h the theoretical values reported i n [26]. T h e additional pressure drop, Ap dd a 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 p r i m i t i v e 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 m i n i m u m section area of cross section of the pipe was reduced b y 75% and the u—velocity was quadrupled. T h e stream function contours for Re = 10, plotted to 1:1 scale, compares well w i t h 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 i n 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 w a l l , i n c l u d i n g the m a x i m u m and the m i n i m u m values, m a x i m u m stream function values, etc., show a monotonic improvement w i t h the finer grid. T h e separation and the re-attachment points as obtained from the stream function results and from the shear stress results are i n total agreement w i t h 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 function contours clearly reveal the effects of the stenosis and instantaneous shear values along the w a l l of the pipe indicate the initiation 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 wall at various moments during a complete cycle showed a very small amount of oscillation near the m i n i m u m value. Chapter 5.2 5: Conclusions 115 Suggestions for Further Developments Some specific recommendations for further studies based o n the work i n this thesis are: • perform the analysis w i t h the transformations V = f( ,y) x to m a p the curved boundary onto a straight line £ = £o, a n d n = tanh(i^x) to m a p the infinite domain to the finite domain. • perform the analysis for multiple stenoses i n 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 a n d compare a n d contrast it w i t h the present study. Bibliography [1] C a i l k , D . A . , and N a g h d i , P . M . , Axisymmetric side a Slender Surface of Revolution, Motion of a Viscous Fluid In- Transactions of A S M E : Journal of A p p l i e d 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 M i s r a , J . C . , Flow in Arteries in the Presence of Stenosis, J o u r n a l of Biomechanics, V o l . 19, 1986, pp. 907-918. [4] Daly, B a r t J . , A Numerical Femoral Arteries, Study of Pulsatile Flow Through Stenosed Canine 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 & F l u i d s , V o l . 15, 1987, pp. 275-280. [6] Deshpande, Through M . D . , Giddens, D . P . , and M a b o n , R . F . , Steady Modeled Vascular Stenosis, Laminar Journal of Biomechanics, V o l . 9, Flow 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, A t l a n t a , G A 30332, U S A , 1977. [8] Forrester, J . H . , and Young, D . F . , Flow Through a Converging-Diverging and its Implications in Occlusive Vascular Disease—I: Theoretical J o u r n a l of Biomechanics, V o l . 3, 1970, pp. 297-305. 116 Tube Development, Bibliography 117 [9] Forrester, J . H . , and Young, D . F . , Flow Through a Converging-Diverging and its Implications imental Results in Occlusive Vascular and Their Implications, Disease—77; Theoretical Tube and Exper- Journal of Biomechanics, V o l . 3, 1970, pp. 307-316. [10] Greenspan, D . , Numerical an Orifice for Arbitrary Studies Reynolds of Viscous, Incompressible Flow Through Number, International Journal for Numerical M e t h o d s 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 & F l u i d s , V o l . 13, 1985, pp. 239-269. [12] Lohner, R . , Finite Elements in CFD: What Lies Ahead, International Journal for N u m e r i c a l Methods i n Engineering, V o l . 24, 1987, pp. 1741-1756. [13] M a c D o n a l d , D . A . , Pulsatile Flow in Catheterized Artery, Journal of Biomechan- ics, V o l . 19, 1986, pp. 239-249. [14] M c L a c h a l a n , N . W . , Ordinary and Physical Non-linear Differential Equations in Engineering Sciences, Oxford University Press, 1950. [15] M o o r t y , S., A Parametric Study of Rigid Body-Viscous Flow Interaction, M.A.Sc. Thesis, T h e University of B r i t i s h C o l u m b i a , Vancouver, C a n a d a , 1987. [16] M o o r t y , S., A Numerical Study of Low Reynolds Number Fluid-Structure Inter- action, J o u r n a l of F l u i d s and Structures, V o l . 3, 1989, pp. 37-60. [17] Nayfeh, A . H . , and M o o k , D . T . , Nonlinear Oscillations, W i l e y International, New Y o r k , 1979. [18] Olson, M . D . , Variational-Finite Element Methods for Two-Dimensional and A x isymmetric Navier-Stokes Equations, Finite Elements in Fluids, 72, editors Gallagher, R . H . , et al., J o h n Wiley, L o n d o n , 1975. V o l . 1, pp. 57- 118 Bibliography [19] P a t t a n i , P . G . , Nonlinear Analysis of Rigid Body-Viscous Flow Interaction, Thesis, T h e University of B r i t i s h C o l u m b i a , Vancouver, Canada, [20] P a t t a n i , P . G . , and Olson, M . D . , Periodic Interaction, Solution Ph. D. 1986. of Rigid Body-Viscous Flow International Journal for Numerical Methods i n F l u i d s , V o l . 7, 1987, pp. 653-695. [21] P h i l p o t , E . , Yoganathan, A . P . , Sung, H . - W . , W o o , Y . - R . , Franch, R . H . , Sahn, D . J . , and V a l d e z - C r u z , L . , In-Vitro Pulmonary Artery Pulsatile Flow Visualization Studies in a Model, Transactions of A S M E : Journal of Biomechanical E n - gineering, V o l . 107, 1985, pp. 368-375. [22] Porenta, G . , Y o u n g D . F . , and Rogge, T . R . , A Finite-Element Flow in Arteries Including Taper, Branches Model of Blood and Obstructions, Transactions of A S M E : Journal of Biomechanical Engineering, V o l . 108, 1986, pp. 161-167. [23] P o z r i k i d i s , C , A Study of Peristaltic Flow, Journal of F l u i d Mechanics, V o l . 180, 1987, pp. 515-527. [24] R i c e , J . G . , and Schnipke, R . J . , An Equal That Does Not Exhibit Spurious Order Velocity Pressure Formulation Pressure Modes, Computer Methods i n A p p l i e d Mechanics and Engineering, V o l . 58, 1986, pp. 135-149. [25] R o b i c h a u d , M . P . , and Tanguy, P . A . , Finite Fluid Problems Element Solution by a Preconditioned of Three-Dim- ensional Incompressible Conjugate Residual Method, International Journal for Numerical Methods i n Engineering, V o 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 Through Planar and Axisymmetric Element Abrupt Analysis Expansions, of Laminar Flows Computers & F l u i d s , 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 Experimental 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 Methods for the Navier-Stokes Equations, Element Solution 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 with Inter-Element Penalties, Solution of Navier-Stokes Problems 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 Model with Varying Degrees of Pulmonic in a Pulmonary Artery Stenosis, Journal of Biomechanics, Vol. 19, 1986, pp. 129-146. [35] Young, D.F., Fluid Mechanics of Arterial Engineering, Vol. 101, 1979, pp. 157-175. Stenoses, Journal of Biomechanical 120 Bibliography [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 tions with Finite of Navier-Stokes Equa- Elements, Computer Methods in Applied Mechanics and Engi- neering, Vol. 57, 1986, pp. 223-237. APPENDIX A Newton-Raphson Procedure T h e modified Newton-Raphson iteration procedure adopted i n this study is outlined here. T h e 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 i n the form /i(xi,x2, ...,xn) = 0 f2(xx,x2, ... ,xn) = 0 fn[xX,X2, . . . ,Xn) (A.l) = 0. E x p a n d i n g 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) + Axi (A.2) where the derivatives are evaluated at ( x i 0 , x 2 o, • • • , xn0) and A x i = Xi — Xio A x 2 = x2 Ax n — X 2 Q (A.3) — Xji XJIQ. Similarly, carrying out the expansion for a l l the equations i n ( A . l ) and rearranging, we obtain 121 Appendix A: Newton-Raphson rfdfA fdfA \dx ) \dxj 2 0 fdfA (dfA \dx ) 0 2 (dfA \dxj \dx ) 2 fdfA \dx ) n Q Q ' Axi ' (dfA \dx ) 0 fdfA 0 _ Q VSa J 122 Procedure n Ax2 0 _ (dfA "* \dx ) n ' < . Axn , —fi(xio,x2o,.. -f2(xw,X20, > == • • • , Xno) < • —fn( x10, #20) • •• > (A.4) > ^no) > Q Equations (A.4) can be written i n the form [T]{Ax} = {—/} where the matrices denote the corresponding terms i n equation ( A . 4 ) . T h e algorithm for obtaining the solution vector {x} is outlined below: 1. Guess the solution vector {x}0. 2. C o m p u t e [T]. 3. C o m p u t e {—/}. 4. O b t a i n {Ax} f r o m [T]{Ax} = {-/}. 5. O b t a i n {x}x = {x}0 + {Ax}. 6. If {x}i is close to {x}0 w i t h i n the specified tolerance go to 8, otherwise set {x}0 = {x}x. 7. If [T] needs to be modified i n this iteration go to 2, otherwise go to 3. 8. Set {x} = {x}i and stop. APPENDIX B Details of Tangent Stiffness Matrix and Nonlinear Load Vector The left hand side of equation (3.12) can be re-written as [T] = K 0 0 K -M .0 0" M + 'A B C" V A 0 Z 0 A. K. where K , and M are the stiffness and mass matrices. Sub-matrices A to £ , using index notation, are given as A u 4- fi x A u 4- FH A v A = Re SfmjAm ' £x ijm B = f) y Do m i ' Re fix imj D U I m. ' cy r>v imi m 8 X- B v 2 • sx rvu. I cx fiu + Au &\mjA vm 0' + 6ij m An Bu 8 X B u + 8 y B v +6?- B v i sy /~iv Re 123 0 8ijm C771 &imjCm 2 0 + 8fmjCm + 8\jmC vm 0 Appendix B: Details Cx V = £ = of Tangent Stiffness Matrix DU SX i cx RK ijm JDU I CI) and Nonlinear Dt) U TDV ^imj m B m + Load Vector ijm- m D fiimjBm + 0~ ii B V V m m 0 RE rA A T h e incremental solution vector is {Ax} = < A B and the load vector is AC K 0 0 (A) 0 K M B 0 -M K C {-f} = \ 6? {A]At + B?B£/2 8f (A]Al + BJBJI/2 + C ? C £ / 2 ) + ^ ( A J A J + B ? B £ / 2 + C*Cj;/2) jk jk + C?C2/2) + 5? ifc (A?Aj£ + B V B £ / 2 + CvC^/2) ) 0 ^•*(AJBJ + AJfBJ) + S* (A?Bl + AJJBy) ik -Re ( Sf (A]Bt jk + AlBJ) + Sf (A]Bl + AJBJ) jk 0 6f (A]Ct + AtC]) + 8y {A)Cl + A j q ) jk jk ^• f c (AJCj; + AZC?) + ^ f c ( A V C J + Aj£CJ) 0 APPENDIX C Determination of the Stream Function Values Consider a n eight-noded isoparametric element as shown i n the figure (3.1). The shape functions iV,- and M,- for the element are presented i n Chapter 3. T h e value of the stream function $ and the coordinates (x, y) at a point can be obtained i n terms of the n o d a l values \ 1 / , - and nodal coordinates (x,-,y,-). *= * = 1, 2,..., 8 x = NiXi V = Nan T h e following procedure is adopted to obtain the value of accurately at each node of the plotting grid which lies w i t h i n the finite element. 1. Assume a quadratic interpolation for £ and n of the form £ = Ai + A2x + A3y + A4x 2 + A5xy + A6y 2 + ATx 2y + A8y 2x 7] = B1 + B2x + B3y + B4x 2 + B5xy + B6y 2 + B7x 2y + B8y 2x (C.l) T h e (x, y) and (£, n) values at each of the eight nodes on the element are known and hence eight equations each i n A 1 , A 2 , . . . , A 8 a n d Bx, B2,..., tained and solved for. Thus for any point (x0,y0) equation ( C . l ) , corresponding (^o?7?©) B% are ob- i n the finite element, using coordinates can be calculated. 2. Determine the m a x i m u m and m i n i m u m value of x and y coordinates of the finite element. T h e approximate area occupied b y this element is a rectangle of 125 Appendix C: Determination of the Stream Function Values length xmax - xmin along the x - a x i s and height ymax - ymin 126 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 w i t h coordinates ( x p , y ) p w i t h i n the finite element as outlined below: (a) Take the i n i t i a l guess value of (£P,77P) to be that obtained from equation ( C . l ) , i.e., iv = M + A2xp + A3yp + A j x , + A5xpyp + A6y 2p + A7x\yp + Asy 2pxp np = B1 + B2xp + B 3y p + B4xl + B&y 2v + B7x 2pyp + B8y 2xp 2 + B5xpyp (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. dNj J = dNj ds dNj ds (d) Determine the (a/, y') coordinates of the point using the calculated value of shape functions. .7 = 1,2,..., 8 (e) Determine the error i n the (x, y) coordinates. Ax = xp — x' Ay = y - y' p Determine the error i n the (£,77) coordinates. Appendix C: Determination of the Stream Function Vaiues 127 (g) T h e new value of (£, n) is given as ? = V' = + V P + ^ A??. W i t h these new values of (£, 77), repeat steps (b) to (f) u n t i l desired convergence is achieved. 4. Once convergence has occurred i n step 3, the stream function value can be obtained using \I> = N^i. APPENDIX D Determination of Stress Components For axisymmetric flow, the three non-trivial components of the stress tensor i n their non-dimensional form are given as n 0-y 'xy du = dv -P + 2 dy _ (dv du s \dx dy. (D.l) Consider an eight-noded isoparametric element as shown i n the figure (3.1). Using the shape functions AT,- and M t - presented i n Chapter 3, equations ( D . l ) are discretized as -PjMj + -PjMj + dNi dx 2Ui 2vi + Ui ' dNj dx dNi dy * = 1,2,...,8; ; = 1,2,...,4. (D.2) dNi' dy Stress components crx, ery and rxy at some point (x, y) G Q can be obtained i n terms of the n o d a l values of the p r i m i t i v e variables (u,v,p). T h e 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 w i t h i n which point (x0, yQ) lies x . 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. 1 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 rj = B1 + B2x + B3y + B4x 2 + Bhxy + B6y 2 + B7x 2y + B&y 2x (D.3) T h e (x, y) and (£, r?) values at each of the eight nodes on the element are known and hence eight equations each i n A\, A2,..., As and B\, B2,..., B$ are obtained and solved for. 3. Determine (£, n) value of each point of the plotting g r i d w i t h coordinates ( x 0 , y0) w i t h i n the finite element as outlined below: (a) Take the i n i t i a l guess value of ( £ 0 , n0) to be that obtained from equation (D.3), i.e., to = Ax + A2x0 + A3y0 + A4x 20 + A5x0y0 r]0 = Bi + B2x0 + B3y0 + B4x\ + B5x0y0 + A6y 20 + A7x 20y0 + + B6yl + B y x dNi (b) Calculate the shape function derivatives - — - and o£ for these values of £ 0 , 77 using the equations given 0 + B7x 2Qy0 A8y 2x0 8 2 0 d N -r-^- for i = 1 , 2 , . . . ,8 or) i n Chapter 3. (c) Calculate the Jacobian J. dNi dN r ds dN 1 ds X% dt dN dt Vi n Xi V i l (d) Determine the (x', y') coordinates of the point using the calculated value of shape functions. x' = NiXi y' = i = l,2,...,8 Niyi (e) Determine the error i n the (x,y) coordinates. Ax = x — x' 0 Ay = y0- y' Appendix D: Determination of Stress Components 130 (f) Determine the error i n the (£,?/) coordinates. J }^J J Ay An (g) T h e 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 ox dNi oy . 4. Once convergence has occurred m step 3, compute —— and —— using K dN dNi N dx dN 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 w a l l Appendix D: Determination of Stress Components 131 To obtain stresses i n any direction other than the global axes, we use the Mohr's circle construction. A s shown i n the figure ( D . l ) , for equilibrium crx, o~y and rxy on unit surface area must balance and crxi [crx + o~y) (crx I cr i x = <?y< = T x , y l = yv 2 ((Tx + 0-y) 1 —— - _(^- ( 7 i/) acting on y/2 surface area. Thus Txiyi - (j ) — Vy 2— — x y (CTX — (Jy) . cos 20 + rxy s i n 29 — — cos 29 + rxy s i n 29 s i n 2 g + acting T ^ c o s 2 9
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Finite element analysis of periodic viscous flow in...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Finite element analysis of periodic viscous flow in a constricted pipe Singh, Rajesh Kumar 1989
pdf
Page Metadata
Item Metadata
Title | Finite element analysis of periodic viscous flow in a constricted pipe |
Creator |
Singh, Rajesh Kumar |
Publisher | University of British Columbia |
Date Issued | 1989 |
Description | The finite element method is extended to analyze axisymmetric periodic viscous flow in 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 wall of the artery. The artery geometry is approximated as an axisymmetric channel. Blood 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 primitive variable representation of the Navier-Stokes equations. Eight-noded isoparametric elements with quadratic interpolation for velocities and bilinear for pressure 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 R[sub w], and the steady Reynolds number R[sub e]. The pulsatile flow is modeled as a linear combination of steady, cosinusoidal and sinusoidal components. The magnitude of each component is determined by minimizing the error in approximation through Galerkin's procedure. Results are obtained for various values of R[sub w] and R[sub e] for laminar flow. Results are also presented for limiting cases and, wherever possible, numerical results thus obtained are compared with analytical and experimental results published in the literature. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-08-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0062609 |
URI | http://hdl.handle.net/2429/28063 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
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
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062609/manifest