Pulsatile Blood Flow Simulations in Computed Tomography(CT) Scan-Based and Idealized Geometries of Human Aorta

(1)

Lappeenranta University of Technology Faculty of Technology

Department of Mathematics and Physics

Pulsatile Blood Flow Simulations in Computed Tomography(CT)

Scan-Based and Idealized Geometries of Human Aorta

The topic of this Master's Thesis was approved by the faculty council of the Faculty of Technology on May 26 2010.

The examiners of the thesis were Professor Heikki Haario and Docent Payman Jalali.

The thesis was supervised by Docent Payman Jalali.

Lappeenranta, August 23, 2010

Zuned Mansuri

Korpisuonkatu 17 C 32, Lappeenranta 53850, Finland.

Phone: +358466303055 zuned.mansuri@lut.

(2)

Abstract

Lappeenranta University of Technology Department of Mathematics and Physics

Zuned Mansuri

Pulsatile Blood Flow Simulations in Computed Tomography(CT) Scan-Based and Idealized Geometries of Human Aorta

Master's Thesis 2010

43 pages, 33gures,5 tables, 0appendices

Examiners: Professor Heikki Haario Docent Payman Jalali

Keywords: Medical imaging, Arterial disease, Hemodynamics, Computational Fluid Dynamics, Wall Shear Stress

Blood ow in human aorta is an unsteady and complex phenomenon. The complex patterns are related to the geometrical features like curvature, bends, and branching and pulsatile nature of ow from left ventricle of heart. The aim of this work was to understand the eect of aorta geometry on the ow dynamics. To achieve this, 3D realistic and idealized models of descending aorta were reconstructed from Computed Tomography (CT) images of a female patient. The geometries were reconstructed using medical image processing code. The blood ow in aorta was assumed to be laminar and incompressible and the blood was assumed to be Newtonian uid. A time dependent pulsatile and parabolic boundary condition was deployed at inlet. Steady and unsteady blood ow simulations were performed in real and idealized geometries of descending aorta using a Finite Volume Method (FVM) code. Analysis of Wall Shear Stress (WSS) distribution, pressure distribution, and axial velocity proles were carried out in both geometries at steady and unsteady state conditions. The results obtained in thesis work reveal that the idealization of geometry underestimates the values of WSS especially near the region with sudden change of diameter. However, the resultant pressure and velocity in idealized geometry are close to those in real geometry.

(3)

Acknowledgements

I express my warm and sincere thanks to the supervisor of the thesis, Dr. Payman Jalali for his close supervision and valuable guidance. I am grateful to him for giving me an opportunity to work in research project related to numerical simulation of blood ow in large arteries. I thank Prof. Heikki Haario for examining the thesis. Their guidance, valuable suggestions and encouragement helped me to successfully complete this work.

I am very thankful to Dr. Tuomo Kauranne for his continuous encouragement, advises, and comments on the work. I also thank him for arranging computational resources at CSC-IT Center for Science, Finland. I especially acknowledge my gratitude to Paritosh Vasava for his contributions to this work. He gave innovative ideas, inspirations, comments, and advises both in this work and in his friendly company during my entire stay in Finland.

The thesis work was funded by the Academy of Finland under the grant No.123938. The entire thesis work was carried out at Department of Energy Technology in Lappeenratna University of Technology. I am thankful to everyone at Department of Energy Tech- nology for their support and encouragement.I would like to thank Chief Radiologist Matti Sauna-Aho, South Karelia Central Hospital at Lappeenranta for providing CT scans used for this work. I would also like to thank Dr. Pertti Kolari for administrative assistance for acquisition of CT images. I am thankful to Department of Mechanical Engineering for providing medical imaging software 3D- Doctor used for processing CT images. I would also like to thank Department of Mathematics for providing nancial support during course of my stay in Lappeenranta.

My wishes and love to my best friends Ashvin, Sandeep and Arjun who shared my happiness and sorrows. Together we have faced every challenges during the studies and my entire stay in Finland. These experiences have helped me grow stronger and condent. I am highly indebted to my beloved friend Hemal for her constant inspiration, trust and belief in my work.

I am forever grateful to my loving parents and my entire family for their sacrices and eorts for upbringing me as a person. I deeply express my sincere love, feelings and respect to them. I would not have been here without their encouragement and support.

Above all, I thank Allah for the many showers of blessings upon my life and may it continue in future.

Lappeenranta, August 23, 2010 Zuned Mansuri

(4)

VOCABULARY

WSS Wall Shear Stress

CFD Computational Fluid Dynamics CT Computed Tomography

CAT Computed Axial Tomography 2D Two Dimensional

3D Three Dimensional

MRI Magnetic Resonance Imaging ROI Region of Interest

DICOM Digital Imaging and Communications in Medicine CAD Computer Aided Design

PDE Partial Dierential Equation FDM Finite Dierence Method FVM Finite Volume Method FEM Finite Element Method UDF User dened Function

(7)

NOTATIONS

General

T Threshold

x, y Coordinate values (X_avg, Y_avg) Centroid

h Cross-sectional distance between any two CT slices d Curve length between any two CT slices

V Velocity of blood at entrance plane

D Diameter of entrance plane in descending aorta geometry

R_e Reynolds number

u x-component of velocity v y-component of velocity z z-component of velocity

t Time

Du

Dt Total or substantive derivative ofu with respect to time

p Static pressure

∇ ·p Pressure gradient

f Body force acting on the uid u_x Velocity in x-direction

u_max Maximum velocity

m Number of rows in matrix representation of an image n Number of columns in matrix representation of an image

(8)

Greek Symbols

ρ Density of blood (kg/m³) µ₀ Viscosity of blood(P a)

δ_i Small element of uid in i^th direction τ_i Stress component in i^th direction λ Volumetric deformation

µ Linear deformation

(9)

List of Figures

1.1 Schematic demonstration of human aorta . . . 1

2.1 Gray-level histogram of a threshold . . . 6

2.2 (a) Cross sectional CT scan image (b) Cross sectional CT scan image after contrast adjustment . . . 8

2.3 CT scan image after segmentation process . . . 9

2.4 3D surface model of descending aorta . . . 9

2.5 7^th degree polynomial t to radius for rst 215 CT slices . . . 11

2.6 Distances between two consecutive CT slices . . . 12

2.7 2D axisymmetric geometry of descending aorta . . . 12

3.1 a boundary layer mesh in 2D axisymmetric model . . . 14

3.2 Wall shear stress in a 2D axisymmetric geometry from dierent mesh sizes 15 3.3 Mesh in cross sectional inow plane of a 3D descending aorta geometry . 16 3.4 Boundary conditions for (a) 3D (b) 2D axisymmetric geometries of descending aorta . . . 20

3.5 Parabolic velocity proles at inlet of (a) 3D (b) 2D axisymmetric geometries of descending aorta . . . 21

3.6 Velocity pulse function . . . 21

4.1 Dierent time steps in a cardiac cycle of inlet velocity prole . . . 23

4.2 WSS distribution in a 2D axisymmetric geometry in steady state simulation . . . 24

4.3 WSS(P a) distribution in a 3D geometry in steady state simulation . . . 24

4.4 Average WSS distribution in a 2D axisymmetric geometry at dierent time steps . . . 26

4.5 WSS(P a) distribution in a 3D geometry at dierent time steps . . . 26

4.6 Distribution of pressure(mmHg)in (a) 2D axisymmetric (b) 3D geometries of the descending aorta in steady state simulation . . . 28

4.7 Pressure (mmHg) distribution in a 2D axisymmetric geometry at dierent times steps . . . 28

(10)

4.8 Pressure (mmHg) distribution on the wall of a 3D geometry at dierent time steps . . . 30 4.9 Velocity(m/s)distribution in a 2D axisymmetric geometry in steady state

simulation . . . 31 4.10 (a) axial cross-sections in a 3D geometry (b) velocity (m/s) proles at

dierent cross-sections of a 3D geometry normal to the ow direction in steady state simulation. . . 32 4.11 Velocity (m/s) distribution in a 2D axisymmetric geometry at dierent

times steps . . . 32 4.12 (a) Axial cross-sections in the 3D geometry (b) Velocity (m/s) proles at

dierent cross sections at dierent time steps . . . 34 4.13 Comparison of 2D and 3D geometries for average WSS distribution in

steady state simulation . . . 35 4.14 Comparison of 2D and 3D geometries for average WSS distribution at

time t=0.15,0.25,0.35,0.45,0.55,0.75 s. . . 35 4.15 Comparison of 2D and 3D geometries for average pressure distribution in

steady state simulation . . . 36 4.16 Comparison 2D and 3D geometries for average pressure distribution at

time t=0.15,0.25,0.35,0.45,0.55,0.75 s. . . 37 4.17 Comparison of 2D and 3D geometries for pressure drop at dierent time

steps t=0.15,0.25,0.35,0.45,0.55,0.75 s. . . 38 4.18 Comparison of 2D and 3D geometries for average velocity distribution in

steady state simulation . . . 38 4.19 Comparison of 2D and 3D geometries for average velocity distribution at

time t=0.15,0.25,0.35,0.45,0.55,0.75 s. . . 39

(11)

List of Tables

4.1 Average WSS in a 2D axisymmetric geometry at dierent time steps of a cardiac cycle . . . 25 4.2 Average WSS in a 3D geometry at dierent time steps . . . 27 4.3 Pressure values at the inlet and the outlet in a 2D axisymmetric geometry

at dierent time steps . . . 29 4.4 Pressure values at the inlet and the outlet in a 3D geometry at dierent

time steps . . . 29 4.5 Maximum velocity at dierent times steps in a 2D axisymmetric geometry 33

(12)

1 Introduction

Cardiovascular diseases are a class of diseases that involve the heart or blood vessels (arteries and veins). They are one of major causes of human death in the modern world. This happens due to our relatively limited understanding towards the cause of the diseases. The diagnosis and treatment techniques may require specic techniques which dier from person to person. The cardiovascular system includes the heart and the blood vessels. The arteries are elastic blood vessels that carry oxygenated blood and vital nutrients from the heart and distribute them to dierent organs of the body.

The aorta is the main systemic artery and the largest artery of the human body. It originates from the left ventricle of the heart and extends down to the abdomen, where it bifurcates into two small branches called common iliac arteries as shown in Figure 1.1. The oxygenated fresh blood from the heart via left ventricle enters the ascending aorta. The common carotid arteries supply blood to the upper part of the body, whereas descending aorta delivers blood to the lower abdomen part of the body. Comparatively, the aorta has the largest geometric tapering as it moves away from the ventricle towards abdomen [1].The walls of the arteries consist of three layers: The inner layer or tunica intima, the middle layer or tunica media and the outer most layer or tunica externa.

The arteries possess good elastic property which supports expansion and contraction of the wall to drive the blood away. The innermost layer is smoother than the outer layer, which helps the blood to ow smoothly.

Figure 1.1: Schematic demonstration of human aorta

(13)

In normal conditions the blood ow does not face any obstacles on the way and the ow is regular and smooth. But sometimes due to the presence of cholesterol, calcium and other substances in blood, a fatty substance called plaque develops near the inner wall of the artery. As the time passes this plaque grows resulting in hardening the arterial wall and narrowing the blood passage. This can cause severe diseases such as the development of atherosclerosis. As a result the arterial wall loses its elastic property which limits the area of blood ow. Therefore the ow turns to be abnormal in the reduced cross sectional area of the artery.

The presence of unusual hemodynamic conditions in the arteries often creates abnormal biological responses [2]. Earlier observations show that most of the aneurysms (a balloon like bulge) are due to atherosclerosis which develops near the aortic arch and the abdominal aorta [3]. The normal arterial blood ow is laminar with secondary ows generated at curves and branches. It is known that early lesions of atherosclerosis occur in the regions of complex vascular geometry, such as bends, tapering and branching [4].These are regions where major ow disturbances occur such as secondary or reversed ow, and ow separation. The skewing of velocity in such regions could cause oscillating direction of wall shear stress (WSS) which can create pockets leading to atherosclerotic diseases [2]. Thus the analysis and study of uid dynamic quantities such as WSS may play an important role in association with other factors responsible for the progress of arterial diseases. This was the main reason to conduct this study.

On the other hand, it has been revealed that Computational Fluid Dynamics (CFD) may be used in analysis and prediction of uid ow, heat transfer, mass transfer, chemical reactions, and biological processes by solving mathematical equations which governs these processes. CFD is widely used in the industrial and academic applications. It is one of the branches in uid mechanics that uses numerical methods and algorithms to analyze and solve problems. Computer based simulations are run using the capacity of several processors at a time. Since a reliable ow simulation requires an accurate denition of vascular geometry, it is necessary to construct an exact model of the real geometry. An approximate model will give some results, but they may dier from the real ones. However, it is important to realize the dierences that the geometry simplication can apply on the results of simulations.

A real 3D model of vessel geometry can be reconstructed from Computed Tomography (CT) or Computed Axial Tomography (CAT) scan images of the body of a patient.

These are medical imaging methods to obtain detailed pictures of the structures inside the body. Modern CT scanners are very fast and they can generate a huge stack of images covering large sections of human body within few seconds. These images are two-dimensional gray scale images and are used for diagnosis and treatment pur-

(14)

poses. Nowadays there are many commercial tools available which can process medical images and reconstruct a real model of the vessel. These programs use image processing techniques like image contrast, image segmentation, and threshold selection for reconstructing 3D models which is then simulated. Previous numerical studies have shown dependence of ow on vessel geometry. Since realistic geometries have a complex structure, a simplied version of vessel geometry can be constructed using information available from the real data but assuming it to be planar with circular cross sections.

This would also reduce the complexity and save our eort. Observations from such geometries may result in better understanding the inuence of curvature, tapering, and irregularities on the ow proles.

In this thesis work, we have focused mainly on the descending aorta. A CFD method was used to investigate the ow problem. Blood ow simulations were performed in both a real 3D human descending aorta reconstructed from CT scan images and in an idealized model of it. Simulations were done in both steady and unsteady conditions.

The working domain and the boundary conditions of the geometries were dened in pre- processor software ANSYS Gambit. The nite volume analysis was performed using the nite volume method(FVM) package of ANSYS Fluent V6. The results of cross sectional velocities, pressure, and shear stresses from both geometries were then compared.

