In [1]:
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 matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
In this file we want to compare our code with the example which we solved in the class for Module 4.
Do you remember that example? Please review the example and then come back here.
It is expected you not to have any problems with applying Dirichlet boundary for temperature because we have already done it for five examples. You just need to fix your parameters.
In [2]:
nx=80 # the number of nodes in lattice direction
dt=1
dx=1
x1=0.
x2=1.0 # the length
alpha=1.22e-3 # heat diffusion coefficient
D=1 # the dimension of the problem
mstep=40000 # the number of time step
omega=1./(0.5+(alpha*D)) #Chapman-Enskog dt=1 and dx=1
Tleft=100.0 # left wall temperature
Tright=0.0 # right wall temperature
k=3 # k=0,1,2 <==== c(2)===c(0)====> c(1)
In [3]:
x=numpy.linspace(x1,x2,nx+1)
w=numpy.zeros(k) # witghting factor
T=numpy.zeros(nx+1) # Temperature matrix
f= numpy.zeros((k, nx+1)) # distribution function
feq= numpy.zeros((k, nx+1)) # Equilibrum distribution function
In [4]:
w[0]=0.0
w[1]=0.5
w[2]=0.5 # k=3
In [5]:
#================== Initial value
T[:]=0.0 #initial temperature
T[0]=100.0
for i in range (k): #k=0,1,2
f[i,:]=w[i]*T[:]
In [6]:
# Main loop : comprised two parts :collision and streaming
#=====================
for n in range(mstep) :
# collision process
for i in range(nx+1):
T[i]=f[0,i]+f[1,i]+f[2,i]
feq[0:k,i]=w[0:k]*T[i]
f[0:k,i]=(1.-omega)*f[0:k,i]+omega*feq[0:k,i]
# streaming process
for i in range(0,nx):
f[1,nx-i]=f[1,nx-i-1]
f[2,i]=f[2,i+1]
# Boundary condition
f[1,0]=Tleft-f[2,0] # constant temperature at left T=1,0
f[2,nx]=Tright-f[1,nx] # constant temperature at right T=0.0
T[0]=100.0
T[nx]=0.0
# end of the main loop
In [7]:
# Finite difference #
Tf=numpy.zeros(nx+1) # Temperature finite difference
To=numpy.zeros(nx+1) # Temperature storaage for previous time
Tf[0]=100.0
Tf[nx]=0.0
nx=50
dx=1./nx
dt=0.163
nt=100 # 2*alpha*dt/dx^2 << 1 === > for stability
for n in range(nt) :
To[:]=Tf[:]
for i in range(1,nx):
Tf[i]=To[i]+ ( (dt*alpha/(dx**2.)) *(To[i+1]-2*To[i]+To[i-1]) )
In [8]:
plt.figure(figsize=(15,6), dpi=100)
plt.xlabel('x', fontsize=10) #x label
plt.ylabel('T', fontsize=10) #y label
plt.plot(x,T, 'r-',label=' Lattice Boltzmann Method');
plt.plot(x,Tf, 'bo',label=' Finite Difference Method ');
plt.legend();