Mise en œuvre du champ électrique et magnétique dans une simulation 2D de particules libres

Question

Je cherche à introduire l'électrodynamique dans ma simulation de gaz parfaits, qui est simplement un groupe de particules à l'intérieur d'un rectangle. Modifier il semble que j'utilisais les mauvaises équations, et c'était la cause de mon problème initial. J'utilise l'environnement Python 3.5 IDLE, numpy pour le vecteur mathématique et pygame pour la visualisation. Les équations que je utilise sont: electric field et magnetic field

Jusqu'à présent, si bon, mon code est en cours d'exécution. Mais pour recréer la physique de manière réaliste, les champs doivent être étendus avec la vitesse c. En outre, lors du calcul de la force, je devrais utiliser les données du temps retardé, mais il est réécrit à chaque nouvelle image. Puis-je stocker ces données ou est-ce lourd à la mémoire? Quelle méthode suggérez-vous?

Le code utilise la classe Ball dont les instances sont les particules. Il a une méthode associée qui calcule E(r) et B(r). Pour calculer la force agissant sur la particule i, je les ajoute tous ensemble à la position i et calculer la force agissant sur i.

Voici le code, où «masse» sert de rayon de la particule ainsi que de charge et la dimension z est simplement ajouté pour le produit croisé. Le code correspondant est celui qui est affecté par la déclaration if electrodynamics == true:

import pygame 
import time 
import random as r 
import numpy as np 

pygame.init() 

display_width = 800 
display_height = 600 

electrodynamics = True 
show_field = False 

Ball_num = 7 
c = 10 
red = (255,0,0) 
white = (255,255,255) 
gameDisplay = pygame.display.set_mode((display_width,display_height)) 
pygame.display.set_caption('bouncy') 
clock = pygame.time.Clock() 

def normalized(a): 
    b = np.linalg.norm(a) 
    if not(b == 0): 
     return a/b 
    else: 
     return np.zeros((1,3), dtype = float) 

def tot_EM_field_at_charge(charges, charge): 

    EM = np.array([[0.,0.,0.],[0.,0.,0.]], dtype = float) 

    for q in charges: 
     EM = EM + q.EM_field(charge.position) 

    return EM 

def Force_on_bally(field, charge): 
    c = np.cross(charge.velocity, field[1]) 
    force = charge.mass*(field[0] + np.cross(charge.velocity, field[1])) 
    return force 

class arrow(object): 

    def __init__(self, length, x, y, charges): 


     self.length = length 
     self.position = np.array([x,y,0]) 

     self.field = tot_EM_field_at_charge(charges, self) 
     self.field_mag= np.linalg.norm(self.field[0]) 
     if self.field_mag == 0:    
      self.field_direction = np.zeros(3,) 
     else:  
      self.field_direction = self.field[0]/self.field_mag 

     self.position_end() 
     self.color() 
     self.show() 
    def position_end(self): 

     self.position_2 = self.position + self.field_direction * self.length 

     return self.position_2 

    def color(self): 
     self.color = self.field_mag 
     if self.color < 0.05: 
      self.color = (46,120,255) 
     elif self.color < 0.1: 
      self.color = (147,145,252) 
     elif self.color < 0.3: 
      self.color = (249,23,28) 
     elif self.color < 0.6: 
      self.color =(251,139,33) 
     elif self.color < 1: 
      self.color = (255,255,127) 
     else: 
      self.color = (255,255,255) 


     return self.color 
    def show(self): 
     pygame.draw.line(gameDisplay, self.color,(int(self.position[0]), int(self.position[1])), (int(self.position_2[0]), int(self.position_2[1]))) 








class Ball(object): 
    def __init__(self, x, y, m, c, v_x, v_y): 
     self.position = np.array([x,y,0], dtype = float) 
     self.position_2 = np.array([x,y,0], dtype = float) 
     self.velocity = np.array([v_x, v_y,0], dtype = float) 
     self.velocity_2 = np.array([v_x, v_y,0], dtype = float) 
     self.acceleration = np.array([0.,0.,0.], dtype = float) 
     self.mass = m 
     self.color = c 

    def acceleration_compute(self,force): 
     a = force/self.mass 
     self.acceleration += a 

    def move(self): 
     self.velocity += self.acceleration 
     self.position += self.velocity 
     self.acceleration *= 0 

    def show(self): 
     pygame.draw.circle(gameDisplay, self.color, [int(self.position[0]), int(self.position[1])], self.mass) 

    def Edgelord(self): 
     if ((self.position[0] + self.velocity[0] >= display_width-self.mass) and self.velocity[0] > 0): 
      self.velocity[0] *= -1 
      self.position[0] = display_width - self.mass + self.velocity[0] 


     elif ((self.position[0] + self.velocity[0] - self.mass <= 0) and self.velocity[0] < 0): 

      self.velocity[0] *= -1 
      self.position[0] = self.mass + self.velocity[0] 


     elif ((self.position[1] + self.velocity[1] >= display_height - self.mass) and self.velocity[1] > 0): 

      self.velocity[1] *= -1 
      self.position[1] = display_height - self.mass + self.velocity[1] 


     elif ((self.position[1] + self.velocity[1] - self.mass <= 0) and self.velocity[1] < 0): 

      self.position[1] = self.mass -self.velocity[1] 
      self.velocity[1] *= -1 

    def EM_field(self, R): 
     radius = np.linalg.norm(R - self.position) 
     if radius != 0: 
      unitradius = (R - self.position)/radius 
     else: 
      unitradius = np.zeros(3,) 

     if np.linalg.norm(radius) != 0 and np.dot(unitradius, self.velocity)!=1: 
      charge  = self.mass/(1 - np.dot(unitradius, self.velocity) ** 3) 


      if radius < self.mass: 
       radius = self.mass 

      radius2  = radius ** 2 

      velocity_in_c = self.velocity/c 

      oneMinusV2 = 1 - np.dot(velocity_in_c, velocity_in_c) 
      uMinusV  = unitradius - velocity_in_c 
      aCrossUmV = np.cross(uMinusV, self.acceleration); 
      Eleft  = (oneMinusV2 * (unitradius - velocity_in_c))/radius2 
      Eright  = np.cross(unitradius, aCrossUmV)/(radius*c**2) 
      E   = charge * (Eleft - Eright) 
      #E = np.zeros(3,) 
      B   = np.cross(unitradius/c, ((charge*c**2) * (Eleft - Eright))) 

      EM_field = np.array([E,B], dtype = float) 
     else: 
      EM_field = np.zeros((2,3), dtype = float) 

     return EM_field 

ballys = [] 

for i in range(Ball_num): 
    #ballys.insert(i, Ball(r.randrange(300,display_width - 5, 10),r.randrange(200,display_height/2,1) , r.randrange(5,10,1),(r.randint(1,255),r.randint(1,255),r.randint(1,255)), r.randint(-200,200)/1000, r.randint(-200,200)/1000)) 
    ballys.insert(i, Ball(200 + i*50, 220 + i*20 , 10,(r.randint(1,255),r.randint(1,255),r.randint(1,255)),0, 0)) 
#ballys.append(Ball(300 + 50, 300, 10,(r.randint(1,255),r.randint(1,255),r.randint(1,255)),10, 0)) 

up = np.zeros(3,) 
down = np.zeros(3,) 
right = np.array([0.,0.,0.]) 
left = np.array([0.,0.,0.]) 
grav = np.array([0.,0.1,0.]) 
repulsion = np.array([0.,0.,0.]) 

crashed = False 

while not crashed : 

    for event in pygame.event.get(): 
     if event.type == pygame.QUIT: 
      pygame.quit() 
      quit() 

     if event.type == pygame.KEYDOWN: 
      if event.key == pygame.K_LEFT: 
       left = np.array([-0.1, 0.,0.]) 
      if event.key == pygame.K_RIGHT: 
       right = np.array([0.1,0.,0.]) 
      if event.key == pygame.K_DOWN: 
       down = np.array([0.,0.1,0.]) 
      if event.key == pygame.K_UP: 
       up = np.array([0.,-0.1,0.]) 

     if event.type == pygame.KEYUP: 
      if event.key == pygame.K_LEFT or event.key == pygame.K_RIGHT or event.key == pygame.K_UP or event.key == pygame.K_DOWN : 
       right = np.array([0.,0.,0.]) 
       left = np.array([0.,0.,0.]) 
       down = np.zeros(3,) 
       up = np.zeros(3,) 

    gameDisplay.fill(white) 

    if show_field == True: 
     for i in range(display_width//20): 
      for j in range(display_height//20): 
       arry = arrow(8, 10 + i*20, 10 + j*20 , ballys) 



    if electrodynamics == True: 
     for bally in ballys: 
      bally.acceleration_compute(Force_on_bally(tot_EM_field_at_charge(ballys, bally), bally)) 

    for i, bally in enumerate(ballys): 

     #if electrodynamics == True: 
      # bally.acceleration_compute(Force_on_bally(tot_EM_field_at_charge(ballys, bally), bally)) 
     bally.Edgelord() 

     bally.acceleration_compute(up) 
     bally.acceleration_compute(down) 
     bally.acceleration_compute(right) 
     bally.acceleration_compute(left) 

     #ballys[i].acceleration_compute(grav * ballys[i].mass) 

     for bally2 in ballys[i+1:]: 

      #checks collisions 
      if np.linalg.norm(bally.position - bally2.position) <= bally.mass + bally2.mass : 

       bally.velocity_2 = (bally.mass * bally.velocity + bally2.mass * bally2.velocity + bally2.mass *(bally2.velocity - bally.velocity))/ (bally.mass + bally2.mass) 
       bally2.velocity_2 = (bally.mass * bally2.velocity + bally.mass * bally.velocity + bally.mass *(bally.velocity - bally2.velocity))/ (bally2.mass + bally.mass) 

       #prevents balls getting stuck in each other and assignes new velocitys 
       if not(np.linalg.norm(bally.position + bally.velocity_2 - (bally2.position + bally2.velocity_2)) <= bally.mass + bally2.mass): 
         bally.velocity = bally.velocity_2 
         bally2.velocity = bally2.velocity_2 

     bally.Edgelord() 

     bally.move() 
     bally.show() 

    pygame.display.update() 
    clock.tick(60) 

pygame.quit() 
quit() 

Answer 1

En regardant votre code, je supposais le problème venait des calculs 1/r et 1/r². En effet, si vous divisez par une valeur vraiment petite, vous finissez par obtenir des valeurs vraiment importantes, ce qui peut faire exploser votre simulation.

Une petite valeur de r se produit lorsque vos particules ont une vélocité suffisante pour se rapprocher l'une de l'autre. Vous résolvez alors la collision entre eux, mais le mal est déjà fait. Vous calculez les deux forces, les intégrez pour la nouvelle vitesse, puis résolvez la collision et les vitesses inverses si nécessaire, mais les vitesses sont déjà énormes.

Pour valider mon hypothèse, je viens d'essayer d'empêcher le radius utilisé dans le calcul trop grande, par exemple:

if radius < ball.mass: 
    radius = ball.mass 

J'ai couru la simulation avec 11 particules (qui explosait dans comme 1 ou 2 minutes sans cela). La simulation a fonctionné pendant environ 15 minutes et n'a pas encore explosé (et ne le sera pas). BTW, pour s'assurer que le problème ne venait pas de ce que vous avez observé dans votre édition, j'ai commenté la ligne ballys[i].acceleration_compute(grav * ballys[i].mass). Il semble que la simulation explose toujours sans (et est également stable en vous assurant que le rayon ne soit pas trop petit)

EDIT: Here is the code I ran and from which I observed a correct behaviour