A Gentle Introduction to Lattice Gas Automaton for Simulation of Fluid Flow with Python. [The FHP model for Navier-Stokes Equations]
Modeling fluid flow has been an area of interest for many researchers and scientists in the field of engineering and applied physics. In the past decades, many techniques have been developed to model the kinematics or the dynamics of the fluid. The common approach is to use the famous Navier-Stokes equations, which fully describe the behavior of fluids. These equations would then be solved using numerical techniques like Finite Differences.
Although these methods provide a very good approximation to the system, they have a high computational cost and they are generally hard to work with. In 1986, Uriel Frisch et al invented a very intelligent model that was able to simulate the fluid flow. Their model, which they called Lattice Gas Automaton (LGA) was able to reproduce the Navier-Stokes equations without any numerical approximations. While most of the existing models for fluid flow at that time started with the macroscopic variables and the partial differential equations to describe the continuous system, LGA followed the opposite direction. The model simulates the macroscopic behavior of the system by starting with simple rules of the microscopic system which get coarse-grained to reproduce the general behavior of the fluid flow. The model can be seen as a small universe that simplifies continuous complex equations into discrete simple rules. LGA is considered to be one of the best models for fluid flow until this moment.
In this post, I’m going to gently explain and describe how LGA works, and then I will implement a full simulation using Python.
LGA is a simple 2D lattice that is divided into LxL cells. Each of these cells has a hexagonal shape. The choice of the hexagonal shape is important since hexagonal cells are usually used to simulate a Navier-Stokes equation of fluid flow. They are preferable because of the isotropy of the momentum advection tensor. Here is how the model looks:
Each one of these cells (hexagons) is used to represent one particle of the fluid. This means that one particle can have up to six surrounding neighbors (one neighbor at each side of the hexagon). The model basically simulates the motion of these particles and their collisions with each other and with any obstacles in the system. The particles move in discrete-time steps with constant velocities where each particle has its own vector (direction of motion). Particles don’t change their direction (keep moving in their chosen vector) until collision happens. Each cell can have one particle at a time which means that if two particles go towards the same cell at the same time, collision is going to happen. This restriction is called the exclusion principle; it ensures that six Boolean variables at each lattice site are always enough to represent the microdynamics.
Since particles move with a constant velocity, then at each time step, they move one cell along their velocity vector (move to one of their neighboring cells). The main physical rules that the model uses to simulate the fluid flow are the conservation of mass and momentum. These physical rules are applied by adding restrictions to the particles’ interaction. The restrictions are a set of configurations of collisions that can change the direction of the particles. In other words, for this model to work, there are only a few configurations that can lead to a non-trivial collision (i.e. a collision in which the directions of motion have changed). These configurations are the following:
- When two particles collide head to head, both of them are deflected by 60 degrees so that the output of the collision is still a zero momentum configuration with two particles.
- When three particles collide with an angle of 120 degrees between each other, they bounce back to their original direction. This ensures that the momentum is conserved.
To model these behaviors mathematically, each cell (hexagon) in the lattice is going to be represented using the following equation:
The equation refers to the number of particles entering cell “r” at time step “t” and going in the direction “i”. Since the model makes the restriction that each cell can have up to one particle only, then “n” can only be 0 or 1. 0 means that no particle is coming to the cell towards that direction and 1 means there is a particle coming towards that direction. Since each cell is a hexagon and it has six sides, then there are six directions for particles to move towards (i=0, 1, …6). This means that each cell can be fully described using six equations.
The evolution equation for each cell is then described using the following formula:
Δr is just the distance that the particle traveled. Since this is a discrete mode and the particles travel one cell per time step, then r+Δr is just the neighboring cell to the cell r. Cᵢ means that the particle is going in the direction i (to the neighboring cell of cell r that is in the direction i where i is one of the six directions of the hexagon). Δt is just the time difference. So this equation is just saying that if there is a particle entering cell r at time t and going towards the direction i, then at time t+Δt that particle will enter cell r+Δr and will still keep going in the direction i. This is assuming that no collision happens at any of the two cells. If collision happens, then that particle can change its direction and move to another cell.
Using that equation, we can now list the cases when collision can happen at a specific site “r”:
- Head-on collision will happen if nᵢ =1 and nᵢ₊₃ = 1 where i is any direction. This is because each cell is a hexagon with six sides and so any two facing sides (head to head sides) are separated by three other sides (so i and i+3). nᵢ =1 means that a particle is entering the cell and going towards direction i and nᵢ₊₃ = 1 means that a particle is entering the cell and going towards direction i+3. This means that a head-on collision should happen. For this to be true, no other particles should enter from the other directions (e.g. all other nᵢ =0). If another particle enters the same site from another direction, then this is going to be a three particle collision instead of a head-on collision. To mathematically describe that, we use the following equation:
If the quantity Dᵢ =1, then head-on collision is going to happen across i and i+3 and no particles will end up going in either of these two directions (because particles will get deflected by 60 degrees as described above). The equation simply multiplied nᵢ by nᵢ₊₃ if the result is 1 then two particles are going in opposite directions (potential head-on collision). The remainder of the equation just checks if a third particle is coming to the cell from other directions (if that’s the case, then Dᵢ should be equal to 0 indicating that head-on collision isn’t happening).
Similarly, for three-body collisions, nᵢ, nᵢ₊₂, nᵢ₊₄ all have to be equal to 1 and all other nᵢ have to equal 0. This is represented using the following equation:
Similar to Dᵢ, Tᵢ is either 0 or 1 and indicates whether three-body collision happens.
Now, let’s talk about the result of collisions. For the three-body collisions, the result is easy as particles just bounce back to where they came from. However, for two-body collisions, particles change directions and there are two possible directions to take (see fig 2 above). For example, a particle can end up going in the direction i because of a head-on collision at nᵢ₊₁ or nᵢ₊₄ or a head-on collision between nᵢ₋₁ and nᵢ₊₂. Each has a 50% probability. Then, we can introduce a random variable q where q=1 means particles deflect to the left and q=0 means particles deflect to the right. Then the number of particles that end up in the direction i due to head-on collision is:
Not that this quantity can only be 0 or 1 because q is a random variable that can only be 0 or 1.
Now, we can finally put all these together to update our evolution equation to be the following:
where Ωᵢ is called the collision term. Its used to tell whether a collision (either head-on or three-body will happen). It’s given by the following formula:
As you might have noticed, Ωᵢ can be 0, 1, or -1. The first term -Dᵢ just checks if head-on collision happens. This term is negtaive because if it happens and no other collision happens then Ωᵢ will be equal to -1 and when added to nᵢ(r, t) the result should be zero which indicates that no particle will remain the in direction i due to the head-on collision. The second two terms just check if a particle will end up in the direction i due to a head-on collision in the other directions. The last two terms check if a three-body collision will happen which should remove or add the particle to the direction i.
Now, we have enough equations to fully model the system. These equations are enough to describe the flow of any fluids on a the 2D grid lattice. Let’s go ahead and start building that model in Python.
Let’s start by building the 2D lattice. I’m going to add walls at the boundary so this system is going to be as fluid in a box:
def initialize():
vectorized_object = np.vectorize(lambda obj: node(False))
config = np.zeros((rows, columns)).astype(object)
config[:] = vectorized_object(config[:])
#adding the walls:
vectorized_object_wall = np.vectorize(lambda obj: node(True))
config[0, :] = vectorized_object_wall(config[0, :])
config[rows-1, :] = vectorized_object_wall(config[rows-1, :])
config[:, 0] = vectorized_object_wall(config[:, 0])
config[:, columns-1] = vectorized_object_wall(config[:, columns-1])
return config# a function to observe and visualize the grid:def observe(config):
vectorized_x = np.vectorize(lambda obj: obj.occupied) view_grid = vectorized_x(config)
plt.figure(figsize = (10, 8))
pylab.cla()
pylab.imshow(view_grid, vmin = 0, vmax = 1, cmap = pylab.cm.binary)
Now, I will create a class for the cells:
class node:
def __init__(self, wall):
self.ni_s = [0, 0, 0, 0, 0, 0]
self.wall = wall
if not wall:
self.occupied = 0
else:
self.occupied = 1The above class will represent each cell (hexagon) in the lattice. Note that Ni_s has six values to represent the six directions. If a node is set to be a wall, then it’s value will be 1 (while 0 means empty cell).
Now, I write a function to get the neighbors of a cell given its coordinates in the lattice:
def get_neighbors(coordinates, config):
row, column = coordinates
neighbors = []
neighbors.append([row, (column+1)%columns]) #right
neighbors.append([(row-1)%rows, (column+1)%columns]) # top right
neighbors.append([(row-1)%rows, (column-1)%columns]) #top left
neighbors.append([row, (column-1)%columns]) #left
neighbors.append([(row+1)%rows, (column-1)%columns]) #bottom left
neighbors.append([(row+1)%rows, (column+1)%columns]) # bottom right
return neighborsNote that I’m using a usual square lattice but to make it hexagonal, I’m only using six surrounding cells for each center cell.
I will write a function to fill the grid randomly with N particles:
#getting all possible coordinates and choosing N radnomly:i_coords, j_coords = np.meshgrid(range(rows), range(columns), indexing='ij')
coordinate_grid = np.array([i_coords, j_coords])
all_coordinates = []
for i in range(1, rows-1):
for j in range(1, columns-1):
all_coordinates.append(list(coordinate_grid[:, i, j]))
def fill_random(N, config, d = False):
samples = random.sample(all_coordinates, N)
for row,column in samples:
config[row, column].occupied = 1
our_neighbors = get_neighbors([row, column], config)
while True:
if not d:
direction = random.choice([i for i in range(6)])
else:
direction = d-1
n_row, n_col = our_neighbors[direction]
if not config[n_row,n_col].wall:
if config[n_row,n_col].occupied == 0:
config[n_row,n_col].ni_s[direction]=1
break
return config
Let’s try the code:
rows = 300
columns = 300
my_grid = initialize()
my_grid = fill_random(1000, my_grid)
observe(my_grid)
That looks good. Now, let’s code down the functions:
def D_i(n, i):
"""
The head-on collision function
"""
return n[i%6]*n[(i+3)%6]*(1-n[(i+1)%6])*(1-n[(i+2)%6])*(1-n[(i+4)%6])def T_i(n, i):
"""
The three-body collision function
"""
return n[i]*n[(i+2)%6]*n[(i+4)%6]*(1-n[(i+1)%6])*(1-n[(i+3)%6])*(1-n[(i+5)%6])
def collision_factor(coords, config, i):
"""
The collision factor Ω
"""
q = random.choice([0, 1])
n = config[coords[0], coords[1]].ni_s
coll = -D_i(n, i)+q*D_i(n, (i-1)%6) + (1-q)*D_i(n, (i+1)%6) - T_i(n, i)+T_i(n, (i+3)%6)
return coll
Now, let’s get to the big part. I write a function to use all of the above and simulate the system:
def update(config):
next_config = initialize()
for row in range(rows):
for column in range(columns): # for each cell in the grid
if not config[row, column].wall: # if it's not a wall #get all of its neighbors
our_neighbors = get_neighbors([row, column], config)
bounce_direc = [] # make this list to store if it's going to bounce off a wall
for i in range(6):# for each neighbor # start by letting the particle that is coming along i (from i+3) get in
if config[row, column].ni_s[i] >0: # then update so that the coming particle can enter the cell
next_config[row, column].occupied = 1
th_row, th_col = our_neighbors[(i+3)%6]
# the cell that sent that particle
if not next_config[th_row, th_col].wall:
next_config[th_row, th_col].occupied = 0 # update it to 0 since the particle moved already
next_config[row, column].ni_s[i] = 0#then for each neghibor, check if something will pop inside them from this cell and change their n_i
other_row, other_col = our_neighbors[i]
#calculate its value using the FHP equation
next_config[other_row,other_col].ni_s[i] = config[row, column].ni_s[i]+\
collision_factor([row, column], config, i)
# check if any particle is hitting the wall:
if next_config[other_row,other_col].wall and next_config[other_row,other_col].ni_s[i]>0:
next_config[other_row,other_col].ni_s[i] = 0
bounce_direc.append((i+3)%6) #store their bounce off directions
# if any particles were hiting the wall, send them back:
if bounce_direc:
next_config[row, column].ni_s[bounce_direc[0]] = 1return next_config
In the function above, I’m basically doing the following steps:
- For each cell in the lattice, if that cell is not a wall:
- Get all the neighbors of the cell
- For each neighbor, see if there is a particle coming in from that neighbor.
- If yes, let that particle in and remove it from the neighbor
- Then, for each of these neighbors, check if a particle will go to them from the current cell using the evolution equation (update their nᵢs)
- If yes, update their nᵢ
- Check if any particle is about to hit the wall, and make it bounce back by updating its neighboring cell’s nᵢ
Finally, let’s run the simulation:
def turn_state(config):
"""
A function to visualize the states
"""
view_grid = np.copy(config)
for i in range(len(config)):
for j in range(len(config)):
if type(config[i,j]) is not int:
view_grid[i, j] = config[i,j].occupied
view_grid = view_grid.astype(int)
return view_gridmy_grid = initialize()
my_grid = fill_random(1000, my_grid)
observe(my_grid)states = []
states.append(turn_state(my_grid))
next_grid = update(my_grid)
states.append(turn_state(next_grid))states.append(turn_state(next_grid))
for i in tqdm(range(200)):
next_grid = update(next_grid)
states.append(turn_state(next_grid))build_animation(states, "video.mp4")
Here is the output:
Looking awsome!!
Usually, to simulate a fluid flow, part of this wall can be removed and replaced with a periodic boundary. We can include that using the following small edits to the code:
def initialize_2():
vectorized_object = np.vectorize(lambda obj: node(False))
config = np.zeros((rows, columns)).astype(object)
config[:] = vectorized_object(config[:])
return configdef update_periodic(config):
#next_config = config
next_config = initialize_2()
for row in range(rows):
for column in range(columns):
our_neighbors = get_neighbors([row, column], config)
for i in range(6):
other_row, other_col = our_neighbors[i]
#calculate its value using the FHP equation
next_config[other_row,other_col].ni_s[i] = config[row, column].ni_s[i]+\
collision_factor([row, column], config, i)if config[row, column].ni_s[i] >0:# then update so that the coming particle can enter the cell
next_config[row, column].occupied = 1
th_row, th_col = our_neighbors[i]
next_config[th_row, th_col].occupied = 0
next_config[row, column].ni_s[i] = 0return next_config
The four walls have been replaced by periodic boundaries.
The above was created using only 1000 particles. Now, you can increase the number of particles and initialize the particles to a specific direction according to the fluid you are trying to simulate.
Hope this article was useful. Here is the full code.
