r5.05.02 Algorithmes d’intégration directe de l’opérateur DYNA_VIBRA / BASE_CALCUL=’PHYS’#
Résumé:
Ce document décrit les schémas d’intégration temporelle qui sont utilisés pour résoudre de façon directe des problèmes de dynamique en mécanique linéaire transitoire . Les schémas de NEWMARK et WILSON- \(\theta\) sont détaillés, ainsi que les schémas «différences centrées à pas constant» et «pas de temps adaptatif».
Méthodes d’intégration temporelle directe implicite et explicite d’un problème dynamique#
On suppose que la structure étudiée a un comportement linéaire et que les équations régissant son équilibre dynamique ont été discrétisées par éléments finis. On obtient un système discret d’équations différentielles du second ordre qu’il s’agit d’intégrer dans le temps.
De façon générale, ces équations prennent la forme suivante:
\(M.{\ddot{X}}_{t}+C.{\dot{X}}_{t}\text{}+K.{X}_{t}={R}_{t}\)
\(M\) est la matrice de masse du système,
\(C\) est la matrice d’amortissement visqueux du système,
\(K\) est la matrice de rigidité élastique du système,
\(R\) est le vecteur des forces externes appliquées au système visqueux.
Le système est du second ordre.
Deux classes de méthodes d’intégration peuvent être distinguées pour intégrer pas à pas les équations d’équilibre: ce sont les méthodes d’intégration explicite et implicite.
Voyons ce qui les distingue en examinant l’intégration temporelle du système linéaire suivant:
\(M.{\ddot{X}}_{t}+C.{\dot{X}}_{t}\text{}+K.{X}_{t}={R}_{t}\)
Ce système différentiel du second ordre peut se ramener à un système du premier ordre:
où:
\(u=\left[\begin{array}{c}X\\ \dot{X}\end{array}\right]u=\left[\begin{array}{c}\dot{X}\\ \ddot{X}\end{array}\right]A=\left[\begin{array}{cc}I& 0\\ 0& M\end{array}\right]B=\left[\begin{array}{cc}0& I\\ -K& -C\end{array}\right]F=\left[\begin{array}{c}0\\ R\end{array}\right]\)
Pour intégrer cette équation différentielle, on utilise une discrétisation \({t}_{i}\) de l’intervalle d’étude ainsi qu’une formule de différences finies pour exprimer la dérivée \(\dot{u}\) .
On appelle méthodes d’intégration explicite les méthodes où, dans (4894) écrite au temps \({t}_{i}\) , seule la dérivée \(\dot{u}\) fait intervenir la variable \(u\) 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 \(K\) . Si de plus, on réalise un « mass-lumping » de façon à rendre la matrice \(M\) diagonale, la détermination de \({u}_{i+1}\) 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\) dans (4894) à un instant postérieur à \({t}_{i}\) , généralement \({t}_{i+1}\) , de façon à déterminer les variables du problème à \({t}_{i+1}\) . Leur détermination passe donc par la résolution d’un système faisant intervenir l’opérateur \(K\) .
Deux notions concernant les schémas d’intégration sont importantes: la consistance 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 de 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 (cf. [bib2]).
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é du schéma.
Certains algorithmes implicites ont la particularité d’être inconditionnellement stables, suivant le choix de certains paramètres, ce qui fait leur intérêt et permet d’intégrer le phénomène dynamique avec un pas de temps arbitrairement grand.
Le schéma de Wilson- \(\theta\) et le schéma de NEWMARK peuvent être explicites pour certains choix de leurs paramètres. Ils sont employés pour leurs propriétés de stabilité inconditionnelle, propres aux schémas implicites. Ils seront donc classés ici dans la catégorie des schémas implicites et on verra sous quelles conditions ils donnent les propriétés de stabilité souhaitées.
Deux schémas d’intégration explicites ont été aussi introduits. Il s’agit des schémas DIFF_CENTRE et ADAPT qui sont tous deux basés sur la méthode des différences centrées. Ils sont conditionnellement stables et nécessite pour être performants un matrice de masse diagonalisée. La stabilité conditionnelle conduit à un contrôle du pas de temps qui, exploité dans le cas du schéma ADAPT, permet une adaptation du pas de temps en fonction de la rapidité des phénomènes modélisés.
Le schéma WILSON-\(\theta\) [bib1]#
Présentation du schéma#
On supposera dans ce qui suit que les solides sont élastiques linéaires. Cette méthode part de l’hypothèse que l’accélération est linéaire entre \(t\) et \(t+\theta .\Delta t\) :
En intégrant (4900) en fonction de la variable \(t\) , on obtient:
On écrit les équations d’équilibre au temps \(t+\theta .\Delta t\) avec \(\theta \ge 1\) :
en exprimant \({\dot{X}}_{t+\theta .\Delta t }\) et \({\ddot{X}}_{t+\theta .\Delta t }\) en fonction de \({\ddot{X}}_{t+\theta .\Delta t }\) et de \({X}_{t}\) , \({\dot{X}}_{t}\) et \({\ddot{X}}_{t}\) par le système (4901), (4902), et en remplaçant dans (4903), il vient:
\(\tilde{K}.{X}_{t+\theta .\Delta t }=\tilde{R}\)
où
\(\tilde{K}=K+\frac{3}{(\theta .\Delta t)}.C+\frac{6}{{(\theta .\Delta t)}^{2}}.M\)
\(\tilde{R}={R}_{t}+\theta .({R}_{t+\Delta t}-{R}_{t})+M.({a}_{0}.{X}_{t}+{a}_{2}.{\dot{X}}_{t}+2.{\ddot{X}}_{t})+C.({a}_{1}.{X}_{t}+2.{\dot{X}}_{t}+{a}_{3}.{\ddot{X}}_{t})\)
\(\begin{array}{ccccccc}{a}_{0}=\frac{6}{{(\theta .\Delta t )}^{2}}& & {a}_{1}=\frac{3}{(\theta .\Delta t )}& & {a}_{2}=2.{a}_{1}& & {a}_{3}=\frac{\theta .\Delta t }{2}\end{array}\)
On remonte aux déplacements, vitesses et accélérations au pas \(t+\Delta t\) par les relations:
\(\begin{array}{c}{\ddot{X}}_{t+\Delta t}={a}_{4}.({X}_{t+\theta .\Delta t}-{X}_{t})+{a}_{5}.{\dot{X}}_{t}+{a}_{6}.{\ddot{X}}_{t}\\ {\dot{X}}_{t+\Delta t}={\dot{X}}_{t}+{a}_{7}.({\ddot{X}}_{t+\Delta t}+{\ddot{X}}_{t})\\ {X}_{t+\Delta t}={X}_{t}+\Delta t.{\dot{X}}_{t}+{a}_{8}.({\ddot{X}}_{t+\Delta t}+2.{\ddot{X}}_{t})\end{array}\)
\(\begin{array}{ccccccccc}{a}_{4}=\frac{{a}_{0}}{\theta}& & {a}_{5}=\frac{-{a}_{2}}{\theta}& & {a}_{6}=1-\frac{3}{\theta}& & {a}_{7}=\frac{\Delta t }{2}& & {a}_{8}=\frac{\Delta {t}^{2}}{6}\end{array}\)
Algorithme complet de la méthode WILSON-\(\theta\)#
a) initialisation:
conditions initiales \({X}_{0},{\dot{X}}_{0}\) et \({\ddot{X}}_{0}\)
choix de \(\Delta t\) et \(\theta\) et calcul des coefficients \({a}_{1}\) ,… \({a}_{8}\) (cf. ci-dessus)
assembler les matrices de rigidité \(K\) et de masse \(M\)
former la matrice de rigidité effective \(\tilde{K}=K+{a}_{0}.M+{a}_{1}.C\)
factoriser \(\tilde{K}\)
b) à chaque pas de temps:
calculer le chargement effectif \(\tilde{R}\) |
|
\(\tilde{R}={R}_{t}+\theta .({R}_{t+\Delta t}-{R}_{t})+M.({a}_{0}.{X}_{t}+{a}_{2}.{\dot{X}}_{t}+2.{\ddot{X}}_{t})+C.({a}_{1}.{X}_{t}+2.{\dot{X}}_{t}+{a}_{3}.{\ddot{X}}_{t})\) |
|
résoudre \(\tilde{K}.{X}_{t+\theta .\Delta t}=\tilde{R}\) |
|
calculer les déplacements au temps \(t+\Delta t\) \(\begin{array}{c}{\ddot{X}}_{t+\Delta t}={a}_{4}.({X}_{t+\theta .\Delta t}-{X}_{t})+{a}_{5}.{\dot{X}}_{t}+{a}_{6}.{\ddot{X}}_{t}\\ {\dot{X}}_{t+\Delta t}={\dot{X}}_{t}+{a}_{7}.({\ddot{X}}_{t+\Delta t}+{\ddot{X}}_{t})\\ {X}_{t+\Delta t}={X}_{t}+\Delta t.{\dot{X}}_{t}+{a}_{8}.({\ddot{X}}_{t+\Delta t}+2.{\ddot{X}}_{t})\end{array}\) |
|
calcul du pas de temps suivant: retour au début |
Condition de stabilité du schéma WILSON-\(\theta\)#
La méthode est inconditionnellement stable pour WILSON- \(\theta\) > 1.37, une valeur couramment employée pour \(\theta\) étant 1.4. De plus, la méthode présente de la dissipation numérique pour \(\theta >1\) , d’autant plus importante que \(\theta\) augmente; mais il peut aussi engendrer des amplifications oscillantes parasites en hautes fréquences.
Le mot clef facteur WILSON permet de spécifier l’emploi de cet algorithme et le choix de la valeur de \(\theta\) . Par défaut, la valeur de \(\theta\) est prise à 1.4.
Pour la valeur de \(\theta\) égale à 1, le schéma de Wilson- \(\theta\) est égal au schéma dit d’« accélération linéaire », c’est-à-dire le schéma de Newmark avec \(\gamma =1/2;\beta =1/6\) qui est conditionnellement stable.
Le schéma de Newmark [bib1], [bib2]#
Présentation du schéma#
Newmark a introduit deux paramètres \(\gamma\) et \(\beta\) pour le calcul des positions et vitesses au pas \(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}=\text{}{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.{\ddot{X}}_{t+\Delta t}+C.{\dot{X}}_{t+\Delta t}+\text{}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 en paramétrisation « D-form »:
\(\tilde{K}.{X}_{t+\Delta t}=\tilde{R}\text{où :}\tilde{K}=K+{a}_{0}.M+{a}_{1}.C\)
\(\tilde{R}={R}_{t+\Delta t}+C.\left\lbrace {a}_{1}.{X}_{t}+{a}_{4}.{\dot{X}}_{t}+{a}_{5}.{\ddot{X}}_{t}\right\rbrace +M.({a}_{0}.{X}_{t}+{a}_{2}.{\dot{X}}_{t}+{a}_{3}.{\ddot{X}}_{t})\)
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}_{\mathrm{²7}}=\gamma .\Delta t\end{array}\)
Algorithme complet du schéma de Newmark#
a) initialisation:
conditions initiales \({X}_{0}\) , \({\dot{X}}_{0}\) et \({\ddot{X}}_{0}\)
choix de \(\Delta t\) et calcul des coefficients \({a}_{1}\) ,… \({a}_{8}\) (cf ci-dessus)
assembler les matrices de raideur \(K\) de masse \(M\) et d’amortissement \(\mathrm{C}\)
former la matrice de rigidité effective \(\tilde{K}=K+{a}_{0}.M+{a}_{1}.C\)
factoriser \(\tilde{K}\)
b) à chaque pas de temps:
calculer le chargement effectif \(\tilde{R}\) |
|
2 |
\(\tilde{R}={R}_{t+\Delta t}+M.({a}_{0}.{X}_{t}+{a}_{2}.{\dot{X}}_{t}+{a}_{3}.{\ddot{X}}_{t})+C.\left\lbrace {a}_{1}.{X}_{t}+{a}_{4}.{\dot{X}}_{t}+{a}_{5}.{\ddot{X}}_{t}\right\rbrace\) |
résoudre \(\tilde{K}.{X}_{t+\theta .\Delta t}=\tilde{R}\) |
|
calculer les vitesses et accélérations au temps \(t+\Delta t\) \(\begin{array}{c}{\ddot{X}}_{t+\Delta t}={a}_{0}.({X}_{t+\Delta t}-{X}_{t})-{a}_{2}.{\dot{X}}_{t}-{a}_{3}.{\ddot{X}}_{t}\\ {\dot{X}}_{t+\Delta t}={\dot{X}}_{t}+{a}_{6}.{\ddot{X}}_{t}+{a}_{7}.{\ddot{X}}_{t+\Delta t}\end{array}\) |
|
5 |
calcul du pas de temps suivant: retour au début |
Conditions de stabilité du schéma de Newmark#
La méthode de Newmark est inconditionnellement stable si:
\(\gamma >0.5\) et \(\beta >\frac{{(2\gamma +1)}^{2}}{4}\)
On introduit un amortissement numérique positif si \(\gamma >0.5\) et négatif si \(\gamma <0.5\) . Lorsque \(\gamma =0.5\) et \(\beta =0\) , la formule de Newmark se réduit au schéma des différences centrées. Une combinaison très souvent employée est \(\gamma =0.5\) et \(\beta =\frac{1}{4}\) , car elle conduit à un schéma d’ordre deux, inconditionnellement stable sans amortissement numérique.
Ce schéma d’intégration est utilisé 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 intégré dans l’opérateur DYNA_VIBRA / BASE_CALCUL=’PHYS’. Le mot-clé facteur NEWMARK permet de spécifier l’emploi de ce schéma d’intégration et le choix de la valeur de \(\beta\) et \(\gamma\) . Par défaut, la valeur de \(\beta\) est prise à \(0.25\) et la valeur de \(\gamma\) est prise à \(0.5\) .
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.
Cependant, si le contenu de la réponse réside dans un ensemble de modes propres, dont le plus haut a une fréquence propre \({F}_{\max}\) , on devra encore respecter un critère sur le pas de temps de la forme:
\(\Delta t<\frac{0.1}{{F}_{\max}}\text{à}\frac{0.01}{{F}_{\max}}\)
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 ces modes élevés.
On peut voir sur le graphe ci-après la diminution d’amplitude d’un système à un degré de liberté, sans amortissement, quand on l’intègre par différentes méthodes (Wilson- \(\theta\) et Newmark \(\gamma =\frac{1}{2},\beta =\frac{1}{4}\) )
On vérifie ici que l’algorithme de NEwmark avec ces paramètres ne présente aucun amortissement numérique.
Par contre, les algorithmes implicites ont également un effet assez sensible d’élongation des périodes propres contenues dans la réponse de la structure qui conduit à un déphasage de la solution calculée. Le graphe ci-dessous présente des pourcentages d’élongation de la période propre d’un système à un degré de liberté sans amortissement.
Sur ces deux graphes, on constate que pour garantir une 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 correctement dans l’analyse numérique.
Schéma des différences centrées à pas constant#
Principe#
Le schéma explicite des différences centrées à pas constant s’écrit:
\(\begin{array}{c}{\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}}({t}_{n},{X}_{n},{\dot{X}}_{n})+o(\Delta {t}^{2})\end{array}\)
en paramétrisation en accélérations ( on ne propose pas la paramétrisation en déplacements), avec les notations suivantes:
.
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. On peut aussi l’écrire de la manière suivante avec les valeurs sur trois pas successifs:
\(\begin{array}{c}{\dot{\mathrm{X}}}_{n+\frac{1}{2}}={\dot{\mathrm{X}}}_{n-\frac{1}{2}}+\Delta t{\ddot{\mathrm{X}}}_{n}\left({t}_{n},{\mathrm{X}}_{n},{\dot{\mathrm{X}}}_{n}\right)+o\left(\Delta {t}^{2}\right)\\ {\mathrm{X}}_{n+1}={\mathrm{X}}_{n}+\Delta t{\dot{\mathrm{X}}}_{n+\frac{1}{2}}\left({t}_{n},{\mathrm{X}}_{n},{\dot{\mathrm{X}}}_{n}\right)+o\left(\Delta {t}^{2}\right)\end{array}\phantom{\rule{4em}{0ex}}\Rightarrow \phantom{\rule{4em}{0ex}}\begin{array}{c}{\dot{\mathrm{X}}}_{n+1}={\dot{\mathrm{X}}}_{n}+\frac{1}{2}\Delta t\left({\ddot{\mathrm{X}}}_{n}+{\ddot{\mathrm{X}}}_{n+1}\right)\\ {\mathrm{X}}_{n+1}={\mathrm{X}}_{n}+\Delta t\left({\dot{\mathrm{X}}}_{n}+\frac{1}{2}\Delta t{\ddot{\mathrm{X}}}_{n}\right)\end{array}\)
L’accélération en \({t}_{n}\) n’est pas immédiatement calculable car la vitesse n’est connue qu’au demi-pas de temps précédent (en \({t}_{n-1/2}\) ), ce qui pose problème pour évaluer les termes d’amortissement. Pour contourner cette difficulté, on calcule l’accélération en \({t}_{n}\) par l‘approximation suivante de l’équation du mouvement:
\({\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n})\stackrel{~}{=}{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n-})={M}^{-1}(F({t}_{n})-\mathrm{K.}{X}_{n}-C.{\dot{X}}_{n-})\)
ce qui constitue une approximation valable si l’amortissement est suffisamment faible (\({\dot{X}}_{n}={\dot{X}}_{n-1/2}+o(1)\) ). Le schéma perd sa précision d’ordre deux si l’amortissement de la structure est important.
D’autres méthodes d’approximation de l’accélération peuvent être envisagées. Celle choisie s’est révélée un bon compromis entre la simplicité et la stabilité, comme l’étude décrite dans la référence [bib4] sur la précision et la stabilité de plusieurs méthodes.
Les champs sont archivés aux instants \({t}_{n},{t}_{;+1},\mathrm{...}\) , la vitesse étant approchée à ces instants par la formule suivante:
\({\dot{X}}_{n+1}={\dot{X}}_{n+1/2}+\frac{\Delta t}{2}{\ddot{X}}_{n+1}({t}_{n+1},{X}_{n+1},{\dot{X}}_{n+1/2})\)
Le schéma des différences centrées «raccourcit» les périodes (distorsion de fréquence). Il conserve le moment cinétique. L’énergie totale du système est oscillante au cours des pas de temps.
Stabilité#
Le schéma des différences centrées est conditionnellement stable . Dans le cas d’un système sans amortissement [bib2], le schéma est stable pour un pas de temps vérifiant \(\Delta t <\frac{2}{{\omega}_{\max}}\) où \({\omega}_{\max}\) est la plus grande pulsation propre du système, soit \(\Delta t <\frac{{T}_{\min}}{\pi}\) . Il faut un minimum de pas de temps pour décrire la plus petite période du système \({T}_{\min}\) .
La valeur limite pour le pas de temps diminue lentement lorsque l’amortissement augmente [bib4]. Le critère de stabilité exige alors:
\(\Delta t<\frac{2}{{\omega}_{\max}}\left(\sqrt{1+{\xi}^{2}}-\xi \right)\)
où \(\xi\) désigne l’amortissement réduit.
Par exemple, pour un amortissement de 0,5%, la condition devient \(\Delta t <\frac{{T}_{\min}}{5}\) .
Algorithme#
En résumé, le schéma tel qu’il est introduit dans le Code_Aster se présente de la façon suivante:
0 |
initialisation: \(\Delta t,{X}_{0},{\dot{X}}_{0}\) donnés \({\ddot{X}}_{0}={M}^{-1}(F(t=0))-K.{X}_{0}-C.{\dot{X}}_{0}\) \({\dot{X}}_{-1/2}={\dot{X}}_{0}-\frac{\Delta t}{2}{\ddot{X}}_{0}\) |
1 |
A chaque pas de temps \({X}_{n},{\dot{X}}_{n-1/2},{\ddot{X}}_{n}\) connus \(\begin{array}{c}{\dot{X}}_{n+1/2}={\dot{X}}_{n+1/2}+\Delta t{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n+1/2})\\ {X}_{n+1}={X}_{n}+\Delta t{\dot{X}}_{n+1/2}\\ {\ddot{X}}_{n+1}={M}^{-1}(F({t}_{n})-\mathrm{K.}{X}_{n+1}-C.{\dot{X}}_{n-1/2})\\ {\dot{X}}_{n+1}={\dot{X}}_{n+1/2}+\frac{\Delta t}{2}{\ddot{X}}_{n+1}\end{array}\) |
2 |
archivage éventuel de \({X}_{n+1},{\dot{X}}_{n+1},{\ddot{X}}_{n+1}\) puis retour à l’étape 1) pour le pas suivant. |
Matrice de masse diagonale#
Le calcul de l’accélération nécessite l’inversion de la matrice de masse. Ce schéma explicite devient plus performant si on utilise une matrice de masse concentrée (« mass lumping ») de telle sorte qu’elle soit diagonale. L’inversion ne nécessite alors plus de factorisation et est immédiate.
C’est pourquoi le schéma de différences centrées n’est licite qu’avec des matrices de masse construites de façon diagonale, par l’option ‘MASS_MECA_DIAG’ de l’opérateur CALC_MATR_ELEM.
Vérification du pas de temps#
On a vu que le schéma de différences centrées est stable à condition que le pas de temps, en l’absence d’amortissement, soit inférieur à une valeur limite, égale à \(\Delta t <\frac{2}{{\omega}_{\max}}\) . Dans la pratique on emploie un pas de temps qui vaut de 5% à 20% du pas de temps critique. Il a donc été introduit un test sur le pas de temps qui vérifie que:
\(\Delta t <0,05\frac{2\pi }{\underset{1\le i\le \text{nddl}}{\max}(\sqrt{\frac{{k}_{ii}}{{m}_{ii}}})}\)
où \({k}_{ii}\) et \({m}_{ii}\) sont les termes diagonaux des matrices de raideur et de masse.
Si cette condition n’est pas vérifiée, l’utilisateur est arrêté avec un message lui indiquant le pas de temps maximal qu’il peut utiliser.
Calcul de l’accélération#
Le calcul de l’accélération se fait ainsi:
pour chaque degré de liberté, on teste si le terme diagonal de la matrice de masse correspondant est nul.
s’il n’est pas nul, le terme d’accélération est calculé selon la formule: \({\ddot{X}}_{n+1}={M}^{-1}(F({t}_{n})-K.{X}_{n+1}-\mathrm{C.}{\dot{X}}_{n-})\)
s’il est nul, le terme d’accélération n’est pas calculé. C’est le cas pour des degrés de liberté dits de Lagrange. S’ils correspondent à des degrés de liberté bloqués, il est licite de ne pas tenir compte de la ligne en question et de ne pas calculer son accélération. Dans le cas où le degré de liberté de Lagrange a été introduit pour définir une liaison entre deux degrés liberté, cela n’a plus de sens. Le schéma est donc alors inutilisable et un test arrête l’exécution avec un message explicite.
Schéma à pas de temps adaptatif#
Principe#
Les méthodes de calcul explicite sont particulièrement indiquées dans la simulation de phénomènes rapides, tels que lapropagation des ondes dans les solides. En revanche, elles conviennent moins bien à des phénomènes plus lents puisque la condition de stabilité du schéma impose un pas de temps de l’ordre de la plus petite période propre du système.
Le schéma adaptatif, basé sur le schéma de différences centrées, a été développé pour permettre le calcul de réponses transitoires dans lesquelles les phénomènes «rapides» et «lents». Par exemple lors d’un impact, dans un premier temps des ondes à haute fréquence se propagent et se dissipent dans la structure. Ensuite, la structure ne répond plus que sur ses modes de basses fréquences, les hautes fréquences s’étant amorties. L’idée est donc d’adapter le pas de temps au fur et à mesure en fonction des phénomènes mis en jeu, en fixant un critère de précision sur la solution.
Schéma#
Le schéma explicite des différences centrées à pas variables s’écrit:
\(\begin{array}{c}{\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}}({t}_{n},{X}_{n},{\dot{X}}_{n})+o(\Delta {t}^{2})\end{array}\)
avec les notations suivantes:
.
L’on constate que le pas de temps varie. Il est indicé: \(\Delta {t}_{n}\) .
Cela a pour conséquence que le schéma n’est plus rigoureusement du second ordre, puisqu’il n’est plus «centré». Plus \(\Delta {t}_{n-1}\) et \(\Delta {t}_{n}\) sont différents, plus l’ordre du schéma est proche de un. De fortes variations du pas de temps s’accompagnent donc d’une baisse de précision. La formule de la vitesse employée conduit à de bons résultats lorsque le pas de temps diminue mais fait baisser la limite de stabilité lorsque le pas de temps augmente. C’est pourquoi on le contraint à ne s’accroître que très progressivement.
Enfin, on utilise les mêmes approximations que pour les différences centrées en ce qui concerne le calcul des accélérations et des vitesses aux pas de temps «entiers»:
l’accélération est estimée par \({\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n})\stackrel{~}{=}{\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n-})\) et \({\ddot{X}}_{n}({t}_{n},{X}_{n},{\dot{X}}_{n-})={M}^{-1}.(F({t}_{n})-\mathrm{K.}{X}_{n}-C{\dot{X}}_{n-})\) ;
et la vitesse stockée est évaluée par \({\dot{X}}_{n+1}={\dot{X}}_{n+}+\frac{\Delta t }{2}{\ddot{X}}_{n+1}\) .
Comme pour le schéma de différences centrées, dont il est inspiré, le schéma à pas adaptatif nécessite l’inversion de la matrice de masse. C’est pourquoi on exige la diagonalisation de la matrice de masse ainsi que les mêmes restrictions sur les degrés de liberté de Lagrange que pour le schéma à différences centrées.
Estimation du pas de temps en fonction de la précision exigée#
Pour définir un critère sur le pas de temps en fonction de la précision exigée sur la solution, on introduit la notion de fréquence apparente perturbée [bib4]:
\({f}_{{\text{AP}}_{n}}=\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{x}-{\ddot{X}}_{x-1}}{{X}_{x}-{X}_{x-1}}\mid }\)
Cette grandeur peut s’interpréter comme la «fréquence instantanée» du système. C’est en effet une approximation de la pente locale de la courbe force/déplacement. Elle est liée à l’erreur due à la troncature dans les développements limités du schéma. Elle permet en outre de tenir compte des forces extérieures et de leurs fluctuations en fréquence.
Dans le cas d’un système à plusieurs degrés de liberté, il est nécessaire de calculer une fréquence apparente pour chaque degré de liberté. On emploie alors le maximum sur toutes les fréquences calculées pour déterminer le pas de temps.
Si le dénominateur tend vers zéro, le fréquence apparente peut devenir très grande et perdre sa signification physique. On obtient alors un raffinement injustifié du pas de temps lorsque la vitesse s’annule. Dans le cas d’oscillations sinusoïdales, c’est le cas deux fois par période. On modifie alors le critère en introduisant la condition suivante:
\(\frac{\mid {X}_{x}-{X}_{x-1}\mid }{\Delta t }\le {\dot{X}}_{\min}\Rightarrow {f}_{{\text{AP}}_{n}}=\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{x}-{\ddot{X}}_{x-1}}{{\dot{X}}_{\min}\Delta t }\mid }\)
On obtient un intermédiaire entre la fréquence apparente perturbée et l’erreur de troncature. La valeur de \({\dot{X}}_{\min}\) n’est pas facile à déterminer a priori et une valeur mal choisie peut conduire à une modération artificielle de la fréquence apparente.
On propose deux méthodes.
Influences des nœuds voisins#
Dans le cas d’un système à plusieurs degrés de libertés, on peut se servir de l’information donnée par les \(1\le j\le \text{nv}\) nœuds voisins du nœud \(i\) :
\({f}_{{\text{AP}}_{n}}=\underset{\text{DX},\text{DY},\text{DZ},\text{DRX},\text{DRY},\text{DRZ}}{\max}(\underset{1\le i\le \text{nb}\text{noeud}}{\max}\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{x}^{i}-{\ddot{X}}_{x-1}^{i}}{{b}_{n}^{i}}\mid })\)
où \({b}_{n}^{i}=\Delta {t}_{n}\max({10}^{-15}{\text{ms}}^{-1},{\dot{X}}_{n+}^{i},\frac{1}{100}\underset{1\le j\le \text{nv}}{\max}({\dot{X}}_{n+}^{j}))\)
Cette méthode requiert le recensement des nœuds voisins et l’estimation des «vitesses» selon chaque type de degré de liberté (translation DX, DY, DZ, et éventuellement rotation DRX, DRY et DRZ) pour ces nœuds voisins.
La méthode programmée simplifie quelque peu cette formule et consiste, pour un degré de liberté donné, \(i\) , à faire à partir de cette position une recherche ascendante et une recherche descendante sur les degrés de liberté dans leur ordre de numérotation défini par NUME_DDL. Les deux premiers degrés de liberté, \(k\) et \(l\) , de même nature trouvés respectivement avant et après le degré de liberté \(i\) sont considérés comme les «voisins». Pour limiter le coût de cette technique, la recherche est faite une fois pour toute en début du calcul transitoire et les «voisins» sont enregistrés dans deux tableaux d’entiers.
L’emploi de cette méthode est amorcé par le mot clef VITE_MIN=‘NORM’.
Utilisation de l’information au temps précédent#
On peut aussi s’appuyer sur l’information apportée par les pas de temps précédents pour estimer la vitesse minimale. On l’estime alors par la formule suivante:
\({\dot{X}}_{\min}=\underset{k<n}{\max}(\frac{\mid {\dot{X}}_{k}^{i}\mid }{100},{10}^{-15}{\text{ms}}^{-1})\)
On a alors:
\({f}_{{\text{AP}}_{n}}=\underset{\text{DX},\text{DY},\text{DZ},\text{DRX},\text{DRY},\text{DRZ}}{\max}(\underset{1\le i\le \text{nb}\text{noeud}}{\max}\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{x}^{i}-{\ddot{X}}_{x-1}^{i}}{{b}_{n}^{i}}\mid })\)
avec \({b}_{n}^{i}=\Delta {t}_{n}\max({10}^{\text{-15}}{\text{ms}}^{\text{-1}},{\dot{X}}_{n+}^{i},\frac{1}{100}\underset{k<n}{\max}({\dot{X}}_{k}^{j}))\)
Cette méthode est enclenchée par le mot clef VITE_MIN=‘MAXI’.
Cette méthode ne peut pas être employée si la vitesse varie trop au cours du calcul, car dans ce cas on aurait à chaque pas:
\(\frac{\mid {X}_{n}-{X}_{n-1}\mid }{\Delta t }\le {\dot{X}}_{\min}^{i}\)
Choix du nombre de pas par période apparente, \(N\)#
Les calculs d’erreur et les critères de stabilité établis pour un système à un seul degré de liberté (voir [bib4]) ont permis d’estimer le nombre de pas \(N\) nécessaires par période apparente pour obtenir une bonne précision. Ces essais ont montré qu’unminimum de 20 pas par période est nécessaire. Ce nombre est paramétrable par l‘utilisateur dans le fichier de commande grâce au mot clef ‘NB_POIN_PERIODE’. Sa valeur par défaut établie à 50 mène à une précision sur l’intégration temporelle de l’ordre de 1 à 2%.
Le pas de temps initial sert comme pas de temps maximal dans l’absolu: \(\Delta {t}_{\text{mac}}=\Delta {t}_{\text{initial}}\) . Pondéré par un coefficient paramétrable par PAS_LIMI_RELA, il sert de pas de temps minimal: \(\Delta {t}_{\min}=\text{PLR}\ast \Delta {t}_{\text{initial}}\)
Heuristique d’évolution du pas de temps#
On définit un indicateur, dit d’«erreur», sur le choix du pas de temps:
\(\text{erreur}=\Delta {t}_{n}{\mathit{Nf}}_{{\text{AP}}_{n}}\)
Il faut que cet indicateur soit inférieur à 1 pour espérer garantir une bonne intégration temporelle de la plus petite période propre. Cependant le schéma adaptatif doit concomitamment éviter l’emploi d’un pas de temps trop petit, qui provoquerait alors un surcoût de calcul, voire l’apparition de «bruits» parasites.
En fonction de l’indicateur, l’algorithme va accroître ou décroître le pas de temps. On définit pour cela deux coefficients, \(\text{CDP}\) , le coefficient de raffinement du pas de temps (mot clef ‘COEF_DIVI_PAS’, valeur par défaut: 1,334) et \(\text{CMP}\) , le coefficient d‘amplification du pas temps (mot clef ‘COEF_MULT_PAS’, valeur par défaut: 1,1).
Lors de cette quête du pas de temps optimal, on définit un nombre d’itération maximal de réduction du pas de temps, \({\text{iter}}_{\max}\) , pour éviter au pas de temps d’évoluer de façon trop brutale, ce qui est préjudiciable à l’ordre du schéma, et pour ne pas lancer une optimisation trop coûteuse.
si l’indicateur d’erreur est supérieur à sa valeur limite, que l’on a pas dépassé le nombre limite de raffinement pour un pas de temps et que le pas de temps reste plus grand que sa valeur minimale fixée a priori , on raffine le pas de temps: \(\Delta {t}_{n}>\frac{1}{{\mathit{Nf}}_{{\text{AP}}_{n}}},\text{iter}<{\text{iter}}_{\max}\text{et}\Delta {t}_{n}>\Delta {t}_{\min}\Rightarrow \frac{\Delta {t}_{n}}{\text{CDP}}\to \Delta {t}_{n}\) ,
si l’indicateur montre que depuis cinq pas consécutifs le pas de temps se révèle trop fin, i.e. \(\Delta {t}_{n}<\frac{0,75}{{\mathit{Nf}}_{{\text{AP}}_{n}}}\) , alors \(\min\left(\Delta {t}_{\min},\text{CMP}\Delta {t}_{n}\right)\to \Delta {t}_{n}\)
Algorithme#
l’algorithme a été programmé dans Code_Aster selon l’organigramme suivant:
0 |
Initialisation: \({X}_{0},{\dot{X}}_{0}\) donnés \({\ddot{X}}_{0}={M}^{\text{-1}}(F(t=0))-K{X}_{0}-C{\dot{X}}_{0}\) \({\dot{X}}_{-}={\dot{X}}_{0}-\frac{\Delta t }{2}{\ddot{X}}_{0}\) récupération des paramètres d’intégration: \(\Delta {t}_{\text{initial}}\) CMP coefficient d’amplification du pas de temps CDP coefficient de réduction du pas de temps PLR limite au raffinement telle que \(\Delta t \ge \mathrm{PLR}\Delta {t}_{\text{initial}}\) N nombre de pas de temps par période apparente iter max nombre maximal de réductions du pas de temps |
1 |
à chaque pas de temps: \({X}_{n},{\dot{X}}_{n-},{\ddot{X}}_{n}\) connus \({t}_{n+1}={t}_{n}+\Delta {t}_{n}\) 1.0 iter=0 1.1: intégration temporelle \({\dot{X}}_{n+\frac{1}{2}}={\dot{X}}_{n-\frac{1}{2}}+\frac{\Delta {t}_{n-1}+\Delta {t}_{n}}{2}{\ddot{X}}_{n}\) \({X}_{n+1}={X}_{n}+\Delta {t}_{n}{\dot{X}}_{n+\frac{1}{2}}\) \({\ddot{X}}_{n+1}={M}^{\text{-1}}.(F({t}_{n+1})-K.{X}_{n+1}-\mathrm{C.}{\dot{X}}_{n+\frac{1}{2}})\) \({\dot{X}}_{n+1}={\dot{X}}_{n+\frac{1}{2}}+\frac{\Delta {t}_{n}}{2}{\ddot{X}}_{n+1}\) 1.2 calcul de la fréquence apparente et de l’erreur sur le pas de temps \({f}_{{\text{AP}}_{n}}=\underset{\text{DX},\text{DY},\text{DZ},\text{DRX},\text{DRY},\text{DRZ}}{\max}(\underset{1\le i\le \text{nb}\text{noeud}}{\max}\frac{1}{2\pi }\sqrt{\mid \frac{{\ddot{X}}_{x}^{i}-{\ddot{X}}_{x-1}^{i}}{{b}_{n}^{i}}\mid })\) \(\text{erreur}=\Delta {t}_{n}{\mathrm{Nf}}_{{\text{AP}}_{n}}\) 1.2 test sur la pertinence du pas de temps * si \(\text{erreur}>1\) et \(\text{iter}<{\text{iter}}_{\max}\) alors \(\Delta {t}_{n}/\text{CDP}\to \Delta {t}_{n}\) mais si \(\Delta {t}_{n}<\Delta {t}_{\min}\) arrêt du calcul avec message d’erreur \(\text{iter}+1\to \text{iter}\) et retour en 1.1 * si \(\text{erreur}>1\) et \(\text{iter}>{\text{iter}}_{\max}\) alors émission d’une alarme et passage au point 2. * si \(\text{erreur}<1\) passage au point 2 avec si \(\text{erreur}<0,75\) depuis 5 pas consécutifs: amplification du pas de temps \(\Delta {t}_{n}=\min\left(\Delta {t}_{\max},\text{CMP}\Delta {t}_{n}\right)\) |
2 |
acceptation de la solution: archivage éventuel de \({\mathrm{X}}_{n+1},{\dot{\mathrm{X}}}_{n+1},{\ddot{\mathrm{X}}}_{n+1}\) puis \(n+1\to n\) : retour en 1 pour le pas de temps suivant |
Conclusion#
L’opérateur DYNA_VIBRA / BASE_CALCUL=’PHYS’ permet le choix entre plusieurs méthodes d’intégration temporelle. Dans leur paramétrage par défaut, les schémas de Wilson et de Newmark sont des schémas implicites inconditionnellement stables. Ils nécessitent donc une inversion de système linéaire à chaque pas de temps mais offrent en contrepartie un choix du pas temps qui n’est restreint que par la finesse avec laquelle on souhaite décrire l’évolution temporelle des phénomènes modélisés.
Les schémas DIFF_CENTRE et ADAPT sont explicites ce qui leur évite, dans le cas d’une matrice de masse diagonale, une inversion de matrice coûteuse. Mais la stabilité conditionnelle de ce type de schéma conduit généralement à l’utilisation de petits pas de temps, conditionnés par la plus petite période propre du système. Il n’est donc pas garanti que les schémas explicites soient systématiquement plus rapides. Cela dépend des phénomènes simulés. Si la physique de ces phénomènes nécessite une discrétisation temporelle fine, le pas de temps employé est naturellement dans l’intervalle de stabilité. Dans le cas contraire, les contraintes de stabilité numérique entraîne une inflation dans le nombre de pas de temps très coûteuse.
Le schéma ADAPT met à profit une information sur le contenu fréquentiel de la réponse pour adapter le pas de temps. La discrétisation du temps n’est donc plus imposé par la plus petite période propre du système mais par sa réponse. Cela peut être un avantage lorsque la fréquence de la réponse évolue dans le temps, comme dans le cas des impacts.
Bibliographie#
[bib3] S. COURTIER-ARNOUX et COLL : « Notes du cours de fonction : Dynamique des Structures » - note HP-61/94/189.
[bib5] A.C. LEGER: «Introduction des schémas explicites «différences centrées» et «pas de temps adaptatif» dans l’opérateur DYNA_VIBRA / BASE_CALCUL=’PHYS’ du Code_Aster .