from IPython.display import display, Image
from IPython.display import display_pretty, display_html, display_jpeg, display_png, display_json, display_latex, display_svg
from IPython.display import Latex
from IPython.display import HTML
Image(filename='logo3.jpg', width=100)
Project 2.1: Derivation of incompressible Navier - Stokes equations
Project 2.2: Time discretization of Navier-Stokes Equations in 2D
Project 2.3: Space discratization and update equations
Project 2.4: Project 2.4: Boundary condition, linear equation and visualization of N-S equations
Project 2.5: Full Code
In following notebook I would like to discuss three very important aspect of numerical solution. Especially visualization is something I did not mentioned before, mainly because it is independent of rest of code.
For $U$ and $V$ we should use Dirichlet boundary conditions, which are the simplest ones, which are fixing boundary point value. Of course in case of North and South boundaries for $U$ and West and East for $V$ we shall average two neighboring points in the following manner:
In case of implementation we would like to define initial conditions before entering time loop as set of 4x2 values representing 4 directions and two vector fields in following way:
uN = np.array(x*0)
vN = np.array(avg(x)*0)
uS = np.array(x*0)
vS = np.array(avg(x)*0)
uW = np.array(avg(y)*0+1.0)
vW = np.array(y*0)
uE = np.array(avg(y)*0+1.0)
vE = np.array(y*0)
It would mean we have constant West -> East flow
Those arrays we are going to concatenate with U and V interior arrays in first step of each iteration according to formula 1.2 as follows:
Ue = op.cncat((uW, U, uE))
Ue = op.cncat((2 * op.vec2col(uS) - Ue[:, 0], Ue, 2 * op.vec2col(uN)-Ue[:, -1]), dim='h')
Ve = op.cncat((vS.T, V, vN.T), dim='h')
Ve = op.cncat((2 * op.vec2row(vW) - Ve[1, :], Ve, op.vec2row(vE) - Ve[:, -1]))
We would like to have constant stream through boundaries so Neumann boundary conditions seems like a perfect match. Namely we would like to have zero flow, which means following equation holds:
For example for y-axis we have: $$\frac {P_{i,j+1} - P_{i,j}}{h} = 0 \quad Thus \quad P_{i,j+1} = P_{i,j}$$
In case of Laplacians there is one decisive parameter for inserting boundary conditions, namely: $$ -L_0 = \frac{1}{h^2}\begin{bmatrix} a_{11} & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & a_{11} \end{bmatrix} \quad(1.5) $$
$a_{11}$ | Boundary condition |
---|---|
1 | Neumann |
2 | Dirichlet |
3 | Dirichlet between points |
For example if we act with this operator on matrix of coordinates $X$ you should see why:
Usually boudary conditions have a form: $$L_u X = Constant$$ $$$$ $$\frac {a_{11}X_{11} + X_{21}}{h^2} = C$$ $$$$ $$a_{11} = -1 => C = \frac {X_{21} - X_{11}}{h} = \frac {dX}{dx}$$
Which means pressure Laplacian is equal to:
Lp = sparse.kron(sparse.eye(ny), Lap(nx, hx, boudary_cond='Neumann')) + \
sparse.kron(Lap(ny, hy, boudary_cond='Neumann'), sparse.eye(nx))
We are going to solve two types of linear equations:
For each velocity component and Poisson equation for pressure:
Here Python flexibility plays important role, to solve those linear expression we might use library associated with scipy, useful for solving sparse matrix equations:
from scipy.sparse.linalg import spsolve
As I mentioned before we want to obtain optimum speed of calculation, so we are going to use vector form of $U$ and $V$. It means our implementation of equations should look like this:
def generate_laplacians(ubx='Neumann', uby="Neumann", vbx="Neumann", vby="Neumann", pbx="Neumann", pby="Neumann",
qbx="Dirichlet_between", qby="Dirichlet_between"):
"""Returns laplacians for vector field (U, V) and scalar fields P and stream function Q
Arguments are just names of boundary conditions
"""
Lu = sparse.eye(nx*(ny-1))+dt / Re*(sparse.kron(sparse.eye(ny), Lap(nx-1, hx, boudary_cond=ubx)) +
sparse.kron(Lap(ny, hy, boudary_cond=uby), sparse.eye(nx-1)))
Lv = sparse.eye(nx*(ny-1))+dt / Re*(sparse.kron(sparse.eye(ny-1), Lap(nx, hx, boudary_cond=vbx)) +
sparse.kron(Lap(ny-1, hy, boudary_cond=vby), sparse.eye(nx)))
Lp = sparse.kron(sparse.eye(ny), Lap(nx, hx, boudary_cond=pbx)) + \
sparse.kron(Lap(ny, hy, boudary_cond=pby), sparse.eye(nx))
Lq = sparse.kron(sparse.eye(ny-1), Lap(nx-1, hx, boudary_cond=qbx)) + \
sparse.kron(Lap(ny-1, hy, boudary_cond=qby), sparse.eye(nx-1))
return Lu, Lv, Lp, Lq
#velocity:
rhs = np.ravel(U)
U = spsolve(Lu, rhs, permc_spec='COLAMD') #colamd is one of preparation method used in spsolve
U = np.reshape(U, (nx-1, ny))
rhs = np.ravel(V)
V = spsolve(Lv, rhs, permc_spec='COLAMD')
V = np.reshape(V, (nx, ny-1))
#pressure:
rhs = np.ravel((np.diff(op.cncat((uW, U, uE)), axis=0)/hx + np.diff(op.cncat((vS.T, V, vN.T), dim='h'))/hy)) #gradient
P = spsolve(Lp, -rhs, permc_spec='COLAMD')
P = np.reshape(P[:], (nx, ny)) # 2.2a
U = U - np.diff(P, axis=0)/hx # 2.2b
V = V - np.diff(P)/hy
Where in formulas 2.2a 2.2b $\Delta t$ canceled out.
Real - time visualization is possible due to incompressible flow argument. It means divergence of vector field $V$, $U$ is 0 so in analogy with electrodynamics we can introduce vector potential function call stream function $Q$ which is defined as follows:
$$\frac {\partial Q}{\partial x} = u \quad \frac {\partial Q}{\partial y} = -v$$
There is an operation called orthogonal gradient which is defined to match equations above, using it we find:
$$(\nabla Q)^o = \bigg[\frac{\partial Q}{\partial x},-\frac{\partial Q}{\partial y}\bigg] = {\bf U}$$
Applying curl on both sides yells:
$$\nabla \times (\nabla Q)^o = -\Delta Q^o = {\nabla \times{ \bf U} = \frac {\partial u}{\partial y} - \frac {\partial v}{\partial x}}$$
Which means calculating stream function Q turns into Poisson Equation. We are assuming Q as zero at the boundaries so final implementation has a following form:
rhs = np.diff(U)/hy - np.diff(V, axis=0)/hx
rhs = np.ravel(rhs)
q = spsolve(Lq, rhs, permc_spec='COLAMD')
Q = np.zeros((nx+1, ny+1))
Q[1:-1, 1:-1] = np.reshape(q, (nx-1, ny-1))
It can’t be that hard! For the last time we connect interior points of velocity field using quiver plot which combined with counter plot of stream function and colored contour plot of pressure should result in beautiful animation.
rhs = np.diff(U)/hy - np.diff(V, axis=0)/hx
rhs = np.ravel(rhs)
q = spsolve(Lq, rhs, permc_spec='COLAMD')
Q = np.zeros((nx+1, ny+1))
Q[1:-1, 1:-1] = np.reshape(q, (nx-1, ny-1))
[x_avg, y_avg] = np.meshgrid(avg(x), avg(y))
cxf = plt.contourf(x_avg, y_avg, P.T)
cx = plt.contour(X, Y, Q.T)
Ue = op.cncat((op.vec2col(uS), avg(op.cncat((uW, U, uE)).T).T, op.vec2col(uN)), dim='h')
Ve = np.concatenate((op.vec2row(vW), avg(op.cncat((op.vec2col(vS), V, op.vec2col(vN)), dim='h')), op.vec2row(vE)))
#Len = (Ue**2 + Ve**2 + eps)**0.5 #normalizing factor
#ax = plt.quiver(x, y, (Ue/Len).T, (Ve/Len).T, scale=50, headwidth=4, headlength=3, headaxislength=3)#with normalisation
#st = plt.streamplot(x, y, (Ue/Len).T, (Ve/Len).T) #stream function
ax = plt.quiver(x, y, Ue.T, Ve.T, scale=50, headwidth=4, headlength=3, headaxislength=3)
tx = fig.add_subplot(111)
for obs in obstacles: #obstacles
tx.add_patch(obs.get_patch()[0])
position=fig.add_axes([0.88,0.15,0.02,0.7])
cbar = plt.colorbar(cxf, cax = position)
plt.draw()
plt.clf()
You can ask me: "Hey obstacles? What obstacles" I decided to take into conideration ractangular obstacles. Thus I write a class obstacle in a following way:
class Obstacle():
def __init__(self, lx, ly, sx, sy, boundary_cond='Dirichlet'):
self.sx = sx
self.sy = sy
self.lx = lx
self.ly = ly
self.x0 = self.sx - self.lx/2.0
self.y0 = self.sy - self.ly/2.0
self.mesh = []
self.bc = boundary_cond
self.polygon = plt.Rectangle((self.x0, self.y0), self.lx, self.ly, fill=True)
self.patches = [self.polygon]
def set_points(self, hx, hy, nx, ny):
map = np.zeros((nx, ny))
map[round(self.x0/hx):round((self.x0+self.lx)/hx), round(self.y0/hy):round((self.y0+self.ly)/hy)] = 1
self.mesh = map
def get_points(self):
return self.mesh
def insert_obj(self, U, V): #set velocity to zero within an object
U[self.mesh[:-1, :]==1] = 0
V[self.mesh[:, :-1]==1] = 0
def get_patch(self):
return self.patches
Example of obstacle class usage:
obstacles = []
obstacles.append(Obstacle(0.1, 0.3, 0.2, 0.4))
for obs in obstacles:
obs.set_points(hx, hy, nx, ny)
#inside a loop:
for obs in obstacles:
obs.insert_obj(U, V)
It should produce following picture:
from IPython.display import YouTubeVideo
YouTubeVideo("jV8QKn6Xcuc")
from IPython.display import HTML
HTML('<iframe src=http://by2pie.wordpress.com width=900 height=500></iframe>')
[1]Seibold Benjamin, MIT: A compact and fast Matlab code solving the incompressible Navier-Stokes equations on ractangular domains,
[2]Matyka Maciej, Wroclaw University: Solution to two-dimensional Incompressible Navier-Stokes Equations with
SIMPLE, SIMPLER and Vorticity-Stream Function Approaches.
Driven-Lid Cavity Problem: Solution and Visualization.