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) |
Partager