# Décomposition de domaine pour l'étude double-échelle de milieux hétérogènes

Willy Leclerc, Philippe Karamian

► To cite this version:

Willy Leclerc, Philippe Karamian. Décomposition de domaine pour l'étude double-échelle de milieux hétérogènes. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717082

HAL Id: hal-01717082

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

Submitted on 25 Feb 2018

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
# Décomposition de domaine pour l'étude double-échelle de milieux hétérogènes

Willy LECLERC<fnref n="1" />, Philippe KARAMIAN<fnref n="1" />

**Résumé** — La communication porte sur la mise en place de la décomposition de domaine pour l'étude double-échelle de milieux hétérogènes. Nous nous intéressons à l'étude des propriétés élastiques effectives de composites dont les renforts sont fortement enchevêtrés et forment un réseau géométrique complexe. Dans un tel contexte, une étude éléments finis basée sur une modélisation numérique viable s'avère difficile à concevoir sans un coût de calcul prohibitif. L'objectif consiste à accélerer l'estimation des propriétés grâce à l'outil de décomposition de domaine. Dans ce but, des adaptations des méthodes du complément de Schur et FETI au cadre de l'homogénéisation double-échelle sont proposées.

**Mots clés** — décomposition de domaine, éléments finis, volume élémentaire représentatif, matériau composite, calcul parallèle.

## 1 Introduction

La microstructure des composites renforcés en fibres aléatoirement réparties présente un fort enchevètement des hétérogénéités. La géométrie du réseau de fibres en est alors complexe et difficile à appréhender dans un cadre éléments finis. Les maillages sont en effet difficile à générer dans un tel contexte, et s'avèrent souvent mal structurés. Un tel inconvénient de conception devient rédhibitoire dans un cadre automatisé pour lequel un échantillonnage entier de motifs représentatifs du matériau doit être mis en place. C'est pourquoi, nous avons opté pour une conception de maillages par pixelisation en 2D ou voxélisation en 3D [1, 2, 3]. Le concept est basé sur une approximation de la géométrie à l'aide d'une grille d'éléments quadrangulaires ou hexaédriques. Il permet de générer des microstructures complexes avec une grande fiabilité. Cependant, une telle approche est particulièrement coûteuse lorsque la finesse de la microstructure est requise ou que le motif représentatif à concevoir est de grandes dimensions. Différentes techniques s'avèrent efficaces pour améliorer le concept parmi lesquelles nous retrouvons les méthodes de lissage [4]. Nous avons choisi le modèle d'approximation d'ordre $n$ basé sur le concept de raffinement de maillage adaptatif [5, 6]. Celui-ci se fonde sur une génération en deux temps du réseau d'hétérogénéités. Tout d'abord les renforts sont approximés selon une grille grossière d'éléments de largeur égale à une dimension caractéristique de l'hétérogénéité considérée. Dans une deuxième étape, un processus de raffinement local permet d'affiner la géométrie en zone utile. L'approche est efficace pour des motifs pour lesquels le rapport d'échelle entre celui-ci et le renfort est faible (inférieure à 10) mais devient inefficace au-delà.

La notion de VER est d'une grande importance dans le domaine des milieux hétérogènes et nécessite une détermination précise laquelle dépend de la configuration matérielle [7, 8]. Ainsi, certains types de matériaux exigent un motif aux larges dimensions pour que celui-ci soit considéré comme VER. Dans ce contexte, l'approche par voxélisation via le modèle d'approximation d'ordre $n$ perd toute son efficacité. Afin de palier au problème, nous considérons l'outil de décomposition de domaine. L'idée est de subdiviser le VER initial en motifs de plus petites tailles pour lesquels l'approche considérée demeure efficace. Dans l'optique du calcul des propriétés effectives du matériau par le concept double-échelle et périodique décrit par [9], les méthodes classiques que sont le complément de Schur [10] et FETI [11] doivent être légèrement modifiées. Il s'agit d'une part de faire apparaître des termes liés aux coefficients homogénéisés, d'autre part de prendre en compte l'hypothèse de périodicité. Le présent papier a pour but de présenter la mise en place de l'outil de décomposition de domaine dans le cadre spécifique des

