Fluid-Structure Interaction Analysis of Carotid Artery Blood Flow A Patient-Specific Investigation of the Impact of Arterial Wall Deformation on Hemodynamics Master’s thesis in Biomedical Engineering SIMON NORDENSTRÖM DEPARTMENT OF MECHANICS AND MARITIME SCIENCES CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2022 www.chalmers.se www.chalmers.se Master’s thesis 2022 Fluid-Structure Interaction Analysis of Carotid Artery Blood Flow A Patient-Specific Investigation of the Impact of Arterial Wall Deformation on Hemodynamics SIMON NORDENSTRÖM Department of Mechanics and Maritime Sciences Division of Fluid Mechanics Chalmers University of Technology Gothenburg, Sweden 2022 Fluid-Structure Interaction Analysis of Carotid Artery Blood Flow A Patient-Specific Investigation of the Impact of Arterial Wall Deformation on Hemodynamics SIMON NORDENSTRÖM © SIMON NORDENSTRÖM, 2022. Supervisor: Dr. rer. nat. Jan Brüning, Institute of Computer-Assisted Cardiovas- cular Medicine, Charité – Universitätsmedizin Berlin Supervisor: Prof. Leonid Goubergrits, Institute of Computer-Assisted Cardiovascu- lar Medicine, Charité – Universitätsmedizin Berlin Examiner: Prof. Henrik Ström, Department of Mechanics and Maritime Sciences, Chalmers University of Technology Master’s Thesis 2022:11 Department of Mechanics and Maritime Sciences Division of Fluid Mechanics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Time-averaged wall shear stress in a patient-specific carotid artery simulated using a Mooney-Rivlin material model. The dark blue colored region is associated with a high risk of atherosclerotic development. Typeset in LATEX Printed by Chalmers Reproservice Gothenburg, Sweden 2022 iv Fluid-Structure Interaction Analysis of Carotid Artery Blood Flow A Patient-Specific Investigation of the Impact of Arterial Wall Deformation on Hemodynamics SIMON NORDENSTRÖM Department of Mechanics and Maritime Sciences Chalmers University of Technology Abstract Cardiovascular disease, predominantly caused by atherosclerosis, is the leading cause of death worldwide. Therefore, scientific methods capable of accurately evaluating risk factors of cardiovascular disease on a patient-specific basis are highly sought- after. Typically, numerical approaches for simulating blood flow assume arterial walls to be rigid structures. However, to obtain more physiologically realistic simu- lations, the interplay between wall deformation and blood flow, the so-called fluid- structure interaction (FSI), must be taken into consideration. In this project, a workflow is established for performing patient-specific FSI simulations of blood flow and arterial wall deformation in the carotid artery – the main source of blood sup- ply to the brain. The workflow is primarily intended to aid in the execution of FSI analyses in the STAR-CCM+ simulation software. Furthermore, the impact of FSI on hemodynamic parameters such as flow velocity, wall shear stress, and oscillatory shear index, is quantified and analyzed for a patient-specific geometry. It is shown that a rigid wall assumption results in underestimations of the areas associated with an increased risk of atherosclerosis, namely regions of low time-averaged wall shear stress and high oscillatory shear index. Moreover, the rigid wall assumption leads to overestimations of the flow velocities during the systolic phase of the cardiac cy- cle. FSI analyses are performed using three different material models: an isotropic linear elastic model, a Neo-Hookean model, and a Mooney-Rivlin model. For the specific boundary conditions and parameter values employed, marginal differences are observed in the hemodynamics between the three material models. Keywords: FSI, Carotid artery, Simulation, Hemodynamics, Workflow, STAR-CCM+. v Acknowledgements I would like to thank my two thesis supervisors, Dr. Jan Brüning and Prof. Leonid Goubergrits, for the project conceptualization and data acquisition, along with their expertise and valuable input throughout the project. I would also like to express my gratitude to all my colleagues at the Institute of Computer-assisted Cardiovascular Medicine for making the thesis an enjoyable experience. Simon Nordenström, Berlin, May 2022 vii List of Acronyms CAS Carotid Artery Stenting CCA Common Carotid Artery CEA Carotid Endarterectomy CFD Computational Fluid Dynamics ECA External Carotid Artery FEM Finite Element Method FSI Fluid-Structure Interaction FVM Finite Volume Method GCI Grid Convergence Index HGO Holzapfel-Gasser-Ogden ICA Internal Carotid Artery ICM Institute of Computer-Assisted Cardiovascular Medicine ILE Isotropic Linear Elastic LHS Left-Hand Side MRI Magnetic Resonance Imaging NO Nitric Oxide OSI Oscillatory Shear Index PCMRI Phase Contrast MRI PDE Partial Differential Equation PISO Pressure-Implicit with Splitting of Operators RHS Right-Hand Side SCL Space Conservation Law SIMPLE Semi-Implicit Method for Pressure Linked Equations TAWSS Time-Averaged Wall Shear Stress WSS Wall Shear Stress ix Contents 1 Introduction 1 1.1 Aim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Theory 5 2.1 Carotid artery anatomy, physiology and pathology . . . . . . . . . . . 5 2.2 State of the art in carotid artery FSI simulations . . . . . . . . . . . 8 2.3 Blood rheology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.4 Fluid flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.5 Vessel wall mechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5.1 Isotropic linear elastic model . . . . . . . . . . . . . . . . . . . 12 2.5.2 Isotropic non-linear hyperelastic models . . . . . . . . . . . . . 13 2.5.3 Anisotropic non-linear hyperelastic model . . . . . . . . . . . 14 2.6 Fluid-structure interaction . . . . . . . . . . . . . . . . . . . . . . . . 15 2.7 Hemodynamic parameters . . . . . . . . . . . . . . . . . . . . . . . . 16 2.8 Fourier decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.9 Richardson’s extrapolation . . . . . . . . . . . . . . . . . . . . . . . . 17 2.10 Analytical displacements in cylinders . . . . . . . . . . . . . . . . . . 17 3 Methods 19 3.1 Solution strategy for the fluid flow . . . . . . . . . . . . . . . . . . . . 19 3.2 Solution strategy for the solid displacements . . . . . . . . . . . . . . 21 3.3 Fluid-structure coupling . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.4 Preprocessing of the carotid artery geometry . . . . . . . . . . . . . . 24 3.5 Meshing of the domains . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.6 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.7 Simulation operations and evaluation . . . . . . . . . . . . . . . . . . 28 3.8 Models and parameter values . . . . . . . . . . . . . . . . . . . . . . 29 4 Results 31 4.1 Overall workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2 Mesh sensitivity study . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.3 Hemodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.4 Wall deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5 Discussion 39 5.1 Workflow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 xi Contents 5.2 Hemodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.2.1 General results . . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.2.2 Impact of FSI and material models . . . . . . . . . . . . . . . 40 5.3 Wall deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 5.4 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.5 Ethical and societal considerations . . . . . . . . . . . . . . . . . . . 44 6 Conclusion 45 Bibliography 47 xii List of Figures 2.1 Anatomy of the carotid arteries. . . . . . . . . . . . . . . . . . . . . . 6 2.2 Generic diagram of a healthy elastic artery . . . . . . . . . . . . . . . 6 3.1 Carotid artery geometry. . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Preprocessing pipeline. . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Superior view of the fluid and solid meshes. . . . . . . . . . . . . . . 27 3.4 Boundary conditions imposed at the inlet and outlets . . . . . . . . . 28 3.5 Points selected for evaluative purposes. . . . . . . . . . . . . . . . . . 29 4.1 Final workflow for performing FSI analyses of carotid arteries. . . . . 32 4.2 CFD mesh independence analysis of TAWSS and OSI risk zone areas 33 4.3 TAWSS from CFD and FSI . . . . . . . . . . . . . . . . . . . . . . . 34 4.4 WSS at peak systole . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.5 OSI from CFD and FSI . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.6 Velocity magnitudes at peak systole . . . . . . . . . . . . . . . . . . . 35 4.7 Velocity magnitude over time . . . . . . . . . . . . . . . . . . . . . . 35 4.8 Magnitude of wall deformation at peak systole . . . . . . . . . . . . . 36 4.9 Magnitude of wall deformation over time . . . . . . . . . . . . . . . . 37 xiii List of Figures xiv List of Tables 2.1 Overview of presented material models. The material parameters are defined later in this section. The complexity described in the last column is in relation to the other material models, and may be interpreted as the relative level of expected accuracy for large deformations. Of note, in addition to obtaining values of the material parameters, the HGO model also requires specifications of the fiber orientations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.1 Fourier coefficients of the boundary conditions. . . . . . . . . . . . . . 28 4.1 Richardson’s extrapolation of fluid domain mesh values. Note that the computation time presented in the last column is for one FSI time step, whereas the rest of the columns correspond to CFD values. . . . 33 4.2 Deviations in the estimated TAWSS and OSI risk areas . . . . . . . . 36 4.3 Deviations in the probed velocities at peak systole . . . . . . . . . . . 36 4.4 Deviations in the wall deformation magnitudes at peak systole. . . . . 37 xv List of Tables xvi 1 Introduction Atherosclerosis, the leading cause of death in the developed world, is the primary cause of most cardiovascular diseases and is signified by so-called atheromatous plaques on the inner walls of the blood vessels [1]. These plaques are build-ups of bodily substances such as cholesterol, calcium, and macrophages (a type of white blood cell). The presence of such plaques leads to a partial or complete obstruction of the blood flow through the blood vessel. Especially severe occlusions of blood vessels may occur when a plaque ruptures, leaving fragments of it called emboli traveling through the cardiovascular system, which may lodge into smaller blood vessels. For instance, approximately 80 % of all strokes are ischemic, caused by emboli that occlude blood vessels of the brain and consequently cut off the oxygen and nutrient supply to parts of the brain [2]. The main arteries supplying the brain with blood are the carotid arteries, which have an estimated prevalence of plaque corresponding to 20 % of the adult population worldwide [3]. The high prevalence of carotid artery disease (i.e., narrowing of the carotid artery due to the presence of plaque) is thought to be partly due to the complex bifurcated geometry of the carotid arteries, which favors atherosclerotic progression [4]. Atherosclerosis is mainly an asymptomatic disease, presenting with symptoms only in the late, potentially al- ready life-threatening, stages of the disease progression. Moreover, atherosclerosis typically develops slowly over years or decades, with old age being one of the pri- mary risk factors for severe disease. With the marked demographic shift towards an older population in nearly every country across the globe, the demand for scientific tools and methods capable of analyzing atherosclerosis is ever-increasing. Simulations of blood flow are of great importance in cardiovascular research, as it may be used to identify and quantify risk factors for numerous cardiovascular dis- eases, including atherosclerosis. Hemodynamic parameters such as wall shear stress (WSS) and oscillatory shear index (OSI) on arterial walls have been proved to be associated with atherosclerotic progression. Regions with low WSS and oscillatory flow pose a risk for plaque development, whereas high WSS tends to lead to destabi- lization and plaque rupture [5], [6], [7]. Traditionally, blood flow has been simulated using computational fluid mechanics (CFD) approaches, treating arterial walls as rigid, fixed bodies. Although highly useful, approaches where arterial compliance is neglected use a major simplification of reality. Arteries expand during systole (contraction phase of the cardiac cycle) and store blood, which is released during diastole (relaxation phase of the cardiac cycle). Arterial compliance therefore has an important function in blood pressure stabilization. In consequence, stiff arteries are risk factors for and consequences of several cardiovascular diseases [8], [9]. 1 1. Introduction In order to obtain physiologically realistic simulations of arterial blood flow, fluid- structure interaction (FSI) must be taken into consideration. During systole, the increased pressure causes arteries to expand; at the same time, the induced geomet- rical changes of the arterial walls affect the way in which the blood flows through the arteries. Therefore, there is a constant mutual influence between the blood flow and the arterial wall deformation. Large arteries such as the carotid arteries generally expand their lumen diameter by 10 % or more during systole, altering the hemody- namics [10], [11]. To perform FSI simulations, a material model must be used to define the behavior of the arterial walls in the structural domain, while fundamental equations for the conservation of mass and momentum must be solved for the fluid domain. Moreover, a coupling must exist between the structural and fluid domains. Although a few previous studies in the field have compared differences in hemody- namic results between FSI and CFD in carotid artery blood flow, studies comparing hemodynamics and deformations between different material models in this context seem to be scarce [12]. This project is performed at the Institute of Computer-Assisted Cardiovascular Medicine (ICM) at Charité - Universitätsmedizin Berlin using the Star-CCM+ [13] software. Research groups at the ICM have great experience with modeling and simulation of arterial blood flow using CFD. However, their simulations have not taken FSI into account. The methodology and results from this thesis will serve to enable more accurate simulations of hemodynamics and arterial wall deformation at the ICM in the future. 1.1 Aim The main aim of the project is to establish a workflow for patient-specific FSI sim- ulations of blood flow in the carotid arteries. As a secondary aim, the study aims to partially fill the current gap in research regarding how different material models compare to each other, and to CFD, when studying the effects on hemodynamics and structural deformations in the carotid arteries using FSI. The issues that will be addressed in this study include the following: • How can the structural and fluid domains be solved separately? • How can a coupling between the two domains be established? • How can a workflow be developed that instructs the user on how to perform FSI simulations of carotid artery blood flow, that are also customizable to the needs of the user and the patient-specific conditions? 1.2 Limitations The following limitations will be made: • Only one, healthy, subject will be analyzed. Hence, no statistical analyses will be made. 2 1. Introduction • The patient-specific carotid artery geometry will not be generated in this project, but will be provided by the ICM. • The fluid flow is assumed to be laminar. This assumption has been made in close to all previous studies of carotid artery blood flow [12]. • Evaluation against patient-specific experimental results will not be performed, due to the complicated nature of conducting such experiments in vivo. • The material models used will be restricted to those available in Star-CCM+. If not specified otherwise, the limitations above were chosen in order to secure a sufficient amount of time available for fulfilling the aim of the project. 3 1. Introduction 4 2 Theory In this chapter, the theoretical underpinnings of this work are introduced. First, a medical introduction is presented in Section 2.1. In Section 2.2, the state of the art in performing FSI simulations of the carotid artery is reviewed. Thereafter, the theory required to simulate the blood flow is described in Sections 2.3 & 2.4. This is followed by a description of various models used for the deformation behavior of the vessel wall, in Section 2.5. The conditions which need to be satisfied for the FSI to be modeled are described in Section 2.6. Following this, the hemodynamic parameters selected for evaluation are presented in Section 2.7. In Section 2.8, Fourier series are described, as they will later be used to decompose the inlet and outlet waveforms. A method for estimating the uncertainty arising from a finite number of mesh cells is introduced in Section 2.9. Finally, analytical formulas for the deformation of ideal cylinders are presented in Section 2.10, as comparisons to running simulations. 2.1 Carotid artery anatomy, physiology and pathol- ogy The carotid arteries are undoubtedly one of the most important components of the cardiovascular system. Their most vital functions include the regulation of blood pressure and the insurance of a sufficient oxygen and nutrient supply to the brain. To the clinician, carotid arteries are of great interest partly due to the frequency and especially severe consequences of atheromatous plaques building up in the region. To researchers in fluid mechanics, it is the complex hemodynamics resulting from the intricate geometry of the carotid arteries that is of interest. In Fig. 2.1a, the anatomical origins of the common carotid arteries (CCAs) are shown. Blood is pumped from the left ventricle of the heart (not shown) into the as- cending aorta, which further subdivides into smaller arteries to deliver blood to the rest of the body. As can be seen, a difference exists between the left and right CCA: the left CCA branches out directly from the aortic arch, whereas the right CCA stems from the brachiocephalic (innominate) artery. The brachiocephalic artery originates from the aortic arch, and apart from the right CCA also branches off into the right subclavian artery, which primarily supplies the right upper limb with blood. The subclavian arteries have branches called the vertebral arteries that to- gether with the carotid arteries are the only suppliers of blood to the brain. Each of the two CCAs bifurcate into an internal carotid artery (ICA) and an external carotid artery (ECA), as can be seen in Fig. 2.1b. The ICA supplies blood to the 5 2. Theory brain, while the ECA is responsible for blood supply to the face and neck. (a) (b) Figure 2.1: Anatomy of the carotid arteries. (a) Origins of the right and left com- mon carotid arteries from the brachiocephalic artery and aortic arch, respectively [14]. (b) Common carotid artery bifurcating into the internal and external carotid arteries [15]. Figure 2.2: Generic diagram of a healthy elastic artery [16]. Three layers (tunicas) are apparent, namely the intima, the media and the adventitia. Like most other elastic arteries, i.e., large arteries in close proximity to the heart, the carotid arteries are composed of three distinct layers, as visualized in Fig. 2.2. The innermost layer, the intima, consists of a single-cell layer of endothelial cells which control and regulate inflammatory responses and the exchange of substances be- tween the blood flowing through the lumen and the cells of the surrounding tissues. In the endothelium, synthesis of nitric oxide (NO) takes place, which plays a crucial 6 2. Theory role in maintaining proper cardiovascular function [17]. It is believed that high WSS stimulates NO to vasodilate (expand) the artery by relaxing smooth muscle cells, increasing blood flow. This is one of many complex processes that, albeit protects against atherosclerosis, may increase the likelihood of plaque destabilization [18]. By the same token, low WSS (among other things) inhibits NO production, leading to favorable environments for plaque formation. The following layer, the tunica media, consists of smooth muscle cells, elastin and collagen fibrils, which together are ar- ranged in helical shapes. Elastin provides flexibility to the blood vessel and ensures that the original shape can be restored after pressure waves have passed through it. Collagen gives sturdiness to the vessel and prevents rupture from high pressures. The outermost layer, the adventitia, is collagen-rich but elastin-poor. The collagen and other constituents of the adventitia are, as in the media, structured in helices. Baroreceptors are present in the adventitia layer in the so-called carotid sinus – typically a dilated part of the ICA at the carotid bifurcation. These baroreceptors are sensitive to increases in pressure, and respond by stretching. In turn, increased baroreceptor stretch triggers a signal to the cardiovascular control center of the brain, that subsequently sends signals to the heart and other parts of the body to activate mechanisms for lowering the blood pressure. More specifically, the carotid sinus sends a signal through the glossopharyngeal nerve (cranial nerve IX) to the nucleus solitarius in the medulla oblongata in the brainstem, which then stimulates parasympathetic responses (e.g., decreases heart rate through increased vagal tone of the sinoatrial node of the heart; causes vasodilation) and inhibits sympathetic re- sponses [19]. Similarly, at low blood pressures, the baroreceptors get less stretched and therefore trigger sympathetic responses while inhibiting parasympathetic influ- ence, which causes an increase in blood pressure. These blood pressure stabilization mechanisms are called the baroreflex. Atherosclerosis has been shown to disturb the effectiveness of the baroreflex. This is thought to occur partially by increases in the neurotransmitter serotonin travelling through the vasa vasorum (i.e., the small blood vessels that supply the arteries themselves with blood, see Fig. 2.2) of the carotid artery, which then negatively impacts the nerve signals to and from the artery [20]. In addition to baroreceptors, important chemoreceptors are also present in the carotid artery. The chemoreceptors are located in the so-called carotid body, sit- uated close to the carotid apex, i.e., the “tip” between the ICA and ECA. The chemoreceptors monitor the pH as well as the partial pressure of oxygen and carbon dioxide in the blood. When, for instance, a decrease in partial pressure of oxygen is detected, suitable bodily responses including increased alveolar ventilation (i.e., augmented uptake of carbon dioxide from the blood to the lungs and release of oxy- gen in the opposite direction) are triggered [21]. The chemoreceptor response follows the same pathways to the cardiovascular control center in the brain as the baroreflex. Two main surgical procedures exist for treating carotid artery disease: carotid en- darterectomy (CEA) and carotid artery stenting (CAS). In CEA, the ICA is cut open, and the plaque inside the artery is removed. In cases where CEA is consid- 7 2. Theory ered too risky, CAS is instead performed. CAS is a procedure in which the plaque is not removed, but instead a stent is inserted at the location of the plaque using a catheter, which then opens up the blocked passage. As CAS and CEA themselves are associated with risks such as stroke, it is of great importance to evaluate po- tential risk factors. CFD and FSI studies have proved useful in evaluating both the risks associated with hemodynamics, and post-operative results [22] [23] [24]. 2.2 State of the art in carotid artery FSI simula- tions Fluid-structure interaction (FSI) simulations of carotid artery blood flow and wall deformations have advanced considerably in the past decades. Current state of the art in the field allows for simulations of patient-specific carotid artery geometries and boundary conditions, utilizing non-Newtonian rheology models as well as phys- iologically realistic constitutive models for the arterial walls. The effect of linear elastic arterial wall compliance on the hemodynamics in a healthy patient-specific carotid artery has been studied by Lopes et al. [25]. FSI simulations with various values of Young’s modulus (also called the elastic modulus), and estima- tions of hemodynamic parameters and displacements were subsequently compared to those obtained through CFD. The values of the Young’s moduli used ranged from 1.1 MPa to 2.2 MPa, as an attempt to simulate the impact of arterial wall stiffening (decreased elasticity) through atherosclerosis. A similar study has been conducted for stenosed carotid arteries, where instead the influence of the elastic modulus of the plaque on the fractional flow reserve was investigated [26]. Ahmadpour et al. [27] also employed a linear elastic model, aiming at quantifying several hemodynamic parameters of idealized stenoses of varying degrees of severity. To avoid shortcom- ings of the linear elastic assumption, an isotropic hyperelastic 5-parameter Mooney- Rivlin material model for the arterial wall was used by Zouggari et al. [28] in an FSI study that analyzed the hemodynamic differences between a stenosed and healthy carotid artery. In the study, hemodynamic parameters such as wall shear stress and velocities estimated at the plaque and at locations downstream of the plaque were compared for FSI and CFD approaches. Perhaps the to date most physiologically realistic simulations of carotid artery hemodynamics were carried out by Kamenskiy et al. [24]. In their study, Kamenskiy et al. performed FSI investigations of the impact of CAS on patient-specific geometries and hemodynamics in a cohort of 15 patients. The influential hyperelastic anisotropic Holzapfel-Gasser-Ogden (HGO) model [29] was employed, combined with patient-specific CCA angiocatheter pres- sures and ICA and ECA velocities obtained from Doppler ultrasound. A multitude of material models have been used in different FSI studies on carotid arteries, but single studies combining and investigating the impact of material model choice seem particularly difficult to come across. Possibly the only study of this kind was per- formed by Kumar et al. in 2021 [30], in which hemodynamic parameters simulated with a linear elastic model and a Mooney-Rivlin model were shown to differ from each other. 8 2. Theory State-of-the-art FSI studies in carotid arteries all share some common features, judg- ing from a recent review paper on the topic [12]. First, they assume a laminar flow regime. Second, they use a strongly coupled, two-way approach for solving the FSI problems. Typically, FEM-based solvers are used both for the solid and fluid domain, with ADINA [31] arguably being the most commonly employed software. Although STAR-CCM+ is a common software for CFD studies, carotid artery FSI simulations appear to have hardly been carried out in the software, potentially sig- naling a need for a workflow that facilitates such simulations. Third, they tend to lack validation of the results against experimental studies, which is one of the important aspects for coming FSI studies to explore. 2.3 Blood rheology Newtonian fluids like water have a constant viscosity, µ, meaning that the resistance to fluid flow is independent of the so-called shear rate γ̇. In layman’s terms, the viscosity may be understood as the thickness of the liquid, and the study of it is part of a branch in fluid mechanics called rheology. More precisely, it represents the internal friction that exists between layers of the fluid flowing at different ve- locities. Blood plasma, which makes up approximately 55 % of the blood volume, mainly consists of water and can therefore be considered Newtonian [32]. However, the remaining part of the blood has cellular components, most notably erythro- cytes, which show a non-Newtonian behavior. At high shear rates, i.e., when the difference in flow velocity between adjacent layers is large, the erythrocytes tend to deform which decreases the viscosity and facilitates flow. At low shear rates, in contrast, the erythrocytes aggregate and form so-called Rouleaux which results in a viscosity increase [33]. In other words, blood exhibits a shear-thinning behavior. Various models are available for estimating the shear rate dependent viscosities of non-Newtonian fluids. Common ones include the Power Law, Casson, and Carreau- Yasuda models. In the Carreau-Yasuda model [34] [35], the viscosity is given by µ(γ̇) = µ∞ + (µ0 − µ∞) (1 + (λγ̇)a)(n−1)/a , (2.1) where µ∞ is the viscosity at infinite shear rate, µ0 the viscosity at zero shear rate, λ the relaxation time, a a shear-thinning parameter, and n the power constant. 2.4 Fluid flow The governing equations for practically every type of fluid flow are the Navier-Stokes continuity and momentum equations [36]. The continuity equation is a statement of conservation of mass and is expressed as ∂ρ ∂t +∇ · (ρv) = 0, (2.2) 9 2. Theory where ρ and v are the fluid density and velocity vector field, respectively. ∇ is the mathematical nabla operator. For incompressible flows, the density is constant and the transient term in Eq. (2.2) is zero, i.e., ∇ · v = 0. (2.3) In other words, in incompressible flows, the divergence of the velocity is zero. The Navier-Stokes momentum equation describes the evolution of the velocity field of fluid flow, and can be regarded as being a formulation of conservation of linear momentum: ∂(ρv) ∂t +∇ · (ρv⊗ v) = ∇ · σ + fb, (2.4) where ⊗ denotes the outer product, σ is the stress tensor, and fb are body forces per unit volume caused by e.g., gravity or fictitious forces. The stress tensor can be decomposed into a volumetric and a deviatoric term. The volumetric stress is caused by pressure, p, while the deviatoric stress is due to viscosity. Mathematically, σ = −pI + τ , (2.5) where τ is the viscous stress tensor and I is the identity matrix. This viscous stress tensor, in turn, depends on the non-constant fluid viscosity, µ(γ̇), and the rate-of-strain tensor, ε. For incompressible flow, the relation is τ (γ̇) = 2µ(γ̇)ε. (2.6) For infinitesimal deformations, the rate-of-strain tensor can be approximated to only contain first order derivatives of velocity, and is then given by ε = 1 2 ( ∇v + (∇v)T ) . (2.7) From the rate-of-strain tensor, the shear rate may be calculated as γ̇ = √ (2ε : ε), (2.8) where the colon : denotes the double dot product. The divergence of the viscous stress tensor can be calculated from Eq. (2.6) and (2.7). Assuming an incompressible flow and therefore a divergence free velocity field according to Eq. (2.3), we obtain ∇ · τ = µ(γ̇)∇2v. (2.9) Finally, inserting Eq. (2.9) and (2.5) into the original Navier-Stokes momentum equation (2.4), a more explicit form for incompressible flow is obtained: ∂(ρv) ∂t +∇ · (ρv⊗ v) = −∇ · (pI) + µ(γ̇)∇2v + fb. (2.10) The first term on the left-hand side (LHS) in Eq. (2.10) represents the transient behavior of the flow, i.e., how the velocity changes over time. The second term on the LHS is due to convection, which can be thought of as the inclination of a particle to move along with the bulk velocity of its surrounding particles. The first 10 2. Theory term on the right-hand side (RHS) describes the fact that fluids tend to move in the direction of the largest (negative) gradient of pressure. The second term on the RHS is a diffusion term, which is due to viscosity. The third term on the RHS, as already noted, signifies the impact of body forces such as gravity on the flow. It should be noted that the Navier-Stokes momentum equation, Eq. (2.10), consists of vector fields, and therefore makes up a separate equation in each of the x, y and z directions. 2.5 Vessel wall mechanics The deformation of the solid wall can be described by an equation for the conserva- tion of linear momentum, similar to the Navier-Stokes momentum equation [36]: ρs ∂2u ∂t2 −∇ · σs − fs = 0, (2.11) where ρs, u, σs and fs are the solid density, displacement, stress tensor and body force per unit volume, respectively. Apart from Eq. (2.11), the mass contained in the deformed and undeformed configurations must also be equal in order to ensure conservation of mass in the solid throughout the deformation process. In order to solve Eq. (2.11), an explicit relationship between the strain and stress in the solid must be established. Numerous material models of various degrees of complexity exist that provide constitutive laws relating strain and stress. An overview of the material models that will be discussed in this thesis can be seen in Tab. 2.1. Of particular interest is to estimate the product of the variation of strain and the stress, representing the variation of strain energy per unit volume caused by stress, δW : δW = δes : σs, (2.12) where σs is the Cauchy stress tensor and es is the Euler-Almansi strain tensor – the strain and stress tensors associated with the deformed solid configuration. Table 2.1: Overview of presented material models. The material parameters are defined later in this section. The complexity described in the last column is in relation to the other material models, and may be interpreted as the relative level of expected accuracy for large deformations. Of note, in addition to obtaining values of the material parameters, the HGO model also requires specifications of the fiber orientations. Directionality Elasticity Material parameters Relative complexity ILE Isotropic Linear elastic G,K Low Neo-Hookean Isotropic Hyperelastic c10, K Medium Mooney-Rivlin Isotropic Hyperelastic C01, C10, K (for 2-term) High HGO Anisotropic Hyperelastic c, k1, k2 (per layer) Very high 11 2. Theory 2.5.1 Isotropic linear elastic model The simplest material models assume the solid to be an isotropic (i.e., having an identical response to strain in every direction) and linear elastic (i.e., having a linear strain-stress relationship) material [36]. For linear elasticity, the stress tensor σs can be related to the strain tensor εs through a stiffness tensor D: σs = Dεs. (2.13) The infinitesimal strain is defined as εs = 1 2  ∂u ∂X + [ ∂u ∂X ]T  , (2.14) with X being the position of a point in the solid in the undeformed configuration. The deformation u maps a point X to its deformed position x, such that x = X+u. Observe that in deformation theory the actual strain is used, whereas fluid deforma- tion is analyzed using the rate-of-strain, compare Eq. (2.14) and (2.7). Of note is also that for infinitesimal strain, the deformed and undeformed configurations and formulations for the strain and stress are identical. Due to the isotropic assumption, the stiffness tensor written in Voigt notation (omitting symmetrical terms) is given by D =  c11 c12 c12 0 0 0 c11 c12 0 0 0 c11 0 0 0 G 0 0 G 0 G  , (2.15) where the coefficients expressed in terms of the material stiffness quantities shear modulus G and bulk modulus K are c11 = 4 3G+K c12 = −2 3G+K. (2.16) The shear modulus of a material is defined as the ratio between shear stress and shear strain, whereas the bulk modulus quantifies the amount of resistance to compression by pressure. The shear and bulk modulus may also be expressed in terms of the Poisson ratio, ν, and Young’s modulus, E, according to G = E 2(1 + ν) K = E 3(1− 2ν) . (2.17) The Poisson ratio is a measure of how much transverse strain a material experiences in relation to the axial strain. It is typically a value between 0 and 0.5, where the former indicates that axial strains do not result in any transverse strain, and the 12 2. Theory latter that the material is fully incompressible and that the volume remains constant under deformation. Young’s modulus is the ratio of tensile (i.e., uniaxial stretching) stress to tensile strain. For completeness, the stress and strain tensors in Voigt notation are written σs =  σ11 σ22 σ33 σ12 σ23 σ13  ; εs =  ε11 ε22 ε33 2ε12 2ε23 2ε13  , (2.18) where the 2εij, i 6= j elements are called engineering shear strains. The factor 2 is necessary as the total and not the averaged shear strains must be used in Voigt notation in order to obtain the same results as with tensor notation. 2.5.2 Isotropic non-linear hyperelastic models To take non-linear stress-strain behavior into account, hyperelastic models can be used [36]. In these models, a so-called strain energy density function Ψ is determined. For hyperelastic material models, the strain-stress relationship is normally analyzed in the undeformed rather than in the deformed configuration. Eq. (2.12) in the undeformed configuration is written as δW = δE : S, (2.19) where E is the so-called Green-Lagrange strain tensor and S is the so-called 2nd Piola-Kirchhoff stress tensor. The Green-Lagrange strain tensor is given by E = 1 2 ( FTF− I ) . (2.20) The 2nd Piola-Kirchhoff stress tensor is obtained through differentiating the strain- energy density function: S = 2∂Ψ ∂C , (2.21) where C is the right Cauchy strain tensor defined as C = FT F, (2.22) and where F, in turn, is the deformation gradient given by F = ∂x ∂X = I + ∂u ∂X . (2.23) Now, what differs the various hyperelastic models from each other is the strain- energy density function employed. The simplest hyperelastic model is the Neo- Hookean model, fundamentally an extension of Hooke’s law to hyperelastic solids. The Neo-Hookean strain-energy density function is 13 2. Theory ΨNeo−Hookean = c10 ( I1 − 3 ) + K 2 (J − 1)2, (2.24) where c10 = G/2 with G being the initial shear modulus, I1 = tr(C) = tr(J−2/3C) is the first invariant of the deviatoric (i.e., isochoric) right Cauchy stress tensor C, K is the initial bulk modulus and J = det(F) = λ1λ2λ3 is the volume ratio with λi denoting the principal stretches. For incompressible solids, the volume is constant and therefore J = 1. As the shear and bulk moduli for a specific material typically are known at least to an approximate degree (i.e., can be computed from the Poisson’s ratio and Young’s modulus using Eq. (2.17)), the Neo-Hookean model can be used without requiring extensive strain-stress experiments. To obtain more material-specific results, a Mooney-Rivlin strain-energy density function, for example, can be used, given as ΨMooney−Rivlin = N∑ i,j=0 [ Cij ( I1 − 3 )i ( I2 − 3 )j ] + K 2 (J − 1)2, (2.25) where N determines the number of terms, Cij are material-specific parameters, and I2 = 1 2 ( (tr(C))2 − tr ( C2)) is the second invariant of the deviatoric right Cauchy stress tensor. For N = 1, a 2-term Mooney-Rivlin model is obtained. Likewise, N = 2 and N = 3 give 5-term and 9-term Mooney-Rivlin models, respectively. Relating the 2-term material-specific parameters to the shear modulus, we get C01 + C10 = µ 2 . (2.26) 2.5.3 Anisotropic non-linear hyperelastic model A material model developed specifically for arterial walls is the Holzapfel-Gasser- Ogden (HGO) model [29]. The HGO model takes the multilayered nature of the arterial walls into account. Normally, the mechanical response of the innermost layer, the intima, is considered negligible in comparison to the media and adventi- tia layers in healthy patients. However, in atherosclerotic patients, the mechanical response of the intima layer may be highly impactful due to the build-up of plaques on the endothelial walls. The media is often considered to be the most important layer for modeling the mechanical behavior of arteries. The adventitia is mechani- cally active mostly at high pressures, as then the collagen helices are stretched to a straightened position, markedly increasing the stiffness of the artery. In the HGO model, the deviatoric part of the strain-energy density function is decomposed into an isotropic part, Ψiso, and an anisotropic part, Ψaniso, according to Ψ ( C, a01, a02 ) = Ψiso (C) + Ψaniso ( C, a01, a02 ) , (2.27) where the helical collagen structures are modeling with two different collagen types (collagen-I and collagen-III) with reference vectors a01 and a02, respectively. The isotropic part is independent of these reference vectors, as it is assumed that collagen is only active at high pressure and produces an anisotropic response. The strain- energy density function in Eq. (2.27) may (approximately) be expressed in terms of tensor invariants as 14 2. Theory Ψ ( C,A1,A2 ) = Ψiso ( I1 ) + Ψaniso ( I4, I6 ) , (2.28) where Ai = a0i ⊗ a0i, i = 1, 2 (2.29) describe the helical structures of the collagen, and the 4th and 6th invariants I4 and I6 are defined as I4 ( C, a01 ) = C : A1; I6 ( C, a02 ) = C : A2. (2.30) I4 and I6 are the squares of the stretches in the a01 and a02 directions, respectively. The isotropic part is modeled with a Neo-Hookean potential, with the arterial wall assumed to be incompressible (compare Eq. (2.24) with J = 1) Ψiso ( I1 ) = c 2 ( I1 − 3 ) , (2.31) where c is a stress-like material parameter equal to the initial shear modulus in the pure Neo-Hookean strain-energy density function. The anisotropic strain-energy density function is assumed to be given by Ψaniso ( I4, I6 ) = k1 2k2 ∑ i=4,6 { e [ k2(Ii−1)2 ] − 1 } , (2.32) where k1 and k2 are stress-like and dimensionless parameters, respectively. The exponent in Eq. (2.32) describes the collagenous behavior of increased stiffening at high pressures. The strain-energy density function suggested in Eqs. (2.28), (2.31) and (2.32) can be used for each layer of the artery, only the values of the material parameters must be adjusted. 2.6 Fluid-structure interaction The boundary conditions at the interface, Γ, between the solid and the fluid are relatively simple. Assuming the widely accepted no-slip condition, the displacement of the solid, ds Γ, must simply equal the displacement of the fluid, df Γ, at the interface. Mathematically df Γ = ds Γ. (2.33) This also holds true for the rate of change of the displacement. That is, a relation between the fluid velocity, vf Γ, and the rate of change of displacement of the solid at the interface can be obtained through n̂ · vf Γ = n̂ · ∂ds Γ ∂t , (2.34) where n̂ is the interface normal vector. Moreover, there must be a balance of trac- tions on the interface, which can be formulated as 15 2. Theory n̂ · σs Γ = n̂ · σf Γ, (2.35) where σs Γ is the solid and σf Γ is the fluid stress at the fluid-structure interface, respectively. Eq. (2.35) describes the dynamics of the fluid-structure interface, while Eqs. (2.33) and (2.34) describe its kinematics. 2.7 Hemodynamic parameters To characterize and quantify the blood flow, numerous hemodynamic parameters have been defined. Apart from flow velocity, some of the most common ones include the wall shear stress (WSS), the time-averaged wall shear stress (TAWSS) and the oscillatory shear index (OSI) [12]. The WSS, τw, is defined as τw = µ(γ̇)∂v ∂n̂ , (2.36) where n̂ is the normal vector of the endothelial surface. The WSS is the force per unit area exerted by the fluid flow on the endothelial wall. The WSS is evaluated at specific instances, and as an alternative the TAWSS can be used which averages the WSS over a cardiac cycle: TAWSS = 1 T ˆ T 0 |τw|dt, (2.37) where T is the period of the cardiac cycle. The OSI is a non-dimensional hemody- namic parameter that describes the amount of flow reversal. The value ranges from 0.0, signifying uni-directional flow, and 0.5, where the flow is fully oscillatory, i.e., experiences a full reversal of flow direction during the duration of a cardiac cycle. The definition for the OSI is given as OSI = 1 2 1− ∣∣∣´ T 0 τwdt ∣∣∣´ T 0 |τw|dt  . (2.38) 2.8 Fourier decomposition Fourier series can be used to approximate periodic functions [37]. In this thesis, pressure and mass flow rate waveforms imposed as boundary conditions represent the functions we wish to describe as Fourier series. A function f(x) is decomposed into a series of orthogonal harmonic sinusoidal components with corresponding weights. Mathematically, a Fourier series FN(x) of order N is given by FN(x) = a0 2 + N∑ n=1 ( an cos (2π P nx ) + bn sin (2π P nx )) , (2.39) where P is the period of the function and a0, an and bn are weights which can be calculated from 16 2. Theory a0 = 2 P ˆ P f(x)dx, (2.40) an = 2 P ˆ P f(x) · cos (2π P nx ) dx, (2.41) bn = 2 P ˆ P f(x) · sin (2π P nx ) dx. (2.42) 2.9 Richardson’s extrapolation Richardson’s extrapolation is a method for estimating the asymptotic value of a converging series from only a few data points [38]. In the context of simulations performed with meshes (grids), Richardson’s extrapolation can be used to estimate the value, f0, of a variable of interest for a mesh with infinitesimally small cells (i.e., with infinite grid resolution). Given the values of the variable of interest calculated using a coarse (f3), medium (f2) and fine (f1) mesh, the true value is assumed to be given by f0 = f1 + f1 − f2 rp − 1 , (2.43) where p represents the order of convergence, and r = b2 b1 = b3 b2 is a constant grid refinement index, with bi being representative grid cell sizes of mesh i. The order of convergence can be obtained through p = ln ( f3−f2 f2−f1 ) ln r . (2.44) The grid convergence index (GCI) is a conservative estimate of the maximal error of the values one could expect to obtain from a specific mesh relative to the asymptotic value [38]. The GCI for grid i is given by GCIi = Fs · ( fi−fi+1 fi ) rp − 1 , i ∈ {1, 2}, (2.45) where Fs is a factor of safety commonly set to 1.25 for three-grid convergence anal- yses [38]. 2.10 Analytical displacements in cylinders The Law of Laplace is commonly used in cardiovascular physiology, as it provides a simple, although highly idealized, relation between the average wall stress (σw) and pressure (P ) in a cylinder, given its radius (r) and wall thickness (tc) [39]: σw = P · r tc . (2.46) 17 2. Theory Assuming an isotropic linear elastic material model for the cylinder and considering displacements in the radial direction (∆r) only, we get the following expression ∆r = P · r2 E · tc , (2.47) where E is the elastic modulus. However, in the Law of Laplace in (2.46), a thin- walled cylinder is assumed (tc � rc). For thick-walled cylinders with Poisson’s ratio ν, the following expression can be derived [40]: ∆r = Pr2 i r E (r2 o − r2 i ) [ (1− ν) + (1 + ν) r2 0 r2 ] , (2.48) where ri and ro are the radii at the inner and outer part of the wall, respectively, and r ∈ [ri, ro] is the radius at the point for which the displacement is calculated. 18 3 Methods This chapter deals with how to set up and perform the desired simulations in STAR- CCM+. Given the theoretical frameworks presented in Chapter 2, we first investi- gate how the fluid and solid domains can be solved separately in Sections 3.1 & 3.2, respectively. Then, measures necessary to couple the two domains are described in Section 3.3. Sections 3.4 & 3.5 cover the preprocessing and meshing steps performed, respectively. Following this, Section 3.6 focuses on boundary condition acquisition and presentation. In Section 3.7, the types of simulations and evaluations conducted are described. Finally, model and parameter specifications are given in Section 3.8. 3.1 Solution strategy for the fluid flow A method for solving the Navier-Stokes equations and other partial differential equa- tions (PDEs) is the Finite Volume Method (FVM) [36]. In FVM, the object under investigation is divided into a grid composed of small cells (finite volumes). For each cell, the Navier-Stokes momentum equation in Eq. (2.10) is first integrated over the entire cell volume. Then, Gauß’s divergence theorem is employed, which relates the divergence of a vector field F in a volume V to the flow across its closed boundary A, ˚ V (∇ · F)dV = ‹ A F · dA, (3.1) where dA = n̂dA, with n̂ being the outward facing normal vectors of the boundary. Gauß’s divergence theorem is used to convert all volume integrals in the Navier- Stokes equations containing divergence operators to surface integrals. I.e., combining Eq. (3.1) and (2.10) integrated over a cell, we obtain d dt ˚ V ρvdV + ‹ A ρv⊗(vr−vg)·dA = − ‹ A pI·dA+ ‹ A µ(γ̇)∇v ·dA+ ˚ V fbdV, (3.2) which is referred to as a weak form of the Navier-Stokes momentum equation (com- pared to the strong form in Eq. (2.10)). The convection term has been adjusted in Eq. (3.2) to account for the movement of the fluid cells themselves due to the fluid-structure interaction. This is accomplished through the so-called arbitrary Lagrangian-Eulerian (ALE) formulation used in Eq. (3.2), where vr is the fluid velocity relative to the moving grid, and vg is the grid velocity. Now, instead of computing exact and complicated integrals, approximations can be made. Since 19 3. Methods each cell volume is relatively small, the volume integrals can be approximated with the values of the integrand at the cell center multiplied by the volume of that cell. In a similar fashion, the surface integrals can be calculated taking integrand values at each cell face center multiplied by the area of that face, and then summing over all faces of a cell. This results in the following discretized equation: ∂ ∂t (ρvV )c+ ∑ f (ρv⊗(vr−vg)·A)f = − ∑ f (pI·A)f + ∑ f [µ(γ̇)∇v·A]f +(fbV )c, (3.3) where the subscripts f and c represent the cell faces and cell center, respectively. In this work, body terms such as gravity are assumed to be negligible in comparison to the acceleration of blood flow in the carotid artery, and the last term in Eq. (3.3) is therefore assumed to be zero. Cell face values are computed using a second-order upwind differencing scheme. The upwind scheme assigns values to a face between two cells by taking the direction of fluid flow into account, and uses values stored upstream (i.e., towards the inlet) in the carotid artery. Gradients were approximated by a hybrid Gauss-least squares method in combination with a Venkatakrishnan gradient limiter, see [36] and [41] for more details. To obtain a fully discretized form of the Navier-Stokes momentum equation, the transient term in Eq. (3.3) must also be discretized. A backwards (implicit) Euler scheme was used, providing a first-order temporal discretization of the form d dt (ρvV )n+1 c = (ρvV )n+1 c − (ρvV )n c ∆t , (3.4) where n + 1 is the current time step, n is the previous time step, and ∆t is the time step size. Having discretized the Navier-Stokes equations, we now have four equations per cell (i.e., the continuity equation (2.3) and the momentum equation Eq. (3.3) in the x, y and z direction), and four unknowns (pressure and the velocity in each direction). However, there is no explicit equation for the pressure. Moreover, the continuity equation is a restriction to the momentum equation in the sense that the number of allowed velocity values are reduced. Therefore, an algorithm called SIMPLE (Semi-Implicit Method for Pressure Linked Equations), in STAR-CCM+ referred to as the Implicit Unsteady solver, is used [36]. First, with the aim of obtaining an equation for the pressure, the Navier-Stokes momentum equation, Eq. (3.3), may be rewritten in the form MV = −∇P (3.5) whereM is a matrix with known element values calculated from Eq. (3.3), V is a vector containing the unknown velocities of all cells in all three directions, and P are the unknown pressures in the cells. We then decompose the M matrix into a diagonal part A and a velocity dependent part H(V) such that: MV = AV−H(V). (3.6) Inserting Eq. (3.6) into Eq. (3.5), multiplying both sides from the left with the inverse of the A matrix, and finally taking the divergence of both sides and using 20 3. Methods that the divergence of the velocity is zero, results in the following Poisson equation for the pressure: ∇ · ( A−1∇P ) = ∇ · ( A−1H(V) ) . (3.7) Since A is diagonal, its inverse is simply the reciprocals of its diagonal elements. In the first step of the SIMPLE algorithm, the (so-called uncorrected) velocities are calculated from the momentum equation (3.5) using an initial guess for the pressures. Since these velocities likely do not fulfill the continuity equation, the pressures are calculated from equation (3.7) using the uncorrected velocities. The velocities are then corrected using the new pressures values. However, since the velocities have been updated, they now do not necessarily fulfill the momentum equation (3.5). Therefore, the values of theM matrix are updated to suit the corrected velocities, and the process is repeated from the first step using the updated pressure field as the initial guess in the momentum equation. Iteratively, pressure and velocity values that satisfy both the momentum and the continuity equations for each time step will be obtained. A similar algorithm to SIMPLE is the PISO (Pressure-Implicit with Splitting of Operators) algorithm. The difference between the two algorithms is that PISO performs the pressure correction two times per iteration instead of one, speeding up convergence. However, this can lead to instability for large so-called Courant numbers (i.e., for large time step sizes and fluid velocities relative to the size of the cells) [36]. The SIMPLE algorithm is therefore chosen as a means to reduce the risk of diverging solutions. Under-relaxation factors are used for both the velocity and pressure to avoid exceedingly large changes between consecutive iterations that might lead to instability. The pressure, for instance, is updated as follows: Pn+1 = (1− ωp)Pn + ωpP̃n+1, (3.8) where ωp is the under-relaxation factor for pressure, Pn+1 is the pressure used as input for the current iteration, Pn is the input pressure from the previous iteration, and P̃n+1 is the pressure obtained from solving the Poisson equation (3.7) in the previous iteration, and which would be used as the input pressure to the current iteration in case of no under-relaxation (ωp = 1). 3.2 Solution strategy for the solid displacements While the FVM may be considered commonly used for fluid flow simulations, a similar method called the Finite Element Method (FEM) is sometimes preferred for solid deformations [36]. FEM is similar to the FVM in that it first discretizes the space into a finite number of small parts, called elements. Another resemblance between the two methods is that the PDEs in question are integrated over the entire volume of the element, resulting in weak forms of these PDEs. A key difference, however, is that the FEM multiplies the PDEs by a test function and thereby utilizes a principle of virtual work for solving the PDEs. The governing equation for the displacement of a solid, Eq. (2.11), written in the weak (Lagrangian) formulation used in the FEM is then 21 3. Methods ˚ V δu · ( ρs ∂2u ∂t2 −∇ · σs − fs ) dV = 0, (3.9) where the test function δu is a virtual displacement satisfying the displacement boundary conditions (assumed to be 0 at the surface), and which should minimize the value of the LHS. By manipulating Eq. (3.9), the following equation can be obtained: ˚ V δu · ρs ∂2u ∂t2 dV − ˚ V δu · fsdV + ˚ V δWdV − ‹ A δu · (σs · n̂)dA = 0, (3.10) where δW (obtained through integration by parts) is the strain energy variation per unit volume caused by stress, expressed either in the deformed, Eq. (2.12), or undeformed, Eq. (2.19), formulation. Gravity and other body forces are ignored also in the solid domain, reducing the second term in Eq. (3.10) to zero. Values for the unknown variables, i.e., the displacement u, are calculated and stored in nodes (i.e., vertices) of the elements. To approximate the unknowns at each point in space inside an element, and not just at the nodes, basis functions are used. Basis functions are polynomials that interpolate the values, typically linearly or quadratically, between nodes. All integrals in Eq. (3.10) are discretized and evaluated using Gaussian quadrature, which approximates an integral by evaluating the integrand at specific integration points [42]. The Gaussian quadrature of a function f(X) integrated over a volume or surface D, reads ˆ D f (X) dD = N∑ i=1 wif (X (ξi)) det Ji (3.11) where N is the number of integration points, wi are weights, Ji is the Jacobian associated with the change of integration variables, and ξi is a vector containing the integration point locations in two dimensions for a surface integral, and in three di- mensions for a volume integral. The integration point locations and weights depend on the number of integration points used, but not on the integrand. For example, for N = 2, the integration points locations are ± 1√ 3 in each dimensional direction, and the weights are all equal to 1. Discretized versions of Eq. (3.10) using Eq. (3.11) can then be combined for all nodes into a system of linear equations, which can be converted into a matrix equation of the form AU = r, (3.12) where A is a stiffness matrix containing known coefficients, U is a vector containing the unknown displacements for all nodes in each of the three-dimensional directions, and r is a residual (force) vector. The stiffness matrix is symmetric and its number of rows is three times the total number of nodes in the system. The basis functions that approximate the displacements in general only have local support, meaning that they are non-zero only close to a specific node, which results in a sparse stiffness matrix. Eq. (3.12) can then be solved for U in each iteration by using MUMPS, a sparse direct solver [43]. The stiffness matrix is updated in each inner iteration 22 3. Methods using a full-Newton approach, i.e., by iteratively estimating the zeros of functions by performing a Taylor expansion and using the inverse of the Jacobian matrix to update the estimate. 3.3 Fluid-structure coupling Having solution strategies at hand for the fluid and solid domain problems separately, these now have to be coupled to allow for a complete simulation which takes into account how the fluid flow affects the solid deformations, and vice versa. As the fluid and solid problems are solved separately with distinct solvers, a so-called partitioned approach is used for the FSI. The fluid and solid domain are assumed to be strongly coupled, and therefore a two-way information exchange across the fluid-structure interface is needed. The traction, i.e. the pressure and wall shear stress, from the fluid domain is mapped to the solid domain across the interface. This corresponds to the dynamic boundary condition, i.e., Eq. (2.35), and leads to deformations in the solid. These solid deformations, in turn, are mapped back to the fluid domain across the interface, corresponding to the kinematic boundary conditions, Eqs. (2.33) and (2.34). Since the solid deformation only gives information on how to displace the fluid cell nodes at the interface, the displacements must then be propagated throughout the entire fluid domain. The calculation of each node displacement is done using B-splines [44]. This process involves the use of control points at which the displacements are known, and a construction of an interpolation field which is used to specify the displacements of the nodes of the mesh. The interpolation field is made up of a linear combination of the values of the control points and specific basis functions at those points. Mathematically, the interpolation field S as a function of the [0, 1]-normalized u, v and w coordinates to the nodes of the mesh is given by S(u, v, w) = l∑ i=0 m∑ j=0 n∑ k=0 Pi,j,kNi,p(u)Nj,p(v)Nk,p(w), (3.13) where Pi,j,k are the (l+1)(m+1)(n+1) number of three-dimensional control points and Ni,p(u) is the i:th B-spline basis function of order p in the u-direction, and sim- ilar for the basis functions in the v and w direction. The B-splines are polynomials recursively defined through the so-called Cox-de Boor formula, see e.g., [45]. As with the basis functions used in FEM, the B-spline basis functions only have local support. This deformation, or more specifically morphing, of the fluid domain is normally taken into consideration in the flow formulations by rewriting the Navier- Stokes momentum equation in ALE form, Eq. (3.2), as described in that section. However, to avoid instability due to undesirably high velocity, mass flux and pres- sure adjustments for the flow, grid flux terms were ignored in the corrections in the SIMPLE algorithm. (See the doctoral thesis from Prather [46] for detailed math- ematical reasoning to substantiate this claim.) To ensure mass conservation while the fluid and solid meshes deform, a Space Conservation Law (SCL) is employed [36]. The SCL is much the same as the compressible continuity equation, Eq. (2.2), stated in integral form for a cell with boundaries moving with velocity vb, 23 3. Methods ∂ ∂t ˚ V dV − ‹ A vb · dA = 0. (3.14) Eq. (3.14) is easily discretized by using the same tactics as those in Section 3.1 for discretizing the Navier-Stokes momentum equation. To stabilize the FSI calculations and to speed up convergence, the fluid traction resulting from movements of the interface were predicted in each iteration. Essentially, the displaced volume of fluid per unit area is estimated and treated as an added mass which contributes to a force correction. This stabilization method is known as Boundary Interface Added Mass in STAR-CCM+ [36]. 3.4 Preprocessing of the carotid artery geometry A patient-specific carotid artery geometry of a healthy subject was used, see Fig. 3.1. The geometry file in STL format – a popular format for storing information about three-dimensional objects – had been obtained from segmentation from MRI (mag- netic resonance imaging) images of the patient and was provided by the ICM. It was assumed that the geometry depicted the inner wall of the carotid artery at diastolic pressure, i.e., during the period of the cardiac cycle when the ventricles of the heart are refilling with blood and preparing to eject it out to the rest of the cardiovascular system. The patient-specific diameters of the CCA, ICA and ECA at the ends were approximately 6.2 mm, 4.2 mm, and 3.0 mm, respectively. The carotid sinus is clearly visible in the frontal view, Fig 3.1a, as the dilated area of the ICA at the bifurcation. The preprocessing procedure is outlined in Fig. 3.2. First, the geometry was im- ported into STAR-CCM+ as an STL file, Fig. 3.2a. Second, the (preliminary) inlet (CCA) and outlets (ICA and ECA) of the carotid artery were constructed using an operation in STAR-CCM+ called Fill Holes, Fig. 3.2b. Third, since the geometry solely consisted of the (infinitesimally thin) inner arterial wall, also the outer wall as well as the space in-between the two had to be created. To this end, the inner wall was extruded (using a Surface Extruder in STAR-CCM+) with a distance of 0.6 mm in the direction of the surface normal. The arterial wall thickness can vary greatly from patient to patient, and the value of 0.6 mm was chosen as the wall thickness is thought to generally lie between 0.5 mm to 1 mm and be approximately 8 % to 10 % of the CCA diameter (in this case equal to 6.2 mm) [47] [48]. Fourth, due to poor quality in the surface mesh of the extruded surface (resulting from complex geometry and the imported STL file being imperfect) a further step using a Surface Wrapper and Surface Remesher was used to ensure closed and manifold surfaces of high enough quality to be used for meshing the volume of the arterial wall, Fig. 3.2d. Fifth, the inlets and outlets were extended with a distance heuristically set to 3 mm, as a measure to reduce non-physical behavior close to the inlet and outlets due to the imposing of boundary conditions, Fig. 3.2e. Last, the fluid domain (includ- ing the extensions) was wrapped to improve the quality of its surface mesh, Fig. 3.2f. 24 3. Methods (a) (b) (c) Figure 3.1: Carotid artery geometry. (a) Rotated frontal view. The CCA, ICA and ECA correspond to the parts of the artery to the left, upper-right, and lower-right, respectively. (b) Rotated medial view. (c) Superior view. (a) Import geometry (b) Enclose fluid (c) Generate wall (d) Repair wall (e) Extend ends (f) Repair fluid Figure 3.2: Preprocessing pipeline. (a) Import the carotid artery geometry. (b) Enclose the fluid surface using the Fill Holes operation in STAR-CCM+. (c) Gener- ate the arterial wall by using a Surface Extruder with a specified offset, in this paper equal to 0.6 mm. (d) Wrap the arterial wall surface to improve the mesh quality using a Surface Wrapper. (e) Extend the CCA inlet and outlets of the ICA and ECA with a Surface Extruder (here, a 3 mm offset was used). (f) Wrap the fluid domain with a Surface Wrapper. 25 3. Methods 3.5 Meshing of the domains Automated volume meshes were generated for the fluid and solid domains. For the fluid domain, polyhedral cells were used to form the mesh. In order to accurately capture the gradients of the fluid velocity and pressure close to the fluid-structure interface, prism layers were used. A mesh independence analysis was performed to ensure that the number of cells used was enough to capture the behavior of the fluid to a sufficient degree. The area of low TAWSS (< 0.4 Pa) and high OSI (> 0.2) were calculated for five different meshes, varying from very coarse (66 610 cells) to very fine (1 282 825 cells), see Fig. 4.2. Richardson’s extrapolation (see Section 2.9) was then used to estimate the values at an infinite amount of cells, i.e., the true values, as well as the GCI for two short-listed meshes. A nearly constant grid refinement ratio, r ≈ 1.5, was employed between each mesh refinement from the first (very coarse) to the fourth (second-to-last in fineness) mesh, enabling Richardson’s extrapolation using three of these values according to Eqs. (2.43)–(2.45). The representative grid cell sizes bi of mesh i were calculated as [49]: bi = ( 1 Ni · V ) 1 3 , (3.15) where Ni is the total number of grid cells in mesh i, and V is the total volume of the carotid artery. The final mesh chosen for the fluid domain based on the grid independence analysis had a target base size of 0.125 mm, and the minimum allowed size of a cell side was set to 20 % of this base size. Simulations were also performed where the number of prism layers ranged from three to five, and their total thick- ness from 0.1 mm to 0.2 mm, to investigate whether any noticeable differences in the TAWSS and OSI risk zones were obtained. For the final mesh, five prism layers of a total thickness of 0.2 mm were used. The final number of cells for the entire fluid domain were 903 232, of which 772 403 made up the inside of the carotid artery geometry, and 130 829 were part of the extensions at the inlet and outlets. The final mesh of the fluid domain can be seen in Fig. 3.3a. For the solid domain, the mesh was composed of tetrahedral cells. The base size and minimum element size were heuristically selected as to result in at least a three layer thick mesh across the entire solid domain (assumed to be enough to enable realistic calculations of the wall deformation, see e.g., [11]), but aiming to minimize the number of four-layered areas (to reduce computational demand in the simulations). From generating meshes while varying these parameters, a final target base size of 0.24 mm and a minimum cell size of 50 % were chosen for the solid mesh. The final number of elements for the solid domain were 486 558. The final mesh of the solid domain can be seen in Fig. 3.3b. 26 3. Methods (a) Fluid mesh (b) Solid mesh Figure 3.3: Superior view of the fluid and solid meshes. (a) Fluid mesh consisting of polyhedral cells and five prism layers. The 3 mm extensions at the inlet and outlets are not shown. (b) Solid mesh consisting of tetrahedral elements. 3.6 Boundary conditions The boundary conditions specify the values of e.g., pressures and velocities (flow rates) at the arterial boundaries, namely, the inlet of the CCA and the outlets of the ICA and ECA. Patient-specific carotid artery pressures can be obtained non- invasively by means of applanation tonometry, or invasively using an angiocatheter [12]. Velocities may be estimated using phase contrast magnetic resonance imaging (PCMRI) or Doppler Ultrasound [12]. Flow rates can be calculated from velocity profiles by integrating over the inside of the artery (i.e., the arterial lumen). In this study, the patient had undergone no such procedures, and therefore patient-specific boundary conditions were not available. In order to get as realistic boundary con- ditions as possible, albeit not pertaining to the patient under analysis, pressure and velocity waveforms were taken from another patient-specific FSI study on carotid artery blood flow. The ICA and ECA flow rate curves and the CCA pressure curve in [50] were digitized into comma-separated values files using the software Graph Grabber v.2.0.2 [51], and were subsequently approximated by Fourier series accord- ing to Eqs. (2.39)–(2.42) in MATLAB R2021a. The degrees of the Fourier series were selected as three for the pressure waveform and four for the two flow rate wave- forms, as a trade-off between model complexity and accurately capturing the general trends of the waveforms. A cardiac cycle period of 1 s, i.e., 60 beats per minute, was assumed. Therefore, to avoid discontinuities between heart beats, the fundamental angular frequency was set to 2π for all Fourier series. The flow rates were converted from volumetric flow rates (ml/s) to mass flow rates (kg/s). In STAR-CCM+, mass flow rate outlets are not explicitly available for use as boundary conditions. For that reason, the mass flow rates were multiplied by minus one and specified as mass flow inlets at the ECA and ICA outlets. The pressure waveform was subtracted by its initial value (80.21 mmHg) at all points, as it was assumed that approximately diastolic pressure was already present in the acquired carotid artery geometry. See Fig. 3.4 and Tab. 3.1 for visualizations and Fourier coefficients of the final boundary conditions used. The solid boundaries surrounding the fluid inlet and outlets were fixed and therefore constrained not to move. 27 3. Methods Table 3.1: Fourier coefficients of the boundary conditions. The respective Fourier series are defined in Eq. (2.39). Fourier coefficients CCA pressure (mmHg) ECA flow rate (ml/s) ICA flow rate (ml/s) a0 23.88 3.0 11.116 a1 -3.869 -0.04154 -0.4895 a2 -6.037 -0.4128 -0.8089 a3 -2.028 -0.4288 -0.4301 a4 N/A 0.09458 0.3791 b1 11.04 0.7380 1.503 b2 1.573 0.2995 -0.01052 b3 -1.764 -0.2191 -0.6115 b4 N/A -0.3069 -0.3223 (a) Pressure inlet (b) Mass flow rate outlets Figure 3.4: Boundary conditions imposed at the inlet and outlets in each cardiac cycle. a) Pressure waveform at the CCA inlet. The pressure is expressed relative to 80.21 mmHg. b) Mass flow rates at the ICA and ECA outlets. Note that the volumetric flow rates in Tab. 3.1 have been converted to mass flow rates here. 3.7 Simulation operations and evaluation Simulations were run in Star-CCM+ for a pure CFD calculation, as well as FSI with the isotropic linear elastic model, and the two isotropic non-linear hyperelastic Neo-Hookean and (2-term) Mooney-Rivlin models. The HGO model was not used, as it is not available in STAR-CCM+. It is nevertheless included in the theory sec- tion, as part of providing a comprehensive workflow. Each FSI simulation was first run for two time steps, using only CFD to solve the fluid domain. The velocity and pressure fields calculated with CFD were then used as inputs, i.e., initial guesses, for that FSI simulation. The reason for this was that performing FSI simulations without initial velocity and pressure fields from a CFD simulation at times resulted in elevated pressures in the order of ±200 mmHg relative to diastolic pressure in 28 3. Methods the initial inner iterations. These elevated pressures, in turn, would cause large deformations of the fluid and solid meshes, resulting in non-positive cell or element volumes and thus an aborted simulation. The evaluated quantities were WSS, TAWSS, OSI, velocities, and wall displace- ments. The area of TAWSS < 0.4 Pa and OSI > 0.2 were compared between the FSI models and CFD in the peak systolic and end diastolic (defined here as the end of the cardiac cycle) geometric configurations. Three points were chosen to probe the wall displacements, visualized in Fig. 3.5. S1 is close to the carotid bifurca- tion point, where large deformations can be expected. S2 and S3 are part of the CCA wall. The analytical expression for thick-walled cylinders in Eq. (2.48) was employed as a comparison. The analytical displacements were calculated at r = ro, using approximate radii at the three points, where the hydraulic radius for rectan- gular ducts was used at S1 (amounting to ≈ 5.4 mm) due to the rather non-circular cross-section. Velocities were also probed in three different points – in the CCA, ICA, and ECA, see Fig. 3.5 – as well as averaged over the entire fluid domain. All aforementioned quantities were also evaluated visually at peak systole. Figure 3.5: Points selected for evaluative purposes. The points with prefix S belong to the solid wall and are used to evaluate the magnitude of wall deformation, while for the rest, fluid velocities are investigated. 3.8 Models and parameter values Both the arterial wall and the fluid were assumed to be nearly incompressible [29] [11]. The Carreau-Yasuda model was used to model the non-Newtonian rheology of the fluid. The rheological parameter values were taken from [52] and converted to the formulation used in Eq. (2.1), giving: µ∞ = 0.0035 Pa · s, µ0 = 0.016 Pa · s, λ = 8.2 s, a = 0.64 and n = 0.22. An under-relaxation factor (of the form in Eq. (3.8)) of 0.7 was used for the viscosity in order to improve convergence. For velocities and pressures, under-relaxation factors of 0.8 and 0.2 were used, respectively. The flow was assumed to be laminar, as has been proved to be the case in healthy carotid arteries due to low Reynolds numbers [12]. Densities of 1060 kg/m3 and 1300 kg/m3 29 3. Methods were used for the blood and vessel wall, respectively [46] [53]. For the isotropic linear elastic model described in Section 2.5.1, a Young’s modulus E = 0.9 MPa and a Poisson ratio ν = 0.49 were assumed [54]. Identical values were used for the Neo-Hookean model, Eq. (2.24), but converted to bulk and shear mod- uli according to Eq. (2.17). For the 2-term Mooney-Rivlin model, the ratio of the C10 and C01 parameters in Eq. (2.25) were computed from [55], and the parameter values were scaled to satisfy Eq. (2.26) in order to ensure consistency with the other models used. This resulted in C01 = 57 kPa and C10 = 94 kPa. A time step size of 5 ms was used as a trade-off between accurate temporal res- olution and computational time needed for running the simulations, as well as to avoid numerical instability from too small time step sizes. 50 inner iterations (i.e., iterations within time steps) were used to reach convergent solutions in each time step. With an assumed heart rate of 1 Hz = 60 beats per minute, this corresponded to 200 time steps and 10 000 inner iterations per cardiac cycle. In order to mini- mize the influence of initial transient behavior on the solutions, two whole cardiac cycles were simulated. All values presented in the results are therefore taken from the second cardiac cycle. The simulations were run using 60 cores in parallel of an AMD EPYC 7702 64-core processor with 256 gigabytes of RAM. 30 4 Results The results of this study consist of 1) A workflow developed for guiding patient- specific FSI simulations of the carotid artery, 2) Mesh sensitivity information and GCIs from Richardson’s extrapolation, 3) Visualizations and quantifications of hemo- dynamic parameters using CFD and FSI with the isotropic linear elastic, Neo- Hookean, and Mooney-Rivlin models, and 4) Visualizations and quantifications of the carotid wall deformation using the three material models. 4.1 Overall workflow The steps to be taken for solving and setting up FSI analyses of patient-specific carotid arteries, which have been described throughout this thesis, are summarized in a workflow in Fig. 4.1. The procedure starts with the preprocessing of the ar- terial geometry, as described in Section 3.4. Then, the fluid and solid continua have to be set up, e.g., by specifying the fluid rheology (see Section 2.3), material model (Section 2.5), parameter values (Section 3.8), and the assumptions and solu- tion strategies (e.g., incompressible flow and the SIMPLE algorithm, Chapter 3) to be employed. Next, the two domains have to be meshed (Section 3.5), and suitable boundary conditions imposed (Section 3.6). To enable FSI, the motion specifications of the fluid and solid domains then have to be selected, where the fluid domain must be morphed using e.g., B-splines, and an interface between the two domains be es- tablished where fluid traction is mapped to the solid domain and solid displacements are mapped to the fluid domain (Sections 3.1, 3.2 & 3.3). Some hemodynamic pa- rameters such as WSS are already available in STAR-CCM+, whereas others such as OSI can be readily implemented directly from, or by manipulating, the expressions from their mathematical definitions (Section 2.7). Indispensable stabilization tech- niques, which prevent uncontrollable mesh deformations. are the boundary interface added mass method and the omittance of gridflux terms in the transport equations (Section 3.3). A further action to increase stability, albeit only for the very first time steps of a simulation, is to calculate initial velocity and pressure fields with CFD, and provide these as initial guesses for the starting fields in FSI (Section 3.7). Having all these settings in place, a simulation can be started. In case of instabil- ity, usually manifesting itself visually through unreasonably large deformations (or non-positive mesh cell volumes, in which case the simulation will crash) or through diverging residuals, additional stabilization methods may be availed of. Load step- ping, Rayleigh damping, and morphing from zero were not used nor have they been mentioned in this study; for more information on these and the other stabilization 31 4. Results methods listed, the reader is referred to e.g., [36]. In case the simulation appears to have been successful, e.g., through residuals converging to satisfactory low values, and results which seem physically and physiologically realistic, the results may be postprocessed. The postprocessing steps may vary greatly according to which types of results are relevant in the specific use case, but should generally include a mesh sensitivity analysis and an evaluation of the hemodynamic parameters of interest. Figure 4.1: Final workflow for performing FSI analyses of carotid arteries. 4.2 Mesh sensitivity study In Fig. 4.2, the computed values of the variables of interest using five fluid meshes of various degrees of fineness are visualized. As can be seen, the values appear to be close to the asymptotes for the third grid (with cell count = 225 653) and for the finer meshes. In Tab. 4.1, the results of Richardson’s extrapolation using the second to fourth meshes are shown. The representative base cell sizes in Tab. 4.1 have a nearly constant refinement ratio equal to 1.5. Although the third grid seems to be close to the asymptotic values in Fig. 4.2, its corresponding GCIs are 2.0 % and 6.1 %, for the TAWSS and OSI risk areas, respectively. Therefore, the third grid was deemed insufficient in resolution. The GCIs of the grid with 772 403 cells are 0.4 % and 0.9 % for the TAWSS and OSI risk areas, respectively. Considering both the 32 4. Results factor of safety in Eq. (4.1), the sub-percentage GCIs, and the relatively low impact an increased number of fluid domain cells has on the overall computation time in FSI simulations, this mesh was selected. Furthermore, varying the number of prism layers and their total thickness resulted in marginal differences in both the TAWSS and OSI area. Figure 4.2: CFD mesh independence analysis of TAWSS and OSI risk zone areas. The number of cells used for the fluid domain throughout the study were chosen to 772 403 (not including the inlet and outlet extensions). Table 4.1: Richardson’s extrapolation of fluid domain mesh values. Note that the computation time presented in the last column is for one FSI time step, whereas the rest of the columns correspond to CFD values. Cell size (mm) Number of cells (−) TAWSS < 0.4 Pa (−) OSI>0.2 (−) Time for one FSI time step (s) b1 0.134 N1 772403 GCI1 0.4 % GCI1 0.9 % 240 b2 0.202 N2 225653 GCI2 2.0 % GCI2 6.1 % 205 b3 0.304 N3 66610 196 4.3 Hemodynamics The hemodynamic parameters are only visualized for CFD and the Neo-Hookean model, as the differences between the various material models are close to impossible to discern visually. The CFD simulation took approximately 1.5 hours to run, 33 4. Results whereas each FSI simulation had a computational time of approximately 25-30 hours. The TAWSS and OSI are depicted in Fig. 4.3 and Fig. 4.5, respectively, and the deviations in the area of the risk regions relative to CFD are displayed in Tab. 4.2. The WSS at peak systole is also presented visually in Fig. 4.4. The velocity magnitudes in a plane cross-section are visualized at peak systole in Fig. 4.6. The velocities probed throughout the cardiac cycle in the three points selected in Fig. 3.5, as well as the average velocity across the entire fluid domain, are plotted in Fig. 4.7. The corresponding deviations from CFD for the three material models can be seen in Tab. 4.3. (a) CFD (b) FSI Figure 4.3: TAWSS computed from simulations using a) CFD, and b) the Neo- Hookean model. (a) CFD (b) FSI Figure 4.4: WSS at peak systole (t = 0.24 s) computed from simulations using a) CFD, and b) the Neo-Hookean model. (a) CFD (b) FSI Figure 4.5: OSI computed from simulations using a) CFD, and b) the Neo-Hookean model. 34 4. Results (a) CFD (b) FSI Figure 4.6: Velocity magnitudes at peak systole (t = 0.24 s) together with line integral convolutions of the velocity computed from simulations using a) CFD, and b) the Neo-Hookean model. (a) ICA (b) CCA (c) ECA (d) Whole lumen Figure 4.7: Velocity magnitude over time in three discrete points visualized in Fig. 3.5, as well as averaged over the entire fluid domain. The velocities are calculated in the a) ICA point, b) CCA point, c) ECA point, d) whole fluid domain. 35 4. Results Table 4.2: Deviations in the estimated TAWSS and OSI risk areas for the various material models compared to the CFD values as baseline. TAWSS < 0.4 Pa (−) OSI > 0.2 (−) End Diastole Peak Systole End Diastole Peak Systole CFD 0 % 0 % 0 % 0 % ILE +2.14 % +4.92 % +6.90 % +10.07 % Neo-Hookean +2.31 % +5.15 % +7.11 % +10.24 % Mooney-Rivlin +2.12 % +4.70 % +6.64 % +9.59 % Table 4.3: Deviations in the probed velocities at peak systole for the various material models compared to the CFD values as baseline. The ICA, CCA, and ECA points can be seen in Fig. 3.5. In the whole artery column, the mean is compared. ICA CCA ECA Whole artery CFD 0 % 0 % 0 % 0 % ILE -5.57 % -4.03 % -0.38 % -4.21 % Neo-Hookean -5.73 % -4.16 % -0.36 % -4.34 % Mooney-Rivlin -5.39 % -3.88 % -0.45 % -4.05 % 4.4 Wall deformation (a) ILE (b) Neo-Hookean (c) Mooney-Rivlin Figure 4.8: Magnitude of wall deformation at peak systole (t = 0.24 s) computed using the a) isotropic linear elastic model, b) the Neo-Hookean model, and c) the Mooney-Rivlin model. 36 4. Results (a) Point S1 (b) Point S2 (c) Point S3 Figure 4.9: Magnitude of wall deformation in selected points (displayed in Fig. 3.5) from simulations using the three material models, as well as the analytical expression for thick-walled cylinders in Eq. (2.48). Table 4.4: Deviations in the wall deformation magnitudes at peak systole between the various material models, with the Mooney-Rivlin values as reference. The three selected points can be seen in Fig. 3.5. S1 S2 S3 Analytical -52.24 % -39.43 % +43.67 % ILE +4.97 % +5.86 % -8.17 % Neo-Hookean +5.56 % +0.78 % +1.72 % Mooney-Rivlin 0 % 0 % 0 % 37 4. Results 38 5 Discussion This chapter is dedicated to a discussion of the results of the thesis, as well as of the work as a whole. First, the strengths and weaknesses of the developed workflow are discussed in Section 5.1. In Section 5.2, the results in terms of hemodynamics are scrutinized, with a special focus on the differences between CFD and FSI. Thereafter, the various material models are compared to each other in Section 5.3, by examining the obtained wall deformations. Limitations of the study are then discussed in Section 5.4. In the last part of the chapter, Section 5.5, potential ethical and societal considerations associated with this type of study are highlighted. 5.1 Workflow The general workflow depicted in Fig. 4.1, together with the in-depth explanations of the various steps provided in Chapter 2 & 3, enable users to set up and per- form patient-specific FSI simulations of carotid artery blood flow in STAR-CCM+. The simulations can be adapted to the patient-specific conditions by 1) Utilizing a patient-specific carotid geometry, 2) Specifying the wall thickness in the preprocess- ing step in which the wall is generated (extruded), 3) Supplying patient-specific pa- rameter values for the material model used, and 4) Imposing patient-specific bound- ary conditions, for instance pressure waveforms obtained with an angiocatheter, and flow rate waveforms through PCMRI. Moreover, the analyses can be tailored to the needs of the user in practically every step of the workflow, for example by specifying which assumptions should be made, implementing more hemodynamic parameters, and incorporating additional stabilization if needed. However, the contributed workflow was only tried out for one patient-specific geom- etry (as well as an idealized carotid geometry generated in the built-in CAD tool in STAR-CCM+, but not described herein). Therefore, no guarantee can be provided that the workflow is exhaustive enough to work for all patient-specific geometries, assumptions, and boundary conditions. Nonetheless, the workflow will most likely be able to aid in the process of performing FSI studies not only of the carotid artery, but also to some extent FSI analyses in general and using other software. 39 5. Discussion 5.2 Hemodynamics 5.2.1 General results In Fig. 4.3, it can be seen that low TAWSS mainly is concentrated in the carotid sinus, the predominant location of the baroreceptors. This suggests that this area is especially prone to atheromatous plaque development, which is in agreement with the current knowledge in the field (see Chapter 1 and Section 2.1). Moreover, the highest TAWSS is seen around the carotid apex, the main site of the chemoreceptors. This is to be expected, as the blood must decelerate abruptly before the apex, to subsequently get redirected to the ICA and ECA. Similar results have been reported in several other studies, for instance in [25] and [56]. Comparing Fig. 4.4 to Fig. 4.3, it is evident that the WSS on the whole is much larger at peak systole than the average WSS over one heart cycle. Clearly, this is to be expected, as the mass flow rates are close to their maximum at peak systole, resulting in high velocities. These high velocities, in turn, signify large WSS due to high velocity gradients in the direction normal to the wall. In Fig. 4.5, high OSI appears to be close to the regions of low TAWSS in Fig. 4.3. Although the magnitude and partly also the spatial distribution of elevated OSI are likely to be highly dependent on the patient-specific geometry, it has been shown that high OSI typically can be seen proximal (closer to the bifurcation) in the ICA for healthy subjects, and slightly distal (closer to the outlet) in patients with an ICA stenosis [57]. In Fig. 4.6, the velocity appears to be zero at the wall, satisfying the imposed no-slip condition. Moreover, a slow-moving recirculation region is clearly visible in the carotid sinus, as well as in the proximal ECA. The velocity magnitudes are within a physiologically reasonable range for the carotid artery, see e.g., [25], [56], and [28]. 5.2.2 Impact of FSI and material models Visually, it is difficult to appreciate any major differences in the TAWSS between CFD (Fig. 4.3a) and FSI (Fig. 4.3b). Upon close inspection, it can be noticed that the TAWSS typically seems to be slightly higher in CFD compared to FSI. Comparing the WSS at peak systole between CFD (Fig. 4.4a) and FSI (Fig. 4.4b), this trend is easier to spot. Quantitatively, Tab. 4.2 elucidates an underestimation of the area of low TAWSS in CFD by approximately 2 % and 5 % compared to FSI, in the end diastolic and peak systolic geometric configuration, respectively. The increased 3 percentage points difference in the peak systolic configuration is likely almost exclusively due to the increased area of the carotid artery at this phase in FSI. These results partly oppose those of Lopes et al. in [25], where CFD was shown to overestimate the area of low TAWSS. Intuitively, however, as FSI results in lower velocities and velocity gradients, the TAWSS should also decrease and thus lead to an increased region with low TAWSS compared to CFD. This speaks in favor of the TAWSS trend obtained in this study, especially as Lopes et al. (somewhat contradictory, one might argue) also concluded an overestimation of the TAWSS magnitudes in CFD, despite their reported increased area of low TAWSS in CFD. 40 5. Discussion Regarding the OSI, the influence of a deformable wall has been shown to be more complicated than for TAWSS: in [58], the OSI increased in numerous regions, while in many others, a decrease was observed. The same holds true observing the dif- ferences in OSI between CFD and FSI in Fig. 4.5a & 4.5b. For instance, at the proximal part of the ECA, the FSI simulation appears to have higher OSI magni- tudes and larger areas of elevated OSI than CFD. Slightly above this region, on the other hand, a contour of low – but higher than in FSI – OSI is apparent in CFD. In Tab. 4.2, it can be seen that with FSI, the area of the elevated OSI risk zone is around 7 % larger in the end diastolic configuration compared to CFD. In the end systolic configuration, an approximately 10 % higher area is estimated with FSI. As was the case for the TAWSS, also here the difference between end diastolic and peak systolic areas are approximately 3 percentage points, plausibly indicating that the pure areal increase between the two phases in FSI is roughly 3 % (around the risk areas). Despite the complex nature of the OSI, the fact that the area of risk increases in FSI compared to CFD might be explainable by considering that gener- ally lower WSS is present in FSI, and that likely the WSS vectors are more prone to oscillate (i.e., change direction) when their magnitudes are low. It should be noted that the TAWSS and OSI calculated in the peak systolic configuration were calcu- lated at peak systole of the second cardiac cycle, thus also including values from the first cardiac cycle (following the first peak systole) due to the time-averaging over one period. However, this likely has only minimal, if any, effect on the computed hemodynamic parameters. The first cardiac cycle appeared to have identical values to the second cycle, apart from for a few time steps at the very beginnings of the simulations. Two interesting observations can be made looking at the CFD and FSI velocity scene in Fig. 4.6a & 4.6b, respectively. First, the velocity magnitudes are generally lower in FSI. This is expected, while an increased lumen volume on average must have lower fluid velocities for the same mass flow rate, otherwise mass would not be conserved. Second, the recirculation region by the carotid sinus is larger in FSI, due to the wall deformation. In Fig. 4.7, it can be seen that the CFD and FSI velocities typically differ more when the pressure (see Fig. 3.4a), and hence wall deformation, is large. As can be anticipated, at the beginning and end of the heart cycle when the geometry is close to identical in both cases, no difference between the velocities is present. On average, throughout the whole fluid domain, FSI results in velocities approximately 4 % lower than in CFD (see Tab. 4.3). In the selected point in the ECA, almost no difference in velocity is obtained with FSI, which is in agreement with Lopes et al. in [25]. Regarding the impact of the various material models on the hemodynamic param- eters, the differences between them appear to be minimal. For instance, looking at the velocity over time plots in Fig. 4.7, the three material models are all close to perfectly aligned in all plots, consistently deviating from the CFD predictions. However, a small but seemingly systematic difference can be seen between the three material models when looking at the quantitative deviations in Tab. 4.2 & 4.3: The deviations from CFD are largest for the Neo-Hookean model, followed by ILE 41 5. Discussion and last Mooney-Rivlin. This is true for all three of the quantified hemodynamic parameters: the velocities, the area of low TAWSS, and the area of high OSI. The only exception is the velocities in the ECA (Tab. 4.3), where the opposite order seems to be the case, i.e., the Mooney-Rivlin model has the highest deviation from CFD, followed by ILE and last the Neo-Hookean model. Given the roughly 20 times longer computation time for running an FSI simula- tion compared to CFD, whether the added accuracy and greater resemblance to physiological reality obtained with FSI is needed should be carefully assessed in each specific use case. In cases where the lowest possible errors in estimated hemo- dynamic parameters are warranted, or where there is interest in investigating the structural behavior of the arterial walls, plaques, or medical devices inserted inside the body, FSI could well be employed. In other cases – especially as other steps of the process, such as image segmentation and boundary condition acquisition, may induce considerably larger inaccuracies in the analyses than a rigid wall assumption – CFD may be preferred due to its higher stability, ease of implementation, and lower computational load. 5.3 Wall deformation Looking at the deformation scenes in Fig. 4.8, it is evident that the largest defor- mation takes place at the bifurcation around the S1 point (defined in Fig. 3.5). In- tuitively, passages with larger cross-sectional areas and lower wall curvature require higher wall stress to prevent them from bursting. This can roughly be explained by the Law of Laplace. Looking at Eq. (2.46), the wall stress in a thin-walled cylinder is proportional to its radius. Moreover, the displacement is proportional to the radius squared, as can be seen in Eq. (2.47). The (hydraulic) radius around the bifurcation is much larger than elsewhere in the artery, namely approximately 5.4 mm compared to 3.1 mm in the CCA, where the S2 and S3 points are positioned. This corresponds to ( rs1 rcca )2 ≈ 3, meaning that according to the Law of Laplace, approximately three times the displacement magnitude in the CCA can be expected around the bifurca- tion. However, the Law of Laplace, as well as the slightly more advanced expression for thick-walled cylinders stated in Eq. (2.48), use major simplifications of reality, as can be seen in Fig. 4.9. In all three points, the analytical predictions are mas- sively off, either underestimating or overestimating the displacement magnitude by roughly ±50% according to Tab. 4.4. The analytical expression may therefore be considered useful for quick but likely highly imprecise estimations of carotid artery deformation. For other blood vessels for which the geometry is less complex (i.e., more cylinder-like), the analytical expression is likely more accurate. In Fig. 4.8, the wall deformation at peak systole looks nearly identical for all three material models. Upon close inspection, it can be seen that the Mooney-Rivlin model in Fig. 4.8c predicts lower displacement magnitudes around the bifurcation than the other two models. This is confirmed in Fig. 4.9a and Tab. 4.4, quantifying the over-estimations of the ILE and Neo-Hookean models in S1 to 4.97 % and 5.56 %, respectively. In the two selected points in the CCA, on the other hand, the Neo- 42 5. Discussion Hookean model agrees almost perfectly with the Mooney-Rivlin predictions, whereas the ILE model deviations are larger. It should be noted that the displacements in the S2 and S3 points are fairly small compared to the S1 point, making the absolute deviations in these points much smaller and conceivably of less interest than the over-estimations in S1. Whereas the difference in magnitude around the bifurcation point can be seen rather easily in Fig. 4.8, the differences in the CCA can only be discerned by carefully observing the spatial extents of the contour lines. The fact that the results from the various material models only vary slightly, likely suggests that the deformations are small enough to still be within the range for which the stress and strain are highly linearly related. The Neo-Hookean model is typically accurate for strains less than 20 % and for uniaxial stretching, whereas the Mooney-Rivlin model is suitable also for larger strains and for biaxial stretching [36]. Partly because the edges around the inlet and outlets are fixed, it is possible that only a marginal amount of longitudinal stretching occurs, effectively making the stress uniaxial. The maximal displacement in our case is approximately 0.55 mm at peak systole, and assuming a radius at the bifurcation of 5.4 mm, this gives (highly simplified) a maximal strain of roughly 10 % in the radial direction, potentially indicating that no large differences between the models should be expected. 5.4 Limitations Apart from the limitations already set at the offset of this project and listed in Section 1.2, additional limitations are present in this thesis. As mentioned in Section 3.6, the boundary conditions were computed from patient- specific measurements from another simulation study of the carotid artery [50]. This, on the one hand, means that physiologically realistic pressure and flow rate wave- forms were used. On the other hand, the boundary conditions are not specific to the patient of which this study was performed. If the areas or radii of the ICA and ECA in the study from which the boundary conditions were taken were available, the flow rates could have been scaled by compa