0

Why is my LBM code ignoring the source term?

I am trying to solve the one-phase stefan problem of freezing using a D2Q5 lattice with an iterative enthalpy method.

Any advice would be appreciated. The code is below.

#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[i]
            index = i
 
    s[n]=index/nx    
#==============================

My code does not produce the expected results when compared to the analytical solution. It apparently solves the heat conduction equation only and ignores the liquid fraction (Fl) dependent source term in the collision step.

Gwen
  • 11

0 Answers0