Modélisation chute avec frottements

[laura12345s]

Bonjour à tous, je cherche à résoudre un problème depuis quelques temps néanmoins, malgré de nombreux essais je reste bloqué à un endroit.

Nous avons une un propulseur à une altitude z=65km qui tombe verticalement vers le sol avec une vitesse initiale de 1020 m/s, je modélise la courbe de sa vitesse et de son altitude au cours du temps par morceaux :

- pour z< 40km : seule la force d'interaction gravitationnelle intervient
- pour z< 40 km : la force de frottement f de l'air intervient avec f= k v^2

Grâce à un programme python je modélise ça mais le problème est le suivant, quand faut-il rallumer les propulseurs avec une force constante F afin que:

- v (z=0m) = 0m/s
- v (z=14m) < 9m/s

La meilleure piste que j'ai est par une analyse énergétique en négligeant les forces de frottements par la suite:

∆Ec= 1/2m (vb^2 - va^2) = ∑Wab(F)= W(F) - W(P) = (Zb-Za)(F-P)

Après j'obtiens un système d'inconnues (Za, Va et F), avec :

- Zb=14
- Zb'=0
- Vb<9
- Vb'=0
- Za
- Va
- F
- P=mg

En le résolvant, j'ai :

- Va= 0.00889108 * sqrt(Za* (m * g - F))
- F=0.0357143 * (6.9494e6 - m * Vb ^ 2)

Je modifie donc mon programme afin de finaliser la modélisation lorsque la première condition sur la vitesse est remplie, or celle ci n'est jamais atteinte ...
Voici la partie du programme en question :



def E(n, a): ##a = vittesse à 14m (choisie entre 0 et 9)
d = 1 ## delta t
g = 9.81 ## force de gravité
v = 1020 ## vitesse initiale
V = [v] ## création de la liste vitesse
t = 0 ##t initial
T = [0] ## création de la liste temps
z = 65000 ## z initiale
m = 25.3e3 ##masse système
f = 0.0357143 * (6.9494e6 - m * a ** 2)
Z = [z]

for i in range(0, n + 1):

while fabs(v - l) > 2:
A = -7.2e-5
v = v + (A * (v ** 2) + g) * d ##approximation équation differentielle
z = z - fabs(v) * d
t = t + d
V.append(v)
T.append(t)
Z.append(z)
l = 0.00889108 * sqrt(z * (m * g - f))

else:
v = v + (g-(f / m)) * d
z = z - fabs(v) * d
t = t + d
V.append(v)
T.append(t)
Z.append(z)

plt.plot(T, Z, 'y')
plt.plot(T, V, 'p')
plt.grid()
plt.show()


Je vous remercie d'avance de votre aide, c'est tout l'intérêt de mon grand oral de résoudre ça et c'est dans 4 jours ...
-----

[jiherve]

bonsoir
c'est quoi "l" où est il initialisé ?
JR

[RomVi]

Bonjour

La gravité exerce une force m.g sur le mobile, le propulseur en exerce une autre de direction opposée, il suffit de faire la soustraction et de continuer ton algorithme comme tu l'a déjà fait.
Ensuite il suffit de faire une boucle en faisant varier la durée d'allumage pour faire tendre le résultat vers la valeur souhaitée.

[laura12345s]

Bonsoir, l est initialisé dans la partie supérieure du programme mais je n'ai pas tout envoyé pour ne pas surcharger, sinon l est calculée dans la boucle while

[A voir en vidéo sur Futura]

[laura12345s]

Merci pour votre réponse,

D'accord mais je me sers de la force calculée dans mon système, dans ma boucle else j'ai déjà mis que la force exercée et F avec "v = v + (g-(f / m)) * d " ? Si oui comment savoir à partir de quand je l'applique?
Pour ce qui est de la durée d'allumage, je veux seulement un atterrissage en douceur donc la valeur souhaitée, vous voulez dire une altitude de 0m ?

[RomVi]

On peut prendre une vitesse nulle à 0m comme objectif.
Tu as déjà effectué un algorithme qui permet de calculer l'évolution de la vitesse en prenant en compte l'accélération et les frottements, donc la vitesse au moment de toucher le sol.
Il faut ajouter une variable, que tu initialises à une valeur arbitraire de hauteur d'allumage "ha" (disons à 500m du sol) pour mettre en route le propulseur, ce qui permet de passer la valeur de g de 9.81 à une valeur négative (g - la poussée du propulseur). L'évolution de la vitesse sera modifiée, le mobile décélère.

Bien sur il serait étonnant d'atteindre l'objectif souhaité (vitesse nulle à 0m) du 1er coup : Si la vitesse nulle est atteinte pour une hauteur différente de 0 (la différence entre h=0 et h obtenue est appelée "résidu") on modifie la valeur de ha en conséquence et on relance le programme du début.
A chaque nouveau résultat on modifie la variable pour atteindre l'objectif souhaité. On peut faire évoluer cette variable par dichotomie, ou par une méthode plus évolué, en établissant une fonction de variation de ha en fonction du résidu (méthode de Newton par exemple, si tu l'a vue).

