Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

Calcul scientifique Python Discussion :

Problèmes rebonds écoulement de sphère (3D)


Sujet :

Calcul scientifique Python

  1. #1
    Futur Membre du Club
    Homme Profil pro
    Étudiant
    Inscrit en
    décembre 2018
    Messages
    6
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Yvelines (Île de France)

    Informations professionnelles :
    Activité : Étudiant

    Informations forums :
    Inscription : décembre 2018
    Messages : 6
    Points : 6
    Points
    6
    Par défaut Problèmes rebonds écoulement de sphère (3D)
    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)

  2. #2
    Membre expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    804
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 804
    Points : 1 597
    Points
    1 597
    Par défaut
    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).

  3. #3
    Futur Membre du Club
    Homme Profil pro
    Étudiant
    Inscrit en
    décembre 2018
    Messages
    6
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Yvelines (Île de France)

    Informations professionnelles :
    Activité : Étudiant

    Informations forums :
    Inscription : décembre 2018
    Messages : 6
    Points : 6
    Points
    6
    Par défaut
    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

  4. #4
    Membre expérimenté

    Homme Profil pro
    Ingénieur calcul scientifique
    Inscrit en
    mars 2013
    Messages
    804
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Alpes Maritimes (Provence Alpes Côte d'Azur)

    Informations professionnelles :
    Activité : Ingénieur calcul scientifique

    Informations forums :
    Inscription : mars 2013
    Messages : 804
    Points : 1 597
    Points
    1 597
    Par défaut
    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à.

Discussions similaires

  1. Réponses: 1
    Dernier message: 17/03/2016, 20h00
  2. [plot3] cône et sphère : problème de tracer
    Par christophe_halgand dans le forum MATLAB
    Réponses: 5
    Dernier message: 25/06/2009, 21h46
  3. [c++] Hook "anti-rebond" et problème touches ALT
    Par jambono dans le forum Windows
    Réponses: 4
    Dernier message: 07/12/2006, 00h25
  4. Problème : ma sphère est toute trouée !
    Par julio26 dans le forum OpenGL
    Réponses: 2
    Dernier message: 14/05/2005, 14h47
  5. Problème avec le rendu de sphère
    Par Francky033 dans le forum DirectX
    Réponses: 10
    Dernier message: 29/08/2003, 23h00

Partager

Partager
  • Envoyer la discussion sur Viadeo
  • Envoyer la discussion sur Twitter
  • Envoyer la discussion sur Google
  • Envoyer la discussion sur Facebook
  • Envoyer la discussion sur Digg
  • Envoyer la discussion sur Delicious
  • Envoyer la discussion sur MySpace
  • Envoyer la discussion sur Yahoo