Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 by Akash Dhruv

Navier-Stokes, Learn to Fly !!

1. Introduction

One of the most important application of numerical methods lies in solving the non-linear navier stokes equation (1a) and (1b). Introduced in 1822 [1] navier stokes equation describe the physics of a fulid flow accurately, but their correct analytical solution is not possible. Though some analytical solutions do exist for special flows, a general solution of these equations is hindered by the presence of the non linear term $\small u_j {\partial u_i}/{\partial x_j}$.

With the increase in computation power and capabilities over the decades, it is possible to numerically solve these equations to pratical accuracy. This has contributed greately towards efficient and cost effective design and analysis.

$$ {\small \frac{\partial u_i}{\partial t} +u_j \frac{\partial u_i}{\partial x_j}=-\frac {1}{\rho} \frac{\partial p}{\partial x_i} + \nu\bigtriangledown^2u_i} \quad \quad (1a) $$

$$ {\small \frac {\partial u_i}{\partial x_i} =0 } \quad \quad (1b) $$

$ u $- velocity, $\nu$- dynamic viscosity, $p$- Pressure

This notebook will present a numerical code to solve navier stokes equation for a Lid-Driven Cavity. Lid-Driven Cavity is a popular computation problem, where the top wall moves with a constant velocity parallel to itself, and rest of the walls are stationary [2-4]. This problem has been highly researched in both two dimensional and three dimensional domains and is a good example to validate a given solver [2]. The equations will be solved on a staggered grid, where unlike a collocated grid, velocities are located on cell faces, and pressure is located on the cell centre. These types of grid avoid odd-even decoupling between pressure and velocity in comparison to collocated grid [5], which makes the soultion free of checkboard patterns. Odd-even decoupling occurs when neighbouring cells don't share same ancestors. Further discussion on different grids can be found in ref [2-4]. Towards the end, an immersed boundary method is also presented. Immersed boundary methods offer advantage of using a cartesian grid over curvilinear grids when assessing the flow over complex geometries (multi-element airfoils, flying snake cross-section,etc) and moving boundaries (rotating cylinder). This is done in many different ways and thus there are all flavors of immersed boundary methods available. In this notebook, a direct forcing approach will be considered [6]. Direct forcing approach is a brute force technique wherein the $u$ and $v$ velocity points just outside the immersed boundary are determined and its wall velocity is imposed on them. This notebook could serve as a good starting point to understand the basic formulation of Navier-Stokes Solver.

2. Conservative form of Navier Stokes Equation

Equation 1(a) and (b) represent what is popularly called, the non-conservative form of navier stokes equations. Let us perform some substituion using the following steps,

$$ \small {\dfrac{\partial u_i}{\partial t} +u_j \dfrac{\partial u_i}{\partial x_j}=-\dfrac {1}{\rho} \dfrac{\partial p}{\partial x_i} + \nu\bigtriangledown^2u_i} $$$$\small \implies {\dfrac{\partial u_i}{\partial t} +\dfrac{\partial u_i u_j}{\partial x_j}-u_i \dfrac{\partial u_j}{\partial x_j}=-\dfrac {1}{\rho} \dfrac{\partial p}{\partial x_i} + \nu\bigtriangledown^2u_i} $$$$\small \implies {\dfrac{\partial u_i}{\partial t} +\dfrac{\partial u_i u_j}{\partial x_j}=-\dfrac {1}{\rho} \dfrac{\partial p}{\partial x_i} + \nu\bigtriangledown^2u_i}, ({Continuity \implies \dfrac{\partial u_j}{\partial x_j}=0}) $$$$\small Taking, P=\dfrac {p}{\rho} $$$$\small \implies {\dfrac{\partial u_i}{\partial t} +\dfrac{\partial u_i u_j}{\partial x_j}=-\dfrac{\partial P}{\partial x_i} + \nu\bigtriangledown^2u_i} \quad\quad(2)$$

