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 pylab import *


    

In [2]:

    

nx=50 # the number of nodes in x direction lattice direction 
ny=50 # the number of nodes in y direction lattice direction 


alpha=0.25*ny # heat diffusion coefficient 
D=3.0 # the dimension of the problem 
mstep=2000 # 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
k=9 # k=0,1,2,3,4,5,6,8,9


    

In [3]:

    

x=numpy.linspace(0,1,nx+1)
y=numpy.linspace(0,1,ny+1)
w=numpy.ones(k)      # witghting factor
T=numpy.ones((ny+1,nx+1) )   # Temperature matrix
f= numpy.ones((k, ny+1,nx+1))   # distribution function


    

In [4]:

    

##================ Initial boundary condition
w[0]=4./9 #4./9
w[1:5]=1./9.
w[5:9]=1./36.


    

In [5]:

    

##================== Initial value

T[0:ny+1,0:nx+1]=0.0
T[0:ny+1,0]=1.0
T[0:ny+1,nx]=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]


    

In [6]:

    

##   Main loop  : comprised two parts :collision and streaming
##=====================
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]    
    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(1,nx):                # bottom and top Boundary
  for l in range(k):
   f[l,0,i]=f[l,1,i]
   f[l,ny,i]=f[l,ny-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]
 T[ny,1:nx]=T[ny-1,1:nx]
#==============================


    

In [7]:

    

#  Finite Difference Method 

#=============================
Tf=numpy.ones((ny+1,nx+1) )   # finite difference Temperature matrix
To=numpy.ones((ny+1,nx+1) ) 

mstep=20000
dx=1./nx
dy=1./ny
alpha=0.25
dt=1e-4


Tf[0:nx+1,0:ny+1]=0.0
Tf[0:ny+1,0]=1.0
Tf[0:ny+1,nx]=0.0


for n in range (mstep):
 To[0:nx+1,0:ny+1]=Tf[0:nx+1,0:ny+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])
   Tf[j,i]=To[j,i]+ (dt*alpha*(Txx+Tyy-4.*To[j,i])/(dx**2))
#
 Tf[0,1:nx]=Tf[1,1:nx]
 Tf[ny,1:nx]=Tf[ny-1,1:nx]  
##


    

In [8]:

    

CS1 = plt.contour(x,y,Tf)
title('Finite Difference Method')
clabel(CS1,inline='True')
plt.xlabel('$x$')
plt.ylabel('$y$')
show()


    

In [9]:

    

CS2 = plt.contour(x,y,T)
title('Lattice Boltzmann Method')
clabel(CS2,inline='True')
plt.xlabel('$x$')
plt.ylabel('$y$')
show()


    

In [10]:

    

T1=T[ny/2,0:nx+1]
T2=Tf[ny/2,0:nx+1]
plt.xlim(-0.001,1.1)
plt.ylim(-0.001,1.1)
plt.plot(x,T1, 'r-'); 
plt.plot(x,T2, 'bo'); 
plt.legend(['LBM','FDM']); 
plt.xlabel('$x$')
plt.ylabel('$y$')
show()