Travail sur le développement de nouveaux éléments finis solide-coques linéaires

ENSAM Rapport de Bachelor Arbeit Travail sur le d´eveloppement de nouveaux ´el´ements finis solide-coques lin´eaires ´ El`eve : Cl´ement RACINET Tuteur : Dr. Hocine CHALAL 28 mai 2021 Attestation d’authenticit´ e Je, soussign´e Cl´ementa RACINET, ´el`eve-ing´enieur du cursus Arts et M´etiers / KIT, atteste que ce travail est personnel, cite en r´ef´erence toutes les sources utilis´ees et ne comporte pas de plagiat. Ich versichere wahrheitsgem¨aß, die Arbeit selbstst¨andig verfasst, alle benutzten Hilfsmittel vollst¨andig und genau angegeben und alles kenntlich gemacht zu haben, was aus Arbeiten anderer unver¨andert oder mit Ab¨anderungen entnommen wurde sowie die Satzung des KIT und der Arts et M´etiers zur Sicherung guter wissenschaftlicher Praxis in der jeweils g¨ultigen Fassung beachtet zu haben. Le 28 mai 2021 1 NOTICE BIBLIOGRAPHIQUE ANNEE : 2020 TYPE DE DOCUMENT : Rapport de Bachelorarbeit CAMPUS DE RATTACHEMENT : Arts et M´etiers Metz AUTEUR : RACINET Cl´ement TITRE : Travail sur le d´eveloppement de nouveaux ´el´ements finis solides-coques lin´eaires. ENCADREMENT : CHALAL Hocine ENTREPRISE PARTENAIRE : Lem3 NOMBRES DE PAGES : 30 NOMBRE DE REFERENCES BIBLIOGRAPHIQUES : 16 MOTS CLES : ´ El´ements fini ; Lin´eaire ; Solide-Coque ; ´ Elasticit´e Lin´eaire ; Grandes d´eformations ; FGM R´esum´e Mon ´etude `a pour but d’aider `a d´evelopper des nouveaux ´el´ements de type Solide-Coques lin´eaires. Ces derniers allient les avantages de ´el´ements solides et les avantages des coques. Pour des simulations de pi`eces de faible ´epaisseur, ce type d’´el´ements peut permettre d’obtenir des r´esultats optimums `a moindre coˆut. De plus, ces ´el´ements ´evitent les probl`emes de verrouillages qui peuvent apparaˆıtre avec les ´el´ements solides ou les ´el´ements coques. Nous pr´esenterons les r´esultats de nos simulations avec ces nouveaux ´el´ements SHB8 et proposerons certaines pistes pour les am´eliorer dans le cas de simulation de mat´eriaux FGM. 2 Table des mati`eres 1 Introduction 4 2 Le principe de la m´ethode des ´el´ements finis 4 2.1 Discr´etisation de probl`emes complexes . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Discr´etisation des calculs dans le cas de la m´ecanique ´elastique . . . . . . . . . . 5 2.2.1 Rappel des ´equations de la m´ecanique ´elastique . . . . . . . . . . . . . . 6 2.2.2 Num´erisation de l’int´egrale et quadrature de Gauss . . . . . . . . . . . . 6 2.2.3 D´efinition et choix d’un ´el´ement fini . . . . . . . . . . . . . . . . . . . . . 6 2.3 Erreurs et limites du calcul par ´el´ements finis : le ph´enom`ene de verrouillage . . 8 2.3.1 Exemple de verrouillage . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4 Solutions pour des ph´enom`enes de verrouillages . . . . . . . . . . . . . . . . . . 8 3 Pr´esentation des outils et mod`eles utilis´es 10 3.1 Lestypesd’´el´ements ................................. 10 3.1.1 Les´el´ements2D................................ 10 3.1.2 Les ´el´ements 3D lin´eaires C3D8 . . . . . . . . . . . . . . . . . . . . . . . 10 3.1.3 Les ´el´ements 3D lin´eaires C3D8R . . . . . . . . . . . . . . . . . . . . . . 11 3.1.4 Les ´el´ements solide-coque . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Lesmod`elesutilis´es.................................. 13 3.2.1 PoutreenFlexion............................... 13 3.2.2 Plaque annulaire fendue . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2.3 La coque h´emisph´erique soumise `a des forces de traction et de compression radiales .................................... 15 3.2.4 La coque cylindrique `a extr´emit´e ouverte soumise `a des forces de traction radiales .................................... 16 4 Mise en ´evidence du probl`eme de verrouillage r´esolu pour les ´el´ements solidecoques 17 4.1 Le cas de la coque h´emisph´erique soumise `a des forces de traction et de compressionradiales...................................... 17 4.2 Le cas de la plaque annulaire fendue . . . . . . . . . . . . . . . . . . . . . . . . 18 4.3 Le cas de la poutre en flexion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.4 Analysedesr´esultats................................. 20 5 Simulation des mat´eriaux FGM 20 5.1 D´efinition du module de Young pour la stabilisation . . . . . . . . . . . . . . . . 21 5.2 Simulations de FGM et comparaisons avec les mod`eles de r´ef´erences . . . . . . . 22 5.2.1 Le cas de la coque cylindrique `a extr´emit´e ouverte soumise `a des forces detractionradiales.............................. 22 5.2.2 Le cas de la poutre en flexion . . . . . . . . . . . . . . . . . . . . . . . . 24 5.2.3 Analyse des r´esultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.3 R´eduction de l’erreur identifi´ee . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.3.1 R´esultats obtenu avec les nouvelles fonctions . . . . . . . . . . . . . . . . 26 6 Conclusion 28 7 Bibliographie 29 3 1 Introduction Il existe plusieurs mani`ere de simuler les tˆoles fines par ´el´ements finis : nous pouvons utiliser des ´el´ements 3D (solides) ou des ´el´ements 2D (de type coque). Cependant ces deux types d’´el´ements poss`edent certains avantages et d´efauts que nous d´etailleront par la suite. L’id´ee est d’allier les avantages des deux types d’´el´ements afin de cr´eer un nouveau type d’´el´ement qui est l’´el´ement Solide-Coque. Nous testerons cet ´el´ements en lin´eaire sur plusieurs mod`eles et nous simulerons notamment des mat´eriaux `a gradient de fonction (FGM). 2 Le principe de la m´ethode des ´el´ements finis La m´ethode des ´el´ements finis (MEF) ou FEM (Finite Element Method) est une m´ethode utilis´ee pour r´esoudre num´eriquement les probl`emes aux ´equations aux d´eriv´ees partielles. La m´ethode telle qu’on la connaˆıt aujourd’hui est relativement jeune puisqu’elle date de 1956 lorsque Turner et Clough entreprirent de d´ecomposer un domaine surfacique quelconque en petits domaines triangulaires `a d´eplacement lin´eaire. [16] [10] Le but de cette m´ethode qui fait partie de la famille des math´ematiques appliqu´ees est la r´esolution approch´ee d’un probl`eme complexe, o`u les solutions du probl`eme v´erifient des conditions d’existences plus faibles que celles des solutions du probl`eme de d´epart et o`u une discr´etisation permet de trouver une solution approch´ee. Comme de nombreuses autres m´ethodes num´eriques, outre l’algorithme de r´esolution en soi, se posent les questions de qualit´e de la discr´etisation : — Existence des solutions — Unicit´e de la solution — Stabilit´e — Convergence — Mesure et minimisation de l’erreur 2.1 Discr´etisation de probl`emes complexes Le concept fondamental de la m´ethode des ´el´ements finis r´eside en sa capacit´e de discr´etisation d’une topologie plus grande et plus complexe. Cette discr´etisation est r´ealis´ee par une ´etape de maillage. Il est d’usage d’utiliser un maillage carr´e ou triangulaire afin discr´etiser la topologie ´etudi´ee. Cependant, rien n’empˆeche l’utilisation d’´el´ements plus exotiques, bien qu’ils soient en g´en´eral plus complexes. De plus, il n’est pas n´ecessaire d’utiliser un maillage r´egulier. Si voulu, il est possible de cr´eer un maillage plus fin `a un endroit pr´ecis afin d’augmenter la pr´ecision et la fiabilit´e de la solution. En effet, en r`egle g´en´erale, plus un maillage est fin, plus l’erreur est faible et donc plus le mod`ele est proche de la r´ealit´e (voire figure 1 ). Comme le montre la figure 2, on appelle traditionnellement h la plus grande dimension d’un ´el´ement et p le degr´e du polynˆome d´ecrivant le cˆot´e ou l’arˆete (p = 1 pour des 4 Figure 1 – Convergence au maillage d’une poutre en flexion. cˆot´es/arˆetes droits, p = 2 pour des cˆot´es/arˆetes pr´esentant une courbure). Il existe donc deux grands types de raffinements de maillages : le raffinement de type h qui consiste `a diminuer la taille globale des ´el´ements ce qui am`ene `a cr´eer un plus grand nombre d’´el´ements et un raffinement de type p qui consiste `a augmenter le nombre de points d’extrapolation des cˆot´es des ´el´ements afin d’obtenir une meilleure extrapolation. Ces deux convergences am´eliorent la pr´ecision et la fiabilit´e des r´esultats mais complexifient le calcul, ce qui se traduit par un temps de calcul plus long. Figure 2 – Raffinement de maillage par diminution de la taille d’´el´ements (type h, `a gauche) ou par augmentation du degr´e des ´el´ements (type p, `a droite).[12] 2.2 Discr´etisation des calculs dans le cas de la m´ecanique ´elastique Apr`es discr´etisation du syst`eme sous forme d’´el´ements, une seconde discr´etisation peut ˆetre n´ecessaire afin de r´esoudre les ´equations de mani`ere num´erique (une r´esolution 5 analytique ´etant trop complexe ou impossible). 2.2.1 Rappel des ´equations de la m´ecanique ´elastique En notant [K] le tenseur de rigidit´e du syst`eme discr´etis´e, [B] le tenseur des d´eriv´ees des fonctions de forme d’interpolation, [C] le tenseur d’´elasticit´e du mat´eriau, [σ]le tenseur des contraintes et [ε]le tenseur d´eformation, on peut appliquer la loi de comportement du mat´eriau (loi de Hook) : [σ] = [C][ε] (1) On ´ecrit que le travail des forces externes appliqu´ees `a la structure pour un d´eplacement virtuel δU est ´egal au travail interne des contraintes pour ce mˆeme d´eplacement : [δU]T[F]=[δU]TZV [B]T[C][B]dV[U] (2) Comme cette relation est vraie pour tout d´eplacement virtuel, on en d´eduit : [U] = [K]−1[F] (3) Avec [K] sous sa forme plus g´en´erale : [K] = ZV [B]T[C][B]dV(4) 2.2.2 Num´erisation de l’int´egrale et quadrature de Gauss Afin de calculer l’int´egrale ais´ement avec un ordinateur, il est n´ecessaire de discr´etiser l’int´egrale. Pour ce faire, nous utilisons la quadrature de Gauss. Comme le montre la Figure 3, nous d´efinissons une fonction de pond´eration $continue strictement positive, qui peut assurer l’int´egrabilit´e de la fonction `a int´egrer. Les ωisont appel´es les coefficients de quadrature (ou poids). Les points xi sont appel´es les nœuds de la quadrature. Nous pouvons donc `a pr´esent calculer le tenseur de rigidit´e [K] `a condition de d´efinir les points d’int´egrations. Ces points d´ependent donc des ´el´ements choisis. 2.2.3 D´efinition et choix d’un ´el´ement fini On appelle ”´el´ement fini”, un triplet (ph´enom`ene, mod´elisation, type de maille). Le ph´enom`ene est le domaine physique qui r´egit l’´el´ement : m´ecanique, acoustique, thermique. On peut associer `a ce ph´enom`ene une mod´elisation. Par exemple un ´el´ement m´ecanique peut ˆetre mod´elis´e comme ´elastique dans le cas des petites perturbations, ´elastique avec des grandes d´eformations, etc ... A chaque mod´elisation, on y associe un type de maille, c’est-`a-dire une g´eom´etrie (ligne, triangle, rectangle, t´etra`edre, prisme, etc ...) `a laquelle on associe des points d’int´egrations (afin de pouvoir effectuer la quadrature de Gauss). Ainsi, il existe de nombreux types d’´el´ements (on en d´enombrait d´ej`a plus de 500 en 2003) comme le montrent les Figures 4 , 5. [8] 6 Figure 3 – Illustration de la quadrature de Gauss Figure 4 – Familles d’´el´ements finis couramment utilis´ees en m´ecanique [1] Figure 5 – Types d’´el´ements avec leurs points d’int´egrations [1] 7 2.3 Erreurs et limites du calcul par ´el´ements finis : le ph´enom`ene de verrouillage On appelle verrouillage d’un ´el´ement une insuffisance de l’approximation polynomiale qui empˆeche l’´el´ement de repr´esenter correctement certains modes de d´eformation. L’´el´ement se verrouille ou se bloque lorsque ces modes sont sollicit´es et les contraintes peuvent devenir tr`es fausses. 2.3.1 Exemple de verrouillage Nous nous pla¸cons dans une situation d’´el´ements massifs bidimensionnels en contraintes planes, o`u ´el´ement r´eel et ´el´ement de r´ef´erence sont confondus. Supposons que le champ vrai `a repr´esenter soit celui d’une flexion pure d´efini par : u=xy v =−1 2x2(5) Les champs de d´eformation et de contrainte associ´es sont : εxx =∂u ∂x =y εyy =∂v ∂y = 0 γxy =∂u ∂y +∂v ∂x = 0 (6) σxx =Ey 1−ν2σyy =νEy 1−ν2σxy = 0 (7) Si nous consid´erons un ´el´ement `a quatre noeuds, la repr´esentation ´el´ements finis de ce champ est : uEF =uh=xy vEF =vh=−1 2(8) Auxquels correspondent : εhxx =y εhyy = 0 γxy =x(9) σhxx =Ey 1−ν2σhyy =νEy 1−ν2σxy =Ex 2(1 + ν)(10) On remarque que le cisaillement est faux. Il en va de mˆeme pour l’´energie de d´eformation qui diff`ere de l’´energie vraie dans un rapport 1 + 1−ν2 2. 2.4 Solutions pour des ph´enom`enes de verrouillages L’int´egration de Gauss, quand elle s’effectue avec un nombre limit´e de points, inf´erieur `a celui fix´e par la th´eorie permet d’obtenir un r´esultat exact et peut apporter un gain de comportement appr´eciable en supprimant les contraintes qui peuvent peser sur la repr´esentation des d´eformations, en particulier celle du cisaillement.(voire la figure 6) Par ailleurs, il est clair que les temps de calcul des termes d’une matrice de raideur sont proportionnels au nombre de points retenus pour l’int´egration. 8 3.2.3 La coque h´emisph´erique soumise `a des forces de traction et de compression radiales Figure 12 – Mod`ele et r´esultat attendu pour la coque h´emisph´erique soumise `a des forces de traction et de compression radiales. [9] [12] 15 3.2.4 La coque cylindrique `a extr´emit´e ouverte soumise `a des forces de traction radiales Figure 13 – Mod`ele et r´esultat attendu pour la coque cylindrique `a extr´emit´e ouverte soumise `a des forces de traction radiales. [9] 16 4 Mise en ´evidence du probl`eme de verrouillage r´esolu pour les ´el´ements solide-coques Nous allons pr´esenter les r´esultats des simulations r´ealis´ees qui montreront qu’avec un seul ´el´ement dans l’´epaisseur, nous pouvons supprimer le ph´enom`ene de verrouillage qui peut apparaitre. Pour garantir la fiabilit´e des ´el´ements nous pr´esenteront les r´esultats pour diverses simulations. 4.1 Le cas de la coque h´emisph´erique soumise `a des forces de traction et de compression radiales Figure 14 – Mise en ´evidence du ph´enom`ene de verrouillage dans le cas de la coque h´emisph´erique soumise `a des forces de traction et de compression radiales pour les ´el´ements C3D8R. 17 4.2 Le cas de la plaque annulaire fendue Figure 15 – Comparaison des r´eultats pour la plaque annulaire fendue avec les ´el´ements SHB8 et les ´el´ements C3D8R Figure 16 – Mise en ´evidence du ph´enom`ene de verrouillage dans le cas de la plaque annulaire fendue pour les ´el´ements C3D8R. 18 4.3 Le cas de la poutre en flexion Figure 17 – Mise en ´evidence du ph´enom`ene de non verrouillage dans le cas de la poutre en flexion pour les ´el´ements SHB8. Figure 18 – Mise en ´evidence du ph´enom`ene de non verrouillage dans le cas de la poutre en flexion pour les ´el´ements SHB8. 19 4.4 Analyse des r´esultats Dans tous les cas pr´esent´es ci-dessus, l’´ecart entre la solution attendue et la solution obtenue avec des ´el´ements SHB8 est de moins de 1%. Dans les deux premiers cas, (le cas de la plaque annulaire fendue et le cas de la coque h´emisph´erique soumise `a des forces de traction et de compression radiales), nous observons une convergence de la solution r´ealis´ee avec des ´el´ements C3D8R vers une solution erron´ee. Il y a bien ph´enom`ene de verrouillage. Cependant, dans le cas de la poutre en flexion, nous avions d´emontr´e plus tˆot que la solution obtenue avec les ´el´ement C3D8R converge bien vers la solution th´eorique. Nous remarquons ´egalement que les ´el´ements SHB8 donnent la solution exacte. On peut donc en conclure que les nouveaux ´el´ements SHB8 permettent dans les cas ´etudi´es d’´eviter certains ph´enom`enes de verrouillages. 5 Simulation des mat´eriaux FGM Un mat´eriau `a gradient de propri´et´es fonctionnelles (FGM en anglais) est un mat´eriau dont les propri´et´es physiques ´evoluent selon une variation spatiale progressive de la composition ou de la microstructure. Initialement, l’id´ee fut ´emise en 1984 par M. Niino et al. (National Aerospace Laboratory of Japan) [Niino, 1987] et pr´esent´ee comme une tentative pour r´esoudre les probl`emes pos´es par la pr´eparation de mat´eriaux de barri`ere thermique destin´es aux structures spatiales et aux r´eacteurs de fusion. Il existe de nombreuses techniques permettant d’obtenir des FGM plus ou moins exotiques. Les techniques conventionnelles telles que la forge permettent d’obtenir un gradient de carbone entre l’´epaisseur et le cœur de la pi`ece. Cependant, il existe ´egalement de nombreuses autres m´ethodes plus complexes (impression 3D, s´edimentation, coul´ee par centrifugation,frittage de poudres, etc...) permettant d’allier, par exemple, des aciers et des c´eramiques. [15] [11] L’avantage d’avoir un grand nombre de points d’int´egrations dans l’´epaisseur est de pouvoir ´egalement simuler des mat´eriaux FGM avec un seul et unique ´el´ement dans l’´epaisseur. En effet, on d´efinit un module de Young qui d´epend de l’´epaisseur pour pouvoir simuler un mat´eriau FGM. Ceci ´evite d’avoir `a mettre une multitude d’´el´ements dans l’´epaisseur avec chacun un module de Young diff´erent ce qui augmente d’autant plus la complexit´e du calcul. Pour nos tests, nous utilisons une fonction qui suit une loi puissance entre un m´etal et une c´eramique : E(z)=(EC−EM)fc+EMfc=z t+1 2n z −t 2;t 2(11) avec n un r´eel positif. Ainsi, pour n= 0, le mat´eriaux est totalement c´eramique et pour n→ ∞ le mat´eriaux est totalement m´etallique. 20 Figure 19 – Distribution de la fraction volumique de la phase c´eramique en fonction de l’exposant n. [5] 5.1 D´efinition du module de Young pour la stabilisation Comme dit dans la partie 2.3, la stabilisation des ´el´ements lin´eaires d´epend du module de Youg. Or pour les mat´eriaux FGM, le module de Young est une fonction. Il faut donc d´efinir lequel des modules de Young il nous faut utiliser. Nous d´efinissons donc 4 modules de Young diff´erents : — Le module de Young de l’´el´ement m´etallique pur : Em — Le module de Young de l’´el´ement c´eramique pur : Ec — La moyenne arithm´etique du module de Young de l’´el´ement c´eramique pur et de l’´el´ement m´etallique pur : EM+EC 2 — La moyenne int´egrale du module de Young d´efinie par la fonction E(z) et qui vaut E(z) = Ec−EM 1+n+EM Pour cela, nous allons nous focaliser sur le mod`ele de poutre en flexion. Nous obtenons les r´esultats pr´esent´es dans la Table 1. On remarque que l’´ecart `a l’erreur moyenne est tr`es faible voire n´egligeable compar´e `a l’erreur moyenne (il y a un facteur 10). De plus, une simulation a ´et´e r´ealis´ee avec un module de Young compl`etement aberrant (1 MPa) et le calcul n’a pas abouti. On peut donc en d´eduire que l’on peut utiliser n’importe quel module de Young compris entre le module de Young M´etallique et le module de Young C´eramique. 21 Valeur de n Ec Em moyenne arithm´etique moyenne int´egrale C´eramique (n=0) 0.46 0.40 0.34 0.46 0.2 2.13 2.07 2.01 2.11 0.5 2.81 2.76 2.69 2.78 1 1.09 1.12 1.03 1.03 2 1.89 1.97 2.04 1.99 5 5.29 5.43 5.40 5.36 m´etal (n=100) 0.74 0.71 0.88 0.51 Table 1 –´ Ecart relatif en pourcentage par rapport `a la solution de r´ef´erence. 5.2 Simulations de FGM et comparaisons avec les mod`eles de r´ef´erences 5.2.1 Le cas de la coque cylindrique `a extr´emit´e ouverte soumise `a des forces de traction radiales Figure 20 – Comparaison des ´el´ements SHB8 et des courbes de r´ef´erences pour les simulations de FGM au point A. 22 Figure 21 – Comparaison des ´el´ements SHB8 et des courbes de r´ef´erences pour les simulations de FGM au point B . Figure 22 – Comparaison des ´el´ements SHB8 et des courbes de r´ef´erences pour les simulations de FGM au point C. 23 Figure 23 – Comparaison des ´el´ements SHB8 et des courbes de r´ef´erences pour les simulations de FGM dans le cas de la poutre en flexion 5.2.2 Le cas de la poutre en flexion 5.2.3 Analyse des r´esultats On remarque que pour tous les cas simul´es, les erreurs sont relativement faibles. Pour les cas o`u le mat´eriau est pure (N= 0 et N→ ∞), les erreurs sont tr`es faibles (moins de 1%). Pour le cas N= 1, l’erreur reste faible (de l’ordre de 1%). Pour les cas N= 0.2,N= 0.5 et N= 2, l’erreur est plus grande mais reste correcte (aux alentours de 2%). Cependant, pour N= 5, l’erreur augmente jusqu’`a `a 5%. 5.3 R´eduction de l’erreur identifi´ee Nous avions d´efini le module de Young par la fonction suivante E(z)=(EC−EM)z t+1 2n +EM(12) Dans le fichier Input de Abaqus (.inp) qui est le fichier qui contient le mod`ele `a calculer, nous intervertissons les modules de Young de l’acier et de la c´eramique et nous rempla¸cons le npar 1/n. En posant x=z t+1 2, nous obtenons donc cette nouvelle fonction pour l’´equation du module de Young : E1(z) = (EM−EC)z t+1 2 1 n +EC(13) 24