Cardiovascular Fluid Mechanics

exercises and computational models

8W090

F.N. van de Vosse (2014)

Eindhoven University of Technology

Department of Biomedical Engineering
Preface 1

Preface
In this manuscript exercises and computational models related to the theory that is
treated in the lecture notes are given. The code MATLAB will be used to compute,
display and analyze the computational results of these models.
Chapter 1

General introduction

1.1 Sinusoidal waves

The periodic signals (f (t) = f (t + T )) like pressure and flow that appear in cardio-
vascular fluid mechanics can be expressed as the sum of a series of sinusoidal waves
sin(ωk t) and cos(ωk t) with ωk = kω = 2πk/T (k = 1, 2, ...) using Fourier transforms.
An important property of sinusoidal waves is that a wave with frequency ωk of any
amplitude Mk and phase angle φk can be written as the sum of a cosine and a sine
wave with the same frequency and zero phase angle:
Mk cos(ωk t − φk ) = Ak cos(ωk t) + Bk sin(ωk t) (1.1)
with Ak = Mk cos(φk ) and Bk = Mk sin(φk ).
Problem 1.1.1 : Prove relation (1.1) with the aid of the goniometrical relation:
cos(x + y) = cos(x) cos(y) − sin(x) sin(y) (1.2)
⇐

1.2 Fourier transforms

Fourier’s theorem states that any periodic function f (t) can be written as the sum
of an infinite series of terms:
∞
A0 X
f (t) = + (Ak cos(ωk t) + Bk sin(ωk t)) (1.3)
2 k=1
with
T
Z /2
2
Ak = f (t) cos(ωk t)dt k = 0, 1, ...
T
−T /2

T
Z /2 (1.4)
2
Bk = f (t) sin(ωk t)dt k = 1, 2, ...
T
−T /2

ωk = kω and ω = 2π/T

2
General introduction 3

With the aid of (1.1) equation (1.3) can also be written as:
∞
A0 X
f (t) = + Mk cos(ωk t − φk ) (1.5)
2 k=1

with Mk = (A2k + Bk2 ) and φk = arctan(Bk /Ak ) for k = 1, 2, ....

p

In practice the Fourier coefficients Ak and Bk are not computed exactly from (1.4).
Mostly the function f (t) is only defined (or measured) at a finite number of time
steps and a numerical integration is performed according to:
N
2 X
Ak = f (tn )cos(ωk tn ) k = 0, 1, 2, ..., N/2
N n=1
(1.6)
N
2 X
Bk = f (tn )sin(ωk tn ) k = 1, 2, ..., N/2
N n=1

with tn = −T /2 + nT /N . The coefficients Ak and Bk are called the discrete Fourier

coefficients. Note that the maximum frequency ωk in the spectrum of f (t) is deter-
mined by one half of the number of samples N in the time domain.

Problem 1.2.1 : The aortic flow is defined by a periodic function q(t) of time with
N=32 samples 0.0375 [s] apart from each other according to:

t=0:0.0375:1.2;
q=[ 3 3 4 29 327 486 466 378 287 217 166 116 62 -22 12 4 -5 4 9 -2 1 -1 -1 -5 0 0 1
1 1 2 2 2 3];

Check if the following program computes the Fourier coefficients as defined in (1.5) in
a proper way.
p11.m :
% p11.m: Fourier methods 1

clear

% Define the aortic flow as a function of time with N (even)

% discrete intervals

N=32;
T=1.2;
t=0:T/N:T;
q=[ 3 3 4 29 327 486 466 378 ...
287 217 166 116 62 -22 12 4 ...
-5 4 9 -2 1 -1 -1 -5 ...
0 0 1 1 1 2 2 2 3];

% Plot this function

subplot(211)
plot(t,q);
axis([0 1.2 -100 500])
4 Cardiovascular Fluid Mechanics - computational models

legend(’original function’)
xlabel(’time [s]’)
ylabel(’flow [ml/s]’)
title(’Aortic Flow’)
grid on
hold on

% Define the fundamental harmonic omega and the number of

% harmonics

om = 2*pi/T;
nh = N/2;

% Compute the Fourier coefficients

A = zeros(1,nh);
B = zeros(1,nh);
A0 = 2*sum(q(1:N))/N;
for ih = 1:nh
omega(ih) = ih*om;
A(ih) = 0;
B(ih) = 0;
for n=1:N
A(ih) = A(ih) + 2*q(n)*cos(omega(ih)*t(n))/N;
B(ih) = B(ih) + 2*q(n)*sin(omega(ih)*t(n))/N;
end
end

% Display the Fourier coefficients

subplot(212)
hold on

bar([0 1:nh],[A0/2 A],’b’);

bar([0 1:nh],[0 B],’r’);

legend(’A_0/2 A_k’,’B_k’)
xlabel(’Harmonic k’)
ylabel(’A_k and B_k [ml/s]’)
title(’Fourier coefficients’);
grid on

% Check the result

qr = zeros(size(t)) + A0/2;
for ih = 1:nh
qr = qr + A(ih)*cos(omega(ih)*t)+B(ih)*sin(omega(ih)*t);
end

subplot(211)
plot(t,qr,’o’);
axis([0 1.2 -100 500])
xlabel(’time [s]’)
ylabel(’flow [ml/s]’)
legend(’original function’,’reconstructed’)
title(’Aortic Flow’)
General introduction 5

grid on

1.3 Complex exponential forms

More compact notations (and manipulations) of sinusoidal waves can be obtained
using the relation:

eix = cos(x) + i sin(x) (1.7)

Problem 1.3.1 : Show that from this it can be derived that

1 ix
cos(x) = Re[eix ] = 2 (e + e−ix )
(1.8)
1 ix
sin(x) = Im[eix ] = 2i (e − e−ix )

and that, equation (1.1) also can be written as

h i
Mk cos(ωk t − φk ) = Re Mk ei(ωk t−φk ) (1.9)

Problem 1.3.2 : If we define

1
F (k) = (Ak − iBk ) (1.10)
2
show that
  " #
∞ ∞
iωk t  A0
2F (k)eiωk t
X X
f (t) = Re  F (k)e = + Re (1.11)
k=−∞
2 k=1

and
T /2
1
Z
F (k) = f (t)e−iωk t dt (1.12)
T
−T /2

