Hello everyone,
I’m currently programming a simulation of a floating crane on water. My code includes physical principles. I’ve been tearing my hair out over it for several weeks. If anyone could help me spot what’s wrong with my code, I’d be really grateful.
All the best,
Jacques
#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()
. You do not even use the