Python - Simulation de mouvement de particules

Chose que tu veux faire

Simulez le mouvement des atomes comme suit en utilisant le potentiel de Leonard Jones

idéal

image.png http://www.engineering-eye.com/AUTODYN/function/index.html

Potentiel et force interatomiques

Potentiel de Leonard Jones

Leonard Jones Potentiel lorsque la distance entre deux particules est $ r_ {ij} $:

\begin{eqnarray}
U(r_{ij})&=&4ε((\frac{σ}{r_{ij}})^{12}-(\frac{σ}{r_{ij}})^{6})\\
&=&4ε(\frac{σ}{r_{ij}})^{6}((\frac{σ}{r_{ij}})^{6}-1)
\end{eqnarray}

Pour les particules $ N $

\begin{eqnarray}
U(\vec{x_1},...,\vec{x_N})&=&4ε\sum_{i=1}^{N}\sum_{j=i+1}^{N}((\frac{σ}{r_{ij}})^{12}-(\frac{σ}{r_{ij}})^{6})\\
\end{eqnarray}

La force agissant sur l'atome $ i $ sous le potentiel de Leonard Jones est due au gradient à la distance interatomique $ r $.

\begin{eqnarray}
\vec{F}_i&=&-∇U(\vec{x_1},...,\vec{x_N})\\
&=&24ε\sum_{j=1,j≠i}^{N}\frac{1}{r_{ij}^2}(\frac{σ}{r_{ij}})^{6}(1-2(\frac{σ}{r_{ij}})^{6})\vec{r}_{ij}\\
\end{eqnarray}

image.png

https://ja.wikipedia.org/wiki/%E3%83%AC%E3%83%8A%E3%83%BC%E3%83%89-%E3%82%B8%E3%83%A7%E3%83%BC%E3%83%B3%E3%82%BA%E3%83%BB%E3%83%9D%E3%83%86%E3%83%B3%E3%82%B7%E3%83%A3%E3%83%AB

Distance de coupure et division des cellules

Pour réduire la quantité de calcul, un rayon de coupure de $ r_ {cut} $ est fourni, et les forces des particules (sur la grille) plus éloignées de $ r_ {cut} $ sont ignorées car leur contribution est faible.

Divisez l'espace en cubes de $ r_ {cut} × r_ {cut} × r_ {cut} $ et gérez les particules appartenant à chaque grille.

Pour les particules dans une grille spécifique, seule la force des particules dans les grilles adjacentes est calculée.

cell.png

Calcul de la position / vitesse des particules

L'équation cinétique de l'équation différentielle du second ordre est transformée en équation différentielle simultanée du premier ordre et calculée par la méthode Lungekutter.

\begin{eqnarray}
m\frac{d^2\vec{r}}{dt^2}=\vec{F}\\
\end{eqnarray}

↓

\begin{eqnarray}
\frac{d\vec{r}}{dt}&=&\vec{f}(\vec{r},\vec{v},t)=\vec{v}\\
\frac{d\vec{v}}{dt}&=&\vec{g}(\vec{r},\vec{v},t)=\frac{\vec{F}(\vec{r},\vec{v},t)}{m}\\
\end{eqnarray}

Mise à jour par la méthode Rungekutta

\begin{eqnarray}
\vec{k_1} &=& dt×\vec{f}(\vec{r},\vec{v},t)=dt×\vec{v}\\
\vec{l_1} &=& dt×\vec{g}(\vec{r},\vec{v},t)=dt×\frac{\vec{F}(\vec{r},\vec{v},t)}{m}\\
\vec{k_2} &=& dt×\vec{f}(\vec{r}+\frac{\vec{k_1}}{2},\vec{v}+\frac{\vec{l_1}}{2},t+\frac{dt}{2})=dt×(\vec{v}+\frac{\vec{l_1}}{2})\\
\vec{l_2} &=& dt×\vec{g}(\vec{r}+\frac{\vec{k_1}}{2},\vec{v}+\frac{\vec{l_1}}{2},t+\frac{dt}{2})=dt×\frac{\vec{F}(\vec{r}+\frac{\vec{k_1}}{2},\vec{v}+\frac{\vec{l_1}}{2},t+\frac{dt}{2})}{m}\\
\vec{k_3} &=& dt×\vec{f}(\vec{r}+\frac{\vec{k_2}}{2},\vec{v}+\frac{\vec{l_2}}{2},t+\frac{dt}{2})=dt×(\vec{v}+\frac{\vec{l_2}}{2})\\
\vec{l_3} &=& dt×\vec{g}(\vec{r}+\frac{\vec{k_2}}{2},\vec{v}+\frac{\vec{l_2}}{2},t+\frac{dt}{2})=dt×\frac{\vec{F}(\vec{r}+\frac{\vec{k_2}}{2},\vec{v}+\frac{\vec{l_2}}{2},t+\frac{dt}{2})}{m}\\
\vec{k_4} &=& dt×\vec{f}(\vec{r}+\vec{k_3},\vec{v}+\vec{l_3},t+dt)=dt×(\vec{v}+\vec{l_3})\\
\vec{l_4} &=& dt×\vec{g}(\vec{r}+\vec{k_3},\vec{v}+\vec{l_3},t+dt)=dt×\frac{\vec{F}(\vec{r}+\vec{k_3},\vec{v}+\vec{l_3},t+dt)}{m}\\
\vec{r}(t+dt)&=&\vec{r}(t)+\frac{\vec{k_1}+2\vec{k_2}+2\vec{k_3}+\vec{k_4}}{6}\\
\vec{v}(t+dt)&=&\vec{v}(t)+\frac{\vec{l_1}+2\vec{l_2}+2\vec{l_3}+\vec{l_4}}{6}\\
\end{eqnarray}

La position $ \ vec {r} (t + dt) $ à $ t + dt $ ・ La vitesse $ \ vec {v} (t + dt) $ est calculée séquentiellement par la formule ci-dessus.

la mise en oeuvre

Classe de base

import numpy as np
from abc import abstractmethod
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation
import random
import math
import datetime

class Particle():
    def __init__(self, pos = [0.0, 0.0, 0.0], vel = [0.0, 0.0, 0.0], force = [0.0, 0.0, 0.0], mass = 1.0, type = -1):
        self.pos = np.array(pos)
        self.vel = np.array(vel)
        self.force = np.array(force)
        self.mass = mass
        self.type = type

    #Leonard Jones Force agissant entre des particules sous potentiel
    @staticmethod
    def force(pos1, pos2, eps = 1.0, sigma = 1.0):
        #return np.array([0, 0, 0])
        r2 = ((pos2 - pos1)**2).sum()
        if abs(r2) < 0.00001:  
            return np.array([0.0, 0.0, 0.0])
        sig6_div_r6 = (sigma**2 / r2)**3
        #print("force : ", 24 * eps * (1/r2) * sig6_div_r6 * (1 - 2 * sig6_div_r6) * (pos2 - pos1))
        return 24 * eps * (1/r2) * sig6_div_r6 * (1 - 2 * sig6_div_r6) * (pos2 - pos1)
    
class Calculater():
    #Position de la particule après dt(pos)Et vitesse(vel)Est mis à jour par la méthode Rungekutta
    @staticmethod
    def runge_kutta(particle, adjuscent_particles, ps, dt, t, eps, sigma):
        k1 = dt * particle.vel
        l1 = 0
        for p in adjuscent_particles:
            l1 += dt * Particle.force(pos1 = particle.pos, pos2 = p.pos, eps = eps, sigma = sigma) / particle.mass
        l1 += dt * ps.force(particle.pos, particle.vel, particle, t)
        
        k2 = dt * (particle.vel + k1 / 2)
        l2 = 0
        for p in adjuscent_particles:
            l2 = dt * Particle.force(particle.pos + l1/2, p.pos, eps, sigma) / particle.mass
        l2 += dt * ps.force(particle.pos + l1/2, particle.vel + k1/2, particle, t + dt/2)
        
        k3 = dt * (particle.vel + k2 / 2)
        l3 = 0
        for p in adjuscent_particles:
            l3 = dt * Particle.force(particle.pos + l2/2, p.pos, eps, sigma) / particle.mass
        l3 += dt * ps.force(particle.pos + l2/2, particle.vel + k2/2, particle, t + dt/2)
        
        k4 = dt * (particle.vel + k3)
        l4 = 0
        for p in adjuscent_particles:
            l4 = dt * Particle.force(particle.pos + l3, p.pos, eps, sigma) / particle.mass
        l4 += dt * ps.force(particle.pos + l3, particle.vel + k3, particle, t + dt)
        
        particle.pos += (k1 + 2*k2 + 2*k3 + k4)/6
        particle.vel += (l1 + 2*l2 + 2*l3 + l4)/6
        
    #Position de la particule après dt(pos)Et vitesse(vel)Est mis à jour par la méthode Euler
    @staticmethod
    def euler(particle, adjuscent_particles, ps, dt, t, eps, sigma):
        f = 0
        for p in adjuscent_particles:
            f += dt * Particle.force(pos1 = particle.pos, pos2 = p.pos, eps = eps, sigma = sigma) / particle.mass
        f += dt * ps.force(particle.pos, particle.vel, particle, t)
        particle.pos += dt * particle.vel
        particle.vel += dt * f
    
class ParticleSystem():
    # 
    def __init__(self, cutoff_r, NX, NY, NZ, e = 0.9, mass = 1, eps = 1, sigma = 1):
        self.cutoff_r = cutoff_r
        self.NX = NX
        self.NY = NY
        self.NZ = NZ
        self.X_MAX = cutoff_r * NX
        self.Y_MAX = cutoff_r * NY
        self.Z_MAX = cutoff_r * NZ
        self.e = e
        self.mass = mass
        self.eps = eps
        self.sigma = sigma
        self.time = 0.0
        self.prev_timestamp = 0.0
        self.fps = 0.1
        self.particles = []
        self.snapshots = []
        self.fig = plt.figure()
        self.ax = self.fig.add_subplot(111, projection='3d')
        self.gif_output_mode = False
        self.init_particles()
        self.update_cells()
        self.setup_particle_type_dict()
    
    #Trouvez la force agissant sur toutes les particules et la position / vitesse après dt.
    def update(self, dt, calc_method):
        self.time += dt
        #update_forces_on_particle()
        self.update_pos_and_vel(dt, calc_method)
        self.update_cells()
    
    '''
    #Calculez la force agissant sur toutes les particules.
    #Aboli car il est peu probable qu'il soit utilisé avec la méthode Rungekutta.
    def update_forces_on_particle():
        for particle in self.particles:
            particle.force = 0
            #Calculer la force interparticulaire agissant sur la particule
            for p in get_particles_adjuscent_cell(particle):
                f = Particle.force(particle.pos, p.pos)
                particle.force += f
            particle.force += self.force(particle)
    '''
    
    #Calcule la position et le temps après dt de toutes les particules_Calculez avec méthode.
    def update_pos_and_vel(self, dt, calc_method):
        for particle in self.particles:
            calc_method(particle = particle, adjuscent_particles = self.get_particles_adjuscent_cell(particle), ps = self, dt = dt, t = self.time, eps = self.eps, sigma = self.sigma)
            #Faites-le rebondir aux deux extrémités de l'espace.(Coefficient de rebond e)
            if particle.pos[0] < 0:
                #particle.pos[0] += 2*(0 - particle.pos[0])
                particle.pos[0] = 0.0
                particle.vel[0] = -self.e * particle.vel[0]
            if particle.pos[0] > self.X_MAX:
                #particle.pos[0] -= 2*(particle.pos[0] - self.X_MAX)
                particle.pos[0] = self.X_MAX - 0.0001
                particle.vel[0] = -self.e * particle.vel[0]
            if particle.pos[1] < 0:
                #particle.pos[1] += 2*(0 - particle.pos[1])
                particle.pos[1] = 0.0
                particle.vel[1] = -self.e * particle.vel[1]
            if particle.pos[1] > self.Y_MAX:
                #particle.pos[1] -= 2*(particle.pos[1] - self.Y_MAX)
                particle.pos[1] = self.Y_MAX - 0.0001
                particle.vel[1] = -self.e * particle.vel[1]
            if particle.pos[2] < 0:
                #particle.pos[2] += 2*(0 - particle.pos[2])
                particle.pos[2] = 0.0
                particle.vel[2] = -self.e * particle.vel[2]
            if particle.pos[2] > self.Z_MAX:
                particle.pos[2] = self.Z_MAX - 0.0001
                #particle.pos[2] -= 2*(particle.pos[2] - self.Z_MAX)
                particle.vel[2] = -self.e * particle.vel[2]
        if self.gif_output_mode is True:
            self.take_snapshot(save_to_snapshots = True)
    
    #Obtenez la liste des particules contenues dans la cellule contenant la cellule de particule et les 26 cellules adjacentes
    def get_particles_adjuscent_cell(self, particle):
        particles = []
        index = self.get_cellindex_from_pos(particle.pos)
        for i in range(index[0] - 1, index[0] + 1):
            for j in range(index[1] - 1, index[1] + 1):
                for k in range(index[2] - 1, index[2] + 1):
                    try:
                        for p in self.cell[i][j][k] if self.cell[i][j][k] is not None else []:
                            if p is not particle: 
                                particles.append(p)
                    except IndexError:
                        continue
        return particles
                    
    #Calculez l'indice de cellule à partir de la position des particules.
    def get_cellindex_from_pos(self, pos):
        return [int(pos[0] / self.cutoff_r), int(pos[1] / self.cutoff_r), int(pos[2] / self.cutoff_r)]
    
    #Mettre à jour la cellule à laquelle appartient chaque particule
    def update_cells(self):
        self.cell = [[[ [] for i in range(self.NX)] for j in range(self.NY)] for k in range(self.NZ)]
        for p in self.particles:
            try:
                index = self.get_cellindex_from_pos(p.pos)
                self.cell[index[0]][index[1]][index[2]].append(p)
            except Exception as e:
                print("Exception : ", e)
                print("pos : ", p.pos)
                #input()
    
    #Obtenir la liste des particules
    def get_particles(self):
        return self.particles
    
    #Obtenir du temps
    def get_time(self):
        return self.time
    
    #Prenez un instantané de l'état des particules. Il n'est pas affiché.
    def take_snapshot(self, save_to_snapshots = False):
        if self.time - self.prev_timestamp > self.fps:
            self.ax.set_title("Time : {}".format(self.time))
            self.ax.set_xlim(0, self.X_MAX)
            self.ax.set_ylim(0, self.Y_MAX)
            self.ax.set_zlim(0, self.Z_MAX)
            scats = []
            for type, particles in self.particle_dict.items():
                x_list = []
                y_list = []
                z_list = []
                for p in particles[0]:
                    x_list.append(p.pos[0])
                    y_list.append(p.pos[1])
                    z_list.append(p.pos[2])
                scats.append(self.ax.scatter(x_list, y_list, z_list, c=particles[1]))
            if save_to_snapshots is True:
                self.snapshots.append(scats)
                print(len(self.snapshots))
            self.prev_timestamp = self.time
    
    #Afficher l'état des particules
    def show_snapshot(self):
        if self.time - self.prev_timestamp > 0.1:
            self.fig = plt.figure()
            self.ax = self.fig.add_subplot(111, projection='3d')
            self.take_snapshot()
            #plt.savefig('box_gravity_time-{}.png'.format(self.time))
            plt.show()
            self.prev_timestamp = self.time
    
    #Mode d'écriture sur un fichier gif
    def start_output_gif(self, fps = 0.1):
        self.gif_output_mode = True
        self.snapshots = []
        self.take_snapshot(save_to_snapshots = True)
    
    #mode d'écriture du fichier gif désactivé. Écrire dans un fichier
    def stop_output_gif(self, filename = "hoge.gif"):
        print("stop_output_gif : ", len(self.snapshots))
        self.gif_output_mode = False
        ani = animation.ArtistAnimation(self.fig, self.snapshots)
        ani.save(filename, writer='imagemagick')
        self.snapshots = []
    
    #Faites un dictionnaire pour chaque type de particule.
    def setup_particle_type_dict(self):
        self.particle_dict = {}
        for p in self.particles:
            if str(p.type) not in self.particle_dict:
                #Attribuer une couleur aléatoire à chaque type de particule
                self.particle_dict[str(p.type)] = [[], tuple([random.random() for _ in range(3)])]
            self.particle_dict[str(p.type)][0].append(p)
    
    #Calculer l'énergie cinétique de toutes les particules
    def get_kinetic_energy(self):
        return  sum([(p.mass*np.array(p.vel)**2).sum()/2 for p in self.particles])
    
    #Initialisation des particules
    @abstractmethod    
    def init_particles(self):
        raise NotImplementedError()
    
    #Le pouvoir de fonctionner en tant que système(Gravité etc.)
    @abstractmethod
    def force(self, pos, vel, particle, t):
        raise NotImplementedError()

Système d'armée de particules en forme de boîte sous gravité

#Système d'armée de particules en forme de boîte sous gravité
class BoxGravitySystem(ParticleSystem):
    def __init__(self, cutoff_r, NX, NY, NZ, e, mass, eps, sigma):
        super().__init__(cutoff_r, NX = NX, NY = NY, NZ = NZ, e = e, mass = mass, eps = eps, sigma = sigma)
        
    #Génère une armée de particules 10x10x10 en forme de boîte en bas. La vitesse initiale est réglée pour aller en diagonale vers le haut.
    def init_particles(self):
        for x in np.linspace(0.1*self.X_MAX, 0.2*self.X_MAX, 10):
            for y in np.linspace(0.1*self.Y_MAX, 0.2*self.Y_MAX, 10):
                for z in np.linspace(0.1*self.Z_MAX, 0.2*self.Z_MAX, 10):
                    self.particles.append(Particle(pos = [x, y, z], vel=[0.1*self.X_MAX, 0.05*self.Y_MAX, 0.5*self.Z_MAX], mass = self.mass))
    
    #la gravité
    def force(self, pos, vel, particle, t):
        return np.array([0.0, 0.0, -particle.mass * 9.8])

if __name__ == '__main__':
    cutoff_r = 1.0
    e = 0.1
    mass = 1.0
    eps = 1
    sigma = 0.5
    dt = 0.001
    system = BoxGravitySystem(cutoff_r = cutoff_r, NX = 100, NY = 100, NZ = 100, e = e, mass = mass, eps = eps, sigma = sigma)
    time = 0
    system.start_output_gif(fps = 0.1)
    while time <= 5:
        system.update(dt = dt, calc_method = Calculater.runge_kutta)
        time = system.get_time()
        print("Time : ", time)
        #system.show_snapshot()
    system.stop_output_gif(filename = "BoxGravitySystem_cutoffr-{}_mass-{}_eps-{}_sigma-{}_dt-{}.gif".format(cutoff_r, mass, eps, sigma, dt))

Grande masse (pas d'interaction interatomique)

gravity_mass-50.gif

Un système pour enfoncer des balles dans le mur

#Un système pour enfoncer des balles dans le mur
class BulletWallSystem(ParticleSystem):
    def __init__(self, cutoff_r, NX, NY, NZ, e, mass, eps, sigma):
        super().__init__(cutoff_r, NX = NX, NY = NY, NZ = NZ, e = e, mass = mass, eps = eps, sigma = sigma)
    
    #Sphère avec vitesse initiale(balle)Et mettre en place un mur fixe
    def init_particles(self):
        #Balle(balle)produire
        bullet_center = [0.35*self.X_MAX, 0.5*self.Y_MAX, 0.5*self.Z_MAX]
        for i in range(200):
            r = 0.05 * self.X_MAX * (random.random() - 0.5) * 2
            phi = 2 * np.pi * random.random()
            theta = np.pi * random.random()
            self.particles.append(Particle(pos = [bullet_center[0] + r*np.sin(theta)*np.cos(phi), bullet_center[1] + r*np.sin(theta)*np.sin(phi), bullet_center[2] + r*np.cos(theta)], vel=[0.2*self.X_MAX, 0, 0], mass = self.mass, type = 1))
        #Génération de mur
        for x in np.linspace(0.49*self.X_MAX, 0.50*self.X_MAX, 2):
            for y in np.linspace(0.3*self.Y_MAX, 0.7*self.Y_MAX, 30):
                for z in np.linspace(0.3*self.Z_MAX, 0.7*self.Z_MAX, 30):
                    self.particles.append(Particle(pos = [x, y, z], vel=[0.0, 0.0, 0.0], mass = self.mass, type = 2))
                    
    #Pas de force externe
    def force(self, pos, vel, particle, t):
        return np.array([0.0, 0.0, 0])

if __name__ == '__main__':
    cutoff_r = 1.0
    e = 0.1
    mass = 1.0
    eps = 1
    sigma = 0.5
    dt = 0.001
    system = BulletWallSystem(cutoff_r = cutoff_r, NX = 100, NY = 100, NZ = 100, e = e, mass = mass, eps = eps, sigma = sigma)
    time = 0
    system.start_output_gif(fps = 0.1)
    while time <= 5:
        system.update(dt = dt, calc_method = Calculater.runge_kutta)
        time = system.get_time()
        print("Time : ", time)
        #system.show_snapshot()
    system.stop_output_gif(filename = "BoxGravitySystem_cutoffr-{}_mass-{}_eps-{}_sigma-{}_dt-{}.gif".format(cutoff_r, mass, eps, sigma, dt))

cutoffr：1.0、ε：0.001、σ：0.1、dt：0.002 BulletWallSystem_cutoffr-1.0_mass-1.0_eps-0.001_sigma-0.1_dt-0.002_2.gif

cutoffr：1.0、ε：0.02、sigma：0.5 BulletWallSystem_cutoffr-1.0_mass-1.0_eps-0.02_sigma-0.5.gif

cutoffr：1.0、ε：0.02、sigma：1.3、dt-0.001 BulletWallSystem_eps-0.02_sigma-1.3_dt-0.001.gif

cutoffr：1.0、ε：0.3、σ：1.3、dt：0.01 BulletWallSystem_eps-0.3_sigma-1.3_dt-0.01.gif

cutoffr：1.0、ε：0.1、σ：1.2、dt：0.01 BulletWallSystem_eps-0.1_sigma-1.2_dt-0.01_2.gif

cutoffr：1.0、ε：0.1、σ：1.5、dt：0.001 BulletWallSystem_eps-0.1_sigma-1.5_dt-0.001.gif

cutoffr：1.0、ε:1、σ：1.5、dt：0.001 BulletWallSystem_eps-1_sigma-1.5_dt-0.001.gif

Un système qui reçoit une force qui tourne comme une tornade

#Un système qui reçoit une force qui tourne comme une tornade
class TornadeSystem(ParticleSystem):
    def __init__(self, cutoff_r, NX, NY, NZ, e, mass, eps, sigma):
        super().__init__(cutoff_r, NX = NX, NY = NY, NZ = NZ, e = e, mass = mass, eps = eps, sigma = sigma)
        
    #Créez 3000 particules à des positions aléatoires
    def init_particles(self):
        for _ in range(3000):
            x = self.X_MAX * random.random()
            y = self.Y_MAX * random.random()
            z = self.Z_MAX * random.random()
            self.particles.append(Particle(pos = [x, y, z], vel=[0.0, 0.0, 0.0], mass = self.mass, type = 1))
        
    #Une force qui tourne autour de l'axe z(rot(r)∝[0,0,1]Pouvoir de devenir) × z
    def force(self, pos, vel, particle, t):
        return 0.05 * pos[2] * np.array([-(pos[1] - self.Y_MAX/2), pos[0] - self.X_MAX/2 , 0.0])# - 0.5 * np.array(vel)

if __name__ == '__main__':
    cutoff_r = 1.0
    e = 0.1
    mass = 1.0
    eps = 1
    sigma = 0.5
    dt = 0.001
    system = TornadeSystem(cutoff_r = cutoff_r, NX = 100, NY = 100, NZ = 100, e = e, mass = mass, eps = eps, sigma = sigma)
    time = 0
    system.start_output_gif(fps = 0.1)
    while time <= 5:
        system.update(dt = dt, calc_method = Calculater.runge_kutta)
        time = system.get_time()
        print("Time : ", time)
        #system.show_snapshot()
    system.stop_output_gif(filename = "BoxGravitySystem_cutoffr-{}_mass-{}_eps-{}_sigma-{}_dt-{}.gif".format(cutoff_r, mass, eps, sigma, dt))

TornadeSystem_cutoffr-1.0_mass-1.0_eps-0.1_sigma-0.1_dt-0.01_bak.gif

TornadeSystem_cutoffr-1.0_mass-1.0_eps-1_sigma-0.5_dt-0.005.gif

Travail en échec TornadeSystem_cutoffr-1.0_mass-1.0_eps-0.001_sigma-0.1_dt-0.01_bak2.gif

Résumé

L'ajustement des paramètres ne fonctionne pas

idéal

image.png

réalité

BulletWallSystem_eps-0.1_sigma-1.2_dt-0.01_2.gif

Recommended Posts

Python - Simulation de mouvement de particules
Simulation d'escalator