r3.03.05 Élimination des conditions aux limites dualisées

Résumé :

La prise en compte des conditions aux limites dans c ode_aster peut être réalisée soit par une technique directe d’élimination dans des cas simples (AFFE_CHAR_CINE), soit grâce à une technique de dualisation (AFFE_CHAR_MECA), et à l’introduction de multiplicateurs de Lagrange, pour les cas les plus généraux.

Néanmoins, cette approche présente deux inconvénients. D’une part, l’ajout de multiplicateurs de Lagrange augmente le nombre de degré de liberté du problème global, et donc la taille du système à résoudre. Ce point peut devenir pénalisant dès lors que le nombre de conditions aux limites augmente (nombreuses interfaces, liaisons cinématiques, impositions de mouvements rigides, etc.). D’autre part, la technique particulière de dualisation conduit à la perte du caractère semi-défini positif de la matrice de raideur. Ce changement de nature peut conduire à des défauts de robustesse, voir à l’échec de la résolution dans le cas des solveurs itératifs.

On présente dans ce document une méthode permettant d’éliminer a posteriori les multiplicateurs de Lagrange pour prendre en compte les conditions aux limites affines. Cette approche ne concerne pas les multiplicateurs de Lagrange liés à la modélisation du contact.

L’élimination des conditions aux limites est réalisée soit par l’utilisation du mot-clef ELIM_LAGR dans SOLVEUR [U4.50.01], soit directement avec l’opérateur ELIM_LAGR [U4.55.03] qui construit un nouveau conceptmatr_asse_elim.

Introduction

La prise en compte des conditions aux limites dans Code_Aster peut être réalisée soit par une technique directe d’élimination dans des cas simples (AFFE_CHAR_CINE), soit grâce à une technique de dualisation (AFFE_CHAR_MECA), et à l’introduction de multiplicateurs de Lagrange, pour les cas les plus généraux. Cette dernière approche est très générale et permet de traiter sous un même formalisme tous les types de conditions aux limites.

Néanmoins, cette approche présente deux inconvénients. D’une part, l’ajout de multiplicateurs de Lagrange augmente le nombre de degré de liberté du problème global, et donc la taille du système à résoudre. Ce point peut devenir pénalisant dès lors que le nombre de conditions aux limites augmente (nombreuses interfaces, liaisons cinématiques, impositions de mouvements rigides, etc.). D’autre part, la technique particulière de dualisation conduit à la perte du caractère semi-défini positif de la matrice de raideur. La recherche de la solution n’est alors plus une recherche de minimum, mais une recherche de point selle []. Ce changement de nature peut conduire à des défauts de robustesse, voir à l’échec de la résolution dans le cas des solveurs itératifs.

On présente dans ce document une méthode alternative à l’introduction de multiplicateurs de Lagrange pour prendre en compte les conditions aux limites affines. Cette approche ne se substitue pas à l’introduction des multiplicateurs de Lagrange, en particulier pour les études présentant du contact. En revanche, dans le cas des études où les conditions aux limites sont figées, cette approche permet des gains de robustesse, et de performance dans le cas d’un nombre important de conditions aux limites. On rappelle dans un premier temps la technique générale qui préside à l’élimination des conditions aux limites, puis on présente la démarche retenue pour éliminer en pratique les conditions aux limites, ainsi que les différents algorithmes associés. Enfin, on présente la technique d’élimination dans le cas d’un problème aux valeurs propres.

Principe de la résolution par élimination

Définition du problème à résoudre

Le système à résoudre peut se mettre sous la forme générale suivante:

(320)\[\left[K\right]\lbrace u\rbrace =\lbrace f\rbrace\]

sous la contrainte

(321)\[\left[B\right]\lbrace u\rbrace =\lbrace {u}_{0}\rbrace\]

Avec l’introduction de multiplicateurs de Lagrange, le problème devient

(322)\[\begin{split}\lbrace \begin{array}{c}\left[K\right]\lbrace u\rbrace +{\left[B\right]}^{T}\lbrace \lambda \rbrace =\lbrace f\rbrace \\ \left[B\right]\lbrace u\rbrace =\lbrace {u}_{0}\rbrace \end{array}\end{split}\]

et se met sous la forme générale

(323)\[\begin{split}\left[\begin{array}{cc}K& {B}^{T}\\ B& 0\end{array}\right]\left\lbrace \begin{array}{c}u\\ \lambda \end{array}\right\rbrace =\left\lbrace \begin{array}{c}f\\ {u}_{0}\end{array}\right\rbrace\end{split}\]

Dans le cas de la double dualisation (AFFE_CHAR_MECA), cette matrice est augmentée, et le système à résoudre devient

(324)\[\begin{split}\left[\begin{array}{ccc}K& {B}^{T}& {B}^{T}\\ B& -\alpha {I}_{d}& \alpha {I}_{d}\\ B& \alpha {I}_{d}& -\alpha {I}_{d}\end{array}\right]\left\lbrace \begin{array}{c}u\\ {\lambda}_{1}\\ {\lambda}_{2}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}f\\ {u}_{0}\\ {u}_{0}\end{array}\right\rbrace\end{split}\]

après avoir posé naturellement \({\lambda}_{1}={\lambda}_{2}=\lambda /2\) .

Pour rendre son caractère semi-défini positif au problème générique \(\mathit{Ku}=f\) , dans le cas où la matrice \(K\) présente une topologie découlant de la mise sous la forme () ou (), on va chercher à décomposer la solution \(u\) du problème sous la forme \(u={u}_{p}+{u}_{g}\) , où

• \({u}_{p}\) est la solution particulière du problème,

• \({u}_{g}\) est la solution générale du problème, associée aux conditions aux limites \([B]\lbrace {u}_{g}\rbrace =\lbrace 0\rbrace\) .

Recherche de la solution particulière

La recherche de la solution particulière consiste simplement à trouver \({u}_{p}\) tel que \([B]\lbrace {u}_{p}\rbrace =\lbrace {u}_{0}\rbrace\) . Pour calculer \({u}_{p}\) , on considère la décomposition sur l’image et le noyau de \(B\) , soit \({u}_{p}={u}_{p-i}+{u}_{p-n}\) , avec \({u}_{p-i}\) appartenant à l’image de \(B\) et \({u}_{p-n}\) appartenant au noyau de \(B\) .

Ce problème peut être résolu au sens des moindres carrés, et dans ce cas, on considère la solution de () de norme \({L}_{2}\) minimale, soit \({u}_{p}={u}_{p-i}\) . \(\lbrace {u}_{p-i}\rbrace\) est la solution du problème de minimisation:

(325)\[\lbrace {u}_{p-i}\rbrace =\underset{\lbrace {u}_{p}\rbrace }{\text{ArgMin}}\left({\Vert [B]\lbrace {u}_{p}\rbrace -\lbrace {u}_{0}\rbrace \Vert }^{2}\right)\]

Pour calculer efficacement \({u}_{p-i}\) , sous réserve que \(B\) soit de rang plein (rang égal à \(\mathit{Nb}\) , donc pas de contraintes redondantes), on exploite la propriété suivante:

(326)\[\left([B]{[B]}^{T}\right){\left([B]{[B]}^{T}\right)}^{-1}={\left[{I}_{d}\right]}_{\mathit{NbxNb}}\]

et donc

(327)\[\left([B]{[B]}^{T}\right){\left([B]{[B]}^{T}\right)}^{-1}\lbrace {u}_{0}\rbrace =\lbrace {u}_{0}\rbrace\]

La solution recherchée est naturellement donnée par:

(328)\[\lbrace {u}_{p}\rbrace =\lbrace {u}_{p-i}\rbrace ={[B]}^{T}{\left([B]{[B]}^{T}\right)}^{-1}\lbrace {u}_{0}\rbrace\]

En pratique, on effectue une factorisation de Cholesky de la matrice symétrique \({\mathit{BB}}^{T}\) , à l’aide de la bibliothèque MUMPS, puis on résout

(329)\[[B]{[B]}^{T}\lbrace y\rbrace =\lbrace {u}_{0}\rbrace\]

Enfin, on obtient \(\lbrace {u}_{p}\rbrace\) :

\(\lbrace {u}_{p}\rbrace =\lbrace {u}_{p-i}\rbrace ={[B]}^{T}\lbrace y\rbrace\)

Recherche de la solution générale

La solution générale \({u}_{g}\) vérifiant \([B]\lbrace {u}_{g}\rbrace =\lbrace 0\rbrace\) , elle appartient au noyau de \(B\) . La recherche de la solution générale passe donc par une étape de construction d’une base \(Z\) du noyau de \(B\) . On peut ensuite projeter le problème sur le noyau, et résoudre ce problème, mieux posé et de taille réduite.

Obtention d’une base du noyau

Formellement, la construction du noyau de \(B\) peut se faire de plusieurs façons:

• Décomposition en valeurs singulières de \(B\) , soit \(B={\mathit{USV}}^{T}\) . \(\mathit{Ker}(B)\) est l’espace engendré par les colonnes de \({V}^{T}\) associées à une valeur singulière nulle.

• Décomposition \(\mathit{QR}\) de \({B}^{T}\) , soit \(B={R}^{T}{Q}^{T}\) . \(\mathit{Ker}(B)\) est l’espace engendré par les colonnes de \(Q\) associées à une valeur nulle sur la diagonale de \({R}^{T}\) .

• Décomposition \(\mathit{LU}\) de \({B}^{T}\) . Le calcul du noyau n’est pas direct, et s’obtient à partir de sous-blocs de \(L\) .

C’est cette troisième méthode qui est mise en œuvre dans code_aster. On utilise la bibliothèque SuperLU pour construire la factorisation \(\mathit{LU}\) de \({B}^{T}\) :

\[ \begin{align}\begin{aligned}:label: eq-11\\ {P}_{r}{B}^{T}{P}_{c}^{T}=\left[L\right]\left[U\right]=\left[\begin{array}{cc}{[{L}_{1}]}_{{N}_{b}\times {N}_{b}}& [:ref:\\Les matrices :math:`{P}_{r}` et :math:`{P}_{c}` désignent respectivement les matrices de permutation imposées aux lignes et aux colonnes de :math:`{B}^{T}` pour la construction de la factorisation :math:`\mathit{LU}` . Les blocs :math:`{L}_{1}` et :math:`{L}_{3}` sont triangulaires inférieurs et le bloc :math:`{U}_{1}` est triangulaire supérieur.\\Remarque:\\En pratique, :math:`{P}_{c}=\mathit{Id}` , il n’y a pas de permutation des colonnes de :math:`{B}^{T}` .\\On en déduit l’expression suivante pour :math:`B` :\end{aligned}\end{align} \]

(330)\[B={U}^{T}{L}^{T}{P}_{r}=\left[\begin{array}{cc}{[{U}_{1}]}_{{N}_{b}\times {N}_{b}}^{T}& {[:ref:\]

Il est ensuite facile de voir qu’une base de \(\mathit{Ker}(B)\) est donnée par:

(331)\[\begin{split}Z={P}_{r}^{-1}\left[\begin{array}{c}-{[{L}_{1}]}^{-T}{[{L}_{2}]}^{T}\\ {[{I}_{d}]}_{(N-{N}_{b})\times (N-{N}_{b})}\end{array}\right]\end{split}\]

On fait ici l’hypothèse que \(B\) est de rang plein. Si tel n’est pas le cas, alors il suffit de construire \({L}_{1}\) à partir du bloc de \(L\) associé aux termes non nuls de la diagonale de \(U\) .

Projection et résolution

Une fois que la base du noyau est calculée, la résolution est directe. On pose \({u}_{g}=Z\overline{u}\) , et la solution du problème est donc directement donnée par

(332)\[{[Z]}^{T}[K][Z]\lbrace \overline{u}\rbrace ={Z}^{T}\lbrace f\rbrace\]

puisque, par construction,

(333)\[{[Z]}^{T}{[B]}^{T}\lbrace \lambda \rbrace =\lbrace 0\rbrace\]

Dans certains cas, il est nécessaire d’accéder aux multiplicateurs de Lagrange. Il faut donc les reconstruire, à partir des solutions \({u}_{p}\) et \({u}_{g}\) . Pour ce faire, on doit résoudre la première ligne du système (), soit

(334)\[{[B]}^{T}\lbrace \lambda \rbrace =\lbrace f\rbrace -[K]\left([Z]\lbrace \overline{u}\rbrace +\lbrace {u}_{p}\rbrace \right)\]

Si \(B\) est de rang plein, la solution s’écrit formellement:

(335)\[\lbrace \lambda \rbrace ={\left([B]{[B]}^{T}\right)}^{-1}[B]\lbrace f\rbrace -[K]\left([T]\lbrace \overline{u}\rbrace +\lbrace {u}_{p}\rbrace \right)\]

et se calcule en pratique en réutilisant la factorisationde \(B{B}^{T}\) calculée pour la recherche de la solution particulière. De façon analogue, si \(B\) n’est pas de rang plein, il est possible de supprimer les contraintes redondantes.

Modes propres et élimination

Cette section fait écho à la section §7 de la documentation R3.03.01, et justifie la méthode d’élimination pour le calcul des modes propres. On s’intéresse ici au problème général suivant:

(336)\[\textrm{Trouver les couples} \left({\omega}_{0}^{2},\lbrace \varphi \rbrace \right)\]

Les matrices \({M}_{c}\) et \({K}_{c}\) doivent prendre en compte les conditions aux limites, et sont ici supposées obtenues à partir de fonctions de formes \({u}_{N}\) vérifiant \({\mathit{Bu}}_{N}=0\) . Or, dans le cas présent, le problème assemblé n’est pas le problème (), mais un problème de taille plus importante, pour lequel les fonctions de formes ne vérifient pas les contraintes cinématiques. Les matrices \({M}_{b}\) et \({K}_{b}\) sont de dimension \((N-\mathit{Nb})\times (N-\mathit{Nb})\) .

En partant des matrices assemblées \(M\) et \(K\) , de taille \(N\times N\) , construites sur la base de fonctions de formes \(u\) ne vérifiant pas \(\mathit{Bu}=0\) , il existe deux alternatives pour résoudre le problème (). La première consiste à augmenter la matrice \(K\) pour faire apparaître les contraintes cinématiques, en ajoutant des inconnues sous la forme de multiplicateurs de Lagrange. C’est l’approche retenue jusqu’à présent dans le code. Le problème à résoudre s’écrit, avec la technique de double dualisation:

(337)\[\begin{split}\left(\left[\begin{array}{ccc}K& {B}^{T}& {B}^{T}\\ B& -\mathit{Id}& \mathit{Id}\\ B& \mathit{Id}& -\mathit{Id}\end{array}\right]-{\omega}_{0}^{2}\left[\begin{array}{ccc}M& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\right)\left\lbrace \begin{array}{c}\varphi \\ {\lambda}_{1}\\ {\lambda}_{2}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}0\\ 0\\ 0\end{array}\right\rbrace\end{split}\]

L’approche par ajouts de multiplicateurs de Lagrange, dans le cas de la recherche des modes propres, n’est pas satisfaisante, puisqu’en plus de la perte des propriétés de positivité de l’opérateur, on ajoute un nombre important de degrés de liberté, et on élargit donc le spectre du problème. Cet élargissement de spectre pose des problèmes numériques, qui obligent à un travail important dans ce cas.

La seconde approche, celle proposée ici, est plus naturelle, et se rapproche du problème initial. Les matrices de masse et de raideur sont certes assemblées sur la base de fonctions de formes ne vérifiant pas les conditions limites, mais la recherche des modes propres et des valeurs propres peut se faire dans un sous espace adapté. Il suffit de construire une base du sous espace des vecteurs \(v\) vérifiant \(\mathit{Bv}=0\) . Ce sous espace est naturellement le noyau de \(B\) . Il suffit alors de rechercher les modes du problème non contraint projeté dans le noyau de \(B\) . Soit \(Z\) une base du noyau, on cherche alors les couples \(\left({\omega}_{0}^{2},\lbrace \psi \rbrace \right)\) qui vérifient

(338)\[\left[{Z}^{T}\left(K-{\omega}_{0}^{2}M\right)Z\right]\left\lbrace \begin{array}{c}\psi \end{array}\right\rbrace =\left\lbrace \begin{array}{c}0\end{array}\right\rbrace\]

On identifie alors \({M}_{b}={Z}^{T}MZ\) et \({K}_{b}={Z}^{T}KZ\) .

Après résolution du problème modal réduit, on peut calculer la valeur des multiplicateurs de Lagrange :math:`lambda ` par la relation :

(339)\[\left\lbrace \lambda \right\rbrace = - \left(BB^T\right)^{-1} B\left(K-{\omega}_{0}^{2}M\right)Z \left\lbrace \psi \right\rbrace\]

Bibliographie

[1] Benzi, M. and Golub, G.H., and Liesen, J.:Numerical solution of saddle point problems, Acta Numerica, Vol. 14, pp1 - 137 (2005).