This results in conservative form of Navier-Stokes equation, which is widely used when conservation of energy is an important concern [5]. For the problem discussed in this notebook and typically for incompressible flows, using either conservative or non-conservative is acceptable.

3. Numerical code on a staggered grid

Let us start working on the code. First we import all the neccesary libraries.



In [30]:

    
import numpy as np
import matplotlib.pyplot as plt
from math import *
from IPython.display import display
from IPython.display import Image

All the parameters to generate the grid and neccessary global constant are defined



In [31]:

    
# Defining Grid Parameters

Lx=1. #m
Ly=1. #m 

Nx=50
Ny=50

dx=Lx/(Nx-1)
dy=Ly/(Ny-1)


# Defining Global Constants

mu=0.01
rho=1. #kg/m^3
Re=100
U=(Re*mu)/Lx


w=1.1 # Relaxation factor for poisson solver

size=6 # Figure size, will be used for plotting

As discussed earlier, on staggered grids the locations of $\small u$, $\small v$ and $\small P$ are different. While $\small u$ and $\small v$ velocities are located on cell faces, the pressure points are on cell centres as seen in the figure below.



In [32]:

    
display(Image(filename='../Project_Akash/image12.jpeg'))

Let's implement this in the code. Starting with defining mesh points



In [33]:

    
# Staggered Grid

xstart=0
ystart=0

x=np.zeros((Nx,1),dtype=float)
y=np.zeros((Ny,1),dtype=float)

for i in range(Nx):
    x[i]=xstart+i*dx

for j in range(Ny):
    y[j]=ystart+j*dy

Next define $\small u$, $\small v$, and $\small P$ points



In [34]:

    
# Pressure points
xp=np.zeros((Nx+1,1),dtype=float)
yp=np.zeros((Ny+1,1),dtype=float)
p=np.zeros((Nx+1,Ny+1),dtype=float)

# u-velocity points
xu=np.zeros((Nx,1),dtype=float)
yu=np.zeros((Ny+1,1),dtype=float)
u=np.zeros((Nx,Ny+1),dtype=float)
ut=np.zeros((Nx,Ny+1),dtype=float)

# v-velocity points
xv=np.zeros((Nx+1,1),dtype=float)
yv=np.zeros((Ny,1),dtype=float)
v=np.zeros((Nx+1,Ny),dtype=float)
vt=np.zeros((Nx+1,Ny),dtype=float)

# Assigning grid values

# Pressure points
for i in range(1,Nx):
    for j in range(1,Ny):
        xp[i]=(x[i-1]+x[i])/2
        yp[j]=(y[j-1]+y[j])/2

xp[0]=xp[1]-dx
yp[0]=yp[1]-dy

xp[-1]=xp[-2]+dx
yp[-1]=yp[-2]+dy

# u-velocity points
for i in range(Nx):
    for j in range(1,Ny):
        xu[i]=x[i]
        yu[j]=(y[j-1]+y[j])/2

xu[0]=xu[1]-dx
yu[0]=yu[1]-dy

xu[-1]=xu[-2]+dx
yu[-1]=yu[-2]+dy

# v-velocity points
for i in range(1,Nx):
    for j in range(Ny):
        yv[j]=y[j]
        xv[i]=(x[i-1]+x[i])/2;

xv[0]=xv[1]-dx;
yv[0]=yv[1]-dy;

xv[-1]=xv[-2]+dx;
yv[-1]=yv[-2]+dy;

[X,Y]=np.meshgrid(x,y)



In [47]:

    
display(Image(filename='../Project_Akash/newgrid.png'))

So, we have created our staggered grid. Notice the points outside the mesh. These points will act as ghost cells when boundary condition is imposed. The next task is to discretize the navier stokes equation on this staggered grid.

Starting with conservative form of navier stokes equation, described in equation 2 and using taylor series expansion, the discrete experession is formed. Consider the staggered grid below.



