F01.justanswer.com



# SIMULATION OF BLOOD FLOW IN BLOOD VESSEL (NELECTING SHEAR WALL STRESS)import numpy as npimport matplotlib.pyplot as pltfrom matplotlib import animationfrom scipy import interpolatefrom numpy import wherefrom math import sinLNWDT=2; FNT=15plt.rcParams['lines.linewidth'] = LNWDTplt.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.5if (init_func==0): def p(x): p = np.zeros_like(x) p[where((x <= 2.5) & (x >= 7.5))] = p0 return pelif(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 relationdef A(p):A = np.zeros_like(p)A = A0 + C*(p-p0) return A# function defining the flux termdef 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 ################################################################ constantsC = 1.0rho = 10.0A0 = 5.0p0 = 10.0# discretizetmin, tmax = 0.0, 10.0 # start and stop time of simulationxmin, xmax = 0.0, 10.0 # start and end of spatial domainNx = 100 # number of spatial pointsCFL = 0.9 # courant number, need CFL<=1 for stabilityx = np.linspace(xmin, xmax, Nx+1) # discretization of spacep = p_init(x)c = np.sqrt(A0/(rho*C))Zc = rho*c/A0q = (np.max(p)-np.min(p))/Zcv = np.max(q/A(p))c += v qmin = np.min(q)q = np.ones_like(p)*qmindx = float((xmax-xmin)/Nx) # spatial step sizedt = CFL/c*dx # stable time step calculated from stability requirementNt = int((tmax-tmin)/dt) # number of time stepstime = np.linspace(tmin, tmax, Nt) # discretization of time# choose solversolvers = [macCormack]# solvers = [lax_friedrich]# solvers = [lax_friedrich, macCormack]# solve using all specified solverssolutions = [] # for storage of the solutionsfor 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 animatefig, 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 plotlines_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 framedef 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 sequentiallydef 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 figureax_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 figureplt.show() ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download