<footnote n="1">¹ LMNO, Université de Caen Basse-Normandie, France, willy.leclerc@unicaen.fr</footnote>
milieux hétérogènes fibreux et du calcul de leurs propriétés effectives par la méthode double-échelle. Il est structuré comme suit, tout d'abord nous décrivons la modélisation numérique par l'approche par voxelisation. Dans un second temps, nous décrivons la méthode d'évaluation des propriétés. Enfin, dans un troisième et dernier temps, nous nous intéressons aux méthodes de décomposition de domaine et aux modifications réalisées pour les adapter au cadre double-échelle.

# 2 Modélisation numérique

## 2.1 Construction de la géométrie

Nous générons des cellules unitaires du composite renforcé en fibres courtes en suivant plusieurs étapes. Tout d'abord, nous considérons des motifs carrés ou cubiques selon la dimension d'étude choisie (2D ou 3D). Les dimensions du motif représentatif dépendent du rapport d'échelle souhaité entre ces dernières et la longueur moyenne de chaque hétérogénéité. Il est important de savoir à ce stade que ce choix est crucial afin de garantir la représentativité du VER final [8]. Nous générons par la suite une microstructure selon un choix prédéfini de paramètres morphologiques, à savoir la longueur, le diamètre, l'orientation ou encore la courbure et la tortuosité de chaque fibre. Il s'agit de tirer aléatoirement chacune de ces variables selon des lois de distributions données. Cependant, la mise en place d'une telle conception s'avère délicate si l'interpénétrabilité des fibres n'est pas autorisée. Ainsi, les méthodes d'adsorption séquentielle aléatoire [12, 13, 14, 15], lesquelles se basent sur une telle hypothèse, exhibent un coût de calcul parfois prohibitif, en particulier pour des fractions volumiques d'hétérogénéités élevées. Dans le présent travail, nous avons préféré concevoir des microstructures pour lesquelles le recouvrement des fibres est possible. L'avantage de cette approche est double, d'une part elle permet de conserver l'efficacité d'une génération aléatoire. D'autre part, la distribution de chaque paramètre morphologique est ainsi respectée sans génération de biais. Dans l'hypothèse où la courbure de la fibre n'est pas prise en compte, nous générons chaque fibre par un cylindre parfait en 3D et une coupe de cylindre dans le plan en 2D. Enfin, chaque hétérogénéité est reproduite par périodicité si elle est en contact avec un ou plusieurs bords en 2D ou faces en 3D du motif unitaire.

## 2.2 Modèle d'approximation d'ordre n et partitionnement

La génération de maillages sur une géométrie complexe avec recouvrement possible des hétérogénéités présente plusieurs inconvénients. Ainsi, si une reconstruction des frontières du réseau de fibres est toujours possible, l'opération peut s'avérer particulièrement coûteuse. D'autre part, les maillages triangulés autour de ces mêmes interfaces sont souvent peu structurés et leurs défauts impactent directement sur la résolution éléments finis. Nous avons préféré opter pour une approche par pixelisation ou voxelisation du domaine. Il s'agit d'approximer la géométrie réelle du réseau d'hétérogénéités selon une grille d'éléments quadrangulaires en 2D ou hexaédriques en 3D [1, 2, 3]. Un tel procédé est parfaitement fiable et automatisable pour n'importe quel type de microstructures. Cependant, il s'avère prohibitif lorsque la géométrie exigée est fine et nécessite un nombre très important d'éléments quadrangulaires ou hexaédriques. Différents concepts permettent d'améliorer l'efficacité de cette approche. Notre choix s'est porté sur le modèle d'approximation d'ordre n décrit par [5, 6] basé sur la notion de raffinement de maillage adaptatif (RMA). Dans cette approche, les éléments sont tout d'abord choisis de taille équivalente au diamètre de chaque fibre puis localement raffinés selon une puissance de 2 afin d'approcher au mieux la géométrie de la microstructure. La puissance n est appelée ordre d'approximation et est définie comme suit,

$$n = \log_2 \left( \frac{\text{Diamètre de la fibre}}{\text{Largeur du voxel}} \right) \quad (1)$$

Le modèle peut être amélioré en vue de privilégier le gain de temps à la précision géométrique. De plus, un tel concept est parfaitement fiable et adapté à une subdivision régulière du domaine en vue de calculs parallélisés. La Figure 1 illustre une microstructure 3D (a) générée aléatoirement et composée de 40 fibres de longueur 0,2. Le réseau est approximé à l'aide du modèle d'approximation d'ordre 0 (b), c'est à dire le concept pour lequel la taille des voxels est égale au diamètre des fibres. Nous observons
que malgré l'aspect irrégulier des hétérogénéités, les paramètres morphologiques sont tous respectés et permettent d'obtenir un réseau visuellement proche de celui initialement mis en place.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Exemple de VER dont la géométrie initiale composée de 40 hétérogénéités est approximée à l'ordre 0.</figcaption></figure>

## 2.3 Efficacité de la méthode

L'efficacité du modèle d'approximation géométrique est étudiée pour différentes valeurs de l'ordre $n$. Une comparaison est effectuée pour les premiers ordres ($n=0, 1$ et $2$) avec une voxelisation classique sans RMA à taille de voxel identique. Nous considérons des VERs pour lesquels les fibres sont aléatoirement distribuées et orientées, et de longueur 0,2. Les simulations numériques sont effectuées à l'aide d'un processeur Intel(R) Xeon(R) W3670 @ 3.20 GHz. Il est remarquable que le raffinement local est particulièrement efficace lorsque l'ordre de la méthode est importante et mène à une géométrie fine et bien approximée. Plus la densité de fibres est importante et plus l'apport en efficacité est réel. Cependant, dans le cadre des tests ici effectués, si l'ordre est égal à 0, l'approche par RMA n'est guère meilleure voire pire que celle sans RMA. Cette observation est liée aux dimensions des éléments utilisés lesquelles génèrent alors un faible stockage mémoire.

<table><caption>Tableau 1 – Comparaison en temps de calcul CPU de la méthode de voxelisation avec et sans RMA</caption><thead><tr><th>Modèle</th><th>1 fibre</th><th>25 fibres</th><th>100 fibres</th><th>500 fibres</th></tr></thead><tbody><tr><td>Voxelisation classique (ordre 0)</td><td>00h00'03"7</td><td>00h01'14"8</td><td>00h05'05"9</td><td>00h26'22"3</td></tr><tr><td>Voxelisation classique (ordre 1)</td><td>00h00'39"4</td><td>00h10'36"9</td><td>00h40'16"8</td><td>03h33'01"5</td></tr><tr><td>Voxelisation classique (ordre 2)</td><td>00h28'56"4</td><td>10h31'00"7</td><td>&gt;1 jour</td><td>—</td></tr><tr><td>Voxelisation avec RMA (ordre 0)</td><td>00h00'03"8</td><td>00h01'15"2</td><td>00h05'04"9</td><td>00h31'19"2</td></tr><tr><td>Voxelisation avec RMA (ordre 1)</td><td>00h00'05"5</td><td>00h01'58"4</td><td>00h08'20"1</td><td>01h01'01"4</td></tr><tr><td>Voxelisation avec RMA (ordre 2)</td><td>00h00'13"9</td><td>00h05'49"8</td><td>00h25'23"4</td><td>03h20'32"2</td></tr></tbody></table>

# 3 Evaluation des propriétés

## 3.1 Homogénéisation double-échelle

Les propriétés élastiques du composite sont estimées à l'aide de la méthode d'homogénéisation multi-échelle introduite par [9]. Il s'agit d'évaluer par éléments finis le tenseur de complaisance du milieu
hétérogène modélisé par un VER. L'approche consiste à prendre en compte l'influence d'hétérogénéités microscopiques voire nanoscopiques à l'échelle macroscopique. Elle s'appuie sur des outils asymptotiques et une hypothèse de périodicité à l'échelle des inclusions. Dans le présent travail, nous nous limitons à deux échelles, l'échelle microscopique y liée au VER et l'échelle macroscopique x liée au matériau. Dans cette approche, l'équation fondamentale décrivant la loi de comportement du matériau est développée selon un paramètre ε lequel décrit le rapport d'échelle entre le composite et son motif représentatif. Ainsi, dans l'hypothèse d'un milieu continu élastique linéaire, nous avons,

$$ egin{cases} -\operatorname{div} \sigma(\mathbf{u}^{\varepsilon}) = \mathbf{f} & (+C.L.) \\ \sigma_{ij}(\mathbf{u}^{\varepsilon}) = C_{ijkl}^{\varepsilon} e_{kl}(\mathbf{u}^{\varepsilon}) \end{cases} \qquad (2) $$

où $\sigma$ est le tenseur des contraintes, $e$ le tenseur de déformations, $C_{ijkl}^{\varepsilon}$ le tenseur de rigidité local et $\mathbf{f}$ un chargement volumique. En outre, les initiales $C.L.$ font référence aux conditions aux limites. Le déplacement $\mathbf{u}^{\varepsilon}$ dépend de $\varepsilon$ comme suit,

$$ \mathbf{u}^{\varepsilon}(x,y) = \mathbf{u}^{0}(x) + \varepsilon\mathbf{u}^{1}(x,y) + \varepsilon^{2}\mathbf{u}^{2}(x,y) + o(\varepsilon^{2}) \qquad (3) $$

Deux nouvelles formulations de l'équation 2 sont obtenues après injection du développement de $\mathbf{u}^{\varepsilon}(x,y)$ dans celle-ci. La première dite microscopique correspond aux termes $\varepsilon^{-1}$, et la seconde dite macroscopique correspond aux termes $\varepsilon^{0}$. Une étude variationnelle effectuée sur la première nous permet alors de réécrire la seconde sous une forme analogue à l'équation 3. D'où, par analogie, nous tirons ensuite une expression du tenseur de complaisance homogénéisé $C_{ijkl}^{hom}$ comme suit,

$$ C_{ijkl}^{hom} = \frac{1}{|Y|} \int_{Y} C_{ijpq}(y) \left[ \delta_{p}^{k} \delta_{q}^{l} + e_{pqy} (\kappa^{kl}(y)) \right] dY \qquad (4) $$

où $Y$ désigne le VER et $\kappa^{kl}(y)$ est une fonction $Y$-périodique dépendant uniquement de $y$ solution du problème local. Ainsi, la formulation des propriétés homogénéisées est la somme de la moyenne intégrale des propriétés et d'un terme correctif lié à l'influence de la microstructure à l'échelle macroscopique.

## 3.2 Discrétisation et calcul éléments finis

Le tenseur de complaisance homogénéisé $C_{ijkl}^{hom}$ est obtenu par résolution d'une nouvelle formulation variationnelle dont la discrétisation par la méthode de Galerkin mène aux systèmes matriciels suivants :

$$ \underbrace{\begin{bmatrix} K & KE^t \\ KE & E \end{bmatrix}}_{\mathcal{K}} \underbrace{\begin{bmatrix} \omega^{kl} \\ P_{ij}^{kl} \end{bmatrix}}_{\mathbf{u}^{kl}} = \underbrace{\begin{bmatrix} \mathbf{0}^{kl} \\ \delta_i^k \delta_j^l \end{bmatrix}}_{\mathbf{f}^{kl}} \qquad (5) $$

où $K$ est la matrice de rigidité, $E$ la moyenne intégrale du tenseur de complaisance sans terme correctif et $KE$ une matrice de couplage entre $K$ et $E$. $\omega^{kl}$ est le déplacement local à l'échelle de la microstructure. La résolution s'effectue pour les 6 $kl$ problèmes et mène à une estimation directe de $P_{ij}^{kl}$ et $\omega^{kl}$. On montre [9] que pour un choix judicieux du second membre, $P_{ij}^{kl}$ n'est autre que le tenseur de souplesse homogénéisé, en d'autres termes, l'inverse de $C_{ijkl}^{hom}$.

## 3.3 Evaluation asymptotique

