r4.06.02 Calcul modal par sous-structuration dynamique classique#
Résumé:
Ce rapport présente les bases théoriques des méthodes de calcul par synthèse modale. Nous débutons par la description de la transformation de RITZ et des méthodes de recombinaison modale qui en découlent. Puis, nous présentons la synthèse modale qui utilise les techniques classiques de sous‑structuration et de recombinaison modale.
Nous présentons, tout d’abord, les techniques de sous-structuration dynamique classique de CRAIG‑BAMPTON et de MAC‑NEAL.
Les outils de calcul par synthèse modale implémentés dans Code_Aster , sont ensuite passés en revue. On détaillera le calcul modal par sous-structuration puis le calcul de réponse harmonique et enfin le calcul de réponse temporel.
On présente ensuite la condensation dynamique par macro-éléments dynamiques en sous-structuration statique.
La dernière section de ce document présente le principe du calcul de résidu en effort, permettant d’évaluer la qualité d’un modèle réduit pour un calcul donné. On détaille également la méthodologie développée pour l’indicateur de qualité pour le calcul modal pour les modèle généralisés.
Synthèse modale#
La sous-structuration dynamique consiste à déterminer le comportement d’une structure à partir des caractéristiques vibratoires de chacun de ses composants ([bib3] et [bib4]). Les méthodes implémentées dans Code_Aster , utilisent simultanément les techniques classiques de recombinaison modale et de sous-structuration dynamique.
Ces méthodes, bien que différentes de celle des éléments finis, adoptent une démarche sensiblement comparable mais introduisent toutefois une approximation supplémentaire. Elles font apparaître trois étapes essentielles :
Étape 1: étude numérique de chaque composant par la détermination de leurs caractéristiques vibratoires. Le travail consiste à identifier des modes propres et des déformées statiques par des techniques classiques de la mécanique vibratoire. Si on assimile chaque sous-structure à un super‑élément, cette étape est similaire au calcul élémentaire.
Étape 2: raccordement des sous-structures. On utilise les caractéristiques vibratoires déterminées précédemment pour chaque composant, et on tient compte de leur liaisonnement. Ce travail constitue l’étape de sous-structuration proprement dite. Il s’apparente à un assemblage.
Étape 3: la résolution et une phase de remontée permettent d’obtenir la solution recherchée dans le repère physique de la structure globale.
Transformation de RITZ#
La transformation de RITZ fait l’objet de la documentation de référence [R5.06.01]. Nous rappelons ici son principe. Pour le problème de la détermination numérique des modes propres réels du système non amorti associé à la structure, que nous désignerons par modes propres , on se ramène à la résolution du problème de minimisation suivant:
Soit \(\delta\) déplacement virtuel, on cherche: \({\mathit{Min}}_{\delta}\frac{1}{2}{\delta}^{T}.(\mathrm{K}-{\omega}^{2}\mathrm{M}).\delta\)
dont la solution :math:`` vérifie:
La méthode de RITZ consiste à chercher la solution de l’équation de minimisation précédente sur un sous‑espace de l’espace des solutions. Considérons la matrice
contenant les vecteurs de la base du sous-espace en question, organisés en colonnes. Restreinte à cet espace de dimension réduite, l’équation de minimisation prend la forme suivante:
\({\mathrm{Min}}_{p=\Phi \delta }\frac{1}{2}{p}^{T}.(K-{\omega}^{2}M).p\)
Soit la solution recherchée:
qui vérifie:
où \(\eta\) est le vecteur des déplacements généralisés, \(\stackrel{ˉ}{K}={\Phi}^{T}K\Phi\) et \(\stackrel{ˉ}{M}={\Phi}^{T}M\Phi\) sont appelées respectivement matrices de raideur et de masse généralisées.
Après avoir résolu le système [éq 2.1-3], l’obtention des modes propres dans la base physique se fait à l’aide de la relation [éq 2.1-2]. La transformation de RITZ permet donc de remplacer le problème aux valeurs propres initial [éq 2.1-1] par un problème de même nature [éq 2.1-3], mais de dimension réduite. Les nouvelles matrices de rigidité et de masse généralisées restent symétriques.
Cependant, cette transformation doit être utilisée avec prudence. En effet, la nouvelle base étant incomplète, une approximation est faite au niveau de la projection: on parle de troncature modale. La précision du résultat final dépend alors du choix des vecteurs de base et l’erreur relative due à cette diminution du nombre d’inconnues doit être estimée.
Recombinaison modale#
Une utilisation classique de la transformation de RITZ, est l’analyse dynamique par recombinaison modale. Elle est couramment utilisée pour le calcul de la réponse d’une structure à une excitation basse fréquence. Nous nous limiterons ici au calcul de la réponse à une excitation d’une structure conservative. Dans ce cas, la méthode des éléments finis nous permet de nous ramener à l’équation différentielle matricielle suivante:
Si l’on applique la transformation de RITZ, avec comme base incomplète de projection, les premiers modes propres de la structure, la relation [éq. ] devient:
où \({\stackrel{ˉ}{f}}_{\mathrm{ext}}={\Phi}^{T}{f}_{\mathrm{ext}}\) est le vecteur des forces généralisées.
Les modes propres sont orthogonaux relativement aux matrices de masse et de rigidité. L’équation différentielle [éq 2.2-2] fait donc apparaître des matrices diagonales: le système est alors constitué d’équations découplées. Chacune d’entre elles est l’équation d’un oscillateur à un degré de liberté de type masse-ressort qui fait apparaître les masse, rigidité et force généralisées relatives au mode \(j\) (respectivement: \({m}_{j}\) , \({k}_{j}\) , \({f}_{j}\) ), voir Figure-a.
Si l’on considère la transformation de RITZ [éq ], au niveau d’un degré de liberté, on a:
\({q}_{i}=\underset{j}{\Sigma}{\Phi}_{ij}{\eta}_{j}\)
où:
\({q}_{i}\) est la \(i\) ème coordonnée du vecteur \(q\) ,
\({\eta}_{j}\) est la \(j\) ème coordonnée du vecteur \(\eta\) ,
\({\Phi}_{ij}\) est la composante de la \(i\) ème ligne et de la \(j\) ème colonne de la matrice \(\Phi\) .
Il apparaît donc que la réponse de la structure s’exprime comme la recombinaison pondérée de réponses d’oscillateurs à un degré de liberté découplés. La transformation de RITZ permet, dans ce cas, de définir un schéma équivalent de la structure, qui fait apparaître les oscillateurs à un degré de liberté associés aux modes propres identifiés. Leur raideur et leur masse sont les rigidités généralisées (\({k}_{j}\) ) et les masses généralisées (\({m}_{j}\) ) des modes correspondants.
Figure -a: Principe de l’analyse dynamique par recombinaison modale.
Réservée essentiellement aux études que l’on qualifie abusivement de basses fréquences [1] , la recombinaison modale consiste à utiliser les propriétés d’orthogonalité des modes propres d’une structure pour simplifier l’étude de sa réponse vibratoire. Outre l’intérêt de diminuer l’ordre du problème numérique à résoudre, la transformation de RITZ en base modale, dans ce cas, permet également de découpler les équations différentielles et de dégager une interprétation physique du résultat obtenu. Selon la fréquence d’excitation, on utilisera une base modale plus ou moins tronquée. Il faut cependant estimer l’erreur de troncature pour s’assurer de la validité du résultat.
Synthèse modale#
De façon générale, les méthodes de synthèse modale consistent à utiliser simultanément la sous‑structuration dynamique (découpage en sous-structures) et la recombinaison modale au niveau de chaque sous-structure. Souvent confondue, par abus de langage, avec la sous-structuration dynamique, la synthèse modale n’est qu’un cas particulier de celle-ci.
La sous-structuration dynamique consiste à considérer le déplacement d’une sous-structure dans le mouvement d’ensemble, comme sa réponse aux forces de liaison qui la relient aux autres composants et aux forces extérieures qui lui sont appliquées.
La synthèse modale signifie que l’on calcule ce mouvement, au niveau de chaque sous-structure, par recombinaison modale. On utilise donc une base de projection qui caractérise chaque sous‑structure. En effet, si la structure globale est trop importante pour être soumise à un calcul modal, les dimensions des sous-structures permettent d’effectuer ce travail. La synthèse modale impose d’étudier d’abord séparément chaque composant, afin de déterminer leur base de projection.
Dans la suite de ce chapitre, nous présentons les types de modes et déformées statiques utilisés dans les méthodes de synthèse modale à l’aide de l’exemple simple suivant:
Figure -a : Principe de la sous-structuration dynamique.
Le vecteur des degrés de liberté de la sous-structure est caractérisé par un exposant qui définit le numéro de la sous-structure, et un indice qui permet de distinguer les degrés de liberté internes (indice \(i\) ), des degrés de liberté de frontière (indice \(j\) ).
\({q}^{k}=\left\lbrace \begin{array}{c}{q}_{i}^{k}\\ {q}_{j}^{k}\end{array}\right\rbrace\)
On est amené, pour étudier la sous-structure \(k\) , à définir une impédance au niveau des degrés de liberté de liaison. Dans le cadre des développements réalisés dans Code_Aster , elle est soit nulle, soit infinie.
Les vecteurs de base utilisés dans les méthodes implémentées dans Code_Aster sont:
les modes normaux,
les modes contraints,
les modes d’attache,
les modes d’interfaces,
que l’on définit ci-dessous.
Les modes normaux#
Les modes propres ou modes normaux sont avantageusement utilisés comme base de projection des sous-structures pour plusieurs raisons:
ils peuvent être calculés ou/et obtenus expérimentalement,
ils offrent des propriétés intéressantes d’orthogonalité par rapport aux matrices de masse et de rigidité de la sous-structure,
ils sont associés à des fréquences propres naturelles de la structure.
Ils peuvent être de deux types selon la condition donnée aux interfaces de liaison, cf. Figures-b et -c:
modes propres à interfaces bloquées,
modes propres à interfaces libres.
Figure -b: Modes propres à interfaces bloquées.
Figure -c: Modes propres à interfaces libres.
Notons que dans le cas d’une sous-structure libre, les modes de corps rigide (ou modes d’ensemble) existants font partie de la base de transformation.
Les déformées statiques#
On définit un mode d’interface à chaque degré de liberté de liaison de chaque sous-structure. Selon les cas, il peut s’agir de modes contraints ou de modes d’attache.
Les modes contraints sont des déformées statiques que l’on joint aux modes normaux à interfaces bloquées pour corriger les effets dus à leurs conditions aux limites. Un mode contraint est défini par la déformée statique obtenue en imposant un déplacement unité sur un degré de liberté de liaison, les autres degrés de liberté de liaison étant bloqués.
Figure -d: Modes contraints.
Les modes d’attache sont des déformées statiques que l’on joint aux modes normaux à interfaces libres pour diminuer l’effet de la troncature modale. Un mode d’attache est défini par la déformée statique obtenue en imposant une force unitaire à un degré de liberté de liaison, les autres degrés de liberté de liaison étant libres.
Figure -e: Modes d’attache.
Dans le cas d’une sous-structure possédant des modes de corps rigide (ici, la sous-structure2), sa matrice de rigidité n’est plus inversible et il n’est pas possible de calculer ses modes d’attache. Il est alors nécessaire de bloquer certains degrés de liberté pour rendre la structure isostatique, voire hyperstatique.
Les déformées harmoniques#
Une troisième base de projection a été introduite dans le cadre des développements concernant le calcul de réponse harmonique par sous-structuration dynamique classique. Il s’agit de la base de Craig-Bampton harmonique.
La base de Craig-Bampton harmonique est constituée de modes propres à interfaces bloquées et de modes contraints harmoniques [bib8]. Ces derniers sont joints aux modes normaux à interfaces bloquées pour corriger les effets dus à leurs conditions aux limites. Un mode contraint harmonique est défini par la réponse de la sous-structure non amortie à un déplacement harmonique, d’amplitude unité et de fréquence donnée, imposé sur un degré de liberté de liaison, les autres degrés de liberté de liaison étant bloqués.
Figure -f: Modes contraints harmoniques.
Cette base modale est plus particulièrement appropriée aux problèmes d’interactions fluide‑structure, pour lesquels les chargements statiques ne sont pas applicables dans Code_Aster (non prise en compte de l’effet de masse ajoutée). Elle peut être utilisée pour tout type de calcul (modal, harmonique et transitoire).
Les déformées d’interface réduits ou « modes de couplage »#
On définit un mode d’interface pour l’ensemble des degrés de liberté de liaison de chaque sous-structure.
Les modes d’interface sont des déformées dynamique que l’on joint aux modes normaux à interfaces bloquées ou libres pour corriger les effets dus à leurs conditions aux limites.
Concrètement, le calcul des modes contraints consiste à effectuer une condensation du comportement de la structure sur les degrés de liberté d’interface et construire ainsi un opérateur d’interface (complément de Schur) qui exprime le comportement de la structure réduite à l’interface. Par une analyse spectrale, on peut alors représenter les déformations de l’interface par combinaison linéaire des modes propres de cet opérateur d’interface.
Figure -g: Mode d’interface.
Conditions de liaison entre sous-structures#
Introduction#
Considérons le problème de deux sous-structures \({S}^{1}\) et \({S}^{2}\) liaisonnées de manière rigide non soumises à des efforts extérieurs. Pour que le mouvement de la structure complète soit continu, il faut imposer l’égalité des déplacements des deux composants à l’interface et la loi d’action-réaction:
où:
\({u}^{1}(M)\) représente le champ de déplacements de la sous-structure \({S}^{1}\) ,
\({u}^{2}(M)\) représente le champ de déplacements de la sous-structure \({S}^{2}\) ,
\({f}_{L}^{1}(M)\) représente le champ des forces de liaison appliquées à la sous-structure \({S}^{1}\) ,
\({f}_{L}^{2}(M)\) représente le champ des forces de liaison appliquées à la sous-structure \({S}^{2}\) .
Suivant la nature des maillages à l’interface, le problème précédent peut se discrétiser différemment.
Cas du maillage d’interface compatible#
Dans cette partie, nous nous limitons au cas de maillages compatibles. Cela signifie qu’ils vérifient les propriétés suivantes:
la restriction du maillage de chacune des sous-structures \({S}^{1}\) et \({S}^{2}\) à l’interface s’appuient strictement sur les mêmes nœuds en leur intersection,
les éléments finis associés à ces mailles de liaison sont de même nature (linéaire, quadratique) de part et d’autre de la frontière.
Dès lors, la condition [éq ] est strictement équivalente à la formulation ci-dessous:
où:
\({q}_{{S}^{1}\cap {S}^{2}}^{k}\) |
est le vecteur des degrés de liberté aux nœuds d’interface \({S}^{1}\cap {S}^{2}\) de la sous‑structure \(k\) , |
\({f}_{{S}^{1}\cap {S}^{2}}^{k}\) |
est le vecteur des forces de liaison aux nœuds d’interface \({S}^{1}\cap {S}^{2}\) de la sous‑structure \(k\) . |
En effet, les maillages des deux sous-structures
et
coïncidant sur l’interface, les fonctions de forme associées aux éléments finis sont les mêmes à l’interface. Il suffit donc d’imposer l’égalité aux nœuds des interfaces de liaison de chaque sous-structure pour imposer l’égalité sur tout le domaine de liaison.
Introduisons les matrices d’extraction des degrés de liberté d’interface \({B}_{{S}^{1}\cap {S}^{2}}^{k}\) :
En utilisant l’équation de projection [éq. ], la condition de continuité des déplacements [éq ] et la formulation ci-dessus appliquées aux deux sous-structures, on obtient:
\({B}_{{S}^{1}\cap {S}^{2}}^{1}{\Phi}^{1}{\eta}^{1}={B}_{{S}^{1}\cap {S}^{2}}^{k}{\Phi}^{1}{\eta}^{1}\)
Soit:
où:
\({L}_{{S}^{1}\cap {S}^{2}}^{1}\) est la matrice de liaison de \({S}^{1}\) associée à l’interface \({S}^{1}\cap {S}^{2}\) , elle exprime la continuité de déplacement entre les deux sous-structures à partir des déplacements généralisés à l’interface,
\({L}_{{S}^{1}\cap {S}^{2}}^{2}\) est la matrice de liaison de \({S}^{2}\) associée à l’interface \({S}^{1}\cap {S}^{2}\) .
Dans le cas du problème aux valeurs propres de la structure globale, muni de ses conditions aux limites, on peut écrire:
Les matrices de masse et de rigidité généralisées, le vecteur des degrés de liberté généralisés et la matrice de liaison qui apparaissent ici, sont définis sur la structure globale. Ils prennent une forme particulière à chaque méthode (Craig-Bampton, MacNeal, …) qui sera explicitée par la suite. Le vecteur des multiplicateurs de Lagrange \(\lambda\) s’interprète par les actions de liaison auxquelles sont soumises les interfaces.
Il s’agit donc d’un problème classique de recherche de valeurs propres, auquel est associée une équation de contrainte linéaire. Dans Code_Aster , ce type de problème est résolu par double dualisation des conditions aux limites [bib6].
Ainsi, on peut démontrer que ce système est également solution du problème de minimisation de la fonctionnelle suivante, dite intégrale d’action :
\(f={\int}_{a}^{b}F(\eta ,\dot{\eta},{\lambda}_{1}{\lambda}_{2})\cdot \mathrm{dt}\)
où \(\mathrm{F}(\eta ,\dot{\eta},{\lambda}_{1}{\lambda}_{2})=\frac{1}{2}{\eta}^{T}\stackrel{ˉ}{\mathrm{K}}\eta +\frac{1}{2}\dot{\eta}\stackrel{ˉ}{\mathrm{M}}\dot{\eta}+{({\lambda}_{1}+{\lambda}_{2})}^{T}\mathrm{L}\eta -\frac{1}{2}{({\lambda}_{1}-{\lambda}_{2})}^{2}\) est le hamiltonien.
Les variables de cette fonctionnelle sont les coordonnées généralisées \(\eta\) , les vitesses associées et les multiplicateurs de Lagrange \({\lambda}_{1}\) et \({\lambda}_{2}\) (en nombre égal à 2 fois le nombre d’équations de liaison). Le dernier terme de la fonctionnelle impose l’égalité des coefficients de Lagrange.
L’extremum est atteint pour les valeurs des variables qui annulent les dérivées de \(f\) , quels que soient \(a\) et \(b\) réels:
\(-{\lambda}_{1}+L\eta +{\lambda}_{2}=0\)
\({L}^{T}({\lambda}_{1}+{\lambda}_{2})+(\stackrel{ˉ}{K}-\omega \mathrm{²}\stackrel{ˉ}{M})\eta =0\)
\({\lambda}_{1}+L\eta -{\lambda}_{2}=0\)
La double dualisation aboutit donc à un problème matriciel symétrique réel. On démontre [R3.03.01] qu’elle permet de rendre les algorithmes de triangulation de matrice inconditionnellement stables.
Cas du maillage d’interface incompatible#
Dans cette partie, nous considérons le cas de maillages incompatibles. Cela signifie qu’ils ne vérifient à priori pas les propriétés suivantes:
la restriction du maillage de chacune des sous-structures \({S}^{1}\) et \({S}^{2}\) à l’interface s’appuient strictement sur les mêmes nœuds en leur intersection,
les éléments finis associés à ces mailles de liaison sont de même nature (linéaire, quadratique) de part et d’autre de la frontière.
La méthode choisie pour gérer l’incompatibilité de maillage est une méthode dérivée de celle du PPCM («Plus Petit Maillage Commun»). Cette méthode est utilisée dans le cadre de plusieurs types de problèmes dans Code_Aster : gestion de maillage incompatible, projection de champ sur des maillages de natures différentes… On définit par conséquent une interface «maître» sur laquelle on projette les nœuds de l’interface «esclave».
L’interpolation des champs aux nœuds de l’interface esclave (éléments entourés de \({S}_{2}\) sur la figure ci-dessous) se fait sur l’élément de l’interface maître sur lequel est projeté chaque nœud, à l’aide des fonctions de formes de cet élément.
Figure -a:
On définit \({\Sigma}_{e}\) comme l’interface restreinte à un élément et \({n}_{e}\) le nombre de ddl de l’élément.
En choisissant l’interface de la sous-structure \({S}_{1}\) comme interface maître, la relation () se réécrit:
\({u}^{1}(M)-\sum_{j=1}^{{n}_{e}}({\int}_{{\Sigma}_{e}}{E}_{j}^{2}(M)\xi (M)\text{dM}){q}_{j}^{2}=0\)
où \(\xi (M)\) sont les fonctions de formes de l’élément sur lequel le nœud est projeté et \({E}_{j}^{2}\) la restriction des fonctions de formes de l’élément \(j\) de l’interface de \({S}_{2}\) .
Cette méthode permet de coupler des sous-structures possédant des maillages de nature différentes (linéaire avec quadratique,…) et des densités différentes. On retiendra toutefois que la règle de bonne utilisation consiste à choisir pour maître le maillage de densité la plus faible.
Compte tenu des conditions définies aux paragraphes précédents, l’expression liant les degrés de libertés doit être écrite de nouveau. En prenant l’interface de la sous-structure \({S}_{1}\) comme interface maître, on a:
\({q}_{j}^{1}=\sum_{l}^{{n}_{e}}{\alpha}_{l}{q}_{l}^{2}\)
Ou encore:
\({q}^{1}=A{q}^{2}\)
où \(A\) est une matrice d’observation entre les ddl des deux interfaces.
Les relations de continuité des déplacements s’expriment, avec les coordonnées généralisées, par les relations matricielles suivantes:
\({L}_{j}^{1}{\eta}^{1}={\tilde{L}}_{j}^{2}{\eta}^{2}\)
\({\tilde{L}}_{j}^{k}=A{B}^{k}{\Phi}^{k}\)
On peut écrire cette relation sous forme matricielle:
\(\left[{L}_{j}^{1}-{\tilde{L}}_{j}^{2}\right]\cdot \left[\begin{array}{c}{\eta}^{1}\\ {\eta}^{2}\end{array}\right]=0\)
Les termes des matrices de liaisons ne sont plus égaux à zéro ou un. Ils ont issus des interpolations des déplacements projetés sur les éléments.
Enfin, en utilisant la double dualisation, le système matriciel obtenu à résoudre, dans le cas du problème modal sur la structure globale est:
\(\Rightarrow \left[\begin{array}{ccc}-\text{Id}& \tilde{A}L& \text{Id}\\ {\tilde{A}}^{T}{L}^{T}& \stackrel{ˉ}{K}& {\tilde{A}}^{T}{L}^{T}\\ \text{Id}& \tilde{A}L& -\text{Id}\end{array}\right]\left\lbrace \begin{array}{c}{\lambda}_{1}\\ \eta \\ {\lambda}_{2}\end{array}\right\rbrace -{\omega}^{2}\left[\begin{array}{ccc}0& 0& 0\\ 0& \stackrel{ˉ}{M}& 0\\ 0& 0& 0\end{array}\right]\left\lbrace \begin{array}{c}{\lambda}_{1}\\ \eta \\ {\lambda}_{2}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}0\\ 0\\ 0\end{array}\right\rbrace\)
Où \(\tilde{A}\) désigne la matrice identité \(\text{Id}\) si la sous-structure est celle dont l’interface est maître et \(A\) , la matrice d’observation, s’il s’agit de celle dont l’interface est esclave.
La double dualisation aboutit donc à un problème matriciel symétrique réel même pour le cas des maillages incompatibles.
Conclusion#
Cette méthode permet donc de traiter la connexion d’interfaces correspondant à des types de base modale différents sans coût de gestion d’une élimination toujours délicate. D’autre part, elle est relativement simple. L’inconvénient majeur de cette formulation est de conduire à des systèmes assemblés finaux de dimensions plus importantes que dans la cas de l’élimination. En effet, le couplage des équations matricielles a été fait en introduisant un nombre de degrés de liberté supplémentaires égal à deux fois le nombre d’équations de liaison. Cet accroissement des dimensions des matrices peut donc être très important. Notons que les degrés de liberté de Lagrange introduits sont, dans ce cas, les forces appliquées aux interfaces pour assurer la liaison entre les deux sous-structures.
Élimination des contraintes linéaires#
La technique de la double dualisation permet d’obtenir une uniformité dans la prise en compte des contraintes. Cette approche est en effet utilisable aussi bien pour les contraintes linéaires que non linéaires, sur les déplacements ou les contraintes. Elle autorise donc, avec une grande souplesse, la réalisation d’une grande variété de calculs.
Néanmoins, dans le cadre de l’analyse modale en général, et des techniques de sous-structuration en particulier, la présence des doubles multiplicateurs de Lagrange peut constituer un inconvénient, à la fois liée à l’augmentation de la taille du problème avec le nombre de contrainte, mais aussi liée au conditionnement numérique du système ainsi obtenu. L’augmentation de la taille du problème est particulièrement nette dans le cas des approches par sous-structurations, où, pour le problème final, la grande majorité des équations à résoudre concerne les contraintes. Pour contourner cette difficulté, on a mis en œuvre une technique d’élimination des contraintes cinématiques sur la base d’une décomposition QR des matrices de liaison.
Cas de contraintes linéairement indépendantes#
Les contraintes linéaires associées à un problème de dynamique peuvent s’écrire sous la forme du système linéaire:
\(q\) désignant l’ensemble des degrés de libertés du modèle complet (c’est-à-dire incluant toutes les sous structures, dans notre cas). L’idée de l’élimination consiste à construire une base du noyau de \(L\) , et à résoudre le problème projeté dans cette base. La construction de cette base se fait simplement à partir de la décomposition QR de \({L}^{T}\) .
Pour une matrice \(L\) de taille \(n\times N\) , où \(n\) est le nombre de contraintes, et \(N\) le nombre de degrés de liberté, cette décomposition s’écrit:
\({Q}_{2}\) est alors directement la base recherchée, ainsi qu’il est montré dans la suite.
Cas de contraintes redondantes#
Cette décomposition permet également d’analyser ce qui se passe en présence de contraintes redondantes (nœuds contraints plusieurs fois, sous-structure «en étoile», etc.). Lorsque les contraintes ne sont pas toutes indépendantes, on peut écrire \(L\) sous la forme d’une combinaison linéaire de \({n}_{I}<n\) contraintes linéairement indépendantes \({L}_{I}\) , soit:
\({\left[L\right]}_{n\times N}={\left[A\right]}_{n\times {n}_{I}}\cdot {\left[{L}_{I}\right]}_{{n}_{I}\times N}\)
La décomposition de \({L}_{I}\) s’écrit directement:
et la décomposition de \(L\) devient:
\({Q}_{K}\) une base de \(\mathrm{Ker}(L)\) . En effet, pour n’importe quel élément \(q\) du noyau de \(L\) , il existe un vecteur \(y\) tel que \(q=[{Q}_{K}]\cdot y\) . Dans ces conditions, on a:
\(\begin{array}{ccc}[L]\cdot q& =& \left[\begin{array}{c}A\cdot {R}_{I}^{T}\\ 0\end{array}\right]\cdot \left[\begin{array}{c}{Q}_{I}^{T}\\ {Q}_{K}^{T}\end{array}\right]\cdot q\\ & =& \left[\begin{array}{c}A\cdot {R}_{I}^{T}\\ 0\end{array}\right]\cdot \left[\begin{array}{c}{Q}_{I}^{T}\\ {Q}_{K}^{T}\end{array}\right]\cdot \left[{Q}_{K}\right]\cdot y\\ & =& [A\cdot {R}_{I}^{T}]\cdot [{Q}_{I}^{T}]\cdot [{Q}_{K}]\cdot y\\ & =& 0\end{array}\)
puisque \(I(L)\) et \(\mathrm{Ker}(L)\) sont orthogonaux. Lorsque les contraintes de \(L\) sont linéairement indépendantes, on a directement \({Q}_{K}={Q}_{2}\) . Dans le cas contraire, il faut construire \({Q}_{K}\) à partir de \({Q}_{2}\) et des colonnes \(n-{n}_{I}\) dernières colonnes de \({Q}_{1}\) .
Pour identifier le rang \({n}_{I}\) de la matrice \(L\) , on utilise la matrice triangulaire supérieure \({R}_{1}\) , issue de la relation (). D’après la relation (), elle doit en effet se mettre sous la forme:
\({\left[{R}_{1}\right]}_{n\times n}=\left[\begin{array}{c}{[{R}_{I}]}_{{n}_{I}\times {n}_{I}}\cdot {[A]}_{n\times {n}_{I}}^{T}\\ {[0]}_{(n-{n}_{I})\times n}\end{array}\right]\)
Le rang \({n}_{I}\) est le nombre de termes diagonaux non nuls de \({Q}_{1}\) . En conséquence, les vecteurs de \({Q}_{1}\) appartenant au noyau de \(L\) sont dont ceux correspondant aux termes diagonaux nuls de \(R\) . En pratique, la recherche de \({n}_{I}\) est dépendante de la valeur fixée comme seuil pour le zéro numérique.
Cas d’un nombre important de contraintes#
Les principes exposés dans les paragraphes précédents peuvent être implémentés numérique tels quels et conserver des performances correctes si les matrices associées aux contraintes sont de tailles limitées, soit au maximum de l’ordre du millier de contraintes, faisant intervenir de l’ordre de la dizaine de millier de degrés de libertés. Dans le cas contraire, une des méthodes consiste à construire ce sous espace de manière itérative.
La démarche globale consiste à modifier successivement les lignes de la matrice \(Q\) , pour toutes les contraintes de \(L\) . La matrice \(Q\) de départ est la matrice identité. Pour illustrer la démarche, considérons un étape particulière en choisissant un ensemble de contraintes issues de \(L\) . On note \({Q}_{k}\) la matrice prenant en compte l’élimination des \(k\) premières contraintes. La construction de \({Q}_{k+1}\) se fait de la façon suivante:
Extraire la ligne \({L}_{k+1}\) de \(L\) ,
Déterminer le nombre d’éléments non nuls de \({L}_{k+1}\) . Ces éléments non nuls sont associés aux degrés de liberté impliqués dans la relation linéaire du type de ().
Si un seul élément de \({L}_{k+1}\) est non nul, on annule le terme diagonal de \({Q}_{k+1}\) associé à ce degré de liberté.
Si plusieurs éléments non nuls sont détectés, on cherche, à partir des \({L}_{k}\) lignes précédentes, quels sont les degrés de liberté déjà impliqués dans une relation linéaire. On construit la sous-matrice de \({\mathrm{Q}}_{k+1}\) permettant d’éliminer les nouvelles contraintes en tenant compte de la présence des contraintes précédentes.
Pour ce dernier point, on note \({q}_{i}\) les degrés de liberté impliqués dans les \({L}_{k}\) relations linéaires précédentes, et \({q}_{l}\) les autres degrés de liberté. Le sous-ensemble de relation linéaire associé à \({L}_{k+1}\) s’écrit donc \([{L}_{i}{L}_{l}]\left\lbrace \begin{array}{c}{q}_{i}\\ {q}_{l}\end{array}\right\rbrace =0\) , et on cherche le bloc de \({Q}_{k+1}\) qui vérifie \([{L}_{i}{L}_{l}]\left[\begin{array}{cc}{Q}_{ii}& {Q}_{il}\\ {Q}_{li}& {Q}_{ll}\end{array}\right]=0\) . Le bloc \({Q}_{ii}\) a été construit préalablement, et ne doit pas être modifié pour ne pas impacter les \({L}_{k}\) relations précédentes. On construit donc \({Q}_{il}\) , \({Q}_{li}\) et \({Q}_{ll}\) tels que:
\({Q}_{il}=0\) ,
\({Q}_{li}=-{({L}_{l})}^{+}{L}_{i}{Q}_{ii}\) , où \({({L}_{l})}^{+}\) est le pseudo-inverse de \({Q}_{l}\) , permettant de résoudre \({L}_{i}{Q}_{ii}+{L}_{l}{Q}_{li}=0\) ,
\({Q}_{ll}\) est construit par décomposition \(QR\) de \({L}_{l}={Q}_{ll}{R}_{ll}\) . Après décomposition, les colonnes de \({Q}_{ll}\) associées aux des termes diagonaux non nuls de \({R}_{lŀ}\) sont mises à zéros.
Plusieurs choix existent pour la construction successive des \({L}_{k}\) , et consistent soit à extraire un nombre de lignes fixe de \(\mathrm{L}\) pour chaque pas, soit à faire une recherche des blocs de contraintes disjointes dans \(L\) .
Méthodes de sous-structuration dynamique classique#
Introduction#
Après avoir étudié séparément les différentes étapes de la sous-structuration, et les techniques qu’elles mettent en jeu, il paraît intéressant de présenter les principales méthodes de sous‑structuration dynamique: la méthode de Craig-Bampton et celle de Mac Neal et la méthode des modes d’interfaces réduits.
Craig-Bampton utilise, comme base de projection des sous-structures, des modes contraints et des modes normaux à interfaces fixes [bib1].
Par ailleurs, MacNeal utilise, comme base de projection des sous-structures, des modes d’attache et des modes normaux à interfaces libres [bib2].
On abordera une méthode dérivée de ces méthodes nommée méthodes des modes d’interface réduits qui utilise comme base de projection des sous-structures, des modes d’interface et des modes normaux à interfaces bloquées.
Méthode de Craig-Bampton#
La présentation suivante ne fait intervenir que deux sous-structures \({S}^{1}\) et \({S}^{2}\) , mais elle est généralisable à un nombre quelconque de composants. Après avoir étudié séparément chaque sous‑structure, leurs bases de projection (modes normaux à interfaces fixes et modes contraints) sont connues. Pour chacune d’entre elles (identifiée par l’exposant \(k\) ), on établit une partition des degrés de liberté, distinguant le vecteur des degrés de liberté internes \({q}_{i}^{k}\) et le vecteur des degrés de liberté de liaison \({q}_{j}^{k}\) :
\({q}^{k}=\left\lbrace \begin{array}{c}{q}_{i}^{k}\\ {q}_{j}^{k}\end{array}\right\rbrace\)
Soit:
\({\varphi}^{k}\) la matrice des vecteurs propres de la sous-structure \({S}^{k}\) ,
\({\psi}^{k}\) la matrice des modes contraints de la sous-structure \({S}^{k}\) .
La base de projection de \({S}^{k}\) est caractérisée par la matrice:
\({\Phi}^{k}=[{\varphi}^{k}{\psi}^{k}]\)
La transformation de RITZ (équation [éq ]), nous permet d’écrire:
\({\eta}_{i}^{k}\) est le vecteur des degrés de liberté généralisés associés aux modes propres de \({S}^{k}\) ,
\({\eta}_{j}^{k}\) est le vecteur des degrés de liberté généralisés associés aux modes contraints de \({S}^{k}\) .
Or les modes normaux sont déterminés avec des interfaces fixes, et chaque mode contraint est obtenu en imposant un déplacement unitaire à un degré de liberté de liaison, les autres étant bloqués.
Les coordonnées généralisées relatives aux déformées statiques sont alors les valeurs des degrés de liberté de liaison:
\({q}_{j}^{k}={\eta}_{j}^{k}\)
Intéressons-nous à la contribution de la sous-structure \({S}^{k}\) d’un point de vue énergétique. Les énergies cinétique et de déformation sont:
\({T}^{k}=\frac{1}{2}{q}^{{k}^{T}}.{M}^{k}.{\dot{q}}^{k}=\frac{1}{2}{({\Phi}^{k}{\dot{\eta}}^{k})}^{T}.{M}^{k}.{\Phi}^{k}{\dot{\eta}}^{k}=\frac{1}{2}{\dot{\eta}}^{k}.{\stackrel{ˉ}{M}}^{k}.{\dot{\eta}}^{k}\)
\({U}^{k}=\frac{1}{2}{q}^{{k}^{T}}.{K}^{k}.{\dot{q}}^{k}=\frac{1}{2}{({\Phi}^{k}{\dot{\eta}}^{k})}^{T}.{K}^{k}.{\Phi}^{k}{\dot{\eta}}^{k}=\frac{1}{2}{\dot{\eta}}^{k}.{\stackrel{ˉ}{K}}^{k}.{\dot{\eta}}^{k}\)
Ces expressions font apparaître les projections des matrices de masse et de rigidité sur la base de la sous-structure. Ces matrices, dites généralisées, vérifient un certain nombre de propriétés:
du fait de l’orthogonalité des modes normaux par rapport aux matrices de rigidité et de masse, le bloc supérieur gauche de ces matrices est diagonal. Par la suite, nous considérerons que ces modes sont normés par rapport à la matrice de masse,
on peut également démontrer que les modes contraints sont orthogonaux aux modes normaux par rapport à la matrice de rigidité [bib4].
Les matrices de rigidité et de masse généralisées ont donc la forme suivante:
où \({\lambda}^{k}\) est la matrice des rigidités généralisées associées aux modes propres de \({S}^{k}\) .
Le choix de la base de projection à interfaces bloquées conduit donc à un couplage des modes normaux et des déformées statiques par la matrice de masse.
Cette méthode est intéressante si l’on envisage une étude numérique des sous-structures. En effet, les modes normaux à interfaces fixes et les modes contraints se prêtent bien au calcul. En revanche, leur détermination expérimentale est délicate, car il est difficile de réaliser un encastrement expérimental de bonne qualité.
On montre en outre dans [bib 1] que cette méthode, est d’ordre 2 en \(\omega /{\omega}_{m}\) où \({\omega}_{m}\) est la plus grande pulsation propre identifiée.
Méthode de Mac Neal#
Il serait possible de présenter cette méthode de façon similaire à celle de Craig-Bampton, l’unique différence résidant dans l’utilisation des modes normaux à interfaces libres et des modes d’attache. Cependant, il paraît intéressant d’adopter une démarche légèrement différente, qui permet d’aboutir à un critère de troncature, et de faire apparaître les modes d’attache comme la contribution statique des modes propres non identifiés.
Comme précédemment, on considère le problème d’une structure à 2 composants. La méthode est généralisable à un nombre quelconque de sous-structures. Dans la suite nous identifierons toute grandeur associée à la sous-structure \({S}^{k}\) par l’exposant \(k\) .
Nous nous limitons, ici, au calcul modal de la structure globale, donc les forces extérieures sont nulles, afin d’alléger les développements théoriques. Aussi, le déplacement de \({S}^{k}\) vérifie l’équation d’équilibre dynamique suivante:
Nous allons considérer que la base de projection de \({S}^{k}\) est constituée de tous ses modes propres à interfaces libres. Ainsi, la dimension du problème projeté est égale à la dimension du problème issu de la modélisation éléments finis. Nous supposons, en outre, qu’un certain nombre de ces modes a été déterminé, les autres étant inconnus:
\({q}^{k}=[{\varphi}_{1}^{k}{\varphi}_{2}^{k}]\left\lbrace \begin{array}{c}{\eta}_{1}^{k}\\ {\eta}_{2}^{k}\end{array}\right\rbrace ={\Phi}^{k}{\eta}^{k}\)
où:
\({\varphi}_{1}^{k}\) |
est la matrice des vecteurs modaux identifiés de la sous-structure \({S}^{k}\) , |
\({\varphi}_{2}^{k}\) |
est la matrice des vecteurs modaux non identifiés de la sous-structure \({S}^{k}\) , |
\({\eta}_{1}^{k}\) |
est le vecteur des degrés de liberté généralisés associés aux modes propres identifiés de \({S}^{k}\) , |
\({\eta}_{2}^{k}\) |
est le vecteur des degrés de liberté généralisés associés aux modes propres non identifiés de \({S}^{k}\) . |
L’équation [éq ] devient, avec les coordonnées généralisées définies précédemment:
\(({\Phi}^{{k}^{T}}.{K}^{k}.{\Phi}^{k}-{\omega}^{2}{\Phi}^{{k}^{T}}.{M}^{k}.{\Phi}^{k}){\eta}^{k}={\Phi}^{{k}^{T}}{f}_{L}^{k}\)
soit:
Les modes propres sont orthogonaux par rapport aux matrices de masse et de rigidité et nous choisissons de les normer par rapport à la matrice de masse. On a donc:
Considérons maintenant l’ensemble des deux sous-structures. Chacune d’entre elles vérifie les équations [éq ] et [éq ]. L’ensemble de ces équations dynamiques constitue le système suivant:
En regroupant les modes identifiés et les modes non identifiés:
L’équation traduisant la transformation de RITZ devient:
Avec ces notations, le système d’équations dynamiques [éq ] devient:
Ce système d’équations traduit le comportement dynamique des sous-structures séparément. Il ne représente pas les mouvements d’ensemble de la structure. Pour cela, il faut lui adjoindre les conditions de liaison entre les deux composants.
Les équations entre sous-structures qui assurent leur liaisonnement dérivent des équations [éq ] et [éq ], ainsi que de l’organisation de la base que nous avons choisie [éq ]:
Les équations [éq ] et [éq ] nous permettent d’exprimer les coordonnées généralisées relatives aux modes non identifiés:
A partir des équations [éq ] et [éq ], on obtient donc:
Or, d’après la formule [éq ], on sait que l’on peut écrire les matrices de liaison, sous la forme:
On a donc, d’après [éq ] et [éq ]:
On voit apparaître la matrice de flexibilité dynamique résiduelle associée aux modes non identifiés:
Ce terme montre que les efforts de liaison aux degrés de liberté d’interface sont non nuls si on ne tient compte dans la modélisation que des modes connus.
Dès lors, le problème matriciel [éq ] peut se ramener au système équivalent suivant:
Ce problème a pour inconnues les coordonnées généralisées associées aux modes propres identifiés et les forces de liaison appliquées à la première sous-structure. Le terme de flexibilité résiduelle, qui apparaît également, n’est pas connu. On propose ci-dessous une modélisation qui fait intervenir les modes d’attache de MacNeal.
Il y a deux cas selon que l’on prend en compte ou non la matrice de flexibilité résiduelle.
Premier cas#
On néglige la matrice de flexibilité dynamique résiduelle associée aux modes non identifiés:
\({R}_{e}(\omega )=0\)
La méthode de sous-structuration dynamique résultante est donc très simple. Elle a le désavantage de s’appuyer sur une méthode de recombinaison modale très sensible aux effets de troncature.
Deuxième cas#
On utilise un développement limité: \({R}_{e}(\omega )={R}_{e}(0)+O({\omega}^{2}/{\omega}_{m}^{2})\) .
Adoptons les notations suivantes:
soit:
\(n\) |
le nombre de modes de la base complète, |
\(m\) |
le nombre de modes non identifiés, |
\(\Phi\) |
la matrice des \(n\) modes normaux à interfaces libres, |
\(\Psi\) |
la matrice des \(n\) modes d’attache (définie pour tous les degrés de liberté du système, pas seulement pour les degrés de liberté d’interface entre les 2 sous-structures), |
\({\lambda}_{i}={\omega}_{i}^{2}\) |
la matrice diagonale des \(n\) valeurs propres. |
La matrice complète des modes d’attache vérifie: \(K\Psi =\mathrm{Id}\text{}\iff \text{}\Psi ={K}^{-1}\) .
La base modale complète des modes normaux à interfaces libres constitue une base orthonormale. La matrice de raideur, exprimée dans cette base s’écrit:
\(\stackrel{ˉ}{K}={\Phi}^{T}.K.\Phi =\lambda\)
De la même manière:
\({\stackrel{ˉ}{K}}^{-1}={\Phi}^{T}.{K}^{-1}.\Phi ={\lambda}^{-1}\)
On en déduit donc une nouvelle expression de la matrice complète des modes d’attache:
\(\Psi ={K}^{-1}=\Phi .{\lambda}^{-1}.{\Phi}^{T}\)
Comme les modes propres sont orthogonaux deux à deux et la matrice des valeurs propres est diagonale, la matrice complète des modes d’attache prend la forme finale:
\(\Psi =\sum_{i=1}^{n}{\varphi}_{i}{\lambda}_{i}^{-1}{\varphi}_{i}^{T}\)
Considérons maintenant la matrice de flexibilité dynamique résiduelle, issue de la méthode de MacNeal:
\({\mathrm{R}}_{e}(\omega )={\Phi}_{2}{({\lambda}_{2}-{\omega}^{2}\mathrm{Id})}^{-1}{\Phi}_{2}^{T}=\sum_{i=m+1}^{n}{\varphi}_{i}{({\lambda}_{i}-{\omega}^{2})}^{-1}{\varphi}_{i}^{T}\)
Lorsque le nombre de modes identifiés est suffisamment important, la contribution dynamique devient négligeable devant la contribution statique:
\(\frac{{\omega}^{2}}{{\omega}_{i}^{2}}\ll 1\Rightarrow {R}_{e}(\omega )\approx {R}_{e}(0)=\sum_{i=m+1}^{n}{\varphi}_{i}{\lambda}_{i}^{-1}{\varphi}_{i}^{T}\)
où \({R}_{e}(0)\) est la matrice de flexibilité statique résiduelle.
Pour mettre en évidence l’effet de troncature modale, nous pouvons approcher cette matrice par son développement à l’ordre 1:
\({R}_{e}(\omega )\approx \sum_{i=m+1}^{n}{\varphi}_{i}{\lambda}_{i}^{-1}(1+\frac{{\omega}^{2}}{{\omega}_{i}^{2}}){\varphi}_{i}^{T}\)
Cette matrice peut être calculée en fonction de la matrice des modes d’attache:
\({R}_{e}(0)=\sum_{i=1}^{n}{\varphi}_{i}{\lambda}_{i}^{-1}{\varphi}_{i}^{T}-\sum_{i=1}^{m}{\varphi}_{i}{\lambda}_{i}^{-1}{\varphi}_{i}^{T}\Rightarrow {R}_{e}(0)=\Psi -\sum_{i=1}^{m}{\varphi}_{i}{\lambda}_{i}^{-1}{\varphi}_{i}^{T}\)
Le second terme de cette formulation est calculable de façon exacte, puisqu’il ne fait intervenir que les modes (à interfaces libres) identifiés.
Enfin, notons que dans la méthode de MacNeal, seule la contribution des modes d’attache aux nœuds d’interface est nécessaire.
La résolution du système [éq 3.3-4] nous permet de déterminer les fréquences propres de la structure globale et les coordonnées généralisées des modes propres. La remontée à l’expression des modes propres dans les bases physiques des sous-structures se fait par la relation suivante:
\({q}_{1}^{k}={\Phi}_{1}^{k}{\eta}_{1}^{k}+{R}_{e}^{k}(0){B}_{1}^{{k}^{T}}{f}_{{L}_{1}}^{T}\)
La méthode de MacNeal aboutit donc, dans le cas du calcul des modes de la structure globale, à un problème aux valeurs propres de taille réduite. Les matrices de masse et de rigidité issues de la sous-structuration sont symétriques. Deux méthodes sont en réalité proposées, selon que l’on prend en compte ou non la flexibilité résiduelle. La littérature sur le sujet tend à montrer que l’usage de la flexibilité résiduelle est indispensable pour obtenir des résultats fiables [bib2].
Méthode des modes d’interface (méthode dite réduite)#
Introduction#
La méthode des modes d’interface est basée sur les travaux de différents auteurs, on trouvera une justification théorique complète dans [bib12]. Nous ne développerons pas tous ces aspects dans ce document car l’utilisation qui en est faite dans Code_Aster est relativement simplifiée: il s’agit juste de diminuer la taille du système réduit par l’utilisation des propriétés spectrales de l’opérateur dynamique condensé à l’interface.
Définition des modes d’interface#
Après avoir étudié chaque sous‑structure séparément pour construire les bases de modes normaux, il s’agit de considérer le problème couplé pour construire les bases de modes d’interface. La procédure actuellement présente dans Code_Aster nécessite de considérer le problème couplé par la sous-structuration dynamique. On utilise donc la méthode de Craig-Bampton sans les modes normaux à interfaces bloquées, ce qui équivaut donc à une méthode de sous-structuration statique par condensation de Guyan. On obtient le système couplé suivant:
\({\stackrel{ˉ}{\mathrm{K}}}^{k}={\psi}^{{k}^{T}}.{\mathrm{K}}^{k}.{\psi}^{k}\text{}{\stackrel{ˉ}{\mathrm{M}}}^{k}={\psi}^{{k}^{T}}.{\mathrm{M}}^{k}.{\psi}^{k}\)
\(\psi\) désigne ici les modes contraints. Le calcul des modes propres de la structure globale condensée à ses interfaces par les méthodes de sous-structuration dynamique consiste à résoudre un problème aux valeurs propres matriciel suivant:
\((\stackrel{ˉ}{K}-{\omega}^{2}\stackrel{ˉ}{M}){\eta}_{R}+{L}^{T}\lambda =0\)
\(L{\eta}_{R}=0\)
Il va de soit que lors de ce calcul des modes d’interface, cela n’a pas de sens de calculer autant de modes que de degrés de liberté d’interface car cela équivaudrait à revenir sur le système de modes contraints. Une erreur de troncature de la base de projection est donc introduite à ce niveau, mais qui reste physique et maîtrisable.
Les modes propres obtenus ne sont définis que sur les degrés de libertés d’interface, il est donc nécessaire d’étendre ces modes sur chacune des sous-structures couplées. On obtient ainsi:
\({\psi}_{R}^{k}={\psi}^{k}{\eta}_{R}\)
Calcul des modes d’interfaces#
Dans la pratique, le calcul des modes d’interfaces avec la méthode présentée au paragraphe 3.4.2 n’est pas réaliste. En effet, l’opération de calcul des modes de contraintes \(\psi\) est très couteuse en mémoire et en temps de calcul, d’autant plus qu’on ne considère ensuite qu’une sous famille de \(N\) vecteurs, \(N\) étant très faible devant le nombre de degrés de liberté de l’interface. L’approche adoptée pour la commande MODE_STATIQUE [U4.52.14] consiste à construire les modes d’interfaces à l’aide d’un pré-conditionneur. Si on note \(i\) les degrés de liberté d’interface, et \(c\) les degrés de liberté complémentaire, en partitionnant la matrice de raideur \(\mathrm{K}\) , les modes de contrainte \(\psi\) vérifient
\(\left\lbrace \begin{array}{c}{\psi}_{i}\\ {\psi}_{\mathrm{c}}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}\mathit{Id}\\ -{\mathrm{K}}_{\mathrm{cc}}^{-1}{\mathrm{K}}_{\mathrm{ci}}\end{array}\right\rbrace\)
Plutôt que de construire directement les mode de contrainte sur ce principe, on assemble construit des matrices de masse \(\tilde{{\mathrm{M}}_{ii}}\) et de de raideur \(\tilde{{\mathrm{K}}_{ii}}\) définies uniquement à l’interface. Ces matrices sont construites sur la base d’un modèle de poutres d’Euler Bernoulli en reprenant la connectivité des nœuds de l’interface. Chaque paire de nœuds connectée de l’interface est alors reliée par une poutre cylindrique dont les propriétés géométriques sont ajustées par rapport à la plus courte distance entre deux nœuds connectés de l’interface. Le diamètre de la poutre est alors fixé à 20% de cette distance minimale. Le matériau utilisé pour les poutre est un acier. On calcul ensuite une famille de mode \(\tilde{{\phi }_{i}}\) de ce modèle d’interface. Il convient ici de calculer un nombre de mode \(M\) assez sensiblement supérieur à \(N\) , où \(N\) est le nombre de modes d’interface que l’utilisateur souhaite retenir.
Pour calculer ces modes, on réalise une première opération qui permet de déterminer le nombre de sous parties disjointes d’une interface donnée. Pour chaque partie disjointe d’interface, on calcule les 6 modes de corps rigide. Si on note \(d\) le nombre de parties disjointes, alors le nombre total \(M\) de mode calculé vaut:
\(M=3(E(\frac{N}{d})+2E(\frac{6}{d})+1)\) si \(E(\frac{N}{d})>7\)
\(M=3(6+E(\frac{N}{d})+2)\)
où \(E(.)\) correspond à la partie entière de la quantité considérée.
Ces \(M\) modes \(\tilde{{\phi }_{i}}\) sont ensuite relevés statiquement sur le modèle de la sous structure, et forment alors la base sur laquelle on recherche les modes d’interface:
\(\tilde{\psi}=\left\lbrace \begin{array}{c}\tilde{{\psi}_{i}}\\ \tilde{{\psi}_{\mathrm{c}}}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}\tilde{{\phi }_{i}}\\ -{\mathrm{K}}_{\mathrm{cc}}^{-1}{\mathrm{K}}_{\mathrm{ci}}\tilde{{\phi }_{i}}\end{array}\right\rbrace\)
On projette ensuite les matrices de masse et de raideur du macro élément sur ces vecteurs,
\(\stackrel{ˉ}{\mathrm{K}}={\tilde{\psi}}^{T}.\mathrm{K}.\tilde{\psi}\text{}\stackrel{ˉ}{\mathrm{M}}={\tilde{\psi}}^{T}.\mathrm{M}.\tilde{\psi}\)
et on en cherche les modes propres:
\((\stackrel{ˉ}{\mathrm{K}}-{\omega}^{2}\stackrel{ˉ}{\mathrm{M}}){\eta}_{R}=0\)
Les modes d’interface retenus pour la construction du macro éléments sont finalement les \(N\) premiers vecteurs définis par:
\(\psi =\tilde{\psi}{\eta}_{\mathrm{R}}\)
Le calcul des matrices \(\tilde{{\mathrm{M}}_{ii}}\) et \(\tilde{{\mathrm{K}}_{ii}}\) , et des modes associés \(\tilde{{\phi }_{i}}\) est rapide et peu couteux en mémoire, et permet de limiter l’opération de relèvement aux seuls vecteurs d’intérêt pour construire les modes d’interface. Naturellement, si on calcule autant de mode \(\tilde{{\phi }_{i}}\) que de DDL d’interface, alors le sous espace décrit est identique à celui décrit par les modes de contrainte
Création des macro-éléments réduits#
Pour chacune des sous-structures, on dispose des bases de modes normaux et des bases des modes d’interface. La création des macro-éléments est alors identique au cas de la méthode de Craig-Bampton. On établit une partition des degrés de liberté, distinguant le vecteur des degrés de liberté internes \({q}_{i}^{k}\) et le vecteur des degrés de liberté de liaison \({q}_{j}^{k}\) :
\({q}^{k}=\left\lbrace \begin{array}{c}{q}_{i}^{k}\\ {q}_{j}^{k}\end{array}\right\rbrace\)
Soit:
\({\phi }^{k}\) la matrice des vecteurs propres de la sous-structure \({S}^{k}\) ,
\({\psi}_{R}^{k}\) la matrice des modes d’interface de la sous-structure \({S}^{k}\) .
La base de projection de \({S}^{k}\) est caractérisée par la matrice:
\({\Phi}^{k}=[{\phi }^{k}{\psi}_{R}^{k}]\)
La transformation de RITZ (équation [éq ]), nous permet d’écrire:
\({\eta}_{i}^{k}\) est le vecteur des degrés de liberté généralisés associés aux modes propres de \({S}^{k}\) ,
\({\eta}_{j}^{k}\) est le vecteur des degrés de liberté généralisés associés aux modes d’interface de \({S}^{k}\) .
Intéressons nous à la contribution du composant \({S}^{k}\) d’un point de vue énergétique. Les énergies cinétique et de déformation sont:
\({T}^{k}=\frac{1}{2}{q}^{{k}^{T}}.{M}^{k}.{\dot{q}}^{k}=\frac{1}{2}{({\Phi}^{k}{\dot{\eta}}^{k})}^{T}.{M}^{k}.{\Phi}^{k}{\dot{\eta}}^{k}=\frac{1}{2}{\dot{\eta}}^{k}.{\stackrel{ˉ}{M}}^{k}.{\dot{\eta}}^{k}\)
\({U}^{k}=\frac{1}{2}{q}^{{k}^{T}}.{K}^{k}.{\dot{q}}^{k}=\frac{1}{2}{({\Phi}^{k}{\dot{\eta}}^{k})}^{T}.{K}^{k}.{\Phi}^{k}{\dot{\eta}}^{k}=\frac{1}{2}{\dot{\eta}}^{k}.{\stackrel{ˉ}{K}}^{k}.{\dot{\eta}}^{k}\)
Ces expressions font apparaître les projections des matrices de masse et de rigidité sur la base de la sous-structure. Ces matrices, dites généralisées, vérifient un certain nombre de propriétés du fait de l’orthogonalité des modes normaux par rapport aux matrices de rigidité et de masse, le bloc supérieur gauche de ces matrices est diagonal. De plus, nous considérerons que ces modes sont normés par rapport à la matrice de masse.
Les matrices de rigidité et de masse généralisées ont donc la forme suivante:
où \({\lambda}^{k}\) est la matrice des rigidités généralisées associées aux modes propres de \({S}^{k}\) .
Mise en œuvre dans Code_Aster#
Étude des sous-structures séparément#
La base de projection de chaque sous-structure est composée de modes propres dynamiques et de déformées statiques.
Les modes propres dynamiques de la sous-structure sont calculés avec l’opérateur classique du Code_Aster: CALC_MODES [U4.52.02]. Dans le cas de la sous‑structuration de Craig-Bampton, les interfaces de liaison doivent être bloquées. Ceci est réalisé avec l’opérateur AFFE_CHAR_MECA [U4.44.01].
L’opérateur DEFI_INTERF_DYNA [U4.64.01] permet de définir les interfaces de connexion de la sous‑structure. En particulier, on précise le type de l’interface, qui peut être soit « CRAIGB » (Craig‑Bampton), soit « MNEAL » (MacNeal), soit enfin « AUCUN » (pas de modes statiques calculés).
L’opérateur DEFI_BASE_MODALE [U4.64.02] permet de calculer la base de projection complète de la sous‑structure. Ainsi, les modes dynamiques calculés précédemment sont recopiés. Par ailleurs, les déformées statiques sont calculées en fonction du type défini dans l’opérateur DEFI_INTERF_DYNA [U4.64.01]. Si le type est « CRAIGB », on calcule les modes contraints des interfaces de la sous‑structure. Si le type est « MNEAL », on calcule les modes d’attache des interfaces de la sous‑structure. Si le type est « AUCUN », on ne calcule pas de déformée statique, ce qui correspond à une base de type MacNeal sans correction statique.
L’opérateur MACR_ELEM_DYNA [U4.65.01] calcule les matrices généralisées de rigidité et de masse de la sous‑structure, ainsi que les matrices de liaison.
Assemblage#
Le modèle de la structure complète est déterminé par l’opérateur DEFI_MODELE_GENE [U4.65.02]. En particulier, chaque sous-structure est définie par le macro-élément qui lui correspond (issu de MACR_ELEM_DYNA) et les angles de rotation qui permettent de l’orienter. Les liaisons entre sous‑structures sont définies par la donnée des noms des deux sous-structures impliquées et de ceux des deux interfaces en vis à vis. Dans le cas de la méthode des modes d’interface, il est nécessaire de spécifier l’option ‘REDUIT’ afin de permettre le couplage entre les macro-éléments.
La numérotation du problème généralisé complet est réalisée par l’opérateur NUME_DDL_GENE [U4.65.03]. Les matrices de masse et de raideur généralisées de la structure complète sont assemblées en fonction de cette numérotation avec l’opérateur ASSE_MATR_GENE [U4.65.04].
Calcul modal par sous structuration dynamique classique#
Introduction#
Pour les problèmes modaux, le système étudié n’est soumis à aucune force extérieure. On ne tient pas compte de la dissipation dans le solide. La méthode de calcul modal par sous-structuration, programmée dans Code_Aster , qui permet de remplacer le problème global par un problème simplifié, procède en quatre temps. Tout d’abord, des modes propres et des déformées statiques sont calculés sur chacune des sous-structures composant le système. Ensuite, le problème global est projeté sur ces champs, et on tient compte des couplages entre les sous-structures, au niveau de leurs interfaces. On peut alors résoudre classiquement le problème réduit obtenu. Finalement, il ne reste plus qu’à en déduire la solution d’ensemble par reconstitution.
Équations dynamiques vérifiées par les sous-structures séparément#
Nous allons considérer une structure \(S\) composée de \({N}_{S}\) sous-structures notées \({S}^{k}\) . Nous supposons que chaque sous-structure est modélisée en éléments finis. Nous avons vu que dans un calcul par sous-structuration dynamique, le comportement vibratoire des sous-structures résulte des forces extérieures qui leur sont appliquées, et des forces de liaison qu’exercent sur elles les autres sous-structures. Ainsi, au niveau de la sous-structure \({S}^{k}\) , et dans le cas du calcul modal nous pouvons écrire:
Le champ des forces de liaison s’écrit en notation complexe:
\({f}_{L}^{k}(t)=\left\lbrace {f}_{L}^{k}\right\rbrace {e}^{j\omegat }\)
Le champ des déplacements s’écrit:
\({q}^{k}(t)=\left\lbrace {q}^{k}\right\rbrace {e}^{j\omegat }\)
Les champs de vitesse et d’accélération s’écrivent:
\({\dot{q}}^{k}(t)=j\omega \left\lbrace {q}^{k}\right\rbrace {e}^{j\omegat }\)
\({\ddot{q}}^{k}(t)=-{\omega}^{2}\left\lbrace {q}^{k}\right\rbrace {e}^{j\omegat }\)
Finalement la sous-structure \({S}^{k}\) vérifie l’équation suivante:
La méthode de synthèse modale consiste à rechercher le champ de déplacement inconnu, issu de la modélisation éléments finis, sur un espace approprié, de dimension réduite (transformation de Ritz). Nous avons vu que pour chaque sous-structure, cet espace est composé de modes propres dynamiques et de déformées statiques:
\({\phi }^{k}\) sont les vecteurs modaux associés aux modes propres dynamiques de \({S}^{k}\) ,
\({\psi}^{k}\) sont les vecteurs modaux associés aux déformées statiques de \({S}^{k}\) ,
\({\eta}_{i}^{k}\) est le vecteur des coordonnées généralisées associées aux modes propres de \({S}^{k}\) ,
\({\eta}_{j}^{k}\) est le vecteur des coordonnées généralisées associées aux déformées statiques de \({S}^{k}\) ,
\({\eta}^{k}\) est le vecteur des coordonnées généralisées de \({S}^{k}\) .
L’équation [éq ] est projetée sur la base de \({S}^{k}\) en tenant compte de [éq ]. Ceci nous permet d’écrire:
En supposant que les modes propres dynamiques et les déformées statiques sont organisés comme le montre la formule [éq ] et en considérant que les vecteurs propres associés aux modes dynamiques sont normés par rapport à la masse modale unitaire, les matrices de masse et de rigidité généralisées prennent la forme suivante:
\({\stackrel{ˉ}{M}}^{k}=\left[\begin{array}{cc}\mathrm{Id}& {\varphi}^{{k}^{T}}{M}^{k}{\psi}^{k}\\ {\psi}^{{k}^{T}}{M}^{k}{\varphi}^{k}& {\psi}^{{k}^{T}}{M}^{k}{\psi}^{k}\end{array}\right]{\stackrel{ˉ}{K}}^{k}=\left[\begin{array}{cc}{\lambda}^{k}& {\varphi}^{{k}^{T}}{K}^{k}{\psi}^{k}\\ {\psi}^{{k}^{T}}{K}^{k}{\varphi}^{k}& {\psi}^{{k}^{T}}{K}^{k}{\psi}^{k}\end{array}\right]\)
où:
\(\mathrm{Id}\) est la matrice Identité,
\({\lambda}^{k}\) est la matrice diagonale des carrés des pulsations propres de la base.
On démontre, dans le cas de la méthode de Craig-Bampton, que les modes normaux et les modes contraints sont orthogonaux vis-à-vis de la matrice de rigidité dont les termes extra-diagonaux sont, dès lors, nuls. Cependant, cette propriété n’est pas utilisée dans l’algorithme programmé dans Code_Aster .
Équations dynamiques vérifiées par la structure globale#
Les équations dynamiques que vérifie la structure globale sont:
\((\left[\begin{array}{ccccc}{\stackrel{ˉ}{\stackrel{ˉ}{K}}}^{1}& & & & \\ & {}^{...}& & & \\ & & {\stackrel{ˉ}{\stackrel{ˉ}{K}}}^{k}& & \\ & & & {}^{...}& \\ & & & & {\stackrel{ˉ}{\stackrel{ˉ}{K}}}^{{N}_{s}}\end{array}\right]-{\omega}^{2}\left[\begin{array}{ccccc}{\stackrel{ˉ}{\stackrel{ˉ}{M}}}^{1}& & & & \\ & {}^{...}& & & \\ & & {\stackrel{ˉ}{\stackrel{ˉ}{M}}}^{k}& & \\ & & & {}^{...}& \\ & & & & {\stackrel{ˉ}{\stackrel{ˉ}{M}}}^{{N}_{s}}\end{array}\right])\left\lbrace \begin{array}{c}{\eta}^{1}\\ {}^{...}\\ {\eta}^{k}\\ {}^{...}\\ {\eta}^{{N}_{s}}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{\stackrel{ˉ}{f}}_{L}^{1}\\ {}^{...}\\ {\stackrel{ˉ}{f}}_{L}^{k}\\ {}^{...}\\ {\stackrel{ˉ}{f}}_{L}^{{N}_{s}}\end{array}\right\rbrace\)
Auxquelles, il faut ajouter les équations de liaison (d’après [éq ]):
\(\forall k,l\text{}{L}_{{S}^{k}\cap {S}^{l}}^{k}\left\lbrace {\eta}^{k}\right\rbrace ={L}_{{S}^{k}\cap {S}^{l}}^{l}\left\lbrace {\eta}^{l}\right\rbrace\)
Ce système est résolu par double dualisation des conditions aux limites [R3.03.01]. Sa formulation finale fait donc intervenir le vecteur des multiplicateurs de Lagrange \(\lambda\) et peut s’écrire sous la forme condensée:
Le problème défini par l’équation [éq ] est symétrique. D’autre part, sa dimension est déterminée par le nombre de modes pris en considération (modes dynamiques et déformées statiques). On est donc amené à résoudre un problème modal classique, de taille réduite, auquel est associée une équation de contrainte linéaire. Sa résolution ne pose donc pas de problème.
Mise en œuvre dans Code_Aster#
Étude des sous-structures séparément#
Les modes propres dynamiques sont calculés avec l’opérateur CALC_MODES [U4.52.02]. Les conditions aux interfaces de liaison sont appliquées avec l’opérateur AFFE_CHAR_MECA [U4.44.01]. L’opérateur DEFI_INTERF_DYNA [U4.64.01] permet de définir les interfaces de connexion de la sous-structure. L’opérateur DEFI_BASE_MODALE [U4.64.02] permet de calculer la base de projection complète de la sous-structure.
L’opérateur MACR_ELEM_DYNA [U4.65.01] calcule les matrices généralisées de rigidité, de masse et éventuellement d’amortissement de la sous-structure, ainsi que les matrices de liaison.
Assemblage et résolution#
Le modèle de la structure complète est défini par l’opérateur DEFI_MODELE_GENE [U4.65.02]. Sa numérotation est réalisée par l’opérateur NUME_DDL_GENE [U4.65.03]. Les matrices de masse, de rigidité et éventuellement d’amortissement généralisées de la structure complète sont assemblées en fonction de cette numérotation avec l’opérateur ASSE_MATR_GENE [U4.65.04].
Le calcul des modes propres de la structure complète est réalisé par l’opérateur CALC_MODES [U4.52.02].
Restitution sur base physique#
La restitution des résultats sur base physique est identique au cas du calcul modal. Elle fait intervenir l’opérateur REST_GENE_PHYS [U4.63.31].
Pour diminuer la durée des traitements graphiques lors des visualisations, il est possible de créer un maillage grossier par l’opérateur DEFI_SQUELETTE [U4.24.01]. Ce maillage, ignoré pendant le calcul, sert de support aux visualisations.
Réponse harmonique par sous structuration dynamique classique#
Introduction#
Pour les problèmes harmoniques, le système étudié est soumis à une force spatialement quelconque, mais sinusoïdale dans le temps. La forme du chargement, la fréquence d’excitation et les propriétés modales, jouent chacune un rôle essentiel. Il faut également tenir compte de la dissipation dans le solide, que l’on peut traduire par l’introduction d’une matrice d’amortissement. La méthode de calcul harmonique par sous-structuration, programmée dans Code_Aster , qui permet de remplacer le problème global par un problème simplifié, procède en quatre temps.
Tout d’abord, des modes propres et des déformées statiques sont calculés sur chacune des sous-structures composant le système. Ensuite, le problème global est projeté sur ces champs, et on tient compte des couplages entre les sous-structures, au niveau de leurs interfaces. On peut alors résoudre classiquement le problème réduit obtenu. Finalement, il ne reste plus qu’à en déduire la solution d’ensemble par reconstitution.
Équations dynamiques vérifiées par les sous-structures séparément#
Nous allons considérer une structure \(S\) composée de \({N}_{S}\) sous-structures notées \({S}^{k}\) . Nous supposons que chaque sous-structure est modélisée en éléments finis. Nous avons vu que dans un calcul par sous-structuration dynamique, le comportement vibratoire des sous-structures résulte des forces extérieures qui lui sont appliquées, et des forces de liaison qu’exercent sur elles les autres sous-structures. Ainsi, au niveau de la sous-structure \({S}^{k}\) , nous pouvons écrire:
Dans un problème harmonique, on impose un chargement dynamique, spatialement quelconque, mais sinusoïdal de pulsation \(\omega\) dans le temps. On s’intéresse alors à la réponse stabilisée du système, sans tenir compte de la partie transitoire.
Le champ des forces extérieures s’écrit:
\({f}_{\mathrm{ext}}^{k}(t)=\left\lbrace {f}_{\mathrm{ext}}^{k}\right\rbrace {e}^{j\omega t}\)
Le champ des forces de liaison s’écrit:
\({f}_{L}^{k}(t)=\left\lbrace {f}_{L}^{k}\right\rbrace {e}^{j\omega t}\)
Le champ des déplacements s’écrit:
\({q}^{k}(t)=\left\lbrace {q}^{k}\right\rbrace {e}^{j\omegat }\)
Les champs de vitesse et d’accélération s’écrivent:
\({\dot{q}}^{k}(t)=j\omega \left\lbrace {q}^{k}\right\rbrace {e}^{j\omegat }\)
\({\ddot{q}}^{k}(t)=-{\omega}^{2}\left\lbrace {q}^{k}\right\rbrace {e}^{j\omegat }\)
Finalement la sous-structure \({S}^{k}\) vérifie l’équation du mouvement suivante:
La méthode de synthèse modale consiste à rechercher le champ de déplacement inconnu, issu de la modélisation éléments finis, sur un espace approprié, de dimension réduite (transformation de Ritz). Nous avons vu que pour chaque sous-structure, cet espace est composé de modes propres dynamiques et de déformées statiques:
\({\varphi}^{k}\) sont les vecteurs modaux associés aux modes propres dynamiques de \({S}^{k}\) ,
\({\psi}^{k}\) sont les vecteurs modaux associés aux déformées statiques de \({S}^{k}\) ,
\({\eta}_{i}^{k}\) est le vecteur des coordonnées généralisées associées aux modes propres de \({S}^{k}\) ,
\({\eta}_{j}^{k}\) est le vecteur des coordonnées généralisées associées aux déformées statiques de \({S}^{k}\) ,
\({\eta}^{k}\) est le vecteur des coordonnées généralisées de \({S}^{k}\) .
L’équation [éq ] est projetée sur la base de \({S}^{k}\) en tenant compte de [éq ]. Ceci nous permet d’écrire:
En supposant que les modes propres dynamiques et les déformées statiques sont organisés comme le montre la formule [éq ] et en considérant que les vecteurs propres associés aux modes dynamiques sont normés par rapport à la masse modale unitaire, les matrices de masse et de rigidité généralisées prennent la forme suivante:
\({\stackrel{ˉ}{M}}^{k}=\left[\begin{array}{cc}\mathrm{Id}& {\varphi}^{{k}^{T}}{M}^{k}{\psi}^{k}\\ {\psi}^{{k}^{T}}{M}^{k}{\varphi}^{k}& {\psi}^{{k}^{T}}{M}^{k}{\psi}^{k}\end{array}\right]{\stackrel{ˉ}{K}}^{k}=\left[\begin{array}{cc}{\lambda}^{k}& {\varphi}^{{k}^{T}}{K}^{k}{\psi}^{k}\\ {\psi}^{{k}^{T}}{K}^{k}{\varphi}^{k}& {\psi}^{{k}^{T}}{K}^{k}{\psi}^{k}\end{array}\right]\)
On démontre, dans le cas de la méthode de Craig-Bampton, que les modes normaux et les modes contraints sont orthogonaux vis à vis de la matrice de rigidité dont les termes extra-diagonaux sont, dès lors, nuls. Cependant, cette propriété n’est pas utilisée dans l’algorithme programmé dans Code_Aster .
Nous ne considérons, comme type de dissipation, que l’amortissement visqueux (c’est le seul qui est supporté par les outils de sous-structuration dans Code_Aster ). Deux méthodes sont utilisables pour prendre en compte cet amortissement:
l’amortissement de Rayleigh appliqué au niveau élémentaire qui consiste à supposer que la matrice d’amortissement élémentaire \({C}_{e}\) associée à chaque élément fini du modèle est une combinaison linéaire des matrices de masse et de rigidité élémentaires \({K}_{e}\) et \({M}_{e}\) :
\({C}_{e}={\alpha}_{e}{K}_{e}+{\beta }_{e}{M}_{e}\)
La matrice d’amortissement est alors assemblée \({C}^{k}\) puis projetée sur la base [éq ]:
\({\stackrel{ˉ}{C}}^{k}=\left[\begin{array}{cc}{\varphi}^{{k}^{T}}{C}^{k}{\varphi}^{k}& {\varphi}^{{k}^{T}}{C}^{k}{\psi}^{k}\\ {\psi}^{{k}^{T}}{C}^{k}{\varphi}^{k}& {\psi}^{{k}^{T}}{C}^{k}{\psi}^{k}\end{array}\right]\)
l’amortissement proportionnel appliqué aux modes propres dynamiques de chaque sous‑structure. La matrice résultante est donc une diagonale incomplète (on ne sait pas associer d’amortissement proportionnel aux déformées statiques):
\({\stackrel{ˉ}{C}}^{k}=\left[\begin{array}{cc}{\xi}^{k}& 0\\ 0& 0\end{array}\right]\)
Équations dynamiques vérifiées par la structure globale#
Les équations dynamiques que vérifie la structure globale sont:
\(\begin{array}{}\left[\begin{array}{ccccc}{\stackrel{ˉ}{M}}^{1}& & & & \\ & \mathrm{...}& & & \\ & & {\stackrel{ˉ}{M}}^{k}& & \\ & & & \mathrm{...}& \\ & & & & {\stackrel{ˉ}{M}}^{{N}_{S}}\end{array}\right]\left\lbrace \begin{array}{c}{\ddot{\eta}}^{1}\\ \mathrm{...}\\ {\ddot{\eta}}^{k}\\ \mathrm{...}\\ {\ddot{\eta}}^{{N}_{S}}\end{array}\right\rbrace +\left[\begin{array}{ccccc}{\stackrel{ˉ}{C}}^{1}& & & & \\ & \mathrm{...}& & & \\ & & {\stackrel{ˉ}{C}}^{k}& & \\ & & & \mathrm{...}& \\ & & & & {\stackrel{ˉ}{C}}^{{N}_{S}}\end{array}\right]\left\lbrace \begin{array}{c}{\dot{\eta}}^{1}\\ \mathrm{...}\\ {\dot{\eta}}^{k}\\ \mathrm{...}\\ {\dot{\eta}}^{{N}_{S}}\end{array}\right\rbrace +\left[\begin{array}{ccccc}{\stackrel{ˉ}{K}}^{1}& & & & \\ & \mathrm{...}& & & \\ & & {\stackrel{ˉ}{K}}^{k}& & \\ & & & \mathrm{...}& \\ & & & & {\stackrel{ˉ}{K}}^{{N}_{S}}\end{array}\right]\left\lbrace \begin{array}{c}{\eta}^{1}\\ \mathrm{...}\\ {\eta}^{k}\\ \mathrm{...}\\ {\eta}^{{N}_{S}}\end{array}\right\rbrace \\ =\left\lbrace \begin{array}{c}{\stackrel{ˉ}{f}}_{\mathrm{ext}}^{1}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{\mathrm{ext}}^{k}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{\mathrm{ext}}^{{N}_{S}}\end{array}\right\rbrace +\left\lbrace \begin{array}{c}{\stackrel{ˉ}{f}}_{L}^{1}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{L}^{k}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{L}^{{N}_{S}}\end{array}\right\rbrace \end{array}\)
soit en harmonique:
\(\begin{array}{}(\left[\begin{array}{ccccc}{\stackrel{ˉ}{\stackrel{ˉ}{K}}}^{1}& & & & \\ & {}^{...}& & & \\ & & {\stackrel{ˉ}{\stackrel{ˉ}{K}}}^{k}& & \\ & & & {}^{...}& \\ & & & & {\stackrel{ˉ}{\stackrel{ˉ}{K}}}^{{N}_{s}}\end{array}\right]+j\omega \left[\begin{array}{ccccc}{\stackrel{ˉ}{\stackrel{ˉ}{C}}}^{1}& & & & \\ & {}^{...}& & & \\ & & {\stackrel{ˉ}{\stackrel{ˉ}{C}}}^{k}& & \\ & & & {}^{...}& \\ & & & & {\stackrel{ˉ}{\stackrel{ˉ}{C}}}^{{N}_{s}}\end{array}\right]-{\omega}^{2}\left[\begin{array}{ccccc}{\stackrel{ˉ}{\stackrel{ˉ}{M}}}^{1}& & & & \\ & {}^{...}& & & \\ & & {\stackrel{ˉ}{\stackrel{ˉ}{M}}}^{k}& & \\ & & & {}^{...}& \\ & & & & {\stackrel{ˉ}{\stackrel{ˉ}{M}}}^{{N}_{s}}\end{array}\right])\left\lbrace \begin{array}{c}{\eta}^{1}\\ {}^{...}\\ {\eta}^{k}\\ {}^{...}\\ {\eta}^{{N}_{s}}\end{array}\right\rbrace \\ =\left\lbrace \begin{array}{c}{\stackrel{ˉ}{f}}_{\mathrm{ext}}^{1}\\ {}^{...}\\ {\stackrel{ˉ}{f}}_{\mathrm{ext}}^{k}\\ {}^{...}\\ {\stackrel{ˉ}{f}}_{\mathrm{ext}}^{{N}_{s}}\end{array}\right\rbrace +\left\lbrace \begin{array}{c}{\stackrel{ˉ}{f}}_{L}^{1}\\ {}^{...}\\ {\stackrel{ˉ}{f}}_{L}^{k}\\ {}^{...}\\ {\stackrel{ˉ}{f}}_{L}^{{N}_{s}}\end{array}\right\rbrace \end{array}\)
Auxquelles il faut ajouter les équations de liaison (d’après [éq ]):
\(\forall k,l\text{}{L}_{{S}^{k}\cap {S}^{l}}^{k}\left\lbrace {\eta}^{k}\right\rbrace ={L}_{{S}^{k}\cap {S}^{l}}^{l}\left\lbrace {\eta}^{l}\right\rbrace\)
Ce système est résolu par double dualisation des conditions aux limites [R3.03.01]. Sa formulation finale fait donc intervenir le vecteur des multiplicateurs de Lagrange \(\lambda\) et peut s’écrire sous la forme condensée:
Le problème défini par l’équation [éq ] est symétrique. D’autre part, sa dimension est déterminée par le nombre de modes pris en considération (modes dynamiques et déformées statiques). On est donc amené à résoudre un problème harmonique classique, de taille réduite, auquel est associée une équation de contrainte linéaire. Sa résolution ne pose donc pas de problème.
Mise en œuvre dans Code_Aster#
Étude des sous-structures séparément#
Les paramètres \({\alpha}_{e}\) et \({\beta }_{e}\) de l’amortissement de Rayleigh sont introduits, le cas échéant, par l’opérateur DEFI_MATERIAU [U4.43.01].
Les traitements des sous-structures sont identiques au cas du calcul modal.
L’opérateur MACR_ELEM_DYNA [U4.65.01] calcule les matrices généralisées de rigidité, de masse et éventuellement d’amortissement de la sous-structure, ainsi que les matrices de liaison. L’amortissement de Rayleigh est pris en compte en complétant l’opérande MATR_AMOR. L’amortissement proportionnel est introduit par l’opérande AMOR_REDUIT.
Le chargement harmonique est défini, au niveau de la sous-structure, par les opérateurs AFFE_CHAR_MECA [U4.44.01] (application de la force sur le maillage), CALC_VECT_ELEM [U4.61.02] (calcul des vecteurs élémentaires associés) et ASSE_VECTEUR [U4.61.23] (assemblage du vecteur de chargement sur le maillage de la sous-structure).
Assemblage et résolution#
Comme dans le cas du calcul modal, le modèle de la structure complète est défini par l’opérateur DEFI_MODELE_GENE [U4.65.02]. Sa numérotation est réalisée par l’opérateur NUME_DDL_GENE [U4.65.03]. Les matrices de masse, de rigidité et éventuellement d’amortissement généralisées de la structure complète sont assemblées en fonction de cette numérotation avec l’opérateur ASSE_MATR_GENE [U4.65.04].
Les chargements sont projetés sur les bases des sous-structures auxquelles ils sont appliqués, puis assemblés à partir de la numérotation issue de NUME_DDL_GENE [U4.65.03] par l’opérateur ASSE_VECT_GENE [U4.65.05].
Le calcul de la réponse harmonique de la structure complète est réalisé par l’opérateur DYNA_LINE_HARM [U4.53.11].
Restitution sur base physique#
La restitution des résultats sur base physique est identique au cas du calcul modal. Elle fait intervenir l’opérateur REST_GENE_PHYS [U4.63.31] et éventuellement l’opérateur DEFI_SQUELETTE [U4.24.01] (création d’un maillage « squelette »).
Conclusion#
La méthode de calcul de réponse harmonique par synthèse modale disponible dans Code_Aster s’appuie sur celle de sous-structuration modale, également programmée. Elle consiste à exprimer l’ensemble des équations dans un espace de dimension réduite, constitué de modes des différentes sous-structures, par une méthode de Rayleigh-Ritz. La définition de ces champs est celle utilisée pour la sous-structuration modale et comprend des modes normaux ainsi que d’autres statiques ou harmoniques. La procédure employée se traduit par une projection des matrices et du second membre sur l’espace restreint.
Réponse transitoire par sous structuration dynamique classique#
Introduction#
L’objet de cette documentation de référence est de présenter les bases théoriques des deux méthodes de calcul de réponse transitoire par sous-structuration dynamique disponibles dans Code_Aster . La première consiste à réaliser un calcul transitoire par sous-structuration pour lequel les équations du problème sont projetées sur les bases associées à chaque sous-structure. La difficulté réside dans la double dualisation des conditions aux limites qui aboutit à une matrice de masse singulière. Pour utiliser les schémas d’intégration explicites (qui nécessitent l’inversion de la matrice de masse), il faut donc modifier le traitement des interfaces dans l’opérateur de calcul de la réponse transitoire. La seconde méthode consiste à déterminer les modes propres de la structure complète par sous‑structuration et à projeter sur cette base les équations du problème transitoire. L’étape de restitution sur base physique des résultats généralisés finaux doit donc tenir compte de cette double projection.
L’opérateur de calcul de réponse transitoire qui reçoit la sous-structuration est l’opérateur DYNA_TRAN_MODAL [U4.53.21]. S’appuyant sur des méthodes de recombinaison modale, il a été conçu pour résoudre des problèmes transitoires en coordonnées généralisées et il est très efficace pour les problèmes de grande taille dont il permet de réduire le nombre de degrés de liberté. D’autre part, il supporte la prise en compte de non-linéarités localisées (aux nœuds) que l’on souhaite généraliser au cas de la sous-structuration.
Dans ce rapport, nous présentons les deux méthodes de calcul transitoire par sous‑structuration disponible dans Code_Aster , ainsi que leur mise en œuvre.
Calcul transitoire par projection sur les bases des sous-structures#
Équations dynamiques vérifiées par les sous-structures séparément#
Soit une structure \(S\) composée de \({N}_{S}\) sous-structures notées \({S}^{k}\) . Nous supposons que chaque sous-structure est modélisée en éléments finis. Le comportement vibratoire des sous-structures résulte des forces extérieures qui lui sont appliquées et des forces de liaison qu’exercent sur elles les autres sous-structures. Ainsi, pour \({S}^{k}\) , nous avons:
Le champ de déplacement inconnu, issu de la modélisation éléments finis, est recherché sur un espace approprié, de dimension réduite (transformation de Ritz) selon l’expression:
La transformation de Ritz [éq ], appliquée au déplacement réel et virtuel conduit à l’équation dynamique transitoire de la sous-structure [éq ], permet d’écrire:
Le problème défini par l’équation [éq ] est symétrique. D’autre part, sa dimension est déterminée par le nombre de modes pris en considération (modes dynamiques et déformées statiques). On est donc amené à résoudre un problème transitoire classique mais de taille réduite.
Équations dynamiques vérifiées par la structure globale#
L’écriture matricielle de l’équation dynamique vérifiée par la structure globale, s’écrit simplement à partir des équations dynamiques vérifiées par chaque sous-structure:
\(\begin{array}{}\left[\begin{array}{ccccc}{\stackrel{ˉ}{M}}^{1}& & & & \\ & \mathrm{...}& & & \\ & & {\stackrel{ˉ}{M}}^{k}& & \\ & & & \mathrm{...}& \\ & & & & {\stackrel{ˉ}{M}}^{{N}_{S}}\end{array}\right]\left\lbrace \begin{array}{c}{\ddot{\eta}}^{1}\\ \mathrm{...}\\ {\ddot{\eta}}^{k}\\ \mathrm{...}\\ {\ddot{\eta}}^{{N}_{S}}\end{array}\right\rbrace +\left[\begin{array}{ccccc}{\stackrel{ˉ}{C}}^{1}& & & & \\ & \mathrm{...}& & & \\ & & {\stackrel{ˉ}{C}}^{k}& & \\ & & & \mathrm{...}& \\ & & & & {\stackrel{ˉ}{C}}^{{N}_{S}}\end{array}\right]\left\lbrace \begin{array}{c}{\dot{\eta}}^{1}\\ \mathrm{...}\\ {\dot{\eta}}^{k}\\ \mathrm{...}\\ {\dot{\eta}}^{{N}_{S}}\end{array}\right\rbrace +\left[\begin{array}{ccccc}{\stackrel{ˉ}{K}}^{1}& & & & \\ & \mathrm{...}& & & \\ & & {\stackrel{ˉ}{K}}^{k}& & \\ & & & \mathrm{...}& \\ & & & & {\stackrel{ˉ}{K}}^{{N}_{S}}\end{array}\right]\left\lbrace \begin{array}{c}{\eta}^{1}\\ \mathrm{...}\\ {\eta}^{k}\\ \mathrm{...}\\ {\eta}^{{N}_{S}}\end{array}\right\rbrace \\ =\left\lbrace \begin{array}{c}{\stackrel{ˉ}{f}}_{\mathrm{ext}}^{1}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{\mathrm{ext}}^{k}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{\mathrm{ext}}^{{N}_{S}}\end{array}\right\rbrace +\left\lbrace \begin{array}{c}{\stackrel{ˉ}{f}}_{L}^{1}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{L}^{k}\\ \mathrm{...}\\ {\stackrel{ˉ}{f}}_{L}^{{N}_{S}}\end{array}\right\rbrace \end{array}\)
Auxquelles, il faut ajouter les équations de liaison:
\(\forall k,l\text{}{\mathrm{L}}_{{S}^{k}\cap {S}^{l}}^{k}\left\lbrace {\eta}^{k}\right\rbrace ={\mathrm{L}}_{{S}^{k}\cap {S}^{l}}^{l}\left\lbrace {\eta}^{l}\right\rbrace \text{}{\mathrm{f}}_{{L}_{{S}^{k}\cap {S}^{L}}}^{k}=-{\mathrm{f}}_{{L}_{{S}^{k}\cap {S}^{L}}}^{l}\)
Ce système peut s’écrire sous la forme condensée:
Double dualisation des conditions aux limites#
Le problème condensé, donné ci-dessus, se présente sous la forme d’un système transitoire auquel est associée une équation de contrainte linéaire (en force et déplacement). Dans Code_Aster , ce type de problème est classique et il est résolu par double dualisation des conditions aux limites [R3.03.01], c’est à dire par l’introduction de variables auxiliaires encore appelées multiplicateurs de Lagrange pour dualiser les conditions aux limites. Après introduction des multiplicateurs de Lagrange, le système matriciel se met sous la forme:
où \({\lambda}_{1}\) et \({\lambda}_{2}\) sont les multiplicateurs de Lagrange.
On le constate, l’introduction des multiplicateurs de Lagrange rend singulière la matrice de masse. Dès lors, l’utilisation des schéma d’intégration explicites développés dans l’opérateur DYNA_TRAN_MODAL [U4.53.21] de Code_Aster est impossible car ils nécessitent l’inversion de la matrice de masse.
Pour rendre la matrice non singulière, il suffit de dualiser avec les mêmes multiplicateurs de Lagrange, la condition sur la dérivée seconde des équations de liaison.
Ainsi la condition de continuité des déplacements [éq ] est modifiée par le système équivalent:
où \(\eta °\) et \(\dot{\eta}°\) sont les déplacements et vitesses généralisés à l’instant initial.
Le système matriciel qui en découle est de la forme [bib10]:
On constate que la matrice de masse a la même forme que la matrice de raideur. Elles sont donc inversibles. Ce système est donc parfaitement équivalent à l’équation [éq ] (il vérifie à tout instant les conditions de liaison) et il peut être traité, sous cette forme, par l’opérateur DYNA_TRAN_MODAL.
Traitement de la matrice d’amortissement#
On constate que la condition de continuité des déplacements, formulée dans l’équation [éq], se traduit par une équation du second ordre non amortie. Lors de la résolution par pas de temps d’un problème transitoire, toute erreur numérique risque de s’auto-entretenir, diminuant ainsi la stabilité de l’algorithme. Pour optimiser l’amortissement de l’erreur numérique, il suffit de dualiser la condition sur la dérivée première des équations de liaison avec les mêmes multiplicateurs de Lagrange multipliés par 2 (de manière à rendre cet amortissement critique):
Le système matriciel qui en découle est de la forme [bib10]:
\(\left[\begin{array}{ccc}-\text{Id}& \mathrm{L}& \text{Id}\\ {\mathrm{L}}^{\mathrm{T}}& \stackrel{ˉ}{\mathrm{M}}& {\mathrm{L}}^{\mathrm{T}}\\ \text{Id}& \mathrm{L}& -\text{Id}\end{array}\right]\left\lbrace \begin{array}{c}\ddot{{\lambda}_{1}}\\ \ddot{\eta}\\ \ddot{{\lambda}_{2}}\end{array}\right\rbrace +\left[\begin{array}{ccc}-2\text{Id}& 2\mathrm{L}& 2\text{Id}\\ 2{\mathrm{L}}^{\mathrm{T}}& \stackrel{ˉ}{\mathrm{C}}& 2{\mathrm{L}}^{\mathrm{T}}\\ 2\text{Id}& 2\mathrm{L}& -2\text{Id}\end{array}\right]\left\lbrace \begin{array}{c}\dot{{\lambda}_{1}}\\ \dot{\eta}\\ \dot{{\lambda}_{2}}\end{array}\right\rbrace +\left[\begin{array}{ccc}-\text{Id}& \mathrm{L}& \text{Id}\\ {\mathrm{L}}^{\mathrm{T}}& \stackrel{ˉ}{\mathrm{K}}& {\mathrm{L}}^{\mathrm{T}}\\ \text{Id}& \mathrm{L}& -\text{Id}\end{array}\right]\left\lbrace \begin{array}{c}{\lambda}_{1}\\ \eta \\ {\lambda}_{2}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}0\\ {\stackrel{ˉ}{\mathrm{f}}}_{\mathit{ext}}\\ 0\end{array}\right\rbrace\)
On constate donc que le traitement de l’erreur numérique sur les équations de liaison conduit à modifier la matrice d’amortissement du problème transitoire. Cette modification est tout à fait comparable à celle qui est réalisée sur la matrice de masse.
Par contre, nous n’avons pas souhaité généraliser ce traitement au cas de la résolution des problèmes transitoires non amortis. Cela nous aurait conduit à créer une matrice d’amortissement temporaire. On aurait pu craindre d’augmenter les temps de calcul, sans réel bénéfice. D’ailleurs, il est tout à fait possible à l’utilisateur de définir une matrice d’amortissement dont les coefficients sont nuls. La modification de celle-ci étant automatique, le système transitoire non amorti sera effectivement résolu, tout en optimisant le traitement de toute erreur numérique intervenant sur les équations de liaison.
Traitement des conditions initiales#
Considérons une sous-structure \({S}^{k}\) caractérisée par sa base de projection \({\Phi}^{k}\) composée de modes normaux et de déformées statiques. On suppose qu’initialement la sous-structure \({S}^{k}\) est soumise à un champ de déplacement ou de vitesse (cela ne modifie en rien la démonstration) noté :math:`` . La transformation de Ritz nous permet d’écrire:
\({\mathrm{q}}_{0}^{k}={\Phi}^{k}{\eta}_{0}^{k}\)
où \({\eta}_{0}^{k}\) est le vecteur des déplacements (ou vitesses) généralisé(e)s de \({S}^{k}\) à déterminer.
Le vecteur des déplacements (ou vitesses) généralisé(e)s initiaux(ales) est déterminé comme suit:
\({q}_{0}^{k}={\Phi}^{k}{\eta}_{0}^{k}\Rightarrow {\Phi}^{{k}^{T}}{q}_{0}^{k}=({\Phi}^{{k}^{T}}{\Phi}^{k})\cdot {\eta}_{0}^{k}\)
\(\mathrm{\Rightarrow }{\eta}_{0}^{k}={({\Phi}^{{k}^{T}}{\Phi}^{k})}^{-1}{\Phi}^{{k}^{T}}{\mathrm{q}}_{0}^{k}\)
Calcul transitoire sur une base modale globale calculée par sous-structuration#
La seconde méthode développée consiste à résoudre le problème transitoire sur la base modale de la structure complète calculée par sous-structuration.
Calcul des modes propres de la structure complète par sous‑structuration#
Chaque sous-structure \({S}^{k}\) est représentée par une base de projection, composée de modes propres dynamiques et de déformées statiques, que nous avons notée \({\Phi}^{k}\) . La base de projection de la structure complète qui en résulte est notée \(\Phi\) .
La base modale de la structure complète est calculée par sous-structuration. Chaque mode obtenu est donc combinaison linéaire des vecteurs des bases de projection des sous-structures :
où:
\({\Phi}_{p}\) |
est la matrice des modes propres de la structure complète, |
\(\alpha\) |
est la matrice des coordonnées modales généralisées de la structure. |
La projection des matrices et des vecteurs constitutifs du problème transitoire, sur la base des modes propres de la structure complète calculés par synthèse modale permet de déterminer:
|
\(\stackrel{ˉ}{\stackrel{ˉ}{M}}={\alpha}^{T}\stackrel{ˉ}{M}\alpha\) |
|
\(\stackrel{ˉ}{\stackrel{ˉ}{K}}={\alpha}^{T}\stackrel{ˉ}{C}\alpha\) |
|
\(\stackrel{ˉ}{\stackrel{ˉ}{C}}={\alpha}^{T}\stackrel{ˉ}{C}\alpha\) |
|
\({\stackrel{ˉ}{\stackrel{ˉ}{f}}}_{\mathrm{ext}}={\alpha}^{T}{\stackrel{ˉ}{f}}_{\mathrm{ext}}\) |
Du fait de l’orthogonalité des modes propres de la structure calculés par synthèse modale, par rapport aux matrices \(M\) et \(K\) , les matrices de masse et de rigidité généralisées obtenues ci‑dessus sont diagonales:
Équation dynamique vérifiée par la structure globale#
La structure complète est soumise aux forces extérieures qui lui sont appliquées. Ainsi, nous pouvons écrire:
Le champ de déplacement inconnu, issu de la modélisation éléments finis, est remplacé par sa projection sur la base des modes propres de la structure, selon la formule:
où \(\eta\) est le vecteur des coordonnées généralisées de la structure.
La transformation de Ritz [éq ], appliquée à l’équation dynamique transitoire de la structure [éq ], permet d’écrire:
L’étape de restitution sur base physique nécessite de tenir compte de la double projection (sur la base modale de la structure complète, puis sur les bases de projection des sous‑structures - cf. éq ).
Le problème défini par l’équation [éq ] est de forme tout à fait classique. On est amené à résoudre un problème transitoire symétrique dont la dimension est déterminée par le nombre de modes calculés par sous-structuration et dont les matrices de masse et de rigidité sont diagonales.
Notons enfin que le traitement des conditions initiales est identique au cas du calcul transitoire par projection sur les bases des sous-structures (cf. § 6.2 ).
Étude comparative des deux méthodes développées#
Les bases théoriques, associées aux deux méthodologies mises en œuvre dans Code_Aster pour réaliser un calcul de réponse transitoire en utilisant les techniques de sous-structuration, ont été présentées dans les chapitres précédents. Nous précisons, ici, leurs particularités essentielles.
La première méthodologie consiste à faire un calcul de réponse transitoire par sous‑structuration. L’équation vérifiée par la structure complète est alors projetée sur les bases des sous-structures. La précision de cette méthode est donc directement déterminée par l’étendue de ces bases. Ces dernières peuvent être enrichies sans aboutir à des temps de calcul prohibitifs car les sous-structures sont, en principe, de tailles relativement réduites. Quoi qu’il en soit, il est difficile d’estimer l’effet de troncature modale à la seule connaissance des modes des bases des sous-structures. D’autre part, les bases de projection des sous-structures sont composées de modes qui ne sont pas tous orthogonaux entre eux (modes propres et déformées statiques). Les matrices de masse et de rigidité généralisées constitutives du problème transitoire final sont donc non-diagonales. Globalement, leur largeur de bande peut être déterminée à partir du nombre de déformées statiques des bases de projection des sous-structures. La durée d’intégration dans l’opérateur de calcul transitoire DYNA_TRAN_MODAL sera donc d’autant plus longue qu’il y aura de degrés de liberté d’interface. D’autre part, le pas de temps d’intégration maximal admissible par les schéma d’intégration explicites est déterminé à partir de la fréquence maximale de la base de projection. Dans le cas d’un calcul transitoire par sous-structuration, cette fréquence résulte, en principe, des modes statiques dont les termes diagonaux sont élevés dans la matrice de rigidité généralisée et faibles dans la matrice de masse généralisée. Dès lors, le pas de temps d’intégration ne peut pas être déterminé a priori. L’expérience montre qu’il est très faible, au regard des fréquences propres des bases des sous-structures et que l’utilisation du schéma d’intégration à pas de temps adaptatif de DYNA_TRAN_MODAL est très profitable.
La deuxième méthodologie consiste à faire un calcul transitoire sur la base modale de la structure complète obtenue par sous-structuration. On sait que l’étape consistant à calculer les modes propres de la structure peut être coûteuse en terme de temps de calcul. Ceci est d’autant plus vrai lorsqu’on considère des forces non linéaires car la base de projection doit alors être suffisamment étendue pour bien représenter la dynamique du système. D’autre part, la base modale sur laquelle est calculée la réponse transitoire est généralement de dimension inférieure à celle déterminée par les vecteurs propres des sous-structures (modes propres et déformées statiques). La double projection revient donc à introduire une fréquence de coupure. On doit donc s’attendre à ce que cette méthode soit moins précise que la précédente. Cependant, le calcul des modes propres permet d’estimer l’effet de troncature modale. D’autre part, il peut permettre de valider les modèles des sous-stuctures si on dispose de résultats expérimentaux. Finalement, l’intérêt essentiel de cette méthode est que les matrices de masse et de rigidité utilisées dans le calcul transitoire sont diagonales. L’intégration numérique est donc très rapide.
Mise en œuvre dans Code_Aster#
Étude des sous-structures séparément#
Si on souhaite introduire un amortissement de Rayleigh, les paramètres
et
de cet amortissement sont définis, par l’opérateur DEFI_MATERIAU [U4.43.01].
Les traitements des sous-structures sont identiques au cas du calcul modal et harmonique. Les modes propres dynamiques sont calculés avec l’opérateur CALC_MODES [U4.52.02]. Les conditions aux interfaces de liaison sont appliquées avec l’opérateur AFFE_CHAR_MECA [U4.44.01].
L’opérateur DEFI_INTERF_DYNA [U4.64.01] permet de définir les interfaces de connexion de la sous‑structure. L’opérateur DEFI_BASE_MODALE [U4.64.02] permet de calculer la base de projection complète de la sous-structure (recopie des modes propres et calcul des déformées statiques).
L’opérateur MACR_ELEM_DYNA [U4.65.01] calcule les matrices généralisées de raideur, de masse et éventuellement d’amortissement de la sous-structure, ainsi que les matrices de liaison. L’amortissement de Rayleigh est pris en compte en complétant l’opérande MATR_AMOR. L’amortissement proportionnel est introduit par l’opérande AMOR_REDUIT.
Le chargement transitoire est défini, au niveau de la sous-structure, par les opérateurs AFFE_CHAR_MECA [U4.44.01] (application de la force sur le maillage), CALC_VECT_ELEM [U4.61.02] (calcul des vecteurs élémentaires associés) et ASSE_VECTEUR [U4.61.23] (assemblage du vecteur de chargement sur le maillage de la sous-structure).
L’opérateur CREA_CHAMP [U4.72.04] qui permet d’affecter un champ sur les nœuds d’un modèle permet de décrire le champ de déplacement initial ou/et le champ de vitesse initial de la sous-structure.
Assemblage du modèle généralisé#
Comme dans le cas du calcul modal et harmonique, le modèle de la structure complète est défini par l’opérateur DEFI_MODELE_GENE [U4.65.02]. Sa numérotation est réalisée par l’opérateur NUME_DDL_GENE [U4.65.03]. Les matrices de masse, de raideur et éventuellement d’amortissement généralisées de la structure complète sont assemblées en fonction de cette numérotation avec l’opérateur ASSE_MATR_GENE [U4.65.04].
Les chargements sont projetés sur les bases des sous-structures auxquelles ils sont appliqués, puis assemblés à partir de la numérotation issue de NUME_DDL_GENE [U4.65.03] par l’opérateur ASSE_VECT_GENE [U4.65.05].
Les déplacements généralisés initiaux et les vitesses généralisées initiales pour chaque sous-structure, sont calculés par l’opérateur ASSE_VECT_GENE [U4.65.05]. Cet opérateur réalise également l’assemblage de ces vecteurs en fonction de la numérotation issue de NUME_DDL_GENE [U4.65.03].
Dans le cas d’un calcul transitoire projeté sur les « bases » des sous-structures, les matrices et vecteurs assemblés généralisés obtenus à l’issu de cette étape sont directement utilisés pour le calcul transitoire. Dans le cas d’un calcul sur la base modale de la structure complète calculée par sous‑structuration, il faut réaliser des opérations spécifiques qui sont présentées au § 5.3 .
Calcul de la base modale de la structure complète et projection#
Ce chapitre est spécifique au calcul transitoire sur base modale calculée par sous‑structuration.
La base modale de la structure complète est calculée avec l’opérateur classique de Code_Aster: CALC_MODES [U4.52.02]. On définit une numérotation du problème généralisé final avec l’opérateur NUME_DDL_GENE [U4.65.03]. Les matrices de masse, de raideur et éventuellement d’amortissement généralisées sont projetées sur la base des modes propres de la structure avec l’opérateur PROJ_MATR_BASE [U4.63.12]. Les vecteurs généralisés correspondant aux chargements extérieurs sont projetés sur la base des modes propres de la structure avec l’opérateur PROJ_VECT_BASE [U4.63.13].
Résolution et restitution sur base physique#
Le calcul de la réponse transitoire de la structure complète est réalisé par l’opérateur DYNA_TRAN_MODAL [U4.53.21].
La restitution des résultats sur base physique fait intervenir l’opérateur REST_GENE_PHYS [U4.63.31] ; elle est identique au cas du calcul modal et harmonique. On peut utiliser l’opérateur DEFI_SQUELETTE [U4.24.01] pour créer un maillage « squelette ». Plus grossier que le maillage de calcul, il permet de réduire les durées des traitements graphiques.
Conclusion#
Nous avons présenté, dans ce rapport, les travaux réalisés pour introduire, dans Code_Aster , le calcul de réponse transitoire linéaire par sous-structuration dynamique. Les méthodes qui ont été choisies consistent, pour la première d’entre elles, à projeter les équations transitoires sur les «bases» de chaque sous-structure, composées de modes propres dynamiques et de déformées statiques et pour la seconde, à calculer les modes propres de la structure complète par sous-structuration et à y projeter les équations transitoires.
Nous avons débuté par un exposé des bases théoriques sur lesquelles s’appuient la première méthode de sous-structuration transitoire pour aboutir à la formulation matricielle du problème final. En particulier, un traitement original de l’équation de continuité des déplacements pour rendre la matrice de masse inversible et pour assurer une stabilité optimale de l’algorithme d’intégration, nous a conduits à modifier la forme des matrices de masse et d’amortissement du problème transitoire.
Pour la seconde méthode, la difficulté essentielle consiste à restituer les résultats obtenus en coordonnées généralisées sur la base physique. En effet, il faut tenir compte de la double projection : sur la base modale de la structure complète d’une part, et sur les bases des sous‑structures d’autre part.
Les développements réalisés se sont traduits par des modifications des opérateurs DYNA_TRAN_MODAL [U4.53.21] et REST_GENE_PHYS [U4.63.31]. Leur syntaxe a été très peu modifiée, de manière à ce que leur utilisation soit identique lors d’un calcul par sous‑structuration et d’un calcul direct par recombinaison modale.
Condensation dynamique par macro-éléments dynamiques en sous-structuration statique#
Introduction#
L’objectif de la condensation dynamique de modèle par sous-structuration statique est de réduire le nombre total de degrés de liberté d’un modèle constitué de plusieurs sous-domaines en condensant la résolution des sous-domaines de comportement linéaire sur les degrés de liberté de leurs interfaces avec des sous-domaines de comportement non linéaire.
L’idée est alors de représenter chacun de ces sous-domaines par un macro-élément de type statique pour pouvoir l’utiliser dans un modèle mixte comprenant également des éléments finis pour les parties non condensées de comportement non linéaire afin de procéder à des analyses non linéaires au moyen des opérateurs de calcul non linéaires STAT_NON_LINE [U4.51.03] ou DYNA_NON_LINE [U4.53.01].
Dans la condensation de chaque sous-domaine, on veut également intégrer la partie dynamique de son comportement obtenue par analyse modale et pour cela, on va donc utiliser des macro-éléments dynamiques comme les macro-éléments statiques nécessaires en sous-structuration statique.
Le principe de cette conversion et leur application est décrit dans les documents [bib13].
Problème classique et principe de la méthode#
Figure 7.2-a: Sous-domaines du problème classique de condensation.
Sur la figure 7.2-a ci-dessus sont schématisés les sous-domaines éventuels que l’on peut rencontrer dans le cadre d’un problème classique de condensation par sous-structuration statique:
les sous-domaines de type \(\mathit{E1}\) : sous-domaines de structure de comportement linéaire constant condensés et modélisés par macro-élément dynamique,
les sous-domaines de type \(\mathit{E2}\) : sous-domaines de structure de comportement potentiellement non linéaire non condensés et modélisés par des éléments finis,
les sous-domaines de type \(\mathit{E3}\) : sous-domaine de sol représentés par un macro-élément issu de la conversion d’une matrice complexe impédance de sol.
Le principe de base de la condensation par macro-élément consiste à partitionner pour chaque sous-domaine de type \(\mathit{E1}\) les degrés de liberté \(j\) de leurs interfaces notées \(\Sigma\) avec les sous-domaines de type \(\mathit{E2}\) ainsi que leurs autres degrés de liberté \(i\) dits «internes».
Les matrices de rigidité \(K\) et de masse \(M\) se décomposent ainsi sur ces degrés de liberté:
\(K=\left[\begin{array}{cc}{K}_{ii}& {K}_{ij}\\ {K}_{ji}& {K}_{jj}\end{array}\right]\) \(M=\left[\begin{array}{cc}{M}_{ii}& {M}_{ij}\\ {M}_{ji}& {M}_{jj}\end{array}\right]\) \(\begin{array}{}i\text{ddl internes}\\ j\text{ddl d'interfaces}\end{array}\)
Le principe de sous-structuration statique revient pour la matrice \(K\) par exemple à la condenser uniquement sur les degrés de liberté \(j\) des interfaces \(\Sigma\) en corrigeant le terme initial \({K}_{jj}\) d’un terme complémentaire dit «complément de Schur» selon le principe de la condensation dite de Guyan [bib14]:
Le calcul du nouveau terme \(\stackrel{ˉ}{{K}_{jj}}\) revient à une projection de la matrice \(K\) sur une base de modes statiques \(\Psi\) de type «contraint» obtenus par réponse statique à des déplacements unitaires imposés en chaque degré de liberté \(j\) des interfaces \(\Sigma\) :
\(\Psi =\left[\begin{array}{}{\Psi}_{ij}\\ I\end{array}\right]=\left[\begin{array}{}-{K}_{ii}^{-1}{K}_{ij}\\ I\end{array}\right]\)
La condensation sur les degrés de liberté \(j\) des interfaces \(\Sigma\) de la matrice de masse \(M\) et des charges des sous-domaines de type \(\mathit{E1}\) se fera de manière analogue par le produit matriciel de la matrice et des vecteurs de chargement sur la base de modes statiques \(\Psi\) .
On constate alors que les termes \({\Phi}^{T}K\Phi\) et \({\Psi}^{\mathrm{T.}}M.\Psi\) peuvent être directement obtenus pour chaque sous-domaine de type \(\mathit{E1}\) au moyen de l’opérateur de calcul d’un macro-élément dynamique MACR_ELEM_DYNA [U4.65.01] et que les termes de charges projetés sont obtenus \({\Psi}^{T}.F\) au moyen de l’opérateur PROJ_VECT_BASE [U4.63.13].
Méthodes de condensation par base complète ou réduite#
Selon le mode de calcul des modes de représentation des mouvements des interfaces, deux méthodes sont actuellement utilisées:
La méthode de calcul par base complète#
Le mouvement \(\mathrm{U}\) d’un point quelconque d’un sous-domaine de type E1 s’exprime à partir de la décomposition sur les modes statiques \(\Psi\) et sur un complément apporté par la partie dynamique du mouvement exprimée par des modes propres \(\Phi\) :
Avec la méthode par base complète, sur l’interface \(\Sigma\) , les modes propres \(\Phi\) ont une contribution nulle et les modes statiques \(\Psi\) ont une valeur unitaire et donc le déplacement \(U\) d’un point de l’interface s’exprime:
Dans ce cas, les degrés de liberté physiques \({X}_{\Sigma}\) de l’interface \(\Sigma\) se confondent avec les coordonnées généralisées \({q}_{\Sigma}\) associées aux modes statiques \(\Psi\) .
Les coordonnées généralisées \(q\) associées aux modes propres \(\Phi\) étant nulles sur l’interface, leur résolution sera découplée de celle des coordonnées généralisées \({q}_{\Sigma}\) . Dans l’assemblage des matrices du problème condensé, les termes de projection sur les modes propres \({\Phi}^{T}.K.\Phi\) et \({\Phi}^{T}.M.\Phi\) seront juxtaposés. Pour pouvoir «mélanger» des degrés de liberté physiques avec des coordonnées généralisées, on utilise l’artifice transparent pour l’utilisateur de faire passer ces dernières pour des degrés de liberté physiques supplémentaires pris parmi les degrés de liberté internes à un domaine de type \(\mathit{E1}\) .
Cet alias est actif dans le modèle où les domaines de type \(\mathit{E1}\) sont condensés par des macro-éléments, mais pour avoir le véritable déplacement des degrés de liberté internes à ces domaines, il s’obtient en post-traitement après résolution du problème global assemblé au moyen de l’opérateur REST_COND_TRAN [U4.63.33] en appliquant la relation initiale:
Au niveau de l’assemblage d’un sous-domaine de type \(\mathit{E1}\) , le stockage initial représenté sur la figure 7.3.1-a comprend des termes à la fois sur des degrés de liberté internes et sur les degrés de liberté d’interface avec un profil de remplissage non complet.
Figure 7.3.1-a: Assemblage initial des sous-domaines \(\mathit{E1}\) avant condensation.
Après condensation par la méthode complète, la plupart des degrés de liberté internes disparaissent, mais une toute petite partie est consacrée pour représenter les coordonnées généralisées associées aux modes propres et les termes du macro-élément statique s’ajoutent aux degrés de liberté d’interface avec un profil de remplissage complet cette fois-ci. Le stockage résultant est représenté sur la figure 7.3.1-b ci-dessous.
Figure 7.3.1-b: Assemblage des sous-domaines \(\mathit{E1}\) avec méthode complète.
L’assemblage des termes du macro-élément statique s’effectue par simple sommation sur les degrés de liberté d’interface et ne nécessite pas en général de liaison supplémentaire sur l’interface quand on prend en compte tous les modes statiques contraints de l’interface. Il n’en est pas ainsi dans le cas particulier de la condensation sur un seul nœud d’interface où là il faut ajouter une relation de liaison, solide sur les degrés de liberté d’interface.
Dans le cas d’un sous-domaines de type \(\mathit{E3}\) , le post-traitement par REST_COND_TRAN est inutile car le déplacement n’a pas de composante sur des modes dynamiques et le déplacement à l’interface est définitivement obtenu dès la résolution du problème global.
La méthode de calcul par réduction modale#
Avec la méthode par réduction modale, le mouvement \(U\) d’un point quelconque d’un sous-domaine de type \(\mathit{E1}\) s’exprime toujours à partir de la décomposition sur les modes statiques \(\Psi\) et sur un complément apporté par la partie dynamique du mouvement exprimée par des modes propres \(\Phi\) :
Mais cette fois-ci, sur l’interface \(\Sigma\) , les modes propres \(\Phi\) ont une contribution non nulle car les modes statiques \(\Psi\) ont une valeur non unitaire et sont en nombre généralement très inférieur au nombre de degrés de liberté physiques \({X}_{\Sigma}\) de l’interface \(\Sigma\) qui, dans ce cas, ne se confondent pas avec les coordonnées généralisées \({q}_{\Sigma}\) associées aux modes statiques \(\Psi\) . Il faut donc distinguer dans l’assemblage et la résolution les termes associés aux degrés de liberté physiques et ceux associés aux coordonnées généralisées. De la même façon que pour les coordonnées généralisées \(q\) associées aux modes propres \(\Phi\) , on utilise l’artifice transparent pour l’utilisateur de faire passer les coordonnées généralisées \({q}_{\Sigma}\) associées aux modes statiques \(\Psi\) pour des degrés de liberté physiques supplémentaires pris aussi parmi des degrés de liberté internes à un domaine de type \(\mathit{E1}\) .
Mais cette fois-ci, la résolution des degrés de liberté physiques \({X}_{\Sigma}\) de l’interface \(\Sigma\) ne sera pas découplée de celle des coordonnées généralisées \({\mathrm{q}}_{\Sigma}\) . Dans l’assemblage des matrices du problème condensé, les termes de projection sur les modes propres \({\Phi}^{T}.K.\Phi\) et \({\Phi}^{T}.M.\Phi\) seront juxtaposés. Mais il faudra alors ajouter une liaison et donc des degrés de liberté de Lagrange entre les degrés de liberté physiques \({X}_{\Sigma}\) de l’interface \(\Sigma\) et les coordonnées généralisées \({q}_{\Sigma}\) . Cette liaison est exprimée par la relation:
Cette relation s’exprima par le mot clé répétable LIAISON_INTERF del’opérateur AFFE_CHAR_MECA [U4.44.01], qui permet de relier chaque sous-domaine de type \(\mathit{E1}\) à un sous domaine de type \(\mathit{E2}\) (cf. figure 7.2-a) en générant un cas particulier de LIAISON_DDLavec les coefficients de la relation linéaire précédente.
Ensuite, pour avoir le véritable déplacement des degrés de liberté internes aux domaines de type \(\mathit{E1}\) , on procède de la même façon en post-traitement après résolution du problème global assemblé au moyen de l’opérateur REST_COND_TRAN [U4.63.33] en appliquant la relation initiale:
Après condensation par la méthode de réduction modale, la plupart des degrés de liberté internes disparaissent, mais une toute petite partie est consacrée pour représenter les coordonnées généralisées associées aux modes propres, une autre petite partie est consacrée aux coordonnées généralisées associées aux modes propres d’interface. Cette fois-ci, les termes du macro-élément statique ne s’ajoutent pas aux degrés de liberté d’interface qui conservent alors un profil de remplissage non complet. Mais on obtient alors un bloc quadrangulaire de degrés de liberté de Lagrange entre les degrés de liberté physiques \({X}_{\Sigma}\) de l’interface \(\Sigma\) et les coordonnées généralisées \({q}_{\Sigma}\) .
Le stockage résultant est représenté sur la figure 7.3.2-a ci-dessous.
Figure 7.3.2-a: Stockage des sous-domaines \(\mathit{E1}\) avec méthode de réduction modale.
On remarquera qu’il n’y a pas de choix de calcul imposé pour la base de réduction de modes d’interface \(\Psi\) . Un choix possible sera celui utilisé en interaction sol-structure [bib15] en considérant une base de modes propres du sous-domaine de type \(\mathit{E1}\) sur tapis de ressorts représentatifs de la rigidité des domaines adjacents de type \(\mathit{E2}\) ou \(\mathit{E3}\) .
Indications pour le contrôle a posteriori de la qualité du modèle réduit#
Un des principes permettant de valider la qualité d’un modèle réduit consiste à estimer, a posteriori, l’écart à l’équilibre qui peut exister après calcul.
Cas général pour un calcul sur un modèle réduit#
Pour caractériser cet écart à l’équilibre, on utilise le principe du résidu en effort. Ce résidu en effort \({F}_{r}\) est calculé
pour les réponses temporelles:
pour les réponses harmoniques :
où \({K}_{e}\) est la matrice de raideur de la structure, \({D}_{v}\) la matrice associée à la dissipation visqueuse, et \(M\) la matrice de masse. On note par ailleurs \({q}_{r}\) le vecteur des degrés de liberté réduit, et \({T}_{r}\) la base de réduction associée. Le terme d’efforts extérieur est défini par le produit d’une matrice de localisation des efforts \(B\) et d’un vecteur \(u\) précisant l’évolution temporelle de l’excitation. Cette écriture permet de séparer les composantes spatiales et temporelles d’un effort souvent noté \(f(t)\) . Concernant la réponse harmonique, \({K}_{h}\) est la matrice représentant le comportement hystérétique, \({K}_{v}\) la matrice représentant le comportement viscoélastique.
On peut également adopter la même démarche pour le calcul des modes propres complexes, détaillée dans la documentation sur la construction des modèles réduits en dynamique [U2.06.04]
Si la solution obtenue est exacte, alors le résidu \({F}_{r}\) est nul, par définition. Dans le cas contraire, on peut construire une mesure de l’écart à l’équilibre à partir de la réponse statique de la structure à ce chargement. On note \({R}_{r}\) le résidu en déplacement associé à chaque mode complexe calculé. Par définition, on a
Une mesure de l’écart à l’équilibre est donnée par l’énergie potentielle élastique de la structure, soit
Un exemple d’utilisation de cette technique est donnée dans la documentation [U2.06.04]
Cas particulier du calcul des modes d’un modèle généralisé#
Le même principe peut être mis en œuvre, mais de façon un peu plus précise, de sorte à faire apparaître des écarts particuliers aux différents équilibres.
Présentation#
Dans le cas d’un calcul par sous structuration avec utilisation de modes d’interfaces, il est important de pouvoir déterminer plus précisément l’origine d’un calcul imprécis. Les sources d’imprécisions sont:
la mauvaise représentation des déplacements aux interfaces, pour chaque sous structure,
le comportement dynamique des parties internes, pour chaque sous structure,
les déplacements différentiels (décollement) pouvant apparaître aux interfaces entre les sous
structures.
La commande CALC_CORR_SSD [U4.52.16] permet, sur la base des calculs des travaux associés à ces différents efforts, de construire des termes de corrections permettant d’améliorer séparément les comportements d’interfaces, et les comportements dynamiques internes de chacune des sous structures.
NB :
Il faudra par ailleurs veiller, dans le processus d’enrichissement, à orthogonaliser séparément les vecteurs à interfaces fixes, et les vecteurs à interface libre, ainsi qu’il est réalisé dans le cas test SDLS122 [V3.02.122], avant de les concaténer et de construire les bases modales.
Calcul des indicateurs et corrections associées#
Pour détailler la construction des différents indicateurs et termes correctifs construits dans la commande CALC_CORR_SSD, considérons un exemple simple constitué de l’assemblage de 2 sous structure le long d’une seule interface. Cette approche se généralise facilement, mais conduit à des notations lourdes.
Définition du problème couplé#
Pour chaque sous structure, on introduit les matrices de masse et de raideur. On note \(c\) les degrés de liberté associés aux parties internes («complémentaires»), et \(i\) les degrés de liberté associés aux interfaces. On introduit aussi la base de réduction associée, composée à la fois de modes à interfaces fixes \({\Phi}_{c}\) , et de modes d’interfaces \({\psi}_{i}\) , relevés statiquement sur le reste de la structure, de sorte que les déplacements correspondants sur les degrés de liberté complémentaires soient \({\psi}_{c}=-{K}_{\mathit{cc}}^{-1}{K}_{\mathit{ci}}{\psi}_{i}\)
Les matrices de masse et de raideur réduites s’écrivent
Du fait de l’introduction de modes d’interfaces, la relation de compatibilité des déplacements d’interfaces entre les sous structures 1 et 2 ne peut plus s’écrire \({q}_{i2}={q}_{i1}\) , mais devient, si la sous structure 2 est esclave :
On définit ensuite une base du noyau des contraintes pour éliminer les relations d’interface, sur la base des éléments présentes dans la section 2.5
Le modèle complet est alors assemblé, puis projeté sur une base du noyau des contraintes:
On note \(\eta\) les degrés de liberté généralisés associé au problème ainsi obtenu. On a donc, si on note \(q\) les degrés de liberté physiques, \(q={T}_{r}{q}_{r}={T}_{r}Q\eta\) .
Écart à l’équilibre#
Si on note \({Z}_{0}={K}_{g}-{\omega}_{0}{M}_{g}\) l’impédance dynamique associée au mode propre de pulsation \({\omega}_{0}\) , et \({\eta}_{0}\) le mode propre réduit associé, on a
Dans ces conditions, la relation permettant de quantifier l’écart à l’équilibre de ce mode particulier devient :
Si la base de projection est de bonne qualité, alors l’écart à l’équilibre doit être faible, et, idéalement, on doit même vérifier
le solveur devant assurer
Pour s’assurer d’avoir un modèle de bonne qualité, on doit donc assurer la nullité de chacun des termes, indépendamment les uns des autres. On doit donc vérifier la nullité des travaux des efforts associé aux charges résiduels dans chaque sous structure, soit
On doit également vérifier la nullité des travaux virtuels associés aux efforts résiduels à l’interface :
Enfin, on doit également vérifier la compatibilité des déplacements, et éviter les déplacements différentiels d’interfaces (décollements)
\(C\) Est l’opérateur de projection des déplacements de l’interface maître sur l’interface esclave. Dans le cas d’interfaces compatibles, \(C\) est la matrice identité.
En pratique, pour être plus simplement exploitables, les différents travaux calculés dans CALC_CORR_SSD sont divisés par la pulsation propre du mode considéré. Cette remise à l’échelle permet de comparer les niveaux d’erreurs entre les différents modes.
Corrections associées aux efforts résiduels#
Chacun de ces terme permet de définir des corrections associées. Ainsi, pour chaque sous structure \(k\) , on construit un sous espace dédié à l’enrichissement basé sur le principe de la correction statique (cf. [U2.06.04]) :
Avec,
Et
Dans le cas où la matrice de raideur serait singulière pour la sous structure \(k\) , on régularise le problème en construisant une matrice shiftée en masse \(K-{\alpha}^{2}M\) , où \(\alpha\) est une valeur faible devant les pulsations propres des premiers modes flexibles de la structure.
Corrections associées aux déplacement différentiels d’interfaces#
Dans le cas où les modes couplés calculés n’assurent pas parfaitement la continuité des déplacements d’interfaces, il est possible d’enrichir les déplacements de part et d’autre de l’interface. Il faut cependant ici considérer deux cas, selon que, pour la sous structure \(k\) , l’interface considérée soit maître ou esclave. Si l’interface est une interface esclave, on se contente alors de projeter le déplacement de l’interface maître sur l’interface esclave, et de relever statiquement ce déplacement sur le reste de la sous structure. La base de réduction complète s’écrit alors
Dans le cas où l’interface considérée est maître, il faut alors extrapoler les mouvements de l’interface esclave. L’approche retenue consiste à construire une base des modes de l’interface maître, et à en rechercher la meilleure sélection qui permette de reconstruire les déplacements de l’interface esclave. Cette combinaison de modes d’interface est alors relevée statiquement sur les degrés de liberté complémentaires de la sous structure considérée. Si on note \(R\) l’opérateur d’extrapolation ainsi définit, la base complète devient, pour une interface maître :
La construction de \(R\) est réalisée comme suit:
On calcule les modes \({\Gamma}_{i\mathit{Mast}}\) de l’interface de la sous structure maître
On cherche la meilleure combinaison linéaire de ces modes qui permette de représenter les déplacements de l’interface esclave. On cherche dont \({\beta}_{0}\) qui réalise
Si la norme de l’erreur est supérieure à 1e-6, alors on calcule un nombre de modes \({\Gamma}_{i\mathit{Mast}}\) supérieur, et on itère tant que la norme d’erreur n’est pas satisfaisante.
\(R\) s’écrit alors formellement
Ces corrections sont calculées pour chaque interface de chaque sous structure. On sépare ensuite chaque famille de vecteurs eu deux familles pour être orthogonalisées :
vecteurs dont les déplacement d’interfaces sont nuls
vecteurs dont les déplacements d’interface sont non nuls
NB :
Il faudra par ailleurs veiller, dans le processus d’enrichissement, à orthogonaliser séparément les vecteurs à interfaces fixes, et les vecteurs à interface libre, ainsi qu’il est réalisé dans le cas test SDLS122 [V3.02.122]
Bibliographie#
ROY, J. CRAIG & M. C. BAMPTON: « Coupling of Substructures for Dynamic Analysis » AIAA Journal, (July 1968), Vol. 6, N° 7, p. 1313-1319.
MAC NEAL: « A hybrid method of component mode synthesis », Computers and Structures, (1971), Vol. 1, p. 581-601.
IMBERT: « Analyse des Structures par Éléments Finis » 1979, Cepadues Édition.
RICHARD: « Méthodes de sous-structuration dynamique en éléments finis ». Rapport EDF HP-61/90.149.
RICHARD: « Méthodes de sous-structuration dans le Code_Aster « . Rapport EDF HP‑61/92.149.
PELLET: « Code de Mécanique Aster - Manuel de référence : Les éléments finis dans Aster « Clé : R3.03.01 « Dualisation des conditions aux limites ».
JACQUART: « Code de Mécanique Aster - Manuel de référence : Dynamique en base modale ». Clé : R5.06.01 « Méthodes de RITZ en dynamique linéaire et non-linéaire ».
KERBER: « Cahier des charges: mise en œuvre de la sous-structuration harmonique » - Rapport D.E.R. HP-61/93.053.
KERBER: « Sous-structuration harmonique dans le Code_Aster « - Rapport D.E.R. HP‑61/93.104.
VARE: « Cahier des charges de l’implémentation du calcul transitoire linéaire par sous‑structuration dynamique dans le Code_Aster « - Rapport D.E.R. HP‑61/94.135/B
VARE: « Documentations Utilisateur et Validation des opérateurs de calcul transitoire par sous-structuration » - Rapport D.E.R. HP-61/94.208/A.
BOURQUIN & F. D’HENNEZEL : “Numerical study of an intrinsic component mode synthesis method” – Mathematical Modelling and Numerical Analysis, 1992.
NICOLAS: Spécification du développement de la traduction MACR_ELEM_STAT / MACR_ELEM_DYNA. Compte rendu EDF/AMA-05.149.
Document Aster [U2.07.02]:Notice d’utilisation de la sous-structuration statique.
Document Aster [U2.06.07]: Interaction sol-structure en analyse sismique avec l’interface
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
4 |
P.RichardEDF-R&D/AMV |
Texte initial |
5 |
G.Rousseau, C.VaréEDF-R&D/AMV |
|
8.4 |
C.BodelEDF-R&D/AMA |
|
10.1 |
G.DEVESA EDF-R&D/AMA |
Ajout du §7 sur la condensation dynamique par sous-structuration statique |
10.2 |
M.CORUS EDF-R&D/AMA |
Ajout d’une section sur l’élimination des contraintes suite aux développement autour des modèles généralisés. |
10.4 |
F.VOLDOIRE EDF-R&D/AMA |
Corrections apportées aux §5.3.3, 5.2, 6.2.2, 6.2.5, 7. |
11.1 |
M.CORUS EDF-R&D/AMA |
Ajout d’une section sur les indicateurs de qualité pour les modèles réduits suite aux développement autour des modèles généralisés. Compléments de la section 2 sur le calcul des modes d’interfaces. |