Collision sans frottement entre N disques

1. Introduction

On considère N disques en mouvement sans frottement sur un plan horizontal (problème bi-dimensionnel). La collision entre les disques se fait sans frottement, avec un coefficient de restitution variable. Les équations obtenues sont valables aussi pour deux sphères en collision (sans frottement) dans l'espace à trois dimensions (modèle des sphères dures).

Lors d'un choc sans frottement, la vitesse de glissement entre les deux disques n'intervient pas dans les équations. C'est pourquoi il n'est pas nécessaire de tenir compte de la rotation des solides.

Le cas de la collision de deux disques est traité dans Collision sans frottement entre deux disques. La détermination de la position des disques au moment de la collision est aussi traitée dans cette page.

On s'intéresse ici au cas d'une collision simultanée entre N_d disques. En pratique, les collisions sont considérées comme simultanées si leur différence d'instant est en valeur absolue inférieure à une tolérance ε qui doit être choisie en fonction du problème traité.

On considère un groupe de N_d disques. Chaque disque de ce groupe est soit en collision avec au moins deux autres disques, soit en collision avec un disque qui fait partie du groupe (pour une des deux raisons).

Le groupe comporte N_p paires de disques qui sont en collision. Chaque disque intervient dans au moins une paire.

2. Mise en équation

Les disques du groupe de collision sont indexés de i=0 à i=N_d-1. Dans le programme de simulation, les disques ont un indice global qu'il conviendra de convertir en indice du groupe (et vice-versa) au moyen de deux tables associatives.

Les vitesses ${\vec{V}}_{i}$ des disques avant la collision sont connue, ainsi que les positions de leur centre ${\vec{x}}_{i}$ . Il faut calculer les vitesses $\vec{V}'_{i}$ après la collision. On introduit pour cela les impulsions ${\vec{p}}_{i, j}$ échangées entre disques pour chaque paire du groupe. Les vitesses et les impulsions constituent les N_d+N_p vecteurs inconnus. Dans le cas bi-dimensionnel des disques, cela fait 2(N_d+N_p) inconnues scalaires. Pour des sphères à trois dimensions, on a 3(N_d+N_p) inconnues.

Voyons à présent les équations à résoudre. Chaque disque du groupe reçoit une impulsion des disques qui sont en paire avec lui (ensemble noté P(i)) :

m_{i} (\vec{V}'_{i} - {\vec{V}}_{i}) = \sum_{j \in P (i)} {\vec{p}}_{i, j} (1)

Pour chaque paire, on a

{\vec{p}}_{i, j} = - {\vec{p}}_{j, i} (2)

ce qui assure la conservation globale de la quantité de mouvement pour le groupe. Notons toutefois qu'une seule impulsion inconnue est introduite pour chaque paire. Pour la paire (i,j) (dans cet ordre), on introduit l'impulsion ${\vec{p}}_{i, j}$ qui est reçue par le premier disque de la paire.

Les équations suivantes expriment l'absence de frottement, c'est-à-dire le fait que les impulsions n'ont pas de composantes tangentielles :

{\vec{p}}_{i, j} \cdot {\vec{T}}_{i, j} = 0 (3)

Le vecteur ${\vec{T}}_{i, j}$ est un vecteur tangent aux disques au point de contact. Pour le problème des sphères en trois dimensions, il y aura deux vecteurs tangents et donc deux équations scalaires pour chaque paire.

Les composantes normales des vitesses relatives sont transformées avec un coefficient de restitution e :

(\vec{V}'_{i} - \vec{V}'_{j}) \cdot {\vec{N}}_{i, j} = - e ({\vec{V}}_{i} - {\vec{V}}_{j}) \cdot {\vec{N}}_{i, j} (4)

Le vecteur normal aux surfaces lors de la collision est simplement le vecteur position relative des deux centres des sphères :

{\vec{N}}_{i, j} = {\vec{x}}_{i} - {\vec{x}}_{j} (5)

Le vecteur tangent (2D) ou les deux vecteurs tangents (3D) sont construits à partir du vecteur normal.

Le décompte des équations scalaires pour un problème 2D (disques) se fait comme suit : 2N_d équations scalaires (1) plus N_p équation scalaires (3) plus N_p équations scalaires (4), ce qui fait bien au total 2(N_d+N_p) équations scalaires.

La forme suivante reprend ces équations en introduisant les notations x et y pour les axes et en plaçant les termes constants à droite :

\begin{matrix} m_{i} v'_{i, x} - \sum_{j \in P (i)} p_{i, j, x} = m_{i} v_{i, x} \\ m_{i} v'_{i, y} - \sum_{j \in P (i)} p_{i, j, y} = m_{i} v_{i, y} \\ p_{i, j, x} T_{i, j, x} + p_{i, j, y} T_{i, j, y} = 0 \\ v'_{i, x} N_{i, j, x} - v'_{j, x} N_{i, j, x} + v'_{i, y} N_{i, j, y} - v'_{j, y} N_{i, j, y} = - e (v_{i, x} N_{i, j, x} - v_{j, x} N_{i, j, x} + v_{i, y} N_{i, j, y} - v_{j, y} N_{i, j, y}) \end{matrix}

Les deux premières équations doivent être écrites pour chaque disque; les deux dernières pour chaque paire. On obtient ainsi un système linéaire de la forme AU=B. La matrice colonne U contient les composantes des vitesses et les composantes des impulsions. Il faut à présent adopter une convention pour l'ordre de ces inconnues. Le programme de détection des collisions fournit une liste de paires de disques entrant en collisions. Le plus simple est donc de parcourir cette liste par indice croissant, ce qui fournit l'ordre pour les impulsions des paires. Le disques sont classés dans l'ordre dans lequel ils apparaissent dans les paires. Par exemple, si les paires sont (2,1), (2,3) et (1,3) (collision à 3 disques et 3 paires), la matrice U est :

U = {(p_{2, 1, x} p_{2, 1, y} p_{2, 3, x} p_{2, 3, y} p_{1, 3, x} p_{1, 3, y} v_{2, x} v_{2, y} v_{1, x} v_{1, y} v_{3, x} v_{3, y})}^{T}

L'implémentation devra bien sûr mémoriser l'ordre des vitesses des disques (sous forme de table associative) pour mettre à jour leur vitesse à la fin du calcul. Les matrices A et B sont construites lors du parcours des paires. La résolution du système linéaire est faite avec la méthode d'élimination de Gauss-Jordan (avec pivot sur les lignes et les colonnes).

3. Implémentation python

L'implémentation python présentée ci-dessous est prévue pour être intégrée dans un programme de simulation complet.

import numpy
import numpy.linalg
import math
from pylab import *
from matplotlib.patches import Ellipse,Arrow

On commence par définir une classe pour représenter un disque de rayon r et de masse m :

class Disque:
    def __init__(self,r,m):
        self.r = r
        self.m = m
    def set_position(self,x,y):
        self.x = x
        self.y = y
    def set_vitesse(self,vx,vy):
        self.vx = vx
        self.vy = vy
    def deplacer(self,t):
        self.x += self.vx*t
        self.y += self.vy*t
    def energie(self):
        return self.m*(self.vx*self.vx+self.vy*self.vy)
    def tracer(self,color):
        return Ellipse([self.x,self.y],2*self.r,2*self.r,color=color,fill=False)
    def tracer_vitesse(self,color):
        return Arrow(self.x,self.y,self.vx,self.vy,width=0.2,color=color)

Le programme de simulation stocke les disques sous forme d'une liste d'objets de la classe Disque. La fonction de détection des collisions fournit, pour chaque groupe de collision, la liste des paires de disques (sous forme d'indice sur la liste). On suppose que les disques sont placés dans leur position de collision. La fonction de résolution des collisions a pour arguments la liste des disques, la liste des paires du groupe de collision, et le coefficient de restitution :

Les principales variables utilisées dans cette fonction sont :

liste_disques : liste des objets disques de la simulation (liste globale)
liste_paires : liste des paires en collision, chacune sous la forme (i,j) paire d'indices.
liste_disques_collision : liste des indices des disques du groupe de collision, construite à partir de la liste des paires.
dict_ordre_disques : dictionnaire dont les clés sont les indices des disques. Pour chaque disque, ce dictionnaire fournit sa position dans la liste liste_disques_collision, qui est la position de la vitesse du disque dans la matrice des inconnues U.
dict_paires : dictionnaire dont les clés sont les indices des disques. Pour chaque disque, ce dictionnaire fournit la liste des paires auquel le disque appartient, sous la forme [p,signe], où p est l'indice de la paire et signe (-1 ou +1) la position dans la paire.

def collisions(liste_disques,liste_paires,e):
        Np = len(liste_paires)
        liste_disques_collision = []
        dict_paires = {}
        dict_ordre_disques = {}
        for p in range(Np):
            paire = liste_paires[p]
            i = paire[0]
            j = paire[1]
            if not(i in liste_disques_collision):
                dict_ordre_disques[i] = len(liste_disques_collision)
                liste_disques_collision.append(i)
                dict_paires[i] = [[p,1]]
            else:
                dict_paires[i].append([p,1])
            if not(j in liste_disques_collision):
                dict_ordre_disques[j] = len(liste_disques_collision)
                liste_disques_collision.append(j)
                dict_paires[j] = [[p,-1]]
            else:
                dict_paires[j].append([p,-1])
        Nd = len(liste_disques_collision)
        Ni = (Nd+Np)*2
        A = numpy.zeros((Ni,Ni))
        B = numpy.zeros((Ni,1))
        for i in range(Nd):
            indice = liste_disques_collision[i]
            d1 = liste_disques[indice]
            B[2*i,0] =  d1.m*d1.vx
            B[2*i+1,0] = d1.m*d1.vy
            A[2*i,2*Np+2*i] = d1.m
            A[2*i+1,2*Np+2*i+1] = d1.m
            for p in dict_paires[liste_disques_collision[i]]:
                ip = p[0]
                signe = p[1]
                A[2*i,2*ip] = signe
                A[2*i+1,2*ip+1] = signe
        for ip in range(Np):
            paire = liste_paires[ip]
            d1 = liste_disques[paire[0]]
            d2 = liste_disques[paire[1]]
            nx = d1.x-d2.x
            ny = d1.y-d2.y
            tx = ny
            ty = -nx
            i1 = dict_ordre_disques[paire[0]]
            i2 = dict_ordre_disques[paire[1]]
            eq1 = 2*Nd+2*ip
            eq2 = 2*Nd+2*ip+1
            B[eq1,0] = 0
            B[eq2,0] = -e*((d1.vx-d2.vx)*nx+(d1.vy-d2.vy)*ny)
            A[eq1,2*ip] = tx
            A[eq1,2*ip+1] = ty
            A[eq2,2*Np+2*i1] = nx
            A[eq2,2*Np+2*i1+1] = ny  
            A[eq2,2*Np+2*i2] = -nx
            A[eq2,2*Np+2*i2+1] = -ny
        s = numpy.linalg.solve(A,B)
        for i in range(Nd):
            indice = liste_disques_collision[i]
            d1 = liste_disques[indice]
            i1 = dict_ordre_disques[indice]
            d1.set_vitesse(s[2*Np+2*i1],s[2*Np+2*i1+1])

4. Exemple

On reprends la fonction de recherche de contact entre deux disques définies dans Collision sans frottement entre deux disques :

def contact(d1,d2):
    Dvx = d1.vx-d2.vx
    Dvy = d1.vy-d2.vy
    Dx = d1.x-d2.x
    Dy = d1.y-d2.y
    a = Dvx*Dvx+Dvy*Dvy;
    b = Dx*Dvx+Dy*Dvy;
    dmin = d1.r+d2.r 
    c = Dx*Dx+Dy*Dy-dmin*dmin;
    delta = b*b-a*c
    if delta <=0 or a==0:
        return False
    rd = math.sqrt(delta)
    t = min((-b+rd)/a,(-b-rd)/a)
    if t<0:
        return False 
    d1.deplacer(t)  
    d2.deplacer(t)
    return True

On considère comme exemple une collision entre 3 disques. Deux disques sont initialement immobile et l'un contre l'autre. Le troisième vient frapper un des deux disques avec un paramètre d'impact variable.

d1 = Disque(1,1)
d2 = Disque(1,1)
d3 = Disque(1,1)
b = 0.0 # parametre d'impact
d1.set_position(0,0)
d1.set_vitesse(0,0)
d2.set_position(2,0)
d2.set_vitesse(0,0)
d3.set_position(-10,b)
d3.set_vitesse(1,0)
liste_disques = [d1,d2,d3]
fig = figure(figsize=(12,6))
ax1 = fig.add_subplot(121,aspect='equal')
ax2 = fig.add_subplot(122,aspect='equal')
ax1.set_xlim(-4,4)
ax1.set_ylim(-4,4)
ax2.set_xlim(-4,4)
ax2.set_ylim(-4,4)
contact(d3,d1)
ax1.add_artist(d1.tracer('b'))
ax1.add_artist(d1.tracer_vitesse('b'))
ax1.add_artist(d2.tracer('r'))
ax1.add_artist(d2.tracer_vitesse('r'))
ax1.add_artist(d3.tracer('g'))
ax1.add_artist(d3.tracer_vitesse('g'))
liste_paires = [[0,2],[0,1]]
e=1.0
collisions(liste_disques,liste_paires,e)
ax2.add_artist(d1.tracer('b'))
ax2.add_artist(d1.tracer_vitesse('b'))
ax2.add_artist(d2.tracer('r'))
ax2.add_artist(d2.tracer_vitesse('r'))
ax2.add_artist(d3.tracer('g'))
ax2.add_artist(d3.tracer_vitesse('g'))

plotA.pdf