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,

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

What is

? $w(q)$

continuous:

$w(1)=0$ differentiable:

${w}^{\prime}(1)=0$ monotonic

unbounded near zero:

$\underset{q\to {0}^{+}}{lim}w(q)=+\mathrm{\infty}$ smooth

There are a lot of functions that satisfy these properties, but I'll arbitrarily go with

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

where

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

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.

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.

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.

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

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

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 (

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

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

, which is trivially satisfied. That means the submatrix $0=0$ isn't invertible, but it is ${\mathbf{S}}_{ee}$ pseudo-invertible. In this case just replaceby the Moore-Penrose inverse ${\mathbf{S}}_{ee}^{-1}$ . ${\mathbf{S}}_{ee}^{+}$