Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Yves Dubief, 2015. NSF for support via NSF-CBET award #1258697.

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)

Steady Laminar Channel Flow

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.

Solution

Here

Taylor Series Expansion

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} $$

Solution

Here

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$.

Solution

Here

Discretization

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}\

\end{matrix}\right]

\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

Numerical solution

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]:
[<matplotlib.lines.Line2D at 0x10a5a0a58>]

Error Quantification

Visual verification is subjective. You should always rely on a statistical approach, such the $L^2$ or $L^\infty$.

$$ 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)


L2 error= 6.2815e-13
Linfinity error= 6.4837e-14

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 [ ]: