r5.06.04 Algorithmes d’intégration temporelle de l’opérateur DYNA_VIBRA / BASE_CALCUL=’GENE’#
Résumé
Ce document décrit les schémas d’intégration temporelle qui sont utilisés pour résoudre dans l’espace des modes des problèmes de dynamique transitoire en mécanique linéaire, avec, pour certains schémas, la prise en compte possible de non linéarités localisées de type chocs, frottements, relations effort-déplacement et effort-vitesse, et l’utilisation possible de la sous-structuration. Les schémas de NEWMARK (implicite, limité ici au linéaire), DEVOGELAERE, et quatre schémas explicites à pas de temps adaptatif, ADAPT_ORDRE1 et ADAPT_ORDRE2RUNGE_KUTTA_32 et RUNGE_KUTTA_54 sont présentés.
Méthodes d’intégration temporelle d’un problème dynamique#
Introduction#
On suppose que les équations régissant l’équilibre dynamique des solides ont été discrétisées par éléments finis. On obtient un système discret d’équations qu’il s’agit d’intégrer dans le temps. Pour cela on choisit une discrétisation \(\lbrace {t}_{i},i\in \mathrm{\mathbb{N}}\rbrace\) de l’intervalle de temps de l’étude \([0,T]\) et on écrit l’équilibre de la structure aux instants \({t}_{i}\) .
De façon générale ces équations prennent la forme suivante:
\(\mathrm{M}.{\ddot{X}}_{t}+\mathrm{C}.{\dot{X}}_{t}+\mathrm{K}.{X}_{t}={\mathrm{R}}_{\text{ext}}(t)+{\mathrm{R}}_{\text{nl}}({X}_{t},{\dot{X}}_{t},{\ddot{X}}_{t})\) éq 2.1 -1
où
\(\mathrm{M}\) est la matrice de masse du système,
\(\mathrm{K}\) est la matrice de rigidité du système,
\(\mathrm{C}\) est la matrice d’amortissement du système,
\({\mathrm{R}}_{\text{ext}}(t)\) est le vecteur des forces extérieures,
\({\mathrm{R}}_{\text{nl}}({X}_{t},{\dot{X}}_{t},{\ddot{X}}_{t})\) est le vecteur des forces non linéaires.
La matrice d’amortissement \(\mathrm{C}\) est en général difficile à évaluer car l’amortissement est souvent fonction de la fréquence. Il est toutefois fréquent de simplifier la prise en compte de l’amortissement en employant le modèle d’amortissement proportionnel, ou modèle de Rayleigh.
Les méthodes de réduction de Rayleigh-Ritz sont présentées dans le document [R5.06.01], Méthodes de Ritz en dynamique linéaire et non linéaire.
Dans le cas où le terme \({\mathrm{R}}_{\text{nl}}({X}_{t},{\dot{X}}_{t},{\ddot{X}}_{t})\) n’est pas nul, la technique des pseudo-forces consiste à projeter sur la base du système linéaire et à maintenir les forces non linéaires au second membre. La technique des pseudo-forces est toujours associée à un schéma d’intégration explicite. De ce fait la prise en compte des non linéarités n’est disponible que pour des schémas explicites. L’ajout des non linéarités ne modifie pas la forme des équations.
Dans la méthode de Ritz, le champ de déplacement \({X}_{t}\) est remplacé par sa projection sur la base modale telle que \({X}_{t}=\Phi {\eta}_{t}\) où \({\eta}_{t}\) est le vecteur des coordonnées généralisées et \(\Phi\) est la base modale, le plus souvent réduite.
Le système dynamique projeté prend la forme suivante, avec \({\eta}_{t}\in {ℝ}^{p}\) :
\({\Phi}^{t}.\mathrm{M}.\Phi .{\ddot{\eta}}_{t}+{\Phi}^{t}.\mathrm{C}.\Phi .{\dot{\eta}}_{t}+{\Phi}^{t}.\mathrm{K}.\Phi .{\eta}_{t}={\Phi}^{t}.{\mathrm{R}}_{\text{ext}}(t)+{\Phi}^{t}.{\mathrm{R}}_{\text{nl}}(\Phi .{\eta}_{t},\Phi .{\dot{\eta}}_{t},\Phi .{\ddot{\eta}}_{t})\) éq 2.1 -2
Lorsque l’hypothèse de Basile ne s’applique pas (amortissement non proportionnel), la matrice d’amortissement projetée n’est pas diagonale. L’intégration du système couplé se fait alors obligatoirement avec l’un des trois schémas suivants: le schéma implicite NEWMARK, le schéma explicite EULER ou les schémas explicites ADAPT_ORDRE1 et ADAPT_ORDRE2.
L’équation obtenue en \({\eta}_{t}\) est de la même forme que l’équation en \({X}_{t}\) . De ce fait, dans la suite du document, on utilisera la notation \({X}_{t}\) aussi bien pour le déplacement en coordonnées généralisées que pour le déplacement dans l’espace physique. Dans le cas de l’opérateur DYNA_VIBRA / BASE_CALCUL=’GENE’, il s’agit de coordonnées généralisées.
Deux classes de méthode peuvent être distinguées dans l’intégration pas à pas des équations d’équilibre, les méthodes d’intégration explicite et les méthodes d’intégration implicite.
Soit le système dynamique linéaire suivant à intégrer dans le temps:
\(M.{\ddot{X}}_{t}+C.{\dot{X}}_{t}+K.{X}_{t}={R}_{\text{ext}}(t)\) éq 2.1 -3
Ce système différentiel du second ordre peut être ramené à un système du premier ordre:
\(N.{\dot{U}}_{t}=H.{U}_{t}+{F}_{t}\) éq 2.1 -4
où \({U}_{t}=(\begin{array}{c}{X}_{t}\\ {\dot{X}}_{t}\end{array})\) , \(N=(\begin{array}{cc}\mathrm{Id}& 0\\ 0& M\end{array})\) , \(H=(\begin{array}{cc}0& \mathrm{Id}\\ -K& -C\end{array})\) , \({F}_{t}=(\begin{array}{c}0\\ {R}_{t}\end{array})\)
Pour intégrer cet système d’équations différentielles, on utilise une discrétisation \(\left\lbrace {t}_{i},i\in \mathbb{N}\right\rbrace\) , ainsi qu’une formule de différences finies pour exprimer la dérivée \({\dot{U}}_{t}\) .
On appellera méthodes d’intégration explicite les méthodes où seule la dérivée \({\dot{U}}_{t}\) fait intervenir des inconnues au temps \({t}_{i+1}\) . De cette façon la détermination des grandeurs cherchées à l’instant \({t}_{i+1}\) ne résulte pas d’une inversion de système faisant intervenir l’opérateur \(H\) . Si de plus, on réalise un «mass-lumping» de façon à rendre la matrice \(M\) diagonale, la détermination de \({U}_{t}\) est particulièrement simple. Ce sont là les principales caractéristiques des méthodes d’intégration explicite.
Les méthodes implicites ou semi-implicites font intervenir la discrétisation de \({U}_{t}\) à un instant postérieur à \({t}_{i}\) , généralement \({t}_{i+1}\) . La détermination des variables passe donc par la résolution d’un système faisant intervenir l’opérateur \(H\) .
Deux notions sont importantes: la consistance, ou l’ordre du schéma d’intégration, et la stabilité.
Les approximations utilisées pour obtenir les opérateurs différentiels définissent la consistance, ou l’ordre du schéma d’intégration. On peut en effet considérer que l’approximation avec laquelle on obtient le déplacement à chaque pas de temps est liée à l’ordre d’approximation des dérivées premières et secondes par rapport au temps.
L’étude de stabilité d’un schéma consiste à analyser la propagation des perturbations numériques au cours du temps. Un schéma stable conserve une solution finie, malgré les perturbations, alors qu’un schéma instable conduit à une explosion numérique ou divergence de la solution.
Pour réaliser une étude de stabilité d’un schéma d’intégration, on met ce dernier sous la forme d’un système récursif linéaire et on détermine les caractéristiques propres de ce système. Si toutes les valeurs propres de l’opérateur de récursivité sont plus petites que 1 en module, le schéma est stable. Sinon il est instable.
Les schémas d’intégration explicite sont généralement conditionnellement stables, ce qui signifie que le pas de temps doit être suffisamment petit pour assurer la stabilité de l’intégration temporelle.
Certains algorithmes implicites ont la propriété d’être inconditionnellement stables, ce qui fait leur intérêt et permet d’utiliser un pas de temps arbitrairement grand.
Les schémas retenus pour l’opérateur DYNA_VIBRA / BASE_CALCUL=’GENE’ sont un schéma implicite, NEWMARK, et quatre schémas explicites, EULER, DEVOGELAERE,ADAPT_ORDRE1 et ADAPT_ORDRE2 (à pas de temps adaptatif). Le choix se fait par le mot-clef METHODE:‘EULER’, ‘DEVOGE’, ‘NEWMARK’, ‘ADAPT_ORDRE1’ ou ‘ADAPT_ORDRE2’.
Méthodes d’intégration implicite#
Introduction#
Les méthodes implicites font intervenir la résolution d’un système matriciel avec l’opérateur précédemment défini. Si les solides sont supposés élastiques linéaires, cela se traduit par la résolution d’un système linéaire à chaque pas de temps.
L’avantage de ces méthodes est leur stabilité inconditionnelle, qui leur permet d’intégrer les équations de la dynamique avec un pas de temps relativement important tout en représentant correctement le comportement des modes les plus bas en fréquence de la structure.
Un version implicite de la méthode de NEWMARK, qui a été programmée dans DYNA_VIBRA / BASE_CALCUL=’GENE’ pour les problèmes linéaires.
La méthode de Newmark [bib1]#
Présentation du schéma#
Newmark [bib1] a introduit deux paramètres \(\gamma\) et \(\beta\) pour le calcul des déplacements et vitesses au pas de temps \(t+\Delta t\) :
\({\dot{X}}_{t+\Delta t}={\dot{X}}_{t}+\Delta t\left[(1-\gamma )\ddot{{X}_{t}}+\gamma {\ddot{X}}_{t+\Delta t}\right]\)
\({X}_{t+\Delta t}={X}_{t}+\Delta t\dot{{X}_{t}}+\Delta {t}^{2}\left[(\frac{1}{2}-\beta )\ddot{{X}_{t}}+\beta {\ddot{X}}_{t+\Delta t}\right]\)
Considérons les équations d’équilibre au temps \(t+\Delta t\) :
\(M\cdot {\ddot{X}}_{t+\Delta t}+C{\dot{X}}_{t+\Delta t}+K{X}_{t+\Delta t}={R}_{t+\Delta t}\)
Reportons les relations précédentes en éliminant \({\dot{X}}_{t+\Delta t}\) et \({\ddot{X}}_{t+\Delta t}\) , il vient:
\(\tilde{K}\cdot {X}_{t+\Delta t}={\tilde{R}}_{t+\Delta t}\) où \(\tilde{K}=K+{a}_{0}M+{a}_{1}C\)
\(\tilde{R}={R}_{t+\Delta t}+C\cdot \left\lbrace {a}_{1}{X}_{t}+{a}_{4}{\dot{X}}_{t}+{a}_{5}{\ddot{X}}_{t}\right\rbrace +M\cdot \left\lbrace {a}_{0}{X}_{t}+{a}_{2}{\dot{X}}_{t}+{a}_{3}{\ddot{X}}_{t}\right\rbrace\)
avec: \(\begin{array}{cccc}{a}_{0}=\frac{1}{(\beta .\Delta {t}^{2})}& {a}_{1}=\frac{\gamma}{(\beta .\Delta t)}& {a}_{2}=\frac{1}{(\beta .\Delta t)}& {a}_{3}=\frac{1}{2\beta }-1\\ {a}_{4}=\frac{\gamma}{\beta}-1& {a}_{5}=\frac{\Delta t}{2}(\frac{\gamma}{\beta}-2)& {a}_{6}=\Delta t.(1-\gamma )& {a}_{7}=\gamma .\Delta t\end{array}\)
Algorithme complet de la méthode de Newmark#
a) initialisation
conditions initiales \({\mathrm{X}}_{0}\) , \(\dot{{\mathrm{X}}_{0}}\) et \(\ddot{{\mathrm{X}}_{0}}\)
choix de \(\Delta t\) \(\gamma ` et :math:\)beta ` et calcul des coefficients \({a}_{\mathrm{1,.}\mathrm{..}},{a}_{8}\) (cf ci-dessus)
assembler les matrices de raideur \(\mathrm{K}\) , de masse \(\mathrm{M}\) et d’amortissement \(\mathrm{C}\)
former la matrice de rigidité effective \(\stackrel{~}{\mathrm{K}}=\mathrm{K}+{a}_{0}\mathrm{M}+{a}_{1}\mathrm{C}\)
factoriser \(\stackrel{~}{\mathrm{K}}\)
b) à chaque pas de temps
calculer le chargement effectif \(\stackrel{~}{\mathrm{R}}\) :
\(\stackrel{~}{\mathrm{R}}={\mathrm{R}}_{t+\Delta t}+\mathrm{C}\cdot \left\lbrace {a}_{1}{\mathrm{X}}_{t}+{a}_{4}{\dot{\mathrm{X}}}_{t}+{a}_{5}{\ddot{\mathrm{X}}}_{t}\right\rbrace +\mathrm{M}\cdot \left\lbrace {a}_{0}{\mathrm{X}}_{t}+{a}_{2}{\dot{\mathrm{X}}}_{t}+{a}_{3}{\ddot{\mathrm{X}}}_{t}\right\rbrace\)
résoudre \(\stackrel{~}{\mathrm{K}}\cdot {\mathrm{X}}_{t+\Delta t}={\stackrel{~}{\mathrm{R}}}_{t+\Delta t}\)
calculer les vitesses et accélérations au temps \(t+\Delta t\)
\({\ddot{\mathrm{X}}}_{t+\Delta t}={a}_{0}\left({\mathrm{X}}_{t+\Delta t}-{\mathrm{X}}_{t}\right)-{a}_{2}{\dot{\mathrm{X}}}_{t}-{a}_{3}{\ddot{\mathrm{X}}}_{t}\)
\({\dot{\mathrm{X}}}_{t+\Delta t}={\dot{\mathrm{X}}}_{t}+{a}_{6}{\ddot{\mathrm{X}}}_{t}+{a}_{7}{\ddot{\mathrm{X}}}_{t+\Delta t}\)
calcul du pas de temps suivant: retour en b)1)
Conditions de stabilité du schéma de Newmark#
La méthode de Newmark utilisée de façon assez répandue dans le domaine de la mécanique, car il permet de choisir l’ordre de l’intégration, d’introduire ou non de l’amortissement numérique, et possède une très bonne précision.
Il est inconditionnellement stable si: \(\gamma >0.5\) et \(\beta >\frac{{(2\gamma +1)}^{2}}{4}\)
On introduit un amortissement numérique positif si \(\gamma >\frac{1}{2}\) et négatif si \(\gamma <\frac{1}{2}\) .
Lorsque \(\gamma =\frac{1}{2}\) et \(\beta =0\) , la formule de Newmark se réduit au schéma des différences centrées. C’est donc alors un schéma explicite.
Une combinaison très souvent employée est \(\gamma =\frac{1}{2}\) et \(\beta =\frac{1}{4}\) , car elle conduit à un schéma d’ordre deux, inconditionnellement stable sans amortissement numérique. C’est le choix qui a été fait dans l’opérateur DYNA_VIBRA / BASE_CALCUL=’GENE’. Le schéma de Newmark de cet opérateur est donc implicite.
Emploi#
Dans DYNA_VIBRA/BASE_CALCUL=’GENE’, ce schéma permet l’intégration de problèmes non-linéaires. Dans le cadre de la sous-structuration dynamique, il permet d’employer une base modale calculée par sous-structuration mais il ne supporte pas le calcul direct sur les bases modales des sous-structures.
Amortissement numérique des schémas implicites#
L’avantage numérique des schémas d’intégration implicite directs réside dans le fait que le pas de temps peut être substantiellement grand par rapport à la plus petite période propre du système sans risquer de causer une instabilité des résultats.
Pour des modes de période propre de l’ordre du pas de temps ou inférieure au pas de temps, les algorithmes d’intégration introduisent un fort amortissement qui contribue à effacer la contribution de des modes élevés (cf [R5.05.02]).
Il n’y pas d’amortissement numérique dans le cas particulier de l’algorithme de NEWMARK avec \(\beta =\frac{1}{4}\) et \(\gamma =\frac{1}{2}\) .
En revanche, les algorithmes implicites on un effet sensible d’allongement des périodes de la réponse de la structure. On constate que pour garantir une bonne précision sur l’amplitude et la phase des déplacements calculés, il faut respecter un critère voisin de:
\(\Delta t<\frac{0.1}{{F}_{\max}}\text{à}\frac{0.01}{{F}_{\max}}\)
où \({F}_{\max}\) est la plus haute fréquence du mouvement que l’on souhaite capturer.
Algorithme TRBDF2#
Présentation#
Cette section traite le schéma d’intégration temporelle TRBDF2 pour Trapezoidal Rule (TR) Backward Differentiation Formula Order 2 (BDF2) et de son utilisation dans le cas de dynamique linéaire et non-linéaire des structures. Nous suivons dans ce qui suit la présentation de [GARDNER2021].
TRBDF2 est un schéma d’intégration temporel en deux étapes. Il a initialement été proposé pour la simulation numérique des microprocesseurs [BANK1985], considérés dans la littérature comme des problèmes raides. Pour illustrer son fonctionnement, on considère la forme canonique:
Ainsi pour intégrer cette équation entre \(t={t}_{n}\) et \(t+\Delta t={t}_{n+1}\) , on enchaîne les deux séquences suivantes:
Tout d’abord, en étape de prédiction, on utilise la règle du trapèze (TR) pour passer de \({t}_{n}\) à \({t}_{n+\gamma }=t+\gamma \Delta t\) , où \(\gamma \in ]0,1[\) . Cela revient à résoudre l’équation suivante où \({f}_{n}=f({y}_{n},{t}_{n})\) :
On utilise ensuite, en étape de correction, le schéma Backward Differentiation Formula d’ordre deux (BDF2, dit aussi schéma d’Euler implicite à trois points) pour passer de \(t={t}_{n}\) et \(t+\Delta t={t}_{n+1}\) en utilisant la solution précédente à \({t}_{n+\gamma }\) :
Ces étapes sont représentées sur la Figure 2.2.3.1-1 .
Figure 2.2.3.1-1Représentation du schéma TRBDF2
Algorithme complet#
a) initialisation
conditions initiales \({X}_{0}\) , \(\dot{{X}_{0}}\) et \(\ddot{{X}_{0}}\)
assembler les matrices de raideur \(K\) , de masse \(M\) et d’amortissement \(C\)
calculer les paramètres \({a}_{1}=\frac{4}{{\left(\gamma \Delta t\right)}^{2}}\) \({a}_{2}=\frac{2}{\gamma \Delta t}\) \({\gamma}_{1}=\frac{1}{\gamma (2-\gamma )}\) \({\gamma}_{2}=\frac{-{(1-\gamma )}^{2}}{\gamma (2-\gamma )}\) \({\gamma}_{3}=\frac{1-\gamma }{2-\gamma }\) . On choisit \(\gamma =2-\sqrt{(2)}\) , ce qui produit la même matrice de rigidité effective pour les deux étapes
former la matrice de rigidité effective \(\stackrel{~}{\mathrm{K}}=\left(\mathrm{K}+{a}_{2}\mathrm{C}+{a}_{1}\mathrm{M}\right)\)
factoriser \(\tilde{K}\)
b) à chaque pas de temps
# 1 ère étape par la méthode du trapèze \(t\to t+\gamma \Delta t\)
calculer une prédiction de l’accélération au temps \(t+\gamma \Delta t\) \({\stackrel{~}{\ddot{\mathrm{X}}}}_{t+\gamma \Delta t}={\mathrm{M}}^{-1}\left({F}^{\mathit{ext}}(t+\gamma \Delta t)-\mathrm{C}{\dot{\mathrm{X}}}_{t}-\mathrm{K}{\mathrm{X}}_{t}\right)\)
calculer le chargement effectif \({\stackrel{~}{\mathrm{R}}}_{t+\gamma \Delta t}\) :
\({\stackrel{~}{\mathrm{R}}}_{t+\gamma \Delta t}=\gamma \Delta t\left(\left({a}_{1}\mathrm{M}+{a}_{2}\mathrm{C}\right){\dot{\mathrm{X}}}_{t}+\frac{1}{2}{a}_{2}\mathrm{M}\left({\ddot{\mathrm{X}}}_{t}+{\stackrel{~}{\ddot{\mathrm{X}}}}_{t+\gamma \Delta t}\right)\right)\)
résoudre \(\stackrel{~}{\mathrm{K}}\cdot \delta {\mathrm{X}}_{t+\gamma \Delta t}={\stackrel{~}{\mathrm{R}}}_{t+\gamma \Delta t}\)
calculer les déplacements, vitesses au temps \(t+\gamma \Delta t\)
\({\mathrm{X}}_{t+\gamma \Delta t}={\mathrm{X}}_{t}+\delta {\mathrm{X}}_{t+\gamma \Delta t}\)
\({\dot{\mathrm{X}}}_{t+\gamma \Delta t}={\dot{\mathrm{X}}}_{t}+2\left(\frac{1}{\gamma \Delta t}\delta {\mathrm{X}}_{t+\gamma \Delta t}-{\dot{\mathrm{X}}}_{t}\right)\)
# 2èmeétape par la méthode BDF2 \(t\to t+\Delta t\)
calculer une prédiction de l’accélération au temps \(t+\Delta t\) \({\stackrel{~}{\ddot{\mathrm{X}}}}_{t+\Delta t}={\mathrm{M}}^{-1}\left({F}^{\mathit{ext}}(t+\Delta t)-\mathrm{C}{\dot{\mathrm{X}}}_{t}-\mathrm{K}{\mathrm{X}}_{t}\right)\)
calculer le chargement effectif \({\stackrel{~}{\mathrm{R}}}_{t+\Delta t}\) :
\(\begin{array}{cc}{\stackrel{~}{\mathrm{R}}}_{t+\Delta t}& =-{a}_{2}\mathrm{M}\left({\mathrm{X}}_{t}-{\gamma}_{1}{\mathrm{X}}_{t+\gamma \Delta t}-{\gamma}_{2}{\mathrm{X}}_{t}-{\gamma}_{3}{\dot{\mathrm{X}}}_{t}\right)\\ & -\left({a}_{1}\mathrm{M}+{a}_{2}\mathrm{C}\right)\left({\dot{\mathrm{X}}}_{t}-{\gamma}_{1}{\dot{\mathrm{X}}}_{t+\gamma \Delta t}-{\gamma}_{2}{\dot{\mathrm{X}}}_{t}-{\gamma}_{3}{\stackrel{~}{\ddot{\mathrm{X}}}}_{t}\right)\end{array}\)
résoudre \(\stackrel{~}{\mathrm{K}}\cdot \delta {\mathrm{X}}_{t+\Delta t}={\stackrel{~}{\mathrm{R}}}_{t+\Delta t}\)
calculer les déplacements, vitesses et accélérations au temps \(t+\Delta t\)
\({\mathrm{X}}_{t+\Delta t}={\mathrm{X}}_{t}+\delta {\mathrm{X}}_{t+\Delta t}\)
\({\dot{\mathrm{X}}}_{t+\Delta t}={\dot{\mathrm{X}}}_{t}+\frac{1}{{\gamma}_{3}\Delta t}\left(\delta {\mathrm{X}}_{t+\Delta t}-{\mathrm{X}}_{t}+{\gamma}_{1}{\mathrm{X}}_{t+\gamma \Delta t}+{\gamma}_{2}{\mathrm{X}}_{t}+{\gamma}_{3}{\dot{\mathrm{X}}}_{t}\right)\)
\({\ddot{\mathrm{X}}}_{t+\Delta t}={\mathrm{M}}^{-1}\left({F}^{\mathit{ext}}(t+\Delta t)-\mathrm{C}{\dot{\mathrm{X}}}_{t+\Delta t}-\mathrm{K}{\mathrm{X}}_{t+\Delta t}\right)\)
calcul du pas de temps suivant: retour en b)1)
Les détails d’obtention des ces expressions sont dans [TARDIEU2021].
Conditions de stabilité#
Le schéma TRBDF2 est d’ordre deux, implicite et donc inconditionnellement stable, car \(\gamma \ge 1/2\) . Il apporte une distorsion de fréquence croissante avec la fréquence. Par ailleurs, il est A-stable et ayant fixé \(\gamma =2-\sqrt{(2)}\) , son rayon spectral est nul, et donc ce schéma est L-stable.
Méthodes d’intégration explicites#
Introduction#
Plusieurs méthodes d’intégration explicite sont présentées: le schéma des différences centrées d’ordre deux, un schéma de Devogelaere-Fu d’ordre quatre et des schémas à pas de temps adaptatif ADAPT et Runge-Kutta. Ces trois méthodes sont disponibles dans l’opérateur DYNA_VIBRA / BASE_CALCUL=’GENE’. Les schémas sont présentés en ne considérant que des forces linéaires. Cependant la prise en compte des forces non linéaires s’en déduit facilement avec la technique des pseudo-forces.
Méthode de Devogelaere-Fu#
Présentation#
Pour présenter l’algorithme de Devogelaere-Fu, abrégé en DEVOGE dans Code_Aster , on met le problème dynamique sous la forme:
\(\mathrm{M}.{\ddot{X}}_{t}+\mathrm{C}.{\dot{X}}_{t}=\mathrm{G}\left(t,{X}_{t}\right)\) où la matrice \(C\) est supposée diagonale.
Les déplacements et les vitesses sont calculés comme suit:
initialisation
\({X}_{-\frac{1}{2}}={X}_{0}-\frac{\Delta t}{2}{\dot{X}}_{0}+\frac{\Delta {t}^{2}}{8}\left({\mathrm{M}}^{-1}.\mathrm{G}\left({t}_{0},{X}_{0}\right)-{\mathrm{M}}^{-1}.\mathrm{C}.{\dot{X}}_{0}\right)\)
\({\dot{X}}_{-\frac{1}{2}}={\left(4I-\Delta t.{\mathrm{M}}^{-1}.\mathrm{C}\right)}^{-1}\left({\left(4I+\Delta t.{\mathrm{M}}^{-1}.\mathrm{C}\right)}^{}.{\dot{X}}_{0}-\Delta t\left(\mathrm{G}\left({t}_{-\frac{1}{2}},{X}_{-\frac{1}{2}}\right)+\mathrm{G}\left({t}_{0},{X}_{0}\right)\right)\right)\)
à chaque pas de temps
\({X}_{n+\frac{1}{2}}={X}_{n}+\frac{\Delta t}{2}{\dot{X}}_{n}+\frac{\Delta {t}^{2}}{24}(4{M}^{-1}.G({t}_{n},{X}_{n})-{M}^{-1}.G({t}_{n-\frac{1}{2}},{X}_{n-\frac{1}{2}})-{M}^{-1}.C(4{\dot{X}}_{n}-{\dot{X}}_{n-\frac{1}{2}}))\)
\({\dot{X}}_{n+\frac{1}{2}}=4{(4I+\Delta t.{M}^{-1}.C)}^{-1}({\dot{X}}_{n}+\frac{\Delta t}{4}(G({t}_{n},{X}_{n})+G({t}_{n+\frac{1}{2}},{X}_{n+\frac{1}{2}})-{M}^{-1}.C.{\dot{X}}_{n}))\)
\({X}_{n+1}={X}_{n}+\frac{\Delta t}{2}{\dot{X}}_{n}+\frac{\Delta {t}_{2}}{6}(4{M}^{-1}.G({t}_{n},{X}_{n})+2{M}^{-1}.G({t}_{n-\frac{1}{2}},{X}_{n-\frac{1}{2}})-{M}^{-1}.C({\dot{X}}_{n}+2{\dot{X}}_{n+\frac{1}{2}}))\)
\({\dot{X}}_{n+1}=6{\left(6I+\Delta t.{\mathrm{M}}^{-1}.\mathrm{C}\right)}^{-1}.\left({\dot{X}}_{n}+\frac{\Delta t}{6}\left(\mathrm{G}\left({t}_{n+1},{X}_{n+1}\right)+4\mathrm{G}\left({t}_{n+\frac{1}{2}},{X}_{n+\frac{1}{2}}\right)+\mathrm{G}\left({t}_{n},{X}_{n}\right)-{\mathrm{M}}^{-1}.\mathrm{C}.\left(4{\dot{X}}_{n+\frac{1}{2}}+{\dot{X}}_{n}\right)\right)\right)\)
Ordre et stabilité du schéma#
Le schéma est d’ordre quatre, les approximations dans l’écriture des dérivées temporelles étant en \(o\left(\Delta {t}^{4}\right)\) . Il possède donc une excellente aptitude à l’intégration de solutions régulières. Son intérêt est en revanche moins manifeste si les fonctions à intégrer présentent des discontinuités (chocs, frottement, etc.)
On peut montrer que pour un système linéaire non amorti le pas de temps garantissant la stabilité vaut: \(\Delta t<\frac{2\sqrt{2}}{{\omega}_{\max}}\) .
Emploi#
Cette méthode est coûteuse en temps de calcul car elle nécessite deux fois l’évaluation du vecteur des forces internes \(G\) , opération particulièrement lourde. Par conséquent elle est peu utilisée en mécanique pour l’intégration directe. Par contre elle est employée par le CEA [bib2] dans le cas des systèmes projetés sur base modale.
Ce schéma permet la prise en compte de non linéarités localisées de type chocs et frottements.
Dans le cadre de la sous-structuration dynamique, il permet d’employer une base modale calculée par sous-structuration mais il ne supporte pas le calcul direct sur les bases modales des sous-structures.
Schémas d’intégration à pas de temps adaptatif#
Introduction : intérêt d’un pas de temps adaptatif#
Réaliser l’intégration temporelle du transitoire d’une structure dans une phase non linéaire pose toujours des problèmes quant au choix du pas temps. L’estimation de l’erreur est rarement accessible lors de l’intégration.
Les schémas d’intégration explicite obligent à respecter un pas de temps maximal pour ne pas diverger. Dans le cas de comportement non linéaire, ce pas ne peut être déterminé a priori et peut changer à chaque itération. Lorsque la rigidité varie très fortement, un pas de temps constant et très fin pour conserver la stabilité du schéma conduit à un nombre d’itérations très grand et à un temps de calcul considérable.
Plusieurs algorithmes d’intégration à pas de temps adaptatif ont donc été développé: ADAPT_ORDRE1, ADAPT_ORDRE2, RUNGE_KUTTA_32 et RUNGE_KUTTA54. Les deux premiers s’appuient sur le schéma de différences centrées, d’ordre deux et sur le schéma d’Euler, d’ordre un. Les deux derniers, ce sont des schémas de la famille de Runge-Kutta avec contrôle de l’erreur. Dans la suite de ce chapitre on ne représentera pas l’algorithme du schéma adaptatif d’ordre un car il est calqué sur le schéma d’Euler. La gestion de l’adaptation du pas de temps est, quant à elle, la même que pour le schéma adaptatif d’ordre deux.
On peut remarquer que ce type de schéma a aussi été programmé dans DYNA_VIBRA/ BASE_CALCUL=’PHYS’ (cf.[R5.05.02]).
Les schémas à pas adaptatifs ADAPT_ORDRE1 et ADAPT_ORDRE2#
Schéma des différences centrées à pas constant#
On présente d’abord le schéma des différences centrées à pas constant sur lequel le schéma ADAPT_ORDRE2 se base. Il s’écrit ainsi:
\({\dot{X}}_{n+\frac{1}{2}}={\dot{X}}_{n-\frac{1}{2}}+\Delta t.{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n})+o(\Delta {t}^{2})\)
\({X}_{n+1}={X}_{n}+\Delta t{\dot{X}}_{n+\frac{1}{2}}+o(\Delta {t}^{2})\)
On constate que la vitesse est exprimée à des indices demi entiers de la discrétisation en temps alors que les déplacements et accélérations sont exprimés aux indices entiers. Écrit de cette façon le schéma est d’ordre deux. Toutefois l’accélération n’est pas immédiatement calculable car la vitesse n’est connue qu’au demi pas précédent. Pour contourner cette difficulté, on peut utiliser plusieurs approximations de la vitesse au pas de temps entier.
méthode 1: supposer que \({\ddot{X}}_{n}\left({X}_{n},{\dot{X}}_{n},{t}_{n}\right)\equiv {\ddot{X}}_{x}\left({X}_{n},{\dot{X}}_{n-\frac{1}{2}},{t}_{n}\right)\) ce qui constitue une approximation valide si l’amortissement est suffisamment faible \({\dot{X}}_{n}={\dot{X}}_{n-\frac{1}{2}}+o(1)\) . Si l’amortissement est important, le schéma perd alors sa précision d’ordre deux.
méthode 2: utiliser une approximation d’ordre un pour la vitesse: ce qui permet de \({\dot{X}}_{n}={\dot{X}}_{n-\frac{1}{2}}+\frac{\Delta t}{2}{\ddot{X}}_{n-1}+o(\Delta t)\) conserver l’ordre deux du schéma.
méthode 3: utiliser un schéma de type prédicteur/correcteur:
prédicteur: \(\lbrace \begin{array}{c}{\dot{X}}_{{n}^{p}}={\dot{X}}_{n-\frac{1}{2}}+\gamma \Delta t{\ddot{X}}_{n-1}\\ {\ddot{X}}_{{n}^{p}}=\ddot{X}({t}_{n},{X}_{n},{\dot{X}}_{{n}^{p}})\end{array}\)
correcteur: \(\lbrace \begin{array}{c}{\dot{X}}_{n}={\dot{X}}_{n-\frac{1}{2}}+\frac{\Delta t}{2}(\beta {\ddot{X}}_{{n}^{p}}+(1-\beta ){\ddot{X}}_{n-1})\\ {\ddot{X}}_{n}=\ddot{X}({t}_{n},{X}_{n},{\dot{X}}_{n})\end{array}\)
où \(\alpha ` et :math:\)beta ` sont deux paramètre à choisir. Park et Underwood [bib6] rapportent qu’effectuer des itérations supplémentaires n’améliore pas de façon sensible la stabilité du schéma.
Adaptation du schéma au pas de temps variable#
Lorsque le pas de temps varie, les expressions du paragraphe précédent ne sont plus valables, l’accélération \({\ddot{X}}_{n}\) n’étant plus nécessairement exprimée au centre de l’intervalle \(\left[{\dot{X}}_{n-\frac{1}{2}},{\dot{X}}_{n+\frac{1}{2}}\right]\) , comme on le voit sur le schéma ci-dessous.
Pour tenir compte de ceci, la vitesse est calculée ainsi:
\({\dot{X}}_{n+\frac{1}{2}}={\dot{X}}_{n-\frac{1}{2}}+\frac{\Delta {t}_{n-1}+\Delta {t}_{n}}{2}{\ddot{X}}_{n}\)
Le schéma ADAPT_ORDRE2 complet s’écrit alors ainsi:
estimation de \({\dot{X}}_{x}\) selon les méthodes 1, 2 ou 3 |
|
\({\dot{X}}_{n+\frac{1}{2}}={\dot{X}}_{n-\frac{1}{2}}+\frac{\Delta {t}_{n-1}+\Delta {t}_{n}}{2}{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n})+o(\Delta {t}^{2})\) |
|
\({X}_{n+1}={X}_{n}+\Delta t{\dot{X}}_{n+\frac{1}{2}}+o(\Delta {t}^{2})\) |
L’ordre du schéma n’est plus rigoureusement égal à deux, le schéma ayant perdu son caractère centré.
Plus \(\Delta {t}_{n}\) et \(\Delta {t}_{n+1}\) sont différents, plus l’ordre du schéma tend vers un. De fortes variations du pas de temps conduisent donc à une perte de précision.
Il est possible de trouver des expressions plus complexes, qui utilisent la vitesse ou l’accélération à l’itération précédente [bib7]. Toutefois la formule présentée ici donne des résultats satisfaisants lorsque le pas de temps diminue mais elle fait baisser la limite de stabilité lorsque le pas de temps augmente. Le remède est de réguler le pas afin qu’il ne s’accroisse que lentement.
Stabilité et précision du schéma#
Pour étudier le schéma, on s’est contenté de l’analyse d’un système à un seul degré de liberté, libre et linéaire, de pulsation propre \(\omega\) et d’amortissement réduit \(\xi\)
\(\ddot{\xi}+2\omega \xi \dot{\xi}+{\omega}^{2}\xi =0\)
La solution approchée, en utilisant le schéma à pas de constant, s’obtient par la relation de récurrence suivante:
\(\mathrm{A}.{Y}_{n}+\mathrm{B}.{Y}_{n-1}=0\)
avec \({Y}_{n}=\left[\begin{array}{c}{\ddot{x}}_{n}\\ {\dot{x}}_{n+\frac{1}{2}}\\ {x}_{n+1}\end{array}\right]\)
\(A\) et \(B\) sont deux matrices qui dépendent de la méthode choisie pour calculer la contribution du terme d’amortissement. On cherche une solution de la forme: \({Y}_{n}=\lambda {Y}_{n-1}\) .
\(\lambda ` est une valeur propre de :math:`{A}^{-1}.B\) et peut s’écriresous la forme suivante:
\({\lambda}_{c}=\exp({\omega}_{c}\Delta t(-{\xi}_{c}\pm i\sqrt{1-{\xi}_{{c}^{2}}}))\) où \({\omega}_{c}\) et \({\xi}_{c}\) sont la pulsation et l’amortissement réduit calculés par l’algorithme.
On peut les comparer à la solution exacte, \({\lambda}_{e}=\exp(\omega \Delta t(-\xi \pm i\sqrt{1-{\xi}^{2}}))\) , ce qui permet d’évaluer l’erreur sur la pulsation et l’erreur sur l’amortissement: \(\frac{\mid {\omega}_{c}-\omega \mid }{\omega}\) et . \(\frac{\mid {\xi}_{c}-\xi \mid }{\xi}\) .
On a étudié [bib 8] et [bib10] les propriétés du schéma selon la méthode employée pour estimer la vitesse aux pas entiers. Il a été empiriquement trouvé que la méthode 3 est à la fois plus précise et plus stable que les méthodes 1 et 2. La méthode, sans surcoût de calcul, permet d’accroître l’ordre du schéma et donne dans la plupart des cas une meilleure précision, sauf en cas d’amortissement faible. Elle est cependant moins stable que la méthode 1. C’est la méthode 2 qui a été finalement retenue dans le schéma ADAPT. Ces études ont permis en outre d’estimer le nombre de points par période nécessaire pour garantir une intégration stable. 20 est une valeur qui donne une bonne marge de sécurité. C’est la valeur choisie par défaut.
Critères d’adaptation du pas de temps#
Les développements précédents permettent de quantifier les erreurs introduites lors du calcul d’un système libre et linéaire. Ces critères ne permettent cependant pas d’adapter le pas de temps. Ils sont en effet délicats à mettre en œuvre dans les cas non linéaires et ne tiennent pas compte des variations de l’excitation.
Un autre approche consiste à étudier l’erreur locale introduite par le schéma à l’aide de développements limités.
La solution exacte d’un système à un degré de liberté vérifie:
\(\begin{array}{c}\lbrace \begin{array}{c}X(t+\frac{\Delta t}{2})=X(t)+\frac{\Delta t}{2}\dot{X}(t)+\frac{\Delta {t}^{2}}{8}\ddot{X}(t)+\frac{\Delta {t}^{3}}{48}\stackrel{\cdots }{X}(t)+o(\Delta {t}^{3})\\ X(t-\frac{\Delta t}{2})=X(t)-\frac{\Delta t}{2}\dot{X}(t)+\frac{\Delta {t}^{2}}{8}\ddot{X}(t)-\frac{\Delta {t}^{3}}{48}\stackrel{\cdots }{X}(t)+o(\Delta {t}^{3})\end{array}\\ \mathrm{\Rightarrow }X(t+\frac{\Delta t}{2})=X(t-\frac{\Delta t}{2})+\Delta t\dot{X}(t)+\frac{\Delta {t}^{3}}{24}\stackrel{\cdots }{X}(t)+o(\Delta {t}^{3})\end{array}\)
La formule d’intégration des différences conduit donc à une erreur de troncature valant:
\({E}_{n}=\frac{\Delta {t}^{3}}{24}\stackrel{\cdots }{X}({t}_{n})\equiv \frac{\Delta {t}^{2}}{12}({\ddot{X}}_{n}-{\ddot{X}}_{n-1})\)
On peut normer cette erreur pour obtenir une erreur relative:
\({e}_{n}=\frac{\Delta {t}^{2}}{12}\frac{\mid {\ddot{X}}_{n}-{\ddot{X}}_{n-1}\mid }{{X}_{n}},{X}_{n}\mathrm{\ne }0\)
Park et Underwood [bib7] ont interprété cette erreur en définissant une «pulsation apparente»:
\({\omega}_{{A}_{n}^{2}}=\frac{{\ddot{X}}_{n}}{{X}_{n}}\)
Appliquée au schéma à différence centrée, cette définition permet d’interpréter l’erreur relative \({e}_{n}\) comme une variation de la pulsation apparente:
\({e}_{n}\equiv \frac{\Delta {t}^{2}}{12}\mid {\omega}_{{A}_{n}^{2}}-{\omega}_{{A}_{n-1}^{2}}\mid\)
De nombreux algorithmes utilisent un critère d’adaptation du pas de temps fondé sur l’erreur de troncature ([bib9], [bib11]). Cependant dans le cas d’un schéma conditionnellement stable, cette méthode ni de s’assurer de la stabilité de l’intégration, ni de garantir une précision pour le calcul des transitoires.
D’autres méthodes utilisent une approximation de la pulsation propre instantanée du système [bib12], à l’aide des matrices de masse et de raideur. Elles ont le défaut de ne pas s’adapter aux forces extérieures et à leurs fluctuations en fréquence.
Il est donc utile de trouver un critère qui tienne compte des deux approches. C’est pourquoi Park et Underwood ont introduit la notion de «fréquence apparente perturbée»:
\({f}_{{\text{AP}}_{n}}=\frac{1}{\mathrm{2p}}\sqrt{\mid \frac{{\ddot{X}}_{n}-{\ddot{X}}_{n-1}}{{X}_{n}-{X}_{n-1}}\mid }\)
Cette grandeur s’interprète comme la fréquence «instantanée» du système.
Dans le cas d’un système à plusieurs degrés de liberté, il faut calculer une fréquence apparente pour chaque degré de liberté et prendre le maximum. Le pas de temps peut être alors choisi pour respecter un minimum de points par période apparente.
Si le dénominateur de l’expression de la fréquence apparente tend vers zéro, celle-ci peut devenir très grande et ne plus avoir de signification. Ceci mène à un raffinement injustifié lorsque la vitesse s’annule. Pour y remédier on ajoute un critère du type:
\(\frac{\mid {X}_{n}-{X}_{n-1}\mid }{\mathrm{Dt}}<{\dot{X}}_{\min}\Rightarrow {f}_{{\text{AP}}_{n}}=\frac{1}{\mathrm{2p}}\sqrt{\mid \frac{{\ddot{X}}_{n}-{\ddot{X}}_{n-1}}{{\dot{X}}_{\min}\mathrm{Dt}}\mid }\)
C’est un intermédiaire entre la fréquence apparente perturbée et l’erreur de troncature. La valeur adéquate de \({f}_{{\text{AP}}_{n}}=\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{n}-{\ddot{X}}_{n-1}}{{X}_{n}-{X}_{n-1}}\mid }\) est difficile à choisir a priori et une valeur inadaptée entraîne une diminution artificielle de la fréquence apparente. Dans le cas d’un système à plusieurs degrés de liberté on contourne cette difficulté en employant les degrés de liberté «voisins»:
\({f}_{{\text{AP}}_{n}}=\underset{1\le i\le \text{nb}\text{dll}}{\max}(\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{n}^{i}-{\ddot{X}}_{n-1}^{i}}{{b}_{n}^{i}}\mid })\)
\({b}_{n}^{i}=\max(\underset{\mid i-j\mid \le \text{lb}}{\max}\mid {X}_{n}^{j}-{X}_{\text{n-1}}^{j}\mid ,{\dot{X}}_{\min}\Delta t )\)
où \(\text{lb}\) désigne par exemple la largeur de bande de la matrice de raideur et où peut être \({\dot{X}}_{\min}\) choisi très petit.
Cette méthode se révèle très efficace dans le cas où les \({X}_{n}^{i}\) désignent des composantes physiques (déplacements). Dans le cas d’une projection sur base modale il n’est pas pertinent d’employer les composantes voisines pour calculer la fréquence apparente. Dans ce cas il vaut mieux revenir au premier critère et utiliser une des deux méthodes suivantes, spécifiées par le mot clef VITE_MIN:
si VITE_MIN=‘NORM’ alors c’est un paramètre variable égal à \(\frac{\Vert {\dot{X}}_{n}\Vert }{100}\) \(\Vert {\dot{X}}_{n}\Vert =\sqrt{\sum_{1\le i\le \text{nbddl}}{({\dot{X}}_{n}^{i})}^{2}}\) . Cette méthode donne de bons résultats lorsque le nombre de degrés de liberté est grand et est inapplicable au cas à un seul degré de liberté. Elle n’est plus indiquée si l’ordre de grandeur de la vitesse est très différent d’un degré de liberté à un autre.
si VITE_MIN=‘MAXI’ alors c’est un paramètre variable et différent pour chaque degré de liberté, \({\dot{X}}_{\min}^{j}=\underset{1\le m\le n}{\max}\frac{\mid {\dot{X}}_{m}^{j}\mid }{1001}\) . Cette méthode a l’avantage de fonctionner quel que soit le nombre de degré de liberté du système mais elle ne peut pas être utilisée si l’ordre de grandeur de la vitesse varie trop au cours du calcul car, dans ce cas, on obtiendrait systématiquement: \(\frac{\mid {X}_{n}^{j}-{X}_{n-1}^{j}\mid }{\Delta t }\le {\dot{X}}_{\min}^{j}\) .
Algorithme du schéma des différences centrées à pas adaptatif#
Les règles évoquées plus haut permettent de fixer un nombre de pas de temps souhaité par période de la réponse en fonction de la précision voulue, \(N\) . Il est ajustable par le mot clef NB_POIN_PERIODE. Le pas de temps \(\Delta {t}_{n}\) doit alors être inférieur à \(\frac{1}{{\mathit{Nf}}_{{\text{AP}}_{n}}}\) . Le mot clef PAS donne le pas de temps initial, \(\Delta {t}_{\text{ini}}\) , et le mot clef PAS_MAXI le pas de temps maximal à ne pas dépasser, \(\Delta {t}_{\max}\) . L’algorithme est décrit schématiquement ci-dessous:
initialisation: \({X}_{0}\) et \({\dot{X}}_{0}\) donnés
\(\Delta {t}_{-1}=0\) , \(\Delta {t}_{0}=\Delta {t}_{\text{ini}}\) et \({\ddot{X}}_{0}=\ddot{X}\left({t}_{0},{X}_{0},{\dot{X}}_{0}\right)\)
initialisation de \({\dot{X}}_{\min}\)
A chaque pas de temps
initialisation de la recherche du pas de temps: \({N}_{\text{iter}}=0\)
calcul de \({\dot{X}}_{n+\frac{1}{2}}\) puis de \({X}_{n+1}\) :
estimation de la vitesse (méthode 2): \({\dot{X}}_{n}={\dot{X}}_{n-\frac{1}{2}}+\frac{\Delta t}{2}{\ddot{X}}_{n-1}\)
vitesse au mi pas: \({\dot{X}}_{n+\frac{1}{2}}={\dot{X}}_{n-\frac{1}{2}}+\frac{\Delta {t}_{n-1}+\Delta {t}_{n}}{2}{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n})\)
déplacement: \({X}_{n+1}={X}_{n}+\Delta t{\dot{X}}_{n+\frac{1}{2}}\)
calcul de l’accélération \({\ddot{X}}_{n+1}\)
calcul de la fréquence apparente:
\(\frac{\mid {X}_{n}-{X}_{n-1}\mid }{\Delta t}\ge {\dot{X}}_{\min}\mathrm{\Rightarrow }{f}_{{\text{AP}}_{n}}=\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{n}-{\ddot{X}}_{n-1}}{{X}_{n}-{X}_{n-1}}\mid }\)
\(\frac{\mid {X}_{n}-{X}_{n-1}\mid }{\Delta t}<{\dot{X}}_{\min}\mathrm{\Rightarrow }{f}_{{\text{AP}}_{n}}=\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{n}-{\ddot{X}}_{n-1}}{{\dot{X}}_{\min}\Delta t}\mid }\)
vérification de l’adéquation entre le pas de temps et la fréquence apparente:
calcul de l’indicateur \(\text{err}=\Delta {t}_{n}{\mathit{Nf}}_{{\text{AP}}_{n}}\)
si \(\text{err}\ge 1\) et \({N}_{\text{itez}}<{{N}_{\text{iter}}}_{\max}\) alors réduction du pas de temps et nouvelle itération de recherche du pas: \(\Delta {t}_{n}\leftarrow 0,75\Delta {t}_{n}\) , \({N}_{\text{iter}}\leftarrow {N}_{\text{iter}}+1\) , retour en 2)
si \(\text{err}\le 0,75\) depuis plus de 5 pas de temps consécutifs, alors accroissement du pas de temps \({N}_{\text{iter}}\leftarrow {N}_{\text{iter}}+1\)
archivage de la solution, calcul éventuel de \({\dot{X}}_{\min}\) et retour en 1) pour l’itération suivante.
Commentaires sur les paramètres de l’algorithme#
Le fait de fixer une borne supérieure \({{N}_{\text{iter}}}_{\max}\) par le mot clef NMAX_ITER_PAS au nombre de réductions du pas de temps permet de s’assurer de la convergence de l’algorithme dans les cas difficiles (par exemple en cas de discontinuité dans les forces extérieures).
Lorsque l’indicateur \(\text{err}\) est supérieur à 1, le pas de temps est multiplié par un facteur fixe (0,75 par défaut mais il peut être modifié par l’utilisateur grâce à l’opérande COEF_DIVI_PAS). Il aurait été possible d’écrire directement: \(\Delta {t}_{n}\to \frac{1}{\text{err}}\Delta {t}_{n}={\mathit{Nf}}_{{\text{AP}}_{n}}\) , ce qui plus intuitif. Mais cette stratégie conduit à un raffinement excessif, la fréquence apparente calculée étant souvent largement supérieure à la fréquence réelle, lorsque l’erreur est grande. Toutefois, en un seul pas de temps, \(\Delta {t}_{n}\) peut être considérablement réduit (facteur \({0,75}^{{N}_{\mathit{iter}}}\) ).
En revanche l’augmentation du pas de temps est toujours beaucoup plus lente (coefficient de multiplication par défaut de 1,1 définissable par COEF_MULT_PAS) et n’a lieu que si l’indicateur est inférieur à 1 pendant cinq pas de temps consécutif. Ces restrictions se justifient par les risques de perte de stabilité ou de précision du schéma lorsque le pas de temps varie trop vite. Un coefficient de 1, 2 ou 1, 3 peut permettre un calcul plus rapide mais expose parfois à des risques d’erreur.
En résumé les valeurs par défaut ont été validées par de nombreux tests et donnent en général satisfaction en termes de précision et de stabilité [bib8].
Performance de l’algorithme#
A précision égale, le nombre d’itérations effectuées par un schéma adaptatif est au moins cinq fois plus faible qu’avec un pas constant dans les phénomènes qui justifient l’utilisation d’un pas variable par l’aspect irrégulier de leur évolution (chocs, excitations discontinues, etc.).
Les études empiriques ont montré que le pas de temps adaptatif permet dans les cas favorables de gagner un facteur deux ou trois en temps de calcul. Ce schéma permet en outre de maîtriser la précision de l’intégration par le moyen du contrôle du nombre de points par période de la réponse.
Dans le cas des systèmes très amortis, les gains peuvent être encore plus importants (calculs cinq à dix fois plus rapides).
En revanche lorsque le pas de temps «idéal» est à peu près constant, l’utilisation du schéma adaptatif se révèle inutile.
Il permet bien sûr la prise en compte des non linéarités localisées de type chocs ou frottements.
En sous-structuration dynamique, il est compatible aussi bien avec l’analyse transitoire sur la base modale restituée sur le système entier ou le calcul transitoire sur les bases distinctes des sous-structures.
Schémas à pas de temps adaptatifs de la famille de Runge-Kutta#
Les méthodes d’intégration explicites à un pas de type Runge-Kutta cherchent à déterminer une solution approchée du problème de Cauchy suivant:
\(\lbrace \begin{array}{c}\dot{y}(t)=f(t,y(t))t\in [{t}_{0,}T]\\ y({t}_{0})={y}_{0}{y}_{0}\in ℝ\end{array}\)
Pour cela, on subdivise l’intervalle \([{t}_{0}:T]:{t}_{0}<{t}_{1}<{t}_{2}<\cdots <{t}_{N}\) et on pose \(\Delta t={t}_{n+1}-{t}_{n}\) . A chaque instant de l’intervalle la solution du problème de Cauchy est donnée par l’expression intégrale:
\(y({t}_{n+1})=y({t}_{n})+{\int}_{{t}_{n}}^{{t}_{n}+\Delta t}f(t,y(t))\)
La solution approchée \({y}_{n}\approx y({t}_{n})\) proposée par les schémas récursifs de Runge-Kutta prennent la forme:
\({y}_{n+1}={y}_{n}+\Delta t\Psi ({t}_{n},{y}_{n},\Delta t)\)
avec:
\(\Psi :[{t}_{0,}T]\times {ℝ}^{d}\times ℝ\to {ℝ}^{d}\)
où on identifie facilement le terme \(\Delta t\Psi ({t}_{n},{y}_{n},\Delta t)\) comme l’approximation du terme intégrale \({\int}_{{t}_{n}}^{{t}_{n}+\Delta t}f(t,y(t))\) .
Ainsi, une méthode de Runge-Kutta à \(s\) étages permettant d’approcher la solution du problème de Cauchy est donnée par la méthode de quadrature suivante:
\(\lbrace \begin{array}{c}{k}_{1}=f({t}_{n},{y}_{n})\\ {k}_{2}=f({t}_{n}+{c}_{2}\cdot \Delta t,{y}_{n}+{a}_{21}\cdot {k}_{1})\\ ⋮\\ {k}_{s}=f({t}_{n}+{c}_{s}\cdot \Delta t,{y}_{n}+{a}_{s1}\cdot {k}_{1}+{a}_{s2}\cdot {k}_{2}+\cdots +{a}_{s\phantom{\rule{2em}{0ex}}s-1}\cdot {k}_{s-1})\\ {y}_{n+1}={y}_{n}+\Delta t\cdot ({b}_{1}\cdot {k}_{1}+{b}_{2}\cdot {k}_{2}+\cdots +{b}_{s}\cdot {k}_{s})\end{array}\)
Un schéma d’intégration de Runge-Kutta est donc entièrement défini par les coefficients \({b}_{i}\) , \({c}_{i}\) et \({a}_{ij}\) avec, par ailleurs, \({c}_{i}={\sum}_{j}{a}_{ij}\) . Dans le cas d’un schéma de Runge-Kutta explicite, les coefficients doivent également vérifier la condition \({a}_{ij}=0\forall j\ge i\) .
En pratique, l’ensemble des coefficients d’un schéma de Runge-Kutta est présenté sous la forme d’un tableau de Butcher comme l’illustre la figure ci-dessous:
\(\begin{array}{c}{c}_{1}\\ {c}_{2}\\ {c}_{3}\\ ⋮\\ {c}_{s}\end{array}\) |
\(\begin{array}{cccc}& & & \\ {a}_{21}& & & \\ {a}_{31}& {a}_{32}& & \\ \cdots & \cdots & \ddots & \\ {a}_{\mathit{s1}}& {a}_{\mathit{s2}}& \cdots & {a}_{ss-1}\end{array}\) |
\({y}_{n+1}\) |
\(\begin{array}{cccc}{b}_{1}& {b}_{2}& \cdots & {b}_{s}\end{array}\) |
De plus, on dira qu’une méthode est d’ordre \(p\) si \(\forall n,0\le n\le N,{e}_{n}=\mathrm{{\rm O}}(\Delta {t}^{p+1}),\Delta t\mathrm{\to }0\)
Schémas de Runge-Kutta emboîtés pour le contrôle du pas de temps adaptatif#
Lors de la résolution numérique du problème de Cauchy décrit plus haut, le choix du pas de temps \(\Delta t\) est déterminant sur l’ordre de grandeur de l’erreur globale. Ainsi, le contrôle de l’erreur commise \({e}_{n}\) lors du calcul de la solution \({y}_{n}\approx y({t}_{n})\) permet de déterminer un choix de pas de temps dit «optimal» de manière à garantir une majoration de l’erreur par une tolérance fournie par l’utilisateur.
Pour ce faire, une procédure classique consiste à employer deux méthodes de Runge-Kutta dites emboîtées. La première méthode d’ordre \(p\) à \(s\) étages sert à calculer la solution approchée \({y}_{n+1}\) , alors que la seconde méthode d’ordre \(\widehat{p}<p\) sert à estimer l’erreur \({e}_{n}=\Vert {y}_{n+1}-{\widehat{y}}_{n+1}\Vert\) pour le contrôle du pas de temps. En général, on a \(\widehat{p}=p-1\) et on note la méthode \(\mathit{RK}\phantom{\rule{0.5em}{0ex}}p(\widehat{p})\) .
L’avantage de cette approche est que l’approximation d’ordre la plus faible \(\widehat{p}\) utilise les mêmes évaluations de \(f\) et donc les mêmes coefficients \({a}_{ij}\) .
Dans le cadre de l’opérateur DYNA_VIBRA / BASE_CALCUL=’GENE’, deux schémas d’intégration explicites de la famille Runge-Kutta sont disponibles:
Le schéma ’RUNGE_KUTTA_32’: il s’agit ici du schéma de Bogacki-Shampine3(2). Ce schéma est d’ordre trois avec 4 étages. Il intègre également une approximation d’ordre deux permettant le contrôle de l’erreur. Bien qu’il y ait 4 étages, il utilise réellement 3 étages car il a la propriété FSAL (First Same As Last). Le tableau de Butcher associé est présenté ci-dessous:
\(\begin{array}{c}0\\ 1/2\\ 3/4\\ 1\end{array}\) |
\(\begin{array}{ccc}& & \\ 1/2& & \\ 0& 3/4& \\ 2/9& 1/3& 4/9\end{array}\) |
\({y}_{n+1}\) |
\(\begin{array}{cccc}2/9& 1/3& 4/9& 0\end{array}\) |
\({\widehat{y}}_{n+1}\) |
\(\begin{array}{cccc}7/24& 1/4& 1/3& 1/8\end{array}\) |
Le schéma ’RUNGE_KUTTA_54’: Il s’agit ici du schéma de Dormand-Prince connu également sous le nom de DOPRI5(4). Ce schéma est d’ordre 5 avec 7 étages. Avec une approximation d’ordre 4 permettant le contrôle de l’erreur. Comme le schéma précédent il a la propriété FSAL et il n’utilise réellement que 6 étages. Le tableau de Butcher associé est présenté ci-dessous:
c0 |
0 |
|||||||
c1 |
1/5 |
1/5 |
||||||
c2 |
3/10 |
3/40 |
9/40 |
|||||
c3 |
4/5 |
44/45 |
-56/15 |
32/9 |
||||
c4 |
8/9 |
19372/6561 |
−25360/2187 |
64448/6561 |
−212/729 |
|||
c5 |
1 |
9017/3168 |
−355/33 |
46732/5247 |
49/176 |
−5103/18656 |
||
c6 |
1 |
35/384 |
0 |
500/1113 |
125/192 |
−2187/6784 |
11/84 |
|
b |
\({\widehat{y}}_{n+1}\) |
5179/57600 |
0 |
7571/16695 |
393/640 |
−92097/339200 |
187/2100 |
1/40 |
d |
\({y}_{n+1}\) |
35/384 |
0 |
500/1113 |
125/192 |
−2187/6784 |
11/84 |
0 |
Ces deux schémas programmés dans l’opérateur DYNA_VIBRA / BASE_CALCUL=’GENE’considèrent un vecteur d’état \({y}_{n}\) comme étant la concaténation du vecteur de déplacements et des vitesses, à savoir:
\({y}_{n}=\left\lbrace \begin{array}{c}{X}_{n}\\ {\dot{X}}_{n}\end{array}\right\rbrace\)
Par ailleurs, la norme utilisée pour le contrôle de l’erreur relative est donnée par [bib14]
\(\mathit{err}=\frac{1}{d}{\sum}_{k=1}^{d}\sqrt{{(\frac{{y}_{n+1}^{k}-{\stackrel{ˆ}{y}}_{n+1}^{k}}{{\mathit{sc}}^{k}})}^{2}}\)
où \(d\) est la dimension du vecteur d’état \(y\) , \({y}_{n+1}^{k}\) et \({\stackrel{ˆ}{y}}_{n+1}^{k}\) les \(k\) -ièmes composantes des vecteurs \({y}_{n+1}\) et \({\stackrel{ˆ}{y}}_{n+1}\) respectivement. Enfin, \({\mathit{sc}}^{k}\) est donné par:
\({\mathit{sc}}^{k}=\mathit{MAX}(\mid {y}_{n}^{k}\mid ,\mid {y}_{n+1}^{k}\mid )+\alpha\)
où \(\alpha\) est un paramètre de régularisation (valeur par défaut \(0,001\) ). Ainsi, l’algorithme contrôle la majoration de l’erreur relative par l’expression \(\mathit{err}\le \mathit{tol}\) où \(\mathit{tol}\) est une tolérance relative donné par l’utilisateur.
Enfin, l’expression du pas de temps optimal fonction de l’erreur commise est:
\(\Delta {t}_{\mathit{opt}}=0,9\cdot \Delta {t}_{n}\cdot {\left(\frac{\mathit{tol}}{\mathit{err}}\right)}^{\frac{1}{p+1}}\)
avec \(p=5\) pour le schéma ’RUNGE_KUTTA_54’ et \(p=3\) pour le schéma ’RUNGE_KUTTA_32’. Pour une meilleure compréhension, les deux algorithmes sont présentés dans le tableau ci-dessous.
Initialisation: \({X}_{0}\) , \({\dot{X}}_{0}\) , \({t}_{0}\) et \(\Delta {t}_{0}\) donnés:
\(n=0\) , \(\Delta t=\Delta {t}_{0}\) et \({\ddot{X}}_{0}=\ddot{X}({t}_{0},{X}_{0},{\dot{X}}_{0})\)
Tant que la condition \({t}_{n}<T\) est satisfaite
Calcul de l’état suivant \({y}_{n+1}\) et de l’erreur relative \(\mathit{err}\) par une méthode de Runge-Kutta
Si la condition \(\mathit{err}\le \mathit{tol}\) est satisfaite (le pas de temps est accepté):
Archivage de \({X}_{0}\) , \({\dot{X}}_{0}\) et \({\ddot{X}}_{0}\)
Passage à l’état suivant:
\({y}_{n}:={y}_{n+1}\)
\({t}_{n}:={t}_{n}+\Delta t\)
\(n:=n+1\)
Calcul de \(\Delta {t}_{\mathit{opt}}\) sous la contrainte \(0,2\cdot \Delta {t}_{n}\le \Delta {t}_{\mathit{opt}}\le 5\cdot \Delta {t}_{n}\) (afin d’éviter des changements brutaux)
Sélection du pas de temps \(\Delta t:=min(\Delta {t}_{\mathit{opt}},{t}_{n}-T)\)
On peut remarquer plusieurs choses, particulière à la méthode ’RUNGE_KUTTA_54’ :
les étages 5 et 6, font le calcul à \(t+\mathit{dt}\) .
que les valeurs de \(d\) correspondent aux coefficients du 6ème étage de la méthode.
l’erreur est calculé à l’aide \((b-d)\chi\) , et le pas est adapté au besoin. On peut donc prendre soit \(b\) soit \(d\) pour calculer la solution à \(t+\mathit{dt}\) .
\(d\) permet de calculer la solution à l’ordre quatre avec un contrôle de l’erreur. On garde donc l’évaluation réalisée au 6ème étage de la méthode pour des raisons de performance CPU. Il n’y a plus besoin d’utiliser ni \(b\) ni \(d\) pour évaluer la solution à \(t+\mathit{dt}\) .
Pour des raisons de performance CPU et d’erreur numérique, on calcule sous forme rationnelle \(b-d\) dans la programmation, pour l’évaluation de l’erreur. Les coefficients \(b\) et \(d\) ne servent donc plus lors la programmation de cette méthode.
Conclusion#
En guise de conclusion, voici résumés les différentes possibilités d’intégration temporelle qu’offre l’opérateur:
Euler (‘EULER’) explicite modifié pour assurer une stabilité conditionnelle,
Schéma de Newmark (‘NEWMARK’) paramétré de façon à être implicite,
Schéma de Devogelaere-Fu (‘DEVOGE’) d’ordre 4,
Quatre schémas adaptatifs explicites:
‘ADAPT_ORDRE1’s’appuyant sur le schéma d’Euleur,
‘ADAPT_ORDRE2’basé sur le schéma de différences centrées,
deux schémas de la famille Runge-Kutta: ’RUNGE_KUTTA_32’ et RUNGE_KUTTA_54’.
Le schéma par défaut est EULER mais ce n’est pas systématiquement le plus adapté. Le schéma de NEWMARK disponible dans le Code_Aster est implicite et garantit une stabilité inconditionnelle mais n’est utilisable que pour des problèmes purement linéaires.
Le schéma DEVOGE est d’ordre 4 et donc est plus précis mais est il est coûteux en temps de calcul.
Les schémas adaptatifs sont plus particulièrement indiqués pour les problèmes avec des non-linéarités localisées, où le pas de temps «idéal» n’est pas constant au cours du transitoire. C’est donc l’expérience de la modélisation qui permet de choisir le schéma le mieux adapté au problème en fonction du rapport (temps de calcul)/précision.
Bibliographie#
K.-J. BATHE, E.-L. WILSON : «Numerical Methods in Finite Element Analysis», Prentice Hall Inc.
ANTUNE, F. AXISA, H.BUNG, F. DOVEIL, E. de LANGRE:«Méthodes d’analyse en dynamique non linéaire des structure» - IPSI.
J.-R. LEVESQUE, P. LABBE et al.: «Module TRANSIS du code POUX in documentation de référence du code POUX» - Rapport interne EDF.
RAUD, P. RICHARD: «Structure et validation du calcul non linéaire des paliers hydrodynamiques du code CADYRO» – Note HP-61/92/041.
G.-D. HAHN: «A Modified Euler Method for Dynamic Analysis» - International Journal for Numerical Methods in Engineering Vol. 32 (1991), pp 943-955.
K.-C. PARK, P.-G. UNDERWOOD: «A Variable-Step Central Difference Method for Structural Dynamic Analysis – Part II: Implementation and performance evaluation» - Computer Methods in Applied Mechanics and Engineering Vol. 23 (1980) pp 259-279.
K.-C. PARK, P.-G. UNDERWOOD: «A Variable-Step Central Difference Method for Structural Dynamic Analysis – Part I: Theoritical Aspects» - Computer Methods in Applied Mechanics and Engineering Vol. 22 (1980) pp 241-258.
JACQUART, S. GARREAU: «Algorithme d’intégration à pas de temps adaptatif dans le Code_Aster » - Note HP-61/95/023.
O.C.. ZIENKIEWICZ, Y.M XIE: «A Posteriori Local Error Estimation and Adaptative Time‑Stepping Procedure for Dynamic Analysis» Earthquake Engineering and Structural Dynamics Vol. 20(1991) pp. 871-887.
A.-C. LEGER, G. JACQUART: «Algorithmes d’intégration de l’opérateur DYNA_VIBRA / BASE_CALCUL=’GENE’ du Code_Aster : documentation de référence» - Note HP-51/96/072.
R.M. THOMAS, I. GLADWELL: «Variable Order Variable Step Algorithms for Second Order systems» - International Journal for Numerical Methods in Engineering Vol. 26 (1998) pp.39‑53.
P.G. BERGAN, E. MOLLESTAD: «An Automatic Time-Stepping Algorithm for Dynamic Problem» - Computer Methods in Applied Mechanics and Engineering Vol. 49 (1985) pp.299‑318.
GAY, S. GRANGER, T. FRIOU: «Présentation d’une méthode de simulation du couplage fluidélastique en régime non linéaire» - Note HT-32/94/015.
HAIRER, S.P. NØRSETT, S. WANNER: «Solving ordinary differential equations I: Nonstiff problems» - Springer-Verlag, second revised edition, 1993
GARDNER : «Numerical Methods for Ordinary and Partial Differential Equations» - book draft,
BANK RE, COUGHRAN JR WM, FICHTNER W, GROSSE EH, ROSE DJ, SMITH RK «Transient Simulations Of Silicon Devices And Circuits» - IEEE Trans CAD 1985 ; CAD-4(4): 436–51
J.BATHE, M. M. I. BAIG:«On a composite implicit time integration procedure for nonlinear dynamics» - Computers and Structures, 83 (2005) 2513–2524
TARDIEU: «Implantation du schéma d’intégration temporelle TRBDF2 dans code_aster» - Note interne EDF R&D 6125-1722-2021-02622-FR