Modeliranje Toka U Krvnim Žilama
Modeliranje Toka U Krvnim Žilama
Modeliranje Toka U Krvnim Žilama
Presentee par
Xiaofei WANG
DOCTEUR DE
i
Abstracts
The vascular network consists of a very large number of segments with various properties
and thus the pulsatile blood flow inside is very complicated. With the time-domain-
based nonlinear 1D model, this thesis studies the blood flow in networks, focusing on
the numerical computing and several applications.
With assumptions of long wave and axisymmetric velocity profile, the 1D governing
equations of mass and momentum are derived by integrating the continuity and Navier-
Stokes equations along the radius. A Kelvin-Voigt viscoelastic model is adopted for the
constitutive equation of the tube. This leads to a nonlinear hyperbolic-parabolic sys-
tem, which is then solved with four numerical schemes, namely: MacCormack, Taylor-
Galerkin, Monotonic Upwind Scheme for Conservation Law (MUSCL) and local dis-
continuous Galerkin. The schemes are implemented in MATLAB and the numerical
solutions are checked favorably against analytical, semi-analytical solutions and clinical
observations. Among the numerical schemes, comparisons are made in four important
aspects: accuracy, ability to capture shock-like phenomena, computational speed and
implementation complexity. The suitable conditions for the application of each scheme
are discussed. After this, a general purpose C++ code is developed and tested on several
networks: a circle of arteries, a human systemic network with 55 arteries and a mouse
kidney with more than one thousand segments. The time dependent distribution of
pressure in the networks is visualized and the propagation patterns of the waves are well
captured. Good speedup is achieved by parallelizations of the code.
The developed code is applied in three studies. First, the coefficients of fluid fric-
tion and wall viscosity are determined with aides of a well defined experimental setup.
Because both the two factors damp the pulse waves, they are difficult to evaluate sepa-
rately. We estimate them in pairs by fitting the 1D viscoelastic model against pressure
waves measured on the experimental setup. The fitted values of viscoelastic parameters
are consistent with values estimated with other methods. The effect of wall viscosity on
the pulse wave has been shown in the same order of that of fluid viscosity. Second, with
time series of pressure and diameter measured in several locations of the sheep arte-
rial network, the viscoelasticity parameters are estimated. With those values, the pulse
waves in the sheep network are simulated and the effect of viscoelasticity is investigated.
Numerical solutions show that the viscoelasticity damps significantly the high frequency
components of the pulse waves. Third, we simulate the change of blood flow induced
by the axillofemoral and femoral-femoral anastomoses with a severe iliac stenosis. The
influence of the bypassing path is studied.
ii
Acknowledgements
The working environment influences a scientist like the growing environment to a fruit
tree. Without the help of many people, I cannot yield with three years of research the
fruitthis manuscript of thesis.
First of all, I would like to thank my supervisor Pierre-Yves LAGREE. Besides being
an active researcher in multiple directions of fluid mechanics, he is also an excellent
advisor, who really cares the full development of PhD students. He believes that curiosity
is a very important driving factor of scientific discoveries and thus respects my research
interests. My co-supervisor, Jose-Maria FULLANA, also gives me advices in a lot of
aspects, from doing research to making presentations. I benefit from both my supervisors
not only in scientific knowledge, but also in the good working conditions they create for
me. They cover generously my national and international travel expenses for conferences
and collaborations. I feel very honored and lucky to do my PhD thesis under their
supervision.
The collaborators of our team of course deserve my gratitudes because this thesis
directly benefits from their works or ideas. I give a big thank you to Prof. Mami MAS-
TUKAWA and Mr. Shohei NISHI (Doshisha Univerisity, Japan) for the collaborations on
the experimental measurements; to Dr. Jean-Frederic GERBEAU (Director of research,
INRIA Paris-Rocquencourt, France) for his guiding me to this research topic; to Prof.
Ricardo L. ARMENTANO (Favaloro University, Argentina) for his groups measurement
data on sheep; to M.D. Salam ABOU TAAM (Hospital Foch, France) for his medical
knowledge in anastomosis surgeries; to Dr. Sylvie LORTHOIS (Director of research,
CNRS, Institut de Mecanique des Fluides de Toulouse, France) and Dr. Cecile DU-
PLAA (Director of research, INSERM, U1034 Adaptation cardiovasculaire a lischemie)
for the geometrical data of a mouse kidney vascular; to Dr. Olivier DELESTRE (Uni-
versite de Nice Sophia-Antipolis, France) for the discussions on the numerical methods;
and to Dr. Helene COULLON (University of Orleans, France) for her useful advices on
code development.
It has been a great pleasure to work in the same office with Julien PHILIPPI, Gounseti
PARE, Luca MARGHERI and Cansu OZHAN. During my PhD research, I have dis-
cussed a lot with colleagues in the Institute of dAlembert and ICS (Institut du Calcul
et de la Simulation, UPMC). I am very sorry that I cannot list their names here one by
one because the list would be very long. I really enjoy the discussions with them and
my research life is much more colorful because of them.
I want to acknowledge the scholarship of China Scholarship Council (2010-1013) and
the French state funds managed by CALSIMLAB and the ANR within the investisse-
ments dAvenir programme under reference ANR-11-IDEX-0004-02 (2013-2014).
I give special thanks to the two referees of my thesis: Prof. Patrick SEGERS (Ghent
iii
University, Belgium) and Prof. Franck NICOUD (Universite Montpellier 2, France). I
also acknowledge the great gentleness of the other four examiners of the jury committee
besides my two supervisors: Prof. Stephane ZALESKI (Institut dAlembert, UPMC,
France), Dr. Irene VIGNON-CLEMENTEL (INRIA Paris-Rocquencourt, France), Dr.
Sylvie LORTHOIS (CNRS, Institut de Mecanique des Fluides de Toulouse, France), and
Prof. Eleuterio F. TORO (University of Trento, Italy).
Last but not least, my deepest gratitude shall go to my family. Without their support,
I could not go as far as today in pursuing a dream in academia.
iv
Contents
1. Introduction 1
1.1. Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2. Blood flow in vascular networks from a mechanical view . . . . . . . . . . 3
1.3. Models for blood flow in distensible tubes . . . . . . . . . . . . . . . . . . 6
1.4. Problems for one-dimensional models . . . . . . . . . . . . . . . . . . . . . 8
1.5. Research context of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.6. Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2. Mathematical models 16
2.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2. 1D mathematical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.1. Integration along radius . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.2. Profile of axial velocity . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.3. Viscoelastic tube law . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.3. Characteristic structure of the system . . . . . . . . . . . . . . . . . . . . 23
2.4. Initial and boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4.1. Initial conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4.2. Inlet and outlet of the homogeneous hyperbolic part . . . . . . . . 25
2.4.3. Conjunction points . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.5. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3. Numerical Schemes 29
3.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2. Numerical solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.1. Operator splitting . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.2. MacCormack scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2.3. Taylor-Galerkin scheme . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2.4. MUSCL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.2.5. Treatment of the parabolic subproblem . . . . . . . . . . . . . . . 36
3.2.6. Local Discontinuous Galerkin scheme . . . . . . . . . . . . . . . . . 37
3.3. Verification and comparison among schemes . . . . . . . . . . . . . . . . . 39
3.3.1. Propagation in a uniform tube . . . . . . . . . . . . . . . . . . . . 40
3.3.2. Attenuation due to the viscosity of blood . . . . . . . . . . . . . . 43
3.3.3. Diffusion due to the viscosity of the arterial wall . . . . . . . . . . 44
v
3.3.4. Shock-like phenomena due to the nonlinearity . . . . . . . . . . . . 45
3.3.5. Reflection and transmission at a branching point . . . . . . . . . . 46
3.3.6. Application on a full systematic arterial network . . . . . . . . . . 47
3.4. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
II. Applications 67
5. Fluid friction and wall viscosity of the 1D blood flow model: study with an
in-vitro experimental setup 68
5.1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.2. Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2.1. One dimensional model . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2.2. Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.2.3. Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . . 72
5.3. Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.3.1. Youngs modulus . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.3.2. Fluid friction and wall viscosity . . . . . . . . . . . . . . . . . . . . 75
5.4. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.5. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
vi
6.4. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.5. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
vii
1. Introduction
1.1. Background
The most easily noticed feature of blood flow in arteries is perhaps the pulsation phe-
nomenon. It has been studied and used in medical practices for thousands of years.
For example, in traditional Chinese medicine, the practitioners palpate on six locations
of the two wrists. Those positions are mapped to various organs as shown in Fig. 1.1.
The pluses are classified by several characteristic elements, such as depth, rate, regular-
ity, width, smoothness etc. The practitioners determine the health status of individual
organs and of the whole body by examining those characteristics. Several traditional
classic medical books (e.g. Nei Jing, 475221 B.C.) have summarized the empirical
knowledge accumulated by generations of practitioners. But we did not find anatomical
evidence for the mapping between the palpating positions and the organs. Moreover
those methods are described in qualitative manners, thus there are no precise standards
for the practitioners to follow. Although the traditional medicine is controversial in
many aspects, it is still an important part of the healthcare system in several asian
countries. Scrutinizing those methods with modern techniques and medical knowledge
is now an active research field [17, 82, 90]. If strong statistical correlations between
diseases and some characteristics of pulses are established, palpating pulses is obviously
a very advantageous prognosis method. In contrast with this empirical approach, which
dose not investigate the wave propagation mechanism, the deterministic approach relies
on deep understanding of the pulse wave behaviours in the cardiovascular system.
The modern concept of blood circulation was established by Harvey in 17th century.
Since then many talented mathematicians and physicists have been attracted to study
the blood flow in the vascular system. In 1775, L. Euler derived the partial differential
equations describing the conservation of mass and momentum for an inviscid fluid. But
he could not solve the equations and thus did not recognize the wave characteristics
of the system. T. Young gave the first scientific description of waves of blood flow in
1808. G. Hagen (in 1839) and Jean L. M. Poiseuille (in 1838) derived independently
the resistance formula of the laminar flow in a long tube. In 1878, A. I. Moens and
D. Korteweg established the formula of the speed of waves in a distensible vessel. B.
Riemann introduced the method of characteristics and it was applied to the study of
arterial blood flow since the decade of 1950. John R. Womersleys approach is to neglect
the nonlinear convection term, and make analysis on frequency domain [72].
The aforementioned studies were done before the advent of powerful digital comput-
ers. In recent decades, the box of research tools has been enriched a lot due to the
progress of computer science, image analysis and computational fluid dynamics (CFD).
The majority of current simulations are done on computers by solving three-dimensional,
1
Figure 1.1.: Positions to palpate pulses in traditional Chinese medicine (from Ref. [21]).
2. To develop non-invasive diagnosis tools. The most well known quantities of blood
flow used as diagnosis purpose must be the heart rate and the systolic/diastolic
pressure. Those biomarkers can be measured easily and are reliable for some dis-
eases. Recently, more tools based on analysis of pulse waves have been proposed
for early detections of vascular disorders, such as pulse wave velocity (PWV), car-
dioankle vascular index, pulse wave separation method et al. [66]. Development of
2
techniques to estimate parameters from flow rate and pressure makes the applica-
tions on this direction promising [40, 44].
3. To optimize the design of intravascular prostheses and surgery planning. We have
seen some simulations of bypass graft surgery in both coronary [68] and periph-
eral [33, 34, 48, 95] vascular systems. The insertion of a stent changes the local
compliance of the vessel. The pressure reflection and WSS alterations incurred
by the stent can be predicted by computing of an unsteady blood flow field with
moving boundaries [53].
3
Figure 1.2.: Arterial system. From Grays anatomy 36th edition by Williams & Warwick.
4
The stiffness and viscosity of the wall are bigger than those of large arteries. Par-
tially due to the weakened buffering effect (reflection at peripheral sites may be
another reason), the pulse pressure (difference between systolic and diastolic pres-
sure) increases, while the mean pressure drops slightly (see Fig. 1.3). As shown
in Table 1.1, due to many branches in this level, the total cross section increases.
The mean flow in those arteries is smaller than in the large arteries, thus the non-
linear convection effect becomes smaller as well. Because of this, the profile of the
waveform becomes less front-steepened than in the large arteries.
Arteriole. This kind of arteries are characterized by layers of muscle cells which are
regulated by autonomic nervous system and hormones. According to Poiseuilles
theory, the resistance of the a tube to the flow is CL/R4 , where L is the length
of the tube, R is the radius and C is a constant coefficient. From Table 1.1, we
see that the characteristic diameters of large arteries and arterioles are 2cm and
0.05cm respectively. Even though we assume the arterioles are 1/10 in length and
10 times in total cross-sectional area of the large arteries, the resistance of the
former is 25600 times of the latter. The main resistance of the network mainly
resides in this level of vessels and the capillary bed. As shown in Fig. 1.3, the
pressure starts to drop sharply at this level. The muscle cells can contract or relax
to change the resistance significantly. The blood supply is thus regulated according
to the metabolic needs of the organs.
Capillaries. The wall of capillaries is formed by only one layer of endothelial linings.
To maximize the total area of the capillary bed, the diameter of capillaries is so
small that the red blood cells have to deform to pass through them. This structure
facilitates material exchange. Due to the buffering effect of the arterial network,
the pressure pulsation becomes very weak here, and the blood moves steadily.
Venules. Vessels at this level collect blood from the capillary bed and forward
them further to larger veins.
Veins. The pressure at this level drops to almost zero (in lying down position).
Compared with arteries, veins have some different features. First, there are valves
in some veins to keep the flow uni-directional. Even though there is no big re-
sistance in the venous network for the blood to go back to the heart, there lacks
pressure gradient to balance the gravity in the lower limbs. Muscular contractions
around deep veins help to propel the blood upward. Second, the wall may collapse
because there are cases when the internal pressure is smaller than external pres-
sure. Third, the wall of veins is thin and the Youngs modulus is much smaller. The
veins are easy to dilate to accommodate blood. According to the wave theory, the
wave speed is also smaller than that in arteries. In same cases, the skeleton muscle
contraction causes very large acceleration of the blood flow, such that shock-like
phenomenon may form in the veins. We will discuss the possibilities in this thesis
(Chapter 3).
5
Figure 1.3.: Systemic blood pressure. From Anatomy & Physiology, OpenStax College.
If we truncate at the level of very small arteries where the pulsation is not significant,
the arterial system can be considered as a very large network of distensible tubes whose
impedance mainly resides in the outlets of the network. As shown in Table 1.1, for a
20kg-dog, the number of vessels in this truncated network can be the order of magnitude
of 105 . The main features of blood flow in this network can be summarized as: pulsation
in time, wave transmission in space and large number of vessel segments.
6
System Section Mean Number Total
level Diameter(cm) of vessels cross section
1 Aorta (19-4.5) 1 (2.8-0.2)
2 Arteries 4.000 40 5.0
3 Arteries 1.300 500 6.6
4 Arteries 0.450 6103 9.5
5 Arteries 0.150 1.1105 19.4
6 Arterioles 0.050 2.8106 55.0
7 Capillaries 0.008 2.7109 1357.0
8 Venules 0.100 1.0107 785.4
9 Veins 0.280 6.6105 406.4
10 Veins 0.700 4103 154.0
11 Veins 1.800 2100 53.4
12 Veins 4.500 110 17.5
13 Venae Cavae (5-14) 2 (0.2-1.5)
Table 1.1.: Vascular dimensions in a 20-kg dog (reproduced from Milnor 1989)
quite realistic boundary conditions with a limited computational cost [8, 88].
7
works, the problem needs to be computed millions of times [39]. With the 1D
model, it is possible in the near future to have fast, accurate, patient-specific and
clinically applicable solutions for the simulation of full bodies.
Lateral leakage and tapering. The aorta is tapering from the proximal to the distal
end. The tapering wall will increase the pulse pressure by continuous reflection.
It is also observed that there are tiny side branches at the aorta. Some studies
show that the effect of tapering wall is compensated by the blood loss to the side
branches. In net effect, it is more like a flow in a straight tube with no lateral
leakage [52]. However, both the branches and tapering can be modeled by the
1D models. The net effect needs further research. In our model, we just assume
that there is no distributive blood loss, because the arterial wall is considered as
impermeable.
8
Wave length is much longer than vessel radius. In the aorta, the wave speed is
about 5 m/s and the radius is about 10 mm [52]. If we assume the period of
one pulse is 1 second, the ratio between the wave length and the radius is 500.
In smaller arteries, the wave speed increases and the radius decreases, the ratio
becomes even bigger. Thus the long wave assumption is fully justified.
Profile of axial velocity. One crucial step in deriving the 1D model is to prescribe a
proper axial velocity profile. The fluid friction coefficient and the correction factor
of the convection term in the momentum balance equation are dependent on this
profile. The ratio between the transient inertia
p force and the viscous force can
be estimated by the Womersley number R /, where R is the internal radius,
the angular frequency and the kinematic viscosity of blood. Table 1.3 shows
the Womersley number at various arteries calculated using a cardiac frequency of
2 Hz, a kinematic viscosity of blood 3.3 106 m2 /s. From the table, we see that
the inertia force are bigger than viscous force in larger vessels. In this case, the
profile is essentially flat in the central part and there is a thin viscous boundary
layer to match the no-slip boundary condition near the wall. When the Womersley
number is very small in the small vessels, the profile is more close to a parabolic or
Poiseuille flow. We will discuss this in more details in chapter 2 and 5. We usually
assume that the profile is the same within one segment, while it may vary through
the network.
Turbulence. Since the mean Reynolds number (computed from the mean flux
between systolic and diastolic ones) in aorta is over 4000 (see Table 1.3 and [79]),
the flow should be described as turbulent if it is in a straight pipe with the given
steady flux. In unsteady case, we consider, the flow is at least transitional to
turbulence if it is not turbulent. In small arteries, the Reynolds number drops a
lot because both the velocity and the diameter of the vessels become smaller. Given
that the Womersley number is also much smaller in small arteries, laminar flow
is very likely to appear in small arteries. That gives justification to Poiseuilles
theory in estimating the peripheral resistance. In medium-sized arteries, there
may be complex transition conditions between turbulent and laminar flow. The
1D models neglect the velocity components in the circumferential direction and
thus a laminar flow is an implicit assumption. However, the dissipating effect of
turbulence on energy can be lumped into the skin friction term.
Constitutive equation of vessel wall. The constitutive equation of the wall is derived
from thin shell theory. The characteristic thickness of wall is 0.2 mm and the radius
of vessel is 2 mm. Then the ratio between the two is 0.1. The wall shows nonlinear
elasticity properties in high pressure and the inertia of the wall may affect the pulse
wave as well. But in most cases, a linear one is enough for a first approximation.
Some recent studies show that viscoelastic property has a considerable effect on
the pulse waves. We will investigate further on the viscoelasticity in this study.
Gravity of blood. Given that the density of mercury is 13.56 times of that of water,
9
Diameter Total Total Fluid Wave Womersley
resistance cross area speed speed number
Table 1.2.: Changes of characteristics of the arterial tree from aorta to arterioles
the hydrostatic pressure at the foot is about 125 mmHg (1 mmHg 133 Pa) for
a normal adult man with a hight of 170 cm in standing posture. This is indeed
a big number compared with the pressure in arteries (80-120 mmHg). In clinical
applications, the blood pressure is measured when the subject is in supine or siting
posture. In those cases, the pressure difference through the whole body caused by
gravity is much smaller than the pressure at the arterial system. In this study, we
will not consider the gravity. But we note that it can be added as a source term if
necessary.
Curvature and bifurcation. Both the curvature and bifurcation will dissipate energy
of the blood flow. Experimental study shows that the energy loss at bifurcations
is very small [50]. Our experiments on medium-sized tubes show that the energy
dissipation caused by curvature is also negligible.
Boundary conditions. Because the network of vessels is very large, we usually have
to truncate it at some levels. The truncated subnetworks have to be described
properly. Reflection coefficient, structured tree [56], 0D models [23] and some
other kind of generalized methods for specific numerical schemes [87] have been
investigated as the boundary conditions of the 1D models.
Extension to veins. With some modifications, the 1D models can also be used in
simulation of venous flow. There are three special features of venous system that
need to be treated: inflow from venules bed, muscular contraction and valves. The
inflow can be considered as a source term in the mass conservation equation; the
muscular contraction can be described by an external pressure in the constitutive
equation of the vessel wall; the valves, which in fact allow very small reverse flows,
can be modeled by very large resistances when the pressure gradient is unfavorable
to the blood flow to the heart [27].
The 1D model may be used to simulate the blood flow in the arterial tree truncated
at the level of arteries as small as 0.1mm in diameter. From aorta to small branches of
the tree, the changes of characteristics are summarized in Table 1.2.
10
Section Radius (cm) Womersley Reynolds
Ascending Aorta 0.75 14.628 4543
Descending Aorta 0.65 12.677 -
Abdominal Aorta 0.45 8.777 -
Femoral Artery 0.2 3.901 605
Carotid Artery 0.25 4.876 1060
Arteriole 0.0025 0.049 -
Capillary 0.0003 0.006 -
Venule 0.002 0.039 -
Inferior Vena Cava 0.5 9.752 -
Main Pulmonary artery 0.85 16.578 -
11
Beside the poor accuracy of non-invasive measurement techniques on the in vivo
data, lack of patient-specific information on model parameters is another reason.
To know the sensitivity of the output to the uncertainties of each parameter of the
1D models, some sensitive analysis have been done [16, 39, 98]. To estimate the
parameters, we can compare the model output with measurements and minimize
the difference by tuning the parameters [40, 44]. Many conjectured applications
of the 1D models rely on the solution of the reverse problems. But researches on
this direction are only on early stage. In these studies, the computations have to
be run millions of times, which contribute to our motivations to develop a fast
numerical code (Chapter 4).
In one class of applications of 1D models, we simulate the diseased blood flow and
various surgery plans hoping that optimized guidance to surgeons may be given
accordingly. Even though there are a lot of difficulties on setting the parameters of
patient-specific models and usually only qualitatively accurate predictions can be
given, we still can extract a lot of clinically relevant information. For example, Hu-
berts et al. show that the mean pressures and flows after an arteriovenous fistula
surgery can be simulated with 1D models with quite satisfactory accuracy [34].
Another successful example is the simulation of blood flow after femoro-popliteal
bypass surgery [95]. Inspired by this, we simulated a whole systemic arterial net-
work with a femoral stenosis and the blood flow after axillo-femoral, femorofemoral
bypass surgeries. The applications of 1D models in this scenario are discussed
(Chapter 7).
12
diffusive term cannot be treated by a simple central numerical flux and a local auxiliary
variable was introduced to stabilize the discontinuous Galerkin scheme (thus the name
local discontinuous Galerkin). We implemented the four numerical methods in Matlab,
verified them and compared them in performance. The proper application conditions
for each scheme were discussed accordingly.
Chapter 4 treats the development of a parallel code. Noting that a fast code is needed
in many cases, we developed a parallel C++ code with two of the four schemes. This
code was verified and the speedup was tested on multicore computers by computing
networks with various number of arteries. Tested on a normal workstation with 12
cores, the classic systemic network with 55 arteries can be computed in a few seconds
for one cardiac cycle. Compared with a sequential one, the speedup is at least 8. We
also highlighted that an anatomically accurate network of mouse kidney with about
1500 segments (built by Dr. Sylvie LORTHOIS and Dr. Cecile DUPLAA) can be
computed in a few minutes. This code enables us to do sensitivity analysis and to fit
parameters in a reverse way of large networks. In both those cases, the code speed is
crucial because the computations have to be run millions of times. Another highlight of
the code is that the configuration of simulations can be set conveniently through several
comma separated value (CSV) files. The topology of the network is set by a directed
acyclic graph (DAG), which can describe a quite general class of networks with various
branching and convergence types. To record the evolution of interested quantities, the
location and time interval can also be set flexibly. All other configuration options such
as viscoelasticity and meshes are also set by CSV files. The output files can be read by
Python and Gnuplot for post-processing. The C++ code does not depend on any third
party libraries, and all steps of the simulation and data analysis can be done with free
open source tools.
Chapter 5 studies the effect of fluid and wall viscosity on the pulse wave with aids of an
experimental setup built in Doshisha University (Mami MASTUKAWA group, Kyoto,
Japan). In the experimental setup, a distensible polymer tube was connected with a
piston pump. We simultaneously recorded the wall displacement with a laser surface
velocimeter (LSV) and the internal pressure with a pressure sensor. The viscoelasticity
of the tube was calculated from the two time series. This value was compared with the
value fitted from the viscoelastic 1D model. The measurement of viscoelasticity is more
accurate than tensile test because the precision of both the LSV and the pressure sensor
is very high. Previous studies fit the viscoelasticity with an inviscid fluid model, or just
used measurements with tensile test without fitting; thus there is a doubt on the values
of the fluid viscosity and wall viscoelasticity, because both of them are damping factors
of the pulse wave. This research helps clear this doubt.
Chapter 6 shows a study of the effect of viscoelasticity on the pulse waves in the
systemic arterial tree of sheep. The research group led by R. L. ARMENTANO (Favaloro
University, Argentina) recorded the diameter and internal pressure at several locations
of the in vivo arterial tree of a group of sheep. We took advantage of the data and
estimated the viscoelasticity values. With those values, we built a sheep network and
computed the pulse waves with and without viscoelasticity. The damping effect of high
frequency components in the waveforms was illustrated.
13
Chapter 7 simulates the blood flow before and after operations of axillofemeral and
femeral-femeral anastomosis. The influence of stenosis on femoral artery and the choice
of grafts with different elasticity were discussed. This work was done in collaboration
with vascular surgeon S. ABOU TAAM (Hospital Foch, Paris, France).
Chapter 8 summarizes the thesis and makes some perspectives.
14
Part I.
15
2. Mathematical models
2.1. Introduction
The previous chapter gives the assumptions of the 1D model and the justifications for
vascular blood flow, especially in large and medium-sized arteries. This chapter first talks
about the derivation of a 1D viscoelastic model, highlighting the approximations made
on the axial velocity profile and the vessel mechanical law. After this, the hyperbolic-
parabolic structure of the 1D model is discussed. The boundary conditions at the inflow,
outflow and internal conjunction points are treated with the characteristic method.
Figure 2.1.: Illustration of a distensible vessel segment. The flow rate is Q, the internal
pressure P , the internal radius R and thus the cross sectional area A =
R2 . The two perpendicular surfaces at the ends of the vessel do not move
horizontally.
16
and all other dependent variables do not depend on . Further invoking the long wave
assumption, we can reduce the mass conservation equation and Navier-Stokes equations
to the following form,
1 (rvr ) vx
+ = 0, (2.1a)
r r x
vx vx vx 1 P vx
+ vr + vx + = r , (2.1b)
t r x x r r r
1 P
= 0. (2.1c)
r
Those governing equations are accompanied by the no-slip boundary condition
vx (R) = 0,
Before going forward, we recall the Reynolds transport theorem for a function F (r, x)
in the one dimensional domain 0 x R(x, t),
Z R Z R
d F (r, x) R
F (r, x)dr = dr + F (r, x)|r=R .
dx 0 0 x x
If one invokes no-slip condition, vx (R) = 0, the last term vanishes. By comparing this
equation with Eq. (2.2), one gets
Z R Z R
(rvr ) d
2 dr + 2 vx rdr = 0,
0 r dx 0
or Z R
d
2Rvr (R) + 2vx rdr = 0.
dx 0
17
We assume there is no gap between the fluid and the wall, thus vr (R) = R/t. Also
noting that the flow rate Q and cross area A are
Z R
Q= 2vx rdr, A = R2 ,
0
we obtain
Q2
Q A P vx
+ + = 2 r . (2.4)
t x A x r r=R
the profile of the axial velocity. In general, we assume that the velocity profile in one
segment has the same shape (r), or vx (r, x, t) = U (x, t)(r) with U (x, t) as the average
fluid velocity. The friction term can be rewritten as Cf U or Cf Q/A, where Cf is
dependent on (r).
18
Figure 2.2.: Velocity profile in a stiff tube with various Womersley numbers. Reproduced
from [65].
The Womersley number, defined as = R / (it is usually denoted as , but we
have used for the momentum correction coefficient), gives a strong indication on the
profile. By applying Womersleys analysis method on a rigid tube, we can derive the
two coefficients in two asymptotic cases [66],
= 4/3 Cf = 8, << 1
= 1 Cf = 2, >> 1.
However, under physiological conditions the Womersley number is about 10 in large
arteries and about 0.1 in small arteries. The values of the two coefficients should fall
between the two cases listed above (see Fig. 2.2).
Hughes and Lubliner [35] introduced a power law velocity profile,
+2 r
vx = U 1 ( ) . (2.5)
R
Through this formula, and Cf are related in the equation
2
Cf = .
1
By varying the parameter , we can approximate the velocity profile under various
Womersley numbers. Fig. 2.3 plots the profile when is in the range from 2 to 9. We
note that when = 2, the profile is a parabolic or Poiseuille one. As increases, the
profile becomes closer to a flat one. When = 9, one gets = 1.1 and Cf = 22.
The latter value was fitted from cardiovascular blood flow [76] and then was adopted by
some following researches (e.g. [2, 48]).
19
In large arteries, predicted by Womersley theory, the central part of the profile is
essentially flat and there is a viscous boundary layer to match the no-slip boundary
condition. Olufsenp et al. [56] assumed such an essentially flat profile and estimated the
viscous layer by /, where is the kinematic viscosity and is the angular frequency
of the pulsatile flow. Noting that a phase difference between Q and wall shear stress is
predicted by Womersley theory, Bessems et al. [10] added a correction term based on
local Womersley number and P/x to the friction term. There are some other models
also based on the local Womersley number to approximate and the friction term, such
as [38, 101]. However, we do not elaborate them here.
Keeping in mind that should vary in the range between 1.0 and 1.3 in the whole
system, we take a uniform value everywhere the 1D model applies. The justification for
this approximation is twofold. First, in large arteries, the thin boundary layer contributes
very little to the integration of the convection term, thus must be very close to 1.0.
Second, in smaller arteries (where is closer to 4/3), the whole convection term becomes
relatively less important because the flow rate decreases and pressure drops more sharply
(the convection term is even neglected to linearize the system in some studies). Thus
we can set to 1.0 or a slightly bigger value (1.1 for example), and this approximation
sacrifices very little of modeling accuracy. For the friction term, however, we consider
Cf depends on the local Womersley number (or the local radius since and are the
same everywhere). The phase difference between Q and the friction is not considered
because the value is vary small [86] and no evidence of clinical relevance about this small
quantity is found. In one word, can be fixed to a value very close to 1.0 and Cf varies
across the whole vascular system, whose value must be fitted from observations.
We collect the final form of 1D governing equations in the following:
A Q
+ = 0, (2.6a)
t x
Q Q2 A P Q
+ ( ) + = Cf . (2.6b)
t x A x A
If we assume = 1, and replace Q by AU , the governing equations can be expressed in
another form,
A (AU )
+ = 0,
t x
U U 1 P U
+U + = Cf .
t x x A
Both of those two forms of governing equations are often seen in literature. In this thesis,
we choose to use the form in (A, Q, P ) because it has a looser limitation on the velocity
profile.
20
1
0.5
Radius [-]
0 =9 =2
-0.5
-1
0 0.5 1 1.5 2
vx [-]
Figure 2.3.: Velocity profile described by the formula proposed by Hughes and
Lubliner [35] when = 2, 3, . . . , 9.
viscoelastic solids in general), when a fixed stress is loaded, the wall keeps extending
gradually (creeping) after an instantaneous extension. If we denote the strain as and
the stress as , the viscoelastic behaviour can be described as
Z
(t) = J(0)(t) + j( )(t )d,
0
There are several other similar viscoelastic models as presented in [84]. When those
models are integrated into a 1D blood flow model, the difference on the prediction
of wave form is minor [61, 77]. A disadvantage of those models is the complexity of
numerical integration incurred by the convolution when the nonlinear 1D models are
solved.
Another widely used viscoelastic model is Kelvin-Voigt model (see Fig. 2.4), which
writes,
d(t)
(t) = E(t) + , (2.8)
dt
21
Figure 2.4.: Kelvin-Voigt viscoelastic model. The coefficients for the spring and dashpot
are E and respectively.
where is a coefficient for the viscosity of the material. The complex modulus for this
model is
E = = E + i.
We note that when = 0, this reduces to a pure elastic model. This model is simpler
when it is integrated with the 1D fluid model. However, previous studies show that
this model captures the viscoelasticity of arteries quite well [2, 4]. Thus we adopt the
Kelvin-Voigt model in this thesis.
We consider a cross section of a vessel with a thin wall (h << R) as shown in Fig.2.5.
We also assume that the arterial wall is isotropic, homogeneous, incompressible, and
moreover that it deforms axisymmetrically with each circular cross-section independently
of the others. The circumferential strain can be expressed as
R R0
= , (2.9)
(1 2 )R0
where R0 is the reference radius without loading and is the Poisson ratio, which is 0.5
for an incompressible material. By Laplaces law, the transmural difference between the
internal pressure P and the external pressre Pext is balanced with the circumferential
stress in the relation
h
P Pext = . (2.10)
R
Combining Eq. 2.8, 2.9 and 2.10, we can easily write the tube law in terms of A and P
A
P = Pext + ( A A0 ) + s , (2.11)
t
with the stiffness coefficient ,
Eh
= ,
(1 2 )A0
and the viscosity coefficient s ,
h
s = . (2.12)
2(1 2 ) A0 A
As
For convenience, we further define Cv = for reasons which will be clear very soon in
the next section.
22
Figure 2.5.: Illustration of the vessel geometry. F is the circumferential tension, h the
thickness of the wall, R the internal radius, P the internal pressure and Pext
the external pressure.
and
0
S=
A ( A0 )
.
Cf Q
A + x 2
3 A
x
In this equation, U is the conservative variable, F the corresponding flux and S the
source term. Note that the flux (scaled by constant density) consists of two parts, the
23
Q2 3
convective Fc and the diffusive Fv . We recognize A due to the fluid flow, 3 A
2 due
to the elasticity, andCv Q
x due to the viscosity of the wall. In general, the suitable
numerical techniques for the convective and diffusive fluxes are different. Thus it is
common to separate the diffusive term and put it on the right side. Thus we may write
the problem in a convection-diffusion form:
U F
+ =S+D (2.15)
t x
with
0
F = Fc , D= 2 .
Cv xQ2
We consider firstly the homogeneous part and later the non-homogeneous part. Expand-
ing the derivative of the flux, the homogeneous part can be written in a quasi-linear form
U U
+ Jc = 0, (2.16)
t x
where Jc is the Jacobian matrix
!
0 1
Jc = Q2
A2
+ c 2Q
2
A
Actually, A is always positive. Therefore c is real, which is the speed of the pressure
wave with respect to the fluid flow. The matrix Jc has two different eigenvalues
Q
1,2 = c. (2.18)
A
Linear algebra shows Jc must be diagonalizable in the form Jc = RR1 . The columns
of R are the right eigenvectors of Jc . Left multiplying Eq. (2.16) by R1 , one obtains
U U
R1 + R1 RR1 = 0.
t x
By introducing a new vector which satisfies U W = R1 , the previous equation can be
transformed into
W W
+ = 0. (2.19)
t x
W1,2 can be readily obtained by integrating U W = R1 componentwise
Q
W1,2 = 4c. (2.20)
A
24
W = [W1 , W2 ]T is called Riemann invariant vector or characteristics. In time-space
plane, W1,2 are constants along the lines Dt X1,2 (t) = 1,2 . In physiological conditions,
1 > 0 > 2 . The two families of characteristic propagate in opposite directions. The
homogeneous part is a subcritical hyperbolic system. For further use, we get the expres-
sions for A and Q by inverting the relation (2.20),
2
(W1 W2 )4 W1 + W2
A= , Q=A . (2.21)
1024 2
In the non-homogeneous part, the skin friction term dissipates the momentum and the
second order derivative of Q is diffusive. Thus the full system has hyperbolic-parabolic
features. In physiological conditions, the Womersley number is not too big and the artery
is almost uniform, thus the source term will be very small and the system is dominated
by the hyperbolicity feature. If the properties of the artery have sharp variations, large
source terms will be introduced. In this case, we will treat the artery as different segments
connected together.
W1 W1
+ 1 (U ) = 0, (2.22a)
t x
W2 W2
+ 2 (U ) = 0. (2.22b)
t x
Since the two eigenvalues have opposite signs, there is exactly one incoming charac-
teristic at each end of the computational domain. The incoming characteristic carries
25
information from outside of the domain and thus is essential to guarantee the problem
to be well-posed. That is to say, the system must be supplemented by B.C.s in the form
The outgoing characteristic carries information from inside of the domain, which
can be given by the differential equations. Since W1,2 are constants along the lines
Dt X1,2 (t) = 1,2 in time-space plane, we can get W2n+1 (0) and W1n+1 (L) by interpola-
tion in the data of the n-th time step:
The characteristics are then transformed to physical variables by relation (2.21) for
numerical computation.
In reality, we rarely have the explicit expression (2.23) for the incoming characteristics.
Usually, we want to impose B.C. in physical term A, Q or P . At the inlet, if An+1 is
given, one can use the relation (2.20) to deduce:
s
n+1 n+1 n+1
W1 = W2 + 8 A .
2
Qn+1
W1n+1 = W2n+1 + 2 .
An
If P n+1 is given, from the wall relation (2.11) simplified with no viscous effect (s = 0),
we in fact impose: r
n+1 n+1 1 1/2
W1 = W2 + 8 (P n+1 + A0 ).
2
At the outlet, some part of the perturbation of outgoing characteristic W1 may be
reflected,
W2n+1 = W20 Rt (W1n+1 W10 ),
where Rt is the coefficient of reflection. If Rt = 0, the B.C. is nonreflecting. That means
the outgoing characteristic goes out without leaving any effect and that the incoming
characteristic is a constant in time. If there are changes of properties in the downstream
of the vessel, usually a nonzero Rt will be incurred.
26
Figure 2.6.: Illustration of a branch. At the conjunction point, there are six unknowns:
Ap and Qp at the outlet of the parent vessel; Ad1 , Qd1 , Ad2 and Qd2 at the
inlets of the daughter vessels.
In these cases, the derivatives of the corresponding variables in the source terms are
very large or even near a singularity, and then the vessel can be treated as several joined
segments with different properties.
Since all of the conjunction points can be treated with the same method, we consider
a branching point as a sample problem: a parent vessel with two daughter arteries (see
Fig. 2.6). At the branching point, there are then six boundary conditions, An+1 p and
n+1 n+1 n+1 n+1 n+1
Qp for the outlet of the parent artery and Ad1 , Qd1 ,Ad2 and Qd2 for the inlets
of the two daughter arteries. From the physical point of view, we have to preserve the
conservation of mass flux
Qn+1
p Qn+1 n+1
d1 Qd2 = 0, (2.25a)
and conservation of momentum flux
n+1
1 Qn+1
p 1 Qdi 2
( n+1 )2 + Ppn+1 ( n+1 ) Pdn+1 = 0 i = 1, 2. (2.25b)
2 Ap 2 Ad i
i
27
The same principle holds for W2 on the two daughter arteries,
Combining Eqs. (2.25a), (2.25b), (2.25c) and (2.25d), there are 6 Eqs. with 6 unknowns.
This nonlinear algebraic system can be readily solved by Newton-Raphson iterative
method with U n as the initial guess. In our test, the computation converges very fast.
Usually a very few iterations are enough for a satisfactory accuracy.
2.5. Conclusion
In this chapter, we derived a viscoelastic 1D model by integrating a simplified conserva-
tion equation and Navier-stokes equqtions along the radius of a vessel. The treatment of
axial velocity profile and tube law was discussed in detail. The final governing equations
of the 1D model form a hyperbolic-parabolic system. The characteristic structure of the
hyperbolic part was presented. Noticing that the hyperbolic part dominates the whole
problem, we treat the boundary conditions with the characteristic method.
28
3. Numerical Schemes
3.1. Introduction
In case of weak nonlinearity (i.e. small perturbation around the equilibrium state [42,
59]), we can linearize the 1D governing equations and find analytical solutions in fre-
quency domain [53, 91]. But for the full nonlinear system, analytical solutions are not
available yet. Thus several numerical schemes have been proposed and used to solve the
system in time domain. We roughly classify them in:
Finite Difference (FD) [22, 56, 64, 67, 76, 80, 102]
29
the mesh size. But if a diffusive term is added to the governing equations, the term will
be hard to treat by an implicit time marching method (e.g. Crank-Nicolson) in the DG
setting, thus the time step may be very severely limited. Therefore, the performance of
each scheme depends on the main features of the studied problems. In fact, the problems
with different main features arise in a wide range of applications. For instance, no shock
is observed in arteries in normal physiological conditions but shock-like phenomena may
arise in veins [12, 24, 47] or in arteries when the human body suffers from a blunt impact
by accident [36]. For another instance, in some conditions diffusive terms or dispersive
terms may arise as source terms [2] and the proper treatment of these terms will pose dif-
ferent levels of difficulty in each numerical framework. Thus to make a cross comparison
of the numerical schemes is interesting and useful.
In this chapter, Section 3.2 describes the numerical solvers. In particular, a large
amount of details of computation are given because this kind of information is scat-
tered in literature. First the hyperbolic-parabolic problem is split into a hyperbolic
subproblem and a parabolic one. Then the hyperbolic subproblem is integrated by Mac-
Cormack, Taylor-Galerkin and MUSCL schemes and the parabolic subproblem is treated
by a Crank-Nicolson method. At the end of this section, a local discontinuous Galerkin
method is presented as an alternative approach without operator splitting. Section 3.3
shows the analytical solutions and numerical results of the proposed schemes. The sys-
tem is linearized and asymptotic solutions are obtained with different source terms in
the system. The effects of the skin friction and the viscosity of the wall on the pulse wave
are clearly observed. Moreover, a wave with a step jump is computed and the ability of
the four schemes to properly capture the shock-like phenomena is tested. After that, a
simple bifurcation is computed and the numerical reflection and transmission coefficients
are compared with the analytical ones predicted using linearized equations. Finally, a
network with 55 arteries is computed. All the numerical solutions are compared favor-
ably with the analytical, semi-analytical solutions or clinical observations. In the last
section, comparisons among the four schemes are made in four important aspects: accu-
racy, ability to capture shock-like phenomena, computational speed and implementation
complexity. The suitable conditions for the application of each scheme are discussed.
30
integration. Thus we applied a fractional step or operator splitting method. Starting
from Eq. (2.15), the original problem is split into to a hyperbolic subproblem,
U F
+ =S (3.1)
t x
and a parabolic one,
U
= D. (3.2)
t
Let us consider the time intervals (tn , tn+1 ), for n = 0, 1, ..., with tn = nt. In every
time interval, the hyperbolic problem is solved to get a predictor U , which is used as
the initial condition (I.C.) of the second problem. The second step can be viewed as a
corrector. The original problem is approximated by a sequential application of the two
subproblems in a certain order.
From data U n , we may make a prediction U by evolving time t of the hyperbolic
subproblem, and correct it with the evolution over t of the parabolic subproblem,
etH etP
U n U U n+1 ,
where etH (etP ) means to solve the hyperbolic (parabolic) subproblem over t. This
method is called Godunov splitting. If the two subproblems are not commutable, the
splitting error is O(t), see Chapter 17 of reference [41].
There is a 3-stage splitting called Strang splitting, which has a leading error term
O(t2 ),
1 1
e 2 tP etH e 2 tP
U n U U U n+1 .
But in most cases the errors induced by the two splittings are very close. That is because
the coefficient of the term O(t) is much smaller then the coefficient of O(t2 ) [41].
We will see in Section 3.3.3 a test case on the diffusion term. The results show that the
Godunov splitting is sufficient for our problem.
at the time step t + t. The finite difference equations (at the interior grid points) are
then :
1. predictor step
t n
Ui = Uin (F Fin ) + tSin , i = 2, ...N
x i+1
31
2. corrector step
1 t t
Uin+1 = (Uin + Ui ) (F Fi1
)+ S , i = 2, ...N
2 2x i 2 i
where F and S are evaluated as functions of the predicted solution U . Note that
the predictor step applies a forward differencing and the corrector step a backward
differencing. The order of the two kinds of differencing can be reversed. The grid points
x1 and xN +1 represent the boundary conditions.
32
The piecewise linear function space associated with the mesh (Figure 3.1) is given as
This is both the trial function space and the test function space in Galerkin framework.
We further define the inner product
Z L
(U, V ) = U V dx.
0
(An+1 n
h , vi ) = RHS1i ,
(Qn+1 n
h , vi ) = RHS2i ,
the terms F (Uhn ), FLW (Uhn ), SLW (Uhn ), SU (Uhn ) and H(Uhn ) are projected onto the trial
function spaceP and expanded by a group finite element method. That is, for example,
[F (Uhn )]1 = j=N j=2 [F n ] v with [F n ] = [F (U n )] . Finally, the matrix form of Eq. (3.7)
j 1 j j 1 j 1
writes
t2
MAn+1 = MAn + tKT [FnLW ]1 (M1 [Fn ]1 + M2 [Fn ]2 )
2 (3.8)
t2
(K1 [Fn ]1 + K2 [Fn ]2 ) + tM[SnLW ]1 ,
2
where
vj
Mij = (vi , vj ), Kij = (vi , )
x
33
and
X X
vi vi
M1 (Su )ij = (Su(1,1) )k vk , vj , M2 (Su )ij = (Su(1,2) )k vk , vj ,
x x
k k
X X
(1,1) vi vj (1,2) vi vj
K1 (H)ij = Hk vk , , K2 (H)ij = Hk vk , .
x x x x
k k
(,)
The form (Su )k indicates the k-th component of the vector at the position (, ) of the
discretized matrix Su . Please note that the operators M1 etc. are functions of Su and
H, therefore they must be updated in every time step.
3.2.4. MUSCL
In this section, we mainly follow the presentation [19] but with a different temporal
integration method. For finite volume method, the domain is decomposed into finite
volumes or cells with vertex xi as the center of cell [xi1/2 , xi+1/2 ], see Figure 3.2. In
every cell, the conservation law must hold,
Z xi+1/2 Z xi+1/2 Z xi+1/2
U F
dx + dx = Sdx.
xi1/2 t xi1/2 x xi1/2
34
where 1 is the biggest eigenvalue of Jc . Other numerical fluxes with less numerical
diffusivity are possible, such as Harten-Lax-van Leer (HLL) flux [11, 19]. Since Rusanov
flux is more simple and robust, it is adopted in this paper. If U and U+ are equal to the
average values at the cells, the scheme will be of first order accuracy. Reconstructions
of U and U+ from U are necessary for a scheme of higher resolution.
Let us consider the techniques of reconstruction. For a scalar s within the i-th cell, we
denote its slope as Dsi , which can be approximated by (si si1 )/x, (si+1 si )/x
or (si+1 si1 )/2x. Then the values of s at the interfaces associated with this cell can
be recovered as
x x
si1/2+ = si Dsi , si+1/2 = si + Dsi .
2 2
The discretization of derivative in space can achieve a second order accuracy by this
method. But the solution will have nonphysical oscillations. Some examples of oscil-
lations induced by these methods can be found in Chapter 6 of reference [41]. Slope
or flux limiter and non-oscillatory solutions are integral characteristics of FV schemes.
MUSCL (monotonic upwind scheme for conservation law) is one popular slope limited
linear reconstruction technique. To present MUSCL, we first define a slope limiter,
min(x,y) if x, y 0,
minmod(x,y) = max(x,y) if x, y 0,
0 else
35
The numerical fluxes Fi+1/2
and Fi1/2 are given by Rusanov flux with the reconstructed
values at the two sides of the interfaces. Note that this is a scheme with five stencils.
The values at x1 and xN +1 are determined by the aforementioned characteristic method.
One ghost cell at each end of the computational domain is needed and we approximate
the values at these cells by the ones at the neighboring boundary cells.
For the temporal integration, we may apply a 2-step second order Adams-Bashforth
(A-B) scheme,
3 1
Un+1 = Un + t (Un ) (Un1 ) .
2 2
This scheme can be initiated by a forward Euler method. Also, a second order Runge-
Kutta (R-K) approach, namely Heun method is possible [75]. It writes
U = Un + t(Un ),
U = U + t(U ),
Un+1 = (U + U )/2.
Comparing the two methods, we note that (U) has to be computed twice in R-K
in every time step while the A-B method only needs once since (Un1 ) is stored in
the previous step and reused in the current step. Because the boundary conditions are
determined dynamically to compute (U), the R-K also incurs one more resolution of
the nonlinear algebraic equations at conjunction points. For these reasons, we choose
the A-B method for the temporal integration, although the R-K method usually allows
a larger time step for convergence.
36
order accurate both on time and space. Moreover, it is unconditionally stable. It is
natural to set a homogeneous Neumann B.C. for the parabolic subproblem, x Up (0, t) =
x Up (L, t) = 0. The subscript p stands for parabolic. We note that a second-order
implicit FE method can also be applied here. But since this subproblem is linear and is
in 1D, the FE method would be exactly equivalent to this FD method.
Similarly, the flux F and the source term S can also be approximated by the direct
summation of piecewise N -th degree polynomials. The local form of the conservation
law on the k-th element is
Uhk Fhk
+ = Shk . (3.12)
t x
Multiplying both sides of Eq. (3.12) with a test function k , and integrating over one
element give k k
Uh k Fh k k k
, + , = Sh , . (3.13)
t Dk x Dk Dk
At the interface of xk , the values of Uh at the two sides, Uhk1 (xk ) and Uhk (xk ), are
not guaranteed equal. A numerical flux Fk is introduced here. Through the numerical
flux, information is communicated between elements. In practice, the second term is
integrated by part again for convenience of computation. Thus we have
k k xk+1
Uh k Fh k k k
k k
, + , + (Fh + F ) = Sh , . (3.15)
t Dk x Dk xk Dk
37
If we introduce Np nodes within the element Dk (Figure 3.3), the local solution can be
expanded as
Np
X
k
Uh (x, t) = Uhk (xki , t)`ki (x), (3.16)
i=1
where `ki (x) is the Lagrange interpolant associated with the i-th node. For the Galerkin
scheme, Eq. (3.15) must hold for every test function `ki (x). Thus we have Np equations
for Np unknowns. In matrix form, the system can be written as,
k
xk+1
k dU k k k k
= M k Sk ,
M + K F + ` (Fh + F ) (3.17)
dt xk
where
d`kj
Mk(i,j) = `ki , `kj k
= `ki ,
Dk
, K(i,j) ,
dx Dk
and `k is the vector of functions (`k1 , `k2 , ..`kNp )T . The system of equations can be turned
into a semi-discrete form,
xk+1
dUk k k k 1 k k
+ Sk ,
= D F + (M ) ` (Fh F ) (3.18)
dt xk
where
d`kj
k k 1 k
D(i,j) = (M ) K =
(i,j) dr ri
is the local differentiation operator [31]. The computation of Mk and Dk is crucial. We
define an affine mapping from a reference element (1, 1) to Dk ,
1+r
x(r) = xk + (xk+1 xk ).
2
The local operators can be readily computed as
Z 1
k k 1 d`j
M(i,j) = Jk `i `j dr, D(i,j) = Jk ,
1 dr ri
where Jk = (xk+1 xk )/2 and `i , `j are the Lagrange interpolants at the reference ele-
ment. Note that the operators Mk and Dk can be precomputed and stored. Legendre-
Gauss-Lobatto points have to be chosen as the interpolation points to minimize compu-
tation error. For more details, we refer to Chapter 3 of reference [31]. For the temporal
integration, a second order A-B scheme is applied for reasons as discussed in Section 3.2.4.
The scheme previously presented can treat a hyperbolic problem. But in this setting
Crank-Nicolson method is hard to apply, because the values at the interfaces are du-
plicated. We consider the problem formulation of Eq. (2.14), where the flux contains
convective part Fc and diffusive part Fv . For the convective part, Rusanov flux as men-
tioned in Section 3.2.4 is applicable. For the diffusive flux, a straight idea is to use the
38
central flux, (Fv (U ) + Fv (U+ ))/2. But as pointed out by Shu el al. [74], this choice is
inconsistent.
To solve this problem, we rewrite the original equations as
U (Fc Cv q)
+ =S
t x
Q
q =0
x
In semi-discrete form, the equations for one element are
xk+1
dUk k k k 1 k k
+ Sk
= D F + (M ) ` (Fh F )
dt xk
xk+1
qk = Dk Qk (Mk )1 `k (Qkh Q )
xk
39
asymptotic solutions are obtained. In case of larger perturbations, the full nonlinear
system allows shocks. The shock-capturing property of each scheme is tested in this case.
After the tests on a single vessel, a simple bifurcation is computed and the reflection
and transmission coefficients are compared with analytical ones predicted by linearized
system. At the end of this section, a network with 55 arteries is computed and the
numerical solutions are checked against clinical observations reported in literature.
A Q Q A Cf 2 Q
+ = 0, + c20 = Q + Cv 2 (3.19)
t x t x A0 x
q
with c0 = 2 A0 , the Moens-Korteweg celerity. To investigate the propagation phe-
nomena at first, we drop the non-homogeneous part (Cf = 0 and Cv = 0). Then
Eqs. (3.19) become dAlembert equations, which admit the pure wave solution. We as-
sume that the initial condition is at equilibrium and the inflow is prescribed as Q(0, t) =
Qin (t) with
2 Tc
Qin (t) = Qc sin( t)H(t + ), t > 0,
Tc 2
where H(t) is the Heaviside function, Tc the period of the sinusoidal wave and Qc the
amplitude. The solution is c0 A = Q = Qin (x c0 t), which means that the waveform
propagates to the right with a speed of c0 .
We propose a numerical test with parameters of the tube inspired by [71]: L = 250cm,
A0 = 3.2168cm2 , = 1.8734 106 Pa/m, = 1.050 103 kg/m3 , and accordingly c0 =
400cm/s. To impose a small perturbation, we choose Qc = 1ml/s and Tc = 0.4s. In
this case the change ratio of the radius is R/R0 = Qc /(2A0 c0 ) = 0.04%, thus the
perturbation is assured small enough. We take the linearized analytical solution at time
t = 0.4s as reference to compute the errors of the numerical solutions. The normalized
error is defined by ||E|| = ||Qnumerical Qanalytical ||rms /Qc , where || ||rms stands for
the root-mean-square error. To specify the time step, we note that it first should satisfy
the CFL (Courant-Friedrichs-Lewy) condition which writes
N +1h hi i
t 6 nCF L min Qi+1
,
i=0 max( Qi
Ai + ci , Ai+1 + ci+1 )
40
nCF L = 12 [19]. Practice shows that nCF L = 1 for MacCormack scheme [22]. A sharp
estimation of the coefficient nCF L for the DG scheme is challenging. We define an
approximate formula, t = Ct NLc0 , to test the stability. In our test, the approximate
threshold values of Ct for the schemes to become unstable are: 0.5 for MUSCL, 0.56 for
Taylor-Galerkin and 1.0 for MacCormack. The results agree with the report in literature.
For the DG scheme, the time step formula is modified accordingly as t = CPt NLc0 , with
P the degree of the polynomial. For the DG scheme, Ct can not be greater than 0.1 (see
Figure 3.4(b)).
To further test the temporal convergence, we fix the mesh (NT G = NF V = NF D = 800,
NDGP1 = NDGP2 = 100) and plot the numerical errors as a function of Ct (see
Figure 3.4(a)). The errors vary slightly for all of the schemes except MUSCL. For the
convergence of the temporal integration, the MUSCL scheme has to choose a smaller
time step than the value prescribed by the CFL condition. But note this is only a
test in linear case, in practical applications, the coefficient Ct may be much smaller for
convergence (Section 3.3.6).
To test the spatial convergence, we fix Ct = 0.1, and vary the number of mesh nodes
N . The log-log plot of ||E|| against x can be seen in Figure 3.4(c). We have two
main observations. First, all of the schemes converge with an order between 1 and 2 and
the DG scheme converges faster (see Figure 3.4(c)). Second, as shown by Figure 3.4(d)
the differences between the analytical solution and all of the numerical solutions are
hardly discernible with a moderate number of mesh points (NT G = NF V = NF D = 800,
NDGP1 = NDGP2 = 100).
To compare the actual speed and accuracy of the four schemes, we set N and Ct (see
Table 3.1) such that the errors achieve the same order of magnitude (see Figure 3.4(f)).
Except the Taylor-Galerkin scheme, all the schemes have the similar accuracy with very
close running time (see Figure 3.4(e) and 3.4(f)). At this point, the Taylor-Galerkin
scheme shows the worst accuracy and needs to run the longest time. We note that large
global matrices arise in Taylor-Galerkin scheme while the operators in other schemes
are local and have small size. That explains the relative poor performance of Taylor-
Galerkin even though a larger time step is allowed by this scheme. We will see that in
case of a network of real size, the largest number of N is about 100 and Taylor-Galerkin
shows a good balanced property between accuracy and speed (Section 3.3.6).
scheme N Ct
Taylor-Galerkin 800 0.5
MUSCL 800 0.3
MacCormack 1600 0.5
DG-P1 200 0.1
DG-P2 100 0.1
41
3
x 10 3
x 10
TaylorGalerkin
MUSCL
7 MacCormack 1.8
1.6
6
DGP1
DGP2
1.4
5
||E||
1.2
||E||
4
0.8
2
0.6
1
0.4
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
Ct Ct
(a) (b)
1.2
TaylorGalerkin
TaylorGalerkin MacCormack
MUSCL MUSCL
MacCormack 1 DGP1
DGP2
DGP1
2 DGP2
10
0.8
0.6
||E||
Flux
3
10
0.4
0.2
4
10 0
1 0 1 0.2
10 10 10 0 50 100 150 200 250
x distance(cm)
(c) (d)
3
x 10
1.5
1200
1000
800 1
Running time(s)
||E||
600
400 0.5
200
0 0
TaylorGalerkin MUSCL Maccormack DGP1 DGP2 TaylorGalerkin MUSCL Maccormack DGP1 DGP2
(e) (f)
Figure 3.4.: Test on a uniform tube. Top left and right: With a fixed mesh (NT G =
NM U SCL = NF D = 800, NDGP1 = NDGP2 = 100), errors as functions of
coefficient Ct . Middle left: Errors as functions of the sizes of elements (cells).
Middle right: All the numerical42solutions for the pulse wave at time 0.4s are
overlapped, the analytical solution is indicated by cross signs. Bottom left
and right: Running time and error of each scheme for the configuration
shown in Table 3.1.
3.3.2. Attenuation due to the viscosity of blood
We now consider the same linearized Eq. (3.19) with the small term due to skin friction
(Cf 6= 0 and Cv = 0). The main dynamics of the system will be grossly the same
traveling wave but attenuated by viscosity of blood. This behaviour can be predicted by
asymptotic analysis. We have a small non-dimensional parameter f = Tc Cf /A0 , which
is the ratio of the characteristic time of pulse Tc to the characteristic time of attenuation
A0 /Cf . In order to see how the waveform slowly evolves when it propagates to, say
right, we make a change of variables to = f t and = x c0 t (slow time, moving
frame). The two differential operators t and x expand as
= + = f c0
t t t
= = .
x x
The solution has the asymptotic expansion
A = A0 + f A1 + ..., Q = Q0 + f Q1 + ...
Substituting these into the governing equations expressed in new variables and collecting
the terms with the same order of f , one has
A0 Q0 A0 A1 Q1
(c0 + ) + f ( c0 + ) + .. = 0
Q0 A0 Q0 Q1 A1 Q0
(c0 + c20 ) + f ( c0 + c20 + ) + .. = 0.
Tc
We take the first order term in f in the first equation, substitute it in the first order
term in f in the second equation. Then we obtain
Q0 A0 Q0
( + c0 + ) = 0.
Tc
From the terms of the zeroth order in f , which involve derivative in only, the solution
must have the form Q0 = c0 A0 (, ) + ( ). Substituting it into the previous equation
generates terms and ( ). These are secular terms and thus can be set null. So we
Q0
have c0 A0 = Q0 and = 2T1 c Q0 , or
For more on asymptotic analysis of blood flow in large blood vessels, we refer to refer-
ence [100].
In Figure 3.5, we plot the snapshots of the waveform at time 0.2s, 0.4s, 0.6s and
0.8s. In the computation, the inflow is a half sinusoidal flux as described in the previous
subsection and the outflow is nonreflecting. The skin friction coefficient Cf is 40, and
43
the parameter 2A0 c0 /Cf is about 2000cm. The damping rate of the amplitude of the
C x
waveform agrees very well with the analytical prediction, exp( 2Af0 c0 ), which is indicated
by the dashed line. Also note that the errors of different schemes are not the same. The
MUSCL scheme causes the peak of the wave to slightly flatten, while all of the other
schemes are dispersive: we have small oscillations at the foot of the signal.
1.2
TaylorGalerkin
MUSCL
MacCormack
1 DG
0.8
0.6
Flux
0.4
0.2
0.2
0 50 100 150 200 250 300 350
Distance(cm)
Figure 3.5.: Attenuation due to the skin friction. The snapshots are at time 0.2s, 0.4s,
C x
0.6s and 0.8s. The dashed line is exp( 2Af0 c0 ) with 2A0 c0 /Cf ' 2000cm.
The flux is normalized with respect to Qc .
Q0 c2 Tc 2 Q0
= 0 . (3.20)
2 2
The solution of this equation can be given by the convolution
Z +
Q0 (, ) = Q0 (0, )G(, )d
and Q0 (0, ) is the initial state. In the test vessel, the parameters are kept the same
as in the case of attenuation. The coefficient Cv is 0.6275m2 /s and v ' 0.1. This
44
corresponds to =5000Pa s, which is in the range of observed values on animals [4].
To facilitate the calculation of the analytical solution, nonreflecting B.C.s are imposed
at the two ends of the vessel and the I.C. is a half sinusoidal waveform for Q (dashed
line in Figure 3.6) and a constant value for A0 . It is clear that half of the initial wave
propagates to right and at the same time the waveform is spread out due to the diffusive
effect. The analytical solution at time 0.4s (indicated by cross signs) agrees well with
the corresponding numerical solutions.
Another point worthy noticing is the operator splitting errors. In the DG scheme, no
operator splitting error is induced. All of the other numerical schemes adopt operator
splitting method. They produce very accurate solutions as well as DG. Thus it verifies
the a priori judgement that Godunov splitting is sufficient.
1.2
TaylorGalerkin
MUSCL
MacCormack
1 DG
0.8
0.6
Flux
0.4
0.2
0.2
0 50 100 150 200 250 300 350
Distance(cm)
Figure 3.6.: Diffusion due to the viscosity of the wall. The dashed line is the initial
condition. One half of the original waveform propagates to right. The
snapshots are at time 0.2s, 0.4s, 0.6s and 0.8s. The analytical prediction
from the convolution at time 0.4s is indicated by cross signs. The difference
between the different numerical solutions is not discernible. The flux is
normalized with respect to Qc .
Q0 1 Q0
= Q0
2A0
One important consequence of nonlinear hyperbolic system is that shocks may arise
even if the initial condition is very smooth. In normal physiological conditions, shocks are
45
not observed in arterial systems. But in venous system, shock-like phenomena may occur
on muscular veins during walking and running. The intramuscular pressure (equivalent
to Pext in our model) can rise to 20 40 kPa in a few milliseconds [6]. In such situation,
experiments and numerical simulations [24, 47] have shown this critical behaviour. For
some large mammals, for instance giraffes, even in static postures, the gravity-driven
flow in a long inclined vein may develop into shock-like waves, like the roll waves in
a shallow-water channel [12, 13]. For another example, the traumatic rupture of the
aorta is responsible for a significant percentage of traffic death and the rupture may be
well accounted for by the shock-like transition resulted from the blunt impact to the
thorax [36]. For possible applications in these situations, we test all the schemes with a
shock-like wave.
To generate a shock, we impose a step jump signal of flux at the inlet. For a vessel of 1
meter, the numbers of elements for Taylor-Galerkin, MacCormack and MUSCL schemes
are 100, 200 and 800 respectively. The DG scheme uses 25 elements and the order of
polynomial is 2. Figure 3.7 shows that the MUSCL scheme with a flux limiter captures
the shock without nonphysical oscillations, whereas the other numerical schemes cause
spurious oscillations. This verifies the total-variation-diminishing (TVD) property of the
MUSCL scheme. But the MUSCL is very diffusive at the shock, thus a very fine mesh
is required. For the DG scheme, limiters may be introduced as well to eliminate the
oscillations [31]. This remedy will be necessary for DG to be applicable on problems
with shocks. On Figure 3.7(b) we plot a case with some viscosity of the wall. The added
moderate physical diffusive term smoothens the wave and all the schemes give almost
the same result.
Zp1 (Zd1
1
+ Zd1
2
) 2Zp1
R= , T = , (3.21)
Zp1 + (Zd1
1
+ Zd1
2
) Zp1 + (Zd1
1
+ Zd1
2
)
where Zp and Zd are the characteristic impedance of the parent and daughter vessels [29,
59]. In Figure 3.8, for sake of illustration, the configuration of the branching and the time
profiles of pressure at two locations are shown. The amplitude is normalized with respect
to Qc = 1106 m3 /s = 1ml/s. For the parent vessel: = 2.3633106 Pa/m, A0 = 4cm2
and for each of the daughter vessels: = 6.3021 106 Pa/m, A0 = 1.5cm2 . The B.C.s
at the outlets of the daughter vessels are nonreflecting. Thus the reflected pulse wave
is generated at the conjunction point. According to the formula (3.21), R = 0.2603
and T = 1.2603. The pressure profiles at the points A and B agree very well with the
analytical predictions. All of the numerical schemes are compatible with this treatment
46
TaylorGalerkin TaylorGalerkin
MUSCL MUSCL
MacCormack MacCormack
1.2 1.2
DG DG
1 1
0.8 0.8
Flux
Flux
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0.2 0.2
30 35 40 45 50 55 60 65 70 75 80 30 35 40 45 50 55 60 65 70 75 80
Distance(cm) Distance(cm)
(a) (b)
Figure 3.7.: A shock in the system. A step jump signal of flux is imposed at the inlet and
a snapshot is shown. The left figure (a) shows that the MUSCL scheme with
a flux limiter captures the shock without nonphysical oscillations, whereas
the other numerical schemes cause spurious oscillations. The right figure (b)
shows that all the schemes give almost the same result for a system with a
moderate physical diffusive term.
of conjunction point. Note that in healthy arterial system, the related arteries of most
conjunctions are well matched such that there are essentially no reflections (R=0) at
the conjunctions [91, 57]. The purpose of the proposed configuration is just to test the
numerical schemes.
47
Figure 3.8.: Reflection and transmission of pressure wave at a branching point. The time
profiles of pressure at points A and B are plotted. The analytical reflection
and transmission coefficients are 0.2603 and 1.2603 (indicated by the dashed
line).
locations may cause a considerable variation. Nevertheless the inclusion of viscosity term
makes it possible to test the numerical schemes in a more realistic condition.
The peak value of the input flux Qc is 500 ml/s. This value is very close to the peak
flow rate at the root of aortic artery [64]. We choose mini=55 i i
i=1 (L /c0 ) as a reference length,
i i
with L the vessel length and c0 the linearized wave speed of the i-th artery. For a coarsest
i Li /ci
possible mesh, the number of elements (cells) of the i-th artery is Nbase = b mini=55 (L0 i /ci ) c,
i=1 0
where bc is the floor function. We computed the relative change of solutions when the
number of the elements (cells) is doubled. Figure 3.9 shows the relative change of the
solutions when the number of the elements (cells) is changed from 2Nbase to 4Nbase .
The relative change of a quantity (for example flux Q) with two meshes N1 and N2
is defined as ||QN1 QN2 ||rms /(Qmax Qmin ), where || ||rms is the root-mean-square
error as before, Qmax and Qmin are the maximum and minimum values within one heart
beat. Figure 3.9 shows that the changes of flux and pressure are less than 1.5% for
all of of the schemes except DG. Thus we plotted in Figure 3.10 the results computed
with mesh 2Nbase . The DG scheme is not tested in this manner because it is already
converged: results in Figure 3.10 show that there is no discernible difference between the
DG solutions with the others even with the coarsest possible mesh. In this computation,
the order of polynomial of DG is 1, thus the total number of free degrees is 2Nbase , which
Li
is equal to those of the other schemes. Time step is prescribed by t = Ct mini=55
i=1 ( N i ci ).
0
The coefficient Ct and the corresponding real time steps in the computation are shown
in Table 3.2.
In Figure 3.10 we plot the history profiles of flux and pressure at the middle of four
representative arteries. All of the numerical solutions agree very well. The main features
of the pressure and flux profiles reported in literature [64, 71] are observed. The peak
48
3 3
x 10 x 10
12 TaylorGalerkin TaylorGalerkin
MUSCL MUSCL
MacCormack MacCormack
5
10
4.5
4
6
3.5
4
3
2
2.5
5 10 15 20 25 30 35 40 45 50 55 5 10 15 20 25 30 35 40 45 50 55
Artery Artery
(a) (b)
Figure 3.9.: Relative changes of the solutions when the mesh is doubled from 2Nbase to
4Nbase . The left figure shows that the relative changes of all the fluxes are
less than 1.3 %. The right figure shows that the relative changes of all the
pressures are less then 0.6% .
Table 3.2.: Time steps and running time for one heart beat using one processor on a
standard Linux work station with MATLAB.
value of pressure waveform increases as we travel down the system. We can also see the
dicrotic notch at artery 1. At artery 37, a reverse flow is observed (see Figure 3.10(f)),
which agrees with clinical measurement [64].
Both in vivo [32, 64] and in vitro [2] studies show that the models with viscoelastic-
ity predict the pulse waves better. This effect is most pronounced at the peripheral
sites [2, 70]. The predictions by the elastic and viscoelastic models are compared at two
locations, see Figure 3.11. We can clearly see the smoothening effect on the pulse curves
at both sites. The biggest relative difference is observed on the flow rate curve at the
peripheral site (see Figure 3.11(c)). This study confirms again the necessity to consider
the viscoelasticity in the 1D model.
49
Flux Pressure
600 TaylorGalerkin
TaylorGalerkin
MUSCL
MUSCL
MacCormarck 16 1 MacCormarck
1 DG
DG
400
14
kPa
ml/s
200
12
0 10
14 8 MacCormarck 8 MacCormarck
DG DG
15
10
ml/s
kPa
12
6
2 9
7.2 7.4 7.6 7.8 8 7.2 7.4 7.6 7.8 8
seconds seconds
(d) (e)
TaylorGalerkin
TaylorGalerkin
MUSCL
MUSCL
120 37 MacCormarck 16 37 MacCormarck
DG
DG
80 14
kPa
ml/s
40
12
0
10
16
12
ml/s
kPa
13
8
10
4
7.2 7.4 7.6 7.8 8 7.2 7.4 7.6 7.8 8
seconds seconds
(a) (h) (i)
Figure 3.10.: The history profiles of flux and pressure at four locations. Ten heart beats
are computed to secure that steady state is achieved, but only the tenth
heart beat is plotted. The differences between the four numerical schemes
are very small. See Table 3.3 for time steps and running time of each
scheme.
50
elastic elastic
visco visco
120 37 16 37
80
14
kPa
ml/s
40
12
0
10
16 19
elastic elastic
visco visco
54 54
16
12
ml/s
kPa
13
8
10
4
7.2 7.4 7.6 7.8 8 7.2 7.4 7.6 7.8 8
seconds seconds
(c) (d)
Figure 3.11.: The comparison between elastic and viscoelastic models (MUSCL scheme).
The viscoelasticity damps the oscillations of high frequency.
51
3.4. Conclusion
In this chapter, we solved the hyperbolic-parabolic system with four numerical schemes:
MacCormack, Taylor-Galerkin, MUSCL and local discontinuous Galerkin. The imple-
mentations were verified with analytical, semi-analytical or clinical observations in many
cases. At first, a single uniform tube was considered. Under the assumption of small
nonlinearities, we obtained asymptotic solutions of the linearized system with different
source terms. The propagation, attenuation and diffusion of the waveform were illus-
trated by both the numerical and analytical solutions. Moreover, in case of a larger
nonlinearity, the shock capturing property of each scheme was tested. After the test
on a single vessel, a simple bifurcation was computed to check the numerical coupling
of different arteries. Finally, we computed a relatively realistic network with 55 ar-
teries. The check of the numerical solutions in all cases was very favorable for all of
the schemes. We can compare the schemes in four aspects: accuracy, shock-capturing
property, computational speed and implementation complexity.
3. For a network of human size, the speed of computation can be ordered from fast
to slow as: Taylor-Galerkin, MUSCL, MacCormack and local DG. The temporal
integration in the Taylor-Galerkin scheme is more efficient than Adams-Bashforth
2-step method. Thus it allows a larger time step with a comparable accuracy. But
if the number of elements for one artery is too large (larger than 500), Taylor-
Galerkin becomes slower because the sizes of the global matrices increase quadrat-
ically and thus the storing and inverting of matrices become very expensive. The
DG scheme prevents the application of Crank-Nicolson method on the diffusive
term. An explicit method called local DG scheme was adopted in this paper. Even
with a moderate diffusion coefficient (within the range observed in physiological
condition), a very small time step is necessary for stability. To compute one heart
beat, the local DG takes about 9 hours while all other schemes take only 20-90
minutes (using one processor on a standard Linux workstation with Matlab).
52
4. From easiest to hardest, the implementation of the schemes can be ordered: Mac-
Cormack, MUSCL, Taylor-Galerkin and local DG.
53
Table 3.3.: Arterial network
l A0 Cv
ID Name (cm) (cm2 ) (106 Pa/cm) (104 cm2 /s) Rt
1 Ascending aorta 4.0 6.789 0.023 0.352
2 Aortic arch I 2.0 5.011 0.024 0.317
3 Brachiocephalic 3.4 1.535 0.049 0.363
4 R.subclavian I 3.4 0.919 0.069 0.393
5 R.carotid 17.7 0.703 0.085 0.423
6 R.vertebral 14.8 0.181 0.470 0.595 0.906
7 R. subclavian II 42.2 0.833 0.076 0.413
8 R.radius 23.5 0.423 0.192 0.372 0.82
9 R.ulnar I 6.7 0.648 0.134 0.322
10 R.interosseous 7.9 0.118 0.895 0.458 0.956
11 R.ulnar II 17.1 0.589 0.148 0.337 0.893
12 R.int.carotid 17.6 0.458 0.186 0.374 0.784
13 R. ext. carotid 17.7 0.458 0.173 0.349 0.79
14 Aortic arch II 3.9 4.486 0.024 0.306
15 L. carotid 20.8 0.536 0.111 0.484
16 L. int. carotid 17.6 0.350 0.243 0.428 0.784
17 L. ext. carotid 17.7 0.350 0.227 0.399 0.791
18 Thoracic aorta I 5.2 3.941 0.026 0.312
19 L. subclavian I 3.4 0.706 0.088 0.442
20 L. vertebral 14.8 0.129 0.657 0.704 0.906
21 L. subclavian II 42.2 0.650 0.097 0.467
22 L. radius 23.5 0.330 0.247 0.421 0.821
23 L. ulnar I 6.7 0.505 0.172 0.364
24 L. interosseous 7.9 0.093 1.139 0.517 0.956
25 L. ulnar II 17.1 0.461 0.189 0.381 0.893
26 Intercostals 8.0 0.316 0.147 0.491 0.627
27 Thoracic aorta II 10.4 3.604 0.026 0.296
28 Abdominal aorta I 5.3 2.659 0.032 0.311
29 Celiac I 2.0 1.086 0.056 0.346
30 Celiac II 1.0 0.126 0.481 1.016
31 Hepatic 6.6 0.659 0.070 0.340 0.925
32 Gastric 7.1 0.442 0.096 0.381 0.921
33 Splenic 6.3 0.468 0.109 0.444 0.93
34 Sup. mesenteric 5.9 0.782 0.083 0.439 0.934
35 Abdominal aorta II 1.0 2.233 0.034 0.301
36 L. renal 3.2 0.385 0.130 0.481 0.861
37 Abdominal aorta III 1.0 1.981 0.038 0.320
38 R. renal 3.2 0.385 0.130 0.481 0.861
39 Abdominal aorta IV 10.6 1.389 0.051 0.358
40 Inf. mesenteric 5.0 0.118 0.344 0.704 0.918
41 Abdominal aorta V 1.0 1.251 0.049 0.327
42 R. com. iliac 5.9 0.694 0.082 0.405
43 L. com. iliac 5.8 0.694 0.082 0.405
44 L. ext. iliac 14.4 0.730 0.137 0.349
45 L. int. iliac 5.0 0.285 0.531 0.422 0.925
46 L. femoral 44.3 0.409 0.231 0.440
47 L. deep femoral 12.6 0.398 0.223 0.419 0.885
48 L. post. tibial 32.1 0.444 0.383 0.380 0.724
49 L. ant. tibial 34.3 0.123 1.197 0.625 0.716
50 R. ext. iliac 14.5 0.730 0.137 0.349
51 R. int. iliac 5.0 0.285 0.531 0.422 0.925
52 R. femoral 44.4 0.409 0.231 0.440
53 R. deep femoral 12.7 0.398 0.223 0.419 0.888
54 R. post. tibial 32.2 0.442 0.385 0.381 0.724
55 R. ant. tibial 34.4 54
0.122 1.210 0.628 0.716
Data adapted from [71] and [4].
4. Development of a parallel code
4.1. Introduction
Patient-specific modeling of blood flow in vascular networks is a very changing task.
There are billions of vessels in the systemic arterial network alone, and they vary largely
in scale. The diameter of vessels ranges from 30 mm (aorta) to 5 103 mm (capillary),
and the Womersley number ranges from 20 to 103 . The main features of blood flow
in various locations can be drastically different. For instance, in the large and medium-
sized vessels, the blood flow is convection-dominated and significantly pulsatile; in the
arteriole and capillary bed, the blood flow is dominated by fluid viscosity and pulsates
slightly in time. Models with various complexity from 3-dimensional to 0-dimensional
have been proposed and the coupling between them is also a hot research topic (multi-
scale modeling). Among those models, the 1D model plays a very special role due to its
balance between modeling complexity and computing cost.
Currently, simulation of a network with about 50 segments for one cardiac cycle can
be done in a few minutes [2, 16, 48]. Keeping this in mind, we want to develop a code
which will have the following features:
2. Easy to be set for large networks. If we only consider the part of the network
where a 1D model is applicable, the number of vessel segments is well beyond our
simulation capacity. Usually the network of the 1D model is truncated at certain
levels and a simplified model is designed to mimic the downstream sub-networks.
This simplified model may be a tapered long tube [51], a structured network of
linearized 1D models [56] or a 0D model [88]. This is usually a good compromise
between economy and efficacy. However, in many cases, the distributed 1D model
cannot be replaced by lumped models. For instance, to simulate some special cir-
culations (e.g. coronary, renal and cerebral circulations), we usually must include
a larger number of vessels in the 1D model to obtain detailed information on the
55
Figure 4.1.: A network described by Directed Acyclic Graph (DAG). Nodes represent
conjunctions, directed edges represent vessels.
distribution of blood flow. Thus, our code should be easily scaled to such large
networks.
56
Figure 4.2.: Input/Output of the code. DAG.csv specifies topology of the network; ves-
sels.csv specifies physical parameters, meshes etc.; outlet.csv is for outlet
boundary conditions, output.csv specifies the locations of interests to be
recorded; time setting.csv is for the time control of computing and recording.
of the arteries. The properties of each vessel are also set through CSV files, such as
the viscoelasticity, the mesh, the reference cross sectional area etc. In the same way, we
specify the size of discretized time step, the location and time intervals to keep records
of the evolutionary history of the system, etc. As shown by Fig. 4.2, the C++ code
reads those files, makes computations and writes the results to specified destinations
for further processing. The computing is done in parallel to improve speed if there are
multi-segments in the network. For convenience, the parameter files can be prepared by
a scripting language. In appendix A, using a circle of arteries as a computing case, we
show more details of the structure of CSV files. We also briefly talk about how to prepare
the parameter files with Python scripts; how to track errors of parameter files; and how
to make animations of the computing result with Gnuplot as the plotting engine. In
fact, the whole workflow of a simulation, including computing and data analysis, can be
done conveniently only with free open source softwares.
57
Figure 4.3.: Illustration of the hierarchy of artery. The class artery defines a virtual
method stepMarch, which is implemented in the subclasses artery FV
and artery FD with different numerical schemes.
Figure 4.4.: Illustration of data structure of a network. Each conjunction contains one
vector of artery pointers to parent vessels and another one to daughter
vessels.
58
At the interior of each artery, the state evolves independently and there are only
interactions between connected arteries in the conjunctions. Thus this problem is very
suitable for parallelization. We take an coarse-grain strategy to parallelize the computing
with OpenMP. Listing 4.1 shows one simplified segment of the code. From line 3 to 16,
the whole computational domain is divided evenly among different threads. From line
19 to 26, the conjunctions are iterated to update boundary conditions of each artery. If
the conjunction does not have any parent arteries, the method isRoot returns true and
we impose an inflow boundary condition. Similarly, if the conjunction only has parent
arteries, the method isLeaf returns true and outflow boundary conditions are imposed
(specified by reflection coefficients). If both the two checks fail, the conjunction must
locate in the interior of the DAG graph and stepMarch is called to solve the system
of coupled equations. In line 27, there is an barrier to synchronize all the conjunctions.
When all the new boundary values of the arteries have been imposed, the next for
loop starts to update the interior values within each artery. The recordings of quantities
designated by the user are also done in parallel manner. They are not shown here for the
reason of space. By this coarse-grain parallelization, the overhead for spawning threads
are minimized. Its performance is better in general than fine-grained parallel code.
59
(a) (b) (c)
Figure 4.5.: One pulse wave propagating through a circle of arteries. The reflection and
transmission at the branching and converging points are clear.
Figure 4.6.: History profiles at the middle of three arteries of the network with a circle
(see Fig. 4.5(a) for the locations). The reflection and transmission coeffi-
cients meet analytical prediction. See Fig. 4.5 for a series of snapshots of
the wave.
60
4.3.2. A network of human systemic arterial tree
A systemic arterial network with 55 main arteries has been computed in literature.
We adapted the parameters of the arteries from paper [91]. Fig. 4.7 shows a series
of snapshots of pressure distribution map of the systemic network. The top subfigure
in each distribution map shows the heart function and the moment of the snapshot is
indicated by the vertical line. The heart function is a curve of flow rate with a period of
0.8 s (heart rate is 70 beats per minute). The systole stage is a half sinusoidal wave with
a period of 0.4 s. One stroke volume is 120 ml, thus the peak value of flow rate at the root
of aorta is about 471 ml/s. The subfigure in the middle shows the pressure at abdominal
aorta. After the peak, the pressure is damped approximately with an exponential rate,
which meets clinical observations. The systole pressure at the aorta is about 110 mmHg,
and the diastole pressure is about 75 mmHg. This model mimics a health adult. From
series of the snapshots, we can recognize the propagation pattern of a pressure wave in
the systemic arterial network. We notice that the peak pressure at the limbs are higher
than at the aorta as shown by Fig. 4.7(e)4.7(g). This may be explained by the blood
leaks to small side branches which are not considered in this network. This simulation
shows that it is possible to use this code to study more complex scenarios.
61
(a) (b) (c)
62
(a) (b) (c)
Figure 4.8.: One pulse propagating in a kidney of a mouse (mmHg). In order to see the
propagation pattern of the pulse wave in this small network, we make the
period arbitrarily short, i.e. 0.1 s.
63
4.3.4. Speedup of parallelization
The speedup was tested on workstations with 12 cores (Intel(R) Xeon(R) X5650 2.67GHz).
We use the tool time in the Linux environment to count the time run by the code. For
the network of human systemic circulation, we computed 10 cardiac cycles, correspond-
ing to simulating 8 seconds of the network. The time step size was set to 0.01ms, and
therefore there are 0.8 million time steps. The number of mesh points in the ith artery
is set with the formula
Li /ci0
N i = 4b c,
mini=55 i i
i=1 (L /c0 )
where L is the length and c0 is the linearized velocity. We roughly have 4 cells in every
one centimeter of segment. In fact, this is a quite refined discretization in both time
and space. The code runs 20 minutes with one core and 3 minutes with 12 cores. That
is to say, for the classic case of the systemic network with 55 segments, we need 18
seconds for one cardiac cycle. For the case of the mouse kidney, we used a mesh with
over 20 thousand mesh nodes, and we simulated 0.8 seconds with the time step size
set as before. The wall clock time was reduced from 29 minutes with one core to 5
minutes with 12 cores. As shown in Fig. 4.9, the speedup can achieve at about 6 with
12 cores. The difference with the ideal speedup curve may be due to the cache effect
because some conjunctions are shared by two processors. The synchronization overhead
among the processors may be another barrier for the speedup. The serial part of the
code on reading and writing data files takes a few milliseconds, thus their influence
on the speedup is minor. To improve the performance of speedup with more cores,
implementing it with MPI is a possible approach [18].
4.4. Discussion
Simulating blood flow in networks with the 1D model is a very important approach since
it can give distributed information of the blood flow in the network and the computation
cost is relatively low. To do a patient-specific simulation, we need a code which is fast
and can be set easily for networks of various sizes. We developed a code with two
numerical schemes: MacCormack and MUSCL. The MUSCL is shock-capturing and
thus the code can be extended to simulate venous flows. The topology of the network is
accounted by a DAG, which can describe networks with various types of conjunctions.
The computing part of the code is written in C++ with an object oriented approach.
The code is parallelized with OpenMP and then tested with networks originated from
real problems. Compared with a sequential code, the speedup with 12 cores is about
6. In particular, we can compute one cardiac cycle of the classic arterial tree with 55
segments in order of second, which is reported in recent papers to be done in order of
minute [2, 48]. The C++ code reads and writes files in CSV formate. For convenience,
we can use a scripting language (like Python) to set the configurations of the simulations
and make various post-processing of the computing result.
64
Figure 4.9.: The speedup the two computing cases tested on a workstation with 12 cores.
4.5. Conclusion
In this chapter, we developed a parallel code for the simulation of blood flow in networks.
The code was verified and tested on various cases. A good speedup was observed on
workstations with multicores. With this code, one cardiac cycle of the classic arterial
tree can be computed in order of second. It is also very convenient to setup this code to
compute very large networks.
65
Listing 4.1: piece of code parallelized with OpenMP
1 #pragma omp p a r a l l e l
2 {
3 int n u m s t r e a d s=omp get num threads ( ) ;
4 // f o r t h e c o n j u n c t i o n s
5 int i S t a r t c o n , iEnd con ;
6 int chunk con= c e i l ( num conjuncs / f l o a t ( n u m s t r e a d s ) ) ;
7 int myID=omp get thread num ( ) ;
8 i S t a r t c o n=myID chunk con ;
9 i f ( (myID+1) chunk con 1 < num conjuncs 1 )
10 iEnd con =(myID+1) chunk con 1;
11 else
12 iEnd con=num conjuncs 1;
13
14 // f o r t h e a r t e r i e s
15 int i S t a r t a r t , i E n d a r t
16 . . . // v a l u e s s e t as t h e c a s e o f c o n j u n c t i o n s
17
18 while ( currentTime<f i n i s h T i m e ) {
19 f or ( int i=i S t a r t c o n ; i<=iEnd con ; i ++) {
20 i f ( c o n j s [ i ]> i s R o o t ( ) )
21 c o n j s [ i ]>BC inflow ( ) ;
22 e l s e i f ( c o n j s [ i ]> i s L e a f ( ) )
23 c o n j s [ i ]>BC outflow ( ) ;
24 else
25 c o n j s [ i ]>stepMarch ( ) ;
26 }
27 #pragma omp b a r r i e r
28
29 f or ( int i=i S t a r t a r t ; i<=i E n d a r t ; i ++)
30 a r t s [ i ]>stepMarch ( ) ;
31
32 currentTime+=s t e p S i z e ;
33 } // end w h i l e
34 } // pragma omp p a r a l l e l end
66
Part II.
Applications
67
5. Fluid friction and wall viscosity of the 1D
blood flow model: study with an in-vitro
experimental setup
5.1. Introduction
Although the modeling of blood flow has a long history of development, it is still a
very challenging problem until now. Recently the time-domain-based 1D model of blood
flow, which is well balanced between modeling complexity and computational cost, has
attracted more and more attention (see e.g. [3, 55, 72, 86, 99]). It not only predicts quite
well the time-dependent distributions of flow rate and pressure in a network, but also
the mechanical properties of the network may be estimated via analysis of the waves.
This makes the 1D model useful for many applications, such as non-invasive diagnosis,
surgical planning, etc. Moreover, in the context of multi-scale modeling, the 1D model
also plays a very important role along with the 3D and 0D (lumped) models.
The basic time-domain-based 1D model predicts nonlinear waves of flow and pressure.
There are also many damping factors in the system, such as fluid viscosity, wall viscoelas-
ticity, geometrical changes of vessels (e.g. curvature, bifurcation), etc. Previous studies
show that the energy loss at bifurcations is minor [50] and in vessels without drastic
geometrical variations (i.e. no severe aneurysms or stenoses), the fluid viscosity and wall
viscoelasticity seem the most significant damping factors. Comparisons between the 1D
model and in-vivo data [32, 64] suggest that the predictions of a viscoelastic 1D model
is significantly more physiological than those of an elastic one. But the comparisons
were only qualitative or semi-quantitative due to the limited accuracy of non-invasive
measurements and the lack of patient-specific parameter values of the 1D model for each
subject.
Quantitative comparisons can be done with in-vitro experimental setups, which can
be well defined and measured. Reuderink et al. [62] connected a distensible tube to a
piston pump, which ejects fluid in pulse forms into the tube. The experimental mea-
surements were compared with the predictions of several formulations of the 1D model.
In the first formulation, they used an elastic tube law and Poiseuilles theory to ac-
count the fluid viscosity. They found this formulation underestimates the damping of
the waves and predicts shocks, which were not observed in the experiments. In another
formulation, they modeled the fluid viscosity with Womersleys theory and adopted a
viscoelastic tube law in the frequency domain. But to use those models, the nonlinear
The experiments were done in Doshisha University (Kyoto, Japan) with the collaboration of Prof.
Matsukawas group.
68
term (convective term) in the momentum equation of the fluid were neglected. Neverthe-
less, a better match between the predictions of this formulation and the experiments was
found. Bessems et al. [9] built a similar experimental setup and in their simulation they
integrated a 3-component Kelvin viscoelastic model to the 1D fluid model. Compared
with the 1D model with an elastic tube law, the viscoelastic model predicts waveforms
in much better agreement with the experiments. However, to estimate the coefficients
of the 3-component Kelvin model, both the convective and fluid viscosity terms were
neglected. Alastruey et al. [2] did a comparison study on an experimental setup of a
network. They measured the coefficients of a Kelvin-Voigt viscoelastic model with tensile
tests instead of fitting them from the waves. They found that the viscoelasticity damps
significantly the high frequent components of the waves at the peripheral locations of the
network. But for the fluid viscosity term, they adopted a value from literature, which
was fitted from waves of coronary blood flow with elastic model for the wall [76].
In this chapter, we study the friction and wall viscoelasticity with the 1D model and
a similar experimental setup where pulse waves are propagating in one distensible tube.
However, there are three main differences between our study and previous ones:
1. Both of the two damping factors are evaluated. Although there are several theories
to estimate the friction term (see, e.g. [10, 38, 56]), the value is rarely determined
experimentally besides the study of Smith et al. with an elastic model [76]. Since
the fluid viscosity and wall viscoelasticity have similar damping influences on the
pulse waves, it is difficult to evaluate them separately from pulse waves. However,
the viscoelasticity has smoothing effect on the waveforms whereas the fluid friction
does not. We take the advantage of this difference and evaluate both of the two
factors in this study.
69
5.2. Methodology
5.2.1. One dimensional model
Let us remind the distensible vessel with a variable internal radius R(x, t), where x is
the axial distance and t is the time. The 1D model is on the flow rate Q, cross-sectional
area A and internal pressure P . The governing equations describe the conservation of
mass and the balance of momentum respectively,
A Q
+ = 0, (5.1)
t x
Q Q2 A P vx
+ ( ) + = 2 , (5.2)
t x A x r r=R
where vx is the axial velocity, is the fluid density and is the kinematic viscosity of
the fluid. In general, the profile of axial velocity also depends on the radius coordinate
r, v.i.z. vx = vx (r, x, t). If we assume the profile has the same shape (r) in one
segment of vessel, the profile function can be separated as vx = U (x, t)(r), with U the
average velocity. If (r) is known, both and the derivative in the friction term can
be calculated. The friction drag can thus be approximated by Cf U = C pf Q/A. The
profile (r) is strongly dependent on the Womersley number defined by R /, where
the newly introduced quantity is the angular frequency. If and are considered as
constant approximately, only R influences and Cf , whose values should be determined
by experiments for vessels with various diameters. In practice, there are various pairs of
values for and Cf , of which several paris are often seen in literature. In one asymptotic
case when the transient inertial force is very large, the profile is essentially flat, then we
have = 1. If we introduce a very thin viscous boundary layer to match the inviscid
core and a no-slip boundary, the friction term can be estimated (see e.g. [10, 56]). In
another asymptotic case, the fluid viscosity force is dominating and a parabolic profile
appears. In this case, = 4/3 and the friction term is 8Q/A. As fitted by Smith
et al. [76] for coronary blood flow, Cf = 22 and is calculated to be 1.1 according
to the power law profile proposed by Hughes and Lubliner [35]. Some other literature
adopted this value of Cf but assumed = 1 for simplification [2, 48].
To describe the viscoelasticity of the wall, several models have been proposed, see
e.g. [32, 84]. Those models incur various levels of difficulties when they are integrated
with the 1D fluid model and solved numerically. Comparisons between predictions of the
1D model with different viscoelasticity models show that the difference is minor [61, 77].
We choose the two-component Kelvin-Voigt model, which relates the strain and stress
in the equation
d
= E + , (5.3)
dt
where E is the Youngs modulus and is a coefficient for the viscosity.
For a tube with a thin wall, the circumferential strain can be expressed as
R R0
= , (5.4)
(1 2 )R0
70
where R0 is the reference radius without loading and is the Poisson ratio, which is 0.5
for an incompressible material. By Laplaces law, the transmural difference between the
internal pressure P and the external pressure Pext is balanced with the circumferential
stress in the relation
h
P Pext = . (5.5)
R
Combining Eq. 5.3, 5.4 and 5.5, we have
dR
P Pexp = e (R R0 ) + s , (5.6)
dt
with
Eh h
e = , and s = .
(1 2 )A0 (1 2 )A0
Note that the R in the denominators of the two coefficients is approximated by R0 under
the assumption that the perturbations are small.
If we assume Pext is constant and insert Eq. 5.6 into the 1D momentum equation to
remove P , the governing equations can be transformed into
A Q
+ = 0, (5.7)
t x
Q Q2 3 Q 2Q
+ + A 2 = Cf + Cv 2 , (5.8)
t x A 3 A x
where
Eh h
= 2
, and Cv = .
(1 )A0 2(1 2 ) A0
For this formulation of the 1D model, there are detailed discussions on the treatment of
boundary conditions and numerical solvers in [93].
71
h (cm) D (cm) L (cm) LA (cm) LB (cm) (kg/cm3 ) (cm2 /s)
0.2 0.8 192 28.3 168.2 1.050103 1102
Figure 5.1.: Experimental setup. Parameters of the tube and fluid are summarized in
Table 5.1.
Pressure-Wall perturbation
To generate a pulse perturbation, we imposed a sinusoidal wave of one full period for
each time of test, see Fig. 5.2(b). The net volume of fluid injected into the tube was zero,
so the tube would return to the original state finally. Fig. 5.3 shows the time series of
pressure and wall displacement, which were measured simultaneously at one point. We
can see that the wave is reflected several times at each end of the tube and the amplitude
is damped roughly exponentially. After about 6 seconds, the tube restores to the initial
equilibrium state. The energy loss due to the wall viscosity is the work done during the
72
flow rate
flow rate
0.3 0.3
time [s] time [s]
(a) (b)
Figure 5.2.: Output function of the pump for the two tests. Left (a): data fitting vie the
1D model, and right (b): Pressure-Wall perturbation. The actual waveform
may deviate slightly from this ideal one due to experimental uncertainties.
test. Assuming that the test starts at t = t0 , where R = R0 and ends at t = te where
R = Re , we integrate the viscoelastic tube law (5.6) along the path of the moving wall
from t0 to te ,
Z Re Z Re Z Re
dR
P dR = [Pext + e (R R0 )] dR + s dR.
R0 R0 R0 dt
Because the radius at the start and the end of the recorded time interval is the same,
R0 = Re , and the external pressure Pext is constant, the first term on the right hand
side disappears. Thus we obtain
R te dR
P dt
s = R tt0e dRdt .
2
t0 ( dt ) dt
But the relationship between dR/dt and the wall velocity Vw measured by LDV is not
straightforward. In the experiment, the tube was put on a desktop in air and the wall
perturbation was measured from the top (see Fig. 5.1). If we assume the mass center of
each cross section essentially doesnt change due to the inertial force, we have dR/dt =
Vw . However, if the inertial force is neglected, we may assume dR/dt = Vw /2. In reality,
the change rate of the radius should fall between the two values, so we estimated a
range between the two extreme cases. After the viscosity part is subtracted, the elastic
coefficient can then be estimated using linear regression.
Tensile test
We prepared some specimens of the polymer of the vessel wall for a tensile test machine
(Shimadzu EZ test). The specimens were elongated at a rate of 0.5 m/min and then
73
2 0.3
1.5
0.2
1
displacement [mm]
pressure [kPa]
0.5 0.1
0
0.5 0
1
0.1
1.5
2 0.2
0 1 2 3 4 5 6 0 1 2 3 4 5 6
time [s] time [s]
(a) (b)
Figure 5.3.: Simultaneous recording of pressure and displacement of wall at one location:
left (a) pressure, right (b) displacement.
released at the same rate. Fig. 5.4 shows one example where the tension force F is
plotted as a function of elongation L. We applied least square method to fit the curve
against a linear function F = C0 + ESL/L, where C0 is a constant, E is the Youngs
modulus, S is the cross-sectional area of the specimens and L is the original length. The
Youngs modulus was then calculated from the slope of the fitted line. The viscosity was
not estimated with this method.
3.5
experimentaldata
3 elasticmodel
2.5
2
force [N]
1.5
0.5
0.5
0 2 4 6 8 10 12 14
elongation [mm]
Figure 5.4.: Tensile test and linear regression with the linear elastic model.
74
8 exp
0.9E0
7 E0
1.1E0
pressure [kPa] 6
0 1 2 3 4 5
time [s]
Figure 5.5.: Pressure time series at measurement point A. E0 is the best fit for the
Youngs modulus. If E0 is perturbed 10%, the arrival time of each peak
changes significantly. Cf = 22 and = 0.9 kPa s.
5.3. Results
5.3.1. Youngs modulus
We have the best fit between the 1D model and the experiment when E = E0 = 2.08
105 Pa. This value falls in the range estimated with the P-W perturbation method and
is about 8% bigger than the tensile test value. Besides the test error, the variance in the
home-made polymer tubes may also contribute to the difference. Fig. 5.5 compares the
experiment with the 1D model when E has a variance of 10% around E0 . The arrival
time of each peak is significantly later if E decreases and vice versa. Please note that
the final state has a higher pressure than the initial state. That is because this time we
imposed a half sinusoidal wave at the inlet and thus a net volume of about 4.5 cm3 fluid
was injected into the tube.
75
9 exp
Cf = 8
8 Cf = 22
Cf = 33
pressure [kPa] 7
0 1 2 3 4 5
time [s]
Figure 5.6.: Pressure time series at measurement point A. The elastic model predicts
shocks. Increasing the friction term can damp the amplitude effectively, but
the shocks still exist. E = 2.08 105 Pa and = 0.
value predicts the best in amplitude but the shocks still exist.
Table 5.2.: Parameters of fluid friction and wall viscosity and the corresponding NRMS.
where N is the number of data points. We fixed Cf at a series of values ranging from 8
to 33 and fitted by varying its value to minimize the NRMS. Around the estimated
minimum point, ten values of were tested with a step size of 0.1 kPas. As an example,
Fig. 5.7 shows NRMS as a function of when Cf = 22. In this case, 1.0 kPas
is obviously the minimum point. Table 5.2 summarizes the tested values of Cf , the
minimum points of and the corresponding residuals of NRMS. From wave1 to wave7,
the value of Cf increases and that of decreases. The residual achieves a minimum at
wave4 and the worst cases are wave1 and wave7 at the two ends of the table.
76
2.5
2.4
2.3
2.2
NRMS [%]
2.1
1.9
1.8
1.7
0.7 0.8 0.9 1 1.1 1.2 1.3 1.4
[kPa s]
In Fig. 5.8, wave4 is plotted alongside wave1 and wave7. First we notice that the
shocks disappear and the amplitude of the three waves is very close to the experimental
data. However, in the first two seconds, the wave-front of wave7 is more steepened than
others. This difference is more clear in the plot of spectrum (Fig. 5.9), which shows that
the high frequent components of wave7 are underdamped. This is because the damping
effect of wall viscosity is stronger on high frequent waves while that of fluid friction
does not depend on frequency in our model. In the last several seconds, only the main
hormonic is left, thus the difference between the three simulated waves is very small.
With wave4, we can estimate that Cf = 22 and = 1.0 kPa s. The viscoelastic
parameters estimated with various methods are summarized in Table 5.3. The values
estimated by the data fitting with the 1D model fall into the range measured by the
P-W perturbation. Even though the tensile test does not give reliable information on
the wall viscosity, this method estimates a very close Youngs modulus as the 1D model.
77
exp
7 wave1
wave4
6 wave7
pressure [kPa]
0 1 2 3 4 5
time [s]
Figure 5.8.: Pressure time series at measurement point A. E = 2.08 105 Pa. The
values of Cf and for the three simulated waves are shown in Table 5.2.
600
exp
elastic
500 wave1
wave4
400 2.4Hz wave7
|P(f)|
300
200 underdamped
100
0
0 2 4 6 8 10 12 14 16 18 20
frequency
Figure 5.9.: Spectrum of the pressure series (only frequencies less than 20 Hz are shown).
E = 2.08 105 Pa. For the elastic case, Cf = 22 and = 0. The values
of Cf and for the three viscoelastic waves are shown in Table 5.2.
78
Measurement A Measurement B
exp 10 exp
8
elast elast
visco 9 visco
7
8
pressure [kPa]
pressure [kPa]
6
7
5
6
4 5
3 4
2 3
0 1 2 3 4 5 0 1 2 3 4 5
time [s] time [s]
(a) (b)
10
8 exp exp
elast elast
visco 9 visco
7
8
pressure [kPa]
pressure [kPa]
6
7
5
6
4 5
3 4
2 3
0 1 2 3 4 5 0 1 2 3 4 5
time [s] time [s]
(c) (d)
Figure 5.10.: Pressure time series at the two measurement points with two numerical
schemes. Left column: point A, right column: point B. Top row: Mac-
Cormack method, bottom row: MUSCL method. The viscoelastic model
predicts much better than the elastic model at both the measurement
points. The MUSCL method depresses the numerical oscillations when
there are shocks. The parameters are: E = 2.08 105 Pa, Cf = 22,
and =1.0 kPa s (visco).
79
exp
7 0.8Cf0
Cf0
1.2Cf0
pressure [kPa] 6
0 1 2 3 4 5
time [s]
In Fig. 5.10, we plotted the pressure waves at the two measurement points: left column
for point A and right column for point B. The shocks predicted by the elastic model
are very obvious. The MacCormack scheme produces numerical oscillations (top row)
whereas the MUSCL scheme depresses them because of a slope limiter of this scheme
(bottom row). For the viscoelastic model, the shocks disappear and a much better
agreement is achieved at both the locations. If the solution is quite smooth, there is
essentially no difference between the two numerical schemes in accuracy. The good
consistency between the tests at the two locations gives us more confidence on the
experiments and numerical simulations.
We also tested the sensitivity of the model output to the uncertainty of the four
parameters: E, Cf , and . As shown by Fig. 5.5, a perturbation of 10% of E can
change significantly the predicted waves, both on the amplitude and velocity. For Cf
and , an uncertainty of about 20% produces a moderate variance on the predicted wave
(see Fig. 5.11 and 5.12). The sensitivity of the output to Cf and is in the same order.
In contrast, if varies in the range from 1.0 to 1.3, there is no noticeable difference
between the predictions (see Fig. 5.13).
5.4. Discussion
With a well defined and measured experimental setup, we investigated the parameters of
a nonlinear 1D viscoelastic model of blood flow. A small variance of stiffness can change
significantly the mean pressure, pulse pressure and wave velocity. The value of vessel
stiffness can be estimated precisely by the 1D model. It shows again that the 1D model
is a very promising tool of diagnosis of the cardiovascular system.
80
exp
7 0.80
0
1.20
6
pressure [kPa]
0 1 2 3 4 5
time [s]
exp
7 =1
=1.1
=1.3
6
pressure [kPa]
0 1 2 3 4 5
time [s]
81
Previous studies show that wall viscosity damps the waveform considerably (see e.g. [2,
32, 60, 70]). However, its influence on the pulse wave is very difficult to be evaluated
in isolation, since the fluid friction is also a considerable damping factor, which is not
easy to be determined precisely. Previous studies usually make a priori estimations
on the fluid friction or even totally neglect this term. The fluid friction is strongly
dependent on the Womersley number and thus various approximation methods have
been proposed [10, 35, 38, 56]. Fitted with an elastic model for coronary blood flow,
Smith et al. [76] determined the friction coefficient as 22. This value is widely seen
in recent literature (e.g. [2, 48, 95]) along with the value 8, which is derived from a
parabolic (poiseuille) velocity profile (e.g. [78, 102]).
Among several models proposed for the viscoelasticity of the wall, we chose the Kelvin-
Voigt model, which can be readily integrated into the 1D governing equations. Armen-
tano et al. [4] measured the pressure and diameter simultaneously of in-vivo arteries
and fitted the Kelvin-Voigt model to estimate the viscosity coefficient in a range from
3.8 1.3 kPa s to 7.8 1.1 kPa s. Alastruey et al. [2] measured the viscosity of a
polymer by tensile test and found a value of 3.0 0.3 kPa s.
In this study, we solved the nonlinear 1D viscoelastic model with MacCormack and
MUSCL schemes. The elastic model predicts shocks, which are captured by the MUSCL
method without non-physical oscillations. The fluid friction and wall viscosity were fitted
in pairs with the 1D model. We obtained quite good agreement between the 1D model
and experiments. If experimental uncertainties are considered, it can be estimated that
Cf = 22 4 and = 1.0 0.3 kPa s (determined by wave3 and wave5 in Table 5.3).
In our study, the frequency of the main harmonic is 2.4 Hz (see Fig. 5.9) and thus the
Womersley number is about 15.5. The blood flow in the main aorta of an adult human
is also featured with such a value. This confirms that in cases of blood flow with a
similar characteristic Womersley number, the Poiseuille model underestimates the fluid
friction (see e.g. [67]). The widely used value Cf = 22 in large arteries is acceptable.
However, in smaller arteries, the Womersley number can be less than one, so parabolic
velocity profile is more likely to appear, which implies that Cf decreases to 8. Thus
the friction term should vary through the whole cardiovascular system and a smaller
value of Cf should be considered if the Womersley number is much smaller.
The viscoelasticity of the wall damps the high frequent components of the wave, thus
the waveform is not very front-steepened, which has been pointed out by many previous
studies (see e.g. [2, 32]). A perturbation of 20% on wall viscosity introduces moderate
variances on pressure waveform, which is similar with the fluid friction (see Fig. 5.11
and 5.12). The output of the 1D model is not very sensitive to uncertainties of the two
damping factors. Thus it is possible to use general values of those two parameters even
in patient-specific simulations with the 1D model. But on in vivo conditions, the wall
viscosity is much larger as measured by Armentano et al. [4] and the surrounding tissues
of the vessel such as fat may also damp the waves like wall viscosity. Thus neither of
the two viscous factors should be neglected if delicate details of the waves are required.
Moreover, the test on shows that its variance in the range from 1.0 to 1.3 almost does
not influence the waveform (see Fig. 5.13).
82
5.5. Conclusion
We studied the parameters of the nonlinear 1D viscoelastic model with a well defined and
measured experimental setup. The 1D model was solved by two schemes, one of which
is shock-capturing. The fluid friction and wall viscosity were fitted in pairs with data
measured at two locations. The estimated viscoelasticity parameters were in consistent
with values obtained with other methods. The good agreement between the predictions
and the experiments indicates that the nonlinear 1D viscoelastic model can simulate
the pulsatile blood flow very well. We showed that the influence of wall viscosity on
the waveform is in the same order of that of fluid friction. Both of them should not be
neglected if delicate details of waveform are required.
83
6. Effect of viscoelasticity of arterial wall
on waves of blood flow: a study on
network of sheep
6.1. Introduction
The medical value of the pulsatile waveform of the blood flow has attracted humans
attention for a long history. For instance, in the ancient civilization of China [29], there
are many records of diagnosis methods based on the pulse waves. But those methods are
described by qualitative words and there is a lack of experimental observations for the
theories behind the methods. Thus it is hard to make those methods to be accepted by
the major of modern medical community, even though some people believe the methods
have some empirical values.
Starting from the original contribution of Euler in the 18th century, many mathe-
maticians and physicists have improved the math modeling of the pulse waves (for more
details, see e.g. [54, 72]). Simulations of blood flow in networks with the one-dimensional
model have been validated favorably with both in vitro experimental [2, 67, 92] and in
vivo clinical data [20, 56, 63, 64, 78, 81]. Because the 1D model is well balanced between
modeling complexity and computational cost, it has a lot of potential applications, rang-
ing from non-invasive diagnosis, to optimization of surgery planning, to performing as
an important part in multi-scale modeling.
However, challenges still exit for patient-specific simulations. One of them is due to
the complex mechanical properties of the arterial wall. In the 1D model, a constitutive
equation of the tube is necessary to close the system of conservation laws of mass and
momentum. Although the viscoelastic behaviour of the wall has been recognized as a
fundamental factor for a long time, most 1D simulations existing in literature adopted
elastic models for simplicity.
Nevertheless, there are several studies of blood flow in network with viscoelastic 1D
model. The viscoelastic models of arterial wall fall into roughly two classes: Fungs
quasilinear viscoelastic models [28] and spring-dashpot models. The first kind of models
are more generalized. But they are more difficult to handle when incorporated with
the 1D model of blood flow, because they involve a creep function and convolutions
have to be computed [77]. Holenstein et al. [32] proposed a model of the first kind and
fitted the parameters from published data. The viscoelastic model makes considerable
different predictions than the elastic one. Reymond et al. [63, 64] adopted Holenstains
The animal experiments were done in Favaloro University (Buenos Aires, Argentine) by Prof. Armen-
tanos group.
84
model and parameter values in their patient-specific simulations. The comparison be-
tween numerical results and in vivo measurements reveals a considerable impact of the
viscoelasticity on the pulse waves. A comparison between the predictions of the 1D
model with different viscoelastic models of the first kind shows that the difference is mi-
nor [61]. Segers et al. [70] incorporated a frequency dependent viscoelastic model with
the linearized 1D model of blood flow. They found that the influence of viscoelasticity
is comparable with that of the elastic nonlinearity [61]. We note that the parameter
values of the viscoelastic model in those studies are fitted from limited available data in
literature.
The second class of viscoelastic models are in analogy to combinations of springs and
dashpots. One of them is the Kelvin-Voigt model which consists of one spring and one
dashpot connected in parallel. The Kelvin-Voigt model is suitable to describe viscoelastic
solids and it is straightforward to be incorporated with the 1D model of blood flow.
Armentano et al. [4] fitted the Kelvin-Voigt model with the time series of diameter and
pressure which were measured simultaneously. Good agreements between the Kelvin-
Voigt model and measurements were obtained at various locations of the arterial network.
Alastruey et al. [2] also adopted the Kelvin-Voigt model and estimated the parameters
by tensile test. They simulated the pulsatile flow in an in vitro experimental setup and
compared this model with an elastic one. The viscoelastic model agrees much better
with measurements than the elastic one. We note that the vessels in the study were
made of polymers which are actually much less viscous than the real arterial wall.
In this paper, the effect of viscoelasticity on the pulse waves was investigated on
the arterial network of sheep. The simultaneous time series of diameter and pressure at
several arteries were collected from a group of sheep. The viscoelasticity coefficients were
estimated by fitting the measurements against the Kelvin-Voigt model. The pulsatile
blood flow in the network was then computed with a nonlinear 1D viscoelastic model.
We observed the smoothing effect of the wall viscosity on the pulse waveforms.
6.2. Methodology
6.2.1. Data acquisition
The experimental data were obtained from a group of eleven sheep (male Merino, be-
tween 25 and 35 kg). Before surgeries, the animals were anesthetized with sodium
pentobarbital (35 mg/kg). The arterial segments of interest (6 cm long) were separated
from the surrounding tissues. To measure the diameter, two miniature piezoelectric
crystal transducers (5 MHz, 2 mm in diameter) were sutured on opposite sides into the
arterial adventitia. The animals were then sacrificed and the arterial segments of interest
were excised for in vitro test.
The arterial segments were mounted on a test bench where a periodical flow was
generated by an artificial heart (Jarvik Model 5, Kolff Medical Inc., Salt Lake City,
USA). The circulating liquid was an aqueous solution of Tyrode. At each arterial seg-
ment the internal pressure was measured using a solid-state pressure micro-transducer
(Model P2.5, Konigsberg Instruments, Inc., Pasadena, USA), previously calibrated us-
85
ing a mercury manometer at 37 degrees. The arterial diameter signal was calibrated in
millimeters using the 1 mm step calibration option of the sono-micrometer (Model 120,
Triton Technology, San Diego, USA). The transit time of the ultrasonic signal with a ve-
locity of 1,580 m/s was converted to the vessel diameter. The experimental protocol was
conformed to the European Convention for the protection of Vertebrate Animals used
for Experimental and Scientific Purposes. For more details on the animal experiments,
please refer to [85].
The synchronized recording of transmural pressures and diameters was applied on the
following seven anatomical locations as shown in Figure 6.1: Ascending Aorta (AA),
Proximal Descending aorta (PD), Medial Descending aorta (MD), Distal Descending
aorta(DD), Brachiocephalic Trunk (BT), Carotid Artery (CA) and Femoral Artery (FA).
The time series of radius and pressure recorded at AA of one sheep are shown in Fig-
ure 6.2 as an example.
Figure 6.1.: Arterial tree of sheep. Experimental data are collected from eleven sheep
at the following seven locations: Ascending Aorta (AA), Proximal De-
scending aorta (PD), Medial Descending aorta (MD), Distal Descending
aorta (DD), Brachiocephalic Trunk (BT), Carotid Artery (CA) and Femoral
Artery (FA). There are three virtual arteries (VA), which are indicated by
dashed lines, to model the side branches when pulse waves are simulated.
Parameters for all the arteries are shown in Table 6.1.
86
12.4 18
12.2
16
12
pressure [kPa]
radius [mm]
11.8 14
11.6
12
11.4
11.2 10
11
8
0 1 2 3 4 0 1 2 3 4
time [s] time [s]
Figure 6.2.: Time series of radius (left) and transmural pressure (right) recorded at the
ascending aorta for one of the eleven sheep.
87
profile.
p In general, the profile depends on the Womersley number which is defined as
R /, with the angular frequency of the pulse wave and the kinematic viscosity
of the fluid. In practice, Cf usually takes an empirical value fitted from experimental
observations. In this study, we assume Cf = 22 as fitted for the blood flow in large
vessels with a Womersley number of about 10 [76]. Those two PDEs are accompanied
by the constitutive equation (6.4), which can be transformed to
p A
P = Pext + ( A A0 ) + s , (6.7)
t
with
Eh A h
= , s = .
(1 2 )A0 2(1 2 ) A0
Those are the governing equations for the blood flow in one segment.
At the conjunction points between segments, the conservation of mass and momentum
is preserved. The energy loss due to the variation of geometry is neglected. At the inlet
of the network (AA), the flow rate is a cyclic half sinusoidal function in time with a
period of 0.5 s. At the outlets, proper reflection coefficients are imposed. The governing
equations are solved numerically by Monotonic Upstream Scheme for Conservation Laws
(MUSCL) and MacCormack methods, which give essentially the same results. The code
has been favorably validated by analytical results and experimental data, see [92, 93].
6.3. Results
6.3.1. Parameters of the arterial wall
Figures 6.3 and 6.4 show the hysteresis loops and time series at the seven arteries.
The comparison between the measurements and the model shows that the Kelvin-Voigt
model captures quite well the wall viscosity. The modeling error mainly comes from
the nonlinearity of elasticity, which is considerable especially at the central part of the
arterial network. Among the seven arteries, the brachiocephalic trunk has the largest
nonlinearity (Figure 6.4, top row). At the aorta, the nonlinearity decreases from the
proximal part to the distal end (Figure 6.3, top row to bottom row). At the peripheral
arteries, represented by carotid artery and femoral artery, the nonlinearity is negligible
(Figure 6.4, second and bottom rows).
Figure 6.5 presents the statistics of the arterial properties of the group of sheep: the
mean values of reference radius R0 , Youngs modulus E, viscosity coefficient , and the
ratio /E with their corresponding standard deviations. We observe that the variations
of R0 , E and between different sheep are much smaller than between different arteries.
The carotid artery has considerable bigger values of viscoelasticity than other arteries.
Another noticeable pattern is that the ratio /E varies slightly among the different
arteries (Figure 6.6(f)).
88
Ascending Aorta
18 18
exp
KelvinVoigt
16 16
pressure [kPa]
pressure [kPa]
14 14
12 12
10 10
8 8
10 10.5 11 11.5 0 1 2 3
radius [mm] time [s]
14 pressure [kPa] 14
12 12
10 10
8 8
9.4 9.6 9.8 10 10.2 10.4 0 0.5 1 1.5 2 2.5 3
radius [mm] time [s]
16 16
pressure [kPa]
pressure [kPa]
14 14
12 12
10 10
8 8
9.2 9.4 9.6 9.8 10 0 0.5 1 1.5 2 2.5 3
radius [mm] time [s]
14 14
12 12
10 10
8 8
Figure 6.3.: Experimental data and the fitted Kelvin-Voigt model. Left column: hys-
teresis loop. Right column: time series.
Brachiocephalic Trunk
18 exp 18
KelvinVoigt
16 16
pressure [kPa]
pressure [kPa]
14 14
12 12
10 10
8 8
Carotid Artery
exp
18 18
KelvinVoigt
16 16
pressure [kPa]
pressure [kPa]
14 14
12 12
10 10
8 8
3.53 3.54 3.55 3.56 0 0.5 1 1.5 2
radius [mm] time [s]
Femoral Artery
18 exp
18
KelvinVoigt
16 16
pressure [kPa]
pressure [kPa]
14 14
12 12
10 10
Figure 6.4.: Experimental data and the fitted Kelvin-Voigt model (continued). Left col-
umn: hysteresis loop. Right column: time series.
90
0.8 0.5
0.7
0.4
0.6
0.3
E [MPa]
0.5
R0 [cm]
0.4 0.2
0.3
0.1
0.2
0.1 0
AA PD MD DD BT CA FA AA PD MD DD BT CA FA
(a) (b)
12
0.035
10
0.03
8
[kPa s]
0.025
/E [s]
0.02
4
2 0.015
0 0.01
AA PD MD DD BT CA FA AA PD MD DD BT CA FA
(c) (d)
Figure 6.5.: Mean values of the reference radius R0 (top left), Youngs modulus E (top
right), viscosity coefficient (bottom left) and relaxation time /E (bot-
tom right) with standard deviations among the group of sheep at the seven
locations of the arterial network.
Figure 6.6 presents the simulated results at three locations: (i) MD as a representative
for central arteries, and (ii) CA and FA for peripheral arteries. It is obvious that the
viscoelastic model predicts smoother waveforms. The high frequency components of the
91
L R0 h E
Artery (cm) (cm) (cm) (MPa) (kPas)
AA 4 0.7125 0.38 0.1475 2.6156
PD 10 0.8841 0.91 0.0811 1.3866
MD 10 0.6422 1.26 0.0603 1.3187
DD 15 0.6204 1.10 0.1363 2.9425
BT 4 0.6729 1.06 0.0656 1.3427
CA 15 0.3051 0.78 0.4025 8.1286
FA 10 0.2118 0.31 0.2213 4.8186
VA1 20 0.2601 1.00 0.2213 1.0000
VA2 20 0.2762 1.00 0.4025 1.0000
VA3 20 1.1441 1.00 0.2213 1.0000
waveform are damped by the wall viscosity, which is most noticeable in CA.
6.4. Discussion
Several previous studies [32, 61, 64, 70, 77] have shown the significant damping effect
of wall viscosity on the pulse waves. However, the simulations are limited by the lack
of exact values of the model parameters, especially for the viscoelasticity of the arterial
network. We estimate the viscoelasticity by evaluating the pressure-diameter relation-
ship from a dataset of direct measurements on several locations of arterial network of
sheep. With the obtained parameter values, we simulate the pulse waves in the sheep
network with the nonlinear 1D model.
There are several models for the viscoelasticity (see e.g. [9, 32, 70, 84]). Fungs quasi-
linear model is more generalized than the spring-dashpot models, but it is more difficult
to be incorporated with the 1D fluid model, thus it is only applicable to limited formu-
lations of the 1D model (e.g. linearized 1D model [32, 70]). We adopt the Kelvin-Voigt
model to fit the pressure-diameter relationship. The results show that this model cap-
tures the viscosity quite well. Valdez-Jasso et al. [85] tested the Kelvin model, which
adopt a stress relaxation constant as an extra parameter. However, their sensitivity
analysis shows that the model prediction depends least on this constant among all the
parameters. Thus, even though the Kelvin-Voigt doesnt include this constant, the va-
lidity is hardly influenced. Moreover, in contrast to nonlinear optimization methods to
estimate the model parameters in [85], we use the linear regression method which is very
fast and the global optimization is readily guaranteed. We also note that the model
needs further improvements to account for the nonlinearity of elasticity, especially for
the central arteries.
By examining the parameter values among the different arteries, we can clearly see that
smaller arteries tend to be stiffer, which has been pointed out by previous studies [84, 85].
However, we want to stress the ratio between and E in this study. By substituting a
92
Medial Descending aorta
6
40 elast elast
visco 5 visco
35
flow rate [ml/s]
pressure [kPa]
30
4
25
20 3
15
10 2
5
1
3.5 3.75 4 3.5 3.75 4
time [s] time [s]
Carotid Artery
pressure [kPa]
1.6 4.5
1.4 4
1.2 3.5
3
1
2.5
0.8 2
0.6 1.5
0.4 1
3.5 3.75 4 3.5 3.75 4
time [s] time [s]
Femoral Artery
elast 6 elast
1.8
visco visco
1.6
flow rate [ml/s]
5
pressure [kPa]
1.4
1.2 4
1 3
0.8
0.6 2
0.4
1
3.5 3.75 4 3.5 3.75 4
time [s] time [s]
Figure 6.6.: Time series of pressure and flow rate at MD , CA and FA. The viscoelastic
model predicts a smoother waveform than the elastic model.
93
harmonic perturbation (, ) = (0 , 0 )eit into the Kelvin-Voigt equation (6.1), we have
0 = (E + i)0 ,
where G(i) = E + i is the complex modulus. The amplitude of G and the phase
shift are s 2
tr tr
|G| = E 1 + , = arctan ,
tf tf
where tr = /E is the viscoelastic relaxation time and tf = 1/ is the typical forcing
time. If a transit perturbation is imposed on the arterial wall, it will restore to equi-
librium but with a phase lag, which is characterized by the viscoelastic relaxation time.
Figure 6.6(f) shows that the viscoelastic relaxation time seems a biological constant.
Re-writing the phase shift as
= arctan (tr ) ,
we can see that when increases, approaches /2. It is obvious that for bigger values
of , the amplitude of G is bigger. This indicates that the higher the frequency of the
waves, the stronger the damping effect of the viscosity. Since the wavefronts are more
steepened toward the peripheral part of the arterial tree due to the advection effect
of blood flow, the damping effect is more significant in this part. This may be a key
protective factor of the micro-circulations.
In a stiffer vascular network, more pulsation energy of high frequency tend to be
damped in micro-circulations, especially in the brain and kidney [52]. This is considered
to be related with the diseases of these organs, which usually coincide with atheroscle-
rosis. The arterial wall is mainly composed of elastin and muscular fibers, and this
composition varies throughout the whole network, from the aorta to the peripheral ar-
teries. The elastin is more related with the elasticity modulus and the muscular fibers
with the viscoelasticity. The smaller arteries usually have more muscular fibers than
large arteries, and this may also be explained by the need of a stronger damping factor
of pulsations right before the micro-circulations.
6.5. Conclusion
We estimated the viscoelasticity of arterial network of sheep by evaluating the pressure-
diameter relationship with a dataset of direct measurements on several locations of the
network. Quite good agreements between Kelvin-Voigt model and measurements were
achieved through a linear regression method. The obtained parameter values were used in
a nonlinear 1D model to simulate the pulse waves in the arterial network. We have shown
the damping effect of the wall viscosity on the high frequency waves, especially at the
peripheral arteries. This was explained by the nearly constant value of the viscoelastic
relaxation time, which is defined by the ratio between the viscosity coefficient and the
Youngs modulus.
94
7. Blood flow induced by axillofemoral and
femoral-femoral anastomoses with a
severe iliac stenosis: a numerical study
7.1. Introduction
There are some medical situations when the blood flow to the infra aortic arteries needs
to be restored. These situations are caused by severe stenoses or total occlusion of the
aorta or iliac arteries. This causes gluteal, pelvic or leg ischemia, leading to symptoms
from effort pain to tissue loss which may lead to amputation if blood flow is not restored
in time. There are different ways to treat such vascular problems. In case of a simple
artery stenosis or a short artery occlusion, endovascular surgery may be a good option
to improve the blood flow. In other cases, conventional surgery has to be practiced.
The best conventional solution is aorto-femoral or aortobifemoral bypass, well known for
their long term patency, and for restoring complete flow to lower legs. When aorta can
not be the bypass origin, because of major calcifications or because of a fragile patient
condition, another way to restore flow is an axillo-femoral bypass between one axillary
artery and one or two femoral arteries. Surgeons have to realize a proximal anastomosis
between axillary artery and the beginning of prosthetic bypass and a distal anastomosis
between the end of the bypass and the femoral artery. The bypass can be a polyesther
or a polytetrafluoroetylen graft. Clinical data show that bypass patency is correlated
with mechanical properties of the bypass. The mechanisms of this correlation may be
alterations of local or global blood flow. Some studies show that peak flow velocity is
a predictive marker of long term bypass patency. Compliance mismatch contribute to
neointima hyperplasia through wall shear stress modifications at both anastomosis sites.
Another point is that the femoral outflow is correlated to the axillary inflow and to the
hemodynamic bypass conditions. Thus, there are a lot of factors influencing both bypass
patency and femoral outflow.
Aorto-iliac occlusive disease blocks blood supply to the influenced leg. Aorto-femoral
bypass grafting is usually first considered when there is no endovascular option. But
when patients have prohibitive risks to undertake aorto-femoral bypass, Axillo-femoral
bypass grafting (AxFBG) and femoral-femoral bypass grafting (FFBG) become alter-
native options to relieve the ischemia legs [49, 58]. AxFBG is a surgery of tunneling
the axillary artery to the femoral artery by a subcutaneous vascular. Similarly, FFBG
is a bypassing between the healthy femoral artery and the diseased one. In practice,
polyethylene terephthalate (Darcon) and expanded polytetrafluoroethylene (ePTFE) are
used as prosthetic vessels.
95
To simulate the blood flow in a large arterial network, a practical and effective model
is the one-dimensional model. In this model, the hemodynamic quantities are described
as mean values on the cross section of arteries. By solving this model, we can get the
flow rate, the cross-sectional area, the pressure and some other quantities which can
be derived further. In this chapter we apply a vascular 1D simulation to investigate
the hemodynamic changes created by such a bypass, and the clinical relevance of the
hemodynamics.
7.2. Methodology
7.2.1. The 1D model of blood flow
For pulsatile flow in a distensible vessel, the 1D model describes the distributions of the
cross sectional area A, the flow rate Q and the internal pressure P . If we denote the axial
distance of the vessel by x and the time by t, the governing equations for conservation
of mass and balance of momentum write
A Q
+ = 0, (7.1)
t x
Q Q2 A P Q
+ + = Cf , (7.2)
t x A x A
where is the blood density and Cf is a coefficient of fluid friction drag. The value
of Cf is proportional to the kinematic viscosity of the blood and the shear strain
at the interface between the blood and the wall. The latter is strongly related with
the Womersley number. For the blood flow in large arteries, which is featured by a
Womersley number of about 10, Cf has an empirical value of 22.
Those equations are complemented by a constitutive equation of the arterial wall,
which relates A and P
P = P (A).
If the arterial wall is simplified as a thin shell with homogenous elastic material, the
transmural pressure P Pext has a linear relation with the elongation of the internal
radius, p
P Pext = ( A A0 ), (7.3)
where
Eh
= ,
(1 2 )A0
with A0 the reference cross-sectional area, E the Youngs modulus, h the thickness of
the wall and the Poissons ration. For incompressible materials like the arterial wall,
= 1/2. At the conjunction points, the mass and momentum are conserved and the
energy dissipations due to discontinuities are neglected as shown to be reasonable in
literature [50].
96
Figure 7.1.: Illustration of AxFBG (left) and FFBG (right) anastomoses to treat a severe
stenosis at the right iliac artery.
7.2.2. Simulations
We adopted the data in literature [71] for the whole systemic network with 55 main
arteries of a general adult human (see Fig. 7.1). To model a severe stenosis in the
right external iliac artery, we chose a 5-cm segment and set the diameter to be 5% of the
original one, which means that the lumen loss is 97.5%. For the AxFBG, the anastomosis
at the proximal end is made at the right subclavian II. The distance from the root of
the right vertebral to the anastomosis point is 12.2 cm. The distal anastomosis is made
right downstream of the stenosis at the right external iliac artery. The length of the
AxFBG bypass graft is set to 40 cm. For the FFBG approach, the proximal anastomosis
is located at the left external iliac with a distance of 4.4 cm from the root of left internal
iliac. For this surgery, the bypass graft is 20 cm in length. The clinically used prothesis
vessels are 0.4 in radius R and about 0.05 cm in thickness of the wall h. If we adopt these
two values and set the Youngs modus E to 9 MPa, the compliance is 1.78 104 /mmHg,
which is defined as the change ratio of diameter per unit of pressure. The compliance
can be calculated with the formula R/(P R) = R/(Eh). The value of compliance is 1.6
for Darcon and 1.9 for ePTFE [69]. Since the elasticity of the grafts in the simulation is
between the two values, we can consider it as an approximation of the real conditions.
All of the parameters of the two bypass grafts are summarized in Table 7.1.
97
L R0 E h Compliance
grafts (cm) (cm) (MPa) (cm) (104 /mmHg)
A-F 40 0.4 9 0.05 1.79
F-F 20 0.4 9 0.05 1.79
450 normal
400 patho
350 AxFBG
flow rate [ml/s]
300 FFBG
250
200
150
100
50
0
7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9
time [s]
Figure 7.2.: The flow rate at the middle of the ascending aorta.
The governing equations were discretized by the MUSCL (Monotonic Upwind Scheme
for Conservation Law) scheme. The code has been validated by analytical solutions
and experiments [92, 93]. Each cell of the mesh was about 0.15 cm in length. The
time step size was set to 0.04 ms. If we double the number of cells, the change ratio
of numerical results were less than 1.3% in both pressure and flow rate. Thus the
numerical convergence was checked. The cardiac output was modeled by a cyclic half
sinusoidal function with a period of 0.8 s. The total volume of one stroke was 120 ml
and accordingly the peak value of flow rate was 471 ml/s. We plotted in Fig. 7.2 the
flow rate at the middle of the ascending aorta, which is very close to the cardiac output
function we imposed at the inlet of the network. Ten heart beats were run to secure
that a stable state was achieved and only the tenth heart beat was plotted.
7.3. Results
Fig. 7.3 shows the flow rate and pressure at the two bypass grafts. The pressure does not
have big difference at the two grafts (Fig. 7.3(b)). In contrast, the flow rate has a signif-
icant bigger pulse amplitude at the AxFBG graft than at the FFBG graft (Fig. 7.3(a)).
To investigate the blood flow at the host and diseased arteries, as shown in Fig. 7.4
we plotted at the right subclavian, left iliac and right femoral arteries. The flow rate
and pressure are plotted in the four conditions: normal, pathological, and after the two
anastomoses. With the severe stenosis, both the flow rate and pressure drop essentially
to zero at the influenced part of the leg (Fig. 7.4(e) and 7.4(f)). While the stenosis
98
Bypass grafts
35 AxFBG 120
30 FFBG
pressure [mmHg]
flow rate [ml/s]
25 110
20
15 100
10 90
5
0 80
-5
7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9
time [s] time [s]
has negligible effect on the flow rate at the subclavian and left iliac arteries (Fig. 7.4(a)
and 7.4(c)), the pressure increases at those sites (Fig. 7.4(b) and 7.4(d)).
After the AxFBG surgery, the flow rate at the upstream of the anastomosis site in-
creases about 50 ml/s at the peak than the normal condition (red line, Fig. 7.4(a)).
The reduction of blood flow at the downstream of anastomosis sees quite small at the
host arm. The pressure drops significantly, even lower than the normal one (red line,
Fig. 7.4(b)). The influence of the FFBG surgery on the left iliac artery is shown in
Fig. 7.4(c) and 7.4(d). Both the flow rate and pressure at the down stream of the anas-
tomosis site are very close to the normal level. At the upstream, the pressure is also
close to the normal level, but the flow rate increases a lot.
At the diseased leg, after the bypass surgery, both the flow rate and pressure restore to
almost the normal level, see Fig 7.4(e) and 7.4(f). But the waveforms for the two surgeries
are lagged in phase than the normal one, especially the one for the AxFBG. Taken as a
whole, the blood flow restored by the FFBG is better than by the AxFBG. This can be
explained by the quasi-symmetry between the two iliac arteries in the network and the
shorter length of the FFBG graft.
7.4. Discussions
Severe stenosis in iliac arteries blocks blood supply to the whole influenced leg. In
medical practice, there are some conditions when endovascular and aorto-femoral bypass
surgeries cannot be practiced due to many risk factors. Those conditions often happen
to old patients who usually have some other complications. AxFBG and FFBG are the
last options for a vascular surgeon if the arteries in one side of armies or legs of the
patient are healthy. Clinical observations show that the success rate of the anastomosis
surgeries are correlated with the hemodynamics.
We simulated the blood flow of the whole arterial network with a severe stenosis in
the right external iliac artery. Both AxFBG and FFBG can restore blood flow to the
99
R. subclavian
normal 130
80 patho
pressure [mmHg]
AxFBG-up 120
flow rate [ml/s]
60 AxFBG-down 110
40 100
90
20
80
0 70
7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9
time [s] time [s]
(a) (b)
L. iliac
70 normal
60 patho 120
pressure [mmHg]
FFBG-up
flow rate [ml/s]
50 FFBG-down 110
40
100
30
20 90
10 80
0
7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9
time [s] time [s]
(c) (d)
R. femoral
20 140
normal
patho 120
pressure [mmHg]
15 AxFBG
flow rate [ml/s]
100
FFBG
80
10
60
5 40
20
0 0
7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9 7.2 7.3 7.4 7.5 7.6 7.7 7.8 7.9
time [s] time [s]
(e) (f)
Figure 7.4.: Flow rate and pressure at the host and diseased arteries before and after the
anastomoses.
100
diseased leg. However, the FFBG has a better result because the phase lag of the restored
flow is smaller compared with the AxFBG. The proximity of the two iliac arteries gives
a significant advantage to FFBG over AxFBG. A shorter prothesis graft introduces
less resistance and that also seems an important factor to take into account when the
bypassing path is designed.
However, for a patient with severe stenosis in one leg, there is a big possibility that
the other leg also suffers from vascular diseases. This factor has to be considered in
the future research. Moreover the flow velocity and wall shear stress (WSS) are seen
as indicators of long term patency of the prosthesis. But the flow velocity predicted by
the 1D model is the average velocity over the cross section and the velocity measured in
clinics is the centerline velocity. To make them comparable, an exact velocity profile has
to be determined. The estimation of WSS is also related with the velocity profile. The
predictions of those two quantities by the 1D model have to validated by experiments.
7.5. Conclusion
Blood flow in the whole arterial network with a severe iliac stenosis is simulated by the
1D model and two anastomosis surgeries are evaluated. In the ideal case when the other
leg does not suffer vascular diseases, the FFBG surgery is more advantageous over the
AxFBG surgery due to the proximity between the host artery and the receptor artery.
Again, we stress that the 1D model can be a useful tool in evaluations and optimizations
of vascular surgeries.
101
8. Summary and Perspectives
While the study on blood flow has a very long history, it is still a very active research
topic with urgent needs on clinics and a lot of challenging problems. Modeling the blood
flow and properly interpreting the results need knowledge of multidisciplinary fields
including physics, applied mathematics, computer science and medicine. 1D model of
blood flow is well balanced between complexity and accuracy. It can predict the time-
dependent distributions of pressure and flow rate at vascular networks. In this thesis,
we first studied the numerical methods for the 1D model and developed a fast parallel
code which can compute large networks. Moreover, with data collected from in vitro
experimental setup and in vivo arterial networks, we studied the coefficient of fluid
friction and wall viscoelasticity. A preliminary study of the hemodynamics induced by
two anastomoses has been done with the 1D model.
8.1. Summary
This section summarizes the main achievements and conclusions of the thesis. Because
the heart output is cyclic, the blood flow in the vascular system is pulsatile. Due to
the distensibility of the blood vessels, the pulsation propagates in the network with a
characteristic speed of 10 m/s. During the propagation, the waves are reflected and re-
reflected at many points of discontinuities of the network, which makes the waveforms
very complicated. One challenge of computer simulation of the blood flow is due to the
large number of vessel segments. While usually the vascular beds are modeled with 0D
models, this approach can not provide distributive information of the blood flow. Thus
one objective of this thesis is to develop a tool, which can compute very large networks
as fast as possible.
First we recalled the derivation of the 1D model, taking the approach of integrat-
ing along the radius of an axisymmetric flow. The axial velocity profile and tube law
have been given special attention, because they are related with two important damp-
ing parameters of the 1D model. If the viscoelasticity of the vessel wall is described
by a Kelvin-Voigt model, a hyperbolic-parabolic system of governing equations for the
pulsatile blood flow can be derived. Many numerical schemes have been proposed to
discretize the governing equations. We implemented in MATLAB four schemes, namely
MacCormack, Taylor-Galerkin, MUSCL (Monotonic Upstream Scheme for Conservation
Laws) and local discontinuous Galerkin. Those schemes belong to various numerical
frameworks and they were cross compared with tests on several cases. The verifications
of the implementations were done by the linearized model and asymptotic solutions were
computed if there are various source terms. The tests show that the schemes have dif-
ferent features. The MacCormack is very simple to implement and quite accurate (of
102
second order in both time and space) when the nonlinearity of the system is not very
large. The Taylor-Galerkin is well balanced between accuracy and speed in the case
of a moderate nonlinearity. The MUSCL can capture shocks without non-physical os-
cillations, which is a unique feature among the four schemes. The local discontinuous
Galerkin has very small numerical dissipation and dispersion. But if the wall viscosity
is large, a very small time step size is necessitated, which makes the scheme very slow.
To improve the computation speed, a parallel C++ code has been developed. The
topology of the network is described by a directed acyclic graph (DAG). To simulate a
specific network, all the configurations are set by comma-separated-value files (CSV). If
the network is complicated, the parameter files can be prepared by scripts of Python
(other higher level languages are also possible). We computed three cases: a circle of
arteries, a systemic arterial tree with 55 vessels, and an anatomically accurate network
of a mouse kidney with over one thousand segments of vessels. The numerical results
were visualized and the propagation pattern of the pulse waves can be seen clearly. Tests
on multi-core workstations show that good speedup is achieved by parallelizations of the
code for big networks.
The fluid friction and wall viscosity are two damping parameters of the 1D viscoelastic
model of blood flow, but their values are difficult to evaluate separately. For the fluid
friction, previous studies either assume a rough estimation or adopt a value in literature,
which is actually fitted with an elastic model (wall viscosity is neglected). We fitted the
two factors in pairs against pressure waves in a well-defined experimental setup. The
study helps determine the coefficient of fluid friction and the influence of wall viscosity
on the pulse waves.
To evaluate the viscoelasticity of in vivo arteries, we fitted the Kelvin-Voigt model
against the time series of pressure and internal radius measured at several locations of
sheep arterial network. With the estimated parameter values, a network was simulated
with elastic and viscoelastic 1D models. The comparison between the numerical results
shows the damping effect of the viscoelasticity on the high frequency components of
the waveform. A simple theory was proposed to explain the damping mechanism with
a characteristic relaxation time, which seems from the experimental data to be almost
constant throughout the network.
Finally, we applied the developed code to simulate a vascular disease and evaluated the
surgery options from the point of view of hemodynamics. For a patient with severe iliac
stenosis, if there are risk factors preventing intravascular angioplasty and aortofemoral
anastomosis, axillo-femoral and femoral-femoral anastomoses become the only options
of surgeons. The simulation produced results compatible with data in literature. We
evaluated the axillo-femoral and femoral-femoral anastomoses on several aspects: the
rebuilding of blood supply to the diseased leg, the influence on the host arteries, and the
blood flow in the bypass prothesis etc.
103
8.2. Perspectives
The long term objective is to develop a tool with the 1D model which can be used in
clinics for diagnosis and surgery evaluation of vascular diseases. From the work of this
thesis, we can progress in the following directions in the near future.
Coupling the venous system. Compared with the modeling of arterial blood flow,
the venous blood flow attracts less attention. Even though the veins have some
distinct features such as collapsibility, unidirectional valves and considerable ex-
ternal pressures at some locations, the venous flow still can be described by the
1D model with some modifications. In the pulmonary circulation, the pulsations
propagate to the venous system. In the systemic circulation, there are also consid-
erable pressure variations in the venous system due to activities or posture changes.
Coupling the two systems both described by the 1D model can help investigate
the interaction dynamics and further may achieve a closed loop of model for the
whole cardiovascular system.
Estimation of the parameters of the model. One challenge for patient-specific sim-
ulations is to properly set the parameters of the model. Uncertainty qualification
methods such as Monte-Carlo analysis can help determine the parameters whose
uncertainties influence significantly the output of the model. Thus special attention
is needed to improve the measurement techniques of those sensitive parameters.
The data assimilation techniques estimate the model parameters from the model
outputs which are usually easier to measure with various non-invasive techniques
such as MIR, ultrasound imaging etc. Those techniques help make the non-invasive
diagnosis closer to real clinical applications. The fast code developed in this the-
sis can be a convenient tool for numerical experiments of the researches in those
directions.
104
Appendix A.
Simulation and data analysis
A.1. Parameter files and running simulations.
We take a circle of arteries (shown in Fig. A.1) as an example to show how to set the
parameter files. In this network, there are four nodes, labeled from 0 to 3 in red, and
four directed edges, also labeled from 0 to 3. The DAG.csv file looks like the list A.1
Each row is a record of an artery, where the first number stands for the starting node
and the second number the end node. The number of row should be the same as that of
the corresponding artery. By default, the node 0 is the inlet of the network. The inlet
could be specified in flow rate, pressure or cross-sectional area. For instance, to specify
the the flow rate, we shall use the key words <time > and <Q Input> to indicate the
columns of time and flux. Cross-sectional area and pressure are indicated by <A Input>
and <P Input>. The outlets are specified by the number of leaf nodes and the reflection
coefficient. In this example, the outlet is node 3 and the reflection coefficient is 0.0 as
shown in A.3.
105
Listing A.2: Inlet.csv
<time> <Q Input >(cm3/ s )
0 . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 e+00 0 . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 e+00
4 . 0 0 0 0 6 6 6 6 7 7 7 7 7 9 7 e 05 3 . 1 4 1 6 4 4 9 6 2 6 6 0 3 2 1 e 04
8 . 0 0 0 1 3 3 3 3 5 5 5 5 5 9 4 e 05 6 . 2 8 3 2 8 9 6 1 5 2 4 2 3 7 9 e 04
1 . 2 0 0 0 2 0 0 0 0 3 3 3 3 3 9 e 04 9 . 4 2 4 9 3 3 6 4 7 6 6 7 9 4 2 e 04
1 . 6 0 0 0 2 6 6 6 7 1 1 1 1 1 9 e 04 1 . 2 5 6 6 5 7 6 7 4 9 8 5 8 8 4 e 03
2 . 0 0 0 0 3 3 3 3 3 8 8 8 8 9 8 e 04 1 . 5 7 0 8 2 1 8 6 1 1 7 3 6 9 9 e 03
...
Similarly, the file for the vessel properties looks like Fig. A.2. The key works are
self-explanatory.
The simulation and recording time are set in the list A.4
In list A.5, we specify the recording location: the first column is the artery and the
second column is the mesh point.
Now all the parameter files are ready and we are going to run the simulations in a
unix-like environment. Before running the simulation, we can set the number of threads
explicitly in the terminal by the command:
export OMP NUM THREADS=4
If we put all those files in a subfolder named Param, and create a subdirectory to put
the result, for instance data, we can initiate the computing by the command:
bloodflow -i param -o data -v
The option -v (verbose) is for debugging of the parameters file. With this option, the
configurations of the simulation will be shown like Fig. A.3.
106
Listing A.4: time setup.csv
< t f i n a l >( s ) , < t s t e p >( s ) , < r e c o r d s t a r t >,< r e c o r d e n d >
2.4 ,4 e 05 ,0 ,2.4
defined methods, all the parameters can be written in a compatible form with the C++
code. List A.6 illustrates the main steps to do this. The file fig data.csv stores processes
data with Python scripts. The first two columns are the coordinates of the 2D network,
and the following columns are the pressure at corresponding points. We give in List A.7
a bash script which makes an animation with Gnuplot as the plotting engine.
107
Figure A.3.: Checking the configurations of the simulation.
a r t s =[ ] # a l i s t o f a r t e r i e s
. . . #g e n e r a t e each a r t e r y , s e t t h e paramters ,
. . . # and append i t t o t h e l i s t a r t s .
c o n j s =[ ] # a l i s t o f c o n j u n c t i o n s
. . . #g e n e r a t e each c o n j u n c t i o n , s e t t h e p a r a m t e r s
. . . # and t h e n append i t t o t h e l i s t c o n j s .
c o n j s t S=timeSetup ( )
. . . #s e t time c o n t r o l parameters , i n f l o w B.C. , e t c .
# g e n e r a t e a network
network=network (ARTS=a r t s ,CONS=c o n j s )
network . tS=tS
#w r i t e a l l t h e paramter f i l e s t o t h e d e s t i n a t i o n d i r e c t o r y .
network . writeParam ( d e s t i n a t i o n d i r e c t o r y )
108
Listing A.7: make animation.sh
#! / b i n / b a s h
function p l o t 1 s n a p s h o t ( ) {
Num snap=$1
p r i n t f v name %02 i $1
g n u p l o t << EOF
set t e r m i n a l g i f
d a t a F i l e= f i g d a t a . c s v
set s t y l e data l i n e s
unset x t i c s
unset y t i c s
unset z t i c s
unset b o r d e r
set o d e s t i n a t i o n d i r / anim $ {name } . g i f
unset key
splot [ ] [ ] [ 6 : 1 6 ] dataFile u 1:2:3 w l l s 1 , \
u 1 : 2 : $ ( ( $ {Num snap }+2))
EOF
}
for i in { 1 . . 4 0 }
do
plot 1snapshot $i
done
c o n v e r t d e l a y 15 l o o p 0 d e s t i n a t i o n d i r / . g i f a n i m a t i o n . g i f
109
A.3. Method of fitting the Kelvin-Voigt model with the data
measured on sheep.
Given the Kelvin-Voigt model for the arterial wall,
Eh Eh 1 h dR
P = 2
2
+ 2
, (A.1)
(1 )R0 (1 ) R (1 )R0 Rdt
with P the transmural pressure, R the internal radius, E the Young modulus, h the
thickness of the arterial wall and the Poisson constant. The time series of P and R
are measured in experiments, and the parameters of the model are wanted. Written in
matrix form, the problem can be expressed as
P (t) = AC,
where A is a N 3 matrix [(1, ..., 1)T , 1/R, dR/(Rdt)], with N the number of data points,
C the 3 1 coefficient vector [Eh/((1 2 )R0 ), Eh/(1 2 ), h/((1 2 )R0 )]T . And
the objective cost function is
v
uN
1u X
J(C) = t ((Pmodel )i Pi )2 ),
N
i
with Pmodel the pressure predicted by the model. We assume that the columns of the
data matrix are independent in the linear space, the errors of the measurement data are
independent and identically distributed. According to the theory of least square method,
C is given by,
C = (AT A)1 AT P.
In the data matrix, the derivative of R is evaluated by a spectral numerical method.
Given a times series R(t) with period T , it can be expanded in Fourier series
X 2i
kt
R= R(k)e T ,
k=
where R is given by Z T
1 2i
R(k) = R(t)e T
kt
dt.
T t=0
For the derivative, one has
dR X 2i 2i kt
= R ke T .
dt T
k=
In the computation,we take advantage of the Discrete Fourier Transform (DFT). The
measurement noises are filtered out through a loop in the calculation. The pseudocode
is listed below:
110
Step 1: Evaluate the DFT of R (assume N as an even number without loss of
generality)
N
2
1 X 2i
R = Re N
nk
.
N
n= N
2
+1
Step 2: If ||Rn || < , set Rn = 0, where is a given criterion whose value should
be optimized through the loop.
2ni
Step 3: Multiply Rn by a factor T to get DR
d
N
2
dR X
d 2i nk
= DRe N
dt
n= N
2
+1
Step 5: Solve the least square problem and evaluate the objective function.
Step 6: Change and return back to Step 2 until the value of objective function
is not decreasing.
111
Bibliography
[1] Global status report on noncommunicable diseases 2010. Technical report, World
Health Organization, 2011.
[2] J. Alastruey, A.W. Khir, K.S. Matthys, P. Segers, S.J. Sherwin, P.R. Verdonck,
K.H. Parker, and J. Peiro. Pulse wave propagation in a model human arterial net-
work: Assessment of 1-d visco-elastic simulations against in vitro measurements.
Journal of biomechanics, 44(12):22502258, 2011.
[3] J. Alastruey, K.H. Parker, J. Peiro, S.M. Byrd, and S.J. Sherwin. Modelling
the circle of willis to assess the effects of anatomical variations and occlusions on
cerebral flows. Journal of Biomechanics, 40(8):17941805, 2007.
[4] R.L. Armentano, J.G. Barra, J. Levenson, A. Simon, and R.H. Pichel. Arterial
wall mechanics in conscious dogs assessment of viscous, inertial, and elastic moduli
to characterize aortic wall behavior. Circulation Research, 76(3):468478, 1995.
[5] K. Azer and C.S. Peskin. A one-dimensional model of blood flow in arteries with
friction and convection based on the womersley velocity profile. Cardiovascular
Engineering, 7(2):5173, 2007.
[6] R.E. Ballard, D.E. Watenpaugh, G.A. Breit, G. Murthy, D.C. Holley, and A.R.
Hargens. Leg intramuscular pressures during locomotion in humans. Journal of
Applied Physiology, 84(6):19761981, 1998.
[7] A.C.L. Barnard, W.A. Hunt, W.P. Timlake, and E. Varley. A theory of fluid flow
in compliant tubes. Biophysical Journal, 6(6):717724, 1966.
[8] C. Bertoglio, P. Moireau, and J.-F. Gerbeau. Sequential parameter estimation for
fluidstructure problems: Application to hemodynamics. International Journal
for Numerical Methods in Biomedical Engineering, 28(4):434455, 2012.
[9] D. Bessems, C.G. Giannopapa, M. Rutten, and F.N. van de Vosse. Experimen-
tal validation of a time-domain-based wave propagation model of blood flow in
viscoelastic vessels. Journal of biomechanics, 41(2):284291, 2008.
[10] D. Bessems, M. Rutten, and F. Van De Vosse. A wave propagation model of blood
flow in large vessels using an approximate velocity profile function. Journal of
Fluid Mechanics, 580:145168, 2007.
112
[11] F. Bouchut. Nonlinear Stability of Finite Volume Methods for Hyperbolic Conser-
vation Laws: And Well-Balanced Schemes for Sources. Springer, Basel switzerland,
2004.
[12] B.S. Brook, S.A.E.G. Falle, and T.J. Pedley. Numerical solutions for unsteady
gravity-driven flows in collapsible tubes: evolution and roll-wave instability of a
steady state. Journal of Fluid Mechanics, 396:223256, 1999.
[13] B.S. Brook and T.J. Pedley. A model for time-dependent flow in (giraffe jugular)
veins: uniform tube properties. Journal of biomechanics, 35(1):95107, 2002.
[15] N. Cavallini, V. Caleffi, and V. Coscia. Finite volume and weno scheme in one-
dimensional vascular system modelling. Computers and Mathematics with Appli-
cations, 56(9):23822397, 2008.
[17] Y. Chen, L. Zhang, D. Zhang, and D. Zhang. Wrist pulse signal diagnosis using
modified gaussian models and fuzzy c-means classification. Medical engineering &
physics, 31(10):12831289, 2009.
[18] H. Coullon, J.-M. Fullana, P.-Y. Lagree, S. Limet, and X. Wang. Blood flow
arterial network simulation with the implicit parallelism library skelgis. In ICCS
2014, 2014.
[19] O. Delestre and P.-Y. Lagree. A well balanced finite volume scheme for blood
flow simulations. International Journal for Numerical Methods in Fluids, page
doi: 10.1002/fld.3736, 2012.
[20] K. DeVault, P.A Gremaud, V. Novak, M.S. Olufsen, G. Vernieres, and P. Zhao.
Blood flow in the circle of willis: Modeling and calibration. Multiscale Modeling
& Simulation, 7(2):888909, 2008.
[22] D. Elad, D. Katz, E. Kimmel, and S. Einav. Numerical schemes for unsteady fluid
flow through collapsible tubes. Journal of biomedical engineering, 13(1):1018,
1991.
113
[23] M. Fernandez, V. Milisic, and A. Quarteroni. Analysis of a geometrical multiscale
blood flow model based on the coupling of odes and hyperbolic pdes. Multiscale
Modeling & Simulation, 4(1):215236, 2005.
[24] P. Flaud, P. Guesdon, and J.-M. Fullana. Experiments of draining and filling
processes in a collapsible tube at high external pressure. The European Physical
Journal Applied Physics, 57(03):31101p131101p13, 2012.
[27] J.-M. Fullana and S. Zaleski. A branched one-dimensional model of vessel networks.
Journal of Fluid Mechanics, 621(1):183204, 2009.
[29] Y.C. Fung. Biomechanics: circulation. Springer Verlag, New York, US, 1997.
[30] J.-F. Gerbeau, M. Vidrascu, and P. Frey. Fluidstructure interaction in blood flows
on geometries based on medical imaging. Computers and Structures, 83(2):155
165, 2005.
[31] J.S. Hesthaven and T. Warburton. Nodal discontinuous Galerkin methods: algo-
rithms, analysis, and applications, volume 54. Springer-Verlag New York Inc, New
York, US, 2008.
[32] R. Holenstein, P. Niederer, and M. Anliker. A viscoelastic model for use in pre-
dicting arterial pulse waves. Journal of biomechanical engineering, 102(4):318325,
1980.
[33] W. Huberts, A.S. Bode, W. Kroon, R.N. Planken, J.H.M. Tordoir, F.N. Van de
Vosse, and E.M.H. Bosboom. A pulse wave propagation model to support decision-
making in vascular access planning in the clinic. Medical engineering & physics,
34(2):233248, 2012.
[35] T.J.R. Hughes and J. Lubliner. On the one-dimensional theory of blood flow in
the larger vessels. Mathematical Biosciences, 18(1):161170, 1973.
114
[36] Y. Kivity and R. Collins. Nonlinear wave propagation in viscoelastic tubes: appli-
cation to aortic rupture. Journal of Biomechanics, 7(1):6776, 1974.
[37] D. Kuzmin. Slope limiting for discontinuous galerkin approximations with a pos-
sibly non-orthogonal taylor basis. International Journal for Numerical Methods in
Fluids, 71(9):11781190, 2013.
[38] P.-Y. Lagree. An inverse technique to deduce the elasticity of a large artery. EPJ
Applied Physics, 9(2):153164, 2000.
[39] C.A.D. Leguy, E.M.H. Bosboom, A.S.Z. Belloum, A.P.G. Hoeks, and F.N. van de
Vosse. Global sensitivity analysis of a wave propagation model for arm arteries.
Medical engineering & physics, 33(8):10081016, 2011.
[40] C.A.D. Leguy, E.M.H. Bosboom, H. Gelderblom, A.P.G. Hoeks, and F.N. van de
Vosse. Estimation of distributed arterial mechanical properties using a wave prop-
agation model in a reverse way. Medical engineering & physics, 32(9):957967,
2010.
[41] R.J. LeVeque. Finite volume methods for hyperbolic problems, volume 31. Cam-
bridge University Press, Cambridge UK, 2002.
[43] S.C. Ling and H.B. Atabek. A nonlinear analysis of pulsatile flow in arteries.
Journal of Fluid Mechanics, 55(03):493511, 1972.
[46] A.C.I. Malossi, P.J. Blanco, and S. Deparis. A two-level time step technique for
the partitioned solution of one-dimensional arterial networks. Computer Methods
in Applied Mechanics and Engineering, 237:212226, 2012.
[49] D. Martin and S.G. Katz. Axillofemoral bypass for aortoiliac occlusive disease.
The American journal of surgery, 180(2):100103, 2000.
115
[50] K.S. Matthys, J. Alastruey, J. Peiro, A.W. Khir, P. Segers, P.R. Verdonck, K.H.
Parker, and S.J. Sherwin. Pulse wave propagation in a model human arterial
network: Assessment of 1-d numerical simulations against in vitro measurements.
Journal of Biomechanics, 40(15):34763486, 2007.
[51] J.P. Mynard and P. Nithiarasu. A 1d arterial blood flow model incorporating ven-
tricular pressure, aortic valve and regional coronary flow using the locally conserva-
tive galerkin (lcg) method. Communications in Numerical Methods in Engineering,
24(5):367417, 2008.
[53] F. Nicoud, H. Vernhet, and M. Dauzat. A numerical assessment of wall shear stress
changes after endovascular stenting. Journal of Biomechanics, 38(10):20192027,
2005.
[54] M.S. Olufsen. Modeling the arterial system with reference to an anestesia simula-
tor. PhD thesis, Roskilde Universitetscenter,, 1998.
[55] M.S. Olufsen, N.A. Hill, G.D.A. Vaughan, C. Sainsbury, and M. Johnson. Rar-
efaction and blood pressure in systemic and pulmonary arteries. Journal of fluid
mechanics, 705:280305, 2012.
[56] M.S. Olufsen, C.S. Peskin, W.Y. Kim, E.M. Pedersen, A. Nadim, and J. Larsen.
Numerical simulation and experimental validation of blood flow in arteries
with structured-tree outflow conditions. Annals of Biomedical Engineering,
28(11):12811299, 2000.
[57] G.L. Papageorgiou, B.N. Jones, V.J. Redding, and N. Hudson. The area ratio of
normal arterial junctions and its implications in pulse wave reflections. Cardiovas-
cular research, 24(6):478484, 1990.
[58] M.A. Passman, L.M. Taylor Jr, G.L. Moneta, J.M. Edwards, R.A. Yeager, D.B.
McConnell, and J.M. Porter. Comparison of axillofemoral and aortofemoral bypass
for aortoiliac occlusive disease. Journal of vascular surgery, 23(2):263271, 1996.
[59] T.J. Pedley. The Fluid Mechanics of Large Blood Vessels. Cambridge University
Press, Cambridge UK, 1980.
[60] R. Raghu and C.A. Taylor. Verification of a one-dimensional finite element method
for modeling blood flow in the cardiovascular system incorporating a viscoelastic
wall model. Finite Elements in Analysis and Design, 47(6):586592, 2011.
[61] R. Raghu, I.E. Vignon-Clementel, C.A. Figueroa, and C.A. Taylor. Compara-
tive study of viscoelastic arterial wall models in nonlinear one-dimensional fi-
nite element simulations of blood flow. Journal of biomechanical engineering,
133(8):081003, 2011.
116
[62] P.J. Reuderink, H.W. Hoogstraten, P. Sipkema, B. Hillen, and N. Westerhof. Lin-
ear and nonlinear one-dimensional models of pulse wave transmission at high wom-
ersley numbers. Journal of biomechanics, 22(8):819827, 1989.
[65] M. Saito. One-dimensional modeling of pulse wave for a human artery model.
Technical report, Universite Pierre et Marie-Curie, 2010.
[68] S. Sankaran, M.E. Moghadam, A.M. Kahn, E.E. Tseng, J.M. Guccione, and A.L.
Marsden. Patient-specific multiscale modeling of blood flow for coronary artery
bypass graft surgery. Annals of biomedical engineering, 40(10):22282242, 2012.
[69] S. Sarkar, H.J. Salacinski, G. Hamilton, and A.M. Seifalian. The mechanical
properties of infrainguinal vascular bypass grafts: their role in influencing patency.
European journal of vascular and endovascular surgery, 31(6):627636, 2006.
[73] Y. Shi, P. Lawford, and R. Hose. Review of zero-d and 1-d models of blood flow
in the cardiovascular system. Biomed. Eng. Online, 10(1):33, 2011.
117
[74] C.W. Shu. Different formulations of the discontinuous galerkin method for the
viscous terms. Advances in Scientific Computing, pages 144155, 2001.
[76] N.P. Smith, A.J. Pullan, and P.J. Hunter. An anatomically based model of tran-
sient coronary blood flow in the heart. SIAM Journal on Applied mathematics,
62(3):9901018, 2002.
[77] B.N. Steele, D. Valdez-Jasso, M.A. Haider, and M.S. Olufsen. Predicting arterial
flow and pressure dynamics using a 1d fluid dynamics model with a viscoelastic
wall. SIAM Journal on Applied Mathematics, 71(4):11231143, 2011.
[78] B.N. Steele, J. Wan, J.P. Ku, T.J.R. Hughes, and C.A. Taylor. In vivo validation of
a one-dimensional finite-element method for predicting blood flow in cardiovascular
bypass grafts. IEEE Transactions on Biomedical Engineering, 50(6):649656, 2003.
[79] P.D. Stein and H.N. Sabbah. Turbulent blood flow in the ascending aorta of
humans with normal and diseased aortic valves. Circulation research, 39(1):5865,
1976.
[80] N. Stergiopulos, D.F. Young, and T.R. Rogge. Computer simulation of arterial
flow with applications to arterial and aortic stenoses. Journal of Biomechanics,
25(12):14771488, 1992.
[81] J.C. Stettler, P. Niederer, and M. Anliker. Theoretical analysis of arterial hemody-
namics including the influence of bifurcations. Annals of biomedical engineering,
9(2):145164, 1981.
[82] A.C.Y. Tang, J.W.Y. Chung, and T.K.S. Wong. Validation of a novel tradi-
tional chinese medicine pulse diagnostic model using an artificial neural network.
Evidence-Based Complementary and Alternative Medicine, 2012, 2011.
[83] E.F. Toro. Riemann solvers and numerical methods for fluid dynamics, volume 16.
Springer, 1999.
[84] D. Valdez-Jasso, D. Bia, Y. Zocalo, R.L. Armentano, M.A. Haider, and M.S.
Olufsen. Linear and nonlinear viscoelastic modeling of aorta and carotid pressure
area dynamics under in vivo and ex vivo conditions. Annals of Biomedical Engi-
neering, 39(5):14381456, 2011.
[85] D. Valdez-Jasso, M.A. Haider, H.T. Banks, D.B. Santana, Y.Z. German, R.L.
Armentano, and M.S. Olufsen. Analysis of viscoelastic wall properties in ovine
arteries. IEEE Transactions on Biomedical Engineering, 56(2):210219, 2009.
[86] F.N. van de Vosse and N. Stergiopulos. Pulse wave propagation in the arterial
tree. Annual Review of Fluid Mechanics, 43(1), 2011.
118
[87] I.E. Vignon and C.A. Taylor. Outflow boundary conditions for one-dimensional
finite element modeling of blood flow and pressure waves in arteries. Wave Motion,
39(4):361374, 2004.
[88] I.E. Vignon-Clementel, C.A. Figueroa, K.E. Jansen, and C.A. Taylor. Outflow
boundary conditions for 3d simulations of non-periodic blood flow and pressure
fields in deformable arteries. Computer methods in biomechanics and biomedical
engineering, 13(5):625640, 2010.
[89] J. Wan, B. Steele, S.A. Spicer, S. Strohband, G.R. Feijo, T.J.R. Hughes, and
C.A. Taylor. A one-dimensional finite element method for simulation-based med-
ical planning for cardiovascular disease. Computer Methods in Biomechanics and
Biomedical Engineering, 5(3):195206, 2002.
[90] H. Wang and Y. Cheng. A quantitative system for pulse diagnosis in traditional
chinese medicine. In Engineering in Medicine and Biology Society, 2005. IEEE-
EMBS 2005. 27th Annual International Conference of the, pages 56765679. IEEE,
2006.
[91] J.J. Wang and K.H. Parker. Wave propagation in a model of the arterial circula-
tion. Journal of biomechanics, 37(4):457470, 2004.
[92] X. Wang, O. Delestre, J.-M. Fullana, M. Saito, Y. Ikenaga, M. Matsukawa, and
P.-Y. Lagree. Comparing different numerical methods for solving arterial 1d flows
in networks. Computer Methods in Biomechanics and Biomedical Engineering,
15(sup1):6162, 2012.
[93] X. Wang, J.-M. Fullana, and P.-Y. Lagree. Verification and comparison of four
numerical schemes for a 1d viscoelastic blood flow model. Computer Methods in
Biomechanics and Biomedical Engineering, In press, 2014.
[94] M. Wibmer. One-dimensional simulation of arterial blood flow with applications.
PhD thesis, Vienna University of Technology, Vienna, Austria, 2004.
[95] M. Willemet, V. Lacroix, and E. Marchandise. Validation of a 1d patient-specific
model of the arterial hemodynamics in bypassed lower-limbs: Simulations against
in vivo measurements. Medical engineering & physics, 2013.
[96] N. Xiao, J. Alastruey, and C. Alberto Figueroa. A systematic comparison between
1-d and 3-d hemodynamics in compliant arterial models. International journal for
numerical methods in biomedical engineering, 30(2):204231, 2014.
[97] N. Xiao, J.D. Humphrey, and C.A. Figueroa. Multi-scale computational model of
three-dimensional hemodynamics within a deformable full-body arterial network.
Journal of computational physics, 244:2240, 2013.
[98] D. Xiu and S.J. Sherwin. Parametric uncertainty analysis of pulse wave propaga-
tion in a model of a human arterial network. Journal of Computational Physics,
226(2):13851407, 2007.
119
[99] Y. Yamamoto, M. Saito, Y. Ikenaga, M. Matsukawa, Y. Watanabe, M. Furuya,
and T. Asada. Experimental study on the pulse wave propagation in a human
artery model. Japanese Journal of Applied Physics, 50(7), 2011.
[100] S. Yomosa. Solitary waves in large blood vessels. Journal of the Physical Society
of Japan, 56:506520, 1987.
[102] M. Zagzoule and J.-P. Marc-Vergnes. A global mathematical model of the cerebral
circulation in man. Journal of Biomechanics, 19(12):10151022, 1986.
120