again with ωk = kω.

⇐
This complex exponential form has great advantages and simplifies calculations and
finding of solutions to differential equations. The two reciprocal relations (1.11) and
(1.12) between a periodic function f (t) and its complex frequency spectrum F (k)
are called the inverse Fourier transform and the Fourier transform respectively.
6 Cardiovascular Fluid Mechanics - computational models

Problem 1.3.3 : Consider again the aortic flow defined by the periodic function q(t)
of time with N=32 samples 0.0375 [s] apart from each other according to:

t=0:0.1:1.2;
q=[ 3 3 4 29 327 486 466 378 287 217 166 116 62 -22 12 4 -5 4 9 -2 1 -1 -1 -5 0 0 1
1 1 2 2 2 3];

Check if the following program computes the Fourier coefficients as in (1.11) and com-
pute the complex frequency spectrum F (k).
p12.m :
% p12.m: Fourier methods 2

clear all

% Define the aortic flow as a function of time with N (even)

% discrete intervals

N=32;
T=1.2;
t=0:T/N:T;
q=[ 3 3 4 29 327 486 466 378 ...
287 217 166 116 62 -22 12 4 ...
-5 4 9 -2 1 -1 -1 -5 ...
0 0 1 1 1 2 2 2 3];

% Plot this function

subplot(211)
plot(t,q);
axis([0 1.2 -100 500])
legend(’original function’)
xlabel(’time [s]’)
ylabel(’flow [ml/s]’)
title(’Aortic Flow’)
grid on
hold on

% Define the fundamental harmonic omega and the number of

% harmonics

om = 2*pi/T;
nh = N/2;

% Compute the Fourier coefficients

A = zeros(1,nh);
B = zeros(1,nh);
A0 = 2*sum(q(1:N))/N;
for ih = 1:nh
omega(ih) = ih*om;
A(ih) = 0;
B(ih) = 0;
for n=1:N
A(ih) = A(ih) + 2*q(n)*cos(omega(ih)*t(n))/N;
General introduction 7

B(ih) = B(ih) + 2*q(n)*sin(omega(ih)*t(n))/N;

end
end

% compute the frequency spectrum

F = [ A0 (A - i*B)/2 ];

% Check the result

qr = zeros(size(t)) + A0/2;
for ih = 1:nh
qr = qr + real( 2*F(ih+1)*exp(i*omega(ih)*t) );
end

subplot(212)
plot(t,q,’-’,t,qr,’o’);
axis([0 1.2 -100 500])
xlabel(’time [s]’)
ylabel(’flow [ml/s]’)
legend(’original function’,’reconstructed’)
title(’Aortic Flow’)
grid on

⇐
Problem 1.3.4 : Show that the same result can be obtained by using the MATLAB
fft-function as given in the next program.
p13.m :
% p13.m: Fourier methods 3

clear

% Define the aortic flow as a function of time with N (even)

% discrete intervals

N=32;
T=1.2;
t=0:T/N:T;
q=[ 3 3 4 29 327 486 466 378 ...
287 217 166 116 62 -22 12 4 ...
-5 4 9 -2 1 -1 -1 -5 ...
0 0 1 1 1 2 2 2 3];

% Plot this function

subplot(211)
plot(t,q)
axis([0 1.2 -100 500])
legend(’original function’)
xlabel(’time [s]’)
ylabel(’flow [ml/s]’)
title(’Aortic Flow’)
grid on

% Define the fundamental harmonic omega and the number of

8 Cardiovascular Fluid Mechanics - computational models

% harmonics nh

om = 2*pi/T;
nh = N/2;

% Compute the frequency spectrum

F = fft(q(1:N))/N;

% Check the result

qr = zeros(size(t)) + real(F(1));
for ih = 1:nh
omega(ih) = ih*om;
qr = qr + real( 2*F(ih+1)*exp(i*omega(ih)*t) );
end

subplot(212)
plot(t,q,’-’,t,qr,’o’);
axis([0 1.2 -100 500])
xlabel(’time [s]’)
ylabel(’flow [ml/s]’)
legend(’original function’,’reconstracted (FFT)’)
title(’Aortic Flow’)
grid on

⇐
General introduction 9

1.4 The windkessel model

Once we have the complex frequency spectrum F (k) of the function f (t), we can
use this to solve for instance the differential equation that represents the windkessel
model relating aortic pressure to aortic flow.
∂pa pa
qa (t) = C + (1.13)
∂t R
If we now use the Fourier transformed equation for each harmonic k
1
Qa (k) = iωk CPa (k) + Pa (k) k = 0, 1, ... (1.14)
R
or:
R
Pa (k) = Qa (k) (1.15)
1 + iωk RC
the aortic pressure pa (t) can be obtained by the inverse Fourier transform:
∞
" #
Pa (0)
2Pa (k)eiωk t
X
pa (t) = + Re (1.16)
2 k=1

Problem 1.4.1 Use the function qa (t) and its transform Qa (k) as flow-input for equa-
tion (1.13) and compute the corresponding pressure pa (t). Play with the parameters C
and R to obtain a physiological pressure pulse. Remind that the model we use here is
not very sophisticated.
p14.m :
% p14.m: Fourier methods for windkessel model

clear

% Define the aortic flow as a function of time with N (even)

% discrete intervals

N=32;
T=1.2;
t=0:T/N:T;
qa=[ 3 3 4 29 327 486 466 378 ...
287 217 166 116 62 -22 12 4 ...
-5 4 9 -2 1 -1 -1 -5 ...
0 0 1 1 1 2 2 2 3];

% Plot this function

subplot(211)
plot(t,qa);
axis([0 1.2 -100 500])
xlabel(’time [s]’)
ylabel(’flow [ml/s]’)
title(’Aortic Flow’)
grid on

% Define the fundamental harmonic and omega

10 Cardiovascular Fluid Mechanics - computational models

om = 2*pi/T;
nh = N/2;

% Compute the frequency spectrum

QA = fft(qa(1:N))/N;

% Compute PA or even better pa = R/(1+i*omega*RC)*qa

R = 0.1; % Periferal resistance [kPa s/ml]

C = 5; % Arterial compliance [ml/kPa]

