Simulation of Blood Flow in Blood Vessel (Nelecting Shear Wall Stress)

# SIMULATION OF BLOOD FLOW IN BLOOD VESSEL (NELECTING SHEAR WALL STRESS)

import numpy as np

import matplotlib.pyplot as plt

from matplotlib import animation

from scipy import interpolate

from numpy import where

from math import sin

LNWDT=2; FNT=15

plt.rcParams['lines.linewidth'] = LNWDT

plt.rcParams['font.size'] = FNT

# function defining the initial condition

# p(x) = p0 for x < 2.5 and x > 7.5

# p(x) = p0 + 0.25p0sin^2(pi*x-2.5/7.5-2.5) for 2.5 < x < 7.5

if (init_func==0):

def p(x):

p = np.zeros_like(x)

p[where((x <= 2.5) & (x >= 7.5))] = p0

return p

elif(init_func==1):

def p(x):

p = np.zeros_like(x)

x_left = 2.5

x_right = 7.5

xm = (x_right-x_left)/2.0

p = where((x>x_left) & (x<x_right), p0 + p0*0.25*np.sin(np.pi*(x-x_left)/(x_right-x_left))**2,f)

return f

# function defining the constitutive relation

def A(p):

A = np.zeros_like(p)

A = A0 + C*(p-p0)

return A

# function defining the flux term

def F(u):

# for Burger’s F = 0.5u^2 {but for blood flow need to define matrices??}

F = np.zeros_like(u)

F = 0.5u^2

return F

# function defining the macCormack scheme

def macCormack(u, t):

"""method that solves u(n+1), for the scalar conservation equation with source term:

du/dt + dF/dx = RHS,

where F = 0.5u^2 for the burger equation

with use of the MacCormack scheme

Args:

u(array): an array containg the previous solution of u, u(n). (RHS)

t(float): an array

Returns:

u[1:-1](array): the solution of the interior nodes for the next timestep, u(n+1).

"""

up = u.copy()

up[:-1] = u[:-1] - (dt/dx)*(F(u[1:]) - F(u[:-1])) + dt*RHS(t-0.5*dt, x[:-1])

u[1:] = .5*(u[1:] + up[1:] - (dt/dx)*(F(up[1:]) - F(up[:-1])) + dt*RHS(t-0.5*dt, x[1:]))

return u[1:-1]

### Main program ###############################################################

# constants

C = 1.0

rho = 10.0

A0 = 5.0

p0 = 10.0

# discretize

tmin, tmax = 0.0, 10.0 # start and stop time of simulation

xmin, xmax = 0.0, 10.0 # start and end of spatial domain

Nx = 100 # number of spatial points

CFL = 0.9 # courant number, need CFL<=1 for stability

x = np.linspace(xmin, xmax, Nx+1) # discretization of space

p = p_init(x)

c = np.sqrt(A0/(rho*C))

Zc = rho*c/A0

q = (np.max(p)-np.min(p))/Zc

v = np.max(q/A(p))

c += v

qmin = np.min(q)

q = np.ones_like(p)*qmin

dx = float((xmax-xmin)/Nx) # spatial step size

dt = CFL/c*dx # stable time step calculated from stability requirement

Nt = int((tmax-tmin)/dt) # number of time steps

time = np.linspace(tmin, tmax, Nt) # discretization of time

# choose solver

solvers = [macCormack]

# solvers = [lax_friedrich]

# solvers = [lax_friedrich, macCormack]

# solve using all specified solvers

solutions = [] # for storage of the solutions

for solver in solvers:

u = np.zeros((2,len(p))) # initialize u

u[0,:] = p # initialize u

u[1,:] = q # initialize u

un = np.zeros((len(time),2,len(x))) # for storage of the numerical solution

for i, t in enumerate(time[1:]):

u[:,1:-1] = solver(u[:,:]) # calculate numerical solution of interior

# enforce boundary conditions:

##########################

### fill in lines here ###

##########################

un[i,:,:] = u[:,:] # storing the solution for plotting

solutions.append(un) # adding the solution to the list of solutions

### animation ##################################################################

# first set up the figure, the axis, and the plot element we want to animate

fig, ax_l = plt.subplots()

ax_r = ax_l.twinx()

ax_l.set_ylim(p0*0.9, 1.25*p0*1.05)

ax_r.set_ylim(qmin*0.3, qmin*1.6)

ax_r.set_xlim(xmin, xmax)

# empty lines for initial plot

lines_l = [ax_l.plot([], [], '-', label='Pressure, '+solver.func_name)[0] for solver in solvers]

lines_r = [ax_r.plot([], [], ':', label='Volume flow, '+solver.func_name)[0] for solver in solvers]

# initialization function: plot the background of each frame

def init_dbl():

for line_l in lines_l:

line_l.set_data([], [])

for line_r in lines_r:

line_r.set_data([], [])

return lines_l, lines_r

# animation function, this is called sequentially

def animate_dbl(i):

for j, solution in enumerate(solutions):

lines_l[j].set_data(x, solution[i][0])

lines_r[j].set_data(x, solution[i][1])

return lines_l, lines_r

# annotate figure

ax_l.legend(lines_l, [l.get_label() for l in lines_l], loc=3, frameon=False)

ax_r.legend(lines_r, [l.get_label() for l in lines_r], loc=2, frameon=False)

ax_l.set_xlabel('x')

ax_l.set_ylabel('Pressure')

ax_r.set_ylabel('Volume flow')

# call the animator

anim = animation.FuncAnimation(fig, animate_dbl, init_func=init_dbl, frames=Nt, interval=50, blit=False)

# show the figure

plt.show()