# Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoï contraint

Julien Yvonnet, David Ryckelynck, Philippe Lorong, Francisco Chinesta

► To cite this version:

Julien Yvonnet, David Ryckelynck, Philippe Lorong, Francisco Chinesta. Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoï contraint. Revue Européenne des Éléments Finis, 2003, 12 (4), pp.487-509. 10.3166/reef.12.487-509. hal-01513837

HAL Id: hal-01513837

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

Submitted on 25 Apr 2017

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
# Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoï contraint

## Méthode des éléments C-naturels

Julien Yvonnet — David Ryckelynck — Philippe Lorong — Francisco Chinesta

Laboratoire de Mécanique des Systèmes et des Procédés.

UMR CNRS ENSAM – ESEM, 151 boulevard de l'Hôpital, F-75013 Paris, France

www.paris.ensam.fr

RÉSUMÉ. La méthode des éléments naturels (MEN) est une nouvelle méthode dite « sans maillage » basée sur l'interpolant sibsonien pour résoudre des équations aux dérivées partielles. La méthode construit une interpolation strictement nodale ce qui simplifie l'application des conditions aux limites. Cependant des difficultés apparaissent pour des domaines non convexes : la linéarité n'est plus assurée sur les bords et des influences non désirées peuvent apparaître entre nœuds présents sur des frontières très proches l'une de l'autre dans des cas fortement concaves comme les fissures. Les solutions proposées jusqu'à présent ne permettent pas de traiter ce type de cas. Nous proposons une démarche où les fonctions de forme sont calculées en se basant sur un diagramme de Voronoï modifié, dit « contraint », permettant de traiter des problèmes de géométrie quelconque tout en garantissant les propriétés de la méthode des éléments naturels.

ABSTRACT. The natural elements method (NEM) is a new technique considered as a "meshless method" based on Sibson co-ordinates for the solution of partial differential equations. The NE shape functions are strictly interpolant which makes easy the imposition of essential boundary conditions. However, issues occur over non-convex boundaries : interpolant is not stricly linear over the whole boundary and interaction between nodes over close boundaries, like cracks, can also occur. Solutions proposed so far fail in cases like these. We propose a methodology to compute the shape functions by mean of a modified constrained Voronoï diagram. Respect of all main properties of the natural elements methods by this way without regard on the gometry of the domain is discussed.

MOTS-CLÉS : méthode des éléments C-naturels, domaines non convexes, diagramme de Voronoï contraint, critère de visibilité, conditions aux limites.

KEYWORDS: C-NEM, non-convex bodies, constrained Voronoï diagram, visibility criterion, essential boundary conditions.
# 1. Introduction

Les méthodes sans maillage (meshless ou mesh-free methods) ont été récemment introduites et suscitent un intérêt grandissant. Elles constituent une nouvelle approche permettant d'offrir une alternative à la méthode des éléments finis ou à celle des volumes finis. Les approches sans maillages évitent en particulier les fréquents et délicats remaillages en cours de simulation pour la résolution de problèmes tels que ceux rencontrés en transformations finies et, plus particulièrement, en présence de propagation de fissures ou de phénomènes présentant des interfaces mobiles. Le principal attrait de ces méthodes est de construire une fonction d'approximation basée uniquement sur une description nodale du domaine. Une structure de maillage sous-jacente d'éléments n'est donc plus nécessaire. Il est ainsi possible de résoudre un grand nombre de problèmes sans remaillage, pour des cas où l'application de la méthode des éléments finis (MEF) se révélait délicate. En effet, dans ces méthodes la qualité de l'approximation n'est pas liée à un critère géométrique sur les éléments tel que le critère d'angle minimal dans la MEF. De ce fait, l'adaptativité devient aisée, puisqu'il suffit de rajouter des nœuds dans certaines zones suivant un critère sans effectuer de remaillage. Plusieurs méthodes utilisant ce principe ont été développées. Nous citons de manière non exhaustive les méthodes suivantes : Smoothed Particle Hydrodynamics (Gingold, 1977), The Diffuse Elements Method (Nayrolles et al., 1992), Element-Free Galerkin (Belytschko, Lu et Gu (1994)), Reproductive Kernel Particle Method, (Liu, Jun et Zhang (1995)), hp-clouds (Duarte et Oden (1996)), Partition of Unity Method (Melenk et Babuška (1996)) et la Méthode des Eléments Naturels (NEM) (Sukumar, Moran et Belytschko (1998)). La Méthode des Eléments Naturels possède la particularité unique, parmi toutes les méthodes citées, de construire une fonction d'interpolation nodale. Les fonctions de forme sont des interpolants de Sibson (Sibson, 1980) qui permettent de garantir par construction : l'interpolation nodale, la propriété de partition de l'unité et la consistance linéaire. Les interpolants utilisés dans les autres méthodes de type *moindres carrés mobiles* (MLS) ou *reproductive kernel (RK)* nécessitent, quant à elles, des calculs supplémentaires coûteux pour satisfaire ces propriétés. La vérification de ces dernières permet d'imposer les conditions aux limites de façon exacte sans traitement particulier. La méthode possède cependant, dans sa version initiale (Braun et Sambridge, 1995), l'inconvénient de ne conserver ses propriétés que sur les domaines convexes. En effet, l'interpolant sibsonien se base sur la construction du diagramme de Voronoï d'un ensemble de nœuds, cet ensemble définissant toujours un domaine convexe. Lorsque le domaine est faiblement non convexe, il est facile d'introduire des variantes de la méthode originelle (Braun et Sambridge, 1995), pour lever cette difficulté. Ces approches se révèlent cependant toutes insuffisantes dans les cas fortement non convexes comme les fissures ou les problèmes d'auto-contact.

