1/19/2017

CFD code implementation

Hi every body, This is my trial to implement CFD code for a cavity using Python, numpy and matplotlib to solve continuity equation, momentum equation and divergence. you can find the code below. 





Here is the code , Python 2.7 and Ipython notebook have been used. 

import numpy as np
import matplotlib.pyplot as plt

g = 50

u = np.zeros((g, g+1))
un = np.zeros((g, g+1))
uc = np.zeros((g, g))

v = np.zeros((g+1, g))
vn = np.zeros((g+1, g))
vc = np.zeros((g, g))

p = np.zeros((g+1, g+1))
pn = np.zeros((g+1, g+1))
pc = np.zeros((g, g))

m = np.zeros((g+1, g+1))

dx = 1.0/(g-1)
dy = 1.0/(g-1)

dt = 0.001
delta = 4.5
error = 1.0
Re = 100.

errors = []



step = 1
step2 = 1

#initializing u, v, p
for i in range(g+1):
    u[g-2][i] = 1.0
    u[g-1][i] = 1.0

for i in range(g+1):
    for j in range(g+1):
        p[i][j] = 1.
       
       
#solving equations
while(error > 0.9):
    for i in range(1, g-1):
        for j in range(1, g):
            un[i][j]= u[i][j] - dt*(  (u[i+1][j]*u[i+1][j]-u[i-1][j]*u[i-1][j])/2.0/dx +0.25*( (u[i][j]+u[i][j+1])*(v[i][j]+v[i+1][j])-(u[i][j]+u[i][j-1])*(v[i+1][j-1]+v[i][j-1]) )/dy)- dt/dx*(p[i+1][j]-p[i][j]) + dt*1.0/Re*( (u[i+1][j]-2.0*u[i][j]+u[i-1][j])/dx/dx +(u[i][j+1]-2.0*u[i][j]+u[i][j-1])/dy/dy)
       

    for j in range(1, g):
        un[0][j] = 0.0
        un[g-1][j] = 0.0

    for i in range(g):
        un[i][0] = -un[i][1]
        un[i][g] = 2-un[i][g-1]
   
# Solve v - Momentum
    for i in range(1, g):
        for j in range(1, g-1):
            vn[i][j] = v[i][j] - dt* ( 0.25*( (u[i][j]+u[i][j+1])*(v[i][j]+v[i+1][j])-(u[i-1][j]+u[i-1][j+1])*(v[i][j]+v[i-1][j]) )/dx +(v[i][j+1]*v[i][j+1]-v[i][j-1]*v[i][j-1])/2.0/dy ) - dt/dy*(p[i][j+1]-p[i][j]) + dt*1.0/Re*( (v[i+1][j]-2.0*v[i][j]+v[i-1][j])/dx/dx+(v[i][j+1]-2.0*v[i][j]+v[i][j-1])/dy/dy )

#v - boundary conditions.
    for j in range(1, g-1):
        vn[0][j] = -vn[1][j]
        vn[g][j]= -vn[g-1][j]

    for i in range(g+1):
        vn[i][0] = 0.0
        vn[i][g-1] = 0.0
   
#Solve Continuity equation

    for i in range(1, g):
        for j in range(1, g):
            pn[i][j] = p[i][j]-dt*delta*(  ( un[i][j]-un[i-1][j] )/dx + ( vn[i][j]-vn[i][j-1] ) /dy  )

# Pressure boundary conditions
    for i in range(1, g):
        pn[i][0] = pn[i][1]
        pn[i][g] = pn[i][g-1]


    for j in range(g+1):
        pn[0][j] = pn[1][j]
        pn[g][j] = pn[g-1][j]

#Display Errors
    error = 0.0

    for i in range(1, g):
        for j in range(1, g):
            m[i][j] = (( un[i][j]-un[i-1][j] )/dx + ( vn[i][j]-vn[i][j-1] )/dy  )
            error = error + abs(m[i][j])
    if(step%5 == 1):
        errors.append(error)
       
    if(step%500 == 1):
       
        print str(step2) + "  - Error: " + str(error) + " for step : " + str(step)
        step2 +=1
        
#iterating u
    for i in range(g):
        for j in range(g+1):
            u[i][j] = un[i][j]
   
#iterating v
    for i in range(g+1):
        for j in range(g):
            v[i][j] = vn[i][j]
   
#iterating p
    for i in range(g+1):
        for j in range(g+1):
            p[i][j] = pn[i][j]

    step += 1

f, ((plt1,plt2),(plt3, plt4)) = plt.subplots(2, 2)

f = plt1.pcolor(u.transpose())
f = plt2.pcolor(v.transpose())
f = plt3.contourf(p.transpose())
f = plt4.plot(errors)
plt4.invert_yaxis()

plt1.set_title("u- velocity contour")
plt2.set_title("v- velocity contour")
plt3.set_title("pressure contour")
plt4.set_title("divergence ..")
plt.show()
f, ((plt1,plt2),(plt3, plt4)) = plt.subplots(2, 2)

f = plt1.pcolor(u.transpose(), shading = "flat")
f = plt2.pcolor(v.transpose(), shading = "flat")
f = plt3.pcolor(p.transpose(), shading = "flat")
f = plt4.plot(errors)

plt4.invert_yaxis()

plt1.set_title("u- velocity contour")
plt2.set_title("v- velocity contour")
plt3.set_title("pressure contour")
plt4.set_title("divergence ..")
plt.show()

ليست هناك تعليقات :

إرسال تعليق