# Une méthode de potentiels numériques explicites pour le calcul multi-échelles des composites non linéaires

Julien Yvonnet, David Gonzalez, Qi-Chang He

► To cite this version:

Julien Yvonnet, David Gonzalez, Qi-Chang He. Une méthode de potentiels numériques explicites pour le calcul multi-échelles des composites non linéaires. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01413984

HAL Id: hal-01413984

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

Submitted on 12 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
# Une méthode de potentiels numériques explicites pour le calcul multi échelles des composites non linéaires

J. Yvonnet¹, D. Gonzalez², Q.-C. He¹

¹ Université Paris-Est, Laboratoire Modélisation et Simulation Multi Echelle MSME FRE
CNRS 3160, 5 Bd Descartes, 77454 Marne-la-Vallée cedex 2, France
julien.yvonnet@univ-paris-est.fr

² I3A, University of saragoza,
Maria de Luna, 3, E-50015 Zaragoza, Spain

**Résumé** – Nous proposons ici une nouvelle méthode multi-échelle pour calculer la réponse des structures hétérogènes non linéaires. Par rapport aux approches existantes basées sur des calculs simultanés aux deux échelles, l'approche découple les calculs microscopiques et macroscopiques. Pour cela, la fonction densité d'énergie effective est calculée en des points discrets de l'espace du chargement macroscopique. Des méthodes d'interpolation pour les grandes dimensions permettent ensuite d'accéder aux contraintes et au tenseur élastique tangent.

**Mots clés** – Méthodes multi-échelles, matériaux non linéaires, potentiels explicites numériques

## 1 Introduction

Le calcul de structures non linéaires est un problème fondamental pour les sciences de l'ingénieur. Cependant lorsque le matériau constituant la structure est hétérogène, il est nécessaire de réaliser une homogénéisation du comportement non linéaire macroscopique. Les difficultés par rapport au cas linéaire sont nombreuses : la forme générale de la loi de comportement n'est pas connue, et l'homogénéisation non linéaire nécessite de résoudre des équations aux dérivées partielles non linéaires avec des coefficients aléatoires ou périodiques pour des conditions aux limites associés à un nombre théoriquement infini de cas de chargement macroscopiques. Un très grand nombre de travaux analytiques et numériques ont été consacrés à ce problème au cours des dernières décennies, mais les limitations restent nombreuses.

