European Journal of Mechanics / B Fluids 77 (2019) 273–280

Contents lists available at ScienceDirect

European Journal of Mechanics / B Fluids

journal homepage: www.elsevier.com/locate/ejmflu

The protocol for using elastic wall model in modeling blood flow
within human artery
∗ ∗
Marcin Nowak a , , Bartłomiej Melka a , , Marek Rojczyk a , Maria Gracka a ,
∗
Andrzej J. Nowak a , Adam Golda b , Wojciech P. Adamczyk a , , Benjamin Isaac c ,
∗
Ryszard A. Białecki a , Ziemowit Ostrowski a ,
a
Silesian University of Technology, Institute of Thermal Technology, Biomedical Engineering Lab, Konarskiego 22, Gliwice 44–100, Poland
b
Gliwice Medical Centre, Department of Cardiology at 4th Municipal Hospital, Kosciuszki 29, 44–100 Gliwice, Poland
c
Institute for Clean and Secure Energy, University of Utah, UT, USA

highlights

• Numerical model of female aorta with congenital disease.

• Application of the iterative fluid–structure interaction technique for predicting elastic nature of the human artery.
• Developed model allows to predict vessel wall displacement in stable way during the entire heart cycle.

article info a b s t r a c t

Article history: Medical diagnostic tools will play a major role in the future for an effective patient treatment and
Received 5 June 2018 reduction their mortality, related to the cardiovascular system diseases (CVDs). There is an urgent need
Received in revised form 20 February 2019 for developing diagnostic procedure to be robust, reliable, accurate and efficient, in the framework of
Accepted 26 March 2019
a paradigm shift. Application of numerical techniques is seen as a perspective tool for such purpose.
Available online 4 April 2019
Nevertheless, existing models need constant improvement in development robust, multi-scale models.
Keywords: This paper elaborates on the development of numerical model for modeling blood flow in the aorta
FSI section. The deformation of the blood vessel was modeled as two-way fluid–structure interaction (FSI)
Blood flow using ANSYS package. Numerical results have shown that the developed model predicts deformations
CFD of the vessels and describes their impact on the pressure, pressure drop and wall shear stresses
Non-Newtonian flow distributions. Differences between rigid and deformed walls were checked based on pressure drop
Elastic blood vessels value. For movable walls, these values were higher both for systole and diastole, which is caused by the
local wall compression during aforementioned moments of the cycle. The significant backflow observed
during the heart cycle is connected with the deformed walls resulting in temporal blood accumulation.
The maximum total deformation of the vessel walls achieved 2.35 mm, and the difference between
the maximum and minimum blood volume was equal 5.2%.
© 2019 Elsevier Masson SAS. All rights reserved.

1. Introduction and its pathological changes, such as atherosclerosis, aneurysm,

strokes, and heart attacks. CVD treatment is one of the biggest
In 2015, cardiovascular diseases (CVDs) caused 31% of deaths expenses on national healthcare and it encompasses cases of both
[1] and they remain the leading cause of patients mortality. Each under- and over-treatment.
year, CVDs cause 45% of all deaths in Europe [2,3]. The main Thus, there is room for the enhancement of diagnostic proce-
causes of these deaths are diseases of the cardiovascular system dures leading to a better outcome of patient treatment, without
increasing its costs. The present paper belongs to the vascular
∗ Corresponding authors. hemodynamics, focusing on flow within the thoracic aorta with
E-mail addresses: marcin.nowak@polsl.pl (M. Nowak), its major branches. The displacement of the vessel walls due to
Bartlomiej.Melka@polsl.pl (B. Melka), wojciech.adamczyk@polsl.pl forces acting on its wall, and resulting from pressure and shear
(W.P. Adamczyk), Ziemowit.Ostrowski@polsl.pl (Z. Ostrowski).
associated with the blood flow are taken into account. Blood
URLs: http://www.itc.polsl.pl/bmelka (B. Melka),
http://www.itc.polsl.pl/wadamczyk (W.P. Adamczyk), vessels are elastic and their diameter changes within the cardiac
http://www.itc.polsl.pl/ostrowski (Z. Ostrowski). cycle profile. This feature smooths the blood flow rate within the

https://doi.org/10.1016/j.euromechflu.2019.03.009
0997-7546/© 2019 Elsevier Masson SAS. All rights reserved.
274 M. Nowak, B. Melka, M. Rojczyk et al. / European Journal of Mechanics / B Fluids 77 (2019) 273–280

