u2.06.13 Conseils généraux d’utilisation de l’opérateur DYNA_NON_LINE#

Résumé :

Ce document présente l’utilisation de méthodes de résolution transitoires (implicites ou explicites) pour la simulation numérique de problèmes dynamiques non linéaires sur base physique.

L’opérateur généraliste de référence pour ce type de calculs se nomme DYNA_NON_LINE et son usage correct sera facilité par le respect de quelques règles de bonne pratique décrites dans ce document.

Ces conseils d’utilisation couvrent:

la définition correcte du modèle au sens dynamique (dont les conditions initiales et aux limites),
la définition de la discrétisation dont le choix du schéma en temps ([R5.05.05], et aussi voir Bibliographie ),
le choix des modèles d’amortissement,
quelques conseils de post-traitement.

Étant donné la grande diversité des problèmes non linéaires, l’utilisateur pourra très utilement compléter sa lecture avec d’autres références plus spécifiques:

[U2.06.03] : sur la modélisation de l’amortissement,
[U2.06.05] : pour l’interaction sol‑structure (linéaire et non‑linéaire),
[U2.06.09] : pour le mono et multi-appuis en calcul sismique des équipements notamment,
[U2.06.10] : sur les spécificités des études de type génie civil sous chargement sismique,
[U2.06.11] : pour l’utilisation de modèles fluide‑structures couplés avec DYNA_NON_LINE,
[U2.04.07] : usage de DYNA_NON_LINE pour résoudre des problèmes en évolution lente mais fortement non linéaires et qui ont du mal à converger avec STAT_NON_LINE (voir bibliographie),
[U2.07.04] : pour la dynamique transitoire non linéaire d’un modèle partiellement réduit par condensation dynamique avec DYNA_NON_LINE,
[U2.06.32] : pour la modélisation des machines tournantes.

La lecture de la documentation [U2.04.01], qui donne des conseils d’utilisation pour l’opérateur STAT_NON_LINE, est aussi fortement recommandée car va surtout approfondir ici les spécificités liées à la dynamique. Tous les aspects communs à STAT_NON_LINE et DYNA_NON_LINE et qui sont détaillés dans la documentation [U2.04.01], comme le choix des paramètres de l’algorithme de Newton, restent valables en dynamiques et ne sont donc pas repris ici.

Définition d’un problème adapté à la dynamique transitoire#

Dans ce chapitre, nous allons proposer des conseils pour la mise au point de la modélisation d’un problème de mécanique dont on souhaite faire la simulation numérique transitoire non linéaire avec Code_Aster. Le chapitre suivant traitera de la discrétisation en espace et en temps, alors que dans ce chapitre seul le modèle continu sera abordé.

Les non-linéarités en dynamique transitoire avec DYNA_NON_LINE sont : * les grands déplacements et grandes transformations, * les relations constitutives non linéaires (plasticité, endommagement, viscoplasticité…), * les non-linéarités d’interface : contact-frottant, collisions…

Modèle#

La première étape est la définition d’un modèle cohérent avec l’hypothèse d’évolution transitoire.

Ainsi, tout système mécanique ne devant pas avoir de mode à énergie cinétique nulle, il convient de s’assurer que la masse volumique est définie en tout point du modèle continu. De même, si l’on veut introduire des éléments discrets pour simuler des corps solides, il faut leur associer une masse.

De même, pour les éléments finis de type coques ou plaques, l’utilisateur peut avoir à s’assurer que tous les degrés de liberté, en particulier les rotations, ont un terme de masse associé. Pour vérifier cette condition, on conseille de spécifier l’option INER_ROTA=”OUI” lors de l’appel à la commande [AFFE_CARA_ELEM] pour tous les modèles COQUE.

Certaines modélisations n’ayant pas de matrice de masse (ce calcul n’a pas été programmé) ne peuvent être utilisées directement avec DYNA_NON_LINE, comme les modèles non locaux. Cette difficulté peut être contournée en leur superposant un modèle classique permettant de représenter l’énergie cinétique (matrice de masse) : les deux modèles s’appuyant sur les mêmes nœuds.

L’usage de certains artifices courants en quasistatique, comme des zones très raides (matériau fictif ayant un très grand module de Young) pour tenir compte de renforts que l’on ne souhaite pas représenter finement géométriquement, par exemple, peut engendrer des perturbations en dynamique. En effet, ce matériau très raide peut générer des oscillations hautes fréquences et des propagations d’ondes dont la célérité est non physique. De plus, avec un schéma en temps explicite, ces zones très raides risquent de faire chuter la valeur du pas de temps critique (condition de stabilité CFL, cf. [R5.05.05]).

En plus de la contribution inertielle, le système aura une contribution dissipative, donc amortissante. En non linéaire, de la dissipation peut être due à la relation de comportement (plasticité…), au frottement dans les liaisons …

Habituellement, la méconnaissance détaillée de tous les mécanismes dissipatifs dans le système est compensée par une représentation simplifiée qui permet de définir un amortissement global. Dans Code_Aster, on dispose de deux modèles d’amortissement visqueux globaux (que l’on peut associer dans la modélisation à des amortisseurs sur des éléments discrets, cf. commande [AFFE_CARA_ELEM] (mot-clé facteur DISCRET ou DISCRET_2D, etc…) : le modèle de Rayleigh et l’amortissement modal.

En pratique, dans les deux cas, il est indispensable d’avoir effectué un calcul modal préalable. En effet, l’amortissement modal est défini sur les modes qui sont alors des arguments de DYNA_NON_LINE. Pour l’amortissement de Rayleigh, la seule manière d’appréhender simplement son sens physique, c’est de le recaler sur des valeurs d’amortissement modal (qui peuvent venir de résultats expérimentaux). Nous reviendrons plus en détail sur l’amortissement de Rayleigh dans le chapitre suivant Discrétisation du problème continu.

D’une manière générale, il faut garder à l’esprit que tous les modèles d’amortissement global ont été conçus pour les cas linéaires et qu’ils servent alors à prendre en compte certains phénomènes dissipatifs comme le frottement ou des liaisons imparfaites. Donc, si on rajoute au modèle la prise en compte de ces non-linéarités, alors qu’on garde les valeurs d’amortissement global recalées sur une approche linéaire, cela peut conduire à avoir trop de dissipation dans le système non linéaire en particulier en basses fréquences.

Idéalement, la prise en compte de toutes les non-linéarités, associée à une discrétisation fine du système, devrait permettre de ne pas avoir à utiliser en plus un amortissement global forfaitaire. Dans la réalité, on est obligé de négliger certains aspects non linéaires pour des raisons de coût CPU et/ou de méconnaissance sur les mécanismes mis en œuvre, et dans ce cas, l’amortissement global a son rôle à jouer. La documentation [U2.06.03] présente plus en détails son usage.

Conditions aux limites et excitations#

La première décision à prendre pour résoudre le transitoire est le choix du référentiel de configuration (le champ de déplacements) : le cas général est un référentiel galiléen (repère absolu). Dans des cas fréquents où le système est entraîné dans un mouvement de translation rigide, on pourra résoudre le transitoire en référentiel relatif. Le calcul se fera en base encastrée, voir Mono-appui. Si le système est sollicité par divers appuis excités différemment, on est dans le cas Multi-appuis et on devra résoudre le transitoire en référentiel absolu (sauf si le comportement est totalement en HPP élastique linéaire).

En quasistatique, il est courant d’imposer des conditions aux limites dont les évolutions en temps sont des fonctions simples à définir, comme des rampes ou des fonctions continues affines par morceaux.

Il est fréquent qu’en dynamique les excitations soient décrites par le produit d’une fonction de l’espace et d’une rampe ou fonction du temps (ou signal en temps). On utilisera ainsi le mot-clé facteur : EXCIT.

Il est recommandé de définir ce signal en temps de sorte de partir et de finir à zéro. Il peut être utile aussi que la liste d’instants choisie pour le transitoire avec DYNA_NON_LINE dure plus longtemps que le signal immposé de sorte de laisser le système étudié revenir au repos par effet de l’amortissement qui a été choisi.

En dynamique, si l’on veut éviter les oscillations parasites de la solution, il faut absolument imposer des quantités suffisamment régulières en temps : donc au moins continûment dérivables, par morceaux. Pour arriver à cela simplement, on peut choisir de se définir des fonctions du temps polynomiales à la place des fonctions affines. D’éventuelles irrégularités dans les chargements imposés peuvent être partiellement compensées par l’utilisation d’un modèle d’amortissement adapté (ou en particulier grâce à un schéma d’intégration en temps dissipatif de type HHT, voir Schéma d’intégration en temps). Néanmoins, des oscillations excessives en non linéaire peuvent compromettre irrémédiablement la suite des calculs car la solution (dont les variables internes) dépend de son historique.

Sur le schéma Lissage d’une évolution imposée pour l’adapter à un calcul transitoire on montre une régularisation possible d’un courbe affine (en pointillés) par une courbe adaptée au calcul dynamique. On peut noter la nécessité, sur cet exemple, de devoir faire commencer l’évolution dynamique en un temps négatif, si l’on veut garder des valeurs comparables pour les instants positifs (en dehors des zones de régularisation).

../../../../_images/Shape138.gif — Fig. 54 Lissage d’une évolution imposée pour l’adapter à un calcul transitoire#

Dans tous les cas où l’on peut calculer le spectre des signaux imposés, la règle de Shannon pour définir l’échantillonnage n’est pas assez stricte pour les calculs non linéaires. En pratique, on conseille d’avoir entre 4 et 10 pas de temps sur la période la plus courte de tous les signaux imposés.

Toute discrétisation insuffisante en temps risque de se traduire par des sur-oscillations hautes fréquences de la solution. Il existe alors un palliatif à cela en choisissant certains schémas en temps, mais au risque d’introduire trop d’amortissement dans le système.

Mono-appui#

Les chargements de type mono-appui (entraînement par un mouvement de translation rigide), qui sont couramment employés lors des calculs sismiques, se caractérisent par des particularités qu’il convient de rappeler ici.

Tout d’abord, la résolution, selon l’hypothèse mono-appui, se fait nativement dans le repère relatif lié au support (radier…). En post-traitement, il faut donc recombiner avec le mouvement d’entraînement pour obtenir les champs cinématiques absolus. Pour l’accélération, c’est en général immédiat car l’accélération d’entraînement est la donnée d’entrée connue sous la forme des accélérogrammes employés.

En revanche, pour reconstruire le déplacement absolu, il faut disposer du déplacement d’entraînement. Dans le cas classique où l’on ne dispose que d’accélérogrammes, on est obligé d’intégrer deux fois en temps pour calculer le déplacement correspondant (cette intégration peut se faire avec l’opérateur [CALC_FONCTION], mot-clé INTEGRE). L’utilisateur doit alors être très attentif aux problèmes de l’intégration numérique de signaux échantillonnés tels des constantes d’intégration qui peuvent induire de la dérive ou des oscillations à très basse fréquence. L’opérateur [CALC_FONCTION] dispose d’options pour corriger ces dérives (mot-clé CORR_DEPL et CORR_ACCE), mais on recommande très fortement de bien contrôler la qualité du résultat intégré final car, dans certains cas, ces corrections sont insuffisantes, y compris la correction par filtrage passe‑haut de [CALC_FONCTION] (avec CORR_ACCE, METHODE = 'FILTRAGE'). Une correction spécifique complémentaire est alors indispensable. Le passage par la FFT pour faire l’intégration en fréquence peut alors constituer une bonne alternative permettant d’aboutir à des signaux intégrés plus corrects ou au moins plus faciles à corriger (juste un décalage d’une constante, par exemple, mais sans dérive).

Multi-appuis#

Parmi les types de chargements qui peuvent demander une adaptation à DYNA_NON_LINE, on peut évoquer le cas du multi-appuis (en accélérations imposées sur les supports) pour les études sismiques.

Cette méthode, originellement développée pour les calculs transitoires en hypothèse linéaire, repose sur la définition de modes statiques élastiques aux appuis. En non linéaire, tant que ces modes restent pertinents, c’est-à-dire que les non-linéarités sont suffisamment faibles pour ne pas demander leur réactualisation, l’approche multi-appuis usuelle reste acceptable.

En revanche, dans le cas général, on doit remettre en cause l’utilisation de modes statiques élastiques calculés initialement. Dans ce cas, pour éviter leur usage, la sollicitation multi-appuis doit être prise en compte en imposant aux appuis les déplacements correspondants à l’intégration des accélérogrammes qui sont des données d’entrée de la méthode multi-appuis linéaire usuelle.

La méthode est alors rigoureuse, quel que soit le type de non-linéarité, mais l’intégration des accélérogrammes doit être menée avec précaution, pour ne pas aboutir à des signaux incorrects.

En dehors de l’échantillonnage, il est indispensable, entre autres, de vérifier l’absence de dérive du signal (une correction peut être faire avec l’option CORR_ACCE de [CALC_FONCTION]). Dans Code_Aster , on peut effectuer cette intégration facilement en se servant de la FFT qui est disponible dans [CALC_FONCTION]. C’est aussi l’occasion d’analyser le contenu fréquentiel du signal. On a alors la relation, dans le domaine fréquentiel: \(U=-A/{\omega}^{2}\) .

On peut aussi se servir de l’option INTEGRE de [CALC_FONCTION]. Pour plus d’informations sur la mise en œuvre des sollicitations de type mono ou multi-appuis, la lecture des documentations [U2.06.09] et [U2.06.10] est toute indiquée.

Contact et collisions#

Les chargements issus du contact par pénalisation peuvent aussi perturber la solution par des oscillations hautes fréquences, qui sont liées à la valeur du coefficient de pénalisation. Ces artefacts peuvent être réduits avec l’usage d’un schéma d’intégration en temps dissipatif en hautes fréquences.

On peut également se doter d’une liste d’instants de calcul qui gère les difficultés par découpage, cf. la commande [DEFI_LIST_INST] avec le mot-clé facteur ECHEC= _F(EVENEMENT = 'DELTA_GRANDEUR', ...), de type « event-driven », en testant si l’incrément d’une composante d’un champ (de vitesse par exemple au point de contact) dépasse un seuil fixé.

Si l’on veut utiliser un schéma d’intégration en temps explicite, un coefficient de pénalisation trop grand va entraîner une chute du pas de temps critique. Comparé à la statique, il peut s’avérer obligatoire de baisser le coefficient de pénalisation, avec l’inconvénient d’augmenter l’interpénétration des solides lors du contact.

On pourra consulter les cas-test : * sdnv100, impact d’une poutre sur une paroi rigide [V5.03.100], * sdnv105, balancement d’un bloc sur une table [V5.03.105], * sdnv301, collision hertzienne de deux billes élastiques [V5.03.301].

Remarque: Pour les schémas en temps explicites ():ref:u2.06.13-Schema_integration_temps), si l’on impose des conditions aux limites en déplacement qui évoluent au cours du temps, il faut tenir compte du fait que ces conditions sont en fait imposées en accélération en explicite (car c’est l’inconnue primale). Cela signifie que l’on doit entrer dans DYNA_NON_LINE la dérivée seconde du signal en déplacement que l’on veut imposer. Cette évolution du déplacement imposé doit donc être dérivable au moins deux fois en temps

Conditions initiales#

Dans la même logique que pour les conditions aux limites, on recommande d’éviter les conditions initiales singulières, comme ce que peuvent générer des chargements imposés ayant des évolutions en temps de type Dirac ou Heaviside.

On privilégiera les conditions régulières donc avec des dérivées (première et si possible seconde) en temps en \({t}_{0}\) nulles.

Si l’on veut partir d’un état précontraint (non vierge), il est tout à fait possible et recommandé d’enchaîner un calcul dynamique à un calcul quasistatique. Dans ce cas, l’état initial (mot-clé facteur ETAT_INIT) doit être alors parfaitement équilibré (déplacements, contraintes, variables internes, efforts appliqués à cet instant) et, sous réserve de régularité dans les conditions imposées au moment du passage statique vers dynamique, l’excitation dynamique devant démarrer avec un signal nul, la solution transitoire ne devra pas présenter d’oscillations non physiques. Cette méthode d’enchaînement de la quasistatique vers la dynamique permet la prise en compte aisée de toutes les précharges, comme par exemple la pesanteur. Par exemple, dans une simulation de lâcher, avec un déplacement initial imposé non nul, alors il faut renseigner les contraintes et variables internes correspondant à cet instant.

Si l’on avait voulu prendre en compte cet état initial précontraint avec uniquement une résolution dynamique, il aurait fallu introduire un fort amortissement initial et attendre que les oscillations se soient dissipées pour avoir l’état statique précontraint. Ensuite, on peut continuer le calcul avec l’amortissement physique. On est obligé de procéder ainsi lorsqu’on utilise un code de calcul qui ne traite pas la quasistatique non linéaire ( cf. Europlexus, voir [bib2]).

Dans certaines situations, on doit enchaîner une simlation en statique, avec des conditions aux limites bloquées puis en dynamique sous l’action d’un champ d’ondes incidentes avec des conditions aux limites de frontières absorbantes et libérées. On consultera [U2.04.08].

Lors de l’initialisation du schéma d’intégration pour le calcul transitoire non linéaire, on cherche donc à inverser la matrice de masse. Si elle est singulière, alors un message avertit l’utilisateur et on met arbitrairement l’accélération initiale à zéro. On peut craindre alors la présence d’oscillations parasites au début du transitoire ; on peut les réduire en adoptant un pas de temps plus petit au tout début du transitoire. D’une manière générale la non-inversibilité de la matrice de masse doit amener l’utilisateur à vérifier son modèle, sauf si c’est volontaire.

Lors de poursuites, deux aspects sont sensibles.

D’une part, il vaut mieux éviter d’avoir de trop brusques variations de pas de temps: la cohérence de la solution au passage de la poursuite en serait affectée. Un fort changement de pas de temps peut être vu comme un filtre. Un grand pas de temps constitue un filtre passe-bas d’une solution initialement calculée avec un pas de temps fin.

D’autre part, si l’on veut changer de type de schéma en temps, certaines règles sont à respecter. Si l’on veut passer d’une méthode d’intégration temporelle implicite (quasistatique ou dynamique) à un schéma explicite, la poursuite sera mathématiquement valide car l’état initial sera équilibré (au résidu près). En revanche, la bascule inverse directe introduit une erreur car l’état initial, venant d’un calcul explicite, ne respectera pas l’équilibre au sens implicite. En effet, cet opérateur va chercher à revérifier l’équilibre: on résout l’équilibre en inversant la matrice de masse, ce qui donne l’accélération initiale. Si ce champ est non nul, cela traduit la non-vérification, au sens statique, de l’équilibre, à l’instant initial. En pratique, cette imprécision peut engendrer des oscillations de la solution (phénomène d’overshooting). D’une manière plus subtile, la bascule du \(\theta\) -schéma vers un schéma du second ordre comme l’accélération moyenne NEWMARK ou HHT complet va introduire une petite erreur sur les termes d’accélération, ce qui peut perturber la solution numérique. Si, malgré tout, on veut faire ces bascules qui induisent des erreurs, on peut minimiser les imprécisions en choisissant de faire ces bascules lors de phases où l’excitation donc la solution évolue peu. En résumé :

Tableau 25 Enchaînements de calculs transitoires#
Origine	Destination	Possible
STAT_NON_LINE	DYNA_NON_LINE	oui
DYNA_NON_LINE avec schéma en temps implicite	DYNA_NON_LINE	oui
DYNA_NON_LINE avec schéma en temps explicite	DYNA_NON_LINE avec schéma en temps explicite	oui
DYNA_NON_LINE avec schéma en temps explicite	DYNA_NON_LINE avec schéma en temps implicite	non

Discrétisation du problème continu#

En complément des conseils précédemment donnés pour le modèle continu, ce chapitre va lister les aspects les plus importants à respecter pour obtenir un modèle discrétisé adapté à [DYNA_NON_LINE].

Maillage#

Comme préalable au calcul transitoire, il est fortement recommandé de mener un calcul modal (avec [CALC_MODES]), afin d’obtenir des informations modales qui vont permettre de qualifier la qualité du modèle en dynamique et d’ajuster certains paramètres. L’objectif n’étant pas de rentrer dans les détails de l’analyse modale, on peut néanmoins rappeler quelques règles.

En général, on peut se définir une fréquence de coupure pour le problème à étudier, et donc une troncature modale associée. La bonne représentation de tous les modes de cette base tronquée peut donner des indications sur les tailles de mailles à employer, en plus des considérations déjà prises en compte pour les calculs quasistatiques. En gros, une dizaine de mailles par longueur d’onde la plus petite est suffisante (à ajuster suivant la richesse des éléments, bien sûr).

L’analyse modale va aussi permettre de vérifier que le modèle est exempt de problèmes comme des contributions non définies à l’inertie ou à la raideur.

Enfin, elle est indispensable pour l’utilisation de l’amortissement modal dans DYNA_NON_LINE ou pour recaler l’amortissement de Rayleigh, comme on va le voir dans ce qui suit ( cf. Modélisation de l’amortissement ).

D’une manière bien plus marquée que pour les calculs quasistatiques, la résolution dynamique s’accommodera assez mal de maillages présentant de brutales variations de tailles d’éléments. En effet, ces zones peuvent s’assimiler à des interfaces qui vont perturber la propagation des ondes. On peut alors voir apparaître des ondes réfléchies qui se superposent aux trains d’ondes « physiques ». De même, si l’on veut une bonne représentation des ondes, le maillage doit être fin sur tout le trajet des ondes: on ne peut se contenter de raffiner que dans certaines zones non-linéaires. Si une onde part d’une zone maillée finement pour aller vers une zone maillée plus grossièrement, elle va subir un filtre et l’onde réfléchie par le bord opposé risque d’être fortement perturbée, voire de disparaître. Si l’on ne s’intéresse qu’à des durées courtes, donc avant tout retour d’onde sur la zone maillée finement, alors cette erreur sur l’onde réfléchie n’est pas pénalisante. En revanche, si l’on veut calculer des solutions sur des durées plus longues, l’erreur commise pourra être non négligeable.

Enfin, contrairement au cas quasistatique, le temps a un sens physique et les échelles de temps que l’on veut analyser sont fortement couplées aux échelles en espace du problème discrétisé. Ainsi, le pas de temps est lié à la taille de maille, ce que l’on peut percevoir immédiatement avec la notion de condition de stabilité CFL pour les schémas explicites.

Schéma d’intégration en temps#

Le temps ayant un sens physique en dynamique, la qualité de sa discrétisation en est d’autant plus sensible. Un schéma d’intégration a pour rôle d’exprimer par une expression approchée deux relations entre les grandeurs (DEPL, VITE, ACCE) à l’instant précédent et celles à l’instant suivant. Les équations du mouvement constituent le troisième paquet d’équations qui permet de fermer le système à résoudre. Selon la grandeur principale (DEPL, VITE, ACCE) que l’on choisit comme pilote pour exprimer ces équations on adopte un type de FORMULATION : DEPLACEMENT, VITESSE, ACCELERATION.

On peut énoncer quelques règles pour le choix du pas de temps :

l’évolution des chargements imposés doit être échantillonnée de manière suffisamment fine (entre 5 à 10 pas de temps par période la plus courte des signaux considérés),
le comportement modal de la structure doit être bien représenté (comme ci-dessus, on doit avoir entre 5 et 10 pas de temps par période la plus faible des modes considérés jusqu’à la fréquence de coupure), en particulier, la distorsion de fréquence causée par le schéma d’intégration en temps doit être limitée : le choix de prendre entre 5 et 10 pas de temps par période la plus faible des modes considérés jusqu’à la fréquence de coupure assure une distorsion de fréquence inférieure au pourcent.

Étant donné le caractère basse fréquence, au mieux moyenne fréquence, de la plupart des problèmes que l’on peut aborder ici, ces deux règles ne sont, en général, pas très pénalisantes.

Le choix du schéma d’intégration temporelle et de ses paramètres se fait avec le mot-clé facteur : SCHEMA_TEMPS =_F(…).

Les cas-test SDLD31 [V2.01.031] et SDLD32 [V2.01.032] offrent une comparaison des schémas d’intégration temporelle proposés dans DYNA_NON_LINE et dans [DYNA_VIBRA] sur des systèmes linéaires simples sous excitation harmonique.

Schéma d’intégration explicite#

Avec un schéma d’intégration explicite, conditionnellement stable, dédié aux problèmes de « dynamique rapide » en présence de transitoires à contenu « hautes fréquences », il faut de plus respecter la condition de stabilité de Courant (dite « CFL » cf. [R5.05.05] et Bibliographie) sous peine de divergence numérique (« explosion » de l’énergie cinétique). Dans Code_Aster, ces schémas sont traités en FORMULATION ='ACCELERATION' (on veut résoudre en inversant la matrice de masse). La matrice de masse, qui doit être inversible, est préférentiellement choisie sous forme diagonale (lumpée) pour des raisons de performance ; cela s’obtient par le mot-clé : MASS_DIAG = 'OUI'. La matrice de masse lumpée permet de corriger une part de la dérive fréquentielle sur de longues durées provenant de l’erreur en temps induite par le schéma d’intégration. Pour les éléments de structure, par exemple DKT/DKTG/DST, il est nécessaire de préciser dans AFFE_CARA_ELEM, sous le mot-clé facteur COQUE, le mot-clé simple INER_ROTA = 'OUI', de sorte d’avoir des termes d’inertie non nuls sur toute la diagonale de la matrice de masse, même s’ils sont faibles, [R3.07.03].

Cependant cette option n’étant pas disponible pour tous les éléments finis, ni avec tous les éléments discrets, il faut alors résoudre avec la matrice de masse consistante, comme en implicite. Pour les éléments de modélisation HM, la matrice de masse a des termes nuls pour certains degrés de liberté spécifiques, par exemple la pression. Dans ce cas aussi, on doit résoudre avec la matrice de masse consistante, à laquelle on ajoute une contribution faible de la matrice de raideur avec le mot-clé COEF_MASS_SHIFT.

On dispose de deux schémas explicites :

différences centrées (DIFF_CENT) qui est non dissipatif,
Tchamwa-Wielgosz (TCHAMWA) qui est dissipatif, d’une manière comparable à HHT.

On préconise de commencer par utiliser un schéma non dissipatif.

Le calcul automatique de la condition de Courant sur base modale est exact et sera toujours valable quel que soit le type d’élément fini utilisé.

Pour le schéma d’intégration explicite de différences centrées DIFF_CENT, le pas de temps critique vaut \(2/{\omega}_{max}\) avec \({\omega}_{max}\) qui est la plus haute pulsation propre du système discrétisé. La présence de raideurs de choc associées à la modélisation pénalisée de contact frottant ou non accroît significativement cette plus haute pulsation propre. En présence d’amortissement physique (coefficient d’amortissement réduit \(\xi\)) le pas de temps critique a pour expression : \(2 \left(\sqrt{1+\xi^2}-\xi\right)/{\omega}_{max}\).

On peut calculer cette pulsation avec CALC_MODES en choisissant l’option 'PLUS_GRANDE' (ou avec l’option PLUS_PETITE en inversant les rôles de la matrice de masse et de raideur).

Pour le schéma TCHAMWA d’intégration explicite dissipatif de Tchamwa-Wielgosz (voir [bib11]), le pas de temps critique est légèrement plus faible que celui du schéma DIFF_CENT et décroît quand on augmente l’amortissement numérique lié au schéma (paramètre PHI, dont la valeur par défaut est \(1.05\)).

La condition de stabilité de Courant peut aussi être approchée, tout au moins sur des modèles avec éléments finis massifs, par \(\Delta t={l}_{\min}/c\) avec \({l}_{\min}\) qui est la plus petite longueur du modèle discrétisé et \(c\) la célérité des ondes de traction au point considéré.

L’opérateur DYNA_NON_LINE se sert de cette formule pour donner une approximation de la condition de stabilité de Courant. Il existe cependant certaines limitations :

on ne sait pas calculer automatiquement la condition de stabilité de Courant associée à la présence de ressorts discrets (ce n’est pas programmé),
on ne corrige pas la formule pour les éléments de structures (coques, plaques et poutres).

Il peut donc exister des cas où la valeur renvoyée n’est pas un minorant de la vraie condition de Courant. Donc, en cas de divergence, il convient de diminuer le pas de temps.

De plus, le calcul de la condition de stabilité de Courant n’est pas réactualisé et ne se fait qu’au début du calcul, se basant sur la célérité des ondes élastiques pour l’état initial. Si le module élastique baisse (endommagement), la condition de stabilité de Courant initiale peut devenir trop sévère. Il n’y a pas de risque de divergence (sauf pour des matériaux dont le module élastique pourrait augmenter), mais le temps CPU pourrait être un peu diminué en réactualisant la condition de stabilité de Courant (comme c’est fait dans les codes dédiés à la dynamique rapide, comme Europlexus [bib2]).

Pour la plupart des structures, la condition de stabilité de Courant est très pénalisante : la célérité des ondes étant souvent de l’ordre de quelques milliers de \(m/s\) , on arrive à des pas de temps de moins de \({10}^{-5}s\) , pour des dimensions d’éléments finis choisis pour des structures habituelles. Cela impose donc un nombre total d’instants très élevé pour décrire le transitoire, donc cela produit un volume considérable de résultats.

Si l’on résout le problème sur base modale réduite (mot-clé facteur PROJ_MODAL =_F()), alors l’inconvénient du très faible pas de temps critique pour un schéma explicite disparaît. En effet, la pas de temps limite sera directement proportionnel à la plus petite période propre de la base modale tronquée. On a la relation: \(\Delta t=2/\omega\) avec \(\omega\) qui est la plus haute pulsation propre du système. Plus la base modale sera tronquée, plus le pas de temps critique associé sera grand.

De plus, le calcul de la condition de Courant est alors immédiat car on connaît explicitement toutes les pulsations de la base, donc la plus haute en particulier.

Avec un schéma explicite, on préconise de se placer légèrement en dessous de la condition de Courant: entre 0,5 et 0,7 fois la condition de stabilité de Courant. Il n’y aura pas de sub-division du pas de temps pour cause de non convergence à l’échelle de l’équilibre global : l’utilisateur reste seul maître du pas de temps tout au long du calcul. Par contre, via le mot-clé facteur COMPORTEMENT, cf. [U4.51.11], l’utilisateur indique ses exigences pour l’intégration de la loi de comportement non linéaire, qui peuvent conduire au découpage automatique défini dans la commande [DEFI_LIST_INST].

Remarques

si l’on impose des conditions aux limites en déplacement qui évoluent au cours du temps, il faut tenir compte du fait que ces conditions sont en fait imposées en accélération en explicite (car c’est l’inconnue primale). Cela signifie que l’on doit entrer dans DYNA_NON_LINE la dérivée seconde du signal en déplacement que l’on veut imposer. Cette évolution du déplacement imposé doit donc être dérivable au moins deux fois en temps
L’usage des éléments quadratiques n’est pas recommandé en explicite : en effet, des oscillations parasites peuvent apparaître sur les champs solutions (déplacement ou contraintes). Ce phénomène peut aussi être amplifié par l’utilisation conjointe de la matrice de masse consistante
Pour la vérification de l’hypothèse de contraintes planes, si l’on utilise la méthode de de Borst, le processus itératif global associé est alors bloqué à la première itération (contrairement au cas implicite). Cela peut entraîner une légère imprécision des résultats, tempérée par le fait que le pas de temps étant par ailleurs très petit, la non-planéité des contraintes ne peut qu’être très modérée. Cependant, il est tout à fait possible d’utiliser la version locale de l’algorithme de de Borst par le biais du mot-clé ITER_CPLAN_MAXI sous le mot-clé facteur COMPORTEMENT, cf. [Comportements non-linéaires]. En effet, ces itérations locales ont bien lieu en explicite, avec un surcoût de calcul associé.

Schéma d’intégration implicite#

