February 4, 2024
Jake Bobowski
A toy model to simulate the settling of particles suspended in a solution due to gravity. The simulation works as follows:
The block of code below takes approximately 2 minutes to run.
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, clear_output
import time
import warnings
with warnings.catch_warnings():
warnings.simplefilter("ignore") # Suppress warning messages that appear when evaluating instances of ln(0) in np.log(count).
N = int(1e4) # Number of particles
num_steps = int(1e2) # Number of interations/steps
x = np.random.rand(N) # Step the initial coordinates of the N particles.
y = np.random.rand(N)
x_initial = x.copy() # Copy in the initial coordinates.
y_initial = y.copy()
fig, (ax, ax1) = plt.subplots(1, 2, figsize=(15, 7)) # Set up the plots to be dynamically updated while the loop runs.
ax.set_aspect('equal')
h = np.linspace(0.05, 1, 19) # We will determine the particle density as a function of the height h from the bottom of the y-range.
bin_edges = np.linspace(0.05, 1, 20) # We will count the particles between for 0.05 < y < 0.1, 0.1 < y < 0.15, ..., 0.95 < y < 1.
for step in range(num_steps):
dx = np.random.binomial(1, 0.5, N) # For each step, the probabilty of a particle scattering left or right is the same.
dy = np.random.binomial(1, 0.25, N) # There is a greater probability that a particle moves down because of gravity.
x_step = abs(np.random.normal(0, 1, N) * 3e-2) # Randomize the sizes of the steps in the x and y directions for each particle and each iteration of the loop.
y_step = abs(np.random.normal(0, 1, N) * 3e-2)
x += np.where(dx == 0, -x_step, x_step) # Update the x coordinates.
x = np.where(x < 0, x + 1, x) # If x < 0, add 1. A particle that moves past the left boundary will appear near the right boundary.
x = np.where(x > 1, x - 1, x) # If x > 1, subtract 1. A particle that moves past the right boundary will appear near the left boundary.
y += np.where(dy == 0, -y_step, y_step) # Update the y coordinates.
y = np.clip(y, 0, 1) # Particles that go beyond the top boundary (y > 1) have their y-coordinate reset to y = 1.
# Particles that go beyond the lower boundary (y < 0) have their y-coordinate reset to y = 0.
ax.clear()
ax.scatter(x_initial, y_initial, s=0.5, color='cyan', label='Initial Positions') # Plot the initial positions of the particles.
ax.scatter(x, y, s=0.5, color='orange', label='Dynamic Positions') # Dynamically plot the new positions.
ax.legend(loc='upper left')
ax.set_xlabel("x position")
ax.set_ylabel("y position")
display(fig)
clear_output(wait=True)
time.sleep(1e-4)
count, _ = np.histogram(y, bins=bin_edges) # Determine the particle cound has a function of the height h.
ax1.clear()
ax1.scatter(h, np.log(count)) # Plot the ln of the particle count as a function of h. The very bottom of the particle chamber (0 < y < 0.05) is excluded.
ax1.set_ylim(0, None) # Set the minimum y-range to be zero
ax1.annotate(f'Step: {step}', xy=(0.02, 0.05), xycoords='axes fraction', fontsize=12, color='black')
ax1.set_xlabel("height")
ax1.set_ylabel("ln(counts)")