*Institute of Industrial Science, The University of Tokyo, Tokyo, Japan **Department of Mechanical Engineering, School of Engineering, The University of Tokyo, Tokyo, Japan

Abstract

Blood vessels are able to expand in order to let more blood flow through them, as well as contract to help control the flow of blood. Fluid-Structure Interaction (FSI) modeling simulates this process of expansion and contraction for more accurate results. The results of two analyses, namely the Fluid and FSI, are presented to compare a difference between the rigid and deformable walls, respectively. In the FSI analysis, strong coupling method (direct method) is used to solve the fundamental descretized equations to satisfy the geometrical compatibility and the equilibrium conditions on the interface between the fluid, representing the blood flow, and the hyper-elastic structure, representing the vascular tissue. A high order Mooney-Rivlin material is used to model the vascular tissue. Automatic mesh update method is adjusting the fluid mesh under the deformation of the structure domain along the time. Peripheral vessels network is applied as outflow boundary condition to the 3D Fluid and FSI simulations. The peripheral network is modeled as a binary symmetric tree attached to the outlet(s) of the 3D model. The small arteries are modeled using 1D scheme, while the arterioles and capilaries are modeled using the structure tree impedance (STI) model. The coupling of the 3D Fluid-Structure finite element model and the multi-scale model (1D - 0D) facilitates the use of more realistic boundary conditions leading to more accurate information on the hemodynamic factors, e.g. wall shear stress (WSS). As, for example, low density lipoprotein tends to deposit in the areas of blood vessels with low WSS. The degeneration of blood vessel walls is often initiated by hemodynamic forces, thus it is valuable to obtain the information on these factors for accurate prevention and prognosis of diseases such as atherosclerosis or middle cerebral aneurysm.

1. Introduction

Atherosclerosis often occurs in the carotid artery. The rupture of an atherosclerosis plaque can result in a stroke in the cerebral or coronary arteries. A cerebral aneurysm, defined as a focal dilatation of the arterial wall, occurs usually at the branching points of arteries supplying blood to the brain. Genetics and risk factors are partially responsible, but also hemodynamic factors, e.g. pressure, flowrate, shear stress, have been proved to play a significant role in forming these vascular diseases, namely the atherosclerosis [1, 2] and aneurysms [3, 4], which are responsible for significant morbidity and mortality.

FSI simulations, in addition to the hemodynamic factors, provide also data on the vessel wall dynamics in the patient-specific models, such as the deformations, principal strains and stresses. The traditional FSI techniques are based on the arbitrary Lagrangian-Eulerian (ALE) method. In this study, the strongly coupled ALE method is used to address the added-mass effect, which is typical for loosely coupled time-stepping algorithms [5, 6]. An example of strongly coupled ALE method used for the coupling of an elastic structure with an incompressible fluid can be found in [7].

There are many models to describe the outflow boundary conditions, such as Windkessel models or lumped parameter models. However, the problem with these models is that a previous knowledge of certain parameters, such as the resistance and compliance, is required to conduct the simulations, therefore, a multi scale model is used, developed based on [8], which does not bear this limitation.

2. Methods

With 3D FSI simulations only, the area of interest is isolated and therefore taken out of the context of the entire circulatory system. This leads to the lack of more realistic boundary conditions. Hence the need to carry out the computer simulations using multi-scale model, here by coupling the 3D FSI model with the 1D-STI model. The 3D FSI domain is coupled to the inflow of the reduced-order 1D-0D simulation using the explicit approach where the data are exchanged between the two models without reformulating the numerical method for the 3D algorithm. The data trasferred from the 3D model to the 1D model is the flowrate and the data returned to the 3D model is the pressure given to the fluid domain assuring proper deformations of the structure domain.

2.1. 3D FSI model

Both the fluid and the structure domains act on each other since the deformation of the structure is determined by both the external and the fluid forces and consequently it affects the fluid movement, as depicted in the Fig. 1.

2.1.1. Governing Equations

Fluid Domain
The blood flowing through a large blood vessel is considered laminar incompressible Newtonian and is governed by the continuity and Navier-Stokes equations.

(1)

(2)

(3)

(4)

The Eqs. 1-4 are continuity equation, Navier-Stokes equation, and Neumann and Dirichlet boundary conditions, respectively, where v is the velocity vector, p is the pressure component, ρ is the fluid density, μ is the viscosity, n is the normal vector, t̅ are the tractions prescribed on Neumann boundary Γ_σ and v̅ are the velocities prescribed on the Dirichlet boundary Γ_u and σ is the Cauchy stress tensor defined as follows