Avec un schéma d’intégration implicite, qui permet de respecter au mieux les conditions d’équilibre dynamique non linéaire, dédié aux problèmes de « dynamique lente » en présence de transitoires à contenu « basses fréquences », dans une phase de mise au point de la simulation transitoire, il est recommandé de faire un comparatif de résultats obtenus pour deux choix de valeurs de pas de temps. En implicite, les schémas proposés sont inconditionnellement stables, mais cela ne signifie pas qu’on peut prendre un pas de temps quelconque ! Un pas trop grand n’amènera pas une divergence, mais l’erreur sur la solution comme la distorsion de fréquence, sera évidemment importante.

Ces schémas sont accessibles en FORMULATION="DEPLACEMENT", FORMULATION ="VITESSE", FORMULATION ="ACCELERATION". La FORMULATION ="DEPLACEMENT" est tout à fait naturelle et privilégiée dans le cas de non-linéarités matérielles, car les forces internes résultent du traitement de la relation de comportement exprimée à partir des déformations elles-mêmes obtenues avec le champ de déplacement. Il n’est pas possible d’utiliser une matrice de masse lumpée avec ces schémas.

Il est recommandé d’assurer un échantillonnage uniforme et régulier du pas de temps dans la mesure du possible. Cependant, en présence de non-linéarités de comportement, il n’est pas possible ni souhaitable d’éviter le redécoupage temporel en cas de non convergence de l’algorithme de Newton, piloté par les directives de subdivision du pas que l’on aura définies dans la commande [DEFI_LIST_INST].

On peut classer les schémas implicites en deux catégories:

accélération moyenne (NEWMARK) d’ordre 2 et qui n’apporte pas de dissipationnumérique: à utiliser en premier,
HHT complet (MODI_EQUI = 'OUI', option par défaut) qui reste d’ordre 2, contrairement au cas de l’accélération moyenne modifiée. Si l’on désire plus d’amortissement en moyenne fréquence, alors le schéma d’accélération moyenne modifié peut être employé avec MODI_EQUI = 'NON', et le schéma devient d’ordre 1. Le schéma HHT complet (MODI_EQUI = 'OUI') est spécifiquement développé pour introduire un amortissement numérique haute fréquence et donc ne pas perturber la réponse physique basse fréquence. L’amortissement est directement piloté par le paramètre ALPHA du schéma,

Si l’on observe des oscillations hautes fréquences dans la solution numérique (en gros, des oscillations dont la période est de l’ordre de quelques pas de temps), on préfèrera choisir le schéma HHT complet, pour commencer avec une valeur de l’ordre de \(–0,1\) pour le paramètre ALPHA. Une valeur de \(–0,3\) constitue une limite haute encore utilisable. Il est fortement recommandé de mener une étude paramétrique sur la valeur du paramètre ALPHA afin de choisir la valeur la plus proche de zéro qui permette d’obtenir une solution correcte sans sur‑amortissement numérique.

On se reportera au paragraphe §5.4 de [R5.05.05].

Les schémas implicites sont à utiliser, prioritairement, avec une formulation en déplacements : FORMULATION = 'DEPLACEMENT'.

Synthèse pour le choix du schéma d’intégration#

Le tableau suivant résume les points essentiels.

Tableau 26 **Tableau : schémas d’intégration.**#
nom	`NEWMARK`	`HHT`	`DIFF_CENT`	`TCHAMWA`
type	implicite	implicite	explicite	explicite
applications	simulations de « longue durée » pour des excitations en basse fréquence	idem	simulations de courte durée	idem
FORMULATION	`"DEPLACEMENT"`	`"DEPLACEMENT"`	`"ACCELERATION"`	`"ACCELERATION"`
matrice de masse	cohérente	cohérente	diagonale préférentiellement	diagonale préférentiellement
stabilité	inconditionnelle	inconditionnelle	conditionnelle	conditionnelle
ordre	ordre 2 sauf si \(\gamma \neq 1/2\)	ordre 2	ordre 2	ordre 1
dissipation numérique	non si \(\gamma=1/2' et :math:\)beta=1/4”	oui en HF si \(\alpha<0\)	non	oui en HF si \(\varphi>1\)
distorsion de fréquence	en \(-{(\omega_{max} \Delta t)}^2\)	en \(-(\omega_{max} \Delta t)^2\)	en \(+(\omega_{max} \Delta t)^2\)
erreur d’amplitude	en \(\omicron (\Delta t )^2\)		en \(\omicron (\Delta t )^4\)
choix du pas de temps	limiter la distorsion de fréquence	limiter la distorsion de fréquence	garantir la stabilité	garantir la stabilité

Pour finir, il convient de signaler que les résultats analytiques sur les caractéristiques (convergence, distorsion de fréquence, erreur, amortisement numérique, overshooting…) des schémas en temps sont obtenus pour un cadre linéaire. Les démonstrations en régime non linéaire sont très rares et se cantonnent à des cas particuliers. En pratique, certaines caractéristiques des schémas peuvent se dégrader en non linéaire. Cela peut expliquer pourquoi il n’est pas forcément indispensable d’utiliser le schéma qui, en linéaire aura des performances exceptionnelles (ordre 4…), mais qu’il vaut mieux privilégier des schémas plus simples et plus robustes, en particulier avec de la dissipation numérique haute fréquence.

De même, en non linéaire, l’évaluation précise de la condition de stabilité de Courant d’un schéma explicite demande une réactualisation de son calcul. En effet, la condition de stabilité de Courant calculée initialement peut ne pas être conservative (par exemple si certains éléments voient leur taille diminuer, ou si des chocs se produisent, avec une modélisation par pénalisation). On ne réactualise pas ce calcul et en cas de divergence, il est recommandé de diminuer le pas de temps et de relancer le calcul.

Les résultats d’analyse linéaire sur les schémas constituent cependant une solide base pour leur analyse ( cf. [R5.05.05] et Bibliographie), tout en sachant que les non-linéarités peuvent perturber le comportement des schémas.

Le pas de temps à imposer pourra être judicieusement borné :

en valeur supérieure, par le pas de temps à respecter pour bien discrétiser les évolutions des chargements imposés et pour bien représenter la plus haute fréquence propre du système dont on veut tenir compte,
en valeur inférieure, par la condition de stabilité de Courant, au sens qu’elle a de plus court temps pour lequel une information peut passer d’un nœud de la maille à un autre.

Entre ces bornes, il est indispensable de mener une étude paramétrique pour s’assurer de la bonne convergence de la solution numérique.

La borne inférieure donne une information cruciale quant aux sous-divisions maximales du pas de temps qu’il faut autoriser lorsque l’algorithme ne converge pas. Laisser le pas de temps se sous-diviser jusqu’à aller nettement sous la condition de stabilité de Courant peut ne servir à rien car la solution calculée risque alors de subir une pollution numérique qui ne facilitera pas la convergence.

Toujours en lien avec les itérations pour la vérification de l’équilibre à chaque pas, on peut remarquer que, dans la plupart des cas, si le pas de temps est suffisamment fin, le nombre maximal d‘itérations à convergence reste modéré: souvent de l’ordre de 10, alors qu’en quasistatique, on peut couramment dépasser ces valeurs.

Donc, l’idée consiste à dire que le pas de temps est d’un bon ordre de grandeur si le nombre d’itérations à convergence reste modéré. Si ce nombre augmente, on peut tenter de réduire légèrement le pas, en respectant toujours les bornes définies ci-dessus. Il existe néanmoins des cas ou l’on peut avoir, ponctuellement, besoin d’autoriser plus d’itérations lors de quelques pas.

Pour finir, si l’on est obligé de changer de schéma en temps pour utiliser un schéma dissipatif, comme HHT, il est indispensable de mener une étude paramétrique sur cet amortissement. En effet, le risque d’introduire une dissipation trop grande n’est pas négligeable, surtout avec le schéma d’accélération moyenne modifiée. Le paragraphe suivant reviendra sur ce point.

Conditions initiales#

On privilégiera des conditions initiales régulières : déplacements et vitesses dérivables en \(t=t_0\), état de contraintes en équilibre à \(t=t_0\) par exemple résultant d’un calcul quasistatique avec [STAT_NON_LINE]. En effet, de nombreux schémas en temps évaluent un prédicteur de l’accélération initiale à partir de l’état initial, en inversant la matrice de masse (si c’est faisable) : si cette accélération initiale est excessive, la perturbation peut se propager sur le transitoire, surtout si l’amortissement est faible et engendrer des oscillations parasite de la solution (« overshooting »). Si la matrice de masse n’est pas inversible, voir la section Modification de la matrice de masse .

Modélisation de l’amortissement#

L’ordre d’introduction et d’utilisation de la dissipation dans une modélisation discrétisée par éléments finis d’une structure est le suivant :

dissipation intrinsèque liée aux relations de comportement non linéaire irréversible des matériaux constitutifs affectés sur les mailles suports des éléments finis, cf. mot-clé facteur COMPORTEMENT=_F(RELATION='...') cf. [Comportements Non Linéaires]. Ces modèles de comportement peuvent dissiper de l’énergie au cours de cycles en déformation-contraintes (par exemple en élastoplasticité) et/ou en vitesses de déformation-contraintes (par exemple en viscoélasticité non linéaire, en viscoplasticité),
dissipation intrinsèque liée aux relations de comportement non linéaire irréversible des liaisons sur des interfaces (collision, frottement…) : cf. mot-clé facteur COMPORTEMENT avec comme relations : DIS_CHOC, DIS_CONTACT, CHOC_ENDO, CHOC_ENDO_PENA, à l’aide des paramètres AMOR_NOR et AMOR_TAN des matériaux correspondants, cf. [R5.03.17] et [AFFE_CARA_ELEM], qui dissipent de l’énergie au cours de cycles en sauts de déplacement-efforts (par exemple en élastoplasticité) et/ou en vitesses de sauts de déplacement-efforts. Des forces d’amortissement proportionnelles à la vitesse sur ces interfaces sont ajoutées aux équations du mouvement, en fonction de l’état d’activation ou non du contact frottant, sans passer par une contribution à la matrice d’amortissement,
dissipation liée à l’amortissement visqueux linéaire affecté à des éléments discrets, cf. commande [AFFE_CARA_ELEM] (mot-clé facteur DISCRET ou DISCRET_2D, mot-clé CARA = ’A_T_D_N’, etc…) ,
dissipation globale de type amortissement structural (Rayleigh ou modal, dit aussi de Wilson & Penzien, que l’on peut utiliser simultanément),
dissipation liée aux frontières absorbantes (condition de Sommerfeld modélisée avec des éléments paraxiaux d’ordre 0) d’un domaine élastique ou d’un domaine fluide parfait (3D_ABSO, D_PLAN_ABSO, 3D_FLUI_ABSO, FLUI_ABSO, 2D_FLUI_ABSO, AXIS_FLUI_ABSO), cf. [r4.02.05 Éléments de frontière absorbante],
dissipation numérique apporté par le schéma en temps.

Idéalement, les premières catégories devraient être suffisantes, mais en pratique, pour des raisons de simplification de la modélisation, il est souvent indispensable d’ajouter de l’amortissement structural, l’amortissement apporté par le schéma en temps étant le dernier recours.

Nous n’aborderons ici que l’usage de l’amortissement structural et celui lié au schéma d’intégration temporelle (pour plus d’informations sur ce point, le lecteur pourra se reporter à la documentation [U2.06.03]).

La dissipation globale de type amortissement structural est prise en compte par l’amortissement modal ou l’amortissement de Rayleigh.

La résolution transitoire avec DYNA_NON_LINE utilise :

la matrice d’amortissement globale assemblant par sommation toutes les contributions : l’amortissement global de Rayleigh, l’amortissement visqueux linéaire affecté à des éléments discrets, les éléments finis de frontières absorbantes (paraxiaux d’ordre 0) solide et fluide, la matrice constante définiee pour les sous-structures condensées (macroéléments dynamiques), et la matrice que DYNA_NON_LINE assemble à partir d’une matrice élémentaire d’amortissement construite en dehors de l’opérateur pour pouvoir ajouter (cumuler) éventuellement un amortissement d’une autre nature, par exemple IMPE_MECA en formulation U-P d’interaction fluide parfait - structure. C’est l’objet du mot-clé facultatif : MATR_ELEM_AMOR de la fournir à DYNA_NON_LINE. Voir le cas-test [zzzz420a],
le vecteur de forces visqueuses assemblant toutes les contributions : amortissement modal.

Il est possible d’utiliser à la fois l’amortissement global de Rayleigh, ainsi que l’amortissement modal et aussi l’amortissement numérique du schéma en temps. Dans ces cas de combinaisons de plusieurs sources d’amortissement (sans compter la dissipation matériau due à la relation de comportement choisie ou le contact‑frottement), l’utilisateur doit vérifier la cohérence globale de la dissipation et donc de ne pas risquer de compter plusieurs fois l’amortissement sur une même zone. Par exemple, si on souhaite combiner de l’amortissement de global (de Rayleigh ou modal) et de la dissipation matériau, en général les règlements imposent de baisser la valeur d’amortissement global comparativement à un calcul purement élastique.

On notera, cf. [R4.02.05], que l’on ne doit pas employer d’éléments paraxiaux dans le cas MULTI_APPUI.

Autre exemple: le cas d’un bâtiment sur appuis parasismiques. Si on veut représenter la dissipation des appuis par un modèle local (loi non linéaire ou amortisseurs) et l’amortissement du bâtiment par de l’amortissement modal, il faudra vérifier que l’on n’introduit pas d’amortissement modal sur les modes qui font principalement travailler les appuis, comme les modes de corps rigide horizontaux (deux translations et la rotation d’axe vertical).

Rappelons enfin que plus on va multiplier les sources de dissipation, plus leur maîtrise et leur interprétation physique seront ardues.

Amortissement structural#

En dynamique transitoire non linéaire, Code_Aster propose deux types de modélisations de ce type :

l’amortissement modal (de Wilson & Penzien), défini de manière uniforme sur l’ensemble du modèle éléments finis, après avoir établi une base modale réduite sur une bande de fréquence représentative donnée, intervenant sous la forme de forces visqueuses assemblées, linéaires du vecteur vitesse,
l’amortissement global de Rayleigh, défini selon l’affectation des matériaux sur le maillage, intervenant sous la forme d’une matrice assemblée effectuant la relation linéaire entre les forces visqueuses et le vecteur vitesse.

Pour un système linéaire à deux degrés de liberté ces deux modélisations donnent les mêmes forces visqueuses si la calibration de l’amortissement de Rayleigh est effectuée sur les deux fréquences propres. Dans les autres cas, les forces visqueuses diffèrent.

Si on résout les équations du mouvement en repère entraîné, en considérant un mouvement d’entraînement de translation rigide, les forces d’inertie d’entraînement étant définies avec la commande [:ref:CALC_CHAR_SEISME<U4.63.01>] (MONO_APPUI= “OUI”) , alors l’amortissement est construit avec le vecteur vitesse relative. On rappelle qu’en dynamique non linéaire, cf. [r5.05.05 Algorithme non linéaire dynamique], si les mouvements imposés aux appuis ne sont pas un mouvement de corps rigide (cas MULTI_APPUI), alors on doit résoudre les équations du mouvement en repère absolu (galiléen).

Amortissement de Rayleigh#

Cette modélisation utilise une matrice globale d’amortissement visqueux \({\bf C}\) comme étant l’assemblage sur le modèle éléments finis d’une combinaison linéaire des matrices de rigidité et de masse. Cette matrice globale d’amortissement ne se diagonalise sur la base des modes propres naturels que si la calibration des paramètres de Rayleigh est identique dans tout le modèle éléments finis. Les paramètres d’amortissement de Rayleigh sont fournis par le mot-clé facteur ’ELAS’ dans l’opérateur [DEFI_MATERIAU].

(80)#\[{\bf C} = \alpha {\bf K} + \beta {\bf M}\]

Trois cas d’identification sont présentés ici pour illustrer les effets induits par cette modélisation :

amortissement proportionnel aux caractéristiques d’inertie : \(\alpha =0\) \(\beta ={\beta}_{i}\) .

Ce cas a été très utilisé en résolution transitoire en formulation ACCELERATION : si la matrice de masse est diagonale, celle d’amortissement l’est encore et le gain en place mémoire est évident. Le coefficient \(\beta\) peut être identifié à l’amortissement réduit expérimental \({\xi}_{i}\) du mode propre \(({\rho}_{i},{\omega}_{i})\) qui participe le plus à la réponse d’où \({\beta}_{i}=2{\xi}_{i}{\omega}_{i}\) . Pour toute autre pulsation on obtient un amortissement modal réduit \(\xi ={\beta}_{i}\frac{{\omega}_{i}}{\omega}\) . Les réponses sur les modes d’ordre élevé \(\omega \gg {\omega}_{i}\) seront très peu amorties et celles sur les modes basses fréquences \(\omega <{\omega}_{i}\) seront trop amortis.

../../../../_images/100000000000013C0000007DB8575CA3B4D0F32F.png — Fig. 55 Allure de l’amortissement proportionnel à la masse#

Cet amortissement proportionnel à la masse est le seul amortissement de Rayleigh qui soit facilement utilisable avec un schéma d’intégration explicite. Il est illustré sur la figure Allure de l’amortissement proportionnel à la masse. En effet, on peut montrer que l’introduction d’un terme proportionnel à la raideur entraîne une baisse du pas de temps critique (condition CFL).

amortissement proportionnel aux caractéristiques de rigidité: \(\alpha ={\alpha}_{j}\) \(\beta =0\) .

Le coefficient peut être identifié, comme précédemment à partir de l’amortissement modal \({\xi}_{j}\) associé au mode \(({\rho}_{j},{\omega}_{j})\) , d’où \({\alpha}_{j}=2\frac{{\xi}_{j}}{{\omega}_{j}}\) . Pour toute autre pulsation on obtient un amortissement modal réduit \(\xi ={\alpha}_{j}\frac{\omega}{{\omega}_{j}}\) . Les réponses sur les modes élevés \(\omega \gg {\omega}_{j}\) sont donc très amorties. On l’illustre sur la figure Allure de l’amortissement proportionnel à la raideur

../../../../_images/10000000000001210000007AD759DF3D0B349A54.png — Fig. 56 Allure de l’amortissement proportionnel à la raideur#

amortissement proportionnel complet : \(\alpha ={\alpha}_{j}\) \(\beta ={\beta}_{i}\).

A partir d’une identification sur deux fréquences propres différentes \(({\rho}_{i},{\omega}_{i})\) et \(({\rho}_{j},{\omega}_{j})\) , nous obtiendrons pour toute autre pulsation un amortissement modal réduit \(\xi ={\alpha}_{j}\frac{\omega}{{\omega}_{j}}+{\beta}_{i}\frac{{\omega}_{i}}{\omega}\) . Dans l’intervalle \([{\omega}_{i},{\omega}_{j}]\) la variation de l’amortissement réduit est faible et en dehors on retrouve la combinaison des inconvénients précédents : les réponses sur les modes extérieurs à l’intervalle sont trop amorties.

../../../../_images/100000000000013D000000775C7AF42F4EDB560B.png — Fig. 57 Allure de l’amortissement de Rayleigh complet#

L’amortissement de Rayleigh complet illustré sur la figure Allure de l’amortissement de Rayleigh complet permet d’avoir une valeur d’amortissement presque constante sur un plateau de fréquence donné, toujours supérieure à \(\sqrt{\alpha \beta}\), ce qui permet de contrôler son effet sur une plage fréquentielle définie en cohérence avec le problème considéré. Cependant, lorsque la raideur s’endommage, les réponses dynamiques transitoires présentent une fréquence centrale qui baisse. L’amortissement de Rayleigh perçu est alors légèrement plus fort en basses fréquences et plus faible en hautes fréquences.

Application à la structure

Les coefficients d’amortissement visqueux global de Rayleigh sont définis, au niveau des caractéristiques du matériau (commande [DEFI_MATERIAU]), par les paramètres AMOR_ALPHA et AMOR_BETA dans le mot-clé facteur ELAS. Il n’y a pas de mot-clé spécifique dans DYNA_NON_LINE pour activer la prise en compte de cet amortissement: si on ne veut pas l’utiliser, il faut retirer toutes les occurrences des mot-clés AMOR_ALPHA et AMOR_BETA de tous les matériaux du modèle éléments finis étudié.

Les valeurs à imposer pour obtenir l’amortissement souhaité \(\xi\) dans l’intervalle des fréquences propres \({f}_{1}\) et \({f}_{2}\) se déduisent des équations suivantes:

(81)#\[\alpha =\frac{\xi}{\pi ({f}_{1}+{f}_{2})}\]

(82)#\[\beta =\frac{4\pi \xi {f}_{1}{f}_{2}}{{f}_{1}+{f}_{2}}\]

où \({f}_{1}\) et \({f}_{2}\) sont les deux fréquences propres bornant l’intervalle d’étude considéré. Dans le cadre de ce document, on cherche des solutions basses fréquences, donc les fréquences \({f}_{1}\) et \({f}_{2}\) seront associées aux premières fréquences du modèle, dont les modes sont cohérents avec le chargement imposé.

Pour donner des ordres de grandeur, l’amortissement modal pour les structures en acier et généralement de l’ordre de quelques %, alors que pour des structures en béton, de type génie civil, on peut monter jusqu’à 5, voire 7 % d’amortissement, pour des calculs globalement linéaires, censés représenter des dissipations d’énergie provenant de dégradations ou fissurations diffuses.

Pour les éléments discrets, les paramètres d’amortissement de Rayleigh n’étant pas définissables dans la commande [AFFE_CARA_ELEM], si on utilise l’opérateur DYNA_NON_LINE, l’amortissement de Rayleigh ne tient donc pas compte de la contribution des éléments discrets. Par contre, la matrice d’amortissement globale assemblée pour la résolution tiendra compte de leur contribution et de celle d’amortissement de Rayleigh provenant des éléments finis (volumiques, surfaciques ou de type poutre…).

Dans le cas des modèles couplés en hydromécanique THM, pour l’analyse d’ouvrages géomécaniques (barrages en terre, colonnes de sol…), cf. [u2.04.05 Notice d’utilisation du modèle THM], DYNA_NON_LINE ne peut construire directement la matrice d’amortissement de Rayleigh. Voir [u2.04.05 Notice d’utilisation du modèle THM] et [U2.04.08]. Un exemple de calcul dynamique non linéaire sur un problème hydromécanique d’une colonne de sol est fourni par le cas-test [wdnp101a]. Cette difficulté peut être contournée en superposant un modèle mécanique classique sur le même maillage permettant de représenter l’énergie cinétique (matrice de masse) et l’énergie élastique (matrice de raideur), les deux modèles s’appuyant sur les mêmes nœuds. Puis de construire la matrice élémentaire d’amortissement avec les paramètres \(\alpha ={\alpha}_{j}\) \(\beta ={\beta}_{i}\) désirés, en utilisant la commande [CALC_MATR_ELEM] avec OPTION=”AMOR_MECA” en donnant les matrices élémentaires de masse et de rigidité (MASS_MECA et RIGI_MECA), puis en spécifiant le concept obtenu dans DYNA_NON_LINE avec le mot-clé MATR_ELEM_AMOR.

La commande DYNA_NON_LINE propose une variante de l’amortissement de Rayleigh utilisant la matrice de raideur tangente au lieu de la matrice de raideur initiale élastique. C’est l’objet du mot-clé facultatif : AMOR_RAYL_RIGI=("TANGENTE"/ "ELASTIQUE"). La valeur par défaut étant ELASTIQUE. On aura alors une matrice d’amortissement de Rayleigh fixe pendant tout le calcul, résultant des valeurs choisies de caractéristiques des matériaux (paramètres AMOR_ALPHA et AMOR_BETA). Bien entendu, la matrice d’amortissement assemblée totale peut être malgré tout variable en présence d’autres sources d’amortissement, comme ce qui provient des relations sur des interfaces : DIS_CHOC, DIS_CONTACT, CHOC_ENDO, CHOC_ENDO_PENA…

Avec la valeur TANGENTE du mot-clé AMOR_RAYL_RIGI, la matrice de raideur réactualisée sera la même que celle qui est utilisée pour le calcul des efforts internes. Le choix de l’amortissement de Rayleigh avec matrice TANGENTE (construite avec les paramètres AMOR_ALPHA et AMOR_BETA des matériaux) introduit une non-linéarité sur les forces d’amortissement et cela peut engendrer des découpes de pas de temps supplémentaires par rapport au calcul de l’amortissement avec la matrice de raideur ELASTIQUE. En effet, les perturbations hautes fréquences parasites sont alors moins amorties qu’avec cette dernière. Si de plus, le matériau possède un comportement non linéaire irréversible avec endommagement adoucissant, alors la matrice de raideur TANGENTE peut avoir des valeurs propres négatives, donc produire, si on la choisit pour construire la matrice d’amortissement de Rayleigh, des instabilités de flottement en hautes fréquences parasites.

Si le modèle contient des variables de commandes (température, irradiation, etc.) alors la matrice de raideur élastique n’est pas constante au cours du temps, on doit alors adopter la matrice de raideur réactualisée TANGENTE pour construire l’amortissement de Rayleigh.

Pour les éléments de joints, on note que la matrice de raideur élastique n’existe pas, seule la matrice de raideur réactualisée TANGENTE est disponible.

Amortissement dû au schéma en temps#

La documentation [R5.05.05] et surtout la note [bib6] présentent cet aspect. On va ici se borner à en rappeler les grandes tendances.

Sur un système à un degré de liberté linéaire (masse ressort, de pulsation propre \(\omega\) ), on peut obtenir la caractérisation suivante de l’amortissement induit par le schéma implicite (accélération moyenne, accélération moyenne modifiée et HHT complet), en fonction du pas de temps et pour différentes valeurs du paramètre ALPHA. On en a une illustration sur la figure Comparaison de l’amortissement dû au schéma en temps.

../../../../_images/101531D000004348000030FDDA3833C6CE87FA8E1.svg — Fig. 58 Comparaison de l’amortissement dû au schéma en temps#

On retrouve bien que seul le schéma d’accélération moyenne ne dissipe pas.

Quand on compare les deux autres schémas, on peut remarquer que:

seul le schéma HHT complet ne perturbe pas le domaine basse fréquence,
pour une même valeur du paramètre ALPHA l’accélération moyenne modifiée introduit beaucoup plus de dissipation que le schéma HHT.

Enfin, il convient de remarquer que la valeur d’amortissement équivalente est dépendante de la pulsation \(\omega\) , et donc va dépendre de l’élément fini considéré. Sur un problème complexe, l’amortissement dû au schéma ne sera donc pas homogène: plus l’élément sera «raide», plus il verra d’amortissement.

De même, si l’on diminue le pas de temps, l’amortissement va baisser.

Afin de mettre en exergue l’influence de l’amortissement haute fréquence des schémas implicites, on va présenter quelques évolutions de l’accélération solution d’un problème linéaire simple de tuyauterie sous séisme.

Sur le graphe Comparaison locale du comportement des schémas en temps, on a zoomé sur une partie de la réponse en accélération pour comparer différentes méthodes de résolution du problème transitoire. La solution de référence en pointillés verts est obtenue par calcul sur base modale (DYNA_VIBRA). La troncature de la base modale filtre naturellement toute perturbation haute fréquence.

../../../../_images/101C9FEA00005F160000455A691D38C3AAB71180.svg — Fig. 59 Comparaison locale du comportement des schémas en temps#

On compare cette solution de référence à un calcul avec le schéma d’accélération moyenne modifié (courbe rouge) qui oscille fortement malgré l’ajustement du paramètre ALPHA.

Enfin, on trace la réponse obtenue avec un schéma HHT complet (courbe noire) qui donne un résultat très proche de la solution de référence: les perturbations hautes fréquences (liées au pas de temps) sont fortement atténuées.

Sur le graphe Comparaison globale du comportement des schémas en temps, on observe les courbes de réponse sur un intervalle de temps plus grand. Cela permet de saisir l’influence du schéma sur la réponse «physique»: donc basse ou moyenne fréquence. Les différentes réponses correspondent à différents schémas d’intégration et valeurs du pas de temps. La solution de référence (le problème étant linéaire) est obtenue par superposition modale (courbe en traits mixtes épais marrons, nommée « modal coupure \(100\mathit{Hz}\) » et obtenue avec DYNA_VIBRA).

../../../../_images/101E016200005CB5000044F0C052965149D713CC.svg — Fig. 60 Comparaison globale du comportement des schémas en temps#

On constate que:

la solution non dissipative (HHT avec ALPHA=0, donc un schéma de type accélération moyenne) avec un pas de temps « grand » surévalue l’amplitude et introduit un déphasage,
les solutions obtenues avec des schémas de type accélération moyenne modifiée (MODI_EQUI='NON'), quel que soit le pas de temps, ont une trop forte dissipation,
les solutions obtenues avec le schéma HHT complet (MODI_EQUI='OUI') avec un pas de temps fin, permettent de bien retrouver la solution de référence.

Pour conclure sur cette partie, étant données les allures de l’amortissement de Rayleigh, ainsi que celui dû au schéma, on peut construire un amortissement total relativement varié.

Pour choisir de manière optimale la valeur du paramètre ALPHA du schéma HHT, il est indispensable de mener une étude paramétrique, afin de trouver la valeur la plus proche de zéro adaptée au problème considéré (on peut conseiller de commencer avec une valeur de -0,1, par exemple). Pour valider ce choix, on peut se baser sur les deux types d’analyses présentées ci-dessus: maîtrise des oscillations en hautes fréquences sur les accélérations ( cf. figure ) et absence de sur‑amortissement numérique sur la réponse en déplacement ( cf. figure ). Sur ce dernier aspect, on peut comparer à la réponse en déplacement obtenue avec un schéma non dissipatif de type NEWMARK.

Avec le schéma HHT complet, l’amortissement à basse et moyenne fréquence étant faible, il convient de ne pas négliger d’introduire dans le modèle une modélisation ad hoc de l’amortissement réel (global et/ou dû aux non‑linéarités).

L’allure de l’amortissement apporté par le schéma de Tchamwa est qualitativement proche de celle de l’accélération moyenne modifiée.

Modification de la matrice de masse#

Dans certains cas spécifiques, comme les calculs sur des modèles mixtes hydromécaniques (HM ou THM, cf. [u2.04.05 Notice d’utilisation du modèle THM]) ou fluide-structure, cf. [r4.02.02 Éléments pour le couplage interaction fluide-structure linéaire avec fluide inerte] et [U2.04.08], les cas ayant un très mauvais conditionnement de la matrice de masse, il peut être utile de modifier la matrice de masse afin de la rendre inversible. C’est aussi utile pour établir une initialisation de l’accélération au démarrage du schéma d’intégration temporelle qui réduise l’overshooting. C’est même indispensable, si ce n’est pas le cas, pour les schémas en temps explicites.

Cette modification de la matrice de masse peut se faire avec le mot-clé COEF_MASS_SHIFT ( cf. documentation [R5.05.05]) dans le mot-clé SCHEMA_TEMPS. Cela permet d’introduire un « shift » (ou décalage) de la matrice de masse \(M\) qui devient: \({\bf M} '={\bf M} +\mathrm{coef.} {\bf K}\). Ce coefficient positif mathrm{coef.} a pour dimensions l’inverse d’une pulsation au carré. Ce décalage de la matrice de masse ne modifie pas les modes.

Cette fonctionnalité permet en outre d’améliorer la convergence des calculs implicites ou d’augmenter la valeur du pas de temps critique en explicite. Il faut cependant signaler que cela revient à perturber le système mécanique, en particulier en décalant les fréquences propres du système vers le bas (et en diminuant la célérité des ondes) : l’erreur dépasse 5% sur les fréquences propres dès que celles-ci sont au-dessus de \((\sqrt{\mathrm{coef.}})^{-1}\). Par exemple, une pulsation propre d’origine \(\omega\) devient alors \(\omega '\) telle que : \(\omega {'}^{2}=\frac{{\omega}^{2}}{1+\mathit{coef}.{\omega}^{2}}\).

On constate d’après cette expression que \(\omega '\) tend vers une valeur maximale \({\mathit{coef}}^{-0.5}\) : le système modifié introduit donc une fréquence de coupure pilotée par la valeur du coefficient de « shift ». Ainsi, par exemple, pour une valeur maximale en pratique pour \(\mathrm{coef}={10}^{-6}\) , la pulsation propre maximale corrigée du système \(\omega '\) vaut \(1000r{\mathrm{d.s}}^{-1}\) , et une valeur pour une fréquence propre d’origine de \(30\mathrm{Hz}\) est corrigée par cette méthode à \(29,481\mathit{Hz}\) . Les fréquences propres situées sous la nouvelle fréquence de coupure sont donc légèrement altérées.

Dans la pratique, il convient de ne pas utiliser de valeur de « shift » trop grande, sauf à se limiter au calcul de solutions à très basse fréquence, voire de type quasistatique. Dans ce dernier cas, l’algorithme transitoire non linéaire peut alors être vu comme un solveur régularisé pour le quasistatique, ce qui peut être utile lorsque la résolution classique statique ne converge pas.

Dans le même ordre d’idée d’utiliser un solveur transitoire pour obtenir une réponse statique dont la convergence avec une méthode implicite habituelle est très difficile, on peut aussi augmenter fortement la masse du système (par exemple en jouant sur la masse volumique) et utiliser un schéma en temps explicite avec un fort amortissement structurel. Dans ce cas on combine la capacité de convergence de la méthode explicite avec un pas de temps critique de stabilité pas trop faible, du fait de l’augmentation de la masse (ce pas de temps de stabilité augmentera comme la racine carrée du coefficient d’augmentation de la masse totale: si on multiplie la masse totale par 100, le pas de temps limite sera multiplié par 10 car la célérité des ondes sera dix fois plus faible). Ce type de méthode a déjà été utilisé pour des cas délicats comme la simulation d’excavation de galeries souterraines, où le sol environnant peut devenir instable.

Instabilité et analyse modale réactualisée#

Au cours de la résolution du problème transitoire il est possible d’utiliser des outils d’analyse aux valeurs propres sur les opérateurs globaux réactualisés. On peut mener deux types d’analyse.

D’une part, le calcul des fréquences propres et modes vibratoires avec la matrice de raideur réactualisée. Cela correspond au mot-clé facteur MODE_VIBR (cette option ayant besoin implicitement de la matrice de masse, elle n’est pas disponible dans [STAT_NON_LINE]). La matrice de raideur peut alors être la raideur élastique, sécante ou tangente à choisir avec le mot-clé MATR_RIGI. Grâce à ce mot-clé facteur on peut suivre l’influence des non-linéarités sur le comportement vibratoire d’une structure. Un exemple d’application serait le cas de structures en béton armé pour lesquelles l’endommagement fait réduire les fréquences propres (avec la raideur tangente ou sécante). Avec une structure élastoplastique, en ayant choisi pour ce faire la raideur tangente, on observe une chute des fréquences propres lors des régimes comportant une incursion plastique, tandis que l’on retrouve les fréquences propres initiales lors des régimes élastiques.

Le graphe Représentation schématique des opérateurs de raideur présente les choix possibles pour la matrice de raideur (pour un matériau endommageant):

../../../../_images/100021B0000069F0000040B3A4BE5E1BF02D329E1.svg — Fig. 61 Représentation schématique des opérateurs de raideur#

D’autre part, on peut, grâce au mot-clé facteur CRIT_STAB=_F(TYPE='FLAMBEMENT'), mener une analyse de stabilité de l’opérateur de raideur. Dans le cas des petites perturbations et où l’on sait calculer la raideur géométrique, alors cette option s’assimile à une analyse de flambement au sens d’Euler sur la matrice de raideur actualisée. Dans les autres cas, quand on ne sait pas calculer la raideur géométrique, alors on bascule sur la recherche de singularité de l’opérateur de raideur seul. Dans tous les cas on obtient des valeurs propres qui vont évoluer au cours du calcul.

Dans le cas du flambement d’Euler, la valeur propre est directement le coefficient multiplicateur du chargement qui permet d’obtenir la charge critique.

Dans le deuxième cas, l’interprétation est moins aisée (les valeurs propres ne sont pas adimensionnelles). Si l’on constate qu’une valeur propre change de signe, cela signifie que la solution calculée a passé une bifurcation et donc que l’on a perdu l’unicité de la solution.

Dans tous les cas, il s’agit d’analyse de stabilité au sens « statique » et, d’ailleurs, cette option est disponible dans STAT_NON_LINE. Actuellement, il n’existe pas d’opérateur permettant de mener une analyse de stabilité au sens dynamique (flottement, divergence) : donc, par exemple, en calculant l’amortissement du système pour détecter quand il devient négatif.

Comme ces opérations demandent un certain coût CPU (comparable à un appel à CALC_MODES sur les quelques premières fréquences, avec le mot-clé NMAX_FREQ, à chaque instant demandé), on a introduit la possibilité de ne calculer ces valeurs propres que pour la liste d’instants d’archivage si elle existe. Pour réduire encore le temps CPU, il est aussi possible de faire plusieurs poursuites et de ne demander le calcul des valeurs propres que sur certains intervalles de temps ou instants particuliers.

On peut aussi récupérer les modes propres (respectivement de flambage) associés à l’exécution du mot-clé MODE_VIBR. Le concept resultat produit par DYNA_NON_LINE est alors enrichi d’une table_container. qui porte le nom 'ANALYSE_MODALE'.

Pour la récupérer dans le concept résultat produit par DYNA_NON_LINE, on fait appel à la commande [RECU_TABLE]. Dans cette table, on pourra récupérer les valeurs des paramètres : instant de calcul considéré, numéro de mode, fréquence (respectivement charge critique).

On pourra ensuite imprimer avec [IMPR_TABLE] cette table, en filtrant selon ces paramètres.

En utilisant la commande [EXTR_TABLE], en spécifiant un TYPE_RESU='MODE_MECA', en filtrant selon les paramètres voulus, on obtient un concept résultat, dont on a choisi le nom, contenant les champs de déplacements des modes calculés, qu’on pourra ensuite visualiser. On pourra consulter le cas-test sdnv106 [V5.03.06].

Archivage et post-traitement#

Le nombre de pas de temps pouvant être très grand, il est fortement recommandé d’utiliser les fonctionnalités d’archivage sur une liste d’instant plus grossière (mot-clé facteur ARCHIVAGE de DYNA_NON_LINE) sous peine d’avoir des bases et fichiers de sorties énormes. Le pas d’archivage peut aller de quelques pas en schéma implicite à 10 à 100 pas en schéma explicite. Attention cependant au risque de rater d’éventuelles valeurs extrêmes dans l’analyse des grandeurs d’intérêt de la solution calculée.

En complément, si l’on a besoin de suivre précisément au cours du temps l’évolution de quelques grandeurs d’intérêt sur quelques GROUP_NO, il existe l’observation (mot-clé facteur OBSERVATION de DYNA_NON_LINE) qui vient compléter l’archivage.

Comme on a pu le voir au paragraphe § Amortissement dû au schéma en temps , si l’on veut analyser les réponses en vitesse ou accélérations au cours du temps, on peut obtenir des courbes assez chahutées. Ces oscillations hautes fréquences peuvent être le signe d’une discrétisation en temps insuffisante (ou d’une irrégularité en temps trop grande) du problème. Il est aussi possible, en utilisant un schéma en temps dissipatif type HHT complet de lisser ces perturbations. Un compromis reste à trouver entre ce lissage et une trop forte dissipation de la réponse. D’une manière générale, il faut bien intégrer que les valeurs instantanées des quantités les moins lissées comme l’accélération sont à manipuler avec précaution. Il vaut mieux chercher à mener son analyse sur des quantités intégrées plus physiquement pertinentes en dynamique comme l’énergie.

En complément, toutes les méthodes d’analyse venant de la quasistatique pour quantifier la qualité d’un résultat (dont les différentes normes résidu en équilibre) sont disponibles et pertinentes avec DYNA_NON_LINE.

Enfin, on peut aussi rappeler que si l’on souhaite mener des analyses fréquentielles, il convient de bien veiller à respecter le domaine de validité de la FFT, par exemple la causalité du signal à traiter. Pour cela, ce signal en temps doit partir et finir à zéro, avec des dérivées nulles. Sans le respect de ces hypothèses de base, l’utilisateur risque d’obtenir des résultats fréquentiels imprécis.

Bilan énergétique#

Pour toute analyse transitoire, il est très utile de pouvoir disposer du bilan énergétique global, sur la durée du transitoire, ce qui permet d’analyser la réponse du système, en particulier la contribution de la dissipation d’énergie par amortissement, et aussi de contrôler la qualité de la solution numérique. La documentation [R4.09.01] présente toutes les fonctionnalités énergétiques disponibles avec le mot-clé facteur ENERGIE=_F(CALCUL='OUI') de l’opérateur DYNA_NON_LINE.

On utilisera pour cela le mot-clé facteur ENERGIE=_F(). Il produit la table de nom PARA_CALC. Le bilan d’énergie peut être extrait de cette table à l’aide de la commande [RECU_TABLE] puis imprimé avec [IMPR_TABLE].

On trouvera dans la table de bilan des énergies au cours du transitoire les grandeurs cumulées à chaque instant :

ENER_CIN, pour l’énergie cinétique, toujours positive ou nulle,
ENER_TOT, pour le travail cumulé de déformation (somme de l’énergie stockée et de l’énergie dissipée par les modèles de comportement irréversibles, incluant le cas des frontières absorbantes),
TRAV_AMOR, pour le travail cumulé des forces d’amortissement issues de la matrice d’amortissement assemblée (Rayleigh, éléments discrets) ou d’amortissement modal,
TRAV_LIAI, pour le travail cumulé incluant le cas des forces internes des matériaux d’interface (DIS_CONTACT et DIS_CHOC, CHOC_ENDO, CHOC_ENDO_PENA`, JOINT_MECA_FROT), cf. R4.09.01, § 1.3.2,
TRAV_EXT, pour le travail cumulé incluant les forces d’excitation au second membre (et des forces d’inertie d’entraînement).

Le reliquat correspond au travail dissipé par le schéma d’intégration numérique.

Fonctionnalités disponibles dans STAT_NON_LINE et non dans DYNA_NON_LINE#

Les méthodes de type recherche linéaire (mixte ou pas) ne sont pas autorisées en dynamique. Cependant, ce manque est à relativiser, sachant que les tentatives d’applications de ces méthodes sur des études de structures en béton armé en dynamique n’ont pas mis en avant d’apport significatif sur la convergence, contrairement à ce que l’on observe en quasistatique. Signalons, néanmoins, qu’aucun argument théorique n’interdirait l’usage de ces méthodes en dynamique.

Enfin, les techniques de pilotage disponibles dans STAT_NON_LINE (longueur d’arc, par exemple) sont interdites en dynamique car elles n’ont alors pas de sens : le temps a, en dynamique, un sens physique.

Spécificités des problèmes fluide-structure couplés#

Il est possible d’utiliser un modèle vibro-acoustique couplé dans DYNA_NON_LINE. Ce modèle se base sur une approche \((u,p,\varphi )\) avec les hypothèses suivantes :

le fluide est de type acoustique linéaire,
la structure doit être considérée comme étant en petites perturbations ou en lagrangien réactualisé.

On peut aussi prendre en compte les surfaces libres soumises à la pesanteur (ondes de ballottement).

La documentation [U2.06.11] présente en détail la mise en œuvre d’un modèle fluide-structure non linéaire couplé pour un calcul de réservoir. La note [bib5] analyse le domaine d’application du modèle fluide-structure couplé.

Optimisation des performances#

En quasistatique, il n’est pas rare de devoir effectuer plus de 10 itérations pour avoir la convergence au sens du résidu en équilibre. En dynamique implicite cette valeur de 10 itérations constitue, en général, une bonne valeur de départ pour le paramètre ITER_GLOB_MAXI du mot-clé facteur CONVERGENCE dans DYNA_NON_LINE. Si l’on ne peut converger en moins de 10 à 20 itérations, il est alors préférable de diminuer le pas de temps plutôt que d’augmenter le nombre maximal d’itérations autorisé.

En explicite, il n’y a pas d’itérations pour l’équilibre, le coût de calcul de chaque pas de temps sera donc constant, quel que soit le niveau de non-linéarité (hormis, éventuellement, la vérification locale du comportement).

L’utilisation, même courante, des schémas d’intégration explicites semble donc très séduisante au vu du temps CPU qui reste maîtrisé. Il faut cependant tempérer cet optimisme en gardant bien à l’esprit que l’on se prive du garde-fou qu’est la vérification précise de l’équilibre et que, par conséquent, la qualité de la solution explicite obtenue doit être analysée avec plus de précautions. L’algorithme explicite ne divergera pas (si l’on respecte la condition de stabilité de Courant), mais la solution obtenue n’est pas garantie par un critère de vérification de l’équilibre. En particulier une étude paramétrique sur le pas de temps est indispensable car l’allure de la solution peut varier fortement au fur et à mesure du transitoire.

De plus, code_aster n’est pas un code optimisé pour les calculs explicites et ses performances en explicite sont modestes, comparées aux codes spécialisés comme Europlexus [bib2].

Réduction de modèle et condensation dynamique par sous-structuration#

Une solution pour diminuer le temps de calcul est de projeter le problème sur une base réduite (base modale ou base de Ritz). On diminue alors grandement le nombre de degrés de liberté et ce type d’approche est disponible dans DYNA_NON_LINE (la résolution gagne aussi à être explicite car la condition de stabilité de Courant sur base modale réduite est peu pénalisante, puisque cette base possède une fréquence de coupure basse). Pour résumer, ce type d’approche est particulièrement adaptée aux problèmes où les non-linéarités restent modérées et localisées. Dès que les non-linéarités deviennent fortes, on peut se poser la question de la réactualisation de la base réduite initiale qui perd de sa cohérence avec la solution courante. Le surcoût de calcul dû au recalcul de la base et aux reprojections vient alors diminuer l’intérêt de ce type de méthode. Le mot-clé facteur PROJ_MODAL permet de faire ce type de calcul avec un schéma d’intégration explicite, les degrés de liberté étant généralisés. Un exemple de calcul est fourni par le cas-test SDNV107a [V5.03.107]. Un autre exemple de calcul dynamique non linéaire sur un problème hydromécanique d’une colonne de sol est fourni par le cas-test WDNP101a [V5.03.107].

On préconise aussi d’utiliser cette démarche réduction /condensation de manière « mixte », pour des problèmes « basses fréquences », où le modèle est constitué de plusieurs sous-domaines en condensant la résolution des sous-domaines de comportement linéaire sur les degrés de liberté de leurs interfaces avec des sous-domaines de comportement non linéaire, traités eux en éléments finis. Chacun de ces sous-domaines linéaires est condensé sur un macro-élément dynamique MACR_ELEM_DYNA [U4.65.01] défini sur chaque interface. Au moyen de l’opérateur - [U4.63.33] après résolution du problème global « mixte » on peut reconstituer les champs dans chaque sous-domaine. Un exemple de calcul dynamique non linéaire par réduction /condensation est fourni par le cas-test FORMA12g [V2.08.012] avec un schéma d’intégration implicite. La méthodologie est présentée en détail dans [U2.07.04].

Parallélisme#

Le parallélisme peut apporter un gain d’autant plus notable que le système à résoudre comportera un grand nombre de DDL et que le nombre de pas de temps sera faible. En effet, l’algorithme de résolution impose des communications importantes à chaque pas (on n’est pas dans une stratégie de type multi-domaine en temps et en espace). En pratique, sur des problèmes comportant quelques centaines de milliers de DDL, le speed-up reste bon jusque 16 à 32 processeurs. Dans la plupart des cas, le solveur itératif PETSC apportera un gain comparativement au solveur MUMPS. En plus du parallélisme au niveau solveur, le parallélisme des étapes élémentaires (résolution du comportement aux points de Gauss) apportera un gain supplémentaire important si l’intégration comportement est coûteuse. Dans tous les cas, il est très utile de consulter dans le fichier de message .mess les mesures de temps CPU par étape, ce qui permet d’identifier les étapes les plus coûteuse de la résolution transitoire pour adapter le parallélisme en conséquence. La documentation [U2.08.06] donne tous les détails pour l’utilisation du parallélisme.

Conclusion#

Ce document présente quelques règles générales pour faciliter l’utilisation de méthodes dynamiques transitoires pour la simulation de systèmes non linéaires : voir [bib3], [bib6], [bib8].

La première étape est l’adaptation du modèle aux méthodes d’analyse dynamique. Il s’agit principalement de s’assurer de la bonne régularité des conditions imposées, du bon échantillonnage des excitations, de la définition correcte de la masse volumique et de l’amortissement global (Rayleigh, modal) ; il est recommandé de commencer par une analyse transitoire linéaire.

Ensuite, il est recommandé de commencer par utiliser un schéma implicite (DYNA_NON_LINE avec un schéma en temps de type NEWMARK pour les problèmes relativement réguliers ou HHT, et THETA_SCHEMA pour les problèmes avec collisions). En effet, les schéma en temps implicites sont les plus performants et les plus généralistes dans Code_Aster en dynamique dite « lente ».

Enfin, pour certaines applications comme la dynamique rapide, le calcul sur base modale ou certains cas de calculs en évolution lente ( cf. [U2.04.07]), l’utilisateur a la possibilité d’utiliser des schémas en temps explicites. Les performances en temps CPU de Code_Aster en explicite sont assez faibles, si l’on compare à un code dédié comme Europlexus [bib2]. De plus, toutes les fonctionnalités disponibles en implicite ne le sont pas en explicite (comme, par exemple, pour le contact où seule la pénalisation est autorisée).

Afin de qualifier la qualité de la solution numérique obtenue, il est indispensable de mener certaines études paramétriques :

comme pour les calculs quasistatiques, en jouant sur la discrétisation spatiale,
en testant différents paramétrages de l’amortissement,
en testant différents pas de temps,
en testant différents schémas en temps, d’abord non dissipatifs et si besoin avec dissipation numérique en hautes fréquences.

Pour valider ce choix de paramètres, il est pertinent d’analyser les réponses en déplacements et en accélérations, ainsi que de se servir des fonctionnalités de bilan d’énergie.

Bibliographie#

[bib1]

Courant, K. Friedrichs & H. Lewy, Über die partiellen Differenzengleichungen der mathematischen Physik, Mathematische Annalen, Vol. 100, No. 1, pp. 32–74, 1928.

[bib2] (1,2,3,4)

Europlexus. Code de calcul par éléments finis en dynamique rapide. Manuel de Référence. CCR/CEA/EDF R&D ( http://europlexus.jrc.it/).

[bib3]

Géradin, M. Hogge & G. Robert, Time Integration of Linear and Nonlinear Dynamic Problems, Finite Element Handbook, Part 4, 1, 4, ED. W. Pilkey, Mc Graw-Hill, 1987.

[bib4]

Greffet, Simulation couplée fluide-structure appliquée aux problèmes d’instabilité non linéaire sous écoulement, Thèse de doctorat, LMT, ENS-Cachan, 2001.

[bib5]

Greffet, Voies d’amélioration de la formulation couplée fluide-structure dans Code_Aster, Note EDF R&D HT-62/02/023/A, 2002.

[bib6] (1,2)

Greffet, Vers de nouvelles méthodes numériques pour l’intégration temporelle dans le Code_Aster. Note EDF R&D HT-62/04/016/A, 2004.

[bib7]

Greffet, Evaluation des méthodes transitoires pour les calculs d’excavation, Note EDF R&D H-T62-2007-02878-FR, 2008.

[bib8]

H.M. Hilber, T.J.R. Hughes & R.M. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthquake Eng. Struct. Dyn., Vol. 5, 283-292, 1977.

[bib9]

T.J.R. Hughes & T. Belytschko, Nonlinear finite element analysis, Zace Services ltd – ICE Division, 2000.

[bib10]

M.N. Newmark, A Method of Computation for Structural Dynamics, Proc. ASCE 85, EM3, 1959.

[bib11]

Tchamwa, C. Wielgosz, Une nouvelle méthode explicite d’intégration directe précise et à dissipation numérique contrôlable, Actes du 13e Congrès Français de Mécanique, Vol. 1, Poitiers, pp. 251-254, septembre 1997.

u2.06.13 Conseils généraux d’utilisation de l’opérateur DYNA_NON_LINE#

Définition d’un problème adapté à la dynamique transitoire#

Modèle#

Conditions aux limites et excitations#

Mono-appui#

Multi-appuis#

Contact et collisions#

Conditions initiales#

Discrétisation du problème continu#

Maillage#

Schéma d’intégration en temps#

Schéma d’intégration explicite#

Schéma d’intégration implicite#

Synthèse pour le choix du schéma d’intégration#

Conditions initiales#

Modélisation de l’amortissement#

Amortissement structural#

Amortissement modal#

Amortissement de Rayleigh#

Amortissement dû au schéma en temps#

Modification de la matrice de masse#

Instabilité et analyse modale réactualisée#

Archivage et post-traitement#

Bilan énergétique#

Fonctionnalités disponibles dans STAT_NON_LINE et non dans DYNA_NON_LINE#

Spécificités des problèmes fluide-structure couplés#

Optimisation des performances#

Réduction de modèle et condensation dynamique par sous-structuration#

Parallélisme#

Conclusion#

Bibliographie#