I N C O M P R E S S I B L E V I S C O U S F L O W I N T U B E S W I T H O C C L U S I O N S By Huaxiong Huang B.Sc. (Mathematics) , Fudan University Shanghai, People's Repubhc of C h i n a , 1984 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S I N S T I T U T E O F A P P L I E D M A T H E M A T I C S D E P A R T M E N T O F M A T H E M A T I C S We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A October 1991 (c) Huaxiong Huang, 1991 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Date npr m\ DE-6 (2/88) A B S T R A C T Viscous, incompressible flow in tubes with partial occlusion is investigated using nu-merical and experimental procedures. The study is related to the problem of atheroscle-rosis, one of the most common diseases of the circulatory system. One of the computational difficulties in solving the incompressible Navier-Stokes equa-tions is the lack of pressure or vorticity boundary conditions. A finite difference approach, referred to as the interior constraint (IC) method, is proposed to resolve this difficulty. As a general numerical method, it is formulated for both the stream function-vorticity and primit ive (physical) variable formulations. The procedure is explained using one-dimensional model with extensive numerical tests presented for two-dimensional cases including flow in a driven cavity, and flow over a backward facing step. Results are obtained with second-order accuracy. Next, the IC method is applied to flow in a tube with an occlusion, which is used as the model for blood flow in stenosed arteries in the study of the pathology of atherosclerosis. Numerical results are obtained for both steady and pulsatile flows. Results are compared with those of S I M P L E , one of the commercially available numerical algorithms. The pulsatile flow study revealed several interesting new features. It suggests that the high shear stress is not likely to initiate atherosclerosis lesions. The recirculation region, which is a prominent feature of the unsteady flow, is more likely to cause the init iat ion and development of the disease. Exper imental measurements for steady flow complement the numerical study and show qualitative agreement. Table of Contents A B S T R A C T ii List of Tables viii List of Figures x N O M E N C L A T U R E xv M E D I C A L PHRASES xvii A C K N O W L E D G E M E N T S xix 1 I N T R O D U C T I O N 1 1.1 Prel iminary Remarks 1 1.2 F l u i d Dynamics in Stenosed Arteries 4 1.3 Flow Model 6 1.4 Computational Methods and Limitations 7 1.5 Scope of Present Investigation 12 2 I N I T I A L A N D T W O - P O I N T B O U N D A R Y - V A L U E P R O B L E M S 15 2.1 First -Order System 15 2.1.1 Initial-value problems (IVPs) 15 2.1.2 Two-point boundary-value problem for a two-equation system . . 17 2.1.3- One-dimensional Navier-Stokes type equations 18 2.1.4 One-dimensional Poisson-type approach 19 2.2 Second-Order System 21 2.2.1 Woods method 22 2.2.2 Integral constraint method 23 2.2.3 Influence matrix method 23 2.2.4 Interior constraint method 24 3 IC M E T H O D F O R T W O - D I M E N S I O N A L F L O W P R O B L E M S 28 3.1 Flow in Driven Cavity 29 3.1.1 Arti f ic ial compressibility approach 29 3.1.2 Pressure Poisson equation approach 34 3.1.3 Stream function-vorticity formulation 39 3.2 Flow over a Backward Facing Step 41 4 S O L U T I O N P R O C E D U R E A N D N U M E R I C A L R E S U L T S 47 4.1 Solution Procedure for Two-Dimensional Problems 47 4.1.1 Stream function-vorticity formulation 47 4.1.2 Pr imit ive variable formulation 48 4.2 Test on One-Dimensional Model Problems 49 4.2.1 One-dimensional initial-value problem 49 4.2.2 One-dimensional stream function-vorticity model 51 4.3 Test on Two-Dimensional Flow Problems 54 4.3.1 The driven cavity flow 54 4.3.2 Flow over a backward facing step 65 5 N U M E R I C A L S I M U L A T I O N O F F L O W I N S T E N O S E D A R T E R I E S 72 5.1 Stream Function-Vorticity Formula 72 5.1.1 Governing equations 72 5.1.2 Boundary conditions 73 5.1.3 Computational domain 74 5.1.4 Representation of constriction 75 5.1.5 Equations and boundary conditions in transform domain 76 5.1.6 Discretization 77 5.2 Pr imit ive Variable Formulation 80 5.2.1 Governing equations 80 5.2.2 Boundary conditions 81 5.2.3 Computational domain 83 5.2.4 Governing equations and boundary conditions in transform domain 83 5.2.5 Discretization 84 5.3 S I M P L E Method 87 5.3.1 Governing equations 87 5.3.2 Representation of constriction 87 5.3.3 Boundary condition and discretization . 88 5.4 Unsteady Flow Calculations 91 5.4.1 Governing equations and boundary conditions 91 5.4.2 T ime discretization 91 5.5 Prel iminary Computation 92 5.5.1 Convergence rate analysis 92 5.5.2 Location of the inlet and outlet boundaries 92 5.5.3 Effect of the grid size 93 6 E X P E R I M E N T A L I N V E S T I G A T I O N 97 6.1 Stenosed Model 98 6.2 Test Facihty 98 6.3 Flow Measurement System - Instrumentation 101 6.3.1 Measurement of velocity 101 6.3.2 Traversing mechanism 104 6.3.3 Pressure transducer 105 6.3.4 Flowmeter 106 6.4 Test Procedure 106 6.5 D a t a Processor 108 7 R E S U L T S F O R S T E A D Y S T A T E F L O W 115 7.1 Comment on Data Presentation 115 7.2 Effects of F low Parameters on Flow Fie ld . 116 7.2.1 Streamline and vorticity distribution 117 7.2.2 Separation and reattachment points 129 7.2.3 Pressure drop across stenosis 131 7.2.4 W a l l vorticity 136 7.2.5 Flow velocity 139 7.3 Comparison with different approaches 140 7.3.1 Pressure drop 140 7.3.2 W a l l vorticity 145 7.3.3 Locations of separation and reattachment points 146 7.4 Experimental Results 149 7.4.1 Pressure distribution 149 7.4.2 Velocity profile 150 7.5 Discussion 159 8 R E S U L T S O F P U L S A T I L E F L O W S I M U L A T I O N 165 8.1 Streamline and Vortic ity 168 8.2 Separation and Reattachment Points 187 8.3 W a l l Vort ic i ty Distribution 188 8 4 Discussion 195 9 C O N C L U S I O N 207 Bibliography 211 A F O R M U L A T I O N O F S I M P L E M E T H O D 219 B A D I M E T H O D A N D L I N E G A U S S - S E I D E L I T E R A T I O N 222 B . l Steady-State Solution • • 222 B.2 Unsteady Solution • • • 224 C W A L L S H E A R S T R E S S A N D V O R T I C I T Y 225 List of Tables 4.1 Convergence rate of I V P 53 4.2 Comparison of solution of I V P 53 4.3 Convergence rate of the 1-d second-order system 53 4.4 Comparison of stream function and vorticity at the center of primary vor-tex for different Reynolds numbers 58 4.5 Comparison of u-velocity along the vertical line through the cavity center. 59 4.6 Comparison of u-velocity along the horizontal line through the cavity center. 59 4.7 Comparison of vorticity along the moving boundary 60 4.8 Convergence rates of solution using the IC method for the driven cavity flow 60 4.9 Convergence rates of sohition using the IC method for the driven cavity flow 60 4.10 Comparison of the stream function and vorticity at the center of the pr i -mary vortex for the driven cavity flow 61 4.11 Comparison of the stream function and vorticity at the center of the pr i -mary vortex for the driven cavity flow 61 4.12 Comparison of the stream function and vorticity at the center of the pr i -mary vortex for the driven cavity flow 62 4.13 Results on separation and reattachment points for the flow over a backward facing step 67 vni 4.14 Characteristics of the primary recirculation region for the flow over a back-ward facing step 67 4.15 Comparison of two algorithms for the driven cavity flow (Re=100, grid size 61 X 61) 68 5.1 Source terms of the stream function-vorticity formulation 94 5.2 Source terms of the primitive variables formulation 94 5.3 Source terms of the conservation form for the S I M P L E method 94 5.4 Convergence rate of the steady state flow 95 5.5 G r i d dependence test on peak vorticity and separation length 95 5.6 Convergence rate of the unsteady flow 95 5.7 Test on entrance and outflow lengths 96 5.8 G r i d dependence test on peak vorticity and separation length 96 7.1 Geometry of stenosis models used in computation 163 7.2 Separation Reynolds numbers for five stenosis models 163 7.3 Comparison of peak wall vorticity for model M 2 163 7.4 Physiological quantities in the animal experiments 164 7.5 Peak wall shear stress {dyn-cm~^) corresponding to the animal experiments. 164 7.6 Calculated pressure drop (mmHg) corresponding to the animal experiments. 164 List of Figures 1.1 A diagrammatic representation of the major branches of the arterial tree.. (After Pedley, T . J . [37]) 14 2.1 Illustration of the 1-d computational domain 27 3.1 (a) Geometry of driven cavity; (b) Illustration of computational domain. 45 3.2 Schematic diagram of flow over a backward facing step 46 3.3 Computational domain for flow over a backward facing step 46 4.1 Stream function contours for the driven cavity flow. . 63 4.2 Vort ic i ty contours for the driven cavity flow 64 4.3 Stream function contours for flow over a backward facing step 69 4.4 Vort ic i ty contours for flow over a backward facing step 70 4.5 Vort ic i ty at the lower wall for flow over a backward facing step 71 4.6 Vort ic i ty at the upper wall for flow over a backward facing step 71 5.1 Geometry of the flow in a tube with a constriction and its boundary con-ditions using the stream function-vorticity formulation 79 5.2 niustration of computational domain using the stream function-vorticity formulation 79 5.3 Geometry of the flow in a tube with a constriction and its boundary con-ditions using the primitive variable formulation 86 5.4 l l lutration of the computational domain using the primitive variable for-mulation 86 5.5 Illustration of the artificial viscosity method using the S I M P L E algorithm. 90 6.1 Stenosis model used in experimental measurements 110 6.2 A schematic diagram of the water tunnel used in experimental study. . . I l l 6.3 L D A system with the traversing gear 112 6.4 Cahbration curve for the Barocel pressure transducer 113 6.5 Cal ibration curve for the orifice flowmeter 113 6.6 A schematic diagram showing instrumentation, data acquisition, and pro-cessing systems 114 7.1 Streamline and vorticity contours for model M l 120 7.2 Streamline and vorticity contours for model M 2 122 7.3- Streamline and vorticity contours for model M 3 124 7.4 Streamline and vorticity contours for model M 4 126 7.5 StreamHne and vorticity contours for model M 5 128 7.6 Separation and reattachment points v.s. Re 130 7.7 Comparison of pressure at the wall and axis for M 3 at Re = 500 132 7.8 Pressure distribution along the wall 134 7.9 Pressure drop across the stenosis 135 7.10 Vort ic i ty distribution along the wall for model M 2 137 7.11 Vort ic i ty distribution along the wafl at Re=100 137 7.12 Peak wall vorticity as a function of Re 138 7.13 Velocity profiles for model M l at Re = 500 141 7.14 Velocity profiles for model M 3 at Re = 500 142 7.15 Velocity profiles for model M5 at Re = 500 143 7.16 Centerhne velocity at Re = 500 144 7.17 Centerhne velocity for model M l 144 7.18 Comparison of pressure drop across the stenosis as obtained by several investigations 147 7.19 Comparison of peak wall vorticity for model M 6 148 7.20 A comparison of separation and reattachment points as predicted by sev-eral investigators 148 7.21 A comparison between numerically and experimentally obtained pressure distribution along the wall 153 7.22 Variat ion of the axial velocity profiles: a comparison between numerical and experimental results 155 7.23 Variat ion of the axial turbulent intensity as affected by the Reynolds number. 157 7.24 Comparison of the present experimental data for axial velocity profile with - those by Ahmed and Giddens [35] 158 8.1 Comparison of wall vorticity distribution 167 8.2 Stream function (top) and vorticity contour (bottom) of steady flow for MO at Re = 100 171 8.3 Stream function contour in a cycle for MO at iîe = 100, Rw = 10,e = l . . 173 8.4 Vortic ity contour in a cycle, MO, iîe = 100, -Ru; = 10, e = 1 175 8.5 Stream function contour in a cycle for MO at Re — 100, Rw = 10, e = 0.1. 177 8.6 Stream function contour in a cycle for MO hi Re = 100, Rw = 3.34, e = 1. 179 8.7 Stream function contour in a cycle for MO at Re = 600, Rw = 10, e = 1. . 181 8.8 Stream function contour in a cycle model M 5 at Re = 100, Rw = 10, e = 1.186 8.9 T ime history of the separation and reattachment points for MO at Re = 100,/îu; = 10, e = 0.1 191 8.10 Time history of the separation and reattachment points for MO at Re = 100, Rw = 10, e = 0.5 192 8.11 T ime history of the separation and reattachment points for MO at Re = 100, Rw = \{),e=l 192 8.12 T ime history of the separation and reattachment points for MO at Re = 100, Rw = 10, e = 1.5 193 8.13 T ime history of the separation and reattachment points for MO at Re = 100, = 7.5, e = 1 193 8.14 T ime history of the separation and reattachment points for MO at Re = 100, = 3.34, e = 1 194 8.15 T ime history of the separation and reattachment points for M 5 at Re = 100, itlu; = 10, e = 1 194 8.16 T ime history of the peak wall vorticity for MO at Re = 100, Rw = 10, e = 1.198 8.17 Effect of e on time history of the peak wall vorticity for MO at Re = 100, Rw = lO 199 8.18 Effect of Rw on time history of the peak wall vorticity for MO at Re = 100, e = 1. . 199 8.19 Effect of Re on time history of the peak wall vorticity for MO at Rw = 10, e = 1 200 8.20 Effect of H on time history of the peak wall vorticity at Re = 100, Rw = 10, e = l 200 8.21 Vort ic i ty distribution along the wall for MO at Re = 100, Rw = 10, e = 0.1.201 8.22 Vort ic i ty distribution along the wall for MO at Re = 100, Rw = 10, e = 1. 201 8.23 Vort ic i ty distribution along the wall for MO at Re = 100, Rw = 10, e = 1.5.202 8.24 Vortic ity distribution along the wall for MO at Re = 100, R.w = 3.34, e = 1.202 8.25 Vort ic i ty distribution along the wall for M 5 at Re = 100, Rw = 10, e = I. 203 8.26 Vort ic i ty distribution along the wall for MO at i?e = 600, /Zu; = 10, e = 1. 203 x m 8.27 Effect of e on vorticity distribution along the wall for MO at Re = 100, Rw = lO 204 8.28 Effect of Rw on vorticity distribution along the wall for MO at iîe = 100, e = 1 204 8.29 Effect of H on vorticity distribution along the wal l at i2e = 100, Rw = 10, € = 1 205 8.30 Effect of Re on vorticity distribution along the wall for MO at Rw = 10, e = l 205 8.31 Effect of grid size on vorticity distribution along the wall for M 5 at Re = 100, i?u; = 10, e = 1 206 A . l Control volume and notation for the S I M P L E method 221 N O M E N C L A T U R E e Ratio of unsteady and steady flux ( Vort ic i ty 77,^ Variables in curvilinear coordinate system A Wavelength of laser beam 1/ Kinematic viscosity p F l u i d density r Shear stress n Physical domain Çî'^ Computational domain dQ, Physical boundary dCl'^ Computational boundary (jj Frequency tp Stream function D Derivative operator V Divergence of velocity Dh Diff'erence operator Do Centered difference operator Forward difference operator D- Backward difference operator Difference operator in ^ direction Difference operator in r] direction D-^ Difference operator in x direction Dl^ Difference operator in y direction F Flow rate / Focus length of lens Doppler frequency H Stenosis height h Mesh size Hs Channel height K Height of a backward facing step in a channel L Stenosis length Li Separation length P Peclet number P Static pressure R Tube radius Re - Reynolds number Rw Frequency Reynolds number S^p Source term in stream function equation Se Source term in vorticity equation Su Source term in x-momentum equation S, Source term in y-momentum equation Sd Source term in continuity equation T Characteristic time, or time period in pulsatile flow t Nondimensional time variable U Average velocity u Velocity component in x-direction V Velocity component in y-direction X Space variable in axial direction y Space variable in radial direction M E D I C A L P H R A S E S Atherosclerosis Atherogenesis Coarctation Endothel ium Hemodynamics Histology Hypertension Hypotension Induration in vivo Lesion Lumen Morphology Nosology Pathology Platelet Sclerosis A degenerative disease of arteries Pathology of atherosclerosis Confinement to a narrowing space A single layer of thin flattened cells that lines internal body cavities Haemodynamics, fluid dynamics of blood flow Tissue structure or organization Abnormally high arterial blood pressure Abnormally low arterial blood pressure The act or process of growing hard or the state of having growing hard In the l iving body of a plant or animal Injury, or an abnormal change in structure of an organ or part due to injury or disease The cavity or passageway of a tubular organ A branch of biology that deals with the form and structure of animals and plants A classification or list of diseases or treatise comprising such a classification The study of diseases, their essential nature, causes, and development, and the structure and functional changes produced by them A minute flattened body Pathological hardening of tissue produced by overgrowth of xvn fibrous tissue and other changes Stenosis A narrowing or constriction of the diameter of any passage, tube, or orifice Thrombus A clot of blood formed within a blood vessel and remaining attached to its place of origin x v u i A C K N O W L E D G E M E N T S Firs t , I would like to express my gratitude to my supervisors, Professors B . R. Sey-mour, and V . J . M o d i , who remained a stimulating source of guidance and encouragement throughout the course of this study. Working with them has been a privilege and a plea-sure. Next , I would also like to express my sincere appreciation to Professor R. Bal iga in M c G i l l University for the interest shown in this work. His enthusiasm and valuable suggestions contributed much to the progress of this research. Thanks are due to the technicians of the Department of Mechanical Engineering. Without their expertise and ingenuity, the experimental part of this study would not have been possible. Last , but most of a l l , 1 would like to thank my wife, Helen, for her patience and loving support. Chapter 1 I N T R O D U C T I O N 1.1 Prel iminary Remarks Atherosderosis, defined by the World Health Organization in its Technical Bulletin in 1958, is 'a variable conabination of changes of the intima of arteries consisting of local ac-cumulations of lipid, complex carbohydrates, blood and blood products, fibrous deposits, and calcium deposits associated with medical changes'. It is the largest single cause of death in the western world [1]. Despite its comparatively short existence as a noso-logical entity (Marchand coined the term atherosclerosis in 1904), the process does not belong to our time only. Fallopius [2] described 'degeneration of arteries into bone'; and Cowper (quoted by Long [3]), percipiently noted that 'the passage of blood is impeded into thickened arteries'. However, the pathology of atherosclerosis is still not completely understood, despite continuous efforts over many years. Among various hypotheses and theories, it is clear that hemodynamic factors make important contributions to both the initiation and localization of atherosclerosis [4,5]. Extensive topographical data relating the distribution of atherosclerosis lesions have been suntimarized by Woolf [1], and interpreted in a variety of fascinating ways. Texon [6,7] has been a strong protagonist of the view that regions where there appears to be a predilection for atherosclerotic lesions to occur are the segmental zones of diminished lateral pressure. The common pathological mechanism at work in all these anatomical situations, as held by Texon, is a suction action in the vessel wall that is produced by a localized decrease in the pressure. In his view, such zones are characterized by curva-ture, branching, tapering, etc.. Caro and his associates [8] maintain that the maximum reduction of pressure from the stagnation value that can occur through the 'Bernoulh mechanism' is approximately 15 mmHg (assuming a maximum aortic blood velocity of 200 cm • s~^). In their view, it is difficult to explain how a reduction of the positive intraluminal arterial pressure from 100 to 85 mmHg could exert any inwardly directed suction force on the endothelium. Measurement of shear stress in vivo is a matter of some considerable difficulty. Fry [9 has succeeded in designing an ingenious system in which the histological changes in the endothelium, which are associated with increases in blood velocity, can be studied and quantified. The investigation suggests a 'critical' value of time-smoothed shear stress to which Ihe endothelial ceU 'yield' developing altered structural and chemical characteris-tics. The critical value of the stress is around 380 dyn • cm'"^ and the yielding process has a relatively short time-constant (probably less than 1 hour). Importantly, a shear stress level lower than the critical value for cell erosion can also produce more subtle alterations in the endothelial cells which might promote atherogenesis. Whether the shear stress of this magnitude can occur in the normally functioning vascular bed is unknown. One of the contradictory observations is that at points of branching, the flow divider edge, which is chronically exposed to high shear stresses, tends to be spared, at least in the early stages of plaque formation [5]. Fry [10] pursued this apparent paradox and argued that certain structural and functional changes in the arterial intima occur as a result of the mechanical stress exerted on the vessel lining by the adjacent blood flow. These resp.">nses are modulated by the magnitude and, no less importantly, the stability of the stress pattern as well as the duration of exposure to the stress. The flow divider edge experience a unidirectional shear stress although its value is relatively high. Cajo, Fitz-Gerald, and Schroter [11] found that fatty streaks and small fatty plaques tend to occur in those regions in which wall shear rate is expected to be low. The development of such lesions is inhibited or at least retarded in those areas where wall shear rates are expected to be high. It was suggested by Caro et al. [8] that this was due to the effect of shear rate on mass transport of lipid-rich macromolecules out of the vessel wall into the lumen. It has been suggested that the 'dwell time' of a platelet in the vicinity of an endothelial surface is a function of the local blood velocity gradient [12]. Areas of low shear rate and flow separation have a particularly long 'dweU time' suggesting the likelihood of platelet-endothelium interaction in the low shear stress region. A dynamic factor that might influence platelet adhesion without necessarily producing a direct effect on the arterial endothelium is stagnation, which can occur downstream of side-branch openings or distal to stenoses. A study by Fox and Hugh [13] suggests that microthrombus formation corresponds with such stagnation regions. Muller-Mohnssen, Kratzer, and Baldauf [14] proposed that at stagnation points flow has a thrombogenic effect by increasing the chances of collision between platelets, other elements of blood and the underlying endothelium. Although the blood flow is generally laminar under normal physiological conditions, unstable and turbulent flow characters may be expected where the velocity is locally increased, for example, due to narrowing of the arterial lumen. It has been held re-sponsible for a number of diverse pathological changes such as the formation of intimai microthrombi [15], modification of endothelial permeability [16] and the stimulation of intimai smooth-muscle prohferation [17,18]. Therefore, the increased colhsion between the formed elements of the blood as well as between the elements and the vessel waU due to the turbulent effect may accelerate the atherogenesis process although it may not responsible for the initiation of the atherosclerosis. One should be very careful in assessing the theories or hypotheses mentioned above since there is no 'direct ' evidence for the involvement of f luid dynamics i n the atheroscle-rotic process [19]. Although it is known that the abdominal aorta, the carotid and coronary arteries, the peripheral arteries, such as the iliacs and femorals (Figure 1.1), as well as regions of arterial branching and sharp curvature, have the greatest predilection for the disease [5], it can not be explained by a single factor. Theories are inconsistent, sometimes contradictory due to our l imited knowledge of the detailed pattern of the disease and fundamental understanding of the associated fluid mechanics. For example, the hemodynamic mechanism has been thought to be associated with increased lateral pressure [20] - [22]; reduced pressure [6,7]; reduced shear stress [8] and increased shear [5,9]. Hence, any effort at acquiring more information and understanding about the flow field is useful particularly even if it results from a greatly simplified, but consistent flow model." 1.2 F luid Dynamics in Stenosed Arteries M a n y researchers have contributed to the study of flow in constricted tubes. Young and Tsai [23] studied some of the flow characteristics in models of arterial stenosis un-der steady and pulsatile flow conditions. Their experiments yielded a description of the extent of the separated flow region, and a measure of the pressure loss across the con-striction. Lee and Fung [24] obtained numerical solutions to the Navier-Stokes equations for flow through a constricted tube. The steady flow field was calculated in the Reynolds number range of 0 to 25. The separation was predicted to occur at a very low Reynolds number. Calculations were not extended to higher Reynolds numbers due to instabihty of the numerical procedure. This was, perhaps, caused by the improper treatment of the boundary conditions when the stream function-vorticity formulation was used. Morgan and Young [25] employed momentum and energy integral equations, based on the work of Forrester and Young [26], to find approximate solution to the axial velocity component of flow in the region of stenosis. This was extended to higher Reynolds number covering the physiological range by Fukushima et al. [27]. Agreement with the experimental data of Young and Tsa i [23] for separation and reattachment points, and pressure drop across the stenoses was reasonably good although the axial momentum equation was simphfled by neglecting the second derivatives of velocity i n the axial direction. Deshpande et al. [28 solved numerically the steady Navier-Stokes equations for the case of a symmetric con-stricted tube using the stream function-vorticity formulation. The results were obtained for several stenosis-models w i th different area reductions, and the Reynolds number i n the practical region. Reasonable agreement wi th the previous reported data was ob-tained. However, the discretization of the constriction geometry and related boundary conditions were not very accurate. The study of blood flow in stenosed arteries may provide the necessary information to verify theories of localization and the init iation of atherosclerosis. Regardless of the exact init iat ing factors, the development of localized arterial stenosis may lead to disor-dered blood flow within and downstream of the constricted regions. Thus , the abil i ty to describe flow through a partial occlusion adds to the insight needed to solve the puzzle of the pathology of atherosclerosis. There is also a widespread interest in determining whether the disordered flow patterns can be used to detect localized arterial disease in its early stages, particularly before it becomes chnicafly significant [29,30]. Beyond that, introduction of a stenosed flow condition by implementing a coarctation in an artery is quite common in animal experirrients [31,32]. It would be useful to combine the fluid dy-namic factors and the results from in vivo experiments attempt to arrive at a reasonable conclusion. Information on pulsatile flow in stenosed arteries is rather l imited. Daly [33] presented results on pulsatile flow through a stenosed canine femoral artery. A Lagrangian-Eulerian procedure was employed to solve the Navier-Stokes equation. The relationship between flow quantities such as time averaged peak wall shear stress, pressure drop across the stenosis, and the stenosis height was discussed. Unfortunately, the case of severe stenosis was not covered and the information was related to a specific flow situation. O 'Br ien and Ehr l i ch [34] have also considered the pulsatile flow in constricted arteries by solving the unsteady Navier-Stokes equations in the stream function-vorticity form. The stenosis was handled by using the conformai mapping technique. Only a mi ld degree of area reduction was investigated. Concurrently, experimental investigations wi th sophisticated tools and techniques have been in progress. A Laser Doppler Anemometer ( L D A ) has been used to mea-sure the velocity field in a constricted tube [35,36]. Al though data on the velocity were quite accurate away from the sohd waU, the L D A system accuracy suffers near the wall . The separation length data, obtained through flow visualization by Young and Tsai [23], are stiU considered to be the most rehable one. Considerable effort has been directed towards to understanding of physical aspects of the problem of flow in constricted arteries. However, in general the results are often either inaccurate or incomplete; and the numerical methods used have certain limitations arising from the form of the governing equations and boundary conditions. There is a considerable scope for innovative contribution here. 1.3 Flow M o d e l Blood , a fluid with suspension of red and white ceUs in plasma, behaves as a non-Newtonian medium in small blood vessels but can be treated as Newtonian in larger arteries [37]. The flow is highly pulsatile and almost entirely laminar in a healthy cir-culatory system, although the peak Reynolds number can be of the order of 10,000. The flow is usually three-dimensional in the distensible passageways, with an average Reynolds number under 1,000 [37]. Thus the blood flow in the arterial system can be A general simulation of the actual circulatory system is indeed difficult, and perhaps beyond reality. However, through judicious simphfication retaining essential features to capture physics of the problem, useful information can be extracted. Thus one may ap-proach the problem in steps of increasing complexities to obtain some appreciation of the complex phenomenon. The blood vessel can be considered as rigid, as an approximation, since the induration of arteries usually accompanies atherogenesis [4]. Although the flow is generally three-dimensional, the two-dimensional simulation will certainly reflect some characters of the general three-dimensional flow. Although the in vivo circulatory system is unsteady, one may begin by studying the steady flow. Each, a small step, attempts shine light on a piece of a great puzzle. In the present study, fluid dynamics of a stenosed artery is simulated as the in-compressible, two-dimensional flow through a rigid tube with a rather simple geometry constriction. Both steady and unsteady situations are considered. The fluid is treated as Newtonian since the focus is on the flow through large arteries. The full Navier-Stokes equations have to be utilized because of an extensive separated flow region. 1.4 Computational Methods and Limitations The equations of motion for an incompressible fluid can be written as: vt + V • v v - -I- V P = 0; (1.1) V v = 0; (1.2) where v is the velocity vector; p the pressure; and Re the Reynolds number. The pressure is a mysterious quantity in incompressible flow from the computational point of view, as there is no exphcit equation for pressure. This arises because the assumption of incompressibility is a singular limit of the compressible flow equations. For compressible flow, the conservation of mass implies that: Pt + v - V / ' + Z ' V - v = 0; (1.3) where pressure p and density p are related through the equation of state p = pip). (1.4) If p in (1.3) is replaced by p, using equation (1.4), an explicit equation for p is obtained as: cpt + cv • "^p + V - v = 0; (1.5) where c = In the limit of incompressible flow, c = 0, and the pressure disappears from the continuity equation. There are two approaches to resolving the problem of absence of exphcit equation for the pressure. Either an equation can be derived from the Navier-Stokes equations with the incompressibility constraint (see, for example, Harlow and Welch [39]), or the pressure may be ehminated from the Navier-Stokes equations to obtain a vorticity transportation equation. However, there is no boundary conditions available for either pressure or vor-ticity. Although these two approaches are related, the problem appears quite differently, and it would be useful to discuss them separately. The finite difference methods for obtaining an exphcit pressure equation can be clas-sified into three major categories: the methods of artificial compressibility [40]; pressure correction [41]; and pressure Poisson equation [39]. In the artificial compressibility method, Chorin [40] used the discretized continuity equation at the boundary to provide the condition for pressure in his original work. The continuity equation is then satisfied everywhere. Patankar and Spalding [41] derived an approximate equation from the continuity and momentum equations for pressure during the iterative procedure by using the idea of staggered grids in Mark -and -Ce l l method [39]. The algorithm, known as S I M P L E , is widely used in practise. The pressure Poisson equation derived by taking the divergence of the momentum equation (1.1) with a substitution of continuity equation (1.2) gives another way to obtain an equation for pressure. However, the precise equivalence of the pressure Poisson equation and the continuity equation is a concern here. Since the equation for pressure is originally impl ic i t , the major difficulty faced by the various methods is the lack of a pressure boundary condition. It is often overlooked that the boundary condition problem and the problem of determining the correct exphcit pressure equation are closely related. As discussed by Gresho and Sani [42], the correct boundary condition for the pressure Poisson equation is simply the momentum equa-tion. It is, however, not quite correct to apply the momentum equation as the pressure boundary condition in the projection method [43]. The projection method, originally pro-posed by Chor in [44], sets up an auxihary velocity field v* and solves the Navier-Stokes equations in two fractional steps. The pressure operator V P = ^ ( V - v ) (1.6) projects V* on a divergence free velocity field v . It, therefore, transforms the problem for a pressure boundary condition into that of a boundary condition for the auxihary velocity v*. It is not correct to apply the momentum equation as the boundary condition if the boundary value of v* is not handled properly. Al fr ink [45] has offered an exphcit description of the required compatibihty condition to ensure full equivalence between the pressure Poisson equation and the continuity equation. It is important that the equiva-lence between the derived equation for the pressure wi th related boundary condition and the incompressible constraint is maintained. The influence matr ix method by Kleiser and Schumann [46] successfully provides the correct pressure boundary condition by ensuring the divergence-free condition in solution of the pressure Poisson equation. Although these methods have been successfully applied in solutions of flow problems, there have been several disadvantages. First of a l l , discretization of the continuity equa-tion at the boundary in the first approach leads to an inefficiency as wefl as inaccuracy of the solution. Secondly, the equation for the pressure correction is an approximation of the original continuity equation, and the accurate discretization near the boundary is difficult, especially for complex geometry. The available variations of the th ird approach are either too complicated for practical use or fail to ensure the divergence-free condition. Another approach, for two-dimensional flows, is to solve the vorticity transport equa-tion by ehminating the pressure terms from the momentum equations, and consider the vorticity, (, as a primary unknown: v ' ^ = - C ; (1.7) i - V ^ C = Ct + rPyC-^.Cv (1-8) This stream function-vorticity method has been widely used for studying two-dimensional flows. However, the condition at the sohd boundary is in terms of the stream function V': where n is the outward unit vector normal to the boundary and s is measured along the boundary. There is no boundary value available for the vorticity (, while boundary condition for stream function ip is over-specified. Methods proposed to overcome this problem can be classified into four types: con-structing a biharmonic function; using an implic it vorticity condition; using integral constraints; and the influence matr ix method. In the first approach the biharmonic equation for ip is derived by ehminating ( from equations (1.7) and (1.8), i - V ' ' V = v V t + 0 , V ' V ' x - î/'x V ' 0 v (1-10) This is a nonhnear fourth-order equation with enough boundary conditions for ip. The problem of vorticity boundary conditions is completely avoided. But the discretization of such an equation gives a relatively full algebraic system and the solution procedures are sometimes unstable [47]. Schreiber and Kefler [48] have presented a stable algorithm to solve such a system by decomposition. B u t it demands relatively more computer effort compared to that involved in simply solving the coupled system of and (. A n alternative is to make use of the second or higher-order form of the vorticity condition at the boundary, as was first suggested by Woods [49]. The method produces an explicit vorticity condition by combining a Taylor expansion and the stream function equation (1.7) at the boundary. Accurate results were obtained by G h i a et al. [50], who employed a third-order form of vorticity condition in their mult igrid procedure. The th ird approach, first proposed by Quartapelle [51], and later reviewed by Dennis and QuartapeUe [52], provides integral constraints for the vorticity and stream function. The method makes use of Green's theorem to transform a local implic it boundary condi-tion into a global integral constraint for the vorticity, and has been applied successfully to study of flow around a circular cylinder and other problems. Finafly, the influence matrix procedure suggested by Kleiser and Schumann [46] was later adapted in the spec-tral methods [53]. It employs the idea of decomposition and source functions to obtain a boundary condition for the vorticity. Both the integral constraint and influence matrix methods are mathematically simple and elegant. However, they can become computa-tionally demanding for problems i n two or three dimensions, and those wi th complex geometries. The computational difficulties involved i n solution of problems involving primitive variables or the stream function and vorticity variables are essentially similar. They can be resolved by using the approach proposed in this thesis. It is referred to as the method of necessary constraints ( N C ) . B y carefully examining the governing discretized equations for incompressible flows, it was observed that there exists a need to impose necessary constraints to solve the Navier-Stokes equations numerically i n either the pr im-itive variables, or stream function and vorticity variables. The constraints used to obtain a discretization system, or to ensure equivalence of the pressure Poisson equation to the continuity equation can be provided in several ways. In general, they can be divided into two categories: boundary constraints ( B C ) by providing the constraints at the boundary; and interior constraints (IC) by giving the constraints inside the domain. Compared with previous methods, the IC method wi l l provide the boundary conditions for pressure and vorticity in a much simpler way. The computational process is considerably simplified since no boundary value is needed. This is particularly important when the geometry of the boundary is comphcated. 1.5 Scope of Present Investigation Because our level of understanding of physiological flow problems, such as flow in stenosed arteries, is relatively poor, numerical solutions to simphfied problems are helpful in pro-viding a foundation for approaching the real situation. The lack of rehable experimental methods to determine the wallshear stress and the size of the recirculation region means that numerical simulation is a useful tool in obtaining such flow information. Numerical simulation using the ful l Navier-Stokes equations, however, encounters a computational difficulty of finding a consistent pressure or vorticity boundary condition at a solid boundary. Exist ing techniques usually result in complicated algorithms, which contribute to inefficiency and inaccuracy. The original objectives of this investigation were (i) to devise and test a simple compu-tational method for the incompressible Navier-Stokes equations for general flow problems which resolves the boundary condition difficulty, (ii) to compute the flow in simphfied models of stenosed arteries under steady and unsteady flow conditions, and (iii) verify some of the computational results by experimental measurement. The first objective was met successfully by proposing a new finite difference method for both primit ive variable and stream function-vorticity formulations. Encouraging computational results were ob-tained'for both steady and unsteady flow. The th ird objective was achieved for a l imited range of parameters. Facial External carotid Carotid bifurcation Inferior Level of inguinal ligament Femoral Saphenous Tibial Figure 1.1 A diagrammatic representation of the major branches of the arterial tree. (After Pedley, T . J . [37]) Chapter 2 I N I T I A L A N D T W O - P O I N T B O U N D A R Y - V A L U E P R O B L E M S This chapter reviews the basic mathematical background to finite difference approxima-tions in initial and two-point boundary-value problems. These are related to fluid flow problems in the primitive, and stream function-vorticity variables. The results are then apphed to two-dimensional flow problems. 2,1 First -Order System A system of first-order ordinary differential equations is considered in the beginning, which is closely related to the Navier-Stokes equations written in the primitive variable form. 2.1.1 Initial-value problems (IVPs) Consider the initial-value problem: D u = f, X > 0; (2.1) u(0) = a; (2.2) where D denotes the derivative respect to a;; u is a known vector (ui, uj, • • •, u„)^; f is a given vector (/i , /z, • • •, / „ )^ ; and a = («i , «2 , • • •, a„)-^ is also specified. This system of equations can be solved by several schemes of varying order. There is a particular interest in solving it by centered differences, because that is the usually employed for the discretization of pressure derivatives in seeking solution to the Navier-Stokes equations in pr imit ive variable formulation. The difference operators Do, D+, and D _ are defined as: Dofi = fi+i - fi-u D+fi = fi+i-fu D.fi = fi-fi.i; (2.3) where / , = / ( x , ) . Dh is used as a general difference operator. It can be one of the three kinds defined above without specification. When a centered difference formula is used to approximate the first-order derivative, the system (2.1), (2.2) becomes: ^ Z ) o U , = f , - , x, = ( z - l ) / i , ^ = l , 2 , • • • ; (2.4) u i = a ; (2.5) where h is the uniform mesh size. The system (2.4), (2.5) is not determined, and a second condition resulting i n a determined algebraic system has to be provided by applying the differential equation (2.1) at x = 0 using a one-sided difference scheme. This is a standard approach in the literature. However, if one checks this approach carefully, it is observed that a second condition is needed because the centered difference does not connect the neighboring points. For example, if / = 0, U2i+i = Cg U2i = Ce satisfy equations (2.4), (2.5) wi th Co = ct and Cg be an arbitrary constant. On these two sets of grid points (even and odd), the solution is separated. A connection between these two sets has to be provided to obtain a correct solution, and this is done in the standard approach by providing a second ini t ia l condition at X = 0. Actual ly , this connection need not be at x = 0, but could be anywhere in the domain. Here the condition which provides this connection is denoted as the necessary constraint. This observation can be generalized in the following conjecture: CLonjecture 2.1 The system (2.1), (2.2) can he discretized by a centered difference scheme if there is a necessary constraint between the neighboring points. The constraint can come from discretizing equation (2.1) at any location x (e.g., x = Q) by a one-sided difference scheme. The scheme is second-order accurate if the constraint also has second-order accuracy. 2.1,2 Two-point boundary-value problem for a two-equation system Consider a system of two equations with separated boundary conditions which is related to the Navier-Stokes type equations. The system is: Du = {, 0<x<l; (2.6) u\0) = a, u'{l) = /3; (2.7) where u = {u^,u^)'^ and f = (/^/^)^. In general, equations (2.6), (2.7) has a unique solution when they are coupled, i.e. when f is a function of both the dependent variables and u^. For example, if = and = 0, solution is determined by the boundary conditions (2.7). The centered difference approximation of system (2.6), (2.7) is: ]-DoU,=ï„ x, = ii-l)h, ^ = 1,•••,7V; (2.8) n ul = a, u], = 0. (2.9) Similarly as in the initial-value problem, two necessary constraints are required for solv-ing (2.8), (2.9). However, the necessary constraints needed to connect the neighboring points are not necessarily derived from both the equations since they are coupled. The constraints can be obtained from one equation only for both the dependent variables. This is summarized in the following: Conjecture 2.2 The system (2.8), (2.9) becomes a determined system when constraints for each variable and are given. The constraints can be provided for both the variables either by applying: • the discretized equation for v? at boundary; or • the discretized equation for using a one-sided difference scheme at two points inside the domain. 2.1.3 One-dimensional Navier-Stokes type equations The two-equation system consider here is a direct analogue of the two-dimensional Navier-Stokes equations [43]: -Dp + g + D'^u = 0\ (2.10) Du = f, 0 < x < l ; (2.11) u(0) = a , u{l) = 0. (2.12) These equations are derived from the two-dimensional Navier-Stokesflequations by as-suming periodic variations infithe y-direction. QHere it is assumed that the solution of the corresponding equation for fithe other velocity component v is known.JlThus g = —Pu and / = f{x) is a given function, ^ is a constant. There is no boundary condi-tion given for p. This means that p is determined only within a constant. A compatibihty condition, the integral constraint /3 = f^ fdx -\- a , must also be satisfied. The solution (ue,Pe) of equations (2.10) - (2.12) can be found by integration: Ue = a -f / f{x)dx; = po + / [9 + D^u)dx; Jo Jo where po is a constant. Hence equations (2.10) and (2.11) have a unique solution if po is given. Consider now the determined system: - Dp + g + D^u = Q; (2.13) Du = / , 0 < X < 1; (2.14) u(0) = a, f\dx = j. (2.15) ^0 The integral relation wiU determine the constant pq. The discretized form of this system by centered difference is: - ^DoPi + Qi + j^D+D.Ui = 0; (2.16) ^ D o u . = /.-, i = 2 , . . . , i V - l ; (2.17) N uo = oc, ^ a , j j , = 7 . (2.18) 1 As in the two-point BVPs, the difference system is indeterminate unless two con-straints for tx,- and pi are provided. The following conjecture proposes a way to impose the necessary constraints: Conjecture 2.3 The system (2.16)-(2.18) becomes determinate when the constraints for each variable pi and u,- are provided by applying: • the discretization of differential equation (2.13) on the boundary; or • the discretization of differential equation (2.13) using a one-sided difference scheme at two points inside the domain. 2.1 A One-dimensional Poisson-type approach Instead of solving the system of equations (2.10)-(2.12), one can obtain a new equation for p by taking the derivative of equation (2.10) and using equation (2.11), D^p =Dg + D^f. (2.19) Condition (2.19) now replaces equation (2.11) in the system (2.10)-(2.12). This is fre-quently used in the numerical simulation of incompressible flow problems in two dimen-sions, and is known as the pressure Poisson equation approach. The equivalence of equations (2.11) and (2.19) is now ensured if and only if Du = / is imposed at some point in the domain. However, the enforcement of Du = / , referred to as the necessary constraintior the Poisson type approach, need not.be on the boundary, but can be apphed anywhere inside the domain. Then the equivalent system to (2.10)-(2.12) is: - Dp -I- ^ -f D^u = 0; (2.20) D^p = Dg-\-D^f, 0<x<l; (2.21) u{0) = a, uil) = l3; (2.22) Du = / , some Xc €[0,1], / pdx = r (2.23) Jo The integral relation determines the integration constant for p. The corresponding dis-crete form is: - ^DoPi gi + ^D+D.ui = 0; (2.24) D^'pi = Dgi + D''fi, i = 2, • • •, iV - 1 ; (2.25) uo = a , UN = /3] (2.26) 1 ^ —DoUi^fi, some i e [l,iV], Ç a . p . = 7. (2.27) Note that equations (2.25) are obtained from the discretized form of equation (2.20) and not directly from equation (2.21). This point will be discussed further for the two-dimensional case. The type of approach is determined by the chosen value ic , and the corresponding index i. The discrete system (2.24)-(2.27) is not determined unless necessary constraints are imposed on u and p. This is summarized in the following conjecture: Conjecture 2.4 The system (2.24)-(2.27) becomes determinate and equivalent discrete form of (2.10)-(2.12) when the constraints for p and u can be provided by applying either: • the discretization of differential equation (2.20) on the boundary; or • the discretization of differential equation (2.20) using a one-sided difference scheme at two points inside the domain. The constraints needed to obtain a determined discretized system can be divided into two types: boundary constraints ( B C ) and interior constraints ( IC). Although this is quite obvious now for one-dimensional model problems, it is not as obvious in higher-dimensional cases. Some of the methods proposed in the literature fail to provide such constraints, resulting in inaccurate solutions, or divergence (details wi l l be given in the discussion of two-dimensional flow problems). Methods discussed in the Hterature fit the first category, the B C method. The IC method, which provides the necessary constraints without introducing extra information on the boundary, is favoured in this thesis. 2.2 Second-Order System In this section, the two-point boundary-value problem for a system of two second-order ordinary differential equations is discussed. Although it can be reduced to a first-order system, the second-order system is discussed directly because of its relation to the Navier-Stokes equations in stream function-vorticity form in the two-dimensional case. The model problem is a system of two second-order equations for dependent variables 0 and C, which has been used by Dennis and Quartapelle [52]: D\^g{xA^)D<: = f{x,rP); (2.28) = 0 < a ; < l ; (2.29) with boundary conditions for one variable tf; only: m = V'(l) = 0, Dm = Di^il) = 0; (2.30) where D and denote first and second order derivatives with respect to x. The functions / and g are given functions of ip and x. No exphcit boundary conditions are given for the variable (, while the boundary conditions for are overspecified. The differential system (2.28)-(2.30) is determined in the sense that four boundary constraints are given for two second-order differential equations. When a standard finite difference approach is used, the system is discretized as: ^ D + D _ C . + §DoC.- = /,•; (2.31) ^D+D.tl>i = - C i , X2<Xi<XN-i; (2.32) 01 = V-ivr = 0, DhtPi = D^rpN = 0- (2.33) The difference system (2.31)-(2.33) cannot be solved in the usual sense, because no con-ditions are available for ( on the physical boundary, xi(= 0) and xjv(= !)• However, if the discretization of the differential equation for C at the physical boundary is considered as the boundary condition for C, it overdetermines the system. Here the physical domain and the physical boundary are those for the original differential equations. 2.2.1 Woods method From Taylor expansion and the condition (2.33): V'2 = E T ^ V i + oih^y, = E ^ ^ V J V + oih'). (2.34) If D'^Tp and D^ip are replaced by —( and —D(^ using (2.29) and substituted into equation (2.34), boundary conditions are obtained for the variable C as: Ci = - ^ ^ 2 - ^C2; CN = - ^ ^ N - i - \CN-I. (2.35) There are two steps involved here: a condition is imposed on the boundary; and the vari-able t/) is expanded on the boundary. This may lead to complications in two-dimensional cases with irregular boundary configurations. 2.2.2 Integral constraint method Integrating equations (2.29) with weight functions w,(a;), gives f\uidx = - f\iD^i;dx, i = l,2. (2.36) Jo Jo With the weight functions (JJI{X) = 1 and u)2{x) — x, the two equations for ^ are: r Cdx = 0, C <dx = 0. (2.37) Jo Jo The relations in (2.37) are discretized and two global constraints for Ci, CAT obtained. In two-dimensional situations, the area integral can be transformed to a boundary integral by Green's theorem and a similar but more complicated formula can be obtained. 2.2.3 Influence matrix method The linear system obtained by setting f{x,il)) = f{x) and g{x,ip) — g{x) is used to illustrate the influence matrix approach. Equation (2.28) can be decomposed into three systems: D'Ca+9-DCa = -L D'^P, = -Ca; (2.38) MO) = V'a(l) = 0, Ca(0) = Ca(l) = 0; (2.39) D'c+ + g-DC+=o, D V + = -C+; (2.40) V'+(0) = v+ ( i ) = o, C+(o) = o, C+(i) = i ; (2.41) and D^C. 5 . DC- = 0, D V - = - C - ; (2.42) V'-(0) = V-(1) = 0, C-(0) = 1, C-(1) = 0. (2.43) With a+ and a_ as two constants to be determined, the superposition C = Co + <^ +C+ + or_C- and xp — ipa + Q^ +V'-i- + a-V*- are solutions of the system of equations (2.28), (2.29). Then two linear constraints for a+ and a_ , known as the influence matrix, are obtained by satisfying the boundary conditions (2.30). The system (2.28), (2.29) can be solved if the values of a+ and «_ are determined. For the nonhnear case, the system is linearized first, and an iterative procedure is then used. A similar procedure can be employed for two-dimensional problems, however with increased computational effort. 2.2.4 Interior constraint method It can be seen that the basic idea behind the methods discussed is to obtain a constraint between the boundary value of C and the over-specified boundary condition for tj;. It is, however, not necessary to obtain such a constraint for the boundary value of C- From the mathematical point of view, C is only an intermediate variable for obtaining a solution for V*- A fourth-order equation for rp can be obtained by eliminating C from the system of equations (2.28), (2.29) as D V + g{x, rk)D^i, = -fix, (2.44) Equation (2.44) with the boundary conditions (2.30) is a determined system. The equiv-alent discretized system, ^D^D.D+D.^i -f ^D^D.Do^i = -fi, (2.45) can be obtained directly from the discretized equations (2.31), (2.32). A relatively full matrix system results from this discretization and special techniques are needed for its solution in the two dimensional case [48] since the usual procedure is unstable [47]. Checking the procedure carefully, one finds that there is no need to calculate the boundary values Ci and (N. AS an intermediate variable, only interior values C» {i = 2 , N — 1) are used. Thus, even if one prefers to solve the system of equations (2.31), (2.32) instead of biharmonic equation (2.45), it is not necessary to calculate the intermediate variable, or 'constraint' on the boundary. The 'constraint'd is needed at the interior points. Here the approach is referred to as the interior constraint (IC) method. The basic idea of the IC method is the use of the boundary conditions (2.33) as equations for the discretized variable rpi at the grid point next to the boundary without using any information from d at that grid point. The Ci variable at the same grid then can be obtained from equation (2.32) instead of equation (2.31). The system (2.31), (2.32) with the boundary conditions (2.33) is now determined by this arrangement from: l^D+D.Ci + ^DoCi = /.•; (2.46) ^D+D-^i = - C o < < (2.47) Dhi^x = DK^N = 0; (2.48) C2 = -^D+D.rP,, CN-I = - • ^ D + D - V ' T V - i , (2.49) as equations at interior points; and ^1 = V'iv = 0; (2.50) are the only boundary conditions needed. If the C value at a boundary point is required, as in some physical problems, it can be obtained by applying equation (2.31) after the interior value is computed. In the IC method, a computational domain is set up inside the physical demain. First, a mesh system with uniform grid size h and nodes a;,- = (i — 1)^, z = 1, • • •, iV is set up on the physical domain, 0 < x < 1. The domain from X 2 to X A T - I is now defined as the computational domain fi"^, i.e. one grid step into the physical domain fî, and the computational boundary dCf^ is also moved one grid step in to X2 and XN^I. This is shown in Figure 2.1. Cl'^ is related to a particular mesh, and approaches 17 when the mesh size approaches zero. In sununary, the IC method produces a determined system by providing the necessary constraints on dCf^. Besides the advantages of simplicity and consistency in the discretiza-tion, the IC method is expected to introduce less error compared to other methods when there is a large variation in the solution near the boundary. Physical Boundary Figure 2.1 Illustration of the 1-d computational domain. Chapter 3 IC M E T H O D F O R T W O - D I M E N S I O N A L F L O W P R O B L E M S The mathematical statements proposed in the previous chapter can be generahzed to two-and three-dimensions, though the analysis is more complicated. The theoretical analysis is not repeated here since it is conceptually similar to the one-dimensional case. Instead, details of the formulation for two-dimensional flow problems are provided to illustrate the concept of necessary constraints. Because of the nature of the solution procedure involved in the finite difference approximation of the Navier-Stokes equations, an exphcit equation for pressure is always required. Thus three approaches will be used: the artificial compressibility, the pressure Poisson equation and the stream function-vorticity formu-lation. The artificial compressibility method introduced by Chorin [40] is quite simple, and hence suitable for explaining the concept. The pressure Poisson equation method has been adapted in many schemes, although the boundary condition for pressure remains a major difficulty. The stream function-vorticity formulation is much simpler for two-dimensional flow problems compared to the primitive variable formulation, and wiU be used in the present study concerning numerical simulation of flow in the arterial system. It wiU be illustrated clearly that pressure and vorticity boundary conditions have to be provided under the guideline of the concept of necessary constraints. It is neither an ar-bitrary choice between several equations, nor a lucky guess. Although there are different ways to impose the pressure and vorticity boundary conditions, they are acceptable only when the necessary constraints are provided. 3.1 Flow in Driven Cavity Although the method is suitable for the general Navier-Stokes equations, for simphcity its application is illustrated through the two dimensional driven cavity problem. The geometry of the problem is illustrated in Figure 3.1a. This problem is recognized as a standard test for the efficiency and accuracy of the numerical simulation of the Navier-Stokes equations, since it contains the general features such as separation and the non-linear nature of many flow problems. Accurate results for the flow variables have been published by many authors in the last ten years. For example, the results of Ghia et al [50] are in good agreement with the studies by Schreiber and Keller [48] for Reynolds numbers up to 5000; both the sets of results are considered as standards for the driven cavity flow problem. The discretized domains are defined as foUows: If h is the uniform grid size, the grids (x,-,2/,) are defined as x,- = {i — l)h,yj = (j — l)h,i,j = 1,---,N. Then the physical domain is Q,{{xi,yj),i,j = 1, • • • ,N}, its boundary is 50, while the computational domain is Cl'^{{xi,yj),i,j = 2, • • • ,iV—1} with its boundary as (Figure 3.16). 3.1.1 Artificial compressibility approach The steady-state Navier-Stokes and divergence equations are written in the form of the artificial compressibility method [40]: vt-I-V • v v - V ^ v - I - V P = 0; (3.1) ept + V • V = 0, x e fî; (3.2) dv v = f(x) or ^ = g(x), x € 5 f i ; (3.3) where e is a small parameter. For the driven cavity problem, the domain is fl : {(x, y), 0 < x,y < 1}. In Cartesian coordinates, the system becomes: Ut + UUx + VUy - ^ ( " x x + "yy) + Px = 0, (3.4) Ut + UVJ: + VVy - 4~(Uxx + Vyy) + Py = 0, (3.5) rte ept + «X + Uy = 0, 0<x,y< 1, (3.6) u = u = 0, at X = 0,1 and y = 0, u = l , i ; = 0, a< y = l . (3.7) The no-slip condition for the velocity is applied at the wall, hence the velocity is zero on all walls except at the top wall which is a moving boundary with a uniform tangential velocity. The notation for the difference operators are introduced here for later use. They are: D^fij = fi^ij - / . . i j ; D^fij = - and Dlfij = fi,j - where fi,j = f{xi,yj). Similar definitions are given for Dofij., D\fij, and D^fij. D^ and are used for general difference operators. Constraints on physical boundary (BC method) As in the one-dimensional model problem, the two momentum equations and the divergence equation are discretized on 17, ''^ '-^ ' ' 2h " 2h ^ ( P M - P " 7 ' ) + ^Dlu,,i + ^ ^ o ^ ^ ' . J = 0, (3.10) 2.< i,j < N - l , = ui, , -= 0, j = 2,---,N-l, (3.11) UNj = VN,j = 0, j = 2 , . . . , iV - 1, (3.12) u , M = y i , i = 0, i = 2,---,N-l, (3.13) Ui,N = 1, Vi,N = 0, i = 2, • • •, iV - 1, (3.14) where aU the superscripts n for the advanced time level are omitted. The constraint needed as the boundary conditions for pressure can be provided by discretizing either the momentum equations or the continuity equation on dCl. Here the formulation for the continuity equation approach is given where a one-sided difference is used for derivatives of the velocity components: -Pîf) + ih(^>''^ + ^ovij) = 0, (3.15) -^^{PNJ - P^NJ) + ^ m ^ N J + D^VNj) = 0, (3.16) i = 2 , . . . , i V - l ; jj{Pi,N - PIN'] + ^{D>i,N + Dlvi,N) = 0, (3.18) i = 2 , . - - , 7 V - l ; Dlmj = - ( 3 u i , , - 4 u 2 j + u3j ) , (3.19) DluN-i,j = {3uN,j - 4u iv - i , i + UN-2,j), (3.20) j = 2 , . . . , i V - l ; Dlvi,i = _(3u.-1 - 4Î;.-,2 + u.-,3), (3.21) Dlvi,N = (3f,-,iv - 4u.-,Ar-i + Vi^N-2), (3.22) i = 2 , - - . , i V - l . This approach was used by Chorin [40] when he proposed the method of artificial compressibility. It did not, however, give good results while it converges very slowly. The reason is simple. One expects to obtain the necessary constraints for both u and p by applying the continuity equation on dO,. But the constraint for p provided by the continuity equation is weak since the time derivative of p nearly disappears at the end of the iteration when the steady state is reached. One can expect that the necessary constraint can be obtained by applying the normal momentum equation on ôfi, which is often used in the literature. It will not be discussed further in this thesis. Constraints inside the domain (IC method) As m the one-dimensional model problem, the constraints for pressure and velocity, when their first derivatives are discretized by centered differences, can be provided by replacing the centered difference in one of the equations on dO,'^ by one-sided derivatives. The one-sided derivatives are used for the first-order derivatives of pressure in both momentum equations on dQ'^. The discretized equations on 0,'^ are: - ^^2(^1^- + ^+^-)« . - . . + à ^ K " ' = 0' (3-23) - YT^i^l^- + DlDDvij ^Dlplj' = 0, (3.24) ' ( K . - P 7 J ' ) + ^{D'o^iJ + mvij) = 0, (3.25) 2 < i j <N-1, = ^ i J = 0, J = 2, • • •, - 1, (3.26) UN,j = VN,: = 0, i = 2, • • •, iV - 1, (3.27) «,M = f.M = 0, i = 2, •. •, - 1, (3.28) Ui,N = 1, Vi,N = 0, i = 2, • • •, iV - 1, (3.29) where all the superscripts n for the advanced time level are omitted. D%pij and D^pij are discretized by central differences, except at the computational boundaries where one-sided differences are used, i.e., DlP2J = -(3p2,;-4p3,i+P4,i), (3.30) ^hPN-lJ = {^PN-l,j — ^PN-2,j + PN-3,j), (3.31) j = 2 , . . . , A r - i ; Dlpi,2 = - ( 3 K 2 - 4 K 3 + K 4 ) , (3.32) ^hPi^-i = (3p.-,Ar-i - 4p,-,iv_2+p,-,iv_3), (3.33) i = 2,---,N-1. Thus no pressure boundary condition is needed on ôfi, while the boundary conditions for the velocity components are given on dQ. Central differences are used for convective terms in both the methods. These are replaced by second-order upwind differences for high Reynolds number calculations. As in the one-dimensional model problems, large errors can be introduced by applying the constraint on dVt while the second approach, the IC method, is expected to give better results when providing the constraint inside fi, i.e. on dO,'^ since a boundary layer exists near the sohd wall. 3.1.2 Pressure Poisson equation approach The pressure Poisson equation is obtained by taking the divergence of the Navier-Stokes equations, - v ^ P = Pt + + uVr -\- vVy - -^V^V -h 2{v^uy - u^Vy), (3.34) Ke where X> = V ' v. Using equation (3.34) to solve for pressure, two major steps are involved: (i) some terms containing V should be omitted to obtain an equation which is not identical to the momentum equation; (ii) the new equation has to be equivalent to the continuity equation, X> = 0. There are several ways to obtain the new equation from equation (3.34), each with a corresponding compatibihty condition. Orszag et al. [43] divided them into two general types: eUiptic and parabohc. For example, the approach originally proposed by Harlow and Welch [39], ensures a divergence-free flow field at the end of the iteration procedure. This approach also reduces the possibility of instabilities in the calculation [47]. However, in the iterative procedure, the equation obtained by this approach is ^ X/\V^-V^-^) + = 0, {x,y) e fi, (3.35) Re^ ^ ' A t instead of — = 0, ix,y)en (3.36) according to the analysis of Orszag et al. [43]. The scheme is eUiptic and does not automatically ensure a divergence-free flow at each time step, except in the steady-state solution. If the other terms in equation (3.34) are also dropped to obtain the pressure Poisson equation, the velocity will not be divergence-free even for the steady-state solution. For example, we obtain the following equation by dropping the term - v ^ P = Vt + V^ uV^ + vVy -\- 2{v^uy - u^Vy). (3.37) Equation (3.37) is different from equation (3.34), indicating V'2? = 0, ( x , y ) 6 n , (3.38) which is Laplace's equation. To ensure a divergence-free flow field, the necessary con-straint or the compatibility condition for equation (3.37) has to be imposed on the bound-ary as X> = 0, {x,y)eda. (3.39) A boundary condition for pressure also needs to be imposed in the calculation. It can come either from the momentum equation or the continuity equation. For the finite difference approach, however, there exists another problem which is often overlooked. For a nonstaggered grid, as pointed out by Abdallah [56], the usual discretization of the system won't satisfy the discretization compatibility condition which is derived from Green's identity. Thus the divergence-free condition is not always satis-fied. The discretization compatibility condition can be interpreted as the condition to ensure the equivalence of the continuity equation and equation (3.34) in the discretized form. After examining the numerical approach by Thompson [54], Aubert and Deville [55] came to a similar conclusion that the discrete, second-order, centered operators do not satisfy the identity = div • grad. Failure to satisfy this discretization compatibility condition wiU lead to oscillations in the iteration procedure or a solution which is not divergence-free. To solve this problem, a consistent discretization should be used. The consistent discretization for the pressure Poisson equation will confirm the identity = div{grad) and satisfy the divergence-free condition by providing the proper boundary condition for the pressure. Let The consistent discretization has to be performed for the momentum and the pressure Poisson equations as follows: Mr,- = 0, MIJ = 0; (3.40) ^{D'oMfj + DlMl) = 0, 2<i,j<N-l. (3.41) Again some terms involving the divergence of the velocity at the advanced time level in equation (3.41) must be dropped. In the following, procedures for obtaining the pressure Poisson equation, necessary constraints and boundary condition for pressure are discussed in detail. As proposed in the previous section for the artificial compressibility method, there are two possible approaches, known as the B C and IC methods. The standard driven cavity problem is used in the discussion for simplicity and conciseness. Constraint on physical boundary (BC method) From the previous discussion about the consistent discretization, the discretized pres-sure Poisson equation on i = 2, iV — 1 and j = 2,N — I requires the discretized momen-tum equations on dCt. Thus to ensure the divergence-free condition ai i = 2, N — I and J = 2, iV — 1, the momentum equation normal to dfl has to be imposed to provide the boundary condition for pressure. Thus the constraint of I> = 0 cannot be imposed on dCl, and the pressure Poisson equation has to ensure a divergence-free flow field inside Ù. The approach of Harlow and Welch [39] who dropped has to be used. Let The discretization form inside 90 is tlien written as: = ^ " . - . i , = (3.42) ^{D-oMfj + DlMl) = 0; (3.43) 2<iJ<N-l. The boundary conditions are: = uij = Uivj = Uivj = 0, J = 2, • • •, iV - 1; (3.44) Ui,i = u.-.i = = 0, u,-,Ar = 1, i = 2, • • • , i V - l ; (3.45) = Mlj = 0, J = 2, • • •, iV - 1; (3.46) = iV^;,, = 1, z = 2 , - - - , i V - l . (3.47) with the necessary constraints: It should be mentioned that the approach introduced here is not consistent in the time discretization [43]. The method introduces error of order 0(1) in the unsteady solution, although the steady-state solution is divergence-free. However, this is the only known approach which ensures a divergence-free velocity in the steady state, in the context of the B C method, for the direct approach to the Navier-Stokes equations. Thompson [54] proposed a procedure, which imposed the artificial compressibility equation on the boundary, to provide the boundary condition for pressure, by ensuring the divergence free condition P = 0 on the boundary. However, this procedure may generate a velocity field which is not divergence-free [55], as the usual discretization of equation (3.34) does not satisfy the discretization compatibihty condition. If a consistent discretization is performed, replacement of the artificial compressibility condition by the normal momentum equation wiU lead to a solution which does not satisfy the divergence-free condition near dCt. Constraint on computational boundary (IC method) As discussed in the section on the artificial compressibility approach, the necessary constraint can be apphed to the pressure boundary condition on dCt^. By providing the necessary constraint on dO,'^, there is freedom to choose the pressure Poisson equation. One can obtain the equation for pressure by dropping either T> or if a divergence-free condition is satisfied on dCl'^. Here the case where the T) is dropped is discussed. The other Ccise is similar, and not presented explicitly. The approach of Harlow and Welch [39] is partially followed, i.e., the necessary con-straint to ensure a divergence-free flow field in fi'^ is = 0 on dCi.'^. The boundary condition for pressure also can be imposed on 50*^ . The simplest way is to apply the ar-tificial compressibility equation on ôfî'^ , then the discretized system for the driven cavity flow in n*^ is: 2 < i,j <N -I; ^t: = Hi = T7- . - i ' (3-48) ^{DlMl^ -f DlMl) = 0, (3.49) 3 < iJ <N -2; ui,j = vij = UNJ = VNj = 0, J = 2, • • •, iV - 1; (3.50) •"i.i = u.M = u.-,jv = 0, u,-,Ar = l , Î = 2 , - - - , i V - 1 ; (3.51) ^ ( K i - P 7 J ' ) + + D'oVij) = 0, (3.52) z- = 2, i = 2 , . - - , i V - l , i = 2, i = 2 , . - - , i V - l . The advantage of the IC method is that the necessary constraint V = 0 can be provided on 50*^ , where the momentum equations are already imposed. One has the freedom to choose the equations which will ensure a divergence-free field and provide the pressure boundary condition, without the necessity of imposing the momentum equations. 3.1.3 Stream function-vorticity formulation The Navier-Stokes equations for steady, incompressible viscous flow can be written in terms of the vorticity C and stream function rp in the form: v V = - C ; (3.53) = V'vCx-V'xCv (3.54) The boundary conditions generally associated with (3.53), (3.54) are i/., = /(5) and ^ = 9{s). (3.55) The system (3.53)-(3.55) is discretized on Q,'^ instead of fi. The physical boundary con-ditions are: ^'(O, y) = 0, ^ '^ (O, y) = 0, 0 < y < 1; (3.56) V'(l,2/) = 0,V'x(l,y) = 0, 0 < y < l ; (3.57) ^^(1,0) = 0,V'y(x,0) = 0, 0 < x < l ; (3.58) rl){x,l) = 0,iPy{x,\) = \, 0 < x < l . (3.59) A second-order central difference approach leads to two discretized equations for C and inside the computational domain: ^ 2 ( ^ 1 ^ - + DlDlXij = ^{DUijD'oCij - D^oi>i,jDlCij); (3.60) ^{DIDZ DlD'.)^ij = -C. - j ; (3.61) i = 3,---,N-2; j = 3 , . - - , 7 V - 2 ; and boundary conditions on the computational domain as: for ; = 2,.--, C2,i = -•^(V -Sj + V'2,i-1 -1- ^2,j+l - 4V'2.i); (3.62) CN-I,J = -J^i'^N-2,} + IpN-U-l + TpN-lJ+l -- 4V ' jv-i,i); (3.63) = 1 , (3.64) = 1, -rpN-2,j, (3.65) i V - 1 ; 0,2 ^jCV'i.a + A-1,2 + ^i+i,2 - 4V'i,2); (3.66) C,i,N-\ = - ^ ( V ' t , i V - 2 + ^i-l,N-l + ^i+l,N-l --4V ' .-,jv-i); (3.67) (3.68) = 1. h jWi,N-2 - 21 (3.69) for é = 2, • • •, A'' — 1. The boundary conditions for C are derived fr^m (3.53) while those for xp are a combination of the two equations in (3.55) using a three point difference for the Neumann condition. A l l the equations and boundary conditions have second-order accuracy , i.e. (9(^^). A second-order upwind scheme for the convective terms is usually employed for the computation of high Reynolds number flows. The details are omitted here for conciseness. As mentioned in the case of the one-dimensional model problem, the usual way to solve the difference analogue of system (3.53)-(3.55) is to derive a boundary condition for C by discretizing the differential equations in Î). The values of C and are determined, except for the values of ( on dQ,. 3.2 Flow over a Backward Facing Step The flow past a backward facing step is another interesting problem. It was chosen by the organizers of a G A M M workshop as the standard test-case, for comparing the performance of a great number of codes for solving the incompressible Navier-Stokes equations [57], because of its simplicity in geometry. It is also amenable to the use of different methods. For example, the driven cavity problem (another standard test case) has two singularities at the boundary which makes apphcation of the flnite element and spectral methods difflcult. Furthermore, it is challenging to set-up precise experiments for driven cavity flows. The flow around a cylinder, although the most interesting of all, is too difficulty because it is unbounded and quickly becomes unsteady for Reynolds number as low as Re = 30. The flow over a backward facing step proves to be a good test-case because of its interesting features. AU codes produce a vortex at the same separation point (Figure 3.2), and the prediction of the position of reattachment then presents a challenge for the performance of any numerical algorithm. As pointed by Morrison and Napolitano [58], most methods face convergence difficulties when Re > 500. This is because the flow structure becomes more complicated physicaUy when Re increases, as shown clearly by the careful experiments of Armaly et al. [59]. Another challenge presented by this problem is the downstream flow condition. Proper outflow conditions must he provided to ensure the convergence of the solution procedure. Morgan et al. [57] have discussed this point in detail. The physical domain is n{(a;,j/),0 < x < lb,y < 1}. The step height = 0.5 (channel height H^ = 1). This particular dimension used by Kim and Moin [60], and Morrison and Napohtano [58], was used here for the purpose of comparing results. Since the method used in the driven cavity flow can be applied exactly except for the type of the boundary condition with the printiitive variable formulation, it is unnecessary to repeat the discussion here. The focus is on the stream function-vorticity formulation since the computational domain fi*^ was set up differently because of the change in outflow conditions. If A x and A y are the uniform grid sizes in the x and y directions, the mesh system has nodes (x,-, y^), x,- = {i - l )Ax, yj = (j - l )Ay, i = 1, • • •, iV, j = 1, • • •, Af. The boundary condition for this problem is the no-slip condition for the velocity at the waU, i.e. ^-(0, y) = /(y) = 0, Vx(0, y) = 0, ^<y<\^ (3-70) V'(x, 0) = 0, 0j,(x, 0) = 0, 0 < x < 15, (3.71) V'(x,l) = ^, Vy(:^,l) = 0, 0 < x < 1 5 . (3.72) A parabohc velocity profile is specified at the inlet, V'(0,î/) = /(y) = 2 4 ( - ^ + ^ - | - h ^ ) , V .(0,y) = 0, ^ (3-73) and a fully developed homogeneous Neumann condition at the outflow, V',(15,y) = Cx(15,y) = 0, 0 < y < 1. (3.74) As pointed out by Armaly et al. [59] from their experiments, the inlet flow profile usually has an effect on the separation length. The downstream condition also can be given by a Dirichlet condition for moderate Re {Re < 100); this gives identical results to the Neumann condition. As shown in Figure 3.3, the computational domain is defined as f2'={(zi, yj),X2 < Xi < ^MtVi ^ Vi ^ VN-i) because there is a vorticity condition at the outflow. The system (3.53)-(3.55) is discretized on fi*^ instead of 0. The discretized system, using a second-order central diff'erence approach, consists of equations for ^ and inside the computational domain: i^^^lDl + -^,DlDy.)A,, = - C , „ (3.76) i = 3 , - - - , M - l , ; = 3 , - - - , A r - 2 ; with the boundaxy conditions on ôfi"^ as: C2,i = —^{i'3,j + ^2,j-i + ^2,j+i + fj-^M; (3.77) CN-1,J = —^{fl>N-2,j + 0iV-l , i - l + V'AT-lJ+l - ^IpN-lj); (3.78) i^2j = \^3j-^f:; (3.79) V'jv-i.j = ^tpN-2,j; (3.80) f o r ; = 2 , . - - , i V - l C.-,2 = - ^ ( V ' . - , 3 + V'.-i,2 + V'.+i,2-4^.-.2); (3.81) Ci,N-i = -^ j (V' . - , iV-2 + V'i-i.iV-i + V'.+i.iv-i - 4V',-,Ar-i); (3.82) V'v,2 = ^V'.-.a; (3.83) V'.-.A^-i = j^i,N-2; (3.84) for Î = 2, • • •, M. All the equations and boundary conditions have second-order accuracy, i.e. 0{h^), where h = max[Ax,Ay]. A second-order upwind scheme for the convective terms is usually employed for the computation of high Reynolds number flows. U=1 0V777777777777777777777777Z?y Computational Domain Physical Domain Figure 3.1 (a) Geometry of driven cavity; (b) Illustration of computational domain. V777777777777777777777777Z. \////////////////////////. Figure 3.2 Schematic diagram of flow over a backward facing step. Physical Domain Figure 3.3 Computational domain for flow over a backward facing step. Chapter 4 S O L U T I O N P R O C E D U R E A N D N U M E R I C A L R E S U L T S Numerical experiments were performed wi th uniform grids using the methods described in the previous chapters. To begin wi th , the one-dimensional model problems with con-straints at different locations ( B C and IC methods) are considered, followed by two-dimensional examples of flow in a driven cavity and over a backward facing step. For the two-dimensional illustrations, the numerical accuracy is assessed through comparison wi th results reported by earlier researchers. 4.1 Solution Procedure for Two-Dimensional Problems In this section, a detailed algorithm for solution of difference systems by the IC method is developed for both stream function-vorticity and primitive variable formulations. 4.1.1 Stream function-vorticity formulation The governing equations (one parabolic in ( and the other elliptic in tp) can be solved by various time discretization schemes. In the present case, they are solved by a standard A l -ternating Direction Implicit (ADI) procedure [61]. For the flow problem at high Reynolds number, the convective terms involved are approximated by the difference scheme pro-posed by M a [62]. For example, the single-step upwind scheme for is given by: where u , j = DgV' i j - The truncation error of this scheme is T.E. = 0{h\ai,,h). (4.3) Here the parameter a , j is introduced to ehminate the highly oscillating behavior due to high Reynolds number. The values of a , j are chosen to be constants of order 0{h) so that the truncation errors remain at the second-order, 0{h'^). The overall solution procedure involves the following sequential steps: • construct the in i t ia l values of tp and C; • solve equations (3.60) and (3.61) for the x-sweep to obtain the advanced values of C and xp; • solve equations (3.60) and (3.61) for the y-sweep to obtain the advanced values of C and tp; • repeat steps 2 and 3 unt i l the convergence criterion is achieved. 4.1.2 P r i m i t i v e v a r i a b l e f o r m u l a t i o n As before, various difference schemes can be used in conjunction with the IC method. For the system of equations (3.8)-(3.14)', the standard A D I solution procedure [61] can also be used to solve the momentum equations. The higher-order upwind scheme [63 (based on a four-point difference scheme) is introduced for the convective terms, i.e. where the parameter q controls the size of the modification. Similar schemes are associ-ated with other convective terms. The overall solution procedure involves the following sequential steps: • construct the initial values of u, u, and p; • solve equations (3.8) and (3.9) for the x-sweep to obtain the advanced values of u and v; • solve equations (3.8) and (3.9) for the y-sweep to obtain the advanced values of u and v; • solve equation (3.10) for the advanced value oi p; • repeat steps 2 to 4 until the steady state is reached. 4.2 Test on One-Dimensional M o d e l Problems Two one-dimensional test-problems were chosen to assess the performance of the IC method. One is an initial-value problem, while the other is a second order two-point boundary-value problem. 4.2.1 One-dimensional initial-value problem Consider a system with one dependent variable: Du = - A u , X > 0; (4.6) u(0) = l ; (4.7) with the exact solution u(x) = e~-^ .^ The primary goal of this test is to check the accu-racy of the numerical approximations for u compared with the analytical solution. The numerical calculations were carried out using the second-order finite difference formulae (2.4), (2.5). To determine the convergence rates, or the order of truncation error a, the root mean square error (Ec) as a function of grid size h was examined by computing E,= Fa-Fa Fa h (4.8) Fa represents value of the exact solutions of u; Fc their computed solution; || • is the discrete L^-noTm on the discretized grid points, and 1 ' ' (4.9) t=i where / is the number of grid points at which the values are calculated. If the accuracy of the numerical solution Fc is of order a, then Fc=:Fa + COih''), (4.10) where C is a constant depending on the exact solution and the grid number (/) used. For two grid sizes hi and h2, the ratio of the root mean square error is: Eel ho h-2 •jc2 « 2 "-2 where Ed and Ec2 are evaluated at the same grid points. The convergence rate is, logEci - logEc2 (4.11) a = (4.12) loghx — logh2 Table 4.1 gives the convergence rate for both formulations in Conjecture 2.1. It clearly shows that a w 2.0 for both formulations, indicating second-order accuracy The computed results, Ubc obtained by B C method, and Uic obtained by IC approach are shown in Table 4.2 for /i = 0.1 and h = 0.05 with the exact solution Uexact- The so-lutions are quite accurate for a smooth solution (A = 1), but with a boundary-layer near X = 0 (A = 5), the boundary constraint method contains larger errors than those associ-ated with the interior constraint approach. This is because the errors in the constraint equations are determined not only by the grid size h, but the constant C in equation (4.10) as well. C depends on the gradient of the exact solution, which is large inside the boundary layer. 4.2.2 One-dimensional stream function-vorticity model As a simple example consider the exact solution of equations (2.28)-(2.30), which models the stream function-vorticity equation in one-dimension: V'(x) = \x'-^x' + \x'; (4.13) Cix) = l - i x - h S x ^ ; (4.14) which satisfy the equations D V ( x ) = C(^); (4.15) D^Ci^) = 30x; (4.16) where f{x) = 30x and g{x) = 0 in equations (2.28), (2.29), the one-dimensional stream function-vorticity illustration. As before, the primary goal of this test is to assess accuracy of the numerical approx-imations for C and rp compared with the analytical solutions. The numerical calculations were carried out using the second-order finite difference formulae (2.31), (2.32) with the boundary conditions (2.33). The computed results f/o.i and UQ.OS are shown in Table 4.3 for h = 0.1 and h = 0.05 with the exact solution Uexact- The convergence rate a ol ip and C obtained by two different grid sizes at the overlapping points of the computed results is predicted as « 2.0. In this Ccise, the IC method gave approximately the same results as the integral constraint approach [52]. However, the present formulation is much simpler and can be easily extended to two-dimensional case. Table 4.1 Convergence rate of IVP a (L - norm) Grid Number ^ Boundary Constraint Interior Constraint 1 1.8111 2 . 3 3 3 4 G11-G21 5 2.2775 1 . 9 0 1 8 Table 4.2 Comparison of solution of IVP X Location (x) Grid (h) U U U exact be ic 0.1 0.606781 0.607076 1 0.5 0.606531 0.05 0.606672 0.606655 0.1 0.000000 0.090186 5 0.5 0.082085 0.05 0.101345 0.084190 Table 4.3 Convergence rate of the 1-d second-order system. * U U U a a 0.1 0.05 exact global (0<x<1 ) (x=0.5 ) (X = 0.5) V 1.88 1.88 0.0338 0.0376 0.0391 C 1.94 1.90 -0.5700 -0.6103 -0.6250 * Globlal rate is calculated by the L ^ norm. 4.3 Test on Two-Dimensional Flow Problems Two standard problems, flow in a driven cavity and flow over a backward facing step, were used as the test cases for the IC method. The convergence rate as well as a comparison of vorticity values and separation lengths are presented. 4.3.1 T h e driven cavity flow Consider now the viscous flow equations for a driven cavity. This problem is recognized as a standard test for the efficiency and accuracy of the numerical simulation of the Navier-Stokes equations, since it contains the general characters such as recirculation and nonhnear nature of many flow problems. Accurate results for the flow variables have been published by many authors over the past ten years. For example, the results of Ghia et al. [50] are in good agreement with the studies by Schreiber and Keller [48] for Reynolds numbers up to 5000; and both results are considered standards for the driven cavity flow problem. However, most of the accurate solutions are obtained by the stream function-vorticity formulation. The results using the primitive variable formulation are not as satisfactory, especially for the pressure. Kim and Moin [60] utihzed the projection method to solve the flow problem in primitive variables, and their results are frequently referenced. Results using stream function-vorticity formulation The computations in this case were carried out for Re = 100, 1000, 5000, and 10000 with a 129 x 129 grid to compare with the results of Ghia et al. [50] and Schreiber and Keller [48]. For Re = 1000 and above, the upwind scheme is used for stability. Table 4.4 clearly shows that the present results for the primary vortex and minimum stream function have at least two or three digits in agreement with most of the results obtained by previous investigators. The velocity components along the vertical and horizontal lines through the cavity center are presented in Tables 4.5 and 4.6 for Re = 1000, 5000, and 10000. They show good agreement up to Re = 5000. It can be seen from Table 4.4 that the vorticity results start to show some discrepancy between the present prediction and those by Ghia et al. [50], and others [48] for Re > 5000. It indicates that the accuracy at high Reynolds number is still an open question because the corner singularities result in a large gradient of vorticity near the wall. Table 4.7 presents the vorticity value at the moving boundary for Re = 100, 1000, 5000, and 10000. It can be seen that the difference starts to appear at Re = 1000 since the wall vorticity is highly sensitive to the grid size and difficult to compute accurately. The calculations were also done for 61x61, 81x81, and 101x101 grids at Re = 100 and 1000 to assess the convergence rate. An analytical solution in this case is not available. Thus the convergence rate is calculated using three different grid sizes hi, hi, and hs. From equation (4.10), it is quite straight forward that Fc3 — Fci h^ — hf Ed, Fc2, and Fc3 are numerical solutions for grid sizes hi, hi, and ^ 3, respectively, and evaluated on the same grid point. Table 4.8 presents the convergence rates for the sequence of solutions that approach the actual solution. The vorticity value predicts second-order rate of convergence since « 2. The convergence rate of stream function value oty,, however, is not in agreement with the theoretical result. This may due to the effect of the corner singularities in the cavity flow. The streamline plots in Figure 4.1 show the development of a central, nearly circular vortex, with bottom secondary vortices. For Re = 5000, a third secondary vortex, near the left top corner, is present, while for Re = 10000, a tertiary vortex in the right bottom region, appears. This is consistent with the results by Ghia et ai [50]. The contour plots in Figure 4.2 show a center region of essentially constant vorticity surrounded by a near circular, thin region of highly oscillatory vorticity. Mesh frequency oscillations near the downstream top corner are evident at Re = 5000, earlier than those observed by Ghia et al [50] since a coarser grid is used in the present computation. This is expected as the mesh with centered differences is not sufficiently fine to resolve the boundary-layer. Results for primitive variables formulation A test on the convergence rate for the IC method was carried out for Re = 100 and Re = 1,000. The purpose of this numerical experiment is to confirm the convergence rates for pressure and velocity obtained by the present method using a primitive variable formulation. The convergence rates are evaluated for u, v, and p as a„, and ap from numerical results obtained by using three different grid sizes, hi, h2, and /13. Table 4.9 provides the computed convergence rates for Re = 100 and 1,000 obtained from grid sizes 61 X 61, 81 x 81 and 101 x 101. Once again, they are close to 2, the theoretical prediction of the second-order accuracy. As pointed before, the focus is on the comparison of two general approaches (BC and IC methods). The numerical results are presented for the artificial compressibility method first, and comparison is based on the primary vortex and minimum stream function values for Re = 1,000. Table 4.10 presents the data, with different grid sizes, for both B C and IC methods. It clearly shows that accurate results are obtained using the IC method while the B C approach does not necessarily give accurate solutions in this case. The IC method is also much more efficient than the B C method, it usually took 4 to 10 times less iterations than the B C method, depending on the grid-size and the Reynolds number. The computation for the IC method was also carried out for Re = 100, and 1000 with a 129 X 129 grid for comparison with the previous results. For the convective terms, the higher-order upwind scheme discussed earher was apphed to the case of i?e = 1,000, while the central difference formula was used for Re = 100. The compressibility parameter e = 1 and a time-step At = 0.005, required to stabilize the difference scheme from iteration to iteration, were used. For Re = 1,000, the integration time for the momentum equations was from t = 0 to < = 70, taking about 14,000 iterations to reach the steady state. The accuracy of the solutions was compared with the results of Ghia et al. [50], Schreiber and Keller [48] and the present stream function-vorticity calculations (Table 4.11). Table 4.12 compares the minimum stream function and primary vortex using the pres-sure Poisson equation method and the results of the artificial compressibility approach for IC and B C techniques at Re = 100. In general. Both the IC and B C methods of pressure Poisson equation approach are quite efficient and the results are in excellent agreement with the results of the artificial compressibility approach by the IC method, i.e. identical in the first three digits. However, there are situations where IC method gives more rehable results while the B C method is not necessarily as accurate (for example, in the artificial compressibility approach). Table 4.4 Comparison of stream function and vorticity at the center of primary vortex for different Reynolds numbers. Minimum Stream Function and Primary Vorticity Re Present Ghia et al. Schreiber & Keller V C V C C 100 -0.10336 -3.16500 -0.10342 -3.16646 -0.10330 -3.18200 (129x129) (129X129) (121 X 121) 1000 -0.11668 -2.02904 -0.11793 -2.04968 -0.11603 -2.02600 (129 X 129) (129 X 129) (141 X 141) 5000 -0.11621 -1.83445 -0.11897 -1.86016 -(129 X 129) (257 X 257) 10000 -0.11286 -1.74198 -0.11973 -1.88082 -0.10284 -1.62200 (129 X 129) (257 X 257) (181 X 181) Secondary Vortices (Bottom Left) Re Present Ghia et al. Schreiber & Keller V C C V C 100 1.39(-6) 0.0141 1.75(-6) 0.0156 2.05(-6) 0.0080 1000 2.15(^) 0.3566 2.31 (-4) 0.3618 2.17(^) 0.3020 5000 1.30(-3) 1.3330 1.36(-3) 1.5306 - -10000 1.57(-3) 1.8611 1.52(-3) 2.0856 1.12(-3) 1.0670 Secondary Vortices (Bottom Right) Re Present Ghia et al. Schreiber & Keller V C C V C 100 9.96(-6) 0.0328 1.25(-5) 0.0331 1.32(-5) 0.0255 1000 1.63(-3) 1.0304 1.75(-3) 1.1547 1.70(-3) 0.9990 5000 2.68(-3) 2.5311 3.08(-3) 2.6635 - -10000 2.66(-3) 2.8615 3.42(-3) 4.0531 2.96(-3) 3.0310 Table 4.5 Comparisons of u-velocity along the vertical line through the cavity center. Grid Location Re=1000 Re: =5000 Re=10000 Point y Present Ghia et al Present Ghia et al Present Ghia et a! 129 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 126 0.9766 0.6590 0.6593 0.4719 0.4822 0.4442 0.4722 125 0.9688 0.5749 0.5749 0.4499 0.4612 0.4510 0.4778 124 0.9609 0.5113 0.5112 0.4491 0.4599 0.4540 0.4807 123 0.9531 0.4663 0.4660 0.4502 0.4604 0.4515 0.4780 110 0.8516 0.3324 0.3330 0.3306 0.3356 0.3198 0.3464 95 0.7344 0.1063 0.1872 0.1941 0.2009 0.1906 0.2067 80 0.6172 0.0563 0.0570 0.0775 0.0818 0.0797 0.0834 65 0.5000 -0.0613 -0.0608 -0.0295 -0.0304 -0.0223 -0.0311 59 0.4531 -0.1068 -0.1065 -0.0711 -0.0740 -0.0621 -0.0754 37 0.2813 -0.2777 -0.2781 -0.2223 -0.2286 -0.2078 -0.2319 23 0.1719 -0.3826 -0.3829 -0.3212 -0.3305 -0.3035 -0.3271 14 0.1016 -0.2953 -0.2973 -0.3996 -0.4044 -0.3684 -0.3800 10 0.0703 -0.2199 -0.2222 -0.4246 -0.4364 -0.4121 -0.4166 9 0.0625 -0.1998 -0.2020 -0.4154 -0.4290 -0.4245 -0.4254 8 0.0547 -0.1791 -0.1811 -0.3971 -0.4117 -0.4312 -0.4274 1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Table 4.6 Comparisons of v-velocity along the horizontal line through the cavity center. Grid Location Re=1000 Re =5000 Re= 10000 Point X Present Ghia et al Present Ghia et al Present Ghia et al 129 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 125 0.9686 -0.2273 -0.2139 -0.5105 -0.4977 -0.5481 -0.5430 124 0.9609 -0.2910 -0.2767 -0.5489 -0.5507 -0.5121 -0.5299 123 0.9531 -0.3514 -0.3371 -0.5434 -0.5541 -0.4732 -0.4910 122 0.9453 -0.4054 -0.3919 -0.5152 -0.5288 -0.4474 -0.4586 117 0.9063 -0.5192 -0.5155 -0.4098 -0.4144 -0.3963 -0.4150 111 0.8594 -0.4215 -0.4267 -0.3568 -0.3621 -0.3396 -0.3674 104 0.8047 -0.3159 -0.3197 -0.2928 -0.3002 -0.2778 -0.3072 65 0.5000 0.0265 0.0253 0.0147 0.0095 0.0151 0.0081 31 0.2344 0.3222 0.3224 0.2690 0.2728 0.2535 0.2722 30 0.2266 0.3306 0.3308 0.2768 0.2807 0.2608 0.2800 21 0.1563 0.3707 0.3710 0.3488 0.3537 0.3276 0.3507 13 0.0938 0.3261 0.3263 0.4175 0.4295 0.3934 0.4149 11 0.0781 0.3032 0.3035 0.4194 0.4365 0.4116 0.4312 10 0.0703 0.2897 0.2901 0.4135 0.4333 0.4177 0.4373 9 0.0625 0.2742 0.2749 0.4023 0.4245 0.4190 0.4398 1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Table 4.7 Comparison of vorticity along the moving boundary. X Re =100 Re= :1000 Re= 10000 Present Ghia et al. Present Ghia et al. Present Ghia et ai 0.0000 - - - — — — 0.0625 40.6632 40.0110 75.7930 75.5980 230.326 209.452 0.1250 22.7228 22.5378 51.7570 51.0557 145.067 145.073 0.1875 16.3789 16.2862 40.7787 40.5437 136.339 127.928 0.2500 12.7773 12.7844 32.2132 32.2953 124..917 116.275 0.3125 10.3216 10.4199 25.2903 25.4341 91.2275 90.0231 0.3750 8.53564 8.69628 20.2328 20.2666 68.3263 67.1400 0.4375 7.24853 7.43218 16.9738 16.8350 58.0016 53.5905 0.5000 6.40695 6.57451 15.1925 14.8901 53.8154 46.8271 0.5625 6.01805 6.13973 14.5205 14.0928 52.8937 44.3287 0.6250 6.12944 6.18946 14.6492 14.1374 53.9932 44.6303 0.6875 6.82795 6.82674 15.3685 14.8061 56.3608 46.8672 0.7500 8.27014 8.22110 16.6630 16.0458 59.2103 50.3792 0.8125 10.8271 10.7414 19.1392 18.3120 61.1194 54.3725 0.8750 15.8090 15.6591 25.6367 23.8707 59.6298 57.7756 0.9375 31.4365 30.7923 48.2273 42.1124 73.9035 66.0352 1.0000 - - - - - -Table 4.8 Convergence rates of solution using the IC method for the driven cavity flow. Grid Number Re a G61-G81-G101 100 2.799 2.078 G61-G81-G101 1 0 0 0 3.566 2.060 Table 4.9 Convergence rates of the solution using the IC method for the driven cavity flow. Grid Number Re a a a U V p G61-G81-G101 100 1.922 2.034 1.963 G61-G81-G101 1000 2.035 2.035 1.879 Table 4.10 Comparison of the stream function and vorticity at the center of the primary vortex for the driven cavity flow. Minimum Stream Function and Primary Vorticity: BC Method IC Method Fte Grid number 6 1 - 0 . 0 9 7 - 1 . 7 1 3 - 0 . 1 0 7 - 1 . 9 1 1 1 0 0 0 81 - 0 . 1 0 3 - 1 . 8 1 6 - 0 . 1 1 2 - 1 . 9 7 8 101 - 0 . 1 0 6 - 1 . 8 7 4 - 0 . 1 1 5 - 2 . 0 1 0 Note: one standard data set is that of Kim & Moin [60], obtained by using the fractional step method in the primitive variable form. The values are -0.116 and -2.026 for the minimum stream function and primary vorticity, respectively. Table 4.11 Comparison of the stream function and vorticity at the center of the primary vortex for the driven cavity flow. Minimun Stream Function and Primary Vorticity Pr imi t i ve Ghia et al Schreiber & Keller Stream"^ Re V C V C V C V ' C -0.103 -3.163 -0.103 -3.166 -0.103 -3.166 -0.103 -3.165 100 (129x129) (129x129) (121x121) (129x129) -0.116 -2.033 -0.118 -2.050 -0.116 -2.026 -0.117 -2.029 1000 (129x129) (129x129) (141x141) (129x129) Note:* primitive-results using the primitive variables formulation; stream-results using the stream function-vorticity formulation. Table 4.12 Comparison of the stream function and vorticity at the center of the primary vortex for the driven cavity flow. Minimun Stream Function and Primary Vorticity Artificial Compressibility Pressure Poisson Equation Re BC Method IC Method BC Method IC Method V C V C V C V C -0.102 -3.109 -0.103 -3.166 -0.103 -3.163 -0.103 -3.163 100 (129x129) (129x129) (129x129) (129x129) Re=100 ^—^ Re=1000 Figure 4.1 Stream function contours for the driven cavity flow. Figure 4.2 Vorticity contours for the driven cavity flow. 4.3.2 Flow over a backward facing step Accurate results were also obtained for flow over a backward facing step. This problem provides an excellent test-case for the accuracy of the numerical method because of the dependence of the reattachment length on the performance of the numerical scheme. The step size / i , is half of the height of the channel Ha. Li is the horizontal distance between the reattachment point of the primary vortex at the lower waU and the step (where the separation occurs). There is also a secondary vortex at the top wall when the Reynolds number increases. is the separation point of the secondary vortex, while L3 is the reattachment point of the secondary vortex at the top wall (Figure 3.2). The location of the primary vortex center is denoted by (ic, Vc)-, while the values of stream function and vorticity are tpc and Co respectively. The computation was carried out on a uniform 101 x 51 grid over a 15 x 1 domain. Figures 4.3 and 4.4 present contours of the stream function and vorticity for Re = 200, 400, 600, and 800. It can be seen that the size of the recirculation region increases with the Reynolds number. At Re = 600 and above, a secondary vortex develops at the top waU. Second-order accuracy was obtained by setting the parameter a,j = Ay in the single step upwind scheme. Table 4.13 presents comparison of the separation and reattachment points Li between the present calculation and the results of Armely et al. [59], Morrison and Napohtano [58], and Kim and Moin [60]. The agreement is reasonably good for both the primitive variables and stream function-vorticity methods. Table 4.14 summarizes the quantitative features of the flow field by presenting the vortex size, location, and values of stream function and vorticity. The vorticity values at the lower and upper boundaries aie shown in Figures 4.5 and 4.6. Agreement between the results by stream function-vorticity and by primitive variables approaches is good for low Reynolds number {Re = 200, 400). The discrepancy at Re = 600, and 800 may be caused by too coarse the grid size. The study shows that accurate solutions can be obtained for both: flow in a driven cavity; and flow over a backward facing step; by the IC naethod. The IC method is also quite efficient in terms of C P U time. Table 4.15 presents a comparison between the IC method and the SIMPLE method for the driven cavity flow with Re = 100, 61 x 61 grid. It shows that the IC method is almost 5 times faster, and more accurate, for both primitive variable and stream function-vorticity formulations by predicting results closer to those in the literature. Table 4.13 Results on separation and reattachment points for the flow over a backward facing step. Re Armaly et ai. Morison & Napolitano P r i m i t i v e Stream S l -a I-, L 2 L 3 •-2 ^ 3 •-3 200 — 5.3 5.3 5.5 400 8.6 - - 8.6 8.0 10.4 8.5 7.7 10.6 8.9 8.4 10.2 600 11.5 - - 10.7 8.7 16.2 10.4 8.4 16.3 11.1 9.1 16.2 800 14.3 - 12.2 9.7 21.0 11.8 9.3 20.9 12.6 10.1 21.0 Table 4.14 Characteristics of the primary recirculation region for the flow over a backward facing step. Re 200 400 600 800 X c 1.050 1.800 2.700 3.450 V c 0.300 0.300 0.300 0.300 p ¥ c -0.03241 -0.03347 -0.03371 -0.03373 Co -2.18306 -2.14078 -2.25827 -2.25337 5.45 8.89 11.14 12.64 x c 1.050 1.800 2.550 3.300 y c 0.300 0.300 0.300 0.300 s ¥ c -0.02565 -0.02631 -0.02635 -0.02624 C c -2.3309 -2.26507 -2.26361 -2.31358 L 1 5.32 8.51 10.43 11.76 Note: P-primitive variable formulation; S-stream function vorticity formulation; Subscript "c' stands for the center of the primary vortex. Table 4.15 Comparison of two algorithms for the driven cavity flow. (Re= 100, grid size 61x61). C P U (second) Iteration mm c . prime IC Method 312.5 4 9 0 - 0 . 1 0 2 5 - 3 . 1 3 4 9 SIMPLE Method 1907.5 250 -0.0850 - 3 . 0 1 5 5 Note: convergence criterion is Div V= 5 x 10 ; computation was performed on a Sparc 1. Chapter 4. SOLUTION PROCEDURE AND NUMERICAL RESULTS Re=200 15.0 0.0 12.5 15.0 Re=600 15.0 Re=800 0 0 2.5 5.0 7.5 10.0 12 5 15.0 Figure 4.3 Stream function contours for flow over a backward facing step. 0 4 8 12 16 X Figure 4.5 Vorticity at the lower wall for flow over a backward facing step. 0 4 8 12 16 X Figure 4.6 Vorticity at the upper wall for flow over a backward facing step. Chapter 5 N U M E R I C A L S I M U L A T I O N O F F L O W I N S T E N O S E D A R T E R I E S 5.1 Stream Function-Vorticity Formula The stream function-vorticity form of the Navier-Stokes equations has a major advantage over the primitive variable form in two-dimensional incompressible flow problems as it satisfies the continuity equation automatically. It also decouples the pressure calculation from the velocity calculation. Thus it eliminates two computational difficulties, i.e. to find the proper boundary condition for the pressure and an exphcit pressure equation which satisfies the incompressibility constraint. However, the vorticity boundary condi-tion at a sohd boundary introduces another difficulty in the calculation. The interior constraint (IC) method is used to resolve this difficulty. 5.1.1 Governing equations For the axisymmetric problems under investigation here, the governing equations can be written, in the stream function-vorticity formulation, using a cyhndrical coordinates. The coupled equations for the stream function (T/») and vorticity {() in a nondimensional form are: •^xx + rpyy =-y( + S^; (5.1) ^ 0 + ^(V'vCx - 0xCv) - Ye'^Cxx + Cyy) = (5.2) where Re = 2UR/i', Rw = R{\lvT)^ are the Reynolds number and frequency Reynolds number (or 'Womersley' parameter) [64], respectively. R is the characteristic upstream radius of the tube; u the kinematic viscosity, U the average velocity, and T the charac-teristic time, tp and C are defined as: u = -xpy, V = - - 0 ^ ; (5-3) y y c = ux - " y ; (5.4) where u and v are axial and transverse components of velocity, respectively. The expres-sions for the source terms 5"^ and are given in Table 5.1. 5.1.2 Boundary conditions This section discusses the problem of finding suitable boundary conditions. Consider the flow in a symmetrical tube with centerhne AB at the bottom as shown in Figure 5.1. The correct boundary values come from physical conditions which are discussed in detail. Upstream boundary The inflow or upstream boundary is labeled as FA in Figure 5.1. A specific profile for velocity is given, which is integrated to set the inflow boundary condition for the stream function as ip = a/Ay^ + b/2y^. Here a and b are constants whose values depend on the type of inflow profile desired (parabolic or uniform) and the radius of the tube, t; = 0 provides another boundary condition for tp as tp^ = 0. Centerhne The centerhne is labeled as AB in Figure 5.1. Axial symmetry is assumed, and a direct consequence of this assumption is zero velocity gradient in radial direction at the axis. This implies that on AB, xp = 0,( = 0. Although the centerhne stream function value is arbitrary, zero is assigned without loss of generality. Sohd boundary The solid boundary is labeled as CF in Figure 5.1. The no-slip velocity condition gives xp = a/A -\- b/2,'ipy = 0. The value of the stream function along the top boundary, including the constricted part, is a given constant which is the flow rate determined from the inlet condition. This condition causes the trouble, since the stream function boundary condition is overspecified, while no boundary condition for vorticity is provided. Most of the approaches in the hterature to resolve this difficulty can be classified as variations of the Woods method [24,34,38],. As mentioned previously, the IC method presents a simple and fundamental way to approach this problem. Here it is adapted to deal with the boundary condition for (. Downstream boundary The outflow region is labeled as BC in Figure 5.1. The outflow condition is yet another challenging problem in computational fluid dynamics since the physical condition at infinity has to be imposed at a finite downstream location. The most common type for tube flow is the fully developed condition, i.e. ipx = 0, Cx = 0, as used by O'Brien and Erhlich [34]. The Dirichlet type condition for ip and by specifying the velocity profile, usually introduces convergence difficulties, and it is only useful for very low Reynolds number flows. Mueller [38] used extrapolation in the axial direction for the'values of stream function and vorticity. This condition aUows the flow to continue in the general direction dictated by the two interior points adjacent to the outflow. In summary, the boundary conditions for tj; and ( for a parabohc velocity profile at the inlet are: . 0 = y 2 ( i _ i j / 2 ) ^ ^ ^ ^ o a t FA; • xp = 0.5, rpy = 0 at CF; • ^ = 0, C = 0 at AB] and • V x = 0, Cx = 0 at BC. 5.1.3 Computational domain The IC method represents a novel approach which resolves the difficulty of determining the correct boundary condition for the vorticity. No extra boundary condition is needed for vorticity at the sohd boundary. One can use the physical boundary condition (no-slip at the wall) directly by setting up a computational domain, inside the physical domain, where the calculations are completed. Arrangement of the computational domain is quite straight forward for the flow in a constricted tube. The procedure can be illustrate the procedure as foUows. Consider a finite difference mesh system covering the physical domain as shown in Figure 5.2. The mesh system is generated by the boundary-fitted coordinate technique [54] which is discussed in detail in the next section. The computational domain fi'^ is bounded by the computational boundary dO,'^, A'B'C'D'E'F'. dVt'^ is defined to be one grid step inside the physical boundary 5fi except at the outflow region and the centerhne. No boundary value of C is needed at dQ since it does not appear in the computation. 5.1.4 - Representation of constriction The irregular shape of the constriction is handled by using the boundary-fitted coordi-nate technique [54]. Other methods are available, for example Deshpande et al. [28] used a rectangular grid system to approximate the irregular boundary by interpolation. This procedure will introduce some errors. Lee and Fung [24] as well as O'Brien and Ehrlich [34] used conformai mapping techniques to map the bell-shaped constriction into a rectangular domain. The technique has an advantage as it satisfies the orthogonality condition, but is difficult, if not impossible, to obtain a suitable transformation for a gen-eral shape or an arbitrary boundary. A more general boundary-fitted coordinate system, which can be generated by solving a second-order partial differential equation numeri-cally, is used in this research. Although not a new idea, the method was introduced in the study of flow through constricted arteries only recently [55,65]. The following system of equations with appropriate boundary conditions transform a two-dimensional region with an arbitrary constriction in the {x,y) plane {A'B'C'D'E'F', Figure 5.2) into a square domain iCv)'-1 J 2 1 {ax^( - 2/3x^r, + JXnr,) = Rx^ + Sxr,; (5.5) ( a y « - 2^y^, + 7y, , ) = Ry^ + Syr,] (5.6) J 2 ' x{0,r,) = 0, xil,r,) = x,, x,(e,0) = 0, x (^ , l ) = ^; (5.7) y(0,7,) = 77, yil,rj) = r,, y(C,0) = 0, yi^,l) = H{0; (5.8) where: îxx + (yy = R{x, y), Vxx + Vyy = S{x, y); a = xl + y^, 7 = x| + y|, /3 = Xr,y^ - x^yr,, J = x^yr, - x„yf ; /r(^) is a given function. Usually, the system is solved numerically. 5.1.5 Equations and boundary conditions in transform domain The governing equations for the flow in the (V», C) form are now written in the transform plane (C, T}) as: 1 J2 [arP^ç - 2^rP^, + 7V',,) - {R^P^ + SiP,) = -yC + S^; (5.9) ^ C t + ^iu^'Q + v^'Cn) + Yij-2«i - 2/5C '^' + 7C..) - (RCi + SCr,)] = 5c; (5.10) where: u^" = i(uy„ - u x j , u^" = j(t;y^ - ux{). The physical boundary conditions areas follows: upstream (^ = 0)-^' = y^(l —l/2y2), V'x = 0; centerhne (77 = O)-0 = = 0; sohd boundary (7/ = 1 ) -0 = 0.5, Xr,xp^+x^ipr, = 0, outflow {i = l)-î/r,V'« - y^V-r, = 0, y„774 - y^ T/,, = 0. 5.1.6 Discretization The discretization of equations (5.5), (5.6) for generating the grid was accompHshed by the usual centered difference for both first- and second-order derivatives to obtain the second-order accuracy. It can be written as: t = l , - - - , M , j = l , . - - , i V ; where D^fij = fi+ij - fij, Difij = fij - D^fij = fi+ij - and /^j = f(^i, rjj). Similar definitions are given for D\, Dl, and DQ. Equations (5.11), (5.12) were solved by an iterative procedure with the associated boundary conditions as: xij = 0, XM,j = XI, yij = yuj = Tj, j = l,---,N; (5.13) Dlxi,i = 0, Xi,N = C, yi,i = 0, yi,N = H{^i,r^), i = l,---,M; (5.14) where D^f^^j = 3/jv,j — 4//v_i,j -f / /v-2,i is the one-sided derivative. A similar definition is given to D^. Discretization of the governing equations is much simpler and more natural using the IC method. Equation (5.10) was discretized on the computational domain, Cf^, while equation (5.9) was needed inside Cl'^ as - ( ^ ^ o V ' . - . i + ^D'oM = -VijCiJ + S^,j, (5.15) i = 2 , . - - , M - l , i = 2 , - - - , i V - l ; Re At + y , / 2 A r ° ^ " ^ ^ 2A7,^°^"^^ lij^i^DiDiOj - ^^DiDlC. ^ i^DlDlU z = 3 , - - - , M - l , i = 3 , . - - , A r _ 2 . The boundary conditions are: V'l.i = y L ( l - \yh\ DlvijDirPi, - DiyijDlrPij = 0, DlyujDii^Mj - DivMjDli^Mj = 0, DlyujDiCuj - DiyMjD^Mj = 0, > = 2 , - - . , i V - l ; Al = Ci,i = 0, rpi,N = 0.5, Dlxi,NDii^i,N + Dixi^NDl^N = 0, i = 2 , - - - , M - l . Centered differences are used for all derivatives for second-order accuracy except Z)^ and D\ on dQP. They are second-order one-sided differences. The higher-order upwind scheme [63] wiU be used for relatively high Reynolds number calculations {Re > 600). y=1 (wall): 14) =0.5, v = 0 . F I y x=0 (inlet) : H'x=0. x=xl (ouflow) C . = 0 . y=0 (centerline) : ijj = C=0. B Figure 5.1 Geometry of the flow in a tube with a constriction and its boundary conditions using the stream function-vorticity formulation. ^ A < \ 0 ' \ \ \ / 1 r c A A ' BB' F' E' D' c A ' BB' Figure 5.2. lUustration of computational domain using the stream function-vorticity formulation. 5.2 Primitive Variable Formulation The primitive (or physical) variable formulation for the fluid motion has some advantages over the stream function-vorticity procedure. It is a direct description of physical quan-tities, and is essentially the same in three- and two-dimensional problems. This makes it easier to extend any numerical schemes from lower to higher dimensions. However, the direct simulation of the primitive variable form presents a computational difficulty. There is no physical value available for pressure at a solid boundary. The IC method gives a direct algorithm for resolving this difficulty. It eliminates the use of pressure at a solid boundary, and provides a simpler way to attack the problem with a complex geometry. This nonstaggered grid approach makes it easy to implement other techniques, such as the boundary-fitted coordinates. As an example, the flow in stenosed arteries is formulated and computation is conducted. 5.2.1 Governing equations The set of governing equations consists of two momentum equations and one continuity equation in the cyhndrical coordinates: Ut + UU^ + VUy - J-(UXX + Uyy) = -p^ S'y, (5.17) rte He 2Rw^ 2 Vt + 4- VVy - -^("XX + Vyy) = - Py + Sy, (5.18) ept + u^ + Vy = Sd. (5.19) u and V are axial and tangential velocity components. The expression for the source terms Su, Sv, and Sd are given in Table 5.2. The artificial compressibility method is used, hence the time derivative in (5.19) doesn't have any physical meaning. The solution of system (5.17)-(5.19) converges to the steady state solution of the incompressible Navier-Stokes equations when the pressure term in the continuity equation (5.19) is zero. 5.2.2 Boundary conditions The physical boundary conditions for the flow in constricted arteries are discussed in this section. Upstream boundary The upstream or inlet boundary is FA in Figure 5.3. A natural choice is to specify the velocity profile: u = ay^ — b;v = 0. a and b are constants depending on the type of the profile chosen. If the fully developed profile is selected, then a = —2 and 6 = 2. Centerhne The centerline is the line of symmetry, AB, in Figure 5.3. The radial gradient of the axial velocity is zero, while the radial velocity is zero at the centerline: Uy = 0,v = 0. Sohd boundary The solid boundary is the line C F in Figure 5.3. The no-slip and no flow conditions are specified at the boundary: u = 0;v = 0. Downstream boundary The outflow or exit boundary is represented by BC in Figure 5.3. This is one of the troublesome parts in the computation of flow in tubes or chan-nels. The flow eventually becomes fully developed downstream of the disturbance. The redevelopment length depends on the Reynolds number. It is usuaUy quite long, and it would be unrealistic to cover this entire length with a computational mesh. Thus, a computational downstream boundary has to be setup where the condition for veloc-ity is provided. There are three kinds of possible boundary conditions, as discussed at the G A M M workshop on flow over a backward facing step [57]. The first kind is the Dirichlet type, i.e. a velocity profile (usually parabolic) is given: u = 2(1 — y^); t; = 0; as in the case of the upstream boundary. This approach is not quite natural and intro-duces instabihty when the Reynolds number becomes large {Re > 100) if the flow at a downstream location is far from fully developed. The second kind is the homogeneous Neumann type, i.e. the fully developed condition implies the axial derivatives of both velocity components are zero: Ux = 0; = 0. This condition is quite popular in the hterature. It is beheved to be more natural than that described in the first kind. It gives convergent solutions with a reasonable accuracy for large Reynolds number even when apphed at a downstream boundary where the flow is far from fully developed. This is because the governing equation is dominated by convective terms which eliminate the effect of the downstream boundary. The last kind is called the natural condition. The stress is assumed to vanish at the downstream boundary: p = Au^x/Re; Ux + Vy = 0. As one can guess from the name, this approach is quite natural and also very popular in the literature. It has been adapted by the G A M M workshop [57], and also by Aubert and Deville [55]. It should be mentioned here that these three kinds of boundary conditions are all correct when the fully developed flow field is established. The problem arises during"computation when the conditions are apphed at a boundary where the flow is not fully developed. There is no absolute standard to suggest preference for the particular set. One should choose the set which introduces the least error in the calculations. In summary, the boundary conditions for u and v are given as: • u = 2(1-y^),v = 0 at FA; • Uy = 0, u = 0 at AB; • u = 0, u = 0 at CF; and • Ux = 0, Ux = 0 at BC. The second kind of the Neumann type, 'fully developed' downstream boundary condition is chosen here for calculation. 5.2.3 Computational domain The setup of the computational domain is rather straight forward. Consider a finite difference mesh system covering the physical domain as shown in Figure 5.4 The mesh system is generated by the boundary-fitted coordinate technique as discussed in detail in the previous section. The computational domain f2'= is bounded by the boundary dCt'', A'B'C'D'E'F\ which is defined one grid step inside the physical boundary dQh- No boundary value for p is needed at dÇlh. since it does not appear in the computation. Thus, the boundary conditions are simply the physical conditions. 5.2.4 Governing equations and boundary conditions in transform domain The flow equations for stenosed arteries are solved in a transform (^ ,17) plane. The relationship between the physical plane (x,t/) and the transform plane (^ ,7/) is exactly the same as discussed in the previous section for the stream function-vorticity approach. The governing equations in the transform plane are: — — U f + (^^"«4 -h u^^u,,) Rt ^ - 2/5"î. + 7t^..) - + Su^)\ -t- -iyr^pi - y^p,) = (5.20) 2Rw^ Re - - ^^""ir, + iv^n) - (R^i + Svr,)] + j{-Xr,pi + x^p^) = 5„; (5.21) + -j{yr,Ui - y^Ur, - Xr,V^ + X^Vr,) = Sd] (5.22) where w^" = jiuyr, - vx^), u ^ " = j{vy^ - uy^). The physical boundary conditions are: • upstream = 0) : u = 2(1 — y^), v = 0; • centerline rj = 0 : Xj,u^ + x^Ur, = 0, u = 0; • solid boundary J7 = l : u = 0, u = 0; • outflow ^ = 1 : yr,u^ - y^Ur, = 0, yr,v^ - y^Vr, = 0. 5.2.5 Discretization The grid generation procedure is the same as in the previous section. The discretization of the governing equation is performed in the Cf^, and the boundary conditions are apphed on dCf^, by providing the interior constraints: Re At + 2 A e ^ ° " ' ' ^ + 2Ar,^°"' '^ i ? e 4 A r ° " ' ' ^ + 2A,7 ° '''' = " i j i A k ? ^ " ^ " " ^ ' ' ^ ^ ' ^ ' - ^ + ^ ^ ^ M ^ 2 P . J ) + 5.,,-; (5.23) Re At ^ 2 A r ° "^^2Ar,^°"^ '^ i i e ^ A ^ ° + 2 A / ° " ' ' ^ ^ ^Re Jl^Ae^^-""" 2AiAn ' ^ Av^^^^-""''^ -{D',y,jDip,j - DiyijDlp,j) + S^^j, (5.24) with J.-,i4ACA7/' i = 2 , - - - , M - l ; i = 2 , - - - , i V - l . Special care was taken in the discretization of the pressure term to provide the necessary constraint needed on dVl'^. JD^ and D^ take the form of one-sided differences instead of the usual centered differences inside il'^. The discretized boundary conditions are: uij = 2(1 - ylj, vij = 0, DhMjDi^Mj - DiyMjDluMj = 0, DlyuMvMj - DlyujDlvuj = 0, j = l,---tN; Dlxi,iDiui,i - D^x,-,iD>,, 1 = 0, Ui.i = 0, u , - , = 0, u.-./v = 0, i = 2 , - - - , M - l . y=1 (wall) :u=v=0, x=0 (inlet) : u=2(1-y2) , v=0.. X=Xl (ouflow): B y=0 (centerline) :u=v=0, Figure 5.3 Geometry of the flow in a tube with a constriction and its boundary conditions using the primitive variable formulation. 5.3 S I M P L E M e t h o d 5.3.1 Governing equations The conservation forms of the momentum and continuity equations are appropriate here since the SIMPLE algorithm is basically a finite volume approach. Hence, the governing equation are written as: ^ ^ u t + {uu)x + {vu)y - ^[(|^«x)x + ( J^"v)v] = -P^ + Su-, (5.25) ^ t ^ t + {uv)x + {vv)y - ^[(^v.). + (;|^t'v)v] = ~Py + S.; (5.26) hyu). + {yv)y] := Sd. (5.27) Here u and v are axial and tangential velocity components, respectively. The expression for the source terms S^, Sv, and Sd are given in Table 5.3. 5.3.2 Representation of constriction It is possible to use the boundary-fitted coordinate technique to transfer the computation in an irregular geometry onto that on a regular domain (a square in the two-dimensional case) for a staggered grid. However, this has several disadvantages: (i) the formulation is much more complicated than on a nonstaggered grid; and (ii) relatively more computer storage is needed [66]. Hence it is not used here. By developing new schemes on a fully non-staggered grid, it is possible to implement the boundary-fitted coordinate technique in a much easier way. Thus the constriction here is handled by a simple approximation method. Patankar 67] introduced a so-called artificial viscosity method when he apphed the SIMPLE al-gorithm to a heat transfer problems involving different media, e.g. fluid and sohd walls. By setting up a mesh covering the whole domain of interest and assigning an artificial viscosity to the sohd walls, the computation can be performed through the whole domain including the sohd waU. The no-slip velocity condition at the sohd wall is a direct con-sequence of the large artificial viscosity of the solid wall. The same idea can be used to approximate the constriction here. The irregular geometry of the boundary is approxi-mated by assigning the viscosity inside the wall, in the constriction area, to be very large (Figure 5.5). Then the computation is performed in the straight tube. It should be mentioned here that a relatively smooth boundary contour was obtained using the boundary-fitted coordinate technique. The boundary approximated by the approach of artificial viscosity used here is not smooth, but represents a stepwise approx-imation. 5.3.3 - Boundary condition and discretization The boundary conditions are the physical ones, the same as in the previous section. There is no boundary value needed for the pressure correction [41]. The discretization of the governing equations is carried out on a staggered grid by using the finite volume method, similar to the M A C method [39]. The continuity equation is discretized on the control volumes, while the momentum equations are discretized between the control volumes. Instead of solving a Poisson equation for pressure as in the original M A C method, a pressure correction equation is derived from the continuity equation [41]. The discretized momentum equations can be written in a symbolic form as: acljUij = aeljUi+ij + awljUi_ij + anljUij+i + asl^Uij_i + slj, (5.28) = «e-j'^.+iJ + au ;V^u._i j -|- anl-Vij+i + a3V_,u,-,_,_i -f slj, (5.29) where the subscripts c, e, w, n, and s stand for the values at the center and neighboring east, west, north, and south control volumes, respectively, u and v are velocity compo-nents. In an iterative procedure, the velocity components and the pressure are assumed to be corrected in each step as u""*"^ = u" + u', v""*"^ — + v' and — p" + p' where ' values are the corrections. A pressure correction equation was derived from the conti-nuity equation (5.27) using a velocity pressure relation simplified from equations (5.28), (5.29) as Details of this discretization are given by Patankar [68]. Expressions for the coefficients and source terms are given in Appendix A . It is apparent that if the viscosity in the constricted region is set to an arbitrary large value, the entire region behaves as a sohd. For example, if a control volume with normal viscosity is on the west side of the control volume with large artificial viscosity, the magnitude of all the other coefficients will be small compared to acij and ae,-,j. Then, from equation (5.28), = Since the velocity inside the wall (u,+ij) is zero, a zero velocity was obtained from the calculation at the waU (u,j), which is the no-slip condition. Figure 5.5 Illustration of the artificial viscosity method using the SIMPLE algorithm. 5.4 Unsteady Flow Calculations The study of hemodynamics in arteries would require apphcation of the IC method to periodic flows. Highly complex nature of the in vivo conditions makes precise description the flow field virtually impossible. However, by using a judiciously simphfied model, useful information relevant to the problem of atherosclerosis can be obtained. 5.4.1 Governing equations and boundary conditions Since only the axisymmetric flow was explored, the stream function-vorticity formulation was used for the computation. The governing equations are the same as in the previous section, but the boundary conditions are now given as functions of time t. They are: • mv,t) = fi{v,t), MO,V,t) = /2(7/,t) at FA; • rPi^, 1,f) = fsit),1,t)MC,l,t) + xd(,l,t)xPr,U, 1,0 = 0 at CF; • tp{(, 0, t) = 0, CU, 0,t) = 0 at AB; and • y,(l,T},t)xPi{l,r,, t)+y^{l,rj, t)0,(l,rj, t) = y,{l, rj, t)Q{l,77, t)+y^{l,ri, t)Ul,t) = 0 at BC; where fi, /a , and /a are given. 5.4.2 T i m e discretization Usually, implicit schemes are preferred for time marching solutions of the Navier-Stokes equations because of their stability. The ADI method [61] aUows the construction of an efficient implicit scheme which leads to the solution of a tridiagonal algebraic system. Several contributions have been reported on the subject (Appendix B). Since the Navier-Stokes equations are nonhnear, the second-order accuracy in time cannot be obtained by applying the A D I method directly. Thus , an iterative procedure, using Crank-Nicholson for time discretization. Ois needed to ensure a second-order accurate solution both in time and space. Details are presented in Appendix B . 5.5 Preliminary Computation In this section, results of the prehminary computations to determine: the convergence rate; and the effect of location of the inlet and outlet boundaries, grid size, and different schemes are discussed. 5.5.1 Convergence rate analysis To ensure accuracy of the scheme, the first important test is the convergence rate analysis. As before, a typical problem was solved using four different grid sizes. The Reynolds number was chosen as 100 since a fine resolution can be obtained for a moderate stenosis (maximum stenosis height H = 0.5, stenosis length I = 4) on a relatively small domain (1 X 20). The grids are 21 x 51, 41 x 101, 61 x 151, and 81 x 201. The calculated convergence rates are listed in Table 5.4, and were found to be close to 2, theoretical predictions for second-order schemes. Values of the maximum wall vorticity, separation and reattachment points are given in Table 5.5 for later comparison. Two convergence rates for the unsteady flow were also calculated (Table 5.6). It can be seen that the second-order accuracy in time was achieved. 5.5.2 Location of the inlet and outlet boundaries Specification of the domain of computation is an important consideration in any nu-merical scheme. Obviously, too large a domain would involve unnecessarily extensive computational effort. On the other hand, too smaU a domain can affect accuracy of the solution in the region of interest. Experience suggested that the downstream location effect to be larger than that due to the inlet if it is not properly defined. Thus, a rela-tively short inlet length (i.n) and a long outlet region (Lout), as suggested by Aubert and Deville [55], appeared reasonable. Computation was carried out to assess both the inlet and outlet location effects. Two sensitive flow quantities, the maximum wall vorticity and the separation length, were chosen as the test standards. The results for typical Reynolds numbers and stenosis geometry are listed in Table 5.7. It can be seen that for the domain chosen {Lin = 8 and Lout = 56, in units of R, the upstream radius), any further increase in the inlet and outlet lengths have limited effect, even when the distance after the reattachment point is quite short (for Re = 1000), and the flow is far from fully developed. 5.5 .3 Effect of the grid size Mesh size is yet another important computational parameter which has to be determined to get relatively accurate results. Generally speaking, any desired accuracy can be ob-tained by reducing the grid size appropriately. However, there is a limit to the grid size because of the round-off error and computer time. Therefore, a compromise is needed. Several tests were carried out for a typical Reynolds number and stenosis geometry. For the second-order scheme, a grid size of 0.04 x 0.02 gave quite accurate results underpre-dicting separation length by 1.5% (Table 5.8). The first-order upwind scheme results in a relatively larger error with a 0.04 x 0.02 grid, under-predicting the separation length and peak waU vorticity by a maximum of 21.7%. The solution improved with a 0.02 x 0.01 grid giving a maximum error of 12.6%. However, the computational time increased by 4 times since the convergence is slow for first-order schemes. 7 Table 5.2 Source terms of primitive variables formulation. Su S ^ S d 2 1 2 1 V V 7 e 7 l e y ^ 7 ^ 7 Table 5.3 Source terms of SIMPLE algorithm. Su 2 1 1 y^^ Re Re ^ 0 Grid Number « X « Y « ; 21x51^1x101-61x151 1.7839 1.8781 1.9268 1.9419 21x51-41x101-81x201 1.7812 1.8811 1.9282 1.9456 Table 5.5 Grid dependence test on peak vorticity and separation length. Re Grid Cmax L s L r 100 21 x51 70.377 0.664 5.621 41 x 101 70.401 0.684 5.679 61 X 151 70.345 0.689 5.690 81 x201 70.317 0.691 5.694 Table 5.6 Convergence rates of the unsteady flow. Time Step at . "s min (Global in t, pointwisein space) (Global in space, pointwise in t) 0.05-0.025-0.0125 2.0084 2.1622 * Note, convergence rate for vorticity at the throat is calculated in a complete time cycle using discrete norm; + convergence rate for stream function and vorticity is calculated at an instantaneous time in the computation domain using space discrete norm. Re Inlet outlet Çmax L s L r length length 8 8 70.377 0.664 5.621 100 16 8 70.377 0.664 5.621 8 16 70.377 0.664 5.626 8 56 148.846 0.476 28.210 1000 16 56 148.637 0.477 28.238 8 80 148.638 0.478 28.233 Table 5.8 Grid dependence test on peak vorticity and separation length. Re Grid Scheme Cmax L s L r 2 6 x 8 1 1 69.200 0.757 4.672 100 51 X 161 1 69.615 0.740 5.115 2 6 x 8 1 2 69.970 0.681 5.611 Note: (i) scheme 1 represents the first-order upwind procedure, scheme 2 is the second-order center difference method; (ii). accurate solution for this problem was listed in Table 5.5 where a fine mesh (81x201) was used. Chapter 6 E X P E R I M E N T A L I N V E S T I G A T I O N A n experimental program was used in conjunction with Ctthe computational results to investigate the flow field Oover a common l imited range of parameters. Of course, the experimental results do not necessarily carry the stamp of truth. They can involve error depending upon the methodology and sensitivity of the instruments. However, wi th modern instrumentation and refined techniques, one can accept experimental information for the study in hand with a good measure of confidence. Experiments were carried out with steady flow for a simple stenosis model over a l imited range of the Reynolds number of interest in the physiological problem under study. The agreement between the numerical simulation and the experimental measurements, and with results by other investigators (when available) tend to substantiate both the methodology and the resulting data. A n interesting feature of the experimental investigation is to provide information which is very difficult to obtain numerically. For instance, finding the effect of distur-bances on the fiow quantities is indeed a challenging task for the numerical simulation at low Reynolds numbers (less than 1000). In fact, such disturbances in the flow exist at most locations of the arterial system, in the forms of branches and curvature. Complete information on this subject would require a thorough investigation by analyzing the full spectrum of velocity fluctuations. Although this is beyond the scope of the present thesis, some useful information concerning as the velocity profile, pressure recovery, and shear rate (shear stress) at the wall , help in appreciating limitations of the numerical results. W i t h this as background, this chapter discusses the experimental methodology used. The results wil l be compared with those of the numerical simulation in the next chapter. 6.1 Stenosed M o d e l The experimental model was constructed from plexigleis. The symmetric constriction followed a cosine curve with a 75% area reduction. The geometry of the model is shown in Figure 6.1. The stenosed model was machined to an accuracy of 10~^ cm. The cosine contour was purposely selected to facilitate comparison as it has been used before in both a theoretical [28] and in an experimental [35] study. It is one of the simplest representations of an arterial stenosis as it provides the smooth variation, which is found in many physiological cases. The same model also is used in the theoretical calculations. 6.2 Test Facility The facihty for the steady state measurements is the water tunnel designed and fabricated in the Department of Mechanical Engineering. It was originally designed by Aminzadeh [69] for his experimental study of hydraulic performance of the aortic prosthetic heart valves using enlarged models, and later modified extensively by Akutsu [70] to require-ments for his investigation in the same general area. The main criterion governing the design was the Reynolds number, which is based on the size of the model and the vis-cosity of the working fluid. Thus, the flowrate is treated as the parameter for reference. The flowrate range in which a steady flow situation can be estabhshed is from 10 cm^/s to 300 cm? I s as governed by the characteristics of the power unit . W i t h water as the working fluid and the model size of 5.08 cm in diameter, the Reynolds number range is around 200 to 8,000. The tunnel is shown schematically in Figure 6.2. It consists of a test circuit including three subassembhes: the test chambers, the fluid return system, and the power unit consisting of a pump and a drive motor. The large test-section is built of four plexiglas walls 2.44 m long, 1.905 cm thick and wide enough to produce an inside cross-section of 20.32 X 20.32 cm. The relatively long length was purposely chosen to ensure sufficient room for installation of flow distribution and straightening devices, and to permit the positioning of the model with ease. A vent, 10.16 cm in diameter and 30.5 cm high, located on the downstream end of the test section provided for fluid expansion as well as an escape route for the air bubbles. It also serves as an effective check against the over pressurization of the test section, particularly near the model location, irrespective of the pump's operating condition. There are five ports to access inside of the section: through each end; v ia two hatches; and through a porthole; located such that one's arm can reach, position, and adjust the model anywhere in the tunnel. In addition, several smaller portholes which could take 1.6 cm N - C plugs were drilled and tapped on the top of the plastic "box". Furthermore, 2.54 cm portholes are also provided on the side face of the "box". These openings were used to mount the model and convey pressure conducting lines. They were sealed off with plugs with "0" - r ings when not in use. Two glass plates, 63.5 X 13.97 X 1.27 cm, recess-mounted in the sides of the test section provided optically flat, homogeneous, and thermally stable wafls for inspection and photography. A drain positioned at the bottom of the "box" facilitated complete emptying and flush cleaning of the tunnel. The tunnel is provided with an alternate channel of flow consisting of a P V C settling chamber, a plexiglas test-section and support blocks. The settling chamber is made from a 25.4 cm diameter P V C tube with a leng*-.h of 104 cm. It is connected to the plexiglas tube test-section carrying the stenosis model. The system is constructed in a way which allows for easy changing of the stenosed model. In the present study, the inner diameter of the upstream and downstream tubes was selected to be 5.08 cm to allow good spatial resolution for measurement of veloc-ity profiles. The entrance length (the length prior to the stenosis) is 50 cm (10 times the diameter). The stenosis model with connecting tubes was installed inside the large test-section, described earlier, for velocity measurements using an L D A system and for pressure measurements using a capacitance type transducer (Barocel). Of critical impor-tance was the transition of the flow from the pump outlet to the test-section. For the large test-section where flow from the pump has to expand from a 7.62 cm diameter tube to the 20.32 x 20.32 cm box, the jet-like flow has to be diffused and spread evenly across the test cross-section. This was achieved by the following arrangement: • several sections of honeycombs to straighten the flow through turbulent exchange and laminar damping, • brass screens of different pore size with or without fiberglass wool in between. Located between the end of the test-section and the power drive unit is the return section essentially comprising a heat exchanger, P V C pipes and elbows with connecting flanges, flow control valves, and a radiator hose. A copper pipe, 350 cm long with 7.62 cm diameter in conjunction with a 244 cm long and 15.24 cm diameter P V C plastic sewer pipe formed an annular single-pass heat exchanger. With the coolant supplied by a water main it was possible to maintain the temperature of the working fluid within 0.2°C. P V C elbows and sections of the radiator hose provided relatively easy, anti-corrosion and vibration free connections between the test-section and the heat changer. The power unit consists of a centrifugal pump (Aurora type G A P B , 200 gal/min, 7.6 m head, 1750 rpm) driven by a three horsepower variable speed D.C motor. The pump impeUer and housing are of cast brass to guard against possible corrosion. The motor is energized by a three phase grid, the voltage being adjusted through an auto-transformer (Variac model 4711) and rectified by selenium diodes. No further smoothing of the D.C. output was reqmred. Flow rate in the tunnel was monitored using a sharp-edge orifice plate mounted up-stream of the pump inlet. The plate location was selected so as to make its reading rela-tively free from upstream and downstream disturbances in the form of elbows, change in section at the pump inlet, pump suction, etc.. A filter (10 micron) in the bypass circuit across the pump is used to minimize dirt contamination of the fluid. 6.3 Flow Measurement System - Instrumentation The physical quantities, such as the pressure distribution, shear stress along the wafl, and the scale of the separation region are of fundamental importance for the flow in stenosed arteries. They are known to have some correlation with the occurrence and the growth of the stenoses. Experimental evidence of these quantities will assist in the understanding flow field in the presence of stenosis. 6.3.1 Measurement of velocity Charting of the velocity field around the stenosis is not an easy task. Recognizing that the flow field is complicated because of separation and the irregular boundary contour, the flow measuring system should: • be able to monitor stagnant and negative flow; and • not perturb the local flow, i.e., it should be noninvasive. Hence, an L D A system appeared to be ideal to that end. A single component LDA system was employed for the measurement of the axial velocity component. L D A is a noncontact optical system for the measurement of fluid flow. It utihzes scattered light from particles in the flow field to measure local velocity of the scattering element. This relatively new technique is based on the invention of the gas laser in the early sixties, with its unique property of spatial and temporal coherence . Most current laser anemometers measure either the rate of change of frequency of the scattered light or the time flight between spatial regions of high intensity in the measuring volume. The former is called Laser Doppler Anemometer (LDA) and the later is called the Laser Transit Time Anemometer (LTA). L D A , being better suited for real time measurements, was preferred. Split laser beams from a single laser source are made to cross the flow field to form a measuring volume, the region of interference fringe pattern. This is accomplished using a beam splitter and a focusing lens. Fluid particles in the measuring volume at a given instant scatter light in all directions with frequency modulation. As a cross-section through the measuring volume consists of alternate light and dark regions, a particle passing through the fringe system emits light pulses at a frequency dependent upon its velocity. This is referred to as the beat or doppler frequency. The standard laser anemometer system has a 180 degree direction ambiguity. In other words, it can not distinguish between forward or backward flow. This would be a serious limitation in the present study. In a refined system the problem is overcome by a frequency shift between the split beams using a Bragg cell thus causing asymmetry in the interference fringe pattern. The use of an L D A for velocity measurement involves careful consideration of the system elements such as: • optics and optical arrangement; • seeding of the flow; • size and type of laser; • type of signal processor; • data acquisition and reduction system. These were considered in conjunction with anticipated characteristics of the flow field: • near zero and/or negative velocity; • complex anatomical shape; • continuous data. After considerable dehberation the following components were selected for the L D A sys-tem: • Ion Laser Technology Ar-Ion laser model 5490 AWC — 00 with 100 mW output, wavelength 5.30 x 10"^; • OEI L D A Transmitter Optic Module consisting of - LD-0- 0102 Optical adapter, - LD — 0 — 210 Polarization rotator module, - LD — 0 — 310 Beam splitter module, - LD — 0 — 420 Double bragg cefl module, - Z/JD - 0 - 610 Lens module / = 310 mm; • DISA 55 X 34 P M Receiving optics; • TSI Photomultiplier model 962 with model 965 power supply; • DISA 55L20 Signal processor consisting of - 55L30 Preamplifier, — 55X37 Frequency tracker, - 55L40 Meter unit. The laser unit is employed with an optics package which is in the forward scatter mode. The beam was split into two parallel beams at a distance of 48mm, then focused by a lens (/ = focal length = 310 mm). At the intersection of the beams an approximately ellipsoidal samphng volume is produced. The dimensions of the sampling volume were 1.3 X 0.3 mm for this investigation. Lax liquid solution was used as seeding particles for the measurement. Using a frequency tracker, the value of the shifted frequency fo is measured and hence the velocity component normal to the fringes can be determined by the following relationship, where A is the wavelength (= 5.30 x 10~^ m for this arrangement), 6 is the beam inter-section angle (9.13°), and V is the component of velocity in the plane of the beams and perpendicular to the bisector of the included angle 0. 6.3.2 Traversing mechanism In order to measure the velocity at different locations in the fluid field, the L D A system was supported on a platform with three degrees of freedom which permitted, within limits, any desired spatial positioning of the measuring volume. The traversing gear consisted of two subassemblies: a base traverse mechanism for rough positioning and a micrometer controlled X — Y translation stage for finer movement in a horizontal plane. The base traverse mechanisn^ was designed and fabricated in the Department of Me-chanical Engineering. It consisted of three platforms, two of which rode on a pair of hardened steel rods with linear bearings permitting movement in a horizontal plane, while the third rides on a modified mechanical jack permitting vertical movement. The maximum travel distance available was 86.6 x 23.7 cm along the horizontal plane and 22.8 cm in the vertical direction. All platform components were made from heavy alumimum plates with 13 mm thickness to minimize static deflection and vibration problems. The fine adjustment horizontal translation stage, supported by the base, was made by A.W. Becker GmbH of Germany. It is free to move 10 cm in each direction with an accuracy of 0.001 mm. Figure 6.3 shows the L D A system with its traverse mechanism. 6.3 .3 P r e s s u r e t r a n s d u c e r Because the pressure (both its absolute value and difference at different locations) is very small (order of 10~^ psi, 70 N/m^), a highly sensitive instrument was required for its measurement is demanded. This was accomplished using a "Barocel Modular Pressure Transducing System" developed by Datametrics Inc. of Waltham, Massachusetts. The type 550 — 5 Barocel sensor is designed to operate with fluid over the pressure range of 0 — 10 psi. The unit is a high precision, stable capacitive voltage divider with a variable element in the form of a thin prestressed steel diaphragm which deflects proportionally to the magnitude of the applied pressure. To isolate the external pressure medium from the sensor diaphragm-capacitance system, the unit uses highly sensitive metalhc bellows. The volume between the bellows, isolator and sensor diaphragm is filled with degassed silicon oil which serves both as a pressure transmitting fluid and as a dielectric. The pressure signal from the external liquid medium is transmitted by the bellows to the sihcon oil which in turn deflects the diaphragm to produce the required change in capacitance. An A . C . carrier voltage at 10 kHz is apphed to the stationary capacitor plates, and a bridge circuit determines an output voltage dependent on the ratio of the capacitance of the diaphragm to each of the stationary plates. The carrier voltage is therefore modulated according to the input pressure. The unit sensitivity is 10"* psi (0.007 N/m^) provided the pressure sensor is fully isolated from external sources of vibration and noise. It was imperative to ensure removal of aU tra^ ces of air pockets from the pressure conducting lines for satisfactory operation. The Barocel was accurately calibrated for steady pressure. The calibration plot is shown in Figure 6.4. 6.3 .4 Flowmeter Although the built-in sharp-edge orifice flowmeter in the water tunnel monitored the average flow rate, it did not give satisfactory results because the present study focused on relatively small flow rates. A more sensitive sharp edged orifice flow meter thus was carefully designed and installed into the flow system to meet the requirement. Before its final installation, the orifice plate was calibrated according to the standard A S M E procedure. Pressure difference across the orifice plate was recorded using a Barocel pressure transducer in conjunction with a Datametrics Manometer model 1018. The calibration plot for the orifice plate is shown in Figure 6.5. 6.4 Test Procedure The success of any experiment depends, to some extent, on the recognition of the ca-pacity and limitation of the equipment used. Performance specifications for instruments normally provide this information. However, they must be checked through actual op-eration. Furthermore, at times, it was necessary to push the equipment to its utmost capability to collect vital data. Hence extensive preliminary tests were carried out be-fore embarking upon an experiment to ensure accuracy (which implies reliability through repeatability) of results. Figure 6.6 shows schematically the arrangement of instrumentation, data acquisition and processing system used for the steady state measurements of pressure and velocity. Pressure distribution along the test-section was measured using a Barocel pressure transducer (Datametrics type 550-5) in conjunction with an electronic manometer (Data-metrics type 1018B). The Barocel is a differential pressure transducer wi th one side con-nected to the location of interest through its system of distribution values, while the other side is exposed to an upstream location. A c t u a l pressure tap locations are shown in Figure 6.1. A l l air bubbles from pressure conducting lines were carefully removed as they affect the results. This was achieved by opening al l the distributing valves to allow the conducting fluid to flow through the pressure conducting line wi th the tunnel run-ning for about 30 minutes. This also ensured a uniform concentration of water solution throughout the test system. As pressure levels were often very low, especiaUy at low Reynolds number, any change in viscosity or density of the working fluid resulted in a zero shift. To account for this, al l measurements were preceded by the monitoring of pressure at no flow condition. After removing aU traces of air bubbles, the distribution valve was opened for the measurement of the zero flow pressure at that tap location. Depending on the viscosity of the solution used, normally a few minutes of settling time was necessary for the pressure to be stabilized. The pressure was recorded and the valve closed. The procedure was repeated at afl the pressure tap locations. Then the tunnel was activated to estabhsh a desired steady flow rate as given by the orifice meter, and a new pressure signal measured. Velocity data were recorded using the L D A system described previously. The circu-lar tube with a constriction was placed inside the large square test-section of the water tunnel with the same simulation fluid surrounding it to minimize the effect of changes in refractive index between the tunnel material (plexiglas) and the solution. To further minimize the refraction effect due to curvature of the circular test-section, velocity mea-surements were conducted in the horizontal plane through the center of the test-section. Only the velocity component in the direction of flow was measured in this experiment. W i t h an appropriate flow rate estabhshed set, the laser beam was first focused inside the test tube wall, giving a constant frequency unaffected by the flow field, thus confirming the setting and working of the L D A system. Next, the measuring volume was moved to the interface between the wall and the fiow field using the traverse mechanism. It was possible to identify the interface accurately by monitoring a sudden change in the L D A signal. Having estabhshed the wall boundary, the measuring volume was progressively moved across the flow field. The signal from the L D A system was filtered prior to the data analysis to ehminate the high frequency noise. The mean velocity value was noted using a digital voltmeter with an appropriate time constant (DISA type 55£>31) and RMS value was recorded using an RMS voltmeter (TSI model 1060), and also calculated from the instantaneous velocity using an IBM personal computer. The average flow rate given by the orifice meter was also recorded. 6.5 Data Processor The analog velocity (from the frequency tracker) and pressure (from the electronic manometer) signals were converted into digital signals by an Analog to Digital Con-verter (ADC) before introducing to the IBM P C . Since the tracker frequency range was high (= 10 kHz) in all measurements, the frequency of the signal fluctuations at the measured Reynolds numbers was at least an order of magnitude lower (< 0.5 kHz), and the frequency of the environmental noise was much higher (> 10^° mHz). Thus there was no detectable error. Also, the environmental noise had no effects on the signal because its intensity was much lower than the laser beam (lower by a factor of 10 )^. In the calculation of mean and time dependent velocities, an average of 800 data block was employed, and usually 20 cycles were used to ensure the reliability. The velocity was expressed as u{^,t) = U{^) + u'{^,t), (6.2) where u is the instantaneous velocity, U is the average velocity, and u' is the velocity fluctuation. The mean pressure was calculated in a similar way. Usually 10 cycles were enough, however, 20 cycles were used when the low frequency fluctuations were large. Centerline Joint Figure 6.1 Stenosis model used in experimental measurements. c 5 Figure 6.2 A schematic diagram of the water tunnel used in experimental study. TE Laser Transmtter Optics Photo Multiplier Recieving Optic Support Beam a X - Y Translation Stage (fine pitch) X - Y - Z Translation Stage 1.6 Volt 1.4 ! 1.2 l i 1 [ 0.8 I 0.6 -0.4 -0.2 -0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 P (psi) Figure 6.4 Calibration curve for the Barocel pressure transducer. 1 F 0.1 b P (psi) 0.01 t 0.001 10 , 100 Q (ml/sec) Figure 6.5 Cahbration curve for the orifice flowmeter. Figure 6.6 A schematic diagram showing instrumentation, data acquisition, and pro-cessing systems. R E S U L T S F O R S T E A D Y S T A T E F L O W As can be expected, the amount of information obtained through a planned variation of system parameters is rather extensive. For conciseness, only results of direct relevance useful in estabhshing trends are presented here. The general approach to data reduction is discussed first, followed by numerical results for model stenoses in the physiological range of the Reynolds number. Results from the literature, when available, are also presented to facihtate comparison wi th those obtained through different approaches discussed earlier. F inal ly , the information is analyzed i n terms of its physiological significance. 7.1 Comment on Data Presentation Before proceeding with the presentation and discussion of results, a comment concerning the manner in which data in this general area of biofluiddynamics has been presented would be appropriate. The knowledge of complex fluid dynamics associated with stenosed arteries is indeed important for understanding the pathology of atherosclerosis. However, the related fluid dynamical parameters are rather numerous thus making it virtuaUy impossible to cover the entire spectrum of variables. As mentioned earlier, the flow quantities which are believed to have a close relationship with the atherosclerosis are the pressure drop (or pressure gradient) across the stenosis area, the maximum waU shear stress, and the separation parameters. Information on these quantities, and their dependence on the geometry of the stenosis, as well as the Reynolds number, wi l l be useful in both pathology and chnical diagnosis. For example, it is interesting to verify from the computation whether the 10 mmHg pressure drop across the stenosis, which causes significant changes in blood pressure level either as hypertension or hypotension, can exist under physiological condition. Obviously, the pressure and velocity distributions for the flow field are necessary to this end. Thus the focus is on the results, presented in a nondimensional form, which reveal influence of the constriction on the fluid flow field. 7.2 Effects of Flow Parameters on Flow Field The main purpose of this numerical investigation is to obtain a relatively complete knowl-edge of how the flow field is affected by the Reynolds number and the stenosis size, i.e., the relative height H and length L (in units of arterial radius). There are five basic stenosis models used in the computation as described in Table 7.1. The area reductions of M l , M2, and M3 are 44%, 56%, and 75% with same stenosis length. The models M4 and M5 have longer and shorter stenosis length compared to Af 3, with area reductions of 75%. As recommended by the World Health Organization [71], three grades of luminal narrowing were recognized as: • no stenosis: no reduction in the diameter (area) of the arterial lumen; • moderate stenosis: more than half the diameter of the original lumen remaining (area reduction less than 75%); • severe stenosis: less than half the original lumen remaining (area reduction more than 75%). It is further defined here that degree of stenoses with area reduction less than 54% is mild. The stenosis models selected in this study belongs to three groups, mild (Ml), moderate (M2), and severe (M3, M4, and M5). 7.2.1 Streamline and vorticity distribution A description of various steady laminar flow patterns which may be encountered is, perhaps, best rendered by a display of streamline and vorticity contours. Figures 7.1 to 7.5, present contours of constant stream function and vorticity for five stenosis models with Re = 100, 500, and 1000. The contours shown in Figure 7.1 is for a mild stenosis (Ml). The two sets of contours are presented in the same diagram: bottom half showing the vorticity and the top half the streamline plots. It can be seen that no separation occurs for Re=100 since no recirculation region is present and the vorticity value is positive in the whole region. The flow is mildly disturbed as it passes through the narrowing. Although asymmetry about the minimum constriction plane (x = 0) is not clearly noticeable in the streamline patterns, it is readily apparent in the vorticity contour. For this mild stenosis model, separation did not occur until the Reynolds number became relatively high. The size of the recirculation region remained very small even when reached a value of 1000. Note, it is not visible in Figure 7.16 for Re = 500, and an incipient tiny vortex can be anticipated in Figure 7.1c for Re = 1000. An increase in the Reynolds number tends to make the flow more asymmetric (in terms of the recovery length for fully developed flow after the constriction) about the x = 0 plane. The eff'ect of degree of stenosis (characterized by height, H) can be assessed through a comparison of Figures 7.1 - 7.3. Figure 7.2 shows contours of the streamlines and constant vorticity for a case of moderate stenosis, with larger area reduction (M2), for Re = 100, 500, and 1000. Corresponding results for severe stenosis (M3) are shown in Figure 7.3. It is apparent that now the separation starts relatively earlier, and the size of the recirculation region is significant even for Re = 100 in model M3 (Figure 7.3a). As the Reynolds number increases, the recirculation region extends a considerable distance downstream (Figure 7.36, 7.3c). The vorticity becomes more asymmetric and a high shear region develops in the constricted area as Re and H increase. Figures 7.4 and 7.5 present the contours of streamline and constant vorticity for stenoses models M4 and M5. The streamline and vorticity patterns as well as the size of the separation region are similar to those observed in Figures 7.3, for a longer stenosis length (Figure 7.4) and a shorter (Figure 7.5) stenosis lengths. The eifect of stenosis length, L, is not as significant as the its height, H. Hence the classification of severity of stenoses on their height [71] is reasonable since the stenosis length plays a less important role. (a) iîe = 100 \ V = 0.5 / V = o -10 (b) Re = 500 (c) Re = 1000 \ V = 0.5 Figure 7.1 Streamline and vorticity contours for model M l . (a) iîe = 100 (b) Re = 500 \ V = 0.5 y v = o (c) Re = 1000 Figure 7.2 Streamline and vorticity contours for model M2. (a) Re = 100 -10 (b) Re = 500 (c) Re = 1000 Figure 7.3 Streamline and vorticity contours for model M3. (a) Re = 100 -10 (b) Re = 500 (c) Re = 1000 Figure 7.4 Streamline and vorticity contours for model M4. (a) Re = 100 (b) Re = 500 \ C = 0 y C = 166.8 10 —r-15 (c) Re = 1000 = \ V = 0.5 E Z Z I Z I ^ y rP = 0 -5 1 0 1 1 5 10 15 2 0 •5 0 5 10 15 2 0 Figure 7.5 Streamline and vorticity contours for model M5. 7.2.2 Separation and reattachment points The locations of separation and reattachment points i , measured in the axial direction starting from the minimum constriction area, are of great physiological interest since they carry the information of the separation region. It is beheved that the low shear stress in the separation region is related to the mass transportation across the arterial wall [8] and long 'dwell time' for platelet endothelium interaction [12]. The occurrence of stagnation at the reattachment point is also believed to have a thrombogenic effect because of the collision between elements in blood and the arterial wall [14]. Figure 7.6 illustrates the effect of Reynolds number and stenosis size on these parameters. It can be seen that the separation length increases almost linearly with the Reynolds number for all the stenosis models. It is clearly demonstrated that the locations of the separation and reattachment points is sensitive to the stenosis height as well as the length of the stenosis. This is because that the occurrence of separation in the boundary layer de-pends on the shape of the boundary. However, the separation length measured between separation and reattachment points in the axial direction, highly depending on stenosis height (comparing M l , M2, and M3), is insensitive to the stenosis length (comparing MS, M4, and M5). In terms of the Reynolds number, the separation occurs quite early for the severe stenosis (M3, M4, and M5), though its locations are slightly different for the three models with the same area reduction and late for mild and moderate stenoses (Ml and M2) (Table 7.2). In summary, the severity of the stenosis (determined by H) has significant effect on the general features such as the separation Reynolds number and the size of the separation region, while the shape of the stenosis (characterized by L) has less effect. 7.2.3 Pressure drop across stenosis The pressure difference in the radial direction is very small, even in the constricted area (Figure 7.7). However, the varaition across the constriction in the axial direction, at given Reynolds numbers, is significant. The flow experiences an increased resistance due to the presence of the stenosis. Although the computational model is a highly simphfied representation of the vascular system, it clearly shows a tendency for a stenosis to reduce the supply of blood to certain regions because of the increased pressure drop. Otherwise extra work has to be done by the heart to provide the same amount of blood for the circulatory system, and proximal hypertension or distal hypotension is apparent. Thus, the pressure losses through such constrictions are of physiological interest. Figure 7.8, present the pressure distribution along the tube wall for all the stenosis models at Re = 100, 500 and 1000. It can be seen that the behavior is essentially similar. There is a rapid fall in pressure (with a reverse pressure gradient) until separation occurs. Pressure changes in the separation region are relatively small. The pressure will recover to Hagen-Poiseuille value provided the outfiow region is long enough for the flow to redevelop into the fully developed case. For the mild stenosis at a relatively low Reynolds number, the adverse pressure gradient will not be large enough to cause the flow to separate. Then a smooth transition to Hagen-Poiseuille flow will occur. To set up a standard for measuring the pressure loss. Young and Tsai [23] used the pressure difference between locations 8 diameters before and after the minimum constricted area as a reference. Although it is arbitrary, since the recovery region depends on the geometry of the stenosis and the Reynolds number, it is used here to measure the pressure loss. Its value for the stenosis models is given in Figure 7.9. The pressure drop, as expected, is larger than that of a tube without a constriction. The trends are as expected with larger pressure drops associated with the severe stenoses. (a) Re = 100 5 r-P -10 -15 h -20 h--25 - 1 0 (b) Re = 500 2 10 Poiseuille Flow 20 30 X 40 M l M2 M3 M4 -Q- M5 50 60 Figure 7.8 Pressure distribution along the wall. (c) Re = 1000 2 -10 h -12 Figure 7.8 Pressure distribution along the wall. 7.2.4 Wal l vorticity The wall shear stress is of physiological importance since there exist two contradictory hypotheses about its role in the initiation of atherosclerosis. The nondimensional wall vorticity is of considerable interest due to its relation to the value of wall shear stress. In the present computation, it is identical to the nondimensional shear stress (Appendix C). The experimental measurement of the wall vorticity is quite difficult, and no rehable method seems to be available. In this situation, the numerical simulation provides some insight into the level of the shear stress involved. Figure 7.10 illustrates the vorticity vari-ation along the sohd surface for a moderate stenosis (M2) at several Reynolds numbers. Values for the Reynolds number selected for display are 100, 500, and 1000, which cor-respond to no separation, incipient separation, and extended recirculation region cases, respectively. It can be seen that the location of the peak vorticity occurs just before the minimum constriction plane, and shifts upstream slightly as the Reynolds number in-creases. The vorticity value falls rapidly downstream of the minimum constriction plane, and soon approaches the same order as that of the Hagen-Poiseuille flow. It can also be seen from Figure 7.11, where the wall vorticity distribution is displayed for all the stenosis models at a given Reynolds number, Re=100, that a significant regions of recir-culation exists for severe stenosis (M3, M4, M5) as evidenced by negative values of the vorticity. A tiny separation region, not visible in Figure 7.10, appears for the moderate model (M2) at this Reynolds number with no separation for M l . Figure 7.12 presents the maximum values of wall vorticity for the stenosis models as functions of Reynolds number. The Hagen-Poiseuille value of 4 is shown for comparison. A rapid increase in the wall vorticity with the Reynolds number for severe stenosis is quite apparent. 60 -10 i ' ' 1 i 1 L - 4 - 2 0 2 ^ 4 6 8 10 Figure 7.10 Vorticity distribution along the wall for model M2. 100 Figure 7.11 Vorticity distribution along the wall at Re=100. Figure 7.12 Peak wall vorticity as a function of Re. 7.2.5 Flow velocity The velocity profiles are presented in Figures 7.13 - 7.15 for the three stenosis models at Re = 500. The regions of reverse flow are clearly seen for severe stenoses M3 and M5 due to the sign change in both the velocity components. The axial velocity profile becomes flatter when H increases and L decreases. Variation of the centerline velocity as affected by the stenoses geometry is presented in Figures 7.16 for Re = 500. It is apparent that the centerhne velocity becomes more asymmetric about the minimum constriction plane a; = 0 (as indicated by the recovery length to the fully developed flow) when the severity of the stenosis increases. For a given stenosis model, the same tendency appears when the Reynolds number increases (Figure 7.17). The velocity at and near the centerhne increases because of the reduction of cross area for the center flow. This reduction of center flow area is obviously due to the separation which occured within the constriction region at the downstream location. Severer stenoses cause larger separation region in both axial and radial dimensions, which reduces more center flow area. The flow redevelops after the separation with the boundary layer growing while velocity at and near the centerline decreases gradually. For a larger Reynolds number (Re — 500, 1000), it takes longer length for the flow reaches the fully-developed condition resulting a slow recovery for the centerline velocity. The maximum velocity moves slightly downstream from the throat since the flow continues to accelerate even after the throat. 7.3 Comparison with different approaches To compare (and. confirm) some of the conclusions, a few comparisons of selected quan-tities were made with results reported in the literature. Although it has been shown in previous chapters that the numerical method developed here gives accurate results, a comparative study also was conducted between the three numerical approaches dis-cussed before, i.e. the IC method for a primitive variable formulation, the IC method for a streamfunction-vorticity formulation, and the SIMPLE algorithm. 7.3.1 Pressure drop Experimental measurements of pressure drop have been reported by Young and Tsai [23 for the stenosis model M6 {H = \,L = 8), while numerical results for both models M5 and M6 are given by Deshpande et al. [28]. A comparison between the present numerical results and their data is presented in Figure 7.18. The agreement among the these results is qualitatively good. The SIMPLE algorithm over-predicts the pressure drop shghtly as expected, since the stepwise approximation to the constriction causes more pressure loss compared to the smooth geometry Deshpande et al. [28] also predicted an increase in the pressure loss compared to the results of the IC method for both M5 and M6 models. Considering that the present results were obtained with second-order accuracy for all the variables and first-order upwind scheme was used by Deshpande et al., it is possible that the artificial viscosity introduced by the first-order scheme dissipates more energy, and results in a large pressure drop. Other possibilities such as the effects of mesh size and treatment of boundary also may exist for this slight discrepancy in the results. (a) u-component -0.06 . 0 0.2 0.4 0.6 0.8 1 y (a) u-component -0.2 • ' I 1 1 I 0 0.2 0.4 0.6 0.8 1 y (a) u-component 5 I 7.3.2 Wal l vorticity A comparison also was made of the wall vorticity because of its crucial role in the phys-iological problem. Furthermore, an accurate determination of the wall vorticity presents one of the major difficulties in the numerical simulation. Thus its predicted values can be used to judge the performance of the numerical method. Figure 7.19 presents variation of the maximum wall vorticity as a function of the Reynolds number for the moderate stenosis model M6. Results obtained by both prim-itive variable and stream function-vorticity formulations are quite consistent. The SIM-P L E algorithm gives higher vorticity values due to the stepwise approximation of the stenosis. The results of a finite difference scheme by Deshpande et al. [28] and a finite element approach by Kawai et al. [65] also are given for the same stenosis model. Both agree weU with the present computation for small Reynolds number, while for relatively large Reynolds number, values predicted by Kawai et al. [65] are smaller than the present analysis, and Deshpande et al. [28] give even smaller values. It is difficult to judge whether the present method over-predicted or [28] and [65 under-predicted the wall vorticity, since no reliable experimental data are available. Table 7.3 presents the calculated values for the bell-shaped stenosis obtained by Lee and Fung [24] and the IC method, both have second-order accuracy. The computational results for a similar stenosis (M5) by Deshpande et al. [28] with first-order accuarcy also are listed for comparison. It is apparent that the present results are close to, but slightly lower than, those of Lee and Fung [24]. It appears that the IC method does not over-predict the wall vorticity. As pointed out before, it is likely that the first-order difference scheme used by Deshpande et al. [28] is responsible for the possible under-prediction, when the Re is large. It is well known that the artificial viscosity introduced for stability reasons has a low 'apparent' Reynolds number [72]. It is also possible that the rectangular mesh in [28] introduced some errors at the irregular boundary. 7.3.3 Locations of separation and reattachment points As pointed before, prediction of the separation and reattachment points is of physiological importance, and is another criterion to test performance of different numerical schemes. Figure 7.20 shows variation of the separation and reattachment points as a function of the Reynolds number for the moderate stenosis M6. It can be seen that the IC method and the calculations of Deshpande et al. [28] match rather well, however, both differ significantly from the experimental data [23]. Kawai et al. [65] predict a slightly later separation. The SIMPLE algorithm gave a much later separation than any other method. The agreement between both the IC methods is excellent, and they are also in good agreement with the results by Deshpande et al. [28] at small Reynolds numbers [Re < 500). Thus it appears that the numerical methods underpredict the separation length by giving an earlier reattachment point than that predicted by the experimental measurement [23]. This may be due to the flow instabihty as reported by Young and Tsai [23], or because of the unsteady character of the separation and reattachment points which make their accurate determination, experimental or numerical, quite difiicult. Slightly earlier reattachment was predicted by the IC method than the calculation of Deshpande et al. [28], but the results of the IC method are always in between those of Kawai et al. [65] and Deshpande et al. [28]. As accurate measurements are extremely difficult to obtain, further conclusions based on the existing results would not be appropriate. (a) Model M6 100 -T-Reynolds Number (b) Model M5 100 1000 Reynolds Number Figure 7.18 Comparison of pressure drop across the stenosis as obtained by several investigations. 70 60 -J 50 -40 J IC ( u - v ) — IC ( ^ ^ • SIMPLE i D e s h p a n d e et a l . X Kawai et a l . 30 -20 -10 -Poiseuille flow n r 10 100 Reynolds Number • -!—I—I—I—r 1000 Figure 7.19 Comparison of peak wall vorticity for model M6. 1200 1000 BOO h Separation points Re 600 r 400 200 -0 -Reattachment points • • SIMPLE Kawai et a l . D e s h p a n d e e l a l . ^ Y o u n g & Tsai 6 8 X 10 12 14 Figure 7.20 A comparison of separation and reattachment points as predicted by several investigators. 7.4 Experimental Results The model experiments to complement the numerical results were only carried out for selected quantities and characteristics of the flow field. A typical stenosis model of severe degree M3 = |, Z = 4) was chosen for comparison. Results were obtained for the pressure distribution along the sohd wall, axial velocity profile, and turbulent intensity. The experimental results were compared with the numerical predictions. 7.4.1 Pressure distribution The pressure distribution along the wall as predicted by the numerical IC method, and given by experiments are presented in Figure 7.21. The Reynolds numbers range covered was 250 - 1000. The correlation between the results, considering the complex character of the flow, is quite satisfactory for Re < 500. Even at higher Reynolds number of 1000, the experimental values agree well with those of the IC method in the convergent part of the stenosis and a few diameters downstream. However, for > 4, there is a large discrepancy in the pressure recovery value. Note, the experiments show sudden pressure recovery while the numerical prediction show a smooth trend. This can be expected as the numerical simulation is for purely laminar flows. The experimental flow field, however, contains a certain amount of turbulent fluctuations. The Reynolds number at the minimum constriction is twice as large as that upstream. Thus the presence of turbulent intensity, relatively larger than upstream values, can be anticipated as confirmed by the experimental data. The transition to or from turbulence could cause the pressure jump. It should be mentioned here that the pressure diff'erence is quite small in the experiment, approaching the limit of the instrumentation, especially for low Reynolds numbers {Re < 500). The measured pressure values fluctuated as shown in Figure 7.21. This fluctuation or uncertainty is expected to become smaller for larger Reynolds number (Figures 7.21c and 7.21d) since the absolute values of pressure in the experiments increase. 7.4.2 Velocity profile The velocity profile is of interest since it provides description of the flow field and has a direct relation to the shear rate (shear stress). Figures 7.22 shows axial velocity profile for a stenosis model of an artery as affected by the Reynolds number. In general, the trends depicted by the experimental and numerical results are similar, however, the local details do vary, particularly at the lower and upper ends of the Reynolds number range investigated {Re = 250,1000). The region of reverse flow is clearly seen. The experimental profiles tend to be flatter because of the turbulent fluctuations they contain. Figure 7.23 focuses on the experimental results on turbulent intensity obtained over a wider range of the Reynolds number {Re = 250—1000). As can be expected, the upstream turbulence intensity at x/R = —2 is rather smaU and virtuafly independent of the radial location. This is because the disturbances are smoothed out before the flow enters the tube. Primarily, the velocity fluctuations develop at the minimum constriction plane, x/R = 0, with larger values near the waU because of the shear layer. The fluctuation, progressively becoming smooth through diffusion, spreads along the radial direction, the profile becoming flatter when moving downstream. This general behavior was present at afl the Reynolds numbers used during the test program. The smoothing process is faster when the Reynolds number increases since the turbulent effect helps to mix and dissipate more energy. This turbulence fluctuation is responsible for the flatter velocity profile measured in the experiment compared to the numerical results. It also causes the discrepancy of the pressure recovery values between the experimental measurements and the numerical data as discussed in the previous section. The present measurements are compared with those of Ahmed and Giddens [35] for the same model, in Figure 7.24 for the Reynolds number of 500 and 2000. The agreement is reasonable, although the comparison is limited to a few locations. (a) Re = 250 _j4 i 1 1 I i I I I I -4 - 2 0 2 4 6 8 10 12 X (b) Re = 500 -12 ^ ' ^ ' ' ' ' -4 -2 0 2 4 6 8 10 12 X Figure 7.21 A comparison between numerically and experimentally obtained pressure distribution along the wall. (c) Re = 750 0 -2 -4 P -6 -8 h -10 -12 -4 - 2 (d) /^ e = 1000 0 -2 h -4 P -6 -8 -10 h 4 X C o m p u l a t i o n (IC) E x p e r i m e n t 10 12 C o m p u t a t i o n (IC) • E x p e r i m e n t -12 -4 -2 0 2 4 6 X 10 12 Figure 7.21 A comparison between numerically and experimentally obtained distribution along the wall. (a) Re = 250 x = - 2 x=0 x=2 x=4 x = 6 0.5 y 0 -0.5 -1 (b) Re = 0.5 y 0 -1 •\ - B^ ^ • \ • • • • ' < 1 !• \. • N. • \ • V a \^ • \^ \ « • \ ^ • \ ^ V A " \ m . • ' •\ • • • • / — • / • / « / • / • / B • r 7 •/ • / • / » / * / * / • / . • / \ \ " 0 1 2 3 4 5 Experiment Computation (IC) ) - 2 x=0 x=2 x=4 x = 6 \ \ \ • \ • \ • 1 • V \" 3 • B J i I 1 • • j m ' / / / • ' / ' / / m / u / • ' f 1 j 0 1 2 3 4 5 Experiment Computation (IC) (c) Re = 750 x = - 2 x=0 x=2 x=4 x = 6 > • • • 1 • r r • !•. •\ A A • • • s B • • • • • 1 */ • / ^ • y"^ • • / /• \ \ I I \ \ L 0 1 2 3 4 Experiment Computation (IC) (d) Rt = 1000 x = - 2 x = 0 x = 2 x = 4 x = 6 • Experiment Computation (IC) Chapter 7. RESULTS FOR STEADY STATE FLOW (a) Re = 250 0.5 h y 0 -0.5 -1 - 2 x=0 x= =2 x=4 x = 6 + f :+ + + + + + + + + + -+ -f- + + t + + • + + + + : + + + + + - i -• i -- l + + 4-+ + + + + + ++ + -+ + + : + + + + : + + + : + : + + + + + + 1 1 (b) Re = 500 y 0 0 .2 .4 .6 .8 1 + RMS 0 .2 .4 .6 .8 1 + RMS (c) Re = 750 x = - 2 1 0.5 y 0 -0 .5 -1 + + (d) Re = 1000 y 0 x=0 x=2 x=4 + + + + + + + + x=6 + + ; + : + + I L J ! 0 .2 .4 .6 .8 1 + RMS 0 .2 .4 .6 .8 1 + RMS (a) Re = 500 y 0 -0 .5 h (b) Re = 2000 0.5 h y 0 -0 .5 f u Present 0 1 2 3 4 5 Ahmed &: Giddens 0 1 Z 3 4 5 Present + Ahmed & Giddens Figure 7.24 Comparison of the present experimental data for axial velocity profile with those by Ahmed and Giddens [35 . 7,5 Discussion The effect of the steady state Reynolds number Re and the smooth constriction parame-ters H and L on the variables associated with the flow through a rigid tube was studied in the present chapter. This include quantitative predictions of the wall shear stress distribution, (and its peak value and location), pressure distribution (and its loss across the constriction), and velocity field. As can be expected, a steep stenosis contour, i.e. X —* 0 or —* 1, is more likely to induce separation, large shear stress, and pressure loss. Deshpande et ai [28], Kawai et aL [65], and others have contributed to the under-standing of this problem through numerical results for the steady state. These studies, however, are incomplete since they cover the physiological range of Re only for a mild stenosis model. The Reynolds number range used for the severe stenosis models was low. Furthermore, the method used in [28] to handle the boundary of a smooth constriction was not quite accurate. The present study covers the Reynolds number range of physiological importance for all the stenosis models. Of course, there is an inherent limitation, presented in any steady state calculation: the unsteady terms are absent. Pulsatile effects, which can be important in the flow through large arteries, have been neglected. In spite of this deficiency, the steady state results presented here can play a useful role. Although they cannot be extrapolated directly to the real fiow in the vascular system, the results do suggest the inadequacy of several theories aimed at initiation of atherosclerosis. It is useful to relate the present results to the experimental observations with coarc-tation in the monkey-aorta. Zarins et al. [31,73] demonstrated that when monkeys with a coarctation in the descending thoracic aorta are placed on an atherogenic diet, the de-velopment of lesions is significantly altered. The coarctation passage itself is consistently spared of disease development while the proximal aorta becomes heavily involved with surface lesions. Distal (downstream) to the stenosis, the degree of coarctation appears to affect the disease severity [74 . Since a high wall shear stress has been postulated as a possible initiating factor in atherosclerosis [9], it is interesting to examine the values by scahng the present nondi-mensional results to conditions approximating the monkey-aorta (Table 7.4). Values of the maximum waU shear stress are tabulated in Table 7.5 for the five models. The values for fully developed Poiseuille flow are also listed in Table 7.5. It can be seen that the maximum waU shear stress, which occurs slightly before the minimum constriction plane, is greatly increased by reducing the lumen area. It exceeds the value of 380 dyn • cm~'^ quoted to cause acute endothelial disruption [9] with severe stenosis for Re = 500 - 1000. However, in the animal models studied by Zarins et al. [31] and Bomberger et al. [32], no lesion development occurred in the constriction area. In contrast, extensive lesions were found to develop in either the proximal (upstream) or distal regions. In many of those animals, the degree of stenosis exceeded 75% with no diminution of blood flow, so that the shear stress values for 75% constriction (models M3, M4, and M5 in the present computations) should represent conservative estimates. Roach and Smith [75], arguing in favour of the high shear stress hypothesis, suggested that a relatively increased shear stress level (not necessarily as high as the critical value suggested by Fry [9]) for a long duration will also damage the endothelial cell by changing the permeability. However, the present numerical data, with some experimental evidences 35] and the animal experiments of Bomberger et al. [32,74] (which was dismissed by Roach and Smith [75]), certainly make the role of the high shear stress in atherosclerosis very questionable. While the argument about the significance of high shear stress continues, the hypoth-esis of flow separation, proposed by Caro and associates [8], has found more acceptance. The idea of mass exchange between the arterial wall and blood proposed by Caro et al. [11] is of interest. The interactions between the local flow fields have been supported by considerable evidence [13,14]. However, more quantitative experiments as weU as numerical simulations are needed to arrive at a final conclusion. An ingenious experiment by Flaherty et al. [76] showed that endothelial ceU nuclei were oriented along the predicted streamlines. Largille and Adamson [77] have also shown that the endothelial ceUs as well as the nuclei are oriented along the streamlines. Legg and Cow [78] have demonstrated that endothelial cells become disoriented in the direction immediately distal to the site of coarctations in the rabbit aorta. Cornhill et al. [79] presented the data which imply that endothelial ceUs tend to ahgn their longitudinal axis in the direction of the local flow. It is, therefore, tempting to speculate that the properties of the wall itself may be affected by the local flow. The disorder in the endothelial cell orientation could be associated with reduction of junctional tenacity between cells, while unidirectional high shear stress may elicit formation of tighter intercellular attachment, thereby reducing transjunction permeability [80 . The change in flow direction near the reattachment and separation points cannot be predicted by steady state models. However, it suggests that the extensive recirculation zone in the distal region of a moderate stenosis, at relative high Reynolds numbers, can be related to the findings from the animal experiments [31,32]. On the other hand, it appears that lesions which develop in the proximal region for a severe stenosis cannot be explained by the steady flow computation. Hypertension is often associated with an increase in the extent and severity of atheroscle-rosis, mostly in locations already predisposed to plaque formation [81] - [86]. Lesions also can be associated with elevated blood pressure in some locations usually spared of the disease. Sako [87] found that lesions developed when the pressure difference across an im-plemented coarctation was over 10 mmHg, while there was no significant atherosclerosis lesion for pressure differences below 10 mmHg. On the other hand, effect of hypoten-sion merits further discussion since reduction of blood pressure is considered to be an appropriate therapeutic approach to atherosclerotic disease in hypertensive patients [32 . It was found that for moderate stenoses, both proximal and distal region were involved with atherosclerosis lesions. For severe stenoses with a pressure drop over 10 mmHg, there was no significant effect on lesions appearing proximal to the coarctation, while lesions distal to the coarctation were strikingly less prominent [32]. This is consistent with other experiments [88]. By calculating the pressure drop across the stenosis in the five models investigated numerically corresponding to the experiments [32] (Table 7.6), it is interesting to find that for moderate and mild stenoses at aU Reynolds numbers, the pressure drop across the stenosis is quite small, below 10 mmHg. The pressure drop does exceed 10 mmHg for severe stenoses (M3, M4, and M5) at Re = 1000. The potentiating effect of hypertension on atherogenesis has been attributed physio-logical changes of the arterial wall such as permeability [82]. The effects of hypotension on aortic intimai permeability or intimai fibrocellular composition is not know [32]. In spite of experimental findings, the mechanism of hypotension with 10 mmHg in pressure drop remains unclear. The present steady state study has provided considerable numerical evidence, which can be used in conjunction with the experimental data, to better understanding of the flow across a stenosis and its consequence. However, the numerical results are for steady, axisymmetric laminar ilow, and direct extrapolation to flow in the circulatory system is not possible as the physiological flow is far more complicated. Even so, this simple steady state model is useful and can serve as a small step in improving the understanding of the phenomena. Chapter 7. RESULTS FOR STEADY STATE FLOW 163 Table 7.1 Geometry of stenosis models used in computation. M1 M 2 M 3 M 4 M 5 H 1/4 1/3 1/2 1/2 1/2 L 4 4 4 8 2 Area Reduction (%) 44 56 75 75 75 Table 7.2 Separation Reynolds numbers for five stenosis models M l M 2 M 3 M4 M 5 Re Sep 309.5 113 25 44 16 Table 7.3 Comparison of peak wall vorticity for model M2. Bell -shaped Sinusoidal Re c max Lee & Fung Present Deshpande et al. Present 10 25 100 300 57.2 72.7 54.8 70.6 49.5 59.7 84.9 105.4 49.5 62.5 94.2 147.4 Chapter 7. RESULTS FOR STEADY STATE FLOW 164 Table 7.4 Physiological quantities in the animal experiments. R (cm) p (g/cm^) y (stokes) 0.35 1.06 0.035 Table 7.5 Peak wall shear stress (dyn/cm ^) corresponding to the animal experiments. Re M l M2 M3 M4 M5 Poiseuille 100 8.3 13.2 36.7 28.6 49.9 2.1 500 67.1 110.6 316.9 240.3 423.2 10.6 1000 168.3 279.2 788.1 607.7 1022.7 21.2 Table 7.6 Calculated pressure drop (mmHg) corresponding to the animal experiments. Re M l M2 M3 M4 M5 100 0.09 0.12 0.24 0.22 0.29 500 0.60 1.07 3.80 3.78 5.16 1000 1.60 3.17 14.67 15.02 15.42 Chapter 8 R E S U L T S O F P U L S A T I L E F L O W S I M U L A T I O N Pulsatile nature of the blood flow and some unexplained observations [32] of in vivo animal experiments suggest that the unsteady flow simulation is essential to iUustrate the hemodynantiic role in atherogensis. Although importance of the unsteady flow characters has been recognized, only a few results are available. This is because unsteady experiment are relatively difficult to perform. Corresponding numerical simulation is also equally challenging. As a first step towards simulation of the pulsatile periodic flow, a simple harmonic component was added to a steady mean flow as f{t) = 1 + esin(2Trt) where t is the nondimensional time. The dynamic factor e is the ratio of the unsteady contribution (for example, the flux) to the mean value. There are four important parameters used to describe pulsatile flow in stenosed arteries. They are: the mean Reynolds number Re; frequency Reynolds number Rw; unsteady ratio e; and the stenosis height H. The present study covers the following range of parameters: 10 < Re < 1000, 0 < Rw < 15, 0 < e < 1.5, and 0 < H < l,; A complete periodic wave form for the flux from in vivo measurements would be more realistic, and could be used [33]. Unfortunately, such information is limited only to specific locations. Because the effects of stenosis geometry and mean Reynolds number have been discussed in the previous chapter for steady state flows, the focus in this chapter is on other two parameters, Rw and e. The information needed to fully describe the unsteady flow fleld is considerable, even for one simple case. It is difficult to decide on the type of information to display to illustrate significant features of the flow field [34]. Here time histories of the flow quantities such as the separation and reattachment points, and wall vorticity are selected for presentation at specified locations. The choice was influenced by the indication that the wall shear stress [9,8 and the separation region [89] are important factors in the development of atherosclerosis under pulsatile flow conditions. Figure 8.1 presents numerical results of the present study as well as the wall vorticity measurements [36] for the stenosis model M3 when Re = 600, Rw = 7.5, and e = 0.67. The agreement between the results is good at the minimum constriction plane and on the upstream locations, i.e. x < 0. However, at a downstream position, the numerical analysis predicts an oscillating back-flow, perhaps due to the coarse numerical grid. The focus is on presentation of the detailed flow-fields at certain critical time instances to assess how the wall vorticity and the separation zones are affected by the two pulsatile parameters, Rw and e. (a) t = 0.25 500 400 h 300 r ^ 200 -100 h -100 (b) t = 0.5 200 150 100 -<f 50 --50 L -100 -4 - 2 0 P r e s e n t " A h m e d & G i d d e n s 4 6 X 8 10 12 14 Figure 8.1 Comparison of wall vorticity distribution. 8.1 Streamline and Vorticity The computed time dependent flow field revealed several interesting features of arterial flows. Consider first fithe steady state streamline and vorticity contours for a mild stenosis MO {H = \, L = 2) at Re = 100, Rw = 10, and e = 1 (Figure 8.2). This stenosis model was chosen as a reference since there is no separation in the steady state. This makes it easier to identify effects of the unsteady component on separation as well as other flow phenomena. As shown in Figures 8.3 and 8.4, the instantaneous streamline and vorticity patterns for a simple pulsatile flow are considerably different from those of the steady state flow throughout the cycle. At i = 0 (Figure 8.3a), where the instantaneous Reynolds number equals the mean one during acceleration of the total flux, the streamlines are not very different from those in the steady flow. However, the vorticity contours show a significant difference (Figure 8.4a). Note, separation does not occur at i = 0.25 (Figure 8.3c) when the flow reached its peak forward value, but a small vortex appeared near the wall distal to the stenosis at t = 0.375 (Figure 8.3(i), with separation and reattachment points spaced close to each other. The vortex develops during the deceleration, when the separation and reattachment points move upstream and downstream, respectively. A high shear region develops at the same time at the proximal side of the stenosis, slightly before the throat (Figure 8.4c?). At ^ = 0.5 (Figure 8.3e), when the instantaneous and mean Re again are equal, there is a large separation region distal to the stenosis. A reverse flow also occurs proximal to the stenosis in the upstream region. A second vortex appears in the distal area, when the reattachment point moves to a downstream location. The vortices start to move away from the wall as time increases. At t = 0.625 (Figure 8.3/), the reattachment points (both downstream and upstream) are swept off the wafl. The backward flow occupies the near wafl region where a free shear layer develops. Large recirculation zones occupy both distal and proximal regions to the stenosis, and a vortex ring forms (Figure 8.4/) corresponding to the detached vortex in the stream function contour (Figure 8.3/). At f = 0.75 (Figure 8.3g), where the flow experiences zero flux at its lower peak, the recirculating flow occupies most of the tube volume. The vortices then move back to the wall as the acceleration period starts. The recirculation region disappears in the proximal region while that in the distal region continues to decease in size at < = 0.875 (Figure 8.3^), and the flow becomes complicated near the wall. By the end of the period, t = 1 (Figure 8.3a), the vortex is smoothed out and the cycle starts again. The instantaneous stream function contours for the same stenosis model with a smaller pulsatile component (e = 0.1) are shown in Figure 8.5. In this case, the effect of the unsteady flow is almost negligible, and the flow is similar to the steady case. The vortex near the wall is so small that it cannot be distinguished ftom the stream function contours. No separation occurred on the proximal side during the entire cycle. Figure 8.6 presents the instantaneous stream function contours for the same stenosis and at the same Re but for a reduced value of Rw {Rw = 3.34). The ratio of the pulsatile component to the mean value is 1. Note, a smaller value of the frequency Reynolds number corresponds to a longer period. Compared to the flow with Rw = 10, it can been seen that now the separation starts later and ends earlier in the cycle. The recirculation zone tends to be smaller in both the proximal and distal regions, with the structure in the distal region similar to that for Rw = 10. Generally speaking, the pulsatile effects are reduced with a reduced Rw. Although we focused on the effects of Rw and e, it also is of interest to have some appreciation as to how the mean Reynolds number. Re, and stenosis height, H, affect the flow field with a pulsatile component. Figure 8.7 shows the instantaneous stream function contours for the same stenosis model at a higher Reynolds number {Re = 600), with other parameters the same as in Figure 8.3. It can be seen that separation starts much earlier (t < 0.125, Figure 8.76), and the separation region is larger and flatter. The structure of the vortex downstream is clear during the deceleration period {t = 0.625, 0.75 Figures 8.7/ and 8.751). The recirculation zones persist longer, and there are vortices in both distal and proximal regions even at ^ = 0.875 (Figure 8.7h). The instantaneous streamlines for a more severe stenosis model M5 are shown in Figure 8.8 for the same flow conditions as in Figure 8.3. The structure now is more complicated. At i = 0 (Figures 8.8a), there is no separation near the stenosis, but a vortex is clearly seen in the distal region far from the stenosis. This vortex is smoothed out when the flow begins to accelerate, but now another vortex emerges no the downstream face of the stenosis {t = 0.0625, Figure 8.86). This vortex grows, and the structure becomes more complicated during both the acceleration and deceleration periods (t = 0.125 to t = 0.375, Figures 8.8c - g). When the instantaneous Reynolds number equals the mean value during the deceleration (i = 0.5, Figure 8.8/1), reverse flow occurs in the proximal area near the wall. The vortices grow and move inwards from the wall as the flow decelerates, with a large recirculation region appearing on both sides, of the stenosis {t = 0.625, Figure 8.8/1). The separation points on both sides are stiU visible at this time, move towards the throat, and persist until the lower peak flow is reached {t = 0.75, Figure 8.8n). Again, the vortices move back to the wall when the flow accelerates and the back-flow in the proximal region disappears, while several recirculation zones stiU exist in the distal region at t = 0.875 (Figures 8.8p, 8.8c). The vortices continue to decrease until t = 1, (Figure 8.8a), when Re equals its mean value and another cycle starts. The flow picture with a more severe stenosis is similar to that for a moderate stenosis model. There are, however, two important differences. Firstly, as expected, the structure of the recirculation zones is much more complicated; secondly, the reversed flow persists much longer both in the proximal and distal regions of the stenosis. -10 0 5 10 15 20 Figure 8.2 Stream function (top) and vorticity contour (bottom) of steady flow for MO at Re = 100. (d) t = 0.375 s 1 1 -I 1 1 i •10 -5 0 5 10 15 2 0 (e) t = 0.5 (f) t = 0.625 (g) t = 0.75 (h) t = 0.875 -10 (a) t = 0 and 1 (b) t = 0.125 (c) t = 0.25 (d) t = 0.375 15 2 0 (e) t = 0.5 (f) t = 0.625 (g) t = 0.75 (h) t = 0.875 (a) f = 0 and 1 -10 —r--5 \ V = 0.5: 10 — I — 15 — t 2 0 (b) t = 0.125 •10 -5 \ V = 0.54; y xp = o —r-10 15 —f 20 (c) t = 0.25 r -10 —r--5 \ = 0.55: y 'i> = o —r-10 15 20 (d) t = 0.375 -10 -5 \ V = 0.54: 10 15 — t 2 0 (e) t = 0.5 -10 -5 (f) t = 0.625 -10 (g) t = 0.75 (h) t = 0.875 —r--5 -10 -5 •10 -5 \ V = 0.5-- 1 — 10 \ V = 0.461 / V = o - T — 10 \ V = 0.45-y rP = 0 \ ip = 0.46 m —f— 10 —r-15 — T — 15 10 15 — I — 15 20 —t 2 0 2 0 2 0 {&)t = Q and 1 - 10 —r--5 \ xP = 0.5 / i> = o —Tr-io 15 —t 2 0 (b) t = 0.125 •10 -5 \ = 0.85 10 15 2 0 (c) t = 0.25 -10 -5 (d) t = 0.375 -10 -5 \ V = i y rP = 0 10 15 \ •tp = 0.851 - T -10 2 0 15 2 0 (e) t = 0.5 •10 (f) t = 0.625 1— 10 (g) t = 0.75 (h) t = 0.875 1— •10 - r --5 \ V = 0.5 —r-10 —r-15 \ V = 0.15 •vv. - .———• V V' = o-10 \ V = 0.15: ^ V = 0 —r--5 10 —T— 15 —t 20 —r 20 15 2 0 (a) f = 0 and 1 (b) t = 0.125 -10 -5 ^ V = o —r-10 15 20 (e) t = 0.5 (f) t = 0.625 1— 10 (g) t = 0.75 (h) t = 0.875 —r--5 ^ \ V- = 0.15 10 15 2 0 -10 {a)t = 0 and 1 (b) t = 0.0625 (c) t = 0.125 (d) f = 0.1875 (e) t = 0.25 (f) t = 0.3125 -10 \ V = 0.96" / V' = o 10 20 (g) t = 0.375 (h) t = 0.4375 -10 \ v = o . 6 9 : : ^ ^ V = o 10 2 0 •10 0 10 2 0 (q) t = 0.875 (r) t = 0.9375 -10 0 10 Figure 8.8 Stream function contour in a cycle model M5 at Re = 100, Rw = 10, e = 1. 8.2 Separation and Reattachment Points Although the streamline pictures show much of the details, the attached wall vortex length behavior can be presented more clearly graphically by considering the wall zero vorticity contours as functions of time. The locations of zero wall vorticity (or shear stress) are the stagnation points as well as separation and reattachment points of the attached vortices. Figure 8.9 shows the time history of the separation point (S curve) and reattachment point (R curve) in the {t, x) plane for the moderate stenosis model, MO, where Re = 100, Rw = 10, e = 0.1. The growth and decrease of the distal vortex from its initial appearance Pi{t,x) through the acceleration and deceleration period is clearly shown. When the pulsatile component increases from e = 0.1 to e = 0.5, the reattachment point is swept off the wall for a certain time period and a recirculation zone appears in the proximal region as the flow experiences deceleration (Figure 8.10). The separation point distal to the stenosis stays on wall, and there is no connection between the recirculation regions upstream and downstream. The reattachment point gradually moves upstream until it joins the separation point finally at P2{t,x). For a large pulsatile component of e = 1, separation occurs earlier around the peak flux {t = 0.25) as shown in Figure 8.11. The reattachment point is swept off the wall in the distal regions earlier. The two recirculation regions join together at P2{t, x), and are separated again at P3{t, x) during the latter part of the deceleration and the earlier part of the the acceleration phases, respectively. No significant difference appeared when the ratio of the pulsatile component increased to e = 1.5, where a negative flux exists during the latter part of the deceleration and early part of the acceleration cycles. The structure is basically the same as that for e = 1 (Figure 8.12). When the frequency Reynolds number was reduced to Rw = 7.5 and 3.34, with a pulsatile component ratio e = 1 (Figures 8.13 and 8.14), there was a tendency for the reattachment points on both sides of the stenosis to be swept off the waU more gradually and latter during the deceleration. The appearance of separation points in the earlier part of the acceleration cycle on both sides tends to be more abrupt at smaller Rw. The time history of the separation and reattachment points for a more severe stenosis model (Figure 8.15) was different compared to that for the moderate stenosis model. The flow field is quite complex since more attached vortices are found in the latter part of the deceleration and in the earlier part of the acceleration phases [P4{t,x) and Ps{t,x)]. As observed from the streamline contours, several vortices develop and move inward from the wafl between these two time stages. 8.3 Wal l Vorticity Distribution Because of its physiological interest, the wall vorticity is discussed in some detail here. The time history of the maximum average wall vorticity is shown in Figure 8.16. The maximum time-averaged wall shear stress also occurs within the stenosed region, and its is larger than for the steady state flow. The instantaneous value experiences a sinusoidal time variation due to the sinusoidal variation of the flux. Increase in the peak value with an increase in the pulsatile component ratio e, frequency Reynolds number Rw, stenosis height H, and mean Reynolds number Re are clearly shown Figures 8.17 - 8.20. They show that periodic variation is persisted and the peak value increases as the flow becomes more unsteady, i.e. Rw and e increase, the stenosis becomes severer, and flow with larger mean Reynolds number as well. The distributions of wall vorticity in the axial direction at several time instants are illustrated in Figure 8.21 - 8.25. Figure 8.21 shows that the difference from the steady state solution is not significant for e = 0.1. However, considerable deviation appears when e increases to 1 and 1.5 as shown in Figures 8.22 and 8.23. It is clear from Figure 8.22 that at i = 0.5, when the instantaneous Re equals its mean value, the distribution of wall vorticity is negative except in the constriction area, indicating that two large recirculation regions exist at both proximal and distal sides of the stenosis at that instant. At t = 1 during the acceleration, no negative value is found. A similar situation is observed in the case for e = 1.5, as shown in Figures 8.23. Figure 8.24 shows the instantaneous vorticity distribution for a smaller Rw. Although the vorticity value is negative at t = 0.75, its distribution along the wall behaves similar as the steady flow at i = 0.25, 0.5, and 1. This is because smaller Rw means larger time period, thus flow is less unsteady. For a more severe stenosis, in Figure 8.25, there is no significant qualitative difference except that a large variation was found in the recirculation region distal to the stenosis region 'for t = 0.25. Figure 8.26 compared to Figure 8.22 shows that the mean Re has little effect on the wall vorticity behavior except for the size of the recirculation region. To address the effects of frequency Reynolds number, the ratio of the pulsatile com-ponent to its mean value, mean Reynolds number, and stenosis severity, on wall vorticity, the distribution of wall vorticity is plotted at i = 0.25 when the flow reaches it peak value and starts to decelerate. Figure 8.27 gives the wall vorticity variation as a function of e. It shows that the length of the recirculation region and the peak value vary monotonically for fixed Re and Rw. Variation of the wall vorticity with Rw and H are shown in Figures 8.28 and 8.29, respectively. Note, the size of the recirculation region and peak value of wall vorticity also change mononically with Rw, except that the dependence is not as obvious as for e. Figure 8.30 shows the effect of Re on the wall vorticity at the same time instant oî t = 0.25. As in the case of the steady state flow, the unsteady flow shows that a large Re result in a large recirculation regions and peak vorticity values. The severity of the stenosis produces a larger separation region and peak wall vorticity value. Note, there are two minimum values of vorticity in the distal recirculation region, which is not present for the steady state calculation with the same spatial resolution. It is possible that this second minimum value is introduced by numerical errors. Analysis for this particular stenosis model was carried out using a finer mesh. The vorticity distribution at < = 0 is shown in Figure 8.31. It can be seen that although the grid has some effect on predicting the exact location of the second minimum value, the general trends remain the same. Hence, this second minimum value is not a numerical error and it corresponds to the second vortex. 0.75 t 0.5 0.25 Reverse Flow 1 R 0.5 1.5 X Figure 8.9 Time history of the separation and reattachment points for MO at Re Rw = 10, e = 0.1. R 0.75 h t 0.5 0.25 Reverse Flow R Reverse Flow Pi - 10 10 20 Figure 8.10 Time history of the separation and reattachment points for MO at iîe = 100, Rw = 10, e = 0.5. 0.75 h t 0.5 'r 0.25 h - 10 10 20 Figure 8.11: Time history of the separation and reattachment points for MO at Re = 100, Rw = 10, c = 1. 0.75 h t 0.5 h 0.25 h Figure 8.12 Time history of the separation and reattachment points for MO at Re Rw = 10, e = 1.5. 0.75 ^ t 0.5 r 0.25 ^ Reverse Flow R P. - 10 0 10 X 20 Figure 8.13 Time history of the separation and reattachment points for MO at Re Rw = 7.5, e = 1. 0.75 t 0.5 !-0.25 - 1 0 R^P3 S R Î Reverse Flow i R u p. 10 20 Figure 8.14 Time history of the separation and reattachment points for MO at Re - 100, Rw = 3.34, e = 1. 0.75 h t 0.5 r 0.25 |3 p Reverse ; Flow R / s teady -10 10 20 Figure 8.15 Time history of the separation and reattachment points for M5 at i?e = 100, Rw = 10, e = 1. 8.4 Discussion General effects of the steady Reynolds number and the smooth constriction parameters, H and L, were investigated numerically in the previous chapter. The pulsatile flow field differ substantially from that in the steady case, even at the same instantaneous Re. To obtain complete description of to describe the unsteady flow in a living animal will require an extensive experimental program and computer resources. The results presented here are limited to a simple pulsatile flow model. Nevertheless, some interesting features of the unsteady flow have emerged from the numerical simulations. Effects of the unsteady flow parameters {Rw and e) on the behavior of the attached and detached vortices as well as the distribution of the wall vorticity have been clearly identified. For the mild stenosis model MO, present data confirmed some of the results of O'Brien and Ehrlich 34] for Re = 100. For larger mean Reynolds number and a severe stenosis, model M5, the pulsatile flow field became rather complex. More than one vortex developed distal to the stenosis, both attached and detached. The fiow quantities, such as the vorticity, experienced large oscillations along the wall and several vortices developed, travelling away and towards the wall, and moving downstream during the cycle. From the data available for canine and human arteries [64], the fiow quantities in the larger arteries cannot be separated into a single harmonic wave plus a mean flow. In addition, the height of the instantaneous peak flux above the mean value can be several times the mean flux, i.e. the ratio of the peak instantaneous flux to the mean one e* can be much larger than 1. Although in the present study the case e = 1.5 is included, it is clear that the situation for attached and detached vortices in a natural stenosis is more complicated on both the proximal and distal sides of the stenosis. While it is not appropriate to correlate the present numerical results directly with the flow in arteries, principle features of the arterial flow can be inferred. They strongly suggest that fluid dynamic effects in atherosclerosis cannot be based on the mean flow» . It is interesting to consider the hypotheses of a high shear stress as a major factor in atherosclerosis [9]. Previously, based on the steady state calculation, it was concluded that a high shear stress is not likely to be responsible for initiation of the atherosclerosis. Ahmed and Giddens [35] also came to a similar conclusion based on their experimental measurements. As Woolf [1] pointed out. Fry [10] argued that the shear stress also needs to be oscillating to cause damage to the endothelial cells, which is consistent with the hypothesis of wall permeability. Other studies have also found that atherosclerosis usually occurs at locations where a unidirectional shear stress is not likely to exist [21,75]. Hence it is necessary to examine the distribution of wall shear stress and its time averaged value under pulsatile flow conditions. It was found that the maximum value of the time averaged wall vorticity also occurs in the constriction area, almost at the same location as in the steady state. The time averaged value is always larger than that in the steady flow value at the same mean Re. The value of wall shear stress varied considerably during the cycle, and is not unidirectional at all. Clearly the constricted area again should be the most preferable location for atherosclerotic lesions in circulatory system considering pulsatile effect of the blood flow. Since in the animal experiments lesions were rarely observed in this area, it can be concluded that the high shear stress is unlikely to be responsible for initiation of the atherosclerosis. Another popular hypothesis suggests that existence of a region of reverse flow is re-sponsible for atherosclerosis [8]. Although there is clinical evidences [21] that preferable locations for atherosclerosis are likely to coincide with separation points, this is mostly based on steady flow computation. \s discussed in the previous chapter, the recircula-tion region for steady flow appearing in the distal region is consistent with the animal experiments of Zarins et al. [31] for a mild stenosis (44% area reduction). Numerical results of the present study of the unsteady flow show that the duration and magnitude of the reversal are considerable for a moderate stenosis and low mean Reynolds num-ber (Re=100). Thus, a separation region appears to consistently coincide with a high occurrence of atherosclerosis lesions proximal to the stenosis in the experiments. For a mild stenosis (44% area reduction), the size of the separation region is smaller, and the duration is also shorter during a cycle. 80 60 40 20 0 -20 f - 0 . 1 f = 0 . 5 - ^ £ = 1 . 0 ^ f = 1 . 5 J' — • ^ \ w . • Steady State 0 0.2 0.4 0.6 0.8 1 t Figure 8.17 Effect of e on time history of the peak wall vorticity for MO at Re = 100, Rw = 10. 60 -10 ' ' ' ' 1 ! 0 0.2 0.4 0.6 0.8 1 t Figure 8.18 Effect of Rw on time history of the peak wall vorticity for MO at Re = 100, e = 1. 120 -20 Figure 8.19 Effect of Re on time history of the peak wall vorticity for MO at Rw = 10, e = 1. 300 250 200 100 50 --50 H = 0.5 H = 0.25 / \ / \ / \ / \ / \ - / \ / \ / \ \ \ / \ \ / ^ — - - ^ ^ ^ ^ \ Steady State / \ / / ! 1 1 ' 0.2 0.4 0.6 t 0.8 Figure 8.20 Effect of H on time history of the peak wall vorticity at Re = 100, Rw = 10, e = 1. 60 -20 ' ' '• 1 • ' ' < -4 -2 0 2 4 6 8 10 X Figure 8.22 Vorticity distribution along the waH for MO at Re = 100, Rw ^ 10, 80 -20 t=0.25 — 1=0.5 t=0.75 t=1.0 Figure 8.23 Vorticity distribution along the wall for MO at Re = 100, Rw = 10, -10 Figure 8.24 Vorticity distribution along the wall for MO at Re = 100, Rw = 3.34, 300 250 h 200 h 150 I--50 100 h Figure 8.25 Vorticity distribution along the wall for M5 at Re = 100, Rw = 10, 120 100 h -20 -4 - 2 0 2 4 6 8 10 Figure 8.26 Vorticity distribution along the wall for MO at Re = 600, Rw = 10, 80 -20 ' ^ ' ' ^ ' > ' -4 - 2 0 2 4 6 8 10 X Figure 8.27 Effect of e on vorticity distribution along the wall for MO at Re = 100, Rw = 10. 60 Steady State -10 ' 1 1 i 1 ' i -4 - 2 0 2 4 6 8 10 X Figure 8.28 Effect of Rw on vorticity distribution along the wah for MO at Re = 100, e = 1. 300 250 200 150 100 50 -50 H = l / 4 - I 1 ' 1 \ — H = l / 2 - w 1 ' --jy \ I ^ ^^^^ ^ — • — . / 1 1 \ / 1 1 1 1 -4 - 2 2 4 X 10 Figure 8.29 Effect of E on vorticity distribution along the wall at Re = 100, Rw = 10, e = 1. 120 100 h 80 60 h 40 h 20 h Re=100 Re=600 -20 - 4 . - 2 0 2 4 6 8 10 x Figure 8.30 Effect of Re on vorticity distribution along the wall for MO at Rw e = 1. 300 250 200 h 150 h 100 r 50 -50 -4 - 2 0 28 X 25 g r i d • 60 X 25 g r i d 100 X 50 g r i d 2 4 X 10 Figure 8.31 Effect of grid size on vorticity distribution along the wah for M5 at i?e = 100, Rw = 10, e= 1. C O N C L U S I O N Viscous flow in tubes with occlusions was investigated in the context of blood flow in the large arteries. Numerical methods were developed and simple experiments performed to study hemodynamic effects on atherosclerosis because of the difficulty in carrying out accurate measurements in vivo. Information on the velocity field and pressure distribu-tions around the occlusion in steady flow was accumulated from the experiments. Basic features obtained from the present numerical simulations were confirmed by these exper-imental measurements. Numerical results for steady and unsteady flows for a large range of flow parameters were compared with previous results. The lack of boundary values for pressure and vorticity provides one of the computa-tional difficulties in the numerical solution of the incompressible Navier-Stokes equations. In this study, a new computation scheme, called the interior constraint (IC) method, has been developed to resolve this difficulty for both stream function-vorticity and primitive variable formulations. Accurate second-order results were obtained for standard test flow problems. Most results agree with the previous data in the literature to two to three sig-nificant digits for the driven cavity flow. In spite of the simphcity, the IC method is quite efficient compared to SIMPLE, a commercially available numerical algorithm. Although the IC method demonstrated several interesting practical advantages, such as simphcity, accuracy, and efficiency, and it has been apphed successfully to flows through constricted tubes, it is not appropriate to claim it is theoretically superior to other methods. However, it is expected that for a certain class of problems, e.g. flows with sohd boundaries, the IC method will be an excellent practical method. As was shown in chapter 2, the necessary constraints can be provided in other ways, resulting in several existing numerical schemes classified as BC methods. It can also be shown that the need for the necessary constraints could be avoided using the idea of staggered grids, which forms the basis for another class of numerical schemes such as the M A C method 39] and the SIMPLE algorithm [41]. Both BC methods and methods using staggered grids have been used for decades. The IC method is proposed as an alternative which is simple and accurate for certain classes of problems. It is quite encouraging to find that both primitive variable and stream function-vorticity variable formulations agree with each other well, although the forms of these two sets of equations are different. As an apphcation of the IC method, blood fiow in stenosed arteries has been in-vestigated. It has been shown that agreement with results reported in the literature is achieved. The results of the present numerical simulations, for both steady and unsteady flow problems, show that there is a lack of evidence to support the high shear stress hypothesis of atherosclerosis. The existence of an extensive recirculation region is likely to be one of the hemodynamic factors responsible for the initiation of atherosclerosis lesions. flThe computational results revealed that low lateral pressure D,is unlikely to be the initiating mechanism. However, the exact mechanism is unknown. There remains, however, a considerable body of work to be carried out to complete this study. Firstly, there is stifl room for improvement in the efficiency of the IC method. For example, the nonhnear terms can be linearized by Newton's method [48], which has the advantage of faster convergence rates. However, the problem of solving the resulting large sparse matrix can be more costly than using Gauss-Seidel iteration combined with the Thomas algorithm for the tridiagonal systems. A more practical improvement, easier and not so costly, is to solve the two coupled linearized equations by a two-by-two block tridiagonal matrix technique, or to apply the direct Poisson solver instead of the ADI procedure. An alternative is to use the spectral approach instead of the finite difference method. Secondly, the assumptions of axisymmetric laminar flow of Newtonian fluid in a rigid passageway can be removed to better approximate the more realistic flow in the circula-tory system, (i) A distensible wall can be incorporated by assuming a visco-elastic arterial wall. The moving wall can be approximated by the same grid generation technique used in the present study. A grid system is then generated numerically at each time instants instead of at the beginning of the calculation, (ii) It would be interesting to take turbu-lent effects into account. A low Reynolds number fc — e model [90] can be used to calculate the disturbances which may result for very low Reynolds number, within physiological range. Numerical results combined with the experimental measurements [35] would be useful in physiological practice, (iii) The inclusion of three-dimensional effects would also be interesting, since the flow in living animal is basically three-dimensional, and different features such as secondary flows cannot be revealed by the two-dimensional simulation 37]. Initial investigations would be for the flow in straight tubes with an asymmetrical constriction, flow in curved tubes, and flow in tubes with branches. For the asymmetrical constriction, the experiments by Young and Tsai [23] showed that a difference in flow quantities exists between the asymmetrical and symmetrical constrictions. Calculations in curved tubes without constrictions [37] are consistent with the previous analytical re-sults [91,92]. The inclusion of both curvature and constrictions in tube flows wifl produce better approximation to what happens in the circulatory system. Numerical simulations of flow in branches have been carried out in two-dimensional channels [97,98]. It wifl be more realistic if a three-dimensional tube branch is investigated. Finally, more accurate laboratory measurements are needed to estabhsh complete picture on such a complicated biofluiddynamics problem. Limited measurements were carried out in the present study. With sophisticated measuring techniques, such as the L D A system, accurate results for mean velocity and velocity fluctuations can be obtained [35,36]. However, near-wall measurements remain a diflSculty, especiafly in the constricted area. The wall shear stress and separation and reattachment points cannot be determined accurately. This is caused partially by the different refractive indexes of the materials in the experiments. The focus location of the split laser beams is shifted when they pass through the tube wal l into the working fluid when curvature is present. By choosing a proper combination of the working fluid and material for the stenosis model, it may be possible to obtain the same refractive indexes. The flow velocity inside the constriction region could then be measured. Hopefully, these results can be used with numerical simulations and in vivo experiments to come to a more thorough understanding of atheroslerosis. Bibliography [1] Woolf, N. (1982) Pathology of Atherosclerosis. Butterworth Scientific, London. 2] Fallopins, G. (1575) Lections de partinus similaribus humani corporis, cited by E . Long In Cowdry's Arteriosclerosis (1967). Ed. by Blumenthal, H.T. , 5-20, Spring-field, Charles C. Thomas. [3] Long, E . (1967) In Cowdry's Arteriosclerosis (1967). Ed. by Blumenthal, H.B., 5-20, Springfield, Charles C. Thomas. [4] Mitchell, J .R.A. and Schwartz, C.J. (1965) Arterial Disease, Oxford, Blackwell Scientific Publication. [5] Schwartz, C .J . and Mitchell, J .R.A. (1962) Observations on location of arterial plaques. Circulation Research, X I , 63-73. [6] Texon, M . (1963) The role of vascular dynamics in the development of atheroscle-rosis. In Atherosclerosis and its Origin. Ed. by Scandel, M . and Bourne. G.H. , 167-195. New York, Academic Press. [7] Texon, M . (1974) Atherosclerosis: Its hemodynamic basic and implications. Medical Ginics of North America, 58, 257-268. [8] Caro, C.G. , Fitz-Gerald, J .M. , and Schroter, R.C. (1971) Atheroma and arterial wall shear: Observations, correlations and proposal of a shear dependent mass transfer mechanism for atherogenesis. Proceedings of the Royal Society of London (B), 177, 109-159. 9] Fry, D.L. (1968) Acute vascular endothelial changes associated with increased blood velocity gradients. Circulation Research, 22, 165-197. 10] Fry, D.L. (1973) Responses of the arterial wall to certain physical factors. In CIBA Foundation Symposium, No. 12 (New Series). Atherogenesis: Initiating Factors. Ed. by Porter, R. and Knight, J . , 93-120. Amsterdam, Excerpta Medica, North-Hollana, Elsevier. [11] Caro, C.G. , Fitz-Gerald, J .M. , and Schroter, R.C. (1969) Arterial wall shear and distribution of early atheroma in man. Nature, 223, 1159-1161. [12] Forstorm, F . J . quoted by Fry, D.L. (1978) Haemodynamic injury. In The Throm-botic Process in Atherogenesis. Ed. by Chandler, A .B . , Eurenins, K., Mcmillan, G.C. , Nelson, C.B., Schwartz, C.J. , and Wessler, S., 353-356, New York, Plenum Press. [13] Fox, J.A. and Hugh, A . E . (1966) Localization of atheroma: a theory based on boundary layer separation. British Heart Journal, 23, 388-399. [14] MuUer-Mohnssen, H., Kratzer, M. , and Baldauf, W.(1978) Microthrombus forma-tion in models of coronary arteries caused by stagnation point flow arising at the predilection site of atherosclerosis and thrombosis. In The Role of Fluid Mechanics in Atherogenesis. Ed. by Nerem, R.M. and Cornhill, J .F . , 12-1-12-8. Ohio State University. 15] Murphy, E .A. , Rowsell, H.E. , Downie, H.G. , Robinson, E .A. , and Mustard, J .F. (1962) Incrustation and atherosclerosis:the analogy between early in vivo lesions and the deposits which occur in the extra-corporeal circulations. Canadian Medical Association Journal, 87, 259-274. [16] Scharfstein, H. , Gutstein, W.H. , and Lewis, L.(1963) Changes of boundary layer flow in model systems: Implications for initiation of endothelial injury. Circulation Research, 13, 580-584. [17] Robertson, J.H.(1960) Influence of mechanical factors on structure of peripheral arteries and the localization of atherosclerosis. Journal of Clinical Pathology, 13, 199-204. [18] Wesolowski, S.A., Fries, C.C. , Sabini, A . M . , and Sawyer, P.N. (1965) The signifi-cance of turbulence in hemic systems and in the distribution of the atherosclerotic lesion. Surgery, 57, 155-162. 19] Nerem, R.M. and Cornhill, J .F . (1978) in Fluid mechanics and atherogensis. Ed. by Nerem, R.M. and Cornhfll, J .F . , 1-2, Ohio State University [20] WiUis, G.C. (1954) Localizing factors in atherosclerosis. Canadian Medical Associ-ation Journal, 10, 1-8. [21] Glagov, S., Rowley, D.A., and Kohout, R.I. (1961) Atherosclerosis of the human aorta and its coronary and renal arteries. A consideration of some haemodynamic factors which may be related to the marked differences in atherosclerotic involve-ment of the coronary and renal arteries. Archives of Pathology, 72, 571-578. [22] Duncan, L . E . , Jr., Buck, K., and Lynch, H. (1965) The effect of pressure and stretching in the passage of labelled albumin into canine aortic wall. Journal of Atherosclerosis Research, 5, 69-79. [23] Young, D.F. and Tsai, F Y . (1973) Flow characteristics in models of arterial stenosis I. Steady flow. J. Biomech., 6, 395-410. 24] Lee, J.S. and Fung, F . C . (1970) Flow in locally constriction tubes at low Reynolds numbers. J. Applied Mech., 37, 9-16. [25] Morgan, B . E . and Young, D.F. (1974) An integral method for the analysis of flow in arterial stenosis. Bulletin of Mathematical Biology, 36, 39-53. [26] Forrester, J .H. and Young, D.F. (1970) Flow Through a converging-diverging tube and its implication in occlusive vascular disease 1. theoretical development, 2. the-oretical and experimental results and their implications J. Biomech., 3, 297-316. [27] Fukushima, T . , Azuma, T . , and Matsuzama, T . (1982) Numerical analysis of blood flow in the vertebral artery. J. Biomech. Eng., 104, 143-147. [28] Deshpande, M.D. , Giddens, D.P., and Mabon, R.F. (1976) Steady laminar flow through modelled vascular stenoses. / . Biomech., 9, 165-174. [29] Giddens, D.P., Mabon, R.F. , and Cassanova, R.A. (1976) Measurement of dis-ordered flows distal to subtotal vascular stenosis in the thoracic aortas of dogs. Circulation Research, 39, 112-119. [30] Khahfa, A .M.A . and Giddens, D.P. (1978) Analysis of disorder in pulsatile flows with apphcation to poststenotic blood velocity measurement in dogs. J. Biomech., II , 129-141. 31] Zarins, O.K., Bomberger, R.A., Taylor, K . E . , and Glagov, S. (1980) Arterial steno-sis inhibits regression of diet induced atherosclerosis. Surgery, 88, 86-92. 32] Bomberger, R.A., Zarins, O.K., Taylor, K . E . , and Glagov, S. (1980) Effect of hyper-tension on atherosclerosis and aortic wall composition. J. Surg. Res., 28, 402-409. [33] Daly, B .J . (1976) A numerical study of pulsatile flow through stenosed canine femoral arteries. J. Biomech., 9, 465-475. [34] O'Brien, V . and Ehrhch, L.W. (1985) I. Simple pulsatile flow in an artery with a constriction. J. Biomech., 18, 117-127. [35] Ahmed, S.A. and Giddens, D.P. (1983) Velocity measurement in steady flow through axisymmetric stenoses at moderate Reynolds numbers. J. Biomech., 16, 505-516. 36] Ahmed, S.A. and Giddens, D.P. (1984) Pulsatile poststenotic flow studies with laser doppler anemometry. J. Biomech., 17, 695-705. 37] Pedley, T . J . (1980) The fluid mechanics of large blood vessels. Cambridge Univer-sity Press, Cambridge; New York. 38] Mueller, T . J . (1978) Application of numerical methods to physiological flows, in Numerical Methods in Fluid Dynamics. Ed. by Wirz, H.J. and Smolderen, J . J . , Von Karman Institute for Fluid Dynamics, Rhode-Saint-Genese, Belgium, Hemisphere publishing corporation, 89-154. 39] Harlow, F . H . and Welch, J .E . (1965) Numerical calculation of time-dependent vis-cous incompressible flow of fluid with free surface. Phys. Fluids, 8, 2182-2189. [40] Chorin, A . J . (1967) A numerical method for solving incompressible viscous flow problems. J. Comput. Phys. 2, 12-26. 41] Patankar, S.V. and Spalding, D.B. (1972) A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows. Int. J. Heat and Mass Transfer 15, 1787-1806. [42] Gresho, P.M. and Sani, R.L. (1987) On pressure boundary conditions for incom-pressible Navier-Stokes equations. Int. J. Num. Meth. Fluids, 7, 1111-1145. [43] Orszag, S.A., Isra«U, M . , and Deville, M.O. (1986) Boundary Conditions for In-compressible Flows. J. Scientific Comput. 1, 75-111. [44] Chorin, A . J . (1968) J. Comput. Phys. 22, 745-762. [45] Alfrink, B .J . (1981) in Proceedings, the Second International Conference on Nu-merical Methods in Laminar and Turbulent Flow, Venice, July 13-16, 1981. Ed. by C. Taylor and B.A. Schrefler (Pineridge Press). 46] Kleiser, L. and Schumann, U. (1979) Treatment of incompressibility and boundary conditions in 3-d numerical spectral simulation of plane channel flows, in Proced-ings, Third GAMM-Conference on Numerical Methods in Fluids Mechanics, Braun-schweig. Ed. by E . H . Hirschel, 165-173. [47] Roache, P.J. (1972) Computational Fluids Dynamics (Hermosa, Albuquerque). [48] Schreiber, R. and Keller, H.B. (1983) Driven cavity flows by efficient numerical techniques. J. Comput. Phys. 49, 310-333. [49] Woods, L . C . (1953) A note on the numerical solution of fourth order differential equations. Aeronauticak Quarterly, 5, part 3, 176. ;50] Ghia, U. , Ghia, K .N. , and Shin, C.T. (1982) J. Comput. Phys. 48, 387-411. [51] QuartapeUe, L . (1981) Vorticity conditioning in the computation of two-dimensional viscous flows. Comput. Phys., 40, 453-477. [52] Dennis, S.C.R. and QuartapeUe, L. (1989) Some uses of Green's theorem in solving the Navier-Stokes equations. Int. J. Numer. Methods Fluids, 9, 871-890. 53] Canuto, C. and Sacchi-Landriani, G. (1986) Analysis of the Kleiser-Schumann method. Numer. Math., 50, 217-243. 54] Thompson, J .F . (1978) Computational Fluid Dynamics, (Lecture Series, No. 1978-4, Von Karman Institute for Fluid Dynamics). 55] Aubert, X . and M . Deville, M . (1983) Steady viscous flows by compact differences in boundary-fitted coordinates. J. Comput. Phys. 49, 490-522. [56] AbdaUah, S. (1987) / . Comput. Phys., 70, 182-192. [57] Morgan, K., Periaux, J . , and Thomasset, F . Eds. (1984) Analysis of Laminar Flow over a Backward Facing Step, Vieweg, Braunschweig. [58] Morrison, J .H. and Napohtano, M . (1988) Efficient solutions of two-dimensional incompressible steady viscous flows. Comput. Fluids, 16, No.2, 119-132. 59] Armaly, B .F . , Durst, F. , Pereira, J .C .F . , and Schonung, B. (1983) Experimental and theoretical investigation of backward facing step flow. J. Fluid Mech., 127, 473-496. 60] Kim, J . and Moin, P. (1985) Apphcation of a fractional-step method to incompress-ible Navier-Stokes equations. J. Comput. Phys. 59, 308-323. 61] Peaceman, D.W. and Rachford, H.H., Jr. (1955) The numerical solution of parabolic and elliptic differential equations. SIAM J., 3, 28-41. [62] Ma, Y.W. (1983) Investigation of numerical methods for solving the Navier-Stokes equations. Mathematica Numerica Sinica, 5 (in Chinese). 63] Fletcher, C .A.J . (1988) Computational Techniques for Fluid Dynamics I (Springer-Verlag, Berlin, Heideberg). 64] Womersley, J.R. (1955) Method for calculation of velocity, rate of flow and viscous drag in arteries when pressure gradient is known. J. Physiol., 127, 553-563. [65] Kawai, H., Sawada, T. , and Tanahashi, T . (1987) Numerical simulation of a blood flow in a stenosed tube, in Role of Blood Flow in Atherosclerosis, Proceedings of International Symposium, Hyogo, Japan. [66] Karki, K . C . and Patankar, S.V. (1988) Calculation procedure for viscous incom-pressible flows in complex geometries. Numer. Heat Trans., 14, 295-307. [67] Patankar, S.V. (1978) A numerical method for conduction in composite materials, flow in irregular geometries and conjugate heat transfer. Proc. Int. 6th Heat Tran. Conf., Toronto, Canada. [68] Patankar, S.V. (1980) Numerical Heat Transfer and Fluid Flow, Hemisphere pub-lishing corporation. [69] Aminzadeh, M . (1975) Hydrodynamic performance of an artificial aortic valve im-plant, Ph.D. Thesis, The University of British Columbia, Vancouver, B.C. , Canada. [70] Akutsu, T . (1985) Hydrodynamic performance of mechanical posthetic heart valves, Ph.D. Thesis, The University of British Columbia, Vancouver, B .C. , Canada. [71] World Health Organization (1958) Classification of atherosclerotic lesions. WHO Technical Report Series, No. 143. Geneva; WHO. [72] Peyret, R. and Taylor, T .D . (1982) Computational methods for fluid flow, Springer-Verlag, New York. [73] Zarins, C.K., Bomberger, R.A., and Glagov, S. (1981) Local effects of stenoses: induced flow velocity inhibits atherosclerosis. Circulation, 64, 221-227. [74] Bomberger, R.A., Zarins, C.K., and Glagov, S. (1981) Subcritical arterial stenosis enhances distal atherosclerosis. / . Surg. Res., 3 0 , 205-212. 75] Roach, M.R. and Smith, N.B. (1983) Does high shear stress induced by blood flow lead to atherosclerosis? Persp. Biol. Med., 26, 2, 287-303. [76] Flaherty, T . J . , Pierce, J . E . , and Ferrans, V . J . (1972) Endothehal nuclear patterns in the canine arterial tree with particular reference to hemodynamic events. Circ. Res., 3 0 , 23-33. [77] Langille, B.L. and Adamson, S.L. (1981) Relationship between blood flow direction and endothelial cefl orientation at arterial branch sites in rabbits and mice. Circ. Res., 48, 481-488. [78] Legg, M.J . and Gow, B.S. (1982) Scanning electron microscopy of endothelium around an exp'îrimental stenosis in the rabbit aorta using a new casting material. Atherosclerosis, 42, 299-318. 79] Cornhill, J .F . , Levesque, M.J . , Herderick, E . E . , Nerem, R.M. , Kilman, J.W., and Vasko, J.S. (1980) Quantitative study of rabbit aorta endothelium using vascular casts. Atherosclerosis, 3 5 , 321-337. 80] Glagov, S., Zarins, C.K., Giddens, D.P., and Ku, D.N. (1988) Hemodynamics and atherosclerosis. Atherosclerosis, 112, 1018-1031. 81] Glagov, S., Rowley, D.A., and Kohout, R.I. (1961) Atherosclerosis of human aorta and its coronary and renal arteries. Arch. Pathol. Lab. Med., 71, 558-571. [82] Breterton, K .N . , Day, A . J . and Skinner, S.L. (1977) Hypertension-accelerated atherogenesis in cholesterol-fed rabbits. Atherosclerosis, 27, 79-87. 83] Chobanian, A . J . (1983) The influence of hypertension and other hemodynamic factors in atherosclerosis. Cardiovascular Dis., 26, 177-196. [84] Hollander, W., Prusty, S., Kirkpatrick, B. (1977) Role of hypertension in ischemic heart disease and cerebral vascular disease in the cynomolgus monkey with coarc-tation of the aorta. Circ. Res., 40 (suppl. I), 70-83. [85] Kannel, W.B. , Schwartz, M.J . , and McNamara, P.M. (1969) Blood pressure and risk of coronary heart disease: the Framingham study. Dis. Chest., 56, 43-52. [86] Robertson, W.B. and Strong, J.P. (1968) Atherosclerosis in persons with hyperten-sion and diabetes mehtus. Lab. Invest., 18, 538-555. [87] Sako, Y . (1962) Effects of turbulent blood flow and hypertension on experimental atherosclerosis. Am. Med. Ass., 179, 36-40. 88] Mestel, A . L . , Spain, D.M. , and Turner, H.A. (1964) Atheroma absence distal to subtotal aorta occlusion. Arch. Pathol. 78, 186. [89] Gessner, F .B. (1973) Hemodynamic theories of atherogenesis. Circ. Res., 33, 259-266. [90] Jones, W.P. and Launder, B .E . (1973) The calculation of low-Reynolds-number phenomena with a two-equation model of turbulence. Int. J. Heat Mass Transfer, 16, 1119-1130. [91] Dean, W.R. (1927) The stream-hne motion of fluid in a curved pipe. Phil. Mag., 4, 673-695. 92] Dean, W.R. (1928) Note on the motion of fluid in a curved pipe. Phil. Mag., 5, 208-223. 93] Nazemi, M. , Kleinstreuer, C , and Archie, J.P., Jr. (1990) Pulsatile two-dimensional flow and plaque formation in a carotid artery bifurcation. J. Biomech., 23, No.10, 1031-1037. [94] Douglas, J . , and Gunn, J . E . (1964) Numer. Math., 6, 428-453. [95] Peyret, R.C.R. (1971) Acad. Sci. A, 272, 1274-1277. [96] Daube, O., and Ta Phuoc Loc (1978) J. Mech., 17, 651-678. [97] Agonafer, D., Watkins, C.B., and Cannon, J .N. (J985) Computation of steady flow in a two-dimensional arterial model. J. Biomech. 18, 695-701. [98] Nazemi, M. , Kleinstreuer, C , and Archie, J.P., Jr. (1990) Pulsatile two-dimensional flow and plaque formation in a carotid artery bifurcation. J. Biomech. 23, 1031-1037. F O R M U L A T I O N O F S I M P L E M E T H O D The details of the discretization of the governing equations in flow problems by the SIMPLE method can be found in [68]. It follows the work of Patankar and Spalding 41] and was done on control volumes with a staggered grid arrangement. Here, only the expression for the coefficients and the source terms in the discretized momentum equations are given for reference. For the problems of interest, an advection-diffusion equation is considered: + ( U $ ) . + iv^)y - ^ [ ( | ^ ^ x ) x + ( j ^ $ „ ) „ ] = -P. + 3*, (A.l) With reference to Figure A . l , the discretized equation for an internal $-grid control volume may be cast as follows: (A.2) where ae = DeA{\Pe\) + max[ -F„ O; aw = Du,A{\P^\) + max[Fu„0 an = JD„A(|P„|) + max[-Fn, 0 as = D,A{\P,\) + max[Fa,0], ac = ae + aw + an + as ac°, 0 2Rw y Ax Ay ac = Re At 219 (A.3) (A.4) (A.5) (A.6) (A.7) (A.8) 5 * = s^yAxAy. (A.9) The function can be selected for the desired scheme for the convective term. For example, if the centered difference scheme is used, A(|P|) = 1 — 0.5|P|. P is the Peclet number (or the cell Reynolds number) defined a.s P = D/F. D are defined as: D, = ^y^^ ^ ^y^y, D^ = ^y^, £», = ^ ^ . (A.IO) * SxeRe^ SxyjRe^ " SynRe^ ' Sy^Re F is the flow rate defined as: = Ueî/Ay, Fy, = u^yAy, F„ = v^yAx, F, = v^yAx. ( A . l l ) Here $ is any physical quantity satisfying the advection-diffusion equation, and it can be u and V in the x and y momentum equations, respectively. The subscripts e, w, n, and s indicate the location at east, west, north, and south control volume faces. T Ay A X N < / X X / \ \ \ \ ' X / / • / f y ^ y / \ \ \ \ ' / >• X y \ N \ \ \ \ \ \ f y y y / n \ \ \ \ \ \ \ s \ \ \ s y Ë—t. îontrol Volume Figure A . l Control volume and notation for the SIMPLE method. A D I M E T H O D A N D L I N E G A U S S - S E I D E L I T E R A T I O N The Alternating Direction Implicit (ADI) [61] method is usually preferred in solving the Navier-Stokes equations. It has been used extensively both in its original form [61] and its modified version [94]. B . l Steady-State Solution We consider a two-dimensional advection-diffusipn equation: $t + - - ^ V ' * = 0, ( B . l ) where A and B are functions depending on the variable $. i?e is a constant (Reynolds number). The ADI method is a two-step method which is written as: - ^ij) + iA2 A ° - ; ^ A . . ) $ „ + (B, A ; -l-Ayy)^^;' = 0, (B.3) where Ai,A2,Bi, and B2 are approximations of A and B. The predictor value $,j can be considered as an approximation to the exact solution at time t = {n + | )Ai if A i , A2, and Bi,B2 are suitable approximations. For example, if only first-order accuracy in t is required, the simple expressions Ai = A2 = A'>j and B\ = B2 = 5,", where A^^- = A($"j), and B^- = 5($^j), can be used. To obtain second-order accuracy in the suitable approximations are: A l = CoA^/^ + (1 - Co - Ci)A.';. + CrA^-\ (B.4) 222 A l = (1 - Co + Cx + C2)Al+' + (Co - C i - 2C2)A'^j + C2A^-\ (B.5) Bi = DoB^^' + (1 _ Do - Di)Brj + D^B^-', (B.6) and B i = (1 - £>o + + Z)2)A""'' + ( ^ 0 - A - 2i52)B]5. + P a ^ r ' ' where C o , C i , C 2 and DQ,DI,D2 are arbitrary parameters. Usually, the two steps (B.4) and (B.5) are coupled, and an iterative procedure is required. To satisfy the stabil ity requirement (following Peyret and Taylor [72]), there is some restriction on the cell Reynolds number Rte — max[\Ai\AxRe, jBi jAyiîe] , and the time step size A i . REC < 2 must be satisfied to obtain diagonal dominance for the tridiagonal algebraic systems in (B.4) and (B.5). If this condition does not hold, the l imitat ion on the time step: 2Ax2 2Aj/2 \Ai\Ax-^^' \B2\Ay-^^ wi l l ensure diagonal dominance. The l imitat ion on the time step introduced by the condition of diagonal dominance, although not very restrictive, can be avoided if an upwind non-centered approximation of the first-order derivatives (according to the sign of A i and ^2) are used instead of centered approximations. In this case, second-order accuracy is lost. To recover the accuracy, Peyret [95] proposed a technique consisting of alternating the direction of the non-centered difference at each step for the primitive variable form. It has been used for the stream function-vorticity equation by Daube and T a Phuoc [96]. The single step method proposed by M a [62] was found to be more stable and less oscillatory. A l l of them are first-order in time and space during the transient stage. Second-order accuracy is obtained at the steady state if the time step (or the parameter a in Ma's scheme) is the same order as Ax and A y . B.2 Unsteady Solution For the solution at the transient stage, second-order accuracy cannot be obtained by the ADI scheme even if a centered difference or high-order upwind scheme is used. This is due to the nonhnearity of the convective terms and the implicit nature of the boundary condition at the sohd wall. To obtain second-order accuracy, an iterative method is used in the computation. The advection-diffusion equation (B.l) is discretized as: ^ ( ^ . i - n) + (A"-^' A°-l-Ax^)^p + (B"+^ A ° = O, (B.9) where A "+2 = (A" -h A"+i)/2, = (B " -h B"+i)/2, and = {^^j + ^lf^)/2. A Crank-Nicolson scheme is used for both convective and diffusive terms. The resulting algebraic system is pentadiagonal, and it can be solved iteratively by Gauss-Seidel line relaxation at each time level. W A L L S H E A R S T R E S S A N D V O R T I C I T Y The vorticity ( is defined as Vx — Uy in the cyhndrical coordinate system, where u and v are velocity components along x (axial) and y (radial) directions, respectively. The shear stress is defined as r^y = Ur + If the surface has an angle 6 with the z-direction, the shear stress is T'^^ = {vy — Ux)sin{29) + (Ux + Uy)cos{29). The shear stress r'^^y and the vorticity C are identical at the sohd boundary, which is the inner wall of the arteries. It wiU be shown as follows. In the general curvilinear coordinate system (^ , i/), we have 7/ = 1, and u = u = 0 on the sohd wall. Hence: u = \i-^r,Vi + X^Vr,) = ^X^Vr,, since = = 0. We also have: sin(2e) = 1 + ^ 2 ' 2s cos{29) = 2 ' 3 = tanO = ^. x^ Then the shear stress at the sohd wall with a slope s = tand becomes: = [vy - Ux)sin{26) + {vx + Uy)cos{29) 1/ ^ 25 1, Lui ! 1 + 52 Appendix C. WALL SHEAR STRESS AND VORTICITY The vorticity at sohd wall is: Thus, T^y = C at the sohd wall.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Incompressible viscous flow in tubes with occlusions
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Incompressible viscous flow in tubes with occlusions Huang, Huaxiong 1991
pdf
Page Metadata
Item Metadata
Title | Incompressible viscous flow in tubes with occlusions |
Creator |
Huang, Huaxiong |
Date Issued | 1991 |
Description | Viscous, incompressible flow in tubes with partial occlusion is investigated using numerical and experimental procedures. The study is related to the problem of atherosclerosis, one of the most common diseases of the circulatory system. One of the computational difficulties in solving the incompressible Navier-Stokes equations is the lack of pressure or vorticity boundary conditions. A finite difference approach, referred to as the interior constraint (IC) method, is proposed to resolve this difficulty. As a general numerical method, it is formulated for both the stream function-vorticity and primitive (physical) variable formulations. The procedure is explained using onedimensional model with extensive numerical tests presented for two-dimensional cases including flow in a driven cavity, and flow over a backward facing step. Results are obtained with second-order accuracy. Next, the IC method is applied to flow in a tube with an occlusion, which is used as the model for blood flow in stenosed arteries in the study of the pathology of atherosclerosis. Numerical results are obtained for both steady and pulsatile flows. Results are compared with those of S I M P L E , one of the commercially available numerical algorithms. The pulsatile flow study revealed several interesting new features. It suggests that the high shear stress is not likely to initiate atherosclerosis lesions. The recirculation region, which is a prominent feature of the unsteady flow, is more likely to cause the initiation and development of the disease. Experimental measurements for steady flow complement the numerical study and show qualitative agreement. |
Extent | 8246292 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2008-12-18 |
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.0079819 |
URI | http://hdl.handle.net/2429/3191 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1992-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_1992_spring_huang_huaxiong.pdf [ 7.86MB ]
- Metadata
- JSON: 831-1.0079819.json
- JSON-LD: 831-1.0079819-ld.json
- RDF/XML (Pretty): 831-1.0079819-rdf.xml
- RDF/JSON: 831-1.0079819-rdf.json
- Turtle: 831-1.0079819-turtle.txt
- N-Triples: 831-1.0079819-rdf-ntriples.txt
- Original Record: 831-1.0079819-source.json
- Full Text
- 831-1.0079819-fulltext.txt
- Citation
- 831-1.0079819.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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0079819/manifest