v5.01.121 SDND121 – Système masse-ressort avec chocs sous excitation forcée#
Résumé:
Ce problème correspond à une analyse transitoire par recombinaison modale d’un système discret non linéaire à un degré de liberté. La non-linéarité consiste en un contact avec choc sur un plan rigide. La masse est soumise à une excitation forcée. Le problème est raide, c’est-à-dire que la raideur est très différente entre les phases de vol libre et les phases de contact. Ce problème permet de tester le bilan énergétique pour les différents schémas d’intégration temporelle, ainsi que la cinématique.
Solution de référence#
Méthode de calcul utilisée pour la solution de référence#
Lors d’une phase de vol libre, l’équation du mouvement s’écrit \(m\ddot{x}+{k}_{\mathrm{rappel}}x={F}_{\mathrm{ext}}(t)\) . Compte tenu de la forme sinusoïdale de \({F}_{\mathrm{ext}}(t)\) , cette équation admet une solution analytique. On peut calculer numériquement, avec une précision aussi grande que l’on veut, l’instant \({t}_{in}\) d’entrée dans le contact, vérifiant \(x({t}_{in})={\mathrm{jeu}}_{\mathrm{choc}}\) .
On est alors dans une phase de contact, dont l’équation est \(m\ddot{x}+{c}_{\mathrm{choc}}\dot{x}+({k}_{\mathrm{rappel}}+{k}_{\mathrm{choc}})x={F}_{\mathrm{ext}}(t)\) . Il existe là aussi une solution analytique. Avec les conditions aux limites issues de l’instant de contact, on peut calculer numériquement l’instant \({t}_{\text{out}}\) de sortie du contact.
En procédant ainsi itérativement, on obtient la solution complète du problème.
Remarque : on ne donne pas ici les formules analytiques. Les fichiers contenant les formules et permettant de calculer la solution complète sont joints avec le fichier de commande.
Résultats de référence#
On teste le bilan énergétique, l’adéquation entre les forces de contact et la cinématique, ainsi que la cinématique.
Pour le bilan énergétique, on calcule les énergies cinétique \({E}_{i}^{\mathrm{cin}}=\frac{1}{2}m\dot{x}{}_{i}^{2}\) , potentielle \({E}_{i}^{\mathit{pot}}=\frac{1}{2}{k}_{\mathit{rappel}}{x}_{i}^{2}\) , de choc \({E}_{i}^{\mathit{choc}}=\frac{1}{2}{k}_{\mathit{choc}}{p}_{i}^{2}\) (\(p\) est la pénétration ; cette expression n’est valable que si il n’y a pas d’amortissement de choc), injectée par la force extérieure \({E}_{i}^{\mathrm{inj}}=\sum_{j=1}^{i}{f}_{j}^{\mathrm{ext}}\dot{x}{}_{j+\delta \frac{1}{2}}\Delta t\) (avec \(\delta =1\) pour le schéma d’Euler, \(\delta =0\) pour le schéma des différences centrées). On obtient l’énergie totale \({E}_{i}^{\mathit{tot}}={E}_{i}^{\mathit{cin}}+{E}_{i}^{\mathit{pot}}+{E}_{i}^{\mathit{choc}}\) . On calcule enfin l’erreur globale sur le bilan énergétique par \({\mathit{erreur}}_{\mathit{globale}}^{\mathit{énergie}}=\sqrt{\frac{\sum_{i}{({E}_{i}^{\mathit{tot}}-{E}_{i}^{\mathit{inj}})}^{2}}{\sum_{i}{({E}_{i}^{\mathit{inj}})}^{2}}}\) , qui vaut idéalement 0.
Pour l’adéquation entre les forces de contact et la cinématique, on calcule la grandeur \({\mathit{erreur}}_{\mathit{globale}}^{\mathit{force}}=\sqrt{\frac{\sum_{i}{({F}_{i}^{\mathit{choc}}-{k}_{\mathit{choc}}{p}_{i})}^{2}}{\sum_{i}{({k}_{\mathit{choc}}{p}_{i})}^{2}}}\) , qui vaut idéalement 0.
Pour la cinématique, on compare les instants calculés d’entrée et sortie de contact, aux instants analytiques.
Incertitude sur la solution#
La solution est analytique par morceaux. Les instants d’entrée et sortie de contact sont déterminés numériquement à \({10}^{-9}s\) près.
Modélisation#
Caractéristiques de la modélisation#
Le patin est modélisé par une maille à un nœud (de type POI1), situé au repos en \(O=(0,0,0)\) . L’obstacle est de type “PLAN_Z”. Pour ne pas avoir de choc avec le plan symétrique, on décale suffisamment le centre des deux plans : ORIG_OBST=(-0.5,0.,0.,). Comme le jeu réel entre le patin et le plan de droite doit être de \(1\mathit{mm}\) , on utilise un jeu artificiel JEU=(0.5+jeu_choc).
Figure 3.1-a : Géométrie modélisée
L’intégration temporelle est réalisée sur une durée \(T=4s\) soit avec le schéma d’Euler, soit avec le schéma des différences centrées (“ADAPT_ORDRE2” qu’on force à un pas de temps constant grâce aux paramètres COEF_MULT_PAS=1.0, COEF_DIVI_PAS=1.0). Le problème est raide (le rapport des raideurs entre les phases de vol libre et de contact est de l’ordre de 5000), on est donc amené à utiliser un pas de temps très faible valant \(\Delta t=4\cdot {10}^{-6}s\) .
On réalise aussi un calcul avec “ADAPT_ORDRE2” à pas effectivement variable, dont le but n’est pas la précision, mais simplement la vérification de l’adéquation entre les forces de contact et la cinématique.
Caractéristiques du maillage#
Le maillage est constitué d’un unique nœud et d’une unique maille de type POI1.
Grandeurs testées et résultats#
Bilan énergétique et cinématique#
Dans les tableaux suivants, on donne les valeurs de l’erreur sur l’énergie et de quelques instants d’entrée / sortie de contact.
Schéma d’Euler :
Valeurs |
Référence |
Aster |
Tolérance |
\({\mathrm{erreur}}_{\mathrm{globale}}^{\mathrm{énergie}}\) |
0 J |
0.092 J |
0.1 J |
1ère entrée de contact (s) |
2.4867876 E-02 s |
2.4868 E-02 s |
1.2E-5 s |
1ère sortie de contact (s) |
2.5260518 E-02 s |
2.5264 E-02 s |
1.2E-5 s |
dernière entrée de contact (s) |
3.886525493 E+00 s |
3.886528 E+00 s |
1.2E-5 s |
dernière sortie de contact (s) |
3.886916559 E+00 s |
3.886916 E+00 s |
1.2E-5 s |
Tableau 4.1-1 : Résultats pour le schéma d’Euler
Schéma des différences centrées :
Valeurs |
Référence |
Aster |
Tolérance |
\({\mathrm{erreur}}_{\mathrm{globale}}^{\mathrm{énergie}}\) |
0 J |
0.063 J |
0.1 J |
1ère entrée de contact (s) |
2.4867876 E-02 s |
2.4868 E-02 s |
1.2E-5 s |
1ère sortie de contact (s) |
2.5260518 E-02 s |
2.5264 E-02 s |
1.2E-5 s |
dernière entrée de contact (s) |
3.886525493 E+00 s |
3.886528 E+00 s |
1.2E-5 s |
dernière sortie de contact (s) |
3.886916559 E+00 s |
3.886916 E+00 s |
1.2E-5 s |
Tableau 4.1-2 : Résultats pour le schéma des différences centrées
Adéquation entre force de contact et cinématique#
Valeur |
Référence |
Aster |
Tolérance |
\({\mathrm{erreur}}_{\mathrm{globale}}^{\mathrm{force}}\) |
0 N |
2.22 E-10 N |
1.E-8 N |
Tableau 4.2-1 : Résultats pour l’adéquation entre force de contact et cinématique
Synthèse des résultats#
On observe que les résultats sont très proches de la solution analytique, aussi bien pour le schéma d’Euler que pour le schéma des différences centrées. Les instants de changement de raideur obtenus par ces deux schémas sont identiques.
Les bilans énergétiques présentent une erreur inférieure à \(0,1J\) dans les deux cas. Cette valeur est à relativiser compte-tenu du pas de temps très faible qui a tendance à réduire les écarts, mais un pas de temps aussi faible est nécessaire pour traiter ce problème raide.
Enfin, il y a une adéquation parfaite (à la précision numérique près) entre les forces de contact et la cinématique, ce qui signifie que le traitement du contact par pénalisation est correctement effectué.