In [1]:
%matplotlib inline
# plots graphs within the notebook
%config InlineBackend.figure_format='svg' # not sure what this does, may be default images to svg format
from IPython.display import display,Image, Latex
from __future__ import division
from sympy.interactive import printing
printing.init_printing(use_latex='mathjax')
import time
from IPython.display import display,Image, Latex
from IPython.display import clear_output
#import SchemDraw as schem
#import SchemDraw.elements as e
import matplotlib.pyplot as plt
import numpy as np
import math
import scipy.constants as sc
import sympy as sym
from IPython.core.display import HTML
def header(text):
raw_html = '<h4>' + str(text) + '</h4>'
return raw_html
def box(text):
raw_html = '<div style="border:1px dotted black;padding:2em;">'+str(text)+'</div>'
return HTML(raw_html)
def nobox(text):
raw_html = '<p>'+str(text)+'</p>'
return HTML(raw_html)
def addContent(raw_html):
global htmlContent
htmlContent += raw_html
class PDF(object):
def __init__(self, pdf, size=(200,200)):
self.pdf = pdf
self.size = size
def _repr_html_(self):
return '<iframe src={0} width={1[0]} height={1[1]}></iframe>'.format(self.pdf, self.size)
def _repr_latex_(self):
return r'\includegraphics[width=1.0\textwidth]{{{0}}}'.format(self.pdf)
class ListTable(list):
""" Overridden list class which takes a 2-dimensional list of
the form [[1,2,3],[4,5,6]], and renders an HTML Table in
IPython Notebook. """
def _repr_html_(self):
html = ["<table>"]
for row in self:
html.append("<tr>")
for col in row:
html.append("<td>{0}</td>".format(col))
html.append("</tr>")
html.append("</table>")
return ''.join(html)
font = {'family' : 'serif',
#'color' : 'black',
'weight' : 'normal',
'size' : 12,
}
fontlabel = {'family' : 'serif',
#'color' : 'black',
'weight' : 'normal',
'size' : 16,
}
from matplotlib.ticker import FormatStrFormatter
plt.rc('font', **font)
The objective of this first assignment is to compute the well-known Poiseuille velocity profile in a channel flow. Assuming that the channel length and width are both very large, the governing equations for an incompressible flow, $$ \partial_xu+\partial_yv+\partial_z w = 0 $$ $$ \rho\left(\partial_t u+u\partial_x u+v\partial_y u+w\partial_z u\right) = -\partial_x p+\mu\left(\partial^2_x u +\partial^2_y u+\partial^2_z u\right) $$ $$ \rho\left(\partial_t v+u\partial_x v+v\partial_y v+w\partial_z v\right) = -\partial_y p+\mu\left(\partial^2_x v +\partial^2_y v+\partial^2_z v\right) $$ $$ \rho\left(\partial_t w+u\partial_x w+v\partial_y w+w\partial_z w\right) = -\partial_z p+\mu\left(\partial^2_x w +\partial^2_y w+\partial^2_z w\right) $$ reduce to
$$ 0=-\frac{d P}{dx}+\mu\frac{d^2u}{dx^2} $$
Demonstrate the above statement.
Given an $C^\infty$ function $f$ defined over an interval $]a,b[$, for any $x\in]a,b[$ and any $h$ such that $x+h\in]a,b[$, the Taylor series expansion writes
$$ f(x+h)=f(x)+\sum_{n=1}^\infty\frac{h^n}{n!}\frac{d^n f}{dx^n} $$
Practically, for a small increment $h$, $f(x+h)$ can be estimated as$$ f(x+h)=f(x)+\sum_{n=1}^{N-1}\frac{h^n}{n!}\frac{d^n f}{dx^n}+{\cal O}(h^N) $$
where$$ {\cal O}(h^N)=\sum_{n=N}^{\infty}\frac{h^n}{n!}\frac{d^n f}{dx^n} $$ is the truncation error which contains all terms affected by $h^n$ with $n\geq N$. Since $h$ is small, the large $n$, the smaller the contribution of $(h^n/n!)f^{(n)}(x)$. The above estimation of $f(x+h)$ is said to be of order $N$
Now, let's consider our channel, discretized as shown in the sketch below. This discretization technique is called finite difference method, which is used to solve PDEs and ODEs. Here our ODE is
$$ \mu\frac{d^2U}{dy^2}+G = 0\text{ with }G=-\frac{dP}{dx} $$
In [4]:
PDF('figures/channel.pdf',size=(200,400))
Out[4]:
Demonstrate that you can estimate the derivative of $U$ at node $i+1/2$, using $U_i$ and $U_{i-1}$
To answer this question, you will use the following Taylor series expansion: $$ U(y\pm\Delta y/2)=U(y)+\sum_{n=1}^\infty\frac{(\pm\Delta y/2)^n}{n!}\frac{d^nU}{dy^n} $$
Using the derivative above, show that a numerical scheme for the second derivative of $U$ is $$ \frac{d^2U}{dy^2}= \frac{U_{i-1}-2U_i+U_{i+1}}{\Delta y^2}+{\cal O}(\Delta y^2) $$
The term ${\cal O}(\Delta y^2)$ represents the truncation error, i.e. all the terms from the Taylor series expansion that are affected with $\Delta y^n$ for $n\geq 2$.Following the work above, a second order approximation of our ODE is
$$ \frac{\mu}{\Delta y^2}U_{j-1}-\frac{2\mu}{\Delta y^2}U_{j}+\frac{\mu}{\Delta y^2}U_{j+1}=-G $$
or$$ a_jU_{j-1}+b_jU_{j}+c_jU_{j+1}=-G $$
for $j=0,N-1$ At the boundaries, the velocity is imposed:$$ U_{j=0}=U_\text{lower wall} $$ $$ U_{j=N-1}=U_\text{upper wall} $$
This system of linear equation can be written in a matrix form: $$ \left[\begin{matrix} 1 & 0 & 0 & \cdots & \cdots & 0 \\ a_1 & b_1 & c_1 & \ddots & \cdots & \vdots \\ 0 & \ddots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & \cdots & \cdots & 0& 1 \\ \end{matrix} \right] \left[\begin{matrix} u_0\ u1 \ \vdots \ \vdots \ u{ny-1}\\left[\begin{matrix} u_\text{lower wall}\\ -G \\ \vdots \\ \vdots \\ u_\text{upper wall}\\ \end{matrix}\right] $$
The matrix being tridiagonal, we can solve this system using the efficient Thomas algorithm as follows:
In [5]:
import numpy as np
def TDMAsolver(a, b, c, d):
'''
TDMA solver, a b c d can be NumPy array type or Python list type.
refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
'''
nf = len(a) # number of equations
ac, bc, cc, dc = map(np.array, (a, b, c, d)) # copy the array
for it in range(1, nf):
mc = ac[it]/bc[it-1]
bc[it] = bc[it] - mc*cc[it-1]
dc[it] = dc[it] - mc*dc[it-1]
xc = ac
xc[-1] = dc[-1]/bc[-1]
for il in range(nf-2, -1, -1):
xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]
del bc, cc, dc # delete variables from memory
return xc
Here we will solve the flow inside a channel of half-height $h=1$, a pressure gradient $G=1$ and a viscosity $\mu=1/10$. Both walls are fixed.
Write the code to solve our ODE.
In [18]:
N = 256
mu = 1/10.
h = 1.
G = 1.
U_lower_wall = 0.
U_upper_wall = 0.
U = np.zeros(N)
U_all = np.zeros(N+2)
y_all = np.linspace(-h, +h, N+2)
dy = y_all[1] - y_all[0]
y = np.linspace(-h+dy, +h-dy, N)
#Matrix coefficients
a = mu/dy**2*np.ones(N)
b = -2.*mu/dy**2*np.ones(N)
c = mu/dy**2*np.ones(N)
#Boundary conditions
#Right hand side vector
r = -G*np.ones(N)
r[0] -= a[0]*U_lower_wall
r[N-1] -= c[N-1]*U_upper_wall
U = TDMAsolver(a, b, c, r)
U_all[1:N+1] = U
U_all[0] = U_lower_wall
U_all[N+1] = U_upper_wall
plt.plot(y_all,U_all,lw=2)
U_exact = G/(2*mu)*(h**2-np.power(y_all,2))
plt.plot(y_all,U_exact,'r--',lw=2)
Out[18]:
$$ L^2 = \left[\sum_{i=0}^{N-1}U^2_j\right]^{1/2} $$ $$ L^\infty = \max_{i=0,N-1}\left\vert U_j\right\vert $$
In [19]:
L2 = np.linalg.norm(U_all-U_exact,2)
print("L2 error= %1.4e" %L2)
Linf = np.linalg.norm(U_all-U_exact,np.inf)
print("Linfinity error= %1.4e" %Linf)
Now write a code for a non-Newtonian fluid, using the following viscosity
In [ ]:
def viscosity_carreau(u):
global mu_0,mu_inf,l_relax,n_power
n = len(u)
mu = np.zeros(n+1)
gradu = np.zeros(n+1)
gradu[1:n] = (u[1:n] - u[0:n-1])/dy
gradu[0] = 2.*u[0]/dy
gradu[n] = -2.*u[n-1]/dy
mu = mu_inf +(mu_0 - mu_inf)*np.power(1.+np.power(l_relax*gradu,2),n_power)
return nu
In [ ]: