Simulation physique python

Simulation physique python - Python - Programmation

Marsh Posté le 03-12-2023 à 22:05:56    

Bonjour à tous,
 
Je suis actuellement en train de programmer une simulation d'une grue flottante sur l'eau. Mon code comprend des principes physiques. Je m'arrache les cheveux dessus depuis plusieurs semaines. Si quelqu'un pouvait m'aider à repérer ce qui ne tourne pas rond dans mon code, j'en serai vraiment très reconnaissant.
 
Cordialement,
 
Jacques
 
Voilà mon code:
 
#programme de simulation de l'oscillation d'une grue flottante
 
   
 
import math
 
import matplotlib.pyplot as plt
 
import numpy as np
 
   
 
### Constantes et données
 
   
 
g = 9.81         # gravitation [m/s**2]
 
rho = 1000       # masse volumique de l'eau [kg/m**3]
 
mci = 0.2        # masse du cylindre [kg]
 
Hci = 0.2        # hauteur du cylindre [m]
 
dci = 0.05        #diamètre du ccylindre [m]
 
d = 0.4          # distance finale entre le cylindre et l'axe situé au centre de la grue
 
   
 
### Paramètres du système
 
   
 
mp = 0.540          # masse de la plateforme [kg]
 
m1 = 2.16          # masse de la grue [kg]
 
m2 = 0.2         # masse du grappin [kg]
 
L = 0.6           # longueur de la plateforme [m]
 
l = 0.6           # largeur de la plateforme [m]
 
hp = 0.05          # hauteur de la plateforme [m]
 
hg = 0.6          # hauteur de la grue [m]
 
D = 0.4           # coëfficient d'amortissement
 
v0 = 0.5            # vitesse à laquelle le bras est étandu [m/s]
 
   
 
### Paramètres de la simulation
 
   
 
step = 0.001     # pas (dt) [s]
 
end = 10.0       # durée [s]
 
theta_0 = 0.0    # angle initial [rad]
 
omega_0 = 0.0    # vitesse angulaire initiale [rad/s]
 
   
 
t = np.arange(0, end, step)
 
theta = np.empty_like(t)           # angle [rad]    
 
omega = np.empty_like(t)           # vitesse angulaire [rad/s]
 
a = np.empty_like(t)               # accélération angulaire [rad/s**2]
 
Xg = np.empty_like(t)              # position en x du centre de gravité du systèmpe
 
Xc = np.empty_like(t)              # positions en x du centre de poussée du système
 
distance = np.empty_like(t)               # distance entre le cylindre et l'axe situé au centre de la grue
 
Ca = np.empty_like(t)
 
Cr = np.empty_like(t)
 
Cd = np.empty_like(t)
 
C_tot = np.empty_like(t)           # somme de tous les couples
 
e_flotteur = np.empty_like(t)      # énergie du flotteur
 
e_poussee = np.empty_like(t)       # éenrgie de poussée
 
e_charge = np.empty_like(t)        # énergie de la charge
 
e_cin = np.empty_like(t)           # énergie cinétique du système
 
e_tot = np.empty_like(t)           # énergie totale du système
 
   
 
   
 
def hauteur_volume_imerge():
 
     
 
    return (mp + m1 + m2 + mci)/(rho*L*l)
 
   
 
   
 
def centrep_centreg():
 
    # centre de gravité et de poussée avant d'appliqué la rotation
 
    # hypothèse que la masse de la grue est centrée, la position en x du centre de gravité vaut donc 0
 
    hc = hauteur_volume_imerge()
 
     
 
    Zc = -hc/2                   # coordonné en z du centre de poussée du volume imergé
 
     
 
    Zgplateforme = hp/2-hc       # coordonné en z du centre de gravité de la plateforme
 
    Zggrue = 2*hc                # coordonné en z du centre de gravité de la grue
 
    Zggrappin = hp-hc+dci/2      # coordonné en z du centre de gravité du cylindre, en sacahnt que le cylindre est couché sur la plateforme à sa position intiale
 
    Zg = (mp*Zgplateforme + (m1 + m2)*Zggrue + mci*Zggrappin)/(mp + m1 + m2 + mci)
 
     
 
    return Zg
 
     
 
     
 
def positions_critiques():
 
     
 
    hc = hauteur_volume_imerge()
 
    theta1 = math.atan(2*hc/L)        # angle de soulèvement
 
    theta2 = math.atan(2*(hp-hc)/L)   # angle de submersion
 
     
 
    return theta1, theta2
 
   
 
   
 
def simulation():
 
    """
 
    pre: -
 
    post: exécute une simulation jusqu'à t=end par pas de dt=step.
 
          Remplit les listes thèta, omega, a des angles d'oscillation, vitesses angulaires et accélérations angulaires.
 
    """
 
    # conditions initiales:
 
   
 
    theta[0] = theta_0
 
    omega[0] = omega_0
 
   
 
    # paramètres
 
    hc = hauteur_volume_imerge()
 
    Zg = centrep_centreg()        # centre de gravité initial (position en z)
 
   
 
     
 
    for i in range(len(t)-1):
 
         
 
        dt = step
 
         
 
        # Calcul du centre de poussée, de gravité et la distance entre les 2
 
        Xg[i] = (Zg-hc)* np.sin(theta[i])
 
        Xc[i] = (l**2/(12*hc))*math.sin(theta[i])
 
        dist = Xc[i] - Xg[i] # distance entre les positions en x du centre de gravité et du centre de poussée après rotation
 
         
 
        # variation de la distance par l'extension du bras de la grue (on suppose la vitesse de ce mouvement constante)
 
        # on considère que le cilindre est au bord de la plateforme au temps t = 0
 
         
 
        distance[i] = L/2 + v0*t[i]                  # distance entre le cylindre et l'axe des z (centre de gravité de la plateforme)
 
         
 
        # calcul de l'inertie
 
 
        I = ((mp + m1)*(hp**2 + L**2))/12 + (m2 + mci)*d**2   # on compte la masse du grappin dans avec la masse du cylindre
 
         
 
        # Calculs des couples
 
         
 
        Ca = (mci + m2)*g*distance[i]             # couple appliqué
 
        Cr = (mp + m1)*g*dist              # couple de redressement
 
        Cd = -D*omega[i]                    # frottement
 
        C_tot[i] = Ca + Cr + Cd
 
         
 
        # calcul de l'angle, la vitesse angulaire, l'accélération angulaire
 
        a[i] = C_tot[i]/I
 
        omega[i+1] = omega[i] + a[i] * dt
 
        theta[i+1] = theta[i] + omega[i] * dt
 
         
 
    # calcul des énergies:
 
     
 
    for i in range(len(t)-1):
 
         
 
        e_flotteur[i] = (m1 + m2 + mp)*g*(Xg[i] - Xg[0])
 
        e_poussee[i] = -(m1 + m2 + mp)*g*(Xc[i] - Xc[0])
 
        e_charge[i] = -Ca*theta[i]
 
        e_cin[i] = I*omega[i]**2 / 2
 
        e_tot[i] = e_flotteur[i] + e_poussee[i] + e_charge[i] + e_cin[i]
 
   
 
         
 
def graphiques():
 
     
 
    theta1, theta2 = positions_critiques()
 
     
 
    # graphiques theta, vitesse angulaire, accélération angulaire
 
     
 
    plt.figure(1)
 
    plt.subplot(3,1,1)
 
    plt.plot(t,theta, label="theta" )
 
    plt.axhline(y = theta1, linestyle = '--', color = 'r', label = 'submersion')
 
    plt.axhline(y = theta2, linestyle = '--', color = 'g', label = 'soulèvement')
 
    plt.legend()
 
    plt.subplot(3,1,2)
 
    plt.plot(t,omega, label="omega" )
 
    plt.legend()
 
    plt.subplot(3,1,3)
 
    plt.plot(t,a, label="a" )
 
    plt.legend()
 
    plt.show()
 
     
 
    # graphique des énergies
 
     
 
    plt.figure(2)
 
    plt.plot(t, e_flotteur, label="flotteur" )
 
    plt.plot(t, e_poussee, label="poussee" )
 
    plt.plot(t, e_charge , label="charge" )
 
    plt.plot(t, e_cin, label="cinétique" )
 
    plt.plot(t, e_tot, label="total" )
 
    plt.legend()
 
    plt.show()
 
     
 
   
 
   
 
   
 
# programme principal:
 
   
 
simulation()
 
graphiques()

Reply

Marsh Posté le 03-12-2023 à 22:05:56   

Reply

Marsh Posté le 04-12-2023 à 22:32:54    

Merci d'utiliser la balise [code] pour poster du code parce que là, ça facilite pas la lecture :/
 
Edit : et merci d'être un peu plus précis sur "ce qui ne tourne pas rond" dans le programme. On n'a pas de boule de cristal :o


Message édité par rufo le 04-12-2023 à 22:35:53

---------------
Astres, outil de help-desk GPL : http://sourceforge.net/projects/astres, ICARE, gestion de conf : http://sourceforge.net/projects/icare, Outil Planeta Calandreta : https://framalibre.org/content/planeta-calandreta
Reply

Marsh Posté le 05-12-2023 à 05:56:35    

J'imagine qu'en lançant le machin et en regardant les graphes on en aurait une idée.
Mais ouais va falloir être un peu plus précis pour encourager les gens à regarder ton truc de matheux, parce que là ça ne semble pas très excitant :o


---------------
Réalisation amplis classe D / T      Topic .Net - C# @ Prog
Reply

Marsh Posté le 05-12-2023 à 07:36:52    

Question bête : Python est sensible à la casse pour les noms de variables ?
Parce que appeler une variable D et une autre d, ça sent bon les embrouilles (même avec un langage sensible à la casse, ça peut embrouiller le lecteur) :/ Pareil pour les variables L et l. D'une manière générale, à part pour le variables de boucle, on évite de créer des variables ayant un nom à une seule lettre car c'est rarement explicite :o


---------------
Astres, outil de help-desk GPL : http://sourceforge.net/projects/astres, ICARE, gestion de conf : http://sourceforge.net/projects/icare, Outil Planeta Calandreta : https://framalibre.org/content/planeta-calandreta
Reply

Marsh Posté le 05-12-2023 à 15:50:36    

Plus de précision, car ton code fonctionne là !

Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed