Thermal lattice Boltzmann model ignoring source term -- python code help please

In summary: I am not sure what is causing this error, but it seems like there may be a problem with the initial conditions or boundary conditions being used. It would be helpful to double-check those and possibly adjust them to avoid this error.In summary, the LBM model for phase change solves a 1D 1 phase problem using the D2Q5 lattice. Relevant equations can be found in the code. The model uses thermal diffusion coefficient, lattice dimensions, number of steps, and wall temperatures as input parameters. The main loop consists of collision and streaming processes, with boundary conditions set at the left and right walls. The model also calculates nodal enthalpy and liquid fraction, with convergence being checked before each iteration.
  • #1
Gwen
1
1
LBM model for phase change- relevant equations found here. Also here.

Code:
#Thermal LBM
#solves 1D 1 phase phase-change
#D2Q5 Lattice

nx=100                                   # the number of nodes in x direction lattice direction
ny=5                                    # the number of nodes in y direction lattice direction
alpha=.5/3                     # heat diffusion coefficient                                   # the dimension of the problem
mstep=1000                               # the number of time step
tau=3.*alpha+0.5

Tleft=0.0                                 # left wall temperature
Tright=1.0                                # right wall temperature
k=5 # k=0,1,2,3,4,5,6,8,9

x=numpy.linspace(0,1,nx+1) #start at zero, end at 1, fill with nx+1 even spaced intervals
y=numpy.linspace(0,1,ny+1)
t=np.zeros(mstep)
s=np.zeros(mstep)
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

Hl=1
Hs=0.5
H=numpy.ones((ny+1,nx+1) )                   # Enthalpy matrix
Fl=numpy.ones((ny+1,nx+1) )                  # Liquid fraction matrix (Fl=1 for liquid, Fl=0 for solid)

##================ Initial boundary condition
w[0]=1./3. #0.0
w[1:5]=1./6. #1./4.

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

T[0:ny+1,0:nx+1]=1.0   #temperature in the whole region (including bottom wall)
T[0:ny+1,0]=0        #temperature on the left wall
T[0:ny+1,nx]=1.0       #temperature one node in from the right wall
T[ny,1:nx]=1.0         #temp one node in from the top wall (and one node in from left and right sides)

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
##=====================
for n in range(mstep) :
    t[n]=n  #track the time
    time=t[n]

    epsilon=1e-8
    error=1
    Fl_old=Fl
    while error>epsilon:
        Fl_old_iter=Fl
        T_old_iter=T

# 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.-1/tau)*f[l,j,i]+(1/tau)*feq-w[l]*(Fl[j,i]-Fl_old[j,i])
 
 #streaming process
# ==========================
        for i in range(nx):
            for j in range(ny,0,-1):  #backwards from top to bottom
                f[2,j,i]=f[2,j-1,i]

        for i in range(nx,0,-1):   #backwards from right to left
            for j in range(ny,0,-1):  #backwards from top to bottom
                f[1,j,i]=f[1,j,i-1]

        for i in range(nx,0,-1):   #backwards from right to left
            for j in range(0,ny):     #forward from bottom to second-to-top lattice node
                f[4,j,i]=f[4,j+1,i]

        for i in range(0,nx):      #forward from left to second-to-right lattice node
            for j in range(0,ny):     #forward from bottom to second-to-top lattice node
                f[3,j,i]=f[3,j,i+1]

# Boundary conditions
#  =============================
        for j in range(0,ny+1) :               #left Boundary. Dirichlet boundary condition: constant temperature.
            f[1,j,0]=( Tleft*(w[1]+w[3]) )-f[3,j,0]

        for j in range(0,ny+1):                #right Boundary. adiabatic
            f[3,j,nx]=f[1,j,nx]

        for i in range(0,nx+1):                # bottom and top Boundary
            f[4,ny,i]=f[2,ny,i]                  #adiabatic
 #================================ #calculate temperature
        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           #Dirichlet BC      
        T[0:ny+1,nx]=T[0:ny+1,nx]   #adiabatic BC        
        T[ny,1:nx]=T[ny-1,1:nx]     #adiabatic BC        
        T[0,1:nx]=T[1,1:nx]         #adiabatic BC       
