# Méthode de surface de réponse pour l'optimisation de structures composites multicouches en utilisant les éléments coque-solide

Moez Hannachi, Hakim Naceur, Jean-Louis Batoz, Salim Belouettar

► To cite this version:

Moez Hannachi, Hakim Naceur, Jean-Louis Batoz, Salim Belouettar. Méthode de surface de réponse pour l'optimisation de structures composites multicouches en utilisant les éléments coque-solide. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01492448

HAL Id: hal-01492448

https://hal.science/hal-01492448

Submitted on 20 Mar 2017

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.

Public Domain
# Méthode de surface de réponse pour l'optimisation de structures composites multicouches en utilisant les éléments coque-solide.

M. Hannachi<sup>a</sup>, H. Naceur<sup>a</sup>, J.L. Batoz<sup>a</sup>, S. Belouettar<sup>b</sup>

RÉSUMÉ. Nous présentons dans la première partie de l'article la formulation d'un nouveau modèle E.F. tridimensionnel à huit nœuds de type Coque-Solide baptisé CSH8gZ7, valable pour l'analyse des structures isotropes et composites minces ou épaisses. La deuxième partie du papier traite de l'optimisation des structures stratifiées, où les variables de conception sont les orientations des fibres ou les propriétés du matériau. Pour l'optimisation, une méthode spécifique de construction de surfaces de réponse basée sur l'approximation diffuse est utilisée. Plusieurs applications numériques académiques et industrielles d'optimisation de structures composites multicouches sont conduites pour la validation des modèles.

ABSTRACT. In the first part of the paper we present the general formulation of a new eight-node solid-shell F.E. model called CSH8gZ7 using a 3D formulation, valid for thin or thick isotropic and composite shell structures. The second part of the paper deals with the optimization of multilayer composite structures, where design variables are fibers orientations and material parameters. To solve the optimization problem we adopt a specific response surface method based on diffuse approximation. Several numerical investigations are conducted in this study to assess current procedure capabilities applied to academic and industrial cases.

MOTS-CLÉS: éléments finis coque-solide, structures composites multicouches, optimisation.

KEYWORDS: Solid-Shell elements; Multilayer Composite structures; Optimization.

## 1. Introduction

Les éléments coque-solide sont des modèles qui ne possèdent que des degrés de liberté de translation pour simuler le comportement de structures de type coque ou solide. Ces éléments sont très avantageux, en raison de leur simplicité cinématique, leur capacité à modéliser des structures complexes composées de blocs épais et de régions à parois minces, en soulignant que le traitement laborieux des rotations lors de l'analyse non linéaire géométrique est évité. L'application de ces éléments à l'étude des structures minces cause des phénomènes de blocages numériques

<footnote><sup>a</sup> Université de Technologie de Compiègne, BP 20529 – 60205 Compiègne, France</footnote>

<footnote><sup>b</sup> Centre de Recherche public Henri Tudor, 29 Avenue J.F. Kennedy, L-1855, Luxembourg</footnote>
(BATOZ et al. 1992) qui engendrent des résultats peu probants. A cet effet, nous avons analysé les phénomènes de blocages (DOMISSY, 1997), (HAUPTMANN et al. 1998) ainsi que les moyens d'y remédier. Durant ces cinq dernières années, les méthodes de surface de réponse (RSM) ont gagné beaucoup d'importance dans l'optimisation des structures (NACEUR et al. 2006). La méthode RSM a l'avantage de remplacer un modèle de réponse complexe par un modèle approché évalué en un nombre limité de points disposés selon un plan d'expérience. Dans le présent travail, nous exploitons la RSM basée sur l'Approximation Diffuse (AD) (BELYTSCHKO et al. 1994). afin de concevoir un algorithme de minimisation original pour des modèles de régression (NACEUR et al. 2006).

## 2. Formulation EF de l'élément CSH8gZ7

### 2.1 Cinématique du modèle coque-solide

Nous rappelons brièvement dans cette section la formulation de l'élément coque-solide à huit nœuds. Les vecteurs position $X$ et déplacement $U$ d'un nœud de la surface moyenne sont défini par :

$$X(\xi, \eta, \zeta) = X_0(\xi, \eta) + \zeta X_n(\xi, \eta) = \sum_{i=1}^{4} N_i(\xi, \eta) \left[ \frac{1-\zeta}{2} X_i^{-} + \frac{1-\zeta}{2} X_i^{+} \right] \quad [1]$$

$$U(\xi, \eta, \zeta) = U_0(\xi, \eta) + \zeta U_n(\xi, \eta) = \sum_{i=1}^{4} N_i(\xi, \eta) \left[ \frac{1-\zeta}{2} U_i^{-} + \frac{1-\zeta}{2} U_i^{+} \right] \quad [2]$$

Les $N_i$ étant les fonctions d'interpolation bilinéaire du quadrangle représentant la surface moyenne. $X_i^{-}$, $U_i^{-}$ et $X_i^{+}$, $U_i^{+}$ sont respectivement les vecteurs position et déplacement du ième nœud des surfaces inférieure et supérieures de l'hexaèdre (Figure 1).

Les déformations covariantes infinitésimales sont données en fonction des coordonnées paramétriques comme suit :

$$\begin{align}
\varepsilon_{\xi\xi} &= X_{,\xi}^{T} U_{,\xi} ; \quad \varepsilon_{\eta\eta} = X_{,\eta}^{T} U_{,\eta} ; \quad \gamma_{\xi\eta} = X_{,\xi}^{T} U_{,\eta} + X_{,\eta}^{T} U_{,\xi} \nonumber \\
\varepsilon_{\xi\zeta} &= X_{,\xi}^{T} U_{,\zeta} ; \quad \gamma_{\xi\zeta} = X_{,\xi}^{T} U_{,\zeta} + X_{,\zeta}^{T} U_{,\xi} ; \quad \gamma_{\eta\zeta} = X_{,\eta}^{T} U_{,\zeta} + X_{,\zeta}^{T} U_{,\eta} \tag{3}
\end{align}$$

<figure><img src="image_1.png" /><figcaption>Figure 1 : modélisation d'un composite.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2 : Axes d'orthotropie.</figcaption></figure>
## 2.2 Loi de comportement

Chaque couche est constituée d'un matériau élastique orthotrope (directions {$m_1$; $m_2$; $m_3$}, Figure 2). Les coefficients non nuls de la matrice de comportement (entre contraintes et déformations physiques) sont donnés par les expressions suivantes :

$$H^{1111} = E_1(1 - \nu_{23}\nu_{32})/K, H^{1122} = E_1(\nu_{21} + \nu_{23})/K, H^{2222} = E_2(1 - \nu_{13}\nu_{31})/K,$$

$$H^{1113} = E_1(\nu_{31} + \nu_{21}\nu_{32})/K, H^{2223} = E_2(\nu_{32} + \nu_{12}\nu_{31})/K,$$

$H^{3333} = E_3(1 - \nu_{12}\nu_{21})/K, H^{1212} = G_{12}, H^{1313} = G_{13}, H^{2323} = G_{23}, \nu_{ij}E_j = \nu_{ji}E_i$ pour $(i,j=1,2,3 \text{ et } i \neq j)$. $E_1, E_2, E_3$ représentent les modules d'Young dans les directions d'orthotropie {$m_1$; $m_2$; $m_3$}, et $\nu_{ij}$ et $G_{ij}$ sont respectivement le coefficient de Poisson et le module de cisaillement.

## 2.3 Principe des travaux virtuels et matrice de rigidité

Une séparation des effets membrane/flexion, épaisseur et cisaillement transversal (Batoz, 1992), (Domissy, 1997) est nécessaire au traitement des différents blocages qui peuvent survenir. L'expression du travail virtuel interne est donnée par: $W_{\text{int}} = W_{\text{int}}^{\text{mf}} + W_{\text{int}}^{\text{ez}} + W_{\text{int}}^{\text{c}}$

$$W_{\text{int}}^{\text{mf}} = \int \delta \varepsilon_s^T \tilde{H}_1 \varepsilon_s \, dv = \delta u_n^T K_{mf} u_n ; K_{mf} = \int B_{mf}^T \tilde{H}_1 B_{mf} \, dv$$

$$W_{\text{int}}^{\text{ez}} = \int \delta \varepsilon_z^T \tilde{H}_2 \varepsilon_z \, dv = \delta u_n^T K_{ez} u_n ; K_{ez}^v = \int B_{ez}^T \tilde{H}_2 B_{ez} \, dv \quad [4]$$

$$W_{\text{int}}^{\text{c}} = \int \delta \gamma_s^T \tilde{H}_3 \gamma_s \, dv = \delta u_n^T K_c u_n ; K_c = \int B_c^T \tilde{H}_3 B_c \, dv$$

