# Modelisations multi-échelles du comportement adoucissant des matériaux à matrice cimentaire

Nathan Benkemoun, Martin Hautefeuille, Jean-Baptiste Colliat, Adnan

Ibrahimbegovic

► To cite this version:

Nathan Benkemoun, Martin Hautefeuille, Jean-Baptiste Colliat, Adnan Ibrahimbegovic. Modelisations multi-échelles du comportement adoucissant des matériaux à matrice cimentaire. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01416070

HAL Id: hal-01416070

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

Submitted on 14 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
Modelisations multi-échelles du comportement
adoucissant des matériaux à matrice cimentaire

N. Benkemoun¹, M. Hautefeuille¹,², J.-B. colliat¹, A. Ibrahimbegović¹

¹ LMT-Cachan (ENS Cachan/CNRS/UPMC/PRES UniverSud Paris)
61 av. du Président Wilson
F-94230 Cachan, France
{benkemoun,hautefeuille, colliat,ai} @lmt.ens-cachan.fr

² Institut für Wissenschaftliches Rechnen (TU-Braunschweig)
D-38092 Braunschweig
wire@tu-bs.de

**Résumé** — Cette communication présente un méso modèle pour le béton et, plus généralement, pour les matériaux hétérogènes quasi-fragiles ainsi que son utilisation, comme description fine, au sein d'une méthode multi-échelles intégrée.

**Mots clés** — méthode des Eléments Finis, méthode multi-échelles intégrée, discontinuités fortes, comportement fragile, éléments biphasiques.

# 1 Introduction

La modélisation d'un matériau hétérogène comme le béton est étroitement liée à l'échelle d'observation. D'un point de vue macroscopique, celui-ci est considéré comme homogène, introduisant implicitement la notion de VER [1]. Les modèles macroscopiques proposés comme la plasticité ou l'endommagement s'inscrivent dans le cadre de la thermodynamique [2] et reposent sur des quantités macroscopiques. Cependant la plupart des phénomènes physiques (couplages, transferts de masses ...) tout comme les mécanismes de rupture ont pour origine des échelles plus fines. Ainsi, la possibilité d'établir un lien (numérique) entre les propriétés des constituants initiaux à ces échelles fines et les propriétés mécaniques résultantes à l'échelle macroscopique (par l'intermédiaire d'un modèle élément fini) est un problème classique. Les différentes techniques d'homogénéisation développées à cet effet répondent à ce type de problèmes dans le cas linéaire [3]. Pour les comportements non-linéaires les résultats obtenus sont nettement moins probants et l'outil numérique est nécessaire.

Dans la suite de ce papier, l'échelle la plus fine numériquement envisageable sera l'échelle dite « mésoscopique ». A cette échelle, les matériaux hétérogènes quasi-fragiles comme le béton peuvent être raisonnablement modélisés par deux ou trois phases distinctes et leur réponse mécanique macroscopique dépend à la fois des propriétés de chaque phase mais aussi de leur distribution spatiale. Etre capable de prendre en compte ces différentes phases et de palier aux problèmes de localisation numérique liés au comportement adoucissant apparaît comme nécessaire pour une modélisation réaliste. Le modèle proposé construit à partir d'un treillis tri-dimensionnel générés par une triangulation de Delaunay présente les atouts pour y parvenir. Afin de travailler avec des maillages non-adaptés, une discontinuité faible (discontinuité dans le champ de déformation) est
introduite au sein des éléments barre constituant le treillis. Trois types d'éléments barres apparaissent ainsi : ceux entièrement dans la matrice ou dans les inclusions et ceux « coupés » entre deux phases, présentant ainsi des propriétés élastiques différentes de part et d'autre. De plus, afin de modéliser la fissuration à cette échelle et l'adoucissement résultant sans dépendance au maillage [4], une discontinuité forte (discontinuité dans le champ de déplacement) est également introduite au sein des éléments barres. Elle permet de capturer à la fois les micro-fissures pouvant apparaître dans les différentes phases mais aussi la rupture à l'interface entre deux phases. La stratégie retenue pour enrichir la cinématique du problème élément fini [5] (introduction des deux types de discontinuités) est celle des Modes Incompatibles ([6], [7]). Cette méthode présente l'avantage d'être peu intrusive et de préserver le nombre d'inconnues total.

# 2 Modèle mésoscopique

## 2.1 Enrichissement de la cinématique

### 2.1.1 Discontinuité faible : matériaux hétérogènes

L'introduction d'une discontinuité faible donne la possibilité de travailler avec des maillages non-adaptés pour les matériaux hétérogènes. Ainsi un élément barre peut être « coupés » en deux avec des propriétés élastiques différentes de part et d'autre. D'un point de vue numérique, deux sous-domaines sont définis dans l'élément barre et une interface apparaît entre ces deux parties. C'est sur cette idée que repose le premier enrichissement de la cinématique, concrétisé par la fonction $G_1$

$$G_1 = \begin{cases} -\frac{1}{\theta\ell}, & x \in [0, \theta\ell] \\ \frac{1}{\ell(1-\theta)}, & x \in [\theta\ell, \ell] \end{cases}$$

### 2.1.2 Discontinuité forte : fissuration des matériaux quasi-fragiles

Le comportement adoucissant lié à l'apparition de fissures dans les matériaux quasi-fragiles entraîne une dépendance au maillage sur le plan numérique. Il a été montré par de nombreux auteurs [4] que l'introduction d'une discontinuité forte permet de capturer cette fissuration sans dépendance au maillage. C'est sur cette idée que repose le second enrichissement de la cinématique, concrétisé par la fonction $G_2$

$$G_2 = -\frac{1}{\ell} + \delta_{\Gamma}$$

## 2.2 Cadre théorique : la Méthode des Modes Incompatibles

La formulation variationnelle de Hu-Washizu [8] combinée à la méthode des déformations améliorées fournit un cadre théorique permettant d'enrichir la cinématique des éléments barres : le champ de déformations est repensé comme la somme du gradient symétrique des déplacements $\nabla^s u$ et des déformations améliorées $\tilde{\varepsilon}_1$ et $\tilde{\varepsilon}_2$.

$$\varepsilon = \nabla^s u + \tilde{\varepsilon}_1 + \tilde{\varepsilon}_2$$

En utilisant les fonctions de forme classiques pour un élément barre, les fonctions $G_1$ et $G_2$ introduites précédemment et en tenant compte du critère d'apparition de la discontinuité forte $\Phi$
d'un élément barre, le problème élément fini s'écrit :

$$
\left\{
\begin{array}{ll}
\underset{e=1}{^{n_{elm}}} [\mathbf{f}_e^{\text{int}}(d_e, \alpha_{1,e}, \alpha_{2,e}) - \mathbf{f}_e^{\text{ext}}] = 0 & \\
-\sigma_1(d_e, \alpha_{1,e}, \alpha_{2,e}) + \sigma_2(d_e, \alpha_{1,e}, \alpha_{2,e}) = 0 & \forall e \in [1, n_{elm}] \\
\theta \cdot \sigma_1(d_e, \alpha_{1,e}, \alpha_{2,e}) + (1-\theta) \cdot \sigma_2(d_e, \alpha_{1,e}, \alpha_{2,e}) - t_{\Gamma} = 0 & \forall e \in [1, n_{elm}] \\
\Phi(t_{\Gamma}, \alpha_{2,e}) = 0 & \forall e \in [1, n_{elm}]
\end{array}
\right.
$$

$\alpha_{1,e}$ et $\alpha_{2,e}$ sont les paramètres d'interpolation correspondant respectivement aux discontinu-
tés faible et forte. Ils sont définis au niveau élémentaire et sont donc des quantités locales. $\sigma_1$ et
$\sigma_2$ sont les contraintes dans chaque sous-domaine, $\theta$ est un scalaire qui paramètre la position de
l'interface au sein de la barre.

La fonction critère utilisée pour introduire la rupture d'un élément barre, c'est à dire la discon-
tinuité forte, s'écrit :

$$
\Phi = t_{\Gamma} - (\sigma_{f} - q)
$$

On note $t_Γ$ la contrainte sur la discontinuité et $\sigma_f$ la contrainte avant rupture. Il est important de
noter que ce critère n'est activé qu'en traction. Il n'y a donc pas de rupture en compression à
l'échelle mesoscopique.Le point clef est que la troisième équation du système fournit une valeur
explicite et directe de la contrainte sur la discontinuité $t_Γ$ :

$$
t_{\Gamma} = \theta \cdot \sigma_{1}(d_{e}, \alpha_{1,e}, \alpha_{2,e}) + (1 - \theta) \cdot \sigma_{2}(d_{e}, \alpha_{1,e}, \alpha_{2,e})
$$

De ce fait le système se réécrit :

$$
\left\{
\begin{array}{ll}
\underset{e=1}{^{n_{elm}}} [\mathbf{f}_e^{\text{int}}(d_e, \alpha_{1,e}, \alpha_{2,e}) - \mathbf{f}_e^{\text{ext}}] = 0 & \\
-\sigma_1(d_e, \alpha_{1,e}, \alpha_{2,e}) + \sigma_2(d_e, \alpha_{1,e}, \alpha_{2,e}) = 0 & \forall e \in [1, n_{elm}] \\
\Phi(t_\Gamma, \alpha_{2,e}) = 0 & \forall e \in [1, n_{elm}]
\end{array}
\right.
\qquad (1)
$$

La linéarisation de (1) par rapport à (d<sub>e</sub>, α<sub>1,e</sub>, α<sub>2,e</sub>) amène au système linéaire suivant :

$$
\left(
\begin{array}{ccc}
\mathbf{a} & \mathbf{b} & \mathbf{c} \\
\mathbf{b}^{\mathrm{T}} & d & e \\
-c^{\mathrm{T}} & -e & f(\alpha_{2})
\end{array}
\right)_{n+1}^{(k)}
\left(
\begin{array}{c}
\Delta d \\
\Delta \alpha_{1} \\
\Delta \alpha_{2}
\end{array}
\right)_{n+1}^{(k+1)}
=
\left(
\begin{array}{c}
-\underset{e=1}{^{n_{elm}}}[\mathbf{f}_e^{\text{int}}(d_e, \alpha_{1,e}, \alpha_{2,e}) - \mathbf{f}_e^{\text{ext}}]^{(k)} \\
\mathbf{h}_{1,e} \\
-\Phi
\end{array}
\right)_{n+1}^{(k)}
$$

L'obtention de (Δd, Δα₁ et Δα₂) se fait au travers de deux niveaux. Le niveau global, correspondant à la première équation de (2) et le niveau local. Dans l'esprit de la méthode de résolution partitionnée, Δα₁⁽ᵏ+¹⁾ et Δα₂⁽ᵏ+¹⁾ sont calculés pour une valeur dₙ+₁ fixée. Une fois cette étape terminée, les residus h₁ et Φ sont nuls et α₁ et α₂ sont à jour. On execute alors une double condensation statique sur Δα₁⁽ᵏ+¹⁾ et Δα₂⁽ᵏ+¹⁾, ceci permettant alors de remonter au niveau global. A ce niveau le système à résoudre s'écrit alors :

$$
\hat{\mathbf{K}}^{(k)} \Delta d^{(k+1)} = - \underset{e=1}{^{n_{elm}}} [\mathbf{f}_e^{\text{int}}(d_e, \alpha_{1,e}, \alpha_{2,e}) - \mathbf{f}_e^{\text{ext}}]^{(k)}
$$

où $\hat{\mathbf{K}}^{(k)}$ s'exprime de la façon suivante :

$$
\hat{\mathbf{K}}^{(k)} = \mathbf{a} - \frac{1}{d.f(\alpha_2) + e^2} \cdot \left( \mathbf{b} \ \mathbf{c} \right) \cdot \left( \begin{array}{cc} f(\alpha_2) & -e \\ e & d \end{array} \right) \cdot \left( \begin{array}{c} \mathbf{b}^{\mathrm{T}} \\ -\mathbf{c}^{\mathrm{T}} \end{array} \right)
$$
Ce méso modèle permet d'obtenir des réponses macroscopiques complexes suivant différents trajets de chargement (compression, cisaillement, bi- tri-compression). Néanmoins, avec une approche mono-échelle et monolithic, il est difficilement envisageable de mener des calculs sur des domaines de taille importante. Dans la suite, nous présentons une approche multi-échelles intégrées à même de s'attaquer à des problèmes plus importants.

## 3 Méthode élément finis multi-échelles

Afin de pouvoir employer le méso modèle précédent au cours d'un calcul de structure, une méthode multi-échelles intégrée – suivant l'acceptation introduite par [9] – permettant un calcul simultanée de la structure et de sa méso structure est envisagée. Chaque élément intervenant dans la discrétisation macroscopique peut accueillir une sous description, plus fine, elle même traitée par la méthode des éléments finis. La compatibilité cinématique entre les deux niveaux de représentation – les déplacements micro $u^m$ et macro $u^M$ – est assurée par des multiplicateurs de Lagrange localisés $\lambda$ [10], comme illustré à la Figure 1. Cette stratégie, originalement développée par [11],

<figure><img src="image_2.png" /><figcaption>Figure 1 – Approche multi-échelles intégrée</figcaption></figure>

s'appuie sur une formulation variationelle dans laquelle l'énergie interne macroscopique est calculée en sommant les contributions microscopiques sur chaque sous domaine microscopique $\Omega_e$ :

$$ \Pi^{\text{int}} = \sum_{e} \left( \int_{\Omega_e} \Psi_e(\varepsilon_e^m) \, d\Omega - \int_{\partial\Omega_e} \lambda_e \cdot (u_e^m - u_e^M) \, dS \right) \quad (2) $$

Il découle de cette formulation une équation aux résidus qui de la même manière que pour les modes incompatibles, enrichissant la cinématique de chaque élément barre, peut être explicitement séparé en un résidu macroscopique $r^M$ (3a) et autant de résidus, $r_e^S$ (3b), locaux à chaque sous-domaine possédant une description microscopique.

$$ r^M = \sum_{e=1}^{n_e} F_e^T \cdot \Lambda_e - f^{\text{ext}} \quad (3a) $$

$$ r_e^S = \begin{bmatrix} G_e^T \cdot \Lambda_e \\ E_e^T \cdot U_e^m + G_e \cdot \mu_e - F_e \cdot U_e^M \\ f_e^{\text{int}} \cdot U_e^m - E_e \cdot \Lambda_e \end{bmatrix} \quad (3b) $$

La résolution des déplacements macro discrétisés, $U^M$, est ainsi dissociée de la résolution des grandeurs discrétisées locales sur chaque sous-domaine : les déplacements micro $U_e^m$, les multiplicateurs de Lagrange localisés $\Lambda_e$ et les amplitudes des mouvements de corps rigides $\mu_e$. Ces derniers apparaissent naturellement au cours de la résolution de chaque problème micro. En effet, ces problèmes ne possèdent que des conditions à effort imposé sur leur frontière et ne peuvent
être résolus qu'à un mouvement de corps rigide près qui doit être déterminé. Les matrices $E_e$, $G_e$ et $F_e$ sont des matrices de projection respectivement des déplacements micro, des amplitudes de mouvement de corps rigide et des déplacements macro sur les frontières des sous-domaines micro.

Par une double condensation statique, le système obtenu en employant un schéma de résolution de Newton à (3) peut être partitionné en trois sous-systèmes. Chacun d'eux est dédié à un niveau de résolution : microscopique, aux frontières de chaque sous-domaine micro et macroscopique. La mise en oeuvre de cette stratégie est réalisée dans le cadre du paradigme orienté composant. Chaque échelle de résolution est traitée par un code dédié. Ceux-ci ont la possibilité d'interagir, de communiquer entre eux en réseaux et d'être exécuter en parallèle. L'interaction des composants au sein de l'architecture MuSCAd [12] permet la résolution de problèmes multi-échelles possédant le meso modèle (section 2) comme niveau fin de description. Un exemple est présenté dans la prochaine section.

# 4 Réponse macroscopique en multi-échelles intégrées

L'architecture MuSCAd permet de réaliser des calculs multi-échelles, ici sur une éprouvette de béton de 10 cm de côté, où chaque élément macro possède une description mésoscopique (Fig.2(a)) présentée à la section 2. Les inclusions sont plus rigide que la matrice et élastiques. La matrice est considérée comme quasi-fragile.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Calcul multi-échelles jusqu'à rupture</figcaption></figure>

Grâce à une méthode de longueur d'arc particulière implementée au niveau microscopique [12], il est possible de réaliser un calcul multi-échelles jusqu'à rupture (Fig.2(b)). Il est donc envisageable, en utilisant l'approche proposée, d'intégrer à un calcul de structure une information provenant d'une échelle plus fine où se produisent effectivement les premiers phénomènes de localisation et d'en capturer les effets à une échelle supérieure.

# 5 Conclusion

Cette commnication a présenté un méso modèle original pour les matériaux à matrice cimenteries permettant de prendre en compte à la fois leur caractère hétérogène et leur comportement quasi-fragile. Il permet une modélisation fine des phénomènes complexes qui interviennent au cours des différentes phases du comportement de ces matériaux, y compris en présence d'adoucissement. Pour permettre un calcul de structure faisant intervenir ce modèle à une échelle plus
fine, une méthode multi-échelles intégrée a été également proposée. L'architecture MuSCAd rend
possible le calcul distribué et simultané aux différentes échelles prises en compte et montre des
bons résultats sur des cas de structures simples à méso structures complexes.

Remerciements

Le financement de ces travaux est en parti assuré par l'Université Franco-Allemande (UFA-
DFH).

Références

[1] M. Ostoja-Starzewski. Material spatial randomness : From statistical to representative volume element. *Probabilistic Engineering Mechanics*, 21 :112-132, 2006.

[2] J. Lemaitre et J.L Chaboche. Mécanique des matériaux solides. *Dunod*, 1988.

[3] Z. Hashin and S. Shtrikman. A variational approach to the theory of the elastic behaviour of multiphase materials. *Journal of the mechanics and physics of solids*, 11 :127–140, 1963.

[4] J.C Simo, J. Oliver, and F. Armero. An analysis of strong discontinuity induced by strain softening solutions in rate-independent solids. *Journal of Computational Mechanics*, 12 :277–296, 1993.

[5] G.N. Wells and L.J. Sluys. Application of embedded discontinuities for softening solids. *Engineering fracture mechanics*, 65 :263–281, 2000.

[6] J.C Simo and M. Rifai. A class of mixed assumed strain methods and the method of incompatible modes. *International Journal for Numerical Methods in Engineering*, 29 :1595–1638, 1990.

[7] A. Ibrahimbegovic and E.L Wilson. A modified method of incompatible modes. *Communications in Applied Numerical Methods*, 7 :187–194, 1991.

[8] K. Washizu. Variational methods in elasticity and plasticity. *Pergamon Press, New York*, 3 edition, 1982.

[9] F. Feyel and J.-L. Chaboche. Multi-scale non-linear FE² analysis of composite structures : damage and fiber size effects. *Revue européenne des Éléments Finis : NUMDAM'00 issue*, 10 :449–472, 2001.

[10] K. C. Park, C. A. Felippa, and U. A. Gumaste. A localized version of the method of lagrange multipliers and its applications. *Computational Mechanics*, 24 :476–490, 2000.

[11] A. Ibrahimbegovic and D. Markovic. Strong coupling methods in multi-phase and multi-scale modeling of inelastic behavior of heterogeneous structures. *Computer methods in applied mechanics and engineering*, 192 :3089–3108, 2003.

[12] Martin Hautefeuille. *Numerical Modeling strategy for Heterogeneous Materials : A FE Multi-scale and Component-based Approach*. PhD thesis, Université Technologique de Compiègne – Technische Universität Braunshweig, en cours.