In [36]:

    
display(Image(filename='../Project_Akash/ugrid.jpg'))

The discretized equation in $x$-direction can be written as follows

$$\small \dfrac {u_{i,j}^{n+1}-u_{i,j}^{n}}{\bigtriangleup t}+\dfrac {(u_{i+1/2,j}^{n})^2-(u_{i-1/2,j}^{n})^2}{\bigtriangleup x}+\dfrac {(u_{i,j+1/2}^{n})(v_{i,j+1/2}^{n})-(u_{i,j-1/2}^{n})(v_{i,j-1/2}^{n})}{\bigtriangleup y}=$$$$\small-\dfrac {P_{i+1/2,j}^{n}-P_{i-1/2,j}^{n}}{\bigtriangleup x}+\nu \dfrac{(u_{i+1,j}^{n}+u_{i-1,j}^{n}-2 u_{i,j}^{n})}{{\bigtriangleup x}^2}+\nu \dfrac{(u_{i,j+1}^{n}+u_{i,j-1}^{n}-2 u_{i,j}^{n})}{{\bigtriangleup y}^2} \quad (3) $$

The expression is forward in time and central in space, for both first and the second derivatives of velocities. This type of discretization has been discussed in module 1 through 4.

However, this equation cannot be directly solved. $\small u_{ij}^{n+1}$ obtained from this equation is not divergence free. To insure that the calculated velocities are divergence free, a multi-step method is employed. First, start with calculating a predicted velocity $u_{t_{ij}}^{{n+1}}$ using Hemholtz equation

$$ \small \dfrac {u_{t_{i,j}}^{{n}}-u_{i,j}^{n}}{\bigtriangleup t}+\dfrac {(u_{i+1/2,j}^{n})^2-(u_{i-1/2,j}^{n})^2}{\bigtriangleup x}+\dfrac {(u_{i,j+1/2}^{n})(v_{i,j+1/2}^{n})-(u_{i,j-1/2}^{n})(v_{i,j-1/2}^{n})}{\bigtriangleup y}=$$$$\small \nu \dfrac{(u_{i+1,j}^{n}+u_{i-1,j}^{n}-2 u_{i,j}^{n})}{{\bigtriangleup x}^2}+\nu \dfrac{(u_{i,j+1}^{n}+u_{i,j-1}^{n}-2 u_{i,j}^{n})}{{\bigtriangleup y}^2} \quad (4) $$

$$\small \implies \dfrac {u_{t_{i,j}}^{n}-u_{i,j}^{n}}{\bigtriangleup t}+\dfrac {(u_{e})^2-(u_{w})^2}{\bigtriangleup x}+\dfrac {(u_{n})(v_{n})-(u_{s})(v_{s})}{\bigtriangleup y}=\nu \dfrac{(u_{E}+u_{W}-2 u_{P})}{{\bigtriangleup x}^2}+\nu \dfrac{(u_{N}+u_{S}-2 u_{P})}{{\bigtriangleup y}^2} \quad (5) $$

$$ \small \implies \dfrac {u_{t_{i,j}}^{n}-u_{i,j}^{n}}{\bigtriangleup t}-F1C_{ij}=FV_{ij}\quad (6) $$

Similarly equations can be written in $\small y$-direction by replacing function $\small F1C$ by $\small F2C$



In [37]:

    
# Defining functions for N-S Solver

def F1C(ue,uw,us,un,vs,vn,dx,dy):
    F1C=-((ue**2)-(uw**2))/dx-((un*vn)-(us*vs))/dy
    return F1C

def FV(uP,uE,uW,uN,uS,dx,dy,mu):
    FV=(mu/dx)*(((uE-uP)/dx)-((uP-uW)/dx))+(mu/dy)*(((uN-uP)/dy)-((uP-uS)/dy))
    return FV

def F2C(vn,vs,ve,vw,ue,uw,dx,dy):
    F2C=-((ue*ve)-(uw*vw))/dx-((vn**2)-(vs**2))/dy
    return F2C

Once the predicted velocities are obtained, the corrected velocity is given by,

$$ \small \dfrac {u_{{i,j}}^{n+1}-u_{t_{i,j}}^{n}}{\bigtriangleup t}=-\dfrac {P_{i+1/2,j}^{n}-P_{i-1/2,j}^{n}}{\bigtriangleup x}\quad (7) $$

Taking divergence of this equation, we obtain the poisson equation for the corrected pressure in terms of only the predicted velocity, making the corrected velocity divergence free. This poisson equation is solved using a Weighted Jacobi Method, using a relaxation parameter, $\small \omega$. This method is similar to the Jacobi Method discussed in previous modules but uses $\small \omega$ to included weighted influence of $\small P_{ij}$ at time step $\small n$ to calculate it's value at $\small n+1$. More discussion on this can be found in Ref [5].

This corrected pressure evaluated from the poisson solver is plugged back into equation (7) to calculate the corrected velocity. Lets begin solving !



In [38]:

    
dt=.4*min(dx,dy) # Defining dt taking sigma=0.4

error=1.
maxit=100 # maximum iterations for pressure poisson solver
counter=0

[X,Y]=np.meshgrid(x,y)

# loop for navier-stokes solver

while(error>0.0001):
    un=u.copy()
    vn=v.copy()
    
    # Predictor Step
    ue1=(u[1:-1,1:-1]+u[2:,1:-1])/2
    uw1=(u[1:-1,1:-1]+u[:-2,1:-1])/2
    us1=(u[1:-1,1:-1]+u[1:-1,:-2])/2
    un1=(u[1:-1,1:-1]+u[1:-1,2:])/2
    vs1=(v[1:-2,:-1]+v[2:-1,:-1])/2
    vn1=(v[1:-2,1:]+v[2:-1,1:])/2
        
    G1=(F1C(ue1,uw1,us1,un1,vs1,vn1,dx,dy)+\
        FV(u[1:-1,1:-1],u[2:,1:-1],u[:-2,1:-1],u[1:-1,2:],u[1:-1,:-2],dx,dy,mu))
    
    ut[1:-1,1:-1]=u[1:-1,1:-1]+(dt/1)*G1

    
    ve12=(v[1:-1,1:-1]+v[2:,1:-1])/2
    vw12=(v[1:-1,1:-1]+v[:-2,1:-1])/2
    vs12=(v[1:-1,1:-1]+v[1:-1,:-2])/2
    vn12=(v[1:-1,1:-1]+v[1:-1,2:])/2
    us12=(u[:-1,1:-2]+u[:-1,2:-1])/2
    un12=(u[1:,1:-2]+u[1:,2:-1])/2
        
    G2=(F2C(ve12,vw12,vs12,vn12,us12,un12,dx,dy)+\
        FV(v[1:-1,1:-1],v[2:,1:-1],v[:-2,1:-1],\
        v[1:-1,2:],v[1:-1,:-2],dx,dy,mu))
    
    vt[1:-1,1:-1]=v[1:-1,1:-1]+(dt/1)*G2    
    

    # Boundary Conditions for predicted velocities
    ut[:,0]=2*0-ut[:,1] # Bottom boundary condition
    ut[:,-1]=2*U-ut[:,-2] # Top boundary condition
    
    
    vt[0,:]=0-vt[1,:] # Left boundary condition
    vt[-1,:]=0-vt[-2,:] # Right boundary condition

    
    # Pressure Poisson Solver
    
    for it in range(maxit):
        for i in range(1,Nx):
            for j in range(1,Ny):
                p[i][j]=w*0.25*(p[i][j+1]+p[i][j-1]+p[i+1][j]+p[i-1][j]-\
                        (1*dy/(dt))*(vt[i][j]-vt[i][j-1])-\
                       (1*dx/(dt))*(ut[i][j]-ut[i-1][j]))+(1-w)*p[i][j]    
        
        p[:,0]=p[:,1]
        p[:,-1]=p[:,-2]

        p[0,:]=p[1,:]
        p[-1,:]=p[-2,:]
    
    # Corrector Step
    
    u[1:-1,1:-1]=ut[1:-1,1:-1]-(dt/(1*dx))*(p[2:-1,1:-1]-p[1:-2,1:-1])
    v[1:-1,1:-1]=vt[1:-1,1:-1]-(dt/(1*dy))*(p[1:-1,2:-1]-p[1:-1,1:-2])

    
    # Boundary Conditions for Correct Velocities
    u[:,0]=2*0-u[:,1] # Bottom boundary condition
    u[:,-1]=2*U-u[:,-2] # Top boundary condition
    
    
    v[0,:]=0-v[1,:] # Left boundary condition
    v[-1,:]=0-v[-2,:] # Right boundary condition
    
    error=max(np.sqrt(((u-un)**2).mean()),np.sqrt(((v-vn)**2).mean()))
    
    counter+=1

uu=(u[:,:-1]+u[:,1:])*0.5
vv=(v[:-1,:]+v[1:,:])*0.5

uu=uu.T
vv=vv.T

The values at ghost cells, are determined by using the calculated values at the grid cells and the boundary conditions. The velocity at top boundary is equal to the free-stream velocity and all other boundaries act as a no-slip wall.

Post processing, data visualization and validation: Contours of velocity magnitude and streamlines are plotted and compared to published data.



In [39]:

    
%matplotlib inline

# Contours of velocity magnitude
plt.figure(figsize=(size, (1.-0)/(1.-0)*size))
plt.contourf(X,Y,np.sqrt(uu**2+vv**2))
plt.colorbar()
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Velocity Contours')
plt.show()


# Streamlines
plt.figure(figsize=(size, (1.-0)/(1.-0)*size))
plt.plot([0,0,1,1,0],[0,1,1,0,0],'k')
plt.streamplot(X, Y, uu, vv, density=3, arrowsize=1, arrowstyle='->')
plt.axis([0, 1, 0, 1])
plt.xlabel('x')
plt.ylabel('y')
plt.title('Velocity Streamlines')
plt.show()

print "Time=",counter*dt,"s"









    












    












    



Time= 1.95918367347 s

These results can be validated with those obtained by Ghia and Idris.et.al[2-4]



In [40]:

    
display(Image(filename='../Project_Akash/Ghia.png'))

Results obtained by Ghia.et.al [2], show the formation of primary, secondary and corner vortices at Re=100

Further accuracy can be obtained by using fracitonal-step methods like RK-3 and Adam-Bashforth time marching schemes. Further discussion on these methods can be found in [5].

4. Immersed boundary method - Direct forcing apporach

Cartesian girds are very computationally efficient and their application to parallel programming is quite straight forward [5-6]. Thus, to utilize this property of a cartesian grid, immersed boundary methods are used.

Consider a cylinder in the middle of the domain of a lid driven cavity, with a radius of 0.1 m.



In [41]:

    
r=0.1 #m # Cylinder radius
L=2*r 

dis=np.linspace(0.,pi,20)
xcircle=0.5+r*np.cos(dis)
ycircle=np.sqrt(r**2-(xcircle-0.5)**2)

xcircle=np.concatenate([xcircle,xcircle[::-1]])
ycircle=np.concatenate([ycircle,-ycircle[::-1]])

ycircle=ycircle+0.5

%matplotlib inline

plt.figure(figsize=(size, (1.-0)/(1.-0)*size))
plt.plot(X.T,Y.T,'g')
plt.plot(X,Y,'g')
plt.plot(xcircle,ycircle,'k')
plt.axis('equal')
plt.show()

There are various apporaches of an immersed boundary method. In this notebook, as discussed earlier, direct forcing apporach is used. In this approach, all the $u$ and $v$ velocity points just outside the immersed boundary and its wall velocity is imposed on those points after the predictor step.

First, all the points inside the body are determined using it's equation. The indices of grid location of velocities, if inside the circle are stored in arrays $\small iu$ and $\small iv$.



In [42]:

    
iu=[]
for i in range(len(xu)):
    for j in range(len(yu)):
        if(sqrt((xu[i]-0.5)**2+(yu[j]-0.5)**2)<=r):
            iu.append([i,j])

iu=np.asarray(iu)

iv=[]
for i in range(len(xv)):
    for j in range(len(yv)):
        if(sqrt((xv[i]-0.5)**2+(yv[j]-0.5)**2)<=r):
            iv.append([i,j])
            
iv=np.asarray(iv)

Next the points outside the solid body are determined. For this, a special function where() from numpy library is used, which returns the indices of array elements where a given condition is satisfied.



In [43]:

    
bul=min(iu[:,0])
bur=max(iu[:,0])
bub=min(iu[:,1])
but=max(iu[:,1])

# find u-velocity points just outside the boundary

cux=np.linspace(bul,bur,bur-bul+1) # array to store indices along x-direction

bu=[]

for i in range(len(cux)):
    
    ind=np.where(iu[:,0]==cux[i])
    
    Min=iu[ind[0][0]][1]
    Max=iu[ind[0][0]][1]
    
    # getting maximum and minimum indices in y-dir  for a given indice in x-dir
    
    for j in range(len(ind[0])):
        if(iu[ind[0][j]][1]<Min):
            Min=iu[ind[0][j]][1]
        if(iu[ind[0][j]][1]>Max):
            Max=iu[ind[0][j]][1]
        
    bu.append([cux[i],Min-1])
    bu.append([cux[i],Max+1])

cuy=np.linspace(bub,but,but-bub+1) # array to store indices along y-direction

for i in range(len(cuy)):
    
    ind=np.where(iu[:,1]==cuy[i])
    
    Min=iu[ind[0][0]][0]
    Max=iu[ind[0][0]][0]
    
    # getting maximum and minimum indices in x-dir  for a given indice in y-dir
    
    for j in range(len(ind[0])):
        if(iu[ind[0][j]][0]<Min):
            Min=iu[ind[0][j]][0]
        if(iu[ind[0][j]][0]>Max):
            Max=iu[ind[0][j]][0]
    
    bu.append([Min-1,cuy[i]])
    bu.append([Max+1,cuy[i]])



bu=np.asarray(bu) # array which stores indices of u-velocity points

%matplotlib inline

plt.figure(figsize=(size, (1.-0)/(1.-0)*size))
plt.plot(X.T,Y.T,'g')
plt.plot(X,Y,'g')
plt.plot(xcircle,ycircle,'k')
plt.axis('equal')


for i in range(len(bu)):
    plt.plot(xu[bu[i,0]],yu[bu[i,1]],'.b')
 

bvl=min(iv[:,0])
bvr=max(iv[:,0])
bvb=min(iv[:,1])
bvt=max(iv[:,1])

cvx=np.linspace(bvl,bvr,bvr-bvl+1) # array to store indices along x-direction

bv=[]

for i in range(len(cvx)):
    
    ind=np.where(iv[:,0]==cvx[i])
    
    Min=iv[ind[0][0]][1]
    Max=iv[ind[0][0]][1]
    
    # getting maximum and minimum indices in y-dir  for a given indice in x-dir
    
    for j in range(len(ind[0])):
        if(iv[ind[0][j]][1]<Min):
            Min=iv[ind[0][j]][1]
        if(iv[ind[0][j]][1]>Max):
            Max=iv[ind[0][j]][1]
        
    bv.append([cvx[i],Min-1])
    bv.append([cvx[i],Max+1])

cvy=np.linspace(bvb,bvt,bvt-bvb+1) # array to store indices along y-direction

for i in range(len(cvy)):
    
    ind=np.where(iv[:,1]==cvy[i])
    
    Min=iv[ind[0][0]][0]
    Max=iv[ind[0][0]][0]
    
    # getting maximum and minimum indices in x-dir  for a given indice in y-dir
    
    for j in range(len(ind[0])):
        if(iv[ind[0][j]][0]<Min):
            Min=iv[ind[0][j]][0]
        if(iv[ind[0][j]][0]>Max):
            Max=iv[ind[0][j]][0]
    
    bv.append([Min-1,cvy[i]])
    bv.append([Max+1,cvy[i]])


bv=np.asarray(bv) # array which stores indices of v-velocity points

for i in range(len(bv)):
    plt.plot(xv[bv[i,0]],yv[bv[i,1]],'xk')

plt.show()

The most complex task is done. The indices of the $\small u$ and $\small v$ velocity points just outside the solid boundary are stored in arrays $\small bu$ and $\small bv$ respectively. Next we use the direct forcing approach to force velocity at these points to be zero (stationary cylinder).

The boundary condition at top will have both and $ \small u$ and $ \small v$ velocity in this case, as described in the code segment below.



In [44]:

    
u[:,:]=0
ut[:,:]=0
v[:,:]=0
vt[:,:]=0

error=1.
counter=0 # counter for time loop

dt=.4*min(dx,dy) # Defining dt taking sigma=0.4

maxit=100
[X,Y]=np.meshgrid(x,y)

# loop for navier stokes solver

while(error>0.0001):
    un=u.copy()
    vn=v.copy()
    
    # Predictor Step
    
    ue1=(u[1:-1,1:-1]+u[2:,1:-1])/2
    uw1=(u[1:-1,1:-1]+u[:-2,1:-1])/2
    us1=(u[1:-1,1:-1]+u[1:-1,:-2])/2
    un1=(u[1:-1,1:-1]+u[1:-1,2:])/2
    vs1=(v[1:-2,:-1]+v[2:-1,:-1])/2
    vn1=(v[1:-2,1:]+v[2:-1,1:])/2
        
    G1=(F1C(ue1,uw1,us1,un1,vs1,vn1,dx,dy)+\
        FV(u[1:-1,1:-1],u[2:,1:-1],u[:-2,1:-1],\
           u[1:-1,2:],u[1:-1,:-2],dx,dy,mu))
    
    ut[1:-1,1:-1]=u[1:-1,1:-1]+(dt/1)*G1

    
    ve12=(v[1:-1,1:-1]+v[2:,1:-1])/2
    vw12=(v[1:-1,1:-1]+v[:-2,1:-1])/2
    vs12=(v[1:-1,1:-1]+v[1:-1,:-2])/2
    vn12=(v[1:-1,1:-1]+v[1:-1,2:])/2
    us12=(u[:-1,1:-2]+u[:-1,2:-1])/2
    un12=(u[1:,1:-2]+u[1:,2:-1])/2
        
    G2=(F2C(ve12,vw12,vs12,vn12,us12,un12,dx,dy)+\
        FV(v[1:-1,1:-1],v[2:,1:-1],v[:-2,1:-1],\
        v[1:-1,2:],v[1:-1,:-2],dx,dy,mu))
    
    vt[1:-1,1:-1]=v[1:-1,1:-1]+(dt/1)*G2    
    
    # Direct Forcing 
    
    for i in range(len(bu)):
       ut[bu[i,0]][bu[i,1]]=0
    
    for i in range(len(bv)):
       vt[bv[i,0]][bv[i,1]]=0
    
    # Boundary Condtions for predicted velocities
    
    ut[:,0]=2*0-ut[:,1] # Bottom boundary condition
    ut[:,-1]=2*U-ut[:,-2] # Top boundary condition
    
    
    vt[0,:]=0-vt[1,:] # Left boundary condition
    vt[-1,:]=0-vt[-2,:] # Right boundary condition
    vt[1:-1,-1]=-0.5 # Top boundary condition (v-velocity)

    # Pressure Poisson Solver
    
    for it in range(maxit):
        for i in range(1,Nx):
            for j in range(1,Ny):
                p[i][j]=w*0.25*(p[i][j+1]+p[i][j-1]+p[i+1][j]+p[i-1][j]-\
                                (1*dy/(1*dt))*(vt[i][j]-vt[i][j-1])-\
                       (1*dx/(1*dt))*(ut[i][j]-ut[i-1][j]))+(1-w)*p[i][j]    
        
        p[:,0]=p[:,1]
        p[:,-1]=p[:,-2]

        p[0,:]=p[1,:]
        p[-1,:]=p[-2,:]

    # Corrector Step
    u[1:-1,1:-1]=ut[1:-1,1:-1]-(dt/(1*dx))*(p[2:-1,1:-1]-p[1:-2,1:-1])
    v[1:-1,1:-1]=vt[1:-1,1:-1]-(dt/(1*dy))*(p[1:-1,2:-1]-p[1:-1,1:-2])

    
    # Boundary Condtions for Corrected Velocities
    u[:,0]=2*0-u[:,1] # Bottom boundary condition
    u[:,-1]=2*U-u[:,-2] # Top boundary condition
    
    
    v[0,:]=0-v[1,:] # Left boundary condition
    v[-1,:]=0-v[-2,:] # Right boundary condition
    v[1:-1,-1]=-0.5 # Top boundary condition (v-velocity)
    
    error=max(np.sqrt(((u-un)**2).mean()),np.sqrt(((v-vn)**2).mean()))
    counter+=1
    
    
uu=(u[:,:-1]+u[:,1:])*0.5
vv=(v[:-1,:]+v[1:,:])*0.5

uu=uu.T
vv=vv.T

The results below show how the velocity field curves around the solid obstacle



In [45]:

    
%matplotlib inline
# Contours
plt.figure(figsize=(size+3, (1.-0)/(1.-0)*(size+3)))
plt.contourf(X,Y,np.sqrt(uu**2+vv**2))
plt.colorbar()
plt.quiver(X,Y,uu,vv,units='width')
plt.fill(xcircle,ycircle,'k')
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Velocity Contours')
plt.show()

print "Time=",counter*dt,"s"









    












    



Time= 1.97551020408 s

5. Conclusion

The development of Navier-Stokes solver described in this notebook is just a basic formulation. More complex flows can be simulated by introducing mathematical models to account for these complexities. Turbulence and thermodynamic models have been developed to visualize the eddies formed at high Reynolds numbers. To simulate different types of flow, different boundary conditions can be imposed. Boundary conditions play an important role and are the heart of a CFD code. I encourage the reader to build on this notebook and keep working on developing complex codes.

6. References

Anderson. J.D., "A History of Aerodynamics: And Its Impact on Flying Machines",1997
Ghia.et.al, "High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method.", Journal of Computational Physics, 1982
Idris.M.S, Azwadi.C.S.N, "Numerical Investigation of 2D Lid Driven Cavity Flow Emphasizing Finite Difference Method of Non Uniform Meshing
Vaidehi Ambatipudi., "Simple Solve for Driven Cavity Flow Problem.", Purdue University
Ferziger and Peric., "Computational Methods for Fluid Dynamics.", 2001
Balaras.E., "Modelling complex boundaries using an external force field on fixed cartesian girds in large eddy simulations." , Computers & Fluids



In [46]:

    
from IPython.core.display import HTML
css_file = '../../numericalmoocstyle.css'
HTML(open(css_file, "r").read())









    Out[46]: