[PYTHON] An article that gives you a little understanding of the rigid sphere collision algorithm

I am using LS-DYNA's Corpuscular Particle Method (CPM), but when mixing multiple gases I wanted to consider the effect of the ratio of the number of particles (I'm wondering if I can write it, so I'll omit the content), so I studied about rigid sphere collision, which is quite before the goal. I don't know what the number is (Qiita also has articles written by @yotapoon and @NatsukiLab), but I want to summarize what I understand.

1. What is a rigid sphere / disk collision algorithm?

** This is a simulation method as shown below. ** (The figure is taken from Wikipedia's Kinetic theory of gases) It seems that it is also used in game programming etc. It is also used in the analysis of kinetic theory of gas, and in Wikipedia's Baltzmann distribution article, as a result of particle collision, particles It can be confirmed that the velocity distribution converges to the Maxwell-Boltzmann distribution. My motivation is to calculate the velocity distribution of gas molecules and understand how the kinetic energy of each gas changes when the gas is mixed.


2. Assumptions and mathematical understanding

First, the following two assumptions are made. Especially having a size is the key to realizing random movement.

  1. The sphere / disk has diameter, velocity, and mass.
  2. The sphere / disk makes a constant velocity linear motion = the velocity does not change except for collision.

In addition, this article will explain the source that was uploaded to github called thermosim. This source employs a so-called Time-Step type simulation. (Click here for details. Https://qiita.com/NatsukiLab/items/476e00fea40b86ece31f) In addition, since the same / similar method was described in the article shown in the reference, it seems to be a general method.

2.1. Collision detection

Let the coordinates of the particle * i *. * J * at time $ t $ be $ r (t) $, the velocity be $ v (t) $, the mass be $ m $, and the diameter be $ d $. First, use the formula below to get a list of particles that are ** already ** colliding at the moment.

|r_i(t) - r_j(t)|\leq\frac{1}{2}(d_i+d_j)

The fact that they have already collided is that the particles that came into contact in the current step (time $ t $) are "rewound" before contact, and the velocity vector is updated from the velocity and coordinates at the time of rewinding. I will. (It is guaranteed that there is no contact in the previous step)

Once you know the particles that are colliding, the next step is to calculate when they collided. First, define the relative velocity and relative position vector as follows. Note that the time $ t = 0 $ is the time of the previous step.


Since the particles move in a constant velocity linear motion, the coordinates at the time of collision are shown below, where the time from the previous step to the collision is $ t_c $.

r_i(t_c)=r_i(0)+v_i(t)・ T_c\\
r_j(t_c)=r_j(0)+v_j(t)・ T_c\\

From the above,|r_i(t) - r_j(t)|=\frac{1}{2}(d_i+d_j)Time whent_cIs shown below.(In fact, since the time of the previous step is set to 0,tcIs the same as the time difference from the previous step to the collision.)

t_c=\frac{-r_{ij}・ V_{ij}-\sqrt{(r_{ij}・ V_{ij})^2-v_{ij}^2(r_{ij}^2-0.25(d_i+d_j)^2)}}{v_{ij}^2}

Simply the above|r_i(t) - r_j(t)|=\frac{1}{2}(d_i+d_j)If you square both sides oft_cSince it is a quadratic equation for, it is derived from the equation of the solution. There are two solutions,+Solution(Larger solution)Is the time it takes for the particles to come into contact after penetrating.

2.2 Post-collision speed update

Even before and after the collision, the velocity of the center of gravity system is based on the property that the center of gravity of the object moves linearly at a constant velocity (http://www.wakariyasui.sakura.ne.jp/p/mech/unndouryhzn/jyuusinnunn.html). The relative velocity of and i, j is used to calculate the velocity after the collision. (Although it is explained according to the program here, the article in the bibliography is easier to understand intuitively.)

Since the momentum is saved before and after the collision,


It will be. If the coefficient of restitution is $ e $,


Therefore, from these,




It means the speed of the center of gravity.

The figure below shows the relationship at the time of a complete elastic collision ($ e = 1 $). image.png

3. Implementation (from thermosim.py/def collide in thermosim)

3.1. Collision detection


    from scipy.spatial.distance import squareform,pdist
    # Find colliding particles
    D = squareform(pdist(self.r))
    ind1, ind2 = np.where(D < .5*np.add.outer(self.d, self.d))
    unique = (ind1 < ind2)
    ind1 = ind1[unique]
    ind2 = ind2[unique]

self.r and self.d are np.arrays containing coordinates and diameters, respectively, and squareform (pdist (self.r)) is used to set the current distance between particles, .5 * np.add. Calculate the distance when particles touch with outer (self.d, self.d) . If the number of particles is n, then an array of nxn is returned. ʻUnique = (ind1 <ind2)` is because the information requires only the upper triangular component (excluding the diagonal).

3.2. Calculation of coordinates at the time of collision


            ru = np.dot(dv, dr)/ndv
            ds = ru + sqrt(ru**2 + .25*(d1+d2)**2 - np.dot(dr, dr))
            if np.isnan(ds):

            # Time since collision
            dtc = ds/ndv

            # New collision parameter
            drc = dr - dv*dtc

dtc has the same meaning and the same calculation method as $ -t_c $ in 2.1. drc should be the relative position vector at the time of collision.

Speed update


            # Center of mass velocity
            vcm = (m1*v1 + m2*v2)/(m1+m2)

            # Velocities after collision
            dvf = dv - 2.*drc * np.dot(dv, drc)/np.dot(drc, drc)
            v1f = vcm - dvf * m2/(m1+m2)
            v2f = vcm + dvf * m1/(m1+m2)

See the collision diagram in Chapter 2 for what dvf means.

Update coordinates


            # Backtracked positions
            r1f = r1 + (v1f-v1)*dtc
            r2f = r2 + (v2f-v2)*dtc

Is this dtc okay? I'm thinking a little.


Recommended Posts

An article that gives you a little understanding of the rigid sphere collision algorithm
An article that just tries a little HTTP request with the curl command
A memorandum of filter commands that you might forget in an instant
Let's do clustering that gives a nice bird's-eye view of the text dataset
The story of Django creating a library that might be a little more useful
What is a rational decision that maximizes the chances of encountering an "ideal home"?
A story that reduces the effort of operation / maintenance
[Python] A program that counts the number of valleys
A script that takes a snapshot of an EBS volume
Make a BOT that shortens the URL of Discord
[PyTorch] A little understanding of CrossEntropyLoss with mathematical formulas
LINE Bot that notifies you of the stocks of interest
Generate that shape of the bottom of a PET bottle
A story that analyzed the delivery of Nico Nama.
[Python] A program that compares the positions of kangaroos.
Create a shape on the trajectory of an object
An example of a mechanism that returns a prediction by HTTP from the result of machine learning
A Python script that allows you to check the status of the server from your browser