(15)

2 Medical image processing and reconstruction of or- gans

A medical image is a kind of two-dimensional (2D) monochrome image representing parts or organs of the body. Medical images are used for clinical purposes especially in diagnosis and in examining diseases. Defects and diseases in physiological processes can cause changes in the normal processes which may aect body performance and the health. It is necessary to have a better understanding of the body system so that we can observe and study patterns of behavior of the system. It can be relatively simple to measure the events that occur at the outer surface of the body, but most of the important organs are located within the body. In such situations the use of radiation becomes crucial in tracking and investigating the problems. There are several techniques for medical imaging like CT scans, Magnetic Resonance Imaging (MRI), Ultrasound, Nuclear scans and Angiographic imaging.

CT is one of the most commonly used imaging methods in medical diagnosis. In this technique, a bunch of x-rays is delivered in multiple directions to create 2D cross sectional slices of the human body. Later, three-dimensional (3D) objects can be reconstructed from the set of cross sectional planar images with the help of computers. While conduct- ing diagnosis on this 2D cross sectional images, several factors may aect the quality and result in the loss of correct information content of the images. Many organs of our interest are located within the body inside protective case. It is hard to acquire information about them. The image obtained at one moment can dier from that taken at other moment. For example the rhythmic contractile activity of the heart and the pulsatile blood ow through the arteries create challenges in angiographic imaging. Apart from all these reasons the safety of the patient from radiation is of at most importance.

This leads us to minimize the use of radiation and encourages us to focus more on the image processing steps which can improve the quality of images for accurate diagnosis purposes [5].

Our aim in this work is to reconstruct a descending aorta model from CT scan images to be used in blood ow simulations. For this purpose, data of a female patient was obtained from the South Karelia Central Hospital, Lappeenranta. The images obtained were 2D gray scaled images. In the following section we will describe concepts and techniques for analyzing medical images.

(16)

2.1 Image processing techniques

Mathematically, a monochrome image is a representation of a matrix of real numbers with m rows and n columns. Each element of the matrix is called an element or a pixel. Each pixel has its own address and its address corresponds to its coordinates. A pixel is a sample of the original image which means more samples give a more precise representation of the image. Accurate and clear vision requires some kind of operation to be done on images. Thus we have performed image processing techniques like image contrast, image segmentation, and thresholding over these CT scan images for better visualization.

2.1.1 Image contrast

CT scanned images may be made darker or brighter to properly identify an organ of the body. Images with low contrast may result from poor illumination or due to a wrong setting of lens during image acquisition. In such situations the image is passed through various processing operations such as image sharpening and contrast adjustment. The aim behind contrast adjustment is to increase the range of gray levels to dierentiate objects from their neighbors. Doing so will make the region of our interest clearly visible.

Contrast is dened by the dierence of brightness between the region of interest and the suitably dened background. An image is lacking in contrast if there are no sharp dierences between the black and white colors.

The pixel values that occur below a specied value are displayed black and the pixel values above a specied value are displayed as white, whereas the pixel values between these values are visible with dierent shades of gray. Since the image in our case has 16 bits per pixel, the range of the grayscale values is 0 to 65,536.

2.1.2 Image segmentation

Detection, localization, diagnosis, and monitoring treatment are the most important aspects of the diagnostic medical studies [7]. So, early detection and accurate localization of the diseases can be very useful in improving treatment procedures [5]. This process consists of two steps; the rst one is selecting the region of interest (ROI) and the second one is segmenting the selected region. Accurate quantication of the region of interest depends upon accurate segmentation of ROI in the image. The logic behind selecting the ROI is to restrict the area of data set and concentrate more on the extracted region of our interest. This in fact, will save computation and processing time. Practically,

(17)

selection of a ROI is time consuming because all the slices have to be processed one by one.

Image segmentation plays a vital role in understanding, analyzing and interpreting medical images. It is the process of dividing the image into various parts or objects. It is the initial and essential step before the observation and conclusion steps based on the image. Segmentation is not only important for feature extraction and visualization but also for image measurement and compression. The segmentation of data obtained with anatomical or structural imaging (an imaging technique that examines structural basis of changes caused by disease) is much easier than with functional imaging (an imaging technique used for detecting or measuring changes in metabolism, blood ow, regional chemical composition, and absorption).

2.1.3 Thresholding

In the image segmentation process, thresholding or threshold selection plays a key role due to its natural properties and simplicity of implementation. In this technique a threshold is selected and an image is divided into groups of pixels having values less than the threshold and groups of pixels with values greater than or equal to threshold [8]. Several thresholding methods have been developed based on either image histograms or on local properties such as mean values and standard deviations or local gradients.

The most common approach is global thresholding. When only one threshold is selected based on the image histogram for the entire image then the thresholding is called global.

Suppose that the gray level histogram of an image f(x, y) consists of light objects on dark background such that their pixels have gray levels grouped into two dominant modes. One way to extract these objects from the background is to select a threshold T that separates these modes.

Figure 2.1: Gray-level histogram of a threshold

(18)

Then any pointf(x, y)> T is called object point (region of interest) otherwise the point is called a background point. The threshold image g(x, y) is dened as:

g(x, y) = 1, iff(x, y)> T

= 0, iff(x, y)≤T (2.1)

The result of the thresholding is a binary image, where pixels with value 1 correspond to objects and pixels with value 0 correspond to background. Global thresholding is a simple and fast method. It works well when the image contain objects with a high contrast between the object and the background. But it may not lead to an accurate result when the objects have overlapping intensity levels. Threshold selection becomes dicult when the noise level in an image is high or the image contrast is low.

2.2 Reconstruction of descending aorta geometry

Since CT scans are 2D images, they are rst processed with the help of the image processing software package 3D-DOCTOR. This application is used for image visualization, surface rendering, volume rendering, and measurement of medical images. It extracts object boundaries using 3D image segmentation functions and creates surface rendering for visualization, object measurement, and quantitative analysis [9].

Figure 2.2 (a) presents a cross sectional CT scan image of the human body taken when lying horizontally. The data was obtained of a 55 year old female patient from South Karelia Central Hospital, Lappeenranta. The data consists of 393 slices (CT scan images) taken from the upper chest to the abdominal region of the body. The le format of images was DICOM-Digital Imaging and Communications in Medicine (DCM) and the type of the image was grayscale with 0.5mmslice thickness (the sum of the image slice thickness plus the gap distance between two neighboring slices).These images consist of 16 bits per pixel (65,536 gray levels). In Figure 2.2 (a) it is dicult to distinguish between bones, tissues and other small organs of the body. Therefore, it is necessary to make visual modications by performing some operations on image which would allow us to dierentiate between the various organs present inside the body. After then the right object of our interest is selected to generate a 3D model. Image contrast was set between32787and33127intensity levels which made the aorta clearly visible just above the spinal cord as shown in Figure 2.2 (b).

(19)

Figure 2.2: (a) Cross sectional CT scan image (b) Cross sectional CT scan image after contrast adjustment

The next step is to select the region of our interest (ROI). This will not only save computational and processing time, but also eort, as we might not be interested in processing the entire image unnecessarily. This process includes dening a closed region or boundary around the cross section of the aorta. After then segmentation operations are performed which can be done both locally and globally. A global ROI is the same for all the CT images which is time ecient, but it might increase the risk of errors in nal results. A local ROI is performed individually on all images which are time consuming in case of a large stack of images, but it gives minimum error in nal results. The error here means objects that are wrongly identied.

Now, the upper and lower thresholds are to be selected. In this stage, the object of our interest is separated from the rest of the objects present in the ROI. Threshold selection will highlight the objects that have values in the dened interval range. The selection of threshold must be done precisely as the nal object depends upon the threshold limits. This procedure plays a vital role in obtaining local microscopic details in the nal surface. The maximum and minimum limits were set to 34416 and 33120, and all the images were segmented. As a result, all the boundaries present in the ROI were approximated by a line and were highlighted in blue color as shown in the Figure 2.3.

There may be several objects that have been wrongly identied and are not of interest.

Such objects can be removed by performing some simple operations.

(20)

Figure 2.3: CT scan image after segmentation process

Finally, all the CT images were passed through the segmentation process and the object of our interest was approximated by a closed boundary. Then a three dimensional surface was generated over these contours as shown in Figure 2.4. The rendered surface is a wireframe made up of small surface polygons. The number of polygons can be controlled by a density parameter. A smaller value of this parameter would generate a wireframe with higher number of surface polygons and takes a longer time to process, whereas a larger value generates a lower number of surface polygons and takes a relatively shorter time to process.

Figure 2.4: 3D surface model of descending aorta

The outer surface of the model aorta geometry is rough but it can be smoothed if required. The generated 3D model is a surface but a volume domain cannot be created inside due to limitation of the software package 3D DOCTOR. For a blood ow simula-

(21)

tion the domain is required to be a volume such that it is suitable for CFD pre-processors to dene domain properties. So the model was processed in Computer Aided Design (CAD) program for volume conversion.

2.2.1 Surface to volume conversion

The IGES le format is one of the most commonly used le formats for interchanging 3D CAD models. It consists of points, lines and areas .The surface generated by the 3D DOCTOR program has a triangular surface mesh. It is not advisable to mesh on an already discretized geometry, because this would destroy the model quality further. One way is to export a set of contours which outlines the geometry to CAD programs. But this path fails, as most CAD programs require a guiding curve to reconstruct a surface and this has to be done manually by selecting each prole. It would be inecient if the geometry is not regular or has a large number of contours. Another possible way is to export a set of points and then create surfaces instead of triangles in the CAD environment. Here a set of points describing surface geometry from 3D DOCTOR was imported to the design program Catia v6 which can construct a solid shell over the imported points. Then the surface le was saved in IGES format which is suitable to create a volume domain by the CFD Pre-processor ANSYS Gambit.

2.3 Transformation to an idealized geometry of the descending aorta

Many studies have been carried out on curved vessel geometries to establish the eects of vessel curvature. Blood ow in the aorta is one of such complex ow situations in the cardiovascular system. This is mainly due to curvature, tapering, branching, and irregular geometry. Previous numerical studies have shown that velocity proles and wall shear stresses depend on varying anatomical structures of vessel geometry [10].

From a uid mechanics point of view, it is known that major ow disturbances such as formation of secondary ow, ow separation etc. occur in regions of complex vascular geometry [3]. Since a realistic descending aorta has complex structure, in this work an attempt was made to consider what eect factors such as geometrical features like curvature and tapering might have on blood ow in an idealized or simplied descending aorta. The ideal geometry was assumed to be planar with circular cross section and a monotonically decreasing diameter from inlet to outlet.

To maintain preciseness, the information available from the real data was used to construct an ideal model of the descending aorta. We know from our earlier discussion

(22)

that the region of interest (the cross section of the aorta) in all the CT scan images is approximated by a closed boundary. This entire set of images with DCM le format were converted to the bitmap image format with797×794pixels and image information was obtained with the help of the mathematical programming tool Matlab. From the image aspect of view, a line is a collection of pixels and each pixel has its own address.

A Matlab code was built which could read CT images and extract the information about those pixel values that describe the approximation line of the aorta cross-section. The pixel values were then transferred to Cartesian coordinate values. As the approximated contour is not in a circle form, the centroid were calculated by taking the mean of all x and y coordinate values. Using the formularadius=√

(X_avg −x)²+ (Y_avg −y)² the radius for all contours was calculated. Where (Xavg, Yavg) is the centroid, x and y are coordinate values for dierent pixels.

Figure 2.5: 7^th degree polynomial t to radius for rst 215 CT slices

The values of radius obtained had uctuations which were believed to be due to the real geometrical features of the aorta. A 7^th degree polynomial was tted to the radius data from 0to215 slices as shown in Figure 2.5 and then onwards up to393 slices were assumed to have constant values. The sudden rise in the values from 130 to 200 slices is due to the structural anatomy of the real aorta. In order to preserve the monotonic trend in the radius, some sets of radius values between0to215slices were also assumed to be constant. Now using these radius values, circles can be generated which represent the contours of the real aorta. Since the ideal geometry was assumed to be straight i.e. without bending, the distance considered between two successive slices was the curve length. In 2D view the stack of slices are imposed on each other. In practice, the distance c in Figure 2.6 is the distance between the centers of any two consecutive slices given by c= √

(x−x )²+ (y−y )². The parameter h is the cross-sectional

(23)

distance between the two slices considered at the time of capturing the CT images. The curve lengths between any two successive slices were calculated by d=√

c²+h².

Figure 2.6: Distances between two consecutive CT slices

Now from the calculated values of radius and consecutive distances between two slices, a three dimensional geometry can be constructed which has circular cross sections. Since the geometry looked planar and axisymmetric, the computational domain was considered as the symmetric part of the aorta geometry. With the help of a journal le containing coordinates, a 2D axisymmetric geometry of descending aorta was built automatically in ANSYS Gambit as shown in Figure 2.7.

Figure 2.7: 2D axisymmetric geometry of descending aorta

(24)

3 Computer simulations of aortic blood ow

In this chapter, we present mathematical methods and physical conditions that were applied during the experimental set up of this work. This includes the general algorithm of a CFD method following the grid generation, uid properties, governing equations and boundary conditions, solver set-up and the computational method for solution.

Fluid dynamics is the science of the uid motion. Fluid ow is studied in one of the three ways (1) Theoretical uid dynamics: This is crucial in understanding concepts and predicting the trends. It can obtain a lot of information using simplifying assump- tions. However, it does not always provide sucient information. (2) Experimental uid dynamics: This is costly and a dicult way to achieve exact conditions. Because it is not usually possible to develop the real situations due to the human errors or due to limit of the repetitions. (3) Numerically (computational uid dynamics (CFD)): This provides a qualitative prediction of uid ow with any conditions. With the use of computers, it has potential to provide tremendous amount of data at a fraction of the cost of experiments [11].

Fluid ows are governed by Partial Dierential Equations (PDE), which represents the conservation laws for mass, momentum, and energy. CFD is the art of replacing such PDE systems by a set of algebraic equations and then can be solved with the help of computers. CFD provides numerical approximation to the equations that govern the uid ow. We have used CFD approach for our study. The main steps in this process are:

1) Pre-processing : This stage consists of input to the ow problem. It involves

• Denition of the computational domain

• Grid generation- dividing domain into a nite number of small sub-domains

• Denition of uid properties and specication of appropriate boundary conditions 2) Solving : This stage solves the problem using an approach from dierent numerical methods like Finite Dierence Method(FDM), Finite Volume Method(FVM), Finite Element Method(FEM). The general steps of these numerical algorithms consists of

• Integration of the governing equations

• Discretization of the integral equations into a system of algebraic equations

(25)

3) Post-processing : In this stage the calculated results are visualized and interpreted by means of

• Vector plots and contour plots

• Streamlines and particle tracking

• 2D and 3D surface plots

3.1 Grid generation

Grid or mesh generation is the initial step in solving process where the computational domain is divided into a nite number of smaller sub-domains called elements or cells over which governing equations of the ow are then solved. It is the discrete representation of the geometry of the problem.

In this work, CFD Pre-processor ANSYS Gambit was used for meshing the geometries.

It can create structured and unstructured meshes in both 2D and 3D geometries that can be used for computational purposes. Structured meshes can be used in a relatively simple and regular geometry, while unstructured mesh is used in complex geometries.

Grid sizes have signicant impact on the rate of convergence, solution accuracy and the required processing time. A good solution can depend upon the properties of meshing like adequate denser grid, use of boundary layer mesh, adjacent cell length ratios etc.

So in order to improve the quality of the mesh, essential care must be taken in choosing the schemes.

For a 2D axisymmetric geometry of the descending aorta, we have used unstructured mesh elements for calculating the solution. Since we are interested in the calculation of wall shear stresses, the boundary layer mesh was employed on boundaries of the geometry and triangular elements were used elsewhere. The mesh elements were kept ner close to the boundaries and coarser away from it.

Figure 3.1: a boundary layer mesh in 2D axisymmetric model

(26)

Dierent mesh sizes were selected to ensure that the mesh we use for solution is grid independent. Wall shear stresses from seven dierent mesh cases were compared with the number of elements in it starting from 19,948 to131,824. As shown in Figure 3.2, wall shear stress for a 2D idealized geometry with dierent mesh sizes is compared. Results dier very slightly from each other with a lower to higher mesh sizes. In the mesh cases containing number of elements higher than 80000, no noticeable dierences are found in the solution and are almost equal. So we have used the mesh that contained about 81,188 numbers of computational mesh elements for our nal solution. In addition, ner meshes above certain level do not give signicant dierences in the accuracy of the solution; they only increase the computational time.

Figure 3.2: Wall shear stress in a 2D axisymmetric geometry from dierent mesh sizes For a 3D geometry of descending aorta, while smoothing and volume conversion process the edges of the geometry automatically rounds o. So the edges of inlet and outlet of the geometry were made at in ANSYS Gambit for proper inow and outow boundary conditions. Then as similar to 2D axisymmetric geometry, the unstructured mesh was generated with boundary layer mesh at the walls of the geometry as shown in Figure 3.3.

The height of the boundary was kept monotonically increasing from the wall towards the center of the geometry. From the computational-time point of view, it was not possible to try several dierent mesh cases with a 3D geometry. Therefore, an optimal mesh size with 799560 numbers of tetrahedral computational cells was chosen which gave reasonably good results within available time.

(27)

Figure 3.3: Mesh in cross sectional inow plane of a 3D descending aorta geometry

3.2 Fluid properties

A uid is any substance, which ows because its particles are not rigidly attached to each other. Fluid includes both liquids and gases. In the study of uid behavior, it is essential to have a proper understanding of the properties of the uid. Several properties of uids include temperature, pressure, density, viscosity, etc. This section briey describes the properties of uid that were considered in our ow problem.

The working uid considered in the aorta geometries was blood. It is the principal vehicle and the medium that provides nutrients and removes waste products from body organs. It consists of plasma uid, which contains about 90−95 percentage of water and dissolved materials like proteins, lipids, carbohydrates. The density of the blood was taken to be ρ = 1060kg/m³. The blood was assumed as an incompressible and Newtonian uid. A Newtonian uid is a uid whose stress versus stain rate curve is linear. The constant of proportionality is called viscosity. Blood is approximately four times more viscous than water. It does not exhibit constant viscosity and is non- Newtonian especially in the microcirculatory system. This happens due to presence of the red blood cells, which are small semisolid particles who increase the viscosity of the blood and aect the behavior of the uid [4]. Generally in large and medium size arteries blood behaves as a Newtonian uid and it is used as a common assumption in blood ow analysis [12]. Hence, the blood was assumed a Newtonian uid with constant viscosity µ = 0.0035P a. The typical Reynolds number value was calculated by R_e = ρV D/µ, whereρis the density of blood,V is the velocity at the entrance plane,Dis the diameter of the aorta at the entrance, µ is the viscosity of the blood, and Re was found in the range associated with laminar ow.

(28)

3.3 Governing equations and Boundary conditions

In this section, we will briey derive the governing equations and the boundary conditions that were used to solve the problem. This includes the derivation of Navier-Stokes equation following the discussion of parabolic and pulsatile ow proles for the inlet boundary conditions.

In a uid ow, the governing equations are the system of equations that describes the uid motion. These equations are mostly described by conservation laws of mass, momentum, and energy.

3.3.1 Conservation laws

The three-dimensional mass conservation or continuity equation for a compressible uid is given by equation

∂ρ

∂t +div(ρu) = 0 (3.1)

However, for an incompressible uid, the density ρis constant and hence equation (3.1) becomes

div(ρu) = 0 (3.2)

The compact form of three-dimensional momentum equation in given by equation

ρDu

Dt =−∇ ·p+∇ ·τ_ij +f (3.3)

where ρ^Du_Dt is the acceleration term, −∇ ·p is pressure gradient, τij are shear stresses that becomes viscosity for an incompressible Newtonian uid, and f are the body forces acting on the uid.

3.3.2 Navier-Stokes equations for Newtonian uid

In many uid ows, the viscous stresses can be expressed as a function of local deformation rate or the strain rate. In three-dimensional ows, the local rate of deformation is composed of linear deformation rate and volumetric deformation rate. The linear

(29)

deformation rate has nine components three of which are linear elongating deformation components dened by s_xx = ^∂u_∂x, s_yy = ^∂v_∂y, s_zz = ^∂w_∂z and and six shearing linear deformation components dened by:

s_xy =s_yx= 1 2

(∂u

∂y + ∂v

∂x )

s_yz =s_zy = 1 2

(∂v

∂z +∂w

∂y )

s_xz =s_zx = 1 2

(∂w

∂x +∂u

∂z )

The volumetric deformation rate is given by ^∂u_∂x+^∂v_∂y+^∂w_∂z =div u. In a Newtonian uid viscous stresses are proportional to the rates of deformation. Therefore, compressible uid involves two viscosities: rstµto relate linear deformation and secondlyλto relate volumetric deformation.

τ_xx = 2µ∂u

∂x +λ div u τ_yy = 2µ∂v

∂y +λ div u τzz = 2µ∂w