(5)

where I is the identity tensor and ε is the strain rate tensor defined as follows

(6)

Substituting Eq. 6 into 5 and combining it with Eq. 2 yields to the following form of the Navier-Stokes equation

(7)

Structure Domain
The fundamental ingredients describing the geometrically nonlinear continuum mechanics are kinematic equation, constitutive equation and equilibrium equation. These equations are summarized in the Tonti diagram depicted in Figure 2.
The C in the Figure 2, is the right Cauchy-Green tensor, the F is the deformation gradient tensor, the E is the Green-Lagrange strain tensor, the S is the second Piola-Kirchhoff stress tensor for hyper-elastic material using Mooney-Rivlin model with high-order strain energy density function W represented by the decoupled form [9]

(8)

with the isochoric contribution, W_iso(C), and the volumetric contribution, W_vol(J). The isochoric contribution used is in the following form

(9)

where the constants C₁-C₉ are derived to fit the stress-strain curve of the material to be modeled and the I₁ and I₂ are invariants of the right Cauchy-Green tensor C. The isochoric contribution is expressed as

(10)

where p is a scalar serving as an indeterminate Lagrange multiplier and can be identified as a hydrostatic pressure, J is the determinant of the deformation gradient tensor, J = det F. A material that keeps the volume constant during the motion, such as the blood vessel tissue, is characterized by the incompressibility constraint J = 1.

Substituting the strain tensor E into

(11)

the second Piola-Kirchhoff stress tensor can be written, in terms of the right Cauchy-Green tensor C, as

(12)

The blood vessel tissue can be assumed to be an isotropic material, i.e. its material properties are not dependent on the direction along which they are measured, therefore the dependence of W on the Cauchy-Green tensor C may be expressed by the material axes, thus the W can become a function of the invariants of C only, i.e. tensor C, as

(13)

with

.

with constraint I₃ = 1, because det C = (det F)² = 1. Now, considering Eq. 13, Eq. 12 can be rewritten as

(14)

It can be readily verified that

.

Considering the above derivations, the second Piola-Kirchhoff stress tensor S becomes

(15)

2.1.2. Finite element formulations

Using the Galerkin method and spatial discretization for Eqs. 1 and 7, the equations in matrix form are as follows

(16)

(17)

where M is the generalized mass matrix, the Λ is the generalized matrix of convective term, K_μ is the fluid viscosity matrix, G is the divergence operator matrix, P and V are the time derivatives of the pressure and velocity vectors, respectively, in ALE coordinates, and F is the external force vector.

The finite element discretization of the structure governing equations in total Lagrangian formulation yields the following linearized form of the equilibrium equation

(18)

where M is the mass matrix, K is the tangent stiffness matrix, F is the external force vector, Q^s is composed of the internal force and inertial force of the structure, and ∆Ü and ∆U are the increments of the acceleration and displacement vectors, respectively.

The fluid and the structure are coupled by equating the displacement, velocity and acceleration on the interface common to both the fluid and structure domains. The equilibrium on the interface is enforced by the element assemblage process. The system equations are then constructed in the incremental form as follows

(19)

where φ^fs and U^s are the unknown variables of the coupled system of the following form

.

where the P is the pressure vector of the fluid, V_i^f is the velocity vector of the fluid independent of the structure, V_c is the coupled velocity vector, V_i^s is the velocity vector of the structure independent of the fluid, U_c is the coupled displacement vector, and U_i^s is the displacement vector of the structure independent of the fluid. Here, in the Eq. 19, the F is the external force vector of the coupled system, Q^fs is the equivalent internal force vector, M^fs is the mass matrix of both the fluid and structure domains, C^f consists of the divergence, viscous and convective terms of the fluid, and the K^s is the tangent stiffness matrix of the structure. These coupled matrices and vectors are of the following form

.

.

.

.

where the subscripts c and i denote the coupling and independent parts of the variables, respectively.

The applied FSI solver is a fully implicit, monolithic ALE-FEM approach. The discretization is based on the Q1/P0 finite element pair, where continuous piecewise linear approximation is used for the kinematic variables (displacement, velocity, acceleration) and discontinuous piecewise constant approximation for the pressure, together with predictor-multicorrector algorithm (PMA) [10] based on the Newmark-β method for time stepping. The resulting nonlinear systems are solved using a Newton-Raphson method, while the linear subsystems are treated via Krylov Subspace method, namely generalized minimal residual method with ILU preconditioner.

2.2. Time integration method

