[PYTHON] Artikel, der Ihnen hilft, den Kollisionsalgorithmus für starre Kugeln ein wenig zu verstehen

Ich verwende die [Corpuscular Particle Method (CPM)] von LS-DYNA (http://public.beuth-hochschule.de/~kleinsch/Expl_FEM/2007_CPM_LSDyna.pdf), aber beim Mischen mehrerer Gasarten Ich wollte den Effekt des Verhältnisses der Anzahl der Partikel berücksichtigen (ich frage mich, ob es in Ordnung ist, es zu schreiben, also werde ich den Inhalt weglassen), also habe ich mich mit der Kollision starrer Kugeln befasst, obwohl der Inhalt ziemlich kurz vor dem Ziel liegt. Ich weiß nicht, wie die Nummer lautet (Qiita hat auch Artikel von @yotapoon und @NatsukiLab), aber ich möchte zusammenfassen, was ich verstehe.

1. Was ist ein starrer Kollisionsalgorithmus zwischen Kugel und Scheibe?

** Dies ist eine Methode zur Durchführung der Simulation, wie in der folgenden Abbildung gezeigt. ** (Die Abbildung stammt aus Kinetischer Theorie der Gase von Wikipedia) Es scheint, dass sie auch in der Spielprogrammierung usw. verwendet wird. Es wird auch bei der Analyse der Kinetik von Gasmolekülen und in dem Baltzmann-Verteilungsartikel von Wikipedia als Ergebnis einer Partikelkollision verwendet. Es kann bestätigt werden, dass die Geschwindigkeitsverteilung zur Maxwell-Boltzmann-Verteilung konvergiert. Meine Motivation ist es, die Geschwindigkeitsverteilung von Gasmolekülen zu berechnen und zu verstehen, wie sich die kinetische Energie jedes Gases ändert, wenn das Gas gemischt wird.

2. Annahmen und mathematisches Verständnis

Erstens sind die folgenden beiden Voraussetzungen. Besonders eine Größe ist der Schlüssel zur Realisierung zufälliger Bewegungen.

1. Kugel / Scheibe hat Durchmesser, Geschwindigkeit und Masse
2. Die Kugel / Scheibe macht eine lineare Bewegung mit konstanter Geschwindigkeit = die Geschwindigkeit ändert sich nur bei Kollision

Darüber hinaus wird in diesem Artikel die Quelle erläutert, die auf github [thermosim] hochgeladen wurde (https://github.com/pierrethibault/thermosim). Diese Quelle verwendet eine sogenannte Zeitschritt-Simulation. (Klicken Sie hier für Details. Https://qiita.com/NatsukiLab/items/476e00fea40b86ece31f) Da die gleiche / ähnliche Methode in dem in der Referenz gezeigten Artikel beschrieben wurde, scheint es sich außerdem um eine allgemeine Methode zu handeln.

2.1 Kollisionsurteil

Die Koordinaten der Teilchen * i *. * J * zum Zeitpunkt $ t $ seien $ r (t) $, die Geschwindigkeit sei $ v (t) $, die Masse sei $ m $ und der Durchmesser sei $ d $. Verwenden Sie zunächst die folgende Formel, um eine Liste der Partikel zu erhalten, die derzeit ** bereits ** kollidieren.

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

Die Tatsache, dass sie bereits kollidiert sind, ist, dass die Partikel, die im aktuellen Schritt (Zeit $ t $) in Kontakt gekommen sind, vor dem Kontakt "zurückgespult" werden und der Geschwindigkeitsvektor aus der Geschwindigkeit und den Koordinaten zum Zeitpunkt des Zurückspulens aktualisiert wird. Ich werde. (Es ist garantiert, dass im vorherigen Schritt kein Kontakt besteht.)

Sobald Sie die kollidierenden Partikel kennen, müssen Sie im nächsten Schritt berechnen, wann sie kollidieren. Definieren Sie zunächst die Relativgeschwindigkeit und den relativen Positionsvektor wie folgt. Beachten Sie, dass die Zeit $ t = 0 $ die Zeit des vorherigen Schritts ist.

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

Da sich die Partikel in einer linearen Bewegung mit einer konstanten Geschwindigkeit bewegen, sind die Koordinaten zum Zeitpunkt der Kollision unten gezeigt, unter der Annahme, dass die Zeit vom vorherigen Schritt bis zur Kollision $ t_c $ beträgt.

r_i(t_c)=r_i(0)+v_i(t)・ T._c\\
r_j(t_c)=r_j(0)+v_j(t)・ T._c\\

Von Oben,|r_i(t) - r_j(t)|=\frac{1}{2}(d_i+d_j)Zeit wannt_cWird unten gezeigt.(Da die Zeit des vorherigen Schritts auf 0 gesetzt ist,tcEntspricht dem Zeitunterschied vom vorherigen Schritt zur Kollision.)

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}

Einfach das oben genannte|r_i(t) - r_j(t)|=\frac{1}{2}(d_i+d_j)Auf beiden Seiten quadratischt_cDa es sich um eine quadratische Gleichung für handelt, wird sie aus der Gleichung der Lösung abgeleitet. Es gibt zwei Lösungen:+Lösung(Größere Lösung)Ist die Zeit, die die Partikel benötigen, um nach dem Eindringen in Kontakt zu kommen?

2.2 Geschwindigkeitsaktualisierung nach der Kollision

Auch vor und nach der Kollision basiert die Geschwindigkeit des Schwerpunktsystems auf der Eigenschaft, dass sich der Schwerpunkt des Objekts linear mit konstanter Geschwindigkeit bewegt (http://www.wakariyasui.sakura.ne.jp/p/mech/unndouryhzn/jyuusinnunn.html). Verwenden Sie die relativen Geschwindigkeiten von und i, j, um die Geschwindigkeit nach der Kollision zu berechnen. (Obwohl es gemäß dem Programm hier erklärt wird, ist der Referenzartikel intuitiv leichter zu verstehen.)

Da der Impuls vor und nach der Kollision gespeichert wird,

m_iv_i+m_jv_j=m_iv_i'+m_jv_j'

Es wird sein. Wenn der Abstoßungskoeffizient $ e $ ist

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

Daher aus diesen,

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

Hier,

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

Es bedeutet die Geschwindigkeit des Schwerpunkts.

Die folgende Abbildung zeigt die Beziehung zum Zeitpunkt der vollständigen elastischen Kollision ($ e = 1 $).

• dvf und dv sind Variablen im Quellcode in Kapitel 3. Es ist schwer zu verstehen, ob es sich um einen Code handelt. ..

3. Implementierung (von thermosim.py/def kollidieren in thermosim)

3.1 Kollisionsurteil

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 und self.d sind np.arrays, die Koordinaten bzw. Durchmesser enthalten, und squareform (pdist (self.r)) wird verwendet, um den aktuellen Abstand zwischen Partikeln festzulegen, .5 * np.add. Berechnen Sie die Entfernung, wenn sich die Partikel mit dem Äußeren berühren (self.d, self.d) . Jedes gibt ein Array von nxn zurück, wobei n die Anzahl der Partikel ist. Der Grund, warum unique = (ind1 <ind2) gesetzt ist, ist, dass die Informationen nur die Komponenten des oberen Dreiecks erfordern (mit Ausnahme der Diagonale).

3.2. Berechnung der Koordinaten zum Zeitpunkt der Kollision

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 hat dieselbe Bedeutung und dieselbe Berechnungsmethode wie $ -t_c $ in 2.1. "drc" sollte der relative Positionsvektor zum Zeitpunkt der Kollision sein.

Geschwindigkeitsaktualisierung

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)

Im Kollisionsdiagramm in Kapitel 2 erfahren Sie, was "dvf" bedeutet.

Koordinaten aktualisieren

thermosim.py

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

Ist das dtc okay? Ich denke ein wenig nach.

Verweise

• https://github.com/pierrethibault/thermosim
• Ich wollte umfangreiche numerische Berechnungen für starre Kugeln durchführen, @ NatsukiLab
• Hochgeschwindigkeitssimulation eines großen starren Kugelsystems, @ yotapoon
• Berechnungs- und Beschleunigungsmethode in großem Maßstab in der Simulation der Molekulardynamik starrer Scheiben, Masaharu Isobe, Physical Properties Research, 72 (1), S. 21-41, 1999-04
• Bewegung des Schwerpunkts, wenn das Gesetz der Impulserhaltung gilt
• Kollision zweier Bälle