Compressed Particles

Consider a system with a bunch of spherical particles compressed inside a container with fixed walls.

How might we generate a configuration like this?

Well, the particle/particle and particle/wall contact forces are what drives the motion of the particles. For frictionless contact, we can model these forces as coming from a contact potential, $U_c$$U_c$, that depends on the relative positions between particles, $\mathrm{\Delta }{\mathbf{p}}_{ij}:={\mathbf{p}}_j-{\mathbf{p}}_i$$\Delta \bold{p}_{ij} \coloneqq \bold{p}_j - \bold{p}_i$. If we define the potential as a piecewise function

then, we can ensure particles will only interact with their immediate neighbors.

What is $w(q)$$w(q)$?

$w(q)$$w(q)$ is the function that will actually define the contact forces when two particles get close to each other. We would like $w$$w$ to be

1. continuous: $w(1)=0$$w(1) = 0$​

2. differentiable: $w^{\prime }(1)=0$$w'(1) = 0$​​

3. monotonic

4. unbounded near zero: $\underset{q\to 0^{+}}{lim}w(q)=+\mathrm{\infty }$$\lim_{q\rightarrow0^+} w(q) = +\infty$

5. smooth

There are a lot of functions that satisfy these properties, but I'll arbitrarily go with $w(q):=k(q-1)\log (q)$$w(q) \coloneqq k \, (q - 1) \, \log(q)$, which makes our potential look like

Now that we have our individual contact potential, the potential of the whole system is given by the sum of pairwise potentials:

where ${\mathbf{w}}_k,{\hat{\mathbf{n}}}_k$$\bold{w}_k, \bold{\hat{n}}_k$ are {point, unit normal} pairs, respectively, that define the walls. With this in mind, the packed configuration of particles is given by minimizing $U_{\text{total}}$$U_{\text{total}}$ with respect to the particle positions. Here's an animation of numerically minimizing $U_{\text{total}}$$U_{\text{total}}$ with gradient descent applied to an initially-overlapping random configuration of particles:

After a few iterations, we've got a configuration that looks reasonable, with particles overlapping slightly with their neighbors.

Resistor Networks

Now, what if these particles were conductive, and we wanted to measure the effective resistance between the top and bottom walls of the packed configuration? Below, we'll discuss two approaches for generating a resistor network for the particle system.

Particle-Centered

One possible way to describe the current flowing through these packed particles is by associating a voltage with each particle's center, and then connecting neighboring particles by a resistor:

The actual model for calculating the resistance between particles is left to the reader, but it will inevitably depend on a few factors like:

• conductivity of each particle

• radius of each particle

• quality of contact between particles

Although this approach is simple, it has a downside that becomes clearer if we zoom in on a specific area:

If current flows from the lower red particle, through the blue particle, and into the upper red particle, which "path" seems more reasonable? The particle-centered approach requires that the current flow all the way to the center of the blue particle first and turn back around to reach the upper red particle. This means the current would have to travel a further distance through the blue conductor, resulting in a seemingly-larger resistance than the (more natural) direct path.

Contact-Centered

Another way to cast the resistor network problem is to associate the voltages not with the centers of particles, but with the points of contact between particles. Then, the current is free to move directly through the conductive particles to reach another point of contact, avoiding the pathology described above. The (slightly more cluttered) graphical representation of this approach is given below:

As before, there is still some modeling that needs to be done to generate the resistance values for the resistors in this network.

Calculating Effective Resistance

Now that we've come up with a resistor network, how do we actually calculate the effective resistance as measured between the top and bottom walls? Before we start analyzing the entire resistor network, let's briefly review how an individual resistor works.

Ohm's law says that the current $i$$i$ that flows through the resistor is proportional to the voltage difference across it:

That's useful if we knew what the voltages were, but we don't actually know that information yet. However, Kirchhoff's current law says that (by charge conservation), at each node in the circuit, the amount of current flowing in must match the current flowing out

In practice, it is helpful to adopt a sign convention where we denote incoming current with a positive sign, and outgoing current with a negative sign. That way, we can write Kirchhoff's current law as just a single summation

We can write the current-balance equations for each node in matrix form as

Where $\mathbf{S}$$\bold{S}$ is the "conductance matrix" for the network. The conductance matrix for an individual resistor is given by

The conductance matrix for the whole network is just the sum of the conductance matrices for each resistor. An example C++ implementation of this conductance matrix calculation for a resistor network is outlined below

x
struct Resistor {
  int node_1;
  int node_2;
  double value;
};

...
  
// initialize the circuit
std::vector< Resistor > network = {
  ...
};
  
// create a matrix of zeroes, with as many 
// rows/columns as we have nodes in the circuit
Matrix S = zeros(num_nodes, num_nodes);

// assemble the conductance matrix, by having each resistor
// add its own contribution, as derived above
for (auto [i, j, R] : network) {
  S[i][i] -= 1.0 / R;
  S[i][j] += 1.0 / R;
  S[j][i] += 1.0 / R;
  S[j][j] -= 1.0 / R;
}

Once we have the conductance matrix for the full network, we can create a simpler, equivalent resistor network with respect to some subset of the nodes by static condensation. That involves partitioning the nodes into two sets: kept ($k$$k$) and eliminated ($e$$e$) and rearranging the conductance matrix into a 2x2 block form, with blocks corresponding to those kept and eliminated nodes

From here, the reduced conductance matrix is given by

We can apply this to our resistor network problem by taking the bottom wall to be node $0$$0$, the top wall to be node $1$$1$, and all remaining nodes numbered $2,3,\cdots$$2, 3, \cdots$ . Then, we compute the conductance matrix entries like in the C++ code snippet above and evaluate the equation $(1)$$(1)$ above to get ${\mathbf{S}}_r$$\bold{S}_r$. If we inspect the entries of ${\mathbf{S}}_r$$\bold{S}_r$, they will have the form $\left[\begin{array}{cc}-s& s\\ s& -s\end{array}\right]$$\begin{bmatrix}-s & s\\ s & -s\end{bmatrix}$, and the equivalent resistance is $R=\frac{1}{s}$$R = \frac{1}{s}$.

Note: in this example one of the red particles near the top right isn't touching any other particles, so its current balance equation is just $0=0$$0 = 0$, which is trivially satisfied. That means the submatrix ${\mathbf{S}}_{ee}$$\bold{S}_{ee}$ isn't invertible, but it is pseudo-invertible. In this case just replace ${\mathbf{S}}_{ee}^{-1}$$\bold{S}_{ee}^{-1}$ by the Moore-Penrose inverse ${\mathbf{S}}_{ee}^{+}$$\bold{S}_{ee}^{+}$.