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
I would like to develop a way to represent vector field $\vec{v} = (u,v)$ on staggered grid and as a result find update equations for its components
Since we take into consideration $3$ quantities, 2-compnent velocity vector $u,v$ and pressure let me consider their spacing
Image(filename='spacing.jpg', width=400)
In the picture above I used grid with $nx = 3, ny=2$ blue($U$), black($V$), red($P$) colors are used for interior points and green($U$), gray($V$), orange($P$) ones for boundary. Dimensions of each quantity can be deduced from the picture: $$dim(U_{interior}) = 2 \times 2 = n_x-1 \times n_y $$ $$dim(V_{interior}) = 3 \times 1 = n_x \times n_y-1 $$ $$dim(P_{interior}) = 3 \times 2 = n_x \times n_y $$
I would like to collect them in following table:
Quantity | Resolution of interior points | With boundary condition |
---|---|---|
pressure P | $n_x \times n_y$ | $n_x +2 \times n_y+2$ |
velocity x-component U | $(n_x-1) \times n_y$ | $n_x+1 \times n_y+2$ |
velocity y-component V | $n_x \times (n_y-1)$ | $n_x+2 \times n_y+1$ |
We can find these derrivative in equation <a href= >2.4 in previous notebook:
From equation it is easy to find linear equation for $U^{*}$:
Where we have linear operator $L_u$, yet to be found. Assuming we found operator we can state update equation in a form: $$$$ $$U^{**} = L_u^{-1}U^* \quad(2.5) $$
You may wonder how our Laplacian will look like. First of all we have to recall matrix representation of Laplacian. For seek of simplicity we would like to introduce 3x3 laplacian:
Where $a_{11}$ and $a_{33}$ are yet to be defined during boundary condition discussion.
If we introduce a grid represented by matrix 3x3 A, we should obtain follwing form of Laplace operator $\frac {\partial^2}{\partial x^2}+ \frac {\partial^2}{\partial y^2}$ :
To increase speed of calculation we would like to use vector form of our $U$ vector. It is obtained by command:
# U = internal structure + Ubc = boudary conditions
rhs = np.ravel((U+Ubc))
For row representation of coefficients L would take a form involving kronecker products with eye matrices,namely representation of 2.5 $L_u$, which inverse is capable of acting on rhs.
We should introduce a function generating Laplacian of given size. I am quite sure my function $Lap$ is not perfect and it can be upgraded but hey, it works! Here we have example of a code generating Laplacians $ L_u ->U-$ vector laplacian, $L_v->V-$ vector laplacian with already implemented boundary condition, Dirichlet by default:
from scipy import sparse
import numpy as np
def Lap(size, delta, boudary_cond='Dirichlet'):
"""Function returning Matrix representing Laplace operator"""
"""arguments:"""
"""size - size of matrix"""
"""delta - discretisation step"""
"""boundary_cond - type of boundary_condition """
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
nx = 3
ny = 3
Re = 1
dt = 1
hx = 1
hy = 1
Lu = sparse.eye(nx*(ny-1))+dt / Re*(sparse.kron(sparse.eye(ny), Lap(nx-1, hx)) +
sparse.kron(Lap(ny, hy), sparse.eye(nx-1)))
Lv = sparse.eye(nx*(ny-1))+dt / Re*(sparse.kron(sparse.eye(ny-1), Lap(nx, hx)) +
sparse.kron(Lap(ny-1, hy), sparse.eye(nx)))
print(Lu)
(0, 0) 3.0 (0, 1) -1.0 (0, 2) -1.0 (1, 0) -1.0 (1, 1) 3.0 (1, 3) -1.0 (2, 0) -1.0 (2, 2) 4.0 (2, 3) -1.0 (2, 4) -1.0 (3, 1) -1.0 (3, 2) -1.0 (3, 3) 4.0 (3, 5) -1.0 (4, 2) -1.0 (4, 4) 3.0 (4, 5) -1.0 (5, 3) -1.0 (5, 4) -1.0 (5, 5) 3.0
Since we are interested in values on the edges of the cell we can calculate velocity field using formula: $$ (U_x)_{i,j} \approx \frac{U_{i+\frac{1}{2},j} - U_{i-\frac{1}{2},j}}{h_x} \quad(2.2)$$
Then $U_x$ is "living" inside Cell in the same position as $P_{i,j}$, something similair happends to $(V_y)_{i,j}$
If you recall forumla 2.3 from previous notebook you know we have to find a way for treating nonlinear terms, namely to find a way to store those quantity in correct places.
We need $U^2$ in cell center and $UV$in the corner so interpolating neighboring values should be enough, so we have:
$$$$
import numpy as np
import math as mt
def avg(A, k=1):
"""
Returns horizontally averaged matrix
"""
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
#Copies of the field:
Ue = U
Ve = V
#Treating U^2, V^2 terms
Ua = avg(Ue[:, 1:-1]) #We are leaving boundary conditions untouched thus [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
Next in case of terms $UV$ we should use following formulas:
$$$$ $$\overline{V}^v_{i,j} = \frac{V_{i,j+\frac{1}{2}} + V_{i,j-\frac{1}{2}}}{2} \quad(2.3c)$$
Where subscribts indicate quantity is averager either averaged horizontally $^h$ or vertically $^v$, finally equation 2.3 from previous notebook changes into:
#UV Terms
Ua = avg(Ue.T) #vertically averaged
Ud = np.diff(Ue)/2.0
Va = avg(Ve) #horrizontally averaged
Vd = np.diff(Ve, axis=0)/2.0
UVx = np.diff(Ua * Va, axis=0) / hx
UVy = np.diff(Ua * Va) / hy
Recalling update equation we find:
#Update equation:
U += - dt * (UVy[1:-1, :] + U2x)
V += - dt * (UVx[:, 1:-1] + V2y)
Believe me or not but I tend to be lazy and probably this is the main reason for this class to exist. I decided to create operator class, which instance $op$ is usefull during concatenation and numpy matrix hokus-pokus operations.
class Operators:
"""
Class responsible for shorthand notation
Functions:
cncat - concatenation ->
matrix_tuple -> tuple of matrixes to concatenate
dim -> vertical/horizontal direction, vertical by default
vec2col -> transforming np array generated as (x,) to column vector (x,1)
vec2row -> transforming np array generated as (x,) to row vector (1,x)
"""
def cncat(self, matrix_tuple, dim='v'):
"""Concatenation"""
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)))
Since this is straight forward declaration I would like to stop here and focus on last notebook concerning: Boundary conditions problem, pressure correction, linear equation solver and visualisation.
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.