INVESTIGATION OF THE BLOOD FLOW IN DEFORMABLE VESSELS USING STABILIZED FINITE ELEMENT METHOD INVESTIGAÇÃO DO ESCOAMENTO SANGUÍNEO EM VASOS

The study and simulation of blood flow in the cardiovascular system have many applications such as pathologies studies, surgical planning, and design of medical devices. Several works within this area consider the problem using a rigid wall assumption, while blood velocity and pressure in large arteries are greatly influenced by vessel wall dynamics. The Coupled Momentum Method (CMM) is based on a strong coupling of degrees-of-freedom of the fluid and the solid domains. This coupling is made by considering that the deformation of the wall, in a variational level, become a boundary condition for the fluid domain. As an advantage, description of motion (Eulerian) may be kept the same and a fixed mesh. In this study, blood is considered a Newtonian fluid and the wall a thin-walled linear elastic. This work is focused on using fluid-structure interaction in 3D geometries to evaluate the blood and vessel dynamics in contrast to the rigid wall formulation what is considered only a fluid dynamic problem. For realistic and physiological parameters, CMM has demonstrated to be a good alternative for the simulation of blood flow in large arteries by virtue of representing more realistic phenomena than a rigid wall model. The results were obtained using FEniCS and Python, and are in agreement with the theoretical and numerical solutions from the literature.


INTRODUCTION
Computational Fluid Dynamics (CFD) is an area of mathematics and branch of fluid mechanics. It is used in the design of many systems and can be used for modeling blood flow (MORRIS et al., 2016). Computational Fluid Dynamics (CFD) is an area of computational science engineering and a branch of fluid mechanics used in the design of many physical systems with great potential in medicine. In this context, Navier-Stokes equations (NSE), which describes the motion of viscous fluid particles based on the conservation of mass and momentum, could be used by the modeling of the blood flow (MORRIS et al., 2016). In some applications, the fluid exerts interactions with another medium, and vice versa, into so-called fluid-structure interaction (FSI) problems (BAZILEVS; TAKIZAWA; TEZDUYAR, 2013).
As described in (BAZILEVS; TAKIZAWA; TEZDUYAR, 2013), FSI problems are present in many fields. For the biomechanics area, works (VOSSE; DONGEN, 1998;FIGUEROA et al., 2006;FIGUEROA, 2009;MORRIS et al., 2016;THOMAS;SUMAM, 2016) describe the importance of FSI in this field of study. Those kinds of problems characterize a set of partial differential equations and boundary conditions related to each of the domains (fluid and solid). They form a single coupled set of non-linear equations which make it difficult to search for analytical solutions (BAZILEVS; TAKIZAWA; TEZDUYAR, 2013). Significant advances in research of computational methods have been developed in the last decades, especially applications in problems of fluid mechanics and FSI (KUHL; HULSHOFF; BORST, 2003;FIGUEROA et al., 2006;BAZILEVS;TAKIZAWA;TEZDUYAR, 2013;COURT;FOURNIÉ, 2015).
One of the most used approaches to include deformation of the vessel wall is the Lagrangian-Eulerian-Arbitrary (ALE) formulation to solve the FSI problems (DONEA S. GIULIANI, 1982;HUGHES, 1987). However, as described in (FIGUEROA et al., 2006), this technique encounters significant challenges for three-dimensional blood flow modeling, it needs to update the geometry of the fluid and structure at each time step, resulting in an expensive approach. Also, some physiological applications require the models to be a time interval of relevance. Therefore, the Coupled Momentum Method for Fluid-Structure Interaction (CMM-FSI) proposed by Figueroa et al. (2006), was made thinking in avoid some difficulties found in the ALE method for 3D geometries. The CMM considers the conventional formulation of finite elements for the Navier-Stokes equations in a rigid domain, adding modifications in such a way that the deformation of the membrane is considered. The idea is based on the fact that the diameters of the arteries are small compared to the wavelength of the cardiac pulse. It causes the solid to exhibit a membrane over the bend, requiring no additional degrees of freedom in the wall of the solid. Thus, the linear membrane is reinforced Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 294-01, 294-15, 2020 DOI: 10.21575/25254782rmetg2020vol5n61387 with transverse shear discussed in (COOK et al., 2007;HUGHES, 2012).
There are studies about blood flow on arteries consider a fixed wall approach (VALEN-SENDSTAD et al., 2011;EVJU;VALEN-SENDSTAD;MARDAL, 2013). This work is focused on showing how numerical simulations for blood flow in rigid walls may be impacted in contrast to models that considers the vessel wall dynamics. This investigation will be done by applying the CMM-FSI with physiological parameters and idealized geometries.

