Instabilité d'un système masse-ressort à 2 noeuds avec RK45

Recherche :

Mot : Pseudo : Filtrer
Bas de page
Auteur	Sujet : Instabilité d'un système masse-ressort à 2 noeuds avec RK45

Feitan21

:hello:

Je cherche à modéliser des déformations d'un système avec un certain nombres de noeuds reliés entre eux par des ressorts.

J'ai un résultat un peu inattendu, lorsque mon noeud 0 est en position 0., 0., 0. et mon noeud 1 en position 1., 1., 1. tout fonctionne, et l'effet du ressort est bien cohérent. Le noeud 0 est toujours de masse infinie et le noeud 1 de masse 1 et j'applique une force similaire à la gravité en Y.

Le déplacement du noeud 1 en Y :

Par contre quand j'inverse la position des deux noeuds le ressort ne fait plus effet et mon noeud diverge vers l'infini.

Déplacement du noeud 1 en Y :

Je n'arrive pas à comprendre pourquoi j'obtiens ce résultat, je ne sais pas si c'est une instabilité numérique due au solver, ou si mon calcul de la force du ressort est erronée. J'ai pourtant essayé d'implementer une version standard avec des boucles imbriquées et une version vectorisée. Les deux donnent les mêmes résultats.

Voici mon code :

Code :

import numpy as np
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
import scipy.integrate
X = np.array([[0., 0., 0.], # 0
[1., 1., 1.]]) # 1
edges = np.array([[0, 1]])
W = np.array([[0, 1],
[1, 0]])
d = cdist(X, X, "euclidean" )
mass = 1 * np.ones(X.shape[0])
mass[0] = np.inf
mass[1] = 1
Force = np.zeros((X.shape[0], 3))
Force[1, 1] = 9.81
N = X.shape[0]
dt = 0.1
T = int(1e4)
k = 5
def compute_F_spring_vec(X, W, di):
d_inv = np.reciprocal(cdist(X, X, "euclidean" ))
d_inv[np.arange(np.shape(d_inv)[0]),np.arange(np.shape(d_inv)[0])] = 0
op1 = np.dot(W, X)
op2 = X * np.repeat(np.sum(W, axis=1), 3).reshape(X.shape)
op3 = np.dot(W * di * d_inv, X)
op4 = X * np.repeat(np.sum(W * di * d_inv, axis=1), 3).reshape(X.shape)
F_spring = - k * (-op1 + op2 + op3 - op4)
return F_spring
def compute_F_spring(X, W, d):
F_spring = []
for i in range(X.shape[0]):
F_i = np.zeros(3)
for j in range(X.shape[0]):
if W[i,j]:
F_i += -k * (np.linalg.norm(X[i, :] - X[j, :]) - d[i, j]) * (X[i, :] - X[j, :])/np.linalg.norm(X[i, :] - X[j, :])
F_spring.append(F_i)
return F_spring
def fun(t, X):
X = X.reshape((int(X.shape[0]/3), 3))
F_spring = compute_F_spring_vec(X, W, d)
F = F_spring + Force
a = F / mass[:, np.newaxis]
X = X + a*t
return X.ravel()
solver = scipy.integrate.RK45(fun, 0, X.ravel(), t_bound=100, rtol=0.001, vectorized=True)
Ylast = np.inf
X1 = []
while solver.status == 'running' and solver.t < 10:
error = solver.step()
X1.append([solver.y[4], solver.t, solver.y[3], solver.y[5]])
if (np.abs(X1[-1][0] - Ylast))/np.abs(Ylast) < 1e-6:
break
else:
Ylast = X1[-1][0]

Quelques précisions :

X contient les positions x,y,z de tous les nodes.
Edges la liaison entre les nodes, ce qui me permet de déduire une matrice de connectivité W (1 les noeuds sont connectés, 0 ils ne le sont pas).
Et d la distance euclidienne entre chaque paire de noeuds

Pouvez-vous m'aider ?

Merci

Message édité par Feitan21 le 09-07-2019 à 17:09:07

Publicité

MaybeEijOrNot

but someone at least

Sans lire le code, si tu places ta masse infinie en haut et que tu la soumets à une force de gravité, ça me semble logique que ça merde, puisse t-il y avoir un rapport ?

Message cité 1 fois

---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.

Feitan21

MaybeEijOrNot a écrit :

Sans lire le code, si tu places ta masse infinie en haut et que tu la soumets à une force de gravité, ça me semble logique que ça merde, puisse t-il y avoir un rapport ?

