v2.08.012 FORMA12 – Travaux pratiques de la formation « introduction à l’analyse dynamique linéaire & non-linéaire »#
Résumé :
Ce cas est est une introduction à l’usage de Code_Aster pour l’analyse transitoire dynamique (linéaire et non-linéaire). Il traite d’une plaque épaisse d’acier, sujette à une excitation sinusoïdale sur une table vibrante.
La plaque est modélisée en éléments 3D. On compare l’analyse transitoire sur base modale avec une analyse transitoire directe. Dans la dernière partie on prend en compte la plasticité de l’acier par une une loi non-linéaire.
Les différentes modélisations abordées sont:
modélisation A : analyse modale ;
modélisation B : transitoire direct ;
modélisation C : transitoire sur base modale avec correction statique a priori ;
modélisation D : transitoire sur base modale avec correction statique a posteriori ;
modélisation E : transitoire non-linéaire;
modélisation F : transitoire linéaire avec non-linéarité de contact;
modélisation G : transitoire non-linéaire par sous-structuration dynamique.
Le premier chapitre présente la géométrie et les données communes à toutes les modélisations.
Modélisation A – Analyse modale#
Objectifs#
Pour avoir une idée des propriétés dynamiques du système, on commence par une analyse modale.
La géométrie est dessinée puis maillée sous Salome_Meca. Les premiers modes de la plaque sont ensuite calculés avec Code_Aster .
Géometrie dans SALOME#
On modélise le maillage dans Salome.
Dans le module Geometry, on utilise le menu New Entity > Primitive > Box pour construire une boîte avec les dimensions de la plaque.
Fig. 502 Création de la géométrie via une box#
Pour spécifier dans Code_Aster les conditions aux limites et identifier les nœuds utiles aux post-traitements, on a besoin de nommer quelques éléments de la géométrie et quelques nœuds.
C’est fait grâce à New Entity > Group > Create.
La plaque est fixée sur un des côtés. On se sert d’un point de la face supérieure, dans un angle opposé au bord encastré, comme nœud de post-traitement.
Pour pouvoir préciser la finesse du maillage dans l’épaisseur de la plaque, il peut être utile de créer un groupe correspondant à l’arrête verticale d’un des coins de la plaque.
Fig. 503 Géométrie de la plaque#
Maillage dans SALOME#
L’étape suivante est la création du maillage par le module Mesh. De nombreux mailleurs (algorithmes) sont disponibles dans le logiciel.
Pour obtenir un maillage réglé et prédictible, un choix possible est:
Hexahedron(i,j,k)/Quadrangle(Mapping)/Wire discretisation, sucessivement correspondant aux maillages 3D, 2D et 1D.
Fig. 504 Algotithme pour générer le maillage 3D#
Fig. 505 Algotithme pour générer le maillage 2D#
Fig. 506 Algotithme pour générer le maillage 1D#
Pour éviter un maillage trop fin dans l’épaisseur de la plaque, on conseille de créer un sous-maillage sur l’une des arrêtes de coin.
N.B.: il ne faut pas oublier de «propager» la discrétisation de l’arrête de coin choisie sur les autres en le précisant dans Add. Hypothesis.
Fig. 507 Algotithme pour générer le maillage 1D avec hypothèse de discrétisation de l’arrête#
Comme compromis entre temps de calcul et précision, on pourra créer 40 éléments dans la largeur de la plaque et 2 dans l’épaisseur.
On obtient alors un maillage grossier mais ce défaut peut être compenser par des éléments quadratiques.
Pour passer d’éléments linéaires à des éléments quadratiques on utilise l’option suivante dans le menu meshing: Modification > Convert to/from quadratic .
A la fin de la procédure, le maillage comprend 4725 nœuds et 800 éléments 3D.
Fig. 508 Information de la maillage créé#
Fig. 509 Maillage créé (vue 3D)#
Analyse modale#
On demande des premiers modes de la structure jusqu’à une fréquence de \(50\mathit{Hz}\) . C’est un choix prudent puisque la fréquence d’excitation ne s’élève qu’à \(15\mathit{Hz}\) .
Pour ce calcul très simple, la meilleure solution est le recours à un assistant. Pour ce faire, cliquer «patte de droite» sur l’entrée Current case d’AsterStudy puis sur Add Stage with Assistant. Laissez-vous guider et lancer votre calcul dans l’onglet History View.
A noter qu’il peut être nécessaire d’augmenter la mémoire allouée au calcul dans la case Memoryde la fenêtre Run parameters (voir Fig. 510).
Post-traitement#
Pour les post-traitements simples, il est très commode de cliquer «patte de droite» dans l’entrée qui correspond au fichier MEDde sortie dans la fenêtre Data Files Summary(voir Fig. 511 ).
Vous pouvez parcourir les différents modes en cliquant sur le bouton:
Fig. 510 Run via AsterStudy#
Fig. 511 Paramètre de lancement de calcul#
Voici, à titre d’exemple quelques déformées modales:
Illustration 1: Premier mode (6,3 Hz) 
|
Illustration 2: Deuxième mode (15.4 Hz) 
|
Illustration 3: Troisième mode (38.6 Hz) 
|
Illustration 4: Quatrième Mode (49.0 Hz) 
|
Pour vérifier la qualité du maillage, on peut, bien entendu, faire une étude de convergence en fonction de la finesse du maillage. Mais sur une géométrie aussi simple, il est possible de trouver une bonne approximation par une formule analytique (par exemple dans Formulas for Stress, Strain, and Structural Matrices , Walter D. Pikey – éd. John Wiley & Sons, Inc. 1994).
Dans le cas étudié, on trouve:
\(\mathit{f1}=6.23\mathit{Hz}\)
\(\mathit{f3}=38.16\mathit{Hz}\)
Les résultats donnés par Code_Aster sont donc plutôt précis, malgré la grossièreté du maillage.
Avec un maillage linéaire, les résultats n’auraient pas été aussi bons.
Modélisation B – Analyse transitoire directe#
Problème et esquisse du calcul#
On souhaite maintenant calculer le mouvement de la plaque sur sa table vibrante. On commence par le calcul le plus naturel, le calcul direct sur base physique. On réalise deux calculs : sans et avec amortissement type Rayleigh.
Une question essentielle réside dans la modélisation du chargement sur la plaque. On choisit ici d’utiliser une force de pesanteur d’intensité et de direction correspondant aux réglages de la table vibrante.
Les principales étapes du fichier de commande sont:
début : DEBUT
lecture du maillage: LIRE_MAILLAGE
construction du modèle: AFFE_MODELE
définition des matériaux : DEFI_MATERIAU
affection des matériaux aux différents groupes du maillage: AFFE_MATERIAU
conditions aux limites : AFFE_CHAR_MECA
définition du chargement : AFFE_CHAR_MECA(PESANTEUR=_F(GRAVITE=300.,DIRECTION=(-1,0,1))
construction des matrices de masse et de raideur et du vecteur de second membre : ASSEMBLAGE
définition de la modulation en temps (sinus de période :math:`15times 2pi ` secondes) : FORMULE
analyse transitoire : DYNA_VIBRA(TYPE_CALCUL=”TRAN”, BASE_CALCUL=”PHYS”)
On étudie la réponse de la structure jusqu’à \(0.5s\) avec un pas de temps de \(0.002s\) (à renseigner sous le mot-clé INCREMENT)
récupération du déplacement en un point : RECU_FONCTION
impression de la courbe : IMPR_FONCTION
Résultats#
Voyons la courbe des résultats:
Fig. 512 Déplacement vertical (DZ) en P#
On observe une réponse sinusoïdale à la fréquence d’excitation. Cette réponse sinusoïdale est néanmoins fortement modulée par la réponse du premier mode propre de la plaque.
Prise en compte de l’amortissement#
On se propose maintenant de prendre en compte l’amortissement dans le cadre du calcul réalisé précédemment. Pour ce faire, on considère le modèle d’amortissement proportionnel de Rayleigh, avec les paramètres suivants :
\(\alpha = 0.00011575 s\) : valeur à renseigner dans AMOR_ALPHA dans DEFI_MATERIAU;
\(\beta = 1.142397 {s}^{-1}\) : valeur à renseigner dans*AMOR_BETA* dans DEFI_MATERIAU .
Ces paramètres correspondent à un amortissement réduit cible de 2%, sur la plage de fréquences [5; 50] Hz.
Une fois ces paramètres renseignés dans l’opérateur DEFI_MATERIAU, il est nécessaire d’assembler la matrice d’amortissement via la macro-commande ASSEMBLAGE (option AMOR_MECA). Le calcul transitoire avec amortissement peut ensuite être réalisé via l’opérateur DYNA_VIBRA, en renseignant la matrice d’amortissement assemblée via l’option MATR_AMOR.
On compare alors les déplacements obtenus avec ceux obtenus sans amortissement. L’évolution temporelle de la composante DZ des déplacements avec et sans amortissement de Rayleigh est présentée ci-dessous :
Fig. 513 Comparaison de déplacement vertical entre deux calculs: sans et avec amortissement Rayleigh#
Modélisation C – analyse transitoire sur base modale#
Problème et indices#
L’analyse transitoire peut aussi se faire plus efficacement en projetant le problème sur une base modale. On se propose de comparer cette méthode avec le calcul transitoire direct.
Dans un premier temps, on propose de résoudre le transitoire sur la famille de modes déterminée dans la modélisation A. Pour ce faire, on construit une base modale avec les modes précédents, on projette les matrices et les vecteurs dans cette nouvelle base et on résout le transitoire. Ces étapes sont précisées ci-après:
calcul des modes comme dans la modélisation A
construction d’un vecteur contenant le second membre par CALC_VECT_ELEM et ASSE_VECTEUR
définition d’une base modale par DEFI_BASE_MODALE
création d’une numérotation généralisée par NUME_DDL_GENE
Si la base est orthonormée, on peut l’indiquer par le mot-clé STOCKAGE=’DIAG’
projection des matrices de masse, de raideur et du vecteur de chargement : PROJ_VECT_BASE et PROJ_MATR_BASE
comme le vecteur est un chargement, on indique TYPE_VECT=”FORC”
calcul transitoire sur base modale : DYNA_VIBRA(TYPE_CALCUL=”TRANS”, BASE_CALCUL=”GENE”)
on calcule la réponse jusqu’à \(0.5s\) avec un pas de temps de \(0.002s\)
récupération du déplacement du point P par RECU_FONCTION(RESU_GENE)
Il est opportun de s’interroger sur quelle base s’exprime alors ce déplacement
impression de la courbe déplacement P/temps : IMPR_FONCTION
On souligne que la commande DYNA_LINE assure toutes ces étapes de manière très synthétique. Attention: par défaut, cette commande détermine l’enrichissement statique. Pour ne pas le faire, il faut préciser ENRI_STAT=”NON”.
Résultats des calculs sur base modale sans correction statique#
Pour le déplacement vertical on obtient des résultats très proches du calcul transitoire direct.
Fig. 514 Comparaison des déplacements en DZ en point P entre base physique et base modale jusqu’à 50 Hz#
Mais si on examine les déplacements dans la direction horizontale, ils sont assez différents !!!
Fig. 515 Comparaison des déplacements en DX en point P entre base physique et base modale jusqu’à 50 Hz#
Quoique la base modale ait été calculée au-delà de \(50\mathit{Hz}\) et qu’on ait respecté la règle habituelle recommandant de couvrir deux fois la fréquence maximale de l’excitation, dans la direction horizontale, la méthode par projection modale s’éloigne du résultat direct.
Que s’est-il passé?
La plaque est en fait très raide dans son plan (direction horizontale). Les modes propres ne représentent donc qu’imparfaitement ce mouvement. Pour le restituer correctement, il faudrait prendre bien plus de modes propres.
Calculs sur base modale avec correction statique à priori#
Une méthode se montre toutefois bien plus performante. Il s’agit de la méthode de la correction statique, qui permet de prendre en compte la contribution «statique» de la structure au mouvement. Elle peut être réalisée de 2 manières:
a priori: on intègre un nouveau mode dans la base modale sur laquelle résoudre le transitoire;
a posteriori: on conserve la base modale initiale mais on ajoute une correction à la solution obtenue.
Nous allons examiner ces 2 possibilités.
Tout d’abord la correction statique a priori (modélisation C). Les étapes du calcul transitoire sur base modale avec correction statique à priori sont les suivantes:
On calcule des modes de la plaque jusqu’à 50 Hz via CALC_MODES.
On calcule du pseudo-mode via MODE_STATIQUE/SPEUDO_MODE.
On définit une nouvelle base modale par ajout du mode statique aux modes « dynamiques » par DEFI_BASE_MODALE / RITZ.
Comme ces différents modes ne sont a priori pas orthogonaux entre eux, on spécifie ORTHO=”OUI” pour réaliser une orthonormalisation. On doit également indiquer quelle matrice sera utilisée pour cette opération (orthogonalité dans la norme de la matrice).
Création d’une numérotation généralisée par NUME_DDL_GENE
Si la base est orthonormée, on peut l’indiquer par le mot-clé STOCKAGE=’DIAG’
projection des matrices de masse, de raideur et du vecteur de chargement : PROJ_VECT_BASE et PROJ_MATR_BASE
Comme le vecteur est un chargement, on indique TYPE_VECT=”FORC”
Calcul transitoire sur base modale : DYNA_VIBRA(TYPE_CALCUL=”TRANS”, BASE_CALCUL=”GENE”)
On calcule la réponse jusqu’à \(0.5s\) avec un pas de temps de \(0.002s\)
On a recours au schéma temporel de Devogelaere-Fu par SCHEMA_TEMPS / SCHEMA=”DEVOGE”
construction de la réponse globale de la structure sur base physique par REST_GENE_PHYS
récupération du déplacement du point P RECU_FONCTION(RESULTAT)
impression de la courbe déplacement P/temps : IMPR_FONCTION
On souligne que la commande DYNA_LINE assure toutes ces étapes de manière très synthétique.
La correction statique a posteriori va être traitée dans la modélisation D.
Fig. 516 Comparaison des déplacements en DX en point P entre base physique et base modale jusqu’à 50 Hz avec correction à priori#
Prise en compte de l’amortissement dans des calculs sur base modale avec correction statique à priori#
On se propose de prendre en compte l’amortissement dans le cadre des calculs réalisés précédemment. Dans cette optique, l’amortissement est ici modélisé via un coefficient d’amortissement réduit, égal à :math: xi =1.86%.
Pour chaque cas (i.e., base modale, correction statique a priori, DYNA_LINE), on procède de manière purement analogue aux calculs sans amortissement, en spécifiant le coefficient d’amortissement réduit précité via les mots-clés AMOR_MODAL/AMOR_REDUIT dans le cas l’opérateur DYNA_VIBRA, et via le mot-clé AMORTISSEMENT dans le cas de l’opérateur DYNA_LINE.
L’évolution temporelle de la composante DZ des déplacements obtenus avec et sans amortissement, dans le cas d’un calcul transitoire sur base modale avec correction statique a priori est présentée ci-dessous :
Fig. 517 Comparaison des déplacements en DX en point P entre deux calculs sur base modale jusqu’à 50 Hz avec correction à priori, avec et sans amortissement modal à 1.86%#
Modèle D- correction statique a posteriori#
Procédure de calcul#
Les étapes sont assez proches de la méthode sans correction statique. Mais on ajoute quelques étapes pour prendre en compte les déformations dues aux effets statiques du chargement.
Attention! Il ne faut pas oublier d’inclure la correction lors du retour vers la base physique.
calcul de la base modale jusqu’à 50 Hz ;
définition du chargement de pesanteur :
AFFE_CHAR_MECA(PESANTEUR=_F(GRAVITE=300.,DIRECTION=(-1,0,1))
calcul de la correction statique : MACRO_ELAST_MULT
projection des matrices de masse, de raideur et du vecteur de chargement : PROJ_BASE
donnée de la fonction multiplicative en temps (sinus de période \(15\times 2\pi\) s) : FORMULE
définition d’une liste d’instants (On peut prendre une liste d’instants sur \(0.5s\) avec un pas de temps de \(0.002s\) ) : DEFI_LIST_REEL
On a aussi besoin des dérivées premières et secondes de la fonction multiplicative en temps: formule (si on connaît leur forme analytique ou CALC_FONCTION pour une approximation numérique)
calcul transitoire sur base modale : DYNA_VIBRA(TYPE_CALCUL=”TRANS”, BASE_CALCUL=”GENE”) – donner derrière MODE_CORR la déformée statique et spécifier derrière D_FONC_DT et D_FONC_DT2 les dérivées premières et secondes
récupération du déplacement du point P RECU_FONCTION(RESU_GENE).Ne pas oublier CORR_STAT=”OUI”
impression de la courbe déplacement P/temps :IMPR_FONCTION
Résultats#
Grâce à la correction statique à postériori, les résultats sont bien meilleurs.
Fig. 518 Comparaison des déplacements DX en point P entre base physique et base modale avec correction statique à posteriori#
Modélisation E – dynamique non-linéaire#
L’objectif principal de ce TP est de mettre en place un calcul non-linéaire transitoire avec prise en compte de la plasticité. Nous allons utiliser l’opérateur DYNA_NON_LINE.
Procédure de calcul#
Pour calculer le comportement non-linéaire dynamique, on utilise l’opérateur d’usage général DYNA_NON_LINE. Il est proche de STAT_NON_LINE et son usage en proche. La principale différence tient au schéma d’intégration en temps et à l’influence de l’inertie et de l’amortissement.
Données générales#
Par rapport aux calculs précédents dans DYNA_VIBRA et DYNA_LINE, la différence fondamentale est que les matrices (de rigidité, d’amortissement et de masse) et les seconds membres (chargements) ne sont pas à calculer avant l’opérateur. En effet, la commande DYNA_NON_LINE est dite «monolithique» car la matrice de rigidité varie pendant le transitoire puisqu’on résout un problème non-linéaire (qui dépend du déplacement). Il faut donc lui donner les informations comme vous le faites dans la macro-commande ASSEMBLAGE pour qu’il évalue lui-même ses opérateurs:
CHAM_MATER
MODELE
éventuellement CARA_ELEM si vous avez des éléments de structure ou des discrets (ce qui n’est pas le cas ici)
Les chargements sont les mêmes que dans les exercices précédents. Ils sont donnés par le mot-clef facteur EXCIT comme dans DYNA_VIBRA. Ne pas oublier d’appliquer le sinus dans FONC_MULT pour la gravité.
Choix du comportement non-linéaire#
Le choix du comportement non-linéaire est fait dans le mot-clef facteur COMPORTEMENT de DYNA_NON_LINE. Si vous ne renseignez rien, on va considérer que le calcul est élastique en petites déformations.
Bien entendu, si vous avez un comportement non-linéaire, il faut définir les caractéristiques matériaux idoines dans DEFI_MATERIAU.
Dans cet exercice, on vous propose de la plasticité de von Mises à écrouissage cinématique linéaire:
dans DEFI_MATERIAU, il faut choisir ECRO_CINE. On propose D_SIGM_EPSI = 2GPa (le module de Young divisé par 100) et une limite élastique SY = 200MPa. Il est hautement recommandé de toujours renseigner le mot-clef ELAS, même quand on fait de la plasticité.
dans DYNA_NON_LINE/COMPORTEMENT, on choisit RELATION = “VMIS_CINE_LINE” * et *DEFORMATION =’PETIT’.
Choix du schéma en temps et discrétisation temporelle#
Le choix du schéma en temps est fondamental en dynamique. Par ailleurs, un problème non-linéaire ajoute une contrainte supplémentaire sur la discrétisation temporelle pour avoir une précision suffisante pour l’intégration de la non-linéarité au niveau du point de Gauss (plasticité). Dans code_aster, la majorité des schémas de discrétisation des équations du comportement sont implicites et précis, mais cette précision peut dépendre de cette discrétisation, en particulier parce qu’elle fait quelques hypothèses sur la nature du problème (voir formation avancée de code_aster). Il s’agit donc d’un paramètre numérique supplémentaire sur lequel il faut être très attentif.
Il convient de choisir une discrétisation temporelle qui a trois caractéristiques:
compatible avec le chargement appliqué. Il faut veiller à avoir suffisamment de points de discrétisation temporelle pour bien capturer ce chargement. L’application de la règle de Shannon est une bonne base: au moins une fréquence de discrétisation deux fois supérieure à celle du chargement le plus pénalisant
la discrétisation temporelle doit également être choisie en fonction des caractéristiques du schéma en temps: condition de stabilité (pour les schémas qui en ont une) et précision
La discrétisation temporelle doit être choisie pour intégrer précisément la partie non-linéaire
On propose de prendre un pas de temps de 0,002 s. Le calcul non-linéaire étant potentiellement très long, on recommande de s’arrêter au temps tfin = 0,1s:
LISTR=DEFI_LIST_REEL(DEBUT=0.,INTERVALLE=_F(JUSQU_A=0.1,PAS=0.002,),)
La contrainte supplémentaire du comportement non-linéaire impose de permettre éventuellement à DYNA_NON_LINE de découper le pas de temps lorsqu’il échoue dans l’intégration. Pour cela, on utilise l’opréteur DEFI_LIST_INST:
LIS=DEFI_LIST_INST(METHODE=”MANUEL”,DEFI_LIST=_F(LIST_INST=LISTR,),)
Cette commande dit que le cadencement du calcul se fera sur la base de la liste d’instants LISTR (et, par exemple, on sauvegardera toutes les données du calcul à ces instants),
Mais qu’en cas de problème (échec d’intégration de la non linéaire par exemple), on l’autorise à découper le pas de temps en 5, via * ECHEC=_F(SUBD_NIVEAU=5,SUBD_PAS_MINI=1e-5,SUBD_METHODE=”MANUEL”,),*. La découpe est récursive: il découpe en 5, s’il échoue, il redécoupe en 5, etc. Il s’arrêtera définitivement (échec) dès qu’il atteindra un pas minimum donné par le paramètre SUBD_PAS_MINI=1E-5.
Le résultat de cette commande est à renseigner directement dans le mot clef INCREMENT/LIST_INST de DYNA_NON_LINE. Attention à bien choisir LIS et non LISTR, sinon la découpe du pas de temps en cas d’échec ne sera pas activée.
Dans le mot-clef SCHEMA_TEMPS, vous pouvez sélectionner votre schéma en temps. Par défaut, il s’agit d’un schéma de Newmark en déplacement. Ce schéma est inconditionnellement stable et conserve l’énergie. D’autres schémas sont disponibles, tel que le schéma dissipatif HHT qui permet d’atténuer les perturbations numériques engendrées par les hautes fréquences. On adopte ici le schéma HHT avec le paramètre ALPHA=-0.1.
Post-traitement#
Le calcul non-linéaire est long et produit beaucoup de données par défaut. En effet, en non-linéaire, non-seulement vous aurez accès aux déplacements, vitesses et accélérations (DEPL, VITE et ACCE) mais aussi aux contraintes aux points de Gauss (SIEF_ELGA) et aux variables internes (VARI_ELGA). Sur des transitoires longs, les données produites peuvent être très volumineuses. Il est donc recommandé d’utiliser le mot-clef ARCHIVAGE pour sélectionner les instants auxquels vous voulez vos données.
Comme on veut tester l’effet de la plasticité, on va calculer le champ SIEQ_ELGA dans CALC_CHAMP pour vérifier le niveau des contraintes (critère de von Mises). Par ailleurs, on vérifie également l’hypothèse DEFORMATION=’PETIT’ en calculant le champ EPEQ_ELGA dans CALC_CHAMP.
Enfin, à des fins de comparaisons, on peut extraire déplacements, vitesses et accélérations sur le point P de la structure grâce à la commande RECU_FONCTION:
DEPL_DZ = RECU_FONCTION(GROUP_NO=”P”,NOM_CHAM=”DEPL”,NOM_CMP=”DZ”,RESULTAT=DYNADNL)
Cette fonction est ensuite imprimable par IMPR_FONCTION.
On peut aussi le faire dans le module post-process d’AsterStudy ou dans Paravis.
Analyse des résultats – Schéma de Newmark#
A la lecture du fichier de message «.mess», on observe que le nombre d’itération s’accroît lorsque le matériau entre dans un domaine non-linéaire, a contrario des pas élastiques linéaires, où la convergence se fait en une unique itération.
Sur le transitoire, on a les courbes suivantes pour le déplacement, vitesse et accélération au point P de la structure avec les schémas de Newmark et HHT :
Fig. 519 Comparaison des déplacements entre deux schémas Newmark et HHT#
Fig. 520 Comparaison des vitesses entre deux schémas Newmark et HHT#
Fig. 521 Comparaison des accélérations entre deux schémas Newmark et HHT#
On voit qu’avec le schéma de Newmark l’accélération est perturbée par des effets à fréquence élevée. Le schéma HHT permet de corriger un peu ces problèmes. La suite des résultats correspond aux calculs effectués avec le schéma HHT.
Les courbes suivantes montrent l’effet de la plasticité sur le déplacement, vitesse et accélération du point P.
Fig. 522 Effet de plasticité sur la réponse en déplacement du point P#
Fig. 523 Effet de plasticité sur la réponse en vitesse du point P#
Fig. 524 Effet de plasticité sur la réponse en accélération du point P#
On observe que, comparé au calcul élastique : * la plasticité « amortit » l’amplitude de la réponse ; * la fréquence apparente de la structure diminue légèrement.
Ces observations sont un signe de l’assouplissement de la structure sous l’effet de la plasticité, ce qui peut être vérifié par l’évolution de la première fréquence propre de la structure, obtenue avec le mot-clef MODE_VIBR=_F(LIST_INST=LISTR, NMAX_FREQ=1, MATR_RIGI=”TANGENTE”, OPTION=”PLUS_PETITE”) :
Fig. 525 Evolution de la fréquence du premier mode en transitoire#
La figure suivante montre la valeur maximale de la contrainte équivalente de von Mises au cours du temps, parmi tous les points de gauss de l’ensemble d’éléments de la structure. On observe un plafonnement de la contrainte à 200 MPa pour l’analyse plastique.
Fig. 526 Effet de plasticité sur la contrainte von Mises au cours du transitoire#
On vérifie que l’hypothèse des petites déformations est bien correcte. La figure suivante montre la valeur maximale du deuxième invariant du tenseur des déformations (EPEQ_ELGA/INVA_2) au cours du temps, parmi tous les points de gauss de l’ensemble d’éléments de la structure :
Fig. 527 Effet de plasticité sur la déformation au cours du transitoire#
La valeur maximale de la déformation est inférieure à 0.2%. On est bien en petites déformations.
De manière générale en dynamique il est recommandé de regarder le bilan d’énergie. Pour cela, on peut activer le calcul des quantités d’énergie pendant tout le transitoire grâce au mot-clef ENERGIE dans DYNA_NON_LINE. Il est possible de récupérer le bilan d’énergie a posteriori du calcul dans une table. Il faut utiliser la commande *RECU_TABLE *:
tableNrj = RECU_TABLE(CO=DYNADNL, NOM_TABLE=”PARA_CALC”)
Fig. 528 Bilan énergétique au cours du transitoire#
La valeur maximale de la dissipation du schéma d’intégration temporel HHT atteinte tout au long du transitoire a été inférieure à 20 J. Le schéma HHT adopté a permis donc de corriger partiellement les perturbations à haute fréquence observées avec le schéma de Newmark, tout en gardant une dissipation négligeable.
Modélisation F – dynamique linéaire avec non-linéarités localisées#
Objectif#
Ce cas est une introduction à l’usage de Code_Aster pour l’analyse transitoire dynamique avec non-linéarités localisées. Le comportement de choc considéré est de type élastique. La non-linéarité est introduite par une relation déplacement / vitesse / effort aux nœuds portant les non-linéarités. Le calcul est réalisé en s’appuyant sur l’opérateur de dynamique linéaire de Code_Aster DYNA_VIBRA. La résolution s’appuie sur un schéma d’intégration explicite, la non-linéarité est traitée au second membre des équations du mouvement : * à l’instant i, les déplacements sont calculés afin d’attendre l’équilibre avec les forces calculées à l’instant i-1 ; * les forces à l’instant i sont calculées à partir des déplacements ainsi calculé et elles sont appliquées à l’instant i+1 afin d’obtenir les déplacements à i+1. On calcule la réponse de la plaque jusqu’à un temps final de 1 seconde, on considère une raideur normale des chocs égale à KN=1E8N/m.
Cadre du calcul#
Une plaque épaisse d’acier est déformée avant d’être lâchée. Deux points de contact limitent alors le mouvement de la plaque. Les deux points de contact sont disposés sous la plaque, en vis-à-vis desextrémités du bord libre (Figure 7.1). Quand la plaque est horizontale, le jeu est nul entre la plaque et les points de contact. La géométrie de la plaque et le matériau sont les mêmes que ceux déjà utilisées pour les modélisations a,b,c et d.
Fig. 529 Position des points de choc#
Plus précisément, la plaque est pré-déformé de 0,25m selon la direction Z et ensuite laissée libre de vibrer. Le contact aura lieu à chaque fois que, pendant le mouvement vibratoire, les nœuds d’extrémité de la plaque situés sur la surface inférieure se retrouveront dans leur position initiale avant déformation.
Maillage#
Le maillage de la plaque est le même que celui déjà utilisé pour les modélisations a,b,c,d et e. Il faut cependant créer:
Un groupe de maille CHOC_a2 et CHOC_b2 qui contiennent les nœuds d’extrémité de la plaque positionnés sur la surface inférieure en vis-à-vis des points de contact;
LIBRE: il s’agit d’un groupe de mailles (SEG2) à l’extrémité libre de la plaque, reliant les groupes de nœuds CHOC_a2 et CHOC_b2 .
L’affectation du chargement du problème statique se fait via l’opérateur MECA_STATIQUE.
Les conditions aux limites et le chargement sont imposées via le mot clé EXCIT.
On se sert ensuite de la commande CREA_CHAMP afin de créer le champ de déplacement qui constituera ensuite la condition initiale du calcul dynamique. Il s’agit du champ de déplacement calculé au dernier instant su calcul statique. Attention, il est important d’utiliser, dans la construction du champ, la même numérotation des degrés de libertés qui a été utilisée lors de l’assemblage, sans le chargement de déplacement imposé.
Calcul dynamique#
Les étapes communes sont les suivantes:
on calcule la base modale via la commande CALC_MODES;
on calcule les modes statiques via la commande MODE_STATIQUE. Un effort unitaire est imposé au nœuds de choc b2 et a2 afin de calculer la déformée statique de la plaque dans ces conditions. Le champ de déplacements ainsi obtenu est utilisé pour enrichir la base modale de la plaque afin de prendre en compte sa réponse statique. L’enrichissement est fait par la commande DEFI_BASE_MODALE;
la base construite n’est pas constitué de modes orthogonaux. Il faut alors l’orthogonaliser. Pour ce-faire on s’appuie à nouveau sur la commande DEFI_MODALE;
on projette ensuite les matrices de masse et de raideur sur la nouvelle base via l’opérateur PROJ_BASE, mais il faut également projeter le vecteur déplacement issu du calcul statique.
définition de la géométrie des lieux de choc. La commande Aster qui permet de définir la géométrie des lieux de choc est DEFI_OBSTACLE (documentation utilisateur U4.44.21). Le mot clé TYPE permet de définir de manière plus précise la physique du choc. Par exemple, les mots clé PLAN_Y, PLAN_Z, CERCLE et DISCRET définissent la géométrie des contacts entre une structure mobile et un obstacle indéformable. Les types BI_PLAN_Y, BI_PLAN_Z, BI_CERCLE et BI_CERC_INT permettent de définir les lieux de contact possibles entre deux nœuds mobiles. Dans le cas présent d’une plaque qui rentre en contact avec des points fixes dans l’espace, le type PLAN_Z est choisi.
paramétrage de l’opérateur de dynamique linéaire DYNA_VIBRA (documentation utilisateur U4.53.03). Le paramétrage de DYNA_VIBRA passe par la définition du comportement de choc. Le mot clé comportement doit être présent autant de fois que le nombre des points de contact, donc deux fois dans le cas présent. Pour chaque occurrence du mot clé comportement il faut ensuite indiquer :
quel est le groupe de nœud de la structure auquel le comportement sera affecté (CHOC_a2 à la première occurrence du mot clé comportement, CHOC_b2 à la deuxième)
si le choc a lieu avec frottement ou pas: mot clé FROTTEMENT
la relation de CHOC(DIS_CHOC pour les contacts élastiques)
la raideur selon la direction normale au choc: RIGI_NORM
la géométrie du lieu de choc;
la normale à la direction de choc NORM_OBST;
l’obstacle défini par la commande DEFI_OBSTACLE, mot clé OBSTACLE;
l’origine de l’obstacle via le mot clé ORIG_OBST. Il s’agit de l’origine du repère de choc;
la normale à la direction de choc NORM_OBST;
le jeu, ou la position du plan de choc dans le repère de choc ;
les matrices de masse et de raideur (celles projetées sur la base enrichie et orthogonalisée)
l’état initial (ETAT_INIT) comprenant le déplacement statique projeté sur base modale
le schéma d’intégration (on choisi un schéma type Runge-Kutta, mot clé RUNGE_KUTTA_32), et l’incrément de temps (instant finale et pas. Un pas égal à 0.0001s a été choisi pour l’étude)
les matrices de masse et de raideur (celles projetées sur la base enrichie et orthogonalisée)
In-fine: on reporte les résultats obtenus sur base physique grâce à la commande REST_GENE_PHYS. On choisit pour cela un pas de temps de 0.001s. La définition de la discrétisation en temps pour REST_GENE_PHYS est faite via LIST_DEFI_REEL.
On imprime les résultats dans un fichier .med via la commande IMPR_RESU.
Fig. 530 Géométrie des lieu de choc pour les éléments PLAN_Z.#
RQ 1. Deux plans de choc symétriques par rapport à l’origine du repère de choc sont créés automatiquement. Afin d’éviter qu’il y ait contact de la structure avec le plan qui n’est pas présent dans la modélisation, il est conseillé d’éloigner l’origine du repère de l’obstacle de la structure. Dans le cas présent, il a été choisi de positionner l’origine du repère à 1 m de distance verticale par rapport à la face inférieure de la plaque. De cette manière les plans de chocs, se trouveront à la position initiale, ce veut dire en contact avec CHOC_a2 et CHOC_b2 avant pré-déformation, (un exemple est présenté à la Figure 7.3-2 pour le point a), et le second plan de choc sera assez éloigné de ne pas entrer en contact avec CHOC_a2 et CHOC_b2 à l’état initial.
Fig. 531 Géométrie du lieu de choc en correspondance du point a.#
RQ 2. le choix du pas d’intégration est un élément crucial vis-à-vis de la qualité de la solution lors de l’utilisation des schémas l‘intégration explicites. Une discrétisation temporelle suffisamment raffinée doit être choisie afin d’obtenir une solution convergée. Cependant, les ressources mobilisées par le calcul dépendent du raffinement en temps. Une étude de sensibilité au raffinement temporel est nécessaire afin d’orienter le choix.
Post-traitment et impression des résultats#
On récupère les déplacements et les efforts au cours du temps en correspondance des nœuds de choc grâce à la commande RECU_FONCTION.
On imprime les résultats au point b au format xmgrace grâce à la commande IMPR_RESU.
L’ouverture des fichiers xmrgace se fait par ligne de commande via la commande
xmgrage monfichier
Enfin, on rappelle que grâce à la commande POST_DYNA_MODA_T (documentation Aster U4,84,02), pour chaque non linéarité de choc, n calcule et on imprime dans une table :
les résultats des chocs et les grandeurs associées:
pour chaque choc, on a la valeur de l’instant où la force est maximale, la valeur de la force
maximale et de l’impulsion, la durée du choc, la vitesse d’impact et le nombre de rebonds,
les données globales du choc :
sur l’ensemble des chocs constatés, on précise : le maximum absolu de force de choc, la
valeur moyenne des maxima de forces de choc et l’écart type des extrema de force de
chocs
l’histogramme décrivant les max des forces d’impact :
on a le nombre de classes de l’histogramme, les valeurs de ses abscisses (force min et max
de chaque classe) et la densité de probabilité de la force maximale de chaque classe.
Questions#
Réaliser l’étude en prenant en compte uniquement un mode vibratoire de la plaque et puis enrichir progressivement la solution.
On remarque que tous les modes entre 0 et 500 Hz jouent sur les efforts et les déplacements au point de contact et que, de plus, l’utilisation d’une base modale incomplète peut conduire à erreur dans le calcul de l’effort de choc.
Fig. 532 Déplacement au nœud b calculé avec un seul mode propre#
Fig. 533 Effort au nœud b calculé avec un seul mode propre#
Fig. 534 Déplacement au nœud b calculé avec les modes propres calculés entre 0 et 500 Hz#
Fig. 535 Effort au nœud b calculé avec les modes propres calculés entre 0 et 500 Hz#
Modélisation G – dynamique non-linéaire avec sous-structuration#
Objectif#
Ce TP fait suite à la modélisation E (dynamique non-linéaire). L’objectif de ce TP est de mettre en place un calcul non-linéaire transitoire avec DYNA_NON_LINE en tenant compte de la sous-structuration dynamique.
La sous-structuration dynamique consiste à partitionner la structure en une ou plusieurs sous-structures afin de réduire la taille du modèle numérique et de pouvoir étudier indépendamment chaque sous-structure. Le modèle d’ensemble s’obtient par la suite en assemblant les modèles des différentes sous-structures.
Chaque sous-structure est projetée sur une base de vecteurs judicieusement choisie. Cette projection permet de réduire la taille du modèle de chaque sous-structure et par conséquent la taille du modèle de la structure assemblée.
Le comportement de la sous-structure doit être obligatoirement linéaire. La portion de la plaque dans laquelle le comportement non-linéaire (plasticité) a été observé dans la modélisation E restera représentée en base physique. Le choix de la région du modèle à appliquer la méthode de sous-structuration doit se faire donc judicieusement.
La figure suivante présente les éléments qui ont plastifiés (en rouge) au long du transitoire de la modélisation E, ce qui correspond à la région proche de l’encastrement de la plaque :
Fig. 536 Résultat de la modélisation E avec la zone plastifiée en rouge.#
La moitié du modèle restant linéaire (côté opposé à l’encastrement) sera donc soumise à la procédure de sous-structuration :
Fig. 537 Maillage avec deux GROUP_MA pour la sous-structuration#
Pour la réalisation de cette étude, le maillage forma12g.med a été enrichi avec la création d’un groupe d’éléments (VOL_1) représentant la partie qui sera modélisé en base physique (côté encastrement), un groupe d’éléments (VOL_2) représentant la partie qui sera représentée en base réduite, et un groupe de nœuds (No_Int) à l’interface de ces deux zones du modèle.
Procédure de calcul#
Création du macro-élément dynamique#
La première étape consiste à la création du macro-élément dynamique qui représentera le comportement de la partie condensée de la structure dans le modèle final. Pour cela, la méthode de Craig-Bampton sera appliquée. Cette étape se fait uniquement sur le groupe d’éléments VOL_2 du maillage.
La technique de Craig-Bampton consiste à choisir comme base de projection, des modes normaux et des modes contraints. Les modes normaux sont des modes propres de la sous-structure où on bloque tous les noeuds des interfaces avec les sous-structures adjacentes. Un mode contraint est défini par la déformée statique obtenue en imposant un déplacement unitaire sur un degré de liberté de liaison, les autres degrés de liberté de liaison étant bloqués. La base de projection contient des modes normaux et tous les tous les modes contraints associés à tous les degrés de liberté de liaison.
Pour calculer les modes normaux, on doit donc encastrer le groupe des nœuds de l’interface No_Int avec l’opérateur AFFE_CHAR_MECA. Avant d’assembler les matrices de masse et rigidité nécessaires au calcul des modes normaux, on profite pour créer le chargement appliqué au groupe VOL_2, qui composera également le macro-élément dynamique. On exécute ensuite l’opérateur ASSEMBLAGE avec les options RIGI_MECA, MASS_MECA, AMOR_MECA et CHAR_MECA.
On calcule les modes normaux avec CALC_MODES. Attention, le nombre de modes retenus doit être multiple du nombre de DDL des éléments finis que composent le domaine condensé (3 dans ce cas), mot-clef CALC_FREQ=_F(NMAX_FREQ=30).
Les modes contraints sont calculés avec l’opérateur MODE_STATIQUE, mots-clefs MODE_STAT=(_F(GROUP_NO= »No_Int », AVEC_CMP=(« DX », « DY », « DZ »))).
Ensuite on doit définir l’interface dynamique avec DEFI_INTERF_DYNA, mot-clef TYPE=’CRAIGB’, afin de pouvoir définir la base modale contenant les modes normaux, les modes contraints et l’interface dynamique à partir de l’opérateur DEFI_BASE_MODALE.
On récupère la numérotation de la nouvelle base à l’aide de NUME_DDL_GENE, et on projette le chargement sur cette base avec PROJ_VECT_BASE.
Le macro-élément est construit par l’opérateur MACR_ELEM_DYNA avec les mots-clefs BASE_MODALE et CAS_CHARGE. On utilise ensuite DEFI_MAILLAGE avec le mot-clef DEFI_SUPER_MAILLE pour définir un nouveau maillage à partir du macro-élément dynamique. Ce nouveau maillage doit être ensuite assemblé au maillage de base du TP, forma12g.med, à l’aide de l’opérateur ASSE_MAILLAGE.
Construction du modèle#
Une fois le macro-élément généré et intégré au maillage initial, la deuxième étape de ce TP consiste à construire le modèle de la plaque, qui aura une moitié représentée en base physique, et l’autre par le macro-élément dynamique. La prise en compte du macro-élément dans le modèle se fait à partir du mot-clef AFFE_SOUS_STRUC de l’opérateur AFFE_MODELE.
Une fois le matériau élasto-plastique affecté au groupe de mailles physiques ‘VOL_1’ et les conditions aux limites d’encastrement appliquées au bord de la plaque, il est important, à titre de vérification, de calculer les modes propres du modèle, pour vérifier si le macro-élément a été proprement construit et intégré dans le modèle qui sera utilisé par la suite.
La suite du calcul se déroule comme dans la modélisation E, la seule différence étant que dans l’opérateur DYNA_NON_LINE il faut ajouter le mot-clef SOUS_STRUC avec l’argument CAS_CHARGE pour que le chargement appliqué à la partie du modèle qui a été condensée soit pris en compte.
Une fois le calcul terminé, l’opérateur REST_COND_TRAN permet de récupérer les champs dans la partie condensée du modèle. La suite du post-traitement se déroule de manière identique à celle de la modélisation E.
Résultats#
Les figures suivantes présentent le déplacement, vitesse et accélération au point P du modèle, pour les modélisations E (modèle complet) et G (avec sous-structuration dynamique). On observe que les courbes se superposent.
Fig. 538 Comparaison des déplacements en point P entre deux calculs transitoires non-linéaires direct et par sous-structuration#
Fig. 539 Comparaison des vitesses en point P entre deux calculs transitoires non-linéaires direct et par sous-structuration#
Fig. 540 Comparaison des accélérations en point P entre deux calculs transitoires non-linéaires direct et par sous-structuration#