A. Dutta D. M. Wang .1. M. Tarbell Department of Chemical Engineering, The Pennsylvania State University, University Park, PA 16802
Numerical Analysis of Flow in an Elastic Artery Model Oscillatory and pulsatile flows of Newtonian fluids in straight elastic tubes are simulated numerically with the aid of Ling and Atabek's "local flow" assumption for the nonlinear convective acceleration terms. For the first time, a theoretical assessment of the local flow assumption is presented, and the range of validity of the assumption is estimated by comparison with perturbation solutions of the complete flow problem. Subsequent simulations with the local flow model indicate that the flow field and associated wall shear stress are extremely sensitive to the phase angle between oscillatory pressure and flow waves (impedance phase angle). This phase angle, which is a measure of the wave reflection present in the system, is known to be altered by arterial disease (e.g., hypertension) and vasoactive drugs. Thus, the paper elucidates a mechanism by which subtle changes in systemic hemodynamics (i.e., phase angles) can markedly influence local wall shear stress values.
1
Introduction Blood vessel walls are elastic, and the diameter of large arteries can vary by ± 5 percent over the cardiac cycle [1]. Traditional interest in elastic vessels has focused on the problem of propagation of pressure and flow pulses in the cardiovascular system [2, 3]. Less attention has been paid to the influence of wall motion on the local flow field at a particular axial position in an artery [4], and the potential influence of wave propagation/reflection on local flow fields seems to have been clearly recognized only recently [5], although Womersley offered a few hints in his very last paper [6]. The general problem of determining local flow fields in elastic tubes is complicated by the necessity of specifying upstream and downstream boundary conditions. This is a great limitation in cardiovascular flow modeling because of the complex architecture of the branching vascular network. To get around this restriction, Ling and Atabek [4] introduced an assumption allowing the axial convective acceleration term, which involves the axial velocity gradient, to be expressed in terms of the local axial velocity. With this assumption, the problem can be solved locally without axial boundary conditions. Ling and Atabek [4] did not provide any real justification for their assumption or any assessment of its validity. In Section 3 of this paper we present a theoretical analysis of their assumption and delineate its range of validity. Later, in Section 5, we compare predictions based on the assumption with analytical perturbation solutions of the full equations and boundary conditions and verify its validity over a broad range of physiological flow conditions. We then use Ling and Atabek's assumption to consider the influence of wave reflection on local flow fields in a series of sinusoidal flow simulations in Section 6. We pay particular attention to the effects on wall shear stress. Finally, we present
multiharmonic physiological flow simulations under conditions characteristic of the proximal aorta in Section 7. 2
Formulation
T o determine velocity fields in arteries, we model blood flow using a homogeneous incompressible Newtonian fluid in an isotropic, thin-walled tube with longitudinal constraint when the fluid is subjected to an oscillatory pressure gradient. In addition, we assume that the tube is straight and axisymmetric. This, of course, restricts our results to straight sections of artery not too near a bend or bifurcation. T h e equations of motion and continuity can be simplified by making the long-wavelength approximation, wR_
«U
c
where to is the angular frequency, R is the radius, and c represents the wave speed. Under this assumption, the axial viscous transport terms are negligible and the radial equation of motion simply reduces to dp/dr = 0, indicating that the pressure is independent of radial position [7]. Hence, the governing equations can be expressed as:
dw dw dw — +u — + w — = dt dr dz Id, r dr
1 dp(z,t) dw +v dr2 P dz ^
dw dz
1 dw r dr
„
(2.2)
In addition, a wall property constitutive equation is required to complete the formulation. For a purely elastic model, which is quite accurate for arteries [8], the relation can be expressed as: R = R(P) The boundary conditions in the radial direction are
Contributed by the Bioengineering Division for publication in the JOURNAL OF BIOMECHANICAL ENGINEERINO. Manuscript received by the Bioengineering Division August 6, 1991; revised manuscript received September 30, 1991.
(2.1)
^ =0 u=0 dr
atr = 0
26/Vol. 114, FEBRUARY 1992
(symmetry)
(2.3)
(2.4)
Transactions of the ASME
Copyright © 1992 by ASME Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 10/07/2013 Terms of Use: http://asme.org/terms
w = 0 u-
M ' dt
at r = R(t,z)
(2.5)
(no-slip)
Axial boundary conditions will not be required with Ling and Atabek's assumption. To overcome the difficulties stemming from the moving boundary, a coordinate transformation, £ = r/(R(t,z)), is introduced, and Ling and Atabek's assumption, dw/dz = f(z, t) I w I, is invoked to eliminate the axial velocity gradient terms. The governing equations now become dt
p dz + R2
ar
flow is large or a is low enough that the harmonic profiles have different shapes, then the estimation of dw/dz will have some error. In Section 5.1 we will carry out a calculation based on Ling and Atabek's assumption and compare it with a more general perturbation solution [9]. 4
The Local Flow Method The problem of flow in an elastic artery was formulated in Section 2 (Eqs.' (2.6)-(2.9)). Neglecting the natural taper of an artery, dR/dz can be expressed as
£ d£ / f dR
u\dw
w/du
/dR\ _ /dR\ Kdz) ~ [dp)
u\
/dp\ \dz)
(4.1)
The derivative dR/dp represents the elastic response of the artery and experimental values are available [14]. Substituting Eq. (4.1) into Eq. (2.7) and applying the boundary conditions on u (Eqs. (2.8) and (2.9)), we finally obtain The unknown function, f(z, t), can be determined with the help of the transformed boundary conditions: dw
w
-0
w=0
u=0
at£ = 0
(2.8)
u=— dt
at £ ; 1
(2.9)
\ W
dR dp , = dp dz
[
I 2
r
£MY/£
£\w\di,
£MY/£-
\w\dt,
\
1 dR Jo dt
£lwltf£ (4.2) £lwlc?£
•>o
3
Assessment of Ling and Atabek's Assumption In order to provide a simple assessment of Ling and Atabek's assumption, we consider a restricted version of it, dw/dz = f(z, t)w, which is identical to the original form when the axial velocity does not change direction within the cross-section. We can integrate this simple equation to observe that, w = F(z, t)g(%, t), which indicates that the dependence of w on z and £ is separable. Under what conditions do we expect this to be valid? Let us consider a general situation in which the axial velocity can be represented by a Fourier series: S
Wn(z,H)e"
(3.1)
If the dependence of w„ on z and £ is separable, then w„(z,£)=F„U)g„(£),
(3.2)
and it follows that
f(z,t) ••
dz _ntiN '
dz N
(3.3)
w
2]
F„(z)gn{$)e"*
This is true only when g„(£) = knG(£), where G(£) is an arbitrary function of £ and k„ is a constant. The above argument indicates that all of the harmonics (including the steady part) should have a similar velocity profile, G(£). For oscillatory viscous flow in a tube, according to the linearized solution [6], the steady velocity profile is parabolic while the unsteady velocity profiles are close to flat when the unsteadiness parameter, a, is large (a = R0(o)/v)1/2, where a> is the fundamental angular frequency, R0 represents the mean radius and v denotes the kinematic viscosity). Therefore, when the mean flow is small compared to the oscillatory flow and a is large, Ling and Atabek's assumption is expected to be a good approximation. Also, when a is very low, all harmonics will have a parabolic profile and the assumption again should provide a good approximation. But, if the mean Journal of Biomechanical Engineering
When w is uniformly positive or negative within the crosssection at a particular instant, the expression in parentheses vanishes. Equations (2.6) and (4.2) now describe the simplified flow problem which does not require axial boundary conditions. However, to complete the specification of the problem, it is necessary to prescribe the local tube wall motion, R(t), and pressure gradient, dp/dz(t). In this section, we concentrate on simple oscillatory functions containing a single harmonic. dp
-
-r- (t)=K+kcos(ut) dz R(,t)=R[l+kRcos(wt-)]
(4.3) (4.4)
The mean parameters, K and R, amplitude parameters, k and kR, and the phase angle (hereafter referred to as the "radius phase angle") must be specified along with the frequency, co. As we shall see, the flow rate, Q(t), and the wall shear stress (WSS), r{t), which are computed as part of the solution of the flow problem, are well approximated by simple oscillatory functions of the form Q(.t) =
Q[l+kQCOs(ut-8)]
T(t)=T[l+krC0S(03t-P)]
(4.5) (4.6)
with very small contributions from the second and higher harmonics. The amplitude of second harmonics were always less than 8 percent of the first harmonic for all the cases we have simulated. In our simulations the radius phase angle, , was varied from 0 to 180 deg and the flow phase angle 6 was found to be relatively insensitive to this variation. The phase angle between radius and flow is given by 6 - . Assuming, as we have, that the wall of an artery is purely elastic, then the radius and pressure waveforms are in phase. In this case, the phase angles between radius and flow and between pressure and flow are equivalent. Now, physiological aortic input impedance data [1,11] are typically characterized by pressure-flow phase angles in the range 0 to - 60 deg for the first four harmonics of the waveforms. As we shall see, the phase angles between pressure FEBRUARY 1992, Vol. 114/27
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 10/07/2013 Terms of Use: http://asme.org/terms
and flow determined in our simulations, embrace the range of normal physiological input impedance data, and thus the results should be relevant to physiological flow analysis. Equations (2.6) and (4.2-4.4) were solved numerically using the IMSL subroutine DPDES which is based on a finite difference method employing the Method of Lines. The choice of the incremental step sizes m m
«
a. c
< r > - -10
5.1 Validity of Ling and Atabek's Assumption. Wang and Tarbell [9] have recently solved the nonlinear flow problem in a finite elastic tube driven by a sinusoidal pressure or flow input using perturbation methods. Their approach requires no assumption for dw/dz, but is limited asymptotically to flows in which the mean is small compared to the oscillatory component. We will first compare the predictions of Wang and Tarbell's analysis with the results of the local flow model for the case of zero mean flow where the perturbation method is most rigorous. Since the mean flow rate is an output of the local flow model, it was necessary to adjust the mean pressure
E
-204T , , I - 9 0 - 7 0 - 5 0 - 3 0 - 1 0 10 30 50 70 90 P-Q (degrees) Fig. 1(6) a = 18, k = 40.5 dynes/cm 3 , fl = 1.4 cm, O = 0, Osc. flow rate = 285 cm3/s Fig. 1 Mean pressure gradient versus P-Q phase angle. The solid lines represent the result of perturbation analysis [9] and the stars represent those of the local flow model.
28 / Vol. 114, FEBRUARY 1992 Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 10/07/2013 Terms of Use: http://asme.org/terms
Transactions of the AS ME
2 1;
Mean Press. Grad., Dynes/Cm 3
a = 12 fl = 1.9 cm, k = 18 dynes/cm 3
Fig. 3(a)
-200
-100
0
100
200
300
400
Mean flow rgte, c m * * 3 / s Fig. 2(a) a = 12, k = 18 dynes/cm 3 , Osc. flow rate
115cm 3 /s
50 40
9
Numerical Analysis of Flow in an Elastic Artery Model Oscillatory and pulsatile flows of Newtonian fluids in straight elastic tubes are simulated numerically with the aid of Ling and Atabek's "local flow" assumption for the nonlinear convective acceleration terms. For the first time, a theoretical assessment of the local flow assumption is presented, and the range of validity of the assumption is estimated by comparison with perturbation solutions of the complete flow problem. Subsequent simulations with the local flow model indicate that the flow field and associated wall shear stress are extremely sensitive to the phase angle between oscillatory pressure and flow waves (impedance phase angle). This phase angle, which is a measure of the wave reflection present in the system, is known to be altered by arterial disease (e.g., hypertension) and vasoactive drugs. Thus, the paper elucidates a mechanism by which subtle changes in systemic hemodynamics (i.e., phase angles) can markedly influence local wall shear stress values.
1
Introduction Blood vessel walls are elastic, and the diameter of large arteries can vary by ± 5 percent over the cardiac cycle [1]. Traditional interest in elastic vessels has focused on the problem of propagation of pressure and flow pulses in the cardiovascular system [2, 3]. Less attention has been paid to the influence of wall motion on the local flow field at a particular axial position in an artery [4], and the potential influence of wave propagation/reflection on local flow fields seems to have been clearly recognized only recently [5], although Womersley offered a few hints in his very last paper [6]. The general problem of determining local flow fields in elastic tubes is complicated by the necessity of specifying upstream and downstream boundary conditions. This is a great limitation in cardiovascular flow modeling because of the complex architecture of the branching vascular network. To get around this restriction, Ling and Atabek [4] introduced an assumption allowing the axial convective acceleration term, which involves the axial velocity gradient, to be expressed in terms of the local axial velocity. With this assumption, the problem can be solved locally without axial boundary conditions. Ling and Atabek [4] did not provide any real justification for their assumption or any assessment of its validity. In Section 3 of this paper we present a theoretical analysis of their assumption and delineate its range of validity. Later, in Section 5, we compare predictions based on the assumption with analytical perturbation solutions of the full equations and boundary conditions and verify its validity over a broad range of physiological flow conditions. We then use Ling and Atabek's assumption to consider the influence of wave reflection on local flow fields in a series of sinusoidal flow simulations in Section 6. We pay particular attention to the effects on wall shear stress. Finally, we present
multiharmonic physiological flow simulations under conditions characteristic of the proximal aorta in Section 7. 2
Formulation
T o determine velocity fields in arteries, we model blood flow using a homogeneous incompressible Newtonian fluid in an isotropic, thin-walled tube with longitudinal constraint when the fluid is subjected to an oscillatory pressure gradient. In addition, we assume that the tube is straight and axisymmetric. This, of course, restricts our results to straight sections of artery not too near a bend or bifurcation. T h e equations of motion and continuity can be simplified by making the long-wavelength approximation, wR_
«U
c
where to is the angular frequency, R is the radius, and c represents the wave speed. Under this assumption, the axial viscous transport terms are negligible and the radial equation of motion simply reduces to dp/dr = 0, indicating that the pressure is independent of radial position [7]. Hence, the governing equations can be expressed as:
dw dw dw — +u — + w — = dt dr dz Id, r dr
1 dp(z,t) dw +v dr2 P dz ^
dw dz
1 dw r dr
„
(2.2)
In addition, a wall property constitutive equation is required to complete the formulation. For a purely elastic model, which is quite accurate for arteries [8], the relation can be expressed as: R = R(P) The boundary conditions in the radial direction are
Contributed by the Bioengineering Division for publication in the JOURNAL OF BIOMECHANICAL ENGINEERINO. Manuscript received by the Bioengineering Division August 6, 1991; revised manuscript received September 30, 1991.
(2.1)
^ =0 u=0 dr
atr = 0
26/Vol. 114, FEBRUARY 1992
(symmetry)
(2.3)
(2.4)
Transactions of the ASME
Copyright © 1992 by ASME Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 10/07/2013 Terms of Use: http://asme.org/terms
w = 0 u-
M ' dt
at r = R(t,z)
(2.5)
(no-slip)
Axial boundary conditions will not be required with Ling and Atabek's assumption. To overcome the difficulties stemming from the moving boundary, a coordinate transformation, £ = r/(R(t,z)), is introduced, and Ling and Atabek's assumption, dw/dz = f(z, t) I w I, is invoked to eliminate the axial velocity gradient terms. The governing equations now become dt
p dz + R2
ar
flow is large or a is low enough that the harmonic profiles have different shapes, then the estimation of dw/dz will have some error. In Section 5.1 we will carry out a calculation based on Ling and Atabek's assumption and compare it with a more general perturbation solution [9]. 4
The Local Flow Method The problem of flow in an elastic artery was formulated in Section 2 (Eqs.' (2.6)-(2.9)). Neglecting the natural taper of an artery, dR/dz can be expressed as
£ d£ / f dR
u\dw
w/du
/dR\ _ /dR\ Kdz) ~ [dp)
u\
/dp\ \dz)
(4.1)
The derivative dR/dp represents the elastic response of the artery and experimental values are available [14]. Substituting Eq. (4.1) into Eq. (2.7) and applying the boundary conditions on u (Eqs. (2.8) and (2.9)), we finally obtain The unknown function, f(z, t), can be determined with the help of the transformed boundary conditions: dw
w
-0
w=0
u=0
at£ = 0
(2.8)
u=— dt
at £ ; 1
(2.9)
\ W
dR dp , = dp dz
[
I 2
r
£MY/£
£\w\di,
£MY/£-
\w\dt,
\
1 dR Jo dt
£lwltf£ (4.2) £lwlc?£
•>o
3
Assessment of Ling and Atabek's Assumption In order to provide a simple assessment of Ling and Atabek's assumption, we consider a restricted version of it, dw/dz = f(z, t)w, which is identical to the original form when the axial velocity does not change direction within the cross-section. We can integrate this simple equation to observe that, w = F(z, t)g(%, t), which indicates that the dependence of w on z and £ is separable. Under what conditions do we expect this to be valid? Let us consider a general situation in which the axial velocity can be represented by a Fourier series: S
Wn(z,H)e"
(3.1)
If the dependence of w„ on z and £ is separable, then w„(z,£)=F„U)g„(£),
(3.2)
and it follows that
f(z,t) ••
dz _ntiN '
dz N
(3.3)
w
2]
F„(z)gn{$)e"*
This is true only when g„(£) = knG(£), where G(£) is an arbitrary function of £ and k„ is a constant. The above argument indicates that all of the harmonics (including the steady part) should have a similar velocity profile, G(£). For oscillatory viscous flow in a tube, according to the linearized solution [6], the steady velocity profile is parabolic while the unsteady velocity profiles are close to flat when the unsteadiness parameter, a, is large (a = R0(o)/v)1/2, where a> is the fundamental angular frequency, R0 represents the mean radius and v denotes the kinematic viscosity). Therefore, when the mean flow is small compared to the oscillatory flow and a is large, Ling and Atabek's assumption is expected to be a good approximation. Also, when a is very low, all harmonics will have a parabolic profile and the assumption again should provide a good approximation. But, if the mean Journal of Biomechanical Engineering
When w is uniformly positive or negative within the crosssection at a particular instant, the expression in parentheses vanishes. Equations (2.6) and (4.2) now describe the simplified flow problem which does not require axial boundary conditions. However, to complete the specification of the problem, it is necessary to prescribe the local tube wall motion, R(t), and pressure gradient, dp/dz(t). In this section, we concentrate on simple oscillatory functions containing a single harmonic. dp
-
-r- (t)=K+kcos(ut) dz R(,t)=R[l+kRcos(wt-)]
(4.3) (4.4)
The mean parameters, K and R, amplitude parameters, k and kR, and the phase angle (hereafter referred to as the "radius phase angle") must be specified along with the frequency, co. As we shall see, the flow rate, Q(t), and the wall shear stress (WSS), r{t), which are computed as part of the solution of the flow problem, are well approximated by simple oscillatory functions of the form Q(.t) =
Q[l+kQCOs(ut-8)]
T(t)=T[l+krC0S(03t-P)]
(4.5) (4.6)
with very small contributions from the second and higher harmonics. The amplitude of second harmonics were always less than 8 percent of the first harmonic for all the cases we have simulated. In our simulations the radius phase angle, , was varied from 0 to 180 deg and the flow phase angle 6 was found to be relatively insensitive to this variation. The phase angle between radius and flow is given by 6 - . Assuming, as we have, that the wall of an artery is purely elastic, then the radius and pressure waveforms are in phase. In this case, the phase angles between radius and flow and between pressure and flow are equivalent. Now, physiological aortic input impedance data [1,11] are typically characterized by pressure-flow phase angles in the range 0 to - 60 deg for the first four harmonics of the waveforms. As we shall see, the phase angles between pressure FEBRUARY 1992, Vol. 114/27
Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 10/07/2013 Terms of Use: http://asme.org/terms
and flow determined in our simulations, embrace the range of normal physiological input impedance data, and thus the results should be relevant to physiological flow analysis. Equations (2.6) and (4.2-4.4) were solved numerically using the IMSL subroutine DPDES which is based on a finite difference method employing the Method of Lines. The choice of the incremental step sizes m m
«
a. c
< r > - -10
5.1 Validity of Ling and Atabek's Assumption. Wang and Tarbell [9] have recently solved the nonlinear flow problem in a finite elastic tube driven by a sinusoidal pressure or flow input using perturbation methods. Their approach requires no assumption for dw/dz, but is limited asymptotically to flows in which the mean is small compared to the oscillatory component. We will first compare the predictions of Wang and Tarbell's analysis with the results of the local flow model for the case of zero mean flow where the perturbation method is most rigorous. Since the mean flow rate is an output of the local flow model, it was necessary to adjust the mean pressure
E
-204T , , I - 9 0 - 7 0 - 5 0 - 3 0 - 1 0 10 30 50 70 90 P-Q (degrees) Fig. 1(6) a = 18, k = 40.5 dynes/cm 3 , fl = 1.4 cm, O = 0, Osc. flow rate = 285 cm3/s Fig. 1 Mean pressure gradient versus P-Q phase angle. The solid lines represent the result of perturbation analysis [9] and the stars represent those of the local flow model.
28 / Vol. 114, FEBRUARY 1992 Downloaded From: http://biomechanical.asmedigitalcollection.asme.org/ on 10/07/2013 Terms of Use: http://asme.org/terms
Transactions of the AS ME
2 1;
Mean Press. Grad., Dynes/Cm 3
a = 12 fl = 1.9 cm, k = 18 dynes/cm 3
Fig. 3(a)
-200
-100
0
100
200
300
400
Mean flow rgte, c m * * 3 / s Fig. 2(a) a = 12, k = 18 dynes/cm 3 , Osc. flow rate
115cm 3 /s
50 40
9