heart cycle, allowing for accumulation of blood during systole and forward pressure waves, traveling from the heart, and backward
releasing it during diastole [4]. pressure waves, generated as the reflection along the cardiovas-
When the blood flow within deformable vessels is modeled, cular system and leading to increased systolic blood pressure [17].
the question arises, whether the wall can be treated as a rigid A pulse wave analysis was performed on a repaired aortic coarc-
boundary or the fluid–structure interaction approach should be tation. Even a successful treatment introduce scar tissue or stent,
used. Rigid walls assumption overestimates the wall shear stress which leads to local stiffening. This stiffening affects the pulse
(WSS), sometimes by 50%, and it causes some qualitative and pressure and heart load, as was shown using a simplified geomet-
quantitative differences when compared with more accurate flex- rical model of a 3D flexible tube. In the area where the CoA has
ible wall approach [5]. Moreover, the blood vessels of children are been removed, additional stiffness was introduced by imposing a
more flexible than those of adults [5], whose additional stiffness six times higher Young modulus than in the remaining portion of
comes from arteriosclerosis. Thus, models of blood flow when the arterial wall. The influence of the time step and grid sizes on
applied to younger patients in particular should account for the the amplitude of the traveling wave was investigated.
cyclic changes of the vessels shape, i.e. fluid–structure interaction In [18], the FSI technique was used to model eight mice
(FSI) should be implemented. Our research uses this approach to carotid left and right bifurcations, obtained from four mice, where
simulate a case of congenital aortic stenosis — coarctation. The animal-specific geometries and in vivo measured boundary con-
Coarctation of the aorta (CoA) is the narrowing of the artery, ditions were used. The conclusion drawn from these investi-
which results in differential blood flow, pressure loss across the gations was that the FSI approach is required to capture the
CoA, turbulent jet flows and causes hypertension in the upper influence of the shear stress in the development of atheroscle-
body circulation [6]. According to different sources, CoA accounts rosis.
for 5%–11% of all congenital heart defects [7–9] and is 1.7 times The short review of the subject indicates that the accurate
more common in boys than girls [7]. In rare circumstances, it prediction of the blood vessels wall displacement and its move-
develops later in life. It may occur as an isolated defect or it ment has a significant impact on the predicted flow pattern
and it cannot simply be ignored. The approach presented in this
may develop in association with other lesions. To fix the con-
publication is focused on the application of the FSI to model
dition, surgical treatment, where the narrowing can be simply
blood flow in the main human artery with coarctation. The au-
removed, is required. Depending on the size of the CoA, the
thors’ earlier research [19] indicates that the displacement of
Dacron graft [10] may be implemented and two free ends are
the artery wall can significant affect the flow pattern. These
reconnected. In most cases, angioplasty can be performed instead
simulations were carried out using a single phase model, where
of using intrusive and difficult surgery. This paper is devoted
the vessel walls were treated as rigid, and the response of the
to the development of simulation model that can be used to
downstream cardiovascular system was captured by the lumped
effectively simulate the blood flow within deformable arteries,
parameter model [20]. This model is based on electrical analogy
including cases where the CoA condition appears.
using three-element-Windkessel and is used to calculate proper
Several papers on application the FSI analysis to model the outlet pressure for the boundary conditions.
aortic blood flow have been published. A comparison between There are some investigations focused on the usage of CFD in
the rigid and FSI model based on the modeling of a human left the CoA treatment, noted [19,21–24], but none of them include
coronary artery bifurcation was presented in [11]. It was shown the FSI model of aorta with this disease. The ascending aorta has
that the overall flow field and displacement distribution is consid- the ability of storing blood during systole and deliver it to the
erably affected by omitting the influence of myocardial motion. rest of the circulatory system. However, this ability is disordered
Other applications of FSI approach for modeling pulsatile flow in the CoA-affected aortas [21]. Thus, the FSI model of these
within aortic aneurysm can be found [12]. In reference [13], an geometries should be used. The numerical studies shown that in
FSI analysis was performed on the blood flow through a Dacron the subjects affected by CoA the higher time-averaged wall shear
graft, which replaced the aortic dissection segment. A simplified stress and oscillatory-shear index is observed. The earlier authors’
cylindrical geometry model was used and the shell model was paper [19] showed the influence of the virtual therapy onto the
assumed for the mechanical solver. The results showed that the hemodynamic properties, which can be found as the Dacron graft
significant change in radial displacement in the interface between implementing.
two materials was not observed. The von Mises stress analysis When compared with its stiff wall counterpart, the application
on the junction between graft and aorta demonstrated that ef- of the FSI dramatically increases the computational time. To keep
fects such as the unboundedness of the layers or rupture are this time within acceptable limits, simplified, constant boundary
not probable due to the much lower stresses than the ultimate conditions at outlets, based on data published in [25,26], were
strength of the aorta. In [14] a numerical analysis of blood flow applied. This means that the lumped model, based on electrical
in a flexible abdominal real aorta was performed. The computed analogy, was not used in this study. The attempts of using Wind-
pressure in the FSI model was 15% lower compared to the rigid kessel model during FSI calculations resulted not only in the very
model. Moreover, these results compared well with the available excessive computational cost, but also in serious convergence
in vivo measurements. [15] created an FSI model of the blood problems. The displacement distribution and the changes of the
flow in a real artery, treating blood as a Newtonian fluid. A blood domain volume are discussed in Section 3. As in [19], in this
monolithic solver was used here. This kind of numerical pro- paper the hemodynamics was represented using the single-phase
cedure solves both fluid and solid equations in one matrix in model. The simulation model was implemented using the ANSYS
Arbitrary Lagrangian Eulerian (ALE) formulation, using finite ele- package, which combines both the ANSYS Fluent CFD solver and
ment method. The rigid wall assumption was concluded to show the ANSYS Mechanical APDL solver. The considered geometry can
some unphysiological results; that is, a flux entering through be classified as strongly deformable biological structure where
the inlet section at the moment when the valve is closed. This two-way implicit coupling procedure needs to be used.
abnormality was not observed in the FSI analysis. In [16] the
FSI analysis on similar geometry — human aorta was presented. 2. Material and methods
Large differences between the rigid wall model were reported,
especially in the case of the flow rate, which increased almost Geometry
100% when compared with the FSI. The geometry of an 8-year old female patient with moderate
The FSI analysis of the pulse propagation in arteries was re- thoracic CoA (∼65% area reduction) was used. Fig. 1 illustrates
alized by [17]. The measured blood pressure is the result of the the shape of the aorta, together with the marked branches. The
 M. Nowak, B. Melka, M. Rojczyk et al. / European Journal of Mechanics / B Fluids 77 (2019) 273–280 275