∂z +λ div u τ_xy =τ_yx =µ

(∂u

∂y + ∂v

∂x )

τ_yz =τ_zy =µ (∂v

∂z + ∂w

∂y )

τ_xz =τ_zx=µ (∂w

∂x + ∂u

∂z )

(3.4)

Substituting the values of equation (3.4) in equation (3.3) and rearranging the terms, we get

ρDu

Dt =−∂p

∂x +div(µ grad u) +S_{M x} ρDv

Dt =−∂p

∂y +div(µ grad v) +S_{M y} ρDw

Dt =−∂p

∂z +div(µ grad w) +S_{M z}

(3.5)

(30)

These equations are called Navier- Stokes equations in three dimensions for a compressible Newtonian uid.

In this work, the uid (blood) was assumed incompressible. Therefore, the density of blood remains constant as per equation (3.2) i.e. div u= 0. As we have not considered gravity to be acting on the uid, the term from body forces disappears and the stress term becomes viscosity.

ρDu

Dt =−∇ ·p+∇ ·µ (3.6)

Finally, the above form of Navier-Stokes equation was used by the commercial code ANSYS Fluent V6 to solve the ow problem.

3.3.3 Boundary conditions

Boundary conditions are a required component of a mathematical model of uid ow. In numerical simulations, it is essential to dene the boundary conditions for the selected region of the problem, as one might not be interested in simulating the entire area.

When solving Navier-Stokes equations, initial and boundary conditions must be applied.

Inappropriate boundary conditions might result in non-physical solutions.

Common boundary conditions are classied either in terms of numerical values that have to be set or in terms of physical conditions. For steady state problems, three types of spatial boundary conditions can be specied: Dirichlet, Neumann, and mixed boundary conditions. Whereas physical boundary conditions commonly used in uid problems are solid walls, inlets, symmetry boundaries, pressure boundary conditions, outow boundary conditions, cyclic or periodic conditions, etc. Boundaries within the uid can be dened as walls to bound uid and solid regions.

Inlet is the entrance region where the uid enters the domain, therefore its velocity, pressure, and mass ow rate must be known. If the uid has turbulent characteristics then these must be specied. Symmetric boundaries are applied when the ow is sym- metrical about some plane and this is commonly used to reduce the computational eort in a problem. Pressure gradient drives the uid in conned ows; therefore, pressure boundary conditions are generally used boundary conditions. Outow boundary conditions are employed at ow exits where the details of ow velocity or pressure are known prior to solution of the problem. This condition is appropriate to use when the ow at exit is almost in a fully developed state.

(31)

Figure 3.4: Boundary conditions for (a) 3D (b) 2D axisymmetric geometries of descending aorta

For our simulations, we have used a velocity inlet and outow condition at the inlet and the outlet respectively as shown in Figure 3.4. Assuming the blood ow to be fully developed at inow region, a parabolic velocity prole was imposed both on a 3D geometry of the descending aorta and a 2D axisymmetric idealized version of it for a steady state simulation. For this purpose, two UDF's (User Dened Function) were written in C compiler that demonstrate the parabolic nature of velocity using the relation given in equation (3.7) and (3.8).

u_x =u_max (

1−(y²+z²)¹² (radius)²

)

,for a 3D parabolic velocity prole (3.7)

u_x =u_max (

1− y² (radius)²

)

,for a 2D parabolic velocity prole (3.8)

(32)

Figure 3.5: Parabolic velocity proles at inlet of (a) 3D (b) 2D axisymmetric geometries of descending aorta

A transient simulation requires a time dependent velocity prole. So similar to the steady state simulation case, two more UDF's were written describing a parabolic velocity prole that changes with time as shown in Figure 3.6.

Figure 3.6: Velocity pulse function

The above velocity pulse was estimated by a 9^th degree polynomial and was used as a velocity inlet boundary condition. One cardiac cycle was0.85secondslong. Simulations were run for four cardiac cycles i.e 3.4seconds. During systolic pressure (the maximum pressure in the aorta when the heart is beating and pumping blood throughout the body) i.e. around 0.3^th−0.4^thsecond of time of cycle, the blood velocity is at its peak value. Whereas during diastolic pressure (the lowest pressure in the aorta in the moment between beats when the heart is resting ) the blood velocity falls.

(33)

3.4 Solver set-up

All the simulations were performed in ANSYS Fluent V6. The walls of the descending aorta models were assumed rigid and stationary. A pressure based algorithm was chosen as the solver type. This solver is generally selected for an incompressible uid. As there is no heat transfer in the blood ow process, energy equations were not solved.

Since turbulence was not expected in the aortic blood ow, a laminar model was used throughout the work. In solution methods, the SIMPLE algorithm was selected for pressure-velocity coupling. First Order Upwind scheme was employed as a numerical scheme for discretization of the momentum equations. Convergence criteria for all the ow parameters like continuity, pressure and velocity were set to 10⁻⁵. In transient cases, the time step size was set to 0.01sec with 340 number of total time steps. And maximum 500 iterations were performed per each time step.

3.5 Computational method

Most commonly used approach in CFD programs are based on either the Finite Volume Method (FVM) or the Finite Element Method (FEM). The other numerical methods include the Finite Dierence Method (FDM), the Spectral Method, the Boundary El- ement Method, etc. ANSYS Fluent V6 uses FVM to solve the governing equations of uid ow.

The Finite Volume Method is one of the most versatile discretization techniques used in CFD. It is a numerical method representing partial dierential equations in the form of algebraic equations. In the rst step, this method divides the domain into a nite number of smaller sub-domains called elements or control volumes, where the variable of interest is located at the center of the control volume. Interpolation determines the variable on the faces of the control volumes. Then it integrates the dierential form of the governing equations over each control volume. Finally, FVM uses an iterative method like the Gauss-Seidel method to solve the resulting system of algebraic equations known as discretized equations. One of the important features of FVM is its conservation properties. Since they are based on applying conservation principles over each control volume, global conservation is ensured. Another advantage of FVM is that the method is not limited to a single grid type. It gives freedom to use both structured and unstructured meshes.

(34)

4 Results and discussion

In this chapter, we present the simulation results of blood ow from both geometries of the descending aorta: real three-dimensional and two-dimensional axisymmetric. Veloc- ity proles at several axial cross-sections, pressure distribution on aorta wall, and WSS distribution from both steady state and unsteady state simulations are discussed in the following sub-sections.

An unsteady state simulation is a time dependent process. Depending upon the time and the inow boundary condition, the quantities like WSS, velocity, and pressure change their magnitudes. Thus, ow might behave dierently at dierent time steps. Dierent cases at every 0.1second of a total 3.4seconds simulation time were captured. Obser- vations from six dierent time steps (as shown in Figure 4.1) of the fourth cardiac cycle are presented and discussed in this work.

Figure 4.1: Dierent time steps in a cardiac cycle of inlet velocity prole

4.1 Wall shear stress distribution

WSS is an important factor in determining possible sites for development of arterial disease such as atherosclerosis ([14],[15],[16]). WSS depends upon the viscosity of the uid and the velocity gradient and, is dened as

τ =µdu

dr (4.1)

(35)

