Discrete Elements

The Nature of Mathematical Modeling (draft)

$ \def\%#1{\mathbf{#1}} \def\mat#1{\mathbf{#1}} \def\*#1{\vec{#1}} \def\ve#1{\vec{#1}} \def\ds#1#2{\cfrac{d #1}{d #2}} \def\dd#1#2{\cfrac{d^2 #1}{d #2^2}} \def\ps#1#2{\cfrac{\partial #1}{\partial #2}} \def\pp#1#2{\cfrac{\p^2 #1}{\p #2^2}} \def\p{\partial} \def\ba{\begin{eqnarray*}} \def\ea{\end{eqnarray*}} \def\eps{\epsilon} \def\del{\nabla} \def\disp{\displaystyle} $

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
scale = .5
width = 10*scale
height = 7*scale
marker = 2*height
linew = marker/2.5
font = 3*height
fig,ax = plt.subplots(figsize=(width,height))
ax.axis([0,width,0,height])
ax.set_axis_off()
def line(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    plt.plot([x0,x1],[y0,y1],'k-',linewidth=0.8*linew)
def point(x,y):
    x = scale*x
    y = scale*y
    plt.plot(x,y,'ko',markersize=marker)
def text(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        fontsize=font,color='black')
def itext(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        style='italic',
        fontsize=font,color='black')
def gtext(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        style='italic',
        fontsize=font,color=(.3,.3,.3))
def arrow(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    ax.annotate('',xy=(x1,y1),xytext=(x0,y0),
        arrowprops=dict(color='black',width=marker/7,headwidth=marker,headlength=marker))
gtext(2.5,2,'molecular\ndynamics')
text(5.5,4,'particle\nsystems')
gtext(8.5,6,'partial\ndifferential\nequations')
itext(2.5,0.7,'discrete')
itext(8.5,0.7,'continuous')
itext(0.5,2,'micro')
itext(0.5,6,'macro')
arrow(1.3,1.2,1.3,7)
arrow(1.3,1.2,10,1.2)
arrow(3.5,2.5,4.6,3.5)
arrow(7.4,5.4,6.3,4.4)
plt.show()
  • molecular dynamics
    • start from fundamental physics
    • femtosecond time steps
    • can take days per microseconds
    • Hollingsworth, Scott A., and Ron O. Dror. "Molecular dynamics simulation for all." Neuron 99, no. 6 (2018): 1129-1143.
  • can satisfy spirit but not letter
    • redefine molecule, macromolecular dynamics
  • mesh-free, computation happens where particles are
    • can handle complexity where PDE approximations can break down
    • is massively parallel
    • trillions of particles \cite{Tchipev:19}
    • Tchipev, Nikola, Steffen Seckler, Matthias Heinen, Jadran Vrabec, Fabio Gratl, Martin Horsch, Martin Bernreuther et al. "TweTriS: Twenty trillion-atom simulation." The International Journal of High Performance Computing Applications 33, no. 5 (2019): 838-854.

Particle Systems

SPH (Smoothed Particle Hydrodynamics)

function $f(\ve{r})$

expand integral as sum

$$ \int f(\ve{r})~dr = \sum_i V_i f_i $$

interpolate with weighting function

$$ f(\ve{r}) = \sum_i V_i f_i W(\ve{r}-\ve{r}_i) $$

preserve normalization

$$ \int W(\ve{r})~d\ve{r} = 1 $$

kernel function approximation, Function Fitting

relate volume element to a density and mass, introducing particles

$V_i = m_i/\rho_i$

$$ f(\ve{r}) = \sum_i \frac{m_i}{\rho_i} f_i W(\ve{r}-\ve{r}_i) $$

derivatives go into the sum

$$ \del f(\ve{r}) = \sum_i \frac{m_i}{\rho_i} f_i \del W(\ve{r}-\ve{r}_i) $$

Navier-Stokes

$$ \rho\left(\ps{\ve{v}}{t}+\ve{v}\cdot\del\ve{v}\right) = - \del p + \rho \ve{g} + \mu \del^2\ve{v} $$

quantities relative to particles in Lagrangian (moving frame) vs Eulerian (fixed frame)

Lagrangian derivatives

$$ \ba \frac{Df}{Dt} &=& \frac{f(\ve{x}+\ve{v}dt,t+dt)-f(\ve{x},t)}{dt} \\ &\approx& \frac{f(\ve{x},t)+\ps{f}{t}dt+\ve{v}\cdot \del f dt-f(\ve{x},t)}{dt} \\ &=& \ps{f}{t}+\ve{v}\cdot \del f \ea $$$$ \rho\frac{D\ve{v}}{dt} = - \del p + \rho \ve{g} + \mu \del^2\ve{v} $$

pressure

$$ \del p(\ve{r}_j) = \sum_i \frac{m_i}{\rho_i} p_i \del W(\ve{r_i}-\ve{r}_j) $$

density

$$ \ba \rho(\ve{r}_j) &=& \sum_i \frac{m_i}{\rho_i} \rho_i W(\ve{r}_i-\ve{r}_j) \\ &=& \sum_i m_i W(\ve{r}_i-\ve{r}_j) \ea $$

viscosity, depends on relative velocity

$$ \mu\del^2\ve{v}(\ve{r}_j) = \mu\sum_i \frac{m_i}{\rho_i} (\ve{v}_i-\ve{v}_j) \del^2 W(\ve{r}-ve{r}_i) $$

equation of state relates pressure to density

converted PDE to ODEs for particle motion

interpolate for fields

computation happens where the particles are, not over a grid

Gingold, Robert A., and Joseph J. Monaghan. "Smoothed particle hydrodynamics: theory and application to non-spherical stars." Monthly notices of the royal astronomical society 181, no. 3 (1977): 375-389.

Fraga Filho, Carlos Alberto Dutra, Carlos Alberto Dutra Fraga Filho, and Castro. Smoothed particle hydrodynamics. Springer International Publishing, 2019.

TLA (Three Letter Acronyms)

MPM (Material Point Method)

used for solid mechanics

hybrid particle-grid method

Lagrangian particles represent material properties with internal state variables

Eulerian grid used for computing spatial derivatives in governing equations

alternate Lagrangian and Eulerian steps

used for a range of governing equations

Sulsky, Deborah, Shi-Jian Zhou, and Howard L. Schreyer. "Application of a particle-in-cell method to solid mechanics." Computer physics communications 87, no. 1-2 (1995): 236-252.

Sołowski, Wojciech T., Martin Berzins, William M. Coombs, James E. Guilkey, Matthias Möller, Quoc Anh Tran, Tito Adibaskoro, Seyedmohammadjavad Seyedan, Roel Tielen, and Kenichi Soga. "Material point method: overview and challenges ahead." Advances in Applied Mechanics 54 (2021): 113-204.

DPD (Dissipative Particle Dynamics)

Groot, Robert D., and Patrick B. Warren. "Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation." The Journal of chemical physics 107, no. 11 (1997): 4423-4435.

PPM (Parallel Particle-Mesh)

Sbalzarini, Ivo F., Jens H. Walther, Michael Bergdorf, Simone Elke Hieber, Evangelos M. Kotsalis, and Petros Koumoutsakos. "PPM–A highly efficient parallel particle–mesh library for the simulation of continuum systems." Journal of Computational Physics 215, no. 2 (2006): 566-588.

MCA (Movable Cellular Automata)

Psakhie, S. G., Y. Horie, G. P. Ostermeyer, S. Yu Korostelev, A. Yu Smolin, E. V. Shilko, A. I. Dmitriev, S. Blatnik, M. Špegel, and S. Zavšek. "Movable cellular automata method for simulating materials with mesostructure." Theoretical and applied fracture mechanics 37, no. 1-3 (2001): 311-334.

FPM (Finite Particle Method)

Liu, M. B., W. P. Xie, and G. R. Liu. "Modeling incompressible flows using a finite particle method." Applied mathematical modelling 29, no. 12 (2005): 1252-1270.

VPM (Vortex Particle Method)

Winckelmans, Grégoire S., and Anthony Leonard. "Contributions to vortex particle methods for the computation of three-dimensional incompressible unsteady flows." Journal of computational physics 109, no. 2 (1993): 247-273.

DSMC (Direct Simulation Monte Carlo)

Oran, Elaine S., C. K. Oh, and B. Z. Cybyk. "Direct simulation Monte Carlo: recent advances and applications." Annual Review of Fluid Mechanics 30, no. 1 (1998): 403-441.

PD (Peridynamics)

start from integral form of solid mechanics governing equations

use particle expansion to replace integral with sum

particles interact through bonds within a horizon

fracture represented by breaking bonds

Silling, Stewart A., and Richard B. Lehoucq. "Peridynamic theory of solid mechanics." Advances in applied mechanics 44 (2010): 73-168.

Hattori, Gabriel, Mark Hobbs, and John Orr. "A Review on the Developments of Peridynamics for Reinforced Concrete Structures: G. Hattori et al." Archives of Computational Methods in Engineering 28, no. 7 (2021): 4655-4686.

DEM (Discrete Element Methods)

SPH starts from PDE and expands in particles

In Differential and Difference Equations found PDE from masses and springs

In Finite Elements found masses and springs from FEA

can start from masses and springs

Cundall, Peter A., and Otto DL Strack. "A discrete numerical model for granular assemblies." geotechnique 29, no. 1 (1979): 47-65.

Tavarez, Federico A., and Michael E. Plesha. "Discrete element method for modelling solid and particulate materials." International journal for numerical methods in engineering 70, no. 4 (2007): 379-404.

originally for granular media, contact interactions

become a generic term

generalize to coarse-grained particles

MIPS (Memoryless Isotropic Point Particles)

simplest particles

RISC vs CISC modeling

contrary to earlier results, can model anisotropy and hysteresis through particle geometry

Strand, Erik, Filippos Tourlomousis, and Neil Gershenfeld. "Mesoscale material modeling with memoryless isotropic point particles." Journal of Computational Science 75 (2024): 102198.

Forces

Newton's equation of motion $\ve{F}=m\ve{a}$

add up forces on particle

demlaw.jpg Figure: A MIPS force law (from [Strand:24])

demforces.jpg Figure: MIPS force laws (from [Strand:24])

24.01.MIPS.jpg Figure: COMSOL (top) vs MIPS (bottom) solutions to the Sandia fracture challenge (right) [Strand:24]

Sorting

$O(N^2)$ particle interactions, $O(N)$ sort

lexicographic sort

construct key

place in pre-defined bins using particle position as key

can use a cumulative sum sort to find the bin sizes

Geometry

Maxwell criterion DOF vs constraints

Maxwell, J. Clerk. "L. on the calculation of the equilibrium and stiffness of frames." The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 27, no. 182 (1864): 294-299.

rectangular lattice shear, cross-brace bookcase

sphere packing

2D triangular

3D Hexagonal Close-Packed (HCP), Face-Centered Cubic (FCC)

cannon ball stack

elastic constants not isotropic

can pour in particles

can sample with noise distribution

Integration

Euler

$$ \ba \ve{a}(t) &=& \ve{F}\left(\ve{x}(t)\right)/m \\ \ve{v}(t+h) &=& \ve{v}(t) +h \ve{a}(t) \\ \ve{x}(t+h) &=& \ve{x}(t) +h \ve{v}(t) \ea $$

first order accurate

semi-implicit

$$ \ba \ve{a}(t) &=& \ve{F}\left(\ve{x}(t)\right)/m \\ \ve{v}(t+h) &=& \ve{v}(t) +h \ve{a}(t) \\ \ve{x}(t+h) &=& \ve{x}(t) +h \ve{v}(t+h) \ea $$

same order

is symplectic, conserves phase space volume

position Verlet

Verlet, Loup. "Computer" experiments" on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules." Physical review 159, no. 1 (1967): 98.

Taylor series, expand position forwards

$$ \ve{x}(t+h) = \ve{x}(t)+h\ve{v}(t)+h^2\ve{a}(t)/2+\ldots $$

expand position backwards

$$ \ve{x}(t-h) = \ve{x}(t)-h\ve{v}(t)+h^2\ve{a}(t)/2+\ldots $$

add these

$$ \ve{x}(t+h) = 2\ve{x}(t)-\ve{x}(t-h)+h^2\ve{a}(t) $$

second order accurate

need to keep old position

can derive from semi-implicit

$$ \ba \ve{x}(t+h) &=& \ve{x}(t)+h\left(\ve{v}(t)+h\ve{a}(t)\right) \\ &\approx& \ve{x}(t)+h((\ve{x}(t)-\ve{x}(t-h))/h+h\ve{a}(t)) \\ &=& 2\ve{x}(t)-\ve{x}(t-h)+h^2\ve{a}(t) \ea $$

velocity Verlet

$$ \ba \ve{x}(t+h) &=& \ve{x}(t)+h\ve{v}(t)+\frac{h^2}{2}\ve{a}(t) \\ \ve{v}(t+h) &=& \ve{v}(t)+h\ve{a}(t)+\frac{h^2}{2}\ds{\ve{a}}{t} \\ \ds{\ve{a}}{t} &\approx& (\ve{a}(t+h)-\ve{a}(t))/h \\ \ve{v}(t+h) &=& \ve{v}(t)+\frac{h}{2}(\ve{a}(t+h)+\ve{a}(t)) \ea $$

need to keep old $\ve{a}$

can combine to save memory

$$ \ba \ve{v}(t+h/2) &=& \ve{v}(t)+h\ve{a}(t)/2 \\ \ve{x}(t+h) &=& \ve{x}(t)+h\ve{v}(t+h/2)\\ \ve{a}(t+h) &=& \ve{F}(\ve{x}(t+h))/m\\ \ve{v}(t+h) &=& \ve{v}(t+h/2)+h\ve{a}(t+h)/2 \ea $$

can combine first and last $v$ updates

dissipation proportional to $\ve{v}$

can be internal, inertial, numerical

Stability

literally a stiff ODE: widely varying time scales

numerical stability limited by the speed of sound

linear array of masses and springs speed equal to square root of the elastic force divided by the mass

convert to a time step by dividing lattice pitch by sound velocity

PBD (Position-Based Dynamics)

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
scale = 1
width = 3*scale
height = 2.7*scale
line = width/3
marker = 3*width
font = 6*width
fig,ax = plt.subplots(figsize=(width,height))
ax.axis([0,width,0,height])
ax.set_axis_off()
def wire(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    plt.plot([x0,x1],[y0,y1],'k-',linewidth=line)
def point(x,y):
    x = scale*x
    y = scale*y
    plt.plot(x,y,'ko',markersize=marker)
def text(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        fontsize=font,color='black')
def arrow(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    ax.annotate('',xy=(x1,y1),xytext=(x0,y0),
        arrowprops=dict(color='black',width=marker/5,headwidth=marker,headlength=marker))
def lightarrow(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    ax.annotate('',xy=(x1,y1),xytext=(x0,y0),
        arrowprops=dict(color=(0.7,0.7,0.7),width=marker/5,headwidth=marker,headlength=marker))
def rectangle(x,y,w,h):
    x = scale*x
    y = scale*y
    w = scale*w
    h = scale*h
    #rect = patches.Rectangle((x,y),w,h,facecolor='gray')
    rect = patches.Rectangle((x,y),w,h,facecolor=(0.7,0.7,0.7))
    ax.add_patch(rect)
def circle(x,y,r):
    x = scale*x
    y = scale*y
    plt.plot(x,y,'ko',markersize=r,markeredgewidth=width/2,
        markerfacecolor='white',markeredgecolor='black')
rectangle(0,0,3,1.25)
circle(0.5,2.5,20)
circle(1.5,0.5,20)
circle(1.5,1.42,20)
lightarrow(0.65,2.25,1.38,1.6)
arrow(0.65,2.25,1.41,0.71)
arrow(1.5,0.71,1.5,1.15)
plt.show()

Figure: Position-Based Dynamics surface constraint

Müller, Matthias, Bruno Heidelberger, Marcus Hennix, and John Ratcliff. "Position based dynamics." Journal of Visual Communication and Image Representation 18, no. 2 (2007): 109-118.

Macklin, Miles, Matthias Müller, and Nuttapong Chentanez. "XPBD: position-based simulation of compliant constrained dynamics." In Proceedings of the 9th International Conference on Motion in Games, pp. 49-54. 2016.

widely used for computer graphics, qualitative rather than quantitative

rigid body collision has infinite gradient

scattering can go into the interior, or be unstable

instead, integrate forces to project $\ve{x}$, correct by applying position constraints, revise the velocity

integration can be for particles or rigid bodies

constraints can be stretching, bending, colliding, ...

constraints can be equality $C(\ve{x}) = 0$ or inequality $C(\ve{x}) \ge 0$

solve with constraint projection

$$C(\ve{x}+\Delta\ve{x}) \approx C(\ve{x})+\del C(\ve{x})\cdot \Delta\ve{x} = 0$$

move $\ve{x}$ in direction of most rapid change in the constraint

$$\Delta\ve{x} = \lambda \del C(\ve{x})$$

substitute to find $\lambda$

$$\Delta\ve{x} = -\frac{C(\ve{x})}{|\del C(\ve{x})|^2}\del C(\ve{x})$$

apply iteratively to all particles

CA (Cellular Automata)

Particle systems use discrete representations, still require algorithms for solution

What about using computation as the representation?

Finkel, Leif H., and Gerald M. Edelman. "Stanislaw Ulam: the warmth and brilliance of an eclectic mind." letters in mathematical physics 10, no. 2 (1985): 243-245.

Von Neumann, John, and Arthur Walter Burks. "Theory of self-reproducing automata." (1966).

Shannon, Claude Elwood, and John McCarthy, eds. Automata studies. Vol. 34. Princeton University Press, 2016.

Games, Mathematical. "The fantastic combinations of John Conway’s new solitaire game “life” by Martin Gardner." Scientific American 223, no. 120-123 (1970): 11.

The idea was developed by Ulam {Finkel:85}, von Neumann {vonNeumann:66}, and colleagues {Shannon:56} in the 1950s. A classic example of a CA is Conway's Game of Life {Gardner:70}, in which occupied sites on a grid get born, survive, or die in succeeding generations based on the number of occupied neighboring sites. Any CA has the same elements: a set of connected sites, a discrete set of states that are allowed on the sites, and a rule for how they are updated. We will start by studying a slighly more complicated system that recovers the Navier-Stokes fluid equations, and then will consider the more general question of how cellular automata relate to computation.

LGA (Lattice Gas Automata)

Hydrodynamics was one of the earliest and best-developed application areas of cellular automata. A cellular automata model of a fluid (traditionally called a lattice gas) is specified by the geometry of a lattice, by the discrete states permitted at each site, and by an update rule for how the states change based on their neighbors. Both partial differential equations and molecular dynamics models use real numbers (for the values of the fields, or for the particle positions and velocities). A lattice gas discretizes everything so that just a few bits describe the state of each site on a lattice, and the dynamics reduce to a simple look-up table based on the values of the neighboring sites. Each site can be considered to be a parcel of fluid. The rules for the sites implement a ``cartoon'' version of the underlying microscopic dynamics, but should be viewed as operating on a longer length scale than individual particles. We will see that the conservation laws that the rules satisfy determine the form of the equivalent partial differential equations for a large lattice, and that the details of the rules set the parameter values.

Frisch, Uriel, Brosl Hasslacher, and Yves Pomeau. "Lattice-gas automata for the Navier-Stokes equation." Physical review letters 56, no. 14 (1986): 1505.

A historically important example of a lattice gas is the FHP rule (named after its inventors, Frisch, Hasslacher, and Pomeau {Frisch:86}. This operates in 2D on a triangular lattice, and the state of each site is specified by six bits (Figure ). Each bit represents a particle on one of the six links around the site, given by the unit vectors $\hat\alpha$. On each link a particle can either be present or absent, and all particles have the same unit velocity: in the absence of collisions, they travel one lattice step ahead in one time step. The simple update rule proceeds in two stages, chosen to conserve particle number and momentum (Figure ). First, collisions are handled. At the beginning of the step a particle on a link is considered to be approaching the site, and after the collision step a particle on a link is taken to be leaving the site. If a site has two particles approaching head-on, they scatter. A random choice is made between the two possible outgoing directions that conserve momentum (always choosing one of them would break the symmetry of the lattice). Because of the large number of sites in a lattice gas it is usually adequate to approximate the random decision by simply switching between the two outgoing choices on alternate site updates. If three particles approach symmetrically, they scatter. In all other configurations not shown the particles pass through the collision unchanged. After the collision step there is a transport step, in which each particle moves by one unit in the direction that it is pointing and arrives at the site at the far end of the link. While these rules might appear to be somewhat arbitrary, it will turn out that the details will not matter for the form of the governing equations, just the symmetries and conservation laws.

carule.jpg Figure: Update rules for an FHP lattice gas.

A simpler rule related to FHP is HPP (named after Hardy, de Pazzis, and Pomeau \cite{Hardy:76}), which operates on a square lattice. Each site is specified by four bits, and direction-changing collisions are allowed only when two particles meet head-on (unlike FHP, here there is only one possible choice for the exit directions after scattering). We will see that HPP and FHP, although apparently quite similar, behave very differently.

Let's label time by an integer $T$, the lattice sites by a vector $\ve{X}$, and the lattice directions by a unit vector $\hat\alpha$. If we start an ensemble of equivalent lattices off with the same update rule but different random initial conditions, we can define $f_\alpha(\vec X,T)$ to be the fraction of sites $\ve{X}$ at time $T$ with a particle on link $\hat\alpha$. In the limit of a large ensemble, this fraction becomes the probability to find a particle on that link. Defining this probability will let us make a connection between the lattice gas and partial differential equations.

At each time step, in the absence of collisions, the fraction of particles at site $\vec X$ at time $T$ pointing in direction $\hat\alpha$ will move one step in that direction: $$ f_\alpha(\vec X+\hat\alpha,T+1) = f_\alpha(\vec X,T) ~~~~. $$ Let's introduce new rescaled variables $\vec x = \delta_x \vec X$ and $t = \delta_t T$ in terms of the (small) space step $\delta_x$ and time step $\delta_t$. Substituting in these variables, collision-free transport becomes $$ f_\alpha(\vec x+\delta_x\hat\alpha,t+\delta_t)-f_\alpha(\vec x,t) = 0 ~~~~. $$ If the probability $f_\alpha$ varies slowly compared to $\delta_x$ and $\delta_t$, we can expand equation (\ref{ca-transp}) in $\delta_x$ and $\delta_t$: $$ \ps{f_\alpha(\ve{x},t)}{t}~\delta_t + \hat\alpha\cdot\del f_\alpha(\ve{x},t)~\delta_x + {\cal O}(\delta^2) = 0 ~~~~. $$ Choosing to scale the variables so that $\delta_x = \delta_t$, to first order this becomes $$ \ps{f_\alpha(\vec x,t)}{t} + \hat\alpha\cdot\del f_\alpha(\vec x,t) = 0 ~~~~. $$ This equation says that the time rate of change of the fraction of particles at a point is equal to the difference in the rate at which they arrive and leave the point by straight transport (remember equation (\ref{flux})). If collisions are allowed, the time rate of change of $f_\alpha$ will depend on both the spatial gradient and on a collision term $\Omega_\alpha$ scattering particles in or out from other directions $$ \ps{f_\alpha(\vec x,t)}{t} + \hat\alpha\cdot\del f_\alpha(\vec x,t) = \Omega_\alpha (\ve{x},t) ~~~~. $$ $f_\alpha$ is the distribution function to find a particle. The collision term $\Omega_\alpha$ will in general depend on the distribution function for pairs of particles as well as the one-particle distribution function, and these in turn will depend on the three-particle distribution functions, and so forth. This is called the BBGKY hierarchy of equations (Bogolyubov, Born, Green, Kirkwood, Yvon \cite{Boer:61}). The Boltzmann equation approximates this by assuming that $\Omega_\alpha$ depends only on the single-particle distribution functions $f_\alpha$.

We derived equation (\ref{CA-boltz}) by making use of the fact that particles travel one lattice site in each time step, and then assuming that $f_\alpha$ varies slowly. Now let's add the conservation laws that have been built into the update rules. The total density of particles $\rho$ at a site $\ve{x}$ is just the sum over the probability to find one in each direction $$ \sum_\alpha f_\alpha(\ve{x}) = \rho(\ve{x}) ~~~~, $$ and the momentum density is the sum of the probabilities times their (unit) velocities $$ \sum_\alpha \hat\alpha f_\alpha = \rho\vec v ~~~~. $$ Since our scattering rules conserve particle number (the particles just get reoriented), the number of particles scattering into and out of a site must balance $$ \sum_\alpha \Omega_\alpha(\vec x) = 0 ~~~~. $$ And since the rules conserve momentum, the net momentum change from scattering must vanish $$ \sum_\alpha \hat\alpha \Omega_\alpha = 0 ~~~~. $$ Therefore, summing equation (\ref{CA-boltz}) over directions, $$ \sum_\alpha \left[ \ps{f_\alpha(\vec x,t)}{t} + \hat\alpha\cdot\del f_\alpha(\vec x,t) \right] = \sum_\alpha \Omega_\alpha $$ $$ \ps{}{t}\sum_\alpha f_\alpha + \sum_\alpha \hat\alpha\cdot\del f_\alpha = 0 $$ $$ \ps{\rho}{t} + \del\cdot (\rho\vec v) = 0 ~~~~. $$ This is the familiar equation for the continuity of a fluid, and has arisen here because we've chosen scattering rules that conserve mass. A second equation comes from momentum conservation, multiplying equation (\ref{CA-boltz}) by $\hat\alpha$ and summing over directions $$ \ps{}{t}\sum_\alpha \hat\alpha f_\alpha + \sum_\alpha \hat\alpha (\hat\alpha\cdot\del f_\alpha) = 0 ~~~~. $$ The $i$th component of this vector equation is $$ \ps{}{t}\sum_\alpha \hat\alpha_i f_\alpha + \sum_\alpha \sum_j \hat\alpha_i \hat\alpha_j \ps{f_\alpha}{x_j} = 0 ~~~~. $$ Defining the momentum flux density tensor by $$ \Pi_{ij} \equiv \sum_\alpha \hat\alpha_i \hat\alpha_j f_\alpha ~~~~, $$ this becomes $$ \ps{\rho v_i}{t} + \sum_j \ps{\Pi_{ij}}{x_j} = 0 ~~~~. $$

We now have two equations, (\ref{ca-mass}) and (\ref{ca-tensor}), in three unknowns, $\rho$, $\ve{v}$, and $\Pi$. To eliminate $\Pi$ we can find the continuum form of the momentum flux density tensor by using a Chapman--Enskog expansion \cite{Huang:87}, a standard technique for finding approximate solutions to the Boltzmann equation. We will assume that $f_\alpha$ depends only on $\vec v$ and $\rho$ and their spatial derivatives (and not on time explicitly), and so will do an expansion in all possible scalars that can be formed from them. The lowest-order terms of the deviation from the equilibrium uniform configuration are $$ \ba f_\alpha &=& \frac{\rho}{6} \left( 1 + 2\hat\alpha\cdot\vec v + A \left[ (\hat\alpha\cdot\vec v)^2 - \frac{1}{2}|\vec v|^2 \right] \right.\n\\ & & + \left. B \left[ (\hat\alpha\cdot\del)(\hat\alpha\cdot\vec v) - \frac{1}{2}\del\cdot\vec v\right] + \cdots \right) ~~~~. \ea $$ The terms have been grouped this way to guarantee that the solution satisfies the density and momentum equations (\ref{ca-density}) and (\ref{ca-momentum}) (this can be verified by writing out the components of each term). In this derivation the only features of the FHP rule that we've used are the conservation laws for mass and momentum, and so all rules with these features will have the same form of the momentum flux density tensor (to this order), differing only in the value of the coefficients $A$ and $B$.

The Navier-Stokes governing equation for a $d$-dimensional fluid with bulk as well shear viscosity (ignoring gravitational forces) is $$ \ps{\rho\vec v}{t} + \rho(\vec v\cdot\del)\vec v = - \del p + \eta \del^2\vec v + (\zeta + \frac{\eta}{d}) \del (\del\cdot\vec v) ~~~~, $$ where $p$ is the pressure, $\eta$ is the shear viscosity, $\zeta$ is the bulk viscosity \cite{Batchelor:67}. Using the Chapman-Enskog expansion to evaluate the momentum flux density tensor in equation (\ref{ca-tensor}) and comparing it with the Navier--Stokes equation shows that they agree if $\zeta=0$, $\eta=\rho\nu = -\rho B/8$ ($\nu$ is the kinematic viscosity), $\mu=A/4$, and $p=\rho/2$. Further, in the Boltzmann approximation it is possible (with a rather involved calculation) to find the values of $A$ and $B$ for a given CA rule \cite{Wolfram:86}.

The simple conservation laws built into our lattice gas have led to the full Navier--Stokes equation; the particular rule determines the effective viscosity of the fluid. While the details of this calculation are complicated, there are some simple and important conclusions. The viscosity for the square lattice (HPP model) turns out to depend on direction and hence is not appropriate for most fluids, but the viscosity of the triangular lattice (FHP model) is isotropic. This profound implication of the lattice symmetry was not appreciated in the early days of lattice gas models, and helped point the way towards the realization that a simple lattice gas model could in fact be very general. In 3D the situation is more difficult because there is not a 3D lattice that gives an isotropic viscosity. However, it can still be achieved by using a quasiperidoic tiling that is not translationally periodic \cite{Boghosian:99}, or a cut through a higher-dimensional lattice such as the 4D Face-Centered Hyper-Cubic (FCHC) rule \cite{Frisch:87}. It is not possible to reduce the viscosity in a simple model like FHP (an attribute that is needed for modeling a problem such as air flow), but this can be done by adding more than one particle type or by increasing the size of the neighborhood used for the rule \cite{Dubrulle:91}.

Doolen, G. D. (Ed.). (1990). Lattice Gas Methods for Partial Differential Equations. Santa Fe Institute Studies in the Sciences of Complexity, Proceedings Volume IV. Reading, MA: Addison-Wesley.

This collection includes many of the important articles which work out the connection between lattice gases and hydrodynamics.

Rothman, D. H., & Zaleski, S. (1997). Lattice-Gas Cellular Automata: Simple Models of Complex Hydrodynamics (Collection Alea-Saclay: Monographs and Texts in Statistical Physics). Cambridge University Press.

Reviews the basic theory and extends it to porous media and fluids with multiple components.

LBM (Lattice Boltzmann Method)

LGA issues

statistical, must average

difficult to tune parameters

Lattice Boltzmann Method (LBM)

a particle in LGA represented by bit string, in LBM replaced by a real number for the single-particle distribution function

combine equations \ref{ca-transp} and \ref{CA-boltz}, take as a numerical method for $f_\alpha$

$$ f_\alpha(\vec x+\delta_x\vec\alpha,t+\delta_t) = f_\alpha(\vec x,t) + \Omega_\alpha (\ve{x},t) $$

generalized $\vec\alpha$ from a single lattice step to vector for multiple velocities

equivalent to a finite-difference form of the Boltzmann equation

again density

$$ \sum_\alpha f_\alpha(\ve{x}) = \rho(\ve{x}) $$

momentum density

$$ \sum_\alpha \vec\alpha f_\alpha = \rho\vec v $$

conservation mass

$$ \sum_\alpha \Omega_\alpha(\vec x) = 0 $$

conservation momentum

$$ \sum_\alpha \vec\alpha~\Omega_\alpha = 0 $$

BGK approximation to collision term

$$ \Omega_\alpha = \frac{f_\alpha-f_\alpha^{eq}}{\tau} $$

for $f_\alpha^{eq}$ can take Maxwell-Boltzmann distribution, power series expansion

LBM can be more accurate and faster than SPH, but can also be less robust and versatile

Chen, Hudong, Shiyi Chen, and William H. Matthaeus. "Recovery of the Navier-Stokes equations using a lattice-gas Boltzmann method." Physical review A 45, no. 8 (1992): R5339.

Timm, Krüger, Halim Kusumaatmaja, Alexandr Kuzmin, O. Shardt, G. Silva, and E. Viggen. "The lattice Boltzmann method: principles and practice." Cham, Switzerland: Springer International Publishing AG (2016).

Computing Automata

Simulating a cellular automata is an particularly simple type of computation. Rather than the many kinds of memory, instructions, and processors in a conventional computer, it requires just storage for the bits on the lattice, and an implementation of the local update rule. This is an extreme form of a SIMD (Single Instruction Multiple Data) parallel computer. The update can be performed by using the bits at each site as an index into a look-up table, so the architecture reduces to memory cycling through a look-up table (where the sequence of the memory retrieval determines the lattice geometry). This means that relatively modest hardware can exceed the performance of general-purpose computers for CA problems \cite{Toffoli:91}.

We've seen that a cellular automata computer can simulate fluids. A wide range of other physical systems also can be modeled by cellular automata rules; an important example is the rendering of 3D graphics \cite{Toffoli:97}. Conversely, instead of using a computer to simulate a CA, a CA can be used to simulate a computer. One way to do this is by implementing Boolean logic in CA rules; this was used to first prove their computational universality \cite{Banks:71}. This approach is attractive for hardware scaling because it builds in physical constraints; any technology that can perform the local cell updates can execute the same global programs \cite{Gershenfeld:10}.

Alternatively, since CAs can model physical systems, and physical systems can compute, CAs can model physical systems that compute \cite{Margolus:84}. Consider the billiard-ball CA in Figure~\ref{cagate} (the underlying lattice is not shown). This is similar to a lattice gas: billiard balls move on a lattice with unit velocity, and scatter off of each other and from walls. Two balls colliding generates the {\tt AND} function, and if one of the streams of balls is continuous it generates the {\tt NOT} function of the other input. These two elements are sufficent to build up all of logic \cite{Hill:93}. Memory can be implemented by delays, and wiring by various walls to guide the balls. The balls can be represented by four bits per site (one for each direction), with one extra bit per site needed to represent the walls.

cagate.jpg Figure: Billiard ball logic.

This kind of computing has many interesting features \cite{Fredkin:82}. No information is ever destroyed, which means that it is reversible (it can be run backwards to produce inputs from outputs) \cite{Bennett:88}, and which in turn means that it can be performed (in theory) with arbitrarily little dissipation \cite{Landauer:61}. Reversibility is also essential for designing quantum cellular automata, since quantum evolution is reversible. A quantum CA is much like a classical CA, but it permits the sites to be in a superposition of their possible states \cite{Lloyd:93}. This is a promising architecture for building quantum computers based on short-range interactions \cite{Chuang:00}.

Life in Life

For some, cellular automata are much more than just an amusing alternative to traditional models of computation \cite{Fredkin:90}. Most physical theories are based on real numbers. This means that a finite volume of space contains an infinite amount of information, since its state must be specified with real numbers. But if there is an energetic cost to creating information (as there is in most theories), then this implies an infinite amount of energy in a finite space. This is obviously unacceptable; something must bound the information content of space. While such a notion can arise in quantum field theories, CAs start as discrete theories that do not have this problem, and so in many ways they are more satisfying than differential equations as a way to specify governing equations. There is nothing less basic about them than differential equations; which is more ``fundamental'' depends on whether you are solving a problem with a pencil or a computer.

References

Problems

  1. Redo the Finite Element beam-bending problem with a discrete element method, and compare the results. Write a parallel solver that uses a triangular lattice in 2D with a linear elastic nearest-neighbor force law, and use the distance between particles to visualize the strain.
  2. Using a discrete element method, vary the force law to simulate dropping onto a surface bodies that are rigid, elastic, viscous, and brittle.
  3. Extra credit: Implement a particle system that computes.

(c) Neil Gershenfeld 3/2/26