Nous proposons ici une approche numérique nouvelle permettant le calcul des structures hétérogènes non linéaires [1]. Par rapport aux approches existantes, cette technique offre les avantages suivants : a) le calcul multi échelles est réalisé de manière « séquentielle » et non « simultanée » : il est alors possible, pour une microstructure donnée, d'effectuer des calculs sur un nombre infini de structures et de cas de chargements macroscopiques, sans nouveaux calculs locaux. Il est également possible de considérer plus qu'une échelle de description microstructurale, ce qui aujourd'hui est inconcevable en terme de temps de calculs avec les approches disponibles ; b) les calculs préparatoires permettant l'identification de la loi macroscopique et le calcul de structures sont découpés. Une fois les calculs préliminaires effectués (qui sont d'ailleurs totalement parallèlisables), le temps lié au calcul de structure est du même ordre de grandeur qu'un calcul élément fini classique lorsque la loi macroscopique est donnée sous forme analytique.

L'idée centrale de cette technique est de calculer numériquement des valeurs discrètes d'un potentiel dont la loi macroscopique dérive et d'interpoler ses valeurs dans l'hyperespace des
déformations macroscopiques pour accéder aux tenseurs des contraintes et au tenseur élastique tangent par dérivation de la fonction continue construite. Les détails de la technique sont fournis ci-dessous.

## 2 Rappels d'homogénéisation non linéaire

Soit un volume élémentaire représentatif (v.e.r.) défini dans un domaine $\Omega$ et de bord $\partial\Omega$ associé à une microstructure. On suppose que le v.e.r. est constitué de $r$ phases, chacune définie dans un domaine $\Omega_r$. L'opération de moyenne sur le volume de $\Omega$ et de $\Omega_r$ est notée $\langle \cdot \rangle$ et $\langle \cdot \rangle_r$, respectivement. Soit $\chi^{(r)}(\mathbf{x})$ la fonction caractéristique associée à une phase $r$ telle que $\chi^{(r)}(\mathbf{x}) = 1$ dans $\Omega_r$ et $\chi^{(r)}(\mathbf{x}) = 0$ ailleurs. Dans chaque phase, nous supposons qu'un potentiel de déformation élastique $\omega^{(r)}(\boldsymbol{\epsilon})$, ou fonction densité d'énergie est connu. Nous considérons ici le cas des petites perturbations, le tenseur $\boldsymbol{\epsilon}$ étant le tenseur des déformations linéarisé. Le potentiel local est alors défini par

$$ \omega(\mathbf{x}, \boldsymbol{\epsilon}) = \sum_{r=1}^{N} \chi^{(r)}(\mathbf{x}) \omega^{(r)}(\boldsymbol{\epsilon}) \qquad (1) $$

La loi de comportement locale peut alors être déterminée par :

$$ \boldsymbol{\sigma}(\mathbf{x}, \boldsymbol{\epsilon}) = \frac{\partial \omega(\mathbf{x}, \boldsymbol{\epsilon})}{\partial \boldsymbol{\epsilon}} \qquad (2) $$

où $\boldsymbol{\sigma}$ est le tenseur des contraintes de Cauchy satisfaisant les équations d'équilibre. Nous définissons les grandeurs macroscopiques, ou effectives, comme les moyennes sur le volume du v.e.r., soient $\bar{\boldsymbol{\epsilon}} = \langle \boldsymbol{\epsilon} \rangle$ et $\bar{\boldsymbol{\sigma}} = \langle \boldsymbol{\sigma} \rangle$. Afin de résoudre le problème local, des conditions aux limites doivent être imposées sur $\partial\Omega$. Par exemple, les conditions aux limites périodiques sont données sous la forme :

$$ \mathbf{u} - \bar{\boldsymbol{\epsilon}}\mathbf{x} \text{ périodique sur } \partial\Omega \text{ et } \boldsymbol{\sigma}\mathbf{n} \text{ anti-périodique sur } \partial\Omega. \qquad (3) $$

Pour un état de déformation macroscopique $\bar{\boldsymbol{\epsilon}}$, le problème non linéaire peut être résolu par exemple par un calcul élément fini. La loi de comportement macroscopique est donnée sous la forme

$$ \bar{\boldsymbol{\sigma}}(\mathbf{x}, \boldsymbol{\epsilon}) = \frac{\partial \langle \omega(\mathbf{x}, \boldsymbol{\epsilon}) \rangle}{\partial \boldsymbol{\epsilon}} = \frac{\partial \bar{\omega}(\bar{\boldsymbol{\epsilon}})}{\partial \boldsymbol{\epsilon}} \qquad (4) $$

L'idée centrale de ce travail est de proposer une forme numérique pour $\bar{\omega}$, en vue de calculer efficacement la loi de comportement $\bar{\boldsymbol{\sigma}} = f(\bar{\boldsymbol{\epsilon}})$ dans un calcul de structure classique de type éléments finis.

## 3 Méthode des potentiels explicites numériques

Pour la suite, nous adoptons les notations suivantes : $\bar{\boldsymbol{\epsilon}}_1 = \bar{\boldsymbol{\epsilon}}_{11}$, $\bar{\boldsymbol{\epsilon}}_2 = \bar{\boldsymbol{\epsilon}}_{22}$, $\bar{\boldsymbol{\epsilon}}_3 = \bar{\boldsymbol{\epsilon}}_{33}$, $\bar{\boldsymbol{\epsilon}}_4 = \sqrt{2}\bar{\boldsymbol{\epsilon}}_{23}$, $\bar{\boldsymbol{\epsilon}}_5 = \sqrt{2}\bar{\boldsymbol{\epsilon}}_{13}$, $\bar{\boldsymbol{\epsilon}}_6 = \sqrt{2}\bar{\boldsymbol{\epsilon}}_{12}$. Les mêmes notations sont utilisées pour les composantes des contraintes. Les composantes des contraintes et du tenseur élastique tangent $\mathbf{L}$ sont alors données par :

$$ \bar{\boldsymbol{\sigma}}_{\alpha} = \frac{\partial \bar{\omega}}{\partial \bar{\boldsymbol{\epsilon}}_{\alpha}} \quad \text{et} \quad \bar{L}_{\alpha\beta} = \frac{\partial^2 \bar{\omega}}{\partial \bar{\boldsymbol{\epsilon}}_{\alpha} \partial \bar{\boldsymbol{\epsilon}}_{\beta}} \qquad (5) $$

Dans le cas général, le potentiel de déformation élastique est une fonction scalaire dans un espace à six dimensions :

$$ \bar{\omega} = \bar{\omega}(\bar{\boldsymbol{\epsilon}}_1, \bar{\boldsymbol{\epsilon}}_2, \dots, \bar{\boldsymbol{\epsilon}}_6) \qquad (5) $$
Cependant, on peut être amené à étudier un problème sous certaines hypothèses, par exemple dans le plan $x_1 - x_2$. Dans ce cas, $\bar{\omega}$ dépend seulement de $\bar{\varepsilon}_1, \bar{\varepsilon}_2$ et $\bar{\varepsilon}_6$, et les valeurs de $\alpha$ et $\beta$ sont alors prises dans l'ensemble $\{1,2,6\}$. Nous présentons ici la méthode dans le cadre général, mais celle-ci peut être adaptée facilement pour des problèmes faisant intervenir moins de composantes de déformations macroscopiques.

## 3.1 Représentation discrète du potentiel

Dans un premier temps, nous discrétisons un domaine $\Delta$ dans l'espace des déformations macroscopiques à six dimensions $E$. Une discrétisation non régulière peut être envisagée. Nous présentons ici le cas d'une discrétisation $E$ en une grille régulière, par souci de simplicité. Pour cela, chaque axe $\bar{\varepsilon}_\alpha$ est gradué uniformément. Des valeurs minimales et maximales de $\bar{\varepsilon}_\alpha$ doivent être choisies, restreintes par l'hypothèse des petites déformations. Nous introduisons un ensemble de points $\{\xi_0^\alpha, \xi_1^\alpha, \dots, \xi_m^\alpha\}$ définis sur l'intervalle $[\bar{\varepsilon}_\alpha^{\min}, \bar{\varepsilon}_\alpha^{\max}]$, avec $\bar{\varepsilon}_\alpha^{\min} = \xi_0^\alpha < \xi_1^\alpha < \dots < \xi_{m_\alpha}^\alpha = \bar{\varepsilon}_\alpha^{\max}$.

Le domaine $\Delta = [\bar{\varepsilon}_1^{\min}, \bar{\varepsilon}_1^{\max}] \times [\bar{\varepsilon}_2^{\min}, \bar{\varepsilon}_2^{\max}] \times \dots \times [\bar{\varepsilon}_6^{\min}, \bar{\varepsilon}_6^{\max}]$ de $E$ est ainsi discrétisé en un ensemble fini de points, appelés nœuds. Si l'on considère une grille régulière, nous posons $p = m_1 + 1 = m_2 + 1 = \dots = m_6 + 1$. Pour un problème faisant intervenir $d$ composantes de déformations macroscopiques, le domaine $\Delta$ contient $p^d$ nœuds.

Nous évaluons ensuite $\bar{\omega}$ en chaque nœud de $\Delta$ par la méthode des éléments finis. Une représentation du domaine discrétisé $\Delta$ est illustrée dans la figure 1.

<figure><img src="image_2.png" /><figcaption>Figure 1 - Schématisation de l'espace discrétisé des déformations macroscopiques. Différentes déformées d'un v.e.r. sont représentées en différents nœuds du domaine.</figcaption></figure>

Soient les composantes de déformations macroscopiques $(\xi_{i_1}^1, \xi_{i_2}^1, \dots, \xi_{i_6}^1)$ en un nœud $i_1 i_2 \dots i_6$. Nous désignons par le suite l'évaluation du potentiel effectif en ce point par $\bar{\omega}_{i_1 i_2 \dots i_6}$. L'ensemble des valeurs évaluées dans les différentes nœuds est associé à une hypermatrice $\bar{W} = [\bar{\omega}_{i_1 i_2 \dots i_6}]$. Après avoir calculé les valeurs de $\bar{W}$ en résolvant $p^d$ problèmes éléments finis, il est alors possible d'accéder aux composantes du tenseur des contraintes et du tenseur élastique tangent par interpolation des valeurs dans $\bar{W}$.
## 3.2 Méthodes d'interpolation

L'interpolation doit être $C^2$, afin d'obtenir une représentation continue du tenseur élastique tangent pour le calcul des structures dans le cadre d'une méthode de résolution de Newton-Raphson. Dans [1] nous avons proposé deux techniques d'interpolation. La première technique est basée sur une méthode de spline cubiques multidimensionnelles [2]. Cette première approche est bien adaptée lorsque $\Delta$ est défini sur un espace de faibles dimensions. Pour les dimensions plus grandes, ce type d'interpolation est extrêmement coûteux aussi bien en temps de calcul qu'en taille mémoire. Une seconde technique d'interpolation, bien adaptée aux espaces de grandes dimensions, est décrite ci-dessous.

Nous utilisons une décomposition en produits tensoriels [3] de l'hypermatrice $\overline{W}$ :

$$ \overline{W} \approx \sum_{r=1}^{R} \phi_{1}^{r} \otimes \phi_{2}^{r} \otimes \dots \otimes \phi_{6}^{r} \qquad (5) $$

où $\phi_i^r$ ($i=1,...,6$ ; $r = 1,...,R$) sont des vecteurs à valeurs réelles associés à $\overline{\epsilon}_i$ et $R$ est un entier. Les vecteurs $\phi_i^r$ sont déterminés en résolvant un problème des moindres carrés pour une valeur de $R$ donnée :

$$ \inf_{\phi_j^r} \left\| \overline{W} - \sum_{r=1}^{R} \phi_1^r \otimes \phi_2^r \otimes \dots \otimes \phi_6^r \right\|^2, \quad r = 1, \dots, R, \quad j = 1, \dots, 6 \qquad (5) $$

Pour atteindre une précision donnée, $R$ peut être augmenté jusqu'à ce que

$$ \left\| \overline{W} - \sum_{r=1}^{R} \phi_{1}^{r} \otimes \phi_{2}^{r} \otimes \dots \otimes \phi_{6}^{r} \right\| < \delta \qquad (5) $$

avec $\delta$ un paramètre de tolérance. Le problème (7) étant non linéaire par rapport aux inconnues $\phi_j^r$, une procédure itérative est nécessaire pour le résoudre. Un algorithme efficace et simple à mettre en œuvre est l'algorithme des moindres carrés alternés [4]. Une fois les vecteurs $\phi_j^r$ déterminés, nous obtenons une représentation en variables séparées de $\overline{\omega}$ sous la forme :

$$ \overline{\omega}(\overline{\epsilon}_1, \overline{\epsilon}_2, \dots, \overline{\epsilon}_6) \approx \sum_{r=1}^{R} \tilde{\phi}_1^r(\overline{\epsilon}_1) \tilde{\phi}_2^r(\overline{\epsilon}_2) \dots \tilde{\phi}_6^r(\overline{\epsilon}_6) \qquad (5) $$

où $\tilde{\phi}_j^r(\overline{\epsilon}_j) = \sum_{k=1}^{n} N_k(\overline{\epsilon}_k)\{\phi_j^r\}_k$ sont les valeurs interpolées de $\phi_j^r$. Ici $N_k$ est une fonction cubique spline unidimensionnelle. Les contraintes peuvent alors être approximées par :

$$ \overline{\sigma}_{\alpha}(\overline{\epsilon}_1, \overline{\epsilon}_2, \dots, \overline{\epsilon}_6) \approx \sum_{r=1}^{R} \left\{ \left( \prod_{\beta \neq \alpha} \tilde{\phi}_{\beta}^{r}(\overline{\epsilon}_{\beta}) \right) \frac{\partial \tilde{\phi}_{\alpha}^{r}(\overline{\epsilon}_{\alpha})}{\partial \overline{\epsilon}_{\alpha}} \right\} \qquad (5) $$

Les composantes du tenseur élastique tangent peuvent être obtenues de la même façon, par dérivées seconde du potentiel continu construit :

$$ \left\{ \begin{array}{l} \bar{L}_{\alpha\beta}(\bar{\epsilon}_1, \bar{\epsilon}_2, \dots, \bar{\epsilon}_6) \approx \sum_{r=1}^{R} \left( \left\{ \prod_{\gamma \neq \alpha} \tilde{\phi}_{\gamma}^{r}(\bar{\epsilon}_{\gamma}) \right\} \frac{\partial \tilde{\phi}_{\beta}^{r}(\bar{\epsilon}_{\beta})}{\partial \bar{\epsilon}_{\beta}} \frac{\partial \tilde{\phi}_{\alpha}^{r}(\bar{\epsilon}_{\alpha})}{\partial \bar{\epsilon}_{\alpha}} \right) \quad \text{si } \alpha \neq \beta \\ \\ \bar{L}_{\alpha\beta}(\bar{\epsilon}_1, \bar{\epsilon}_2, \dots, \bar{\epsilon}_6) \approx \sum_{r=1}^{R} \left( \left\{ \prod_{\gamma \neq \alpha} \tilde{\phi}_{\gamma}^{r}(\bar{\epsilon}_{\gamma}) \right\} \frac{\partial^2 \tilde{\phi}_{\alpha}^{r}(\bar{\epsilon}_{\alpha})}{\partial \bar{\epsilon}_{\alpha}^2} \right) \quad \text{si } \alpha = \beta \end{array} \right. \qquad (5) $$
## 4 Exemple numérique

Dans cet exemple, une structure décrite dans la figure 2 est constituée d'un matériau hétérogène, dont le v.e.r. est également représenté. Le potentiel local dans chaque phase du v.e.r. est un potentiel de type loi puissance. Les détails peuvent être trouvés dans [1]. L'hypermatrice $\overline{W}$ est calculée suivant la procédure décrite ci-dessus. Les contraintes et le module tangent élastique approximés, dérivées de la fonction continue interpolée, sont utilisés dans l'algorithme de Newton-Raphson pour déterminer l'équilibre de la structure. Pour évaluer la qualité de la méthode, nous comparons note solution avec un calcul Eléments finis multi-niveaux [5]. Les résultats de la figure 2 montrent que les résultats sont en très bonne adéquation avec le calcul multi-niveaux utilisé comme solution de référence.

<figure><img src="image_3.png" /><figcaption>Figure 2 - Structure hétérogènes, maillage, et volume élémentaire représentatif.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3 - Structure hétérogènes, maillage, et volume élémentaire représentatif.</figcaption></figure>
Références

[1] J. Yvonnet, D. Gonzalez, Q.-C. He. A method using numerically explicit potentials for the homogenization of arbitrarily nonlinear composites, submitted.

[2] C. Habermann, F. Kindermann. Multidimensional spline interpolation : theory and applications. *Comput. Econ.* 30, 153-169, 2007.

[3] H.A.L. Kiers. Toward a standardized notation and terminology in multiway analysis. *J. Chmometr.* 14, 105-122, 2000.

[4] T. Zhang, G.H. Golub. Rank-one approximation. *SIAM J. Matrix Anal. Appl.* 23(2), 534-550, 2001.

[5] F. Feyel. Multiscale $FE^2$ elastoviscoplastic analysis of composite structure. *Compt. Mater. Sci.* 16(1-4), 433-354, 1999.