Lattice Boltzmann simulation
Jishnu Rajendran
Lattice Boltzmann methods are a family of computational methods for simulating the evolution of fluid flows in systems. In this notebook, we implement a 2D square domain with rigid barrier using Lattice Boltzmann method. The simulation is based on the LBM model and the boundary conditions are discussed in the following sections. Instead of solving the Navier–Stokes equations directly, a fluid density on a lattice is simulated with streaming and collision (relaxation) processes. As a versatile model, the dynamics of the fluid can be simulated fairly straight. The LBM model can be made to mimic common fluid behaviour like vapour/liquid coexistence, and so fluid systems such as liquid droplets can be simulated.
1. Lattice Boltzmann simulation of a 2D square domain with rigid barrier
""" Lattice-Boltzmann method for fluid simulation Simple rectangular barrier @author: Jishnu """ import numpy, time, matplotlib.pyplot, matplotlib.animation height = 80 # dimensions of lattice width = 200 viscosity = 0.02 # viscosity omega = 1 / (3*viscosity + 0.5) # parameter for relaxation u0 = 0.1 # initial and in-flow speed f_n = 4.0/9.0 # lattice-Boltzmann weight factors o_n = 1.0/9.0 o_36 = 1.0/36.0 performanceData = True # True if performance data is needed
Here we initialize the dimensions of the domain and the initial conditions of fluid flow. We will choose steady flow and initialise the density and velocity arrays.
# Initialize arrays --steady rightward flow: n0 = f_n * (numpy.ones((height,width)) - 1.5*u0**2) # particle densities along 9 directions nN = o_n * (numpy.ones((height,width)) - 1.5*u0**2) nS = o_n * (numpy.ones((height,width)) - 1.5*u0**2) nE = o_n * (numpy.ones((height,width)) + 3*u0 + 4.5*u0**2 - 1.5*u0**2) nW = o_n * (numpy.ones((height,width)) - 3*u0 + 4.5*u0**2 - 1.5*u0**2) nNE = o_36 * (numpy.ones((height,width)) + 3*u0 + 4.5*u0**2 - 1.5*u0**2) nSE = o_36 * (numpy.ones((height,width)) + 3*u0 + 4.5*u0**2 - 1.5*u0**2) nNW = o_36 * (numpy.ones((height,width)) - 3*u0 + 4.5*u0**2 - 1.5*u0**2) nSW = o_36 * (numpy.ones((height,width)) - 3*u0 + 4.5*u0**2 - 1.5*u0**2) rho = n0 + nN + nS + nE + nW + nNE + nSE + nNW + nSW # macroscopic density ux = (nE + nNE + nSE - nW - nNW - nSW) / rho # macroscopic x velocity uy = (nN + nNE + nNW - nS - nSE - nSW) / rho # macroscopic y velocity
1.1. Effects of a Barrier in a steady flow
We choose a simple rectangle barrier in the domain.
barrier = numpy.zeros((height,width), bool) # True wherever there's a barrier barrier[(height//2)-8:(height//2)+8, (height//2)-4:(height//2)+4] = True # simple linear barrier barrierN = numpy.roll(barrier, 1, axis=0) # sites just north of barriers barrierS = numpy.roll(barrier, -1, axis=0) # sites just south of barriers barrierE = numpy.roll(barrier, 1, axis=1) barrierW = numpy.roll(barrier, -1, axis=1) barrierNE = numpy.roll(barrierN, 1, axis=1) barrierNW = numpy.roll(barrierN, -1, axis=1) barrierSE = numpy.roll(barrierS, 1, axis=1) barrierSW = numpy.roll(barrierS, -1, axis=1)
2. Move all particles by one step along their directions of motion (pbc):
def stream(): global nN, nS, nE, nW, nNE, nNW, nSE, nSW nN = numpy.roll(nN, 1, axis=0) # axis 0 is north-south; + direction is north nNE = numpy.roll(nNE, 1, axis=0) nNW = numpy.roll(nNW, 1, axis=0) nS = numpy.roll(nS, -1, axis=0) nSE = numpy.roll(nSE, -1, axis=0) nSW = numpy.roll(nSW, -1, axis=0) nE = numpy.roll(nE, 1, axis=1) # axis 1 is east-west; + direction is east nNE = numpy.roll(nNE, 1, axis=1) nSE = numpy.roll(nSE, 1, axis=1) nW = numpy.roll(nW, -1, axis=1) nNW = numpy.roll(nNW, -1, axis=1) nSW = numpy.roll(nSW, -1, axis=1) # Using boolean arrays to handle barrier collisions (bounce-back): nN[barrierN] = nS[barrier] nS[barrierS] = nN[barrier] nE[barrierE] = nW[barrier] nW[barrierW] = nE[barrier] nNE[barrierNE] = nSW[barrier] nNW[barrierNW] = nSE[barrier] nSE[barrierSE] = nNW[barrier] nSW[barrierSW] = nNE[barrier]
def collide(): """ Calculates the collision step of the Lattice Boltzmann Method (LBM) algorithm. Updates the macroscopic variables `rho`, `ux`, and `uy` based on the population distributions `n0`, `nN`, `nS`, `nE`, `nW`, `nNE`, `nNW`, `nSE`, and `nSW`. Parameters: None Returns: None """ global rho, ux, uy, n0, nN, nS, nE, nW, nNE, nNW, nSE, nSW rho = n0 + nN + nS + nE + nW + nNE + nSE + nNW + nSW ux = (nE + nNE + nSE - nW - nNW - nSW) / rho uy = (nN + nNE + nNW - nS - nSE - nSW) / rho ux2 = ux * ux uy2 = uy * uy u2 = ux2 + uy2 omu215 = 1 - 1.5*u2 uxuy = ux * uy n0 = (1-omega)*n0 + omega * f_n * rho * omu215 nN = (1-omega)*nN + omega * o_n * rho * (omu215 + 3*uy + 4.5*uy2) nS = (1-omega)*nS + omega * o_n * rho * (omu215 - 3*uy + 4.5*uy2) nE = (1-omega)*nE + omega * o_n * rho * (omu215 + 3*ux + 4.5*ux2) nW = (1-omega)*nW + omega * o_n * rho * (omu215 - 3*ux + 4.5*ux2) nNE = (1-omega)*nNE + omega * o_36 * rho * (omu215 + 3*(ux+uy) + 4.5*(u2+2*uxuy)) nNW = (1-omega)*nNW + omega * o_36 * rho * (omu215 + 3*(-ux+uy) + 4.5*(u2-2*uxuy)) nSE = (1-omega)*nSE + omega * o_36 * rho * (omu215 + 3*(ux-uy) + 4.5*(u2-2*uxuy)) nSW = (1-omega)*nSW + omega * o_36 * rho * (omu215 + 3*(-ux-uy) + 4.5*(u2+2*uxuy)) # Force steady rightward flow at ends # no need to set 0, N, and S component nE[:,0] = o_n * (1 + 3*u0 + 4.5*u0**2 - 1.5*u0**2) nW[:,0] = o_n * (1 - 3*u0 + 4.5*u0**2 - 1.5*u0**2) nNE[:,0] = o_36 * (1 + 3*u0 + 4.5*u0**2 - 1.5*u0**2) nSE[:,0] = o_36 * (1 + 3*u0 + 4.5*u0**2 - 1.5*u0**2) nNW[:,0] = o_36 * (1 - 3*u0 + 4.5*u0**2 - 1.5*u0**2) nSW[:,0] = o_36 * (1 - 3*u0 + 4.5*u0**2 - 1.5*u0**2)
# Compute curl of the velocity field: def curl(ux, uy): """ Calculates the curl of a vector field. Parameters: ux (numpy.ndarray): The x-component of the vector field. uy (numpy.ndarray): The y-component of the vector field. Returns: numpy.ndarray: The curl of the vector field. """ return numpy.roll(uy,-1,axis=1) - numpy.roll(uy,1,axis=1) - numpy.roll(ux,-1,axis=0) + numpy.roll(ux,1,axis=0)
3. Visualization of the simulation
# for animation. theFig = matplotlib.pyplot.figure(figsize=(8,3)) fluidImage = matplotlib.pyplot.imshow(curl(ux, uy), origin='lower', norm=matplotlib.pyplot.Normalize(-.1,.1), cmap=matplotlib.pyplot.get_cmap('jet'), interpolation='none') bImageArray = numpy.zeros((height, width, 4), numpy.uint8) # an RGBA image bImageArray[barrier,3] = 255 # set alpha=255 barrier sites only barrierImage = matplotlib.pyplot.imshow(bImageArray, origin='lower', interpolation='none') # Function called for each successive animation frame: startTime = time.perf_counter() #frameList = open('frameList.txt','w') # file containing list of images def nextFrame(arg): # (arg is the frame number, which we don't need) global startTime if performanceData and (arg%100 == 0) and (arg > 0): endTime = time.perf_counter() print( "%1.1f" % (100/(endTime-startTime)), 'frames per second' ) startTime = endTime #frameName = "frame%04d.png" % arg #matplotlib.pyplot.savefig(frameName) #frameList.write(frameName + '\n') for step in range(15): # adjust number of steps for smooth animation stream() collide() fluidImage.set_array(curl(ux, uy)) return (fluidImage, barrierImage) # return the figure elements to redraw animate = matplotlib.animation.FuncAnimation(theFig, nextFrame, interval=0.5, blit=True) matplotlib.pyplot.show()