Nous proposons dans cet article une nouvelle méthode, baptisée *méthode des éléments naturels contraints* ou *méthode des éléments C-naturels* (C-NEM en anglais), adaptée aux géométries quelconques. La particularité de l'approche est de
calculer les fonctions de forme, ou interpolants sibsoniens, sur la base d'un diagramme de Voronoï modifié, appelé *diagramme de Voronoï contraint*. Ce diagramme est une construction géométrique applicable à un domaine de géométrie quelconque, à condition de connaître une description de sa frontière. Sa construction revient à introduire un critère dit de *visibilité* (Organ *et al.*, 1996) dans l'algorithme de construction du diagramme de Voronoï proposé à l'origine par Bowyer (Bowyer, 1981). Ce diagramme sert alors de structure de donnée unique et fixe pour calculer les fonctions de forme en tout point du domaine. La construction permet, en particulier, d'appliquer les méthodes d'intégration récemment introduites et adaptées à ce type de représentation (Chen *et al.*).

Dans la section 2 nous revenons sur la définition de l'interpolation naturelle, sur la manière dont elle est liée au diagramme de Voronoï et sur ses propriétés principales. La section 3 est un exposé des différents problèmes qui interviennent sur les domaines non-convexes. Notre approche est présentée dans la section 5, après avoir effectué, dans la section 4, une étude des solutions actuellement proposées pour traiter les domaines non convexes. Des illustrations numériques, données en fin de section 5, confirment la vérification du maintien des propriétés de l'interpolation au voisinage d'une forte concavité (la pointe d'une fissure).

## 2. Rappels sur l'interpolant sibsonien

Dans cette section, nous rappelons brièvement la construction des interpolants sibsoniens ainsi que leurs propriétés principales. Des détails sur ces descriptions peuvent être trouvés dans (Sukumar *et al.*, 1998). Nous montrons ensuite les problèmes induits par le caractère convexe de la construction de ces interpolants.

### 2.1. Construction

Le diagramme de Voronoï (ou de Dirichlet) est une construction géométrique unique réalisée à partir d'un ensemble de nœuds distincts dans un espace euclidien de dimension quelconque $E^K$. Pour des raisons évidentes de pratique, seuls des cas tels que $K \le 3$ sont ici considérés. Soit un ensemble $n_1, \dots, n_N$ de nœuds distincts dans $E^K$, on peut définir des zones $T_I$ telles que :

$$T_I = \{ \mathbf{x} \in E^K : d(\mathbf{x}, n_I) < d(\mathbf{x}, n_J), \forall J \neq I \}$$

où d est une distance euclidienne et $T_I$ est la cellule de Voronoï associée au nœud $n_I$. Elle est l'intersection de plusieurs demi-espaces, chacun limités par un plan médiateur d'un segment $n_I n_J$. De ce fait les cellules sont convexes et polyédriques si le nœud est à l'intérieur du domaine convexe formé par l'ensemble des nœuds, ou ouverte si le nœud est sur la frontière $\Gamma$ de ce domaine, comme le montre la figure 1.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** *Diagramme de Voronoï d'un ensemble de nœuds du plan.*</figcaption></figure>

La cellule de Voronoï $T_I$ est une région de l'espace telle que tout point $\mathbf{x}$ à l'intérieur de cette cellule est plus proche du nœud $\mathbf{n}_I$ que de n'importe quel autre nœud $\mathbf{n}_J$. Deux nœuds partageant une face de leurs cellules de Voronoï respectives sont appelés *voisins naturels*.

Pour construire les interpolants sibsoniens, il est nécessaire d'introduire les cellules de Voronoï du second ordre. On peut définir une cellule du second ordre $T_{IJ}$ de la façon suivante :

$$ T_{IJ} = \{ \mathbf{x} \in E^K : d(\mathbf{x}, \mathbf{n}_I) < d(\mathbf{x}, \mathbf{n}_J) < d(\mathbf{x}, \mathbf{n}_K), \forall K \neq I, \forall K \neq J \} $$

Un point quelconque se trouvant dans $T_{IJ}$ a comme voisin le plus proche $\mathbf{n}_I$ et comme second voisin le plus proche $\mathbf{n}_J$. Introduisons un point $\mathbf{x}$ en lequel on désire calculer les fonctions de forme naturelles. Sur la figure 2, la cellule délimitée par le contour *abcde* est la cellule du second ordre associée au point $\mathbf{x}$ et au nœud $\mathbf{n}_I$.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** *Construction des interpolants sibsoniens.*</figcaption></figure>
Définissons $\kappa(\mathbf{x})$ une mesure de Lebesgue (longueur, aire et volume en 1D, 2D et 3D respectivement). Soit $\kappa_I$ la mesure de Lebesgue de la cellule $T_I$ et $\kappa_{x,J}$ la mesure de Lebesgue de la cellule $T_{x,J}$. L'interpolant sibsonien calculé au point $\mathbf{x}$ associé au voisin $\mathbf{n}_I$ est le rapport entre la mesure de Lebesgue $\kappa_x$ de la cellule de Voronoï du permier ordre associée à $\mathbf{x}$ et la mesure de Lebesgue $\kappa_{xI}$ de la cellule du second ordre associée à $\mathbf{x}$ et $\mathbf{n}_I$.

$$ \phi_I(\mathbf{x}) = \frac{\kappa_{xI}}{\kappa_x} $$

Si l'on considère le cas bidimensionnel de la figure 2, la fonction de forme associée au nœud $\mathbf{n}_I$ calculée en $\mathbf{x}$ est donnée par :

$$ \phi_I(\mathbf{x}) = \frac{\text{Aire}(abcde)}{\text{Aire}(af gh\theta)} $$

Le schéma d'interpolation pour un champ $\mathbf{u}$ à valeurs dans un espace $\mathbb{R}^\alpha$ quelconque, $\mathbf{u}(\mathbf{x}) : \Omega \in E^K \to \mathbb{R}^\alpha$, est de la forme :

