In [25]:
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
One Dimensional Heat Conduction
Considering the plate with hot wall on the left boundary and the cold wall on the right, we need to solve a Navier-stokes equation which has only diffusion part. The procedures can be categorized as follow:
1- create dimensions and Lattice nodes and physical parameters
2- calculate collision function
3- streaming the distribution function after particles' collision
4- imply proper boundary condition on the walls for temperature
5- obtain temperature based on summing temperature distribution function
6- repeat the procedure until the criteria of convergence is satisfied
Note:
1- We want to solve one dimensional equation therefore it is reasonable to use D1Q3. It means the lattice velocities are placed on the left and right hand side of the center node . Each point has its weighting coefficient.
2- As we discussed in the introduction summing distribution function $f$ reveals the accumulation or probability of particles defined as density. You know that density is a scalar parameter. For obtaining temperature distribution or any transport phenomena we need to simulate the scalars with defining distribution function. This is the strategy which will be implemented in each problem. However: The equilibrium distribution function for Advection-diffusion is different from diffusion. In fact equilibrium function will be modified by some parameters to involves the effect of velocities to have advection effects ( we will talk about this issue in Lid-driven and natural convection Example )
3- In our numerical class we learned Finite difference method so here we validate our code with making comparison with FDM. For FDM we consider forward time step and central differences for diffusion part. We know that central is second order accuracy and creates stability solution for diffusion problem but not for advection one.
Important Points in LBM:
Proper unit conversions are required to model problems with LBM. Always find all independent dimensionless numbers for given problem before anything else is done. At first we need to convert equations into non-dimensional format then some physical non-dimensional parameters will appear. We use these parameters to make compatible and identical physical conditions between LBM and original format. LBM's equation is similar to Navier-stokes because both of them are transport equations with diffusion and advection parts.
Question 1: What should we do when we have source terms? Answer later because in this example we do not have source term.
Question 2: How can we deal with the boundaries which have Flux? If you remember we learned in Numerical class how can we obtain proper boundary condition when we have flux? The same procedure is used here by considering the fact that the scalar value in each node is equal to the sum of distribution function (which itself has two directions in 1D or 8 directions in 2D).
For our problem now we need to know how to apply constant value for temperature or insulated boundary. What form the distribution function would take? Considering the left boundary we have:
$$\begin{array}{l} {f_{1}^{eq} (0,t)-f_{1}^{} (0,t)+f_{2}^{eq} (0,t)-f_{2}^{} (0,t)=0} \\ {f_{1}^{eq} (0,t)=\omega _{1} .T(0,t)=\omega _{1} .T_{w} } \\ {f_{2}^{eq} (0,t)=\omega _{2} .T(0,t)=\omega _{2} .T_{w} } \\ {x=0\, \, :\, \, \, \, f_{1} (0)+f_{2} (0)=T_{w} } \\ {STREAMING\, :\, \, f_{2} (0)=f_{2} (1)} \\ {f_{1} (0)=T_{w} -f_{2} (0)} \end{array}$$For having insulated boundary the distribution function in one direction needs to be damped by the distribution function of opposite direction right on the boundary: for example if the left wall was adiabatic then we had :
$$\begin{array}{l} {f_{1} (0)=f_{3} (0)} \\ {f_{5} (0)=f_{7} (0)} \\ {f_{8} (0)=f_{6} (0)} \end{array}$$Two dimensional heat conduction equation:
$\frac{\partial T}{\partial t} =\alpha \, \, [\, \, \frac{\partial }{\partial x} (\frac{\partial T}{\partial x} )+\frac{\partial }{\partial y} (\frac{\partial T}{\partial y} )\, \, ]$$\frac{\partial T}{\partial t} =\alpha \, \, [\, \, \frac{\partial }{\partial x} (\frac{\partial T}{\partial x} )+\frac{\partial }{\partial y} (\frac{\partial T}{\partial y} )\, \, ]$
Connection with FD and the Navier-Stokes Equations:
The Lattice Boltzmann method can be viewed as a special finite difference method .The primary reason why LBM can serve as a method for fluid simulations is that the Navier-Stokes equations can be recovered from the discrete equations through the "Chapman-Enskog" procedure, a multi-scaling expansion technique .We are familiar with this equation thanks to the heat transfer course. As we discussed for modelling we need to convert the equation to LBM format and find about the relation between diffusion coefficients (in this example$\alpha $) or as another example we can refer to the kinematic viscosity in fluid problems and relaxation time in LBM. For achieving this goal we need to follow ``Chapman-Enskog'' method. Please review this section from original textbook. By applying Chapmann -Enskog we found that: (For D1Q3 D=2 and for D2Q9 it takes D=3) $$\alpha =\frac{\Delta x^{2} }{D\Delta t} (\frac{1}{\omega } -\frac{1}{2} )$$
In [26]:
nx=100 # the number of nodes in lattice direction
dt=1 # Lattice spacing between two particle in x direction
dx=1 # Time interval between two step
ck=dx/dt # Lattice speed
csq=ck*ck
u=0.1
x1=0. # position left wall
x2=1.0 # position right wall (the length)
alpha=0.25 # Heat diffusion coefficient
D=1 # The dimension of the problem
mstep=int (0.04*nx**2) # The number of time step
omega=1./(0.5+(alpha/(dt*csq))) # Chapman-Enskog dt=1 and dx=1
Tleft=1.0 # left wall temperature
Tright=0.0 # right wall temperature
k=3 # k=0,1,2 <==== c(2)===c(0)====> c(1)
In [27]:
xl=numpy.linspace(x1,x2,nx+1)
xl=xl*100
w=numpy.zeros(k) # weighting factor
cx=numpy.zeros(k) # lattice speeds
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 [28]:
w[0]=0.0 # weighting coefficients
w[1]=0.5
w[2]=0.5 # k=3
cx[0]=0.0 # node speed coefficients
cx[1]=1.0
cx[2]=-1.0
In [29]:
T[:]=0.0 #initial temperature
T[0]=1.0
for i in range (k): #initial distribution function
f[i,:]=w[i]*T[:]
Main loop : comprised two parts :collision and streaming
In [30]:
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]*( 1.+(cx[0:k]*u) )
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]=1.0
T[nx]=0.0
# end of the main loop
Finite difference # Advection-Diffusion
Do you remember what criteria should be used in FDM to guarantee stabilization ? We help you to remember :
2alphadt/dx^2 << 1
In [31]:
dt=0.25 # time step
mstep=1600 # number of iteration
nx=200 # number of nodes in x direction
lx=100. # length of the plate
dx=lx/nx # space between two adjacent node
x=numpy.linspace(0,lx,nx+1) # matrix which store each node position
Tf=numpy.zeros(nx+1) # Temperature finite difference
To=numpy.zeros(nx+1) # Temperature storaage for previous time
Tf[0]=1.0 # Boundary condition left wall
Tf[nx]=0.0 # Boundary condition right wall
In [32]:
for n in range(mstep) :
To[:]=Tf[:]
for i in range(1,nx):
adv=dt*u*(To[i]-To[i-1])/dx
Tf[i]=To[i]+ ( (dt*alpha/(dx**2.)) *(To[i+1]-2*To[i]+To[i-1]) -adv)
In [33]:
plt.figure(figsize=(10,5), dpi=100)
plt.xlabel('x', fontsize=10) #x label
plt.ylabel('T', fontsize=10) #y label
plt.plot(x,Tf, 'r-',label=' Lattice Boltzmann Method');
plt.plot(xl,T, 'b*',label=' Finite Difference Method ');
plt.legend();