#==============================   #calculate nodal enthalpy and liquid fraction
        for i in range(nx+1):
            for j in range(ny+1):
                H[j,i]=0.5*T[j,i]+0.5*Fl[j,i]
        for i in range(nx+1):
            for j in range(ny+1):
                if H[j,i]<=Hs:
                    Fl[j,i]=0
                elif H[j,i]>Hs and H[j,i] < Hl:
                    Fl[j,i]=(H[j,i]-Hs)/(Hl-Hs)
                else:
                    Fl[j,i]=1
#==============================   #convergence? If no, go back
        for i in range(nx+1):
            for j in range(ny+1):
                error_Fl=abs(np.max(np.max((Fl[j,i]-Fl_old_iter[j,i])/Fl[j,i])))
                error_T=abs(np.max(np.max((T[j,i]-T_old_iter[j,i])/T[j,i])))
                error=np.max([error_Fl, error_T])      
#find position of phase change boundary (where Fl<=0.5)  
    Fl_col=Fl[3,:]<=0.5
    max = Fl_col[0]
    index = 0
    for i in range(1,len(Fl_col)):
        if Fl_col[I] >= 0.5:
            max = Fl_col
            index = i[/I][/I]
 
        s[n]=index/nx   
#==============================[/I][/I]
 
Last edited by a moderator:
Technology news on Phys.org
  • #2
I get a division by zero error at line 119
 

1. How does the Thermal Lattice Boltzmann model work?

The Thermal Lattice Boltzmann model is a numerical method used to simulate fluid flow and heat transfer at a mesoscopic level. It works by dividing the fluid domain into a lattice grid and assigning probabilities for fluid particles to move from one grid point to another. These probabilities are based on the fluid's macroscopic properties, such as density and velocity, and are updated at each time step to simulate the flow and transport of heat.

2. What is the purpose of ignoring the source term in the Thermal Lattice Boltzmann model?

The source term in the Thermal Lattice Boltzmann model accounts for any external forces or heat sources acting on the fluid. By ignoring this term, the model becomes simpler and more computationally efficient, making it easier to implement and analyze. However, this also means that the model may not accurately capture all aspects of the fluid flow and heat transfer, so it should be used with caution and for specific applications.

3. Can the Thermal Lattice Boltzmann model be implemented in Python?

Yes, the Thermal Lattice Boltzmann model can be implemented in Python. There are various open-source codes and libraries available that provide the necessary functions and algorithms for simulating fluid flow and heat transfer using the model. Additionally, Python's flexibility and readability make it a popular choice for researchers and engineers working with this model.

4. Are there any limitations to using the Thermal Lattice Boltzmann model?

Like any numerical model, the Thermal Lattice Boltzmann model has its limitations. It is most suitable for simulating low-speed, incompressible flows and may not accurately capture turbulent or compressible flow phenomena. Additionally, the model's accuracy can be affected by the chosen lattice geometry and grid resolution, so careful validation and sensitivity analysis are crucial.

5. How can I validate the results of the Thermal Lattice Boltzmann model?

The results of the Thermal Lattice Boltzmann model can be validated by comparing them to analytical solutions or experimental data. This can help identify any discrepancies or limitations of the model and provide insights for improving its accuracy. Additionally, conducting sensitivity analyses by varying model parameters can also help validate the results and assess the model's robustness.

Similar threads

  • Programming and Computer Science
Replies
1
Views
943
  • Programming and Computer Science
Replies
5
Views
1K
  • Programming and Computer Science
Replies
6
Views
1K
  • Programming and Computer Science
Replies
4
Views
4K
  • Programming and Computer Science
Replies
6
Views
3K
  • Programming and Computer Science
Replies
34
Views
2K
  • MATLAB, Maple, Mathematica, LaTeX
Replies
1
Views
3K
  • Programming and Computer Science
Replies
8
Views
5K
  • Atomic and Condensed Matter
Replies
3
Views
861
  • Engineering and Comp Sci Homework Help
Replies
2
Views
3K
Back
Top