- Email: [email protected]
The assignment of velocity profiles in finite element simulations of pulsatile flow in arteries
Comput.
Bid.
Med.
0 1997 Elsevier
Vol. 27, No. 3, pp. 233-247. 1997 Science Ltd. All rights resewed Printed in Cheat Britain 0010-4S25/97 $17.00+0.00
PII: Soo10425(5q@oooG1
THE ASSIGNMENT OF VELOCITY PROFILES IN FINITE ELEMENT SIMULATIONS OF PULSATILE FLOW IN ARTERIES A. REDAELLI*‘*‘, F. BOSCHETTI’** and F. INZOLI*~~ ’ Dipartimento di Bioingegneria, Politecnico di Milano, P.za L. da Vinci 32,20133 Milano, Italy * CeBITeC, Politecnico di Milan0 and Ospedale S. Raffaele, Milano, Italy ’ Dipartimento di Energetica, Politecnico di Milano, Milano, Italy (Received 3 1 May 1996; revised 20 November 1996) Abstract-In this paper we present a new method for the assignment of pulsatile velocity profiles as input boundary conditions in finite element models of arteries. The method is based on the implementation of the analytical solution for developed pulsatile flow in a rigid straight tube. The analytical solution provides the fluid dynamics of the region upstream from the fluid domain to be investigated by means of the finite element approach. In standard fluid dynamics finite element applications, the inlet developed velocity profiles are achieved assuming velocity boundary conditions to be easily implementable-such as flat or parabolic velocity profiles-applied to a straight tube of appropriate length. The tube is attached to the inflow section of the original fluid domain so that the flow can develop fully. The comparison between the analytical solution and the traditional numerical approach indicates that the analytical solution has some advantages over the numerical one. Moreover, the results suggest that subroutine employment allows a consistent reduction in solving time especially for complex fluid dynamic model, and significantly decreases the storage and memory requirements for computations. 8 1997 Elsevier Science Ltd. Computational fluid dynamics Pulsatile flow FORTRAN User’s subroutine FIDAP
Womersley
Finite element method
1. INTRODUCTION In recent years the computational approach-based on finite element,finite volume or finite difference methods-has been successfullyused to simulate blood flow in physiological and pathological models of arterial districts. Some of these studies were performed according to the simplifying, yet nonphysiological, assumption of steady flow conditions [l-7]. Indeed, the numerical solution of the complete Navier-Stokes equations,including the transientterm, is expensivein terms of computing time and requires high performance computers.It is only in recent years that they have gradually become available and that a number of studiesconcerning unsteady flow have been published [8-161. One of the main problems associatedwith the numerical simulation of pulsatile flow is the definition of the inflow boundary conditions, which can be seteither in terms of pressure gradients or velocity profiles. The latter approach is more widely used since it ensures greater computational stability. Moreover, the pressuregradient is often unknown. A conventional way to obtain pulsatile developed boundary conditions is to attach a straight tube to the inflow section of the original fluid domain and assignflat or parabolic velocity profiles to the tube as inlet boundary conditions. The length of the tube is specifically chosenso that the flow can develop well. The drawback of this approachis that it requires a buffer spaceand therefore the solution of a greater number of equations. In this paper, the analytical solution of the NavierStokes equation with pulsatile boundary conditions is implemented. The analytical solution was calculated by Hale et al. [17] for a homogeneous, incompressible, Newtonian fluid in the hypothesis of a rigid, straight, cylindrical tube. * Author to whom correspondence should be addressed. 233
234
al.
A.REDAELLIU
Hale’s solution is applied to the assignment of pulsatile velocity profiles as input boundary conditions in finite element models of arteries. The analytical solution simulates the fluid dynamics of the region upstream from the fluid domain under investigation. An automated procedure, requiring the setting of very few parameters-the properties of blood and inlet boundary geometry-provides the instantaneous velocity profile needed to perform the finite element simulation. This paper also provides a comparison between the analytical method and the traditional approach which requires the finite element simulation of the upstream region. This study has two main objectives: to easily calculate the correct inlet boundary conditions and to reduce the computational time involved. 2. MATERIAL
AND METHODS
2.1. Theory In the case of a cylindrical, rigid tube (of radius R) and homogeneous, incompressible, Newtonian fluid flow, the Navier-Stokes equation is expressed as:
au,--
(1)
at -
where u,=u,(l;t) is the axial velocity, r is the radial co-ordinate, p=p(t,z) is the pressure throughout the tube, p and p are the viscosity and the density of the Newtonian fluid, respectively. Consider a periodic pressure gradient expressed as a cosinusoidal Fourier series as follows: & - = - Real(pa( w)ej”), dz
(2)
where w is the frequency of the periodic phenomenon. The solution of the Navier-Stokes equation is expressed by:
uJr,t)=Real{
Tejti[
1 - ‘WI},
(3)
where Jo is the Bessel function of order zero and Wo = R(wrlv) “* is the Womersley number. Hence, volumetric flow rate is expressed by: 2 vqwo
J,(v=jWo) Jo(vqwo>
(4)
I)
where J, is the Bessel function of order 1. Using periodic signals, it is very useful to express equations (3) and (4) in the frequency domain as proposed by He et al. [ 181. Moreover, the following dimensionless forms are defined and used: &fQz-
q
?rR2ii
A
=-
44R
&wTv;
Pfi
where ii is the average velocity in the tube. The Fourier transforms (^) of the dimensionless velocity and flow rate are:
U,(r,t)= O(J2jej”=
1
(5)
Simulations
of pulsatile
flow in arteries
235
(6)
where: Re= 5’ P a= Woti (i=O,1,2,...) By combining equations (5) and (6), it is possible to express the non-dimensional velocity U, as a function of the dimensionless flow rate Q for each harmonic in the frequency domain.
(7)
It is therefore necessary to apply the inverse Fourier transform dimensional form for the output array.
and restore the
2.1.1. Algorithm description. A scheme of the algorithm is presented in Fig. 1. The algorithm is parametric with respect to the current geometry and the flow rate. Two files provide these information. The first contains the vessel radius and the co-ordinates of the centre of the inlet boundary surface. The second contains the discretized flow rate waveform. The waveform data are quadratically interpolated to achieve 64 values in the period. At the initial time step, an initialization algorithm (t=ts,a,) is run first. The Fast Fourier transform of the flow rate Q(t) is first calculated, according to Ferziger’s method [19]. The approximated formulae reported by Abramowitz and Stegun [20] for the calculation of Bessel functions of order zero, Jo(r,t), and of order one, J,(t) are used. When substituted in equation (7), they allow the computation of the axial velocity in the frequency domain. The inverse Fourier transform is then applied to determine the axial velocity as a function of time and spatial position. The time discretization obtained is for 64 samples per period. The solution vector is stored in an output database. When the successive time steps are performed the output database is used as an input database (algorithm t,,,
A. REDAELLIet al.
236
are stored at t,,,, is called VELO.DAT. This means that the subroutine performs the previous calculations only nnodes times: once the VELO.DAT file has been created the subroutine only reads the values of axial velocity and time from that file. For each time step, the FIDAP simulation is run once the boundary conditions have been set for all of the inlet nodal points. If the actual time is lower than fend, time is incremented. 3.2. Design testing The velocity profiles were computed by means of the developed software for two different test cases, using the aortic [21] and the femoral [12] waveforms as input flow rate
l
vessel
l
centre coordinates
radius:
l l l
Fourier transform of Q(r) Bessel function J&J) Bessel function Jo(t) Bessel function Jt(f)
Algorithm (t = f,‘,,,) :............................................. ............................................
I............... .. .... ...... . .. ... .. ...... . .... ..... ... ... ......... . .... Interpolate U&J) to achieve the value for the current time
:.......................... ....... ....,... .Algorithm ...... ....... . ..(t,,.,, ........et.....<.. .t&.....
Fig. I. Flow diagram of the algorithm for the calculation of Womersley developed profiles.
Simulations
of pulsatile
flow
in arteries
237
databases (Fig. 3). The aortic flow simulation was performed with a mean Reynolds number of 4160 which corresponds to a mean flow rate of 5 Vmin and a radius of 10 mm. The femoral flow simulation was performed with a mean Reynolds number of 420 which corresponds to a mean flow rate of 0.2 Cumin and a radius of 4 mm. Two different flow rate frequencies were considered: 60 and 90 bpm. The corresponding Womersley numbers were, respectively, 15 and 19, for the aortic flow simulation and 6 and 8 for the femoral artery. Parallel to this, the velocity profiles were calculated by means of FIDAP as the outlet velocity profiles of a rigid, straight, cylindrical tube. This, as previously stated, is the usual
VELO.DAT
No
END simulation Fig. 2. Flow diagram of FIDAP-W0MER.f subroutineinterface.
A. REDMLLI
238
etal.
(a)
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time itir&$,r,
-2 4
Fig. 3. Selected inputs: non-dimensional curves used to simulate (a) the aortic [21] and (b) the femoral [ 121 flow rates. They are non-dimensionalized with respect to (a) the maximum aortic flow rate (Q.3, and (b) the average femoral flow rate (Q.,), respectively. The two profiles are used as inputs for the testing of the simulation models.
procedure to simulate pulsatile, developed velocity profiles: add a straight inlet tube to a model which simulates pulsatile flow in a vascular district, the flow develops and tends asymptotically to the analytical solution after an appropriate number of diameters. Two axi-symmetric finite element models of a straight vessel (aorta and femoral artery) were created using FIDAP, with a vessel length of 40 diameters in both cases. In order to evaluate the number of diameters needed to obtain a fairly well developed flow, three crosssections were also defined through each vessel at a distance of 10, 20 and 30 diameters, respectively, from the inlet. As far as the inlet boundary conditions of the straight vessel are concerned, most finite element codes, including FIDAP, allow the implementation of polynomial formulations. In the present study, the velocity boundary conditions at the inlet section were assumed to be either flat or parabolic, typical of undeveloped and developed steady flows, respectively. Indeed, the proper number of diameters is expected to depend on the inlet velocity profile of the straight tube. Figure 4 gives a qualitative description of the two approaches-analytical and numerical-taken in the present study. Quadrilateral nine-node isoparametric annular elements were employed for the mesh of the finite element models. The total number of elements was 1777 for a total of 2000 nodes. The mesh discretization was 16 in the radial direction, and 30 in the axis of symmetry direction; the grid size was finer close to the wall and towards the inlet section. Further assumptions for the models were: homogeneous and Newtonian fluid (p= 1060 kg/m3, ,z= 3 X 10 -3 kglms), impermeable walls, no slip boundary conditions at
flat and pambolic verOcify pn#los used ot the inletrecth of the addkdjluid domains
boundmy inlet veloci@ conditian
Fig. 4. Different approaches utilized to calculate the inlet boundary velocity profiles for FIDAP fluid domains.
Simulations of pulsatile flow in arteries
239
Table 1. Models and simulations
Aorta $(ZD axi-symmetric) Femoral artery §(2D axi-symmetric) Femoral artery *
Use of W0MER.f
Use of FIDAP with flat input velocity profile
Use of FIDAP with parabolic input velocity profile
x X X
X X -
X X -
Q Frequencies: 60 and 90 bpm. * Frequency of 60 bpm.
the wall, uniform pressure at the outlet section, neglected gravitational effects. To verify the efficacy of the subroutine for 3D models, it was also applied to assessing the velocity profiles in a 3D straight model with femoral artery flow rate conditions and a frequency of 60 bpm. The 3D mesh was mapped; the mesh elements were nine-node isoparametric elements (bricks) and the total number of elements was 18777. The model assumptions and simulation procedure were the same as described for the axi-symmetric models. The full Navier-Stokes equations were solved using Gale&in’s weighted residual approach in conjunction with finite element approximation. A quasi-Newton solver was employed for the solution method. The time integration technique adopted was the implicit backward Euler with a fixed time step equal to l/64 of the period. It was necessary to perform three cycles to reach flow and pressure stability. The computations were carried out on a HP 735. The simulations performed are summarized in Table 1. 3.3. Testing of the results Figures 5.1 and 5.2 refer to the aortic flow simulations (f=60 bpm). They compare the inlet velocity time courses provided by the W0MER.f subroutine with the correspondent velocity time courses obtained by FIDAP from the straight vessel simulation. The velocity time courses refer to the centre-line node (r=O; panels a), to an intermediate node (r=0.5R; panels b) and to a node close to the wall (rk0.9R; panels c). The curves obtained with FIDAP were calculated at a distance from the inlet section of 10, 20, 30 and 40 diameters (#>, respectively, assuming a parabolic (Fig. 5.1) and a flat (Fig. 5.2) inlet profile. In both cases the axial velocities calculated by FIDAP differ from section to section (10 4 to 40 (6), indicating that the flow was still developing towards the pulsatile steady state condition represented in the figures by the analytical solution. The difficulties in achieving the pulsatile steady state condition experienced by the numerical code were due to the high inertial convective acceleration in the central region (panels a, Fig. 5.1 and 5.2). In fact a closer correspondence between the velocity time courses is shown at those points in the proximity of the wall characterized by lower inertia (panel b and c of Fig. 5.1 and 5.2). The same calculations were performed for the femoral artery with a frequency of 60 bpm. In this case too, both a parabolic (Fig. 5.3) and a flat (Fig. 5.4) inlet profile were taken into account for the FIDAP simulations. Unlike the aortic model, the low inertia associated with the femoral flow allows the full development of the velocity for FIDAP simulations after a few diameters. Results from aorta and femoral artery flow simulations, with a frequency of 90 bpm, are given in Fig. 6.1 and 6.2. Only the data at 40 diameters were taken into consideration, since this length corresponded to the velocity profiles closest to the developed ones. The nodes for which the velocity time courses are shown are the same as in Fig. 5.1-5.4. The results are similar to those at 60 bpm and, again, the femoral artery model shows a closer correspondence between the analytical and the FIDAP solutions. The results of each test case are summarized in Table 2, which reports the differences between the analytical values and the FIDAP values at the four tube length considered. Table 2 refers to the three previously defined locations (~0, r=0.5R, r=0,9R). The differences are weighted with respect to the maximum values calculated analytically for
240
A. REDAELLI
et
al.
sorta(/=6obpm) pamtdic
aorta (/I
pwjile
104
v on%)
pwfile
104
. . . . . . . x4 30 4
233 I
0.00 1 0.0
60 bpm)
jlaf
0.2
0.4
0.6 time (I)
0.8
. . . . . . . 34 -3Q4 -a14
1.0
0.2
0.0
0.4
1)
0.6 h
0.8
1.0
(9
1)
v Ws)
0.0 0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
tims w
-0.80 1
C)
5.1. femoral ytrty pambdic
femoral artery (i-
(f= 60 bpm)
j&Jr
profile
60 bpm)
pwJilc
v (Ids)
0.80 i
0.20 0.10 0.00 -0.10
b)
.O
104
. . .._._
204
Simulations of pulsatile flow in arteries
241
that node in the cycle. The mean values (-c SD) over the cycle of the weighted differences are reported for each case. The instantaneous velocity profiles of the aortic and of the femoral arteries with a frequency of 90 bpm are shown in Fig. 7 at five time steps. Again, the W0MER.f results are compared with those obtained with either flat or parabolic inlet velocity profiles from FIDAP. The 3D velocity profiles at the inlet surface of the 3D femoral artery model (frequency=60 bpm) are reported in Fig. 8. As previously stated, this simulation was performed in order to verify the subroutine ability to handle nodes in the 3D domain. 4. DISCUSSION The results presented in the previous section show a good agreement between W0MER.f and FIDAP velocity calculations where an appropriate length for the vessel in FIDAP simulations is adopted. It should, however, be noted that after 40 diameters, the FIDAP simulations give slightly different values than those computed by the W0MER.f. The first reason for this may be that 40 diameters are not sufficient to reach stable solutions. For steady flows, it is possible to define a characteristic length (Lc) equal approximately to 0.058DRe (where D is the diameter and Re is the Reynolds number) at which a centreline velocity is within 1% of its final value [22]. For the vessels under consideration, Lc is equal to 241 diameters for the aorta model and to 24 diameters for the femoral one. Even if Lc refers to steady flows, a dependence on Re is also to be expected for pulsatile flow. Indeed, the centre-line velocity time course of the aorta under both the parabolic and flat inlet profile assumptions (Fig. 5.la and 5.2a), shows great variability from one section to another. The second reason is that the instantaneous flow rates for FIDAP and W0MER.f simulations are not exactly the same. This is evident from Fig. 7b, where in the second instant considered, the velocity profile calculated by the subroutine overestimates the corresponding FIDAP profiles all over the vessel cross-section. The opposite then happens in the third and fourth instants. These behaviours are detectable also in Fig. 6.2. For instance, the fourth instant of Fig. 7b almost corresponds to the negative peak of Fig. 6.2, where the velocities calculated by means of the W0MER.f subroutine are lower than those calculated with FIDAP (Fig. 6.2a, b, c). The differences in flow rates may be due to the following factors: 1. In finite element codes, derived quantities such as flow rate, energy losses and so on, depend on interpolation algorithms adopted across the elements and may induce slight but appreciable differences with respect to the exact solution; 2. The Fast Fourier Transform-direct and inverse-used in the W0MER.f subroutine may introduce approximations into the shape of the inflow curve (direct FFT) and into that of the outflow velocity (inverse FlT). Furthermore, some general features can be pointed out, as follows: the parabolic profile is less likely to reach the Womersley profile than the flat one, especially at high flow rates, as shown in Fig. 5. la and 5.3a; the centre-line velocity is more dependent on the inlet velocity profile in the central region of the vessel due to the higher inertial contents of the flow. This occurrence is particularly evident in Fig. 7a, where the FIDAP velocity profile overestimates or underestimates the W0MER.f profiles in the central region of the vessel, depending on Fig. 5. Comparison of the velocitytimecourses calculated with W0MER.f andwith FIDAP.The cycle frequencyis 60bpm.The four FIDAP curvesare.calculatedat a distancefrom the inlet sectionof 10.20.30and40diameters, respectively. Panels (a),(b) and(c) referto thevelocitytime courses calculated: (a)at thesymmetryaxisof thelumen(r=O); (b) midwaybetweenthe axis and thevesselwall(r=OX); (c) closeto thewall(r=0.9i?).5.1.Aotic modeloutput:FIDAPparabolic inlet profile.5.2.Aortic modeloutput:FIDAP flat inletprofile.5.3.Femoral modeloutput:FIDAP parabolicinletprofile.5.4.Femoralmodeloutput:FIDAP flat inletprofile.
242
A.
REDAELLI
et al.
sons (I= 90 bpm) femoral artery if= 90 bpm) v Ws)
0.0
0.2
0.4
0.6 -
w
0.8
1.0 M
0.20 0.10 0.M) 0.10 -0.20
Cl
Fig. 6. Comparison of the velocity time courses calculated with W0MER.f and with FIDAP. Both flat and parabolic inlet velocity profiles are considered for FIDAP simulations. The cycle frequency is 90 bpm. (a) r=O; (b) r=O.SR; (c) r=0.9R. 6.1. (left) Aortic model output. 6.2. (right) Femoral model output.
whether the inlet velocity profile is parabolic or flat. The considerations below are related to the issue of the computational time. If a length of 40 diameters is used, the 2D FIDAP straight vessel model takes approximately 330 s of CPU time for each cycle to achieve the developed pulsatile velocity profile. This time does not depend on whether the inlet profile is parabolic or flat. The CPU time, in turn, is approximately 470 s when W0MER.f is linked to the the same 2D FIDAP straight model. The difference between 470 s and 330 s gives the time taken by the subroutine W0MER.f to calculate the boundary conditions (140 s). Hence, the use of W0MER.f means approximately 190 s gained (58%) over the “straight vessel” methodology. If, however, a length of only 20 diameters is used for the “straight vessel” methodology (low Reynolds numbers such as in the femoral case), the application of the W0MER.f subroutine apparently gives no benefit. Nevertheless, this result obviously depends on the morphology of the original domain and on the degree of discretization of the finite element mesh. Indeed, in the case of complex fluid domains, such as 3D ones, a consistent time saving is to be expected if the proposed algorithm is employed. In fact, the solution time depends almost quadratically on the number of elements in the mesh, which significantly increases if an inlet straight vessel of convenient length is added to the original fluid domain. As far as the algorithm (&, < f < tend)(Fig. 1) is concerned, the complexity is fin), where n is the number of nodes lying on the inlet boundary surface. This also means a significant reduction in the computational storage and memory requirements.
Calculated differences
(%)
- 25.62 14.0 8.2*5.5 - 25.6k 12.5 9.3*9.3 - 10.4+ 10.7 0.5*3.1 -5.5+4.1 0.323.1
a-F10
(%)
- 20.028.6 5.3k5.8 - 15.4+ 12.1 6.5k9.5 - 0.8k3.2 -0.13.2 - 1.2k3.0 -0.1+3.1
a-F20
r=O.O
(o/o)
- 10.8k6.6 4.1 k5.9 - 16.7a9.7 5.2k9.8 - 0.723.2 -0.lk3.2 -0.7Et3.1 -0.lk3.1
a-F30
(%)
- 14.31t5.5 2.9k6.2 - 8.2+: 13.9 5.Ok9.8 - 0.6i3.2 -0.lk3.2 -0.6k3.1 -0.1 k3.1
a-F40
(%)
- 8.7k8.8 1.3k6.0 - 8.529.7 l.lk9.9 - 1.2k3.6 -0.li3.3 -0.8k3.2 -0.123.1
a-F10
(%)
- 2.12 1.4 0.3k6.4 -4.O+ 11.2 -0.2k10.2 -0.523.3 -0.lk3.3 -0.5k3.2 -0.1r3.2
a-F20
r=OS
(%)
- 1.6k2.1 O.Ok6.4 - 2.82 10.5 -0.4k10.3 - 0.5k3.3 -0.lk3.3 -0.5~3.2 -0.lk3.2
a-F30
(%)
- 1.8+ 1.9 -0.li6.5 - 1.9+ 10.5 -0.4k10.3 - 0.553.3 -0.1+3.3 -0.5i3.2 -0.123.2
a-F40
(%)
r=0.9).
calculated
(%)
Weighted
o.ozt4.9
2.221.3 -0.928.0 -0.31t13.3 - 2.Ozt 13.2 -0.224.6 - 0.0*4.7 0.2zt4.8
a-F40
in the cycle.
o.ort4.9
(%)
0.0+4.9
a-F30
1.6+ 1.9 - l.Ik8.1 1.5~tl3.8 - 2.02 13.2 - 0.2k4.7 0.0~~4.8 - 0.1 i4.8
(%)
r=0.9
5.2e4.7 - 1.4k8.1 1.3* 12.3 - 2.3* 13.2 - 0.1 t4.7 O.Ok4.8 0.024.8
a-F20
analytically
7.8k7.0 - 2.328.2 6.6+ 16.2 -3.4k13.2 3.1 i5.4 - 0.2~4.9 1.1 k4.8 -0.124.9
a-F10
differences (mean + SD) between the analytical solution (a) and the FIDAP simulations, weighted with respect to the maximum values are reported, for each test case, for the four considered tube lengths (FIO, F20, F30, F40) at three nodes of the inlet section (r=O, r=0.5,
Aorta (f=60) parabolic Aorta (f=60) flat Aorta (f=90) parabolic Aorta (f=90) flat Femoral (f= 60) parabolic Femoral (f=60) flat Femoral (f=90) parabolic Femoral (f=90) flat
Test cases
Table 2.
2 3. rl
2 5 g. 2 s 5’
#’
244
A. RED-
et al.
5. CONCLUSIONS
t=O.ws _-‘E -_ -I I-_ ,,--
This paper presents a computer-generated analytical solution of developed pulsatile flow. The analytical solution is used to assign the input velocity boundary conditions to finite element models of arteries since it provides the fluid dynamics of the region upstream from the investigated fluid domain to be investigated by means of the finite element approach. Results indicate that the analytical solution is potentially advantageous with respect to aorta (f= 90bpm) t = 0.2%
t =0.13s
1.00
0.00
mls
t = 0.38s
t =OSls
-
flat
__ - . -
parabolic
femoral artery (f= 90 bpm)
t=O.oos
t = 0.2%
t =0.13s
I
I
0.10
0.00 Ids
t = 0.38s
-
t = 0.51s
flat
. _ _ . _ parabolic -analytic
Fig. 7. Comparison of 2D axi-symmetric velocity profiles obtained with the W0MER.f subroutine and with FIDAP. with flat and parabolic inlet velocity profiles: (a) aorta; (b) femoral artery. The re-sults refer to five instants of the cycle, with a frequency of 90 bpm.
Simulations of pulsatile flow in arteries
245
the numerical one. Its application is much simpler and more immediate: the traditional method provides results which are affected by both the length of the attached tube and the inflow velocity profile. The method presented here, on the other hand, needs only the vessel radius, the co-ordinates of the centre of the inlet section and the inlet flow rate curve as inputs, thanks to its parametric feature. The moderate effort required to implement this method, together with the reduction in computing time it involves, indicates that the algorithm presented in this paper is a powerful tool in the optimal solution of fluid-dynamic problems assiciated with cardiovascular circulation. The subroutine is available upon request by E-mail. 6. SUMMARY In this paper we present a new method for the assignment of pulsatile velocity profiles as input boundary conditions in finite element models of arteries. The method is based on the implementation of the analytical solution of developed pulsatile flow for a rigid straight tube to the finite element approach. The analytical solution provides the fluid dynamics of the region upstream from the fluid domain to be investigated by means of the finite element approach. In usual fluid dynamic finite element applications, the developed velocity profiles are achieved by assuming easily implementable velocity boundary conditions-such as flat or parabolic velocity profiles-applied to a straight tube of appropriate length. The tube is attached to the inflow section of the original fluid domain, so that the flow can develop fully. The proposed method is expected to improve the computational approach both in terms of flexibility and reduction of computing time. As a test case, the algorithm developed is implemented as a subroutine to a finite element
Fig. 8. 3D velocity profiles at the inlet surface of a 3D model of a straight tube with the femoral artery flow rate (frequency =60 bpm). The four surfaces on the top are observed from an elevation angle of 40” with respect to the flow axis, the two surface on the bottom are. observed from an elevation angle of - 40”. The mesh is mapped. The inlet boundary surfaces are discretized with 121 nodes. Their co-ordinates are passed to W0MER.f. For each node, the subroutine calculates the corresponding velocity time c&se on the basis of its distance from the centre of the inlet boundarv surface. These data are stored in the VELO.DAT output file. At each time steo, FIDAP reads the proper velocity values from the output file and assign‘s the velocity boundary conditions as reported in the picture.
246
A.REDAELLI
etal.
commercial code named FIDAP. The results are compared with those obtained with the same finite element package using the traditional approach mentioned above. Two situations, involving the aortic and femoral flow rates in straight tubes of appropriate dimensions, were considered. Results indicate that the analytical solution is potentially advantageous with respect to the numerical one. Its application is simple and immediate. Thanks to its parametric feature, only the vessel radius, the co-ordinates of the centre of the inlet section and the inlet flow rate curve are needed to run FIDAP simulations. It is therefore possible to assign pulsatile inlet velocity boundary conditions to different fluid domains quite easily. Moreover, the results suggest that subroutine employment allows a consistent reduction of the solution time, especially for complex fluid dynamic models, and brings about a significant reduction in the computational storage and memory requirements. In fact, the execution time depends almost quadratically on the number of elements in the model, which significantly increases if a straight inlet vessel of convenient length is added to the original fluid domain. The algorithm presented has a complexity e
Simulations of pulsatile flow in arteries About the Author-&mwA Bosctmrn graduated in Electronic Engineering in 1990 from the Politecnico of Milan, and received a PhD in Biomedical Engineering from the MURST in 1994. She is currently a researcher at CeBITeC at the San Raffaele Hospital in Milan. About the Author-Fmo INZOLI graduated in Mechanical Engineering in 1984 from the Politecnico of Milan, and received a PhD in Energetics from the MURST in 1993. He is currently a researcher at the Department of Energetics of the Politecnico of Milan.
247
Bid.
Med.
0 1997 Elsevier
Vol. 27, No. 3, pp. 233-247. 1997 Science Ltd. All rights resewed Printed in Cheat Britain 0010-4S25/97 $17.00+0.00
PII: Soo10425(5q@oooG1
THE ASSIGNMENT OF VELOCITY PROFILES IN FINITE ELEMENT SIMULATIONS OF PULSATILE FLOW IN ARTERIES A. REDAELLI*‘*‘, F. BOSCHETTI’** and F. INZOLI*~~ ’ Dipartimento di Bioingegneria, Politecnico di Milano, P.za L. da Vinci 32,20133 Milano, Italy * CeBITeC, Politecnico di Milan0 and Ospedale S. Raffaele, Milano, Italy ’ Dipartimento di Energetica, Politecnico di Milano, Milano, Italy (Received 3 1 May 1996; revised 20 November 1996) Abstract-In this paper we present a new method for the assignment of pulsatile velocity profiles as input boundary conditions in finite element models of arteries. The method is based on the implementation of the analytical solution for developed pulsatile flow in a rigid straight tube. The analytical solution provides the fluid dynamics of the region upstream from the fluid domain to be investigated by means of the finite element approach. In standard fluid dynamics finite element applications, the inlet developed velocity profiles are achieved assuming velocity boundary conditions to be easily implementable-such as flat or parabolic velocity profiles-applied to a straight tube of appropriate length. The tube is attached to the inflow section of the original fluid domain so that the flow can develop fully. The comparison between the analytical solution and the traditional numerical approach indicates that the analytical solution has some advantages over the numerical one. Moreover, the results suggest that subroutine employment allows a consistent reduction in solving time especially for complex fluid dynamic model, and significantly decreases the storage and memory requirements for computations. 8 1997 Elsevier Science Ltd. Computational fluid dynamics Pulsatile flow FORTRAN User’s subroutine FIDAP
Womersley
Finite element method
1. INTRODUCTION In recent years the computational approach-based on finite element,finite volume or finite difference methods-has been successfullyused to simulate blood flow in physiological and pathological models of arterial districts. Some of these studies were performed according to the simplifying, yet nonphysiological, assumption of steady flow conditions [l-7]. Indeed, the numerical solution of the complete Navier-Stokes equations,including the transientterm, is expensivein terms of computing time and requires high performance computers.It is only in recent years that they have gradually become available and that a number of studiesconcerning unsteady flow have been published [8-161. One of the main problems associatedwith the numerical simulation of pulsatile flow is the definition of the inflow boundary conditions, which can be seteither in terms of pressure gradients or velocity profiles. The latter approach is more widely used since it ensures greater computational stability. Moreover, the pressuregradient is often unknown. A conventional way to obtain pulsatile developed boundary conditions is to attach a straight tube to the inflow section of the original fluid domain and assignflat or parabolic velocity profiles to the tube as inlet boundary conditions. The length of the tube is specifically chosenso that the flow can develop well. The drawback of this approachis that it requires a buffer spaceand therefore the solution of a greater number of equations. In this paper, the analytical solution of the NavierStokes equation with pulsatile boundary conditions is implemented. The analytical solution was calculated by Hale et al. [17] for a homogeneous, incompressible, Newtonian fluid in the hypothesis of a rigid, straight, cylindrical tube. * Author to whom correspondence should be addressed. 233
234
al.
A.REDAELLIU
Hale’s solution is applied to the assignment of pulsatile velocity profiles as input boundary conditions in finite element models of arteries. The analytical solution simulates the fluid dynamics of the region upstream from the fluid domain under investigation. An automated procedure, requiring the setting of very few parameters-the properties of blood and inlet boundary geometry-provides the instantaneous velocity profile needed to perform the finite element simulation. This paper also provides a comparison between the analytical method and the traditional approach which requires the finite element simulation of the upstream region. This study has two main objectives: to easily calculate the correct inlet boundary conditions and to reduce the computational time involved. 2. MATERIAL
AND METHODS
2.1. Theory In the case of a cylindrical, rigid tube (of radius R) and homogeneous, incompressible, Newtonian fluid flow, the Navier-Stokes equation is expressed as:
au,--
(1)
at -
where u,=u,(l;t) is the axial velocity, r is the radial co-ordinate, p=p(t,z) is the pressure throughout the tube, p and p are the viscosity and the density of the Newtonian fluid, respectively. Consider a periodic pressure gradient expressed as a cosinusoidal Fourier series as follows: & - = - Real(pa( w)ej”), dz
(2)
where w is the frequency of the periodic phenomenon. The solution of the Navier-Stokes equation is expressed by:
uJr,t)=Real{
Tejti[
1 - ‘WI},
(3)
where Jo is the Bessel function of order zero and Wo = R(wrlv) “* is the Womersley number. Hence, volumetric flow rate is expressed by: 2 vqwo
J,(v=jWo) Jo(vqwo>
(4)
I)
where J, is the Bessel function of order 1. Using periodic signals, it is very useful to express equations (3) and (4) in the frequency domain as proposed by He et al. [ 181. Moreover, the following dimensionless forms are defined and used: &fQz-
q
?rR2ii
A
=-
44R
&wTv;
Pfi
where ii is the average velocity in the tube. The Fourier transforms (^) of the dimensionless velocity and flow rate are:
U,(r,t)= O(J2jej”=
1
(5)
Simulations
of pulsatile
flow in arteries
235
(6)
where: Re= 5’ P a= Woti (i=O,1,2,...) By combining equations (5) and (6), it is possible to express the non-dimensional velocity U, as a function of the dimensionless flow rate Q for each harmonic in the frequency domain.
(7)
It is therefore necessary to apply the inverse Fourier transform dimensional form for the output array.
and restore the
2.1.1. Algorithm description. A scheme of the algorithm is presented in Fig. 1. The algorithm is parametric with respect to the current geometry and the flow rate. Two files provide these information. The first contains the vessel radius and the co-ordinates of the centre of the inlet boundary surface. The second contains the discretized flow rate waveform. The waveform data are quadratically interpolated to achieve 64 values in the period. At the initial time step, an initialization algorithm (t=ts,a,) is run first. The Fast Fourier transform of the flow rate Q(t) is first calculated, according to Ferziger’s method [19]. The approximated formulae reported by Abramowitz and Stegun [20] for the calculation of Bessel functions of order zero, Jo(r,t), and of order one, J,(t) are used. When substituted in equation (7), they allow the computation of the axial velocity in the frequency domain. The inverse Fourier transform is then applied to determine the axial velocity as a function of time and spatial position. The time discretization obtained is for 64 samples per period. The solution vector is stored in an output database. When the successive time steps are performed the output database is used as an input database (algorithm t,,,
A. REDAELLIet al.
236
are stored at t,,,, is called VELO.DAT. This means that the subroutine performs the previous calculations only nnodes times: once the VELO.DAT file has been created the subroutine only reads the values of axial velocity and time from that file. For each time step, the FIDAP simulation is run once the boundary conditions have been set for all of the inlet nodal points. If the actual time is lower than fend, time is incremented. 3.2. Design testing The velocity profiles were computed by means of the developed software for two different test cases, using the aortic [21] and the femoral [12] waveforms as input flow rate
l
vessel
l
centre coordinates
radius:
l l l
Fourier transform of Q(r) Bessel function J&J) Bessel function Jo(t) Bessel function Jt(f)
Algorithm (t = f,‘,,,) :............................................. ............................................
I............... .. .... ...... . .. ... .. ...... . .... ..... ... ... ......... . .... Interpolate U&J) to achieve the value for the current time
:.......................... ....... ....,... .Algorithm ...... ....... . ..(t,,.,, ........et.....<.. .t&.....
Fig. I. Flow diagram of the algorithm for the calculation of Womersley developed profiles.
Simulations
of pulsatile
flow
in arteries
237
databases (Fig. 3). The aortic flow simulation was performed with a mean Reynolds number of 4160 which corresponds to a mean flow rate of 5 Vmin and a radius of 10 mm. The femoral flow simulation was performed with a mean Reynolds number of 420 which corresponds to a mean flow rate of 0.2 Cumin and a radius of 4 mm. Two different flow rate frequencies were considered: 60 and 90 bpm. The corresponding Womersley numbers were, respectively, 15 and 19, for the aortic flow simulation and 6 and 8 for the femoral artery. Parallel to this, the velocity profiles were calculated by means of FIDAP as the outlet velocity profiles of a rigid, straight, cylindrical tube. This, as previously stated, is the usual
VELO.DAT
No
END simulation Fig. 2. Flow diagram of FIDAP-W0MER.f subroutineinterface.
A. REDMLLI
238
etal.
(a)
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 time itir&$,r,
-2 4
Fig. 3. Selected inputs: non-dimensional curves used to simulate (a) the aortic [21] and (b) the femoral [ 121 flow rates. They are non-dimensionalized with respect to (a) the maximum aortic flow rate (Q.3, and (b) the average femoral flow rate (Q.,), respectively. The two profiles are used as inputs for the testing of the simulation models.
procedure to simulate pulsatile, developed velocity profiles: add a straight inlet tube to a model which simulates pulsatile flow in a vascular district, the flow develops and tends asymptotically to the analytical solution after an appropriate number of diameters. Two axi-symmetric finite element models of a straight vessel (aorta and femoral artery) were created using FIDAP, with a vessel length of 40 diameters in both cases. In order to evaluate the number of diameters needed to obtain a fairly well developed flow, three crosssections were also defined through each vessel at a distance of 10, 20 and 30 diameters, respectively, from the inlet. As far as the inlet boundary conditions of the straight vessel are concerned, most finite element codes, including FIDAP, allow the implementation of polynomial formulations. In the present study, the velocity boundary conditions at the inlet section were assumed to be either flat or parabolic, typical of undeveloped and developed steady flows, respectively. Indeed, the proper number of diameters is expected to depend on the inlet velocity profile of the straight tube. Figure 4 gives a qualitative description of the two approaches-analytical and numerical-taken in the present study. Quadrilateral nine-node isoparametric annular elements were employed for the mesh of the finite element models. The total number of elements was 1777 for a total of 2000 nodes. The mesh discretization was 16 in the radial direction, and 30 in the axis of symmetry direction; the grid size was finer close to the wall and towards the inlet section. Further assumptions for the models were: homogeneous and Newtonian fluid (p= 1060 kg/m3, ,z= 3 X 10 -3 kglms), impermeable walls, no slip boundary conditions at
flat and pambolic verOcify pn#los used ot the inletrecth of the addkdjluid domains
boundmy inlet veloci@ conditian
Fig. 4. Different approaches utilized to calculate the inlet boundary velocity profiles for FIDAP fluid domains.
Simulations of pulsatile flow in arteries
239
Table 1. Models and simulations
Aorta $(ZD axi-symmetric) Femoral artery §(2D axi-symmetric) Femoral artery *
Use of W0MER.f
Use of FIDAP with flat input velocity profile
Use of FIDAP with parabolic input velocity profile
x X X
X X -
X X -
Q Frequencies: 60 and 90 bpm. * Frequency of 60 bpm.
the wall, uniform pressure at the outlet section, neglected gravitational effects. To verify the efficacy of the subroutine for 3D models, it was also applied to assessing the velocity profiles in a 3D straight model with femoral artery flow rate conditions and a frequency of 60 bpm. The 3D mesh was mapped; the mesh elements were nine-node isoparametric elements (bricks) and the total number of elements was 18777. The model assumptions and simulation procedure were the same as described for the axi-symmetric models. The full Navier-Stokes equations were solved using Gale&in’s weighted residual approach in conjunction with finite element approximation. A quasi-Newton solver was employed for the solution method. The time integration technique adopted was the implicit backward Euler with a fixed time step equal to l/64 of the period. It was necessary to perform three cycles to reach flow and pressure stability. The computations were carried out on a HP 735. The simulations performed are summarized in Table 1. 3.3. Testing of the results Figures 5.1 and 5.2 refer to the aortic flow simulations (f=60 bpm). They compare the inlet velocity time courses provided by the W0MER.f subroutine with the correspondent velocity time courses obtained by FIDAP from the straight vessel simulation. The velocity time courses refer to the centre-line node (r=O; panels a), to an intermediate node (r=0.5R; panels b) and to a node close to the wall (rk0.9R; panels c). The curves obtained with FIDAP were calculated at a distance from the inlet section of 10, 20, 30 and 40 diameters (#>, respectively, assuming a parabolic (Fig. 5.1) and a flat (Fig. 5.2) inlet profile. In both cases the axial velocities calculated by FIDAP differ from section to section (10 4 to 40 (6), indicating that the flow was still developing towards the pulsatile steady state condition represented in the figures by the analytical solution. The difficulties in achieving the pulsatile steady state condition experienced by the numerical code were due to the high inertial convective acceleration in the central region (panels a, Fig. 5.1 and 5.2). In fact a closer correspondence between the velocity time courses is shown at those points in the proximity of the wall characterized by lower inertia (panel b and c of Fig. 5.1 and 5.2). The same calculations were performed for the femoral artery with a frequency of 60 bpm. In this case too, both a parabolic (Fig. 5.3) and a flat (Fig. 5.4) inlet profile were taken into account for the FIDAP simulations. Unlike the aortic model, the low inertia associated with the femoral flow allows the full development of the velocity for FIDAP simulations after a few diameters. Results from aorta and femoral artery flow simulations, with a frequency of 90 bpm, are given in Fig. 6.1 and 6.2. Only the data at 40 diameters were taken into consideration, since this length corresponded to the velocity profiles closest to the developed ones. The nodes for which the velocity time courses are shown are the same as in Fig. 5.1-5.4. The results are similar to those at 60 bpm and, again, the femoral artery model shows a closer correspondence between the analytical and the FIDAP solutions. The results of each test case are summarized in Table 2, which reports the differences between the analytical values and the FIDAP values at the four tube length considered. Table 2 refers to the three previously defined locations (~0, r=0.5R, r=0,9R). The differences are weighted with respect to the maximum values calculated analytically for
240
A. REDAELLI
et
al.
sorta(/=6obpm) pamtdic
aorta (/I
pwjile
104
v on%)
pwfile
104
. . . . . . . x4 30 4
233 I
0.00 1 0.0
60 bpm)
jlaf
0.2
0.4
0.6 time (I)
0.8
. . . . . . . 34 -3Q4 -a14
1.0
0.2
0.0
0.4
1)
0.6 h
0.8
1.0
(9
1)
v Ws)
0.0 0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
tims w
-0.80 1
C)
5.1. femoral ytrty pambdic
femoral artery (i-
(f= 60 bpm)
j&Jr
profile
60 bpm)
pwJilc
v (Ids)
0.80 i
0.20 0.10 0.00 -0.10
b)
.O
104
. . .._._
204
Simulations of pulsatile flow in arteries
241
that node in the cycle. The mean values (-c SD) over the cycle of the weighted differences are reported for each case. The instantaneous velocity profiles of the aortic and of the femoral arteries with a frequency of 90 bpm are shown in Fig. 7 at five time steps. Again, the W0MER.f results are compared with those obtained with either flat or parabolic inlet velocity profiles from FIDAP. The 3D velocity profiles at the inlet surface of the 3D femoral artery model (frequency=60 bpm) are reported in Fig. 8. As previously stated, this simulation was performed in order to verify the subroutine ability to handle nodes in the 3D domain. 4. DISCUSSION The results presented in the previous section show a good agreement between W0MER.f and FIDAP velocity calculations where an appropriate length for the vessel in FIDAP simulations is adopted. It should, however, be noted that after 40 diameters, the FIDAP simulations give slightly different values than those computed by the W0MER.f. The first reason for this may be that 40 diameters are not sufficient to reach stable solutions. For steady flows, it is possible to define a characteristic length (Lc) equal approximately to 0.058DRe (where D is the diameter and Re is the Reynolds number) at which a centreline velocity is within 1% of its final value [22]. For the vessels under consideration, Lc is equal to 241 diameters for the aorta model and to 24 diameters for the femoral one. Even if Lc refers to steady flows, a dependence on Re is also to be expected for pulsatile flow. Indeed, the centre-line velocity time course of the aorta under both the parabolic and flat inlet profile assumptions (Fig. 5.la and 5.2a), shows great variability from one section to another. The second reason is that the instantaneous flow rates for FIDAP and W0MER.f simulations are not exactly the same. This is evident from Fig. 7b, where in the second instant considered, the velocity profile calculated by the subroutine overestimates the corresponding FIDAP profiles all over the vessel cross-section. The opposite then happens in the third and fourth instants. These behaviours are detectable also in Fig. 6.2. For instance, the fourth instant of Fig. 7b almost corresponds to the negative peak of Fig. 6.2, where the velocities calculated by means of the W0MER.f subroutine are lower than those calculated with FIDAP (Fig. 6.2a, b, c). The differences in flow rates may be due to the following factors: 1. In finite element codes, derived quantities such as flow rate, energy losses and so on, depend on interpolation algorithms adopted across the elements and may induce slight but appreciable differences with respect to the exact solution; 2. The Fast Fourier Transform-direct and inverse-used in the W0MER.f subroutine may introduce approximations into the shape of the inflow curve (direct FFT) and into that of the outflow velocity (inverse FlT). Furthermore, some general features can be pointed out, as follows: the parabolic profile is less likely to reach the Womersley profile than the flat one, especially at high flow rates, as shown in Fig. 5. la and 5.3a; the centre-line velocity is more dependent on the inlet velocity profile in the central region of the vessel due to the higher inertial contents of the flow. This occurrence is particularly evident in Fig. 7a, where the FIDAP velocity profile overestimates or underestimates the W0MER.f profiles in the central region of the vessel, depending on Fig. 5. Comparison of the velocitytimecourses calculated with W0MER.f andwith FIDAP.The cycle frequencyis 60bpm.The four FIDAP curvesare.calculatedat a distancefrom the inlet sectionof 10.20.30and40diameters, respectively. Panels (a),(b) and(c) referto thevelocitytime courses calculated: (a)at thesymmetryaxisof thelumen(r=O); (b) midwaybetweenthe axis and thevesselwall(r=OX); (c) closeto thewall(r=0.9i?).5.1.Aotic modeloutput:FIDAPparabolic inlet profile.5.2.Aortic modeloutput:FIDAP flat inletprofile.5.3.Femoral modeloutput:FIDAP parabolicinletprofile.5.4.Femoralmodeloutput:FIDAP flat inletprofile.
242
A.
REDAELLI
et al.
sons (I= 90 bpm) femoral artery if= 90 bpm) v Ws)
0.0
0.2
0.4
0.6 -
w
0.8
1.0 M
0.20 0.10 0.M) 0.10 -0.20
Cl
Fig. 6. Comparison of the velocity time courses calculated with W0MER.f and with FIDAP. Both flat and parabolic inlet velocity profiles are considered for FIDAP simulations. The cycle frequency is 90 bpm. (a) r=O; (b) r=O.SR; (c) r=0.9R. 6.1. (left) Aortic model output. 6.2. (right) Femoral model output.
whether the inlet velocity profile is parabolic or flat. The considerations below are related to the issue of the computational time. If a length of 40 diameters is used, the 2D FIDAP straight vessel model takes approximately 330 s of CPU time for each cycle to achieve the developed pulsatile velocity profile. This time does not depend on whether the inlet profile is parabolic or flat. The CPU time, in turn, is approximately 470 s when W0MER.f is linked to the the same 2D FIDAP straight model. The difference between 470 s and 330 s gives the time taken by the subroutine W0MER.f to calculate the boundary conditions (140 s). Hence, the use of W0MER.f means approximately 190 s gained (58%) over the “straight vessel” methodology. If, however, a length of only 20 diameters is used for the “straight vessel” methodology (low Reynolds numbers such as in the femoral case), the application of the W0MER.f subroutine apparently gives no benefit. Nevertheless, this result obviously depends on the morphology of the original domain and on the degree of discretization of the finite element mesh. Indeed, in the case of complex fluid domains, such as 3D ones, a consistent time saving is to be expected if the proposed algorithm is employed. In fact, the solution time depends almost quadratically on the number of elements in the mesh, which significantly increases if an inlet straight vessel of convenient length is added to the original fluid domain. As far as the algorithm (&, < f < tend)(Fig. 1) is concerned, the complexity is fin), where n is the number of nodes lying on the inlet boundary surface. This also means a significant reduction in the computational storage and memory requirements.
Calculated differences
(%)
- 25.62 14.0 8.2*5.5 - 25.6k 12.5 9.3*9.3 - 10.4+ 10.7 0.5*3.1 -5.5+4.1 0.323.1
a-F10
(%)
- 20.028.6 5.3k5.8 - 15.4+ 12.1 6.5k9.5 - 0.8k3.2 -0.13.2 - 1.2k3.0 -0.1+3.1
a-F20
r=O.O
(o/o)
- 10.8k6.6 4.1 k5.9 - 16.7a9.7 5.2k9.8 - 0.723.2 -0.lk3.2 -0.7Et3.1 -0.lk3.1
a-F30
(%)
- 14.31t5.5 2.9k6.2 - 8.2+: 13.9 5.Ok9.8 - 0.6i3.2 -0.lk3.2 -0.6k3.1 -0.1 k3.1
a-F40
(%)
- 8.7k8.8 1.3k6.0 - 8.529.7 l.lk9.9 - 1.2k3.6 -0.li3.3 -0.8k3.2 -0.123.1
a-F10
(%)
- 2.12 1.4 0.3k6.4 -4.O+ 11.2 -0.2k10.2 -0.523.3 -0.lk3.3 -0.5k3.2 -0.1r3.2
a-F20
r=OS
(%)
- 1.6k2.1 O.Ok6.4 - 2.82 10.5 -0.4k10.3 - 0.5k3.3 -0.lk3.3 -0.5~3.2 -0.lk3.2
a-F30
(%)
- 1.8+ 1.9 -0.li6.5 - 1.9+ 10.5 -0.4k10.3 - 0.553.3 -0.1+3.3 -0.5i3.2 -0.123.2
a-F40
(%)
r=0.9).
calculated
(%)
Weighted
o.ozt4.9
2.221.3 -0.928.0 -0.31t13.3 - 2.Ozt 13.2 -0.224.6 - 0.0*4.7 0.2zt4.8
a-F40
in the cycle.
o.ort4.9
(%)
0.0+4.9
a-F30
1.6+ 1.9 - l.Ik8.1 1.5~tl3.8 - 2.02 13.2 - 0.2k4.7 0.0~~4.8 - 0.1 i4.8
(%)
r=0.9
5.2e4.7 - 1.4k8.1 1.3* 12.3 - 2.3* 13.2 - 0.1 t4.7 O.Ok4.8 0.024.8
a-F20
analytically
7.8k7.0 - 2.328.2 6.6+ 16.2 -3.4k13.2 3.1 i5.4 - 0.2~4.9 1.1 k4.8 -0.124.9
a-F10
differences (mean + SD) between the analytical solution (a) and the FIDAP simulations, weighted with respect to the maximum values are reported, for each test case, for the four considered tube lengths (FIO, F20, F30, F40) at three nodes of the inlet section (r=O, r=0.5,
Aorta (f=60) parabolic Aorta (f=60) flat Aorta (f=90) parabolic Aorta (f=90) flat Femoral (f= 60) parabolic Femoral (f=60) flat Femoral (f=90) parabolic Femoral (f=90) flat
Test cases
Table 2.
2 3. rl
2 5 g. 2 s 5’
#’
244
A. RED-
et al.
5. CONCLUSIONS
t=O.ws _-‘E -_ -I I-_ ,,--
This paper presents a computer-generated analytical solution of developed pulsatile flow. The analytical solution is used to assign the input velocity boundary conditions to finite element models of arteries since it provides the fluid dynamics of the region upstream from the investigated fluid domain to be investigated by means of the finite element approach. Results indicate that the analytical solution is potentially advantageous with respect to aorta (f= 90bpm) t = 0.2%
t =0.13s
1.00
0.00
mls
t = 0.38s
t =OSls
-
flat
__ - . -
parabolic
femoral artery (f= 90 bpm)
t=O.oos
t = 0.2%
t =0.13s
I
I
0.10
0.00 Ids
t = 0.38s
-
t = 0.51s
flat
. _ _ . _ parabolic -analytic
Fig. 7. Comparison of 2D axi-symmetric velocity profiles obtained with the W0MER.f subroutine and with FIDAP. with flat and parabolic inlet velocity profiles: (a) aorta; (b) femoral artery. The re-sults refer to five instants of the cycle, with a frequency of 90 bpm.
Simulations of pulsatile flow in arteries
245
the numerical one. Its application is much simpler and more immediate: the traditional method provides results which are affected by both the length of the attached tube and the inflow velocity profile. The method presented here, on the other hand, needs only the vessel radius, the co-ordinates of the centre of the inlet section and the inlet flow rate curve as inputs, thanks to its parametric feature. The moderate effort required to implement this method, together with the reduction in computing time it involves, indicates that the algorithm presented in this paper is a powerful tool in the optimal solution of fluid-dynamic problems assiciated with cardiovascular circulation. The subroutine is available upon request by E-mail. 6. SUMMARY In this paper we present a new method for the assignment of pulsatile velocity profiles as input boundary conditions in finite element models of arteries. The method is based on the implementation of the analytical solution of developed pulsatile flow for a rigid straight tube to the finite element approach. The analytical solution provides the fluid dynamics of the region upstream from the fluid domain to be investigated by means of the finite element approach. In usual fluid dynamic finite element applications, the developed velocity profiles are achieved by assuming easily implementable velocity boundary conditions-such as flat or parabolic velocity profiles-applied to a straight tube of appropriate length. The tube is attached to the inflow section of the original fluid domain, so that the flow can develop fully. The proposed method is expected to improve the computational approach both in terms of flexibility and reduction of computing time. As a test case, the algorithm developed is implemented as a subroutine to a finite element
Fig. 8. 3D velocity profiles at the inlet surface of a 3D model of a straight tube with the femoral artery flow rate (frequency =60 bpm). The four surfaces on the top are observed from an elevation angle of 40” with respect to the flow axis, the two surface on the bottom are. observed from an elevation angle of - 40”. The mesh is mapped. The inlet boundary surfaces are discretized with 121 nodes. Their co-ordinates are passed to W0MER.f. For each node, the subroutine calculates the corresponding velocity time c&se on the basis of its distance from the centre of the inlet boundarv surface. These data are stored in the VELO.DAT output file. At each time steo, FIDAP reads the proper velocity values from the output file and assign‘s the velocity boundary conditions as reported in the picture.
246
A.REDAELLI
etal.
commercial code named FIDAP. The results are compared with those obtained with the same finite element package using the traditional approach mentioned above. Two situations, involving the aortic and femoral flow rates in straight tubes of appropriate dimensions, were considered. Results indicate that the analytical solution is potentially advantageous with respect to the numerical one. Its application is simple and immediate. Thanks to its parametric feature, only the vessel radius, the co-ordinates of the centre of the inlet section and the inlet flow rate curve are needed to run FIDAP simulations. It is therefore possible to assign pulsatile inlet velocity boundary conditions to different fluid domains quite easily. Moreover, the results suggest that subroutine employment allows a consistent reduction of the solution time, especially for complex fluid dynamic models, and brings about a significant reduction in the computational storage and memory requirements. In fact, the execution time depends almost quadratically on the number of elements in the model, which significantly increases if a straight inlet vessel of convenient length is added to the original fluid domain. The algorithm presented has a complexity e
Simulations of pulsatile flow in arteries About the Author-&mwA Bosctmrn graduated in Electronic Engineering in 1990 from the Politecnico of Milan, and received a PhD in Biomedical Engineering from the MURST in 1994. She is currently a researcher at CeBITeC at the San Raffaele Hospital in Milan. About the Author-Fmo INZOLI graduated in Mechanical Engineering in 1984 from the Politecnico of Milan, and received a PhD in Energetics from the MURST in 1993. He is currently a researcher at the Department of Energetics of the Politecnico of Milan.
247