为什么我的 LBM 代码会忽略源术语？
我正在尝试使用具有 迭代焓法 的 D2Q5 晶格来解决冻结的单相 stefan 问题。
如有任何建议，我们将不胜感激。代码如下。
#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
#==============================
与解析解相比，我的代码没有产生预期的结果。它显然只求解了热传导方程，而忽略了碰撞步骤中与液体分数 (Fl) 相关的源项。