February 6, 2024
Jake Bobowski
A toy model to simulate the diffusion of partlicles that are injected into the centre of container filled with a fluid. Like using an eyedropper to put a drop of food colouring in the middle of a jar filled with water. The simulation works as follows:
The block of code below takes several minutes to run.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # Importing 3D plotting tools
from IPython.display import display, clear_output
import time
import warnings
with warnings.catch_warnings():
warnings.simplefilter("ignore") # Suppress warning messages.
N = int(1e4) # Number of particles
num_steps = int(5e2) # Number of interations/steps
x = np.ones(N)*0.5 # Step the initial coordinates of the N particles to be exactly in the centre of the 1 x 1 cell (x, y) = (0.5, 0.5).
y = np.ones(N)*0.5
z = np.ones(N)*0.5
x_initial = x.copy() # Copy in the initial coordinates.
y_initial = y.copy()
z_initial = z.copy()
fig = plt.figure(figsize=(13, 5)) # Set up the figure for 3D plotting
ax = fig.add_subplot(121, projection='3d') # 3D subplot for particle positions
ax1 = fig.add_subplot(122) # 2D subplot for particle count
ax.set_box_aspect([1,1,1]) # Equal aspect ratio for 3D plot
h = np.linspace(0, 1, 20) # 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, 1, 21) # We will count the particles between 0 < y < 0.05, 0.05 < y < 0.1, ..., 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.5, N) # For each step, the probabilty of a particle scattering up or down is the same.
dz = np.random.binomial(1, 0.5, N) # For each step, the probabilty of a particle scattering up or down is the same.
x_step = abs(np.random.normal(0, 1, N) * 2e-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) * 2e-2)
z_step = abs(np.random.normal(0, 1, N) * 2e-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.where(y < 0, y + 1, y) # If y < 0, add 1. A particle that moves past the bottom boundary will appear near the top boundary.
y = np.where(y > 1, y - 1, y) # If x > 1, subtract 1. A particle that moves past the top boundary will appear near the bottom boundary.
z += np.where(dz == 0, -z_step, z_step) # Update the y coordinates.
z = np.where(z < 0, z + 1, z) # If z < 0, add 1. A particle that moves past the bottom boundary will appear near the top boundary.
z = np.where(z > 1, z - 1, z) # If z > 1, subtract 1. A particle that moves past the top boundary will appear near the bottom boundary.
ax.clear()
ax.scatter(x_initial, y_initial, z_initial, color='cyan', label='Initial Positions') # Plot the initial positions of the particles.
ax.scatter(x, y, z, s = 1, color='orange', label='Dynamic Positions') # Dynamically plot the new positions.
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.set_zlim(0, 1)
ax.legend(loc='upper left')
ax.set_xlabel("x position")
ax.set_ylabel("y position")
ax.set_zlabel("z position")
display(fig)
clear_output(wait=True)
time.sleep(1e-4)
count, _ = np.histogram(z, bins = bin_edges) # Determine the particle count has a function of the height h.
ax1.clear()
ax1.scatter(h, count) # Plot the particle count as a function of h.
ax1.set_xlim(0, 1)
ax1.set_ylim(0, None) # Set the minimum y-range to be zero
ax1.annotate(f'Step: {step}', xy=(0.02, 0.95), xycoords='axes fraction', fontsize=12, color='black')
ax1.set_xlabel("height")
ax1.set_ylabel("counts")