pa = zeros(size(t))+R*real(QA(1));
for ih = 1:nh
omega(ih) = ih*om;
pa = pa + real( (R/(1+i*omega(ih)*R*C))*2*QA(ih+1)*exp(i*omega(ih)*t) );
end
subplot(212)
plot(t,pa);
axis([0 1.2 0 20])
xlabel(’time [s]’)
ylabel(’pressure [kPa]’)
title(’Aortic Pressure’)
grid on

Problem 1.4.2 Check the following program that will optimize R and C, play with
the initial values and try to extend it with other components.
p15.m :
% windk.m
%========
clear
clf

global qao pao t

% load pressure and flow

load lvpq;

% define time axis and plot

nt = length(pao);
t = 0:1/(nt-1):1;
subplot(211)
plot(t,qao);

% define initial values for R and C

R = mean(pao)/mean(qao);
% R = 1/25;
C = 2/R;
General introduction 11

% compute pressure from qao, R and C

pw = windkp(qao, t, R, C);

subplot(212)
plot(t,pao,t,pw,’o’);

% estimate C

guess = [R C];
options(1) = 1;
est = fminsearch(’fwindkp’,guess,options);

R
R = est(1)
C
C = est(2)

pw = windkp(qao, t, R, C);

subplot(212)
plot(t,pao,t,pw,’o’);

subplot(211)
qest=min(qao)+(max(qao)-min(qao))*(pao-pw-min(pao-pw))/max((pao-pw)-min(pao-pw));
plot(t,qest,’o’,t,qao);

%windkp.m
%========

function pw = windkp(q, t, R, C)

% compute pressure from q, R, C

nt = length(t);
fq = fft(q)/nt;

qtot = zeros(size(t))+real(fq(1));
pw = qtot*R;
nh = floor(nt/2);

f1=1/t(nt);
for ih = 1:nh,
f(ih) = ih*f1; % define frequency(1:nh)
omega(ih) = 2*pi*f(ih); % define omega (1:nh)
qq(ih,:) = real(2*fq(ih+1)*exp(i*omega(ih)*t)); % flow (nh:nt)
z = 1/(i*omega(ih)*C+1/R);
fpw(ih+1) = fq(ih+1)*z;
p(ih,:) = real(2*fpw(ih+1)*exp(i*omega(ih)*t)); % flow (nh:nt)
pw = pw + p(ih,:); % total flow
end

%fwindkp.m
%=========

function f = fwindkp(par)
12 Cardiovascular Fluid Mechanics - computational models

global qao pao t

f = windkp(qao, t, par(1) , par(2) ) - pao;

f = sum(f.^2);

⇐
Chapter 2

Newtonian flow in blood vessels

The flow in arteries in many cases can be considered as an instationary boundary

layer flow. In order to obtain more insight in the parameters that determine this
kind of flow, the flow between two parallel plates will be considered first. In a later
stage of this course the same parameters will appear in a more realistic model of
pulsating artery flow.

2.1 Flow between two parallel moving plates

Consider incompressible Newtonian flow between two parallel plates of infinite length
given by the planes y = 0 and y = h (see figure 2.1). The upper plate can be moved
in x-direction with any desired velocity V (t).

V (t)
y=h

y=0
x

Figure 2.1: Two parallel moving plates

Problem 2.1.1 : Show that the velocity v(y, t) is given by the diffusion equation:

∂2v

∂v
− ν =0


∂y 2




 ∂t

(2.1)

 v(0, t) = 0





v(h, t) = V (t)


13
14 Cardiovascular Fluid Mechanics - computational models

If the velocity V (t) is periodic in time, we may assume that also v(y, t) is periodic.
Let V (t) = V̂ cos(ωt) with ω a given angular frequency. We will search for a solution
of the form:

v = v̂eiωt (2.2)
Substitution of (2.2) in (2.1) gives:
iω ∂2
v̂ − 2 v̂ = 0 (2.3)
ν ∂y

Problem 2.1.2 :
Show that substitution of
p iω p iω
y y
v̂(y) = c1 e ν + c2 e− ν (2.4)
yields the solution:
e(1+i)y/δ − e−(1+i)y/δ
r
2ν
v̂ = V̂ (1+i)h/δ with δ= (2.5)
e − e−(1+i)h/δ ω
√ √
Note: i = 12 2(1 + i)

Problem 2.1.3 :
Use the following MATLAB code to investigate the influence of the parameter δ.
p21.m :

%Flow between two parallel moving plates

clear
clg

% Create a ’time axis’ and the velocity of the upper plate

T = 1;
t=0:T/100:T;
omega = 2*pi;
Vhat = 1e0;
V = Vhat * exp(i*omega*t);

% Plot the velocity of the upper plate as a function of time

subplot(121);
plot(t/T,real(V)/Vhat);
xlabel(’t/T [-]’);
ylabel(’V/V0 [-]’);
title(’Velocity of upper plate’);
grid on

% Create an ’y-axis’ and the solution of the diffusion problem

h = 1;
Newtonian flow in blood vessels 15

y=0:h/100:h;
% nu = 1e-1;
nu = input(’nu = ’);
delta = sqrt(2*nu/omega);
ey = y*(1+i)/delta;
eh = h*(1+i)/delta;

vhat = Vhat * (exp(ey)-exp(-ey))./(exp(eh)-exp(-eh));

% Plot the velocity

subplot(122);
for time=0:0.1:1
v=real(vhat * exp(i*omega*time));
plot(y/h,v/Vhat);
hold on;
end
xlabel(’y/h [-]’)
ylabel(’v/V0 [-]’)
str=[’delta =’,num2str(delta)];
text(0.1,0.7,str);
title(’Velocity between plates’)
grid on

⇐
Now we will analyze the flow for more complicated movements of the upper wall.
First we will assume the upper wall velocity to be:

V (t) = V̂ + V̂ cos(ωt) (2.6)

Problem 2.1.4 :
Prove that the solution of the problem will be like:
" ! #
y e(1+i)y/δ − e−(1+i)y/δ iωt
v(y, t) = Re V̂ + V̂ e (2.7)
h e(1+i)h/δ − e−(1+i)h/δ
Adjust the code given above to incorporate the steady component of the upper wall
velocity.
p22.m :
%Flow between two parallel moving plates
clear
clg