$$ \mathbf{u}^{\mathbf{h}}(\mathbf{x}) = \sum_{I=1}^{n} \phi_I(\mathbf{x}) \mathbf{u}_I $$

Où $\mathbf{u}_I$ est la valeur nodale du champ $\mathbf{u}$ en $\mathbf{n}_I$ et $\phi_I(\mathbf{x})$ la fonction de forme associée à ce même nœud. On peut définir le domaine d'influence d'une fonction de forme associée à un nœud de la manière suivante : il s'agit de l'union des cercles (sphères en 3D) vides circonscrits à des triangles (tétraèdres en 3D) dont un des sommets est $\mathbf{n}_I$, comme représentée sur la figure 3.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Support de la fonction de forme associée au nœud $\mathbf{n}_I$.</figcaption></figure>

## 2.2. Propriétés

Nous présentons ici brièvement les propriétés fondamentales des fonctions de forme de la méthode des éléments naturels.
### 2.2.1. Interpolation nodale

Cette propriété peut être exprimée par :

$$ \phi_I(\mathbf{x}_J) = \delta_{IJ} $$

qui implique que comme dans la méthode des éléments finis, une fonction $\phi_I(\mathbf{x})$ vaut 1 au nœud $\mathbf{n}_I$ et zéro en tout nœud $\mathbf{n}_J$ tel que $I \neq J$. Cela entraîne que l'interpolation construite passe par les valeurs nodales et permet l'application directe des conditions aux limites de Dirichlet pour les nœuds concernés. Cette propriété que l'on peut trouver dans la littérature sous le nom de *propriété du delta de Kronecker*, fait défaut dans certaines autres méthodes meshless (SPH, la Méthode des Moindres Carrés Mobiles, ...). De ce fait, imposer les conditions aux limites demande des efforts de calculs supplémentaires. Plusieurs méthodes ont été proposées dans ce but : citons pour information la méthode du multiplicateur de Lagrange (Belytschko et al., 1994), une méthode de pénalité (Zhu et al., 1998), la méthode de transformation (Chen et al., 1996), l'approche basée sur le principe de D'Alembert (Günther et al., 1997), l'introduction de fonction poids singulière (Kaljevic et al., 1997) et la méthode *boundary singular kernel* (Chen et al., 1998). Nous évoquons plus loin une approche couplant méthodes meshless et Elements Finis (Sukumar et al., 2001).

### 2.2.2. Partition de l'unité

Les interpolants sibsonien, qui sont les fonctions de forme utilisées par la Méthode des Eléments Naturels, possèdent la propriété suivante :

$$ \sum_{I=1}^{n} \phi_{I}(\mathbf{x}) = 1 \quad \forall \mathbf{x} \in \Omega $$

Cette propriété peut être assimilée à une consistance linéaire d'ordre zéro ce qui signifie que l'interpolant construit reproduit de manière exacte les fonctions constantes.

### 2.2.3. Consistance linéaire

Toute interpolation permettant de reproduire exactement une fonction variant linéairement vis à vis des variables d'espace possède la propriété appelée consistance linéaire. Sibson (Sibson, 1980) a montré que l'interpolant sibsonien vérifie cette propriété. Celle-ci peut s'écrire de la façon suivante :

$$ \mathbf{x} = \sum_{I=1}^{n} \phi_{I}(\mathbf{x}) \mathbf{x}_{I} \quad \forall \mathbf{x} \in \Omega $$
De nombreuses méthodes sans maillage (SPH, ...) ne vérifient pas, sans introduction de traitements spécifiques, la consistance linéaire. Il est à noter que la méthode proposée étant de type partition de l'unité, il est possible d'enrichir la base d'approximation pour obtenir une consistance plus que linéaire, voire sur une autre base de fonctions que les polynômes.

### 2.2.4. Linéarité stricte de l'approximation sur les bords

Sur la frontière d'un domaine convexe, les fonctions de forme sont strictement linéaires entre deux nœuds voisins. La fonction d'approximation étant une combinaison linéaire des fonctions de forme, celle-ci est également linéaire sur la frontière de ce domaine. La démonstration peut être trouvée dans (Sukumar, 1998).

## 3. Problématique associée aux domaines non convexes

La méthode des éléments naturels offre donc des propriétés séduisantes, tant que les domaines sont convexes. Malheureusement, en calcul des structures, les domaines étudiés présentent le plus souvent des géométries non convexes. Nous présentons dans cette section les problèmes intervenant pour le cas des géométries non convexes.

### 3.1. Non-linéarité de l'approximation sur les bords non convexes

La propriété énoncée en 2.2.4. est due au fait que sur les domaines de géométrie convexe, les cellules de Voronoï associées aux nœuds de la frontière ont, dans le cas 2D, des aires infinies (en 3D des volumes infinis) comme on peut le voir sur la figure 4. Les fonctions de forme étant calculées comme des rapports d'aires (Sukumar, 1998), la contribution des nœuds intérieurs, qui elle est finie, devient négligeable devant celle des nœuds sur la frontière. Pour un domaine non convexe, il est aisé de vérifier, comme le montre la figure 4, que les cellules de Voronoï associées aux nœuds le long de la frontière où le domaine n'est pas convexe, $\Gamma_1$ sur la figure 4, ont des aires finies. La démonstration faite pour les cas convexes reste vérifiée pour $\bar{\Gamma} = \Gamma_2 \cup \Gamma_3 \cup \Gamma_4$ mais est fausse pour $\Gamma_1$: pour un point x localisé sur $\Gamma_1$, il existe des contributions non nulles de la part de nœuds internes. La fonction d'interpolation $u^h(x)$ n'est donc pas linéaire entre 2 nœuds adjacents de $\Gamma_1$.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Domaine non convexe. L'approximation **u**<sup>h</sup>(x) est linéaire sur $\tilde{\Gamma} = \Gamma_2 \cup \Gamma_3 \cup \Gamma_4$ mais pas sur $\Gamma_1$ (Sukumar, 1998).</figcaption></figure>

**3.2. Chevauchement du support des fonctions de forme au travers des fissures**

Dans les cas fortement non convexes tels que les fissures, des influences
parasites peuvent intervenir de manière très importante et introduire des erreurs qui
sont du même ordre de grandeur que la solution recherchée. En effet, tout point à
l'intérieur du domaine d'influence d'un nœud comme décrit figure 3 subit
l'influence de celui-ci. Ainsi un point x (voir figure 5) peut se trouver dans le
domaine d'influence d'un nœud situé de l'autre côté de la fissure. Ces problèmes
n'existent pas dans la méthode des éléments finis, car un nœud ne subit l'influence
que des nœuds auxquels il est connecté, il ne peut donc subir l'influence de nœuds
situés de l'autre côté d'une fissure. Ces défauts sont présentés en détail dans
(Sukumar *et al.*, 2001) et (Cueto *et al.*, 2000).

<figure><img src="image_6.png" /><figcaption>**Figure 5.** *Modélisation d'une fissure* (Sukumar, 1998).</figcaption></figure>
### 3.3. Domaines d'intégration

Dans les méthodes dites de Galerkin l'intégration numérique est habituellement utilisée pour évaluer les intégrales de la formulation variationnelle dite faible. Dans la méthode des éléments finis, la méthode de Gauss est appliquée sur chaque élément. Dans ce cas, le support des fonctions de forme coïncide avec le domaine utilisé pour l'intégration. Deux difficultés apparaissent dans le cadre des méthodes meshless : premièrement, les fonctions de forme sont généralement rationnelles (problème également rencontré pour de nombreux éléments finis), la méthode de Gauss se révèle alors non optimale. Deuxièmement, on ne peut définir les domaines d'intégration comme dans la MEF, puisqu'il n'y a pas de maillage. Dolbow et Belytscko ont montré (Dolbow et al., 1998) que si le domaine d'intégration était choisi sans lien avec le support des fonctions de forme, des erreurs considérables pouvaient être commises. Il ne faut donc pas utiliser la triangulation de Delaunay comme base de l'intégration. (Chen et al., 2001) montrent que l'intégration doit être effectuée pour chaque nœud sur un domaine représentatif entourant le nœud. La cellule de Voronoï est à ce titre un bon choix de domaine représentatif. Pour tenir compte des limites du domaine matériel, le domaine d'intégration est défini comme l'intersection de la cellule de Voronoï et de la frontière du domaine (voir figure 6).

<figure><img src="image_7.png" /><figcaption>**Figure 6. Diagramme de Voronoï et domaines d'intégration.**</figcaption></figure>

Une difficulté survient dans les cas fortement non convexes, comme le voisinage d'une fissure, lorsque le diagramme de Voronoï est construit sans tenir compte de la frontière. La simple intersection des cellules avec la frontière ne définit plus alors, comme l'illustre la figure 7, des volumes élémentaires recouvrant l'ensemble du domaine étudié. La méthode que nous proposons dans la section 5.2. permet de résoudre cette difficulté de manière simple.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Exemple de non recouvrement du domaine pour une construction des domaines d'intégration basée sur l'intersection directe du diagramme de Voronoï avec la frontière.</figcaption></figure>

# 4. Traitement des domaines non convexes : approches actuelles

Dans ce paragraphe nous ne traitons pas de l'approche consistant à décomposer le domaine non convexe en sous-domaines convexes car cette approche n'est applicable que sur des domaines possédant une ou plusieurs symétries (Sukumar *et al.*, 1998).

## 4.1. Interpolants non sibsoniens et couplages EF-EN

(Sukumar *et al.*, 2001) proposent d'utiliser des interpolants différents, non sibsoniens, basés sur des rapports de longueurs plutôt que sur des rapports d'aires. Les propriétés de cet interpolant induisent la linéarité des fonctions de forme sur la frontière que le domaine soit convexe ou non. Ce type d'interpolant a cependant l'inconvénient de présenter une dérivée discontinue sur les cercles de Delaunay. De plus, cet interpolant n'empêche pas les influences non désirées au travers de fissures, c'est-à-dire pour les cas fortement non convexes.

Les mêmes auteurs proposent une seconde approche pour garantir la linéarité sur la frontière afin de permettre une vérification aisée des conditions aux limites. Cette
approche consiste à placer une couche d'éléments finis sur les bords. En transformations finies il est alors nécessaire d'introduire un critère géométrique sur la déformation des éléments, et donc également des procédures de remaillage, ce qui est en contradiction avec les raisons qui poussent à utiliser une méthode sans maillage.

## 4.2. Les formes alpha

(Cueto et al., 2000) adoptent un point de vue tout à fait original et proposent de ne pas changer l'interpolation mais de restreindre les voisins naturels intervenant dans le calcul de l'approximation. En effet, la liste des voisins naturels d'un nœud définit de façon immédiate la zone d'influence de la fonction de forme associée à ce nœud. Or, pour les domaines fortement non convexes, les difficultés proviennent justement du fait de l'étendue indésirable de certains domaines d'influence.

Dans cette optique, la proposition est la suivante : éliminer de la liste des voisins naturels tous les nœuds qui appartiennent à un tétraèdre dont le rayon de la sphère circonscrite dépasse une certaine valeur fixée $\alpha$. Cette méthode permet effectivement de garantir la linéarité le long du bord d'un domaine non convexe et évite les influences de nœuds non désirées.

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Restriction des voisins naturels par la méthode des formes alpha. Défaut de la méthode pour les cas fortement non convexes.</figcaption></figure>

Cueto et ses collaborateurs proposent deux conditions permettant de définir la valeur $\alpha$ maximale : une condition basée sur la distance entre les bords d'une fissure et leur axe médian (Cueto et al., 2001) ainsi qu'une condition basée sur le rayon de courbure local de la frontière. La distance entre deux nœuds consécutifs sur le bord doit être inférieure aux deux paramètres cités. La densité nodale sur le bord leur est inversement proportionnelle. L'application de ces conditions puis la suppression des voisins suivant la méthode des formes $\alpha$ garantissent, sur la frontière, la linéarité de l'interpolation. Cependant, lorsque la distance entre deux surfaces tend vers zéro (auto-contact) ou lorsqu'une surface présente un angle vif ou une pointe de fissure et
que le rayon de courbure local s'annule, la vérification des conditions précédentes devient impossible, la densité locale théorique devant alors être infinie.

### 4.3. Le critère de visibilité (Organ et al., 1996)

Les approches tentant de réaliser une interpolation naturelle à partir de la seule description nodale se sont révélées infructueuses dans le cas de domaines trop fortement non convexes. L'introduction de la connaissance de la frontière en plus du nuage de points semble nécessaire.

Plusieurs démarches visant à modifier le support des fonctions de forme, grâce à la prise en compte de la description de la frontière, ont été développées pour des méthodes sans maillage, voir par exemple les approches *Element Free Galerkin* (EFG) ou *Smoothed Particles Hydrodynamics* (SPH). Ces démarches sont basées sur l'utilisation de critères de visibilité visant à limiter, en un point donné, le nombre de nœuds influents (ou encore nœuds visibles). La liste des nœuds influents en un point étant fonction de l'approche sans maillage utilisée, il est nécessaire d'adapter le traitement en fonction de l'approche concernée.

Krysl et Belytschko ont proposé une façon de définir le domaine d'influence à proximité de fissures et autres domaines non convexes dans le cadre de la méthode EFG. Ce critère est décrit en détails dans (Organ et al., 1996) et testé dans la simulation de propagation de fissures (Fleming et al., 1997), (Belytschko et al., 1997). Dans l'approche de Krysl et Belytschko, un point **x** est dans le domaine d'influence d'un nœud $n_I$ s'il remplit deux conditions. Il doit premièrement être contenu dans la zone géométrique définissant le domaine d'influence d'un nœud en général (par exemple un cercle de rayon donné dans le cas des interpolants de type moindres carrés mobiles). Deuxièmement, le nœud $n_I$ doit être visible pour un observateur placé sur le point **x** en considérant les frontières opaques.

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

<figure><img src="image_11.png" /></figure>
**Figure 9.** *Support de fonction de forme de type moindres carrées mobiles à proximité d'une fissure modifié par le critère de visibilité — Fonction de forme associée au nœud n<sub>1</sub> (Belytschko et al., 1996).*

Cette approche introduit deux types de discontinuités dans les fonctions de forme : une discontinuité attendue en AC (figure 9), qui correspond à la fissure et une discontinuité numérique, en AB. Cette discontinuité introduit des erreurs (Organ et al., 1996). D'autres critères ont alors été élaborés pour supprimer cette discontinuité numérique : les critères de diffraction et de transparence (Organ et al., 1996).

# 5. Méthode des éléments naturels contraints ou méthode des éléments C-naturels (C-NEM)

Nous proposons dans cette section une modification de la méthode des éléments naturels permettant de réaliser une interpolation sur des domaines non convexes ou quelconques et évitant tous les problèmes énoncés dans la section 3. Nous introduisons la notion de diagramme de Voronoï contraint qui servira de support au calcul des fonctions de forme. La construction du diagramme de Voronoï contraint nécessite l'application d'un critère de visibilité analogue à celui présenté dans le paragraphe précédent.

## 5.1. Diagramme de Voronï contraint

Pour construire le diagramme de Voronoï contraint nous introduisons le critère de visibilité défini dans le paragraphe 5.1.2 ci-dessous. Le terme diagramme de Voronoï *borné* ou *contraint* peut être trouvé dans (Klein, 1996). Le terme *contraint* a été choisi par analogie avec la triangulation de Delaunay contrainte (Borgers), car comme nous le verrons, le diagramme de Voronoï contraint est le dual de la triangulation de Delaunay contrainte. C'est ce terme que nous choisirons par la suite. La plupart des définitions données dans les références citées sont restreintes au cas 2D et utilisent pour décrire la frontière la notion de *Planar Straight Graph* qui se réfère à un ensemble de nœuds et de segments définissant le contour fermé d'un polygone. La définition du diagramme de Voronoï contraint que nous proposons est par contre valable pour une dimension $K \le 3$.

### 5.1.1. Définitions préalables

Nous notons $\Gamma$ la frontière du domaine $\Omega$ étudié. $\Omega$ est un ouvert de $\mathbb{R}^K$. $\Gamma$ est défini par la liste de segments en 2D et de faces triangulaires en 3D qui la composent.
Nous notons $I(\mathbf{n})$ l'ensemble des nœuds intérieurs à $\Omega$ et $E(\mathbf{n})$ l'ensemble des nœuds appartenant à $\Gamma$.

### 5.1.2. Définition du critère de visibilité

Un point $\mathbf{x}$ est visible d'un nœud $\mathbf{n}_I$ si et seulement si il vérifie l'une des trois conditions suivantes :

- $\mathbf{n}_I \in I(\mathbf{n})$ et il n'existe aucune intersection entre le segment $x\mathbf{n}_I$ et $\Gamma$

- $\mathbf{x} \in \Omega, \mathbf{n}_I \in E(\mathbf{n})$ et il n'existe aucune intersection entre le segment $x\mathbf{n}_I$ et $\Gamma$

- $\mathbf{x}$ est situé sur une arête ou une face contenant $\mathbf{n}_I, \mathbf{n}_I \in E(\mathbf{n})$

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Critère de visibilité appliqué aux nœuds : nœuds visibles (□) depuis $\mathbf{n}_I$. Cas où $\mathbf{n}_I \in \Omega$ — Cas où $\mathbf{n}_I$ est sur $\Gamma$.</figcaption></figure>

### 5.1.3. Définition du diagramme de Voronoï contraint

Le diagramme de Voronoï contraint est défini par un ensemble de cellules que nous appelons cellules de voronoï contraintes. A chaque nœud $\mathbf{n}_I$ est associée la cellule de voronoï contrainte $T_I^C$ définie par :

$$T_I^C = \{ \mathbf{x} \in \Omega : d(\mathbf{x}, \mathbf{n}_I) < d(\mathbf{x}, \mathbf{n}_J), \forall J \neq I, \mathbf{n}_J \text{ visible depuis } \mathbf{n}_I \}$$

Autrement dit tout point $\mathbf{x}$ contenu dans $\Omega$ à l'intérieur de $T_I^C$ est plus proche du nœud $\mathbf{n}_I$ que de n'importe quel nœud $\mathbf{n}_J$ visible depuis le nœud $\mathbf{n}_I$.

Les cellules associées aux nœuds de $I(\mathbf{n})$ sont convexes et polyédriques. Les cellules associées aux nœuds de $E(\mathbf{n})$ ne sont pas nécessairement convexes. Il est à noter que sur la partie droite de la figure 11 nous avons prolongé les arêtes de voronoï traversant la frontière $\Gamma$.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** A gauche ensemble de nœuds P(n). Au centre diagramme de Voronoï de P(n). A droite diagramme de Voronoï de l'ensemble P(n) contraint par la frontière imposée Γ.</figcaption></figure>

## 5.2. Interpolation C-NEM

Dans la méthode des éléments naturels contraints ou méthode C-NEM, la fonction d'interpolation $\mathbf{u}^h(\mathbf{x})$ est calculée par :

$$ \mathbf{u}^h(\mathbf{x}) = \sum_{I=1}^{V} \phi_I^C(\mathbf{x}) \mathbf{u}_I $$

où : V est le nombre de voisins naturels visibles depuis le point $\mathbf{x}$, $\phi_I^C$ sont les interpolants sibsoniens identiques à ceux utilisés dans la méthode des éléments naturels mais calculés sur la base du diagramme de Voronoï contraint. Les valeurs $\mathbf{u}_I$ sont les inconnues nodales des voisins naturels visibles depuis le point $\mathbf{x}$.

Notons que lorsque $\mathbf{x}$ s'approche d'une arête ou d'une face de $\Gamma$, le contour des cellules de voronoï concernées est prolongé à l'infini à l'extérieur de $\Omega$ ce qui permet de retrouver des contributions en aire (pour le 2D) ou en volume (pour le 3D) infinies. Un exemple de tracé du prolongement des arêtes des cellules de voronoï contraintes est visible sur la partie droite de la figure 11 dans le cas 2D (prolongement des arêtes des cellules). Ces contributions infinies permettent de retrouver la linéarité de l'interpolation sur la frontière du domaine. Il est à noter que dans la pratique, le fait d'avoir de aires infinies ne pose pas de problèmes d'un point de vue numérique car sur les bords on utilise une approximation linéaire justifiée par le fait que l'interpolant sibsonien, restreint par le critère de visibilité, devient linéaire sur n'importe quel arête du bord.

Dans cette méthode, les nœuds qui interviennent pour le calcul des fonctions de forme et qui apportent une contribution non nulle sont d'abord sélectionnés comme étant des voisins naturels du point $\mathbf{x}$, puis subissent une deuxième sélection ne maintenant que ceux qui sont visibles depuis le point $\mathbf{x}$. La figure 12 montre un exemple de calcul de fonction de forme par la méthode C-NEM. Le diagramme de
Voronoï contraint est représenté. Les prolongements externes au domaine, et infinis,
des arêtes de Voronoï n'ont pas été représentés pour des raisons de lisibilité. On peut
noter que seuls les voisins naturels vus depuis le point **x** sont pris en compte pour le
calcul des fonctions de forme. Les nœuds présents sur la frontière opposée de la
fissure n'interviennent pas dans le calcul. Le diagramme de Voronoï contraint
fournit ainsi une structure de données fixe et unique pour le calcul des fonctions de
forme en tout point à l'intérieur du domaine. Il n'est pas nécessaire de traiter le
problème par sous-domaines convexes et on peut calculer les fonctions de forme de
manière continue en se déplaçant à l'intérieur du domaine. Enfin le diagramme de
Voronoï contraint permet de représenter correctement les domaines d'intégration $T_I^C$
décrits dans le paragraphe 3.3. (Comparer les figures 7 et 12)

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Calcul de fonction de forme au voisinage d'une concavité par la méthode C-NEM. A gauche introduction d'un point dans le diagramme de Voronoï contraint. A droite : les nœuds de la frontière opposée ne sont pas retenus par le critère de visibilité et sont ignorés pour le calcul des fonctions de forme.</figcaption></figure>

**5.3. Vérification des propriétés aux abords d'une forte concavité**

Nous montrons ici que la méthode C-NEM vérifie, contrairement à la méthode
des éléments naturels, les principales propriétés énoncées dans la section 2.2.

Ces propriétés sont : la linéarité sur les bords, la Partition de l'Unité et la
consistance linéaire. Nous montrons également que l'interpolation construite est
continue près de la pointe d'une fissure.

5.3.1. *Continuité de l'interpolation*

Sur la figure 13 est tracé le diagramme de Voronoï contraint associé à un nuage
de point décrivant un domaine présentant une fissure. Il est à noter que la frontière $\Gamma$
est nécessaire à la description du domaine, sans cette frontière la fissure ne serait pas
identifiable.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Discrétisation d'un domaine non convexe ; à gauche : maillage et frontière. A droite : diagramme de Voronoï contraint du domaine.</figcaption></figure>

Nous avons vu en 4.3. que, dans le cadre des interpolants de type moindres carrés mobiles avec critère de visibilité, les fonctions de forme étaient discontinues aux alentours d'une fissure. Ceci est dû au fait que le support d'une fonction de forme est de géométrie fixée et qu'il ne dépend pas des nœuds du voisinage. L'application du critère de visibilité conduit, dans ce cas, à une discontinuité comme le montre la traversée du segment AB sur la figure 9.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Support des fonctions de forme naturelles contraintes $\phi_I^C$ aux abords d'une fissure.</figcaption></figure>

Dans le cas de la méthode C-NEM, chaque frontière du support des fonctions de forme associées aux nœuds proches de la concavité passe par la pointe de celle-ci (la forte non convexité du domaine étant décrite par sa frontière, cette dernière possède nécessairement un nœud en pointe de fissure). Chaque support de fonction de forme n'est ainsi tronqué que par la frontière physique du domaine. Les fonctions de forme restent donc continues sur le domaine, et l'interpolation qui en découle également (c'est une combinaison linéaire des fonctions de forme). De cette continuité « naturelle » découle qu'aucun lissage des fonctions de forme, en adoptant des critères tels que la transparence ou la diffraction qui étaient délicats à programmer surtout en 3D, n'est nécessaire.
La figure 15 représente le voisinage de la pointe de la fissure. Nous avons calculé les fonctions de forme des nœuds du voisinage sur le parcours circulaire autour de la pointe. On vérifie que les fonctions de forme sont continues tout le long du parcours. De plus, elles restent infiniment dérivables sauf si l'on passe par les nœuds, lieux où les fonctions de forme naturelles ne sont pas dérivables. Nous avons également vérifié sur ce parcours la partition de l'unité à la précision de la machine (10⁻¹⁴), cette propriété est examinée dans le paragraphe suivant.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** *Parcours circulaire autour de la pointe de la fissure.*</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** *Fonctions de forme associées aux nœuds du voisinage le long du parcours circulaire autour de la pointe de la fissure.*</figcaption></figure>
### 5.3.2. Consistance linéaire

En un point donné du domaine, le fait de réduire l'ensemble des nœuds voisins naturels aux seuls nœuds visibles, ne change pas la nature de l'interpolation. Celle-ci est toujours calculée par l'interpolant sibsonien sur un ensemble donné de nœuds. Tout se passe comme si les nœuds non visibles n'existaient pas. Comme le point est situé à l'intérieur du domaine, l'ensemble des nœuds visibles définit un domaine convexe contenant le point. Sibson (Sibson, 1980) démontre que, dans le cas des interpolants naturels, les propriétés de Partition de l'unité :

$$ \sum_{I=1}^{n} \phi_{I}(\mathbf{x}) = 1 $$

et consistance linéaire :

$$ \mathbf{x} = \sum_{I=1}^{n} \phi_{I}(\mathbf{x}) \mathbf{x}_{I} $$

sont vérifiées pour peu que le point où sont calculées les fonctions de forme soit à l'intérieur du domaine convexe formé par ses voisins naturels. Dans l'approche C-NEM, seule la liste des voisins ayant changée, les deux propriétés restent, en chaque points du domaine, vérifiées par la nouvelle interpolation.

La figure 17 se focalise sur le voisinage de la pointe d'une fissure. Nous avons dans cet exemple appliqué un champ scalaire linéaire du type :

$$ T(y) = \alpha y $$

Nous avons ensuite calculé la fonction d'interpolation $T^h$ telle que :

$$ T^h(\mathbf{x}) = \sum_{I=1}^{V} \phi_I^C(\mathbf{x}_I) T_I $$

le long du parcours décrit sur la figure 17. La figure 18 donne les résultats numériques s'y référant. Sur les trajets AB, CD et EF, nous vérifions bien que la fonction augmente linéairement comme la solution exacte $T(y)$. Sur les trajets BC et DE, $T^h$ est bien constante, comme la solution exacte $T(y)$. Sur cet exemple, la consistance linéaire est vérifiée numériquement à la précision de calcul près ($10^{-14}$). A noter que la discontinuité de la dérivée de $T^h$ n'est due qu'au changement de direction sur le parcours.
<figure><img src="image_19.png" /><figcaption>**Figure 17.** *Parcours carré autour de la fissure.*</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** *Fonction scalaire T<sup>h</sup> calculée le long du parcours carré ABCDEF représenté figure 17.*</figcaption></figure>

### 5.3.3. Linéarité de l'approximation sur les bords

Le test suivant permet de confirmer numériquement la linéarité de l'approximation sur la frontière du domaine par la méthode des éléments naturels contraints. Sur la figure 19, on applique comme en 5.3.2. un champ scalaire linéaire du type :

$$T(y) = \alpha y$$

Nous appliquons cette fois des conditions aux limites en imposant les valeurs T1 et T2 sur les nœuds du bord **n**₁ et **n**₂ respectivement. Le point où l'on évalue T se déplace suivant le trajet ABCD (voir figure 19). Il faut vérifier que lorsque l'on calcule les fonctions de forme en un point de CD, celui-ci ne subit que l'influence des nœuds **n**₁ et **n**₂. Les courbes de la figure 20 illustrent deux cas : pour le premier,
les valeurs nodales imposées sont T1 = 0 et T2 = 0. On constate que sur le trajet CD, T^h calculée est nulle. Dans l'exemple 2, T1 = 50 et T2 = 0. On vérifie que sur le trajet CD, la température T^h est bien linéaire. Cette linéarité a également été vérifiée à la précision de calcul près (10⁻¹⁴).

<figure><img src="image_21.png" /><figcaption>Figure 19. Parcours ABCD.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 20. Vérification numérique de la linéarité sur le bord du domaine non convexe (entre C et D)</figcaption></figure>
# 6. Conclusion

La méthode des éléments naturels contraints ou méthode des éléments C-naturels a été présentée. Celle-ci permet de construire une interpolation sur des domaines quelconques en vérifiant toutes les propriétés de la méthode des éléments naturels : continuité, Partition de l'Unité, consistance linéaire et linéarité sur les bords, même au voisinage de fortes concavités comme les fissures. La méthode décrite se base sur le diagramme de Voronoï contraint qui permet d'éviter de décomposer le domaine en sous-domaines convexes et d'offrir une structure de données unique permettant de calculer les fonctions de forme en tout point d'un domaine de géométrie quelconque. Une description réaliste des domaines d'intégration est de plus obtenue directement. La méthode permet donc de traiter des problèmes sur des domaines de géométries quelconques et de manière efficace et promet d'être un outil essentiel pour l'analyse des phénomènes en transformations finies où les géométries peuvent évoluer de manières complexes et non convexes. Par ailleurs, l'interpolation introduite est, comme pour l'approche NEM, nodale ce qui permet d'appliquer des conditions aux limites par simple substitution des inconnues dans le système linéaire résultant de la discrétisation. L'extension de cette méthode au cas 3D est l'un des travaux de recherche actuel des auteurs.

# 12. Bibliographie

Belytschko T., Kronggauz Y., Fleming M., Organ D., Liu W. K., «Smoothing and accelerated Computations in the Element free Galerkin Method», 1995.

Belytschko T., Kronggauz Y., Organ D., Fleming M., « Meshless Methods : An overview and Recent Developments », Comput. Meth. Appl. Mech. Engng., vol. 139, 1996, p. 3-47.

Belytschko T., Liu W. K., Singer M.« On adaptativity and error Criteria for meshfree methods », New Advances in Adaptive Computational Methods in Mechanics, edited by P. Ladeveze and J. T. Oden.

Belytschko T., Lu Y. Y., Gu L., « Element-Free Galerkin Methods », Int. J. Numer. Meth. Eng., vol. 37, 1994, p. 229-256.

Borgers B., « Generalized Delaunay triangulations of non-convex domains», Comput. Math. Appl., vol. 20, n°7, p. 45-49.

Bowyer A., « Computing the Dirichlet tessellation», The Computer Journal, vol. 24, n°2, 1981, p. 162-166.

Braun J., Sambridge S., « A numerical method for solving partial aquations on highly irregular grids», Nature, vol. 376, 1995, p. 655-660.

Chen J. S., Pan C., Wu C. T. Liu W. K., « Reproductive Kernel Particle Methods for large deformation analysis of non linear structures », Comput. Meth. Appl. Mech. Engng., vol. 139, 1996, p. 195-129.

Chen J.S., Wang H. P., « New boundary conditions treatments un meshfree computation of contact problems », Comput. Meth. Appl. Mech. Engng, 1998.
Chen J-S., Wu C-T., Yoon S., You Y., « A stabilized conforming nodal integration for Galerkin mesh-free methods », *Int. J. Numer. Meth. Engng*, vol. 50, 2001, p. 435-466.

Cueto E., Doblaré M., Gracia L., « Imposing essential boundary conditions in the natural elements method by means of density-scaled alpha-shapes », *Int. J. Numer. Meth. Engng*, vol. 49, 2000, p. 519-546.

Cueto E., Cegoñino J., Calvo B., Doblaré M., « On the imposition of essential boundary conditions in Natural Neighbour Galerkin Methods », *Accepté en communication dans Int. J. Numer. Meth. Engng*, 2002.

Dolbow J., Belytschko T., «Numerical integration of the Galerkin weak form in meshfree methods», *Computational Mechanics*, vol. 23, 1999, p.219-230.

Fleming M., Chu Y. A., Moran B., Belytschko T., «Enriched element-free Galerkin Methods for crack tip-fields», *Int. J. Numer. Meth. Engng*, vol. 40, 1997, p. 1483-1504.

Günther F. C., Liu W. K., « Implementation of boundary conditions for meshless methods », *Comput. Meth. Appl. Mech. Engng.*, submitted, 1997.

Kaljevic I., Saigal S., « An improved Element Free Galerkin Formulation », *Int. J. Numer. Meth. Engng*, vol. 40, 1997, p.2953-2974.

Klein R., Lingas A., « A linear-time randomized algorithm for the bounded Voronoï diagram of a simple polygon », *International Journal of Computational Geometry and Applications*, vol 6, 1996, p. 263 - 278.

Nayrolles B., Touzot G., Villon P., « Generalizing the finite element method : Diffuse approximation and diffuse elements », *Computational Mechanics*, vol. 10, issue 5, 1992, p.307-318.

Organ D., Fleming M. Terry T. Belytschko T., « Continuous meshless approximations for Nonconvex bodies by Diffraction and Transparency », *New Advances in Adaptive Computational Mechanics*, vol. 19, 1996, p. 1-11.

Sambridge M., Braun J., Mc Queen M., « Geophysical parameterization and interpolation of irregular data using natural neighbors », *Geophysi. J. Int.*, vol. 122, 1995, p. 837-857.

Sibson R., « A vector Identity for the Dirichlet tesselations », *Math. Proc. Camb. Phil. Soc.*, vol. 87, 1980, p. 151-155.

Sukumar N., Moran B., Belytschko T., « The Natural Elements Method in solid mechanics », *Int. J. Numer. Meth. Engng*, vol. 43, 1998, p. 839-887.

Sukumar N., Moran B., Semenov A. Y., Belikov V. V., « Natural neighbor Galerkin methods », *Int. J. Numer. Meth. Engng*, vol. 50, 2001, p. 207-219.

Watson D. F., « Computing the n-dimensional Delaunay tessellation with application to Voronoï polytopes », *The Computer Journal*, vol. 24, n°2, 1981, p. 167-172.

Zhu T., Atluri S. N., « A modified collocation method and a penalty formulation for enforcing the essential boundary conditions in the Element-free Galerkin Method », *Computational Mechanics*, vol. 21, 1998, p. 211-222.