SIMULATION OF CARDIAC THERMAL AND FLUID-STRUCTURE INTERACTION PHYSICS USING THE FINITE ELEMENT METHOD: METHODS, DEVELOPMENT AND APPLICATIONS By Tejas Patel A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering – Doctor of Philosophy 2023 ABSTRACT The human heart is a highly complex organ, and its primary function is to pump blood through the arteries, veins and to perfuse all other body tissues and organs, including itself. In the last decade, cardiac simulations have become increasingly crucial to gain clinical insight into cardiac function, treatment, and testing. Nowadays, multi-physics cardiovascular simulations applied to patient- specific modeling can help in the diagnosis of cardiovascular diseases and in studying relevant clinical treatments. Hence, our central objective here is to develop a generalized multi-physics finite element (FE) framework that includes thermal-fluid structure interaction coupling to study cardiac function and treatments. First, we developed a stabilized FE based flow solver with heat transfer to study hemodynamics. A python based open-source FE library (FEniCS) is used from ground-up to custom-build the solver. We benchmark and validate the solver and study convergence for classical test cases at intermediate Reynolds and Peclet number. Second, we utilize the solver to investigate cryoballoon ablation (CBA), which is a minimally invasive surgery that uses freezing or cryoenergy to treat atrial fibrillation (AF). To begin with, we use a patient-specific left atrium (LA) geometry and realistic pulmonary vein (PV) blood flow boundary conditions to validate hemodynamics of the LA chamber. Next, we position a cryoballoon (CB) at the pulmonary vein ostium to simulate incomplete occlusion during cryotherapy and investigate the factors affecting lesion formation. We observe that lesion size is highly sensitive to the CB position and balloon tissue contact area. The threshold gap for lesion formation is 2.4 mm. We also note that as the balloon tissue contact area increases, the surgery is more effective, and the power absorbed across the CB reduces. Third, we extend our development to a fully coupled fluid-structure interaction (FSI) solver with heat transfer using FEniCS. The FSI solver (named vanDANA) that uses the immersed boundary (IB) method is based on the Distributed Lagrange Multiplier based Fictitious Domain method and the interpolation of variables is conducted using the smeared delta-functions. Additionally, the structure can be set as incompressible or compressible. We benchmark our solver and analyze the scalability on HPC. This builds a solid foundation for the future use of this solver. Copyright by TEJAS PATEL 2023 This thesis is dedicated to the patience and all the sacrifices made by my mom and dad, to friends in Michigan, and to my personal growth during this journey. Thank you for always being present - knowingly/unknowingly for me. v PREFACE The content included in Chapter 3, Chapter 4, Appendix A and Appendix B of this thesis has contributed towards the following publication: Patel, Tejas, Chris Li, Farshad Raissi, Ghassan S. Kassab, Tong Gao, and Lik Chuan Lee. "Coupled thermal-hemodynamics computational modeling of cryoballoon ablation for pulmonary vein isolation." Computers in Biology and Medicine (2023): 106766. vi TABLE OF CONTENTS CHAPTER 1 : GENERAL BACKGROUND ................................................................................. 1 1.1 Heart Anatomy and cardiac cycle diagram ..................................................................... 2 1.2 Background of this dissertation ...................................................................................... 5 1.3 Objectives of this dissertation ....................................................................................... 13 CHAPTER 2 : REVIEW OF THE STATE OF ART ..................................................................... 14 2.1 Challenges, complications, and the need for cryotherapy simulations ......................... 15 2.2 Fluid structure interactions: complexities, algorithms, and applications ...................... 17 CHAPTER 3 : DEVELOPMENT OF STABILIZED FINITE ELEMENT MODEL TO INVESTIGATE THERMAL-HEMODYNAMICS IN CARDIAC CHAMBERS ....................... 21 3.1 Methods......................................................................................................................... 22 3.2 Numerical details .......................................................................................................... 29 3.3 Benchmarks................................................................................................................... 33 CHAPTER 4 : PATIENT SPECIFIC COMPUTATIONAL MODELING OF CRYOBALLOON ABLATION FOR PULMONARY VEIN ISOLATION ............................................................... 36 4.1 Patient specific geometry and physical setup ............................................................... 37 4.2 Validation for Left Atrium hemodynamics ................................................................... 40 4.3 LA+CB results .............................................................................................................. 42 4.4 Discussion ..................................................................................................................... 50 4.5 Conclusion .................................................................................................................... 53 CHAPTER 5 : DEVELOPMENT OF IMMERSED BOUNDARY BASED FICTITIOUS DOMAIN METHOD FOR MODELING OF FLUID STRUCTURE INTERACTIONS WITH HEAT TRANSFER ...................................................................................................................... 54 5.1 Methods......................................................................................................................... 55 5.2 vanDANA solver........................................................................................................... 63 5.3 Benchmarks................................................................................................................... 67 5.4 HPC performance.......................................................................................................... 76 CHAPTER 6 : SUMMARY .......................................................................................................... 78 6.1 Limitations .................................................................................................................... 79 6.2 Concluding remarks ...................................................................................................... 80 BIBLIOGRAPHY ......................................................................................................................... 82 APPENDIX A : EFFECTS OF STABILIZATION SCHEMES .................................................. 91 APPENDIX B : CRYOBALLOON POSITIONING................................................................... 92 vii CHAPTER 1 : GENERAL BACKGROUND 1 1.1 Heart Anatomy and cardiac cycle diagram The human heart is a highly complex organ, and its primary function is to pump blood through the arteries, veins and to perfuse all other body tissues and organs, including itself. This process facilitates the exchange of gases, fluids, electrolytes and heat exchange between the cells and the ambient environment [1]. The cardiac anatomy mainly comprises of two atriums and two ventricles separated by atrio-ventricular (AV) valves, and other arteries and veins that deliver blood to and Figure 1.1: a) Heart anatomy with four chambers, two atria and two ventricles, separated by AV valves. The red color on the left and blue color on the right side depicts the flow of oxygenated blood and deoxygenated blood respectively. b) The Wiggers diagram. It demonstrates pressure-volume relationship for different heart chambers in sync with a healthy electrocardiogram (ECG). Figure adopted from Klabunde [1]. away from the heart chambers (see Figure 1.1a). The right heart comprises of the right atrium (RA) and right ventricle (RV) that are connected by a tricuspid AV valve and the left side comprises of the left atrium (LA) and the left ventricle (LV) that are connected by a mitral valve (MV). These valves function with the help of chordae tendineae and papillary muscles. Functionally, the heart 2 can be viewed as two pumps connected in series, one responsible for systemic circulation and the other one for pulmonary circulation. The systemic circulation comprises of the blood vessels inside and outside of all the organs excluding the lungs and the pulmonary circulation denotes blood flow within the lungs that allows the exchange of gases and oxygen in the alveoli. Deoxygenated blood collected from the body in the systemic circulation passes through the inferior and superior vena cava, and then enters the RA and the right ventricle (RV). The deoxygenated blood is then pumped by the RV into the pulmonary artery and lungs for oxygenation. The oxygenated blood from the lungs enters the LA through the pulmonary veins (PV’s) in the left heart. The oxygenated blood enters the LV, which then pumps the blood into the aorta and systemic circulation. The pumping capacity of the heart is usually expressed by cardiac output = amount of blood pumped during each contraction (stroke volume) * heart rate (HR). In Figure 1.1b, the cardiac cycle diagram (aka. The Wiggers Diagram) illustrates the pressure and volume waveforms in the left heart, which is normally divided into two different phases, namely systole and diastole. Systole refers to ventricular contraction and ejection whereas diastole refers to ventricular relaxation and filling. In the left heart, the MV closes when pressure in the LV is higher than that in the LA, marking the beginning of systole. This onset of systole is denoted by the QRS complex which represents ventricular depolarization. Next, with all the valves closed, the intraventricular pressure begins to rise until there is a pressure difference across the LV and aorta. This pressure difference causes the aortic valve to open, and blood is then ejected out of the ventricles. About approximately 180 ms after the QRS complex, a T-wave initiates ventricular repolarization. The ventricle muscle starts to relax, and the intraventricular pressures start to drop. This marks the beginning of diastole and soon after, the aortic valve closes, and the LV volume remains constant (since all valves are closed). This volume is called End-Systolic Volume (ESV). 3 In healthy humans, the ESV is approximately 50 mL. Pressure in the LV continues to drop, and when the LV pressure falls below the LA pressure, the MV opens. This initiates ventricular filling which is mainly divided into two phases: passive filling without atria contraction and filling because of atria contraction. Passive filling occurs before the P-wave in the ECG, and this accounts for the bulk of LV filling. As the LV continues to fill up with blood, they become stiffer (less compliant) and the intraventricular pressures start to rise. The P-wave in the ECG is associated with electrical depolarization of the atria, which results in atrial contraction and an increase in LA pressure as marked by the a-wave in the LA pressure curve in Figure 1.1b. Atria contraction (sometimes known as atrial kick) forces more blood into the LV before the MV closes at the end of diastole. Filling due to atrial contraction accounts for only 10% of ventricular filling. This behavior is reflected in the spectral doppler MV velocity graph in Figure 1.2, wherein the E peak represents rapid filling during LV relaxation and A peak represents filling due to the atrial kick. In heathy human subjects, the average E/A ratio ranges between 0.75 to 1.5. Figure 1.2: Normal Mitral Valve (MV) inflow velocity measured for a healthy subject using spectral doppler echocardiography. Here Deceleration Time (DecT) is defined as the time taken for the maximum velocity at the E peak to drop to zero. Figure adopted from [2]. 4 1.2 Background of this dissertation 1.2.1 Atrial Fibrillation Atrial fibrillation (AFib) is one of the most common chronic cardiac arrythmias where the upper chambers of the heart (atria) beat out of synchronization with the lower chambers i.e., the ventricles. This increases the risk of stroke and heart failure, and the cardiac output reduces by approximately 30%. In healthy humans, the ECG comprises of the P-wave, QRS complex and the T-wave. For patients suffering from AFib however, the P-waves are missing, and the heartbeat is highly erratic, see Figure 1.3. In these patients, electrical depolarization currents arise from multiple sites throughout the atria, resulting to quiver and uncoordinated beating with no discernable P-waves. The symptoms of AFib range from nothing to fatigue, shortness of breath, dizziness, palpitations, or chest pain. Although AFib is not life threatening by itself, it can be seriously uncomfortable, is associated with many cardiovascular diseases and requires medical treatment [3]. There are different classifications of AFib [4]: • Paroxysmal Atrial Fibrillation: Intermittent episodes occur and usually stop within 24 hours, without any treatment. • Persistent Atrial Fibrillation: Last longer than seven days and may or may not cease on its own. • Long-standing Atrial Fibrillation: Last longer than a year and changes the cardiac structure. • Permanent Atrial Fibrillation: Permanent and can cannot be treated by medical practices and can severely diminish quality of life. 5 Depending on the patients’ symptoms, AFib is treated with medicine to control the heart rhythm and using blood thinners to prevent clotting. If medications do not help, the patient often undergo cardioversion or a cardiac ablation treatment. (a) (b) Figure 1.3: Pathway of electrical signals and ECG of (a) a healthy patient vs (b) an AFib patient. Note how irregular electrical impulses are generated at the PV opening in case of an AFib heart. 1.2.2 Cryoballoon ablation In AFib patients, irregular electro-physiological signals usually begin where the PV is attached to the LA. In these patients, it becomes very important to cease the erratic beating and hence the transmission of irregular electrical activity from the PV to the rest of the cardiac muscle in the atria. There are non-invasive (e.g., cardioversion) and invasive AFib treatments that include catheter ablation or the maze procedure. Cardioversion supplies low-energy shocks to regulate a 6 patient’s heartbeat [5]. However, this is only a temporary solution and will not cure the underlying cause. There are two types of minimally invasive procedures for treating AFib, namely, radiofrequency ablation (RFA) and cryoballoon ablation (CBA). Radiofrequency ablation or point- by point focal ablation is a more traditional approach that uses heat to destroy the heart tissue. Cryoballoon ablation, on the other hand, is a recent cryo-energy based technique. It has become one of the most prominent surgical procedures to treat AFib [6]. Medtronic, based out of Minneapolis MN, is the leading producer of devices used for CBA surgery. Since the release of CB2 (second-generation cryoballoon) by Medtronic [7, 8]: Arctic Front Advance in 2012 and Arctic Front Advance Pro in 2020, efficiency of CBA surgery has increased drastically. Figure 1.4: Cryoballoon Ablation catheter. For CBA, vascular access is obtained at the femoral vein and the LA is instrumented via a trans-septal puncture. Figure taken from Patel et al. [9]. Ideally, CBA is conducted stepwise, first for the left PV’s and then for the right PV’s. During CBA, the electrophysiologist positions the balloon and then keep a heavy forward pressure to occlude the PV ostium (see Figure 1.4) [8, 9]. Before ablation, the cryoballoon (CB) is manually flushed back and forth in the PV antrum to eliminate any trapped air-bubbles in the blood stream. An 8- pole circular recording catheter attached to the frontal tip of the CB continuously monitors PV 7 potentials in real time. By exploiting the Joule-Thompson effect, liquid nitrous oxide is injected in the CB and a phase change from liquid-to-gas gradually absorbs heat from the LA. When the temperature in the cardiac tissue falls below 23! C, electrical dormancy is observed, and isolation is confirmed by the loss of PV potentials. A conduction delay (or a conduction block) is normally observed at the catheter [10] and the time taken to achieve this isolation is called time to isolation (TTI), which is normally between 30-60 sec. After TTI, the physiologist ideally holds the CB in place as per the recommended dosing protocol - about 2 to 4 minutes per freeze. Permanent cellular damage is detected at temperatures below −20! C, if held for 60 seconds. On the other hand, if the temperature is less than −40! C (necrosis temperature), the hold duration is not important, and the cardiac tissue becomes irreversibly nonviable. At this stage, irreparable lesion formation is initiated at the balloon-tissue contact surface and for an optimal scar/lesion, temperature monitored at the CryoConsole normally ranges between −40! C to −50! C [10]. However, this is only the return cryogenic gas temperature and not the actual temperature of the cryoballoon. Nadir temperatures at the CB surface can go as low as −60! C to −70! C. The most important factors controlling lesion thickness and depth for a successful CBA surgery is namely, balloon-tissue contact area (particularly no leakage at all), time of freeze and the source of cryo-energy. Overall, CBA surgery forms a contiguous lesion and is less time consuming (lower TTI) compared to RFA, which means that the patient is on general anesthesia for a shorter duration. This makes CBA safer than RFA. Cryoballoon ablation has also led to more predictable surgeries with less reliance on electrophysiologist skills and improved quality of life post-surgery. Due to its obvious advantages over RFA, CBA has become the preferred choice of ablation surgery to treat AFib patients. 8 1.2.3 Fluid Structure Interactions Fluid structure interaction (FSI) is the interaction of single or multiple flexible structures (bodies) in adjacent fluid media. Depending on the physical nature, geometry and the incompressibility of the media, the structural deformation can be insignificant or large. Such interactions between the fluid and structure can be found in many day-to-day applications such as aerodynamic flows, bioengineering, marine applications, and turbomachinery. Almost all bio-driven physics in the human body e.g., hemodynamics of aneurysms, arterial stenosis, blood flow in the heart chambers, circle of Willis, and the flow of urine in kidneys and the bladder involve FSI. Developing numerical models to study FSI is crucial and can be used as a practical tool to probe the physical behavior of biosystems. In problems involving FSI, one needs to solve the Navier-Stokes (NS) equations [11] for fluid flow as well as the structural mechanics equations describing deformation of the flexible body. The fluid and structure equations are coupled in both directions, meaning each medium exerts a force on the other and subsequently influences their temporal behavior. For example, in vortex induced vibrations (VIV) over bluff bodies, a condition is reached when the vortex shedding frequency is close enough to the body’s natural frequency and the body begins to oscillate naturally. The forces resulting from vortex shedding cause the body to oscillate and subsequently the body’s oscillation affects the nature of the flow around it. Hence to build a stable design, one needs to consider vortex shedding, forces acting on the body and the natural frequency and amplitude of the body’s oscillation at various flow speeds. Vortex induced vibrations is a very fundamental FSI phenomenon and is widely studied to investigate vibrations caused due to flow over cables, bridges, heat exchanger tubes, ocean drilling risers etc. Another widely studied FSI is the cyclic pumping of blood through arteries. Arteries are compliant and deformable, and the blood flow exerts pressure on the walls causing them to stretch 9 and deform. However, if the shear stress is too high or cyclic it can lead to endothelial dysfunction that may contribute to the development of atherosclerosis. This is because increased shear stress can cause damage to the endothelial cell lining, leading to the formation of plaque and narrowing of the blood vessel. Simulations assuming that the arterial wall is rigid usually overestimates the true wall shear stress [12]. To address this issue, simulations should consider FSI to accurately estimate the wall shear stress and help identify areas of the arterial system that are at a higher risk of developing endothelial damage. Figure 1.5 shows an idealized flow pattern over plaque in an artery. Here, on the distal side, the shear stress is oscillatory and this compounds to more plaque formation over time. Also in real practice, it is fundamentally imperative to include patient specific geometries for such simulations. Overall, understanding how shear stress affects arterial walls can help identify risk factors associated with diseases and help optimize the design of treatments. Figure 1.5: Ideal distribution of wall shear stress an arterial segment with atherosclerotic plaque formation. Figure adopted from Young Cho et al [13]. 1.2.4 Coding in FEniCS The FEniCS project (https://fenicsproject.org [14]), initiated in 2003, is an open-source python platform to solve partial differential equations (PDE’s) using the finite element method (FEM). The backend of FEniCS mainly comprises of fundamental building blocks: DOLFIN, FIAT, FFC (FEniCS Form Compiler), UFL (Unified Form Language) and UFC (Unified Form Assembly 10 Code). FEniCS includes automated code generation capabilities, which generate C++ code for a variational weak form written in high level UFL language. This eliminates the need for manually coding complex subroutines and allows one to quickly generate code in just a few lines of python. This permits quick prototyping, debugging, and testing of the physical problem at hand. Moreover, DOLFIN also provides a set of preconditioned Krylov solvers and PETSc as a linear algebra backend, that support easy assembly and conversion to matrix-vector form before solving it iteratively. In addition, for experienced programmers, FEniCS offers straightforward compilation of custom-built C++ codes into their high-level coding structure. Due to all these advantages, we choose FEniCS to develop our FEM codes for modeling FSI. Figure 1.6 shows an example FEniCS script for solving a 3D elastic deformation problem in case of a hyper-elastic material. For post-processing, the results can be easily visualized on the opensource ParaView (Kitware Inc.) [15]. 11 from dolfin import * # Create mesh and define function space mesh = UnitCubeMesh(24, 16, 16) V = VectorFunctionSpace(mesh, "Lagrange", 1) # Define functions du = TrialFunction(V) # Incremental displacement v = TestFunction(V) # Test function u = Function(V) # Displacement from previous iteration B = Constant((0.0, -0.5, 0.0)) # Body force per unit volume T = Constant((0.1, 0.0, 0.0)) # Traction force on the boundary # Kinematics d = u.geometric_dimension() I = Identity(d) # Identity tensor F = I + grad(u) # Deformation gradient C = F.T*F # Right Cauchy-Green tensor # Invariants of deformation tensors Ic = tr(C) J = det(F) # Elasticity parameters E, nu = 10.0, 0.3 mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu))) # Stored strain energy density (compressible neo-Hookean model) psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2 # Total potential energy Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds # Compute first variation of Pi (directional derivative about u in the direction of v) F = derivative(Pi, u, v) # Solve variational problem solve(F == 0, u, bcs, form_compiler_parameters=ffc_options) # Save solution in VTK format file = File("displacement.pvd"); file << u; Figure 1.6: Simple FEniCS script to solve a 3D elasticity problem for a hyper-elastic cube. This python code is copied from the FEniCS tutorial: demo hyperelasticity.py. 12 1.3 Objectives of this dissertation This dissertation is organized as follows: 1. Chapter 3: Develop a computational modeling framework in FEniCS to investigate flow and heat transfer to investigate cardiovascular flows. 2. Chapter 4: Simulate patient specific CBA surgery for pulmonary vein isolation. • Validate present computational framework for patient-specific left atrium hemodynamics. • Investigate the factors affecting CBA surgery. • Analyze temperature distribution and efficiency for various cryoballoon positions in the PV antrum. 3. Chapter 5: Develop an immersed boundary (IB) distributed Lagrange multiplier (DLM) based fictitious domain (FD) method in FEniCS to model FSI for biophysical flows. • Create a user-friendly open-source solver and publish/document it on GitHub: vanDANA (https://github.com/patelte8/vanDANA). • Validate and measure HPC-scalability of vanDANA solver (IB-FSI algorithm) using classical benchmark problems. 13 CHAPTER 2 : REVIEW OF THE STATE OF ART 14 2.1 Challenges, complications, and the need for cryotherapy simulations One of the major challenges in a CBA surgery is to achieve complete occlusion at the PV ostium. In case of a leakage at the PV ostium, energy (cold) is convected into the blood stream, and this hinders lesion formation. Another critical challenge is to monitor the real time evolution and distribution of temperature (cold) in the PV antrum. This is very critical in preventing complications like the phrenic nerve and esophageal injury. Phrenic nerve injury (PNI) is not a common occurrence in RFA. Due to the anatomical proximity of the phrenic nerve to the right PV’s and more contiguous lesions in CBA, however, many studies [16, 17, 18] have reported that the risks of PNI for CB2’s are increased up to 8% - 10%. On the other hand, esophageal injury is observed for both RFA and CBA due to proximity of the esophagus to the posterior LA wall (see Figure 2.1). Hence, it is very important for the electrophysiologist to monitor real time temperature in the esophagus by placing an esophageal probe near the CB [19]. If the balloon and tissues are not in proper contact, then TTI could be longer than 60 sec or there could be no-isolation at all. In such conditions, a better and safer strategy is to abort the freeze and reposition the CB at another position on the PV ostium. Considering all the risks involved, optimum CB positioning is very crucial for an effective CBA surgery [20]. Figure 2.1: Figure adopted from Sarairah S.Y. et al. [19]. 15 Today with the growth in computational modeling, these challenges to CBA can be addressed by conducting simulations of the temperature field for patient specific geometries prior to the surgery. Numerical models to simulate cryotherapy can be developed to investigate patient specific LA geometries. Some numerical studies have been conducted to investigate the thermal field in the PV antrum. Getman et. al. [7] studied the relationship between TTI and freeze duration where they modeled the temperature at various tissue depths using COMSOL. Their simulations for both the 28 mm CB2’s: Arctic Front Advance and Arctic Front Advance Pro demonstrated similar TTI’s across various cardiac tissue depths up to 5 mm, where they observed a TTI of about 60 sec for a lesion depth of 3 mm. Similarly, Xia et al. [21] used a more simplistic 2D-axisymetric geometry for the balloon-PV contact surface and they compare their numerical findings with an in-house experimental balloon setup surrounded by a freshly cut porcine myocardium tissue. They conducted a thorough tissue damage analysis and found that if the temperature of the tissue drops below −40! C (necrosis temperature), the duration of hold is unimportant due to direct cell damage [22]. Whereas, for temperatures less than −20! C, they found that it is important to hold the CB for at least 60 sec to produce necrosis of the tissue. They observe that the necrosis −40! C isotherm only reaches a depth of about 2 mm after 240 sec of ablation. At depths greater than 2 mm, the hold time becomes decisive to lesion/scar formation. Both these studies [7, 21] model the cold by solving for the Pennes heat conduction (or energy) equation in time and accounting for convection due to blood flow in the LA antrum. Today, with the advent of CBA surgery and the diverse range of governing factors affecting the surgery, there is a need for patient-specific computational modeling within a LA geometry to simulate the coupling between blood flow and heat transfer during CBA. Our goal mainly lies in 16 investigating how the presence of a CB affects hemodynamics in the LA, and how incomplete occlusion (leakage) affects temperature distribution at the PV ostium and lesion formation. 2.2 Fluid structure interactions: complexities, algorithms, and applications To simulate real-world FSI, requires a robust, accurate, and flexible algorithm that can handle complex 3D moving geometries. In such geometries, the implementation of boundary conditions and mesh generation are major challenges and also the fundamental reason that influence the accuracy of results. Mesh generation is critical especially when the objects are moving or deforming and if not done efficiently, can lead to high computational overheads [23]. Most FSI mesh-based algorithms can be classified into two families: boundary-fitted methods and non- boundary-fitted methods, according to whether the boundary-fitted mesh is used for the solution of the flow field [24, 25]. The boundary fitted method is the traditional approach where the mesh conforms to the boundary of the domain and remeshing is required as the interface deforms with time. The classical boundary fitted method is Arbitrary-Lagrangian Eularian (ALE) method [26]. Besides a fine grid resolution at the FSI interface, the ALE method also requires discretization of the governing equations on grids aligned with the boundary. The other class of methods are the non-boundary fitted methods, where the flow field is computed on a stationary Eulerian grid and remeshing is not required for the fluid grid. This saves a considerable chuck of the computational cost in terms of remeshing the flow grid. Sometimes, remeshing might be required for the structural grid if the extent of deformation is large which is computationally less expensive as the solid grid is less in volume as compared to the fluid grid. The immersed boundary (IB) and fictitious domain (FD) algorithms are examples of non-boundary fitted methods. The IB method typically loses accuracy at the solid-fluid interface due to the numerical difficulty in constructing a conservative mapping function between the Eulerian (fluid) mesh and 17 the Lagrangian (solid) mesh. As a result, IB methods cannot resolve the fine details of boundary flow for accurate shear stress computation. Today, there are two different FD based schemes in the literature: non-body-force-based scheme and the body-force-based scheme. The method used by Farnell et al. [27] for the simulation of a filament in a flowing soap film belongs to the first category. The second, is the IB method proposed by Peskin [28] which is a body force based FD method that has been applied to a wide range of FSI problems (Zhu and Peskin [29], Eggleton and Popel [30]). In this method, the flexible body moves at the same velocity as the local fluid, and then affects the fluid motion through an elastic force that is calculated with the known deformation of the body and is introduced into the fluid momentum equation as a pseudo body force. Here the key idea is to virtually extend the exterior fluid domain into the solid where a distribution of pseudo body forces is employed to enforce the structural movement through Lagrange multipliers. Such methods have been implemented by Kamensky et al. [31], Kapada et al. [32] and Yu [33]. Instead of using the macroscopic NS equations to solve for the fluid flow, one can flexibly combine particle based methods with the immersed boundary method. Particle-based methods model the fluid as a collection of discrete particles rather than as a continuous fluid. Two commonly used particle-based methods are Smoothed Particle Hydrodynamics (SPH) and Lattice Boltzmann Method (LBM). SPH [34] is a Lagrangian method that models the fluid as a collection of particles, where each particle is assigned a mass, velocity, and density. The fluid equations are solved by computing the forces between the particles using a kernel function that smooths out the interaction between neighboring particles. LBM [35, 36, 37, 38] is a mesh-free method that discretizes the fluid domain into a lattice of cells, where each cell represents a small volume of fluid. The fluid equations are solved by tracking the distribution of particles within each cell using a probabilistic approach [39]. LBM is well-suited for simulating complex fluid behavior such as 18 turbulence and multiphase flows. Recently, particle-based methods have become more and more popular because it offers advantages in terms of parallelization (to efficiently model large-scale simulations), simple implementation, and robustness for FSI applications. Models describing FSI in cardiovascular systems In the last two decades, the use of FSI modeling has increased in cardiovascular engineering, especially in the field of arterial blood flow and cardiac hemodynamics. Computational models of FSI are used to model the abdominal aorta [40], carotid bifurcation [41], heart valves [42] and cerebral aneurysm [43]. Torii, Wood et al. [44] has studied the effects of wall compliance on a patient-specific right coronary artery with severe stenosis and found noticeable differences in the instantaneous wall shear stress produced by the FSI and rigid wall models. Apart from the difference in patient geometries, mechanical properties of the vessel wall also vary significantly. The larger arteries are more elastic than the smaller ones. Deformations up to 10% of the vessel radius are common in large arteries. This elasticity gives rise to the Windkessel effect and provides a passive mechanism for smoothing the pulsatile blood flow from the heart. The Windkessel effect is mainly attributed to the distension of large arteries during systole, functioning as a buffering reservoir for blood [45]. Besides arteries, lately FSI models have also been used to simulate heart valve dynamics [46] and cardiac chamber flows [47]. The earliest coupled heart-valve FSI model was developed by Peskin and McQueen’s in the 1970s [48, 49, 50] using the classical IB approach. Chandran and Kim [51] reported a study with MV dynamics in a simplified LV chamber model during diastolic filling using an immersed interface-like approach. One of the main limitations of this coupled models is the simplified representation of the biomechanics of the LV wall. Recently, Gao et al. [47] have developed a MV–LV model which has full FSI and is based on a realistic geometry and an experimental model for soft tissue mechanics. Such complex FSI simulations 19 require a robust multiphysics coupling algorithm that allows large deformations and incorporates turbulence modeling and contact mechanics. Given the complexities and challenges, we aim to develop an open-source robust IB based FSI solver using the delta-functions which can simulate any complex 3D geometry. Moreover, we also intend to benchmark, document, and measure the scalability of our IB-FSI algorithm on HPC. This will enable us to create a solid foundation for the future use of this solver. 20 CHAPTER 3 : DEVELOPMENT OF STABILIZED FINITE ELEMENT MODEL TO INVESTIGATE THERMAL- HEMODYNAMICS IN CARDIAC CHAMBERS 21 3.1 Methods 3.1.1 Thermal-Hemodynamics coupling model Hemodynamics and heat transfer in the LA are described by solving the continuity, incompressible NS and energy conservation equations. To non-dimensionalize the governing equations [11], we introduce a characteristic length 𝐿"# and velocity 𝑉"# to define the non-dimensional control $!" &!" ' parameters (ℂ): Reynolds number 𝑅𝑒 = ' and Prandtl number 𝑃𝑟 = (, where 𝜐 is the kinematic viscosity (m2/s) and 𝛼 is the thermal diffusivity (m2/s) of the fluid. We also define the non-dimensional Peclet number 𝑃𝑒 = 𝑅𝑒 ∙ 𝑃𝑟, which is a measure of the ratio of advective to diffusive transport of temperature. We restrict our modeling to moderate 𝑅𝑒 and assume laminar flow everywhere. Correspondingly, the complete non-dimensional governing equations is given as follows: In the fluid domain Ω: 𝜕𝒖 + 𝒖 ∙ ∇𝒖 − ∇ ∙ 𝛔 − 𝒇 = 𝟎, (3.1) 𝜕𝜏 ∇ ∙ 𝒖 = 𝟎, (3.2) )* ,#* (3.3) )+ + 𝒖 ∙ ∇Τ − -. − 𝑄 = 0, Constitutive model - Newtonian Fluid: 2 𝝈(𝒖, 𝑝) = −𝑝𝑰 + 𝑬(𝒖), (3.4) 𝑅𝑒 / (3.5) 𝑬(𝒖) = 0 (∇𝒖 + (∇𝒖)1 ), On the boundaries: 𝒖 = 𝒖𝛚 𝑎𝑡 𝜕Ω𝒖3 ; 𝛁𝒖 ∙ 𝒏 = 𝒖𝛚 𝑎𝑡 𝜕Ω𝒖5 , (3.6) 22 𝑝 = 𝑝6 𝑎𝑡 𝜕Ω73 ; ∇𝑝 ∙ 𝒏 = 𝑝6 𝑎𝑡 𝜕Ω75 , (3.7) Τ = Τ6 𝑎𝑡 𝜕Ω*3 ; ∇Τ ∙ 𝒏 = Τ6 𝑎𝑡 𝜕Ω*5 . (3.8) In the above equations, 𝒖, 𝑝 and Τ are the primary unknown variables defining the fluid velocity, pressure, and temperature fields, respectively. The blood is assumed to behave as a Newtonian fluid. Correspondingly in Eqs. (3.4) and (3.5), 𝛔 and 𝑬 are, respectively, the stress and strain rate tensors. In Eq. (3.1), 𝒇 is the prescribed body force, 𝜏 is the non-dimensional time and ∇(∙) is the gradient operator, whereas in Eq. (3.3), 𝑄 is a prescribed heat source. The Dirichlet and Neumann conditions are prescribed for the unknown variables 𝒖, 𝑝 and Τ at the boundary in Eqs. (3.6-3.8), with 𝒏 denoting the vector normal to the boundary and (𝒖𝛚 , 𝒖𝛚 , 𝑝6 , 𝑝6 , Τ6 ,Τ6 ) 𝝓 𝝓 denoting the prescribed values. In these equations, the boundary is denoted by 𝜕Ω𝝓 = 𝜕Ω3 ∪ 𝜕Ω5 𝝓 𝝓 such that 𝜕Ω3 ∩ 𝜕Ω5 = { }, where 𝝓 ⊆ { 𝒖, 𝑝, Τ }. 3.1.2 Stabilized Finite element method The FE method is used for spatial discretization of the governing equations. Denoting 𝒗, 𝑞 and 𝑤 as the respective test functions for 𝒖, 𝑝 and Τ, the resultant Galerkin weak formulation is given as follows: 𝝏𝒖 〈 , 𝒗〉 + 〈𝒖 ∙ 𝛁𝒖, 𝒗〉 + 〈𝛔, 𝛁𝒗〉−〈𝛔 ∙ 𝒏, 𝒗〉𝝏𝛀 − 〈𝒇, 𝒗〉 = 𝟎, 𝝏𝝉 (3.9) 〈𝛁 ∙ 𝒖, 𝑞〉 = 0, (3.10) 𝜕Τ 1 〈 , 𝑤〉 + 〈𝒖 ∙ ∇Τ, 𝑤〉 + \〈∇Τ, ∇𝑤〉 − 〈Τ6 , 𝑤〉);% ] − 〈𝑄, 𝑤〉 = 0, 𝜕𝜏 𝑃𝑒 $ (3.11) 23 where {𝒖 ∈ 𝐻/ (Ω), 𝒖 = 𝒖𝛚 𝑎𝑡 𝜕Ω𝒖3 }; `𝑝 ∈ 𝐻/ (Ω), 𝑝 = 𝑝6 𝑎𝑡 𝜕Ω73 a; {Τ ∈ 𝐻/ (Ω), Τ = Τ6 𝑎𝑡 𝜕Ω*3 }; {𝒗 ∈ 𝐻/ (Ω), 𝒗 = 𝟎 𝑎𝑡 𝜕Ω𝒖3 }; `𝑞 ∈ 𝐻/ (Ω), 𝑞 = 0 𝑎𝑡 𝜕Ω73 a and {𝑤 ∈ 𝐻/ (Ω), 𝑤 = 0 𝑎𝑡 𝜕Ω*3 }. In the above equations, 〈∙ ,∙ 〉 and 〈∙ ,∙ 〉𝝏𝛀 denote the inner product of two functions i.e., 〈𝑧/ , 𝑧0 〉 = ∫𝛀 𝑧/ 𝑧0 ∙ 𝑑Ω and 〈𝑧/ , 𝑧0 〉); = ∫); 𝑧/ 𝑧0 ∙ 𝑑s for functions 𝑧/ and 𝑧0 . Because of the high 𝑅𝑒 (~10< ), 𝑃𝑒(~10= ) and velocities, hemodynamics in the heart chambers is dominated by advection. As such, it is necessary to incorporate appropriate stabilization schemes to the Galerkin weak formulation to ensure numerical stability [52]. To do so, we applied the well-established Streamline-Upwind Petrov Galerkin (SUPG) [53, 54] and the crosswind stabilization scheme developed by Codina [55]. The SUPG stabilization scheme is applied by introducing the following to Eqs. (3.9) and (3.11): ∫; 𝛾>?-@ (𝒖) 𝑃(𝒖, 𝒎) ∙ 𝑹𝝓 𝑑Ω, (3.12) where, & 0 A 0 0 # 2 2 o|𝒖|o 4 𝛾>?-@ (𝒖) = 𝛼 jk m + n r + 9 k 0m t , (3.13) ∆𝜏 ℎ ℂℎ 𝑃(𝒖, 𝒎) = 𝒖 ∙ ∇𝒎 . (3.14) In the above, ℎ denotes the cell diameter, ℂ denotes the respective control parameter (i.e., 𝑅𝑒 or 𝑃𝑒) and the coefficient 𝛼 is prescribed with a value of 0.85. Next, we introduce another residual pressure-based stabilization term to Eq. (3.9), namely the pressure stabilizing Petrov Galerkin (PSPG) scheme [19] that is given as: 𝜕𝒖 u 𝛾->-@ (𝒖) ∇𝑞 ∙ k + 𝒖 ∙ ∇𝒖 − ∇ ∙ 𝛔 − 𝒇m 𝑑Ω, (3.15) ; 𝜕𝜏 24 where 𝛾->-@ = 𝛾>?-@ . The crosswind stabilization scheme adopted from Codina [18] is also applied by introducing the following as additional term to the weak formulation into Eqs. (3.9) and (3.11): u v𝛾BC (𝝓) Λ(𝒖, 𝝓)x ∶ 𝛁𝒎 𝑑Ω, (3.16) ; where Λ(𝒖, 𝝓) = 𝚪(𝒖) ∙ ∇𝝓, (3.17) 1 |o𝑹𝝓 o| 𝛾BC (𝝓) = 𝛽# ℎ (3.18) 2 o|∇𝝓|o 𝒖⨂𝒖 𝑰 − |𝒖|𝟐 𝑖𝑓 𝒖 ≠ 𝟎 In the above, 𝚪(𝒖) = } •, 𝝓 is the unknown variable, and 𝛽# = max k0, 𝜁 − 𝟎 𝑖𝑓 𝒖 = 𝟎 0 F|𝛁𝝓|F ℂ I JF𝑹𝝓 FJ m [17]. In both Eqs. (3.12) and (3.16), 𝑹𝝓 and 𝒎 are the residual and test function. The crosswind stabilization scheme works in such a way that a non-linear numerical dissipation proportional to the element residual is added in the crosswind direction in regions where the convective effects are dominant. It should be noted that SUPG, PSPG and crosswind stabilizations are highly consistent because they are all residual based terms and vanish as the residual approaches to zero. Next, we note that when the control parameters (𝑅𝑒, 𝑃𝑒) are large, applying only the SUPG, PSPG and crosswind stabilization schemes is not sufficient to prevent non-physical oscillations. Hence, to further increase stability and accuracy of the flow solver, we also include two additional stabilization terms to Eq. (3.9). Specifically, the Least-Squares on Incompressibility Constraint (LSIC) stabilization is used to improve the accuracy and conditioning of mass balance in the LA. This stabilization scheme is implemented by adding the following to Eq. (3.9): 25 ∫; 𝛾&>LB (∇ ∙ 𝒖) (∇ ∙ 𝒗) 𝑑Ω, (3.19) 0 where 𝛾&>LB = < M.. To address the numerical instability arising from flow reversal at the boundaries with high 𝑅𝑒 that is common in biological flow problems [20], we include the following penalty-based backflow stabilization term as: 𝐾7 u v(𝑰 − 𝒏⨂𝒏) 𝒖x ∙ 𝒗 𝑑𝑠, (3.20) * N;) 7 which is active only on the boundaries ∂Ω3 into Eq. (3.9). The addition of this term is to confine the flow to be in the normal direction at the boundary where pressure is imposed. In the above equation, 𝑰 is the Identity matrix, ⨂ is the cartesian outer product on vectors and 𝐾7 is a penalty factor, which we prescribed a value of 10O . 3.1.3 Time discretization We use the standard first order forward Euler scheme for time-discretization of the inertia term in Eq. (3.9). Because of the high 𝑅𝑒 and 𝑃𝑒, and to avoid any strict time step restriction [56], we use the unconditionally stable Crank Nicolson time-discretization scheme for terms associated with the advection of velocity and temperature, as well as the temperature diffusion and (elliptic) viscous terms. In the non-linear advection, crosswind, SUPG and PSPG stabilization terms, the convective velocities, and constants 𝛾>?-@ , 𝛾BC , 𝛾->-@ are discretized using an explicit time-integration scheme. The convective velocity in the backflow stabilization and LSIC terms, however, are discretized using an implicit time-integration scheme. The complete time discretized weak form with all the stabilization terms is given below: 26 𝒖𝜽Q𝟏 − 𝒖𝜽 〈 , 𝒗〉 + 〈𝒖𝜽 ∙ ∇𝒖𝜽Q𝟏 𝜽Q𝟏 U 𝜽Q𝟏 U 𝒄𝒌 , 𝒗〉 + 〈𝝈v𝒖𝒄𝒌 , 𝑝 x, ∇𝒗〉−〈𝛔v𝒖𝒄𝒌 , 𝑝 x ∙ 𝒏, 𝒗〉𝝏𝛀 ∆𝜏 − 〈𝒇, 𝒗〉 + u 𝛾>?-@ v𝒖𝜽 x v𝑃v𝒖𝜽 , 𝒗x + ∇𝑞x ; 𝒖𝜽Q𝟏 − 𝒖𝜽 ∙k + 𝒖𝜽 ∙ ∇𝒖𝜽Q𝟏 𝜽Q𝟏 U 𝒄𝒌 − ∇ ∙ 𝛔v𝒖𝒄𝒌 , 𝑝 x − 𝒇 𝑑𝑥 m ∆𝜏 + u ‹𝛾BC v𝒖𝜽 x Λv𝒖𝜽 , 𝒖𝜽Q𝟏 xŒ ∶ 𝛁𝒗 𝑑𝑥 ; + u 𝛾&>LB v∇ ∙ 𝒖𝜽Q𝟏 x (∇ ∙ 𝒗) 𝑑𝑥 + 〈10O (𝑰 − 𝒏⨂𝒏) 𝒖𝜽Q𝟏 , 𝒗〉N;* = 𝟎 (3.21) ) ; 〈∇ ∙ 𝒖𝜽Q𝟏 , 𝑞〉 = 𝟎 (3.22) Τ UQ/ − Τ U 1 〈 , 𝑤〉 + 〈𝒖𝜽Q𝟏 ∙ ∇Τ#V , 𝑤〉 + \〈∇Τ#V , ∇𝑤〉 − 〈Τ6 , 𝑤〉);% ] − 〈𝑄, 𝑤〉 ∆𝜏 𝑃𝑒 $ + u 𝛾>?-@ v𝒖𝜽Q𝟏 x 𝑃v𝒖𝜽Q𝟏 , 𝑤x ; 𝜕Τ UQ/ 𝛻 0 Τ#V ∙ k +𝒖 𝜽Q𝟏 ∙ ∇Τ#V − − 𝑄 m 𝑑𝑥 𝜕𝜏 𝑃𝑒 + u ‹𝛾BC vΤ U x Λv𝒖𝜽Q𝟏 , Τ UQ/ xŒ ∶ ∇𝑤 𝑑𝑥 = 0. ; (3.23) 𝒖𝜽,𝟏 Q𝒖𝜽 *.,& Q*. In the above equation, 𝒖𝜽Q𝟏 𝒄𝒌 = 𝟐 and Τ#V = 𝟐 are the Crank-Nicolson velocity and temperature, respectively. Quantities with superscripts 𝜃 + 1 and 𝜃 denote the current and previous timesteps, respectively. The time-discretized equations with the stabilization terms can be solved either using a coupled scheme or a decoupled fractional stepping algorithm. Here, we use the fractional stepping algorithm that is computationally more efficient [57]. In a traditional fractional stepping algorithm (Table 3.1), the computation is split into two separate steps where Eq. (3.21) is decoupled from Eq. (3.22) and solved separately for 𝒖, 𝑝. 27 Table 3.1: Algorithm of the fractional stepping scheme Predictor step : Solve for intermediate velocity (𝒖∗ ) by ignoring the pressure gradient (∇𝑝,) in Eq. (3.9). Projection step : a. Project the pressure field onto a divergence free velocity space (using the intermediate velocity) i.e., 𝛁 ∙ 𝒖∗ = 0. Solve iteratively and minimize pressure residual. b. Perform velocity correction to accommodate for the pressure gradient. In our problem setting, we adopted the Incremental Pressure Correction Scheme (IPCS) which uses pressure gradient from the previous time step at the predictor step, along with semi- implicit formulations for the advection terms as we found that the scheme is robust, accurate and can accommodate larger time steps that reduces the overall simulation time [58]. The two step - IPCS scheme is given as follows: Step 1 - Velocity prediction: 𝒖∗ − 𝒖𝜽 〈 , 𝒗〉 + 〈𝒖𝜽 ∙ ∇𝒖∗𝒄𝒌 , 𝒗〉 + 〈𝝈v𝒖∗𝒄𝒌 , 𝑝U x, ∇𝒗〉−〈𝛔v𝒖∗𝒄𝒌 , 𝑝U x ∙ 𝒏, 𝒗〉𝝏𝛀 − 〈𝒇, 𝒗〉 ∆𝜏 + u 𝛾>?-@ v𝒖𝜽 x 𝑃v𝒖𝜽 , 𝒗x ; 𝒖∗ − 𝒖𝜽 ∙ k + 𝒖𝜽 ∙ ∇𝒖∗𝒄𝒌 − ∇ ∙ 𝛔v𝒖∗𝒄𝒌 , 𝑝U x − 𝒇m 𝑑𝑥 ∆𝜏 + u ‹𝛾BC v𝒖𝜽 x Λv𝒖𝜽 , 𝒖∗ xŒ ∶ ∇𝒗 𝑑𝑥 + u 𝛾&>LB (∇ ∙ 𝒖∗ ) (∇ ∙ 𝒗) 𝑑𝑥 ; ; + 〈10O ∗ (𝑰 − 𝒏⨂𝒏) 𝒖 , 𝒗〉N;* = 𝟎 (3.24) ) Step 2 - Pressure projection & velocity correction: 1 〈∇v𝑝UQ/ − 𝑝U x, ∇𝑞〉 + 〈∇ ∙ 𝒖∗ , 𝑞〉 ∆𝜏 ∗ )∇𝑞 n 𝒖𝜽 − 𝒖𝜽A𝟏 (3.25) + u 𝛾->-@ (𝒖 + 𝒖𝜽 ∙ ∇𝒖𝜽𝒄𝒌 − ∇ ∙ 𝛔v𝒖𝜽𝒄𝒌 , 𝑝U xr 𝑑𝑥 = 0 ; ∆𝜏 28 𝒖𝜽Q𝟏 − 𝒖∗ 〈 , 𝒗〉 + 〈∇v𝑝UQ/ − 𝑝U x, 𝒗 〉 = 0. (3.26) ∆𝜏 𝒖∗ Q 𝒖𝜽 In the above equations, 𝒖∗ is the intermediate velocity and 𝒖∗𝒄𝒌 = 0 is the Crank- Nicolson velocity. We used (P2-P1) discretization for the flow velocity field 𝒖 and pressure field 𝑝, and P2 discretization for the temperature field Τ. Correspondingly, we select 𝜁 = 0.7 in the crosswind stabilization constant 𝛾BC [55]. The system of equations given in Eq. (3.24-2.26) and Eq. (3.23) are solved using the high-level open-source FE library FEniCS [14]. The code can be found in https://github.com/patelte8/vanDANA. 3.2 Numerical details We chose to make use of a modified fractional time-stepping scheme to solve for the pressure- velocity coupling instead of a coupled scheme. The coupled scheme is fully implicit, where Eq’s. (3.9) and (3.10) are solved simultaneously. Compared to the IPCS fractional time-stepping scheme, the coupled scheme is more stable and can accommodate larger time steps. However, the coupled scheme is computationally more expensive because of larger matrix sizes and requires stable pairs of FE discretization that satisfy the inf-sup or Ladyzhenskaya-Babuska-Breezi (LBB) condition [59, 60]. The Taylor-Hood (TH) element is one such pair, in which the velocity field 𝒖 and pressure field 𝑝 are, respectively, discretized using quadratic (P2) and linear (P1) Lagrange-type elements. Equal order (P1-P1) discretization that does not obey the LBB condition can be used only with pressure stabilization PSPG [61, 62] - Eq. (3.15), to circumvent pressure instabilities arising due to the unstable LBB pairs. Here, it is crucial to understand that the SUPG stabilization (Eq. (3.12)) is responsible to stabilize the node-to-node oscillations present due to the high advective nature of the flow and is inherently required regardless of whether TH elements or PSPG stabilization are used. We mainly introduce the PSPG stabilization because we observe that it improves the 29 convergence rate of pressure projection equation (Eq. (3.25)) [62]. On the other hand, the fractional time-stepping scheme theoretically allows one to circumvent the LBB condition and use equal order (P1-P1) discretization without the need for PSPG stabilization. Hence, we can also use P1- P1-P1 discretization for 𝒖, 𝑝 and Τ. Furthermore, we also introduce crosswind stabilization into our weak formulation. The crosswind stabilizing operator is introduced to prevent undershooting/overshooting of the field variables in the presence of sharp gradients, especially in the vicinity of boundary layers [63]. The crosswind stabilization operator adds a non-linear numerical dissipation proportional to the element residual in the direction perpendicular to streamlines in regions where the convective effects are dominant. To better understand the effects of crosswind, LSIC and PSPG stabilizations, we also conducted additional numerical tests (see Appendix). 3.2.1 Implementation in FEniCS FEniCS allows us to write the weak form in the high level UFL format and later assemble it into matrix-vector form to solve it iteratively. However, for a large problem size, using the naive weak form for assembly can be computationally expensive. Hence it is good practice to split the weak form into constant time-independent matrices and vectors wherever possible and perform a pre- assembly before solving for a transient problem. For the discretization in Eq. (3.24), we consider [ the velocity field 𝒖 as piecewise quadratic (P2): 𝒖 = ∑/ 0 𝜓‘Y (𝑥)𝒖𝒊 and the pressure field 𝑝 [ piecewise linear (P1): 𝑝 = ∑/ * 𝜓Y (𝑥)𝑝Y . Here, 𝜓‘ and 𝜓 are quadratic and linear FE basis functions, 𝐾\ and 𝐾7 are the Degrees of Freedom (DOF), and 𝒖𝒊 and 𝑝Y are the nodal values for velocity and pressure respectively. To decompose the weak form, we predefine the following set of constant matrices: 30 [0 [0 2 ∇𝜓 ‘ Y + ∇𝜓‘Y 1 𝐴Y] = “ u 𝜓‘Y ∙ 𝜓‘] ∙ 𝑑𝑥 𝐾Y] = “ u n r : ∇𝜓‘ ∙ 𝑑𝑥 ] ;1 𝑅𝑒 ;1 2 / / [0 [0 𝑂Y] = “ u 𝑰: ∇𝜓‘] ∙ 𝑑𝑥 𝐹Y] = “ u 𝒇 ∙ 𝜓‘] ∙ 𝑑𝑥 / ;1 / ;1 1 2 ∇𝜓‘Y + ∇𝜓‘Y r 𝐵Y] = “ u (𝜓 𝒏) ∙ 𝜓‘] ∙ 𝑑𝑠 𝑀Y] = “ u n 𝒏 ∙ 𝜓‘] ∙ 𝑑𝑠 N;1 𝑅𝑒 N;1 2 [0 𝐿𝑆Y] = “ u 𝛾&>LB v𝛁 ∙ 𝜓‘] x v𝛁 ∙ 𝜓‘] x ∙ 𝑑𝑥 𝐵𝑆Y] = “ 10O u (𝑰 − 𝒏⨂𝒏) 𝜓‘Y ∙ 𝜓‘] ∙ 𝑑𝑠 (3.27) *1 / ;1 N;) Other terms that depend on the evolving solution are assembled separately in the temporal loop. Their corresponding matrices are, v𝒖∗𝒊 + 𝒖𝜽𝒊 x Convective term: 𝒗𝒋 𝑌Y] v𝜓‘, 𝒖𝜽 x 2 SUPG stabilization: 𝒗𝒋 ›𝑈𝑆Y] v𝜓, 𝜓‘, 𝒖𝜽 x 𝒖∗𝒊 − 𝜒] v𝜓, 𝜓‘, 𝒖𝜽 xž Crosswind stabilization: 𝒗𝒋 𝐶𝑊Y] v𝜓‘, 𝒖𝜽 x 𝒖∗𝒊 . (3.28) In the above equation, 𝜒] is a vector of size 𝐾\ . Next, we use all the matrices in Eq. (3.27) and (3.28) to rewrite the weak form: 𝐴Y] ∗ v𝒖∗𝒊 + 𝒖𝜽𝒊 x v𝒖∗𝒊 + 𝒖𝜽𝒊 x 𝒗𝒋 ¡ v𝒖𝒊 − 𝒖𝜽𝒊 x + 𝑌Y] v𝜓‘, 𝒖𝜽 x − 𝑂Y] 𝑝YU + 𝐾Y] + 𝐵Y] 𝑝YU ∆𝜏 2 2 v𝒖∗𝒊 + 𝒖𝜽𝒊 x − 𝑀Y] + 𝐹Y] + 𝑈𝑆Y] v𝜓, 𝜓‘, 𝒖𝜽 x 𝒖∗𝒊 − 𝜒] v𝜓, 𝜓‘, 𝒖𝜽 x 2 + 𝐶𝑊Y] v𝜓‘, 𝒖𝜽 x 𝒖∗𝒊 + 𝐿𝑆Y] 𝒖∗𝒊 + 𝐵𝑆Y] 𝒖∗𝒊 ¢ = 0 (3.29) 31 On further simplification, 𝐴 𝐾 𝑀 𝑌 £¡ Y] + Y] − Y] + Y] ¢ + 𝑈𝑆Y] + 𝐶𝑊Y] + 𝐿𝑆Y] + 𝐵𝑆Y] ¤ 𝒖∗𝒊 ∆𝜏 2 2 2 𝐴Y] 𝐾Y] 𝑀Y] 𝑌 = £¡ − + ¢ − Y] ¤ 𝒖𝜽𝒊 + [𝑂Y] − 𝐵Y] ] 𝑝YU + 𝜒] ∆𝜏 2 2 2 (3.30) An efficient algorithm can be constructed for Eq. (3.30) as follows: • Pre-set sparsity pattern of all the matrices (since mesh is stationary). _23 • Preassemble the following matrix calculations: (𝐿𝑆Y] + 𝐵𝑆Y] ), (𝑂Y] − 𝐵Y] ), ∆+ and _ [23 a23 ‹ ∆+23 − 0 + 0 Œ. Time loop: t < total time: b23 • Assemble matrices: 𝑈𝑆Y] , 𝐶𝑊Y] and − 0 and vector 𝜒] . b23 _ [23 a23 • Evaluate 𝑌Y] ← − 0 + ‹ ∆+23 − 0 + 0 Œ 𝜆] = 𝑌Y] 𝒖𝜽𝒊 + [𝑂Y] − 𝐵Y] ] 𝑝YU + 𝜒] _ • Evaluate 𝑌Y] ← −𝑌Y] + 2 ∆+23 + 𝐶𝑊Y] + 𝑈𝑆Y] + v𝐿𝑆Y] + 𝐵𝑆Y] x • Solve 𝑌Y] 𝒖∗𝒊 = 𝜆] using BICGSTAB with a jacobi preconditioner and appropriate boundary conditions. • Solve Eq. (3.25), (3.26) and (3.23) and repeat. The above optimized algorithm is implemented in FEniCS, and it is highly efficient. The matrices that depend on the evolving solution (Eq. (3.28), 1st step in time loop) are assembled using the naive approach in FEniCS. This is very straight-forward, since one can use the lhs( ) and rhs( ) functions in conjunction with the assemble( ) function call. Note that we use a constant time step ∆𝜏 for our solver. A similar optimization is used for the weak form of energy equation Eq. (3.23) 32 and is not explained here for the sake of simplicity. FEniCS makes it easy to implement the above optimized algorithm in a few lines of high-level python code. The accuracy and convergence of our algorithm has been well tested, and a few benchmark results are presented in the next section. The code’s performance shows good strong scalability up to 84 cores on Michigan State University’s - AMD EPYC 7H12 cluster. 3.3 Benchmarks 3.3.1 2D Lid-driven cavity We solve the standard 2D lid-driven cavity case for 𝑅𝑒 = 1000 and 𝑅𝑒 = 2500. Our domain size is (1,1) and the mesh resolution and time step are 200 x 200 and 0.005 respectively. The mesh is uniform throughout and is composed of triangular elements. The top boundary has a velocity 𝑈 = 1 in the positive x-direction and the other boundaries are no-slip. We use wall boundary condition for pressure at all the four walls. For this setup, we do not solve for the energy equation. Figure 3.1 shows the center line velocity (Ux, Uy) for 𝑅𝑒 = 1000 and 2500 respectively. These are an excellent match with Botella and Peyret [64]. (a.) (b.) Figure 3.1: A comparison of velocity profiles at the center line of the domain with Botella and Peyret [64] for (a) 𝐑𝐞 = 𝟏𝟎𝟎𝟎 (b) 𝐑𝐞 = 𝟐𝟓𝟎𝟎. 33 3.3.2 3D flow over a heated sphere We simulate flow over a 3D heated sphere in a cylindrical domain with dimensions: 𝐿 = 31 and 𝑟 = 10. The 𝑅𝑒 is defined with respect to the sphere diameter and inlet velocity and is set as 𝑅𝑒 = 100, 300 and 1000. The size of the cylindrical domain is considered such that it has negligible effect on the simulation results. We run all the 𝑅𝑒 cases for 𝑃𝑟 = 0.7 and compare our results with Gao et al. [65] and Rodriguez et al. [66]. Note that we do not consider backflow stabilization for this flow problem. The boundary conditions at the inflow and circumferential boundaries are uniform velocity 𝑈 = (1, 0, 0) and at the sphere surface is no-slip. An outflow boundary condition is applied at the outlet. For temperature: at the inflow Τ = 0, at the sphere surface Τ = 1 and for the outlet and circumferential boundaries a Neumann condition is used. We use a non-uniform tetrahedral mesh which is fine at the sphere surface, and it gets coarse as we move downstream and at the far circumferential boundaries. The time step is calculated such that the maximum runtime courant number stays below 0.3. A comparison of our drag coefficient and Nusselt number, with the existing literature is given in Table 3.2. Our results showcase a good match at all the 𝑅𝑒. Table 3.2: Drag coefficient and Nusselt number at Re = 100, 300, 1000 and Pr = 0.7. Cd Nu (Pr = 0.7) Re present Benchmark [65, 66] present Benchmark [66] 100 1.1 1.09 6.91 - 300 0.667 0.659 10.84 10.63 1000 0.468 0.466 17.78 17.4 Next, we also perform a convergence study for the case of 𝑅𝑒 = 1000, and the code is run for a range of mesh sizes. Table 3.3 shows order of convergence and error norm with respect to the benchmark values for all the tetrahedral mesh sizes simulated. Here ℎ is the average cell size of 34 c I the elements close to the sphere and the order of convergence 𝒑, is calculated as 𝑙𝑛 ‹c& Œ /𝑙𝑛 ‹I&Œ. # # We observe that our solver achieves fourth and second order accuracy for velocity and temperature respectively. Table 3.3: Error norm and order of convergence for different mesh sizes. The velocity and temperature are both quadratic (P2). Cd error p Nu error p ℎ 3.75 × 10A0 1.61 × 10A/ --- 1.23 × 10A/ --- 3 × 10A0 6.83 × 10A0 3.84 8.32 × 10A0 1.75 2.25 × 10A0 2.18 × 10A0 3.97 4.83 × 10A0 1.89 1.5 × 10A0 4.30 × 10A< 4.00 2.18 × 10A0 1.96 35 CHAPTER 4 : PATIENT SPECIFIC COMPUTATIONAL MODELING OF CRYOBALLOON ABLATION FOR PULMONARY VEIN ISOLATION 36 4.1 Patient specific geometry and physical setup To investigate CBA, simulations are performed on a patient specific LA geometry reconstructed from the computed tomographic images of a 54-year-old female patient (Patient III in [67]) suffering from paroxysmal atrial fibrillation (Figure 4.1). The model dimensions given in Table 4.1. Table 4.1: LA patient specific model dimensions LSPV LIPV RSPV RIPV MV Area (mm2) 322.5 209.9 188.35 437.6 2061.1 Av. Dia. (mm) 20.26 16.34 15.48 23.6 51.24 Figure 4.1: (a) Schematic of LA geometry with an ellipsoidal CB at the RIPV; the four PV’s and MV serve as inlets and outlet for blood flow in the LA respectively. (b) Close up view of the CB; LAA: LA appendage. Flow rate measurements (at a HR - 63 bpm) of 4 PVs acquired in healthy patients over 1 cardiac cycle [68] are imposed as boundary conditions at the 4 PV ostium serving as inlets to the LA (Figure 4.2a). No-slip boundary condition is applied on the LA walls. Pressure is set as 0 at the )7 outlet of the LA (i.e., MV annulus), whereas for all the other boundaries, we prescribed )𝒏 = 0. A temperature of 37°C is prescribed as Dirichlet boundary condition at the four PV ostia. Neumann 37 boundary condition with zero heat flux (i.e., Τ6 = 0) is applied at the MV outlet. Initial temperature in the problem domain is set as 37°C. The diameter and velocity of the LSPV inlet is taken as the characteristic length and velocity scales, respectively (i.e., 𝐿"# = 20.26 mm and 𝑉"# = 9.3 𝑐𝑚/𝑠). The corresponding control parameters are 𝑅𝑒 ≈ 628 and 𝑃𝑟 = 27 for this problem. Values of the parameters are given in Table 4.2. Table 4.2: Physical properties of blood at 370C. Kinematic Thermal Specific heat Density (kg/m3) viscosity (m2/s) conductivity (W/m0C) (J/kg0C) 1060 3 × 10Ae 0.52 3800 We consider two baseline cases, one without the CB (LA-only case) and the other with a CB to simulate CBA (LA+CB case). In the latter, the cryoballoon (Medtronic CB2) geometry is idealized as an ellipsoid with major/minor diameter = 23/19.85 mm located in the right inferior pulmonary vein (RIPV) [7]. The CB is prescribed a constant minimum temperature of -70°C. At the LA wall, we prescribed a temperature dependent heat source 𝑄 = 𝛽(1 − Τ) to describe heat transfer due to blood perfusion in the LA wall. Here, 𝛽 is the coefficient proportional to the perfusion rate and its quadratic dependence on non-dimensional temperature Τ is given as follows: 6.18 Τ 0 − 7.39 Τ + 2.21 𝑓𝑜𝑟 Τ ≥ 0.725 𝛽=· ». (4.1) 0.1 𝑓𝑜𝑟 Τ < 0.725 The maximum value of 𝛽 at 37°C is taken as 14% (mass fraction of LA) of the rate of coronary perfusion: 0.85 ml/gm/min. At lower temperatures, we assume that perfusion remains constant and reduces to 10% of its original baseline value. While the dependence of vascular perfusion 𝛽 with temperature Τ is not patient-specific, the reduction of β with Τ should hold in general for every patient. Note that the non-dimensional temperature Τ used above is different from the actual 38 temperature Τf#g , which is dimensional. So, Τf#g lies between the temperature ΤhYi = −70°C and Τhfj = 37°C, and is given by: Τf#g = ΤhYi + (Τhfj − ΤhYi ) ∙ Τ. (4.2) (a) (b) Figure 4.2: Prescribed patient-specific blood flow velocity waveforms for the (a) PV inlets and mitral regurgitation (MR) [69]. (b) Perfusion constant 𝜷 as a function of temperature 𝚻𝒂𝒄𝒕 . The behavior of perfusion coefficient 𝛽 given in Eq. (4.1) with actual temperature (Τf#g ) is shown in Figure 4.2b. We note that a similar vascular perfusion-temperature relationship was also employed by Pennes [70] and Getman et. al. [7]. Moreover, we note that Eq. (3.3) reduces to the Pennes heat conduction equation [70] at the LA wall boundary because the conductivity and density of cardiac tissue is very close to that of blood (Table 4.2) and the velocity at the LA wall boundaries is zero here. Therefore, the inclusion of a heat source at the boundary is more physical than prescribing a Neumann boundary condition. We neglect metabolic heat generation of the cardiac tissue and assume it to be negligible when compared to the cryogenic heat flux across the CB. 39 4.2 Validation for Left Atrium hemodynamics To validate the model for LA hemodynamics, a simulation was performed for 25 heart cycles in the LA-only case until the MV flow reaches periodic steady state. In this simulation, the FE tetrahedral mesh has 1.16 million degrees of freedom (DOF) and is graded so that it is fine close to the LA wall. The time step is selected so that the maximum Courant number throughout the simulation runtime is less than 0.12. Note that here, we do not solve for the energy Eq. (3.3). For this case, Figure 4.3a shows the velocity field vectors at the MV plane. Figure 4.3b shows a comparison of our instantaneous MV velocity (averaged over the MV surface area given in Table 4.1) obtained from the LA-only simulation case with that measured from a healthy patient using spectral doppler echocardiography [71]. The simulated MV velocity waveform is close to the measurement, showing an E peak associated with the passive LV filling and an A peak associated with the contraction of the atria (i.e., atrial kick) during LV filling. The model-predicted E/A ratio and deceleration time (DT) are 1.4 and 240 ms, respectively. The greatest difference between the simulated MV velocity waveform and measurement occurs during systole (0.76 – 0.95ms), where the model predicted some (negative) backflow when the MV is closed during systole. This behavior is attributed to the pressure boundary conditions imposed and the assumption of a rigid LA wall geometry. 40 Figure 4.3: Simulation results from LA-only case (a) Velocity vectors at the MV plane. Vectors are scaled in (iv) to illustrate the presence of back/reverse flow. Comparison of (b) simulated MV velocity with measurements from spectral doppler images [71]; and (c) simulated pulmonary vein (PV) pressure waveforms with that measured based on the ultrasound transit time by Smiseth, Otto et. al. [72]. Note: measurements are scaled uniformly based on the prescribed heart rate, velocity, and pressure. Figure 4.3c shows a comparison between the simulated PV pressure waveforms with that measured from a healthy patient based on the ultrasound transit time by Smiseth et. al. [72]. Similar to the MV velocity waveform, the simulated PV pressure waveform is consistent with the measurements during diastole. The largest difference between the model prediction and measurement occurs during systole, where the simulated result shows a slight pressure rise (small peak during systole in Figure 4.3c) that is not found in the measurements. 41 In Table 4.3, we compare the simulated results with the measurement, during diastole and systole using a standard Welch t-Test. We note that the mean for both MV velocity and PV pressure match more closely during diastole as compared to systole. The same is reflected by the higher p- values in case of diastole. Also, despite the deviation of simulated results during systole, the absolute t-value is still less than t-critical which implies a good statistical match. Table 4.3: Comparison of MV velocity and PV pressure for LA-only case during diastole and systole using a standard Welch t-Test. Note: for p-values < 0.05, the comparison is statistically different. Mean Variance t-value (t-critical) p-value MV velocity Diastole present 70.81 1000.79 0.285 (1.98) 0.78 measured [71] 68.83 1098.01 Systole present -10.45 97.59 -2.05 (2.14) 0.06 measured [71] -4.19 8.18 PV pressure Diastole present 0.16 0.09 0.43 (1.99) 0.67 measured [72] 0.13 0.06 Systole Present -0.23 0.013 1.67 (2.16) 0.11 measured [72] -0.31 0.010 4.3 LA+CB results Having validated the model, CBA was simulated in the LA+CB case where the CB was positioned at a fixed position that is close to the lower end of the RIPV ostium (Figure 4.1b). The simulation was performed for 40 heart cycles until a periodic steady state for the power absorbed across the 42 CB is obtained. Because a graded fine mesh is necessary to resolve the gap between the CB and the LV wall, the FE tetrahedral mesh for this case has 1.88 million DOF (Figure 4.1b). The time step for this simulation case is selected such that the maximum Courant number throughout the runtime is less than 0.25. The computing wall time is approximately 46 hours for simulation of the LA+CB case. Figure 4.4a shows the temperature distribution at the RIPV ostium after periodic steady state is reached after about 20 seconds in the LA+CB case. By defining the scar/lesion as the region in the LA wall with temperature less than −20! C [21], the model predicted the lesion area 𝑨𝒍 to be 296.63 𝑚𝑚0 . The maximum lesion width along the flow direction is 6.8 mm. To determine a relationship between the local gap 𝑔 (defined as the minimum distance to the LA wall for any given point on the CB) and the local LA wall temperature 𝑇, we project both the local 𝑔 and 𝑇 onto the CB surface (𝜕Ω Bn ) (Figure 4.4b and c) and plot out their distribution (i.e., 𝑔 – 𝑇 distribution) in Figure 4.4d. The distribution shows that the temperature plateau when the gap 𝑔 is above approximately 6mm. On the other hand, the gap threshold at which the temperature falls below that associated with the lesion formation (i.e., −20! C) is approximately 2.5 mm. By defining the balloon-tissue contact area 𝐴# as the area where the gap 𝑔 < 3mm and an efficiency factor 𝒆 given by the ratio between 𝐴o and 𝐴# , we find that 𝐴𝒄 = 422.8 mm0 and 𝒆 = 0.7 in the LA+CB case. Figure 4.4e shows the rate of energy transfer (power) across the CB with time. At the initiation of CBA, the power is high, then it drops and flattens after about 24s once a periodic steady state is reached. The average value of power absorbed across the CB is approximately 94 W in the LA+CB case. 43 (a) (b) (c) (e) (d) Figure 4.4: Simulation results from LA+CB case (a) Periodic steady state temperature profile at the RIPV ostium. (b) Gap distribution projected onto the CB. (c) Temperature at neighboring RIPV ostium wall (corresponding to the gap) projected onto the CB. (d) Plot for local temperature vs local gap distribution; extracted from Figure 4.4c and 4.4b respectively. (e) Power absorbed across the CB vs time during the cryoablation therapy. Using the LA+CB case as baseline, we then perform simulations with: 1) 12 distinct CB positions where the RIPV ostium is not completely occluded, 2) 4 different PV blood velocities (for healthy patients) at 120%, 90%, 75% and 60% of that given in the Figure 4.2a. This flow rate measurement is applied to all inlets associated with the PVs. 3) Flow reversal at the MV to investigate whether mitral regurgitation (MR) affects the formation of lesion in CBA. To simulate flow reversal at the MV, we use the blood flow rate measurements in a PV of a patient suffering from severe MR [69] (Figure 4.2a). Similar temperature boundary conditions are applied in all the simulations. 44 4.3.1 Effects of CB position on lesion formation Figure 4.5 shows the local 𝑔 – 𝑇 distribution for 3 CB positions out of the 12 simulated CB positions. Here P1, P8 and P11 are notably different in terms of the CB position and therefore, form lesions at different locations at the RIPV ostium. The results show that the temperature fall below that for lesion formation (−20! C) at locations where the gap is less than 2.5 mm at the RIPV ostium irrespective of the CB position. At locations with 𝑔 > 5 − 6 mm, the temperature plateaued at approximately the blood temperature (37! C). CB Gap distribution Temperature Temperature vs gap position (mm) distribution (0C) a. P1 b. P2 Figure 4.5: CB position, P1 to P11 characterized by gap distribution (projected on the CB) and its effect on corresponding temperature distribution at RIPV ostium. The lesion area 𝑨𝒍 , is the total area on either side of the RIPV ostium where 𝚻𝒂𝒄𝒕 < −𝟐𝟎𝟎 𝐂. Also, a plot for local temperature vs gap shows that lesion formation starts if g < 2.5 mm. 45 Figure 4.5 (cont’d) c. P3 d. P4 e. P5 f. P6 46 Figure 4.5 (cont’d) g. P7 h. P8 i. P9 j. P10 47 Figure 4.5 (cont’d) k. P11 Figure 4.6 shows the efficiency factor 𝒆 and average rate of energy transfer (power absorbed in watts) across the CB for all the 12 simulated positions arranged with increasing 𝑨𝒄 . The figure shows that the lesion area 𝑨𝒍 and 𝒆 are increased with increasing 𝑨𝒄 . There is also considerable decrease in power as the efficiency factor increases. These results imply that less energy per unit area is required to maintain the low nadir temperatures of the CB if it is positioned in such a manner that 𝑨𝒄 is high. Figure 4.6: Efficiency factor 𝒆 and average power absorbed across the CB surface (watts) with increasing balloon-tissue contact area 𝑨𝒄 for 12 different CB positions in the RIPV ostium. 48 Figure 4.7 shows the scatter plot between the LA wall temperature Τf#g (see Eq. (4.3)) with the local gap 𝑔 and the balloon-tissue contact area 𝑨𝒄 . The results show that the LA wall temperature is less sensitive to 𝑨𝒄 than 𝑔. The non-linear relationship between Τf#g with 𝑔 and 𝑨𝒄 can be described by the fitted equation: −4268.6 + 1971.47𝑔 − 4.82𝐴# k m ∙ 100 𝑓𝑜𝑟 1.25 𝑚𝑚 ≤ 𝑔 ≤ 9.75 𝑚𝑚 Τf#g = 3666.34𝑔 + 2.44𝐴# (4.3) & 115 𝑚𝑚0 ≤ 𝐴# ≤ 450 𝑚𝑚0 . The coefficient of determination (R2) for the fit is 0.79. Figure 4.7: Non-linear regression fit of the LA wall temperature 𝚻𝒂𝒄𝒕 with the local gap 𝒈 and the balloon-tissue contact area 𝑨𝒄 . 4.3.2 Effects of PV blood velocity and mitral regurgitation on lesion formation Compared to the position of CB, the effects of PV blood velocity and MR are less substantial. Specifically, when the velocity magnitude in the PV is halved, the power absorbed changes by only up to 6%, and the efficiency factor 𝑒 (or lesion area 𝑨𝒍 ) increases anywhere between 5% to 38% depending on the CB position. The changes in 𝑒 and 𝑨𝒍 are larger when 𝑨𝒄 is low. For 49 example, in the case of P3 (Figure 4.8), power drops only by 3.2% and 𝑨𝒍 increases by only 8.5% when the velocity magnitude is halved. In the presence of MR, we observe only a 2% change in lesion size and approximately 3.9% change in power absorbed across the CB. 120 % 60 % P3 (a) (b) Figure 4.8: Temperature distribution and lesion area 𝑨𝒍 at the RIPV ostium for CB position P3 when the PV blood flow velocity magnitude is (a) 120 % (b) 60 % of that given in Figure 4.2a. 4.4 Discussion We investigated 3 main factors that can affect the lesion size for CBA, i.e., catheter balloon position, PV blood velocity and MR. By prescribing PV inlet flow boundary conditions measured in human subjects, our model can predict MV velocity profile that agrees with those measured in humans [71], especially during diastole. The framework also predicts both the E and A peak (E/A ratio = 1.4) and the deceleration time (DT = 240 ms) that are within the normal value range (0.75 < E/A < 1.5 and DT > 140 ms) found in healthy subjects. 4.4.1 Patient-specific model calibration for cryoablation simulations To simulate patient-specific cryotherapy, the governing parameters in Eq. (3.1-3.8) can be made patient-specific by prescribing the 𝑅𝑒, 𝑃𝑟, flow and temperature boundary conditions based on measurements obtained in each patient. Specifically, the blood flow velocity at the PV ostia can be measured using doppler ultrasound or patient specific segmentation can be performed using a 50 4D (3D - time resolved) MRI or CT imaging scan [68]. Moreover, the pressure and temperature in the LA can in principle, be measured using a catheter. In previous studies, Gonzalo et al. [73] and Morales et al. [74] have used in-vivo 4D-CT and 4D-MRI measurements to segment and resolve the LA geometry and flow for their hemodynamic simulations. Using such imaging, it can be very difficult to precisely model areas with less flow for e.g., the LAA, however, to measure flow velocities at major vessels like the PV’s, a set of planes perpendicular to the vessel of interest can be placed to extract the normal projection of the velocity. For medical imaging and segmentation one can utilize tools like MeVisLab [75] or 3D Slicer [76]. Thereafter, these images can also be used to calculate the average diameter of the PV’s and MV. 4.4.2 Factors affecting cryoablation Of the 3 factors that may affect the effectiveness of CBA, namely, CB position, PV inlets velocity and MR, we found that the lesion size is most sensitive to CB position, which in turn, depends on the patient specific PV anatomy. The lesion size depends directly on the size of the balloon-tissue contact area (defined as per area where gap less than 3 mm) and the power absorbed by the CB. Correspondingly, we find that a gap less than 2.5 mm is a threshold to produce a temperature < −20! C for lesion formation. This supports the hypothesis that lesion formation is largely influenced by balloon-tissue contact area 𝐴# and time of freeze [10]. In the LA+CB case that has a high 𝐴# , the average lesion width of 6.8 mm for a 10 cm/s RIPV velocity is comparable with experimental observations in a porcine heart by Parvez et. al. [77]. Their experiments revealed similar lesion widths for a vertical orientation of the CB at 6g contact force and 40 cm/s PV velocity. They noted an increase in lesion size with a decrease in PV velocity. For instance, lesion width increased by 39% from 10.4 mm to 14.5 mm when the PV velocity is 51 halved (from 40 cm/s to 20 cm/s). The same quantitative behavior is observed for other CB positions simulated in Figure 4.5. 4.4.3 Comparison with previous models of cryoablation and LA flow Computational models simulating cryoablation in the LA are limited. Most models [7, 21, 78] are based on idealized geometries as well as boundary conditions and consider only heat transfer across the LA wall to analyze transmural lesion depth. Specifically, a numerical simulation performed by Mussig et al. [78] with the electromagnetic and thermal simulation software CST Darmstadt, ignores the convective effects of blood flow in the PV antrum. The CB geometry in that model, however, is realistic and replicates the model used by Medtronic [7]. For a 23 mm CB set at −50! C, they observed that after 20 sec of ablation the temperature is −50! C and −24! C at a depth of 0.5 and 1 mm respectively. Instead of prescribing a constant CB temperature, Xia et al. [21] applied a constant heat source of −30 W to the CB contact area of a segment of the PV antrum wall. They predicted about 45 s for the temperature to reach −20! C at a depth of 0.5 mm. Also, at depths greater than 2 mm, temperature remains between −20! C to −40! C and does not drop below −40! C, suggesting that the hold time becomes important for lesion/scar formation at that depth. To the best of our knowledge, there is only one computational model [7] that consider blood flow and heat transfer in modeling CBA. That model, however, is based on an idealized geometry of a PV ostium and does not consider flow leakage across the CB. In that study, Getman et al. [7] modelled the temperature at various tissue depths using the non-isothermal flow module in COMSOL. They considered balloon-tissue interactions with the CB set at −60! C. For tissue depths less than 0.5 mm, they predicted TTI (Τf#g < 23! C) and the time required for optimal lesion (Τf#g < −20! C) are 33 and 40 s, respectively. In comparison to these studies, we predicted a shorter time of (~ 24 sec) to steady state lesion formation in the LA+CB case. This difference may 52 be due to (a) A lower nadir CB temperature; (b) Zero thickness of the PV ostium; (c) Difference in boundary conditions at the PV inlets. We also note that accounting for LA wall thickness will increase the power absorbed by the CB and increase the duration required for lesion formation. For power absorbed by the CB during the procedure, the model predicted it to be approximately 165 W for most of the CB positions considered (Figure 4.6). This value is approximately half of that measured in a recent study that uses an ex-vivo calorimetric assessment on a bovine epicardial tissue [79]. That study revealed an ablative power of more than 325 W to generate 2 cm diameter freeze zones in less than 120 s. In the case of excess leakage past the CB, our model also predicts a substantial increase in the ablative power. Similar behavior was observed by Ghosh and McGuire [80]. Data for cooling power, however, is typically not reported for cryoablation probes. Nevertheless, most nitrous oxide probes provide approximately 8.5 W/mm with a range of 50 – 130 W depending on the size/length of the cryocatheter [79, 81]. We also note that for a fixed ablative power, nadir CB temperatures is inversely proportional to the amount of leakage past the CB. 4.5 Conclusion We show that our framework can reproduce measurements of MV velocity and PV pressure waveforms measured in human subjects. The modeling framework predicts that the effectiveness of cryoablation therapy based on the lesion size is very sensitive to the CB position, while PV inlet velocity and MR do not have substantial effects on the lesion size. The framework can be applied in future clinical investigations to optimize cryoablation therapy in patient specific LA-PV geometries. 53 CHAPTER 5 : DEVELOPMENT OF IMMERSED BOUNDARY BASED FICTITIOUS DOMAIN METHOD FOR MODELING OF FLUID STRUCTURE INTERACTIONS WITH HEAT TRANSFER 54 5.1 Methods 5.1.1 Mathematical formulation We use a Eulerian description to solve for the fluid motion and a Lagrangian description to solve for the solid deformation and stresses in the Immersed Boundary (IB) method. We employ the governing continuity, momentum, and energy equations along with kinematic constraints in each domain. A schematic is given in Figure 5.1 and the corresponding boundary value problem (b.v.p.) in dimensional form is given as: For the fluid domain: 𝛺/𝑃(𝑡) r𝒖𝒇 (5.1) 𝜌q rg = 𝛻. 𝝈𝒇 + 𝜌q 𝒇𝒇 𝛻. 𝒖𝑓 = 𝟎 (5.2) r*𝒇 (5.3) 𝜌q 𝐶7q rg = 𝑘q 𝛻 0 Τq + 𝑄𝒇 For the solid domain: 𝑃(𝑡) r𝒖𝒔 (5.4) 𝜌" rg = 𝛻. 𝝈𝒔 + 𝜌" 𝒇𝒔 𝐽𝜌" = 𝜌! (5.5) r*𝒔 (5.6) 𝜌" 𝐶7" rg = 𝑘" 𝛻 0 Τ" + 𝑄𝒔 𝝈𝒔 = −𝑝" 𝑰 + 𝐺(𝑩 − 𝑰) (5.7) Kinematic constraints: 𝑃(𝑡) 𝑑𝒙 = 𝒖𝒔 𝑑𝑡 (5.8) 55 Figure 5.1: A schematic showing the FSI domain and respective boundaries along with the reference and current configuration for the solid. In the above equations we use the constitutive law for a Newtonian fluid (Eq. (3.4)) and a hyper- elastic solid (Eq. (5.6)) and Ω is the overall fluid domain including the fictitious fluid and 𝑃(𝑡) is the solid current configuration. In Eqs. (5.1-5.3), 𝒖𝒇 , 𝑝q , Τ𝒇 , 𝜌q , 𝒇q , 𝝈𝒇 are the fluid velocity, pressure, temperature, density, body force and stress tensor respectively. In Eqs. (5.4-5.7), 𝒖𝒔 , 𝑝" , Τ𝒔 , 𝜌" , 𝒇" , 𝝈𝒔 denotes the solid velocity, pressure, temperature, density, body force and Cauchy stress tensor at the current configuration 𝒙𝜽Q𝟏 . Also, 𝐽 is the determinant of deformation (gradient) tensor 𝐹 which describes the material deformation with respect to the reference )𝒙 configuration and is given as 𝐹 = )𝑿 = 𝛁𝒙. In the solid constitutive law Eq. (5.7), 𝑩 = 𝑭 ∙ 𝑭𝑻 is the finger tensor and G is the shear modulus of the solid material. Additionally, in the energy Eqs. (5.3) and (5.6), 𝑘q , 𝑘" are the thermal conductivity, 𝐶7q , 𝐶7" are the specific heats, 𝑄𝒇, 𝑄𝒔 are the heat sources for the fluid and solid domain respectively. On the boundaries: 𝒖𝒇 = 𝒖𝜞 ; Τ𝒇 = Τ𝜞 𝑜𝑛 𝜞 (5.9) 56 On the FSI interface: )*𝒇 )*𝒔 (5.10) 𝒖𝒇 = 𝒖s ; 𝝈𝒔 . 𝒏 = 𝝈𝒇 . 𝒏 + 𝑭𝒆 ; Τ𝒇 = Τ𝒔 ; 𝑘q )𝒏 = 𝑘" )𝒏 𝑜𝑛 𝜕𝑃 Here, 𝒏 denotes the vector normal to the FSI interface and 𝑭𝒆 is the non-hydrodynamic forces at the fluid-solid interface, like external contact forces. Note that in Eq. (5.5), 𝜌! is the solid density in the reference configuration and if the solid is incompressible i.e., 𝐽 = 1, we also need to compute {6 oi| the solid pressure 𝑝" . For a compressible hyper-elastic solid, however, we can replace 𝑝" by | , where 𝜆! is related to the bulk compressibility of the solid, and compute only the solid deformation. In our model, we implement both the incompressible and compressible formulations for the solid constitutive law. We also note that at high values of 𝜆! ~100 𝐺, the solid becomes close to incompressible (𝐽~1). 5.1.2 Finite element discretization The FEM is used for spatial discretization of the governing equations in both the fluid and solid domain. Let 𝒗, ω be the test functions for the momentum, q be the test function for the continuity, and γ, 𝛾! be the test functions for the energy equations - for the fluid and solid domains respectively; then the Galerkin weak formulation is given as: Momentum: 𝑑𝒖𝒇 𝑑𝒖𝒔 〈𝜌q − 𝜌q 𝒇𝒇 , 𝒗〉;/- + 〈𝝈𝒇 , 𝛻𝒗〉;/- + 〈𝜌" − 𝜌" 𝒇𝒔 , 𝝎〉- 𝑑𝑡 𝑑𝑡 (5.11) + 〈𝝈𝒔 , 𝛻𝝎〉- = u 𝑭𝒆 ∙ 𝝎 𝒅𝒙 )- Continuity: 〈𝑞 , ∇ ∙ 𝒖𝒇 〉;/- = 0 (5.12) Energy: 𝑑Τq 𝑑Τ" 〈𝜌q 𝐶7q − 𝑄q , 𝛾〉;/- + 〈𝑘q 𝛻Τq , 𝛻𝛾〉;/- + 〈𝜌" 𝐶" − 𝑄" , 𝛾" 〉- 𝑑𝑡 𝑑𝑡 (5.13) + 〈𝑘" 𝛻Τ" , 𝛻𝛾" 〉- = 0 57 where `𝒖𝒇 , 𝒖𝒔 , Τ𝒇 , Τ𝒔 , 𝑝" , 𝑝q ∈ 𝐻/ (Ω)a and `𝒗, 𝑞, 𝝎, 𝛾, 𝛾𝑠 ∈ 𝐻/ (Ω)a. In the above equations, 〈∙ ,∙ 〉 denotes the Euclidean inner product of two functions in their respective domains. Fictitious Domain (FD) / Distributed Lagrange Multiplier (DLM) formulation Next, using the fictitious domain (FD) formulation (initially proposed by Yu [33]), we extend the fluid domain into the immersed solid and assume that the interior of the solid is filled with a fictitious fluid (see Figure 5.2). This fictitious fluid is constrained to move with the same velocity as the solid by enforcing a distributed Lagrange multiplier 𝜆 which acts as a pseudo volumetric force. For coupling between both the domains we use the Dirac-delta function kernels to calculate the pseudo forces. Figure 5.2: Fictitious domain formulation where the fluid is extended into the solid domain and 𝒖𝒇 = 𝒖𝒔 in the overlapping region. Using the FD/DLM formulation, we get, Fluid: 𝑑𝒖𝒇 (5.14) 〈𝜌q − 𝜌q 𝒇𝒇 , 𝒗 〉; + 〈𝝈𝒇 , 𝛻𝒗〉; + 〈𝑞 , ∇ ∙ 𝒖𝒇 〉; = 〈𝝀 , 𝒗〉- 𝑑𝑡 Fluid temperature: 𝜕Τ𝒇 1 〈 − 𝒖𝒇 ∙ ∇Τ𝒇 − 𝑄q , 𝛾〉; + 〈∇Τ𝒇 , ∇𝛾〉; = 〈𝜆 1 , 𝛾〉- (5.15) 𝜕𝑡 𝑅𝑒𝑃𝑟 Solid: 𝑑𝒖𝒔 〈v𝜌" − 𝜌q x − (𝜌" 𝒇𝒔 − 𝜌q 𝒇𝒇 ), 𝝎〉- + 〈𝝈𝒔 − 𝝈𝒇 , 𝛻𝝎〉- 𝑑𝑡 (5.16) 𝒆 = u 𝑭 ∙ 𝝎 𝒅𝒙 − 〈𝝀 , 𝝎〉- )- Lagrange multiplier constraint: 〈𝒖𝒇 − 𝒖𝒔 , 𝝃〉- = 0 (5.17) 58 Solid temperature Lagrange r*𝒔 〈(𝜌𝑟𝐶7~ − 1) − (𝑄" − 𝑄q ), 𝛾" 〉- + rg multiplier: V7 A/ (5.18) 〈 ∇Τ" , ∇𝛾" 〉- = −〈𝜆 1 , 𝛾" 〉- -. where 𝜉 and 𝛾" are the test functions on the Lagrange multiplier space. Here, the distributed Lagrange multiplier (DLM), is introduced in the FD formulation to separate the fluid and solid equations from Eq. (5.11) and relax the velocity constraint. Note that, 𝜆 1 is the temperature-based Lagrange multiplier, which acts as a pseudo heat source [82], besides 𝜆 is the momentum based Langrange multiplier which acts as a pseudo body force in the momentum equation. To impose incompressibility constraint for the solid, we solve for solid pressure and explicitly add the following to Eq. (5.16): u (𝐽 − 1). 𝜁𝑑𝑥 = 0 (5.19) - It is worth pointing out that considering the system as a whole, in the pair of Eqs. (5.14, 5.16) and Eqs. (5.15, 5.18), the fictitious Lagrange multiplier terms cancel out. 5.1.3 Numerical scheme To solve for the fluid velocity-pressure coupling in Eq. (5.14), we use the optimized algorithm for incremental pressure correction scheme (IPCS) described in section 3.1.3 and 3.2.1 and also introduce stabilization schemes, namely: SUPG (Eq. (3.12)), PSPG (Eq. (3.15)) and crosswind stabilization (Eq. (3.16)). Detailed information on the need for stabilization schemes is provided in section 3.2. Further, to non-dimensionalize the governing equations, we introduce the following characteristic length 𝐿"# and velocity 𝑉"# scales, which changes the variable unknowns to ××× 𝒖𝒇 = 𝒖𝒇 𝒖 78 7! 𝒙 ×××𝒔 = 𝒔 , ××× , 𝒖 𝑝q = • # , 𝑝Ø" = • # and 𝒙 Ø = . The subsequent non-dimensional solid $!" $!" 8 $!" 8 $!" 𝑳 𝒔𝒄 • B V { @ parameters are 𝜌~ = • ! , 𝐶7~ = B*! , 𝑘~ = V ! , ××× 𝜆! = • $6 # and 𝐺̅ = • # . To summarize, we have 8 *8 8 8 !" 8 $!" 59 $!" &!" $!" a set of four independent control parameters for our IB-FSI formulation: 𝑅𝑒 = • , 𝐹𝑟 = √ƒ& , !" • $!" # † 𝑃𝑟 = ( , 𝐸𝑐 = B , where 𝜈 is the kinematic viscosity given as 𝜈 = • , 𝑔 is the acceleration *8 (*: A*; ) 8 due to gravity, 𝛼 is the thermal diffusivity and we introduce 𝐸𝑐, that accounts for heat generation due to viscous dissipation of the fluid. Another important control parameter that appears in the energy equation is the Peclet number; 𝑃𝑒 = 𝑅𝑒 ∙ 𝑃𝑟. This completes the dimensionless FD/DLM formulation with FE stabilizations. The numerical algorithm in order of solving the equations and its corresponding numerical discretization is given below: 𝒖𝒇 ∗ − 𝒖𝒇 𝜽 𝟑𝒖𝒇 𝜽 − 𝒖𝒇 𝜽A𝟏 Fluid: 〈 + ∙ 𝛁𝒖∗𝒄𝒌 , 𝒗 〉; ∆𝝉 𝟐 𝛁𝒖∗𝒄𝒌 + (𝛁𝒖∗𝒄𝒌 )1 + 〈−𝑝q U 𝑰 + , ∇𝒗 〉; 𝑅𝑒 Ý 𝒈 (5.20) + 〈𝑝q U ∙ 𝒏, 𝒗 〉N; − 〈 0 , 𝒗 〉; 𝐹𝑟 + 〈γ>?-@ (𝒖𝒇 𝜽 ) 𝑃(𝒖𝒇 𝜽 , 𝒗), 𝑅U 〉; + 〈γBC (𝒖𝒇 𝜽 ) Λ(𝒖𝒇 𝜽 , 𝒖𝒇 ∗ ), ∇𝒗 〉; = 〈𝝀𝜽 , 𝒗 〉- 1 〈∇v𝑝q UQ/ − 𝑝q U x, ∇𝑞〉𝛀 + 〈∇ ∙ 𝒖𝒇 ∗ , 𝑞〉𝛀 + 〈𝛾->-@ v𝒖𝒇 ∗ x∇𝑞, 𝑅U 〉𝛀 = 0 ∆𝜏 (5.21) 𝒖𝒇 𝜽Q𝟏 − 𝒖𝒇 ∗ 〈 , 𝒗〉 + 〈∇v𝑝q UQ/ − 𝑝q U x, 𝒗 〉𝛀 = 0. (5.22) ∆𝜏 Fluid Τ𝒇 UQ/ − Τ𝒇 U 1 temperature: 〈 , 𝛾〉𝛀 + 〈𝒖𝒇 𝜽Q𝟏 ∙ ∇Τ#V , 𝛾〉𝛀 + 〈∇Τ#V , ∇𝛾〉𝛀 ∆𝜏 𝑃𝑒 2𝐸𝑐 −〈 〈𝑬(𝒖𝒇 ), ∇𝒖𝒇 〉𝜽Q𝟏 , 𝛾〉; 𝑅𝑒 (5.23) + 〈𝛾>?-@ v𝒖𝒇 𝜽Q𝟏 x 𝑃v𝒖𝒇 𝜽Q𝟏 , 𝛾x, 𝑅U 〉; + 〈𝛾BC vΤ𝒇 U x Λv𝒖𝒇 𝜽Q𝟏 , Τ𝒇 UQ/ x, ∇𝛾〉; = 〈𝜆 1 U , 𝛾〉- 60 ∆𝒙𝒔 𝜽Q𝟏 UQ/ Solid: 〈𝜌𝑟 𝟐 , 𝝎 〉-6 + 〈(∇! 𝝎)1 , −𝑝" 𝑭A𝟏 + 𝐺(𝑭𝑻 − 𝑭A𝟏 )〉-6 ∆𝝉 𝒖𝒇 𝜽Q𝟏 ∆𝒙𝒔 𝜽 =〈 + (𝜌𝑟 − 1) , 𝝎 〉-6 + 〈𝐽UQ/ − 1, 𝜁〉-6 ∆𝝉 ∆𝝉𝟐 (5.24) Ý 𝒈 + 〈(𝜌𝑟 − 1) 0 , 𝝎 〉-6 − 〈𝝀𝜽 , 𝝎〉-6 𝐹𝑟 DLM 𝜽Q𝟏 ∆𝒙𝒔 𝜽Q𝟏 Constraint: 〈𝒖𝒔 − , 𝝃 〉- = 𝟎 ∆𝝉 (5.25) 𝒖𝒔 𝜽Q𝟏 − 𝒖𝒇 𝜽Q𝟏 〈 , 𝝎𝒔 〉- = 〈𝜆UQ/ − 𝜆U , 𝝎𝒔 〉- ∆𝝉 (5.26) Solid Τ𝒔 UQ/ − Τ𝒔 U 2𝐸𝑐 temperature 〈𝜌𝑟𝐶7~ − 1 , 𝛾" 〉- + 〈 〈𝑬(𝒖𝒇 ), ∇𝒖𝒇 〉𝜽Q𝟏 , 𝛾" 〉- ∆𝜏 𝑅𝑒 (5.27) Lagrange 𝑘~ − 1 +〈 ∇Τ",#V , ∇𝛾" 〉- = −〈𝜆 1 UQ/ , 𝛾" 〉- multiplier: 𝑃𝑒 In the above equations, 𝒖∗𝒄𝒌 and Τ#V is the Crank-Nicolson velocity and temperature given as 𝛁𝒖∗𝒄𝒌 = 0.5v𝒖𝒇 ∗ + 𝒖𝒇 𝜽 x and Τ#V = 0.5 (Τ𝒇 UQ/ + Τ𝒇 U ) and ∆𝒙𝒔 is the incremental solid displacement. Also, Eq. (5.22) is solved in the reference configuration 𝑃! and the gradients in both the current and reference configuration are related as 𝛁 = 𝛁! (𝑭A/ ). Note that in the above equations we drop the bar on top of the non-dimensional variables and neglect any external hydrodynamic forces at the fluid-structure interface i.e., 𝑭𝒆 = 𝟎 for the sake of simplicity. Also, one needs to keep in mind that Eq. (5.22) applies to a perfectly incompressible solid and we compute both the solid pressure 𝑝" and displacement 𝒙𝒔 . However, if the solid is compressible, we {6 oi| replace the solid pressure by | and remove 〈𝐽 − 1, 𝜁〉-6 term from Eq. (5.24) and solve only for the incremental solid displacement ∆𝒙𝒔 . To solve Eq. (5.20) and Eq. (5.23), we use the stabilized biconjugate gradient iterative method with hypre-amg preconditioners at an absolute tolerance value of 10-8. To solve the solid problem 61 Eq. (5.24), we use the non-linear Newton’s iteration for a relative tolerance value of 10−6. The ∆𝒙𝒔 from the previous time step is provided as initial guess to the Newton solver and we observe that it takes about 3-4 iterations to converge. 5.1.4 Interpolation functions To interpolate fictitious forces, heat source, velocity, or temperature from the Eulerian mesh to the Lagrangian mesh or vice-versa, we use the 4-point conservative piecewise-delta function [83] given below: 𝑢‡ (𝑋) = ∫ˆ 𝑢‡ (𝑥)𝛿(𝑥 − 𝑋) 𝑑𝑥 ; 𝜆(𝑥) = ∫- 𝜆(𝑋)𝛿(𝑥 − 𝑋) 𝑑𝑥 ; (5.28) Τ𝒔 (𝑋) = ∫ˆ Τ𝒔 (𝑥)𝛿(𝑥 − 𝑋) 𝑑𝑥 ; 𝜆 1 (𝑥) = ∫- 𝜆(𝑋)𝛿(𝑥 − 𝑋) 𝑑𝑥 1 𝑥−𝑋 𝑦−𝑌 𝑧−𝑍 𝛿(𝒙 − 𝑿) = < ∅k m∅k m∅k m ℎ ℎ ℎ ℎ (5.29) 1 ⎧ (3 − 2|𝑟| − ê1 + 4|𝑟| − 4|𝑟|0 ) |𝑟| ≤ 1 ⎫ ⎪ 8 ⎪ ∅(𝑟) = 1 1 ≤ |𝑟| ≤ 2 (5.30) ⎨ (5 − 2|𝑟| − ê−7 + 12|𝑟| − 4|𝑟|0 ) 2 ≤ |𝑟| ⎬ ⎪8 ⎪ ⎩ 0 ⎭ In the IB-FSI solver, the delta functions essentially determine the order of accuracy of the IB method. Due to the inherent nature of the delta functions, it is obligatory that the fluid grid surrounding the solid is structured and for good accuracy, we recommend that the solid mesh size stays between 0.8h - 1.5h at all times where h is the uniform fluid grid size [33, 84]. Here, note that the solid grid can be unstructured and non-uniform. Also, the delta functions utilize piecewise functions ∅ with different support areas and continuity properties which directly dictate their behavior in terms of numerical stability and accuracy [85]. We use the 4-point piecewise function (Eq. (5.30)) since we observe that it provides a good balance between stability and accuracy. 62 5.2 vanDANA solver Our code is open-source and available to use in the form a GitHub repository: https://github.com/patelte8/vanDANA. vanDANA is an efficient FEM Immersed Boundary (IB) based Flow-thermal FSI solver utilizing the FEniCS library and it is a python package with main executable (vanDANA.py) and three submodules (common, utilities, user_inputs). The directory tree for vanDANA code is given below in Figure 5.3. Figure 5.3: Directory tree for vanDANA IB-FSI code. 1. The common module is the heart of the solver and mainly comprises of different variational problems: flow_variational_problem.py, flow_temperature_variational_problem.py, solid_variational_problem.py, lagrange_variational_problem.py. Each of these files contain separate class definitions and subroutines for pre-assembly, run-time assembly, solvers, and post-processing for their respective variational problems. We also include some basic post- processing subroutines for calculation of drag, lift, Nusselt number, Jacobian, and vorticity. 63 Hence, the user need not require making any changes to the common module unless one needs to implement a custom user-defined subroutine. 2. The utilities module deals with the handling of background functions for the code and comprises of the files: utils.py, read.py and write.py. Almost all the functionality for read, write, restart, MPI, counters, memory usage and timing of different modules is initialized and carried out using this module. 3. Next, the user_inputs module provides input to the code mainly in terms of control parameters, time step, meshes, initial and boundary conditions for any given physical problem. Here, we allow user control from two main files namely: user_parameters.py and boundary_initial_conditions.py. A sample example of user parameters.py is given below in Figure 5.4. All parameters are grouped in the form of python dictionaries which makes them easy to append during runtime. restart = False # Restart parameter # Physics of the problem # --------------------------------------------------------------------- problem_physics = dict( solve_temperature = True, solve_FSI = True, compressible_solid = True, viscous_dissipation = False, body_force = False, ) def f_dir(dim): # Body force direction n = -1*tensors.unit_vector(1, dim) return n interpolation_fx = 'phi4' # Delta-function interpolation # FEM stabilization and constants # --------------------------------------------------------------------- stabilization_parameters = dict( Figure 5.4: All control parameters are provided as input to the vanDANA code from user_parameters.py. 64 Figure 5.4 (cont’d) # Navier-stokes SUPG_NS = False, # explicit PSPG_NS = False, # explicit crosswind_NS = False, # implicit # Energy-equation SUPG_HT = False, # explicit crosswind_HT = False # implicit ) alpha = Constant(0.85) # SUPG/PSPG stabilization constant C_cw = Constant(0.7) # Crosswind stabilization constant # Physical parameters # --------------------------------------------------------------------- physical_parameters = dict( g = 9.81, # Gravity (m/s2) # Fluid rho_f = 1, # Density (kg/m3) nu = 1, # Dynamic viscosity (kg/m.s) Spht_f = 1, # Specific heat (J/kg.C) K_f = 1, # Thermal conductivity (W/m.C) # Solid rho_s = 10, # Density (kg/m3) Sm = 0, # Shear modulus (N/m2) Ld = 0, # Compressibility (N/m2) Spht_s = 0.11, # Specific heat (J/kg.C) K_s = 1.2 # Thermal conductivity (W/m.C) ) def calc_non_dimensional_solid_properties(g, rho_f, nu, Spht_f, K_f, rho_s, Sm, Ld, Spht_s, K_s, Lsc, Vsc, T0, Tm, Tsc): rho = rho_s/rho_f Spht = Spht_s/Spht_f K = K_s/K_f Ld = 2000 # Ld/(rho_f*Vsc*Vsc) Sm = 500 # Sm/(rho_f*Vsc*Vsc) return rho, Spht, K, Ld, Sm # Characteristic scales # --------------------------------------------------------------------- characteristic_scales = dict( Lsc = 1, # m Vsc = 1, # m/s 65 Figure 5.4 (cont’d) T0 = -1*52, # lower_temp (C) Tm = 37 # higher_temp (c) ) # Temporal control # --------------------------------------------------------------------- time_control = dict( C_no = 0.35, # Maximum possible Courant number dt = 0.0025, # Time-step T = 100, # Total runtime adjustable_timestep = True ) # FEM degree of variables # --------------------------------------------------------------------- fem_degree = dict( velocity_degree = 2, pressure_degree = 1, temperature_degree = 2, displacement_degree = 2, lagrange_degree = 1 ) # Non-dimensional numbers # --------------------------------------------------------------------- def calc_non_dimensional_numbers(g, rho_f, nu, Spht_f, K_f, rho_s, Sm, Ld, Spht_s, K_s, Lsc, Vsc, T0, Tm, Tsc): Re = rho_f*(Vsc*Lsc)/nu Pr = Spht_f*nu/K_f Ec = (Vsc*Vsc)/(Spht_f*(Tm-T0)) Fr = Vsc/sqrt(g*Lsc) return Re, Pr, Ec, Fr # Enter "True" if you want to post-process data # --------------------------------------------------------------------- post_process = True # File printing / solid-remeshing control # --------------------------------------------------------------------- print_control = dict( a = 40, # for printing variables and restart files b = 20, # for post processing data c = 20, # for simulation_wall_time text file d = 5, # for remeshing solid current-configuration mesh e = 20 # for runtime_tsp_courant_no_stats text file ) 66 In user_parameters.py, we have provided the following important features, namely: • The restart parameter, if "True" will restart the simulation from the last saved time step. • The order of finite-element basis functions for all variables can be specified using the fem_degree dictionary. • In the time_control dictionary if adjustable_timestep is set as True, the code will automatically calculate the time step during run-time by limiting the maximum Courant number to time_control['C_no']. Otherwise, if adjustable_timestep is set as False, it will use time_control['dt'] as a constant time step. • The calc_stream_function flag if True, enables calculation of vorticity and stream function only for a 2D problem setup. This is because the steam function is defined only in 2D, using ∇0 𝜓 = −𝜔, where 𝜔 is the vorticity. • The effect of stabilization schemes can also be tuned using stabilization constants: alpha and C_cw. To execute the code, vanDANA.py is the main executable file that pulls information from all the three modules and runs the main function vanDANA_solver (args). This sets the entire workflow defined in Eqs. (5.20-5.27). Overall, vanDANA code is aimed to be user-friendly and hence to gain more user-control for ease of testing and restarting, we use the argparse library [86] to add keywords from the command line, for e.g.: mpirun.mpich -n 64 python3 vanDANA.py - restart=True -T=20. 5.3 Benchmarks In this section we validate four different FSI benchmarks. In all these benchmark cases, we used second order (quadratic) elements for velocity and solid displacement and linear elements for pressure and Lagrange multiplier. 67 5.3.1 Transverse flow over a flexible beam This problem was introduced by Baaijens [34] and later also studied by Yu [33]. A thin flexible leaflet is located transverse to the flow direction at the center of the 2D channel as shown in Figure 5.5. For velocity boundary conditions, a no-slip at Γ0 , free slip at ΓO and a parabolic time varying 0‰ velocity profile given as 𝑈 = 1.5𝑦(2 − 𝑦)sin ( 1 𝑡) is provided at Γ/ and Γ< , where the origin is defined as the bottom left corner and the total runtime is 𝑇 = 10. For pressure we provide wall boundary conditions on all the sides. The channel length is 10 and height 𝐻 = 1. The density ratio is set as 𝜌~ = 1 and the solid is considered incompressible with stiffness 𝐺 = 10< . Note that since this is an incompressible leaflet i.e., 𝐽 = 1, we also solve for the solid pressure. Figure 5.5: Schematic geometry for transverse flow over a flexible plate. The leaflet is held stationary at its bottom end and its length and height are 0.8 and 0.0212 respectively. All body forces are neglected. The Reynolds number considered here is 100 and the fluid mesh is uniform with a mesh size of ℎq = 1/128. The solid mesh size is ℎ" = 1/100 along the length of the leaflet and it has 3 grid cells along its thickness. The time step is set as 0.005. Figure 5.6 shows our leaflet position in the channel at 8 different time instants over one cycle of oscillation. These results depict and an excellent qualitative match with those of Yu [33]. 68 Figure 5.6: Transverse flow over a flexible leaflet: Re 100, 𝝆𝒓 = 𝟏, 𝑮 = 𝟏𝟎𝟑 . A comparison of present results with Yu [33] for leaflet positions in the channel over one periodic cycle. 5.3.2 Turek and Hron (FSI2) Next, we consider the classical FSI2 benchmark by Turek and Hron [87] - flow over an elastic beam attached behind a rigid cylinder. The schematic is given in Figure 5.7. Here 𝐷, the diameter of the rigid cylinder is defined as the characteristic length scale 𝐿"# and is chosen as 𝐷 = 1. In Figure 5.7, other model dimensions are given as 𝐻 = 4.1, 𝐿 = 11 and 𝑙 = 3.5, ℎ = 0.2. The center of the rigid cylinder is (2,2) and at time 𝑡 = 0, the flexible beam is symmetric with respect to the cylinder center. The boundary conditions for the top and bottom boundary are no-slip, left e?6 Œ(•AŒ) boundary is a constant inlet velocity profile 𝑈 = •# and the right boundary is outflow. Here, 𝑈! = 1 is defined as the characteristic velocity scale 𝑉"# . For pressure, we apply wall boundary 69 condition at all the boundaries, except the at the right boundary the pressure is set as 0. The 𝑅𝑒 is defined with respect to 𝐷 and 𝑈! , and is set as 100. Figure 5.7: A schematic diagram for the FSI2 benchmark [87]: Flow past a flexible beam (compressible) attached behind a rigid cylinder. For the properties of the flexible beam, the non-dimensional shear modulus is 𝐺 = 500, density ratio is 𝜌~ = 10, and the solid is considered compressible with non-dimensional compressibility as 𝜆! = 2000. The fluid mesh is 575 x 215, and the flexible beam has a resolution of 175 x 10. The time step is set such that the average courant number during runtime is around 0.35 and the total non-dimensional run time is T = 100 sec. The total DOF’s for the case setup is 1.15 million. Figure 5.8: Comparison of present y-displacement of trailing tip of the flexible beam with Turek and Hron [87]. Our results showcase an excellent match for the y-displacement of the trailing tip of the flexible beam (Figure 5.8), strouhal number and average drag coefficient (Table 5.1) with results from the 70 Ž literature. The drag coefficient is calculated as 𝐶3 = !.=• 12 million. Second, while vanDANA using FEniCS, both the fluid and solid mesh are split into N-MPI partitions each. However, in such IB-FSI simulations in most cases the solid mesh is smaller in volume than the fluid mesh and such a partitioning strategy retards the scalability of the solver as a whole. Third, the present Lagrange multiplier based heat transfer algorithm is unstable at higher values of 𝜌~ 𝐶7~ − 1 and needs implicit treatment of the unsteady term in Eq. (5.25) [82]. Fourth, solving for physical problem cases with large values of 𝑘~ , can be unstable. However, this can be tackled by using a smaller time step. Last, the IB algorithm in conjunction with the delta-function interpolation kernel requires the fluid cartesian mesh to be uniform in all three directions at least in the region overlapping the solid. In case of an internal flow problem, this uses a larger fluid mesh volume. Here, one needs to judiciously provide physical flow boundary conditions since it can instigate solver divergence. 6.2 Concluding remarks This thesis can mainly be summarized in two parts: • Patient specific CBA modeling: We have developed a computational stabilized FE framework to simulate the dynamic thermal-fluid coupling in a realistic LA geometry during cryoablation therapy. We show that the framework can reproduce measurements of MV velocity and PV pressure waveforms measured in healthy human subjects. The modeling framework predicts that the effectiveness of cryoablation therapy based on the lesion size is very sensitive to the CB position, while PV inlet velocity and MR do not have substantial effects on the lesion size. 80 The framework can be applied in future clinical investigations to optimize cryoablation therapy in patient specific LA-PV geometries. • Development of vanDANA IB-FSI solver: We have developed a fully coupled FE based immersed boundary FSI code using the Distributed Lagrange Multiplier based fictitious domain method and have extended it to deal with heat transfer. We have used FEniCS, and our code is open source, documented, user-friendly and is available to use in the form of a GitHub repository: vanDANA. The vanDANA code is benchmarked against complex FSI problem cases and its scalability has been well tested on HPC. In future, our computational framework can be used to simulate fully coupled FSI physics (including heat transfer) for complex cardiovascular flows. The future direction for this project should be focused on addressing the limitations and improving the computational speed and scalability of the vanDANA solver. Other computational models (or subroutines) that incorporate turbulence modeling, active strain, contact mechanics, and non- Newtonian models for blood flow can be implemented into the vanDANA solver. This will be of use to capture complex flow physics and investigate cardiovascular diseases in patient specific geometries more accurately. 81 BIBLIOGRAPHY [1] R. E. Klabunde, Cardiovascular physiology concepts, Lippincott Williams & Wilkins, 2011. [2] "Assessment of Mitral Valve Function | Radiology Key.," July 3, 2022. [Online]. Available: https://radiologykey.com/assessment-of-mitral-valve-function. [3] V. Markides and R. J. Schilling, "Atrial fibrillation: classification, pathophysiology, mechanisms and drug treatment," Heart 89, no. 8 , pp. 939-943, 2003 . [4] Z. Nesheiwat, A. Goyal and M. Jagtap, Atrial Fibrillation (Nursing), Treasure Island (FL): StatPearls Publishing. PMID: 30252328, 2022 Jul 31. [5] M. Sucu, V. Davutoglu and O. Ozer, "Electrical cardioversion," Annals of Saudi medicine 29, no. 3 , pp. 201-206, 2009. [6] O. M. Wazni, G. Dandamudi, N. Sood, R. Hoyt, J. Tyler, S. Durrani and M. Niebauer, "Cryoballoon ablation as initial therapy for atrial fibrillation," New England Journal of Medicine 384, no. 4 , pp. 316-324, 2021 . [7] M. Getman, E. Wissner, R. Ranjan and J. Lalonde, "Relationship between time-to-isolation and freeze duration: Computational modeling of dosing for Arctic Front Advance and Arctic Front Advance Pro cryoballoons," J Cardiovasc Electrophysiol. 30, p. 2274–2282, 2019. [8] M. Yacoub and R. Sheppard, "Cryoballoon Pulmonary Vein Catheter Ablation of Atrial Fibrillation," StatPearls, 2021. [9] N. Patel, J. Maradey and P. Bhave, "Atrial Fibrillation Ablation: Indications and Techniques," Curr Treat Options Cardio Med. 21, p. 43, 1936. [10] W. Su, R. Kowal, M. Kowalski, A. Metzner and J. Svinarich, "Best practice guide for cryoballoon ablation in atrial fibrillation: The compilation experience of more than 3000 procedures," Heart Rhythm. 12, p. 1658–1666, 2015. [11] T. Patel and A. Lakdawala, "A dual grid, dual level set based cut cell immersed boundary approach for simulation of multi-phase flow," Chemical Engineering Science. 177, p. 180– 194, 2018. [12] A. Athani, N. N. N. Ghazali, I. A. Badruddin, S. Kamangar, A. E. Anqi and A. Algahtani, "Investigation of two-way fluid-structure interaction of blood flow in a patient-specific left coronary artery," Bio-Medical Materials and Engineering 33 no. 1 , pp. 13-30, 2022. [13] Y. I. Cho, D. J. Cho and R. S. Rosenson, "Endothelial shear stress and blood viscosity in peripheral arterial disease," Current atherosclerosis reports 16, pp. 1-10, 2014. 82 [14] A. Logg, K. Mardal and G. Wells, Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Manual, 2012. [15] A. Henderson, ParaView Guide, A Parallel Visualization Application, Kitware Inc., 2007. [16] V. Parikh and M. Kowalski, "Comparison of phrenic nerve injury during atrial fibrillation ablation between different modalities, pathophysiology and management," Journal of atrial fibrillation 8, no. 4 , 2015. [17] D. L. Packer, R. C. Kowal, K. R. Wheelan, J. M. Irwin, J. Champagne, P. G. Guerra and M. Dubuc, "Cryoballoon ablation of pulmonary veins for paroxysmal atrial fibrillation: first results of the North American Arctic Front (STOP AF) pivotal trial," Journal of the American College of Cardiology 61, no. 16, pp. 1713-1723, 2013. [18] R. Casado-Arroyo, G.-B. Chierchia, G. Conte, M. Levinstein, J. Sieira, M. Rodriguez- Mañero and G. D. Giovanni, "Phrenic nerve paralysis during cryoballoon ablation for atrial fibrillation: a comparison between the first-comparison between the first-and second- generation balloon," Heart Rhythm 10, no. 9 , pp. 1318-1324, 2013. [19] S. Y. Sarairah, B. Woodbury, N. Methachittiphan, D. M. Tregoning, A. R. Sridhar and N. Akoum, "Esophageal thermal injury following cryoballoon ablation for atrial fibrillation," Clinical Electrophysiology 6, no. 3, pp. 262-268., 2020. [20] N. Kulkarni, W. Su and R. Wu, "How to prevent, detect and manage complications caused by cryoballoon ablation of atrial fibrillation," Arrhythmia & Electrophysiology Review 7, no. 1, p. 18, 2018. [21] Y. Xia, B. Liu, P. Ye and B. Xu, "Thermal field and tissue damage analysis of cryoballoon ablation for atrial fibrillation," Applied Thermal Engineering. 142, p. 524–529, 2018. [22] A. A. Gage, J. M. Baust and J. G. Baust, "Experimental cryosurgery investigations in vivo," Cryobiology 59, no. 3 , pp. 229-243., 2009. [23] P. Crosetto, S. Deparis, G. Fourestey and A. Quarteroni, "Parallel algorithms for fluid- structure interaction problems in haemodynamics."," SIAM Journal on Scientific Computing 33, no. 4 , pp. 1598-1622., 2011. [24] A. M. Hess, Simulation and Design of Soft Robotic Swimmers with Artificial Muscle, Michigan State University, 2019. [25] Z. Lin, A. Hess, Z. Yu, S. Cai and T. Gao, "A fluid–structure interaction study of soft robotic swimmer using a fictitious domain/active-strain method," Journal of Computational Physics 376 , pp. 1138-1155., 2019. 83 [26] A. W. Bergersen, A. Slyngstad, S. Gjertsen, A. Souche and K. Valen-Sendstad, "turtleFSI: A robust and monolithic FEniCS-based fluid-structure interaction solver," Journal of Open Source Software 5, no. 50, p. 2089, 2020. [27] D. Farnell, T. David and D. Barton, "Numerical simulations of a filament in a flowing soap film," International Journal for Numerical Methods in Fluids 44, no. 3, pp. 313-330., 2004. [28] C. S. Peskin, "Numerical analysis of blood flow in the heart," Journal of computational physics 25, no. 3, pp. 220-252., 1977. [29] L. Zhu and C. S. Peskin, "Simulation of a flapping flexible filament in a flowing soap film by the immersed boundary method," Journal of Computational Physics 179, no. 2, pp. 452- 468., 2002. [30] C. D. Eggleton and A. S. Popel, "Large deformation of red blood cell ghosts in a simple shear flow," Physics of fluids 10, no. 8 , pp. 1834-1845., 1998. [31] D. Kamensky, M.-C. Hsu, D. Schillinger, J. A. Evans, A. Aggarwal, Y. Bazilevs, M. S. Sacks and T. J. Hughes, "An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves," Computer methods in applied mechanics and engineering 284 , pp. 1005-1053, 2015. [32] C. Kadapa, W. G. Dettmer and D. Perić, "A fictitious domain/distributed Lagrange multiplier based fluid–structure interaction scheme with hierarchical B-spline grids," Computer Methods in Applied Mechanics and Engineering 301, pp. 1-27, 2016. [33] Z. Yu, "A DLM/FD method for fluid/flexible-body interactions," Journal of computational physics 207, no. 1 , pp. 1-27., 2005. [34] C. Antoci, M. Gallati and S. Sibilla, "Numerical simulation of fluid–structure interaction by SPH," Computers and Structures 85, p. 879–890, 2007. [35] S. Chen and G. Doolen, "Lattice boltzmann method for fluid flows," Annu. Rev. Fluid Mech. 30, pp. 329-364, 1998. [36] X. He and L. Luo, "Lattice boltzmann model for the incompressible navier–stokes equation," Journal of Statistical Physics 88, pp. 927-944, 1997. [37] H. Huang and X.-Y. Lu, "An ellipsoidal particle in tube Poiseuille flow," Journal of Fluid Mechanics 822, pp. 664-688, 2017. [38] Q. Wei, X. Yuan-Qing, T. Fang-Bao, T.-X. Gao, X.-Y. Tang and W.-H. Zu, "IB-LBM simulation on blood cell sorting with a micro-fence structure," Bio-Medical Materials and Engineering 24, no. 1, pp. 475-481, 2014. 84 [39] S.-G. Cai, Computational fluid-structure interaction with the moving immersed boundary method, Doctoral dissertation, Compiègne, 2016. [40] C. M. Scotti, A. D. Shkolnik, S. C. Muluk and E. A. Finol, "Fluid-structure interaction in abdominal aortic aneurysms: effects of asymmetry and wall thickness," Biomedical engineering 4 , pp. 1-22., 2005. [41] S. Tada and J. M. Tarbell, "A computational study of flow in a compliant carotid bifurcation– stress phase angle correlation with shear stress," Annals of biomedical engineering 33, pp. 1202-1212., 2005. [42] J. H. Lee, A. D. Rygg, E. M. Kolahdouz, S. Rossi, S. M. Rett, N. Duraiswamy, L. N. Scotten, B. A. Craven and B. E. Griffith, "Fluid–structure interaction models of bioprosthetic heart valve dynamics in an experimental pulse duplicator," Annals of biomedical engineering 48, pp. 1475-1490., 2020. [43] T. E. Tezduyar, K. Takizawa, T. Brummer and P. R. Chen, "Space–time Fluid– structure Interaction Modeling of Patient-Specific Cerebral Aneurysms," Int. J. Numer. Methods Biomed. Eng., 27(11), p. pp. 1665–1710. , 2011. [44] T. R, W. NB, H. N, D. AW, W. AR, H. AD, D. J, F. DP, M. J, Y. GZ and T. SA., "Fluid– structure Interaction Analysis of a Patient-Specific Right Coronary Artery with Physiological Velocity and Pressure Waveforms," Commun. Numer. Methods Eng., 25(5), p. 565–580., 2009. [45] L. Taelman, J. Degroote, P. Verdonck, J. Vierendeels and P. Segers, "Modeling hemodynamics in vascular networks using a geometrical multiscale approach: numerical aspects," Annals of biomedical engineering 41 , pp. 1445-1458, 2013. [46] G. Luraghi, F. Migliavacca and J. F. R. Matas, "Study on the accuracy of structural and FSI heart valves simulations," Cardiovascular Engineering and Technology 9, pp. 723-738, 2018. [47] H. Gao, L. Feng, N. Qi, C. Berry, B. E. Griffith and X. Luo, "A coupled mitral valve—left ventricle model with fluid–structure interaction," Medical engineering and physics 47 , pp. 128-136., 2017. [48] C. S. Peskin, "Flow patterns around heart valves: a numerical method," Journal of computational physics 10, no. 2, pp. 252-271, 1972. [49] C. S. Peskin, "Numerical analysis of blood flow in the heart," Journal of computational physics 25, no. 3 , pp. 220-252., 1977. 85 [50] M. D. M., C. S. Peskin and E. L. Yellin, "Fluid dynamics of the mitral valve: physiological aspects of a mathematical model," American Journal of Physiology-Heart and Circulatory Physiology 242, no. 6, pp. H1095-H1110, 1982. [51] K. B. Chandran and H. Kim, "Computational mitral valve evaluation and potential clinical applications," Annals of biomedical engineering 43 , pp. 1348-1362, 2015. [52] E. Oñate, M. Manzan, E. Oñate and M. Manzan, "Stabilization Techniques for Finite Element Analysis of Convection-Diffusion Problems," International Center for Numerical Methods in Engineering, 2000. [53] K. Hansen, A. Arzani and S. C. Shadden, "Finite element modeling of near‐wall mass transport in cardiovascular flows," International journal for numerical methods in biomedical engineering 35, no. e3148, 2019. [54] V. John and P. Knobloch, "A computational comparison of methods diminishing spurious oscillations in finite element solutions of convection-diffusion equations," p. 122–136, 2006. [55] R. Codina, "A discontinuity-capturing crosswind-dissipation for the finite element solution of the convection-diffusion equation," Computer Methods in Applied Mechanics and Engineering. 110, p. 325–342, 1993. [56] W. Kress and P. Lötstedt, "Time step restrictions using semi-explicit methods for the incompressible Navier–Stokes equations," Computer Methods in Applied Mechanics and Engineering. 195 , p. 4433–4447, 2006. [57] P. Ryzhakov and J. Marti, "A Semi-Explicit Multi-Step Method for Solving Incompressible Navier-Stokes Equations," Applied Sciences. 1 , 2018. [58] K. Selim, A. Logg and M. Larson, "An Adaptive Finite Element Splitting Method for the Incompressible Navier-Stokes Equations," Computer Methods in Applied Mechanics and Engineering. 209–212, p. 54–65, 2012. [59] T. Gelhard, G. Lube, M. Olshanskii and J. Starcke, "Stabilized finite element schemes with LBB-stable elements for incompressible flows,," Journal of Computational and Applied Mathematics, vol. 177, p. 243–267., 2005. [60] D. Boffi, F. Brezzi and M. Fortin, Mixed Finite Element Methods and Applications, 44, 2013. [61] T. Tezduyar and S. Sathe, "Stabilization Parameters in SUPG and PSPG Formulations," Journal of Computational and Applied Mechanics. 4, p. 71–88, 2003. 86 [62] J. Peterson, A. Lindsay and F. Kong, "Overview of the incompressible Navier–Stokes simulation capabilities in the MOOSE framework," Advances in Engineering Software. 119, p. 68–92., 2018. [63] S. Haßler, A. M. Ranno and M. Behr, "Finite-element formulation for advection-reaction equations with change of variable and discontinuity capturing," Computer Methods in Applied Mechanics and Engineering 369, no. 113171, 2020 . [64] O. Botella and R. Peyret, "Benchmark spectral results on the lid-driven cavity flow," Computers & Fluids. 27 , p. 421–433, 1998. [65] T. Gao, Y. Tseng and X. Lu, "An improved hybrid Cartesian/immersed boundary method for fluid-solid flows," International Journal for Numerical Methods in Fluids. 55, p. 1189–1211, 2007. [66] I. Rodriguez, O. Lehmkuhl, M. Soria, S. Gómez, M. Domínguez-Pumar and L. Kowalski, "On the boundary layer development and heat transfer from a sphere atmoderate Reynolds numbers," Tenth International Conference on Computational Fluid Dynamics, p. 9–13, 2018. [67] T. Fastl, C. Tobon-Gomez, A. Crozier, J. Whitaker, R. Rajani, K. McCarthy, D. Sanchez- Quintana, S. Ho, M. O’Neill, G. Plank, M. Bishop and S. Niederer, "Personalized computational modeling of left atrial geometry and transmural myofiber architecture," Medical Image Analysis. 47 , p. 180, 2018. [68] J. Lantz, V. Gupta, L. Henriksson, M. Karlsson, A. Persson, C. Carlhäll and T. Ebbers, "Impact of Pulmonary Venous Inflow on Cardiac Flow Simulations: Comparison with In Vivo 4D Flow MRI," Annals of Biomedical Engineering. 47, p. 413, 2019. [69] A. Klein, W. Stewart, J. Bartlett, G. Cohen, F. Kahan, G. Pearce, K. Husbands, A. Bailey, E. Salcedo and D. Cosgrove, "Effects of mitral regurgitation on pulmonary venous flow and left atrial pressure: An intraoperative transesophageal echocardiographic study," J Am Coll Cardiol. 20, p. 1345–1352., (992. [70] H. Pennes, "Analysis of tissue and arterial blood temperatures in the resting human forearm," J Appl Physiol. 1 , p. 93–122., 1948. [71] T. Marinko, J. Dolenc and C. Bilban-Jakopin, "Cardiotoxicity of concomitant radiotherapy and trastuzumab for early breast cancer," Radiology and oncology 48, vol. no. 2, p. 105. , 2014. [72] O. Smiseth, C. Thompson, K. Lohavanichbutr, H. Ling, J. Abel, R. Miyagishima, S. v. Lichtenstein and J. Bowering, "The pulmonary venous systolic flow pulse—its origin and relationship to left atrial pressure," J Am Coll Cardiol. 34, p. 802–809, 1999. 87 [73] A. Gonzalo, M. García-Villalba, L. Rossini, E. Durán, D. Vigneault, P. Martínez-Legazpi and O. Flores, "Non-Newtonian Patient-specific Simulations of Left Atrial Hemodynamics," bioRxiv, pp. 2021-06., 2021. [74] X. Morales, J. Mill, G. Delso, F. Loncaric, A. Doltra, X. Freixa, M. Sitges, B. Bijnens and O. Camara, "4D flow magnetic resonance imaging for left atrial haemodynamic characterization and model calibration," In Statistical At Atlases and Computational Models of the Heart. M&Ms and EMIDEC Challenges: 11th International Workshop, STACOM 2020, Springer International Publishing, pp. 156-165., 2021. [75] J. Egger, M. Gall, J. Wallner, P. Boechat and A. Hann, "HTC Vive MeVisLab integration via OpenVR for medical applications," PLOS ONE 12(3): e0173972., 2017. [76] A. Fedorov, R. Beichel, J. Kalpathy-Cramer, J. Finet, J.-C. Fillion-Robin, S. Pujol and C. Bauer, "3D Slicer as an image computing platform for the Quantitative Imaging Network," Magnetic resonance imaging 30, vol. no. 9 , pp. 1323-1341., 2012. [77] B. Parvez, V. Pathak, C. Schubert and M. Wood, "Comparison of Lesion Sizes Produced by Cryoablation and Open Irrigation Radiofrequency Ablation Catheters," Journal of Cardiovascular Electrophysiology. 19, p. 528–534., 2008. [78] R. Müssig, M. Heinke and J. Hörth, "Cryoballoon model and simulation of catheter ablation for pulmonary vein isolation in atrial fibrillation, Current Directions in Biomedical Engineering. 4," p. 473–475., 2018. [79] J. Baust, A. Robilotto, K. Snyder, R. v. Buskirk and J. Baust, "Evaluation of a new epicardial cryoablation system for the treatment of Cardiac Tachyarrhythmias," Research Article Trends Med. 18 , p. 1–9, 2018. [80] J. Ghosh and M. McGuire, "Atrial flow dynamics as a determinant of tissue temperature during balloon cryoablation," EP Europace. 20, p. f451–f457., 2018. [81] J. Baust, A. Robilotto, P. Guerra, K. Snyder, R. v. Buskirk, M. Dubuc and J. Baust, "Assessment of a novel cryoablation device for the endovascular treatment of cardiac tachyarrhythmias," SAGE Open Med., vol. 6, 2018. [82] Z. Yu, X. Shao and A. Wachs, "A fictitious domain method for particulate flows with heat transfer," Journal of Computational Physics 217, vol. no. 2, pp. 424-452, 2006. [83] Y. Yu, "Numerical Methods for Fluid-Structure Interaction: Analysis and Simulations," PhD diss., Brown University, 2014. [84] X. Wang and L. T. Zhang, "Interpolation functions in the immersed boundary and finite element methods," Computational Mechanics 45, no. 4 , pp. 321-334, 2010. 88 [85] X. Yang, Z. L. Xing Zhang and G.-W. He, "A smoothing technique for discrete delta functions with application to immersed boundary method in moving boundary simulations," Journal of Computational Physics 228, no. 20, pp. 7821-7836, 2009. [86] [Online]. Available: https://docs.python.org/3/library/argparse.html. [87] S. Turek and J. Hron, "Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow," Springer Berlin Heidelberg, 2006. [88] F.-B. Tian, H. Dai, H. Luo, J. F. Doyle and B. Rousseau, "Fluid–structure interaction involving large deformations: 3D simulations and applications to biological systems," Journal of computational physics, no. 258, pp. 451-469, 2014. [89] Z. Yu, Y. Wang and X. Shao, "Numerical simulations of the flapping of a three-dimensional flexible plate in uniform flow," Journal of Sound and Vibration 331, no. 20 , pp. 4448-4463, 2012. [90] É. d. Madron, V. Chetboul and C. Bussadori, Clinical Echocardiography of the Dog and Cat, 2015, p. 1–348.. [91] G. D’Angelo, E. Moro, G. Nicolosi, V. Dall’Aglio, R. Mimo, S. Mangano and D. Zanuttini, "Color Doppler identification of early diastolic turbulence in the left atrium in patients with mitral valve insufficiency: persistence of regurgitation or inertia phenomenon?," Giornale Italiano Di Cardiologia. 20 , p. 700–704, 1990. [92] S. Singh and S. Kumar, "A study on the effect of metabolic heat generation on biological tissue freezing," ScientificWorldJournal, 2013. [93] "Ablation Products for Atrial Fibrillation - CryoConsole Cardiac Cryoablation System | Medtronic, (n.d.).," [Online]. Available: https://www.medtronic.com/us-en/healthcare- professionals/products/cardiac-rhythm/ablation-atrial-fibrillation/cryoconsole-cardiac- cryoablation.html. [94] N. V. C. Gillespie, A. Lin and B. Knight, "Ice formation in the left mainstem bronchus during cryoballoon ablation for the treatment of atrial fibrillation," Heart Rhythm. 13 , p. 814–815, 2016. [95] S. Bordignon, A. Fürnkranz, B. Schmidt and K. Chun, "Remaining ice cap on second- generation cryoballoon after deflation," Circulation: Arrhythmia and Electrophysiology. 5, 2012. 89 [96] A. Odutayo, C. X. Wong, A. J. Hsiao, S. Hopewell, D. G. Altman and C. A. Emdin, "Atrial fibrillation and risks of cardiovascular disease, renal disease, and death: systematic review and meta-analysis," bmj 354, 2016. 90 APPENDIX A : EFFECTS OF STABILIZATION SCHEMES First, we consider the LA+CB case without the LSIC and PSPG stabilization schemes for the momentum equation. The results for both lesion area (Figure A.1a) and power absorbed across the CB are similar to that shown in Figure 4.4a and e. Subsequently, the case without LSIC and PSPG stabilization for the momentum equation is utilized to investigate the isolated effects of crosswind stabilization on the energy equation (Figure A.1b) and the momentum equation (Figure A.1c). When the crosswind stabilization is disabled for the energy equation, localized undershoots (−132! C) /overshoots (216! C) are observed in the temperature field in the vicinity of the CB that is outside the problem’s temperature range (between −70! C to 37! C, see Figure A.1b). On the other hand, when crosswind stabilization is disabled for the momentum equation, we observe an unsteady lesion, and the position of the lesion oscillates in way that is non-physical even though temperature gradients downstream of the CB are better captured without undershooting/overshooting of the temperature (Figure A.1c). Based on these simulations, we conclude that the crosswind stabilizing operator (Eq. (3.16)) is necessary and serves mainly to capture sharp gradients in the field variable. (a) (b) (c) Figure A.1 Effect of crosswind stabilization: (a) LA+CB case without PSPG and LSIC stabilization for the momentum equation (b) LA+CB case without PSPG and LSIC stabilization for the momentum equation and without crosswind stabilization for the energy equation (c) LA+CB case without PSPG, LSIC and crosswind stabilization for the momentum equation. 91 APPENDIX B : CRYOBALLOON POSITIONING Under identical simulation conditions (section 4.1), a total of 12 CB positions (LA+CB case + P1 - P11) are simulated to study the effects of CB position on lesion formation. The exact positions of the CB (center O) with respect to 3 constant points A, B and C on the LA surface are given in Table B.1. Here, we also introduce 𝛼 which is the angle between the CB’s minor axis and OA (see Figure B.1). Note that point A is located centrally at the RIPV opening, and point B and C are selected arbitrary at the top and side of the RIPV where it merges with the LA. Positioning the CB in this manner provides a range of different balloon-tissue contact areas 𝑨𝒄 to analyze incomplete occlusion in case of cryotherapy. Figure B.1: Geometrical positioning of the CB (center O) in the RIPV ostium. Point A is located centrally at the RIPV opening. n is normal to the plane passing through the two major axes of the CB and the angle between n and OA is defined as 𝜶. Point B and C are selected arbitrary at the top and side of the RIPV where it merges with the LA. Table B.1: The distance OA, OB, OC, angle 𝜶 and balloon tissue contact area 𝑨𝒄 for all CB positions. CB position OA (mm) OB (mm) OC (mm) 𝜶 (deg) 𝑨𝒄 (mm2) P1 31.69 16.25 15.36 7.37 252.9 P2 33.95 18.67 16.66 20.16 116.7 P3 30.67 21.88 16.31 14.27 207 P4 30.11 17.43 14.81 9.18 223.4 92 Table B.1 (cont’d) P5 29.67 21.30 15.72 14.30 213.4 P6 28.45 23.40 16.52 16.26 444.8 P7 29.71 19.74 14.94 19.85 193.2 P8 30.14 20.74 15.82 10.94 187 P9 31.14 23.18 16.66 11.79 259 P10 30.52 20.54 16.10 4.21 167.4 P11 31.09 19.67 14.81 14.52 176 93