In [148]:
import numpy
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import *
Before going to module 4 of our Numerical class, we want to check the example to completely sure about our code. In this code we have just changed our boundary conditions as follow:
Tleft temperature: hot
Tright temperature: cold
Ttop temperature: cold
TBottom wall: adiabatic or insulated
In [149]:
nx=100 # the number of nodes in x direction lattice direction
ny=100 # the number of nodes in y direction lattice direction
u=0.1
v=0.4
alpha=1.0 # heat diffusion coefficient
D=3.0 # the dimension of the problem
mstep=200 # the number of time step
omega=1./(0.5+(alpha*D)) #Chapman-Enskog dt=1 and dx=1
Tleft=1.0 # left wall temperature
Tright=0.0 # right wall temperature
Ttop=0.0
k=9 # k=0,1,2,3,4,5,6,8,9
In [150]:
x=numpy.linspace(0,1,nx+1)*nx
y=numpy.linspace(0,1,ny+1)*ny
w=numpy.ones(k) # witghting factor
cx=numpy.ones(k)
cy=numpy.ones(k)
T=numpy.ones((ny+1,nx+1) ) # Temperature matrix
f= numpy.ones((k, ny+1,nx+1)) # distribution function
In [151]:
w[0]=4./9 #4./9
w[1:5]=1./9.
w[5:9]=1./36.
cx[0]=0.0
cx[1]=1.0
cx[2]=0.0
cx[3]=-1.0
cx[4]=0.0
cx[5]=1.0
cx[6]=-1.0
cx[7]=-1.0
cx[8]=1.0
cy[0]=0.0
cy[1]=0.0
cy[2]=1.0
cy[3]=0.0
cy[4]=-1.0
cy[5]=1.0
cy[6]=1.0
cy[7]=-1.0
cy[8]=-1.0
In [152]:
T[0:ny+1,0:nx+1]=0.0
T[0:ny+1,0]=1.0 # Left T=1
T[0:ny+1,nx]=0.0 # right T=0.0
T[ny,0:nx+1]=0.0 # Top T=0.0
for i in range(nx+1):
for j in range(ny+1):
for l in range (k): #k=0,1,2,3,4
f[l,j,i]=w[l]*T[j,i]
Main loop : comprised two parts :collision and streaming
In [153]:
for n in range(mstep) :
# collision process
for i in range(nx+1):
for j in range(ny+1):
for l in range (k):
feq=w[l]*T[j,i]*( 1.+ 3.*(cx[l]*u+cy[l]*v ) )
f[l,j,i]=(1.-omega)*f[l,j,i]+omega*feq
#streaming process
# ==========================
for i in range(nx):
for j in range(ny,0,-1):
f[6,j,i]=f[6,j-1,i+1]
f[2,j,i]=f[2,j-1,i]
for i in range(nx,0,-1):
for j in range(ny,0,-1):
f[5,j,i]=f[5,j-1,i-1]
f[1,j,i]=f[1,j,i-1]
for i in range(nx,0,-1):
for j in range(0,ny):
f[8,j,i]=f[8,j+1,i-1]
f[4,j,i]=f[4,j+1,i]
for i in range(0,nx):
for j in range(0,ny):
f[7,j,i]=f[7,j+1,i+1]
f[3,j,i]=f[3,j,i+1]
# Boundary condition
# =============================
for j in range(0,ny+1) : #left Boundary
f[1,j,0]=( Tleft*(w[1]+w[3]) )-f[3,j,0]
f[5,j,0]=( Tleft*(w[5]+w[7]) )-f[7,j,0]
f[8,j,0]=( Tleft*(w[8]+w[6]) )-f[6,j,0]
for j in range(0,ny+1): #right Boundary
f[3,j,nx]=( Tright*(w[1]+w[3]) )-f[1,j,nx]
f[7,j,nx]=( Tright*(w[5]+w[7]) )-f[5,j,nx]
f[6,j,nx]=( Tright*(w[8]+w[6]) )-f[8,j,nx]
for i in range(0,nx+1): # top Boundary
f[4,ny,0:nx+1]=( Ttop*(w[4]+w[2]) )-f[2,ny,0:nx+1]
f[7,ny,0:nx+1]=( Ttop*(w[7]+w[5]) )-f[5,ny,0:nx+1]
f[8,ny,0:nx+1]=( Ttop*(w[8]+w[6]) )-f[6,ny,0:nx+1]
for i in range(1,nx): # bottom Boundary
for l in range(k):
f[l,0,i]=f[l,1,i]
#================================
for i in range(nx+1):
for j in range(ny+1):
sum=0.0
for l in range (k):
sum=sum+f[l,j,i]
T[j,i]=sum
T[0:ny+1,0]=Tleft
T[0:ny+1,nx]=Tright
T[0,1:nx]=T[1,1:nx] # adiabatic boundary
T[ny,1:nx]=Ttop
#==============================
#==============================
Finite Difference Method
In [154]:
Tf=numpy.ones((ny+1,nx+1) ) # finite difference Temperature matrix
To=numpy.ones((ny+1,nx+1) )
mstep=1000
dx=1.0
dy=1.0
alpha=1.0
dt=0.2
Tf[0:nx+1,0:ny+1]=0.0
Tf[0:ny+1,0]=1.0
Tf[0:ny+1,nx]=0.0
In [155]:
for n in range (mstep):
To[0:ny+1,0:nx+1]=Tf[0:ny+1,0:nx+1]
for j in range (1,ny):
for i in range (1,nx):
Txx=(To[j,i+1]+To[j,i-1])
Tyy=(To[j+1,i]+To[j-1,i])
advx=dt*u*(To[j,i]-To[j,i-1])/dx
advy=dt*v*(To[j,i]-To[j-1,i])/dy
advt=advx+advy
Tf[j,i]=To[j,i]+ (dt*alpha*(Txx+Tyy-4.*To[j,i])/(dx**2))-advt
#
Tf[0:ny+1,0]=Tleft
Tf[0:ny+1,nx]=Tright
Tf[0,1:nx]=Tf[1,1:nx] # adiabatic boundary
Tf[ny,1:nx]=Ttop
##
In [156]:
CS1 = plt.contour(x,y,Tf)
title('Finite Difference Method')
clabel(CS1,inline='True')
plt.xlabel('$x$')
plt.ylabel('$y$')
show()
In [157]:
CS2 = plt.contour(x,y,T)
title('Lattice Boltzmann Method')
clabel(CS2,inline='True')
plt.xlabel('$x$')
plt.ylabel('$y$')
show()
In [158]:
plt.figure(figsize=(7,10), dpi=100)
T1=T[ny/2,0:nx+1]
T2=Tf[ny/2,0:nx+1]
plt.xlim(0,100)
plt.ylim(-0.05,1.05)
plt.plot(x,T1, 'r-');
plt.plot(x,T2, 'bo');
plt.legend(['LBM','FDM']);
plt.xlabel('$x$')
plt.ylabel('$y$')
show()
In [158]: