r5.05.05 Algorithme non linéaire dynamique#
Résumé:
L’opérateur DYNA_NON_LINE [DYNA_NON_LINE] s’emploie pour l’analyse dynamique non linéaire des structures par une intégration directe en temps. Les non-linéarités peuvent provenir du comportement du matériau, des liaisons (contact-frottement), ou de grandes transformations géométriques (grands déplacements et grandes rotations).
L’organisation de DYNA_NON_LINE s’apparente fortement à celle de l’opérateur quasi-statique non linéaire STAT_NON_LINE [r5.03.01]. A priori , toutes les relations de comportement développées dans le cadre de STAT_NON_LINE fonctionnent dans celui de DYNA_NON_LINE.
On présente ici la formulation générale du problème dynamique non linéaire afin de préciser les articulations entre les aspects purement dynamiques et ceux déjà traités dans d’autres opérateurs ou formulations disponibles dans code_aster : gestion des conditions aux limites, des couplages fluide-structure, de l’amortissement, du calcul en repère relatif, puis les propriétés du schéma d’intégration numérique temporel, qui s’élabore indépendamment de toute relation de comportement. On expose comment celui-ci s’articule avec l’algorithme de Newton pour traiter les non-linéarités matérielles et géométriques.
On dispose de trois schémas en temps implicites complémentaires et efficaces:
la famille de Newmark,
l’accélération moyenne modifiée («HHT» dans DYNA_NON_LINE, avec l’option MODI_EQUI=’NON’),
le schéma HHT complet («HHT» dans DYNA_NON_LINE, avec l’option MODI_EQUI=’OUI’), ainsi que sa version « sans overshooting » proposée dans [bib35] (“NOHHT” dans DYNA_NON_LINE).
L’opérateur propose aussi deux schémas explicites:
les différences centrées,
le schéma dissipatif de Tchamwa-Wielgosz.
On donne quelques conseils pour un bon usage que la documentation [Conseils généraux d’utilisation de l’opérateur DYNA_NON_LINE] vient compléter.
Dynamique non linéaire : discrétisation spatiale du problème continu#
Résoudre un problème de dynamique non linéaire nécessite de décrire tout d’abord les équations du problème continu, puis de présenter leur discrétisation spatiale, ici en éléments finis, et enfin de décrire la méthode d’intégration temporelle, associée au traitement des non-linéarités matérielles et géométriques.
Discrétisation du problème de dynamique linéaire#
On note \(u\) le champ de déplacements absolus par rapport à la configuration de référence, et paramétré par l’instant \(t\) , appartenant à l’espace affine des champs admissibles \({V}_{\text{adm}}\) .
La méthode directe consiste à résoudre le problème issu de la discrétisation par éléments finis de la formulation en déplacement.
La discrétisation du travail virtuel des forces d’inertie, dans un champ \(\delta v\in {V}_{\text{adm}}^{0}\) , l’espace vectoriel directeur de \({V}_{\text{adm}}\) , s’écrit:
\({\int}_{\Omega}\rho \ddot{u}.\delta vd\Omega {=}^{t}\delta V.M.\ddot{U}\)
La discrétisation de la variation virtuelle du travail dissipé en viscosité (amortissement apporté par une dépendance des contraintes en fonction des vitesses de déformation) est:
\({\int}_{\Omega}c\dot{u}.\delta vd\Omega {=}^{t}\delta V.C.\dot{U}\)
On précise au r5.05.05-terme-inertie comment l’opérateur d’amortissement \(\mathrm{C}\) est construit dans DYNA_NON_LINE.
La discrétisation du travail virtuel des efforts intérieurs en élasticité linéaire s’écrit:
\({\int}_{\Omega}\varepsilon (u).A.\varepsilon (\delta v)d\Omega {=}^{t}\delta V.K.U\)
Enfin, \(\mathrm{L}\) désigne le second membre issu de la discrétisation du travail virtuel des forces extérieures.
En élasticité linéaire, cela aboutit donc à considérer le système différentio-algébrique hyperbolique suivant, pour les degrés de liberté \(\mathrm{U}\) , avec les conditions initiales:
\(\begin{array}{}\text{Trouver}U\in {ℝ}^{n}:\\ M.\ddot{U}+C.\dot{U}+K.U=L\\ U({t}_{0})={U}_{0}\\ \dot{U}({t}_{0})={\dot{U}}_{0}\end{array}\)
accompagnées de conditions aux limites.
Les conditions initiales sont fournies à l’algorithme par le mot-clé ETAT_INIT (opérandes DEPL et VITE).
Si l’état initial résulte d’une simulation en statique linéaire ou non linéaire, on ne prend pas en compte de vitesse initiale, et le déplacement, ainsi que les variables d’état (contraintes, variables internes), sont extraits du résultat de cette simulation à l’instant de départ considéré.
Le système d’équilibre dynamique devient instable si on peut trouver une pulsation \(\omega\) complexe, dont la partie imaginaire est négative, pour laquelle on peut annuler le déterminant de: \(-{\omega}^{2}\mathrm{M}+i\omega \mathrm{C}+\mathrm{K}\) .
Discrétisation du problème de dynamique non linéaire#
On se place maintenant dans un cadre mécanique non linéaire.
On note le travail virtuel \(\underset{\Omega}{\int}\sigma (\mathrm{u},\dot{\mathrm{u}},t).\varepsilon (\delta \mathrm{v})d\Omega\) de déformation (dit aussi des forces internes) du problème de mécanique non linéaire. Cependant on ajoute fréquemment aux forces internes provenant de la loi de comportement, les forces internes d’amortissement visqueux construites de manière globale, cf. [r5.05.05-amortissement]: \(C\cdot \dot{U}\) , où \(\mathrm{C}\) est la matrice d’amortissement. Le travail virtuel des forces internes s’écrit donc après discrétisation:
\({}^{t}\delta \mathrm{V}.\left(\mathrm{R}(\mathrm{U},\dot{\mathrm{U}},t)+\mathrm{C}.\dot{\mathrm{U}}\right){=}^{t}\delta \mathrm{V}.\left({}^{t}\mathrm{Q}\left(\mathrm{U}\right).\sigma \left(\mathrm{U},\dot{\mathrm{U}},t\right)+\mathrm{C}.\dot{\mathrm{U}}\right)\)
où l’on a volontairement distingué les forces visqueuses linéaires (opérateur \(\mathrm{C}\) ) des autres forces internes. Dans le cas des petits déplacements, l’opérateur de déformation assemblé \({}^{t}\mathrm{Q}\) est constant (et défini sur la configuration initiale confondue avec la déformée).
Le bilan d’énergie mécanique s’écrit:
\(\mathrm{L}.\delta \mathrm{v}={\int}_{\Omega}\rho \ddot{\mathrm{u}}.\delta \mathrm{v}d\Omega +{\int}_{\Omega}\sigma \left(\mathrm{u},\dot{\mathrm{u}},t\right).\epsilon \left(\delta \mathrm{v}\right)d\Omega\)
Le champ de contraintes \(\sigma\) à l’instant \(t\) s’écrit de façon générale \(\sigma (\mathrm{u},\dot{\mathrm{u}},\mathrm{Z},t,\mathrm{H})\) , si l’on note \(Z\) le champ de variables de commande, comme par exemple \(T\) le champ de températures, et \(\mathrm{H}\) l’histoire passée de la structure jusqu’à l’instant \(t\) . Pour les comportements incrémentaux, l’histoire est l’ensemble des états (champs de déplacements, contraintes et variables internes) à l’instant précédent.
Dans le cas linéaire ( cf. r5.05.05-discretisation-dynamique-lineaire), le travail virtuel de déformation devient:
\(\mathrm{R}\left(\mathrm{U},\dot{\mathrm{U}}\right)+\mathrm{C}.\dot{\mathrm{U}}=\mathrm{K}.\mathrm{U}+\mathrm{C}.\dot{\mathrm{U}}\)
où \(\mathrm{K}\) est la matrice de rigidité élastique de la structure.
Terme d’inertie#
On note \({}^{t}\delta \mathrm{V}.\mathrm{M}\left(\mathrm{U},\dot{\mathrm{U}},\ddot{\mathrm{U}}\right)\) après discrétisation, le travail virtuel \({\int}_{\Omega}\rho \ddot{\mathrm{u}}.\delta \mathrm{v}\phantom{\rule{0.5em}{0ex}}d\Omega\) des forces d’inertie du système.
On note \(M\) la matrice de masse du système en petites transformations. En grandes rotations de structure (comme par exemple les poutres, cf. [r5.03.40]), ce travail virtuel est une fonction non linéaire de \(U(t)\) et de ses dérivées temporelles (spécifiquement des degrés de liberté de rotation); on fait donc apparaître le terme d’accélération habituel avec une correction non linéaire: \(\mathrm{M}\left(\mathrm{U},\dot{\mathrm{U}},\ddot{\mathrm{U}}\right)=\mathrm{M}\left(\mathrm{U}\right).\ddot{\mathrm{U}}+{\mathrm{L}}_{\text{GR}}^{\text{iner}}\left(\mathrm{U},\dot{\mathrm{U}},\ddot{\mathrm{U}}\right)\) . Dans les autres cas, \(\mathrm{M}\left(\mathrm{U},\dot{\mathrm{U}},\ddot{\mathrm{U}}\right)=\mathrm{M}\left(\mathrm{U}\right).\ddot{\mathrm{U}}\) , qui peut varier si on réactualise la géométrie, ou qui est constante en petits déplacements.
Remarque: cas de la réactualisation géométrique ( DEFORMATION=”PETIT_REAC” ):
Pour un problème dynamique avec réactualisation géométrique, en ayant admis la conservation de la masse, on constate que la matrice de masse n’est pas modifiée. En effet: en notant respectivement \(det\mathrm{J}(\xi )\) et \(\det{J}_{\mathrm{RG}}(\xi )\) le déterminant de la matrice jacobienne de la transformation élément réel \({\Omega}_{e}\) vers élément de référence \({\Omega}_{r}\) et celui de la transformation élément réel réactualisé \({\Omega}_{\mathrm{RG}}\) vers élément de référence –cf. [r3.01.00], §3.5– le travail virtuel des forces d’inertie s’écrit:
\({\int}_{{\Omega}_{\mathit{RG}}}{\rho}_{\mathit{RG}}\ddot{\mathrm{u}}.\delta \mathrm{v}d{\Omega}_{\mathit{RG}}={\int}_{{\Omega}_{r}}{\rho}_{\mathit{RG}}\ddot{\mathrm{u}}.\delta \mathrm{v}.det{\mathrm{J}}_{\mathit{RG}}(\xi )d{\Omega}_{r}={\int}_{{\Omega}_{r}}\rho \ddot{\mathrm{u}}.\delta \mathrm{v}.det\mathrm{J}(\xi )d{\Omega}_{r}\)
car \(\rho_{\mathrm{RG}} = \rho \cdot \det J_{\mathrm{RG}}^{\mathrm{EL}} = \rho \cdot \frac{\det J}{\det J_{\mathrm{RG}}}\) suite à la conservation de la masse dans la réactualisation (jacobien \({J}_{\mathrm{RG}}^{\mathrm{EL}}\) ).
Amortissement#
Il est loisible d’utiliser des éléments discrets sur lesquels on fait porter un comportement d’amortissement via une matrice agissant sur les degrés de liberté, cf. [u4.42.01], mais l’amortissement peut aussi concerner les modèles massifs ou de structures. L’opérateur d’amortissement \(\mathrm{C}\) de ces derniers peut être dans code_aster défini de deux façons dans DYNA_NON_LINE, cf. aussi [r5.05.04], [u2.06.03]:
1) |
une façon globale sur une base de modes propres \(({\Phi}_{k})\) , dite amortissement modal (selon l’hypothèse de BASILE) établie au préalable sur la structure élastique, exprimés sur la base «physique» de la modélisation par éléments finis. On définit ainsi un coefficient par mode choisi. Le mot-clé: AMOR_MODALde l’opérateur DYNA_NON_LINEpermet de lui fournir la base de modes et les coefficients d’amortissement réduit (selon l’hypothèse de BASILE). En effet, les amortissements sont déterminés expérimentalement par analyse modale sur les résonances.
Les déplacements \(\mathrm{U}\) sont donc projetés sur les modes pour obtenir leurs coordonnées généralisées: \({\eta}^{k}{=}^{t}{\Phi}_{k}.\mathrm{U}\) . La matrice d’amortissement modal est:
\(\mathrm{C}=\left(\mathrm{K}{\Phi}_{k}\right).{C}_{\text{modal}}^{k}{.}^{t}\left(\mathrm{K}{\Phi}_{k}\right)\) avec le scalaire \({C}_{\text{modal}}^{k}=2\frac{{\xi}_{k}}{{\omega}_{k}{.}^{t}{\Phi}_{k}.K.{\Phi}_{k}}\) où \({\xi}_{k}\) est le facteur d’amortissement modal à la pulsation \({\omega}_{k}\) et \({k}_{k}{=}^{t}{\Phi}_{k}.K.{\Phi}_{k}\) est la raideur généralisée du mode \(k\) .Malheureusement cette matrice peut avoir un profil très plein et rendreonéreuse (les matrices \(C\) et \(K\) n’ayant pas le même profil), comme on va le voir au [ |
2) |
une façon globale/locale dite amortissement visqueux proportionnel (selon l’hypothèse de RAYLEIGH) à partir des matrices de raideur élastique \(K\) et de masse \(M\) . Les paramètres sont donnés par le matériau sur les éléments finis du modèle (mots-clé AMOR_ALPHA/AMOR_BETAde la commande DEFI_MATERIAU). La matrice d’amortissement visqueux est \(C=\alpha K+\beta M\) . Elle est diagonalisable sur la base des modes propres réels, ce qui rend possible de faire un calcul transitoire sur base modale en découplant les modes: voir [r5.06.04]. Cette formulation, en linéaire, conduit à un facteur d’amortissement lié à la fréquence \(f\) : \(\xi =\alpha \pi f+\beta /\left(4\pi f\right)\) . Dans le cas non linéaire, cette évaluation n’a plus cours. On ajuste en pratique les coefficients \(\alpha ,\beta\) de telle sorte que l’amortissement \(\xi\) soit quasiment uniforme dans la plage \([{f}_{1,}{f}_{2}]\) de fréquence d’intérêt pour la structure étudiée. D’où ainsi de manière raisonnable: \(\alpha =\frac{\xi}{\pi ({f}_{1}+{f}_{2})}\) et \(\beta =\frac{4\pi \cdot \xi \cdot {f}_{1}\cdot {f}_{2}}{{f}_{1}+{f}_{2}}\) Cette modélisation de l’amortissement est peu réaliste si la structure répond sur une large bande de fréquences, car l’amortissement introduit varie alors beaucoup sur la gamme étudiée. Si la loi de comportement du matériau est non linéaire dissipative, le choix des paramètres d’amortissement concerne la plage où la structure reste quasiment élastique. Au cours d’un calcul transitoire non-linéaire, on va utiliser pour la résolution une matrice assemblée de raideur qui peut être la matrice tangente (réactualisée à chaque pas, par exemple) ou la matrice élastique. Historiquement, pour des raisons pratiques, on construisait l’amortissement de Rayleigh en se servant de la matrice de raideur telle qu’elle était spécifiée par l’utilisateur pour le calcul des efforts internes. Donc si l’utilisateur choisissait la matrice tangente, on avait: \(\mathrm{C}=\alpha {\mathrm{K}}^{\mathrm{T}}+\beta \mathrm{M}\) Ceci rend délicat l’interprétation de l’effet de l’amortissement proportionnel. Notamment en cas d’apparition de valeurs propres négatives de la matrice \({K}^{T}\) (par exemple en cas d’endommagement du matériau), l’amortissement peut devenir négatif et renforcer les instabilités! Sans aller jusqu’à cela, sur un simple exemple élasto-plastique, on observe une baisse de l’amortissement visqueux car la raideur tangente baisse. Cette perte d’amortissement pourra être en partie compensée par la dissipation matériau, mais le niveau de dissipation total sera mal maîtrisé. C’est pour cela que l’on a décidé, en version 10 du code, de donner à l’utilisateur le choix pour la définition de l’amortissement de Rayleigh. Soit on adopte la matrice de raideur déjà employée pour les efforts internes (comme dans les versions précédentes du code), soit on spécifie l’usage de la matrice élastique. Ce deuxième choix est particulièrement conseillé pour les lois adoucissantes et pour les lois globales de type GLRC, cf.[R7.01.32]. |
Liaisons#
En pratique, on peut avoir des conditions de liaison bilatérales ou unilatérales, ou des liaisons de type «impédance» ou «absorbantes», cf. [r4.02.05] .
Liaisons bilatérales
Les liaisons bilatérales s’écrivent sous forme de la relation suivante: \(\mathrm{B}.\mathrm{u}={\mathrm{u}}^{\mathrm{d}}\left(t\right)\) .
Elles sont parfaites: elles ne dissipent pas d’énergie. Ce sont des liaisons «holonomes», où la vitesse \(\dot{U}\) n’intervient pas. Ces liaisons sont en général dualisées par Code_Aster , cf. [r3.03.01].
Liaisons unilatérales
Le système mécanique objet de la simulation par éléments finis peut entrer en contact (liaisons unilatérales) avec un «obstacle», qui est un solide dont on connaît a priori le mouvement, d’où la définition d’un jeu \({d}_{0}(t)\) . Les liaisons unilatérales (par exemple le contact unilatéral) s’écrivent sur la configuration à tout instant \(t\) : \(\mathrm{A}\left(\mathrm{u}\right).\mathrm{u}\left(t\right)\le {\mathrm{d}}_{0}\left(t\right)\) (non-pénétration ou vérification que le jeu effectif reste positif ou nul dans toute configuration). L’opérateur \(\mathrm{A}\left(\mathrm{u}\right)\) peut dépendre de la configuration en présence de grands déplacements par réactualisation au cours du pas de temps, cf. [r5.03.50].
On ne considérera que des liaisons «holonomes», du type \(A(U,t)=0\) , qui ne font intervenir que les valeurs des degrés de liberté \(U(t)\) et le temps explicitement si l’obstacle est mobile. On ne considérera pas de liaisons «non holonomes», par exemple du type roulement sans glissement, qui font intervenir la vitesse explicitement et s’écrivent \({\mathrm{A}}_{1}\left(\mathrm{U},t\right).\dot{\mathrm{U}}+{\mathrm{A}}_{2}\left(\mathrm{U},t\right)=0\) , \({\mathrm{A}}_{2}\) et la dépendance directe en temps n’étant présente que s’il y un obstacle mobile.
Si l’obstacle est immobile, la liaison en tant que telle sera explicitement indépendante du temps (on dit aussi «scléronome»).
Ces «charges de contact» sont définies par l’opérateur DEFI_CONTACT. La présence de liaisons unilatérales exige de définir les vitesses du solide dans un espace fonctionnel particulier afin d’assurer l’existence de solutions du système de la dynamique. En effet, lors des instants (dénombrables) d’impact, les vitesses \(\dot{U}({t}^{-})\) et \(\dot{U}({t}^{+})\) peuvent ne pas coïncider. Il est nécessaire pour garantir ce résultat que les données (chargement, équations de liaisons) vérifient une propriété d’analycité, qui est acquise en pratique avec la discrétisation par éléments finis choisie (voir [bib6]).
On n’introduit pas a priori de relation constitutive des impacts (comportement dissipatif en rebond simple, exprimé à l’aide d’un coefficient de restitution normale \(e\in \left[0,1\right]\) ), exprimée sur les différentiels de vitesses des deux points en impact:
\({\dot{\mathrm{U}}}_{2\text{\_norm}}^{+}\left(t\right)={\dot{\mathrm{U}}}_{1\text{\_norm}}^{+}\left(t\right)-e\left({\dot{\mathrm{U}}}_{2\text{\_norm}}^{-}\left(t\right)-{\dot{\mathrm{U}}}_{1\text{\_norm}}^{-}\left(t\right)\right)\)
Ce type de comportement est introduit d’habitude pour traiter le contact-impact de corps rigides, alors que la modélisation numérique avec solides déformables permet de représenter directement le comportement vibratoire sous le choc et les non-linéarités matérielles. Mais il est possible d’ajouter dans la modélisation des éléments discrets de contact-choc placés sur l’interface en contact, munis de la loi DIS_CHOC, qui apportent une dissipation d’amortissement (à condition de supposer de petits mouvements…).
On peut adjoindre aux liaisons unilatérales un comportement de frottement (critère de COULOMB), qui dissipe de l’énergie en glissement relatif des surfaces en contact.
On fait parfois référence au fait que le coefficient de frottement dynamique est plus faible que celui mesuré en quasi-statique (adhérence), voir [bib27], [bib33]. Cela provient du fait qu’en dynamique des vibrations de haute fréquence sur la réaction normale apparaissent et affaiblissent la valeur du seuil de frottement de COULOMB. On n’aurait donc pas besoin de fournir deux valeurs de coefficients à Code_Aster , puisqu’on modélise les solides déformables en contact-frottement (à condition de pouvoir simuler ces vibrations haute fréquence, qui sont généralement atténuées par les schémas d’intégration temporelle…).
On sait que la dissipation est une condition nécessaire pour l’existence de solutions théoriques au problème de dynamique avec impact [bib32]. L’emploi du schéma HHT, voir [r5.05.05-variante-schema-newmark], qui introduit de la dissipation numérique peut s’avérer nécessaire.
Liaisons dissipatives locales
Des relations spécifiques (comme DIS_CHOC) sont conçues pour traiter certains types de liaisons ponctuelles dissipatives, agissant directement sur les degrés de liberté d’éléments discrets du système, voir [r5.03.17]. Elles constituent une loi de comportement en forces généralisées fonction de déplacements généralisés intégrée comme l’ensemble des forces internes \(R(U,\dot{U},t)\) des structures étudiées.
Liaisons absorbantes
Les liaisons de type «absorbantes», cf. [r4.02.05], permettent de simuler le «filtrage» d’une partie de la réponse dynamique, en empêchant la sortie d’ondes réfléchies et diffractées, et évitent que le champ «incident» ne soit perturbé sur une frontière du modèle: voir le [r5.05.05-equations-mouvement]. Elles introduisent des termes amortissants du type \({A}_{\text{abso}}(\dot{u})\) sur une frontière du solide considéré.
Problème dynamique discrétisé#
La dualisation des conditions aux limites de DIRICHLET \(B.u={u}^{d}(t)\) et des conditions unilatérales conduit après discrétisation à définir les inconnues à tout instant \(t\) : \((U,\lambda ,\mu )\) , où \(\lambda\) représente les «multiplicateurs de LAGRANGE» des conditions aux limites de DIRICHLET [r3.03.01], et \(\mu\) représente les «multiplicateurs de Lagrange» des conditions unilatérales.
Le problème dynamique non linéaire s’écrit, avec les conditions initiales[R3.03.01], [r5.03.50]:
Trouver la trajectoire \(U(t)\) :
\(L\) représente le vecteur des forces externes (chargements mécaniques). Ces forces peuvent dépendre du temps et de l’espace. On suppose que, comme les liaisons, elles dépendent de façon régulière des paramètres, ce qui assure l’existence de la solution du problème (théorème de CAUCHY). On peut considérer des forces «suiveuses» \(L(U,t)\) , par exemple la pression, si on prend en compte les changements de géométrie.
Le vecteur \({}^{t}B.\lambda\) s’interprète comme l’opposé des réactions d’appui aux nœuds correspondants (\(B\) est l’opérateur linéaire exprimant le passage aux degrés de liberté des appuis). Le vecteur \({}^{t}A.\mu\) s’interprète comme les forces nodales dues au contact (\(B\) est l’opérateur linéaire exprimant le passage aux degrés de liberté des zones en contact).
L’analyse de stabilité du système d’équilibre dynamique (4894) est plus complexe qu’en linéaire, mais une condition suffisante de perte de stabilité est la possibilité de trouver une pulsation \(\omega\) à partie imaginaire négative pour laquelle on peut annuler le déterminant de: \(-{\omega}^{2}M+i\omega {C}^{T}+{K}^{T}\) , défini sur les opérateurs tangents, à l’instant considéré.
Conditions initiales#
Les conditions initiales \({U}_{0},{\dot{U}}_{0}\) sont fournies au code par le mot-clé ETAT_INIT (opérandes DEPL et VITE).
Si l’état initial résulte d’une simulation en statique linéaire ou non linéaire, le déplacement, ainsi que les variables d’état (contraintes, variables internes), sont extraits du résultat de cette simulation, et la vitesse initiale est par défaut supposée nulle.
Prise en compte d’un état initial précontraint#
Si le problème dynamique à résoudre «suit» un état mécanique initial, deux principales situations s’offrent à nous:
on souhaite calculer les déplacements \(u\) et les contraintes dynamiques à partir de l’état «vierge» initial, où tous les champs sont nuls;
on souhaite calculer les déplacements \(u\) et les contraintes dynamiques en «différentiel» à partir d’un état pré-chargé, dans le cadre élastique.
Dans la première situation, il convient de faire au préalable le calcul de l’état statique, éventuellement non linéaire (matériau, grandes transformations), préalable à la dynamique.
Ainsi, si la structure a un comportement non linéaire, on doit procéder directement en dynamique non linéaire, après avoir évalué l’état initial par simulation en statique, éventuellement non linéaire (avec l’opérateur STAT_NON_LINE), et le champ de déplacement est évalué depuis le début de l’historique. Il peut être nécessaire de prendre en compte les variations de géométrie. Les informations nécessaires pour décrire l’état initial (résultats d’une simulation précédente via le concept résultat EVOL_NOLI ou champs mécaniques nécessaires: DEPL, SIGM, VARI) sont fournies par le mot-clé ETAT_INIT, par exemple si on se trouve dans le cadre d’un comportement incrémental, voir [u4.51.03].
L’état initial peut être obtenu aussi par une simulation en dynamique «très lente», en ayant soin de mettre une rampe «lente» de dépendance au temps sur les efforts statiques appliqués, ainsi qu’un fort amortissement (physique, cf. [r5.05.05-terme-inertie] ou/et numérique cf. [r5.05.05-variante-schema-newmark]). Cette manière de procéder a l’avantage d’injecter tant dans l’opérateur de la phase de prédiction (4293) de l’algorithme de NEWTON, que dans celui de la phase de correction (4297) des termes \(\widehat{\mathrm{K}}\) venant de la matrice de masse et de celle d’amortissement, pour établir l’état mécanique statique. Cela est précieux dans des situations de contact-frottement, d’endommagement… pour améliorer la convergence.
La seconde situation concerne le cas d’une structure qui a subi un pré-chargement thermo‑mécanique «ordinaire», conduisant à un état d’équilibre élastique linéaire. Si l’on mesure par \(\mathrm{u}\) le déplacement à partir de cet état pré-chargé, qui a engendré un état de contraintes \({\sigma}_{1}\) , alors l’énergie de déformation élastique est complétée par un terme de raideur géométrique:
\({\int}_{\Omega}\epsilon \left(\mathrm{u}\right)\phantom{\rule{0.5em}{0ex}}.\mathrm{A}.\epsilon \left(\delta \mathrm{v}\right)\phantom{\rule{0.5em}{0ex}}d\Omega +{\int}_{\Omega}{\sigma}_{1}.\nabla \mathrm{u}\otimes \nabla \left(\delta \mathrm{v}\right)\phantom{\rule{0.5em}{0ex}}d\Omega {=}^{t}\delta \mathrm{V}.\left(\mathrm{K}+{\mathrm{K}}_{\mathrm{g}}\right)\phantom{\rule{0.5em}{0ex}}.\mathrm{U}\)
On assemble alors simplement les matrices \(\mathrm{K}\) et \({\mathrm{K}}_{\mathrm{g}}\) , et on effectue la résolution en dynamique linéaire, comme au [r5.05.05-discretisation-dynamique-lineaire]. Ce peut être par exemple le cas d’une étude sismique sur un barrage voûte. Avec DYNA_NON_LINE, il est nécessaire de fournir par le mot-clé ETAT_INIT, le champ de contraintes SIGM résultat de l’état pré-chargé et de préciser au niveau du mot-clé DEFORMATION sous COMPORTEMENT la prise en compte des termes non linéaires de déformation.
Il est fréquent que la contribution de \({\mathrm{K}}_{\mathrm{g}}\) soit négligeable: on peut alors se contenter d’une analyse dynamique sur la base d’un état initial totalement vierge. On notera aussi que le calcul de la matrice \({\mathrm{K}}_{\mathrm{g}}\) n’est pas disponible pour tous les éléments finis proposés par Code_Aster , voir les docs. [u3].
Problèmes couplés fluide-structure vibro-acoustique#
On pourra pour plus de détails se reporter aux documents [r4.02.02], [r4.02.04], [r4.02.05].
On considère les petits mouvements en approche eulérienne d’un fluide parfait compressible, éventuellement baignant une paroi d’une structure solide. Le fluide est dit barotrope:
on considère de petites perturbations irrotationnelles autour de l’état initial (hydrostatique): \({\vec{\mathrm{U}}}_{\text{fl}}={\vec{\mathrm{u}}}_{\text{fl}}^{0}+{\vec{\mathrm{u}}}_{\text{fl}}\) , \(P={P}_{0}+p\) et \({\rho}_{f}={\rho}_{0}+\rho\) , l’indice 0 désignant la partie permanente des champs,
la loi de comportement du fluide donne les contraintes «fluctuantes»: \(\sigma =-p\text{Id}=-\rho {c}_{0}^{2}\text{Id}\approx {\rho}_{0}{c}_{0}^{2}\left(div\vec{{\mathrm{U}}_{\text{fl}}}\right)\phantom{\rule{0.5em}{0ex}}\text{Id}\) , car \({\rho}_{0}={\rho}_{f}\left(1+div{\vec{\mathrm{U}}}_{\text{fl}}\right)\) ,
les vitesses fluides dérivent d’un potentiel \({\vec{\mathrm{v}}}_{\text{fl}}={\dot{\vec{\mathrm{u}}}}_{\text{fl}}=\overrightarrow{\nabla \dot{\phi }}\) et sont modélisées à l’aide des champs \(\left(p,\phi \right)\) : pression fluctuante et potentiel de déplacement, qui ne sont pas indépendants car: \(\dot{p}=-{\rho}_{0}{c}_{0}^{2}\Delta \dot{\phi }\) (en combinant équation de continuité et loi de comportement).
On admet que l’on ne considère pas de forces fluctuantes de volume \(\vec{f}\) s’exerçant sur le fluide. L’équilibre dynamique du fluide s’écrit: \(\vec{\nabla p}+{\rho}_{f}{\vec{\ddot{\mathrm{u}}}}_{\text{fl}}=\vec{0}\) (équation d’Euler linéarisée), qui est valable pour un fluide non pesant compressible ou pesant incompressible; par contre, pour un fluide pesant compressible, les approches eulérienne et lagrangienne ne coïncident pas même en petits mouvements: ce cas n’est pas traité par Code_Aster .
L’équilibre dynamique du fluide va être écrit sous forme variationnelle sous l’action d’une pression fluctuante \({p}_{\text{imp}}\) imposée sur une partie de la frontière. En régime harmonique, l’équilibre dynamique du fluide se traduit par la formulation variationnelle de l’équation de Helmholtz, [bib26].
Code_Aster possède une formulation symétrisée, voir [r4.02.02], avec des éléments \((P,\varphi )\) où \(P\) est le vecteur des degrés de liberté de pression et \(\varphi\) est le potentiel pour décrire les perturbations (fluctuations) dans le fluide, sachant que \(\overrightarrow{\nabla p}+{\rho}_{0}\overrightarrow{\nabla \ddot{\varphi}}=\overrightarrow{0}\) . Une équation en \(P\) résout l’équilibre dynamique dans le domaine fluide, celle en \(\varphi\) traduisant l’équation de continuité dérivée combinée à la loi de comportement fluide. Des conditions aux limites en \(P\) et \(\varphi\) complètent le système d’équations pour décrire les évolutions du fluide. Ainsi, à cause de la formulation \((P,\phi )\) , si sur une frontière \(\partial {\Omega}_{\text{f\_p}}\) du domaine fluide, une pression fluctuante est appliquée: \(P={P}_{\text{imp}}(t)\) , on doit aussi y imposer une condition sur \(\varphi\) : \(\varphi (t)={\varphi}_{\text{imp}}(t)\) , vérifiant \({\rho}_{0}{\ddot{\varphi}}_{\text{imp}}(t)=-{P}_{\text{imp}}(t)\) .
Comme on envisage une frontière commune fluide-structure \(\partial {\Omega}_{\text{FS}}\) , où on définit la normale \(\overrightarrow{n}\) sortante du domaine structure vers le fluide, le chargement de paroi du fluide est couplé avec le déplacement de la structure. Les déplacements normaux sont continus: \((\overrightarrow{{u}_{\mathrm{fl}}}+{\overrightarrow{u}}_{\mathrm{fl}}^{0})\cdot \overrightarrow{n}={\overrightarrow{u}}_{\mathrm{st}}\cdot \overrightarrow{n}\) sur \(\partial {O}_{\text{FS}}\) , \({\overrightarrow{u}}_{\text{fl}}^{0}\) désignant la partie permanente du déplacement fluide. En utilisant le potentiel fluide, on a: \(\overrightarrow{\nabla \dot{\varphi}}.\overrightarrow{n}={\overrightarrow{\dot{u}}}_{\text{st}}.\overrightarrow{n}\) sur \(\partial {\Omega}_{\text{FS}}\) . De manière duale, les vecteurs contraintes sont continus: \(-p\overrightarrow{n}=\sigma .\overrightarrow{n}\) sur \(\partial {\Omega}_{\text{FS}}\) . La structure reçoit ainsi le chargement fluctuant du fluide: \({\int}_{\partial {\Omega}_{\text{FS}}}p\partial \overrightarrow{v}.(-\overrightarrow{n})\text{dS}\) .
Si on envisage une surface libre (soumise à une pression constante), traitée aussi en description eulérienne, voir [r4.02.04], on note par \(z\) (et \(Z\) les degrés de liberté associés après discrétisation) l’altitude fluctuante (petite) de la surface libre \(\partial {\Omega}_{\text{SL}}\) et la fluctuation de pression eulérienne dans le fluide vérifie: \({\int}_{\partial {\Omega}_{\text{SL}}}(p-{\rho}_{0}\text{gz})\delta z\text{dS}\) , dans toute altitude virtuelle \(\delta z\) . C’est la seule conséquence de la pesanteur que l’on puisse prendre ne compte dans cette formulation.
Par ailleurs, on peut avoir besoin de prendre en compte une frontière artificielle avec un milieu infini (qui doit traiter la condition de radiation à l’infini): Code_Aster propose des éléments finis de frontière absorbante (éléments paraxiaux ou anéchoïques), voir [r4.02.05]. Comme ils apportent a priori un terme en dérivée troisième du temps, à cause de l’introduction du champ \(\varphi\) , on préfère le traiter à l’aide d’une transformation en un terme non symétrique en \(\dot{P}\) qui sera reporté au second membre, de manière explicite.
Au total on obtient les équations semi-discrètes différentio-algébriques du problème couplé:
accompagné des conditions initiales:
\(U({t}_{0})={U}_{0}\) , \(\dot{U}({t}_{0})={\dot{U}}_{0}\) , \(P({t}_{0})={P}_{0}\) , \(\varphi ({t}_{0})={\varphi}_{0}\) et \(Z({t}_{0})=0\) ,
et des conditions aux limites: \(U(t)={U}_{\text{imp}}\) sur le bord \(\partial {\Omega}_{{S}_{u}}\) de la structure, des éventuelles conditions unilatérales, et \(P(t)={P}_{\text{imp}}(t)\) avec \(\varphi (t)={\varphi}_{\text{imp}}(t)\) , vérifiant \({\rho}_{0}{\ddot{\varphi}}_{\text{imp}}(t)={P}_{\text{imp}}(t)\) sur le bord \(\partial {\Omega}_{\text{f\_p}}\) du fluide.
Remarque :
On note que dans le cas non linéaire, on remplace dans (4895) le terme \(KU\) par les forces internes non linéaires \(R(U,\dot{U},Z,t)\) .
Les différents opérateurs sont:
les matrices \(K\) , \(M\) , \(C\) définies plus haut pour la structure solide,
\({Q}_{\text{fl}}\) est la matrice construite à partir de \({\int}_{{\Omega}_{f}}\frac{1}{{\rho}_{0}{c}_{0}^{2}}p.\text{qd}\Omega\) , qui a le sens physique d’une énergie élastique du fluide,
\({H}_{\text{fl}}\) est construite à partir de \({\int}_{{\Omega}_{f}}-{\rho}_{0}\overrightarrow{\nabla \varphi }.\overrightarrow{\nabla \Psi }d\Omega\) , et décrit le transport de masse fluide,
\({M}_{\text{fl}}\) est construite à partir de \({\int}_{{\Omega}_{f}}\frac{1}{{c}_{0}^{2}}\varphi .\text{qd}\Omega\) , et décrit l’inertie du fluide,
\({M}_{\text{FS}}\) est construite à partir de \({\int}_{\partial {\Omega}_{\text{FS}}}{\rho}_{0}\varphi .\delta \overrightarrow{v}.(-\overrightarrow{n})\text{dS}\) (\(\overrightarrow{n}\) est la normale du domaine structure vers le fluide), et décrit le débit massique eulérien à l’interface fluide-structure,
\({A}_{\text{fa}}\) est construite à partir de \({\int}_{{\Omega}_{f}}{\rho}_{0}p.\psi d\Omega\) , et désigne l’opérateur frontière absorbante, agissant sur \(\dot{P}\) , modifiant l’équation en \(\varphi\) de la paroi absorbante: les éléments (modèle 3D_FLUI_ABSO) absorbants dé-symétrisent le système issu de la formulation \((P,\varphi )\) et on va reporter ce terme au second membre, qui sera noté \({L}^{\text{abso}}\) , par discrétisation temporelle explicite à chaque itération, voir le paragraphe
r5.05.05-schema-integration-temporelle,\({K}_{z}\) est la «raideur» de surface libre, construite à partir de \({\int}_{\partial {\Omega}_{\text{SL}}}{\rho}_{0}\text{gz}.\delta z\text{dS}\) ,
\({M}_{z}\) provient du travail de la pression fluctuante en surface libre, construite à partir de \({\int}_{\partial {\Omega}_{\text{SL}}}{\rho}_{0}z.\psi \text{dS}\) ,
le terme \({L}_{\text{st}}\) contient, entre autres chargements, l’effet de la pression hydrostatique exercée par le fluide sur la structure.
En résumé, la prise en compte d’un domaine fluide en évolution fluctuante barotrope, interagissant avec la structure conduit à considérer dans le système dynamique non linéaire (4894) enrichi sur des degrés de liberté particuliers \((P,\varphi )\) et \(Z\) :
un opérateur de forces d’inertie \(M(U).\ddot{U}\) enrichi par:
\(\mathrm{M}^{\text{fs}} \left( \begin{array}{c} \ddot{\mathrm{U}} \\ \ddot{\mathrm{P}} \\ \ddot{\varphi} \\ \ddot{\mathrm{Z}} \end{array} \right) = \left( \begin{array}{cccc} 0 & 0 & \mathrm{M}_{\text{FS}} & 0 \\ 0 & 0 & \mathrm{M}_{\text{fl}} & 0 \\ {}^{t}\mathrm{M}_{\text{FS}} & {}^{t}\mathrm{M}_{\text{fl}} & \mathrm{H}_{\text{fl}} & \mathrm{M}_{\mathrm{z}} \\ 0 & 0 & {}^{t}\mathrm{M}_{\mathrm{z}} & 0 \end{array} \right) \left( \begin{array}{c} \ddot{\mathrm{U}} \\ \ddot{\mathrm{P}} \\ \ddot{\phi} \\ \ddot{\mathrm{Z}} \end{array} \right)\)
un opérateur de forces intérieures \(R(U,\dot{U})\) enrichi par:
\(K^{\text{fs}} \left( \begin{array}{c} U \\ P \\ \phi \\ Z \end{array} \right) = \left( \begin{array}{cccc} 0 & 0 & 0 & 0 \\ 0 & Q_{\text{fl}} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & K_{z} \end{array} \right) \left( \begin{array}{c} U \\ P \\ \phi \\ Z \end{array} \right)\)
un second membre enrichi par le report au second membre par discrétisation temporelle explicite de \(-{A}_{\mathrm{fa}}\dot{P}\) sur les degrés de liberté \(\varphi\) .
Le fluide doit rester en petits mouvements (hypothèse de base de cette modélisation), mais on peut considérer des grands mouvements de la structure, baignée par le fluide, via une réactualisation de la géométrie \(\mathrm{X}\to \mathrm{X}+\Delta \mathrm{U}\) (via le mot-clé COMPORTEMENT, opérande DEFORMATION:”PETIT_REAC”, valable si on considère de petites rotations) des frontières \(\partial {\Omega}_{\text{FS}}\) , ce qui fait recalculer le terme \({M}_{\text{FS}}\) mais aussi tous les autres puisque le domaine \({\Omega}_{f}\) a évolué; les champs scalaires \((P,\phi )\) sont alors simplement transportés à l’identique sur la géométrie réactualisée. Voir aussi le cas-test FDNV100 [v8.03.100].
Prise en compte de lois de comportement visqueux et amortissement#
Une loi de comportement visqueux, voir [r5.03.08], se traduit comme une loi élastoplastique par une évolution du travail des forces intérieures \(R(U,\dot{U},t)\) . Ainsi, tout en amenant un amortissement «physique» dans l’équilibre dynamique, elle ne produit pas de contribution directe au terme d’amortissement \(C.\dot{U}\) de l’équation d’équilibre dynamique, mais cependant de manière indirecte si on a choisi un amortissement de Rayleigh ( cf. r5.05.05-terme-inertie) via la matrice de raideur tangente du schéma d’intégration.
En effet avec une loi de comportement visqueux, le tenseur des déformations comporte une partie élastique, une partie thermique, une partie anélastique (connue) et une partie visqueuse, déviatorique (déviateur des contraintes noté \(\tilde{\sigma}\) ), par exemple vérifiant:
\(\begin{array}{c}{\epsilon}_{\text{tot}}={\epsilon}_{\mathrm{e}}+{\epsilon}_{\text{th}}+{\epsilon}_{\mathrm{a}}+{\epsilon}_{\mathrm{v}}\\ \sigma =\mathrm{A}\left(T\right).{\epsilon}_{\mathrm{e}}\\ {\dot{\epsilon}}_{v}=g\left({\sigma}_{\text{eq}},\lambda ,T\right)\frac{3}{2}\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}\end{array}\)
Pour la relation de comportement visqueux LEMAITRE, la fonction \(g\) est explicite, mais ce n’est pas toujours le cas. Après discrétisation implicite en temps, l’écoulement visqueux est:
\(\frac{\Delta {\epsilon}_{\mathrm{v}}}{\Delta t}=\frac{3}{2}g\left({\sigma}_{\text{eq}},{\lambda}^{-}+{\left(\Delta {\epsilon}_{v}\right)}_{\text{eq}},T\right)\frac{\stackrel{~}{\sigma}}{{\sigma}_{\text{eq}}}\)
On peut aussi adopter un schéma semi-implicite, qui semble donner de meilleurs résultats:
\(\frac{\Delta {\varepsilon}_{v}}{\mathrm{\Delta t}}=\frac{3}{2}g({({\sigma}^{-}+\frac{\Delta \sigma }{2})}_{\text{eq}},{\lambda}^{-}+\frac{{(\Delta {\varepsilon}_{v})}_{\text{eq}}}{2},{T}^{-}+\frac{\mathrm{\Delta T}}{2})\frac{({\tilde{\sigma}}^{-}+\frac{\Delta \tilde{\sigma}}{2})}{{({\sigma}^{-}+\frac{\Delta \sigma }{2})}_{\text{eq}}}\)
Après résolution d’une équation non linéaire locale par élimination de \(\Delta {\epsilon}_{\mathrm{v}}\) pour calculer \({\sigma}_{\text{eq}}={\left({\sigma}^{-}+\mathit{\Delta \sigma }\right)}_{\text{eq}}\) , on obtient ainsi la contrainte à la fin du pas de temps courant \(\sigma ={\sigma}^{-}+\Delta \sigma\) .
Equations en mouvement « relatif » [r4.05.01]#
Généralités#
Dans de nombreuses applications, notamment en séisme, on souhaite calculer directement le champ de déplacement de la structure déduit du mouvement d’«entraînement» provenant des déplacements imposés \({\mathrm{U}}_{d}\left(t\right)\) des supports de la structure.
On note alors \({u}_{a}\) le déplacement absolu de la structure: \({\mathrm{u}}_{\mathrm{a}}={\mathrm{u}}_{\text{ent}}+\mathrm{u}\) , \({\mathrm{u}}_{\text{ent}}\) étant le déplacement d’entraînement (par exemple, \({\mathrm{u}}_{\text{ent}}\) peut être un champ incident: il est appelé alors «déplacement pseudo-statique»). Dans le cas MONO_APPUI, c’est un mouvement de corps rigide. Dans le cas MULTI_APPUI, \({\mathrm{u}}_{\text{ent}}\) est défini comme étant la solution statique linéaire de la structure soumise aux déplacements imposés aux appuis, voir [r5.05.05-decomposition-relevement-elastostatique].
Et on appelle \(u\) le déplacement «relatif» (appelé ainsi par abus de langage si \({\mathrm{u}}_{\text{ent}}\) n’est pas de corps rigide).
En effet la codification (RCC-M, ASME…) introduit la distinction entre «contraintes primaires» dues au mouvement vibratoire «relatif» et «contraintes secondaires» dues au mouvement vibratoire d’«entraînement». La pertinence de cette distinction disparaît a priori évidemment dès qu’on envisage un comportement non linéaire du matériau.
Si l’on traite le problème d’interaction dynamique avec le sol (qui est un demi-espace infini), en séisme par exemple, le champ de déplacement relatif \(\mathrm{u}\) vérifie: \(\underset{\mathrm{x}\to \infty }{\lim}\mathrm{u}\left(\mathrm{x}\right)=0\) : seul le champ incident \({\mathrm{u}}_{\text{ent}}\) est perceptible à l’infini– c’est la donnée de chargement sismique. On utilise pour définir ce chargement de déplacement imposé la commande AFFE_CHAR_MECA et le mot-clé facteur ONDE_PLANE, sur une frontière donnée du maillage considéré.
Dans ce cas de problème d’interaction avec le sol, on ne connaît pas a priori le déplacement d’«entraînement» directement appliqué à la structure, puisqu’il résulte de la réponse totale couplée: aussi le cas MULTI_APPUI, où le déplacement est nécessairement connu à l’interface, n’est-il pas applicable.
Par contre, ne pouvant pas simuler en éléments finis avec Code_Aster le champ dans tout le demi-espace infini (le sol), on est conduit à placer des frontières «élastiques absorbantes», cf. [r4.02.05], au bord du maillage de sol. Le travail virtuel associé à ces frontières absorbantes, de normale sortante \(n\) , est traité comme un second membre (pour les éléments finis absorbants paraxiaux de Code_Aster , qui sont d’ordre 0), car il est intégré de façon explicite dans le schéma d’intégration temporelle ( cf. [r5.05.05-schema-integration-temporelle]) de DYNA_NON_LINE: on prend en effet la valeur \({\dot{u}}^{-}\) de la vitesse à l’instant précédent. La forme linéaire associée vaut:
Après discrétisation spatiale, on déduit le second membre:
Remarque: cas MONO_APPUI :
Pour un problème dynamique traité en déplacements relatifs (cas MONO_APPUI), seul le premier terme de (4897) subsiste. On pourra se reporter à [r4.02.05].
Remarque: cas d’un problème avec interaction fluide-structure:
Pour une structure subissant des déplacements imposés, en présence d’interaction fluide-structure, où le domaine fluide n’est pas directement lié au «support» imposant le signal d’entraînement, il est possible de résoudre le système d’équilibre dynamique en terme de déplacement «relatif» \(U\) de la structure et de variables fluides \((\mathrm{P},\varphi )\) «absolues» et de cote de surface libre du fluide \(Z\) «absolue», cf. [r5.05.05-problemes-couples-fluide-structure-vibro-acoustique]. En effet, on peut montrer que \(({P}_{\text{ent}},{\varphi}_{\text{ent}})=0\) et \({Z}_{\text{ent}}=0\) , sur la base de (4895). On pourra se reporter à [r4.02.05].
On ne pourra pas considérer un tel type de décomposition champ d’«entraînement»–champ «relatif» en cas de chargements de type pression imposée fluctuante sur une paroi du domaine fluide.
On cherche à exploiter cette distinction dans le système dynamique non linéaire discret (4894) pour simplifier la prise en compte des déplacements imposés \({U}_{d}(t)\) .
Fig. 197 Décomposition du mouvement absolu#
On sépare les conditions de liaison, voir (4894), en un groupe de conditions d’«entraînement» par les supports: \({B}_{s}.{U}_{a}={U}_{d}(t)\) et un groupe de conditions de liaisons notées par \({B}_{L}.{U}_{a}=0\) que l’on souhaite imposer directement au déplacement absolu \({U}_{a}\) de la structure (par exemple des liaisons internes comme «3D_POU»…) auxquelles sont associés leurs paramètres \({\lambda}_{L}\) de LAGRANGE. On envisage donc par la suite ces deux familles: \({B}_{s}\) pour les mouvements d’«entraînement» par les supports \(\Gamma\) et \({B}_{L}\) pour les liaisons bilatérales «ordinaires».
Si les supports sont en nombre fini (ce qui sera le cas après discrétisation de toute manière), on note \({V}_{e}\) l’espace vectoriel des champs de déplacements de la structure «entraînée» \({u}_{\text{ent}}\) , de dimension finie \({N}_{s}\) , que l’on va définir ci-après. On décompose les conditions de DIRICHLET \({U}_{d}(t)\) , sur une base \(({X}_{\text{ks}})\) de déplacements des supports: \({U}_{d}(t)={U}_{d}^{k}(t){X}_{\text{ks}}\) , \(k=1\to {N}_{s}\) parcourant tous les «degrés de liberté entraînés» par les supports.
Décomposition par relèvement élastostatique#
On construit un «relèvement élastostatique», c’est-à-dire une base de \({V}_{e}\) à partir des solutions statiques élastiques linéaires de la structure sous les seuls déplacements imposés de base \(({X}_{\text{ks}})\) des supports (aucun chargement en force imposée). Après discrétisation par éléments finis, cela revient à résoudre \(k=1\to {N}_{s}\) problèmes d’élastostatique (matrice de raideur \(K\) ):
Trouver \(({\Psi}_{k},{\lambda}_{s{\Psi}_{k}},{\lambda}_{L{\Psi}_{K}})\) tels que
On appelle «modes statiques linéaires» ces \({N}_{s}\) solutions \({\Psi}_{k}\) (renseignés via l’opérande «MODE_STAT» de DYNA_NON_LINE). Ils sont calculés au préalable par l’opérateur MODE_STATIQUE [u4.52.14] avec l’option MODE_STAT.
On a nécessairement avec (4898): \({}^{t}{\Psi}_{\ell}.K.{\Psi}_{k}{+}^{t}{X}_{\mathrm{\ell s}}.{\lambda}_{s{\Psi}_{k}}=0\) , \(\forall \ell ,k=1\to {N}_{s}\) .
Remarque:
Il est essentiel que les mêmes liaisons \({B}_{L}\) soient représentées dans la résolution du problème de recherche de modes statiques (4898) que dans le problème dynamique complet.
Le champ de déplacements de la structure «entraînée» \({u}_{\text{ent}}\) , après discrétisation par éléments finis, est donc décrit par les \({U}_{\text{ent}}={U}_{\text{ent}}^{k}(t){\Psi}_{k}\) , vérifiant notamment \({B}_{s}.{U}_{\text{ent}}={U}_{d}^{k}(t){X}_{\text{ks}}\) sur les supports, parcourant le sous-espace discret \({V}_{e}\) , des «degrés de liberté entraînés» par les supports.
Le sous-espacediscret \({V}_{e}\) est donc engendré par la base \(({\Psi}_{k})\) . On a nécessairement par construction:
\({U}_{\text{ent}}^{k}(t)={U}_{d}^{k}(t),\forall k\)
Les degrés de liberté de déplacements de la structure «entraînée» sont donc directement donnés par la valeur des \({U}^{d}(t)\) .
Ayant caractérisé le sous-espace discret \({V}_{e}\) à partir des modes statiques linéaires, après discrétisation, étudions un champ de déplacement \(W\) de la structure, cinématiquement admissible quelconque, mais nul sur les supports: \({B}_{s}.W=0\) et vérifiant les liaisons «ordinaires» \({B}_{L}.W=0\) . En vertu de (4898), on a nécessairement:
\(\left\lbrace \begin{array}{l} {}^{t}W \cdot K \cdot \Psi_k + {}^{t}W \cdot {}^{t}B_s \cdot \lambda_{s\Psi_k} + {}^{t}W \cdot {}^{t}B_L \cdot \lambda_{L\Psi_k} = 0 \\ {}^{t}W \cdot B_s \cdot \Psi_k = {}^{t}W \cdot X_{\text{ks}} \\ {}^{t}W \cdot B_L \cdot \Psi_k = 0 \end{array} \right. \quad \forall W, \text{ tel que } B_s \cdot W = 0\)
D’où simplement le problème statique linéaire:
L’opérateur de raideur élastique \(K\) étant défini positif (ayant éliminé les modes de corps rigides), on constate que tout champ de déplacement absolu \({W}_{a}\) cinématiquement admissible de la structure, après discrétisation, peut s’écrire par décomposition unique:
avec \({V}_{e}\) engendré à partir des modes statiques \(({\Psi}_{k})\) , et \(V\) contenant les champs dits «degrés de liberté actifs» tels que \({B}_{s}.W=0\) (nuls sur les supports). On appelle \(V\) le sous-espace discret des «degrés de liberté actifs».
Cas MONO_APPUI et MULTI_APPUI#
Pour un chargement de type MONO_APPUI, les modes statiques sont simplement les modes de corps rigides de la structure: \(K.{\Psi}_{k}=0\) , vérifiant les liaisons «ordinaires» \({B}_{L}\cdot {\Psi}_{k}=0\) .
Si le chargement est MULTI_APPUI, les modes statiques \(({\Psi}_{k})\) sont quelconques.
Cette décomposition sur deux sous-espaces supplémentaires \({V}_{a}=V\oplus {V}_{e}\) de tout champ cinématiquement admissible bâtie à l’aide de l’opérateur d’élasticité \(K\) est applicable dans toute évolution non linéaire de la structure, y compris avec des chocs…, à condition que les liaisons bilatérales restent les mêmes au cours de l’historique.
Exploitons cette décomposition et projetons maintenant le problème dynamique non linéaire (4894) séparément sur le sous-espace \(V\) puis sur le sous-espace \({V}_{e}\) , en exploitant (4898). Le résultat se simplifie (un peu seulement) car \({B}_{s}.W=0\) (c’est-à-dire les «degrés de liberté actifs» ne travaillent pas dans les réactions des appuis supports \({B}_{s}\) ), \({B}_{L}.{\Psi}_{k}=0\) (c’est-à-dire les modes statiques ne travaillent pas dans les réactions de liaison \({B}_{L}\) ), \({B}_{s}.{\Psi}_{k}={X}_{\text{ks}}\) :
Trouver \({U}_{a}={U}_{e}+U\) , \({\lambda}_{s}\) , \({\lambda}_{L},\mu\) tels que:
On constate que cela est assez compliqué.
Restreignons-nous d’abord au problème dynamique où les opérateurs d’inertie \(M\) et de forces intérieures \(R\) sont linéaires, et en absence de frontières absorbantes: \({L}^{\text{abso}}({U}_{e},{\dot{U}}_{e},\dot{U})=0\) .
On fait l’hypothèse dans Code_Aster que: \(C.{\dot{U}}_{e}=0\) , y compris en multi-appuis (alors que cela n’est exact qu’en mono-appui, où les modes statiques sont les modes rigides \(K.{\Psi}_{k}=0\) sans déformation). Cela revient à négliger la contribution des déplacements d’entraînement de la structure aux forces d’amortissement visqueux.
Les modes statiques \({\Psi}_{k}\) vérifiant (4898), le système (3682) se restreint à:
On constate sur la première de ces équations (3683), que grâce aux hypothèses faites, on peut se restreindre à résoudre un problème dynamique sur le champ de déplacement «relatif» \(U\) , ayant bloqué les degrés de liberté sur les supports (\({B}_{s}.U=0\) ), à condition de fournir au préalable le terme \({}^{t}W.M.{\ddot{U}}_{e}\) , ainsi que les modes statiques \(({\Psi}_{k})\) .
L’objet de l’opérateur CALC_CHAR_SEISME [u4.63.01] est justement de calculer le terme \(-{}^{t}W.M.{\ddot{U}}_{e}\) , transformé en un concept de type «charge», à l’aide de l’opérateur AFFE_CHAR_MECA [u4.25.01]. On peut aussi simplement introduire une charge de «pesanteur» unitaire dans la direction voulue, et amplifiée par le signal temporel d’accélération. L’avantage est d’exploiter directement la donnée en accélérogramme (par exemple produit à partir d’un spectre), sans avoir à l’intégrer deux fois en temps avec les incertitudes que cette opération engendre. On peut produire directement les contraintes dites «primaires» induites par la dynamique en «mouvement relatif».
Cet avantage est perdu si des liaisons unilatérales sont présentes, puisque la condition \(A.(U+{U}_{d}^{k}(t){\Psi}_{k})\le {d}_{0}(t)\) demande la valeur de \({U}_{d}^{k}(t)\) pour être exprimée correctement, excepté si on a la chance que le mouvement d’entraînement n’ait pas de composante sur l’interface où le jeu de la liaison unilatérale est calculé, c’est-à-dire si les butées de ces appuis unilatéraux sont fixes dans le repère absolu!
La seconde équation de (3683) fournit les réactions sur les supports \({\lambda}_{s}\) (le reste étant déterminé par la résolution du problème en déplacement «relatif» \(U\) ); mais on constate là aussi qu’il est nécessaire de connaître \({U}_{d}^{k}(t)={U}_{e}^{k}(t)\) sur les appuis.
C’est la même chose si l’on veut reconstituer la solution complète pour le post-traitement en contraintes avec leur contribution dite «secondaire», liée à \({U}_{e}={U}_{d}^{k}(t){\Psi}_{k}\) . Cette contribution est nulle en mono-appui, puisqu’il s’agit alors un mode rigide.
Maintenant on considère le problème dynamique d’une structure interagissant avec un sol (milieu «infini»), source d’une onde sismique incidente, cf. [r4.05.01] et [r4.02.05]. On est nécessairement dans un cadre «MONO_APPUI», où les modes statiques sont les modes rigides de la structure, d’où: \(K.{\Psi}_{k}=0\) et on admet habituellement que \({t}^{}{\Psi}_{k}.\mathrm{C}.\dot{\mathrm{U}}=0\) (en effet, physiquement, l’amortissement provient des déformations dans le solide). On utilise donc des éléments de frontière absorbante, et le terme \({L}_{\text{abso}}({U}_{e},{\dot{U}}_{e},\dot{U})\) est présent dans[éq2.6-7]. L’onde incidente est fournie directement par le signal \({U}_{d}(\overrightarrow{x},t)\) . On construit donc les termes de[éq 2.6-2], ainsi que le terme \({}^{t}W.M.{\ddot{U}}_{e}\) . Le système[éq2.6-7] devient dans le cas mono-appui:
Revenons maintenant au problème général non linéaire (3682). On constate que la première équation sur les degrés de liberté «actifs» comporte nécessairement un couplage avec le champ \({\mathrm{U}}_{e}={\mathrm{U}}_{d}^{k}(t){\Psi}_{k}\) dans l’opérateur d’inertie comme dans celui de forces internes. On doit donc, conformément au schéma d’intégration temporelle et de résolution non linéaire, développé au [r5.05.05-schema-integration-temporelle], reconstituer à chaque instant de calcul la valeur de déplacement absolu \({U}_{a}={U}_{e}+U\) , les contraintes…
En conclusion, on peut traiter le problème dynamique en mouvement relatif (en admettant que les forces d’amortissement ne dépendent que de lui), selon les conditions suivantes:
situation |
Possibilité de MONO_APPUI |
Possibilité de MULTI_APPUI |
comportement linéaire en petites transformations, sans frontière absorbante, avec des conditions unilatérales telles que les jeux ne soient pas modifiés par le mouvement d’entraînement, que le domaine fluide ne soit pas directement chargé, et avec un chargement d’accélération imposée \({\ddot{U}}_{e}\) , cf. (3683) |
oui |
oui |
comportement quelconque et avec frontière absorbante, avec des conditions unilatérales telles que les jeux ne soient pas modifiés par le mouvement d’entraînement, et avec un chargement d’accélération imposée \({\ddot{U}}_{e}\) |
oui |
NON |
situation |
Possibilité de MONO_APPUI |
Possibilité de MULTI_APPUI |
comportement linéaire en petites transformations, sans frontière absorbante, avec des conditions unilatérales telles que les jeux ne soient pas modifiés par le mouvement d’entraînement, que le domaine fluide ne soit pas directement chargé, et avec un chargement d’accélération imposée \({\ddot{U}}_{e}\) , cf. (3683) |
oui |
oui |
comportement quelconque et avec frontière absorbante, avec des conditions unilatérales telles que les jeux ne soient pas modifiés par le mouvement d’entraînement, et avec un chargement d’accélération imposée \({\ddot{U}}_{e}\) |
oui |
NON |
Dans le cas «NON», on ne peut que traiter le problème dynamique dans le repère absolu.
Dans le cas «NON», on ne peut que traiter le problème dynamique dans le repère absolu.
Schéma d’intégration temporelle : schéma de NEWMARK et méthode de NEWTON#
Le problème mécanique à analyser étant modélisé en éléments finis, selon la démarche décrite au [r5.05.05-dynamique-non-lineaire], on calcule les champs de déplacements, vitesses et accélérations aux nœuds en une suite discrète d’instants de calcul \({t}_{1},{t}_{2}\dots {t}_{i-1},{t}_{i}\dots {t}_{N}\) : \({\left\lbrace \lbrace {t}_{i}\rbrace \right\rbrace }_{1\le i\le N}\) .
L’utilisateur de DYNA_NON_LINE peut actuellement choisir entre cinq schémas temporels à un pas, dont trois implicites: celui de NEWMARK (1959), sa variante dite «accélération moyenne modifiée», ou celui de HILBER-HUGUES-TAYLOR (HHT, 1977): voir le paragraphe [r5.05.05-variante-schema-newmark] et deux explicites: celui des différences centrées (qui est un cas particulier de la famille de Newmark) et le schéma de Tchamwa-Wielgosz (qui introduit une dissipation numérique haute fréquence) ). L’état de la structure étant connu à l’instant \({t}_{i-1}\) , on cherche à calculer son état à l’instant \({t}_{i}\) par une méthode de prédiction-correction.
Schéma de NEWMARK#
On présente ici ce schéma sous sa forme classique ([bib1] et [bib2]) relative à un mouvement de translation ou de petite rotation. Pour les grandes rotations d’éléments de structure [bib3], les formules sont plus compliquées, mais elles permettent de la même façon d’actualiser la vitesse et l’accélération angulaires en fonction de l’accroissement de déplacement, qui est dans ce cas le vecteur-rotation.
On note ci-après par \(\Phi\) la configuration, c’est-à-dire le paramétrage du système par les degrés de liberté des éléments finis: déplacements et rotations \(U\) , pression \(P\) , potentiel \(\varphi\) …
Le schéma de NEWMARK repose sur les développements suivants du vecteur-configuration, fonction du temps, où \(\beta\) et \(\gamma\) sont deux paramètres:
Il s’agit de la paramétrisation en accélération, dite «A-form», où on exprime les déplacements et les vitesses avec les accélérations. L’équation (4900) s’écrit aussi avec (4901), si \(\gamma \ne 0\) :
\(\Phi (t+\Delta t)\approx \Phi (t)+\frac{\Delta t}{\gamma}((\gamma -\beta )\dot{\Phi}(t)+\beta \dot{\Phi}(t+\Delta t))+\frac{\Delta {t}^{2}}{2\gamma }(\gamma -2\beta )\ddot{\Phi}(t)\)
Ces paramètres \(\beta\) et \(\gamma\) sont fournis respectivement via les opérandes BETA et GAMMA du mot-clé SCHEMA_TEMPS (SCHEMA=”NEWMARK”) de la commande DYNA_NON_LINE.
Voir le [r5.05.05-qualites-defauts-schema-newmark] pour les caractéristiques du schéma en fonction des valeurs de ces paramètres.
Les crochets aux seconds membres des équations (4900) et (4901) sont ainsi des moyennes pondérées de \(\ddot{\Phi}(t)\) et de \(\ddot{\Phi}(t+\Delta t)\) .
On remarque qu’on ne peut avoir \(\beta =0\) . Au premier pas de calcul on exploite directement les conditions initiales: \(\Phi (0)\) et \(\dot{\Phi}(0)\) . Il faut aussi \(\ddot{\Phi}(0)\) , sauf si \(\gamma =1\) : voir la remarque faite au [r5.05.05-phase-prediction].
En pratique ces expressions ne sont pas utilisables car on va devoir exprimer les valeurs en déplacements à l’instant \(t+\Delta t\) à partir de celles à l’instant \(t\) . Entre ces deux instants, on note l’incrément et la moyenne des champs:
\(\Delta \Phi =\Phi (t+\Delta t)-\Phi (t)\) ; \(\stackrel{ˉ}{\Phi}=\frac{1}{2}(\Phi (t)+\Phi (t+\Delta t))\)
Ces résultats seront utilisés au [r5.05.05-point-vue-energetique]. L’équation (4900) donne:
Et, d’après (4901):
Les expressions (4902) et (4903) constituent la paramétrisation en déplacement, dite «D-form», où on exprime les accélérations et les vitesses avec les déplacements, ce qui sera avantageux en non linéaire pour exprimer les lois de comportement. Enfin, on peut écrire:
Dans le cas des grandes rotations d’éléments de structure les expressions homologues à (4902) et (4903) sont plus complexes [bib5], mais les relations qui suivent sont assez aisément transposables.
Au cours d’un pas de temps (de \(t\) à \(t+\Delta t\) ), où les valeurs à l’instant \(t\) sont connues, les équations (4903) et (4902) définissent les accroissements de vitesse \(\delta \dot{\Phi}\) et d’accélération \(\delta \ddot{\Phi}\) correspondant à un accroissement de déplacement arbitraire \(\delta \Phi\) à partir de la position à l’instant \(t\) , dont on aura besoin lors des itérations de correction de NEWTON (au sein du pas de temps, voir le paragraphe r5.05.05-phase-correction-methode-newton) :
Plaçons-nous à l’instant \(t={t}_{i-1}\) , et on écrit l’équilibre[éq2.2.4-1] après discrétisation spatiale à l’instant \({t}_{i}={t}_{i-1}+\Delta t\) , éventuellement avec les éléments complémentaires apportés en[éq2.2-6]. On note \({U}_{i}=U({t}_{i})\) les degrés de liberté au nouvel instant \({t}_{i}\) et en exploitant les termes du schéma temporel (4903) et (4902), on aboutit au système non linéaire d’équilibre dynamique:
avec:
Remarque :
On constate que la matrice \(\stackrel{ˆ}{K}\) contribuant à la raideur généralisée du système à résoudre est enrichie de termes provenant de la matrice de masse \(M\) du système (à condition que l’on ait affecté une masse volumique sur tous les éléments finis du modèle) et de la matrice d’amortissement \(C\) *. On verra plus loin que la linéarisation des forces internes* \(R({U}_{i},{\dot{U}}_{i},t)\) contribue aussi à \(\stackrel{ˆ}{K}\) . Ainsi, même si le système mécanique considéré comporte des modes rigides, par exemple dans une étude où un solide est en chute libre, la matrice de masse vient «éviter» que la matrice de «raideur» \(\stackrel{ˆ}{K}\) ne soit pas factorisable. En quelque sorte, ce sont les forces d’inertie du corps solide (dans ses modes rigides) qui assurent l’équilibre avec les forces extérieures, ce qui ne pourrait se faire en quasi-statique.
Cependant, on observe que le pas de temps \(\Delta t\) apparaît aussi. S’il est trop grand, les termes de masse ne seront pas assez importants face à ceux de raideur, et la matrice risque d’être quasi non factorisable («pivots quasi-nuls»).
Remarque :
Les termes avec l’exposant \({}^{\text{fs}}\) apparaissent si le système contient des domaines fluides (voir r5.05.05-problemes-couples-fluide-structure-vibro-acoustique et les termes venant du domaine fluide dans (4895)); on rappelle qu’il n’est pas prévu de chargement fluide \({L}^{\mathrm{fl}}\) .
Remarque :
Contrairement au cas d’amortissement de Rayleigh (voir [r5.05.05-terme-inertie]), la matrice d’amortissement \(C\) n’apparaît pas dans la matrice \(\stackrel{ˆ}{K}\) en présence d’amortissement modal (mot-clé AMOR_MODAL), car dans ce cas, la force \(C.{\dot{U}}_{i-1}\) est construite directement sans stockage de la matrice \(C\) , et est reportée au second membre \(\stackrel{ˆ}{L}({t}_{i})\) .
Le terme \({L}^{\mathrm{abso}}\) dans le second membre de (4289) apparaît en présence de frontières absorbantes d’un fluide ou d’un solide (formulation élastique), voir le paragraphe r5.05.05-equations-mouvement. Il est traité selon un schéma explicite en fonction des champs solutions obtenus à \({t}_{i-1}\) .
Remarque :
Contrairement au \(\theta\) –schéma également utilisé en thermique transitoire, cf. r5.02.01, où les dérivées sont écrites de manière explicite, tandis que l’équation de conservation est écrite à un instant fictif résultant d’une combinaison avec le paramètre \(\theta\) des valeurs aux instants \({t}_{i-1}\) et \({t}_{i}\) , cf. r5.05.05-schemas-explicites, les équilibres dynamiques sont vérifiés sur les instants \({t}_{i-1}\) et \({t}_{i}\) , tandis que les dérivées sont des combinaisons déterminées par le schéma de Newmark (4900) et (4901). En revanche, le schéma HHT complet, voir le paragraphe r5.05.05-schema-hht-methode-acceleration-moyenne-modifiee exploite un décalage, à la manière du \(\theta\) –schéma.
Il reste encore à traiter les termes non linéaires: les forces internes \(R({U}_{i},{\dot{U}}_{i},t)\) et la contribution non linéaire des forces d’inertie en grandes rotations de structures \({L}_{\mathrm{GR}}^{\mathrm{iner}}({U}_{i},{\dot{U}}_{i},{\ddot{U}}_{i})\) ( cf. r5.05.05-amortissement).
Au premier instant \({t}_{1}\) , on voit dans le second membre (4289) que l’on a besoin de \({U}_{0},{\dot{U}}_{0}\) , fournis par les conditions initiales, mais aussi de \({\ddot{U}}_{0}\) à cause du schéma de NEWMARK (sauf si \(\gamma =2\beta =1\) ): voir la remarque faite au paragraphe r5.05.05-phase-prediction.
On utilise donc une méthode de Newton pour résoudre ce problème non linéaire (4286).
Phase de prédiction#
Le système (4287) est non linéaire et est intégré, après une prédiction d’EULER par linéarisation, à l’aide d’une méthode itérative de NEWTON, comme en statique non linéaire [r5.03.01]. Le calcul de cette prédiction peut être légèrement erroné, tant que la phase de correction par itérations de NEWTON [r5.05.05-phase-correction-methode-newton] est capable de corriger à convergence…
Pas de temps général#
À la phase de prédiction, on exploite la solution au pas précédent ou les valeurs de l’état initial, et on note la matrice de raideur tangente initiale du pas ( cf. définition de \(Q\) au [r5.05.05-discretisation-dynamique-non-lineaire]):
Ces termes ( cf. r5.05.05-discretisation-dynamique-non-lineaire) sont évalués sur le pas précédent, le premier terme du second membre de (4290) n’apparaissant qu’en grands déplacements (\(Q\) n’étant alors pas constant).
Si le comportement est linéaire, la matrice \({K}_{i-1}\) est simplement la matrice de rigidité élastique \(K\) de la structure.
On peut aussi décider pour économiser du temps calcul de ne pas réactualiser cette matrice (seulement dans le cas où le pas de temps ne varie pas), et de prendre la matrice de rigidité élastique, voir [u4.51.03], mot-clé Newton, opérande PREDICTION, valeur ‘ELASTIQUE’, plutôt que ‘TANGENTE’. Si on ne spécifie rien, le choix par défaut fait par code_aster est cohérent avec celui fait sur les itérations de correction de NEWTON décrites ci-après au paragraphe r5.05.05-phase-correction-methode-newton.
Par ailleurs, en présence de grandes rotations (éléments de structure: poutres…), on doit aussi dériver le terme non linéaire d’inertie \({L}_{\mathrm{GR}}^{\mathrm{iner}}({U}_{i},{\dot{U}}_{i},{\ddot{U}}_{i})\) :
La matrice \(\stackrel{ˆ}{K}\) ayant été établie par (4288), cette nouvelle matrice \({K}_{i-1}^{\mathrm{M1}}\) est combinée à la matrice \(\stackrel{ˆ}{K}+{K}_{i-1}\) , voir ci-dessous (4293).
On note l’accroissement de chargement \(\Delta \stackrel{ˆ}{L}({t}_{i})\) , \(\stackrel{ˆ}{L}\) étant défini avec ref:éq 3.1-10 <éq 3.1-10>]:
On définit aussi les accroissements de conditions imposées \(\Delta {U}_{d}({t}_{i})\) , et on regroupe dans \(\Delta {L}^{\text{anél}}({t}_{i})\) les dépendances des contraintes en fonction des divers paramètres ou «variables de commande» \(Z\) de la loi de comportement du matériau constitutif: tels que la température…:
\(\Delta {L}^{\text{anél}}({t}_{i})=-{}^{t}Q({U}_{i-1}).\frac{d\sigma }{dZ}{\mid}_{({U}_{i-1},{\dot{U}}_{i-1},{\ddot{U}}_{i-1},Z)}.\Delta Z\)
En présence d’amortissement modal les forces d’amortissement sont reportées au second membre. On ajoute alors à \(\Delta \stackrel{ˆ}{L}({t}_{i})\) le terme correspondant: \(-C.{\dot{U}}_{i-1}\) .
Le terme \(\Delta {L}^{\text{abso}}({t}_{i})=-{A}_{\text{fa}}\Delta \dot{P}({t}_{i-1})+\Delta {L}_{\text{abso}}^{\text{ordre}0}({U}_{\text{ent}}({t}_{i}),{\dot{U}}_{\text{ent}}({t}_{i}),\dot{U}({t}_{i-1}))\) désigne la contribution intégrée en explicite à partir des solutions à \({t}_{i-1}\) (pour éviter d’avoir à traiter une matrice non symétrique, cf. [r4.02.05]) d’une frontière fluide absorbante et d’une frontière élastique absorbante, cf. (4895) et (4897).
On calcule alors des valeurs prédictives pour le pas de temps en cours \((\Delta {U}_{i}^{0},\Delta {\lambda}_{i}^{0},\Delta {\mu}_{i}^{0})\) :
Prédiction
où on définit la prédiction \({U}_{i}^{0}={U}_{i-1}+\Delta {U}_{i}^{0}\) pour le nouvel instant \({t}_{i}={t}_{i-1}+\Delta t\) , sous les conditions de contact établies au pas de temps précédent.
Si on a choisi MATRICE= ‘ELASTIQUE’ dans le mot-clé NEWTON, on ne réévalue pas à chaque pas de temps \(\stackrel{ˆ}{K}+{K}_{i-1}+{K}_{i-1}^{\mathrm{M1}}=\stackrel{ˆ}{K}+{K}_{0}+{K}_{0}^{\mathrm{M1}}\) , ce qui évite le coût de ré-assemblage et d’inversion, mais augmente le nombre d’itérations de correction.
Après l’établissement d’un candidat solution de (4293), sur la base des liaisons déjà en contact au pas précédent, sans assurer de nouvelle vérification du critère de contact, on lance l’algorithme de contraintes actives pour satisfaire les conditions de contact: on corrige ainsi \({U}_{i}^{0},{\lambda}_{i}^{0},{\mu}_{i}^{0}\) [r5.03.50].
Remarque: matrice tangente singulière:
On vérifie sur (4293) que si la matrice tangente \({K}_{i-1}\) est singulière (cas d’un mode rigide, d’un matériau endommagé, ou plateau ductile…), la dynamique est pilotée par les forces d’inertie, et que, la matrice \(\stackrel{ˆ}{K}\) étant en général régulière (cf. remarque faite au [r5.05.05-schema-newmark]), l’on trouve malgré tout bien un prédicteur \({U}_{i}^{0}\) , à condition que la précision dans la matrice \(\stackrel{ˆ}{K}\) ne soit pas perdue (pivot quasi-nul) à cause d’un choix de pas de temps trop grand qui joue en \(1/\beta \Delta {t}^{2}\) dans (4288). C’est par exemple le cas en dynamique de chute libre (cf. la remarque de l’équation 3.1-10).
Premier pas de temps#
Si l’on se trouve au premier pas de temps, d’une étude ou d’une reprise (poursuite), le prédicteur est calculé différemment pour tenir compte de l’état initial (\({U}_{0}\) , \({\dot{U}}_{0}\) , \({\sigma}_{0}\) ):
avec:
et l’accélération \({\ddot{\mathrm{U}}}_{0}^{0}\) évaluée par la résolution préalable du système (on simplifie en supposant les liaisons de contact bloquées, les itérations de NEWTON se chargeant de corriger):
Remarque:
Il serait plus exact de calculer: \(\mathrm{M}{\ddot{\mathrm{U}}}_{0}^{0}+{}^{t}B.{\lambda}_{0}+{}^{t}\mathrm{A}.{\mu}_{0}=\mathrm{L}({t}_{0})+{\mathrm{L}}^{\text{anél}}({t}_{0})-{}^{t}\mathrm{Q}.{\sigma}_{0}-\mathrm{C}{\dot{\mathrm{U}}}_{0}\) , mais pour simplifier sachant que cela n’aura que peu d’influence sur la suite des solutions, on néglige \(\mathrm{C}{\dot{\mathrm{U}}}_{0}\) dans le second membre de (4296).
On doit noter que:
la matrice \(M\) doit être inversible: on devra affecter une masse volumique sur tous les éléments finis du modèle,
il peut être nécessaire d’établir par un calcul statique (éventuellement non linéaire) l’état d’équilibre sous le chargement initial, donc les contraintes \({\sigma}_{0}\) , avant l’imposition des conditions dynamiques initiales \({U}_{0},{\dot{U}}_{0}\) . En effet, sinon, l’accélération \({\ddot{U}}_{0}^{0}\) pourrait être «excessive» et conduire sur une branche non désirée d’équilibre,
en présence de «charges cinématiques», cf. [u4.44.03], ces dernières sont mises à zéro à cette étape de calcul de \({\ddot{U}}_{0}^{0}\) .
Cependant, certains éléments finis ne possèdent pas de masse sur tous les degrés de liberté, par exemple les poutres avec gauchissement POU_D_TGjustement sur les degrés de liberté de gauchissement, ou en couplage fluide/structure. Pour les éléments de couplage vibro-acoustiques [r5.05.05-problemes-couples-fluide-structure-vibro-acoustique] en effet, la matrice de masse du problème (4895) n’est pas inversible. Dans ces cas-là, la matrice \(M\) n’est pas inversible, et on se contente de prendre une accélération initiale nulle sur l’ensemble du modèle et la suite des itérations devra se charger de corriger cette prédiction moins «bonne». On choisira alors un petit pas de temps pour assurer la convergence au moins au début du transitoire.
Le terme \({\mathrm{L}}^{\text{abso}}({t}_{0})\) est évalué à partir de la vitesse initiale \({\dot{U}}_{0}\) ( cf. [r4.02.05]).
Phase de correction par la méthode de NEWTON#
On cherche les valeurs \((\Delta {U}_{i},\Delta {\lambda}_{i},\Delta {\mu}_{i})\) des incréments de déplacements et paramètres de LAGRANGE depuis les valeurs \(({U}_{i-1},{\lambda}_{i-1},{\mu}_{i-1})\) obtenues à l’équilibre précédent (instant \({t}_{i-1}\) ). On prend comme valeurs initiales \((\Delta {U}_{i}^{0},\Delta {\lambda}_{i}^{0},\Delta {\mu}_{i}^{0})\) obtenues à l’issue de la phase de prédiction, avant de commencer les itérations de la méthode de NEWTON.
À chaque itération \(n\) de NEWTON, on note par des \(\delta\) les évolutions conduisant à l’estimation des incréments \(\Delta\) (du pas de temps \(i\) au pas de temps \(i+1\) ) à convergence des itérations: \(\Delta {U}_{i}^{n+1}={U}_{i}^{n}+\delta {U}_{i}^{n+1}-{U}_{i-1}\) et \(\Delta {\lambda}_{i}^{n+1}={\lambda}_{i}^{n}+\delta {\lambda}_{i}^{n+1}-{\lambda}_{i-1}\) (de même pour les \(\mu\) ).
On doit résoudre alors un système permettant de déterminer \((\delta {U}_{i}^{n+1},\delta {\lambda}_{i}^{n+1},\delta {\mu}_{i}^{n+1})\) , incréments des déplacements et des paramètres de LAGRANGE depuis le résultat \(({U}_{i}^{n},{\lambda}_{i}^{n},{\mu}_{i}^{n})\) de l’itération précédente:
Correction (itération ° \(n\) )
avec un second membre appelé «résidu», car il tend vers zéro à la convergence. On note ( cf. équilibre (4894)):
Le terme \(\stackrel{ˆ}{L}({t}_{i})\) est défini comme dans (4289); la matrice \(\stackrel{ˆ}{K}\) est donnée par (4288).
En présence d’amortissement modal, les forces d’amortissement sont reportées au second membre. Selon la valeur donnée au mot-clé AMOR_MODAL, REAC_VITE, on ajoute à \(\stackrel{ˆ}{L}({t}_{i})\) le terme réactualisé: \(-C.{\dot{U}}_{i}^{n}\) , ou non réactualisé \(C.{\dot{U}}_{i-1}\) .
La matrice \({K}_{i}^{n}\) est la matrice de l’application linéaire tangente de la partie «forces internes» du système d’équations non linéaires (4287); elle vaut donc:
En l’absence de forces suiveuses, le dernier terme est nul. Les forces suiveuses peuvent être: la pression exercée sur les bords d’éléments massifs, le chargement de pesanteur pour les éléments de câble, la force centrifuge en grands déplacements, le chargement de pesanteur pour toutes les modélisations THM des milieux poreux non saturés [r7.01.10].
Si on envisage la réactualisation de la géométrie (en grands déplacements), on a plus précisément:
Le premier terme est la contribution du comportement comme en petites transformations, à la différence que cette contribution est ici évaluée en configuration actuelle. Le second terme est la contribution de la géométrie qui n’est pas présente en petites transformations. Dans le cadre de la réactualisation PETIT_REAC, ce terme n’est pas présent dans le calcul de la matrice tangente. Il est pris en compte pour les autres options de non linéarité géométrique (GREEN et SIMO_MIEHE)
La matrice \({K}_{i}^{\mathrm{Mn}}\) est la matrice de l’application linéaire tangente de la partie «forces d’inertie» du système d’équations non linéaires (4287) qui vaut donc:
En pratique on peut utiliser la «vraie» matrice tangente \({K}_{i}^{n}\) , mais cela présente un coût calcul certain (calcul et factorisation), ou se contenter d’une réactualisation de temps en temps: voir en [u4.51.03] le mot-clé facteur Newton, mot-clé MATRICE.
Remarque: matrice tangente singulière:
On vérifie sur (4297) que si la matrice tangente \({K}_{i}^{n}\) est singulière (cas d’un matériau endommagé, ou plateau ductile…), la dynamique est pilotée par les forces d’inertie, et que, la matrice \(\stackrel{ˆ}{K}\) étant régulière, l’on trouve malgré tout bien un correcteur \(\delta {U}_{i}^{n+1}\) (comme si on était par exemple en situation de chute libre). Cependant, si le pas de temps est «trop élevé», la matrice \((\stackrel{ˆ}{K}+{K}_{i}^{n}+{K}_{i}^{\mathrm{Mn}})\) peut être mal conditionnée et le solveur lui trouve un pivot quasi-nul.
Après chaque itération de Newton ayant établi un candidat solution de (4297) sans vérifier le critère de contact, on lance l’algorithme de contraintes actives pour satisfaire les conditions de contact: on corrige ainsi \(\delta {U}_{i}^{n+1}\) . Pour les autres algorithmes de contact-frottement de Code_Aster , tel que celui utilisé par les éléments discrets «à choc» (DIS_CHOC), cf. [bib17], l’itération de Newton (4297) traite toutes les non-linéarités en même temps.
Mise à jour#
Dans le cas des petites rotations, pour les modélisations habituelles (éléments massifs, poutres, plaques, coques, éléments discrets…), la mise à jour à la fin de l’itération \(n+1\) s’appuie sur les formules (4905) et (4286):
Dans le cas des grandes rotations des éléments de structure (poutres…) la mise à jour, nettement plus complexe, est indiquée en [r5.03.40].
Critère d’arrêt#
Le critère de convergence globale de l’algorithme de NEWTON est identique à celui pratiqué dans STAT_NON_LINE, cf. [bib19]. Il traduit la vérification de l’équilibre dynamique.
À l’instant \({t}_{i}\) , on arrête les itérations au rang \(n\) , dès que l’inégalité suivante est satisfaite:
\(\eta\) est une tolérance, introduite en donnée par l’utilisateur (mot-clé CONVERGENCE, opérande RESI_GLOB_RELA), valant par défaut 10–6, et \({\parallel .\parallel }_{\infty}\) est la norme du maximum sur les degrés de liberté.
Le dénominateur de (4303) est une norme du chargement à l’instant \({t}_{i}\) , à laquelle on rapporte le numérateur, qui est une norme des forces non (encore) équilibrées.
Remarque :
Il comporte une contribution des charges anélastiques, utiles pour normaliser l’erreur sur l’équilibre en absence de charge extérieure. Il doit contenir les forces internes dues aux variables de commande. On y ajoute de plus les forces d’inertie, afin que dans les cas rares où aucune autre force n’agit, le critère reste calculable.
Tout comme en statique, on doit accorder de l’importance à une convergence correcte car sinon, les estimations de forces intérieures et de réactions de contact, vérifiant les relations de comportement et de liaison de contact-frottement s’éloignent de l’équilibre, et ce n’est pas une force d’inertie \(M({U}_{i}^{n}).{\ddot{U}}_{i}^{n}\) trop éloignée de la valeur «exacte» qui va suffire à produire au pas de temps suivant une bonne réponse dynamique transitoire.
On pourra aussi employer d’autres critères de convergence comme proposé dans STAT_NON_LINE [u4.51.03].
Qualités et défauts du schéma de NEWMARK#
Propriétés du schéma de NEWMARK#
Ce paragraphe et le suivant reprennent partiellement certaines parties de [bib2]. On pourra aussi consulter [bib22], [bib25] et [bib30].
Généralités#
On définit à l’aide des méthodes d’analyse numérique plusieurs types de propriétés pour un schéma d’intégration temporelle. Voici leur signification:
Convergence: la solution tend vers une limite quand le pas de temps \(\Delta t\) tend vers0;
Précision: taux de convergence quand le pas de temps \(\Delta t\) tend vers 0;
Consistance et ordre du schéma: le résidu est borné par \({(\Delta t)}^{k+1}\) (exemple \(k=\mathrm{2 }\) : pour la règle du trapèze): un polynôme d’ordre \(k+1\) est intégré exactement;
Stabilité: une perturbation finie de l’état initial n’entraîne pas de perturbation exagérément amplifiée («explosion numérique») à l’état ultérieur; il faut que l’amplification (rayon spectral) soit inférieure à 1: pour un schéma de forme générale \({U}_{i}=A{U}_{i-1}+{L}_{i}\) , on doit avoir \(\rho (A)\le 1\) . C’est une condition nécessaire pour ne pas diverger! Si l’amplification \(\rho (A)<1\) il y a atténuation numérique.
L’erreur en temps \(e\) induite par le schéma de NEWMARK (4900), (4901), qui approche un développement de Taylor de la solution, est donnée par, cf. [bib26]:
que l’on peut normaliser par le vecteur position \(\parallel X\parallel\) ou l’amplitude \(\parallel U\parallel\) . On peut aussi remplacer dans (4758) la norme \({L}_{2}\) (\(\parallel \Delta {\ddot{U}}_{i}\parallel\) ) par la norme \({L}_{\infty}\) : \({\mid \Delta {\ddot{U}}_{i}\mid }_{\infty}=\underset{\mathrm{noeuds}k}{\max}(\mid \Delta {\ddot{U}}_{i}^{k}\mid )\) . Code_Aster ne fournit pas, actuellement, l’information de cet estimateur.
L’analyse numérique du schéma de NEWMARK se fait sur le traitement de l’équation de l’oscillateur linéaire à un degré de liberté, sans phénomène dissipatif:
Les schémas temporels ne donnent pas la solution exacte de (4582):
\(x={C}_{1}\cos\omega t+{C}_{2}\sin\omega t\)
avec \(\omega =\sqrt{k/m}\) , \({C}_{1}\) et \({C}_{2}\) dépendant des conditions initiales, mais une solution approchée :
\({x}_{i}=\left[\exp(-\xi \tilde{\omega}{t}_{i})\right]({C}_{1}^{'}\cos\tilde{\omega}{t}_{i}+{C}_{2}^{'}\sin\tilde{\omega}{t}_{i})\)
Deux erreurs sont donc apportées:
d’une part, il s’introduit un amortissement artificiel défini par le taux d’amortissement réduit \(\xi\) , qui est le décrément logarithmique divisé par \(2\pi\) ;
d’autre part, la pulsation \(\omega\) correspondant à la période exacte \(T\) est remplacée par la pulsation \(\tilde{\omega}\) correspondant à une période \(\tilde{T}\) .
Ces erreurs dépendent du rapport \(\frac{\Delta t}{T}\) et du schéma lui-même.
En écrivant l’équilibre du système à 1 degré de liberté (4582), sans effort imposé, et les équations du schéma de NEWMARK (4900), (4901) on obtient:
d’où en terme d’accroissements:
Les propriétés de cette matrice d’amplification –qui doit conserver le caractère harmonique de la réponse, plus précisément celles de son polynôme caractéristique, - servent à caractériser les propriétés du schéma, en régime linéaire, [bib25].
Les propriétés du schéma de NEWMARK sont résumées ci-après.
À cause de la manière choisie d’exprimer le schéma ( cf. (4900) et (4901)), on ne peut pas prendre \(\beta =0\) .
Pour les problèmes linéaires , le schéma est inconditionnellement stable , même en présence d’amortissement physique –-une perturbation n’est pas amplifiée par le schéma– c’est-à-dire stable quelle que soit la taille du pas de temps, si les paramètres satisfont aux inégalités:
Si \(2\gamma >1\) et \(\beta <\frac{1}{4}{(\gamma +\frac{1}{2})}^{2}\) le schéma est conditionnellement stable : le pas de temps doit être choisi inférieur à: \(\Delta t\le \frac{{T}_{\min}}{\pi}{({(\gamma +\frac{1}{2})}^{2}-4\beta )}^{-1/2}\) -cette relation étant valable en absence d’amortissement physique \(\xi =0\) mais aussi en présence d’amortissement physique–- en fonction de la plus petite période \({T}_{\min}\) de vibration du système étudié, voir [bib22], [bib25].
En présence d’amortissement physique \(\xi >0\) , on peut se permettre des \(\Delta t\) légèrement supérieurs.
En présence de contact avec impact, la stabilité du schéma est assurée si \(2\beta \ge \gamma \ge \mathrm{½}\) , [bib31].
Le tableau suivant montre quelques cas particuliers:
\(\beta\) |
\(\gamma\) |
méthode |
type |
propriétés |
\(\frac{1}{12}\) |
\(\frac{1}{2}\) |
«Fox-Goodwin» |
implicite |
Schéma d’ordre 4 en temps; pas de dissipation numérique en absence d’amortissement matériau: donc pas d’atténuation d’amplitude due au schéma. Conditionnellement stable: \(\Delta t\le \sqrt{3/2}/{(\pi {f}_{\max})}^{-1}\) , \({f}_{\max}\) étant la fréquence vibratoire maximale «visée» dans la simulation, pas de dissipation numérique. Erreur en temps: \(1+{(\omega .\Delta t)}^{4}/480+...\) |
\(\frac{1}{6}\) |
\(\frac{1}{2}\) |
accélération «linéaire» |
implicite identique à la \(\theta\) -méthode WILSON, avec \(\theta =1\) |
Schéma d’ordre 2 en temps, pas de dissipation numérique en absence d’amortissement matériau: donc pas d’atténuation d’amplitude due au schéma. Conditionnellement stable: \(\Delta t\le \sqrt{3}/{\left(\pi {f}_{\max}\right)}^{-1}\) , \({f}_{\max}\) étant la fréquence vibratoire maximale «visée» dans la simulation, pas de dissipation numérique Erreur en temps: \(1-{(\omega .\Delta t)}^{2}/24+...\) |
\(\frac{1}{4}\) |
\(\frac{1}{2}\) |
«règle du trapèze» ou accélération moyenne |
implicite |
Schéma d’ordre 2 en temps, pas de dissipation numérique en absence d’amortissement matériau: donc pas d’atténuation d’amplitude due au schéma. Stabilité inconditionnelle en \(\mathrm{\Delta t}\) . Erreur en temps: \(1-{(\omega .\Delta t)}^{2}/12+...\) : les fréquences sont décalées vers le bas, mais une matrice de masse consistante limite ce défaut (puisque provoque l’effet inverse) |
\(\beta\) |
\(\gamma\) |
\(\alpha\) -méthode «accélération moyenne modifiée» |
Implicite
voir [ |
Schéma d’ordre 2 en temps
Erreur en temps: \(1-{(\omega .\Delta t)}^{2}(3{\alpha}^{2}+1)/12+...\)
Dissipation numérique: \(|\alpha |\omega .\Delta t/2+...\) voir [ |
Le schéma de NEWMARK est du second ordre en temps (au pire) si et seulement si \(\gamma =1/2\) . Dès que \(\gamma >1/2\) , le schéma de NEWMARK est d’ordre 1, et introduit une dissipation numérique proportionnelle à \((\gamma -\frac{1}{2})\Delta t\) . Afin d’avoir un amortissement numérique croissant avec la fréquence, il convient de choisir: \(\beta \ge \frac{{(\gamma +1/2)}^{2}}{4}\) , l’égalité étant le choix optimal, cf. le schéma HHT, voir le [r5.05.05-variante-schema-newmark].
Si on choisissait \(\gamma <1/2\) , le schéma de NEWMARK apporterait un amortissement numérique négatif qui amènerait une instabilité.
En utilisant les accroissements (4584), on déduit l’erreur en temps (4758) du schéma, moyennée sur une période \(2\pi /\omega\) , et normalisée par rapport à l’amplitude, pour l’oscillateur à 1d.d.l., cf. [bib26]:
Cette expression peut servir d’erreur de référence au schéma appliqué à la réponse dynamique d’une structure quelconque.
Remarque :
Les propriétés du schéma de Newmark obtenues sur le cas de l’oscillateur à 1 degré de liberté sont généralisables au cas du système dynamique linéaire à plusieurs degrés de liberté (en utilisant une projection sur base modale) pourvu que l’opérateur d’amortissement soit proportionnel aux opérateurs de masse et de raideur réduits.
Remarque complémentaire :
Le schéma de Newmark est connu pour introduire une erreur notable dès qu’une rotation non infinitésimale intervient dans la réponse dynamique du solide. Par ailleurs, le schéma de Newmark peut engendrer aussi une erreur si l’on intervient trop «brutalement» sur des variations de paramètres (pas de temps, \(\beta\) ou \(\gamma\) ) dans la succession de la discrétisation temporelle sur la durée du transitoire.
Schéma d’accélération moyenne ou règle du trapèze#
La «règle du trapèze» (\(\beta =1/4\) , \(\gamma =1/2\) ), inconditionnellement stable , est la plus communément adoptée, associée à une masse consistante (option par défaut MASS_DIAG = “NON”). On montre, par exemple dans [bib30], que la distorsion en fréquence relative apportée par le schéma d’intégration, pour un système amorti (l’amortissement réduit est noté \(\xi =\frac{c}{2\sqrt{\text{km}}}=\frac{c}{\mathrm{2m}\omega }\) ) à 1 degré de liberté en vibration libre, de pulsation propre amortie \({\omega}_{d}=\omega \sqrt{1-{\xi}^{2}}\) :
est:
où \({\omega}_{d}^{\text{Newmark}}\) désigne la valeur approchée de \({\omega}_{d}\) par le schéma (via une analyse spectrale de la réponse).
En vibrations forcées, la dépendance de (3854) en \({\omega}^{2}\Delta {t}^{2}\) reste valable. De manière générale, les fréquences sont décalées vers le bas, et ce décalage croît si le pas de temps augmente et décroît pour un amortissement croissant.
Si l’on souhaite utiliser un schéma où les fréquences sont décalées vers le haut, il convient de choisir l’option MASS_DIAG = “OUI”.
On montre aussi dans [bib30], que la distorsion en amortissement relative apportée par le schéma d’intégration de la «règle du trapèze» n’est pas nulle pour un système comportant un amortissement matériau. Pour de petits pas de temps:
où \({\xi}^{\mathit{Newmark}}\) désigne la valeur de \(\xi\) résultant du schéma.
L’amortissement total diminue, et décroît quand le pas de temps augmente. En vibrations forcées, la dépendance de distorsion en amortissement (3855) devient en \(\omega \Delta t\) , cf. [bib30].
Pour des pas de temps plus importants proches de la pulsation de coupure numérique due au schéma, telle que \({\stackrel{ˆ}{\omega}}_{d}^{\text{Newmark}}\Delta t=\pi\) , cf. [bib30], on obtient l’estimation (cas de vibrations libres):
On peut résumer donc par le fait que la précision est bonne si on choisit \(\Delta t=\frac{1}{10\omega }\) ( cf. (3854)), pour la pulsation \(\omega\) du système la plus élevée que l’on cherche à décrire, cf. [r5.05.05-choix-pas-temps].
Point de vue énergétique#
Dans le cas d’un système conservatif élastique linéaire en petites transformations sans amortissement physique, on peut facilement évaluer les variations des différentes énergies au cours du pas de temps \(\Delta t\) entre \({t}_{i}\) et \({t}_{i+1}\) . Entre ces deux instants, on note l’incrément et la moyenne des champs:
\(\Delta \Phi =\Phi ({t}_{i+1})-\Phi ({t}_{i})\) \(\stackrel{ˉ}{\Phi}=\frac{1}{2}(\Phi ({t}_{i+1})+\Phi ({t}_{i}))\)
Ainsi on a respectivement pour l’énergie cinétique, l’énergie de déformation ; \(M\) et \(K\) étant symétriques, et le potentiel des efforts extérieurs indépendant du temps –le système étant conservatif:
\(\Delta {E}_{\text{cin}}=\frac{1}{2}\Delta \left({}^{t}\dot{\Phi}.\mathrm{M}.\dot{\Phi}\right)={\Delta}^{t}\dot{\Phi}\phantom{\rule{1em}{0ex}}.\mathrm{M}.\overline{\dot{\Phi}}\)
\(\Delta {E}_{\text{déf}}=\frac{1}{2}\Delta \left({}^{t}\Phi .\mathrm{K}.\Phi \right)={\Delta}^{t}\Phi \phantom{\rule{1em}{0ex}}.\mathrm{K}.\overline{\Phi}\)
\(\Delta {V}_{\mathit{ext}}=-\Delta \left({}^{t}\Phi .\mathrm{F}\right)=-{\Delta}^{t}\Phi \phantom{\rule{1em}{0ex}}.\overline{\mathrm{F}}\)
En utilisant l’équilibre aux instants \({t}_{i}\) et \({t}_{i+1}\) . sous l’action des efforts extérieurs, on obtient:
Les variations des diverses énergies, dont l’énergie totale du système –ou hamiltonien: \(\Delta {E}_{\text{tot}}=\Delta {E}_{\text{cin}}+\Delta {E}_{\text{déf}}+\Delta {V}_{\text{ext}}\) , s’expriment, en exploitant (3857) pour éliminer le chargement \(F\) et les termes du schéma (4903) et (4904):
Et on obtient que la variation d’énergie totale du système s’exprime:
On vérifie que lorsque \(\gamma =1/2\) et \(\beta =1/4\) (règle du trapèze), on n’a pas de dissipation numérique d’énergie apportée par le schéma. On remarque qu’un choix différent de \(\gamma\) et/ou \(2\beta \ne \gamma\) amène une dissipation proportionnelle au pas de temps.
On vérifie aussi que si \(\gamma >\phantom{\rule{0.5em}{0ex}}1/2\) et \(2\beta -\gamma >0\) , l’amortissement numérique est positif, cf. [bib25], ce qui est conforme aussi à ce qui a été énoncé au [r5.05.05-generalites].
Propriétés du schéma de NEWMARK pour les problèmes non linéaires#
Pour les problèmes non linéaires (grandes déformations, non-linéarités matériau), le schéma est inconditionnellement stable , si \(\beta \ge \frac{1}{4}\) et avec \(\gamma =1/2\) , cf. preuve en [bib23].
En présence de contact unilatéral , on trouve dans la littérature, cf. [bib24], le conseil de prendre \(\beta =\gamma\) (par exemple avec \(\beta =\gamma =1/2\) ) ce qui assure la compatibilité des vitesses durant la phase où le contact est maintenu entre deux solides.
Cela s’obtient directement à partir des équations (4900) et (4901); le saut \([\text{ . }]\) de vitesse entre les deux points restant en contact sur le pas \(\Delta t\) (\(\left[\Phi \left(t\phantom{\rule{0.5em}{0ex}}\right)\right]=0\) et \(\left[\Phi \left(t\phantom{\rule{0.5em}{0ex}}+\phantom{\rule{0.5em}{0ex}}\Delta t\right)\right]=0\) ) à l’instant \(t+\Delta t\) est en effet:
\(\left[\dot{\Phi}(t+\Delta t)\right]=(1-\gamma /\beta ).\left[\dot{\Phi}(t)\right]+(1-\gamma /(2\beta )).\Delta t.\left[\ddot{\Phi}(t)\right]\)
Si on choisit \(\beta =\gamma\) , alors \([\dot{\Phi}(t+\Delta t)]=o(\Delta t)\)
.Cependant ce choix n’est pas compatible avec celui d’un optimum en terme de dissipation numérique, tel que défini par la \(\alpha\) -méthode ( cf. [r5.05.05-variante-schema-newmark]): il faudrait prendre \(\alpha =-1-\sqrt{2}\) qui donne un amortissement numérique énorme! On n’incite donc pas l’utilisateur à suivre cette recommandation.
Si le contact unilatéral est résolu par une méthode de pénalisation, on n’échappe pas au risque d’instabilité même en introduisant un amortissement numérique.
Choix des pas de temps#
Le pas de temps à choisir doit respecter un certain nombre de conditions.
La première, évidente, est qu’il doit être adapté à l’échantillonnage temporeldes chargements appliqués au système étudié. Accessoirement, il peut être opportun de reconsidérer une modélisation» avec une dépendance temporelle trop «violente des chargements appliqués, en adoucissant des ruptures de pente par exemple.
On conseille, pour des raisons de précision (critère de type «Shannon» sur la fréquence de coupure), un pas de temps à choisir tel que:
où \({T}_{\min}\) désigne la plus petite période de vibration du système que l’on souhaite étudier.
Le pas d’espace des mailles éléments finis choisi intervient aussi: le pas de temps \(\Delta t\) maximal à choisir est de l’ordre de \(h/c\) , où \(c\) est la célérité des ondes de compression élastiques du matériau et \(h\) une taille caractéristique des mailles, si l’on cherche à décrire partiellement des phénomènes en haute fréquence, pour lesquels cependant cette formulation numérique de l’élastodynamique n’est pas véritablement adaptée (il existe d’autres méthodes numériques pour ce faire).
Enfin dans le cas de solides présentant des modes rigides (chute libre par exemple), conformément à la remarque faite au [r5.05.05-schema-newmark], il convient de choisir un pas de temps \(\Delta t\) suffisamment petit pour que les termes de masse soient du même ordre que ceux de raideur (dans la matrice \(\stackrel{ˆ}{K}\) , cf. (4288)). Ainsi, on pourra choisir:
où \(L\) est le diamètre du solide considéré en «chute libre» sous l’accélération de pesanteur \(g\) . Ainsi l’incrément de déplacement \(g\Delta {t}^{2}/2\) lors de ce pas de temps est faible devant les dimensions du solide, ou, d’une autre manière, similaire à un déplacement élastique sous un champ d’action de même ampleur.
Il est possible de faire se succéder plusieurs analyses dynamiques, sur des intervalles de temps successifs, communiquant par des reprises en choisissant comme état initial le résultat du dernier pas étudié auparavant, en utilisant des pas de temps très distincts selon l’idée a priori que l’on a de la réponse du système étudié. On ne possède pas de résultat général avec ce type de choix en terme de convergence…
On sait aussi que le traitement des collisions est sensible au cadencement des pas de temps, par rapport aux instants «réels» de choc: on devra étudier la sensibilité de la réponse obtenue.
Cependant, un pas de temps «trop petit» peut exacerber les oscillations induites par la discontinuité.
On pourra aussi s’orienter vers le choix d’un schéma explicite, d’ordre moins élevé (ordre 1), cf. par exemple [bib31] même si un tel choix conduit à prendre des pas de temps très faibles.
Il est cependant conseillé de maintenir constant le pas de temps durant une phase de réponse dynamique linéaire stationnaire, pour garder les propriétés énoncées précédemment.
Remarques :
On doit noter que Code_Aster ne propose pas aujourd’hui de méthode multi-domaine en temps et espace, qui permettrait de définir un pas de temps par zone dans le solide étudié, ni de critère d’erreur en dynamique.
Dans l’opérateur DYNA_NON_LINE, on peut choisir de résoudre l’équilibre en déplacement ou en vitesse pour les schémas implicites (pour les schémas explicite c’est en accélération). Le mot-clé suivant permet ce choix: FORMULATION = “DEPLACEMENT” ou “VITESSE”.
Une variante du schéma de Newmark : les schémas HHT et d’accélération moyenne modifiée#
Motivation#
En analyse mécanique, on souhaite que les basses fréquences soient reproduites le plus fidèlement possible.
On souhaite par contre que les hautes fréquences soient atténuées par le calcul parce qu’elles peuvent engendrer des instabilités numériques et que les contraintes mécaniques associées sont en général faibles.
Les courbes donnant le facteur d’amortissement \(\xi\) en fonction de \(\frac{\Delta t}{T}\) (période \(T=2\pi /\omega\) ) doivent donc:
partir de l’origine avec une tangente horizontale pour donner un très faible amortissement aux basses fréquences,
être des fonctions croissantes pour atténuer les hautes fréquences et ce d’autant plus qu’elles sont plus élevées.
Pour tenter d’atteindre ces objectifs, Hilber, Hughes et Taylor (HHT) ont proposé dans [bib4] de définir les paramètres \(\beta\) et \(\gamma\) de Newmark en fonction d’un troisième paramètre \(\alpha\) négatif par les relations suivantes, qui sont calquées sur les conditions de stabilité (4901):
Ce choix offre le meilleur compromis sur la précision et l’amortissement en hautes fréquences.
Ce paramètre \(\alpha\), négatif, est fourni via l’opérande ALPHA du mot-clé SCHEMA_TEMPS (SCHEMA=”HHT”) de [u4.53.01].
Schéma HHT et méthode d’accélération moyenne modifiée#
On obtient ainsi [bib4], [bib25]:
Par ailleurs, l’équilibre dynamique (4894), discrétisé en temps à l’instant \({t}_{i}={t}_{i-1}+\Delta t\) est modifié en introduisant un «décalage», contrôlé aussi par le coefficient \(\alpha \le 0\) , sur les forces intérieures:
où les forces extérieures sont évaluées à l’instant \({t}_{i+\alpha }=(1+\alpha ){t}_{i}-\alpha {t}_{i-1}={t}_{i}+\alpha \Delta t\) :
Le système d’équations non linéaires (4287), se récrit donc:
avec, en suivant (4288) et (4289) et \(\beta\) et \(\gamma\) fonction du paramètre \(\alpha\) , cf. (4758):
Les modifications que ce schéma apporte aux phases de prédiction et correction de la méthode de Newton [r5.05.05-phase-prediction] et [r5.05.05-phase-correction-methode-newton] sont les suivantes:
phase de prédiction (4293): mettre la nouvelle expression de \(\stackrel{ˆ}{K}\) (4592), ainsi que \(\left(1+\alpha \right)\,{\mathrm{K}}_{i-1}\) à la place de \({\mathrm{K}}_{i-1}\) dans (4293) et \(\left(\alpha -1\right)\,\mathrm{R}\left({\mathrm{U}}_{i-1},{\dot{\mathrm{U}}}_{i-1},{t}_{i-1}\right)\) à la place de \(-\mathrm{R}\left({\mathrm{U}}_{i-1},{\dot{\mathrm{U}}}_{i-1},{t}_{i-1}\right)\) dans l’expression de \(\Delta \widehat{\mathrm{L}}\left({t}_{i}\right)\) dans (4292);
phase de correction (4297): mettre \(\left(1+\alpha \right)\,{\mathrm{K}}_{i}^{n}\) et \(\alpha {\mathrm{F}}_{i}^{n}\) à la place de \({\mathrm{K}}_{i}^{n}\) et \({\mathrm{F}}_{i}^{\mathrm{n}}\).
- Remarque
L’organisation algorithmique globale est donc inchangée par rapport au schéma à accélération moyenne modifiée: seuls la matrice du système et le second membre sont modifiés, et il convient de stocker les vecteurs obtenus au pas précédents.
Quand \(\alpha =-1\) (d’où \(\beta =1\) , \(\gamma =3/2\) ), le schéma HHT devient alors explicite ( cf. (4590)):
avec:
et
assurant en même temps la plus grande dissipation possible en hautes fréquences.
- Remarque
On peut utiliser les deux variantes suivantes avec le mot-clé facteur SCHEMA_TEMPS(SCHEMA=”HHT”):
- D’une part, en choisissant pour le mot-clé simple MODI_EQUI=’NON’, on utilise l’équation d’équilibre classique où ni les
termes de (4590) à (4593) «décalés» par \(\left(1+\alpha \right)\) ou \(\alpha\) sur l’amortissement et les raideurs, ni l’évaluation à \({t}_{i+\alpha }\) des seconds membres ne sont traités, ce qui fait perdre un ordre sur le schéma (de 2 à 1). Il s’agit donc en réalité de ce qu’on note dans la littérature la «méthode d’accélération moyenne modifiée», qui se limite donc à définir la relation optimale entre les paramètres du schéma de Newmark selon les (4758) à (4583);
D’autre part, en choisissant pour le mot-clé simple MODI_EQUI=’OUI’ (valeur par défaut du code), cela revient à adopter le schéma complet HHT présenté ci-dessus avec décalage des termes d’efforts intérieurs et extérieurs.
Avec le schéma HHT complet (HHT avec MODI_EQUI=’OUI’), on traite les termes venant des déplacements et des forces suiveuses imposés comme dans le cas du schéma d’accélération moyenne (HHT avec MODI_EQUI=’NON’).
Propriétés du schéma d’accélération moyenne modifiée#
Il faut:
pour que les conditions de stabilité du schéma soient remplies, cf. [r5.05.05-proprietes-schema-newmark-problemes-non-lineaires]: le schéma est inconditionnellement stable .
Le schéma «d’accélération moyenne modifiée» (appelé aussi \(\alpha\) -méthode), amène donc \(\beta ={(\gamma +\frac{1}{2})}^{2}/4\) , ce qui est le choix optimal pour apporter de l’amortissement croissant sur les hautes fréquences.
Quand \(\alpha =0\), le schéma «d’accélération moyenne modifiée» redevient la «règle du trapèze» et l’amortissement est nul [Fig. 198].
La valeur \(\alpha =-1\) , soit \(\gamma =\frac{3}{2}\) et \(\beta =1\) produit la plus forte dissipation en haute fréquence, mais détruit beaucoup la précision sur les modes de basse fréquence.
En pratique, dans le schéma «d’accélération moyenne modifiée» originel, on limite à \(\alpha \in \left[-\frac{1}{3},0\right]\) , ce qui assure la monotonie de l’accroissement de l’amortissement en fonction de la fréquence. Le choix \(\alpha =-0.10\) semble être efficace, même s’il est possible de prendre de valeurs de \(\alpha\) plus élevées en valeur absolue.
La figure Fig. 198, extraite de [bib5], donne les variations de l’amortissement numérique \(\xi\) en fonction de \(\frac{\Delta t}{T}\) pour quelques valeurs de \(\alpha\) , où \(T\) est la période du système à 1 d.d.l. Cette figure appelle les commentaires suivants:
la «règle du trapèze» (\(\alpha =0\) ) est séduisante parce qu’elle n’apporte pas d’amortissement numérique pour un système non amorti physiquement, cf. [
r5.05.05-schema-acceleration-moyenne], mais elle peut être instable en non linéaire,quand \(\alpha \ne 0\) , les courbes ne sont pas à tangente horizontale à l’origine. Comme le schéma «d’accélération moyenne modifiée» (\(\alpha\) -méthode) apporte un important amortissement artificiel, certains utilisateurs [bib5] choisissent un pas de temps, puis déterminent la valeur du paramètre \(\alpha\) pour que, dans la plage de fréquences qui les intéresse, le taux d’amortissement numérique soit du même ordre que le taux d’amortissement mécanique réel qui, alors, n’est pas pris en compte,
l’erreur en temps ne dépend pas de \(\alpha\) , cf. [bib30].
Fig. 198 Taux d’amortissement numérique en fonction du pas de temps du schéma d’accélération moyenne modifiée#
On constate que l’on a pour les petits pas de temps \(\frac{\Delta t}{T}\) l’amortissement numérique suivant:
ce qui permet d’estimer dans la gamme de fréquence visée l’amortissement numérique moyen apporté par le schéma d’accélération moyenne modifiée (\(\alpha\) -méthode), qui vient s’ajouter à l’amortissement physique éventuellement déjà présent.
La figure Fig. 199, également extraite de [bib4], donne les variations de l’erreur relative en période en fonction de \(\frac{\Delta t}{T}\) . L’erreur est une fonction croissante de \(\vert \alpha \vert\) .
Fig. 199 Erreur relative sur la période en fonction du pas de temps du schéma d’accélération moyenne modifiée#
Propriétés du schéma HHT#
Il faut:
pour que les conditions de stabilité du schéma soient remplies, cf. [r5.05.05-proprietes-schema-newmark-problemes-non-lineaires]: le schéma est inconditionnellement stable .
Quand \(\alpha =0\) , le schéma HHT redevient la «règle du trapèze» et l’amortissement numérique est nul [Fig. 198]. La valeur \(\alpha =-1\) conduit à écrire l’équilibre au pas précédent et est donc déconseillée.
En pratique, dans le schéma HHT, on se limite à \(\alpha \in \left[-\frac{1}{3},0\right]\) , ce qui assure la monotonie de l’accroissement de l’amortissement en fonction de la fréquence. Le choix \(\alpha =-0.10\) semble être efficace. L’erreur en temps croît moins vite que l’amortissement numérique en fonction de \(\alpha\) [bib5].
La figure [Fig. 199], extraite de [bib25], donne les variations de l’amortissement numérique \(\xi\) en fonction de \(\frac{\Delta t}{T}\) pour quelques valeurs de \(\alpha\) , en les comparant avec le schéma d’accélération moyenne modifiée. On constate que:
la méthode HHT amortit moins les basses fréquences que le schéma d’accélération moyenne modifiée;
de plus, pour une même valeur du paramètre \(\alpha\) , l’amortissement numérique est nettement plus important, en « moyenne » et « haute » fréquence \((\omega \Delta t>1)\) dans le cas des schémas de type accélération moyenne modifiée, que pour le schéma HHT. Par exemple, pour \(\omega \Delta t=2\) , et pour \(\alpha =-0.05\) , on a un rapport 2 sur les valeurs d’amortissement numérique;
on conseille donc avec le schéma HHT d’envisager la prise en compte d’un amortissement physique [
r5.05.05-terme-inertie], qu’on aurait pu ne pas prendre en compte en ayant choisi le schéma d’accélération moyenne modifiée.
Fig. 200 Taux d’amortissement numérique en fonction du pas de temps du schéma HHT, comparaison avec celui produit par le schéma d’accélération moyenne modifiée, cf. [bib29]#
Schéma HHT et overshooting#
Comme mentionné dans le paragraphe r5.03.01-proprietes-schema-hht, le schéma HHT initialement introduit dans [bib4] induit un amortissement dans les hautes fréquences, tout en présentant un amortissement limité dans les basses fréquences. Ceci permet notamment d’éliminer des modes associés à des fréquences élevées, pour lesquels la discrétisation spatiale et/ou temporelle n’est pas suffisante.
Toutefois, certaines composantes de la solution numérique peuvent être faiblement amorties (notamment dans le cas de systèmes raides), ce qui induit l’existence d’oscillations altérant la solution, conduisant au phénomène connu sous le nom d”overshooting. En particulier, Hulbert et Hughes ont démontré l’existence d’overshooting en vitesse dans le cas du schéma HHT (voir [bib37]). De plus, Erlicher et al. ont montré que l’entièreté de la famille des schémas \(G-\alpha\) (dont le schéma HHT fait partie) peut présenter de l’overshooting (voir [bib36]). Par ailleurs, Bauchau et al. [bib38] proposent également une discussion relative aux effets de l’overshooting sur les réponses obtenues avec le schéma HHT.
Ainsi, en dépit de ses propriétés d’amortissement numérique appréciables, le schéma HHT peut être sujet à des phénomènes d’overshooting. Dans ce contexte, Kaiping a proposé dans [bib35] une modification concernant l’ensemble des schémas de type \(G-\alpha\), dans l’optique de limiter les effets d’overshooting. En particulier, une nouvelle version du schéma HHT, intitulée NOHHT (No Overshooting HHT), a pu être proposée, via une modification des coefficients du schéma. On renvoie à [bib35] pour de plus amples détails relatifs à cette modification.
Le schéma NOHHT de [bib35] peut être utilisé dans [u4.53.01], via le mot-clé SCHEMA_TEMPS.
Dans l’optique de comparer les schémas HHT et NOHHT, on considère le cas élémentaire (introduit dans [bib35]) d’un système à un degré de liberté, de masse \(m=1\) kg, de période propre \(T=1\) s, et doté d’un amortissement réduit \(\xi \in [0,1[\). Afin d’étudier l’effet de l’overshooting en vitesse, on calcule l’erreur absolue entre les vitesses obtenues avec la solution analytique du problème et la solution numérique au premier instant suivant l’instant initial, et ce pour plusieurs valeurs de pas de temps \(h > 0\). Les résultats obtenus dans les cas \(\xi=0\) (respectivement \(\xi=0.5\)) sont donnés dans la Figure Fig. 201 (respectivement Fig. 202). Dans le cas du schéma HHT classique, l’erreur en vitesse croît significativement lorsque le pas de temps choisi devient grand. En revanche, cet effet est sensiblement limité dans le cas du schéma NOHHT. Ce dernier schéma permet donc bien de limiter les effets de l’overshooting.
Fig. 201 Étude de l’overshooting du schéma HHT - Évolution de l’erreur absolue en vitesse en fonction de la discrétisation temporelle, dans le cas sans amortissement (\(\xi=0\))#
Fig. 202 Étude de l’overshooting du schéma HHT - Évolution de l’erreur absolue en vitesse en fonction de la discrétisation temporelle, dans le cas avec amortissement (\(\xi=0.5\))#
Les schémas explicites#
code_aster permet l’usage de deux schémas explicites différents: celui des différences centrées et celui de Tchamwa-Wielgosz [bib34], qui, lui, est dissipatif. Ces schémas, contrairement aux schémas implicites, sont conditionnellement stables: il existe un pas de temps limite qu’on ne peut dépasser sous peine de divergence.
Ces schémas sont très bien adaptés aux problèmes de dynamique rapide, où l’on cherche à analyser les phénomènes de propagation d’ondes. Leur coût devient prohibitif, au moins sur base physique, quand on doit faire des simulations numériques sur des temps longs, comme c’est le cas quand on cherche à obtenir une réponse vibratoire ou quand le chargement imposé est de longue durée: à partir de plusieurs secondes, comme, par exemple, en séisme.
Il faut cependant signaler que, comparées à une résolution avec schéma implicite, les approches explicites coûtent bien moins cher sur un pas de temps. En effet, alors qu’en implicite on doit véritablement résoudre le système global assemblé de grande dimension à chaque pas (voire même à chaque itération en non linéaire), en explicite on s’affranchit de cette inversion très coûteuse en utilisant un opérateur diagonal (basé sur la seule matrice de masse, donc qui ne demande pas de réactualisation: avec le mot clé MASS_DIAG = “OUI”). Avec un schéma explicite, on ne peut résoudre qu’en accélération (mot clé FORMULATION = “ACCELERATION”).
Le schéma des différences centrées#
Ce schéma (mot-clé SCHEMA_TEMPS, SCHEMA=”DIFF_CENT”), le plus courant dans les codes de dynamique rapide, est un cas particulier de la famille de Newmark. Il se construit en prenant \(\beta =0\) et \(\gamma =1/2\) . Cette méthode d’intégration en temps, qui n’introduit pas de dissipation numérique, est d’ordre 2.
Son pas de temps de stabilité (condition CFL: Courant Friedrichs et Levy [bib25]) est de \(2/\omega\) , avec \(\omega\) qui est la plus haute pulsation propre de la structure considérée.
- Remarques
Sur base physique, on peut montrer que le pas de temps critique est équivalent au temps minimal que met une onde pour traverser le plus petit élément fini du maillage. En pratique, on approxime cette valeur en calculant le minimum, pour tous les éléments finis de: \(\Delta t=D/C\) avec d qui est le diamètre de la plus grande sphère inscrite et \(c\) qui est la célérité des ondes élastiques de compression dans l’élément considéré.
Cette méthode de calcul de la condition CFL, directement applicable sur les éléments massifs, demande à être corrigée pour les éléments de structure.
Pour les applications industrielles classiques, la condition CFL oblige à prendre des pas de temps très petits (généralement de l’ordre de \({10}^{-5}s\) . pour les structures métalliques).
En revanche, sur base modale, de par la forte troncature généralement faite pour les cas industriels, le pas de temps critique, directement donné par \(2/\omega\) avec \(\omega\) qui est la pulsation de coupure, sera bien plus grand que sur base physique.
La résolution explicite se fait toujours en accélération (l’opérateur à inverser se réduit alors à la matrice de masse seule que l’on recommande de «lumper», donc de la rendre diagonale).
Le schéma de Tchamwa-Wielgosz [bib34]#
Ce schéma récent peut être vu comme le pendant du schéma HHT pour l’explicite. L’objectif recherché est le même: introduire une dissipation numérique haute fréquence, qui ne perturbe pas «trop» la solution basse fréquence. Ce schéma de Tchamwa-Wielgosz [bib34] (mot-clé SCHEMA_TEMPS, SCHEMA=”TCHAMWA”) ne se déduit pas de la famille de Newmark.
Tout comme HHT, ce schéma est paramétré par un réel, nommé \(\phi\) . Si l’on prend \(\phi\) =1, le schéma ne dissipe pas. En pratique il est recommandé de choisir une valeur de \(\phi\) comprise entre 1 et 1,1.
Le pas de temps critique vaut celui du schéma des différences centrées divisé par \(\phi\) .
Performance des schémas explicites dans code_aster#
code_aster n’est pas un code dédié à la dynamique rapide, c’est un code généraliste orienté vers des approches implicites. Il n’est donc pas optimisé pour les résolution explicites, en particulier pour ce qui est du calcul des efforts intérieurs et de la structuration des objets générés, comme la structure de donnée solution.
Il n’est donc pas recommandé d’utiliser un schéma explicite pour des problèmes de grande taille, d’autant plus que l’on veut calculer sur un très grand nombre de pas.
Calculs explicites sur base de Ritz avec sous-domaines dans code_aster#
Les commentaires du paragraphe précédent sont toutefois à nuancer dans la mesure où il est possible de procéder à des calculs explicites sur base de Ritz avec sous-domaines (avec degrés de liberté «généralisés») dans Code_Aster .C’est-à-dire qu’on peut au moyen du mot-clé facteur PROJ_MODAL réaliser le calcul avec un schéma d’intégration en temps explicite sur une base modale (ou de Ritz) de projection préalablement calculée.
Cette option intervient dans le cas où l’on veut condenser dynamiquement une partie du modèle au comportement linéaire, en ne calculant strictement par DYNA_NON_LINE que des sous-domaines au comportement non linéaire. Ceci, afin de réduire la taille du modèle de calcul. Dans ce cas, il est nécessaire de calculer une base modale de Ritz sur l’ensemble du modèle: le sous-domaine non linéaire modélisé pour le calcul faisant appel à DYNA_NON_LINE et les autres sous-domaines supposés linéaires condensés dynamiquement. Cette base doit être orthogonalisée par rapport aux matrices de masse et de rigidité linéaire de l’ensemble du modèle. Elle doit simplement être représentative des mouvements activant l’ensemble du modèle.
Le domaine de calcul modélisé comprend ainsi un sous-domaine non linéaire de calcul \(I\) et une série de sous-domaines linéaires \(E\) «externes» au calcul qui seront condensés dynamiquement.
Les matrices \(M\) , \(C\) , \(K\) , les forces \(F\) peuvent se décomposer en une partie associée au domaine \(I\) et une partie associée au(x) sous-domaine(s) \(E\) . L’équilibre de la structure résolu par DYNA_NON_LINE peut s’écrire alors en terme de degrés de liberté d’éléments finis globaux \(U\) :
Soit \(\lbrace \Phi \rbrace\) une base modale de Ritz orthogonalisée par sous-structuration élémentaire sur l’ensemble des sous-domaines \(I\) et \(E\) . On utilise alors la transformée de Ritz: \(U = \Phi \, q`\) et donc en projetant (4600) par \({}^{t}\Phi\) à gauche on obtient:
Les modes \(\Phi\) réduits à \(I\) sont notés \(\phi\) et les degrés de liberté \(U\) réduits à \(I\) sont notés \(u\) donc (4601) devient pour chaque pas de temps:
On a donc besoin des matrices généralisées de masse \(\tilde{M}=({\tilde{M}}_{I}+{\tilde{M}}_{E})\) diagonale ainsi que des matrices des sous-domaines \(E\) : \({\tilde{K}}_{E}\) et \({\tilde{C}}_{E}\) pleines et du vecteur généralisé \({\tilde{F}}_{E}^{\text{ext}}\) qui participent aux termes complémentaires au problème dynamique standard réduit au sous-domaine \(I\) qui serait résolu sans condensation dynamique en absence de sous-domaines \(E\) par DYNA_NON_LINE:
\(\tilde{M}.\ddot{q}={\phi }^{t}.{F}_{I}^{\text{ext}}–{\phi }^{t}.({F}_{I}^{int}(u,\dot{u})+{C}_{I}.\dot{u})\)
À chaque pas de temps, on évalue et ajoute les termes complémentaires issus de la condensation dynamique, cf. (4602) , à ceux du problème dynamique standard sans condensation dynamique et on résout avec le schéma d’intégration temporelle explicite en coordonnées généralisées \(q\) , ayant exprimé avec le schéma explicite les termes en \(\mathrm{u},\dot{\mathrm{u}}\) sur le sous-domaine \(I\) .
La connaissance de la base modale de Ritz orthogonalisée \(\lbrace \phi \rbrace\) permet ensuite alors d’obtenir à la fois les degrés de liberté \(u\) sur le sous-domaine de résolution \(I\) et \(U\) sur l’ensemble du modèle constitué de ses sous-domaines.
On ne renseignera derrière le mot-clé MODE_MECA que les modes \(\lbrace \phi \rbrace\) obtenus par réduction de la base de Ritz au domaine non linéaire de calcul \(I\) traité par DYNA_NON_LINE, derrière le mot-clé MASSE_GENE la matrice de masse généralisée \(\tilde{M}\) , ainsi que derrière les mots-clé RIGI_GENE AMOR_GENE les opérateurs généralisés des sous-domaines \(E\) : \({\tilde{K}}_{E}\) , \({\tilde{C}}_{E}\) et enfin avec le mot-clé facteur EXCIT_GENE \({\tilde{F}}_{E}\) . Un exemple de calcul est fourni par le cas test SDNV107A [v5.03.107].
La nuance principale par rapport à la restriction d’utilisation d’un schéma explicite pour des problèmes de grande taille consiste à la fois par la possibilité de réduire la taille du modèle de calcul et par l’utilisation d’une condition CFL sur le pas de temps liée à la fréquence maximale de la base de projection, ce qui permet de réduire fortement le nombre de pas de temps de calcul.
Shift de masse#
Il est possible d’entrer un coefficient \(\mathrm{coef}\) derrière l’opérande COEF_MASS_SHIFT dans le mot-clé SCHEMA_TEMPS afin de réaliser un «shift» de la matrice de masse \(M\) qui devient: \(M'=M+\mathrm{coef.}K\) .
L’introduction de ce coefficient permet d’inverser en dynamique avec schéma explicite sur base physique la matrice de masse quand celle-ci a des termes nuls pour certains degrés de liberté spécifiques, par exemple la pression pour les éléments de modélisation HM. Il est à noter que ce shift n’est pas nécessaire en dynamique avec schéma explicite sur base modale (présence du mot-clé PROJ_MODAL) car dans ce cas la matrice à inverser est la matrice de masse généralisée qui est toujours inversible.
L’entrée de ce coefficient permet également d’améliorer fortement la convergence en dynamique avec schéma implicite ou explicite quel que soit le type de modélisation en imposant une fréquence de coupure inversement proportionnelle à la valeur de \(\mathrm{coef}\) .Cette pratique, s’apparentant à la méthode de «selective mass scaling» proposée par Lars Olovsson dans de multiples publications, permet de satisfaire les critères de convergence avec des pas de temps plus importants au prix d’une distorsion de l’ensemble des fréquences propres du système, légère cependant pour les basses fréquences.
Une pulsation propre d’origine \(\omega\) devient alors \(\omega '\) elle 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}\) .
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}\) .
Remarques :
Il convient de signaler que même si on active l’option de «shift» de masse, pour le calcul de l’accélération initiale, c’est toujours la matrice de masse non shiftée qui est utilisée. Le fait qu’elle puisse être non-inversible ne pose néanmoins pas de problème, car dans ce cas, on impose l’accélération initiale à 0 en tout point.
La matrice de rigidité utilisée pour ce «shift» est en pratique celle qui est calculée lors de la phase de prédiction de l’algorithme de Newton.
Exemple#
Il s’agit la figure Fig. 203 du mouvement plan de pendulaison d’une barre extensible \(\mathrm{AB}\) de longueur unité, rotulée en \(A\) à un support fixe, libre en \(B\) et abandonnée sans vitesse dans le champ de pesanteur terrestre à partir d’une position définie par l’angle \({\theta}_{0}\) . On néglige tous les phénomènes de dissipation mécanique.
Comme l’amplitude \({\theta}_{0}\) peut être grande – nous la prendrons de \(90°\) – le point \(B\) subit de grands déplacements et le problème est non linéaire.
La période théorique est:
\(T=1,6744s\)
Le calcul du mouvement du pendule par le schéma HHT (\(\alpha\) -méthode) avec \(\alpha =0\) («règle du trapèze») constitue le cas-test SDNL100.
La figure Fig. 204 représente l’évolution pendant une période de la cote du point \(B\) calculée par le schéma «d’accélération moyenne modifiée» avec trois valeurs de \(\alpha\) . La période est divisée en 40pas de temps égaux.
La courbe en trait plein est relative à \(\alpha =0\) . On n’observe pratiquement aucune erreur.
Fig. 203 Pendule de grande amplitude#
Fig. 204 Schéma «d’accélération moyenne modifiée», 40 pas de 0.0419 s#
La courbe en trait d’axe est relative à \(\alpha =-0,1\) . On observe un taux d’amortissement d’environ 2% alors que la figure Fig. 203, pour \(\frac{\Delta t}{T}=0,025\) , prévoit 0,8%. C’est que cette courbe a été établie en linéaire, alors que le mouvement de notre pendule est non linéaire.
La courbe en pointillés est relative à \(\alpha =-0,3\) . Le taux d’amortissement est d’environ 5,8%, alors que celui prévu par la figure Fig. 198 est d’environ 2,2%. L’écart est encore dû à la non-linéarité du problème.
Enfin, les courbes en trait d’axe et surtout en pointillés révèlent un raccourcissement de la période calculée par rapport à la période théorique.
Conclusion#
Le schéma d’intégration temporelle de Newmark et ses variantes dites «accélération moyenne modifiée» et HHT, accompagnés de la méthode de Newton, permettent de traiter de nombreux types de problèmes de dynamique non linéaire, pour des non-linéarités matérielles ou géométriques.
Le traitement des problèmes dynamiques hautement non linéaires et non réguliers, tels que l’analyse des structures en grands déplacements ou de contact-frottement, est sujet à l’instabilité numérique, même avec des méthodes d’intégration temporelle inconditionnellement stables dans le domaine linéaire. On ne parvient alors à intégrer les équations du mouvement par rapport au temps qu’en introduisant de la dissipation artificielle. Tout l’art est de doser cette dissipation pour que, dans la gamme des fréquences qui présentent un intérêt mécanique, elle soit à peu près équivalente à l’amortissement naturel, sans trop décaler le spectre vibratoire de la structure.
Deux schémas explicites sont aussi proposés. On recommande aux utilisateurs de l’opérateur DYNA_NON_LINE de lire la documentation [u2.06.13] dont le rôle est d’aider à bien utiliser ces méthodes de résolutions transitoire.
Bibliographie#
K.J. BATHE: Finite element procedures in engineering analysis. Prentice-Hall (1982).
AUFAURE: Méthodes directes d’analyse dynamique des structures en non linéaire. Note HI-70/93/124 (27 janvier 1994).
H.M. HILBER, T.J.R. HUGHES and R.L. TAYLOR: Improved numerical dissipation for time integration algorithms in structural dynamics. Earthq. Engng Struct. Dyn. 5, 283-292 (1977).
J.-L. LILIEN: Contraintes et conséquences électromécaniques liées au passage d’une intensité de courant dans les structures en câbles. Thèse, Université de Liège (1983).
BALLARD: Dynamique des systèmes mécaniques discrets avec liaisons unilatérales parfaites. C.R.Acad.Sci. Série IIb, Paris, pp.953-958, 1999.
[bib7] [r5.02.01] Algorithme de thermique linéaire transitoire.
[bib8] [r5.03.01] Algorithme non linéaire quasi-statique (opérateur STAT_NON_LINE).
[bib9] [r5.06.01] Méthodes de RITZ en dynamique non linéaire.
[bib10] [r3.03.01] Dualisation des conditions aux limites.
[bib11] [r5.03.08] Intégration des relations de comportement viscoélastiques dans l’opérateur STAT_NON_LINE.
[bib12] [r5.03.50] Contact unilatéral par des conditions cinématiques.
[bib13] [r4.02.02] éléments vibro-acoustiques.
[bib14] [r4.02.04] Couplage Fluide - Structure avec Surface Libre.
[bib15] [r4.05.01] Réponse sismique par analyse transitoire.
[bib16] [r5.03.40] Modélisation statique et dynamique des poutres en grandes rotations.
[bib18] [r4.02.05] éléments de frontière absorbante.
[bib20] [u4.53.01] Opérateur DYNA_NON_LINE.
[bib21] [r5.05.04] Modélisation de l’amortissement en dynamique linéaire.
SANS: Amortissement physique et amortissement numérique dans les schémas d’intégration numériques des études structurales. Note HT-61/01/028/B, octobre 2002.
T.J.R. HUGHES: A note on the stability of Newmark’s algorithm in non linear structural dynamics, pp. 383-386 (19/08/1975).
MILLARD: Contact, frottement, adhésion: bases et avancées récentes en modélisation et simulation numérique. Cours IPSI 10/2004, ch.3-8.
GREFFET: Vers de nouvelles méthodes numériques pour l’intégration temporelle dans le Code_Aster . Note EDF/AMA HT-62/04/016/A, 12/2004.
NOELS, et al.: Automatic time stepping algorithms for implicit numerical simulations of non linear dynamics, Adv. Eng. Software, vol. 33, pp. 589-603, 2002.
M.O OTHMAN, A.H. ELKHOLY: Measurement of friction coefficient under reciprocating sliding conditions. Eur. J. Mech. A/Solids, n°2, 215-225, 1994.
[bib28] M.OHAYON: Interactions fluides-structures. Éd. Masson, Paris, 1992.
M.GERADIN, D.RIXEN: Théorie des vibrations; application à la dynamique des structures, Éd. Masson, Paris, 1996.
P.PEGON: Alternative characterization of time integration schemes. Comput. Meth. Appl. Mech. Engrg. 190, pp. 2707-2727, 2001.
J.A.C.MARTINS, J.T.ODEN: Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws. Nonlinear Analysis, Theor. Meth. Appl., vol. 11, n°3, pp. 407-428, 1987. Voir aussi: J.T.ODEN, J.A.C.MARTINS: Models and computational methods for dynamic friction phenomena. Comp. Meth. Appl. Mech. Eng., vol. 52, pp. 527-634, 1985.
IBRAHIMBEGOVIC, E.L.WILSON: Unified computational model for static and dynamic frictional contact analysis. Int. J. for Num. Meth. In Eng., vol. 34, pp. 233-247, 1992.
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.
KAIPING : A new family of generalized-\(\alpha\) time integration algorithms without overshoot for structural dynamics, Earthquake Engineering & Structural Dynamics, vol. 37, no. 12. Wiley, pp. 1389–1409, May 09, 2008.
ERLICHER, L. BONAVENTURA, O.S. BURSI : The analysis of the generalized-\(\alpha\) method for non-linear dynamic problems. Computational Mechanics 2002; 28:83–104.
G.M. HULBERT, T.J.R. HUGHES : An error analysis of truncated starting conditions in step-by-step time integration: consequences for structural dynamics. Earthquake Engineering and Structural Dynamics 1987; 15:901–910.
O.A. BAUCHAU, G. DAMILANO, N.J. THERON : Numerical integration of non-linear elastic multi-body systems. International Journal for Numerical Methods in Engineering 1995; 38:2727–2751.