# Technique(s) de raccord 2D-3D pour l'analyse non-intrusive de structures composites stratifiées

Guillaume Guguin, Olivier Allix, Pierre Gosselet

► To cite this version:

Guillaume Guguin, Olivier Allix, Pierre Gosselet. Technique(s) de raccord 2D-3D pour l'analyse non-intrusive de structures composites stratifiées. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717829

HAL Id: hal-01717829

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

Submitted on 26 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
# Technique(s) de raccord 2D-3D pour l'analyse non-intrusive de structures composites stratifiées

Guillaume Guguin<fnref n="1" />, Olivier Allix<fnref n="1" />, Pierre Gosselet<fnref n="1" />

**Résumé** — La simulation de structures industrielles à base de plaques composites requiert la prise en compte de phénomènes critiques arrivant à des échelles très petites devant la taille de la structure. Typiquement le délaminage se développe à l'échelle du micron pour des structures de l'ordre du mètre. Cela implique, pour obtenir une simulation fiable, de manipuler des modèles éléments finis contenant de très nombreux degrés de liberté. Notre objectif est de réduire la taille des modèles en proposant de coupler de manière non-intrusive un modèle plaque à un modèle 3D de façon à créer un modèle hybride, performant et pertinent, sans remettre en cause les principales pratiques des concepteurs.

## 1 Introduction

L'utilisation croissante des matériaux composites dans le milieu industriel pour fabriquer des composants critiques nécessite de nouvelles simulations prenant en compte de plus en plus de phénomènes. Un exemple est l'utilisation de liaisons boulonnées entre des plaques composites : les détériorations intervenant à des échelles variées en interaction (décohésion fibre-matrice, rupture fibre, délaminage), elles peuvent avoir des conséquences importantes sur l'intégrité de la structure complète. Leur prédiction implique de traiter des modèles éléments finis comportant de très nombreux degrés de liberté. Afin de réduire la complexité des calculs sans perdre en précision, on se propose de coupler un modèle plaque (global) classique, à un modèle 3D (local) fin, de façon à créer un modèle hybride, à la fois performant et pertinent. Ce couplage est réalisé de manière non-intrusive de façon à être facilement adaptable aux différents codes élément finis utilisés par les industriels ou en recherche.

En effet, les modélisations plaques sont connues pour donner des résultats corrects loin des zones critiques telles que les bords ou des défauts internes. Sur ces dernières zones, on substitue un modèle 3D (voir la figure 1). La difficulté réside dans la réalisation du couplage entre les deux modèles, le long de l'interface γ entre les deux zones.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Modèle 3D et modèle hybride</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Modèle couplé non-intrusif</figcaption></figure>

Cette étude se place dans le cadre d'une stratégie itérative non-intrusive [5] (développée sur des problèmes purement 3D) qui permet de converger vers le modèle souhaité à partir de calculs conduits sur un modèle global de l'ensemble de la structure $ω_i \cup ω_p$ et sur un modèle fin de la zone d'intérêt $\Omega_I$. Elle autorise la réutilisation des modèles et logiciels des industriels, mais permet également de les remplacer par des outils dédiés.

Cette approche peut être vue comme une correction des stratégies de "submodelling", très classiques

<footnote n="1">¹ LMT-Cachan — ENS-Cachan, CNRS, UPMC, Pres UniverSud Paris, 61 av. du Président Wilson, 94230 Cachan
guguin@lmt.ens-cachan.fr, allix@lmt.ens-cachan.fr, gosselet@lmt.ens-cachan.fr</footnote>
dans l'industrie, où suite à un calcul plaque global, un zoom 3D est calculé sur la zone d'intérêt avec des conditions aux limites déduites du calcul plaque. Ces techniques uniquement descendantes conduisent à des résultats erronés sur les quantités d'intérêts recherchées puisque les effets des phénomènes calculés dans la zone 3D sur le reste de la structure sont négligés. Dans notre approche, ces effets sont pris en compte par une stratégie de correction itérative sur le modèle plaque global (voir la figure 2). Dans cette stratégie, l'étape de descente est porteuse d'information en déplacement et celle de remontée est basée sur un calcul de résidu en effort entre les deux modèles.

Il existe une grande variété d'améliorations apportés aux théories des plaques classiques, souvent par la proposition de nouvelles formes de déplacement [7, 8] ou de contraintes [2] et en déduisant les modèles en découlant. Ces théories d'ordre plus élevé sont plus précises mais requièrent l'utilisation de logiciels ou procédures particuliers. De plus les modèles plaques sont souvent limités à des cas linéaires, sans présence de défauts ou d'endommagement, par conséquent une procédure permettant de reconstruire la solution 3D à partir de la solution plaque est nécessaire.

La reconstitution de quantités 3D à partir du calcul plaque n'est pas aisée. Les hypothèses cinématiques plaque (section droite et rigide) appliquées sur un modèle 3D imposent trop de contraintes à la structure et provoquent des effets de bords qui polluent la solution notamment sur les contraintes normales transverses. Ces dernières jouant un rôle important dans les cas de délaminage il est nécessaire de les représenter correctement. Des techniques améliorées ont été proposées pour obtenir des quantités cinématiques et statiques 3D de qualité mais leur mise en oeuvre n'est pas compatible avec notre objectif de couplage non-intrusif avec un code du commerce où les éléments finis disponibles sont limités et où le post-traitement de dérivées d'ordre supérieur est une source importante d'erreur.

Pour des cas isotropes, on peut appliquer une prédiction des contraintes/correction des déplacements par une méthode analytique, en comparant l'écart entre les contraintes attendues 3D et celles calculées [1]. Cependant cette procédure devient rapidement lourde dans le cas d'une plaque stratifiée (introduction de variables pour chaque couche). D'autre part une approche exacte de la solution plaque peut se révéler lourde et nécessiter des régularités particulières [6] peu compatibles avec le calcul élément fini. Certains auteurs ont aussi développé des stratégies basées sur l'ajout de gauchissements et une approche asymptotique [9, 10] qui via un calcul 1D dans l'épaisseur trouve les coefficients à appliquer au modèle de plaque. Cette procédure sert d'une part à déterminer le meilleur modèle Reissner-Mindlin mais fournit aussi les informations nécessaires à la reconstruction 3D du champs de contraintes et de déplacement.

Le raccord développé ici, contraint par l'objectif de non-intrusivité, repose sur une version étendue du raccord mixte [4] en remplaçant la définition analytique des contraintes dans l'épaisseur par des relèvements numériques qui permettent de s'adapter facilement à n'importe quelle séquence de plis.

Après une mise en place du problème, deux stratégies sont proposées dans la section 2. La première fait intervenir un raccord faible en déplacement entre les deux modèles de manière à transmettre les informations pertinentes du modèle plaque (effet à grande longueur caractéristique). Elle revient à rendre stationnaire un Lagrangien et souffre des risques classiques d'instabilité LBB en cas de discrétisations non-cohérentes entre le domaine plaque et le domaine 3D. La seconde stratégie contourne les problèmes de stabilité en s'appuyant sur un raccord en déplacement pur où les champs plaques sont corrigés par des gauchissements, eux aussi issus du relèvement numérique. Ces deux types de raccords ont été mis en oeuvre par un script pilotant le code_aster®, et un premier exemple est présenté dans la section 3.

## 2 Présentation du raccord plaque/3D

On considère un problème 3D sur le domaine $\Omega = \omega \times [-h, h]$. Avec $f$ les forces de volume, $g$ les efforts imposés sur $\partial_N \Omega$ et $u_D$ le déplacement imposé sur la partie complémentaire du bord $\partial_D \Omega$. Le comportement 3D est donné par le tenseur de Hooke $\mathbf{H}$. On note $\varepsilon(\mathbf{u})$ la partie symétrique du gradient de déplacement et $\sigma$ le tenseur de contrainte de Cauchy.

Soit $\mathcal{U}(\Omega)$ l'espace des déplacements cinématiquement admissible.

$$ \mathcal{U}(\Omega) = \{ \mathbf{u} \in H^1(\Omega), \mathbf{u} = \mathbf{u}_D \text{ sur } \partial_D \Omega \} $$

et on note $\mathcal{U}_0(\Omega)$ l'espace vectoriel associé.
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Illustration du problème 3D</figcaption></figure>

Le problème 3D s'écrit :

$$
\left\{
\begin{array}{l}
\text{Trouver } \underline{u} \in \mathcal{U}(\Omega) \text{ tel que } \forall \underline{u}^* \in \mathcal{U}_0(\Omega) \\
\int_{\Omega} \underline{\sigma} : \underline{\varepsilon}(\underline{u}^*) d\Omega = \int_{\Omega} \underline{f} \cdot \underline{u}^* d\Omega + \int_{\partial_N \Omega} \underline{g} \cdot \underline{u}^* dS = 0 \\
\underline{\sigma} = \mathbb{H} : \underline{\varepsilon}(\underline{u})
\end{array}
\right.
\qquad (1)
$$

On considère maintenant la modélisation plaque, on note avec un ~ les quantités planes.

On suppose une forme de déplacement de type Reissner-Mindlin, la section reste droite mais peut
former un angle avec la normale au plan moyen.

$$
\underline{v} = \tilde{\underline{v}} + \nu_z \underline{e}_z + z \tilde{\underline{\theta}} \wedge \underline{e}_z \quad (2)
$$

avec $\tilde{\underline{v}}$, $\nu_z$ et $\tilde{\underline{\theta}}$ des fonctions définies sur $\omega$.

La dérivation de cette forme de déplacement conduit aux déformations plaques :

$$
\begin{equation}
\begin{aligned}
\underline{\varepsilon}(\underline{v}^P) &= \begin{pmatrix} \underline{\tilde{\varepsilon}} \\ \text{sym} \\ \underline{\tilde{\varepsilon}}_z \end{pmatrix} \quad \text{avec} \quad \begin{cases} \underline{\tilde{\varepsilon}} = \tilde{\underline{\gamma}}(\tilde{\underline{v}}) + z \tilde{\underline{\chi}}(\tilde{\underline{\theta}}) \\ \underline{\tilde{\varepsilon}}_z = \frac{\tilde{\underline{\nabla}}(\underline{v}_z) + \tilde{\underline{\theta}} \wedge \underline{e}_z}{2} \end{cases} \\
\underline{\tilde{\gamma}}(\tilde{\underline{v}}) &= \begin{pmatrix} \underline{v}_{x,x} & \frac{\underline{v}_{x,y} + \underline{v}_{y,x}}{2} \\ \text{sym} & \underline{v}_{y,y} \end{pmatrix}, \qquad \underline{\tilde{\chi}}(\tilde{\underline{\theta}}) = \begin{pmatrix} \underline{\theta}_{y,x} & \frac{\underline{\theta}_{y,y} - \underline{\theta}_{x,x}}{2} \\ \text{sym} & -\underline{\theta}_{x,y} \end{pmatrix}
\end{aligned}
\tag{3}
\end{equation}
$$

En utilisant la notation $\langle \rangle$ pour l'intégration dans l'épaisseur, le tenseur des contraintes plaque devient :

$$
\underline{\tilde{N}} = \langle \underline{\tilde{\sigma}} \rangle \qquad \underline{\tilde{M}} = \langle z \underline{\tilde{\sigma}} \rangle \qquad \underline{\tilde{Q}} = \langle \underline{\tilde{\sigma}}_z \rangle \tag{4}
$$

Une autre hypothèse du modèle plaque est que le terme $\sigma_{zz}$ est négligeable devant les autres composantes du tenseur de contraintes. On peut alors écrire le comportement en terme de comportement plan $\tilde{\mathbb{H}}$ pour les composantes planes et avec $\tilde{\mathbb{B}}$ pour les composantes hors-plan (en 3D, $\tilde{\underline{\sigma}}_z = \tilde{\mathbb{B}} \underline{\varepsilon}_z$).

Les relations de comportement plaque s'écrivent alors :

$$
\begin{align}
\tilde{\underline{N}} &= \langle \tilde{\mathbb{H}} \rangle : \tilde{\underline{\gamma}}(\tilde{\underline{v}}) \nonumber \\
\tilde{\underline{M}} &= \langle z^2 \tilde{\mathbb{H}} \rangle : \tilde{\underline{\chi}}(\tilde{\underline{\theta}}) \nonumber \\
\tilde{\underline{Q}} &= \langle \tilde{\mathbb{B}} \rangle : \tilde{\underline{\varepsilon}}_z(\underline{v}_z, \tilde{\underline{\theta}}) \tag{5}
\end{align}
$$

L'espace des champs cinématiquement admissibles plaque est :

$$
\mathcal{U}^P(\omega) = \{ (\tilde{\underline{v}}, \tilde{\underline{\theta}}, \nu_z) \in H^1(\omega)^5, \text{satisfaisant les CL de Dirichlet} \} \quad (6)
$$

Et le problème plaque s'écrit<fnref n="1" /> :

$$
\left\{
\begin{array}{l}
\text{Trouver } (\tilde{\underline{v}}, \tilde{\underline{\theta}}, \nu_z) \in \mathcal{U}^P(\omega) \text{ tel que } \forall (\tilde{\underline{v}}^*, \tilde{\underline{\theta}}^*, \nu_z^*) \in \mathcal{U}_0^P(\omega) \\
\int_{\Omega} (\tilde{\underline{N}} : \tilde{\underline{\gamma}}(\tilde{\underline{v}}^*) + \tilde{\underline{M}} : \tilde{\underline{\chi}}(\tilde{\underline{\theta}}^*) + \tilde{\underline{Q}} : \tilde{\underline{\varepsilon}}_z(\nu_z^*, \tilde{\underline{\theta}}^*)) d\omega = \int_{\Omega} \underline{f} \cdot \underline{v}^* d\Omega + \int_{\partial_N \Omega} \underline{g} \cdot \underline{v}^* dS \quad (7) \\
\underline{\tilde{N}} = \langle \tilde{\mathbb{H}} \rangle : \tilde{\underline{\gamma}}(\tilde{\underline{v}}), \quad \underline{\tilde{M}} = \langle z^2 \tilde{\mathbb{H}} \rangle : \tilde{\underline{\chi}}(\tilde{\underline{\theta}}), \quad \underline{\tilde{Q}} = \langle \tilde{\mathbb{B}} \rangle : \tilde{\underline{\varepsilon}}_z(\underline{v}_z, \tilde{\underline{\theta}})
\end{array}
\right.
$$

<footnote n="1">1. Pour des raisons de simplicité, on ne développe pas le chargement en termes de quantité plaque.</footnote>
## 2.1 Proposition de couplage par Lagrangien

On définit le problème couplé entre le domaine 3D $\Omega_I$ et le domaine plaque $\omega_P$ le long de l'interface $\gamma \times [-h, h]$. Afin de simplifier la présentation, on traite ici le cas d'une interface droite de normale $\mathbf{e}_x$.

Trouver $\underline{u} \in \mathcal{U}(\Omega_I)$, $\underline{v} \equiv (\tilde{\nu}, \tilde{\Theta}, v_z) \in \mathcal{U}^P(\omega_P)$, tel que $\forall \underline{u}^* \in \mathcal{U}_0(\Omega_I)$, $\underline{v}^* \equiv (\tilde{\nu}^*, \tilde{\Theta}^*, v_z^*) \in \mathcal{U}_0^P(\omega_P)$,
$$ \int_{\Omega} \underline{\sigma} : \underline{\varepsilon}(\underline{u}^*) d\Omega = \int_{\Omega_I} \underline{f} \cdot \underline{u}^* d\Omega + \int_{\partial_I \Omega} \underline{g} \cdot \underline{u}^* dS + \sum_k \int_{\gamma_k} \lambda_k \int_z \underline{u}^* \cdot \underline{\tau}_k \cdot \mathbf{e}_x dz dl $$
$$ \int_{\omega} (\underline{\tilde{N}} : \underline{\tilde{\gamma}}(\tilde{\nu}^*) + \underline{\tilde{M}} : \underline{\tilde{\chi}}(\tilde{\Theta}^*) + \underline{\tilde{Q}} \cdot \underline{\tilde{\varepsilon}}_z^*(v_z^*, \tilde{\Theta}^*)) d\omega = \int_{\Omega} \underline{f} \cdot \underline{v}^* d\Omega + \int_{\partial_N \Omega} \underline{g} \cdot \underline{v}^* dS - \sum_k \int_{\gamma_k} \lambda_k \int_z \underline{v}^* \cdot \underline{\tau}_k \cdot \mathbf{e}_x dz dl $$
$$ \int_{\gamma} \lambda_k^* \int_z (\underline{u} - \underline{v}) \cdot \underline{\tau}_k \cdot \mathbf{e}_x dz dl = 0, \quad \forall k, \forall \lambda_k^* \in \Lambda(\gamma) $$
$$ \underline{\sigma} = \mathbb{H} : \underline{\varepsilon}(\underline{u}), \quad \underline{\tilde{N}} = \langle \tilde{\mathbb{H}} \rangle : \underline{\tilde{\gamma}}(\tilde{\nu}), \quad \underline{\tilde{M}} = \langle z^2 \tilde{\mathbb{H}} \rangle : \underline{\tilde{\chi}}(\tilde{\Theta}), \quad \underline{\tilde{Q}} = \langle \tilde{\mathbb{B}} \rangle : \underline{\tilde{\varepsilon}}_z(v_z, \tilde{\Theta}) \quad (8) $$

($\underline{\tau}_k$) est une famille de tenseurs de contraintes qui caractérise le couplage (elle dépend des matériaux et de la séquence d'empilement). Le raccord est imposé par des multiplicateurs de Lagrange ($\lambda_k$) qui assurent que le saut de déplacement ne génère pas de travail sur l'interface dans les champs ($\underline{\tau}_k \cdot \mathbf{e}_x$).

Si on analyse le terme de couplage dans la partie plaque on a :

$$ \int_{\gamma} \lambda_i \int_{z} \underline{v} \cdot \underline{\tau}_i \cdot \mathbf{e}_x dz dl = \int_{\gamma} \lambda_i (\underline{\tilde{N}}_i \cdot \mathbf{e}_x \cdot \tilde{\nu} + \underline{\tilde{M}}_i \cdot \mathbf{e}_x \cdot \tilde{\Theta} + \underline{\tilde{Q}}_i \cdot \mathbf{e}_x v_z) \quad (9) $$

où $(\underline{\tilde{N}}_i, \underline{\tilde{M}}_i, \underline{\tilde{Q}}_i)$ sont les contraintes plaque associées aux contraintes 3D $\underline{\tau}_i$. On note que pour la normale donnée $\mathbf{e}_x$, seulement 5 composantes sont actives, elles correspondent aux 5 degrés de liberté de plaque.

Si on analyse le terme de couplage dans la partie 3D des équations d'équilibre, on observe que cela implique la condition limite suivante :

$$ \int_{z} \underline{\sigma} \cdot \mathbf{e}_x \cdot \underline{u}^* dz = \sum_{i} \lambda_i \int_{z} \underline{\tau}_i \cdot \mathbf{e}_x \cdot \underline{u}^* dz \quad (10) $$

Donc $\underline{\sigma} \cdot \mathbf{e}_x$ doit être une combinaison linéaire des ($\underline{\tau}_i$). Si on choisit $\underline{u}^*$ le déplacement plaque, alors on observe que cette relation implique que chaque effort généralisé du problème 3D soit une relation linéaire des efforts généralisés associés à ($\underline{\tau}_i \cdot \mathbf{e}_x$) pondérés par les coefficients ($\lambda_i$).

Notons que l'utilisation d'un Lagrangien nécessite une étude précise pour le choix de la discrétisation des multiplicateurs ($\lambda_i$) le long de l'interface, afin d'éviter les instabilités LBB, en particulier pour des discrétisations incompatibles entre la plaque et le domaine 3D.

## 2.2 Construction de ($\underline{\tau}_i$)

L'étude d'un empilement de plis orthotropes ne permet pas de disposer, contrairement à [4] qui traite d'une plaque isotrope, d'une formule analytique simple des contraintes dans l'épaisseur ($\underline{\tau}_i$). Une telle formule étant lourde à obtenir et changeant pour chaque empilement, on propose un relèvement numérique de ($\underline{\tau}_i$). Pour cela on construit 8 cas avec des conditions limites choisies pour activer en particulier un certain effort généralisé (voir 1 et Figure 4-8.). Selon le principe de Saint-Venant, en prenant une plaque suffisamment grande, on évite la pollution de la solution intérieure par les effets de bords, et cette solution ne dépend que du matériau/empilement/maillage.

De chaque chargement $j = 1:8$, on extrait au milieu de la plaque les valeurs du tenseur de contrainte 3D ($\underline{\sigma}_j$) variant dans l'épaisseur. À partir de ce profil de contrainte, en utilisant (4), on calcule les efforts généralisés plaque ($\underline{\tilde{N}}_j, \underline{\tilde{M}}_j, \underline{\tilde{Q}}_j$) que l'on regroupe dans la notation suivante :

$$ (\mathcal{F}_j^i)_{i=1,8} = (N_{jxx}, N_{jxy}, M_{jxx}, M_{jxy}, Q_{jx}, N_{jyy}, M_{jyy}, Q_{jy}) \quad (11) $$

Grâce à la linéarité du problème, un état de contrainte issu d'une solution de St Venant est une combinaison des ($\underline{\sigma}_j$) que l'on a extrait. Les ($\underline{\tau}_i$) sont construit dans le but de découpler les contributions des efforts généralisés :

$$ \sum_i \mathcal{F}_j^i \underline{\tau}_i = \underline{\sigma}_j \Rightarrow [\underline{\tau}_i] = [\mathcal{F}_j^i]^{-1} [\underline{\sigma}_j] \quad (12) $$
<table><caption>Tableau 1 – Construction des solutions de St Venant</caption><thead><tr><th>j</th><th>chargement</th><th>sollicitation principale</th><th>figure</th></tr></thead><tbody><tr><td>1</td><td>σ<sub>xx</sub> = -1 sur ∂Ω<sub>0</sub>, σ<sub>xx</sub> = 1 sur ∂Ω<sub>L</sub></td><td>N<sub>xx</sub></td><td>4</td></tr><tr><td rowspan="2">2</td><td>σ<sub>xy</sub> = -1 sur ∂Ω<sub>L</sub>, σ<sub>xy</sub> = 1 sur ∂Ω<sub>0</sub></td><td rowspan="2">N<sub>xy</sub></td><td rowspan="2">5</td></tr><tr><td>σ<sub>xy</sub> = -1 sur ∂Ω<sub>a</sub>, σ<sub>xy</sub> = 1 sur ∂Ω<sub>-a</sub></td></tr><tr><td>3</td><td>σ<sub>xx</sub> = z sur ∂Ω<sub>0</sub>, σ<sub>xx</sub> = z sur ∂Ω<sub>0</sub></td><td>M<sub>xx</sub></td><td>6</td></tr><tr><td rowspan="2">4</td><td>σ<sub>xy</sub> = -z sur ∂Ω<sub>0</sub>, σ<sub>xy</sub> = z sur ∂Ω<sub>L</sub></td><td rowspan="2">M<sub>xy</sub></td><td rowspan="2">7</td></tr><tr><td>σ<sub>xy</sub> = z sur ∂Ω<sub>-a</sub>, σ<sub>xy</sub> = -z sur ∂Ω<sub>a</sub></td></tr><tr><td>5</td><td>σ<sub>xz</sub> = 1 sur ∂Ω<sub>L</sub></td><td>Q<sub>x</sub></td><td>8</td></tr><tr><td>6</td><td>même problème que 1 avec une rotation de 90°</td><td>N<sub>yy</sub></td><td></td></tr><tr><td>7</td><td>même problème que 3 avec une rotation de 90°</td><td>M<sub>yy</sub></td><td></td></tr><tr><td>8</td><td>même problème que 5 avec une rotation de 90°</td><td>Q<sub>y</sub></td><td></td></tr></tbody></table>

<figure><img src="image_5.png" /></figure>

où les crochets représentent les notations matrice/vecteur. Alors les $(\tau_{\equiv i})$ sont obtenus après inversion d'une matrice $8 \times 8$.

Dans les graphes 9-11 on a représenté un exemple de profil de contrainte que l'on peut obtenir pour le matériau décrit section 3.

<figure><img src="image_6.png" /></figure>

## 2.3 Approche par gauchissement

On se propose de remplacer le raccord par Lagrangien par un raccord où le déplacement 3D $u$ a une forme donné sur $\Gamma: u = v + w$ où $v$ est le déplacement plaque, et $w$ est un gauchissement qui dépend de l'état de contrainte de la plaque ($F|_{\gamma}$).
Le problème couplé s'écrit alors :

$$
\begin{align*}
& \text{Trouver } \underline{u} \in \mathcal{U}(\Omega_I), \underline{v} \equiv (\tilde{\nu}, \tilde{\Theta}, v_z) \in \mathcal{U}^P(\omega_P), \\
& \text{tel que } \forall \underline{u}^* \in \mathcal{U}_0(\Omega_I), \underline{v}^* \equiv (\tilde{\nu}^*, \tilde{\Theta}^*, v_z^*) \in \mathcal{U}_0^P(\omega_P), \\
& \int_{\Omega} \underline{\underline{\sigma}} : \underline{\underline{\varepsilon}}(\underline{u}^*) d\Omega = \int_{\Omega_I} \underline{f} \cdot \underline{u}^* d\Omega + \int_{\partial_I \Omega} \underline{g} \cdot \underline{u}^* dS \\
& \int_{\omega} (\underline{\underline{N}} : \underline{\underline{\gamma}}(\tilde{\nu}^*) + \underline{\underline{M}} : \underline{\underline{\chi}}(\tilde{\Theta}^*) + \underline{\underline{Q}} \cdot \underline{\underline{\varepsilon}}_z(v_z^*, \tilde{\Theta}^*)) d\omega = \int_{\Omega} \underline{f} \cdot \underline{v}^* d\Omega + \int_{\partial_N \Omega} \underline{g} \cdot \underline{v}^* dS \tag{13} \\
& \underline{u}|_{\Gamma} = \underline{v}|_{\gamma} + \sum_k w_k \mathcal{F}_{|\gamma}^k \\
& \langle \underline{\underline{\sigma}} \cdot \underline{\underline{n}} \cdot \underline{\underline{v}}^* \rangle = \underline{\underline{N}} : \underline{\underline{\gamma}}(\tilde{\nu}^*) + \underline{\underline{M}} : \underline{\underline{\chi}}(\tilde{\Theta}^*) + \underline{\underline{Q}} \cdot v_z^* \\
& \underline{\underline{\sigma}} = \mathbb{H} : \underline{\underline{\varepsilon}}(\underline{u}), \quad \underline{\underline{N}} = \langle \underline{\underline{H}} \rangle : \underline{\underline{\gamma}}(\tilde{\nu}), \quad \underline{\underline{M}} = \langle z^2 \underline{\underline{H}} \rangle : \underline{\underline{\chi}}(\tilde{\Theta}), \quad \underline{\underline{Q}} = \langle \underline{\underline{B}} \rangle : \underline{\underline{\varepsilon}}_z(v_z, \tilde{\Theta})
\end{align*}
$$

De manière analogue à l'approche de la section précédente, les fonctions de gauchissement sont déduites des problèmes de Saint-Venant (Figure 4-8). Des 8 expériences numériques, on obtient l'état de contrainte ($\sigma_j$) et le déplacement ($\underline{u}_j$) au centre de la plaque. En supposant que la partie plaque de ce déplacement développe les mêmes efforts généralisés que le déplacement 3D dans les chargements de St Venant, on peut extraire celle-ci :

$$
\begin{align}
v_{i_x} &= \langle \underline{\tau}_1 : (\underline{U}_i \otimes \underline{e}_x) \rangle \notag \\
v_{i_y} &= \langle \underline{\tau}_2 : (\underline{U}_i \otimes \underline{e}_x) \rangle = \langle \underline{\tau}_6 : (\underline{U}_i \otimes \underline{e}_y) \rangle \notag \\
\theta_{i_y} &= \langle \underline{\tau}_3 : (\underline{U}_i \otimes z\underline{e}_x) \rangle = \langle \underline{\tau}_7 : (\underline{U}_i \otimes z\underline{e}_y) \rangle \tag{14} \\
\theta_{i_x} &= -\langle \underline{\tau}_3 : (\underline{U}_i \otimes z\underline{e}_x) \rangle \notag \\
v_{i_z} &= \langle \underline{\tau}_5 : (\underline{U}_i \otimes \underline{e}_x) \rangle = \langle \underline{\tau}_8 : (\underline{U}_i \otimes \underline{e}_y) \rangle \notag \\
\underline{W}_i &= \underline{U}_i - (v_{i_x}\underline{e}_x + v_{i_y}\underline{e}_y + v_{i_z}\underline{e}_z + z\theta_{i_y}\underline{e}_x - z\theta_{i_x}\underline{e}_y) \notag
\end{align}
$$

On pose ainsi $\underline{W}_j = \underline{U}_j - \underline{V}_j^P$ pour chaque chargement de St Venant. Par linéarité entre les efforts et les gauchissements on peut réutiliser la matrice précédemment définie pour découpler les efforts pour les gauchissements

$$ [\underline{w}_i] = [\mathcal{F}_j^i]^{-1} [\underline{W}_j] \quad (15) $$

Les figures 12-14 présentent les trois composantes du gauchissement associé à un chargement de cisaille-ment.

<figure><img src="image_7.png" /><figcaption>Fig. 12 – Profil de $W_x^{Q_x}$ dans Fig. 13 – Profil de $W_y^{Q_x}$ dans Fig. 14 – Profil de $W_z^{Q_x}$ dans l'épaisseur l'épaisseur l'épaisseur</figcaption></figure>

## 3 Exemple d'application

Afin d'illustrer les deux approches, on prend l'exemple d'une plaque carré ($2L = 10, h = 0,4$) composée de 4 plis d'un matériaux orthotrope (voir ci-desous). L'orientation de chaque plis est donnée par la
configuration $(-45/45)_s$ :

$$ 
\begin{aligned}
E_L &= 25 \times 10^6, \quad E_T = E_N = 1 \times 10^6, \\
G_{LT} &= G_{LN} = 0.5 \times 10^6, \quad G_{TN} = 0.2 \times 10^6, \\
\nu_{LT} &= \nu_{TN} = \nu_{LN} = 0.25
\end{aligned}
$$

Un déplacement transverse $U_d$ est imposé à l'extrémité de la plaque. Le comportement plaque est défini par l'homogénéisation de Dhatt and Batoz [3]. Le maillage est effectué avec des éléments DSQ. Le modèle local 3D est maillé avec des éléments quadratique HEXA20, et supporte la localisation des contraintes due à l'encastrement près du bords. Sur l'interface, il y a correspondance noeud à noeud (autrement dit, une bande dans l'épaisseur d'éléments 3D est connectée à 2 éléments plaque). Dans le cas du raccord en Lagrangien, les multiplicateurs sont supposés constants par bande d'éléments 3D (ce qui permet d'éviter les instabilités LBB). La référence est le modèle 3D défini sur toute la structure.

<figure><img src="image_8.png" /><figcaption>Fig. 15 – $\sigma_{xx}$ dans l'épaisseur</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 17 – $\sigma_{xz}$ dans l'épaisseur</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 16 – $\sigma_{zz}$ dans l'épaisseur</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 18 – $\sigma_{yz}$ dans l'épaisseur</figcaption></figure>

Les contraintes des graphes 15-18 sont prises au milieu de l'interface dans l'épaisseur. Pour les deux méthodes, $\sigma_{xz}$ est de l'ordre de grandeur de $\sigma_{xx} h/L$ et $\sigma_{zz}$ de l'ordre de grandeur de $\sigma_{xx} h^2/L^2$. Néanmoins, on observe que le raccord par gauchissement donne des résultats plus proches de la référence. De plus ce dernier est plus simple à mettre en oeuvre (conditions de Dirichlet standard, au lieu de MPC sur chaque ligne de noeuds), et plus robuste puisque non sujette aux instabilités LBB.

On s'intéresse ensuite au résultat de la modélisation hybride en réalisant des itérations entre le calcul plaque et le calcul 3D en injectant le résidu en effort généralisé sur le modèle plaque [5]. On obtient alors un champ de déplacement corrigé, qui permet de relancer le calcul local. Le modèle couplé vers lequel on converge est illustré sur les figures 19 et 20.
<figure><img src="image_12.png" /><figcaption>Fig. 19 – Erreurs relatives au cours des itérations</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 20 – Evolution de $Q_x$ global et local au cours des itérations</figcaption></figure>

# 4 Conclusion

On a développé ici deux approches permettant d'étendre une technique de raccord non-intrusif à des modèles de dimensionnalité différente. Elles reposent sur le relèvement numérique d'un ensemble de fonctions créé *a priori*. Par sa simplicité, sa robustesse et son efficacité, l'approche en déplacement semble surclasser celle en multiplicateurs de Lagrange. Une modélisation hybride plaque/3D est obtenue à convergence de la méthode de couplage non-intrusive.

Les travaux en cours portent sur l'utilisation de modèles endommageables au sein de la zone 3D.

# Références

[1] O. Allix and C. Dupleix-Couderc, *Plate theory as a mean to compute precise 3D solutions including edge effects and related issues*, New Trends in Thin Structures : Formulation, Optimization and Coupled Problems, 1-28, 2010.

[2] Auricchio, F and Sacco, E, Refined first-order shear deformation theory models for composite laminates, 70,381,2003.

[3] Batoz, JL and Dhatt, GS, *Moderation des structures par éléments finis, Volume 2 : poutres et plaques*, Paris : Hermes, 1990.

[4] R. W. McCune, C. G. Armstrong, D. J. Robinson. *Mixed-dimensional coupling in finite element models*, International Journal for Numerical Methods in Engineering, 725-750, 2000.

[5] L. Gendre, O. Allix, P. Gosselet, F. Comte. *Non-intrusive and exact global/local techniques for structural problems with local plasticity*, Computational Mechanics, 233-245, 2009.

[6] P. Ladevèze, *La théorie exacte de la flexion des plaques*, Comptes Rendus de l'Académie des Sciences-Series IIB-Mechanics-Physics-Astronomy, 407-414, 1998.

[7] Reddy, JN, *A general non-linear third-order theory of plates with moderate thickness*, International Journal of Non-Linear Mechanics, 677-686, 1990.

[8] Savoia, M and Laudiero, F and Tralli, A, *A two-dimensional theory for the analysis of laminated plates*, Computational mechanics, 14,1,38–51,1994.

[9] Yu, W., *Mathematical construction of a Reissner–Mindlin plate theory for composite laminates*, International journal of solids and structures, 6680-6699, 2005.

[10] Yu, Wenbin and Hodges, Dewey H and Volovoi, Vitali V, *Asymptotically accurate 3-D recovery from Reissner-like composite plate finite elements*, Computers & structures, 81,7,439–454,2003.