  A rotating superfluid forms an array of quantized vortex lines which
determine its angular velocity. The spasmodic evolution of the array under the
influence of deceleration, dissipation, and pinning forces is thought to be
responsible for the phenomenon of pulsar glitches, sudden jumps in the spin
frequency of rotating neutron stars. We describe and implement an $N$-body
method for simulating the motion of up to 5000 vortices in two dimensions and
present the results of numerical experiments validating the method, including
stability of a vortex ring and dissipative formation of an Abrikosov array.
Vortex avalanches occur routinely in the simulations, when chains of unpinning
events are triggered collectively by vortex-vortex repulsion, consistent with
previous, smaller-scale studies using the Gross-Pitaevskii equation. The
probability density functions of the avalanche sizes and waiting times are
consistent with both exponential and log-normal distributions. We find weak
correlations between glitch sizes and waiting times, consistent with
astronomical data and meta-models of pulsar glitch activity as a
state-dependent Poisson process or a Brownian stress-accumulation process, and
inconsistent with a threshold-triggered stress-release model with a single,
global stress reservoir. The spatial distribution of the effective stress
within the simulation volume is analysed before and after a glitch.
