# Méthode Arlequin et décomposition de domaine :
aspects théoriques et logiciels

Josselin Touzeau, Vincent Chiaruttini, Frédéric Feyel, Hachmi Ben Dhia

► To cite this version:

Josselin Touzeau, Vincent Chiaruttini, Frédéric Feyel, Hachmi Ben Dhia. Méthode Arlequin et décomposition de domaine : aspects théoriques et logiciels. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01413605

HAL Id: hal-01413605

https://hal.science/hal-01413605v1

Submitted on 10 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Méthode Arlequin et décomposition de domaine : aspects théoriques et logiciels

J.Touzeau¹, V.Chiaruttini¹, F.Feyel¹, H.Ben Dhia²

¹ Onera (DMSM/CEMN)
29 avenue de la division Leclerc, 92322 CHATILLON
{josselyn.touzeau,vincent.chiaruttini,frederic.feyel}@onera.fr

² Ecole Centrale Paris (LMSS-MAT)
Grande voie des vignes, 92295 CHATENAY-MALABRY
hachmi.ben-dhia@ecp.fr

**Résumé** — Nous présentons quelques aspects concernant l'utilisation d'une approche de type Arlequin pour la simulation numérique de problèmes multiéchelles. Pour assurer une résolution performante, un solveur de type décomposition de domaine est développé à partir d'une approche standard sans recouvrement. L'adaptation au cas recouvrant est détaillée pour permettre la gestion des déplacements de solides rigides, la prise en compte de maillages non-compatibles et l'utilisation de préconditionneurs. Cette approche a été mise en œuvre au sein du code éléments finis ZéBuloN et les premières expérimentations numériques sont présentées.

**Mots clés** — multiéchelle, décomposition de domaine, méthode Arlequin.

## 1 Introduction

Les objectifs de réduction de poids, dans un contexte de préservation d'une sécurité maximale, conduisent les industriels de l'aéronautique à devoir utiliser des modèles de plus en plus élaborés pour évaluer finement la durée de vie des composants critiques. La complexité des modèles physiques entrant en jeu devient telle, qu'elle conduit nécessairement à développer de nouvelles approches numériques spécifiques à leur intégration dans les codes de calcul. L'objectif des méthodes multiéchelles et multimodèles est de répondre à ce besoin : permette la prise en compte efficace d'un ensemble de phénomènes physiques de natures différentes, interagissant, mais intervenant à différents niveaux d'espace et/ou de temps.

Dans cette voie, plusieurs approches ont été développées récemment pour traiter ce type de problèmes multi-échelles. Par exemple, les décompositions de domaine de type micro-macro, comme la méthode LATIN [1], réalisent un partitionnement de structures en sous-structures et interfaces donnant la possibilité de scinder les champs de déplacement et de contraintes en champs microscopique et macroscopique. L'approche FE² [2] permet d'écrire des problèmes mécaniques non-linéaires à une échelle microscopique pour ensuite, par l'intermédiaire d'équations de localisation et d'homogénéisation, obtenir des relations de comportement pour mener la résolution à l'échelle de la structure. Une autre approche, la méthode Arlequin [3, 4], sur laquelle nous avons basé nos travaux, permet de superposer, dans des zones d'intérêt, des modèles complexes sur une structure possédant un comportement global souvent plus simple, en effectuant une partition d'énergie et en assurant la communication entre les modèles.
La stratégie de résolution proposée dans la suite conduit à condenser le problème Arlequin au niveau des interfaces et à utiliser un solveur itératif de Krylov [5]. Les aspects concernant le traitement d'une interface recouvrante, des incompatibilités de maillage, de la transmission des déplacements de solides rigides et du calcul des compléments de Schur sont présentés.

La section 2 rappelle les principes de l'approche Arlequin, la section 3 décrit les aspects concernant le solveur décomposition de domaine, la section 4 présente la mise en œuvre numérique de la méthode au sein du code ZéBuloN<fnref n="1" /> et les premiers résultats de simulations numériques.

## 2 La méthode Arlequin

L'approche Arlequin repose sur le couplage de différents modèles à l'aide d'un opérateur de collage. Cet opérateur assure la transmission des efforts et la continuité du champ de déplacement au sens faible dans les zones de superposition.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Problème multimodèle constitué de deux domaines.</figcaption></figure>

Pour simplifier, on présente l'approche dans le cadre d'un problème modélisé à l'aide de deux formulations distinctes sur deux domaines ($\Omega_1$ et $\Omega_2$) possédant une intersection non vide ($\Omega_1 \cap \Omega_2 \neq \emptyset$) représenté sur la figure 1. Pour chaque domaine $\Omega_i$ de $\mathbb{R}^n$ ($n = 2$ ou 3) est posé un problème classique d'élasticité linéaire caractérisé par : un déplacement imposé $u_i^d$ sur une partie $\partial_u \Omega_i$ de la frontière $\partial\Omega$, un effort $F_i^d$ imposé sur la partie complémentaire $\partial_F \Omega_i$ et un tenseur d'élasticité $A_i$. Sur chaque domaine $\Omega_i$, la formulation variationnelle du système s'écrit, à l'aide des espaces cinématiquement admissibles usuels $U_i^{ad}$ :

$$ \int_{\Omega_i} \sigma_i : \varepsilon(v_i) \, d\Omega_i = \int_{\partial_F \Omega_i} F_i^d \cdot v_i \, dS_i, \quad \forall v_i \in U_i^{ad0} \quad (1) $$

où $\sigma_i$ représente le tenseur des contraintes sur le domaine $\Omega_i$ et $\varepsilon(u_i) = 1/2 (\nabla u_i + \nabla^T u_i)$ le tenseur des déformations linéarisé.

Pour assurer la liaison entre les systèmes, un terme de couplage est ajouté sur une zone $S_c \subset \Omega_1 \cap \Omega_2$ :

$$ C(\lambda, \mu) = \int_{S_c} K (\nabla \lambda \nabla \mu + 1/l^2 \lambda \mu) \, d\Omega, (\lambda, \mu) \in \mathcal{M}^2 \quad (2) $$

où $\mathcal{M}$ est un espace admissible régulier défini sur $S_c$, $l$ une longueur caractéristique et $K$ un coefficient de l'ordre de grandeur la rigidité des milieux en présence.

Des champs de pondération $\alpha_i$ sont introduits pour assurer la répartition des énergies de même nature entre les différents modèles (en respectant $\alpha_1 + \alpha_2 = 1$ sur $\Omega_1 \cup \Omega_2$).

À l'aide des équations précédentes, il est possible d'établir l'écriture variationnelle la formu-

<footnote n="1">¹Disponible sur Internet à http://www.nwnumerics.com.</footnote>
lation Arlequin :

$$
\begin{align}
\int_{\Omega_1} \alpha_1 \sigma_1 : \varepsilon(v_1) \, d\Omega_1 + C(\lambda, v_1) &= \int_{\partial_F \Omega_1} \alpha_1 F_1^d \cdot v_1 \, dS_1, \quad \forall v_1 \in \mathcal{U}_1^{ad0} \nonumber \\
\int_{\Omega_2} \alpha_2 \sigma_2 : \varepsilon(v_2) \, d\Omega_2 - C(\lambda, v_2) &= \int_{\partial_F \Omega_2} \alpha_2 F_2^d \cdot v_2 \, dS_2, \quad \forall v_2 \in \mathcal{U}_2^{ad0} \tag{3} \\
C(u_1, \mu) - C(u_2, \mu) &= 0, \quad \forall \mu \in \mathcal{M} \nonumber
\end{align}
$$

Pour résoudre numériquement le système variationnel précédent, une discrétisation de type
éléments finis est retenue et conduit à la résolution du système linéaire suivant :

$$
\begin{bmatrix} K_1 & 0 & C_{12}^T \\ 0 & K_2 & -C_{21}^T \\ C_{12} & -C_{21} & 0 \end{bmatrix} \cdot \begin{pmatrix} U_1 \\ U_2 \\ \Lambda \end{pmatrix} = \begin{pmatrix} F_1 \\ F_2 \\ 0 \end{pmatrix} \quad (4)
$$

où $K_1, K_2$ sont les matrices de rigidité pondérées, $C_1, C_2$ les matrices de couplage et $F_1, F_2$ les forces généralisées pondérées.

3 Solveur décomposition de domaine associé à une modélisation Ar-
lequin

Le calcul du système discrétisé (4), et plus encore sa résolution, présentent des difficultés
importantes :

- dans une utilisation générale de la méthode, des maillages incompatibles se superposent au niveau des zones de collage, ce qui impose une attention particulière dans le calcul des blocs de collage $C_{ij}$;

- les champs de pondération doivent être appliqués en veillant à ne pas trop détériorer le conditionnement des blocs $K_i$ ;

- les blocs $K_i$ peuvent posséder une singularité selon le type de condition imposée sur le sous-domaine, il convient alors de déterminer l’élément du noyau qui satisfait au mieux à la physique du problème ;

- la factorisation du système complet, parce qu'il possède de nombreux termes nuls sur la diagonale, requiert une stratégie de pivotement partiel que ne possèdent pas directement les solveurs directs traditionnels efficaces (de type dissection emboîté ou dscpack, par exemple).

Les deux premiers points seront abordés dans la section suivante consacrée à l'architecture de la mise en œuvre efficace de l'approche Arlequin dans le code orienté objet ZéBuloN. Concernant les aspects solveurs, nous avons d'une part développé un outil simple permettant la transmission des mouvements de solide rigide sur les domaines flottants, et d'autre part mené des développements en vue de proposer une démarche résolution itérative basée sur une méthode de type décomposition de domaine.

Pour traiter efficacement le raccord d'un domaine flottant i avec un domaine bien positionné j, une nouvelle condition de liaison a été développée. L'ensemble des champs de déplacement admissibles du domaine flottant T<sup>ad</sup> est défini à partir d'une bijection de Ker K<sub>i</sub>. Puis un problème de minimisation est écrit pour assurer au domaine flottant le mouvement de corps rigide traduisant au mieux le déplacement u(x<sub>p</sub>) d'un ensemble de nœuds x<sub>p</sub> du maillage j au sens des moindres carrés ; soit à déterminer le champ de vecteur torsoriel τ(x) ∈ T<sup>ad</sup> tel que :

$$
\text{Min}_{\tau(x) \in \mathcal{T}^{\text{ad}}} \sum_{k=1}^{p} \| \tau(x) - u(x_p) \|^2 \quad (5)
$$

Dans la démarche de résolution, le traitement de ce type de raccord est mené de la même façon
que la transmission de modes rigides.
La stratégie retenue pour la résolution du système (4) repose sur l'utilisation d'une méthode de type décomposition de domaine. Habituellement, ce type d'approche est retenu dans le cas de domaines non-recouvrants et d'interfaces compatibles. Néanmoins, d'un point de vue algébrique il n'existe pas d'obstacle particulier au traitement par décomposition de domaine du système linéaire produit par une modélisation Arlequin. Le système d'équation (4) peut se décomposer en deux parties :

$$ \begin{bmatrix} K_1 & C_{12}^T \\ C_{12} & 0 \end{bmatrix} \cdot \begin{bmatrix} U_1 \\ \Lambda \end{bmatrix} = \begin{bmatrix} F_1 \\ \Xi \end{bmatrix}, \quad \begin{bmatrix} K_2 & -C_{21}^T \\ -C_{21} & 0 \end{bmatrix} \cdot \begin{bmatrix} U_2 \\ \Lambda \end{bmatrix} = \begin{bmatrix} F_2 \\ -\Xi \end{bmatrix} \quad (6) $$