$$\text{avec : } \langle \varepsilon_s \rangle^V = \langle \varepsilon_{xx} \ \varepsilon_{yy} \ \gamma_{xy} \rangle, \langle \varepsilon_z \rangle = \langle \varepsilon_{xx} \ \varepsilon_{yy}^v \ \varepsilon_{z'z'}^v \rangle, \langle \gamma_s \rangle = \langle \gamma_{xz'} \ \gamma_{yz'} \rangle$$

## 2.4 Méthodes d'élimination de blocages

Pour réduire le blocage en cisaillement nous utilisons la méthode classique ANS (Batoz, 1992). Les déformations de cisaillement transverse $\gamma_{\xi\xi}^{ANS}$ et $\gamma_{\eta\zeta}^{ANS}$ sont calculées suivant le plan de la surface moyenne $\zeta = 0$. Elles sont exprimées en fonction des quantités $\gamma_{\xi}$ et $\gamma_{\eta}$ définies au milieu des côtés de la surface moyenne.

Le blocage trapézoïdal se produit quand les éléments hexaédriques à huit nœuds sont utilisés pour modéliser des coques à simple ou double courbure. Il peut être réduit en faisant une interpolation bilinéaire des contraintes normales transversales définies aux quatre coins de la surface moyenne de l'élément.

Afin d'obtenir une distribution linéaire de la contrainte normale dans la direction épaisseur, nécessaire à éviter le phénomène de blocage en épaisseur, nous enrichissons les champs de déformations normales dans la base locale, en ajoutant un degré de liberté interne $\alpha$. Ce 7ème paramètre sera éliminé par la technique de condensation statique.
### 3. Optimisation:

Le problème d'optimisation consiste à minimiser une fonction $f(x)$, soumise à un ensemble $m$ de contraintes $g_i(x) \le 0, j = 1,...,m$. Avec $f$ appelée fonction objectif, $x_i$ variables de conception et $g_j$ les limitations.

La méthode de surface de réponse que nous utilisons a été présentée dans de précédents travaux des auteurs (NACEUR et al. 2006), elle consiste à remplacer la fonction coût par une approximation $\tilde{f}$ à partir d'un nombre limité d'évaluations de la fonction réelle. Le pseudo problème d'optimisation devient alors :

$$ \text{Minimiser } \tilde{f} \qquad [5] \\ \text{avec } \tilde{g}_i(x) \le 0, \quad j = 1,...,m $$

La fonction $\tilde{f}$ peut être définie en fonction de la base polynomiale P et des coefficients $a$ non linéaires en $x$ telle que : $\tilde{f}(x) = p^T(x) \cdot a(x)$. Ces coefficients $a$ sont déterminés par minimisation d'une erreur $J(a)$ entre la fonction d'approximation et la fonction EF.

$$ J(a) = \sum_{i=1}^{N} w(\|x_i - x\|) (p^T(x_i - x) a - f(x_i))^2 \quad [6] $$

Moyennant des coefficients de pondérations diffus $w(x)$, la minimisation de $J(a)$ permet d'obtenir les coefficients $a$ et par conséquent la fonction $\tilde{f}$. Pour plus de détails nous renvoyons le lecteur à la référence NACEUR et al. 2006.

### 4. Optimisation d'une lame composite

#### 4.1. Analyse et position du problème

La structure est une coque épaisse composite (Figure 3), formée de deux couches d'épaisseurs identiques, la couche inférieure est en alliage d'Aluminium 6066-T651 (E=67GPa, $v=0.34$) et d'une couche externe orthotrope en fibre d'aramide haut module (E$_{11}$=130Gpa, E$_{22}$=E$_{33}$=5Gpa, $v_{12}=v_{13}=0.25$, $v_{23}=0$, G$_{23}$=5Gpa et G$_{12}$=G$_{13}$=10Gpa). La structure est soumise à une charge en A de 0.425kgf/mm et encastrée à son bord BC. La structure est maillée avec 1596 éléments coque-solide, et 3588 nœuds (Figure 3). Dans un premier temps la validation du modèle de maillage est faite en utilisant les paramètres matériaux initiaux et en comparant nos résultats à ceux en utilisant ABAQUS. La figure 7 montre la solution correspondante à la flèche sous charge en comparaison à ABAQUS qui donne 13.676 mm, nous remarquons que notre résultat 13.747 mm est en très bon accord avec celui obtenu par ABAQUS (Figure 5). Cela permet la validation du modèle E.F.
<figure><img src="image_3.png" /><figcaption>Figure 3 Géométrie et chargement.</figcaption></figure>

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

