# Un élément fini $C^1$ pour l'analyse thermomécanique de poutres sandwich

Philippe Vidal, Olivier Polit

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

To cite this version:

Philippe Vidal, Olivier Polit. Un élément fini $C^1$ pour l'analyse thermomécanique de poutres sandwich. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01813027

HAL Id: hal-01813027

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

Submitted on 12 Jun 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_2.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Un élément fini C1 pour l'analyse thermomécanique de poutres sandwich

Philippe Vidal* — Olivier Polit*

*Laboratoire de Mécanique de Paris X (LMpX)
1, chemin Desvallières
F-92410 Ville d'Avray
{vidal,polit}@u-paris10.fr

RÉSUMÉ. Ce travail présente un nouvel élément fini poutre 1D prenant en compte le couplage thermomécanique pour des structures stratifiées. Le but de cette étude est de proposer un EF simple et efficace. Il est basé sur une approche de type couche équivalente avec l'exploitation d'un modèle raffiné pour la partie mécanique. Ceci permet de s'affranchir du calcul des coefficients correcteurs en cisaillement transverse. En ce qui concerne la partie thermique, une approche type "layerwise" est mise en oeuvre. Il faut noter que cet élément assure à la fois la continuité des déplacements et des contraintes de cisaillement aux interfaces entre les couches, ainsi que la continuité de la composante transverse du flux thermique. Les conditions limites sur les faces supérieures et inférieures sont également vérifiées. Les résultats obtenus sont comparés à des solutions exactes. Une étude de convergence (vitesse, précision) est également effectuée.

**ABSTRACT.** This paper deals with a new 1D beam finite element including thermomechanical coupling for laminated structures. A simple and efficient F.E. is presented. It is based on an equivalent single layer approach with a refined model for the mechanical part, avoiding the use of transverse shear correction factor. For the thermal part, a layerwise approach is carried out. It should be noted that both the continuity of displacements and shear stress, and the continuity of the transverse component of the heat flux, at the layer interfaces are satisfied. Boundary conditions on top and bottom faces are also satisfied. Results are compared to exact solutions. A convergence study is carried out.

**MOTS-CLÉS :** couplage thermomécanique, élément fini, poutre multi-couches, contrainte transverse d'ordre élevé, layerwise.

**KEYWORDS:** thermomechanical coupling, multilayered beam, finite element, higher order transverse shear, layerwise approach.
# 1. Introduction

Les composites stratifiés, de plus en plus utilisés en aéronautique, spatial... grâce à leur excellente propriété mécanique, sont soumis à des sollicitations thermiques sévères. Afin de prendre en compte les effets locaux aux interfaces qui conditionnent la rupture, il s'avère donc nécessaire de considérer le couplage thermomécanique dans les calculs de dimensionnement. Dans ce but, nous développons ici un élément fini poutre 1D prenant en compte ces couplages pour des poutres stratifiés. Pour les composites, les théories classiques des stratifiés, ou les théories de déformation du premier ordre, ou d'ordre supérieure ont été étendues aux problèmes thermiques en poutre et plaque [RED 97, TAN 89], mais elles ne permettent pas de prendre en compte la continuité du flux thermique dans l'épaisseur. Ainsi, des approches "couches discrètes" sont apparues pour pallier cet inconvénient : [ALI 99, KAP 03] qui nécessite un nombre élevé d'inconnues. Nous nous situons ici dans une approche de type layerwise pour la variation de températures dans l'épaisseur, introduite par [RED 87] en élasticité. De plus, la continuité du flux à l'interface des couches et la vérification des conditions aux limites sur les faces supérieure et inférieure sont prises en compte dans notre approche en utilisant la méthodologie introduite par [BLA]. Pour la partie mécanique, une théorie de type couche équivalente est utilisée avec l'exploitation d'un modèle raffiné. Un élément fini conforme efficace et sans pathologie numérique a été développé dans [GAN 99], et est réutilisé dans ce travail.

Un premier exemple de poutre sandwich soumis à un chargement thermique est présenté. Ensuite, un calcul couplé thermomécanique est mis en œuvre et les résultats obtenus sont comparés à une solution de référence. Des études de convergence sont effectuées afin de montrer l'efficacité de cet E.F à la fois sur les aspects vitesse et précision en déplacement/température et contrainte/flux.

# 2. Le problème de conduction

Considérons une poutre de section rectangulaire occupant le domaine $B = \Omega \times [-\frac{h}{2} \le z \le \frac{h}{2}] \times [-\frac{b}{2} \le x_2 \le \frac{b}{2}]$. Elle est constituée de N couches pouvant présenter chacune un comportement isotrope transverse (Cf. Figure 1). Le problème à résoudre consiste à chercher la température $\theta(M)$ et le flux $q(M)$ qui vérifient les conditions limites, l'équation de la chaleur, et la loi de Fourier.

## 2.1. Modèle "couche discète"

On écrit le champ de températures continu dans chaque couche sous la forme :

$$ \theta(x_1, x_3 = z) = \sum_{k=1}^{N} \theta^{(k)}(x_1, x_3) \chi_{[z_k, z_{k+1}[} \quad \text{où} \quad \chi_{[z_k, z_{k+1}[} = \begin{cases} 1 & \text{si } z \in [z_k, z_{k+1}[, \text{ 0 sinon.} \\ 0 & \text{sinon.} \end{cases} \tag{1} $$

On néglige la variation de $\theta$ dans la direction $(0, x_2)$. La variation de la température suivant $x_1$ est liée à la discrétisation 1D. Pour la direction $z$, la technique layerwise
est utilisée. On considère une variation quadratique de la température dans ces deux directions. En introduisant les coordonnées adimensionnelles suivant $x_1$ et $z$, notées $\xi$ et $\zeta$ respectivement, on a dans chaque couche notée (k) :

$$ \text{pour } \zeta\epsilon[-1, 1] \text{ et } \xi\epsilon[-1, 1], \zeta = 2\frac{[z - (z_k + z_{k+1})/2]}{z_{k+1} - z_k} $$

$$ \theta^{(k)}(\xi, \zeta) = \frac{1}{2}\zeta(\zeta - 1)\theta_b^{(k)}(\xi) + \frac{1}{2}\zeta(\zeta + 1)\theta_b^{(k+1)}(\xi) + (1 - \zeta^2)\theta_m^{(k)}(\xi) \quad [2] $$

où $\theta_b^{(k)}$, $\theta_m^{(k)}$ sont les températures en bas et au milieu respectivement de la couche (k), conformément à la Figure 1.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** empilement du stratifié et description des inconnues en température dans l'épaisseur</figcaption></figure>

Les inconnues du problème de conduction thermique sont donc les températures $\theta_b^{(k)}$ et $\theta_m^{(k)}$ pour $k = 1, ..., N$, et $\theta_b^{(N+1)}$ (Figure 1). Ainsi, pour un stratifié à N couches, et une approximation quadratique, le nombre d'inconnues s'élèvent à $2N+1$ températures, si l'on suppose une seule subdivision numérique par couche.

## 2.2. Continuité du flux thermique

L'approche développée précédemment assure naturellement la continuité des températures entre les couches du stratifié, mais ne respecte pas celle du flux de température dans l'épaisseur. Pour respecter la physique, nous imposons cette continuité aux interfaces entre les différentes couches, ce qui va permettre de réduire le nombre d'inconnues indépendantes du problème. Ainsi, les N-1 conditions de continuité à l'interface entre la couche (k) et (k + 1) s'écrivent :

$$ [-k_3^{(k)} \theta_{z}^{(k)}]_{z=z_{k+1}} = [-k_3^{(k+1)} \theta_{z}^{(k+1)}]_{z=z_{k+1}}, \quad k = 1, ..., N-1 \quad [3] $$

avec $k_3$ coefficient de conductivité transverse. On prend également en compte les conditions limites en surface haute et basse du composite (température et/ou flux), ce qui n'induit plus que N inconnues $\theta_m^{(1)}, \theta_m^{(2)}, ..., \theta_m^{(N)}$, soit une inconnue par couche.
<table><caption>Tableau 1. températures en $x_1 = L/2$ et flux thermique $q_3$ dans l'épaisseur</caption><thead><tr><th colspan="3">températures</th><th colspan="5">flux q<sub>3</sub></th></tr><tr><th>z</th><th>exact</th><th></th><th>1/1/1 c</th><th>1/1/1</th><th>1/2/1</th><th>2/3/2</th><th>4/6/4</th><th>exact</th></tr></thead><tbody><tr><td>0.00</td><td>0.000</td><td>0.000</td><td>-0.066</td><td>-0.064</td><td>-0.064</td><td>-0.064</td><td>-0.064</td><td>-0.065</td></tr><tr><td>0.10<sup>-</sup></td><td>0.013</td><td>0.013</td><td>-0.064</td><td>-0.065</td><td>-0.065</td><td>-0.065</td><td>-0.065</td><td>-0.066</td></tr><tr><td>0.10<sup>+</sup></td><td></td><td></td><td></td><td></td><td></td><td>-0.065</td><td>-0.065</td><td></td></tr><tr><td>0.20<sup>-</sup></td><td>0.026</td><td>0.026</td><td>-0.063</td><td>-0.067</td><td>-0.067</td><td>-0.067</td><td>-0.067</td><td>-0.068</td></tr><tr><td>0.20<sup>+</sup></td><td></td><td></td><td>-0.063</td><td>-0.000</td><td>-0.059</td><td>-0.064</td><td>-0.066</td><td>-0.068</td></tr><tr><td>0.30<sup>-</sup></td><td>0.027</td><td>0.027</td><td>-0.441</td><td>-0.377</td><td>-0.398</td><td>-0.397</td><td>-0.395</td><td>-0.396</td></tr><tr><td>0.30<sup>+</sup></td><td></td><td></td><td></td><td></td><td></td><td></td><td>-0.394</td><td></td></tr><tr><td>0.40<sup>-</sup></td><td>0.028</td><td>0.028</td><td>-0.820</td><td>-0.753</td><td>-0.738</td><td>-0.730</td><td>-0.733</td><td>-0.733</td></tr><tr><td>0.40<sup>+</sup></td><td></td><td></td><td></td><td></td><td></td><td>-0.724</td><td>-0.731</td><td></td></tr><tr><td>0.50<sup>-</sup></td><td>0.030</td><td>0.030</td><td>-1.198</td><td>-1.129</td><td>-1.077</td><td>-1.094</td><td>-1.087</td><td>-1.089</td></tr><tr><td>0.50<sup>+</sup></td><td></td><td></td><td></td><td></td><td>-1.060</td><td></td><td>-1.087</td><td></td></tr><tr><td>0.60<sup>-</sup></td><td>0.032</td><td>0.032</td><td>-1.576</td><td>-1.505</td><td>-1.481</td><td>-1.463</td><td>-1.469</td><td>-1.472</td></tr><tr><td>0.60<sup>+</sup></td><td></td><td></td><td></td><td></td><td></td><td>-1.457</td><td>-1.469</td><td></td></tr><tr><td>0.70<sup>-</sup></td><td>0.035</td><td>0.036</td><td>-1.955</td><td>-1.881</td><td>-1.903</td><td>-1.899</td><td>-1.887</td><td>-1.891</td></tr><tr><td>0.70<sup>+</sup></td><td></td><td></td><td></td><td></td><td></td><td></td><td>-1.886</td><td></td></tr><tr><td>0.80<sup>-</sup></td><td>0.040</td><td>0.040</td><td>-2.333</td><td>-2.257</td><td>-2.324</td><td>-2.340</td><td>-2.354</td><td>-2.356</td></tr><tr><td>0.80<sup>+</sup></td><td></td><td></td><td>-2.333</td><td>-2.338</td><td>-2.338</td><td>-2.352</td><td>-2.356</td><td>-2.356</td></tr><tr><td>0.90<sup>-</sup></td><td>0.513</td><td>0.514</td><td>-2.400</td><td>-2.401</td><td>-2.401</td><td>-2.386</td><td>-2.390</td><td>-2.390</td></tr><tr><td>0.90<sup>+</sup></td><td></td><td></td><td></td><td></td><td></td><td>-2.386</td><td>-2.390</td><td></td></tr><tr><td>1.00</td><td>1.000</td><td>1.000</td><td>-2.467</td><td>-2.465</td><td>-2.465</td><td>-2.479</td><td>-2.483</td><td>-2.483</td></tr><tr><td>ddls</td><td>45</td><td></td><td>45</td><td>75</td><td>105</td><td>195</td><td>405</td><td>-</td></tr></tbody></table>

## 2.3. Résultats numériques

On se place dans le cadre d'une méthode d'élément fini conforme, avec un E.F de poutre 1D à 3 noeuds. Une technique de condensation classique est employée pour imposer les conditions de continuité sur le flux, et les conditions limites sur les faces supérieures et inférieures. Une poutre sandwich composée d'un coeur et d'une peau isotrope ($k_{skin} = 0.5$W/(K.m), $k_{core} = 50.5$W/(K.m)) d'épaisseur respective $h - 2e$, et $e$ est soumise à un chargement thermique : sur le sommet $\bar{\theta}_t = \theta_c \sin(\pi x_1/L)$ avec $\theta_c = 1$, sur les trois autres cotés, $\bar{\theta} = 0$. De plus, on a $S = L/h = 2$, $e = 0.2h$, et $h = 1$.

Le tableau 1 montre que la température est en excellente corrélation avec la solution exacte. On note que le fait d'imposer des conditions de continuité sur le flux $q_3$ ("c") améliore sensiblement le résultat au niveau des interfaces des couches. L'errur sur la valeur maximale du flux est inférieure à 1 % dans tous les cas, même avec seulement une subdivision par couche. Les différentes subdivisions (1/2/1, 2/3/2, 4/6/4) témoignent d'une convergence prononcée vers la solution exacte dans le cas où l'on impose aucune condition de continuité. Néanmoins, elles nécessitent beaucoup de ddls.
### 3. Le couplage thermomécanique

Le couplage thermomécanique s'effectue par la relation de comportement. Le problème couplé à résoudre se décompose alors en deux étapes : tout d'abord, résoudre le problème thermique, qui constitue ensuite un chargement pour le problème thermomécanique. On ne décrira pas ici les approximations E.F. pour la mécanique, dont les détails sont donnés dans [GAN 99]. Un modèle sinus avec ("sin-c") et sans continuité ("sin") est exploité. On se concentre sur les résultats numériques concernant les contraintes qui sont déterminantes dans le cadre de la conception.

Le problème à résoudre (Cf. [ALI 99]) est le suivant : une poutre stratifiée (0°/90°/0°) est simplement posée, et soumise à un chargement thermique tel que

$$ \theta(x_1, z) = \theta_0 \frac{2z}{h} \sin(\frac{\pi x_1}{L}) $$

$$ \frac{E_T}{E_L} = 25, \frac{G_{LT}}{G_{LT}} = 0.5, \frac{G_{TT}}{G_{TT}} = 0.2, \nu_{LT} = \nu_{TT} = 0.25, \frac{\alpha_T}{\alpha_L} = 1125. $$

Le tableau 2 témoigne de résultats très satisfaisants pour la contrainte de cisaillement transverse avec peu d'éléments (8 éléments), alors que la valeur convergée de la contrainte axiale est obtenue avec 2 éléments, ce qui montre les performances et la convergence rapide de cet E.F.

<table><caption>**Tableau 2.** $\bar{\sigma}_{13}(0, \pm h/6)$, $\bar{\sigma}_{11}(L/2, \pm h/6)$ pour différents nombres de ddls (sin-c) : étude de convergence - S=50</caption><thead><tr><th>N</th><th>Nbre ddls</th><th>σ̄<sub>11</sub>(L/2, ±h/6)</th><th>erreur</th><th>σ̄<sub>13</sub>(0, ±h/6)</th><th>erreur</th></tr></thead><tbody><tr><td>2</td><td>12</td><td>± 374.0</td><td>&lt; 1%</td><td>0.491</td><td>17%</td></tr><tr><td>4</td><td>24</td><td>± 374.1</td><td>&lt; 1%</td><td>0.574</td><td>3.5%</td></tr><tr><td>8</td><td>48</td><td>± 374.1</td><td>&lt; 1%</td><td>0.596</td><td>&lt; 1%</td></tr><tr><td>10</td><td>60</td><td>± 374.1</td><td>&lt; 1%</td><td>0.598</td><td>&lt; 1%</td></tr><tr><td>12</td><td>72</td><td>± 374.1</td><td>&lt; 1%</td><td>0.599</td><td>&lt; 1%</td></tr><tr><td>exact</td><td></td><td>± 371.4</td><td></td><td>0.595</td><td></td></tr></tbody></table>

Le tableau 3 montre une très bonne corrélation pour la contrainte axiale, et pour la contrainte de cisaillement transverse, avec la solution exacte, pour des poutres semi-épaisse à très mince.

<table><caption>**Tableau 3.** $\bar{\sigma}_{13}(0, \pm h/6)$, $\bar{\sigma}_{11}(L/2, \pm h/6)$ pour différentes valeurs de S</caption><thead><tr><th rowspan="2">S</th><th colspan="4">σ̄<sub>11</sub>(L/2, ±h/6)</th><th colspan="4">σ̄<sub>13</sub>(0, ±h/6)</th></tr><tr><th>sin-c</th><th>sin</th><th>exact</th><th>erreur</th><th>sin-c</th><th>sin</th><th>exact</th><th>erreur</th></tr></thead><tbody><tr><td>20</td><td>± 374.0</td><td>374.1</td><td>371.5</td><td>&lt; 1%</td><td>1.509</td><td>1.500</td><td>1.441</td><td>4%</td></tr><tr><td>50</td><td>± 374.1</td><td>374.1</td><td>371.4</td><td>&lt; 1%</td><td>0.6008</td><td>0.6009</td><td>0.5948</td><td>&lt; 1%</td></tr><tr><td>100</td><td>± 374.1</td><td>374.1</td><td>371.4</td><td>&lt; 1%</td><td>0.3005</td><td>0.3004</td><td>0.2987</td><td>&lt; 1%</td></tr></tbody></table>

La répartition de la contrainte de cisaillement transverse dans l'épaisseur, qui est calculée par intégration des équations d'équilibre et adimensionnée (Cf. [ALI 99]), est en très bonne adéquation avec la solution exacte (Cf. figure 2). Il est à noter que les courbes "exactes" sont obtenues à partir du relevé de points sur les graphes de
<figure><img src="image_4.png" /><figcaption>**Figure 2.** $\bar{\sigma}_{13}(0, z)$ pour $S=20$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** $\bar{\sigma}_{11}(L/2, z)$ pour $S=20$</figcaption></figure>

[ALI 99]. Le résultat est également très satisfaisant pour la contrainte axiale (Cf. figure 3).

**4. Conclusions**

L'élément fini développé ici montre de bonnes performances aussi bien au niveau précision que vitesse de convergence, pour des problèmes thermiques et thermomécaniques. Il faut de plus noter qu'il tient compte des phénomènes physiques, puisqu'il assure la continuité du flux transverse, et des contraintes transverses aux interfaces des couches, ainsi que la vérification des conditions limites, à la fois en déplacement/température, et contrainte/flux thermique transverses. Cela permet en outre de réduire le nombre d'inconnues du problème, ce qui le rend finalement simple et efficace.

**5. Bibliographie**

[ALI 99] ALI J., BHASKAR K., VARADAN T., « A new theory for accurate thermal/mechanical flexural analysis of symmetric laminated plates », *Compos. Struct.*, vol. 45, 1999, p. 227-232.

[BLA] BLANC M., TOURATIER M., « A constrained discrete layer model for heat conduction in laminated composites », *soumis à Computers and Structures*.

[GAN 99] GANAPATHI M., PATEL B., POLIT O., TOURATIER M., « A $C^1$ finite element including transverse shear and torsion warping for rectangular sandwich beams », *Int. J. of Numerical Methods in Engng*, vol. 45, 1999, p. 47-75.

[KAP 03] KAPURIA S., DUMIR P., AHMED A., « An efficient higher order zigzag theory for composite and sandwich beams subjected to thermal loading », *Int. J. of Solids and Structures*, vol. 40, 2003, p. 6613-6631.

[RED 87] REDDY J., « A generalisation of two-dimensional theories of laminated composite plates », *Comm. in Applied and Numerical Methods*, vol. 3, 1987, p. 173-180.

[RED 97] REDDY J., *Mechanics of laminated composite plates - theory and analysis*, CRC Press, Boca Raton, FL, 1997.

[TAN 89] TANIGAWA Y., MURAKAMI H., OOTAO Y., « Transient thermal stress analysis of a laminated composite beam », *Int. J. Thermal Stresses*, vol. 12, 1989, p. 25-39.