Figure 4.2: WSS distribution in a 2D axisymmetric geometry in steady state simulation Figure 4.2 shows the distribution of WSS along the wall of a 2D axisymmetric geometry in steady state simulation. Here, we can observe that the WSS is uniform up to 0.12m curve length of the geometry. A sudden climb to the maximum value of WSS about 2.80P a can be noticed in the narrow tapering region between 0.125−0.175m of the geometry. Then for the remaining curve length from 0.175−0.2459m, WSS declines and stabilize towards the outlet.

Figure 4.3: WSS (P a) distribution in a 3D geometry in steady state simulation

(36)

Figure 4.3 shows the distribution of WSS in a 3D geometry in steady state simulation.

Here, WSS is non-uniformly distributed. Dierent contours can be seen with higher WSS at the bending of the geometry. Few large jumps in the contours value were also observed within a small area of region near the curvature bend. The maximum value of WSS was about 8.72P a.

In unsteady simulation, the WSS magnitude and distribution vary during the systole and diastole of the cardiac cycles in a 2D axisymmetric geometry. The average WSS at six dierent time steps is given in Table 4.1. It is observed that the average WSS is higher at the time t= 0.25,0.35s that occurs during the systole of the cycle, when the blood is introduced to the aorta with maximum pressure. While the WSS is low during diastole of the cycle, when the heart is in the resting position.

Table 4.1: Average WSS in a 2D axisymmetric geometry at dierent time steps of a cardiac cycle

Time step Average WSS

(s) (P a)

t=0.15 0.0453 t=0.25 0.8188 t=0.35 1.3585 t=0.45 0.1930 t=0.55 0.7053 t=0.75 0.0603

A graphical representation of the average WSS in 2D axisymmetric geometry along the curve length at dierent time steps is presented in Figure 4.4. The points in the graph in Figure 4.4 represent the average WSS values measured from20equally divided sections.

There is not much variation in WSS at the time t = 0.15, 0.45, and 0.75s along the curve length. This is due to the fact that input ow proles at these times steps were at and zero. At the time t= 0.25and0.55s, a little higher uniformly distributed WSS can be noticed. During the time t = 0.35s, high WSS is observed. The distribution at this time step is uniform up to0.125m. After this point, WSS elevates to the maximum value about 2.63P a and then stabilizes near the outlet.

(37)

Figure 4.4: Average WSS distribution in a 2D axisymmetric geometry at dierent time steps

Figure 4.5: WSS (P a)distribution in a 3D geometry at dierent time steps

(38)

WSS distribution on aorta wall of a 3D geometry at dierent time steps is presented in Figure 4.5. Here the distribution is observed to be non-uniform at all time steps. At time t= 0.15 and 0.75s, the aorta wall experiences very low WSS. The distribution of WSS is fairly uniform in the upper half of the vessel. The lower half region experiences a high WSS with the highest value noticed at the largest curvature bend of the aorta wall. At time t= 0.35s, maximum WSS can be observed. Here the maximum value of WSS reached about 6.65P a. The average values of WSS in 3D geometry at dierent time steps are presented in Table 4.2.

Table 4.2: Average WSS in a 3D geometry at dierent time steps Time step Avergae WSS

(s) (P a)

t=0.15 0.0651 t=0.25 0.9061 t=0.35 1.4393 t=0.45 0.3859 t=0.55 0.5769 t=0.75 0.1022

4.2 Pressure distribution

The distribution of pressure in both geometries of the descending aorta in steady state simulation is shown in Figure 4.6. The pressure is uniform in the beginning of a 2D axisymmetric geometry as shown in Figure 4.6 (a). Near the contracting region, pressure uniformly decreases as blood ow moves ahead passing through the tapering of the geometry. In 3D geometry, the pressure is non-uniformly distributed on the aorta wall.

Lower pressure values were observed at the inner bending in comparison with the outer bending. This pressure dierence may provide some kind of force that may drive the blood through the bend [17].

(39)

Figure 4.6: Distribution of pressure(mmHg)in (a) 2D axisymmetric (b) 3D geometries of the descending aorta in steady state simulation

Figure 4.7: Pressure (mmHg) distribution in a 2D axisymmetric geometry at dierent times steps

(40)

The distribution of pressure in the 2D axisymmetric geometry at dierent time steps is presented in Figure 4.7. It can be observed that the pressure varies with time. The pressure contour changes uniformly from inlet towards outlet in all the cases. At time t = 0.15 and 0.75s, the range of pressure values can be observed to be smaller in comparison with the other times steps that makes it look uniform without any dierences in the contour. At time t = 0.25 and 0.35s, the variation in the pressure values can be noticed with high pressure at the inlet and low towards the outlet. The pressure is reversed at time t= 0.45and 0.55s with higher pressure at the outlets. Pressure at the inlet and the outlet at dierent time instants are presented in Table 4.3.

Table 4.3: Pressure values at the inlet and the outlet in a 2D axisymmetric geometry at dierent time steps

Time step Inlet Pressure Outlet Pressure

(s) (mmHg) (mmHg)

t=0.15 101.6421 101.5695

t=0.25 101.4307 80.0436

t=0.35 101.5461 89.9236

t=0.45 101.8516 119.8359

t=0.55 101.7500 113.8971

t=0.75 101.6485 102.4639

The pressure values at the inlet and the outlet in a 3D geometry at dierent time steps are presented in Table 4.4.

Table 4.4: Pressure values at the inlet and the outlet in a 3D geometry at dierent time steps

Time step Inlet Pressure Outlet Pressure

(s) (mmHg) (mmHg)

t=0.15 101.7145 101.5209

t=0.25 101.4521 80.2589

t=0.35 101.5532 89.8078

t=0.45 101.9515 119.2064

t=0.55 101.8550 113.5056

t=0.75 101.7248 102.3860

Figure 4.8 illustrates the distribution of pressure in 3D geometry at dierent time instants. The variations in the pressure contours can be observed at dierent time steps.

(41)

The pressure contours are uniformly distributed in the upper half region of the geometry, whereas below the tapering region, the distribution of pressure is non-uniform.

At the time t = 0.15 and 0.75s, the pressure is uniform and varying with a small variation throughout the geometry. Large pressure dierence was noticed at time t = 0.25s. During the peak of the systole i.e.t = 0.35s, the geometry experiences high pressure. Reversed pressure can be observed at the time t = 0.45 and 0.55s. The pressure values were observed low near curvature bending of the geometry with few pressure dierence contours at all the time steps. The highest-pressure can be noticed at the outlet during time t = 0.45s of the cycle.

Figure 4.8: Pressure (mmHg) distribution on the wall of a 3D geometry at dierent time steps

(42)

4.3 Velocity distribution

The velocity ow in a 2D axisymmetric geometry in steady state simulation is presented in Figure 4.9. A parabolic velocity prole at the inlet with a maximum value of0.45m/s is carried until the mid of the geometry. This is because of almost uniform diameters and negligibly small dierences of the axisymmetric geometry caused by the smaller tapering coecient. High velocity gradients can be observed near the tapering region of the geometry. The velocities at the walls were zero. After passing through the tapering region, the velocity increases and reaches the maximum value of 0.677m/s close to the outlet.

Figure 4.9: Velocity(m/s) distribution in a 2D axisymmetric geometry in steady state simulation

The velocity proles at several axial planes are presented in Figure 4.10 (b). These pictures describe the distribution of velocity magnitudes in 3D geometry of the descending aorta. From Figure, it is observed that parabolic ow at the inlet begins with a maximum velocity 0.45m/s distributed at the center and zero velocity on the walls. It can be observed that in sections (s4, s5, s6), the ow accumulates on outer curvature of the geometry. This causes the velocity proles to be C-shaped. Similarly, the ow accumulates towards the outer curvature in the sections (s7, s8, s9). The highest velocity was about 0.652m/s near the largest curvature, where there is sudden tapering in the geometry. Later near the outlet, the ow regains parabolic shape.

(43)

Figure 4.10: (a) axial cross-sections in a 3D geometry (b) velocity (m/s) proles at dierent cross-sections of a 3D geometry normal to the ow direction in steady state simulation.

Figure 4.11: Velocity (m/s) distribution in a 2D axisymmetric geometry at dierent times steps

(44)

In unsteady state simulation, dierent velocity proles were applied at dierent time steps of a cardiac cycle. The velocity distribution in a 2D axisymmetric geometry at dierent time steps is presented in Figure 4.11. It can be observed that the velocity is low at time t = 0.15,0.55,0.85s as shown in Figure 4.11. The variation of velocity magnitude in these cases is small which is represented in the shades of color blue. The velocities at the walls were zero. The maximum velocity 0.48m/s is reached at time t = 0.35s that occurs during the systole. In the straighter region of the geometry, the distribution of velocity is uniform with the maximum velocity magnitude accumulated at the center of the ow. Near the tapering region of the geometry, the velocity magnitude rises with high gradients. The maximum velocities at dierent time steps are given in Table 4.5.

Table 4.5: Maximum velocity at dierent times steps in a 2D axisymmetric geometry Time step Maximum Velocity

(s) (m/s)

t=0.15 0.1026

t=0.25 0.2209

t=0.35 1.4818

t=0.45 0.3578

t=0.55 0.1429

t=0.75 0.1304

The distribution of velocity in axial planes in a 3D geometry at dierent times of a cycle is presented in Figure 4.12. The cross sections shown in Figure 4.12(a) were taken normal to the direction of ow in the geometry. During time t = 0.15,0.55,0.75s, the applied inlet velocity was very low, which makes it dicult to distinguish between the various contours of velocity magnitudes and hence appears at with the color blue.

During the time t = 0.25,0.45s, blood ows with the highest velocity. Maximum velocity was noticed at systole time t = 0.35s. At all the time steps, ow is skewed towards the outer curvature of the geometry, which can be observed in section z1. This causes velocity proles to be C-shaped. The velocity magnitudes is maximized towards the outer bending curvatures, which can be observed in sections z1, z2. Later near the outlet, the ow regains parabolic shape.

(45)

Figure 4.12: (a) Axial cross-sections in the 3D geometry (b) Velocity (m/s) proles at dierent cross sections at dierent time steps

4.4 Comparison of 2D and 3D geometries for wall shear stress, pressure and velocity distribution

Here we analyze and compare WSS, pressure, and velocity distribution in both the geometries at steady and unsteady conditions.

Figure 4.13 shows the distribution of the average WSS values in both the geometries in steady state simulation. The points in the graph represent the average values of 20 equally divided sections. It can be observed that WSS in a 2D axisymmetric case begins uniformly, while 3D case starts with a growing trend having comparatively higher WSS values. This happens because the region up to the curve length of0.1min 2D geometry is more regular. The corresponding length in 3D geometry is irregular because of the rough surface that causes non-uniform WSS distribution. At the region near the tapering, WSS jumps to the maximum with somewhat higher values in 3D geometry. After the sharp tapering, WSS declines smoothly in 2D case whereas some uctuations can be observed in 3D case due to its anatomical irregularities. In summary, the idealization of geometry leads to the underestimation of WSS though the trend of WSS variation is qualitatively

(46)

in agreement between the two geometries.

Figure 4.13: Comparison of 2D and 3D geometries for average WSS distribution in steady state simulation

Figure 4.14: Comparison of 2D and 3D geometries for average WSS distribution at time t=0.15,0.25,0.35,0.45,0.55,0.75 s.

(47)

Figures 4.14 illustrates WSS distribution at dierent time steps in 3D (solid lines) and 2D axisymmetric (dotted lines) geometries of the descending aorta. The points in the graph represent the average WSS values of 20 equally divided sections from both the geometries. At the time steps t = 0.15 and 0.75s, WSS along the curve length of both the geometries is low and the overall distribution is uniform. The WSS increases with rising inlet velocity magnitude. At the timet = 0.25s, both geometries experience high WSS. Maximum WSS can be observed during the systole att = 0.35s. An elevation in the magnitude of WSS can be observed in the tapering region 0.125−0.175m length of both the geometries. Sudden reduction of the ow area in this region causes velocity gradients resulting in high WSS. This increment in 2D axisymmetric geometry can be observed linear, whereas 3D geometry has uctuations.

Figure 4.15 shows the distribution of average pressure in both the geometries in steady state simulation. The points in the graph represent the average pressure values of 20 equally divided sections from both geometries. The distribution in 2D axisymmetric geometry begins uniformly with comparatively higher pressure values. The corresponding distribution of pressure in 3D geometry is low with uctuations. After passing a distance of 0.12m, pressure declines smoothly in 2D case whereas the corresponding pressure in 3D geometry has a few oscillations and ends with higher pressure values.

The pressure dierence in 3D and 2D axisymmetric geometries were 27.67mmHg and 38.41mmHg, respectively.

Figure 4.15: Comparison of 2D and 3D geometries for average pressure distribution in steady state simulation

(48)

Figure 4.16 illustrates the average pressure distribution in 3D and 2D axisymmetric geometries at dierent time steps. At the timet = 0.15,0.75s, the pressure distribution in both geometries can be observed to be uniform. During the time t = 0.25,0.35s, pressure decreases in both geometries. While at time t = 0.45,0.55s, reversed pressure was observed. At these time steps, the distribution moves from lower (inlet) to higher- pressure (outlet) values. Highest pressure can be seen at the time t = 0.45s near the outlet of both the geometries. Pressure distribution in 2D axisymmetric geometry can be observed to be little smoother than in 3D geometry. This is because of the rough surface walls in real geometry.

Figure 4.16: Comparison 2D and 3D geometries for average pressure distribution at time t=0.15,0.25,0.35,0.45,0.55,0.75 s.

Figure 4.17 illustrates the pressure drop at dierent time steps of a cycle in 3D and 2D axisymmetric geometries. At time t = 0.15,0.75s, the pressure drop in both the geometries is measured to be small. A noticeable pressure dierences of approximately 12mmHg was observed at timet = 0.35,0.55s in both geometries. High pressure drops of approximately 21mmHg and 17mmHg can be observed at the time t = 0.25,0.35, respectively. Comparatively, the pressure dierences in 2D axisymmetric geometry were observed to be a little higher than in 3D geometry.

Pulsatile Blood Flow Simulations in Computed Tomography(CT) Scan-Based and Idealized Geometries of Human Aorta