# Optimisation topologique de composites pour maximiser la résistance à la fracture

Daicong Da, Julien Yvonnet, Liang Xia

► To cite this version:

Daicong Da, Julien Yvonnet, Liang Xia. Optimisation topologique de composites pour maximiser la résistance à la fracture. 14e Colloque National en Calcul des Structures (CSMA 2019), May 2019, Giens, France. hal-04096020

HAL Id: hal-04096020

https://hal.science/hal-04096020

Submitted on 12 May 2023

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 Attribution 4.0 International License
# Optimisation topologique de composites pour maximiser la résistance à la fracture

D. Da<fnref n="1" />, J. Yvonnet<fnref n="1" />, L. Xia<fnref n="2" />

**Résumé** — Nous proposons dans ce travail [1, 2] une méthodologie pour maximiser la résistance de composites de type matrice/inclusions pour lesquelles l'endommagement prend la forme de microfissuration dans la matrice et les interfaces. Pour cela, une approche par optimisation topologique est adoptée, dans laquelle la fonction objectif est directement définie comme l'énergie de fracturation de la structure composite, et obtenue par un calcul complet de fissuration par champ de phase, depuis l'initiation des fissures jusqu'à la rupture de la structure. La méthode BESO est adoptée pour l'optimisation topologique des inclusions.

**Mots clés** — optimisation topologique, microfissuration, champ de phase.

## 1 Introduction

Les propriétés physiques et mécaniques effectives des matériaux hétérogènes sont déterminées d'une part par leur composition et les fractions volumiques de leurs constituants, mais peuvent être drastiquement modifiées d'autre part par la géométrie des différentes phases et la présence d'interfaces. Par exemple, il a été montré que changer la topologie à fraction volumique constante peut permettre dans le cas linéaire de maximiser des modules élastiques ou d'obtenir des valeurs de ces modules n'existant pas naturellement dans la nature (coefficient de Poisson négatif, coefficients d'expansion thermique négatif, etc.). Récemment, les avancées dans les techniques de fabrication additive/impression 3D (voir par exemple [3]) on rendu possible la fabrication de tels matériaux, ouvrant la voie vers le design de matériaux totalement nouveaux. L'optimisation topologique des microstructures a été étudiée de manière extensive pour les propriétés linéaires (voir par exemple [4]), et plus récemment dans un cadre non linéaire (voir par exemple [5]). L'une des préoccupations majeure en ingénierie est la conception de matériaux plus résistants à la rupture. Dans ce travail, nous explorons la possibilité d'utiliser l'optimisation topologique dans cet objectif dans les composites. Ce sujet est resté relativement peu exploré jusqu'à présent. L'idée est d'utiliser une méthode de simulation numérique robuste (ici la méthode de champs de phase [6, 7] étendue à la fissuration interfaciale [8]) permettant la simulation complète de la rupture dans une structure composite depuis l'initiation des microfissures en vue d'obtenir une courbe de réponse et optimiser une fonction objectif, ici définie comme l'énergie de fissuration, et d'ensuite réaliser une optimisation topologique des inclusions en calculant la sensibilité de cette fonction par rapport à la géométrie des inclusions (voir Fig. 1).

## 2 Modélisation de la microfissuration avec endommagement interfacial par la méthode de champ de phase

Dans cette section nous décrivons brièvement la méthode utilisée pour modéliser la microfissuration dans la structure composite. La méthode se base sur l'approche variationnelle de l'endommagement [6, 7] étendue ici à l'endommagement interfacial, afin de permettre une initiation de fissures à partir des interfaces dans une structure non initialement fissurée. Nous utilisons le cadre proposé dans Nguyen et al. [8]. Les principales équations sont rappelées ci-dessous. L'énergie totale pour le solide en présence

<footnote n="1">¹ Université Paris-Est, Laboratoire Modélisation et Simulation Multi Echelle MSME, CNRS UMR 8208, 5 bd Descartes,
77454 Marne-la-Vallée, France, julien.yvonnet@univ-paris-est.fr</footnote>

<footnote n="2">² State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, Wuhan, China</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – (a) Structure composite en présence de fissures matricielles et interfaciales ; (b) réponse mécanique de la structure jusqu’à sa rupture.</figcaption></figure>

de fissures dans les phases et les interfaces est exprimée par :

$$E = \int_{\Omega} W_{u}^{e}(\varepsilon^{e}(\mathbf{u}, \beta), d) d\Omega + \int_{\Omega} [1 - \beta(\mathbf{x})] g_{c} \gamma_{d}(d) d\Omega + \int_{\Omega} \psi^{l}(\mathbf{w}) \gamma_{\beta}(\beta) d\Omega, \quad (1)$$

où $W_u^e$ est la densité d'énergie élastique, $\gamma_d$ est une densité de fissures et $\gamma_\beta$ une densité d'interfaces (voir [8] pour plus de détails), $g_c$ désigne la ténacité dans une phase et $\psi^l$ est une densité d'énergie dépendant du saut de déplacements aux interfaces et $\varepsilon^e$ est définie telle que :

$$\varepsilon = \varepsilon^e + \bar{\varepsilon}, \quad (2)$$

où $\bar{\varepsilon}$ est la partie de la déformation induite par un saut de déplacements régularisé aux interfaces. Par minimisation de l'énergie ci-dessus sous la contrainte d'irréversibilité du champ d'endommagement $d(\mathbf{x})$, l'évolution de $d(\mathbf{x})$ peut être obtenue par résolution du problème suivant :

$$\left\{
\begin{array}{l}
2(1-d)\mathcal{H} - (1-\beta)\frac{g_c}{l}(d-l^2\nabla^2d) = 0, \quad \text{in } \Omega \\
d(\mathbf{x}) = 1, \quad \text{on } \Gamma^c \\
\nabla d(\mathbf{x}) \cdot \mathbf{n} = 0, \quad \text{on } \partial\Omega,
\end{array}
\right.
\quad (3)$$

où $\beta(\mathbf{x})$ est un champ régularisé dévrivant la présence des interfaces et $\mathcal{H}$ une fonction d'histoire (voir [8]). La forme faible associée à (3) est donnée par :

$$\int_{\Omega} \left\{ (2\mathcal{H} + [1-\beta]\frac{g_c}{l})d\delta d + [1-\beta]g_c l \nabla d \cdot \nabla(\delta d) \right\} d\Omega = \int_{\Omega} 2\mathcal{H}\delta dd\Omega. \quad (4)$$

La variation de l'énergie totale par rapport au champ de déplacements $\mathbf{u}$ fournit la forme faible associée au problème d'équilibre mécanique :

$$\int_{\Omega} \frac{\partial W_{u}^{e}}{\partial \varepsilon^{e}} : \varepsilon^{e}(\delta \mathbf{u}) d\Omega + \int_{\Omega} \frac{\partial \psi^{l}(\mathbf{w})}{\delta \mathbf{w}} \cdot \partial \mathbf{w} \gamma_{\beta}(\beta) d\Omega = \int_{\Omega} \mathbf{f} \cdot \delta \mathbf{u} d\Omega + \int_{\partial \Omega_{F}} \bar{\mathbf{F}} \cdot \delta \mathbf{u} d\Gamma \quad (5)$$

et $\mathbf{w}$ est un saut de déplacements régularisé (voir [8]).

### 3 Optimisation topologique

On considère un domaine $\Omega$ discrétisé en $N_i$ éléments finis. Chaque élément $i$ est associé à une densité fictive utilisée pour la définition de la topologie $\rho_i$. Nous adoptons le cadre d'optimisation topologique de la méthode BESO (Bi-directional evolutionary structural optimization [9]) pour sa simplicité de mise en œuvre. Dans ce contexte, les propriétés de chaque phase sont interpolées par rapport aux densités comme :
$$E_i = \rho_i E_{\text{inc}} + (1 - \rho_i) E_{\text{mat}} \quad (6)$$

où $E_{\text{inc}}$ et $E_{\text{mat}}$ désignent les modules de Young des inclusions et des phases, respectivement. Les coefficients de Poisson et les autres propriétés matérielles sont supposées rester constantes pour des raisons de simplicité. Pour un chargement imposé en déplacement, la fonction objectif d'énergie de fracturation de la structure est équivalente à maximiser le travail durant toute l'histoire du chargement depuis l'initiation jusqu'à la rupture. En pratique, le travail des efforts extérieurs $J$ est évalué numériquement par intégration comme :

$$J \approx \frac{1}{2} \sum_{n=1}^{n_{\text{load}}} \left( \mathbf{f}_{\text{ext}}^{(n)} + \mathbf{f}_{\text{ext}}^{(n-1)} \right)^{\mathrm{T}} \Delta \mathbf{u}^{(n)}. \quad (7)$$

Ci-dessus, $n_{\text{load}}$ est le nombre total d'incréments de déplacements, $\Delta \mathbf{u}^{(n)}$ est le n-ième incrément de déplacement et $\mathbf{f}_{\text{ext}}^{(n)}$ est la force nodale généralisée au n-ième incrément de chargement. Le problème d'optimisation topologique peut être défini comme :

$$
\begin{array}{ll}
\max_{\rho} & J(\rho, \mathbf{u}, \mathbf{d}, \beta) \\
\text{subjected to:} & \mathcal{R} = 0 \\
& V(\rho) = \sum \rho_i v_i = V_{\text{req}} \\
& \rho_i = 0 \text{ or } 1, i = 1, \dots, N_i.
\end{array}
\quad (8)
$$

Ci-dessus, $v_i$ désigne le volume du i-ème élément, $V(\rho)$ et $V_{\text{req}}$ sont le volume total et le volume visé, et $\mathcal{R}$ désigne le résidu associé à l'équilibre mécanique :

$$\mathcal{R} = \mathbf{f}_{\text{ext}} - \mathbf{f}_{\text{int}}. \quad (9)$$

Dans (9), $\mathbf{f}_{\text{int}}$ est défini dans chaque élément comme une force interne et donné en fonction de la variable de design topologique $\rho_i$ et des contraintes par :

$$\mathbf{f}_{\text{int}} = \sum_{i=1}^{N_i} \rho_i \int_{\Omega_i} \mathbf{B}^{\mathrm{T}} \sigma^e d\Omega_i. \quad (10)$$

La sensibilité $\partial \hat{J}/\partial \rho_i$ est donnée par :

$$\frac{\partial \hat{J}}{\partial \rho_i} = -\frac{1}{2} \sum_{n=1}^{n_{\text{load}}} \left\{ \left( \lambda^{(n)} \right)^T \int_{\Omega_i} \mathbf{B}^T (\boldsymbol{\sigma}^e)^{(n)} d\Omega_i + \left( \boldsymbol{\mu}^{(n)} \right)^T \int_{\Omega_i} \mathbf{B}^T (\boldsymbol{\sigma}^e)^{(n-1)} d\Omega_i \right\}. \quad (11)$$

avec (voir [2])

$$\boldsymbol{\lambda}_{\text{F}}^{(n)} = \left( \mathbf{K}_{\text{tan,FF}}^{(n)} \right)^{-1} \mathbf{K}_{\text{tan,FE}}^{(i)} \Delta \mathbf{u}_{\text{E}}^{(n)} \quad (12)$$

et

$$\boldsymbol{\mu}_{\text{F}}^{(n)} = \left( \mathbf{K}_{\text{tan,FF}}^{(n-1)} \right)^{-1} \mathbf{K}_{\text{tan,FE}}^{(n-1)} \Delta \mathbf{u}_{\text{E}}^{(n)}. \quad (13)$$

# 4 Résultats

Nous appliquons la démarche à plusieurs exemple de structures contenant des inclusions dont les interfaces peuvent s'endommager pour participer à l'évolution de la microfissuration. La structure étudiée dans la figure 2 (a) est une plaque renforcée par des inclusions initialement de formes carrées (Fig. 2 (b)). Un chargement en traction est appliqué sur les bords de la structure. Une fissure initiale est présente sur le bord inférieur. On peut voir dans la figure 2 (c) le résultat de l'optimisation topologique des différentes inclusions. Dans les figures 3 (a)-(d) on peut observer l'évolution de la microfissuration au cours du chargement dans la structure optimisée. L'endommagement pouvant s'initier aux interfaces en plus de la matrice, les différentes inclusions peuvent adopter des formes distinctes pour maximiser l'énergie de fracturation. L'évolution de la croissance des microfissures dans la géométrie optimisée peut être observée dans la figure 3.
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – (a) Structure contenant des renforts, (b) géométrie des renforts avant optimisation, (c) après optimisation.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Evolution de la croissance des microfissures pour la géométrie optimisée.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Comparaison entre les réponses des structures non optimisées et optimisées.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Optimisation topologique des inclusions pour une structure périodique, avec contrainte de périodicité sur la forme des inclusions.</figcaption></figure>

On peut comparer dans la Fig. 4 la différence entre la réponse de la structure initiale et optimisée, se traduisant par une augmentation significative de la contrainte maximale à la rupture ainsi que de l'énergie de fracturation.

Dans l'exemple suivant, on considère une structure composite soumise à de la flexion 3 points dont les inclusions sont initialement cylindriques, sans pré-fissure (voir Fig. 5). on a ici forcé la périodicité des inclusions durant l'optimisation. La géométrie des inclusions optimisée, le trajet de fissuration et la comparaison entre les réponses des structures optimisées et non optimisées sont fournies dans la figure 5, avec une augmentation de 82.5 % de l'énergie de fissuration, et montrant le potentiel de l'approche. Pour ces deux exemples, la fraction volumique des inclusions est maintenue constante.

Références

[1] L. Xia, D. Da, and J. Yvonnet. Topology optimization for maximizing the fracture resistance of quasi-brittle composites. *Computer Methods in Applied Mechanics and Engineering*, 332:234–254, 2018.

[2] D. Da, J. Yvonnet, L. Xia, and G. Li. Topology optimization of particle-matrix composites for optimal fracture resistance taking into account interfacial damage. *International Journal for Numerical Methods in Engineering*, 2018.

[3] Z. Quan, Z. Larimore, A. Wu, J. Yu, X. Qin, M. Mirotznik, J. Suhr, J.-H. Byun, Y. Oh, and T.-W. Chou. Microstructural design and additive manufacturing and characterization of 3d orthogonal short carbon fiber/acrylonitrile-butadiene-styrene preform and composite. *Composites Science and Technology*, 126:139–148, 2016.

[4] J. E. Cadman, S. Zhou, Y. Chen, and Q. Li. On design of multi-functional microstructural materials. *Journal of Materials Science*, 48(1):51–66, 2013.

[5] F. Fritzen, L. Xia, M. Leuschner, and P. Breitkopf. Topology optimization of multiscale elastoviscoplastic structures. *International Journal for Numerical Methods in Engineering*, 106(6):430–453, 2016.

[6] G.A. Francfort and J.J. Marigo. Revisiting brittle fracture as an energy minimization problem. *J. Mech. Phys. Solids*, 46(8):1319–1342, 1998.

[7] C. Miehe, M. Hofacker, and F. Welschinger. A phase field model for rate-independent crack propagation : Robust algorithmic implementation based on operator splits. *Computer Methods in Applied Mechanics and Engineering*, 199:2776–2778, 2010.

[8] T.T. Nguyen, J. Yvonnet, Q.-Z. Zhu, M. Bornert, and C. Chateau. A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography. *Computer Methods in Applied Mechanics and Engineering*, (312):567–595, 2016.

[9] X. Huang and Y.M. Xie. Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method. *Finite Elements in Analysis and Design*, 43(14):1039–1049, 2007.