%pylab inline
import pylab as pl
import math
import numpy as np
import matplotlib.pyplot as plt
import csv
# number of the grid points
ymax = 31;
# material properties
rho = 1000.0;
mu = 0.00102;
nu = mu / rho ;
#dimensions of the computational domain
l = 0.5 ;
h = 1 ;
w = 1.0 ;
# define Reynolds number
Re = 2000 ;
# compute the average velocity
ua = (Re*nu)/(2.0*h);
# compute the pressure drop
dp = (12.0*mu*l*ua)/(h*h);
# compute the volume flow rate
Q = (w*dp*h*h*h)/(12.0*mu*l);
# create the grid
dy = h/(ymax-1.0)
y = np.zeros(ymax)
y[0] = 0.0
for i in range(1,ymax):
y[i] = y[i-1]+dy
# compute the analytical solution
u =np.zeros(ymax)
#define the boundary conditions...define them inside the loop
u[0] = 0.0
u[ymax-1] = 0.0
const = dp/(2.0*mu*l )
for j in range(1,ymax-1):
u[j] = const * y[j]*(h-y[j])
print ua
plt.plot(u,y,'x')
plt.ylabel('Position [m]')
plt.xlabel('Velocity [m/s]')
plt.title('1D laminar velocity profile in a channel ')
plt.show()
Populating the interactive namespace from numpy and matplotlib 0.00102