où $\Lambda$ représente les efforts de collage partagés par chacun des domaines et $\Xi$ les efforts créés par les déplacements locaux via l'opérateur de collage.

Si l'on considère (pour simplifier) que les matrices de rigidité $K_i$ sont régulières, sur le domaine 1, en utilisant la première ligne il est possible d'isoler $U_1$ et d'établir :

$$ -C_{12} K_1^{-1} C_{12}^T \Lambda = S_1^p \Lambda = \Xi - C_{12} K_1^{-1} F_1 \quad (7) $$

où $S_i^p$ est le complément de Schur primal du domaine 1. Le complément de Schur dual est obtenu en isolant $\Lambda$ dans (7), soit :

$$ [-C_{12} K_1^{-1} C_{12}^T]^{-1} \Xi = S_1^d \Xi = \Lambda + S_1^d C_{12} K_1^{-1} F_1 \quad (8) $$

Les problèmes condensés primal et dual s'écrivent alors respectivement :

$$ (S_1^p + S_2^p) \Lambda = S^p \Lambda = F_\Lambda \quad (9) $$

$$ (S_1^d - S_2^d) \Xi = S^d \Xi = F_\Xi \quad (10) $$

Comme dans le cas des décompositions de domaine sans recouvrement, la résolution de ces équations est obtenue à l'aide d'un solveur itératif de Krylov auquel il est appliqué les augmentations nécessaires au traitement des éventuels mouvements de corps rigides.

Les bonnes performances du solveur vont être liées à l'association de préconditionnneurs efficaces. Dans le cas d'une approche primale, l'inverse du complément de Schur primal global est approximée par la somme pondérée des inverses des compléments de Schur primaux locaux soit la somme des compléments de Schur duaux (et réciproquement en approche primale) :

$$ (\hat{S}^p)^{-1} = (S_1^p)^{-1} + (S_2^p)^{-1} = S_1^d - S_2^d \quad (11) $$

$$ (\hat{S}^d)^{-1} = (S_1^d)^{-1} - (S_2^d)^{-1} = S_1^p + S_2^p \quad (12) $$

En pratique, l'utilisation des préconditionneurs conduit, dans la démarche itérative, à effectuer de nombreux produits matriciels avec les compléments de Schur. Si les modes rigides sont correctement pris en compte, les termes diagonaux nuls dans le complément de Schur dual (8) rendent sa factorisation problématique. Pour permettre ce traitement, nous avons modifié la numérotation du système de façon à assurer la factorisation des termes d'interface après celle des termes internes.

# 4 Mise en œuvre numérique

L'approche Arlequin a été mise en œuvre au sein du code ZéBuloN co-développé par l'Onera. Les efforts de collage étant de même nature mathématique qu'un champ de déplacement, nous avons choisi de définir un maillage spécifique pour introduire ces inconnues et écrire l'opérateur de couplage C. L'intégration dans le code de calcul a été ensuite effectuée par le développement de deux classes :

- un objet permettant de lier la zone de collage au problème élément fini correspondant en procédant à une recherche des éléments intersectant avec ceux de la zone de collage ;
- un objet élément Arlequin basé sur les éléments de la zone de collage qui permet d'écrire le terme de collage et de supprimer l'énergie excédentaire.

Pour traiter d'éventuels maillages incompatibles, le nombre de points de Gauss sur le maillage de la zone de collage est augmenté de façon significative (par exemple, jusqu'à 81 points de Gauss pour des éléments triangulaires), et la formulation de l'opérateur de collage est menée en s'appluyant sur les techniques mortier. La fonction de pondération est calculée explicitement sur le maillage d'interface à l'aide d'une distance entre les zones où la fonction possède une valeur constante.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Fissuration d'une plaque en traction - interface compatible.</figcaption></figure>

Pour valider la démarche, nous avons commencé par étudier le problème d'une plaque fissurée, considérée comme une structure plane en acier, ayant un comportement linéaire élastique (module d'Young de 200000 MPa, coefficient de poisson de 0,3) et soumise à un chargement de traction constante (cf. figure 2). La discrétisation est menée à l'aide d'une interface compatible, des fonctions d'interpolation linéaires sont utilisées sur le maillage complet (totalisant 5 165 degrés de liberté) alors que des fontions quadratiques sont employées sur le patch (9010 inconnues), la zone de collage utilise un maillage avec des fonctions d'interpolation linéaires (1326 inconnues).

<figure><img src="image_4.png" /><figcaption>Figure 3 – Convergence des normes relatives des résidus pour les différentes approches. À gauche pour le cas d'une interface compatible, à droite, pour une interface non-compatible.</figcaption></figure>

Nous avons comparé les résultats obtenus avec les différentes approche duale/primale et leurs variantes préconditionnées, sur la figure 3. Ces résultats présentent le nombre d'itérations à convergence des différentes approches pour un niveau d'erreur relative GMRes identique de $10^{-8}$. Dans ce cas d'interface compatible, les approches primale et duale non-préconditionnées présentent des comportements très différents, les opérateurs d'interfaces produisant un complément de Schur dual a priori mieux conditionné. L'utilisation des préconditionneurs dans ce cas conduit à des performances comparables pour les deux approches. On peut noter que dans tous les cas les erreurs relatives produites par l'approche Arlequin (en comparaison avec la solution obtenue à l'aide d'un seul maillage) sont en déplacements, en moyenne sur le patch, inférieure à 0,2%.
<figure><img src="image_5.png" /><figcaption>Figure 4 – Fissuration d'une plaque en traction - interface incompatible.</figcaption></figure>

La méthode Arlequin nous permet de traiter des cas de raccord non-compatibles, nous avons étudié le cas d'un plaque fissurée en traction avec un raccord non-compatible (cf. figure 4). La résolution de ce problème avec les différentes démarches (cf. figure 3) illustre la difficulté, dans ce contexte, d'obtenir des préconditionneurs performants, les opérateurs collage influant très négativement sur le conditionnement des compléments de Schur.

# 5 Conclusion et perspectives

Nous avons présenté plusieurs aspects concernant la mise en œuvre d'une approche multi-échelle de type Arlequin associée à un solveur par décomposition de domaine. Pour mener à bien cette démarche plusieurs outils ont été développés pour la transmission des modes rigides, la facturation des opérateurs à termes diagonaux nuls ou la prise en compte des raccords incompatibles.

Les premiers résultats obtenus avec cette démarche sont encourageant, mais il démontre la difficulté de prise en compte des interfaces incompatibles dans le cadre d'un solveur itératif. C'est dans cette direction que nous envisageons de poursuivre nos travaux, à travers le calcul de termes couplage mieux conditionnées, ainsi que dans la recherche de préconditionneurs simplifiés.

# Remerciements

Les auteurs tiennent à remercier Pierre Gosselet pour ses conseils lors de l'adaptation du solveur parallèle de ZéBuloN à la résolution d'un problème Arlequin.

# Références

[1] P. Ladevèze, D. Dureisseix. Une nouvelle stratégie de calcul micro/macro en mécanique des structures. *C. R. Acad. Sci. Paris*, Série II b(327), 1237-1244, 1999.

[2] F. Feyel, J-L. Chaboche. FE² multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials. *Comput. Methods Applied Mech. Engrg*, 183, 309-330, 2000.

[3] H. Ben Dhia. Problèmes mécaniques multi-échelles : la méthode Arlequin. *C. R. Acad. Sci. Paris*, Série II b(326), 899-904, 1998.

[4] H. Ben Dhia, G. Rateau. The Arlequin method as a flexible engineering design tool, *IJNME*, 2005.

[5] P. Gosselet, C. Rey. Non-overlapping domain decomposition methods in structural mechanics. *Archives of Computational Methods in Engineering*, Volume 13, Num 4, 515-572, 2006.