REFLEX

<figure><img src="image_5.png" /><figcaption>Figure 4 : Tête chercheuse composite centrale.</figcaption></figure>

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

ABAQUS

<caption>Figure 5 : Comparaison des déplacements (REFLEX, ABAQUS)</caption>

## 4.2. Optimisation des paramètres

Le problème d'optimisation consiste à trouver l'angle d'orientation optimal des fibres de la couche externe orthotrope en fibre d'aramide ainsi que le module d'Young de la couche en alliage d'aluminium, sachant que sous la charge le déplacement de la structure est imposé à 15mm. La fonction objective choisie est basée sur le critère général de Hill, le problème d'optimisation s'écrit ainsi :

$$ 
\left\{ 
\begin{array}{l} 
\text{Min } J(\theta, E) = \sum_{e=1}^{nelt} \left( F(\sigma_{22} - \sigma_{33})^2 + G(\sigma_{33} - \sigma_{11})^2 + H(\sigma_{11} - \sigma_{22})^2 + \right. \\ 
\left. \qquad + 2L \tau_{23}^2 + 2M \tau_{31}^2 + 2N \tau_{12}^2 \right) \\ 
\text{avec les limitations } \rightarrow \begin{cases} U_A(\theta, E) \geq 15 \text{mm} \\ -22.5^\circ \leq \theta \leq +22.5^\circ \\ 60 \text{GPa} \leq E \leq 79 \text{GPa} \end{cases} 
\end{array} 
\right. 
$$

Le modèle de surface de réponse global basé sur l'AD est construit en utilisant le plan d'expérience basé sur la tête chercheuse centrale composite à 9 paramètres (Figure 4). Cette fonction de réponse est estimée par :

$$ \tilde{J}(\theta, E) = 939800 + 58864.99 \theta + 29763.57 E + 11549.15 \theta E + 146500 \theta^2 + 1431.34 E^2 $$

Basée sur les données du même plan d'expérience, la limitation en déplacement est également approchée par notre modèle basé sur l'AD, elle est donnée par :

$$ \tilde{U}(\theta, E) = 13.08 + 3 \theta - 1.22 E - 0.17 \theta E + 4.58 \theta^2 + 0.14 E^2 $$
<figure><img src="image_7.png" /><figcaption>Figure 6 : Surface de réponse de la fonction objective.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7 : Représentation des iso contours de la solution optimale.</figcaption></figure>

La solution optimale obtenue après dénormalisation des variables correspond à $\langle \theta, E \rangle = \langle 3.4^\circ, 60.12\text{GPa} \rangle$. Le déplacement obtenu pour la structure optimisée est de $U=15.114\text{mm}$ ce qui est très satisfaisant par rapport à la limitation $U \ge 15$.

## 5. Conclusion

Dans cet article, nous avons proposé un nouveau modèle de coque-solide baptisé CSH8gZ7 très performant, valable pour les matériaux isotropes et composites multicouches. Ce modèle E.F. implémenté dans l'environnement REFLEX a été couplé à une méthode spécifique de surface de réponse basée sur l'AD pour l'optimisation de structures composites. La méthodologie proposée a montré une efficacité en terme de précision et capacité d'optimisation de structures composites multicouches. Des travaux sont en cours pour l'extension du modèle au cas de flambement linéaire et pour les problèmes non linéaires de façon générale.

## 6. Bibliographie

BATOZ J.L., DHATT G., "Modélisation des structures par éléments finis", vol. 2, 3, 1992, Paris, Editions Hermès.

BELYTSCHKO T., LU Y.Y., GU L., "Element-Free Galerkin Methods", Int. j. numer. thods eng., vol. 37, 229-256, 1994.

DOMISSY E., "Formulation et évaluation d'éléments finis volumiques modifiés pour l'analyse linéaire et non linéaire des coques", Thèse de doctorat, Université de Technologie de Compiègne, 1997.

HAUPTMANN R., SCHWEIZERHOF K., "A systematic development of solid-shell element formulation for linear and non-linear analysis employing only displacement degrees of freedom", Int. J. Numer. Meth. Engrg., vol.42, 1998, p.49-69.

NACEUR H., GUO Y.Q., BEN-ELECHI S., "Response surface methodology for design of sheet forming parameters to control springback effects", Computers & Structures 84, 1651–1663, 2006.