Fig. 2. Division of the artery model into nine bodies with marked thicknesses
for each of the section given in mm. The upper left corner shows a sample
connector merging Sections 2 and 1.

using 1,250,000 elements. To accurately represent the flow field

while keeping the computational time within reasonable limits,
the upper branches of the aorta were discretized using finer mesh
than the main aorta section. The sensitivity of the numerical
results on changing the mesh resolution was checked using three
different mesh resolutions: coarse (900,000), medium (1,250,000)
and fine (1,700,000). Carried out simulations has shown minor
Fig. 1. Geometry used for numerical simulation with marked branches.
influence of element size on predicted pressure distribution as
well as predicted vessel wall displacement. For each case the
mesh size was decreased as well as the prism boundary layer was
geometry was acquired during Gadolinium-enhanced magnetic densified, until the controlled variables during the cycle achieved
resonance angiography (MRA). During the MRA process, the pa- constant level and do not change via further densification. These
tient was in supine position with held breath [8]. variables were the local and average static pressure values. On
Base geometry was represented as a raw unstructured trian- the solid side, the FEM mesh was tested via imposing constant
gulated surface format (STL). The MRA method does not acquire pressure on the inner side of the wall, and it was densified
the geometry of the solid domain, due to that, the MRA figures until the total deformation stops varying apparently. The domain
need to be used initially to build the vessel’s external shape. To which represents the vessel walls was discretized using 90,000
simulate fluid flow within the blood vessel, the internal part of
elements. To model the cyclic distortion of the vessel, dynamic
the vessel needs to be filled by the fluid, which means that the STL
mesh with activated smoothing and remeshing options was used.
file needs to be converted into the fluid/solid body. The Geomagic
Both generated meshes are shown in Fig. 3.
Design X (3D Systems Corporation) was used for that purpose.
The most common assumption behind the FSI is the constant Mathematical model
thickness of the wall of the blood vessel in the entire domain. The main advantage of using the FSI methodology is the pos-
A thinner wall in the region of smaller vessels produces bigger sibility to capture displacement of the blood vessel by coupling
deformations, leading to higher flow rates. Several different ap- the fluid and the mechanical solver, i.e. ANSYS Fluent and ANSYS
proaches can be used to specify the vascular wall thickness. [14] Mechanical APDL. The fluid–solid interactions can be modeled
has proposed a very simple approach where constant thicknesses using a one-way or two-way coupling approach. In the first, the
were used for the abdominal and iliac artery, equal to 1 mm fluid influences the shape of the blood vessel geometry, while the
and 0.5 mm, respectively. A more universal approach has been changes in the shape or the solid does not affect the fluid motion.
presented in [27] where patient’s age and the diameter of the This assumption can be appropriate if the wall is stiff enough to
vessel were taken into account. In the present study, the thickness keep the low deformation caused by the fluid flow. In this case,
of the blood vessel has been defined as 10% of the local effective the information is passed from the CFD solver to the mechanical
vessel radius [5]. Using the Geomagic Design X code, the source solver and the stress analysis in the solid satisfies the assumption
STL geometry was first fixed and the resulting shape was then of small deformations. Because the Young modulus of the walls of
subdivided into nine segments, each of approximately constant the arteries is low, the vessels undergo substantial deformations
diameter. induced by the blood flow [5], which means that the geometry
The outer surface of each zone was then created using the Au- change enough to affect the flow field and deformations. To
tosurface or Loft operations of Geomagic package [28] depending capture this behavior, two-way coupling needs to be invoked.
on the domain regularity. In this study, we used the implicit iteration scheme, where the
After creating nine solids, they were merged into one by iterative data exchange between the solid and fluid solver is
connecting the adjacent subdomains. When linking segments carried out at every time step until convergence is reached. The
of different wall diameter, connectors with linearly changing data exchange between the coupled solvers is schematically de-
thicknesses have been introduced to guarantee smooth geometry picted in Fig. 4. The stop criterion within each time step was to
transition. The sections of the geometrical model with various keep the normalized change between transferred values (RMS)
vessel thickness and a sample connector are presented in Fig. 2. below 1%. However, this procedure was numerically expensive.
The thicknesses of two zones, 2 and 7, presented in Fig. 2 The calculation time was approximately 123 h, using 12 2.4 GHz
were calculated as an arithmetic mean of contiguous zones while Intel Xeon E5-2620 cores.
avoiding discontinuous changes of vessel thicknesses. In the en- The walls of the arteries consist of three layers: adventitia,
tire computational domain, the wall thickness was within the media and intima (which include endothelium and elastica in-
0.15 mm to 0.85 mm range. The entire domain was discretized terna) [29] each of different thickness and mechanical properties.
276 M. Nowak, B. Melka, M. Rojczyk et al. / European Journal of Mechanics / B Fluids 77 (2019) 273–280

Fig. 4. Flowchart of the 2-way iteratively implicit approach used for FSI.

Fig. 3. Discretized sections of the computational domains: blood vessel (left)

and fluid (right).

To account for the presence of the collagen fibers, anisotropic

constitutional laws should be used. However, at the present stage
of our research, such complexity could not be justified. Thus, a
simple linear elastic model has been employed using averaged
values of Young modulus equal to 1.7 MPa, Poisson ratio 0.45 and Fig. 5. Pulsatile inlet velocity profile with two characteristic points marked —
density 1100 kg/m3 [16]. systole: 0.14 s and diastole: 0.62 s.
The blood is a non-Newtonian fluid, which means the depen-
dence between shear stress and the shear rate is nonlinear [4].
This comes mainly form the presences of the blood cells, primar- The outlets of the blood domain, as well as the inlet, are
ily hematocrit [30]. To account 150 for this behavior, the Carreau considered as the artificial boundary [10], since they do not
viscosity model [30] has been applied in this paper. The model represent any physical boundary, but only the region of the
reproduces the experimental results in the large range of the interest in our investigation. In this study, the pressure boundary
shear rate [30] and takes the form condition was set for each outlet with the constant value of pres-
] m−1 sure equal to zero pascal. However, to represent the condition
η = η∞ + (η0 − η∞ ) · 1 + (λγ̇ )2 2
[
(1) at the artificial boundaries more accurately, the Kirchhoff-based
where η stands for the local non-Newtonian viscosity, γ defines model as Windkessel should be used [10]. These kind of models
the local shear rate, η0 and η∞ shear-rate viscosities are equal to make possible to impose a real boundary condition in the model
0.056 Pa s and 0.00345 Pa s, respectively. λ = 3.313 s stands for of the limited cardiovascular system section. Furthermore, they
characteristic time and m is the power-law index equal to 0.3568. make it possible to avoid the spurious effects such as artificial
The density of blood, which is a mixture of plasma, red blood pressure wave reflections, which presence could be considered
cells, white blood cells and platelets, was averaged to 1051 kg/m3 as the limitation of the simpler 0 Pa condition.
[30,31]. The fluid solver (ANSYS Fluent) calculates the solution For the structural solver, the fixed support condition was set
based on the solving set of continuity and Navier–Stokes partial for every vessel ending cross-section. In this way, these cross-
differential equations in the unstructured mesh. sections are prevented from moving and deforming.
The velocity-inlet boundary condition was implemented into Choosing an adequate time step size is essential to obtain
the solution procedure using User Defined Function (UDF). To convergence of the both solvers and the System Coupling (RMS
accurately represent flow characteristic, the velocity profile was parameter). The step size should be decreased if the forces or
represented by five polynomial functions. The pulsation pro- displacements change significantly from one time step to the
file was created from the 20 measuring points, which were next; however, a value that is too small may deteriorate the
acquired using a cardiac-gated, 2D, respiratory compensated, stability. Three values of time step sizes were set for different
phase-contrast (PC) cine sequence with through-plane velocity time intervals, which are 0.01 s, 0.0075 s, and 0.005 s.
encoding [8]. Fig. 5 illustrates the implemented velocity profile FSI cases, including the flow of the incompressible liquid
with marked measurement points. through the thin structures with low Young modulus, are al-
The cardiac output was provided with the measurement data ways unstable. The FEM solver in the first coupling iteration
and was calculated via time integral using ascending aortic flow deforms the fluid–solid interface due to forces resulting from the
data [8,32]. This value was used for the validation of the numer- fluid field. In the next, the pressure calculated by FVM solver
ical model. may strongly increase or decrease because of the fluid volume
 M. Nowak, B. Melka, M. Rojczyk et al. / European Journal of Mechanics / B Fluids 77 (2019) 273–280 277

Fig. 6. The local values calculated from the definition of the Reynolds number,
located on the inlet (left) and CoA (right) longitudinal sections. The upper
contours show the values during systole (0.14 s) and the lower contours —
show the values during diastole (0.62 s). The location of the sections is shown
in the right-hand picture. Fig. 7. Wall pressure comparison between the systole (left) and diastole (right).

change. In the subsequent coupling iteration, the interface in the of the computational domain. In the dominant portion of the
FEM solver will deform in the opposite direction. The described domain and times the Reynolds number falls below the 2300
instability mechanism, including the alternating changes of the limit, therefore the flow was treated as laminar at this stage of
pressures, forces and displacements during the solution proce- the research.
dure within every single time step, was controlled by using the The no-slip shear boundary condition was set at the fluid and
volume-based solution stabilization [33] and appropriate time solid interface. This condition causes the fluid to stick to the wall
step sizes. and move with the same velocity as the wall [35].
Blood flow in the arteries is characterized by low Reynold
numbers and minor turbulence effect, which is only observed 3. Results and discussion
under rare circumstances [5]. The factor that could result in
transition from laminar to turbulent flow is the increase of the The pressure field on the external walls of the blood domain is
flow velocity, which could appear during physical exercise or presented in Fig. 7 for the two characteristic points of the cardiac
because of the presence of a prosthetic implant or aortic steno- cycle: systole (0.14 s) and diastole (0.62 s) representing the maxi-
sis [10]. The laminar flow assumption could limit the number of mum and minimum velocities. The calculated maximum pressure
solved equations and its validity was checked by calculating the exerted on the blood walls of 4088 Pa occurred on the ascending
estimated Reynolds number ρ ude /µ, where ρ stands for the fluid aorta — the inlet. During diastole, the value of lowest pressure
density, u is the averaged velocity related to the whole object, de equal to 108 Pa was observed at the inlet. The minimum value
is the equivalent diameter and µ is the local dynamic viscosity. of static pressure, during systole, was observed in the region of
The characteristic dimension depends on each case. For the flow CoA. The Virtual Therapy simulation shown in [19] demonstrated
through the pipe, de is assumed to be the pipe diameter. The that removing the narrowing section in the descending aorta
Reynolds number was evaluated for two zones of the aorta: the positively influences the pressure field in the modeled section
ascending aorta, due to the maximum average diameter, and the and helps to stabilize the hyperpressure problem. This can be
CoA, due to the maximum velocities. The characteristic dimension attributed to the reduction of the local pressure drop in the CoA.
was computed by averaging the local diameter (effective lumen The in vivo measurements of the pressure drop across the
diameter) and was equal 0.017 m for the ascending and 0.0064 m stenosis can cause many side effects such as palpitations, chest
for the CoA zone. The local value of the viscosity was calculated pain, shortness of breath, headache, nausea or fatigue [8]. There-
by Fluent from the Carreau model. Fig. 6 depicts the contours of fore, this value can be non-invasively calculated via CFD or FSI
the local Reynolds number calculated at two characteristic points modeling as the difference of the area-weighted pressure average
of the cardiac cycle and on the longitudinal zone sections passing between two transverse sections — before and after CoA. We
through the center of the zones. created the simpler CFD model of the aorta to compare the cal-
Given that the Reynolds number relates to the whole object, culated pressure drop values between the rigid and flexible wall
the representative values for each zone were evaluated as area- assumption. For CFD the pressure drop reached 715 Pa for systole
weighted average within the sections presented in Fig. 6. This and 2 Pa for diastole. In the case of FSI model, the values equals
approach will give some indicative values that will show whether 1320 Pa (systole) and 17 Pa (diastole). The higher values for
the flow might be turbulent or not. the flexible wall model resulted from the local wall compression
During systole, the Reynolds number was equal 3172 and during aforementioned moments of the cycle.
3396, and during diastole 29 and 89, both for the inlet and Fig. 8 illustrates the velocity field and the wall of the blood
CoA zone, respectively. The flow through the pipe is assumed as domain. The highest values were observed in the CoA section and
laminar, when Re < 2300 [34]. It can be seen that the highest within the right common carotid artery. The maximum values
values exceed this limit and the flow is then defined as transi- reached about 2.90 m/s in the CoA region. The supplementary
tional (partially turbulent) [34]. However, these high values occur video of the velocity field shows minor backflow in the region
only during small part of the cardiac cycle and limited regions of the descending aorta and supraaortic vessels. However, we
278 M. Nowak, B. Melka, M. Rojczyk et al. / European Journal of Mechanics / B Fluids 77 (2019) 273–280

Fig. 10. Wall shear stress contour in the location of the maximum values for
the FSI model.

Fig. 8. Velocity field in m/s during the systole (left) and diastole (right).

Fig. 11. Fluid cross-sections through the place of the maximum deformation,
shown on the right-hand side.

observed here. The mesh cross-section for the whole cardiac cycle
Fig. 9. Comparison of the wall shear stress during the systole and diastole.
is presented on the supplementary videos. The highest values
of the pressure were situated in the location of the maximum
instantaneous deformation of the vessel wall.
analyzed the changes of the blood volume during the cycle and The total mass of the blood entering and leaving the compu-
we observed, that the backflows are during the moments when tational domain during the cardiac cycle (0 s−0.7025 s) can be
the volume is growing, so it is related to the blood accumulation calculated using the following simple expression
resulted from the flexible vessel approach. The additional factor ∫ 0.7025
influencing the backflow might be the applied simplified outlet m= ṁ (t ) dt (2)
pressure boundary condition (0 Pa). 0
The wall shear stress is one of the most important variables where ṁ (t ) represents the inlet mass stream for the mass flowing
that affect the deformation of the blood vessel wall. Fig. 9 illus- into the domain, and the sum of the every five outlet mass
trates the calculated wall shear stress at two characteristic points. streams for the mass flowing out of the domain, in the moment t.
For both CFD and FSI model, the maximum values can be seen at The integrals were calculated using the numerical method based
the upper-branch vessels, especially in the branching region (see on the parabolic interpolation between the three subsequent
Fig. 10). points of the function ṁ (t ). The mass inflowing and outflowing
Fig. 11 presents four pressure fields with shown cross-section the domain during cardiac cycle was equal in both cases, approx-
of the numerical mesh in the region of the descending aorta. The imately 0.040 kg/cycle, and the difference between the results
moments of the cardiac cycle were chosen to show the range of was equal 0.028%. The calculated cardiac output, which is the
cross-section location in space. Therefore, the aortic movement mass inflowing and outflowing the domain per minute, was about
and cross-section area changes can be analyzed from Fig. 11. The 3.264 l/min [8], which differs by 0.60% from the measurement
dynamic mesh behavior realized by smoothing and remeshing is data.
 M. Nowak, B. Melka, M. Rojczyk et al. / European Journal of Mechanics / B Fluids 77 (2019) 273–280 279

Fig. 12. Mass flow at the inlet to the artery and sum up flow through the
outlets.

Fig. 13. Distribution of the blood vessel deformation during the moment of the
The changes of the blood flow rate at the inlet and outlets maximum deformation occurrence (right), compared with the vessel wall shape
during one heart cycle is presented in Fig. 12. It can be seen that before the solution procedure (left). For better visualization, the displacement
mass stream ṁin does not directly coincide with the sum of the magnitude is set as five.

outlet streams ṁout at each time instants. This is mainly caused

by the flexibility of the vessel walls, which cause an accumulation
of the blood. The blood volume within the domain was changing of presented approach for modeling of blood flow accounting
between 24.9 cm3 and 26.2 cm3 , respectively. The negative values for the deformation of the vessel wall. The predicted difference
of blood flow rate ṁout resulted from the reversed flow at the between maximum and minimum volume of the blood domain
outlets. approached 4.7%, where the maximal displacement of the vessel
The contour of the total deformation at the maximum distor- wall was 2.35 mm. These two values support the conclusion
tion that appears at the 0.225 s of the cardiac cycle is presented that FSI should be applied to accurately simulate the blood flow
in Fig. 13. The values are calculated as within arteries whose walls are characterized by low stiffness.
The knowledge that we have gained in this study will next be
√
U = Ux2 + Uy2 + Uz2 (3)
used for modeling the blood flow where the simple single-phase
where Ux , Uy and Uz are the component deformations in di- approach will be replaced by a multifluid model, where plasma,
rections x, y, z, respectively, calculated relative to the central red and white blood cells, and the platelets will be modeled
coordinate system. The maximum value of deformation reached as separated inter-penetrating continuum phases and lumped
about 2.35 mm and was located on the descending aorta. This parameter model will account for the influence of the remaining
value is 4.7 times greater than the local vessel wall thickness. portion of the cardiovascular system.
In [15], the maximum deformation was approximately two times Interested readers are also referred to the manuscript repos-
greater. Nevertheless, in the cited work, the geometry concerns itory, which stores an animation of the vessel wall deformation,
a much older patient (30 years old compared with 8 years in the wall pressure, velocity vectors and dynamic mesh operation.
current study) whose aortas’ length and diameters are greater. For
the vessel stenosis, the maximum deformation achieved above 3
Acknowledgments
times higher value than the CoA vessel thickness.

4. Summary This research is supported by National Science Centre (Poland)

within projects No 2014/13/B/ST8/04225 and No 2014/15/D/ST8
This paper is the continuation of our previous work [19,36] /02620. This help is gratefully acknowledged herewith.
and presents an application of a hybrid approach for modeling
blood flow within elastic human thoracic arteries with lesion in
Appendix A. Supplementary data
form of coarctation at its proximal part of descending section. The
displacement of the artery wall and the influence of the moving
wall on the blood flow were modeled using the fluid–structure Supplementary material related to this article can be found
interaction approach using the commercial ANSYS package. Wall online at https://doi.org/10.1016/j.euromechflu.2019.03.009.
deformation was modeled using ANSYS Mechanical APDL while
the pressure and shear forces that act on the vessel walls due References
to blood flow were evaluated using CFD ANSYS Fluent solver.
The iterative data exchanged between the mechanical and fluid [1] WHO | Cardiovascular diseases (CVDs), WHO. URL http://www.who.int/
dynamics portion of the model was accomplished by invoking mediacentre/factsheets/fs317/en/.
implicit two-way coupling between wall movement and blood [2] M. Cierniak-Piotrowska, G. Marciniak, J. Stańczyk, Article in polish:
flow. The blood flow was modeled using a single-phase approach, Statystyka zgonów i umieralności z powodu chorób układu krążenia,
with non-Newtonian fluid properties. To mimic the flow rate gen- Zachorowalność i umieralność na choroby układu krążenia a sytuacja
demograficzna Polski (2015) 46–80.
erated by cardiac cycles, the inlet velocity profile was determined
[3] E. Wilkins, L. Wilson, K. Wickramasinghe, P. Bhatnagar, J. Leal, R. Luengo-
based on the measured mass flow rate. Fernandez, R. Burns, M. Rayner, N. Townsend, European Cardiovascular
The pressure, velocity WSS distribution, and also deforma- Disease Statistics 2017 edition, European Heart Network. URL http://www.
tion of the vessel wall were used to investigate the application ehnheart.org/images/CVD-statistics-report-August-2017.pdf.
280 M. Nowak, B. Melka, M. Rojczyk et al. / European Journal of Mechanics / B Fluids 77 (2019) 273–280

[4] P.D. Morris, A. Narracott, H. von Tengg-Kobligk, D.A. Silva Soto, S. [19] B. Melka, W. Adamczyk, M. Rojczyk, A.J. Nowak, A. Golda, Z. Ostrowski,
Hsiao, A. Lungu, P. Evans, N.W. Bressloff, P.V. Lawford, D.R. Hose, Virtual Therapy Simulation for Patient with Coarctation of Aorta Using
J.P. Gunn, Computational fluid dynamics modelling in cardiovascular CFD Blood Flow Modelling, Springer International Publishing, Cham, 2017,
medicine, Heart (2015) http://dx.doi.org/10.1136/heartjnl-2015-308044, pp. 153–160, http://dx.doi.org/10.1007/978-3-319-47154-9_18.
heartjnl–2015–308044– URL http://heart.bmj.com/content/early/2015/10/ [20] N. Westerhof, F. Bosman, C.J. De Vries, A. Noordergraaf, Analog studies of
28/heartjnl-2015-308044.full. the human systemic arterial tree, J. Biomech. 2 (1969) 121–143.
[5] Y. Bazilevs, K. Takizawa, T.E. Tezduyar, Computational Fluid- [21] J.F. LaDisa, C.A. Taylor, J.A. Feinstein, Aortic coarctation: Recent develop-
Structure Interaction, John Wiley & Sons, Ltd, Chichester, UK, 2013, ments in experimental and computational methods to assess treatments
http://dx.doi.org/10.1002/9781118483565, URL http://doi.wiley.com/ for this simple condition, Progr. Pediatr. Cardiol. 30 (1–2) (2010) 45–49,
101002/9781118483565. http://dx.doi.org/10.1016/j.ppedcard.2010.09.006.
[6] T. Mansi, K. McLeod, M. Pop, K. Rhode, M. Sermesant, A. Young (Eds.), [22] J. Lantz, T. Ebbers, J. Engvall, M. Karlsson, Numerical and experimental
Statistical atlases and computational models of the heart, in: Imaging assessment of turbulent kinetic energy in an aortic coarctation, J. Biomech.
and Modelling Challenges, in: Lecture Notes in Computer Science, vol. 46 (11) (2013) 1851–1858, http://dx.doi.org/10.1016/j.jbiomech.2013.04.
10124, Springer International Publishing, Cham, 2017, http://dx.doi.org/ 028.
10.1007/978-3-319-52718-5, URL http://link.springer.com/101007/978-3- [23] A. Menon, D.C. Wendell, H. Wang, T.J. Eddinger, J.M. Toth, R.J. Dho-
319-52718-5. lakia, P.M. Larsen, E.S. Jensen, J.F. LaDisa, A coupled experimental and
[7] E. Rosenthal, Coarctation of the aorta from fetus to adult: curable con- computational approach to quantify deleterious hemodynamics, vascular
dition or life long disease process?, Heart 91 (11) (2005) 1495–1502, alterations, and mechanisms of long-term morbidity in response to aortic
http://dx.doi.org/10.1136/hrt.2004.057182, URL http://heart.bmj.com/cgi/ coarctation, J. Pharmacol. Toxicol. Methods 65 (1) (2012) 18–28, http:
doi/101136/hrt.2004057182. //dx.doi.org/10.1016/j.vascn.2011.10.003, arXiv:NIHMS150003.
[8] Source of the medical data from, http://www.vascularmodel.com/, [24] U. Gülan, B. Lüthi, M. Holzner, A. Liberzon, A. Tsinober, W. Kinzelbach,
(accessed: 01.03.16). An in vitro investigation of the influence of stenosis severity on the
[9] P.M.S. Syamasundar Rao Patnana, Coarctation of the Aorta: Background, flow in the ascending aorta, Med. Eng. Phys. 36 (9) (2014) 1147–1155,
Pathophysiology, Prognosis. URL https://emedicine.medscape.com/article/ http://dx.doi.org/10.1016/j.medengphy.2014.06.018.
895502-overview. [25] M. Toloui, B. Firoozabadi, M. Saidi, A numerical study of the effects
[10] A.V.L. Formaggia, A. Quarteroni, Cardiovascular Mathematics, 2009. URL of blood rheology and vessel deformability on the hemodynamics of
http://link.springer.com/10.1007/978-88-470-1152-6. carotid bifurcation, Sci. Iran. 19 (1) (2012) 119–126, http://dx.doi.org/10.
[11] M. O. J. M. M. Malve, A. Garcia, Unsteady blood flow and mass transfer 1016/j.scient.2011.12.008, URL http://linkinghub.elsevier.com/retrieve/pii/
of a human left coronary artery bifurcation: FSI vs. CFD, Int. Commun. S102630981100263X.
Heat Mass Transfer 39 (6) (2012) 745–751, http://dx.doi.org/10.1016/j. [26] A. Amindari, L. Saltik, K. Kirkkopru, M. Yacoub, H.C. Yalcin, Assessment
icheatmasstransfer.2012.04.009. of calcified aortic valve leaflet deformations and blood flow dynamics
[12] K.M. Khanafer, J.L. Bull, R. Berguer, Fluid–structure interaction of turbulent using fluid structure interaction modeling, Inform. Med. Unlocked 9 (July)
pulsatile flow within a flexible wall axisymmetric aortic aneurysm model, (2017) 191–199, http://dx.doi.org/10.1016/j.imu.2017.09.001, URL http://
Eur. J. Mech. B/Fluids 28 (1) (2009) 88–102, http://dx.doi.org/10.1016/j. linkinghub.elsevier.com/retrieve/pii/S2352914817300369.
euromechflu.2007.12.003. [27] P. Reymond, P. Crosetto, S. Deparis, A. Quarteroni, N. Stergiopulos,
[13] R. Jayendiran, B. Nour, A. Ruimi, Dacron graft as replacement to dissected Physiological simulation of blood flow in the aorta: Comparison of
aorta: A three-dimensional fluid–structure-interaction analysis, J. Mech. hemodynamic indices as predicted by 3-D FSI, 3-D rigid wall and
Behav. Biomed. Mater. 78 (2017) (2018) 329–341, http://dx.doi.org/10. 1-D models, Med. Eng. Phys. 35 (6) (2013) 784–791, http://dx.doi.
1016/j.jmbbm.2017.11.029, URL http://linkinghub.elsevier.com/retrieve/pii/ org/10.1016/j.medengphy.2012.08.009, URL http://www.sciencedirect.com/
S1751616117305179. science/article/pii/S1350453312002226.
[14] M. Alishahi, M.M. Alishahi, H. Emdad, Numerical simulation of blood [28] 3D Systems, Geomagic. URL http://www.geomagic.com/en/products/
flow in a flexible stenosed abdominal real aorta, Sci. Iran. 18 (6) (2011) design/overview.
1297–1305, http://dx.doi.org/10.1016/j.scient.2011.11.021. [29] H. Gray, H.V. Carter, T.P. Pick, R. Howden, Anatomy : descriptive and
[15] P. Crosetto, P. Reymond, S. Deparis, D. Kontaxakis, N. Stergiopulos, A. Quar- surgical, Bounty (2012).
teroni, Fluid–structure interaction simulation of aortic blood flow, Comput. [30] Y.I. Cho, K.R. Kensey, Effects of the non-Newtonian viscosity of blood on
& Fluids 43 (1) (2011) 46–57, http://dx.doi.org/10.1016/j.compfluid.2010. flows in a diseased arterial vessel. Part 1: Steady flows, Biorheology 28
11.032. (September) (1991) 241–262.
[16] M. Neidlin, S.J. Sonntag, T. Schmitz-Rode, U. Steinseifer, T.A.S. Kaufmann, [31] T. Kenner, The measurement of blood density and its meaning, Basic Res.
Investigation of hemodynamics during cardiopulmonary bypass: A mul- Cardiol. 84 (2) (1989) 111–124, http://dx.doi.org/10.1007/BF01907921.
tiscale multiphysics fluidstructure-interaction study, Med. Eng. Phys. 10 [32] M. Lavdaniti, Invasive and non-invasive methods for cardiac ouput mea-
(0) (2016) http://dx.doi.org/10.1016/j.medengphy.2016.01.003, URL http: surement, Int. J. Caring Sci. 1 (3) (2008) 112–117, http://dx.doi.org/10.
//www.elsevier.com/locate/medengphy. 1016/j.bpg.2011.02.003.
[17] L. Taelman, J. Degroote, A. Swillens, J. Vierendeels, P. Segers, Fluid– [33] ANSYS⃝ R
Academic Research, Release 17.2, Help System, ANSYS, Inc..
structure interaction simulation of pulse propagation in arteries: Numerical [34] F. White, Fluid Mechanics, McGraw-Hill, New York, 2010.
pitfalls and hemodynamic impact of a local stiffening, Internat. J. Engrg. [35] H. Versteeg, W. Malalasekera, Introduction to Computational Fluid
Sci. 77 (2014) 1–13, http://dx.doi.org/10.1016/j.ijengsci.2013.12.002. Dynamics, Vol. 44, 2005, URL http://www.ncbi.nlm.nih.gov/pubmed/
[18] D. De Wilde, B. Trachet, N. Debusschere, F. Iannaccone, A. Swillens, J. 6686412.
Degroote, J. Vierendeels, G.R. De Meyer, P. Segers, Assessment of shear [36] B. Melka, M. Gracka, W. Adamczyk, M. Rojczyk, A. Golda, A.J. Nowak,
stress related parameters in the carotid bifurcation using mouse-specific R.A. Białecki, Z. Ostrowski, Multiphase simulation of blood flow within
FSI simulations, J. Biomech. 49 (11) (2016) 2135–2142, http://dx.doi.org/ main thoracic arteries of 8-year-old child with coarctation of the aorta,
10.1016/j.jbiomech.2015.11.048. Heat Mass Transf. (2017) URL http://link.springer.com/101007/s00231-
017-2136-y.