[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.

v_{ij}=v_i-v_j\\
r_{ij}=r_i(0)-r_j(0)\\

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,

m_iv_i+m_jv_j=m_iv_i'+m_jv_j'

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

e=-\frac{v_i'-v_j'}{v_i-v_j}

Therefore, from these,

v_i'=V-e\frac{m_j}{m_i+m_j}v_{ij}\\
v_j'=V+e\frac{m_i}{m_i+m_j}v_{ij}\\

here,

V=\frac{m_iv_i+m_jv_j}{m_i+m_j}

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 $).

• dvf and dv are variables in the source code in Chapter 3. It's hard to understand if it's a code. ..

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

3.1. Collision detection

thermosim.py

    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

thermosim.py

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

            # 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

thermosim.py

            # 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

thermosim.py

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

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

References

• https://github.com/pierrethibault/thermosim -I wanted to do large-scale numerical calculations of rigid spheres, @ NatsukiLab -High-speed simulation of large-scale rigid sphere system, @ yotapoon -Large-scale calculation and speed-up method in rigid disk molecular dynamics simulation, Masaharu Isobe, Physical Properties Research, 72 (1), pp.21-41, 1999-04 -Movement of the center of gravity when the law of conservation of momentum holds -Collision of two balls