% Create a ’time axis’ and the velocity of the upper plate

T = 1;
t=0:T/100:T;
omega = 2*pi;
V0 = 1e0;
Vhat = 1e0;
V = Vhat + Vhat * exp(i*omega*t);
16 Cardiovascular Fluid Mechanics - computational models

% Plot the velocity of the upper plate as a function of time

subplot(121);
plot(t/T,real(V)/V0);
xlabel(’t/T [-]’);
ylabel(’V/V0 [-]’);
title(’Velocity of upper plate’);
grid on

% Create an ’y-axis’ and the solution of the diffusion problem

h = 1;
y=0:h/100:h;
% nu = 1e-1;
nu = input(’nu = ’);
delta = sqrt(2*nu/omega);
ey = y*(1+i)/delta;
eh = h*(1+i)/delta;

vhat = Vhat * (exp(ey)-exp(-ey))./(exp(eh)-exp(-eh));

% Plot the velocity

subplot(122);
for time=0:0.1:1
v = Vhat*y/h + real(vhat * exp(i*omega*time));
plot(y/h,v/V0);
hold on;
end
xlabel(’y/h [-]’)
ylabel(’v/V0 [-]’)
str=[’delta =’,num2str(delta)];
text(0.1,1.7,str);
title(’Velocity between plates’)
grid on

⇐
A more general code is given below. Here the velocity of the upper plate is described
by 6 harmonics.
p23.m :
%Flow between two parallel moving plates
clear
clg

% Create a ’time axis’ and the velocity of the upper plate

T = 1;
t=0:T/100:T;
lt = length(t);
nhar = 6;
omega = [2*pi:2*pi:nhar*2*pi]’;
V0 = 1e0;

% Construct V from its harmonics

Newtonian flow in blood vessels 17

Vhat0 = 1e0;
Vhat = [1e0 0e0 0e0 0e0 0e0 0e0]’;

V = Vhat0;
for ihar = 1: nhar
V = V + Vhat(ihar) * exp(i*omega(ihar)*t);
end

% Plot the velocity of the upper plate as a function of time

subplot(121);
plot(t/T,real(V)/V0);
xlabel(’t/T [-]’);
ylabel(’V/V0 [-]’);
title(’Velocity of upper plate’);
grid on

% Create an ’y-axis’ and the solution of the diffusion problem

h = 1;
y=0:h/100:h;
ly = length(y);
nu = 1e-1;
nu = input(’nu = ’);

% Compute solution for each harmonic

vhat = zeros(nhar,lt);
for ihar = 1:nhar
delta = sqrt(2*nu/omega(ihar));
ey = y*(1+i)/delta;
eh = h*(1+i)/delta;
vhat(ihar,1:ly) = ...
Vhat(ihar) * (exp(ey)-exp(-ey))./(exp(eh)-exp(-eh));
end

% Plot the velocity

subplot(122);
for time=0:0.1:1
v = Vhat0*y/h;
for ihar = 1:nhar
v = v + real(vhat(ihar,1:ly) * exp(i*omega(ihar)*time));
end
plot(y/h,v/V0);
hold on;
end
xlabel(’y/h [-]’)
ylabel(’v/V0 [-]’)
str=[’delta =’,num2str(delta)];
text(0.1,0.5,str);
title(’Velocity between plates’)
grid on

Problem 2.1.5 : Check the code given above by using simple harmonics. Note that
a higher frequency should give the same results as a lower viscosity.
18 Cardiovascular Fluid Mechanics - computational models

Problem 2.1.6 : Finally we can run the code with the first 10 harmonics of an aorta
flow pulse.

Vhat0 = 0.11;
Vhat = [ 0.1436-0.1421i ...
0.0111-0.1566i ...
-0.0529-0.0884i ...
-0.0481-0.0391i ...
-0.0400-0.0246i ...
-0.0409-0.0096i ...
-0.0302+0.0068i ...
-0.0168+0.0088i ...
-0.0137+0.0060i ...
-0.0113+0.0099i ];

⇐
Newtonian flow in blood vessels 19

2.2 Bessel functions

Bessel functions Jn (s) play an important role in many physical problems in cylin-
drical coordinate systems and are solutions of:
∂ 2 Jn 1 ∂Jn n2
+ + (1 − )Jn = 0 (2.8)
∂s2 s ∂s s2
For integer values of n the Bessel functions are given by:
2k+n
∞
X (−1)k s
Jn (s) = (2.9)
k=0
k!(n + k)! 2

In MATLAB the Bessel functions can be computed with the aid of the function
bessel(n,s).
Problem 2.2.1 : Examine the properties of J0 and J1 graphically for real and complex
arguments. Use the following program as an example:
p41.m :
% Examine the Bessel function J0
clear

% First a rough estimate

smax = input(’maximum value of s=smax+i*smax = ’);
ds=smax/32;
[x,y]=meshgrid(0:ds:smax,0:ds:smax);
s=x+i*y;
J0=bessel(0,s);

figure(1)

subplot(221);
surf(x,y,real(J0));
xlabel(’real(s)’)
ylabel(’imag(s)’)
title(’real(J0)’);
shading interp
subplot(222);
surf(x,y,imag(J0));
xlabel(’real(s)’)
ylabel(’imag(s)’)
title(’imag(J0)’);
shading interp

subplot(223);
surf(x,y,abs(J0));
xlabel(’real(s)’)
ylabel(’imag(s)’)
title(’abs(J0)’);
shading interp
subplot(224);
surf(x,y,unwrap(angle(J0))/pi);
xlabel(’real(s)’)
ylabel(’imag(s)’)
title(’phase(J0)*pi’);
20 Cardiovascular Fluid Mechanics - computational models

shading interp

% Now we use that in our application we have s = i^(3/2)*(real number)

% so s = x+i*x

x = 0:ds:smax;
s = x + i*x;
J0=bessel(0,s);

figure(2)
subplot(221)
plot(x,real(J0));
xlabel(’x’);
ylabel(’real(J0)’);
title(’real(J0)’);
subplot(222)
plot(x,imag(J0));
xlabel(’x’);
ylabel(’imag(J0)’);
title(’imag(J0)’);
subplot(223)
plot(x,abs(J0));
xlabel(’x’);
ylabel(’abs(J0)’);
title(’abs(J0)’);
subplot(224)
plot(x,unwrap(angle(J0))/pi);
xlabel(’x’);
ylabel(’imag(J0)’);
title(’angle(J0)*pi’);

% Finally we will used the argument s= i^(3/2)*alpha*r/a and scale

% with i^(3/2)*alpha

alpha = input(’alpha = ’);

s0 = i^(3/2)*alpha;
r = 0:1/33:1;
s = s0*r;
J0 = bessel(0,s);
J00 = bessel(0,s0);

figure(3)
subplot(221)
plot(r,real(J0/J00));
xlabel(’r’);
ylabel(’real(J0/J00)’);
title(’real(J0/J00)’);
subplot(222)
plot(r,imag(J0/J00));
xlabel(’r’);
ylabel(’imag(J0/J00)’);
title(’imag(J0/J00)’);
subplot(223)
plot(r,abs(J0/J00));
xlabel(’r’);
ylabel(’abs(J0/J00)’);
Newtonian flow in blood vessels 21

title(’abs(J0/J00)’);
subplot(224)
plot(r,unwrap(angle(J0/J00))/pi);
xlabel(’r’);
ylabel(’imag(J0/J00)’);
title(’angle(J0/J00)*pi’);

2.3 Flow in the aorta

The complex Fourier coefficients for the flow in the aorta are given by:

n
0 0.1100
1 0.1436 -0.1421i
2 0.0111 -0.1566i
3 -0.0529 -0.0884i
4 -0.0481 -0.0391i
5 -0.0400 -0.0246i
6 -0.0409 -0.0096i
7 -0.0302 0.0068i
8 -0.0168 0.0088i
9 -0.0137 0.0060i
10 -0.0113 0.0099i

The diameter of the aorta is assumed to be 2cm, the dynamic viscosity η of blood
equals 5 · 10−3 P a · s, the density ρ is 103 kg/m3 . We want to compute the velocity
profiles in the aorta for a heart beat of 1 per second assuming the aorta to be a rigid
tube.

Problem 2.3.1 : Use the following MATLAB code to reconstruct the aortic flow from
its Fourier coefficients and decide how many harmonics must be taken into account to
obtain a representative flow pulse.

⇐
p42.m :
clear;
% Here are the first ten complex Fourier coefficients of the aortic flow
% in [l/s].

Q = [ 0.1100
0.1436 - 0.1421*i
0.0111 - 0.1566*i
-0.0529 - 0.0884*i
-0.0481 - 0.0391*i
-0.0400 - 0.0246*i
-0.0409 - 0.0096*i
22 Cardiovascular Fluid Mechanics - computational models

-0.0302 + 0.0068*i
-0.0168 + 0.0088*i
-0.0137 + 0.0060*i
-0.0113 + 0.0099*i ]; % flow [l/s]

% Set some parameters

nh = 8; % number of harmonics to take into account

f1 = 1; % frequency of the first harmonic 1 Hz
nt = 64; % number of timesteps
t = linspace(0,1/f1,nt); % time axis

% reconstruct the flow from its Fourier coefficients

qtot = zeros(size(t)); % initialize the total flow

% loop over the harmonics

for ih = 0:nh,

% define the frequency

f(ih+1) = ih*f1; % define frequency (1:nh)

omega(ih+1) = 2*pi*f(ih+1); % define omega (1:nh)

%compute the flow

q(ih+1,:) = Q(ih+1)*exp(i*omega(ih+1)*t); % flow (nh:nt)

qtot = qtot + q(ih+1,:); % total flow (1:nt)
end

% plot the total flow

plot(t,real(qtot)*1000);
title(’Aortic flow’);
xlabel(’t/T [-]’);
ylabel(’flow [ml/s]’);
grid

In order to compute the velocity profiles from

ˆ
" #
i ∂p J0 (i3/2 αr/a)
v̂z (r) = 1− (2.10)
ρω ∂z J0 (i3/2 α)
∂p
we need the pressure gradient ∂z . This can be derived from the longitudinal impedance
ZL and flow q according to:
ˆ
∂p
ZL = − /q̂ (2.11)
∂z
and:
ρ 1
ZL = iω 2
(2.12)
πa 1 − F10 (α)
We will first make a MATLAB function to compute the Womersley function F10
using the following program:
Newtonian flow in blood vessels 23

%f10.m calculates the womersley function for a given alpha

function [f] = f10(alpha)

s=alpha*i^(3/2);
if s ~= 0,
h1=2*bessel(1,s);
h2=s*bessel(0,s);
f=h1/h2;
else
f=1+0*i;
end

Problem 2.3.2 : Check this function by comparing its output with the figure given in
the lecture notes.

⇐
The longitudinal impedance ZL of the aorta now can be computed as follows:

% zl.m calculates the relative longitudinal impedance

% for a given alpha

function [z] = zl(alpha)

if alpha ~= 0,
z = (i*alpha*alpha/8)/(1-f10(alpha));
else
z = 1 + 0*i;
end

Problem 2.3.3 : Plot the longitudinal impedance as a function of α and interpret the
result.

⇐
The computation of the pressure gradient now is straight forward and can be per-
formed by a simple extension of the program that was used to reconstruct the aortic
flow:
p43.m :
clear;
% Here are the first ten complex Fourier coefficients of the aortic flow
% in [l/s].

Q = [ 0.1100
0.1436 - 0.1421*i
0.0111 - 0.1566*i
-0.0529 - 0.0884*i
-0.0481 - 0.0391*i
-0.0400 - 0.0246*i
-0.0409 - 0.0096*i
24 Cardiovascular Fluid Mechanics - computational models

-0.0302 + 0.0068*i
-0.0168 + 0.0088*i
-0.0137 + 0.0060*i
-0.0113 + 0.0099*i ]; % flow [l/s]

% Set some parameters

nh = 8; % number of harmonics to take into account

f1 = 1; % frequency of the first harmonic [1/s]
a = 0.01; % radius of the aorta [m]
eta = 5e-3; % dynamic viscosity [Pa.s]
rho = 1e3; % density [kg/m^3]
nu = eta/rho; % kinematic viscosity [m^2/s]
nt = 64; % number of timesteps [-]
t = linspace(0,1/f1,nt); % time axis

% reconstruct the flow from its Fourier coefficients

qtot = zeros(size(t)); % initialize the total flow

pztot = zeros(size(t)); % initialize the total pressure gradient
PZ = zeros(size(Q)); % initialize pressure gradient (Fourier)
Z = zeros(size(Q)); % initialize impedance

% loop over the harmonics

for ih = 0:nh,

% define the frequency

f(ih+1) = ih*f1; % define frequency (1:nh)

omega(ih+1) = 2*pi*f(ih+1); % define omega (1:nh)

% compute the flow

q(ih+1,:) = Q(ih+1)*exp(i*omega(ih+1)*t); % flow (nh:nt)

qtot = qtot + q(ih+1,:); % total flow (1:nt)

% compute alpha and longitudinal impedance Z

alpha(ih+1) = a*sqrt(omega(ih+1)/nu);
Z(ih+1) = 8*eta*zl(alpha(ih+1))/(pi*a^4);

% compute the pressure gradient (divide by 1000 because Q is in l/s)

PZ(ih+1) = - Z(ih+1)*Q(ih+1)/1000;
pz(ih+1,:) = PZ(ih+1)*exp(i*omega(ih+1)*t);
pztot = pztot + pz(ih+1,:);
end

% plot the total flow

subplot(211);
plot(t,real(qtot)*1000);
title(’Aortic flow’);
xlabel(’t/T [-]’);
ylabel(’flow [ml/s]’);
Newtonian flow in blood vessels 25

grid

% plot the total pressure gradient (divide by 1000 to make kPa/m)

subplot(212);
plot(t,real(pztot)/1000);
title(’Aortic pressure gradient’);
xlabel(’t/T [-]’);
ylabel(’pressure gradient [kPa/m]’);
grid

We will extend the program to compute the velocity profiles. The Womersley profiles
can be computed with the function:

% womer.m computes the Womersley profiles for

% given a, omega, rho and nu and for -a<r<a

function [w] = womer(a, omega, rho, nu ,r)

alpha = a*sqrt(omega/nu);

s = i^(3/2)*alpha;
J00 = bessel(0,s);
s = i^(3/2)*alpha*r/a;
J0 = bessel(0,s);
w = (i/(rho*omega))*(1-J0/J00);

Problem 2.3.4 : Use the code as given below to extend the program and check the
Womersley profiles for α = 2, 4, 8, 16.

⇐
p44.m :
clear;
clg;
% Here are the first ten complex Fourier coefficients of the aortic flow
% in [l/s].

Q = [ 0.1100
0.1436 - 0.1421*i
0.0111 - 0.1566*i
-0.0529 - 0.0884*i
-0.0481 - 0.0391*i
-0.0400 - 0.0246*i
-0.0409 - 0.0096*i
-0.0302 + 0.0068*i
-0.0168 + 0.0088*i
-0.0137 + 0.0060*i
-0.0113 + 0.0099*i ]; % flow [l/s]

% Q = [ 0 1 0 0 0 0 0 0 0 0 0 ];
% Set some parameters
26 Cardiovascular Fluid Mechanics - computational models

nh = 10; % number of harmonics to take into account

f1 = 1; % frequency of the first harmonic [1/s]
a = 0.01; % radius of the aorta [m]
eta = 5e-3; % dynamic viscosity [Pa.s]
rho = 1e3; % density [kg/m^3]
nu = eta/rho; % kinematic viscosity [m^2/s]
nt = 64; % number of timesteps
t = linspace(0,1/f1,nt); % time axis
r = -a:a/20:a;

% reconstruct the flow from its Fourier coefficients

qtot = zeros(size(t)); % initialize the total flow

pztot = zeros(size(t)); % initialize the total pressure gradient
PZ = zeros(size(Q)); % initialize pressure gradient (Fourier)
Z = zeros(size(Q)); % initialize impedance
v = zeros(length(r),length(t));
u = zeros(length(r),length(t));

% loop over the harmonics

for ih = 0:nh,

% define the frequency

f(ih+1) = ih*f1; % define frequency (1:nh)

omega(ih+1) = 2*pi*f(ih+1); % define omega (1:nh)

% compute the flow

q(ih+1,:) = Q(ih+1)*exp(i*omega(ih+1)*t); % flow (nh:nt)

qtot = qtot + q(ih+1,:); % total flow (1:nt)

% compute alpha and longitudinal impedance Z

alpha(ih+1) = a*sqrt(omega(ih+1)/nu);
Z(ih+1) = 8*eta*zl(alpha(ih+1))/(pi*a^4);

% compute the pressure gradient (divide by 1000 because q is in [l/s])

PZ(ih+1) = - Z(ih+1)*Q(ih+1)/1000;
pz(ih+1,:) = PZ(ih+1)*exp(i*omega(ih+1)*t);
pztot = pztot + pz(ih+1,:);

% compute the velocity profile

if alpha(ih+1) ~= 0,
w = womer(a, omega(ih+1), rho, nu, r);
else
w = -(a+r).*(a-r)/(4*eta);
end
u = conj(w’)*exp(i*omega(ih+1)*t);
uu = PZ(ih+1)*u;
v = v + uu;
plot(r/a,real(uu(:,[8 16 24 32 40 48 56 64])));
axis([-1,1,-1.0,1.5])
Newtonian flow in blood vessels 27

xlabel(’r/a [-]’);
ylabel(’velocity [m/s]’);
title([’alpha = ’ num2str(alpha(ih+1))]);
grid on
pause

end

% plot the total flow

subplot(221);
plot(t,real(qtot));
title(’Aortic flow’);
xlabel(’t/T [-]’);
ylabel(’flow [l/s]’);
grid

% plot the total pressure gradient

subplot(223);
plot(t,real(pztot)/1000);
title(’Aortic pressure gradient’);
xlabel(’t/T [-]’);
ylabel(’pressure gradient [kPa/m]’);
grid

subplot(122);
plot(r/a,real(v(:,[8 16 24 32 40 48 56 64])));
title(’Velocity profiles’);
xlabel(’r/a [-]’);
ylabel(’v [m/s]’);
grid

% Make a nice plot

figure(2)
subplot(111)
surf(t,r/a,real(v));
shading interp
title(’Velocity profiles in a rigid arota’);
xlabel(’t/T [-]’);
ylabel(’r/a [-]’);
zlabel(’velocity [m/s]’);

Note that conj(w’) is used to transpose w without taking its complex conjugate.

Problem 2.3.5 : Use different columns Q to analyse the properties of profiles that
result from two or more harmonics.

⇐
28 Cardiovascular Fluid Mechanics - computational models
Chapter 3

Wave phenomena in blood

vessels

3.1 Traveling waves

Traveling waves can be represented by the following function:
p(z, t) = p̂ei(ωt−kz) (3.1)

with ω the angular frequency and k = kr + iki the complex wave number. We
consider the following function:
t−0.25 2
p(0, t) = e−( 0.1) (3.2)

Problem 3.1.1 : Analyze the influence of the value of the wave number k by changing
parameters that determine the Moens-Korteweg wave speed in the following program.

⇐
p61.m :
% Travelling wave

clg
clear

% Define pressure pulse at z=0

T = 1;
N=64;
t=0:T/N:T;

p0 = exp(-((t-0.25)/0.1).^2);

plot(t,p0)
xlabel(’t/T [-]’);
ylabel(’p/p0 [-]’);
grid on

% Compute the complex Fourier coefficients

29
30 Cardiovascular Fluid Mechanics - computational models

om0 = 2*pi/T;
P0 = fft(p0)/N;

% Check what number of harmonics is OK

pc = zeros(size(t)) + real(P0(1));
nh = 8;
for ih = 1:nh
omega(ih) = ih*om0;
pc = pc + 2*real( P0(ih+1)*exp(i*omega(ih)*t) );
end
plot(t,p0,t,pc);
xlabel(’t/T [-]’);
ylabel(’p/p0 [-]’);
grid on

% Define the wave number k0 = omega/c0 with c0 the

% Moens-Korteweg wave speed

a0 = 1e-2;
h = a0/10;
E = 4.5e5;
rho = 1e3;
mu = 0.5;

c0 = sqrt(h*E/(2*rho*a0*(1-mu*mu)));
k0 = omega/c0;

% Compute the traveling wave for a number of z-locations

hold off
nz = 16;
dz = 1/nz;
for iz = 1: nz
z(iz)=(iz-1)*dz;
pz(iz,:) = zeros(size(t)) + real(P0(1));
for ih = 1: nh
k(ih) = k0(ih);
pz(iz,:) = pz(iz,:) ...
+ 2*real( P0(ih+1)*exp(i*omega(ih)*t-i*k(ih)*z(iz)) );
end

% Plot
plot(t,pz(iz,:));
xlabel(’t/T [-]’);
ylabel(’p/p0 [-]’);
grid on
hold on

end
Wave phenomena in blood vessels 31

In the next code the attenuation of the wave is introduced by incorporation of a

friction function according to
F10 (ω)
f0 (ω) = iωρ (3.3)
1 − F10 (ω)
and by assuming the tube to be viscoelastic according to

E = Er (1 + ifv ) (3.4)
p62.m :
% Travelling wave

clg
clear

% Define pressure pulse at z=0

T = 1;
N=64;
t=0:T/N:T;

p0 = exp(-((t-0.25)/0.1).^2);

plot(t,p0)
xlabel(’t/T [-]’);
ylabel(’p/p0 [-]’);
grid on

% Compute the complex Fourier coefficients

om0 = 2*pi/T;
P0 = fft(p0)/N;

% Check what number of harmonics is OK

pc = zeros(size(t)) + real(P0(1));
nh = 8;
for ih = 1:nh
omega(ih) = ih*om0;
pc = pc + 2*real( P0(ih+1)*exp(i*omega(ih)*t) );
end
plot(t,p0,t,pc);
xlabel(’t/T [-]’);
ylabel(’p/p0 [-]’);
grid on

% Define the wave number k0 = omega/c0 with c0 the

% Moens-Korteweg wave speed

a0 = 1e-2;
h = a0/10;
E = 4.5e5*(1+0.2*i);
rho = 1e3;
mu = 0.5;
eta = 3.5e-3;
32 Cardiovascular Fluid Mechanics - computational models

nu = eta/rho;

c0 = sqrt(h*E/(2*rho*a0*(1-mu*mu)));
k0 = omega/c0;

% Compute the traveling wave for a number of z-locations

hold off
nz = 16;
dz = 1/nz;
for iz = 1: nz
z(iz)=(iz-1)*dz;

% Compute k

alpha = a0*sqrt(omega(ih)/nu);
k = k0./sqrt(1-f10(alpha));

Zp = 8*eta/(pi*a0^4);
q0 = 1e-5;
pz(iz,:) = zeros(size(t)) + real(P0(1))-z(iz)*Zp*q0;
for ih = 1: nh
pz(iz,:) = pz(iz,:) ...
+ 2*real( P0(ih+1)*exp(i*omega(ih)*t-i*k(ih)*z(iz)) );
end

% Plot
plot(t,pz(iz,:));
xlabel(’t/T [-]’);
ylabel(’p/p0 [-]’);
grid on
hold on

end

Problem 3.1.2 : Investigate the attenuation of the waves by changing the viscous
properties of the fluid and the wall.

⇐
For a slowly tapered tube Pedley (1980) derives:

q(z, t) = Y (z)p(z, t) (3.5)

and
Y (0) 1/2 −ik̄z
p(z, t) = p(0, t){ } e (3.6)
Y (z)

with
Zz
1
k̄ = kdz ′ (3.7)
z
0
Wave phenomena in blood vessels 33

Problem 3.1.3 : Compute the pressure flow and velocity profiles as a function of
z and t in the aorta assuming that the diameter of the aorta changes from 2.5cm to
2.0cm over a distance of 25cm assuming that Young’s modulus changes from 2.5N/m2
to 4.5N/m2 . Take h = a0 (z)/10 and η = 5 · 10−3 P a · s.

⇐
34 Cardiovascular Fluid Mechanics - computational models
Chapter 4

Non-Newtonian flow in blood

vessels

4.1 The Carreau-Yasuda model

In the next MATLAB code the viscosity of blood is given as a function of the shear
rate. We will complete the program with code to fit these data with an appropriate
viscosity model.
First we will plot the viscosity as a function of the shear rate and try to find values
for the Carreau-Yasuda model:
p71.m :
% Viscosity model