Blood Flow Equations
Blood flow in the large vessels can be approximated as an incompressible flow for a Newtonian fluid, described by the NSE. The assumption of a Newtonian fluid requires the viscous stresses to be linear functions of the components of the strain-rate tensor. Incompressible means no change in specific mass over time. Since mass can neither be destroyed nor created, the principle of mass conservation, that expresses local conservation of mass at any point in a continuous medium is given by Eq. (1). Besides mass conservation a fluid also obeys conservation of momentum (i.e., Newton's second law). The equation for momentum balance is given by Eq. (2). Therefore, for a Newtonian fluid, there are two kinds of stresses acting in a fluid particle: internal pressure (σ 1 ) and due to the pressure and viscous stresses (σ 2 ), given by Eq. (3).ρ + ∇.(ρu) = 0, where u is the blood velocity vector, ρ is the blood specific mass, p is the blood pressure, σ is the total stress tensor acting in the blood particle, f is external forces acting in the fluid, I is the d × d identity tensor and µ is the fluid dynamic viscosity. The superscript˙means the time derivative d dt .
Writing the conservation laws, together with constitutive relationships for the fluid, and the initial and boundary conditions for the problem it is possible the define mathematically the Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 294-01, 294-15, 2020 DOI: 10.21575/25254782rmetg2020vol5n61387 4 blood flow problem: where Ω is the problem domain, Γ g are the Dirichlet boundaries, Γ h are the Neumann boundaries and Ω 0 is the initial condition of domain, n is the normal vector, g is the prescribed Dirichlet condition, h is the prescribed Neumann condition and finally, u 0 and p 0 are the initial velocity and pressure respectively.

Discretization of Blood Flow Equations
We start defining the trial solution and weighting function spaces for the semi-discrete formulation of the problem given by Eq. (4)): where Ω represents the spatial domain, H 1 (Ω) represents the usual Sobolev space of functions with square-integrable values and first derivatives in Ω. P k (Ω e ) is the local polynomial approximation space. For this work, a stabilized formulation of the semi-discrete Galerkin finite element formulation is used for solving the NSE. After integrating by parts and rearranging term, we can obtain the SUPG/PSPG/LSIC formulation as follow: Find u h ∈ U h and p h ∈ P h , in such way that: ∀w h ∈ W h and q h ∈ P h . Following the methods proposed by Hughes and colleagues, described in (TAYLOR; HUGHES; ZARINS, 1998), the stabilization parameters τ pspg (PSPG), τ supg (SUPG) and τ lsic (LSIC) are given by: where R is the residual of momentum (Eq. 2) and h e is the local mesh size.

Vessel Wall Equations
The classic elastodynamics equation used to describe the motion of the vessel wall in a domain Ω s is given by: where v is the displacement vector,v the second time derivative of the displacement field, ρ s is the wall density, σ s is the wall stress tensor, which depends on the constitutive relationship, and F v represents the external forces.
The vessel walls are herein approximated using linear elastic and thin-walled structure assumption, so the derivation of σ s is graphically represented in Fig. 1.
The consideration of a thin-walled structure (membrane) model for the vessel wall nodes will enable a strong coupling of degrees-of-freedom of the fluid and solid domains and will result in an expression for the unknown integral traction. This approximation is justified by experimental evidence showing that the vessel wall constitutive behavior can be Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 294-01, 294-15, 2020 DOI: 10.21575/25254782rmetg2020vol5n61387 reasonably assumed as linear within the physiological range of pressures (ZHOU; FUNG, 1997). Within a thin membrane assumption, we can neglect variations across the thickness and using symmetries, it is possible to work with a reduced form of the stress tensorσ s and material constants, given by: where E and ν are Young's modulus and Poisson's ratio coefficients respectively. For a solid, homogeneous plate, the parameter k = 5/6 accounts for a parabolic variation of transverse shear stress through the membrane in these works (FIGUEROA et al., 2004;COOK et al., 2007).

