Bonjour,
Dans le cadre d'un projet d'informatique j'ai eu à coder l'écoulement de sphères en 2D puis en 3D sur python.
Pour cela je résout l'équation différentielle caractérisant le déplacement des sphères une à une jusqu'à ce qu'elles s’immobilisent ( je les considère comme immuables une fois arrêtées) grâce à la méthode d'Euler, je prend en compte la réaction avec le sol et les autres sphères, les frottements de l'air et le poids .
Mon problème se situe lors des contacts inter-sphères et avec le sol, je modélise le sol et les murs comme un ressort avec un fort coefficient de raideur (comme indiqué dans l'énoncé du sujet) et en faisant varier les coefficients de frottement de l'air et de raideur du sol j'obtient deux résultats non concluants : soit les sphères s'immobilise dès qu'elles touchent le sol sans aucun rebond, soit elles rebondissent indéfiniment en gagnant parfois de la hauteur à chaque saut.
J'ai essayé de réduire le pas de résolution de l'équation différentielle mais cela ne résout pas le problème.
Voici mon code :
Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
from math import *
import numpy as np
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from operator import itemgetter
 
def tri3D(L, x, y, z):
    '''Entrée : L : la liste des boules
                x,y,z: les coordonnées de la boule B qui chute
       Sortie : Liste triée par distance croissante avec B'''
    l2 = []
    for k in L:
        a, b, c, d, e = k[:5]
        l2 += [(a, b, c, d, e, sqrt((x - a) ** 2 + (y - b) ** 2 + (z - c) ** 2))]
    X = sorted(l2, key=itemgetter(5))
    return(X)
 
def pesanteur3D(boule): #boule correspond à un quintuplet (coord en x, coord en y,coord en z rayon de la boule, masse de la boule).
    '''Renvoie la composante selon y du poids appliqué à centre de masse de la boule en question'''
    return - boule[4] * 9.81
 
def moule3D(boule,k,largeur):
    '''Entrée : k : le coefficient de raideur du moule/récipient,
                largeur : la largeur du  cube dans lequel les boules tombent
       Sortie : Projection sur Ox, Oy et Oz de la résultante des réaction de support par rapport au récipient'''
    x, y, z, r = boule[:4]
    d1, d2, d3, d4, d5 = (x - r), (y - r), ((largeur - x) - r), (largeur - y - r), (z - r)    
    #d1 (respectivement d2, d3, d4, d5) les distance par rapport au mur de gauche (devant, droite, derrière, sol)
    F1x = abs(k * d1 *(d1 <= 0))
    F2y = abs(k * d2 *(d2 <= 0))
    F3x = -abs(k * d3 *(d3 <= 0))
    F4y = -abs(k * d4 *(d4 <= 0))
    F5z = abs(k * d5 *(d5 <= 0))
    return(F1x + F3x,F2y + F4y,F5z)                             #Composantes des forces
 
def angleplan(x1,y1,x2,y2):
    '''Entrée : (x1, y1) et (x2, y2) les coordonnées de deux points dans le plan (Ox, Oy)
       Sortie : Renvoie l'angle formé par le segment reliant les deux points et l'axe Ox'''
    if x1 == x2:
        if y1 > y2:
            return pi/2
        else :
            return -pi/2
    t = (y1 - y2)/(x1 - x2)
    if x1 > x2:
        return atan(t)
    else :
        if y1 > y2:
            return atan(t) + pi
        else :
            return atan(t) - pi
 
def bille3D(boule, L, k):
    '''Entrée :L la liste des boules déja placées triée par odre croissant de distance avec B,
               k le coefficient de raideur des sphères ( pris identique à celui du sol)
       Sortie : La projection sur 0x, 0y ,Oz de la résultante des réactions de support de la boule en question/ aux autres boules'''
    x, y, z, r = boule[:4]
    pro_x = pro_y = pro_z = 0
    compteur = 0
    while compteur < len(L) and L[compteur][5] <= (r + L[compteur][3]): #Calcul si contact
        i = L[compteur]
        (x1, y1, z1, r1), d, f, thetax = i[:4], i[5], 0, angleplan(x,y,i[0],i[1])
        d1 = sqrt((x - x1) ** 2+(y - y1) ** 2)
        if z > z1:
            thetaz = acos(d1 / d)
        else :
            thetaz = -acos(d1 / d)
        f=abs(k * (d - r - r1))                                         #Norme de Fb
        pro_x += f * cos(thetax) * cos(thetaz)                          #Projection sur Ox, Oy, Oz
        pro_y += f * sin(thetax) * cos(thetaz)
        pro_z += f * sin(thetaz)  
        compteur += 1
    return (pro_x, pro_y, pro_z)
 
def euler3D(boule, L, k, h, c, largeur):
    '''Entrée : boule : la boule que l'on souhaite placer;
                L : liste des boules déja placées; 
                k : coefficient de raideur; 
                h : pas de l'eq diff; 
                c  : coefficient de viscosité de l'air; 
                largeur : coté du cube
       Sortie : liste L correspondant à l'ancienne avec la nouvelle boule'''
    plt.close()
    x, y, z, r, m = boule
    compteur = 0
    vx = vy = vz = t = 0
    L2=tri3D(L, x, y, z)
    ax=affichage3D(L,largeur)
    while t==0 or abs(vy)>10**-9 or abs(vx)>10**-9 or abs(vz)>10**-9:    #Calcul si mouvement
        fgz=pesanteur3D(boule)                                           #Poids
        (fbx, fby, fbz)=bille3D(boule, L2, k)                            #Réactions avec les autres billes 
        (fmx, fmy, fmz)=moule3D(boule, k, largeur)                       #Réaction avec le support
        vxold, vyold, vzold = vx, vy, vz
        vx += h * ((fbx + fmx - c * vx) / m)                             #Calcul grâce à l'eq diff
        vy += h * ((fby + fmy - c * vy) / m)
        vz += h * ((fgz + fbz + fmz - c * vz) / m)
        x += h * vxold                                  
        y += h * vyold
        z += h * vzold
        t += h
        boule = (x, y, z, r, m)
        L2 = tri3D(L2, x, y, z)
        if compteur % 100 == 0:                        #Réaffichage tout les certaines fois pour pas ralentir trop le programme
           ax.scatter(x, y, z, s = 2900)
        plt.pause(0.001)
        compteur += 1
    #ax.scatter(x, y, z, s = 2900)
    #plt.pause(1)
    return (L + [boule])
 
def affichage3D(L, largeur):
    '''Entrée : L la liste des boules présentes dans le récipient, droite le quadruplet modélisant le vase 
       Affiche les boules dan le récipient'''
    plt.close()
    fig = plt.figure(figsize=(10, 10))
    ax = fig.gca(projection='3d')
    x, y = [0, 0, largeur, largeur, 0],[0, largeur, largeur, 0, 0]
    ax.plot(x, y, 0)
    ax.plot(x, y, largeur)
    for i in [[0,0], [largeur,largeur]]:
        for k in [[0,0], [largeur,largeur]]:
            ax.plot(i, k, [0, largeur])
    for b in L:
        ax.scatter(b[0], b[1], b[2], s = 3000, edgecolor='black')
    plt.show()
    return ax
 
def remplissage3D(l, m, n):
    '''Entrée : largeur : coté du cube,
                m : masse des boules,
                n : nombre de boules
       Modélise la chute des n boules de masse m dans un cube de largeur l'''
    L = []
    r = 0.3
    for k in range(n):
        print(k)
        x = random.randint(0,l*10000)/10000
        y = random.randint(0,l*10000)/10000
        z = l+1.5
        boule = (x, y, z, r ,m)
        L = euler3D(boule, L, 200, 0.01, 10, l)
    affichage3D(L, l)
    return(L)