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.

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

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

- The sphere / disk has diameter, velocity, and mass.
- 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.

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,

```
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

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

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

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

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

`thermosim.py`

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

Is this `dtc`

okay? I'm thinking a little.

- 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

Recommended Posts