Discussion :

[Wyxau]

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)

[lg_53]

Bonjour

Elles sont bien jolies vos fonctions, mais comment les appelez vous ? Si vous ne nous dites pas cela, comment voulez vous qu'on fasse tournez votre code ?

Tentez de débugger sur un cas simple : prenez juste une sphère et voyez comment elle réagit avec le sol (vous devriez observez des rebonds successifs, de plus en plus petit).

[Wyxau]

Bonjour,
Excusez moi de ne pas y avoir pensé...
Pour faire chuter des sphères il faut appeler la fonction remplissage3d avec par exemple comme variables : (3, 0.3, n) n étant le nombre de boule que l'on souhaite faire chuter

J'ai essayé avec un cas simple d'une boule et ce qu'il se passe c'est que hormis pour des boules de masse relativement faible il n'y a pas d'amortissement, des fois au contraire elle rebondit de plus en plus haut ce qui je pense est due au fait que le pas de résolution de l'eq diff dans ces cas la est trop faible ce qui fait que la boule a le temps de s'enfoncer trop profondément avant que le programme se rende compte qu'elle est dans le sol, comme le sol est modélisé par un ressort elle repars avec une trop grande impulsion.

Je pense que le problème de mon programme est qu'hormis les frottements de l'air il n'y a pas de perte d’énergie, alors que lors de l'impact au niveau du sol il devrait y avoir une perte d'energie. Le fait que pour des boules de faible masse cela fonctionne se justifie ainsi par le fait que les frottements de l'air ont une plus grande influence lors de la chute de légers objets.

Merci de votre réponse

[lg_53]

Vous avez votre code 2D ?

Car là sur le graphique en 3D, c'est difficile de voir...
D'autant plus que l'affichage de la boule reste lorsqu'elle descend, donc au lieu de voir une boule descendre, on a une trainée de boule. Et tout à la fin le graph disparait pour ne réafficher que la boule dans son état final. Faudrait modifier les données de votre graphique, plutot que de tout effacer et de retracer. Faites une recherche google pour voir comment on fait une animation avec matplotlib, vous aurez des petits exemples, comme celui là par exemple :
https://pythonmatplotlibtips.blogspo...atplotlib.html


De plus s'il y a plusieurs boules, vous ne voulez pas qu'elles bougent toutes en meme temps ?

Bon tout ceci ne résoud pas votre problème de départ, mais va tout de meme nous aider à y voir plus clair, car là ce n'est pas évident.

Oui ça peut être du au pas h. Vous devriez calculer une cfl à chaque itération dans euler3D. Elle dépend généralement de votre pas d'espace, et de la vitesse de la boule. Donc si la vitesse grandit, il est possible que vous violiez cette condition CFL. Là au moins si c'est le cas, vous pouvez afficher un warning par exemple, ou bien provoquer l'arret de la simulation (selon ce que vous voulez). Et si ca continue à tourner, et bien c'est que le problème viendra d'ailleurs et on aura écarté cette source là.