Custom CFD Codes

# CFD Code 1 - This code runs on the staggered grid finite volume method (FVM)

# and uses explicit time (calculates next step from previous). The

# cells are indexed from bottom-left to top-right. The boundary

# conditions include velocity inlet, velocity outlet, and walls; much

# like a 2D wind-tunnel. Inside of the domain, their is a region where

# the velocities are zero, this is meant to act like a flat plate

# although this needs to be improved. The next code should be written

# using more functions to clean up the code and make it more general.

# The next code also should be written utilizing GPU compute for the

# appropriate functions.

import numpy as np

import matplotlib.pyplot as plt

import timeit

import scipy

import scipy.linalg

import scipy.sparse

import scipy.sparse.linalg

from numba import jit

start_time = timeit.default_timer()

def pressure_poisson2(p, b, dx, dy):

pn = np.empty_like(p)

it = 0

err = 1e5

tol = 1e-8

maxit = 100

beta = 1.1

while err > tol and it < maxit:

pn = p.copy()

for i in range(1, nx+1):

for j in range(1, ny+1):

ap = Ap[j, i]

an = An[j, i]

aso = As[j, i]

ae = Ae[j, i]

aw = Aw[j, i]

rhs = b[j, i] - 1.0*(ae*p[j, i+1] + aw*p[j, i-1] + an*p[j+1, i] + aso*p[j-1, i])

p[j, i] = beta*rhs/ap + (1-beta)*p[j, i]

err = np.linalg.norm(p - pn)

it += 1

return p, err

# ======================================================================================= #

# Input Parameters

AR = 4 # Aspect Ratio of the wind tunnel (length/height)

velocity = 10.0 # wind tunnel velocity

nsteps = 100 # number of time steps

nx = 400 # number of cells horizontally

ny = nx//AR # number of cells vertically

h, w = ny//5, nx//50 # height and width of the flat plate expressed as a percentages

lx = 2.0 # length of the wind tunnel

ly = lx/AR # height of the wind tunnel

nu = 0.05 # dynamic viscosity

# ======================================================================================= #

# define some boiler plate

dx, dy = lx/nx, ly/ny # step size of the finite volumes

t = 0.0 # time starts at t0 = 0

# cell centered coordinates

xx = np.linspace(dx/2.0, lx - dx/2.0, nx, endpoint=True)

yy = np.linspace(dy/2.0, ly - dy/2.0, ny, endpoint=True)

xcc, ycc = np.meshgrid(xx, yy)

# x-staggered coordinates

xxs = np.linspace(0, lx, nx+1, endpoint=True)

xu, yu = np.meshgrid(xxs, yy)

# y-staggered coordinates

yys = np.linspace(0, ly, ny+1, endpoint=True)

xv, yv = np.meshgrid(xx, yys)

Ut, Vt = velocity, 0.0 # velocity at top wall

Ub, Vb = 0.0, 0.0 # velocity at bottom wall

Ul, Vl = 0.0, 0.0 # velocity at left wall

Ur, Vr = 0.0, 0.0 # velocity at right wall

print('Reynolds Number:', Ut*lx/nu)

dt = round(0.5 * min(0.25*dx*dx/nu, 0.25*dy*dy/nu, 4.0*nu/Ut/Ut), 5)

print('time-step, dt =', dt, "[sec]")

# ======================================================================================= #

# initialize velocities

u = np.zeros([ny+2, nx+2]) # include ghost cells

v = np.zeros([ny+2, nx+2]) # include ghost cells

u[2:-2, 1:-1] = velocity # initial condition velocity along tunnel

u[ny//2 - h//2:ny//2 + h//2, nx//2-w:nx//2] = 0.0 # velocities are zero inside Flat plate

v[ny//2 - h//2:ny//2 + h//2, nx//2-w:nx//2] = 0.0

ut = np.zeros_like(u)

vt = np.zeros_like(u)

# ======================================================================================= #

# initialize the pressure

p = np.zeros([nx+2, ny+2]) # include ghost cells

# build pressure coefficient matrix

Aw = 1.0/dx/dx*np.ones([ny, nx])

Ae = 1.0/dx/dx*np.ones([ny, nx])

As = 1.0/dy/dy*np.ones([ny, nx])

An = 1.0/dy/dy*np.ones([ny, nx])

Aw[:, 0] = 0.0 # set left wall coefficients

Ae[:, -1] = 0.0 # set right wall coefficients

An[-1, :] = 0.0 # set top wall coefficients

As[0, :] = 0.0 # set bottom wall coefficients

Ap = -(Aw + Ae + An + As)

n = nx*ny

d0 = Ap.reshape(n)

de = Ae.reshape(n)[:-1]

dw = Aw.reshape(n)[1:]

ds = As.reshape(n)[nx:]

dn = An.reshape(n)[:-nx]

A1 = scipy.sparse.diags([d0, de, dw, dn, ds], [0, 1, -1, nx, -nx], format='csr')

# ======================================================================================= #

for n in range(0, nsteps):

u[-1, 1:] = 0 # top wall

u[0, 1:] = 0.0 # bottom wall

u[1:-1, 1] = 0.0 # left wall

u[2:-2, 1] = velocity # left wall inflow

u[1:-1, -1] = 0.0 # right wall

u[2:-2, -1] = velocity # right wall outflow

v[-1, 1:-1] = 0.0 # top wall

v[1, 1:-1] = 0.0 # bottom wall

v[1:, 0] = 0.0 # left wall

v[1:, -1] = 0.0 # right wall

u[ny//2 - h//2:ny//2 + h//2, nx//2-w:nx//2] = 0.0 # velocities are zero inside Flat plate

v[ny//2 - h//2:ny//2 + h//2, nx//2-w:nx//2] = 0.0

# do x-momentum first - u is of size (nx + 2) x (ny + 2) - only need to do the interior points

for i in range(2, nx + 1):

for j in range(1, ny + 1):

ue = 0.5 * (u[j, i + 1] + u[j, i])

uw = 0.5 * (u[j, i] + u[j, i - 1])

un = 0.5 * (u[j + 1, i] + u[j, i])

us = 0.5 * (u[j, i] + u[j - 1, i])

vn = 0.5 * (v[j + 1, i] + v[j + 1, i - 1])

vs = 0.5 * (v[j, i] + v[j, i - 1])

# convection = - d(uu)/dx - d(vu)/dy

convection = - (ue * ue - uw * uw) / dx - (un * vn - us * vs) / dy

# diffusion = d2u/dx2 + d2u/dy2

diffusion = nu * ((u[j, i + 1] - 2.0 * u[j, i] + u[j, i - 1]) / dx / dx + (u[j + 1, i] - 2.0 * u[j, i] + u[j - 1, i]) / dy / dy)

ut[j, i] = u[j, i] + dt * (convection + diffusion)

# do y-momentum - only need to do interior points

for i in range(1, nx + 1):

for j in range(2, ny + 1):

ve = 0.5 * (v[j, i + 1] + v[j, i])

vw = 0.5 * (v[j, i] + v[j, i - 1])

ue = 0.5 * (u[j, i + 1] + u[j - 1, i + 1])

uw = 0.5 * (u[j, i] + u[j - 1, i])

vn = 0.5 * (v[j + 1, i] + v[j, i])

vs = 0.5 * (v[j, i] + v[j - 1, i])

# convection = d(uv)/dx + d(vv)/dy

convection = - (ue * ve - uw * vw) / dx - (vn * vn - vs * vs) / dy

# diffusion = d2u/dx2 + d2u/dy2

diffusion = nu * ((v[j, i + 1] - 2.0 * v[j, i] + v[j, i - 1]) / dx / dx + (v[j + 1, i] - 2.0 * v[j, i] + v[j - 1, i]) / dy / dy)

vt[j, i] = v[j, i] + dt * (convection + diffusion)

# we will only need to fill the interior points. This size is for convenient indexing

divut = np.zeros([ny + 2, nx + 2])

divut[1:-1, 1:-1] = (ut[1:-1, 2:] - ut[1:-1, 1:-1]) / dx + (vt[2:, 1:-1] - vt[1:-1, 1:-1]) / dy

prhs = 1.0 / dt * divut

pt, info = scipy.sparse.linalg.bicg(A1, prhs[1:-1, 1:-1].ravel(), tol=1e-10)

p = np.zeros([ny + 2, nx + 2])

p[1:-1, 1:-1] = pt.reshape([ny, nx])

# time advance

u[1:-1, 2:-1] = ut[1:-1, 2:-1] - dt * (p[1:-1, 2:-1] - p[1:-1, 1:-2]) / dx

v[2:-1, 1:-1] = vt[2:-1, 1:-1] - dt * (p[2:-1, 1:-1] - p[1:-2, 1:-1]) / dy

t += dt

# Display Calculation Status

if round((100*t/(dt*nsteps)), 1) % 55 == 0:

print()

if round((100*t/(dt*nsteps)), 1) % 5 == 0:

print(" -> ", round(100*t/(dt*nsteps)), "%", end="", sep="")

print("\nTotal Time:", t, "[sec]")

# ======================================================================================= #

# Plotting - plot the streamlines and contours of the final velocity field

ucc = 0.5*(u[1:-1, 2:] + u[1:-1, 1:-1])

vcc = 0.5*(v[2:, 1:-1] + v[1:-1, 1:-1])

speed = np.sqrt(ucc*ucc + vcc*vcc)

x = np.linspace(0, lx, nx)

y = np.linspace(0, ly, ny)

xx, yy = np.meshgrid(x, y)

nn = 1

fig1 = plt.figure('1')

plt.contourf(speed[5:-5, 5:-5], levels=100)

plt.colorbar()

fig2 = plt.figure('2')

plt.xlim([xx[0, 0], xx[0, -1]])

plt.ylim([yy[0, 0], yy[-1, 0]])

plt.streamplot(xx, yy, ucc, vcc, color=np.sqrt(ucc*ucc + vcc*vcc), density=10, cmap='plasma', linewidth=0.25)

end_time = timeit.default_timer()

print("Calculation Time:", end_time-start_time, "[sec]")

plt.show()

The following results where obtained using :

1. aspect ratio = 1

2. grid size = 512x512

3. number of steps = 100

4. velocity = 5 m/s