Wall shear stress distribution in the human carotid siphon during pulsatile flow

J. Biondmics Printed

Vol. 21. No. 8. pp. 663-671.

0021-9290/X8 $3.00 + .oO Pergamon Press pk

1988

in Great Britain

WALL SHEAR STRESS DISTRIBUTION IN THE HUMAN CAROTID SIPHON DURING PULSATILE FLOW K. PERKTOLD, H. FLORIAN, D. HILBERT and R. PETER Institute of Mathematics, Technical University Graz, Graz, Austria Abstract-Wall shear stress distribution in the carotid siphon, which is a multiple curved segment of the internal carotid artery, is investigated numerically under physiological flow conditions. The computer simulation of flow through the model segment is based on the time-dependent, three-dimensional Navier-Stokes equations, solved numerically with a finite element method. The study shows the behavior of the wall shear stress-vector field and identities the zones of high and low wall shear stress values during the cardiac cycle.

INTRODUCTION

Since the early seventies it has been accepted that hemodynamic factors are of importance in the initi-

ation and development of atherosclerotic lesions. In many independently performed investigations it has been proved that arterial lesions occur preferentially in certain arterial zones which are related to the local hemodynamics, whereby high and low values of flow parameters velocity and velocity gradient and of flow induced mechanical stress are of relevance. According to Fry’s shear stress hypothesis (Fry, 1968), early lesions are observed in regions with high velocity gradient, i.e. high wall shear stress. According to Caro’s boundary layer diffusion hypothesis (Caro et al., 1971), early atherosclerotic lesions tend to occur in regions with low velocity, i.e. low wall shear stress. Since the publication of these hypotheses, many investigations have been carried out which should explain the hemodynamic influence on arterial disease processes. Pedley (1980) summarizes results up to 1978. A survey on fluid mechanics in atherogenesis up to 1979 is presented by Nerem and Cornhill (1980). In the proceedings of the specialists meeting on ‘Fluid Dynamics as a Localizing Factor for Atherosclerosis’ held at Heidelberg, F.R.G.. June 1982, the status of the knowledge is documented up to 1982 (Schettler et al., 1983). More recent aspects are described in SchmidSchGnbein and Wurzinger (1986). In spite of substantial efforts, no detailed correlation between fluid dynamic parameters or hydrodynamic forces and the development of atherosclerosis has been found to date. One of the requirements for the explanation of the complex mechanism is the knowledge of reliable values for the hemodynamic factors of interest under realistic, physiologically relevant flow conditions. The most important factor among these seems to be the wall shear stress. The determination of this quantity is rather difficult. Conventional experimental techniques determine wall shear stress values by measuring the velocity close to the wall. Relatively new is the Receiced

7 April 1987; in

recisedform

10 NoDember 1987.

application of an electrochemical technique for the direct measurement of wall shear stress (Lutz et al., 1977; Pei et al., 1985). Another approach uses mathematical models to study local flow and stress patterns, whereby the fluid dynamic equations governing the flow in large arteries are solved numerically. In the past years mostly simplified two-dimensional models have been analysed. Useful information is obtained from these implifications for basic studies (Friedman and Ehrlich, 1984; Perktold et al., 1984). A significant disadvantage of using two-dimensional idealizations is that the important secondary motion effect is not taken into account. The study of this effect requires realistic three-dimensional models. For threedimensional viscous flow problems substantial computer efforts are needed. Thus, we have developed an economical finite element method based on the velocity-pressure formulation of the Navier-Stokes equations. The procedure uses a decomposition theorem for a vector field, which allows the uncoupling of pressure and velocity. Here the wall shear stress pattern in the carotid siphon, which is a multiple curved section of the human internal carotid artery, located at the base of the skull, is analysed numerically. This segment is a typical curvature site in the human artery system with high predilection to the development of lesions. MODEL

ASSUMPTIONS

AND MATHEMATICAL

EQUATIONS

The carotid siphon is idealized here by means of a curved tube with three bends and constant internal tube diameter (D = 5.6 mm). The geometrical model is constructed as listed in Table 1. The segment sections are seen from the flow entrance. Figure 1 shows the front view, plan view and side view of the geometrical model in the Cartesian coordinate system and the pulse wave form at the entrance. The Row direction through the model segment is illustrated by arrows. The various crosssections are denoted by Sl, . . , S5, , SSO. The 663

K.

664

PERKTOLD~C

Table 1

Segment section

Plane of Relative Relative curvature, radius of length curvature of turned clockthe axis RJD wise from 0” Angle &ID

65” Bend Bl Straight tube Ll 48” Bend B2 Straight tube L2 Bend B3 145” Straight tube L3

al.

where II is velocity and p is pressure. The relations between the normalized and the physical variables, symbolized by an overbar, are as follows

x=-x,

0

1.6

t=ol,

u=-p

1 0

1 _

0.4 2.2

90”

1.0

212”

P’---yP PUO

1.0

where L denotes the reference length, U, the reference velocity, w angular frequency and p density. The Reynolds number and the Stokes number are defined as

1.5

cross-section of the flow entrance is Sl. The marked point Pl and the generating line M, as indicated give the orientation of the other cross-sections. The computer simulation of pulsatile bloo& flow through the segment is carried out on the basis of the time-dependent three-dimensional Navier-Stokes equations for incompressible Newtonian fluid flow. The assumption that blood is a Newtonian fluid is acceptable for local flow-calculations in large arteries. Using the time-averaged inflow Reynolds number Re and the Stokes number 11the equation system can be written in normalized non-dimensional variables 2

&-$+(u.V)u+Vp-;Au

l_

L

= 0

v.u=o

(1)

Re=-,

PLUO

q = (2LZpW//l)“2,

p

respectively; p is dynamic viscosity. The normalization of the variables is realized through the tube radius as reference length and the time-mean and spatial-mean inflow velocity as reference velocity. The specified boundary conditions for the Navier-Stokes equations are: time-dependent fully developed velocity at the inlet cross-section Sl, belonging to the given pulse wave form, zero-velocity condition at the rigid wall and zero-stress condition at the outflow cross-section. The idealized inflow velocity profiles are calculated in a pre-processor step. The pulse wave form corresponds to the physiological pulse wave presented by Bharadvaj et al. (1982). The zero-stress condition at the outflow boundary can be

FRONT VIEY

B3 PULSE

WAVE

FORM

normalized

PLAN VIEU

(3)

Z

I .oo

.20

.40

.60

T/TP .BO ' v.00

Fig. 1. Model geometry of the human carotid siphon in the Cartesianco-ordinate system and pulse wave form; Bl, B2, B3 denote the three bends, Sl, _ . _ , ,530are cross-sections, M is a generating line, Pl, . . , P50 are intersection points of the line M and the circumferenceof Sl, . . , S50; TP denotes the normalized time per one cardiac cycle.

665

Wall shear stress distribution

expressed mathematically

NUMERICAL

by i

-pn+zG

au =

0

(4)

where n = (n,, n,,, n,)rdenotes the outward pointing unit vector at the outflow boundary. Condition equation (4) results from the restriction of the stress tensor at the outflow boundary section, where vanishing normal and tangential forces are assumed. The stress tensor describes the forces within the flow field.

The definition of an appropriate condition at the outflow boundary is always difficult. It is necessary that the assumed condition does not influence the flow field within the domain of interest. Often it makes sense to assume the zero-stress condition at an artificial outflow boundary, which is located downstream from the actual outflow cross-section to guarantee that the boundary condition does not affect the flow pattern. Such an assumption is applied here for the carotid siphon. The numerical calculation of three-dimensional pulsatile flow requires great computer efforts. Thus, an efficient method has been developed for approximation of the Navier-Stokes equations. The method is based upon the decomposition of a vector field and the finite element method. A vector field can be decomposed into a divergence-free field with zero normal component at the boundary of the flow domain and into an irrotational field. This fact can be applied to the equation of motion (Chorin, 1968). A detailed presentation of the numerical method, developed here is given by Perktold (1985). The mechanical stress, which occurs in the flow field, is described by the stress tensor. The reduction of the stress tensor for incompressible flow and no-slip condition at the wall leads to the equation for the wall shear stress, which is in non-dimensional form

where u, is the velocity in the tangential plane at a point on the wall, n is the unit normal vector.

El I’,

-

Ll

I32

RESULTS

.

The numerical calculation is carried out for the Reynolds number Re = 100,defined with the tube radius, and the Stokes number q = 6.12, i.e. pulse frequency is 80 min - ‘. The finite element subdivision of the computational domain uses eight node isoparametric bricks and results in 6696 elements with 7755 nodes for each velocity component and 6696 pressure nodes. This means a total node number of 29,961 per time step. The numerical results are displayed in non-dimensional units. For presentation of the wall shear stress the axis of the siphon model is straightened and the wall is cut open along the generating line M and is unrolled. By doing this the siphon wall is presented in rectangular form. The length of the rectangle is the length of the segment axis and the width is the circumference of the tube (Fig. 2). An outline of velocity results is shown in Fig. 3 and Fig. 4. Axial and secondary velocities are presented for maximum flow rate (T/TP= 0.22) and minimum flow rate (T/TP= 0.4) at different cross-sections. The axial velocity is displayed at the angle positions - 45”, O”, 45”, 90”. The indicated angle positions at the different cross-sections are measured in the local plane coordinate systems, where the mark 90” corresponds to the generating line M. The axial velocity profiles are presented in the direction of the arrows from - 1 to 1. In curved sections a velocity component occurs in the plane perpendicular to the axis. In the segment curved three times flow is diverted thrice and strong secondary circulation arises (Fig. 4). Secondary circulation infiuences the distribution of axial velocity. Comprehensive numerical results for the flow pattern in the carotid siphon are given in Perktold et al. (1987). Figure 5 shows the vector field of the wall shear stress during the pulse cycle for specified phase angles: T/TP= 0.0, 0.13, 0.22, 0.3, 0.4, 0.5. The stress vectors are indicated by thin lines extending from marked points. The wall shear stress varies with time throughout the pulse cycle and with position. During the period the values are largely different. while the

83

L2

rl

L3

19

25 -

-.

outer

.-.-

wall

17 9

--

outer _-_-.

w311 3o”

f’;

-si

M

SiO

oul_e,wall .-. 520

530

sio

sso

Fig. 2. Vessel wall: for presentation of numerical results the curved segment is straightened and the wall is unrolled, the sections of the segment are indicated.

K.

T/TP

-

SIO

s20 3

3-

2

2

2

2

1

1

1

1

1

1

1

3 ICI+ -1 0

ICI -1

1

3

0

3

c:--. -1 0

1

1

3 i/-1

0

iI+ -1

3

3

2

2

2

2

2

1

1

1

1

1

1

1

3 ELI+ -1 0

c:i -1

1

3

0

3

3

1

3 I""

-‘I

d

a -1

i

2

2

2

2

1

1

;

1 3 ilIl+

1

0

1 c?L -1

0

T/TP

1

-

CL -1

1

3

3

2

2

2

1

1

3 2

E.I+ -1 0

1

3 ICI+ -1 0

1

0

1

-1

0

IczL

1

3j'

a -1

1

0

2

2

1

1

1 D -1

1 450

2 3 c?L -1

0

3

1 1

1

00

2 0 -1

0

3

ic.JL -1 0

1

-450

1

2 Q -1

900

s50

5x0

2 ICI+ -1

450

s30

2

3

00

al.

0.22

Sl

-450

PERKTOLD~~

0

1

3 i.II -1 0

0 -1

1

O

’

0

1

900

0.40

Sl

SIO

s30

s20

s50

s40

1 31

31

3f

3

31

2-

2$

2 1 3 -4 k--\. 0

3

2I

2

2

2

1

1

‘If--+

lLz+‘-%

1 00

3

2-1 4

0

L -1

1

2

0

1

3 i, -1 ,

0

1

21

1 ILL--

1

3

450:

-1

0

1

0

1

3 I-,. -1 0 I

1

‘b -1 A

Fig. 3. Axial velocity profiles at maximum (T/I?? = 0.22) and minimum flow rate (T/W = 0.4) at different cross-sections and angle positions across the vessel diameter; cross-sections and angle positions are indicated. The inset shows the angle positions at the cross-section of flow entrance Sl. The profiles are presented in the direction of the small arrows of the inset.

-450 1

Wall shear stress distribution P5 SS ,0- y, 4*** ‘\ 4.. . . . a*\

0

*‘. . . . . . .‘) . .._..a-. . . _ . . . . .; ‘, , , . . . #.;, t-. . . *.__/ s3cl I -.*,‘--, :. . l . \

0

,.__“,;’ ,“_ -.,:: 1.. -*a 1 ‘(1 : .* * 1’,’ . ,. _ * * 1.. .-__‘,‘* ’

P45 s45 ‘\ *.:I ’ ‘,‘-v_,

.,,--__- N---_-_- * 1,,--__- 2:: ,:. I-_r 1i ,/’ 1 ‘;

.O T/TP

-

0.22

T/?P

-

0.40

’

u

-

u

-

2.0

1.0

T/TP

-

-

T/TP

-

0.22

-

0.40

u-l.0

u

c(

- I.0

H

Fig. 4. Vector field and profiles of secondary velocity at maximum and minimum flow rate: the vectors of secondary velocity (velocity and direction of flow) are symbolized as thin lines extending from the marked points.

structure of the stress pattern is almost the same throughout the main part of the period. A totally different structure is displayed only for the phase angle of minimum flow rate at T/TP= 0.4, where reverse flow mainly occurs near the wall. The hatched area indicates the zones of reversed flow. Figure 6 shows the distribution of the wall shear stress during the period at selected wall points. The alterations of the stress direction and of the stress values at the selected points are displayed. The diagram illustrates the angle segment in which the shear stress acts at the wall throughout the period. The shear stress vector corresponding to the maximum flow rate (T/TP

= 0.22) is denoted by the letter ‘c’, corresponding to the minimum flow rate (T/TP= 0.4) by the letter ‘e’. Figure 7 presents the contour plots of nondimensional wall shear stress values at maximum (A) and minimum flow rate (B), for time-averaged wall shear stress values: T/TP= 0.0-0.4 (C)and T/TP = O.bl.0, i.e. one cycle (D) as well as for steady flow under Reynolds number Re = 100. The plots indicate the sites of elevated stress at the two selected phase angles and the sites of elevated mean stress. The comparison of the contour plots for pulsatile, timeaveraged and steady flow demonstrates that the structure of the stress distribution is stable and that the wall

aqi

oss

0%

woq

El

El

-

y

lu!pualxa

ocs

OB-

OES

mkm

OES

MA.+w

E8

oi7s

oc =

EB

2-I

Zl

ozs

28

za

ozs

ozs

Cl

0

01s

LB

dL/l

ra

01s

- dl/l

ots

1-I ODE'0

(aw

ts

LS

LS

!d

!d

Jeaqs IleM ‘paw!pu!

OOti'O - dl/l

00s.

ssws

pasJaAaJJO sauoz aql saimpu!

aJe sJov24

‘Ip2rnlassah aql beau dog

saug u!qI se pazyoqmh

05s

-_ = z_ zw- 3 :

oss

ayl

qu!od

paqmu

ObS

Ggt=

OES

E”2.3tl ozs

022’0

01s

dl/l LS

‘5 53~

El

El

El -

01065

MIL.3ki

&a

&a

oils

E8

Zl

Zl

Zl

02s

oqs

za

28

ze

ti

000'0

~a

. .

. .

ie

i//p;

. . . ..-q-v. . . .

I:::: if..*.._. . ..-a_.

.._

.

ta

.

I:: f;i

.

IS

;_

_

.-SC

:

:-Lb

... . I z

. . . . . . . , . . . . . . *.... ,... . ,... .

. .-

s..... . . . . . .

- dl/l OLS

. .

LS

i-

\p$

,.:::: i

- dl/l O!S

L-l

O&l'0

cl

!d

6

id

. ____- ----- - _aCrX_r/CaaaCa.* *-..--- -_... . c.. : -. ..._..._..-----_ sz .-_C_._..,.Mr.C .C,., *-..----d > / ,,,,,,,,,,.,,,,...__- - .. .. .. .. .. .. . : o:::::.?.“” : ::::::::;~~‘ I2 <*r., ,~,,,~,‘,,“‘~-------I:‘:::. .__ _ _ _ _ p~$;:>::__-_____ _ _. _. . : r.....\\\....,.. \ . . . . . .._ ~~s~~~~~~,,~~~====:::::::~ . . ..\.\\.... _I_ !d

-

ie WJE paqsley

:a13h 3v!p~w aql Buynp SS9JlS Jeaqs l~ern~o play Joj3ah

MOMu~nru!u!w) p’o = d& se saseqd amy pawlas

669

Wall shear stress distribution

Pi

el % i 25

-

e-

e -

I

’

c

KC

e-

e

!

17

_._e&.

._.__.__~~~

-.-

TVI I

c e

&

c

I

c

c

I

&

9-p

e

!

e

SlO .

25

0 17

c

?%

c

e-

e /

c

I __._

-‘;

.._.-.-

__._.

eFT.-.-.-~/c

c

-.-.-

e p’

‘----r--

c

H

I

R SIO

s20

s40

s30

Rel’tu l/TP: -

!-

I

I Pi

.-.

9

I

I

c

Pi = PI0

:

.02

.06 .lO .14 .I8 .20 .22 .24 .26 .30 .34 .36 .40 .42 .46 .50 .54 .56 .62 .66 .70 .J4 .76 .82 .86 .90 .94 .96

c.:T/TP

=

.22

e:

20.00

-

T/TP = .40

Fig. 6. Vector field of wall shear stress; variation throughout the cardiac cycle. The wall shear stress vectors at different wall points are displayed for selected time phases T/TP as indicated. The wall shear stress vectors for maximum and minimum flow rate are described with ‘c’ and ‘e’, respectively.

zones of high and low stresses do not or only slightly shift when the flow rate alters. This is not true for reversed flow. The observation indicates that the structure of the wall shear stress distribution for forward directed flow is almost independent of the flow rate. As a consequence it results that wall shear stress structure is determined through the geometry of the flow field.

CONCLUSION

Wall shear stress is of considerable physiological importance, however, it is difficult to measure this quantity. Our theoretical model uses the incompressible Navier-Stokes equations for a Newtonian fluid. The equations are solved numerically with a newly developed finite element procedure. The timedependent wall shear stress is calculated from the finite element-velocity field near the wall. The application of computer yields detailed results for the spatial and temporal distribution of the wall shear stress. In this study the arterial wall has been assumed to be not distensible. As shown by Womersley (1957) for

flow in straight tubes the assumption of wall distensibility leads to a slight quantitative improvement of the calculated flow pattern without influence on the essential pattern of flow in the vessel. The investigation shows strong time-varying secondary motion because of the three curves. The secondary velocity interacts with the axial motion and produces highly skewed velocity profiles. Thus, the resulting wall shear stress pattern is rather complex. Theoretical studies of the type presented here should be appropriate for bioengineers and physiologists analysing unresolved questions in the mechanism through which fluid mechanics might act in the initiation and progression of arterial lesions. A comparison of the theoretical results presented here with wall shear stress measurements is not possible, because of the lack of detailed results for pulsatile

flow through

the carotid

siphon.

Acknowledgement-The authors wish to thank Professor Dr T. Kenner, Institute of Physiology, University of Graz, for his assistance concerning medical problems and questions. This investigation is supported by Austrian research fund (Fonds .zur Forderung der wissenschaftlichen Forschung in Osterreich), Project-Nr. P 5954 P, Vienna, Austria.

670

K. PERKTOLD et al.

. i

I

.-

a.

i4

h

.-

0-l

-S-

ul

._

a

G;

(u

s

.

.-

a

g

,: F

’ m

* ._

a

Wall shear stress distribution REFERENCES

Bharadvaj, B. K., Mabon, R. F. and Giddens, D. P. (1982) Steady ftow in a model of the human carotid bifurcation. Part I-Flow visualization. J. Biomechanics 15. 349-362. Caro, C. G., Fitz-Gerald, J. M. and Schroter, R. C. (1971) Atheroma and arterial wall shear: observation, correlation and proposal of a shear dependent mass transfer mechanism of atherogenesis. Proc. R. Sot. 177, 109-159. Chorin, A. J. (1968) Numerical solution of the Navier-Stokes equations. Math. Comput. 22, 741-762. Friedman, M. H. and Ehrlich, L. W. (1984) Numerical simulation of aortic bifurcation flows: the effect of flow divider curvature. .I. Biomechanics 17, 881-888. Fry, D. L. (1968) Acute vascular endothelial changes associated with increased blood velocity gradients. Circulation Res. 22, 165-197. Lutz, R. J., Cannon, J. N., Bischoff, K. B., Dedrick, R. L., Stiles, R. K. and Fry, D. L. (1977) Wall shear distribution in a model canine artery during steady flow. Circulation Res. 41, 391-399. Nerem, R. M. and Cornhill, J. F. (1980) The role of fluid mechanics in atherogenesis. .I. 6iomech. Engng 102, 181-189.

671

Pedley, T. J. (1980) The Fluid Mechanics of Large Blood Vessels, p. 51. Cambridge University Press, Cambridge. Pei, Z.-H., Xi, B.-S. and.Hwang, N.H.C. (1985) Wall shear stress distribution in a model human aortic arch: assessment by an electrochemiql technique. .I. Biomechanics 18, 645-656. Perktold, K., Gruber, K., Kenner, T. and Florian, H. (1984) Calculation of pulsatile flow and particle paths in an aneurysm-model. Basic Res. Cardiol. 79, 253-26 I. Perktold, K. (1985) Numerische tisung des instationiiren dreidimensionalen Nauier-Stokes Problems. Bericht der Math.-Stat. Sektion im Forschungszentrum Graz, Nr. 255. Graz. Perktold, K., Florian, H. gnd Hilbert, D. (1987) Analysis of pulsatile blood Row: a carotid siphon model. J. biomed. Engng 9,46-53. Schettler, G., Nerem, R. M.. Schmid-Schonbein, H., MGrl, H. and Diehm, H. (Eds.) (1983) Fluid Dynamics as a Localizing Factor for Atherosclerosis. Springe;, Berlin. Schmid-Schb;nbein. H. and Wurzineer. L. J. (19861 Transnort phenomena in pulsating post-st&tic vdrtex how in’ arteries. Nouv. Ren. Fr. Hematol. 28, 257-267. Womersley, J. R. (1957) An elastic tube theory of pulse transmission and oscillatory flow in mammalian arteries. WADC technical report TR-614.