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
Using derived in previous notebooks method we would like to fully implement numerical sulution to fluid dynamics!
First let me define all used imports:
import numpy as np
import math as mt
from scipy import sparse
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt
import matplotlib.ticker
np.set_printoptions(threshold=np.inf)
Next step would be introducing usefull classes and functions, which should provide easier to follow code:
class Operators:
def cncat(self, matrix_tuple, dim='v'):
if dim == 'v':
return np.concatenate(matrix_tuple, axis=0)
if dim == 'h':
return np.concatenate(matrix_tuple, axis=1)
def vec2col(self, vec):
return np.reshape(vec, (max(vec.shape), 1))
def vec2row(self, vec):
return np.reshape(vec, (1, max(vec.shape)))
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):
U[self.mesh[:-1, :]==1] = 0
V[self.mesh[:, :-1]==1] = 0
def get_patch(self):
return self.patches
def generate_laplacians(ubx='Neumann', uby="Neumann", vbx="Neumann", vby="Neumann", pbx="Neumann", pby="Neumann",
qbx="Dirichlet_between", qby="Dirichlet_between"):
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
def Lap(size, delta, boudary_cond='Dirichlet'):
Lap_row = np.concatenate((np.array([-1, 2, -1]), np.zeros(size-2)))
Lap_op = Lap_row[:]
for i in range(size-1):
Lap_op = np.concatenate((Lap_op, Lap_row))
deletion = [0] + list(range(-size+1, 0))
for i in deletion:
Lap_op = np.delete(Lap_op, i)
look_up_tab = {'Neumann': 1, 'Dirichlet': 2, 'Dirichlet_between': 3}
Lap_op[0] = Lap_op[-1] = look_up_tab[boudary_cond]
return np.reshape(Lap_op, (size, size)) / delta ** 2
def avg(A, k=1):
if A.shape[0] == 1:
A = np.transpose(A)
if k < 2:
B = (A[1::][:] + A[0:-1][:])/2
else:
B = avg(A, k-1)
if len(A.shape) < 2:
B = np.reshape(B, (1, len(B)))
return B
Time has come to initialaze varabiles, implemented here situation is steady west-east flow with few obstacles:
op = Operators()
Re = 1e1
dt = 1e-4
tf = 4e-0
nx = 40
ny = 40
lx = 1
ly = 1
nsteps = 2
eps = 2.2204e-16
nt = 10
obstacles = []
obstacles.append(Obstacle(0.1, 0.3, 0.2, 0.4))
obstacles.append(Obstacle(0.1, 0.3, 0.2, 0.85))
obstacles.append(Obstacle(0.1, 0.6, 0.4, 0.3))
obstacles.append(Obstacle(0.1, 0.6, 0.6, 0.7))
obstacles.append(Obstacle(0.1, 0.6, 0.8, 0.3))
nt0 = int(tf / dt)
x = np.linspace(0, lx, nx+1)
hx = lx / float(nx)
y = np.linspace(0, ly, ny+1)
hy = ly / float(ny)
[X, Y] = np.meshgrid(x, y)
for obs in obstacles:
obs.set_points(hx, hy, nx, ny)
#initial conditions
U = np.zeros((nx - 1, ny)) #nx - kolumny
V = np.zeros((nx, ny - 1)) #ny - wiersze
#boundary
uN = np.array(x*0)
vN = np.array(avg(x)*0)
uS = np.array(x*0)
vS = np.array(avg(x)*0+0.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)
'''
uN = np.array(x*0 + 4.0)
vN = np.array(avg(x)*0)
uS = np.array(x*0)
vS = np.array(avg(x)*0)
uW = np.array(avg(y)*0)
vW = np.array(y*0)
uE = np.array(avg(y)*0)
vE = np.array(y*0)#gora A
'''
a = op.vec2col(uS[1:-1]*2)
b = op.vec2col(uN[1:-1]*2)
c = op.vec2row(vW[1:-1]*2)
d = op.vec2row(vE[1:-1]*2)
Lu, Lv, Lp, Lq = generate_laplacians(uby='Dirichlet', vby='Dirichlet', pby='Dirichlet')
plt.ion()
plt.show()
fig = plt.figure()
tx = fig.add_subplot(111)
for obs in obstacles:
tx.add_patch(obs.get_patch()[0])
plt.draw()
plt.clf()
#trying to pring Lu
Lu
<1560x1560 sparse matrix of type '<class 'numpy.float64'>' with 7642 stored elements in Compressed Sparse Row format>
After initialisation what is left is a time loop connected with visualisation step:
Usage example, for low Raynolds number:
for i in range(nt0):
for obs in obstacles:
obs.insert_obj(U, V)
Ue = op.cncat((uW, U, uE))
Ue = op.cncat((2 * op.vec2col(uS)-op.vec2col(Ue[:, 0]), Ue, 2 * op.vec2col(uN)-op.vec2col(Ue[:,-1])), dim='h')
Ve = op.cncat((vS.T, V, vN.T), dim='h')
Ve = op.cncat((2 * op.vec2row(vW) - Ve[0, :], Ve, 2 * op.vec2row(vE) - Ve[-1, :]))
Ua = avg(Ue.T)
Ud = np.diff(Ue)/2.0
Va = avg(Ve)
Vd = np.diff(Ve, axis=0)/2.0
UVx = np.diff(Ua * Va, axis=0) / hx
UVy = np.diff(Ua * Va) / hy
Ua = avg(Ue[:, 1:-1])
Ud = np.diff(Ue[ :, 1:-1], axis=0)/2.0
Va = avg(Ve[1:-1][:].T).T
Vd = np.diff(Ve[1:-1][:])/2.0
U2x = np.diff(Ua**2, axis=0) / hx
V2y = np.diff(Va**2) / hy
U += - dt * (UVy[1:-1, :] + U2x)
V += - dt * (UVx[:, 1:-1] + V2y)
#velocity:
rhs = np.ravel(U)
U = spsolve(Lu, rhs, permc_spec='COLAMD')
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))
P = spsolve(Lp, -rhs, permc_spec='COLAMD')
P = np.reshape(P[:], (nx, ny))
U = U - np.diff(P, axis=0)/hx #kolumnowo
V = V - np.diff(P)/hy #wierszowo
if i == 1 or mt.floor(nsteps*i/nt) > mt.floor(nsteps*(i-1)/nt):
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:
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()
from IPython.display import YouTubeVideo
YouTubeVideo("https://youtu.be/gzwRNP38rRo")
For heigher Raynolds number something interesting start to happend:
from IPython.display import YouTubeVideo
YouTubeVideo("https://youtu.be/K_i3bvIbI5U")
Instead of West - East flow we can define flow towards the center:
from IPython.display import HTML
HTML('<iframe src=http://by2pie.wordpress.com width=900 height=500></iframe>')