# 2D Fluid Simulation using FHP LGCA (Lattice Gas Cellular Automata) # Simulates fluid flow in a circular channel. # Particles go out from right side and enter back from left. # Reference: # Lattice Gas Cellular Automata and Lattice Boltzmann Models by Wolf-Gladrow # FB - 20140818 import math import random from PIL import Image imgx = 512; imgy = 512 # image size image = Image.new("RGB", (imgx, imgy)) pixels = image.load() # simulation parameters: tilesX = 32 tilesY = 32 n = 8 # coarse graining tile size is n by n timeSteps = 300 nodesX = tilesX * n nodesY = tilesY * n nodes = [[[0 for x in range(nodesX)] for y in range(nodesY)] for z in range(6)] obstacle = [[0 for x in range(nodesX)] for y in range(nodesY)] # insert a square obstacle in the middle for y in range(nodesY / 4): for x in range(nodesX / 4): obstacle[y + nodesY / 2 - nodesY / 8][x + nodesX / 2 - nodesX / 8] = 1 # fill-up with fluid flowing towards right for y in range(1, nodesY - 1): # do not include top/bottom walls for x in range(nodesX): if obstacle[y][x] != 1: nodes[0][y][x] = 1 for t in range(timeSteps): # run the simulation # HANDLE COLLISIONS # collisions at non-boundary nodes for y in range(1, nodesY - 1): # do not include top/bottom walls for x in range(nodesX): if obstacle[y][x] != 1: cell = [nodes[z][y][x] for z in range(6)] numParticles = sum(cell) # only 2 or 3 symmetric particle collisions implemented here if numParticles == 3: if cell[0] == cell[2] and cell[2] == cell[4]: # invert the cell contents for z in range(6): nodes[z][y][x] = 1 - cell[z] elif numParticles == 2: # find the cell of one of the particles p = cell.index(1) # its diametric opposite must occupied as well if p > 2: pass elif cell[p + 3] == 0: pass else: # randomly rotate the particle pair clockwise or # counterclockwise if random.randint(0, 1) == 0: # counterclockwise nodes[0][y][x] = cell[5] nodes[1][y][x] = cell[0] nodes[2][y][x] = cell[1] nodes[3][y][x] = cell[2] nodes[4][y][x] = cell[3] nodes[5][y][x] = cell[4] else: # clockwise nodes[0][y][x] = cell[1] nodes[1][y][x] = cell[2] nodes[2][y][x] = cell[3] nodes[3][y][x] = cell[4] nodes[4][y][x] = cell[5] nodes[5][y][x] = cell[0] # collisions along top/bottom walls (no-slip) for x in range(nodesX): cell = [nodes[z][0][x] for z in range(6)] nodes[0][0][x] = cell[3] nodes[1][0][x] = cell[4] nodes[2][0][x] = cell[5] nodes[3][0][x] = cell[0] nodes[4][0][x] = cell[1] nodes[5][0][x] = cell[2] cell = [nodes[z][nodesY - 1][x] for z in range(6)] nodes[0][nodesY - 1][x] = cell[3] nodes[1][nodesY - 1][x] = cell[4] nodes[2][nodesY - 1][x] = cell[5] nodes[3][nodesY - 1][x] = cell[0] nodes[4][nodesY - 1][x] = cell[1] nodes[5][nodesY - 1][x] = cell[2] # collisions at obstacle points (no-slip) for y in range(nodesY): for x in range(nodesX): if obstacle[y][x] == 1: cell = [nodes[z][y][x] for z in range(6)] nodes[0][y][x] = cell[3] nodes[1][y][x] = cell[4] nodes[2][y][x] = cell[5] nodes[3][y][x] = cell[0] nodes[4][y][x] = cell[1] nodes[5][y][x] = cell[2] # HANDLE MOVEMENTS nodesNew = [[[0 for x in range(nodesX)] for y in range(nodesY)] for z in range(6)] for y in range(nodesY): for x in range(nodesX): cell = [nodes[z][y][x] for z in range(6)] # propagation in the 0-direction neighbor_y = y if x == nodesX - 1: neighbor_x = 0 else: neighbor_x = x + 1 nodesNew[0][neighbor_y][neighbor_x] = cell[0] # propagation in the 1-direction if y != nodesY - 1: neighbor_y = y + 1 if y % 2 == 1: if x == nodesX - 1: neighbor_x = 1 else: neighbor_x = x + 1 else: neighbor_x = x nodesNew[1][neighbor_y][neighbor_x] = cell[1] # propagation in the 2-direction if y != nodesY - 1: neighbor_y = y + 1 if y % 2 == 0: if x == 0: neighbor_x = nodesX - 1 else: neighbor_x = x - 1 else: neighbor_x = x nodesNew[2][neighbor_y][neighbor_x] = cell[2] # propagation in the 3-direction neighbor_y = y if x == 0: neighbor_x = nodesX - 1 else: neighbor_x = x - 1 nodesNew[3][neighbor_y][neighbor_x] = cell[3] # propagation in the 4-direction if y != 0: neighbor_y = y - 1 if y % 2 == 0: if x == 0: neighbor_x = nodesX - 1 else: neighbor_x = x - 1 else: neighbor_x = x nodesNew[4][neighbor_y][neighbor_x] = cell[4] # propagation in the 5-direction if y != 0: neighbor_y = y - 1 if y % 2 == 1: if x == nodesX - 1: neighbor_x = 0 else: neighbor_x = x + 1 else: neighbor_x = x nodesNew[5][neighbor_y][neighbor_x] = cell[5] nodes = nodesNew print '%' + str(100 * t / timeSteps) # show progress # Create an image from the final state # Calculate average velocity vectors for tiles aveVelocityVectorMag = [[0.0 for x in range(tilesX)] for y in range(tilesY)] aveVelocityVectorAng = [[0.0 for x in range(tilesX)] for y in range(tilesY)] pi2 = math.pi * 2.0 dx = [math.cos(i * pi2 / 6.0) for i in range(6)] dy = [math.sin(i * pi2 / 6.0) for i in range(6)] for ty in range(tilesY): for tx in range(tilesX): vx = 0.0 vy = 0.0 for cy in range(n): for cx in range(n): for z in range(6): if nodes[z][ty * n + cy][tx * n + cx] == 1 \ and obstacle[ty * n + cy][tx * n + cx] == 0: vx += dx[z] vy += dy[z] aveVelocityVectorMag[ty][tx] = math.hypot(vx, vy) / n ** 2.0 aveVelocityVectorAng[ty][tx] = (math.atan2(vy, vx) + pi2) % pi2 for ky in range(imgy): iy = nodesY * ky / imgy jy = tilesY * ky / imgy for kx in range(imgx): ix = nodesX * kx / imgx jx = tilesX * kx / imgx if obstacle[iy][ix] == 1: # paint the obstacle(s) red = 0 grn = 0 blu = 255 else: # use vector magnitude and angle for coloring aveVelVecMag = aveVelocityVectorMag[jy][jx] aveVelVecAng = aveVelocityVectorAng[jy][jx] red = int(aveVelVecMag * 255) grn = int(aveVelVecAng / pi2 * 255) blu = 0 pixels[kx, ky] = (red, grn, blu) image.save("FHP_LGCA_2DFluidSim.png", "PNG")