[Geo77b]

∆Ec= 1/2m (vb^2 - va^2) = ∑Wab(F)= W(F) - W(P) = (Zb-Za)(F-P)

Après j'obtiens un système d'inconnues (Za, Va et F), avec :

- Zb=14
- Zb'=0
- Vb<9
- Vb'=0
- Za
- Va
- F
- P=mg

En le résolvant, j'ai :

- Va= 0.00889108 * sqrt(Za* (m * g - F))
- F=0.0357143 * (6.9494e6 - m * Vb ^ 2)

Si ∆Ec= 1/2m (vb^2 - va^2) = (Zb-Za)(F-P)

Si Za = Zb' = 0 -> Va = 0
-> 1/2m (vb^2) = (Zb)(F-P)
-> F = (1/2m Vb² / Zb) + mg
-> F = (1 / (2*Zb)).(m.Vb² + m.g.2.Zb) = 0.0357143 * (6.9494e6 + m.Vb²)
avec Zb=14

Sauf erreur.

[laura12345s]

Super merci beaucoup, j'essaye ça demain

[laura12345s]

Il y a peut-être une erreur dans le signe mais Za≠0 car Za est l'altitude où l'on allume les propulseurs

[Geo77b]

Il y a peut-être une erreur dans le signe mais Za≠0 car Za est l'altitude où l'on allume les propulseurs

Les propulseurs sont-ils allumés entre Zb et Zb' ?

[laura12345s]

