- Email: [email protected]
Numerical 3D-simulation of pulsatile wall shear stress in an arterial T-bifurcation model
Numerical 3D-simulation of pulsatile wall shear stress in an arterial T-bifurcation model K. Perktold and R. Peter Institute
of Mathematics,
Received
April
1989,
Technical
accepted
July
University
Graz, Steyrergasse
17,A-8010 Graz, Austria
1989
ABSTRACT The structure of pulsatile bloodjow and wall shear stress in a 90” T-bifurcation model is analyxd numerically. The nonlinear Navier-Stokes equations for time-dependent incompressible .hewtonian JluidJow are approximated using a newly developed pressure correction, Jinite element method. The wall shear stress is calculated from the Jinite element velocifv jield. The investigation shows viscousJow phenomena such asjow separation and stagnation and the distribution of high and low wall shear stress during the pulse cycle. Furthermore, the effect of a sharp corner at the bifurcation edge on the wall shear stress is anaivsed. Detailed 1ocalJow investigation is required to examine&id &namic contribution to the development of arterial diseases such as atherosclerosis and thrombosis. Keywords:
Arterial
bifurcation,
pulsatile
blood
flow, numerical
INTRODUCTION Many biomechanical studies concentrate on local flow behaviour and wall shear stress in blood vessels. The goal of such investigations is the determination of the haemodynamic factors in the development of vascular diseases. Clinical observations show that pathological vessel alterations occur preferentially in zones where the flow is locally disturbed and extreme wall shear stresses appear. The detailed analysis of flow phenomena and comparison with the pattern of diseases should contribute to the comprehension of the haemodynamic influence on pathoThe effects of interest mainly logical alterations. concern the occurrence of flow recirculation and flow stagnation zones as well as areas of high velocity gradient where the endothelium can be damaged and platelets can be stimulated. The recirculation zones where stasis results are treated in the ‘lowshear’ theory of Caro et al.‘, while the aspects of elevated shear stresses are expressed in the ‘highshear’ theory of Fry2. The main interest in local flow dynamics has been concentrated on vessel bifurcations and vessel bends, where diseases often occur. Several recent papers are devoted to experimental studies of the flow structure in 90” T-bifurcation models under different assumpthe results of lasertions. In these papers, or hot-film measurements” Doppler measurements” and of experimental flow visualizations5,6 are presented. We examine here the flow structure and wall shear stress under physiological pulsatile flow conditions using computer simulation. The mathematical approach yields detailed information on the local pulsatile arterial flow. Efficient numerical methods to approximate the viscous flow equations and modern computer technology permit the development Correspondcncc
and rrprint
requests
to: Karl Prrktold
f,1990 Butterworth & Co (Publishers) Ltd 0141-5425/90/010002-II $03.00 2 J. Biomcd. Eng. 1990. Vol. 12, January
investigation,
flow pattern,
wall shear
stress
of physiologically realistic models. Amongst the numerical procedures which can be applied, the finite element method has the most advantages and because of its flexibility it is well suited to the analysis of complicated arterial flow. The appropriate application of the finite element method on viscous flow equations can produce accurate numerical results for fluid velocity, an essential requirement for the calculation of realistic wall shear stresses during the pulse cycle. In this study, numerical simulation is performed using a three-dimensional 90” T-junction model. The pulsatile flow pattern and the wall shear stress distribution are analysed corresponding to a physiological velocity pulse waveform occurring in the left coronary artery. The calculation is carried out for a mean inflow Reynolds number Re = 80, whereby the radius of the parent vessel is used in the definition. The flow in the branch is assumed to be 300/, of the entrance flow. Numerical results are presented for non-dimensional axial and secondary flow velocity as well as for wall shear stress in physical quantities.
MATHEMATICAL MODEL Computer simulation of local blood flow phenomena in arteries requires the full Navier4tokes equations; simplifications of the non-linear flow equations cannot be successful because of the occurrence of flow separation and secondary motion. In this study blood is assumed to be an incompressible, homogeneous, Newtonian fluid ; its specific rheological properties are not taken into account. This idealization can be made for the calculation of blood flow in large arteries where the shear rate is relatively high. The vessel wall is assumed to be rigid; an assumption that is admissible for the determination of local flow
phenomena in short artery segments where pulse wave propagation is not treated. In the primitive variables flow velocity and pressure, the pulsatile, incompressible, Newtonian, 3DNavier Stokes equations can be written as
presses the velocity components and the pressure in terms ofinterpolation functions over each element of the subdivided flow domain. Here 3D brick elements are used. The appropriate solution for velocity and pressure in an element (e) takes the form
1)
v.u=o
(2)
the normalized where u = (u, ZJ,x!) 7 and p denote non-dimensional flow velocity and pressure, respectively. The velocity components U, G, UJ are related to the Gartesian coordinates x,y, Z. The normalization is done with respect to a characteristic length I,, a characteristic velocity U,, and the kinematic viscosity I’. Re and q denote the time-mean Reynolds number and the Stokes number, respectively. and are defined as R e = <:li,!
(3a)
V
(3b) The non-dimensional and pressure are
units
length,
time,
velocity
4) physical units. The other mathewhere - indicates matical symbols have standard meanings: w is angular frequency and e is constant fluid density. For the flow problems which are treated here the boundary of the flow domain partitions into two sections with different types of boundary conditions. At the inflow boundary and at the vessel wall the velocity is given where the flow at the entrance is prescribed with time-dependent velocity profiles, and at the rigid wall the no-slip condition is valid %,,,. = s(l
and
UWll = 0
(5)
The outflow is assumed traction-free. The condition of vanishing surface traction can be described mathematically as (6) where II is the normal unit vector at the outflow boundary. The choice of appropriate boundary conditions at the flow entrance and exit often produces problems. It is required that the assumed conditions do not disturb the numerical results in the flow domain of interest. Numerical
method
The Navier-Stokes equations are spatially discretized using the Galerkin finite element method where the differential equations are reformulated with an equivalent variational formulation in appropriate The finite element method exfunction spaces7.
(7) P’(4 =
MP;,
(4
(8)
where .,V,, i = 1,. . ,8 and ;21 are the element interpolation functions for the velocity and the pressure andu:,i= l,... ,8, represent the unknown velocity values in the eight corner nodes of‘ the element (e). The pressure is taken only in the element centre pb; the interpolation of the velocity components is trilinear; the pressure approximation is constant within an element (M = 1). The approximate function for the velocity has inter-element continuity. The pressure is discontinuous. The numerical procedure developed for the viscous incompressible Navier -Stokes equations uses pressure correction. The decomposition of a vector field is applied to the Navier- Stokes equations’. introducing the auxiliary vector field u”+ “2 and applying implicit time integration yield a numerical algorithm resulting in four calculation steps at each time level tN+ ’ where the nodal velocities u” + ’ and the element oressure 6” + ’ are unknown values (At
=
1” +
I _
irr,:
’
vector 1. calculate the auxiliary the linearized system : [C +
field Cr + ‘/’ frorn
A!k‘l ( tf) + Atk”2Jzf +
2. calculate
the pressure
correction
q”+ ’ : (10)
3. calculate the divergence free d lf 1.. u’ !+I _ d’ + I ;L’ + &(;,,- 1Qq” + 1 4. calculate
the updated
pressure
velocity
field
ill, p” + ’ :
n + I = p” + qn + I P
(12)
In Equations (9~.12) c’ is the mass matrix, Cd is the lumped mass matrix, Kl (zz”) is the linearized convection matrix, E;2 is the diffusion matrix, Qis the gradient matrix, QT is the divergence matrix, cac& is the known right-hand side vector containing given boundary conditions and known velocities. The procedure as outlined here was developed for the efficient numerical treatment of 3D time-dependent fluid flow under low and moderate Reynolds numbers. A description of the method has been published by Perktold” and in modified form by Hilbert”. Wall shear stress Shear stress
stresses in the flow field are described with the tensor, For incompressible fluids and no-slip
J. Riomed. Eng. 1990, Vol. 12, Septembrr
3
P&tile
wall shear stress: K. Perktold and R. Peter
I
branch vessel
Figure 1 90” T-bifurcation model. The arrows indicate the flow direction. D, main vessel diameter; d, branch vessel diameter (d = D/2)
I 1
R
0.5
1
T/TP
I u1.i! / ‘.I .I U
0
0:o
0.1
0.r
0.2
0.S
0.3
0.7
0.6
I
R
0
I \\\\\\1\ 0.8
0.9
U -
4.0
T/TP
-
Figure 2 Velocity pulse waveform (by R.M. Nerem, personal communication) and fully developed profiles in a long straight tube. Normalized quantities, T/TP, time per pulse cycle
conditions stress is
at the
rigid
vessel
wall
the
wall
shear
(13) where u, denotes the tangential velocity normal unit vector at the wall. NUMERICAL
RESULTS
For the 90” T-bifurcation
4
J. Biomed.
and II is the
model
Eng. 1990, Vol. 12, January
treated
in this study
the branch tube diameter is half the main tube diameter. The bifurcation edge is slightly rounded. The model geometry is shown in Figure 2. The flow rate in the branch is 30% of the entrance flow rate during the pulse cycle. The applied finite element subdivision of the computational domain (where symmetry is taken into account) creates 5848 eightnode isoparametric 3D brick elements and results in a total node number of 27 934 for the three velocity components and the pressure. The calculation is carried out under physiological flow conditions. The velocity pulse waveform and the corresponding fully developed velocity profiles at various pulse phases in a long straight tube are shown in Figure 2. The pulse wave contour shows transient reversal flow partially during the pulse cycle. The velocity profiles which are calculated for a long straight tube are used as the time-dependent inflow condition in the main tube. The calculation is performed for the time-mean Reynolds number Re = 85 and the Stokes number 7 = 4.36, where the main vessel radius is taken as the characteristic length (L = 0.2 cm). The chosen Stokes number agrees with the pulse frequency of 85 strokes/min. The kinematic (Newtonian) viscosity of blood is assumed to be v = 0.038 cm’s_ ‘. The averaged inlet velocity is 16 cm s- ‘. The presentation of the numerical results concerns the axial and secondary flow velocity (non-dimensional) as well as the vector field and the contour lines of the pulsatile wall shear stress (in physical units). The magnification factors for the displays are indicated in the figures. A comparison of wall shear values for the slightly rounded corner and a sharp corner at the bifurcation edge is also performed. Velocity
field
The flow velocities in the plane of branching (plane of symmetry) and in the plane perpendicular to the branching plane through the axis of the branch at different time phases during the pulse cycle are illustrated in Figure 3. The diagram shows axial velocity profiles at different positions upstream (marked as A, B) and downstream (marked as C, D, E, F), the branching in the main tube and in the branch (marked as G, H, I) and the vector field of flow velocity during the same time phases. The plots illustrate the flow deviation in the main tube from the centre towards the wall and the skewed profiles with high gradients in the branch for forward directed flow ( T/ ‘TP = 0.5 1, 0.77). Flow separation is observed during deceleration flow (T/TP = 0.77) at the upstream (outer) wail of the branch. In Figure 4 axial velocity in the main tube and in the branch is displayed in 3D representation and the position and orientation of the specified flow crosssections in the flow field is indicated. The marked point at each main tube cross-section corksponds to the vessel wall opposite the branching, and the marked point at each branch cross-section corresponds to the branch wall opposite the flow divider. In Figure 5 the secondary flow velocity at different flow cross-sections and at selected pulse phases is displayed. A significant feature of viscous 3D fluid
Pulratile mall shear strew: 6. Perktold and R. Peter
T/TP
-
0.77
T/TP
-
0.51
T/TP .,ttttttrttfI ..tttttttttr1 .,+tttftfrtt~ .,cttttfrttr,L ,(,tffff~~,,:~~~~----“~ ,a t t t 1 f I~II//%~Z:,“ITI,Z_ttfttttl//~~---------ffffj//f///. :::r::,:.‘, ttfllllltt’
-
0.51
,1tttttttt, +ttttttttt‘ ::Itttttttttr .,4ttttttttta .,4tttttttti. .,,tttttttttl ..*tttffffft+
Figure 3 prrpcndicular
Axial
velocity
to the main
profiles
and velocity
vector
vessel axis (non-dimensional
field in the plane
of symmetry
T/TP
and
in the plane
through
the branch
-
0.12
axis
:I
flow in curved tubes is the appearance of secondary motion in planes normal to the direction of the curved main flow motion. Secondary motion in curved tubes arises as a consequence of the nonuniformity of velocity (high centre stream velocity) in real viscous fluid flow. The secondary motion effect is known to occur in curved tubes where the main flow direction changes. However, in real viscous bifurcation flow the direction of the flow stream changes also and thus secondary flow velocity occurs. Figure 5 shows the flow velocity in different flow cross-sections during the pulse cycle with time phases as indicated. The flow motion in the crosssections results from flow diversion because of the flow division and from the secondary motion effect. The portion of the pure secondary motion produced by the curved fluid stream is not obvious. The reason for this is that the direction of the local axis of the curved fluid stream does not correspond to the geometrical axis which contains the normal vector
on the cross-sections of representation. Flow circulation is observed at the downstream cross-sections for forward directed flow (T/77 = 0.51, 0.77) in the main tube and in the branch vessel. The secondary flow effect promotes the development of asymmetrical velocity profiles whose steep velocity gradients occur in the boundary layer at the branching wall. Wall shear stress The wall shear stress vector field at selected pulse phase angles is shown in Figures 6 and 7. The sign of the vectors corresponds to the tangential velocity direction near the wall. For comparison note the different magnification factors. The wall shear stress presentation uses a plane representation field for the vessel wall. The main vessel is cut open opposite the branching; the branch-cut is downstream of the flow divider along a generating line. Both vessels are unrolled. The correspondence of the two walls is
P&&e
wall shear stms: K. P&old
u
-
2.00
and R. P&t
T/TP
-
I
-
\E
T/TP
u
-
3.00
-
Figure 4
3D representation of axial (G, H, I) are indicated in Figure 3
6
J. Biomed.
Eng. 1990, Vol. 12, January
T/TP velocity
-
-
0.00
u -
1.00
-
u -
1.00
-
0.12
0.51 during
u the pulse
cycle.
-
The specified
3.00
-
flow cross-sections
T/TP
-
0.09
T/TP
-
0.39
T/TP
-
0.77
(A, B), (C, D, E, F) and
.::.;. Cl C
..a
*
r
.
.
.
’
,
.
.
.
.
.
.
.
: _ . . ’ .;
’
: :
l
-’
. . . ,.
[.
. .
. .
_.f’ : . .
.
I:
.
. ;’
.
.
*
.
”
1
. . __..
.
’I
‘, :
,
.
.I
. . . .
.,(
T/TP
- 0.12
c .:-,._ ,;:;;_ . B
.L
,::*
“. ‘.,
. .
_.
.
.
-
.
.
--
..C
*
.
.
..*
c
.
.
,.
L$i3&y
-
-._
,: . .
.
.
-c
. .-:e )
..-
L
.
F
0 ‘.
..-..
. . I
-
,.
.
.
”
.
’
-._.*
-..-a .:...c. . .
,
-
0.50
H
’ . ‘,‘a, .
. .
_
.
.
I
,
.’
,‘:’
’
.
,,,
:2
a. . .
”
I
1
:
C
’
. .
.
‘.
:.
U
-;:*
.
’
T/TP
- 0.51
T/TP
- 0.77
:.
..---.
-\
-4
__--_
_.---. . --..*e---
cd-
--*_-
L .
_---
.-
l
:
:.
c
:
-‘.d ,:.
Pulsatile wall shear slress: K. Perktold and R. Peter
T/TP
xv -
-
15.0
[N/mZl
c
TlTP
L
Figure 6 generating
-
20.0
-
T/TP
0.00
-
J. Biomed.
0.39
[N/m21
T/TP --
L
_(
Wall shear stress vector field in physical units Nm ‘. Plane wall representation: line opposite the branching, and the branch vessel is cut open along the generating
Eng. 1990, Vol.
12, J~IUILW~
5.0
0.
12
[N/m21
,. .,
marked with small arrows in Figure 8 (where the wall shear stress contour lines are presented). The plots illustrate the disturbance of the wall shear stress vector field during the pulse cycle resulting from the branching. Figures 6 and 7 show that zones of high wall shear stress are found in the branching region of the main vessel and around the circumference at the entry region of the branch vessel. The highest wall occurs at maximum shear stress r, = 29.1 Nm-’ flow and is located in the entry region of the branch vessel. Figure 8 illustrates the wall shear stress distribution by means of contour lines. The difference in shear stress between the contour lines is indicated. The numbers are the wall shear stress values in units of Nm- *. The locations of maximum and minimum wall shear stress values are marked with symbols. The contour plots show the large variations in wall shear stress magnitude. The density of the contour lines indicates the local wall shear stress gradient. The location of maximum wall shear stress can be found at the branch vessel entrance region near the flow divider and at the upstream side of the branching in the main vessel. The minimum wall shear stress in the main vessel (at forward directed flow) verified at T/TP = 0.0, 0.51,0.71 occurs at the wall opposite the branching, Comparison of wall shear stress for a slightly rounded and for a sharp corner at the bifurcation edge is seen in Figure 9. The shape of the geometrical
8
-
L
-
-
40.0
-
0.51
IN/m21
the main vessel is cut open along the line extending from the flow divider
model under consideration is of substantial importance in wall shear stress. For 2D bifurcation flow the influence of the rounded flow divider is analysed by Friedman and Ehrlich’ ’ . The effect of radius of curvature of the wall at the flow divider in 3D 90” Tbifurcation models has been investigated by Karino et aL5 under steady flow. As reported there were significant differences in the recirculation zone in the rounded and the square T-junction. In Figure 9 the non-dimensional wall shear stress along the two branching walls indicated as Sl and S2 are shown for a slightly rounded and a square divider. Although the roundness is only slight, the influence in the magnitude of the wall shear stress and in the recirculation zone can be observed. In the case of the mode1 with a rounded edge the peak wall shear stress (at maximum flow rate T/TP = 0.5 1) is significantly reduced. The rounded edge reduces the zones in the branch in which reversed flow occurs.
CONCLUSION These studies were carried out for a 3D mode1 of a 90” T-vessel junction assuming Newtonian viscosity for blood and rigid walls; for local studies in short segments of relatively large arteries these simplifications are permitted. Substantial flow phenomena in large arteries are only slightly influenced as a result of these assumptions. The most important flow parameter is the wall
I
Y
-
_
___,-
___._,~~
_-
_
.
.
4
L^
-
A t,
_.-
.
Puisalile mall shear skw:
K. Perktnid and R. Peter
1 T/TP
t-r‘-
T/TP
- 2.O[N/m21
av
14 --?
- 0.00
A-L
- 0.09
- l.O[N/m21
E I.
5.
l
.
1
min 0.0 ln *
max 9.5
min 0.7 in l max 14.7 in 0
Ln *
mln 0.2 Ln
l
.-
2.
mln 0.5 ln l max 5.2 ln 0
max 6.0 Ln *
I
I
1.
T/TP
Ar,
-
- 0.12
T/TP
1.0[N/m21
m 2.
AL
- 0.39
- 2.01N/mzl
i !
.-.
!5.
, 4
L--
mtn 0.2 ln
l
max 3.6 ln *
1
mln 1.0 tn l max 6.2 In 0
\
P v
TITP
27.
mtn 0.0 In
l
max 24.0
10
J. Biomed.
Eng. 1990, Vol.
contour contour
12, January
l
mln 3.9 In t max 7.7 ln 0
max 11.7 ln 0
r
- 0.51
T/TP
- 3.0[N/m21 .l_
LU.
mln 1.4 In + mex 29.1 Ln 0
In *
Figure 8 Wall shear stress values related to thr selected
nr,
mln 0.4 In
plots lines;
I min 0.0 In
l
max
14.6 in 0
A-L -
- 0.77 2.O[N/m21
mln 0.9 ln l max 21.5 Ln 0
in physical units during the pulse cycle. The numbers give the wall shear location of maximum and minimum wall shear stress values are indicated
stress
Pulmltle ml1 Jhear stres.r: h‘. Perktold and R. Prier I . . ..a*. . . . . , . . . . . . , . . . . ..,. . . . . . . . . . . . .,a..... ,..----..-------, - -a-----.C_C
. . . . . .
I-L Lb
\
s2
Sl
SHARP EDGE
,iiAw
% 20 R. g 0.
-Sl
.R -s* go
..flow detachment . flow reattachment separation
T/TP = 0.51 t
2.00
-
ROUNDED EDGE
a Jb
*>1
-
TENDENCY OF
Pulsatile wall shear stress: K. Perktold and R. Peter
shear stress. The results clearly show the periodically alternating directions of wall shear stress. It is a known fact that in addition to the magnitude of the shear stress the alternating mechanical stress disturbs the function of the endothelial cells. An aspect of great importance in the research of mechanisms in atherogenesis is haemodynamic Caro’s ‘low-shear hypothesis”.From observations it was found that atherosclerotic lesions were localized at the outer walls of bifurcations and T-junctions5. Theoretical and experimental studies under physiological flow conditions indicate that sites where lesions occur correspond with the locations of low flow or low flow recirculation. The detailed investigation of the flow and stress pattern as well as of the pattern of diseases indicates a strong correlation between the sites of low flow and low wall shear stress and the preferred locations for the genesis of atherosclerosis. To enable correlations between haemodynamics and arterial lesions, physiologically relevant geometrical conditions and flow conditions are required. From the clinical point of view it is essential to know that arterial diseases appear first at locations of low flow and low wall shear stress and not at locations of high velocity gradients and high wall shear stress. In human vessel surgery it must be taken into account that areas where the flow is iocaily low should be avoided whenever possible. The mathematical approach permits a detailed analysis of the stress field, and the computer simulation of pulsatile flow yields high resolution. The detailed investigation of local patterns of flow separation and wall shear stress shows the topographic distribution of low and high velocity gradients and the non-uniformity of wall shear stress. In the region of branching, the wall shear stress varies dramatically. Zones of low and high wall shear stress are found in close proximity to the leading edge of the flow divider and to the entrance region of the branch. A comparison of our results using three-dimensional model geometry and a physiological pulse waveform with other theoretical or experimental wall shear stress results is not practicable, because no
12
J. Biomed. Eng. 1990, Vol.
12, January
detailed pulsatile wall shear stress results in a 90” Tvessel bifurcation are available in the literature. As far as our results can be compared with experimental results, good qualitative agreement can be observed.
ACKNOWLEDGEMENT This investigation is supported by Research Fund (Fonds zur..Forderung schaftlichen Forschung in Osterreich), P 6794 P, Vienna, Austria.
the Austrian der wissenProject Nr.
REFERENCES 1. Caro
2.
3.
4.
5.
6.
CG, Fitz-Gerald JM, Schroter RC. Atheroma and arterial wall shear: observation, correlation and proposal of a shear dependent mass transfer mechanism of atherogenesis. pr’roc R S0r 1971; 177: 109-59. Fry DL. Acute vascular endothelial changes associated with increased blood velocity gradients. C’irc Res 1968; 22: 165597. Ku DN, Liepsch D. The effect of non-Newtonian viscoelasticity and wall elasticity on Row at a 90” bifurcation. BiorheololCy 1986; 23: 359970. Talukter N, Nerem RM. Flow characteristics in vascular graft models. DZge,\l First Int Conf Mech Med Biol, Aachen. 1978; VII: 281-4. Karino T, Motomiya M, Goldsmith HL. Flow patterns in model and natural branching vessels. In: Schettler G, Ncrem RM, Schmid-Schonbein H, Marl H, Diehm H, eds. Fluid l>ynamics as a Locatiriy Factor for Atherosclerosis. Berlin: Springer-Vcrlag, 1983; 60-70. Liepsch D, Moravec S. Pulsatile flow of non-Newtonian fluid in distensible models of human arteries. Riorheology 1984;21: 571-86.
7. Temam R. .Nuvier-Stokes Equations. North-Holland, Amsterdam: 1979. 8. Girault VG, Raviart PA. Finite Elemenl Methods for NauierStokes Eguations. Berlin: Springer-Verlag, 1980. 9. Perktold K. On .Numrrical Simulation of ‘Three-Dimensional Physiological Flow Problems, Ber Math Stat Sektion, Forschungsges. Joanneum, Graz, Nr.280, 1987. 10. Hilbert D. An efficient Navier-Stokes solver and its application to fluid flow in elastic tubes, In: .Numeriral MethodJ, colloquia Socielatis Jdnos Bo!pai, 50, North Holland, 1987; 423-31. 11. Friedman, MH, Ehrlich LW. Numerical simulation of aortic bifurcation flows: the effect offlow divider curvature. 3 Biomech 1984; 17: 881-8.
of Mathematics,
Received
April
1989,
Technical
accepted
July
University
Graz, Steyrergasse
17,A-8010 Graz, Austria
1989
ABSTRACT The structure of pulsatile bloodjow and wall shear stress in a 90” T-bifurcation model is analyxd numerically. The nonlinear Navier-Stokes equations for time-dependent incompressible .hewtonian JluidJow are approximated using a newly developed pressure correction, Jinite element method. The wall shear stress is calculated from the Jinite element velocifv jield. The investigation shows viscousJow phenomena such asjow separation and stagnation and the distribution of high and low wall shear stress during the pulse cycle. Furthermore, the effect of a sharp corner at the bifurcation edge on the wall shear stress is anaivsed. Detailed 1ocalJow investigation is required to examine&id &namic contribution to the development of arterial diseases such as atherosclerosis and thrombosis. Keywords:
Arterial
bifurcation,
pulsatile
blood
flow, numerical
INTRODUCTION Many biomechanical studies concentrate on local flow behaviour and wall shear stress in blood vessels. The goal of such investigations is the determination of the haemodynamic factors in the development of vascular diseases. Clinical observations show that pathological vessel alterations occur preferentially in zones where the flow is locally disturbed and extreme wall shear stresses appear. The detailed analysis of flow phenomena and comparison with the pattern of diseases should contribute to the comprehension of the haemodynamic influence on pathoThe effects of interest mainly logical alterations. concern the occurrence of flow recirculation and flow stagnation zones as well as areas of high velocity gradient where the endothelium can be damaged and platelets can be stimulated. The recirculation zones where stasis results are treated in the ‘lowshear’ theory of Caro et al.‘, while the aspects of elevated shear stresses are expressed in the ‘highshear’ theory of Fry2. The main interest in local flow dynamics has been concentrated on vessel bifurcations and vessel bends, where diseases often occur. Several recent papers are devoted to experimental studies of the flow structure in 90” T-bifurcation models under different assumpthe results of lasertions. In these papers, or hot-film measurements” Doppler measurements” and of experimental flow visualizations5,6 are presented. We examine here the flow structure and wall shear stress under physiological pulsatile flow conditions using computer simulation. The mathematical approach yields detailed information on the local pulsatile arterial flow. Efficient numerical methods to approximate the viscous flow equations and modern computer technology permit the development Correspondcncc
and rrprint
requests
to: Karl Prrktold
f,1990 Butterworth & Co (Publishers) Ltd 0141-5425/90/010002-II $03.00 2 J. Biomcd. Eng. 1990. Vol. 12, January
investigation,
flow pattern,
wall shear
stress
of physiologically realistic models. Amongst the numerical procedures which can be applied, the finite element method has the most advantages and because of its flexibility it is well suited to the analysis of complicated arterial flow. The appropriate application of the finite element method on viscous flow equations can produce accurate numerical results for fluid velocity, an essential requirement for the calculation of realistic wall shear stresses during the pulse cycle. In this study, numerical simulation is performed using a three-dimensional 90” T-junction model. The pulsatile flow pattern and the wall shear stress distribution are analysed corresponding to a physiological velocity pulse waveform occurring in the left coronary artery. The calculation is carried out for a mean inflow Reynolds number Re = 80, whereby the radius of the parent vessel is used in the definition. The flow in the branch is assumed to be 300/, of the entrance flow. Numerical results are presented for non-dimensional axial and secondary flow velocity as well as for wall shear stress in physical quantities.
MATHEMATICAL MODEL Computer simulation of local blood flow phenomena in arteries requires the full Navier4tokes equations; simplifications of the non-linear flow equations cannot be successful because of the occurrence of flow separation and secondary motion. In this study blood is assumed to be an incompressible, homogeneous, Newtonian fluid ; its specific rheological properties are not taken into account. This idealization can be made for the calculation of blood flow in large arteries where the shear rate is relatively high. The vessel wall is assumed to be rigid; an assumption that is admissible for the determination of local flow
phenomena in short artery segments where pulse wave propagation is not treated. In the primitive variables flow velocity and pressure, the pulsatile, incompressible, Newtonian, 3DNavier Stokes equations can be written as
presses the velocity components and the pressure in terms ofinterpolation functions over each element of the subdivided flow domain. Here 3D brick elements are used. The appropriate solution for velocity and pressure in an element (e) takes the form
1)
v.u=o
(2)
the normalized where u = (u, ZJ,x!) 7 and p denote non-dimensional flow velocity and pressure, respectively. The velocity components U, G, UJ are related to the Gartesian coordinates x,y, Z. The normalization is done with respect to a characteristic length I,, a characteristic velocity U,, and the kinematic viscosity I’. Re and q denote the time-mean Reynolds number and the Stokes number, respectively. and are defined as R e = <:li,!
(3a)
V
(3b) The non-dimensional and pressure are
units
length,
time,
velocity
4) physical units. The other mathewhere - indicates matical symbols have standard meanings: w is angular frequency and e is constant fluid density. For the flow problems which are treated here the boundary of the flow domain partitions into two sections with different types of boundary conditions. At the inflow boundary and at the vessel wall the velocity is given where the flow at the entrance is prescribed with time-dependent velocity profiles, and at the rigid wall the no-slip condition is valid %,,,. = s(l
and
UWll = 0
(5)
The outflow is assumed traction-free. The condition of vanishing surface traction can be described mathematically as (6) where II is the normal unit vector at the outflow boundary. The choice of appropriate boundary conditions at the flow entrance and exit often produces problems. It is required that the assumed conditions do not disturb the numerical results in the flow domain of interest. Numerical
method
The Navier-Stokes equations are spatially discretized using the Galerkin finite element method where the differential equations are reformulated with an equivalent variational formulation in appropriate The finite element method exfunction spaces7.
(7) P’(4 =
MP;,
(4
(8)
where .,V,, i = 1,. . ,8 and ;21 are the element interpolation functions for the velocity and the pressure andu:,i= l,... ,8, represent the unknown velocity values in the eight corner nodes of‘ the element (e). The pressure is taken only in the element centre pb; the interpolation of the velocity components is trilinear; the pressure approximation is constant within an element (M = 1). The approximate function for the velocity has inter-element continuity. The pressure is discontinuous. The numerical procedure developed for the viscous incompressible Navier -Stokes equations uses pressure correction. The decomposition of a vector field is applied to the Navier- Stokes equations’. introducing the auxiliary vector field u”+ “2 and applying implicit time integration yield a numerical algorithm resulting in four calculation steps at each time level tN+ ’ where the nodal velocities u” + ’ and the element oressure 6” + ’ are unknown values (At
=
1” +
I _
irr,:
’
vector 1. calculate the auxiliary the linearized system : [C +
field Cr + ‘/’ frorn
A!k‘l ( tf) + Atk”2Jzf +
2. calculate
the pressure
correction
q”+ ’ : (10)
3. calculate the divergence free d lf 1.. u’ !+I _ d’ + I ;L’ + &(;,,- 1Qq” + 1 4. calculate
the updated
pressure
velocity
field
ill, p” + ’ :
n + I = p” + qn + I P
(12)
In Equations (9~.12) c’ is the mass matrix, Cd is the lumped mass matrix, Kl (zz”) is the linearized convection matrix, E;2 is the diffusion matrix, Qis the gradient matrix, QT is the divergence matrix, cac& is the known right-hand side vector containing given boundary conditions and known velocities. The procedure as outlined here was developed for the efficient numerical treatment of 3D time-dependent fluid flow under low and moderate Reynolds numbers. A description of the method has been published by Perktold” and in modified form by Hilbert”. Wall shear stress Shear stress
stresses in the flow field are described with the tensor, For incompressible fluids and no-slip
J. Riomed. Eng. 1990, Vol. 12, Septembrr
3
P&tile
wall shear stress: K. Perktold and R. Peter
I
branch vessel
Figure 1 90” T-bifurcation model. The arrows indicate the flow direction. D, main vessel diameter; d, branch vessel diameter (d = D/2)
I 1
R
0.5
1
T/TP
I u1.i! / ‘.I .I U
0
0:o
0.1
0.r
0.2
0.S
0.3
0.7
0.6
I
R
0
I \\\\\\1\ 0.8
0.9
U -
4.0
T/TP
-
Figure 2 Velocity pulse waveform (by R.M. Nerem, personal communication) and fully developed profiles in a long straight tube. Normalized quantities, T/TP, time per pulse cycle
conditions stress is
at the
rigid
vessel
wall
the
wall
shear
(13) where u, denotes the tangential velocity normal unit vector at the wall. NUMERICAL
RESULTS
For the 90” T-bifurcation
4
J. Biomed.
and II is the
model
Eng. 1990, Vol. 12, January
treated
in this study
the branch tube diameter is half the main tube diameter. The bifurcation edge is slightly rounded. The model geometry is shown in Figure 2. The flow rate in the branch is 30% of the entrance flow rate during the pulse cycle. The applied finite element subdivision of the computational domain (where symmetry is taken into account) creates 5848 eightnode isoparametric 3D brick elements and results in a total node number of 27 934 for the three velocity components and the pressure. The calculation is carried out under physiological flow conditions. The velocity pulse waveform and the corresponding fully developed velocity profiles at various pulse phases in a long straight tube are shown in Figure 2. The pulse wave contour shows transient reversal flow partially during the pulse cycle. The velocity profiles which are calculated for a long straight tube are used as the time-dependent inflow condition in the main tube. The calculation is performed for the time-mean Reynolds number Re = 85 and the Stokes number 7 = 4.36, where the main vessel radius is taken as the characteristic length (L = 0.2 cm). The chosen Stokes number agrees with the pulse frequency of 85 strokes/min. The kinematic (Newtonian) viscosity of blood is assumed to be v = 0.038 cm’s_ ‘. The averaged inlet velocity is 16 cm s- ‘. The presentation of the numerical results concerns the axial and secondary flow velocity (non-dimensional) as well as the vector field and the contour lines of the pulsatile wall shear stress (in physical units). The magnification factors for the displays are indicated in the figures. A comparison of wall shear values for the slightly rounded corner and a sharp corner at the bifurcation edge is also performed. Velocity
field
The flow velocities in the plane of branching (plane of symmetry) and in the plane perpendicular to the branching plane through the axis of the branch at different time phases during the pulse cycle are illustrated in Figure 3. The diagram shows axial velocity profiles at different positions upstream (marked as A, B) and downstream (marked as C, D, E, F), the branching in the main tube and in the branch (marked as G, H, I) and the vector field of flow velocity during the same time phases. The plots illustrate the flow deviation in the main tube from the centre towards the wall and the skewed profiles with high gradients in the branch for forward directed flow ( T/ ‘TP = 0.5 1, 0.77). Flow separation is observed during deceleration flow (T/TP = 0.77) at the upstream (outer) wail of the branch. In Figure 4 axial velocity in the main tube and in the branch is displayed in 3D representation and the position and orientation of the specified flow crosssections in the flow field is indicated. The marked point at each main tube cross-section corksponds to the vessel wall opposite the branching, and the marked point at each branch cross-section corresponds to the branch wall opposite the flow divider. In Figure 5 the secondary flow velocity at different flow cross-sections and at selected pulse phases is displayed. A significant feature of viscous 3D fluid
Pulratile mall shear strew: 6. Perktold and R. Peter
T/TP
-
0.77
T/TP
-
0.51
T/TP .,ttttttrttfI ..tttttttttr1 .,+tttftfrtt~ .,cttttfrttr,L ,(,tffff~~,,:~~~~----“~ ,a t t t 1 f I~II//%~Z:,“ITI,Z_ttfttttl//~~---------ffffj//f///. :::r::,:.‘, ttfllllltt’
-
0.51
,1tttttttt, +ttttttttt‘ ::Itttttttttr .,4ttttttttta .,4tttttttti. .,,tttttttttl ..*tttffffft+
Figure 3 prrpcndicular
Axial
velocity
to the main
profiles
and velocity
vector
vessel axis (non-dimensional
field in the plane
of symmetry
T/TP
and
in the plane
through
the branch
-
0.12
axis
:I
flow in curved tubes is the appearance of secondary motion in planes normal to the direction of the curved main flow motion. Secondary motion in curved tubes arises as a consequence of the nonuniformity of velocity (high centre stream velocity) in real viscous fluid flow. The secondary motion effect is known to occur in curved tubes where the main flow direction changes. However, in real viscous bifurcation flow the direction of the flow stream changes also and thus secondary flow velocity occurs. Figure 5 shows the flow velocity in different flow cross-sections during the pulse cycle with time phases as indicated. The flow motion in the crosssections results from flow diversion because of the flow division and from the secondary motion effect. The portion of the pure secondary motion produced by the curved fluid stream is not obvious. The reason for this is that the direction of the local axis of the curved fluid stream does not correspond to the geometrical axis which contains the normal vector
on the cross-sections of representation. Flow circulation is observed at the downstream cross-sections for forward directed flow (T/77 = 0.51, 0.77) in the main tube and in the branch vessel. The secondary flow effect promotes the development of asymmetrical velocity profiles whose steep velocity gradients occur in the boundary layer at the branching wall. Wall shear stress The wall shear stress vector field at selected pulse phase angles is shown in Figures 6 and 7. The sign of the vectors corresponds to the tangential velocity direction near the wall. For comparison note the different magnification factors. The wall shear stress presentation uses a plane representation field for the vessel wall. The main vessel is cut open opposite the branching; the branch-cut is downstream of the flow divider along a generating line. Both vessels are unrolled. The correspondence of the two walls is
P&&e
wall shear stms: K. P&old
u
-
2.00
and R. P&t
T/TP
-
I
-
\E
T/TP
u
-
3.00
-
Figure 4
3D representation of axial (G, H, I) are indicated in Figure 3
6
J. Biomed.
Eng. 1990, Vol. 12, January
T/TP velocity
-
-
0.00
u -
1.00
-
u -
1.00
-
0.12
0.51 during
u the pulse
cycle.
-
The specified
3.00
-
flow cross-sections
T/TP
-
0.09
T/TP
-
0.39
T/TP
-
0.77
(A, B), (C, D, E, F) and
.::.;. Cl C
..a
*
r
.
.
.
’
,
.
.
.
.
.
.
.
: _ . . ’ .;
’
: :
l
-’
. . . ,.
[.
. .
. .
_.f’ : . .
.
I:
.
. ;’
.
.
*
.
”
1
. . __..
.
’I
‘, :
,
.
.I
. . . .
.,(
T/TP
- 0.12
c .:-,._ ,;:;;_ . B
.L
,::*
“. ‘.,
. .
_.
.
.
-
.
.
--
..C
*
.
.
..*
c
.
.
,.
L$i3&y
-
-._
,: . .
.
.
-c
. .-:e )
..-
L
.
F
0 ‘.
..-..
. . I
-
,.
.
.
”
.
’
-._.*
-..-a .:...c. . .
,
-
0.50
H
’ . ‘,‘a, .
. .
_
.
.
I
,
.’
,‘:’
’
.
,,,
:2
a. . .
”
I
1
:
C
’
. .
.
‘.
:.
U
-;:*
.
’
T/TP
- 0.51
T/TP
- 0.77
:.
..---.
-\
-4
__--_
_.---. . --..*e---
cd-
--*_-
L .
_---
.-
l
:
:.
c
:
-‘.d ,:.
Pulsatile wall shear slress: K. Perktold and R. Peter
T/TP
xv -
-
15.0
[N/mZl
c
TlTP
L
Figure 6 generating
-
20.0
-
T/TP
0.00
-
J. Biomed.
0.39
[N/m21
T/TP --
L
_(
Wall shear stress vector field in physical units Nm ‘. Plane wall representation: line opposite the branching, and the branch vessel is cut open along the generating
Eng. 1990, Vol.
12, J~IUILW~
5.0
0.
12
[N/m21
,. .,
marked with small arrows in Figure 8 (where the wall shear stress contour lines are presented). The plots illustrate the disturbance of the wall shear stress vector field during the pulse cycle resulting from the branching. Figures 6 and 7 show that zones of high wall shear stress are found in the branching region of the main vessel and around the circumference at the entry region of the branch vessel. The highest wall occurs at maximum shear stress r, = 29.1 Nm-’ flow and is located in the entry region of the branch vessel. Figure 8 illustrates the wall shear stress distribution by means of contour lines. The difference in shear stress between the contour lines is indicated. The numbers are the wall shear stress values in units of Nm- *. The locations of maximum and minimum wall shear stress values are marked with symbols. The contour plots show the large variations in wall shear stress magnitude. The density of the contour lines indicates the local wall shear stress gradient. The location of maximum wall shear stress can be found at the branch vessel entrance region near the flow divider and at the upstream side of the branching in the main vessel. The minimum wall shear stress in the main vessel (at forward directed flow) verified at T/TP = 0.0, 0.51,0.71 occurs at the wall opposite the branching, Comparison of wall shear stress for a slightly rounded and for a sharp corner at the bifurcation edge is seen in Figure 9. The shape of the geometrical
8
-
L
-
-
40.0
-
0.51
IN/m21
the main vessel is cut open along the line extending from the flow divider
model under consideration is of substantial importance in wall shear stress. For 2D bifurcation flow the influence of the rounded flow divider is analysed by Friedman and Ehrlich’ ’ . The effect of radius of curvature of the wall at the flow divider in 3D 90” Tbifurcation models has been investigated by Karino et aL5 under steady flow. As reported there were significant differences in the recirculation zone in the rounded and the square T-junction. In Figure 9 the non-dimensional wall shear stress along the two branching walls indicated as Sl and S2 are shown for a slightly rounded and a square divider. Although the roundness is only slight, the influence in the magnitude of the wall shear stress and in the recirculation zone can be observed. In the case of the mode1 with a rounded edge the peak wall shear stress (at maximum flow rate T/TP = 0.5 1) is significantly reduced. The rounded edge reduces the zones in the branch in which reversed flow occurs.
CONCLUSION These studies were carried out for a 3D mode1 of a 90” T-vessel junction assuming Newtonian viscosity for blood and rigid walls; for local studies in short segments of relatively large arteries these simplifications are permitted. Substantial flow phenomena in large arteries are only slightly influenced as a result of these assumptions. The most important flow parameter is the wall
I
Y
-
_
___,-
___._,~~
_-
_
.
.
4
L^
-
A t,
_.-
.
Puisalile mall shear skw:
K. Perktnid and R. Peter
1 T/TP
t-r‘-
T/TP
- 2.O[N/m21
av
14 --?
- 0.00
A-L
- 0.09
- l.O[N/m21
E I.
5.
l
.
1
min 0.0 ln *
max 9.5
min 0.7 in l max 14.7 in 0
Ln *
mln 0.2 Ln
l
.-
2.
mln 0.5 ln l max 5.2 ln 0
max 6.0 Ln *
I
I
1.
T/TP
Ar,
-
- 0.12
T/TP
1.0[N/m21
m 2.
AL
- 0.39
- 2.01N/mzl
i !
.-.
!5.
, 4
L--
mtn 0.2 ln
l
max 3.6 ln *
1
mln 1.0 tn l max 6.2 In 0
\
P v
TITP
27.
mtn 0.0 In
l
max 24.0
10
J. Biomed.
Eng. 1990, Vol.
contour contour
12, January
l
mln 3.9 In t max 7.7 ln 0
max 11.7 ln 0
r
- 0.51
T/TP
- 3.0[N/m21 .l_
LU.
mln 1.4 In + mex 29.1 Ln 0
In *
Figure 8 Wall shear stress values related to thr selected
nr,
mln 0.4 In
plots lines;
I min 0.0 In
l
max
14.6 in 0
A-L -
- 0.77 2.O[N/m21
mln 0.9 ln l max 21.5 Ln 0
in physical units during the pulse cycle. The numbers give the wall shear location of maximum and minimum wall shear stress values are indicated
stress
Pulmltle ml1 Jhear stres.r: h‘. Perktold and R. Prier I . . ..a*. . . . . , . . . . . . , . . . . ..,. . . . . . . . . . . . .,a..... ,..----..-------, - -a-----.C_C
. . . . . .
I-L Lb
\
s2
Sl
SHARP EDGE
,iiAw
% 20 R. g 0.
-Sl
.R -s* go
..flow detachment . flow reattachment separation
T/TP = 0.51 t
2.00
-
ROUNDED EDGE
a Jb
*>1
-
TENDENCY OF
Pulsatile wall shear stress: K. Perktold and R. Peter
shear stress. The results clearly show the periodically alternating directions of wall shear stress. It is a known fact that in addition to the magnitude of the shear stress the alternating mechanical stress disturbs the function of the endothelial cells. An aspect of great importance in the research of mechanisms in atherogenesis is haemodynamic Caro’s ‘low-shear hypothesis”.From observations it was found that atherosclerotic lesions were localized at the outer walls of bifurcations and T-junctions5. Theoretical and experimental studies under physiological flow conditions indicate that sites where lesions occur correspond with the locations of low flow or low flow recirculation. The detailed investigation of the flow and stress pattern as well as of the pattern of diseases indicates a strong correlation between the sites of low flow and low wall shear stress and the preferred locations for the genesis of atherosclerosis. To enable correlations between haemodynamics and arterial lesions, physiologically relevant geometrical conditions and flow conditions are required. From the clinical point of view it is essential to know that arterial diseases appear first at locations of low flow and low wall shear stress and not at locations of high velocity gradients and high wall shear stress. In human vessel surgery it must be taken into account that areas where the flow is iocaily low should be avoided whenever possible. The mathematical approach permits a detailed analysis of the stress field, and the computer simulation of pulsatile flow yields high resolution. The detailed investigation of local patterns of flow separation and wall shear stress shows the topographic distribution of low and high velocity gradients and the non-uniformity of wall shear stress. In the region of branching, the wall shear stress varies dramatically. Zones of low and high wall shear stress are found in close proximity to the leading edge of the flow divider and to the entrance region of the branch. A comparison of our results using three-dimensional model geometry and a physiological pulse waveform with other theoretical or experimental wall shear stress results is not practicable, because no
12
J. Biomed. Eng. 1990, Vol.
12, January
detailed pulsatile wall shear stress results in a 90” Tvessel bifurcation are available in the literature. As far as our results can be compared with experimental results, good qualitative agreement can be observed.
ACKNOWLEDGEMENT This investigation is supported by Research Fund (Fonds zur..Forderung schaftlichen Forschung in Osterreich), P 6794 P, Vienna, Austria.
the Austrian der wissenProject Nr.
REFERENCES 1. Caro
2.
3.
4.
5.
6.
CG, Fitz-Gerald JM, Schroter RC. Atheroma and arterial wall shear: observation, correlation and proposal of a shear dependent mass transfer mechanism of atherogenesis. pr’roc R S0r 1971; 177: 109-59. Fry DL. Acute vascular endothelial changes associated with increased blood velocity gradients. C’irc Res 1968; 22: 165597. Ku DN, Liepsch D. The effect of non-Newtonian viscoelasticity and wall elasticity on Row at a 90” bifurcation. BiorheololCy 1986; 23: 359970. Talukter N, Nerem RM. Flow characteristics in vascular graft models. DZge,\l First Int Conf Mech Med Biol, Aachen. 1978; VII: 281-4. Karino T, Motomiya M, Goldsmith HL. Flow patterns in model and natural branching vessels. In: Schettler G, Ncrem RM, Schmid-Schonbein H, Marl H, Diehm H, eds. Fluid l>ynamics as a Locatiriy Factor for Atherosclerosis. Berlin: Springer-Vcrlag, 1983; 60-70. Liepsch D, Moravec S. Pulsatile flow of non-Newtonian fluid in distensible models of human arteries. Riorheology 1984;21: 571-86.
7. Temam R. .Nuvier-Stokes Equations. North-Holland, Amsterdam: 1979. 8. Girault VG, Raviart PA. Finite Elemenl Methods for NauierStokes Eguations. Berlin: Springer-Verlag, 1980. 9. Perktold K. On .Numrrical Simulation of ‘Three-Dimensional Physiological Flow Problems, Ber Math Stat Sektion, Forschungsges. Joanneum, Graz, Nr.280, 1987. 10. Hilbert D. An efficient Navier-Stokes solver and its application to fluid flow in elastic tubes, In: .Numeriral MethodJ, colloquia Socielatis Jdnos Bo!pai, 50, North Holland, 1987; 423-31. 11. Friedman, MH, Ehrlich LW. Numerical simulation of aortic bifurcation flows: the effect offlow divider curvature. 3 Biomech 1984; 17: 881-8.