Fick’s First Law: Diffusion to Steady State

ficks_first_law_steady_state.py

import numpy as np, matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt, io, os, math
import PIL.Image as Image
import imageio.v2 as imageio

# Parameters
L = 1.0
D = 1.0
Nx = 201
x = np.linspace(0, L, Nx)
dx = x[1]-x[0]

C_left = 1.0
C_right = 0.0

# Explicit scheme stability: dt <= dx^2/(2D)
dt = 0.45 * dx**2 / D
t_end = 0.6  # nondimensional-ish
n_steps = int(t_end/dt) + 1

# Initial condition
C = np.ones(Nx)*0.5
C[0] = C_left
C[-1] = C_right

alpha = D*dt/dx**2

# Frame settings
fps = 24
duration_s = 7.0
n_frames = int(fps*duration_s)

# Choose frame times (nonlinear spacing for nicer approach)
# Use quadratic in normalized time to show early dynamics more
frame_times = (np.linspace(0, 1, n_frames)**1.2) * t_end
frame_indices = np.clip((frame_times/dt).astype(int), 0, n_steps-1)

# Precompute target indices set
targets = set(frame_indices.tolist())

# Storage
frames = []
flux_t = []
time_t = []

# Steady-state flux expected
J_ss = D*(C_left - C_right)/L  # since C(x) linear, dC/dx = (C_right-C_left)/L, J=-D*dC/dx
J_ss