# Aspects numériques de la comparaison éléments finis-mesures de champs Application à l'identification de propriétés de matériaux

Gustavo Silva, Rodolphe Le Riche, Jérôme Molimard, Alain Vautrin

► To cite this version:

Gustavo Silva, Rodolphe Le Riche, Jérôme Molimard, Alain Vautrin. Aspects numériques de la comparaison éléments finis-mesures de champs Application à l'identification de propriétés de matériaux. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01509272

HAL Id: hal-01509272

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

Submitted on 16 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 CC0 1.0 - Universal - International License
# Aspects numériques de la comparaison éléments finis-mesures de champs

## Application à l'identification de propriétés de matériaux

G.H.C. Silva, R. Le Riche, J. Molimard et A. Vautrin

École Nationale Supérieure des Mines de Saint-Étienne, France
158, cours Fauriel F-42023 Saint-Étienne cedex 2

**RÉSUMÉ.** La comparaison entre des champs simulés par éléments finis (EF) et des données obtenues par mesures de champs devient une étape de plus en plus courante des procédures modernes d'identification. L'application des techniques standards d'interpolation (e.g., définition d'un voisinage et interpolation par fonctions linéaires, cubiques, B-splines) peut résulter en un filtrage implicite de l'information de haute fréquence et introduit des quantités qui ne sont liées ni au modèle de la structure ni aux mesures. Alternativement, l'utilisation du maillage et des fonctions de forme EF n'introduit aucune autre hypothèse que celles de la simulation. En ce sens, nous qualifions l'interpolation EF d'exacte. Cependant, une implémentation directe de l'interpolation EF conduit à une complexité algorithmique quadratique en le nombre d'éléments et de points de mesure, ce qui n'est pas acceptable lorsque de larges champs de données doivent être comparés plusieurs fois. Cet article décrit comment, à travers une pré-sélection des éléments par un maillage virtuel et l'utilisation de boîtes englobantes, la complexité algorithmique de l'interpolation EF peut rester linéaire. Une illustration des performances de l'algorithme est donnée pour une plaque trouée.

**ABSTRACT.** The comparison of finite elements (FE) and experimental data fields have become ever more prevalent in identification. The application of standard interpolation algorithms for this purpose (e.g. through linear, cubic or B-spline functions) may result in implicit filtering of high-frequency information, and creates a dependency of the field to non-physical parameters. Alternatively, using the existing FE mesh and shape functions to determine mesh degrees of freedom at data point coordinates makes no assumptions beyond those already made in the FE model. In this sense, interpolation using element shape functions is exact. However, basic implementations of this technique generally display a quadratic computation complexity with respect to mesh size and number of data points, which is impractical when large data fields must be compared repeatedly. This document aims at improving the technique's efficiency and linearizing its computation cost via a pre-selection strategy of the element containing each data point. With a combination of cross-products, bounding-boxes and a multi-dimensional storage array, the resulting algorithm shows a significant improvement in efficiency when compared to basic strategies. A sample application is given using a plate with a hole.

**MOTS-CLÉS:** Interpolation, éléments finis, mesures de champs

**KEYWORDS:** Interpolation, finite elements, field measurements
# 1. Introduction

Les comparaisons entre résultats de simulations par éléments finis (EF) et champs de données se généra- lisent à de nombreux domaines, dans le but de réaliser des identifications (e.g., en géophysique, (Fukushima et al., 2005), en mécanique pour des conditions limites, (Padmanabhan et al., 2006), ou des propriétés de matériaux (Bledzki et al., 1999, Bruno et al., 2002)), des transferts de champs (e.g., (Beckert, 2000)) ou des visualisations (e.g., (Nikishkov, 2003)). Les solveurs EF fournissent des résultats aux nœuds d'un maillage ou aux points d'intégration et, en général, les données ne coincident pas avec ces points. Pour comparer ces deux champs, il est nécessaire d'en projeter un sur l'autre.

Trois familles de techniques sont alors possibles : l'approximation, l'interpolation<fnref n="1" /> (Mathworks, 2004) et le placement de nœuds pour lequel le maillage EF est contraint à coincider avec les données (Bruno et al., 2002). Cependant, puisque la précision du calcul EF dépend du maillage, le placement de nœuds introduit un couplage malheureux entre la position des données expérimentales et la précision de la simulation. Les approches les plus générales de l'interpolation et de l'approximation (e.g., interpolation linéaire ou cubique à partir d'une triangulation de Delaunay, B-splines) introduisent des paramètres non physiques qui sont, par exemple, susceptibles d'éliminer l'information haute fréquence par un filtrage implicite. Dans cet article, seules les fonctions de forme des EF sont utilisées pour réaliser l'interpolation. Ceci n'introduit pas d'autres hypothèses que celles du modèle EF. L'efficacité de l'implémentation des techniques d'interpolation EF est discutée et une application à une plaque trouée conclut la présentation.

# 2. Techniques possibles de comparaison de champs

## Méthodes générales

La fonction griddata du logiciel Matlab (Mathworks, 2004) fournit quatre méthodes pour estimer la quantité $u(\mathbf{p})$ au point arbitraire $\mathbf{p}$ à partir d'un ensemble de données $u(\mathbf{n})$ : interpolations linéaires et cubiques à partir de triangulations de Delaunay, plus proche voisin et B-spline. Toutes ces méthodes peuvent être résumées par l'équation

$$u(p) = \sum_{i=1}^{\text{nb. de nœuds en } \epsilon(p)} W_i(\vec{x}(p), \vec{x}(\mathbf{n}^{\epsilon(p)})) u(n_i^{\epsilon(p)}). \quad [1]$$

où $W_i$ est une pondération fonction des coordonnées $\vec{x}$ du point d'interpolation et des voisins $\mathbf{n}^{\epsilon(p)}$ ($\epsilon(p)$ est un voisinage de $p$). Le krigeage (Cressie, 1993) est une autre technique d'interpolation populaire en géostatistique car il fournit une mesure de la variance de l'interpolateur. Cependant le krigeage, comme les autres méthodes générales, introduit des nouveaux paramètres sans sens physique direct (portée et échelle du variogramme en krigeage, taille de noyaux, degrés de polynômes, ...). De plus, ces méthodes ont un comportement mal défini sur des frontières internes telles qu'un trou dans une plaque car elles ne les décrivent pas explicitement.

## Projection des données sur l'espace du modèle EF

E. Pagnacco et D. Lemosse (Pagnacco et al., 2006) ont proposé d'approcher les données expérimentales aux nœuds du maillage EF. Pour cela, il faut trouver les valeurs aux nœuds $u(\mathbf{n}^{\text{mesh}})$ compatibles avec le champ EF et qui représentent le mieux possible, au sens des moindres carrés, les données expérimentales $u'(\mathbf{p})$,

$$\min_{u(\mathbf{n}^{\text{mesh}})} \sum_{i=1}^{\text{nb. de points en } \mathbf{p}} \left[ u'(p_i) - \sum_{j=1}^{\text{nb. de nœuds en } e(p_i)} N_j(\vec{\xi}(p_i)) u(n_j^{e(p_i)}) \right]^2, \quad [2]$$

où $e(p)$ est l'élément contenant le point $p$, $N_j$ est la $j$-ième fonction de forme et $\vec{\xi}(p)$ sont les coordonnées locales dans l'élément $e(p)$ du point $p$. Cette méthode a besoin d'un nombre de points expérimentaux supérieur au nombre de nœuds du maillage et que les données expérimentales couvrent la région de simulation

<footnote n="1">1. Contrairement à l'interpolation, l'approximation n'est pas obligée de passer par les points de mesures ou simulés.</footnote>
(pour que la solution de (2) soit unique). Le coût en calculs implique la détermination des coordonnées locales $\xi(p_i)$ de chaque point expérimental et la résolution du problème de moindres carrés linéaire (2). La projection filtre la partie hautes fréquences des expériences (le bruit expérimental).

## Interpolation à partir des fonctions de forme EF

L'interpolation à partir des fonctions de forme EF est réalisée en trois étapes : la détermination de l'élément d'appartenance du point expérimental $p$, $e(p)$, la transformation de ses coordonnées globales $\vec{x}(p)$ en coordonnées locales $\vec{\xi}(p)$ (transformation inverse), et le calcul du degré de liberté $u(p)$ à partir des fonctions de forme,

$$u(p) = \sum_{j=1}^{\text{nb. de nœuds en } e(p)} N_j(\vec{\xi}(p)) u(n_j^{e(p)}). \quad [3]$$

La quasi-intégralité du coût de calcul est pris par les deux premières étapes, les calculs de $e(p)$ et $\vec{\xi}(p)$. Nous revenons sur elles dans les deux parties suivantes. L'interpolation à partir des fonctions de forme EF n'impose aucune restriction sur le nombre et la position des points expérimentaux. De plus, les données expérimentales $u'$ ne sont pas modifiées.

# 3. Détermination de l'élément d'appartenance

## 3.1. Tests d'appartenance à un élément

Les tests d'appartenance d'un point $p$ à un élément peuvent être exacts ou approchés. Un test approché n'a d'intérêt que si il est substantiellement plus rapide que les tests précis. Il sert alors à éliminer, rapidement, des éléments auxquels $p$ n'appartient pas. Le test par boîte englobante (Rassineux, 1994) est le test approché utilisé dans ce travail. Il procède par la création d'un rectangle entourant l'élément fini (min et max appliqués aux coordonnées des nœuds) suivi de la vérification que $p$ est bien entre les min et max. Le test par boîte englobante élimine correctement de nombreux éléments auxquels $p$ n'appartient pas si ceux-ci sont de forme presque rectangulaire. Au contraire, pour des éléments très distordus, il diagnostiquera souvent à tord l'appartenance. Le test par boîte englobante est donc complété d'un test exact pour établir avec certitude l'appartenance ou non de $p$ à l'élément. Ici, les éléments utilisés sont linéaires et le test exact d'appartenance vérifie que $p$ est dans le cône formé par chaque nœud et les faces de l'élément qui en émanent. En 2D, les produits vectoriels ($\text{face}_1 \times \vec{n}_p$) et ($\vec{n}_p \times \text{face}_2$) doivent avoir le même signe.

## 3.2. Recherche dans l'ensemble des éléments : le maillage virtuel

La méthode la plus simple pour trouver l'élément d'appartenance de $p$ parmi les $N_e$ éléments est d'effectuer séquentiellement un test d'appartenance sur tous les éléments jusqu'à l'obtention du premier résultat positif. Dans certains cas, il faudra parcourir une grande partie des $N_e$ éléments avant de trouver celui qui contient $p$ (cf. exemple sur la figure 1(a)).

Le maillage virtuel (dont des idées proches sont déjà décrites dans (Rassineux, 1994)) a pour objectif de restreindre la liste des éléments pouvant potentiellement contenir $p$. Il repose, d'une part, sur une méthode efficace de détermination de l'appartenance de $p$ à une sous-région de l'espace, ce qui est typiquement obtenu par un pavage régulier qui se superpose au maillage. Dans le cas de pavés rectangulaires de tailles $dx_i$ et indicés par $\gamma_i$ selon chaque dimension $i$ (cf. figure 1(b)), on retrouve le pavé à partir des coordonnées de $p$ par une simple opération arithmétique,

$$\gamma_i = \text{floor} \left( \frac{x_i(p)}{dx_i} \right) \quad [4]$$

D'autre part, chaque pavé doit contenir une liste d'éléments finis, plus petite que l'ensemble de tous les éléments, ayant potentiellement une intersection non nulle avec lui. Ainsi, après l'application de la relation (4) pour trouver le pavé contenant $p$, les tests élémentaires d'appartenance ne sont appliqués qu'aux éléments
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Enumération (a) séquentielle, (b) avec un maillage virtuel.</figcaption></figure>

de la liste réduite. Par rapport à une recherche séquentielle, le maillage virtuel possède un coût supplémentaire d'initialisation car il faut, pour chaque pavé, construire la liste des éléments ayant potentiellement une intersection non nulle avec lui (de nouveau, un principe de boîtes englobantes est utilisé). Cependant cette opération est effectuée une seule fois et est largement amortie lorsque plusieurs $p$ sont à repérer dans le maillage. Une discussion détaillée du coût du maillage virtuel en fonction de la taille de ses éléments et de la distorsion du maillage est donnée dans (Silva et al., 2007).

### 3.3. Transformation inverse

La transformation inverse, depuis les coordonnées globales vers les coordonnées locales, $\vec{x}(p) \to \vec{\xi}(p)$ est une opération non triviale (contrairement à la transformation directe) nécessitant la résolution d'un système d'équations non linéaires. Une approche analytique pour les éléments QUAD4 et une approche itérative générale sont présentées en détails dans (Silva et al., 2007). Quelque soit la méthode utilisée pour la transformation inverse, la précision de l'opération peut être vérifiée en retransformant les coordonnées locales dans le repère global, $\vec{x}(p) \to \vec{\xi}(p) \to \vec{x}'(p)$, et en s'assurant que $||\vec{x}(p) - \vec{x}'(p)||$ est proche de 0. L'approche itérative est en fait la minimisation de $||\vec{x}(p) - \vec{x}'(p)||$ en changeant $\vec{\xi}(p)$.

## 4. Performance de la méthode

### 4.1. Complexité algorithmique

Les différentes méthodes de détermination de l'élément d'appartenance peuvent être comparées en termes de nombre moyen d'opérations (+, -, x, /, =, comparaisons, appels aux fonctions mathématiques standards), $\hat{T}_c$. $\hat{T}_c$ dépend du nombre d'éléments $N_e$, du nombre de points à interpoler $N_p$, du type d'éléments (QUAD4 ici), de la distorsion du maillage (à travers $\overline{N}_{cl}$ le nombre moyen de boîtes englobantes qui contiennent un point $p$), et du maillage virtuel (à travers $\overline{N}_{ve}$ et $\overline{N}_k$, le nombre moyen de pavés listant un élément fini et le nombre moyen d'éléments finis par pavé, respectivement). Les $\hat{T}_c$, calculés dans (Silva et al., 2007) sous des hypothèses d'uniformité de distribution des $p$ et des formes d'éléments, pour différentes combinaisons de tests d'appartenance et de recherche dans la liste des éléments sont les suivants :

<table><thead><tr><td>test appartenance</td><td>recherche élément</td><td>complexité T&#770;<sub>c</sub></td></tr></thead><tbody><tr><td>prod. vectoriel</td><td>séquentielle</td><td>18, 75N<sub>e</sub>N<sub>p</sub> + 41, 25N<sub>p</sub></td></tr><tr><td>boîte + prod. vectoriel</td><td>séquentielle</td><td>1, 25N<sub>e</sub>N<sub>p</sub> + 20, 75N<sub>cl</sub>N<sub>p</sub> + 110.5N<sub>p</sub> + 20N<sub>e</sub></td></tr><tr><td>boîte + prod. vectoriel</td><td>maillage virtuel</td><td>(N&#770;<sub>ve</sub> + 37)N<sub>e</sub> + [1, 25N&#770;<sub>k</sub> + 20, 75N&#770;<sub>cl</sub> + 116, 5]N<sub>p</sub> + 7</td></tr></tbody></table>

L'utilisation de boîtes englobantes et de maillage virtuel ajoute des coûts d'initialisation (les termes indépendants de $N_p$), mais les termes quadratiques en $O(N_p \times N_e)$ disparaissent.
## 4.2. Expériences numériques sur la plaque trouée

Des expériences numériques sont réalisées sur une plaque trouée en variant les nombres d'éléments ($N_e$) et de points à estimer ($N_p$). Ces expériences sont conduites au moyen d'une implémentation en C++ des algorithmes discutés, sur un ordinateur portable avec processeur Pentium 4. Le temps total d'interpolation $T_t$ (qui comprend la détermination de l'élément d'appartenance, la transformation inverse et l'application des fonctions de forme EF) est mesuré pour différentes combinaisons de $N_p$ et $N_e$. On remarque en figure 2(a) que $T_t$ augmente linéairement avec $N_p$ quand $N_e$ est maintenu constant, mais que la pente croît avec $N_e$, ce qui confirme la complexité théorique en $O(N_p \times N_e)$ des stratégies sans maillage virtuel. Avec le maillage virtuel (figure 2(b)), la pente reste constante car la complexité croît en $O(AN_e + BN_p)$, où A et B sont des constantes qui dépendent de la distorsion du maillage EF et de la finesse du maillage virtuel. Par contre, il y a une ordonnée à l'origine non nulle due à l'initialisation du maillage virtuel. L'utilisation du maillage virtuel accélère l'interpolation d'un facteur 10 à 60.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Temps d'interpolations mesurés</figcaption></figure>

L'algorithme d'interpolation (avec maillage virtuel) est maintenant utilisé pour calculer les cartes d'erreurs de la Figure 3. Il s'agit d'une plaque composite trouée ($E_{11} = 100$ GPa, $E_{22} = 36$ GPa, $G_{12} = 25$ GPa, $\nu_{12} = 0.54$ GPa) encastrée à gauche et soumise à une traction surfacique de 1555 kPa à droite. Les expériences sont simulées au moyen d'un calcul EF à 30000 éléments. L'erreur en un point $p$ est une distance Euclidienne normalisée des déplacements,

$$J(p) = \sqrt{\left(\frac{u(p)^{ref} - u(p)^{test}}{u_{max}^{ref} - u_{min}^{ref}}\right)^2 + \left(\frac{v(p)^{ref} - v(p)^{test}}{v_{max}^{ref} - v_{min}^{ref}}\right)^2} \quad [5]$$

Chaque carte comporte $N_p = 10000$ points à estimer et est calculée en environ 0.5 sec CPU. Les seules erreurs dont la figure 3 rend compte sont celles de la discrétisation puisque la géométrie et les matériaux sont les mêmes dans le modèle de référence et les modèles tests. Ainsi, lorsque $N_e$ augmente, l'erreur moyenne diminue car la frontière du trou est mieux décrite.

## 5. Bibliographie

* Beckert A., « Coupling fluid (CFD) and structural (FE) models using finite interpolation elements », *Aerospace Science and Technology*, vol. 4, n° 1, p. 13 - 22, January, 2000.

* Bledzki A. K., Kessler A., Rikards R., Chate A., « Determination of elastic constants of glass/epoxy unidirectional laminates by the vibration testing of plates », *Composites Science and Technology*, vol. 59, n° 13, p. 2015-2024, October, 1999.

* Bruno L., Furgiuele F. M., Pagnotta L., Poggialini A., « A full-field approach for the elastic characterization of anisotropic materials », *Optics and Lasers in Engineering*, vol. 37, p. 417-431, April, 2002.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Cartes d'erreurs *J*(*p*) entre une plaque “expérimentale” simulée avec *N*<sub>e</sub> = 30000 et des mo-dèles EF pour lesquels (a) *N*<sub>e</sub> = 700 and (b) *N*<sub>e</sub> = 1500.</figcaption></figure>

Cressie N., *Statistics for spatial data*, 1st edn, Wiley-Interscience, New York, 15 January, 1993.

Fukushima Y., Cayol V., Durand P., « Finding realistic dike models from interferometric synthetic aperture radar data : The February 2000 eruption at Piton de la Fournaise », *Journal of Geophysical Research-Solid Earth*, 23 March, 2005.

Mathworks, *Matlab : The language of technical computing*, 7 edn, Mathworks, Natic, June, 2004.

Nikishkov G., « Generating contours on FEM/BEM higher-order surfaces using Java 3D textures », *Advances in Engineering Software*, vol. 34, n° 8, p. 469-476, August, 2003.

Padmanabhan S., Hubner J. P., Kumar A. V., Ifju P. G., « Load and boundary condition calibration using full-field strain measurement », *Experimental Mechanics*, vol. 46, n° 5, p. 569 - 578, October, 2006.

Pagnacco E., Lemosse D., « A coupled FE based inverse strategy from displacement field measurement subject to an unknown distribution of forces », *Photomechanics*, Clermont-Ferrand, France, 10 July, 2006.

Rassineux A., Maillage automatique tridimensionnel par une méthode frontale pour la méthode des éléments finis, PhD thesis, Université Henri Poincaré, Vandoeuvre les Nancy, France, 1994. (in French).

Silva G. H. C., Le Riche R., Molimard J., Vautrin A., Exact and efficient interpolation using finite element shape func-tions, Technical report, LTDS, St. Etienne, France, 4 January, 2007. available https://hal.archives-ouvertes.fr/hal-00122640/en/.