La taille du VER est une question cruciale de laquelle dépend la représentativité du motif élémentaire et le coût de calcul nécessaire à l'obtention des propriétés homogénéisées. Ainsi, si un VER de grande taille permet d'assurer la question de la représentativité au détriment de l'efficacité, un petit VER permet lui une génération efficace au détriment de la fiabilité. Notre choix est d'évaluer asymptotiquement les propriétés en effectuant un tirage aléatoire de motifs élémentaires. Le processus allie efficacité et fiabilité mais nécessite une étude préalable de la taille optimale. Cette dernière est effectuée à l'aide d'outils statistiques tels la portée intégrale [8]. D'autre part, le nombre de réalisations $n_r$ à effectuer peut être évalué pour une erreur donnée par simple connaissance du paramètre de variance [5, 8, 16]. En effet,
l'erreur absolue $\epsilon_{abs}$ commise sur une propriété E est liée au volume du VER, à l'écart-type de E, $D_E(V)$ et à $n_r$ comme suit,

$$ \epsilon_{abs} = \frac{1.96D_E(V)}{\sqrt{n_r}} \qquad (6) $$

Les propriétés effectives sont alors estimées de la façon suivante. Tout d'abord nous estimons le tenseur de complaisance pour chaque VER d'un échantillonnage donné de motifs représentatifs par la méthode double-échelle décrite précédemment. Puis, dans un second temps, nous évaluons l'espérance du champ des résultats lequel tend asymptotiquement vers le tenseur de complaisance homogénéisé.

# 4 Décomposition de domaine

## 4.1 Subdivision du domaine

Chaque motif représentatif peut être partitionné dans l'optique d'une mise en place parallélisée à la fois de la génération du maillage et de la résolution éléments finis. La subdivision s'effectue automatiquement après la phase de génération de la microstructure et mène à un partitionnement régulier du domaine initial. Un traitement spécifique est alors réalisé afin d'assurer à la fois la périodicité des bords et faces extérieurs et la continuité aux interfaces internes entre sous-domaines. La Figure 2 montre l'exemple d'un VER 3D composé d'une quarantaine d'hétérogénéités et subdivisé en 8 sous-domaines de mêmes dimensions. Il est remarquable que chaque inclusion traversant une ou plusieurs faces intérieures est prise en compte par l'ensemble des sous-domaines visités. Dans l'ensemble de la présente section, nous notons $\Gamma^i$ l'ensemble des noeuds situés sur les frontières intérieures, $\Gamma^e$ l'ensemble des noeuds situés sur les frontières extérieures et $\Gamma = \Gamma^i \cap \Gamma^e$ la réunion des deux précédents ensembles. $\Omega_n$ est le domaine correspondant à la subdivision numéro n du domaine initial $\Omega$.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Exemple de partitionnement à 8 sous-domaines d'un VER 3D</figcaption></figure>

## 4.2 Méthode du complément de Schur

Dans une première approche, nous nous intéressons à la méthode du complément de Schur [10]. Le but est de mettre en place une méthode facilement implémentable et de la modifier afin d'évaluer directement les propriétés du matériau selon la technique double-échelle vue en sous-section 3.1. Le processus usuel consiste à réduire le problème initial à un problème d'interface prenant en compte les frontières intérieures mais aussi extérieures du fait de l'hypothèse de périodicité du VER. La modification principale se résume alors à faire intervenir des termes supplémentaires liés à l'homogénéisation. Le processus initial se compose d'une renumérotation du système matriciel 5. Ainsi, dans l'hypothèse d'une subdivision en N sous-domaines, la matrice renumérotée $\tilde{\mathcal{K}}$ de $\mathcal{K}$ est comme suit,
$$
\tilde{K} = \begin{bmatrix}
1 & \dots & i & \dots & N & \Gamma & C \\
\hline
1 & K_1 & 0 & \dots & \dots & K_{\Gamma 1}^t & KE_1^t \\
\vdots & 0 & \ddots & \ddots & \ddots & \vdots & \vdots \\
i & \vdots & \ddots & K_i & \ddots & K_{\Gamma i}^t & KE_i^t \\
\vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\
N & 0 & \dots & \dots & 0 & K_N & KE_N^t \\
\Gamma & K_{\Gamma 1} & \dots & K_{\Gamma i} & \dots & K_{\Gamma N} & KE_{\Gamma}^t \\
C & KE_1 & \dots & KE_i & \dots & KE_N & E
\end{bmatrix} \quad (7)
$$

La matrice de rigidité *K* est subdivisée en *N* matrices $K_i$ où *i* est le numéro du sous-domaine. Il faut
garder à l'esprit que chaque matrice $K_i$ est seulement composée des contributions strictement intérieures
au domaine $\Omega_i$ correspondant. En effet, les contributions aux interfaces au domaine *i* sont rassemblées
dans une matrice $K_{\Gamma i}$ pour laquelle nous considérons indistinctement noeuds aux interfaces extérieures
et intérieures. De plus, nous supposons que chaque noeud à l'interface constitue une unique contribution
et ce quelque soit le nombre de sous-domaines auxquels celui-ci appartient. $KE_i$ représente une matrice
de couplage analogue à *KE* sur le sous-domaine *i*. La moyenne intégrale $E_i$ du tenseur de complaisance
*C* sur le sous domaine *i* n'est ici pas directement visible. Elle apparaît dans la matrice *E*, comme étant la
moyenne intégrale effectuée sur l'ensemble du VER, du fait que :

$$
E_i = \frac{1}{|\Omega_i|} \int_{\Omega_i} C dY \qquad (8)
$$

et,

$$
E = \sum_i \frac{1}{|\Omega_i|} \int_{\Omega_i} C dY = \sum_i E_i \quad (9)
$$

D'autre part, les second membre $\tilde{\mathbf{f}}^{kl}$ et solution $\tilde{\mathbf{u}}^{kl}$ sont définis par analogie avec $\mathbf{f}^{kl}$ et $\mathbf{u}^{kl}$ tels que
vus en section 3.2. Ainsi, nous obtenons un système matriciel renuméroté $\tilde{\mathcal{K}}\tilde{\mathbf{u}}^{kl} = \tilde{\mathbf{f}}^{kl}$ dont la résolution
s'effectue via le problème d'interface suivant,

