Unsteady simulation of Carreau fluid model for Pulsatile blood flow through a 3D Bifurcating Artery

CONTENT:

1. Abstract.
2. Introduction.
3. Problem Specification.
4. Geometry and Mesh Generation.
5. Physical Model Set-up.
6. Numerical Solution.
7. Simulation and Results.
8. Verification and Validation.
9. Conclusion.

______________________________________________________________________________________________________________________________________________________________________________ 

                                                                                                                                                            Abstract

Knowledge of blood flow dynamics in arteries is essential for the diagnosis and treatment of cardiovascular and vascular diseases. The health of arteries is decided by the deposition of fatty material on their inner walls. Due to the numerous challenges in performing in-vivo experiments, the use of computational fluid dynamics (CFD) to simulate pulsatile flow in the blood artery has become a promising alternative approach. This report aims to describe the calculation of Wall Shear Stress (WSS) around the walls of the bifurcating artery using the Carreau-model to model blood's viscosity. It also aims to compare the velocity profile across different cross-sections of the bifurcating artery using the software like ANSYS Fluent. Besides this, tracing of massless particles injected at specific locations is also observed to deliver more insight into the blood flow in the bifurcating artery. The role of this study is vital to understand the diagnosis and treatment of cardiovascular diseases.

______________________________________________________________________________________________________________________________________________________________________________

INTRODUCTION

Blood flow in human arteries is a particular area of interest because of its importance in the various health care process. It is essential since most people in developed countries die because of heart attacks or some other diseases related to blood flow. Atherosclerosis or heart stroke are two fatal disease examples related to the blood flow. Atherosclerosis affects the flow of blood. This effect, at carotid artery bifurcation, may be vital. Atherosclerosis affects millions of people worldwide, leading to such diseases as heartbreak and stroke. Complex hemodynamics plays a critical role in the development of atherosclerosis, the process of aging, and many other disease processes. Bio-fluid mechanics play a significant role in the cardiovascular system, so it is essential to understand the forces and movement of blood cells and whole blood and the interaction between the blood cells and the vessel wall. According to the study of Liepsch, D., viscosity, and flow behavior change, specifically the creation of vortices and flow disturbances, can be used to show how medication can influence flow behavior. Experiments have demonstrated that hemodynamics may have a strong influence on the creation of aneurysms and varicose veins. Thus, health care practitioners must understand fluid dynamic factors such as flow rate ratio, pressure and velocity gradients, flow behavior, velocity distribution, and shear stress on the wall and blood cells. These mechanical factors are primarily responsible for depositing the cells and lipids, which is a leading cause of atherosclerosis. The interaction between blood cells and of the cells with the vessel may cause the formation of plaques and agglomerations. These deposits are found predominantly at arterial bends and bifurcations where blood flow is disturbed, and secondary flows are developed, and flow separation regions are found.

Currently, in the era of the vast availability of computational fluid dynamics (CFD) tools and machines with high computational efficiency, there is a possibility to perform numerical simulations of the complex physical phenomena occurring inside the human organism, e.g., blood flow in the cardiovascular system. Several fields of this science have already been intensely explored. Hence, a tremendous amount of publications appeared. General blood flow and hydrodynamics (rheology, viscosity, governing equations) and flow inside significant arteries have been broadly discussed. CFD allows for efficient and accurate computations of hemodynamic features of both normal and abnormal situations in the cardiovascular system, in vivo stimulation of coronary artery flow changes. CFD is different from medical imaging visualization. Medical imaging techniques such as coronary angiography or computed tomography angiography provide anatomic alterations of the coronary artery wall due to plaques, thus allowing only assessing the degree of lumen changes such as stenosis or occlusion. In contrast, CFD analysis identifies hemodynamic changes in the coronary artery, even before the plaques are formed at the artery wall or can occlude the vessels. Therefore, to some extent, CFD allows early detection of coronary artery disease and improves the understanding of the progression of plaques, which are considered of paramount importance to clinical treatment.

The complex flow fields which exist in the region of arterial bifurcations are characterized by substantial spatial and temporal variations in wall shear and particle trajectory, creating environments considerably different than those found in simple in vitro systems, such as Couette or channel flow, that are often employed in the study of fluid dynamic effects on cell function. Understanding the interactions between blood flow and the biological behavior of cells in the arterial wall will undoubtedly require more excellent knowledge of the response of cells to flow field phenomena that are representative of those occurring in actual arterial bifurcations. Bio-medical researchers rely on computational fluid dynamics to model and understand the physical mechanisms behind the formation and progression of hemodynamic disorders. Wall shear stress (WSS) exerted on the walls of the blood vessel due to the flow of blood is one of the main pathogenic factors leading to the development of such disorders. The magnitude and distribution of the WSS in a blood vessel can provide an insight into the locations of possible aneurysm growth. Moreover, blockages that build up over time can be predicted by having a qualitative understanding of the flow profile. 

Computational Fluid Dynamics can be used for modeling, and understanding such vital internal flows and insights gained from such studies can help design patient-specific treatments. This project work is devoted to model the three-dimensional internal blood flow in a bifurcating artery. The required computational mesh is generated, and appropriate boundary conditions are set up that was needed for the simulation. The Non-Newtonian behavior of blood flow was modeled using the Carreau model of viscosity. Moreover, a realistic time-varying boundary condition will be implemented at the inlet using User Defined Functions (UDF) in order to mimic the pulsatile nature of blood flow. As part of post-processing results, Wall Shear Stress (WSS) around the walls of the artery was calculated, and plots of velocity contour on different cross-sections of the artery are obtained to compare the velocity distribution throughout the length of the artery. Studying blood velocity distribution in arteries can be proved very useful in designing implantable devices, while it can help choose between surgical and interventional procedures. Furthermore, tracing of massless particles injected at specific locations during the physics set-up is also performed to deliver more insight into the blood flow in the bifurcating artery. 

______________________________________________________________________________________________________________________________________________________________________________

PROBLEM SPECIFICATION

Blood flows through the bifurcating artery from the inlet (to the left in the figure below) and exits from the two outlets (to the right). The diameter of the artery at the inlet is around `6.3 mm`. The diameter of outlet 1 is around `4.5 mm`, and the diameter of outlet 2 is around `3.0 mm`. The density of blood is `1060 kg//m^3`. As blood is a non-newtonian fluid, the coefficient of viscosity of blood is not a constant, but rather, it is a function of velocity gradients. Here we use the Carreau-model to model blood's viscosity. Since blood flow is pulsatile and cyclic, the velocity profile at the inlet is a function of time. The pressure at the outlet is defined to be constant (`"100 mm of Hg"`). More details on boundary conditions will be provided next to the next lesson, "Physical Model Set-up."

______________________________________________________________________________________________________________________________________________________________________________

GEOMETRY and MESH GENERATION

Geometry

This is a model created from a luminal casting of a carotid artery bifurcation and downloaded from the Grabcad community to carry out the simulation. The model was initially designed and developed in Solidworks by Aaron. The geometry consists of one big inlet and two small daughter outlets. The diameter of the mother (inlet) tube is almost uniform for the entire length, while those of the daughter (outlets) tubes lopsidedly from the bifurcation point to the exit. Different views of the final geometry are shown below.

The area specification of complete geometry is written in the following table.

Mesh Generation

The mesh generated is consisted only of tetrahedral elements with curvature as the size function. First, the body sizing is applied to the whole domain with element size as `0.8 mm`. Then again, second body sizing is used by creating a sphere of influence whose center coincides with the geometric center of the body, and radius equals to `12mm` with `0.4mm` as element size. At last, Inflation is applied to the whole geometry with artery walls as the boundary scoping method having a maximum thickness of `0.6 mm`. The generated mesh contains `45240` nodes and `147611` elements. The first two figure shows only the inflation mesh at outlets and last three figures shows the complete mesh.

______________________________________________________________________________________________________________________________________________________________________________ 

PHYSICAL MODEL SET-UP

After generating the mesh, now its important to describe the physical nature of our flow problem; in this section, we are going to define the physics of the described situation. This section will contain three parts. In first, we will describe the Governing equation and physical properties of fluid, i.e., kinematic viscosity (given by Carreau viscosity model) and other parameters, if required. In the second part, we will initialize our initial and boundary conditions. The last part will active the discrete phase in the physical model tab by injecting two massless particles at specific locations in the artery.

Governing Equations and Specification of Physical Properties

Before starting a CFD simulation, it is always good to take a look at the governing equations underlying the physics. In this case, although we have additional complexities such as pulsatile flow and non-Newtonian fluids, the governing equations are the same as any other fluids problem. The most fundamental governing equations are the continuity equation and the Navier-Stokes equations. Given below a quick review of the equations.

Continuity Equation:          `(del rho)/(delt)+ nabla(rhov)=0`

However, as blood can be regarded as an incompressible fluid, the rate of density change is zero. Thus the continuity equation above can be further simplified in the form below:    `nabla*v=0`

The Navier-Stokes Equation:      `rho((delv)/(delt)+v*nablav)=-nablap+munabla^2v+f`     