Le résultat que vous trouvez pour F est plus logique, dans mon système : (A->B) et (A->B') sont deux événements différents, le premier est l'évènement allumage propulseurs à atterrissage au sol et le deuxième jusqu'à z=14m donc oui les propulseurs sont allumés de Zb à Zb'

[Geo77b]

Le résultat que vous trouvez pour F est plus logique, dans mon système : (A->B) et (A->B') sont deux événements différents, le premier est l'évènement allumage propulseurs à atterrissage au sol et le deuxième jusqu'à z=14m donc oui les propulseurs sont allumés de Zb à Zb'

Mais ce ne sont pas les mêmes propulseurs, donc pas le même F.

[laura12345s]

Si, le but c'est que les deux conditions soient remplies : une vitesse inférieure à 9m/s à 14 m d'altitude et une vitesse nulle à 0m d'altitude, avec le moment A d'allumage des propulseur identique.

[Geo77b]

Donc le calcul de F est le même en remplaçant Za par Zb'=0.

[laura12345s]

Exact, merci

[laura12345s]

Voici le code, mais il doit y avoir une erreur en quelque part car les boucles tournent indéfiniment

def E(a):  ##a = vittesse à 14m
d = 1  ## delta t
g = 9.81  ## force de gravité
v = 1020  ## vitesse initiale
V = [v]  ## création de la liste vitesse
t = 0  ##t initial
T = [t]  ## création de la liste temps
z = 65000  ## z initiale
m = 25.3e3  ##masse système
f = 0.0357143 * (6.9494e6 + m * a ** 2)
    Z = [z]
    print(f)

    while z >= 40000:
            z = -0.5 * g * (t ** 2) - V[0] * t + Z[0]
            v = fabs(-g * t - V[0])
            t = t + d
            V.append(v)
            T.append(t)
            Z.append(z)
    else:
            b=z
            x=v
            y=t
            p=15
h=500
while fabs(p) > 10:
                h=h+1
z=b
                v=x
                t=y
                while z >= h:  ## hauteur allumage propulseurs
A = -7.2e-5
v = v + (A * (v ** 2) + g) * d  ## équation differentielle
z = z - fabs(v) * d
                        t = t + d
                        V.append(v)
                        T.append(t)
                        Z.append(z)
                else:
                    while z >= 0:
                            v = v +(g-(f/m) )* d  ## équation differentielle
z = z - fabs(v) * d
                            t = t + d
                            V.append(v)
                            T.append(t)
                            Z.append(z)
                    else:
                            p = v
            else:
                print("allumage à z =", h, "m")

[Geo77b]

les boucles tournent indéfiniment

Pour voir comment les boucles évoluent, vous pouvez afficher des données pendant le calcul.
ex :

[laura12345s]

En fait j'avais déjà affiché pour voir comment les boucles tournaient et elles ne tournaient pas comme il fallait, qu'est-ce qui est équivalent à "end while" en python?

[gts2]

La fin de while c'est le retour à l'indentation précédente.

[laura12345s]

Voilà ce que j'ai essayé mais ça n'aboutit pas

    
    while z >= 40000:
        z = -0.5 * g * (t ** 2) - V[0] * t + Z[0]
        v = fabs(-g * t - V[0])
        t = t + d
        V.append(v)
        T.append(t)
        Z.append(z)
        b = z
        x = v
        y = t
        p = 15
        h = 1500

    while  fabs(p) > 10:
            h = h + 10
            z = b
            v = x
            t = y
            while z >= h:  ## hauteur allumage propulseurs
                A = -7.2e-5
                v = v + (A * (v ** 2) + g) * d  ## équation differentielle
                z = z - fabs(v) * d
                t = t + d
                V.append(v)
                T.append(t)
                Z.append(z)
                print("frott: ""z = ", z, "   ; v = ", v)

            while z >= 0:
                v = v + (g - (f / m)) * d  ## équation differentielle
                z = z - fabs(v) * d
                t = t + d
                V.append(v)
                T.append(t)
                Z.append(z)
                p = v
                print("allum     ""z = ", z, "   ; v = ", v)
    print("allumage à z =", h, "m")

allum     z =  7574.267658692109    ; v =  1170.4157289857694
allum     z =  6406.74478628381    ; v =  1167.522872408299
allum     z =  5242.114770452981    ; v =  1164.6300158308286
allum     z =  4080.377611199623    ; v =  1161.7371592533582
allum     z =  2921.5333085237353    ; v =  1158.8443026758878
allum     z =  1765.581862425318    ; v =  1155.9514460984174
allum     z =  612.523272904371    ; v =  1153.058589520947
allum     z =  -537.6424600391056    ; v =  1150.1657329434765
allum     z =  37702.51785657748    ; v =  1242.7371434225297
allum     z =  36462.67356973242    ; v =  1239.8442868450593
allum     z =  35225.72213946483    ; v =  1236.9514302675889
allum     z =  33991.66356577471    ; v =  1234.0585736901185

[gts2]

Bonjour,

J'ai essayé, c'est assez problématique, peut-être un problème de données, d'où viennent-elles ?

[laura12345s]

C'est bien possible, pour les données sur un article de space X on trouve qu'après avoir ralenti les propulseurs pendant 30 s après la manœuvre, le propulseur entame sa chute libre à une vitesse Mach 3 donc 1020 m/s, le soucis c'est que je ne trouve pas l'altitude à laquelle commence cette chute en tout cas il y a une erreur car les propulseurs sont en fait déployés à 60km d'altitude donc il faut peut être considérer que la manoeuvre est finie lorsqu'il à une altitude de 50km, pour la masse je m'étais basé sur les données à l'atterrissage d'un sujet de bac, en revanche je viens de chercher à nouveau sur internet et la masse avant l'allumage des propulseurs serait de 72,5 tonnes et une masse à l'atterrissage de 23,5 tonnes.

[jiherve]

une idée en passant ,en partant du "bas":
h = 0 v= 0m/s ,h=14 v = 9m/s
on trouve F = (81/28 + g)*m.
Donc F ne dépend plus que de la masse cela devrait simplifier si celle ci est constante.
Dans la vraie vie cela sent le Runge Kutta, reste à écrire les équa diff.

JR

[laura12345s]

Merci, je vais essayer ça, oui c'est sûr vu que la masse varie aussi mais j'ai bien peur que 5min soient vraiment trop courtes pour expliquer tout ça ...

[laura12345s]

Merci à vous tous je vais rester avec cette modélisation:
Avec la solution finale :" allumage à z = 3052 m"

def E(h):  ##Altitude d'allumage propulseur
    g = 9.81  
    v = 1020  ## vitesse initiale
    V = [v]  
    t = 0  ##t initial
    T = [t]  
    z = 50000  ## altitude initiale
    m =72.5e3  ##masse système
    f = (81/28 + g)*m ##(m/28) * (28*g + a ** 2)
    Z = [z]
    print(f)

    while z >= 40000:
        d=1
        z = -0.5 * g * (t ** 2) - V[0] * t + Z[0]
        v = fabs(-g * t - V[0])
        t = t + d
        V.append(v)
        T.append(t)
        Z.append(z)

    while z >= h:
            d=0.001
            A =-6.45e-5
            v = v + (A * (v ** 2) + g) * d
            z = z - fabs(v) * d
            t = t + d
            V.append(v)
            T.append(t)
            Z.append(z)
    print("frott: ""z = ", z, "   ; v = ", v)

    while h>z >= 0:
            m=25.3e3
            v = v + (g - (f / m)) * d
            z = z - fabs(v) * d
            t = t + d
            V.append(v)
            T.append(t)
            Z.append(z)
    print("allum     ""z = ", z, "   ; v = ", v)



    print("allumage à z =", h, "m")

E(3052)
920957.1428571428
frott: z =  3051.6891977229657    ; v =  402.87415155863823
allum     z =  -0.0008474943186998716    ; v =  -0.9704744323826948
allumage à z = 3052 m