$$
\left\{
\begin{array}{l}
\tilde{S} \mathbf{u}_{\Gamma} = \tilde{\mathbf{f}}_{\Gamma}^{kl} \\
\tilde{S} = \tilde{K}_{\Gamma} - \sum_{i} \tilde{K}_{\Gamma i} K_{i}^{-1} \tilde{K}_{\Gamma i}^{t}
\end{array}
\right.
\qquad (10)
$$

où $\tilde{S}$ est une matrice très proche du complément de Schur et dont la particularité est de prendre en compte
des termes supplémentaires liés aux coefficients homogénéisés. Les expressions de $\tilde{K}_{\Gamma i}$, $\tilde{K}_{\Gamma}$, $\tilde{\mathbf{f}}_{\Gamma}^{kl}$ et $\mathbf{u}_{\Gamma}$ sont
comme suit,

$$
\tilde{K}_{\Gamma i} = \begin{bmatrix} K_{\Gamma i} \\ KE_i \end{bmatrix} \qquad \tilde{K}_{\Gamma} = \begin{bmatrix} K_{\Gamma} & KE_{\Gamma}^t \\ KE_{\Gamma} & E \end{bmatrix} \qquad \tilde{\mathbf{f}}_{\Gamma}^{kl} = \begin{bmatrix} \mathbf{0}^{kl} \\ \delta_i^k \delta_j^l \end{bmatrix} \qquad \mathbf{u}_{\Gamma} = \begin{bmatrix} \omega_{\Gamma}^{kl} \\ P_{ij}^{kl} \end{bmatrix} \qquad (11)
$$

La résolution du problème d'interface 10 permet alors une obtention directe de la matrice de complai-
sance homogénéisée. L'avantage de la méthode réside dans la dimension finale du problème à résoudre
laquelle se limite aux contributions aux interfaces et aux 6 termes liés à la technique d'homogénéisation.
Ainsi, le système matriciel $\tilde{\mathcal{K}}\tilde{\mathbf{u}}^{kl} = \tilde{\mathbf{f}}^{kl}$ n'est que partiellement résolu sans que cela n'impacte l'évaluation
des propriétés élastiques.

4.3 Méthode mixte FETI 1 - complément de Schur

L'approche proposée dans le cadre de la méthode du complément de Schur est dite primale. Le terme
s'explique par la manière dont sont prises en compte à la fois les conditions de continuité et celles
de périodicité ; c'est à dire, en associant à chaque noeud d'interface la sommation des contributions de chaque sous-domaine auquel il appartient. Une seconde approche dite duale consiste à prendre en compte les conditions aux limites sous forme de multiplicateurs de Lagrange lesquels peuvent être vus comme les réactions aux noeuds d'interfaces. La méthode FETI mise en place par [11] fournit un cadre à la fois théorique et numérique au concept dual. La méthode consiste à travailler dans l'espace des multiplicateurs de Lagrange contraint par un second niveau de multiplicateurs prenant en compte les modes de corps rigides liés à l'existence de domaines flottants. Cependant, dans la présente approche où la périodicité induit le caractère flottant de l'ensemble des sous-domaines, nous avons préféré opter pour un concept mixte primal/dual. Il s'agit de prendre en compte les conditions de continuité sous forme duale et les conditions aux limites périodiques sous forme primale. Enfin, identiquement à la méthode du complément de Schur, tout noeud situé à l'interface est pris en compte une unique fois quel que soit le nombre de sous-domaines auquel il appartient. Ainsi, le problème d'interface FETI, lié à une renumérotation du système 5 analogue à celle vue en sous-section 4.2, est comme suit,

$$ \begin{bmatrix} F_I & G_I \\ G'_I & \emptyset \end{bmatrix} \begin{bmatrix} \lambda^{kl} \\ \alpha^{kl} \end{bmatrix} = \begin{bmatrix} d^{kl} \\ 0^{kl} \end{bmatrix} \qquad (12) $$

où $F_I$ est le complément de Schur associé au problème dual et prenant en compte les connections aux interfaces et les coefficients homogénéisés. $G_I$ décrit les traces des modes de corps rigides aux interfaces. Enfin, le second membre $d^{kl}$ se réduit à une matrice identité prolongée de 0 du fait que nous ne considérons aucune force additionnelle. La résolution s'effectue par un gradient conjugué préconditionné projeté sur l'espace des multiplicateurs ou par un algorithme de type Uzawa. La solution $\lambda^{kl}$ fournit alors à la fois les multiplicateurs de Lagrange liés aux conditions de périodicité, les déplacements locaux liés aux conditions de continuité et la matrice de souplesse homogénéisée.

# 5 Validation numérique

Nous générons un échantillonnage de VERs 2D pour lesquels les fibres sont aléatoirement distribuées et orientées. Nous fixons la longueur et le diamètre de chaque hétérogénéité respectivement à 1/12 et 1/240. Chaque inclusion est droite sans courbure dans le plan et le recouvrement entre les renforts est autorisé. Nous considérons, en outre, une subdivision régulière et sans recouvrement de chaque VER en sous-domaines de même forme. La densité de fibres est tirée aléatoirement entre 0 et 30 fibres par cellule unitaire ce qui correspond à une fraction surfacique inférieure à 25%. Le maillage est quant à lui généré selon le modèle d'approximation géométrique d'ordre 0. En d'autres termes, nous considérons le modèle pour lequel les éléments sont de taille égale au diamètre des hétérogénéités. L'évaluation des propriétés homogénéisées est effectuée à l'aide de la méthode du complément de Schur présentée en sous section 4.2 par la méthode double-échelle présentée en sous-section 3.1. Chaque fibre est supposée suivre une loi de comportement isotrope transverse. Les moduli de Young longitudinal et transversal sont fixés respectivement à 1050 GPa et 600 GPa. Le module de cisaillement est fixé à 450 GPa. La matrice est une résine époxy de moduli de Young et de cisaillement fixés respectivement à 4,2 et 1,55 GPa. Afin de vérifier la fiabilité de la méthode du complément de Schur dans le présent contexte des composites à fibres courtes aléatoirement dispersées, nous considérons différentes subdivisions en 4, 9, 16 et 36 sous-domaines. Les Figures 3a) et b) illustrent l'influence de la densité de fibres sur les moduli de Young et de cisaillement effectifs du milieu hétérogène supposé isotrope. Il est important de préciser que le maillage de chaque sous-domaine est obtenu après un découpage du maillage du VER initial, de sorte que le nombre total de degrés de liberté soit conservé. Un comparatif des résultats avec ceux obtenus selon le même processus sans décomposition de domaine montre une bonne fiabilité du modèle. Un certain écart est cependant remarquable à forte densité d'hétérogénéités où l'erreur relative atteint une valeur de près de 3% à densité de 30 fibres par cellule unitaire pour le cas du module de Young. Il a été vérifié que l'écart est lié au traitement aux interfaces garantissant la continuité du milieu entre les sous-domaines et non à la méthode mise en place.
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Influence de la densité de fibres sur a) le module de Young homogénéisé et b) le module de cisaillement homogénéisé pour différentes subdivisions des VERs</figcaption></figure>