Discretization of Wall Vessel Equations
In order to relate the solid problem with the lateral boundary of the fluid domain, Ω s needs to be mapped on Γ s . Assuming a thin-walled structure, integrals defined on Γ s h and Ω s can be related with integrals on the lateral boundary of the fluid domain Γ s according to the following expressions (FIGUEROA et al., 2006): The surface traction t f acting on the fluid lateral boundary due to the interaction with the solid is equal and opposite to the surface traction t s acting on the vessel wall due to the fluid, and since the wall is a thin-enough structure, the internal surface traction t s will be felt uniformly through the wall thickness ζ as follows: So it is possible to relate the unknown integral term Γs −w.t f dΓ of Eq. (8) with the weak form for the vessel wall problem given by Eq. (16).
Since a strong coupling of the degrees-of-freedom of the fluid and solid domains is considered, the displacement, velocity and acceleration fields on the fluid-solid interface are identical. Furthermore, since the weak form of the elastodynamics equations has the same differentiability requirements on the functional spaces as the fluid weak form, we can adopt for the vessel wall problem the same type of functional spaces as in the fluid domain, given by Eqs. (5) and (6). By that way, using the body force b s of Eq. (16), the expression of the weak form for the solid domain using the semi-discrete Galerkin finite element formulation is: Find for all w ∈ W h v . The acceleration term (u ,t ) has been written as the time derivative of the velocity (u) rather than (v) since the goal is to express the vessel wall equations in terms of the fluid unknown parameters.

RESULTS
The CMM-FSI formulation starts from a conventional stabilized finite element for the NSE in a rigid domain and modifies it in such a way that the deformation of the wall domain surrounding the fluid is taken into account. In other words, the CMM-FSI we will solve the NSE with a Neumann boundary condition of traction in the fluid-structure interface.
In this method, the acceleration termv will be written as the time derivative of the velocity (u) rather than as the second time derivative of the displacement field (v) since the goal is to express the vessel wall equations in terms of the unknown fluid parameters. Therefore, the wall is fixed by constraining the degrees-of-freedom only in the nodes located at the inlet and outlet rings. Thus, all other nodes of the wall are allowed to move in any direction. Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 294-01, 294-15, 2020 DOI: 10.21575/25254782rmetg2020vol5n61387 The results were performed using real and physiological parameters (FIGUEROA et al., 2006), and can be seen in Table 1. µ f E ρ s ν s k ζ g/cm 3 dyn.s/cm 2 dyn/cm 2 g/cm 3 --cm 1.06 0.04 4.07 × 10 6 1.09 0.5 5/6 0.03 Source: (FIGUEROA et al., 2006).

Validation of the Implemented Fluid Structure Model
The behavior of fluid structure model will be validated in comparison with the blood flow theory by Poiseuille: where p 0 = 1413 dyn/cm 2 , c = −19.76 dyn/cm 3 is the axial pressure drop, which is related to flow rate, viscosity and pipe radius, R is the artery radius, r is the radial distance and V max = 12.5 cm/s is the maximum velocity found at the inlet centerline.
A simple cylindrical model with nominal radius and vessel length of 0.3 cm and 4.2 cm, respectively, is used within a mesh of 24067 elements and 4782 nodes. For the inlet, it is fixed a steady-state flow mapped to a parabolic velocity profile. For outlet, it is prescribed a constant pressure of 1330 dyn/cm 2 . The numerical simulation was performed until the steady-state condition in order to compare with Poiseuille equation (PRITCHARD; MITCHELL, 2016). Once reached, the wall must have such small displacement that does not change the flow profile and pressure decay.  Source: The authors.
Finally, Fig. 3 represents the wall displacement and the blood flow streamlines for the initial and final phases of the simulation. The radial displacement after reaching the steadystate condition must be next to zero, and for this study, a radial displacement of 0.0005 cm was found, as shown in Fig. 3b. Moreover, a maximum radial displacement of 0.0017 cm is found, equivalent to a 1.13% deformation. These results are in agreement with such that a maximum deformation of 5% is obtained with a physiologic range of pressures (MILNOR, 1989).

Rigid Versus Deformable Model
For a more in-depth comparison between rigid and deformable wall models, this section presents a flow study in an idealized model of the carotid artery using pulsatile flow waves. Pulsatile inflow condition is used in reason to observe the wave propagation and difference between pressure spikes phenomena in both models. For this purpose, an idealized carotid Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 294-01, 294-15, 2020 DOI: 10.21575/25254782rmetg2020vol5n61387 artery, which is a simple cylindrical model with nominal radius and vessel length of 0.3 cm and 12.6 cm respectively, is used. A mesh with 38897 elements and 8178 nodes was made to perform the simulations. Figure 4 shows the sections (S1 and S2), which were used for evaluate the comparison. These locations were chosen to be far enough from the inlet and outlet. Thus, the radial deformation is not affected by the presence of the fixed boundaries. For the inlet, a pulsatile periodic flow wave (Fig. 5a) is prescribed, mapped to a parabolic velocity profile. For outlet, it is used pressure 100 mmHg (constant) and a time dependent pressure wave (Fig. 5b).   Figures 6 and 7 shows that the pressure field and flow waves have discernible solution for each case. For the pressure field, Fig. 7b, display that pressures obtained with rigid wall models present a higher pressure pulse value. As described in (VLACHOPOULOS; O'ROURKE; NICHOLS, 2011), this pressure difference of ∆p = 145.18 − 134.53 = 10.65 mmHg is a well-known phenomenon in the cardiovascular system, where stiffer vessels tend to experience higher pressure pulses, showing higher values in systoles (peaks) and lower in diastole (valley). On the other hand, Fig. 7a does not present this expected result. Figueroa et al. (2006) explain that some outflow boundary conditions, as probably this case, generates unrealistic high pressure pulses in deformable models.
The velocity field also presents noticeable differences. First of all, it is possible to notice a phase lag of 0.015 s between the rigid and deformable waves in Fig. 6a  wave propagation phenomena. This phase lag provides a means to estimate the pulse velocity (MILNOR, 1989). Considering that the vessel length is 12.6 cm, it produces a pulse velocity of approximately 840 cm/s, in agreement with the speed range discussed in (VLACHOPOULOS; O'ROURKE; NICHOLS, 2011).
Finally, Fig. 8 shows that the displacement have a totally different behavior for each model, for both dynamics and its magnitude. Moreover, negative values are shown due the referential and the maximum radial displacement is 0.007 cm and it is equivalent to 2.33% deformation. Figure 8: Comparisons between rigid and deformable models for radial displacement overtime in sections S1 and S2 with distinct outlet pressure functions.
Source: The authors.

CONCLUSIONS
This work has successfully implemented a fluid-structure interaction method for simulating blood flow in 3D deformable domains. We have obtained satisfactory results when comparing with a theoretical case and also shown that the assumptions of linear elasticity and thin-walled structure were appropriated for physiological parameters.
The simulation of blood flow using an idealized artery, i.e., a simple pipe channel was performed until the steady-state, and the results were very similar to both rigid case and theoretical values. Thus, an idealized model of the common carotid artery was used to conduct a comparative study between rigid and deformable wall assumptions. For this study, real and physiological blood flow and pressure waves were used to represent a more realistic simulation. The results have demonstrated that for 3D geometric models, each model choice brings a different result, which should be analyzed. Therefore, since the pressure differences of both models are quite large, it is important to notice that even using a model that represents more Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 6, p. 294-01, 294-15, 2020 DOI: 10.21575/25254782rmetg2020vol5n61387 realistic phenomena, the selection of physiological boundary conditions is mandatory when seeking clinical studies.
For future works, the use of the impedance boundary conditions for outflow is an idea, due to the wave propagation phenomena be directly linked to this condition. For the study of pathologies as the stenosis and aneurysm, the use of non-Newtonian rheology models is also an option because they may represent phenomena closer to physiology. Thus, to perform a patientspecific study, using a 3D model of some artery reconstructed from clinical images.