clear

gdot = [0.01 0.02 0.05 ...

0.1 0.2 0.5 ...
1.0 2.0 5.0 ...
10 20 50 ...
100 200 500];
vis = [ 0.1837 0.1478 0.1000 ...
0.0648 0.0457 0.0260 ...
0.0154 0.0109 0.0071 ...
0.0050 0.0038 0.0032 ...
0.0027 0.0024 0.0023 ];

hold off
loglog(gdot,vis,’o’);
xlabel(’shear rate’)
ylabel(’viscosity’)
axis([1e-3 1e3 1e-3 1]);
grid

eta0 = 0.1;
eta1 = 0.01;
lambda = 1.0;
a = 1.0;
n = 1.0;

35
36 Cardiovascular Fluid Mechanics - computational models

eta = caryas(gdot, eta0, eta1, lambda, a, n);

hold on
loglog(gdot,eta);

Here the function caryas is defined by the MATLAB function caryas.m:

% The Carreau-Yasuda model

function eta = caryas(gdot,eta0,eta1,lambda,a,n)

eta = eta1 + (eta0 - eta1)*(1+(lambda*gdot).^a).^((n-1)/a);

Problem 4.1.1 : Implement the code given above and try to find proper values for
η0 , η∞ , λ, a and n.

⇐
Now we will use the MATLAB function fmins to estimate the parameters. The
function fmins finds a minimum of a function of several variables, starting from an
initial estimate. The call est = fmins(’fname’,guess) returns a vector est that
minimizes the function fname(est) starting from an initial guess vector guess. In
order to use fmins we need to write the function fname that has to be minimized.
In the case of the Carreau-Yasuda model this function (stored in fcaryas.m) reads:

% Compute the least square difference between the Carreau-Yasuda model

% and the experimental data. This is the function that is minimized
% by fmins.

function q = fcaryas(par)

global gdot vis

q = caryas(gdot,par(1),par(2),par(3),par(4),par(5)) - vis;
q = sum(q.^2);

Note that the variables gdot and vis are defined as global variables allowing them
to be referenced inside functions without passing them through the parameter list.
p72.m :
% Viscosity model
clear
global gdot vis

gdot = [0.01 0.02 0.05 ...

0.1 0.2 0.5 ...
1.0 2.0 5.0 ...
10 20 50 ...
100 200 500];
vis = [ 0.1837 0.1478 0.1000 ...
0.0648 0.0457 0.0260 ...
Non-Newtonian flow in blood vessels 37

0.0154 0.0109 0.0071 ...

0.0050 0.0038 0.0032 ...
0.0027 0.0024 0.0023 ];

hold off
loglog(gdot,vis,’o’);
xlabel(’shear rate’)
ylabel(’viscosity’)
axis([1e-3 1e3 1e-3 1]);
grid

eta0 = 0.2;
eta1 = 0.002;
lambda = 50.0;
a = 0.5;
n = 0.5;

eta = caryas(gdot, eta0, eta1, lambda, a, n);

hold on
loglog(gdot,eta);

guess = [eta0 eta1 lambda a n]’;

est = fmins(’fcaryas’,guess);
eta = caryas(gdot, est(1), est(2),est(3),est(4),est(5));
hold on
loglog(gdot,eta);

Problem 4.1.2 : Use the program given above to find the best estimation for the
parameters in the Carreau-Yasuda model. Remind that the result of the estimation
strongly depends on the initial guess that is provided.

4.2 The Quemada model

Problem 4.2.1 : Use a similar programs to fit the data with the Quemada viscosity
model.

⇐
38 Cardiovascular Fluid Mechanics - computational models

4.3 Viscoelastic materials

Both viscoelastic fluids and viscoelastic solids can be described by a two-mode
Maxwell model according to:

X − λ1 (t−tn ) − λ1 (t−tn )
τ (t) = δγn G0 e 0 + G1 e 1 (4.1)
n

Problem 4.3.1 : Show that for a viscoelastic solid with a finite relaxation time at
least two modes are needed to describe its response to changes in shear.

Problem 4.3.2 : Use the program given below to simulate the response of a Newto-
nian and a viscoelastic fluid and a Hookean and viscoelastic solid.
p74.m :
%p34.m : compute sigma of a visco-elastic material as a
% function of gamma and gdot using a 2-mode Maxwell
% model

clear
clg

% Set the material parameters

G0=1e5 ; % Shear modulus 0 [Pa]

l0=1e-3; % Time constant 0 [s]
G1=0e5 ; % Shear modulus 1 [Pa]
l1=1e3 ; % Time constant 1 [s]

% Define the shear gamma [-] as a ramp-function

dt=0.002;
N=501;
t=0:dt:(N-1)*dt;
t1 = 0.2;
t2 = 0.4;
for k=1:N
tk = t(k);
if (0 <= tk & tk < t1)
gamma(k)=0;
elseif ( t1 <= tk & tk < t2)
gamma(k)=0.25*(tk-t1)/(t2-t1);
elseif ( t2 <= tk & tk <= 1)
gamma(k)=0.25;
end
end

% plot gamma as a function of time

subplot(311)
plot(t,gamma);
grid on
legend(’Shear’);
Non-Newtonian flow in blood vessels 39

ylabel(’gamma [-]’);

% Compute shear rate gdot [1/s] by differentiation of gamma

gdot(1) = 0;
gdot(2:N) = (gamma(2:N)-gamma(1:N-1))/dt;

% Plot gdot as a function of time

subplot(312)
plot(t,gdot);
grid on
legend(’Shear rate’);
ylabel(’gdot [1/s]’);

% Compute tau

tau=zeros(1,N);

% Loop over the events in gamma

for n=2:N;
tn = t(n);
dgam = gamma(n)-gamma(n-1);
% Add influence of this event to stress at later times
for k=n:N;
tk = t(k);
dtau = dgam*(G0*exp(-(tk-tn)/l0)+G1*exp(-(tk-tn)/l1));
tau(k) = tau(k) + dtau;
end
end

% Plot tau

subplot(313)
plot(t,tau);
grid on;
xlabel(’t [sec]’);
ylabel(’tau [Pa]’);
legend(’Shear stress’)