Références

[1] L. Mishnaevsky Jr, K. Derrien and D. Baptiste *Effect of Microstructure of Particle Reinforced Composites on the Damage Evolution: Probabilistic and Numerical Analysis*, Composites Science and Technology 64, 1805-1818, 2004.

[2] L. Mishnaevsky Jr *Three-dimensional Numerical Testing of Microstructures of Particle Reinforced Composites*, Acta Materialia 52, 4177-4188, 2004.

[3] L. Mishnaevsky Jr *Automatic Voxel-based Generation of 3D Microstructural FE Models and its Application to the Damage Analysis of Composites*, Materials Science and Engineering 407, 11-23, 2005.

[4] S. Boyd and R. Müller *Smooth Surface Meshing for Automated Finite Element Model Generation from 3D Image Data*, Journal of Biomechanics 39, 1287-1295, 2006.

[5] W. Leclerc, P. Karamian, A. Vivet and A. Campbell *Numerical Evaluation of the Effective Elastic Properties of 2D Overlapping Random Fibre Composites*, Technische Mechanik 32, 358-368, 2012.

[6] W. Leclerc, P. Karamian and A. Vivet *An Efficient Stochastic and Double-scale Model to Evaluate the Effective Elastic Properties of 2D Overlapping Random Fibre Composites*, Computational Materials Science 69, 481-493, 2013.

[7] W. J. Drugan and J. R. Willis, *A Micromechanics-based Nonlocal Constitutive Equation and Estimates of Representative Volume Element Size for Elastic Composites*, Journal of the Mechanics and Physics of Solids 44, 497-524, 1996.

[8] T. Kanit, S. Forest, I. Galliet, V. Mounoury and D. Jeulin *Determination of the Size of the Representative Volume Element for Random Composites: Statistical and Numerical Approach*, International Journal of Solids and Structures 40, 3647-3679, 2003.

[9] E. Sanchez-Palencia *Non-homogeneous Media and Vibration Theory*, Springer, Berlin, 1980.

[10] V.I. Agoshkov and V.I. Lebedev *The Poincaré-Steklov's Operators and the Domain Decomposition Methods in Variational Problems*, Computational Processes and Systems, Nauka, 173-227, 1985.

[11] C. Fahrat, F.X. Roux *A Method of Finite Element Tearing and Interconnecting and its Parallel Solution Algorithm*, International Journal of Numerical Methods in Engineering 32, 1205-1227, 1991.

[12] M.D. Rintoul and S. Torquato *Reconstruction of the Structure of Dispersions*, Journal of Colloid and Interface Science 186, 467-476, 1997.

[13] J. Segurado and J.L.Lorca *A Numerical Approximation to the Elastic Properties Sphere-reinforced Composites*, Journal of the Mechanics and Physics of Solids 50, 2107-2121, 2002.

[14] H.J. Böhm, A. Eckschlager and W. Han *Multi-inclusion Unit Cell Models for Metal Matrix Composites with Randomly Oriented Discontinuous Reinforcement*, Computational Materials Science 25, 42-53, 2002.

[15] B.R. Barello and M. Lévesque *Comparison between the Relaxation Spectra obtained from Homogenization Models and Finite Elements Simulation for the Same Composite*, International Journal of Solids and Structures 45, 850-867, 2008.

[16] E. Ghossein and M. Lévesque *A Fully Automated Numerical Tool for a Comprehensive Validation of Homogenization Models and its Application to Spherical Particles Reinforced Composites*, International Journal of Solids and Structures 49, 1387-1398 2012.