From the description in the post and comments, this is basically an $n$-body simulation with (1) collisions and (2) Brownian motion as the driver of motion.
For any $n$-body simulation, you need to evolve the following two differential equations,
\begin{align}
\frac{\mathrm{d}\mathbf{x}_i}{\mathrm{d}t}&=\mathbf{v}_i \\
m\frac{\mathrm{d}\mathbf{v}_i}{\mathrm{d}t}&=\sum_j\mathbf{F}_{ij}
\end{align}
where the $i$ indicates one of the $n$ particles. In many cases, $\mathbf{F}$ is the gravitational force; in other cases, the Lorentz force is used (in which case this is generally called a particle in cell simulation due to added complexities of the background field). Since you are interested in the random motion of the particles in the space, you should be able to use a Langevin equation,
$$m\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}=-\lambda\mathbf{v}+{\boldsymbol\eta}(t)$$
where ${\boldsymbol\eta}(t)$ is the stochastic white noise term. This also can be written as a stochastic differential equation,
$$m\,\mathrm{d}\mathbf{v}=-\lambda\mathbf{v}\mathrm{d}t+\mathrm{d}W$$
where $\mathrm{d}W$ is the Wiener process. Note that I am assuming that the particles don't interact with each other except through collisions, which I'll discuss later; if there are interactions with forces, then those need to be addressed in your force equation.
These can be evolved using pretty simple Euler integration scheme (see also my answer here), so for a single time step, you would have an update like,
for n = 1 to num_particles
x[n] = x[n] + v[n] * dt
v[n] = (1 - lambda * dt / m) * v[n] + dW(dt) / m
where dW(dt)
returns a random normal value with mean 0 and variance dt
.
After each update, you need to check the distance between all of the particles to see if any two are "close enough" to be considered as colliding. This is naively a $\mathcal{O}(n^2)$ comparison, though I imagine more intelligent systems have been developed & could be used. If any two particles are close enough, you need to compute the resulting velocities given your elastic scattering formula.
Likewise, if you are considering your domain to be rigid, then you need to also compute reflections off the walls for particles that are also close to the wall. For periodic boundaries, you would need to check particles that crossed one wall and place it at the other side (e.g., $x_i\sim1.03\to0.03$ when considering a domain of (0,1))
Overall, your simulation outline would be something like,
- Initialize everything ($\mathrm{d}t$, number of particles, positions, velocities, etc)
- While the number of collisions are less than 10:
- Evolve the particles one time step using the Euler update above
- Check for collisions between particles & increment the collision counter as necessary
- Check for collisions with the wall or, alternatively, check for wrapping if using periodic boundaries
- Increment total time by $\mathrm{d}t$ (may not be necessary, but might be useful)
- Output useful details (e.g., current time, iteration, collisions detected, etc) then return to step 2.1
Processing for the mean free path can be done somewhere in the "check for scattering" part, assuming you know the particles' initial positions.