Python physics simulation

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,


#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.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.plot(t,omega, label="omega")
    plt.plot(t,a, label="a")
    # graphique des énergies
    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")
# programme principal:

Hi and welcome.

I don’t know that anyone is going to be willing to review several hundred lines of code, in order to help you to:

so why not “clue us in” and tell us what’s wrong, rather than have us try to find that out.

By “what’s wrong”, I mean: what is the code doing that it should not be doing, or what is it not doing, that it should be doing; throw us bone, least ways :slight_smile:


The results for angle theta in radian, angular velocity in rad/sec and angular acceleration in rad/sec^2 are not consistent with reality. The values are far too large. Here is the result.

Thanks in advance for your help,


1 Like

The only thing that I can see (and this is purely from a code analysis view point, as I know nothing about the math/physics here) is that at code line 107, you have:
Zc = -hc / 2 # coordonné en z du centre de poussée du volume imergé

… but Zc is never used. Maybe this is nothing; maybe it’s everything.

Edit to add: maybe that line is not at 107 in your copy. The system I use does an “auto reformat” before the analysis , and as such, my code line numbers will not match yours. I’m sure you can find that code line, in any case.

Start using smaller known values with no plotting and do asserts to see if the algorithm you defined works or not, when that’s working then move on to plots and so on. It’s always good to first get it to work in a small scale.

1 Like

Hi !

As a mechanical engineer myself, I can tell your calculation for the center of mass of the entire crane is wrong :upside_down_face:. You do not even use the hg variable, which is one of the most important.

Now figuring out the correct formula unfortunately has nothing to do with Python, so you should rather ask for help on physics/engineering forums on this aspect.

Good luck !