*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. FluidStructure 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 hyperelastic structure, representing the vascular tissue. A high order MooneyRivlin 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 FluidStructure finite element model and the multiscale 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 patientspecific models, such as the deformations, principal strains and stresses. The traditional FSI techniques are based on the arbitrary LagrangianEulerian (ALE) method. In this study, the strongly coupled ALE method is used to address the addedmass effect, which is typical for loosely coupled timestepping 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 multiscale model, here by coupling the 3D FSI model with the 1DSTI model. The 3D FSI domain is coupled to the inflow of the reducedorder 1D0D 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.
 
 Figure 1: FluidStructure interaction scheme. 
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 NavierStokes equations.
  (1) 
  (2) 
  (3) 
  (4) 
The Eqs. 14 are continuity equation, NavierStokes 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 NavierStokes 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.
 
 Figure 2: Tonti diagram for the governing equations of the structure domain. 
The
C in the Figure 2, is the right CauchyGreen tensor, the
F is the deformation gradient tensor, the
E is the GreenLagrange strain tensor, the
S is the second PiolaKirchhoff stress tensor for hyperelastic material using MooneyRivlin model with highorder 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_{1}
C_{9} are derived to fit the stressstrain curve of the material to be modeled and the
I_{1} and
I_{2} are invariants of the right CauchyGreen 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 PiolaKirchhoff stress tensor can be written, in
terms of the right CauchyGreen 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 CauchyGreen 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_{3} = 1, because det
C = (det
F)
^{2} = 1. Now, considering Eq. 13, Eq. 12 can be rewritten as
  (14) 
It can be readily verified that
Considering the above derivations, the second PiolaKirchhoff 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 ALEFEM 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 predictormulticorrector algorithm (PMA) [10] based on the Newmarkβ method for time stepping. The resulting nonlinear systems are solved using a NewtonRaphson 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+∆t}F^{(k)} −
^{t}Q^{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. 2022
 Begin NewtonRaphson loop
 Solve the coupled system, Eq. 24
 Form correctors, Eqs. 2628
 End NewtonRaphson 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 NewtonRaphson 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 (1D0D 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 crosssectional area,
A_{0} is the initial area,
p is the pressure,
ρ is the density,
E is Young’s modulus,
h is the wall thickness and
r_{0} 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_{0} according to the formula based on compliance estimates [12,13] as follows
  (33) 
The
k_{1},
k_{2} and
k_{3} are constants of following values and units,
k_{1} = 2.0×10
^{6} kg·sec^{−2}·m^{−1},
k_{2} = −2.253×10
^{3} m^{−1} and
k_{3} = 8.65×10
^{4} kg·sec^{−2}·m^{−1} [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 semianalytical 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 wavepropagation 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 nonNewtonian 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
_{1}C
_{9} 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.
 
 Figure 3: The undeformed configuration of the carotid artery model. 
Table 1: The C_{1}C_{9} material constants.


C_{1}: 
1.7x10^5 
C_{2}: 
0.0 
C_{3}: 
1.0x10^5 


C_{4}: 
0.0 
C_{5}: 
1.0x10^5 
C_{6}: 
0.0 


C_{7}: 
0.0 
C_{8}: 
0.0 
C_{9}: 
0.0 

 
 Figure 4: The boundary conditions applied to the 3D models. 
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
^{−3} m^{2} · s^{−1} 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. Noslip condition is prescribed on the FSI boundary.
The Fig. 5(a) shows the WSS distribution of the 3D FSI model coupled with the 1D0D model. The values are extracted at the peak systolic velocity of 1.4
m·s^{−1} 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 FSI1D0D analysis. It can be observed that, due to the wall expansion under the pressure induced by the 1D0D model, the WSS values of the 3D FSI1D0D model are smaller. As can also be observed in the graphs of the Fig. 6.

(a) 3D FSI1D0D 

(b) 3D Fluid1D0D 
Figure 6: The Flowrate profiles of 3D FSI1D0D (a) and the 3D Fluid1D0D (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.
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.
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 highorder MooneyRivlin hyperelastic material is used. 1D0D model incorporating the hematocrit effect is developed. The feasibility of coupling the strongly coupled FSI analysis with the 1D0D 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 NextGeneration Integrated Simulation of Living Matter, as a part of the Development and Use of the NextGeneration 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 shearstress,” Circulation Research, vol. 53, no. 4, pp. 502–514, 1983.

 [3] 
J.J. Yeung, H.J. Kim, T.A. Abbruzzese, I.E. VignonClementel, M.T. DraneyBlomme, 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, “Addedmass effect in the design of partiotioned algorithms for fluidstructure 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 semiimplicit 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 structuredtree 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.