One thing to notice in the Navier-Stokes equation is that the viscosity coefficient of `mu` is not a constant but rather a function of shear rate. Blood gets less viscous as the shear rate increases (shear thinning). Here, we model the blood viscosity using the Carreau fluids model. Carreau model is used for non-Newtonian shear-thinning fluids. It is characterized by an apparent viscosity `mu` which gradually decreases with increasing shear rate. At low shear rates, the apparent viscosity approaches a Newtonian plateau where the viscosity is independent of shear rate (zero shear viscosity, `mu_0`). Furthermore, a similar plateau at very high shear rates (infinite shear viscosity, `mu_("inf")`). This model is mostly used for food, beverages, and also blood flow applications. The properties of blood as given in the below table. The following correlation represents the Carreau model:

                                                                                      `mu_(eff)(dot gamma)=u_["inf"]+(mu_o-mu_["inf"])(1+(lamda.dot gamma)^2)^((n-1)/2)`

Specifying Initial and Boundary Conditions

Governing equations are the same for every CFD problem, and then there should be something that makes our results unique for every flow problem. This is done by Boundary conditions, which are unique for every CFD problem. 

Wall: The most straightforward boundary condition to determine is the artery wall. We need to define the wall regions of this model and set it to the "wall." From a physical viewpoint, the "wall" condition dictates that the velocity at the wall is zero.

Inlet: As we know, mammalian blood flow is pulsatile and cyclic in nature. Thus the velocity at the inlet is not set to be a constant, but instead, in this case, it is a time-varying periodic profile. Pulsatile flow at the inlet implemented using a user-defined function (UDF), which mimics blood flow in the human body. The pulsatile profile within each period is considered to be a combination of two phases. During the systolic phase, the velocity at the inlet varies in a sinusoidal pattern. The sine wave during the systolic phase has a peak velocity of `0.5m//s` and a minimum velocity of `0.1m//s`. Assuming a heartbeat rate of `120` per minute, the duration of each period is `0.5s`. This model for pulsatile blood flow is proposed by `"Sinnott M. et al., 2006"`. The UDF is shown in the below graph:

To describe the profile more clearly, mathematical description and the code that is directly interpreted in fluent is also given below:

/***********************************************************************/
/* UDFs for specifying time dependant velocity profile boundary condition
*/ 
/***********************************************************************/


#include "udf.h"//file that contains definitions for define functions and fluent operations
#define PI 3.141592654

DEFINE_PROFILE(inlet_velocity,th,i)
{
	face_t f;
	begin_f_loop(f,th)
		double t = (CURRENT_TIME*2-floor(CURRENT_TIME*2))/2; //t is the local time within each period

	{
		if(t <= 0.218)
			F_PROFILE(f,th,i) = 0.5*sin(4*PI*(t+0.0160236));
		else
			F_PROFILE(f,th,i) = 0.1;
	}
	end_f_loop(f,th);
}

Outlet: The systolic pressure of a healthy human is around `"120 mm Hg"` and the diastolic pressure of a healthy human is around `"80 mm Hg"`. Thus, taking the average pressure of the two phases, we use `"100 mm Hg"` (around `13332` Pascal) as the static gauge pressure at the outlets.

Particle Injections

The motivation for tracking the particles here is to set up for creating particle trajectories later on in the last section. The discrete phase model is activated by injecting the particles in the artery, which provides a better insight into the blood flowing in the bifurcating artery.  Two massless particles are injected at specified locations, as shown in the following table.

______________________________________________________________________________________________________________________________________________________________________________ 

NUMERICAL SOLUTIONS

In this section, First reference entities are specified that are used by Ansys Fluent for calculating non-dimensionalization, then solution monitors were set up for monitoring Drag and also to set up calculation activities to export solution data for post-processing later on. The drag is calculated by integrating the shear and pressure at the wall. The drag coefficient is then calculated by non-dimensionalizing the drag. The reference entities used in the non-dimensionalization are defined in the Reference Values panel in Fluent. Note that the "Reference Values" will not change your solution for the velocity and pressure at the cell centers. It will affect only the drag coefficient and any other non-dimensional quantities calculated from the solution. So its imperative to provide proper reference values, and for this particular situation, the reference values used are shown below.

For solution methods, the SIMPLE scheme was used in the simulation for pressure-velocity coupling and second-order discretization scheme for pressure and second-order upwind for momentum. The transient model has been used for simulation with the time step of `0.005s` and `100` timesteps in total; thus, the simulation is performed for flow time up to `0.5s`. The hybrid method is used for the initialization of the solution, while the maximum iterations at each timestep for the solution to get converged is set to `200`.

______________________________________________________________________________________________________________________________________________________________________________ 

SIMULATIONS and RESULTS