Hello, merci pour ta réponse.
Je n'applique la force que sur le node 1, quelque-soit sa position. D'ailleurs si j'applique la gravité sur le noeud de masse infinie j'obtiens comme attendu un non déplacement puisque l'accélération est nulle.

MaybeEijOrNot

but someone at least

Pour inverser les positions tu fais quoi ?

X = np.array([[1., 1., 1.], # 0
[0., 0., 0.]]) # 1

Et :

Code :

mass[0] = 1
mass[1] = np.inf

EDIT : et qu'est-ce que ça fait si tu fais :

Code :

X = np.array([[0., 0., 0.], # 0
[-1., -1., -1.]]) # 1

Message cité 1 fois
Message édité par MaybeEijOrNot le 09-07-2019 à 20:04:32

---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.

Feitan21

MaybeEijOrNot a écrit :

Pour inverser les positions tu fais quoi ?

X = np.array([[1., 1., 1.], # 0
[0., 0., 0.]]) # 1

Et :

Code :

mass[0] = 1
mass[1] = np.inf

?

EDIT : et qu'est-ce que ça fait si tu fais :

Code :

X = np.array([[0., 0., 0.], # 0
[-1., -1., -1.]]) # 1

Hello,

Pour inverser les positions je change juste les coordonnées x,y,z de mon noeud 0 et 1

X = np.array([[1., 1., 1.], # 0
[0., 0., 0.]]) # 1

Je ne touche pas à la masse.

D'ailleurs si je rentre tes paramètres de masse, j'ai bien mon point 1 fixe, et mon point 0 qui est cohérent.

Si je teste

Code :

X = np.array([[0., 0., 0.], # 0
[-1., -1., -1.]]) # 1

J'ai bien convergence, et cette valeur est bien cohérente puisque ma force est positive. Je tombe sur le même résultat qu'avec un noeud 1 placé en 1, 1, 1

Message édité par Feitan21 le 10-07-2019 à 14:43:41

MaybeEijOrNot

but someone at least

Et tu obtiens quoi comme matrice d'accélération dans les deux configurations ?

---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.

Feitan21

Je sors les images de la matrice d'accélération en Y pour le noeud 1 uniquement, le noeud 0 étant fixe dans les deux cas:

Pour la bonne partie, celle qui converge :

Et pour celle qui diverge :

En zoomant on remarque que le début est correct mais l'instabilité augmente au fur et à mesure avec des oscillations :

Et les données brutes, avec la matrice d'accélération en x,y,z pour le noeud 1 dans les deux cas : https://docs.google.com/spreadsheet [...] sp=sharing

MaybeEijOrNot

but someone at least

Si ce n'est pas la masse c'est peut-être la force.

Si tu inverses ta configuration, ta force ne devrait pas être :

Code :

Force[0, 1] = 9.81

:??:

Message cité 1 fois

---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.

Feitan21

MaybeEijOrNot a écrit :

Si ce n'est pas la masse c'est peut-être la force.

Si tu inverses ta configuration, ta force ne devrait pas être :

Code :

Force[0, 1] = 9.81

:??:

Non ma force est exacte, je ne change pas la position des noeuds dans ma liste, juste mes positions. C'est ce qui me permet d'appliquer toujours ma force en [1,1], soit la deuxième dimension du node 1. le node 0 étant toujours fixe. De plus j'ai bien mon noeud 0 qui est fixe, donc je ne pense pas que ce soit un problème d'application de la force.

Après pour une erreur dans le calcul de la force c'est possible, mais je trouve ça étonnant que ça fonctionne dans de nombreuses conditions et qu'il y ait quelques configurations ou ça ne fonctionne pas.

Message édité par Feitan21 le 10-07-2019 à 19:09:43

FORUM HardWare.fr

Programmation

Python

Instabilité d'un système masse-ressort à 2 noeuds avec RK45

Sujets relatifs
XSLT 1.0 - Grouper par liste de noeuds identiques	Système d'achievements/badges et statistiques
Système de log de champs de BD mis à jour	php - script de récupération infos système
Programme simple pour récupérer informations système d'une machine lin	Système de recommandation sur hébergeur mutualisé
Envoi d'emails en masse	Projet de système à développer
Systeme de session	système référentiel / masse vitesse et position relative.
Plus de sujets relatifs à : Instabilité d'un système masse-ressort à 2 noeuds avec RK45

Page générée en 0.062 secondes