[PYTHON] Article qui vous aidera à comprendre un peu l'algorithme de collision de sphères rigides

J'utilise la [méthode des particules corpusculaires (CPM)] de LS-DYNA (http://public.beuth-hochschule.de/~kleinsch/Expl_FEM/2007_CPM_LSDyna.pdf), mais lors du mélange de plusieurs types de gaz Je voulais considérer l'effet du rapport du nombre de particules (je vais omettre le contenu car je me demande s'il est correct de l'écrire), alors j'ai étudié la collision de sphères rigides, bien que ce soit bien avant le but. Je ne sais pas quel est le nombre (Qiita a également des articles écrits par @yotapoon et @NatsukiLab), mais je veux résumer ce que je comprends.

1. Qu'est-ce qu'un algorithme de collision sphère / disque rigide?

** Il s'agit d'une méthode pour effectuer la simulation comme indiqué dans la figure ci-dessous. ** (La figure est tirée de la théorie cinétique des gaz de Wikipedia) Il semble qu'elle soit également utilisée dans la programmation de jeux, etc. Il est également utilisé dans l'analyse de la cinétique des molécules de gaz, et dans l'article de Wikipédia sur la distribution Baltzmann, à la suite d'une collision de particules, On peut confirmer que la distribution de vitesse converge vers la distribution de Maxwell-Boltzmann. Ma motivation est de calculer la distribution de vitesse des molécules de gaz et de comprendre comment l'énergie cinétique de chaque gaz change lorsque le gaz est mélangé.

2. Hypothèses et compréhension mathématique

Premièrement, les deux suivants sont les prérequis. Surtout avoir une taille est la clé pour réaliser un mouvement aléatoire.

1. La sphère / disque a un diamètre, une vitesse et une masse
2. La sphère / disque effectue un mouvement linéaire à vitesse constante = la vitesse ne change pas sauf en cas de collision

De plus, cet article expliquera la source qui a été téléchargée sur github appelée thermosim. Cette source utilise une simulation de type Time-Step. (Cliquez ici pour plus de détails. Https://qiita.com/NatsukiLab/items/476e00fea40b86ece31f) De plus, étant donné que la même méthode / méthode similaire a été décrite dans l'article présenté dans la référence, cela semble être une méthode générale.

2.1. Jugement de collision

Soit les coordonnées des particules * i *. * J * au temps $ t $ $ r (t) $, la vitesse soit $ v (t) $, la masse $ m $ et le diamètre $ d $. Tout d'abord, utilisez la formule ci-dessous pour obtenir une liste des particules qui sont ** déjà ** en collision pour le moment.

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

Le fait qu'ils soient déjà entrés en collision est qu'ils «rembobinent» les particules qui sont entrées en contact dans l'étape courante (temps $ t $) avant le contact, et mettent à jour le vecteur vitesse à partir de la vitesse et des coordonnées au moment du rembobinage. Je vais. (Il est garanti qu'il n'y a pas de contact à l'étape précédente)

Une fois que vous connaissez les particules qui entrent en collision, l'étape suivante consiste à calculer le moment où elles sont entrées en collision. Tout d'abord, définissez la vitesse relative et le vecteur de position relative comme suit. Notez que l'heure $ t = 0 $ est l'heure de l'étape précédente.

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

Puisque les particules se déplacent dans un mouvement linéaire à une vitesse constante, les coordonnées au moment de la collision sont indiquées ci-dessous, en supposant que le temps entre l'étape précédente et la collision est $ 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\\

De ce qui précède,|r_i(t) - r_j(t)|=\frac{1}{2}(d_i+d_j)Heure à laquellet_cEst montré ci-dessous.(En fait, puisque l'heure de l'étape précédente est mise à 0,tcEst le même que le décalage horaire entre l'étape précédente et la 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}

Simplement ce qui précède|r_i(t) - r_j(t)|=\frac{1}{2}(d_i+d_j)Carré des deux côtés det_cPuisqu'il s'agit d'une équation quadratique pour, elle est dérivée de l'équation de la solution. Il existe deux solutions,+Solution(Solution plus grande)C'est le temps qu'il faut aux particules pour entrer en contact après avoir pénétré.

2.2 Mise à jour de la vitesse après la collision

Même avant et après la collision, la vitesse du système du centre de gravité est basée sur la propriété que le centre de gravité de l'objet se déplace linéairement à une vitesse constante (http://www.wakariyasui.sakura.ne.jp/p/mech/unndouryhzn/jyuusinnunn.html). Utilisez les vitesses relatives de et i, j pour calculer la vitesse après la collision. (Bien que cela soit expliqué selon le programme ici, l'article de référence est plus facile à comprendre intuitivement)

Puisque l'élan est sauvegardé avant et après la collision,

m_iv_i+m_jv_j=m_iv_i'+m_jv_j'

Ce sera. Si le coefficient de répulsion est $ e $

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

Par conséquent, à partir de ceux-ci,

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}\\

ici,

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

Cela signifie la vitesse du centre de gravité.

Le chiffre de la relation au moment de la collision élastique parfaite ($ e = 1 $) est le suivant.

• dvf et dv sont des variables dans le code source du chapitre 3. Il est difficile de comprendre si c'est un code. ..

3. Mise en œuvre (à partir de thermosim.py/def collide dans thermosim)

3.1. Jugement de collision

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 et self.d sont des np.arrays contenant respectivement des coordonnées et des diamètres, et squareform (pdist (self.r)) est utilisé pour définir la distance actuelle entre les particules, .5 * np.add. Calculez la distance au contact des particules avec l'extérieur (self.d, self.d) ». Chacun renvoie un tableau de nxn, où n est le nombre de particules. ʻUnique = (ind1 <ind2) est parce que l'information ne nécessite que la composante triangulaire supérieure (à l'exclusion des composantes diagonales).

3.2. Calcul des coordonnées au moment de la 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 a la même signification et la même méthode de calcul que $ -t_c $ dans 2.1. drc doit être le vecteur de position relative au moment de la collision.

Mise à jour de la vitesse

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)

Voir le diagramme de collision dans le chapitre 2 pour ce que signifie dvf.

Mettre à jour les coordonnées

thermosim.py

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

Est-ce que «dtc» est correct? Je réfléchis un peu.

Les références

• https://github.com/pierrethibault/thermosim
• Je voulais faire des calculs numériques à grande échelle pour les sphères rigides, @ NatsukiLab
• Simulation à grande vitesse d'un système de sphères rigides à grande échelle, @ yotapoon
• Méthode de calcul et d'accélération à grande échelle dans la simulation de dynamique moléculaire de disque rigide, Masaharu Isobe, Recherche sur les propriétés physiques, 72 (1), pages 21-41, 1999-04
• Mouvement du centre de gravité lorsque la loi de conservation de l'élan tient
• Collision de deux balles