The PMA is used based on the Newmark-β method to find the solution of the nonlinear dynamic response of a finite element system. Firstly, form the predictor for the displacement vector of the structure independent of the the fluid at the beginning of each time step as follows

(20)

where V^s is vector consisting of the coupled velocity vector and the velocity vector of the structure independent of the fluid, i.e. part of the φ^fs vector, defined as

.

as well as the predictors for the velocity and acceleration vectors at the beginning of the time step

(21)

(22)

Secondly, solve coupled system

(23)

where ^t+∆tF^(k) − ^tQ^fs(k) is the residual at the beginning of the k^th iteration for the time increment ∆t, thus the Eq. 23 can be rewritten as

(24)

where the effective mass matrix M^* is given as

(25)

Thirdly, form the correctors as follows

(26)

(27)

(28)

where the γ and β parameters are chosen to assure the integration accuracy and stability, in the present study as

.

Thus, the iteration procedure for solving coupled system using the PMA can be summarized as follows
  Begin time step loop
   Increment time step t + ∆t
   Form predictors, Eqs. 20-22
   Begin Newton-Raphson loop
    Solve the coupled system, Eq. 24
    Form correctors, Eqs. 26-28
   End Newton-Raphson loop
  End time step loop

2.3. Mesh update

On the interface with a deformable structure, the fluid nodes are attached to the structure. Therefore, mesh control method is used to maintain the mesh quality in the fluid domain. An elasticity Dirichlet boundary value problem is solved to control the mesh motion of interior nodes of the fluid domain to avoid the mesh distortion. The mesh update functions are called either at the end of every Newton-Raphson iteration or only when converged, i.e. at the end of each time iteration, depending on the given problem set with different level of deformations.

2.4. 1D model

The binary structural peripheral vessel network is modeled based on a fractal nature of arterial network, i.e. when a vessel bifurcates into two daughter vessels, the diameter of the daughter vessels is determined from the parent vessel's diameter using the power law

(29)

where ε = 2.76 is representing the optimal value for blood flow [11]. The peripheral network (1D-0D model), based on [8], is modeled as binary symmetric tree attached to the outlet(s) of the 3D model. The small peripheral arteries, with diameter down to 0.1 mm, are modeled using the 1D model. The 1D model is governed by three equations, namely the mass conservation equation 30, momentum equation 31 and the state equation 32.

(30)

(31)

(32)

where q represents the flow rate, A is the cross-sectional area, A₀ is the initial area, p is the pressure, ρ is the density, E is Young's modulus, h is the wall thickness and r₀ is the initial radius. The state equation relates the fluid influence on the vessel wall to its compliant properties.

The varying stiffness of the small arteries depending on their diameter can be accounted for by stating the Young's modulus E as a function of the vessel radius r₀ according to the formula based on compliance estimates [12,13] as follows

(33)

The k₁, k₂ and k₃ are constants of following values and units, k₁ = 2.0×10⁶ kg·sec⁻²·m⁻¹, k₂ = −2.253×10³ m⁻¹ and k₃ = 8.65×10⁴ kg·sec⁻²·m⁻¹ [8].

2.5. 0D model

After the given diameter of 0.1 mm as the threshold boundary for the 1D model, the following arterioles and capillaries are modeled using the STI model, down to the diameter of 5 μm. Using the analogy of electric circuits, the pressure drop within the STI model is related to the flow rate in the frequency domain

(34)

where Z(x,ω) is the impedance of the STI model. The Eq. 34 is used to derive the impedance in the frequency domain at the root, x = 0, of the STI model using the semi-analytical approach based on the linearized form of the governing Eqs. 30 and 31 as follows

(35)

where ω is the angular frequency, L is the length of the vessel, g is variable related to wave-propagation velocity and compliance of the vessel wall, g = cC, where c is the wave propagation speed and C is the compliance [8]. If ω = 0, then the relation 35 becomes

(36)

where μ is dynamic viscosity and l_rr is the length to the radius ratio.

The non-Newtonian effect of the blood is taken into consideration in this model. Hematocrit of the blood decreases with the diameter of the blood vessel once the diameter becomes smaller than 300 μm [14]. The apparent viscosity of the blood is taken as a function of the hematocrit and the diameter of the vessel [15].

3. Simulations and Results

The patient specific data are used to create the 3D model of the common, external, and internal carotid arteries (CCA, ECA, ICA, respectively) in this paper. The initial configuration of the model is depicted in Fig. 3. The resulting mesh configuration consists of 32,090 hexahedral elements for the fluid domain and 16,696 hexahedral elements for the solid domain. The meshes for the structure and fluid domains are conforming. The inner inlet diameter is 6.68 mm, the wall thickness is constant of 1.37 mm and the nine material constants C₁-C₉ from the Eq. 9 are stated in Table 1. The boundary conditions for the fluid domain are explained in Fig. 4. At the inlet, pulsatile flow is prescribed varying in time according to the measured data profile. The outlets of the model are connected to the inlet of the 1D model.

Table 1: The C1-C9 material constants.
	C1:	1.7x10^5	C2:	0.0	C3:	1.0x10^5
	C4:	0.0	C5:	1.0x10^5	C6:	0.0
	C7:	0.0	C8:	0.0	C9:	0.0

According to the Eq. 29, the number of generations of the 1D and 0D models depends on the initial diameters of the two outlets of the 3D model, namely the outlets of the ECA and ICA arteries. The resulting number of generations of the 1D model connected to the ICA is 16 and the consequent number of the STI generations is 13. As the outlet diameter of the ECA is smaller, also the number of the generations of the 1D model gets smaller, i.e. 14. The consequent number of the STI generations connected to the 1D model connected to the ECA outlet of the 3D model is 13. The times step used is 1/1024 sec, the μ of the incompressible Newtonian fluid of the 3D model is 4.06 × 10⁻³ m² · s⁻¹ and the Reynolds maximum and minimum numbers of the pulsatile inflow velocity profile are 1000 and 300, respectively. All the inlet and outlets structure nodes are fixed in all directions, the main point of the interest is the bifurcation area of the 3D model where the atherosclerosis occurs mostly. No-slip condition is prescribed on the FSI boundary.

(a) 3D FSI-1D-0D		(b) 3D Fluid-1D-0D
Figure 5: The WSS at the peak velocity of the 3D FSI-1D-0D (a) and the 3D Fluid-1D-0D (b) simulations.

The Fig. 5(a) shows the WSS distribution of the 3D FSI model coupled with the 1D-0D model. The values are extracted at the peak systolic velocity of 1.4 m·s⁻¹ in the center of the lumen. To compare the effect of the deformable walls, the Fig. 5(b) shows the results of the 3D Fluid analysis under the same simulation conditions as of the 3D FSI-1D-0D analysis. It can be observed that, due to the wall expansion under the pressure induced by the 1D-0D model, the WSS values of the 3D FSI-1D-0D model are smaller. As can also be observed in the graphs of the Fig. 6.

(a) 3D FSI-1D-0D

(b) 3D Fluid-1D-0D

Figure 6: The Flowrate profiles of 3D FSI-1D-0D (a) and the 3D Fluid-1D-0D (b) simulations.

The flowrates of the ICA and ECA at the peak systole in the FSI case are smaller than those in the fluid analysis only. On the other side, at the diastole, the flowrates are higher in case of the FSI analysis. This is the effect of the vessel walls expanding and contracting and thus decreasing and increasing the resulting flowrates, respectively. Therefore, at the diastole where the flowrates are increased due to the contraction of the vessel walls, the WSS distribution would show higher values for the FSI analysis, compared to the fluid analysis only. Consequently, when compared averaged over the whole cardiac cycle during which the increased values of the WSS at the peak systole even out with the decreased values of the WSS at the diastole, the time averaged WSS (TAWSS) distribution does not show significant difference between the two analyses, as shown in Fig. 7.

(a) 3D FSI-1D-0D		(b) 3D Fluid-1D-0D
Figure 7: The TAWSS of 3D FSI-1D-0D (a) and the 3D Fluid-1D-0D (b) simulations.

However, it is reasonable to expect that the wall dynamics would influence the direction of the WSS vector during the cardiac cycle, which is represented by the oscillatory shear index (OSI). Fig. 8 compares the OSI distribution, where significant difference can be observed as expected. It is generally accepted, that the regions with low WSS and high OSI values are the locations with the highest probability for the occurrence of the atherosclerosis.

(a) 3D FSI-1D-0D		(b) 3D Fluid-1D-0D
Figure 8: The OSI of 3D FSI-1D-0D (a) and the 3D Fluid-1D-0D (b) simulations.

FSI analysis provides additional information on the vessel wall dynamics, such as maximal principal strains and stresses, as shown in Fig. 9. The results are extracted at the peak systole. The strains show there is maximum deformation of 1.4% close to the bifurcation, slightly off the area with the highest OSI values. The principal stresses have their maximum values in the same area of the highest OSI values.

(a) principal strains	(b) principal stresses
Figure 9: The maximal principal strains (a) and stresses (b) of the structure domain.

4. Conclusion

Strongly coupled 3D FSI analysis with high-order Mooney-Rivlin hyper-elastic material is used. 1D-0D model incorporating the hematocrit effect is developed. The feasibility of coupling the strongly coupled FSI analysis with the 1D-0D model is shown. Compared to the fluid analysis only, some hemodynamic factors can be overestimated if the vessel wall deformations are not accounted for. The expansion and contraction of the blood vessel walls during the cardiac cycle change the direction of the WSS vector causing the values of the OSI to decrease. FSI analysis provides additional information about the vessel wall dynamics on top of the hemodynamic factors provided by the fluid part of the analysis. The maximal principal strains and stresses have their maximum values around the same region where also the highest values of OSI occur together with the lowest values of the WSS.

5. Acknowledgements

This research is supported by the Research and Development of the Next-Generation Integrated Simulation of Living Matter, as a part of the Development and Use of the Next-Generation Supercomputer Project of the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan.

6. References

	[1]	M.H. Friendman, G.M. Hutchins, C.B. Bargeron, O.J. Deters, and F.F. Mark, "Correlation between intimal thickness and fluid shear in human arteries," Atherosclerosis, vol. 39, no. 3, pp. 425-436, 1981.
	[2]	C.K. Zarins, D.P. Gidens, B.K. Bharadvaj, V.S. Sottiurai, R.F. Mabon, and S. Glagov, "Carotid bifurcation atherosclerosis quantitative correlation of plaque localization with flow velocity profiles and wall shear-stress," Circulation Research, vol. 53, no. 4, pp. 502-514, 1983.
	[3]	J.J. Yeung, H.J. Kim, T.A. Abbruzzese, I.E. Vignon-Clementel, M.T. Draney-Blomme, K.K. Yeung, I. Perkash, R.J. Herfkens, C.A. Taylor, and R.L. Dalman, "Aortic hemodynamic and morphologic adaptation to chronic spinal cord injury," Journal of Vascular Surgery, vol. 44, no. 6, pp. 1254-1265, 2006.
	[4]	J.D. Humphrey and C.A. Taylor, "Intracranial and abdominal aortic aneurysms: similarities, differences, and need for a new class of computational models," Anual Review of Biomedical Engineering, vol. 10, no. 1, pp. 221-246, 2008.
	[5]	P. Causin, J.F. Gerbeau, and F. Nobile, "Added-mass effect in the design of partiotioned algorithms for fluid-structure problems," Comput. Meth. Appl. Mech. Eng., vol. 194, pp. 4506-4527, 2005.
	[6]	C. Foerster, W.A. Wall, and E. Ramm, "Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows," Comput. Meth. Appl. Mech. Eng., vol. 196, pp. 1278-1293, 2007.
	[7]	M.A. Fernandez, J.F. Gerbeau, and C. Grandmont, "A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid," Int. J. Numer. Meth. Eng., vol. 69, pp. 794-821, 2007.
	[8]	M.S. Olufsen, C.S. Peskin, W.Y. Kim, E.M. Perdersen, A. Nadim, and J. Larsen, "Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions," Annals of Biomedical Engineering, vol. 28, pp. 1281-1299, 2000.
	[9]	G.A. Holzapfel, Nonlinear solid mechanics, a continuum approach for engineering, John Wiley, 2000.
	[10]	A. Huerta and W.K. Liu, "Viscous flow with large free surface motion," Comput. Meth. Appl. Mech. Eng., vol. 69, pp. 277-324, 1988.
	[11]	M.S. Pollanen, "Dimensional optimization at different levels at the arterial hierarchy," J. of Theoretical Biology, vol. 159, pp. 267-270, 92.
	[12]	P. Segers, F. Dubois, D. de Wachter, and P. Verdonck, "Role and relevancy of a cardiovascular simulator," Cardiovascular Engineering, vol. 3, pp. 48-56, 1998.
	[13]	N. Stergiopulos, D.F. Young, and T.R. Rogge, "Computer simulation of arterial flow with applications to arterial and aortic stenosis," Journal of Biomechanics, vol. 25, pp. 1477-1488, 1992.
	[14]	R.A. Freitas Jr., Nanomedicine, Volume I: Basic Capabilities, Landes Bioscience, 1999.
	[15]	A.R. Pries, T.W. Secomb, T. Gebner, M.B. Sperandio, J.F. Gross, and P. Gaehtgens, "Resistance to blood flow in microvessels in vivo," Circulation Research, vol. 75, pp. 904-915, 1994.