Lattice Boltzmann Simulation

reference:
https://github.com/pmocz/latticeboltzmann-python/blob/main/latticeboltzmann.py
https://physics.weber.edu/schroeder/fluids/

The lattice Boltzmann method is a simple way to simulate the fluid dynamics. You can find a detail information here.

The code below is mostly based on one developed here. I changed the left boundary condition, the initial flow pattern, and the shape of an object.

For the following example, I put a bar as an object and Karman vortex sheet can be seen.

In [2]:
%config InlineBackend.figure_format='retina'
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import matplotlib.animation as animation
mpl.rcParams['animation.html'   ]='jshtml'
In [51]:
def main(Nx,Ny,Nt,u0,visco):
    # collision timescale
    tau       = (3*visco + 0.5)
    
    # Lattice speeds / weights
    NL   = 9
    idxs = np.arange(NL)
    cxs  = np.array([0, 0, 1, 1, 1, 0,-1,-1,-1])
    cys  = np.array([0, 1, 1, 0,-1,-1,-1, 0, 1])
    weights = np.array([4/9,1/9,1/36,1/9,1/36,1/9,1/36,1/9,1/36]) # sums to 1
    
    # Initial Conditions
    F = np.ones((Ny,Nx,NL)) 
    for i,cx,cy,w in zip(idxs,cxs,cys,weights):
        F[:,:,i]=w*(F[:,:,i]+cx*3*u0+abs(cx)*4.5*u0**2-1.5*u0**2)
           
    # object boundary
    X, Y = np.meshgrid(range(Nx), range(Ny))
    objct = (X - Nx/4)**2 + (Y - Ny/2)**2 < (Ny/16)**2 #cylinder
    objct = (Nx/4<X)&(X<Nx/4+5)&(Ny/3<Y)&(Y<Ny-Ny/3)   #bar
    
    # history
    isav=200
    hst=np.zeros((Ny,Nx,Nt//isav))
    
    # Simulation Main Loop
    for it in range(Nt):
        # Drift
        for i, cx, cy in zip(idxs, cxs, cys):
            F[:,:,i] = np.roll(F[:,:,i], cx, axis=1)
            F[:,:,i] = np.roll(F[:,:,i], cy, axis=0)
        
        # Left boundary    
        for i,cx,cy,w in zip(idxs,cxs,cys,weights):
            if(i!=0 and i!=1 and i!=5): 
                F[:,0,i]=w*(np.ones(Ny)+cx*3*u0+abs(cx)*4.5*u0**2-1.5*u0**2)
       
        # Set reflective boundaries
        bndryF = F[objct,:]
        bndryF = bndryF[:,[0,5,6,7,8,1,2,3,4]]
        
        # Calculate fluid variables
        rho = np.sum(F,2)
        ux  = np.sum(F*cxs,2) / rho
        uy  = np.sum(F*cys,2) / rho
        
        # Apply Collision
        Feq = np.zeros(F.shape)
        for i, cx, cy, w in zip(idxs, cxs, cys, weights):
            Feq[:,:,i] = rho * w * ( 1 + 3*(cx*ux+cy*uy)  + 9*(cx*ux+cy*uy)**2/2 - 3*(ux**2+uy**2)/2 )
        F += -(1.0/tau) * (F - Feq)
        
        # Apply boundary 
        F[objct,:] = bndryF
        
        # store the vorticity into hst
        if (it % isav == 0):
            ux[objct] = 0
            uy[objct] = 0
            vorticity = (np.roll(ux, -1, axis=0) - np.roll(ux, 1, axis=0))\
                      - (np.roll(uy, -1, axis=1) - np.roll(uy, 1, axis=1))
            vorticity[objct] = np.nan
            hst[:,:,it//isav]=vorticity
    
    return locals()
In [52]:
Nx=200
Ny=80
Nt=60000
u0=0.08
visco=0.01
dat=main(Nx,Ny,Nt,u0,visco)
locals().update(dat)
In [61]:
def update_anim(it):
    
    fig.clf() #clear the figure
    
    ax=fig.subplots() #add subplots

    fig.tight_layout() #reduce spacing around the figure
    cmap = plt.cm.bwr
    cmap.set_bad('black')
    
    im1=ax.imshow(hst[1:Ny-1,1:Nx-1,it] ,aspect='equal', origin='lower', cmap='bwr',vmin=-0.04, vmax=0.04)
    
    cb1=fig.colorbar(im1, ax=ax, shrink=0.8,label=r'$\nabla\times{\bf u}$')
    ax.set_xlabel('X')
    ax.set_ylabel('Y')

nt=Nt//isav

fig=plt.figure(figsize=(10,4))
anim=animation.FuncAnimation(fig,update_anim,frames=nt)
#plt.subplots_adjust(wspace=0.1,hspace=0.01, bottom=0,top=1)

plt.close()
anim
#anim.save('sample.gif', writer='imagemagick') #save the animation as a gif file

/Users/mnakanot/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:9: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("bwr"))
  if __name__ == '__main__':
Out[61]: