v6.03.114 FORMA03 - Travaux pratiques de la formation « Utilisation avancée » : charge limite d’une plaque trouée#
Résumé:
Ce test \(\mathrm{2D}\) en contraintes planes quasi-statique permet d’illustrer sur un cas simple les questions relatives à la modélisation élastoplastique ; il met en évidence les effets de structure, de concentration de contraintes, de charge limite.
Il s’agit d’une plaque rectangulaire homogène, trouée en son centre, constituée d’un matériau élastoplastique avec écrouissage isotrope, dont l’état initial est non contraint, qui est soumise à une traction à ses extrémités. On s’intéresse à la solution élastoplastique en charge.
L’objectif du test est de montrer les possibilités de modélisation, l’utilisation de la commandeSTAT_NON_LINE et le post-traitement avec la plate-forme Salome-Meca.
La modélisation A correspond au calcul à force imposée en élasticité. Il illustre l’usage de la commandeSTAT_NON_LINE dans une configuration simplifiée (calcul purement élastique linéaire). Il sert également de référence pour les autres modélisations.
La modélisation B correspond au calcul à force imposée, de référence avec le comportement VMIS_ISOT_TRAC, et illustre l’usage des différents paramètres de la commandeSTAT_NON_LINE, ainsi que les commandes de dépouillement.
La modélisation C explicite le mode opératoire pour effectuer le calcul jusqu’à la charge limite, en utilisant le pilotage du chargement par un déplacement.
La modélisation D est identique à la modélisation C sauf qu’elle utilise un chargement suiveur.
Solution de référence#
Solution élastique#
En élasticité, pour une plaque infinie, comportant un trou de diamètre \(a\) , soumise à un chargement \(P\) selon \(y\) à l’infini, la solution analytique en contraintes planes et coordonnées polaires \((r,\theta )\) est:
En particulier, au bord du trou (\(r=a\) ), on a:
Et le long de l’axe \(x\) :
Numériquement, pour \(P=1\mathit{MPa}\) , et pour une plaque infinie, on a
Point |
Composante |
Calcul |
\(\mathrm{MPa}\) |
\(A\) |
\(\mathit{SIXX}\) |
\({\sigma}_{\theta \theta }(r=a,\theta =\pi /2)\) |
–1 |
\(B\) |
\(\mathit{SIYY}\) |
\({\sigma}_{\theta \theta }(r=a,\theta =0)\) |
3 |
Pour une plaque de dimension finie , les abaques [bib1] permettent d’obtenir le coefficient de concentration de contraintes, et on trouve que pour une traction de \(1\mathrm{MPa}\) , \(\mathrm{SIGYY}\) maximum vaut environ \(3.03\mathit{MPa}\) au point \(B\) .
Solution élastoplastique (charge limite)#
En élastoplasticité, par une approche statique en contraintes planes, on peut obtenir une borne supérieure de la charge limite pour une bande de largeur \(\mathrm{2L}\) finie et de longueur infinie, comportant un trou de largeur \(\mathrm{2a}\) et soumise à une contrainte imposée à l’infini \(p\) :
Ici on obtient comme borne supérieure de la charge limite: \({p}_{\lim}^{-}=0.9\times 270=243\mathit{MPa}\) . (On prend ici \({\sigma}_{y}=270\mathrm{MPa}\) , car la charge limite est identique entre un matériau élastoplastique parfait et un matériau dont la courbe de traction présente une asymptote horizontale à \(270\mathrm{MPa}\) ). Dans ce test (en particulier la modélisation B), on aimerait trouver, par un calcul élastoplastique, une approximation de cette charge limite, sachant que les méthodes analytiques permettent d’en connaître une borne supérieure. Nous prendrons donc la valeur \({p}_{\lim}^{-}\) comme référence.
Références bibliographiques#
Analyse limite des structures fissurées et critères de résistance. F. VOLDOIRE : Note EDF/DER/HI/74/95/26 1995
«Stress concentration factors», Peterson R.E., Wiley, 1974.
Mise en œuvre du TP#
Déroulement du TP#
Il s’agit de mener à bien le calcul élastique en générant la géométrie, le maillage et le fichier de commandes AsterStudy à l’aide de la plate-forme Salome-Meca.
Ce TP permet de:
Mettre en œuvre un calcul non-linéaire standard dans le module AsterStudy: gestion du chargement, des matériaux, du comportement et des paramètres de STAT_NON_LINE ;
Comprendre et mettre en œuvre la notion de pilotage;
Faire des post-traitements «évolués» (en particulier tracer des courbes).
Géométrie#
On créera la face plane du quart supérieur droit de la plaque.
Lancer le module Geometry.
Les principales étapes pour construire cette géométrie sont les suivantes :
Pour définir les contours de la plaque, onpeut, par exemple, utiliser l’outil «Sketcher» (Menu New Entity → Basic→2D Sketch). Il est plus simple de commencerpar le point \(B\) de coordonnées \((10,0)\) . En partant de \(B\) , pour l’arc de cercle, utiliser Element Type(Arc) et Destination(Direction/Perpendicular), et définir le rayon 10 et l’angle et le rayon 90°. On obtient le point \(A\) . Puis utiliser Element Type(Line) et donner les autres points (\(G\) , \(F\) , \(D\) ) par leurs coordonnées absolues. Terminer par Sketch Closure.
On obtient alors un contour fermé (Sketch_1)sur lequel on doit construire une face (Menu New Entity →Build →Face). La géométrie de la plaque est alors complète.
Construire des groupes utiles pour le calcul. Ici on construitles 5groupes des arêtes sur lesquels s’appuieront les conditions aux limites (symétries et chargement) : gauchepour le bord \(\mathit{AG}\) , hautpour le bord \(\mathit{GF}\) etbaspour le bord \(\mathit{BD}\) , droitepour le bord \(\mathit{FD}\) et troupour l’arc \(\mathit{AB}\) . Menu New Entity →Group →Create Group: Sélectionner le type d’entité géométrique (ici la ligne, edge) et sélectionner le bord directement dans la fenêtre graphique, ensuite cliquer sur Add, un numéro d’objet doit alors apparaître. On peut changer le nom du groupe avant de levalider par Apply.
On peut également créer les groupes de nœuds, qui seront utiles pour le post-traitement ou le pilotage (Menu New Entity → Group → Create Group): cinq groupes de sommet \(A\) , \(B\) , \(D\) , \(F\) et \(G\) .
Maillage#
On créera un maillage plan du quart supérieur droit le la plaque, en éléments quadratiques, pour avoir une précision suffisante.
Lancer le module Mesh.
Les principales étapes pour générer le maillage sont les suivantes :
Construire lemaillage (Menu Mesh→ Create Mesh). Sélectionnerla géométrie à mailler Face_1, puis choisir Algorithm → NETGEN 1D-2Den ajoutantHypothesis → NETGEN2D Parameters. Dans cette hypothèse, sélectionnerFineness→ Fine et cocher la case Second Order avant deApply.
Calculer le maillage (Menu Mesh → Compute ) . Une fenêtre des informations de maillage doit apparaître, et on obtient alors un maillage raffiné près du trou avec de grands éléments dans le haut de la plaque.
Pour raffiner ce maillage, cliquer à droite sur le maillage et choisir Edit Mesh, puis éditer les paramètres NETGEN 2D Parameters:
On peut diminuer Max Size en choisissant par exemple \(10\) .
Si on veut raffiner autour du trou, on peut dans l’onglet Local sizes ajouter le groupe trou \(\mathit{AB}\) avec le bouton On edge , et puis il suffit de modifier la valeur associée. En la diminuant (par exemple \(2\) ), le maillage sera raffiné autour du trou.
Calculer lemaillage(MenuMesh→ Compute).
Pour faire passer le maillage de linéaire à quadratique: «Modification -> Convert to/from quadratic».
Créer les groupes demailles correspondants aux groupe géométriques(Menu Mesh→Create Groups from Geometry).Sélectionnertous les groupes géométriques.
Exporter lemaillage au format MED.
Modélisation A#
Caractéristiques de la modélisation#
Calcul élastique d’un quart de la plaque sur un modèle en contraintes planes (C_PLAN). Le chargement est défini dans le § 1.2 . On charge jusqu’à \(p=10\mathit{MPa}\) .
Caractéristiques du maillage#
On utilise un maillage quadratique.
Calcul élastique avec STAT_NON_LINE#
Il s’agit de mener à bien le calcul élastique en générant le fichier de commandes AsterStudy à l’aide de la plate-forme Salome-Meca. La modélisation est C_PLAN. On doit trouver les mêmes résultats qu’en utilisant la commande MECA_STATIQUE.
Réalisation du calcul#
Lancer le module AsterStudy.
Puis en colonne gauche, cliquer surl’onglet Case View.
On définit le fichier de commandes du cas de calcul (cliquer à droite à CurrentCase et choisir Add Stage).
Nota:ajouter des commandes par le menu Commands→ Show All.
Les principales étapes pour la création et le lancement du cas de calcul sont les suivantes:
Lire le maillage au format MED: Commande LIRE_MAILLAGE.
Orienter la normale du bord sur lequel le chargement de traction sera appliqué: Catégorie Mesh / Commande MODI_MAILLAGE / ORIE_PEAU_2D en affectant le groupe haut dans GROUP_MA. On garde le même nom du maillage en utilisant reuse.
Définir les éléments finis utilisés: Commande AFFE_MODELE pouraffecter le phénomène MECANIQUE et la modélisation en contraintes planes 2D (C_PLAN) à tous les éléments.
Définir le matériau: Commande DEFI_MATERIAU.Choisir ELAS et saisir les valeurs du module de Young et du coefficient de Poisson.
Affecter le matériau à tous les éléments: Commande AFFE_MATERIAU.
Affecter les conditions aux limites cinématiques: Commande AFFE_CHAR_CINE/ MECA_IMPO pour la symétrie sur le quart de plaque (les groupes gauche et bas).
Affecterle chargement: Commande AFFE_CHAR_MECA / FORCE_CONTOUR pour la force répartie sur le haut de la plaque. Le plus simple est de définir une charge unitaire (FY=1.0) , que l’on multipliera ensuite par une fonction rampe avec le temps.
Créer une fonction rampe linéaire \(f=t\) pour multiplier le chargement mécanique unitaire: Commande DEFI_FONCTION. Par exemple, elle varie entre \((\mathrm{0.,}0.)\) et \((\mathrm{1000.,}1000.)\)
Créer la discrétisation temporelle à l’aide des commandes DEFI_LIST_REEL et DEFI_LIST_INST.Par exemple, on peut déterminer la fin d’instant à 10s pour correspondre le chargement à \(10\mathit{MPa}\) .
Calculer l’évolution élastique: CommandeSTAT_NON_LINE. On met COMPORTEMENT / RELATION=”ELAS”, la liste d’instant définie précédemment dans INCREMENT, les matériaux dans CHAM_MATER, MODELE et également les conditions aux limites et le chargement (CHARGE + FONC_MULT) dans EXCIT.
Pour lancerle cas de calcul, en colonne gauche, cliquer sur l’onglet History View.
Post-traitementsdes résultats avec Aster_Study#
Pour un calcul non-linéaire, la commande STAT_NON_LINEsort en standard trois champs (selon les options du mot-clef ARCHIVAGE:
Le champ des déplacementsaux nœudsDEPL;
Le champ des contraintesaux pointsde Gauss SIEF_ELGA;
Le champ des variables internesaux points de Gauss VARI_ELGA.
Pour le post-traitement, on propose en plus de calculer à l’aide de CALC_CHAMP avec reuse:
Le champ des contraintesaux nœuds(SIGM_NOEU) par l’option CONTRAINTE.
Les contraintes équivalentes(Von Mises, Tresca, etc.) aux points de Gauss, SIEQ_ELGA par l’option CRITERES.
On propose d’imprimer les résultats au format MED avec IMPR_RESU afin de les visualiser dans Results orParavis.
On propose ensuiteplusieurspost-traitements plusévolués (facutatif):
Extraction de la contrainte SIYYen fonction du déplacement vertical DYpour le point \(G\) :
Extraire le déplacement vertical DYau point \(G\) : commande RECU_FONCTION en utilisant lerésultat, et choisissantle champ DEPLetla composante DY.
Extraire la contrainte SIYYau point \(G\) : commande RECU_FONCTION en utilisant lerésultat, et choisissantle champ SIGM_NOEU etla composante SIYY.
Imprimer la fonction \(\mathit{SIYY}=f(\mathit{DY})\) : Commande IMPR_FONCTION au format XMGRACE (mot-clef COURBE → FONC_Xet FONC_Y).
Extraction de la contrainte SIYY sur l’arête inférieure \(\mathit{BD}\) :
Calculerle champ des contraintes par éléments aux nœuds (SIGM_ELNO): commande CALC_CHAMP / CONTRAINTE.
Extraire une table de contrainte SIYYaux certains points sur l’arête \(\mathit{BD}\) : commande MACR_LIGN_COUPE permet d’extraire dans une tableles composantes d’un champ suivant un chemin donné (mot-clef LIGN_COUPE). Appliquer la commande sur la composante SIYYdu champ SIGM_ELNO, aux10 points sur le chemin \(\mathit{BD}\) (en donnant les coordonnées des points \(B\) et \(D\) ).
Imprimer une courbe à partir d’une table: commande IMPR_TABLE pourimprimer au format XMGRACEla table précédente en filtrant sur le dernier instant (par exemple, FILTRE / NOM_PARA = ‘INST’, VALE= 10). Les paramètres aux axes X et Y de la courbe sont définis par NOM_PARA:abscisse curviligne (ABSC_CURV)et contrainte (SIYY).
Extraction des contraintes sur le b ord du trou \(\mathit{AB}\) :
Extraire la contrainte \({\sigma}_{\theta \theta }\) le long du bord du trou: commande MACR_LIGN_COUPE. Dans le mot-clé LIGN_COUPE, on précise TYPE=ARC et le repère POLAIRE. On définit 10 points sur le bord du trou en donnant les coordonnées du point de départ \(B\) et du centre \(O\) . On filtre au dernier instant INST = 10et produit ainsi une table.
**Remarque: dans le repère* POLAIRE en 2D, La signification des composantes est: DXrayon \(r\) , DYangle \(\theta ` \*.* *Donc* *pour* :math:`{\sigma}_{\theta \theta }\) on a NOM_CMP = SIYY.
Extraire la fonction \({\sigma}_{\theta \theta }=f(s)\) avec l’abscisse curviligne \(s=\theta R\) (\(s\in [0,R\pi /2]\) ): commande RECU_FONCTION.On définit PARA_X avec l’abscisse curviligne ( ABSC_CURV) et PARA_Y avec \({\sigma}_{\theta \theta }\) (SIYY).
En élasticité, on peut comparerce résultat numérique avec la référence analytique (équation ).
On peut créer la formule analytique:commande FORMULE. Elle dépend de la coordonnée curviligne \(S\) (NOM_PARA=”S” / VALE=”p*(1.+2.*cos(2.*S/R))”avec \(p=\mathrm{10MPa}\) et \(R=\mathrm{10mm}\) ).
Créer une liste des valeursréelles de \(S\) de 0 à \({s}_{\max}=\pi R/2\) : Commande DEFI_LIST_REEL.
Interpolation de la formule à partir de la liste de \(S\) : commande CALC_FONC_INTERP.
Commande IMPR_FONCTION pour imprimer au format XMGRACEles courbes de deuxfonctions: celle issue du résultatnumérique \({\sigma}_{\theta \theta }=f(s)\) et celle de la formuleanalytique.
Extraction de la force résultante verticale en haut de la plaque en fonction du déplacement vertical:
Calcul de l’option FORC_NODA avec la commande CALC_CHAMP / FORCE.
Extraire le déplacement verticalau point \(G\) : commande RECU_FONCTION en utilisant lerésultat, et choisissantle champ DEPL, composante DY.
Extraire de la force résultante verticale sur le bord hautde la plaque: commande MACR_LIGN_COUPE.On choisit RESULTATet NOM_CHAM=FORC_NODA. Dansle mot-clef LIGN_COUPE, on définitdes points sur lesquels on souhaite de calculer les résultante: donner les coordonnées de \(G\) et \(F\) , le nombre des points, et RESULTANTE=DY.On produit alorsune table.
Extraire la fonction de la résultante verticale suivant letemps \({\mathit{Resultante}}_{\mathit{DY}}=f(t)\) : commande RECU_FONCTION.
Imprimer lafonction \({\mathit{Resultante}}_{\mathit{DY}}=f(\mathit{DY})\) : commande IMPR_FONCTION au format XMGRACE.
Grandeurs testées et résultats#
On teste la valeur des composantes de contraintes pour le chargement de \(\mathrm{10MPa}\) :
Composante |
Type de référence |
Valeur |
Tolérance |
SIGM_NOEU– \(\mathit{SIYY}\) en \(B\) |
ANALYTIQUE |
\(30\mathrm{MPa}\) |
1.00% |
SIGM_NOEU– \(\mathit{SIXX}\) en \(A\) |
ANALYTIQUE |
\(-10\mathrm{MPa}\) |
2,00% |
Modélisation B#
Caractéristiques de la modélisation#
On fait trois calculs sur un modèle en contraintes planes (C_PLAN) :
Calcul élastique: on charge jusqu’à \(p=10\mathrm{MPa}\) ;
Calcul élasto-plastique: on charge jusqu’à \(p=230\mathrm{MPa}\) ;
Calcul élasto-plastique puis décharge: on charge jusqu’à \(p=230\mathrm{MPa}\) puis on décharge jusqu’à \(p=0\) .
Caractéristiques du maillage#
On utilise le même maillage que la modélisation A (qui comporte 315 TRIA6 et 686 nœuds).
Calcul élasto-plastique avec STAT_NON_LINE#
Il s’agit de mener à bien le calcul élasto-plastique avec écrouissage isotrope donné par une courbe de traction telle que la contrainte uniaxiale tende vers une valeur constante (\(270\mathrm{MPa}\) ).
Il existe donc une charge limite pour cette structure dont une borne supérieure est connue \({p}_{\lim}<243\mathit{MPa}\) . Dans cette modélisation, on ne charge que jusqu’à \(\mathrm{230MPa}\) et on procède à un retour élastique. Le cas de la charge limite sera traité dans le paragraphe suivant.
Préparation du fichier de commande#
Lancer le module AsterStudy.
Puis en colonne gauche, cliquer surl’onglet Case View.
On définit le fichier de commandes du cas de calcul (cliquer à droite à CurrentCase et choisir Add Stage).
Nota:ajouter des commandes par Menu Commands→ Show All.
On définit le fichier de commandes du cas de calcul. Le fichier de commande est très semblable à la modélisation précédente, ci-dessous, en gras, on indique les différences:
Lire le maillage au format MED: Commande LIRE_MAILLAGE.
Orienter la normale du bord sur lequel le chargement de traction sera appliqué: Catégorie Mesh / Commande MODI_MAILLAGE / ORIE_PEAU_2D en affectant le groupe haut dans GROUP_MA. On garde le même nom du maillage en utilisant reuse.
Définir les éléments finis utilisés: Commande AFFE_MODELE pour affecter le phénomène MECANIQUE et la modélisation en contraintes planes 2D ( C_PLAN ) à tous les éléments.
Lire la courbe de traction fournie dans le fichier forma03b.21 : Commande LIRE_FONCTION / NOM_PARA = ‘EPSI’.
D éfinir le matériau : Commande DEFI_MATERIAU / ELAS et TRACTION (affecter la courbe de traction).
Affecter le matériau à tous les éléments: Commande AFFE_MATERIAU.
Affecterles conditions aux limites cinématiqueset le chargement: Commande AFFE_CHAR_CINE/ MECA_IMPO pour la symétrie sur le quart de plaque (les groupes gauche et bas).
Affecter le chargement : Commande AFFE_CHAR_ MECA / FORCE_CONTOUR pour la force répartie sur le haut de la plaque. Le plus simple est de définir une charge unitaire ( FY=1.0 ) , que l’on multipliera ensuite par une fonction rampe avec le temps.
Créer une fonction rampe linéaire \(f=t\) pour multiplier le chargement mécanique unitaire: Commande DEFI_FONCTION. Par exemple, elle varie entre \((\mathrm{0.,}0.)\) et \((\mathrm{1000.,}1000.)\)
Créer la discrétisation temporelle à l’aide des commandes DEFI_LIST_REEL et DEFI_LIST_INST. Par exemple, on peut déterminer 30 pas de temps jusqu’à 300s pour correspondre le chargement maximal à \(300\mathit{MPa}\) . Dans DEFI_LIST_INST, activer la découpe automatique du pas de temps: ECHEC / EVENEMENT=”ERREUR” et ACTION=”DECOUPE”.
Calculer l’évolution élasto-plastique : Commande STAT_NON_LINE. On met COMPORTEMENT / RELATION= “ VMIS_ISOT_TRAC “, la liste d’instant définie précédemment dans INCREMENT, les matériaux dans CHAM_MATER, MODELE et également les conditions aux limites et le chargement (CHARGE + FONC_MULT) dans EXCIT.
Calcul élastique#
Si on indique INST_FIN=10sdans le mot-clef INCREMENT de la commandeSTAT_NIN_LINE, on aura bien appliqué une force de \(10\mathit{MPa}\) , ce qui équivaut strictement au cas élastique.
Voici quelques éléments à vérifier:
Vérifier qu’on retrouve les mêmes résultats que dans la modélisation précédente: déplacements et contraintes.
Vérifier l’indicateur de plasticité (VARI_ELGA) . Il doit être nul partout, idem pour la déformation plastique cumulée (EPSP_ELGA).
Observer le tableau de convergence: nombre d’itérations de Newton.
Faites varier la discrétisation temporelle (nombre de pas de temps).
Comparer la contrainte \({\sigma}_{\theta \theta }\) sur le bord du trou et comparer avec la solution analytique.
Calcul élasto-plastique en charge#
Si on indique INST_FIN=230sdans le mot-clefINCREMENT de la commande STAT_NIN_LINE, on aura bien appliqué une force de \(230\mathit{MPa}\) *.*
On doit avoir plastification de la structure car la contrainte résultante dans une partie de la structure est supérieure à la limite élastique (qui vaut \(\mathrm{200MPa}\) ).
Avec la discrétisation en 30 pas de temps du chargement et l’activation de la découpe du pas de temps, il n’y aura pas de convergence. L’algorithme de Newton échoue. Pour essayer de faire converger, vous pouvez jouer sur plusieurs paramètres:
Augmenter la discrétisation du chargement (attention à ne pas aller trop loin, moins de 100 pas) pour que le calcul ne soit pas trop long !).
Augmenter le nombre d’itérations maximum de Newton qui vaut 10 par défaut ( STAT_NON_LINE / CONVERGENCE / ITER_GLOB_MAXI ).
Activer la recherche linéaire (STAT_NON_LINE / RECH_LINEAIRE ).
Augmenter le nombre de subdivisions possibles du pas de temps dans la gestion de la non-convergence ( DEFI_LIST_INST/ECHEC/SUBD_NIVEAU ).
La combinaison de toutes les techniques précédentes.
Voici une combinaison qui fonctionne dans notre cas:
Discrétisation du chargement en 50 pas;
STAT_NON_LINE / CONVERGENCE / ITER_GLOB_MAXI = 20;
Recherche linéaire activée;
Découpe jusqu’à cinq sous-niveaux ( DEFI_LIST_INST/ECHEC/SUBD_NIVEAU=5 ).
Avec cette combinaison, on obtient la convergence en 358 pas (au lieu des 30 initiaux, à cause des découpes) et plus de 4300 itérations.
Quelques résultats intéressants à observer:
A l’instant final du calcul, pour le chargement maximum, on peut remarquer sur les isovaleurs de déformation plastique cumulée, la localisation des déformations plastiques (variable interne V1) au voisinage de \(B\) . On pourra utiliser la visualisation aux points de Gauss pour visualiser la déformation plastique équivalente cumulée aux lieux où elle est calculée.
Pour un chargement inférieur à \(66,7\mathrm{MPa}\) , il n’y a pas de plastification.
Jusqu’à \(230\mathit{MPa}\) , on est constamment en charge.
La valeur maximum du critère de Von Mises aux points de Gauss est toujours inférieure ou égale à \(270\mathrm{Mpa}\) , ce qui montre que la solution vérifie bien la loi de comportement.
Observons le tableau de convergence à un pas quelconque:
Instant de calcul: 2.297250000000e+02 - Niveau de découpe: 2
Nous sommes à l’instant \(229.725\) , on a découpé deux fois la liste initiale et on converge en 8 itérations de Newton. On remarque:
La recherche linéaire n’a pas été très coûteuse: une itération seulement (par défaut STAT_NON_LINE / RECH_LINEAIRE / ITER_LINE_MAXI=3 ). On voit également qu’il n’y a pas de recherchelinéaire en prédiction.
La convergence est testéesur le critère RESI_GLOB_RELA. Sans changer les valeurs par défaut de la commande STAT_NON_LINE / CONVERGENCE, on doit donc avoir RESI_GLOB_RELA inférieur à \(1.0\times {10}^{-6}\) pour atteindre la convergence.
La matrice tangente n’est calculée qu’en prédiction, sur plusieurs pas de temps, on voit qu’elle est toujours calculée à la première itération. Ceci correspond bien au réglage par défaut de la commande STAT_NON_LINE / NEWTON: REAC_INCR=1 et REAC_ITER = 0.
Il est possible de demander à la commande STAT_NON_LINEd’afficher plus d’informations: suivre la valeur d’un degré de liberté (mot-clef SUIVI_DDL), ou demander à afficher l’endroit où le critère de convergence est le plus mauvais (mot-clef AFFICHAGE/INFO_RESIDU=”OUI”). Ce dernier réglage permet, par exemple, de savoir quel endroit de la structure pilote la convergence.
A la fin du transitoire des informations statistiques sont affichées:
Statistiques sur tout le transitoire.
Nombre de pas de temps : 358
Nombre d’itérations de Newton : 4353
Nombre de factorisations de la matrice : 358
Nombre d’intégrations du comportement : 8715
Nombre de résolutions K.U=F : 4353
Nombre d’itérations de recherche linéaire : 4004
Temps CPU consommé dans le transitoire : 2 m 29 s
dont temps « perdu » dans les découpes : 55.640 s -> la liste d’instant est efficace à 62.8 %
Temps assemblage matrice : 0.790 s
Temps construction second membre : 5.690 s
Temps total factorisation matrice : 1.390 s
Temps total intégration comportement : 2 m 4 s
Temps total résolution K.U=F : 3.030 s
Temps autres opérations : 14.120 s
Nous n’allons pas détailler toutes les informations mais faire quelques remarques:
Pour ce on voit que le poste le plus consommateur est l’intégration de la loi de comportement, bien devant la factorisation et la résolution du système. C’est souvent le cas en 2D, mais c’est surtout lié au fait qu’on utilise une version non complète de Newton. La matrice n’est factorisée qu’une fois par pas de temps (on le voit aussi sur le nombre de factorisations: 358, comme le nombre de pas de temps).
La liste de temps initiale découpé en 50 pas de temps n’était pas la plus efficace: le temps perdu dans le calcul à cause des échecs de convergence et donc de la redécoupe du pas de temps est d’environ un tiers du temps total.
Pour améliorer substantiellement la convergence, il suffit d’activer le Newton complet: STAT_NON_LINE / NEWTON / REAC_ITER=1. Sur 50 pas de temps on obtient les résultats suivants:
Statistiques sur tout le transitoire.
Nombre de pas de temps : 50
Nombre d’itérations de Newton : 152
Nombre de factorisations de la matrice : 152
Nombre d’intégrations du comportement : 202
Nombre de résolutions K.U=F : 152
Temps CPU consommé dans le transitoire : 5.330 s
Temps assemblage matrice : 0.220 s
Temps construction second membre : 0.520 s
Temps total factorisation matrice : 0.650 s
Temps total intégration comportement : 3.280 s
Temps total résolution K.U=F : 0.110 s
Temps autres opérations : 0.550 s
En plus de l’augmentation de la vitesse de convergence (30 fois plus rapide), on observe aucune découpe du pas de temps. On peut découper encore plus grossièrement la liste d’instants.
Calcul élasto-plastique en charge puis décharge#
Nous allons maintenant procéder à la décharge. Pour cela, On définit une nouvelle rampe en forme de chapeau:
\(F=0.\) pour \(t=0\) .
\(F=230.\) pour \(t=230.\) .
\(F=0.\) pour \(t=300.\) .
Puis il faut définir une nouvelle la liste d’instants en rapport: Commande DEFI_LIST_REEL(30 pas de temps jusqu’à \(t=230.\) , puis 10 pas de temps jusqu’à \(t=300.\) ) et de la commande DEFI_LIST_INST en activant la découpe automatique du pas de temps avec ECHEC / EVENEMENT=”ERREUR” et ACTION=”DECOUPE”).
Faites un nouveau calcul (unenouvelle commandeSTAT_NON_LINE) en prenant INST_FIN =300. dans le mot-clef INCREMENT de la commandeSTAT_NON_LINE , et en utilisant la nouvelle rampe et la nouvelle liste d’instant.
Si on reprend la stratégie permettant de mener le calcul jusqu’au bout en plasticité (c’est-à-dire jusqu’à \(p=230\mathit{MPa}\) en utilisant STAT_NON_LINE / NEWTON / REAC_ITER=1 ), en minimisant la découpe temporelle (exercice précédent), on s’aperçoit que cette stratégie doit être améliorée sur la partie décharge plastique car elle ne converge pas avec ces réglages.
On rappelle que la décharge se fait de manière élastique et crée une déformation inélastique lorsque le chargement est nul. Pour que le calcul converge, il est donc nécessaire d’activer la matrice élastique en prédiction (STAT_NON_LINE /NEWTON / PREDICTION=”ELASTIQUE”).
Il s’agit de mener à bien le calcul élasto-plastique avec écrouissage isotrope donné par une courbe de traction telle que la contrainte uniaxiale tende vers une valeur constante (\(270\mathrm{MPa}\) ).
Il existe donc une charge limite pour cette structure dont une borne supérieure est connue: \({p}_{\lim}<243\mathit{MPa}\) . Dans cette modélisation, on va montrer comment mener le calcul au delà de cette charge limite grâce au pilotage.
Grandeurs testées et résultats#
On teste la valeur des composantes de contraintes pour le calcul élastique (chargement de \(\mathrm{10MPa}\) ), on doit trouver la même chose que dans la modélisation A, soit:
Composante |
Type de référence |
Valeur |
Tolérance |
SIGM_NOEU– \(\mathit{SIYY}\) en \(B\) |
AUTRE_ASTER |
Identique A |
|
SIGM_NOEU– \(\mathit{SIXX}\) en \(A\) |
AUTRE_ASTER |
Identique A |
|
On teste la valeur des composantes de contraintes et des variables internes pour le calcul élasto-plastique (à l’instant correspond au chargement de \(\mathrm{230MPa}\) ):
Composante |
Type de référence |
Tolérance |
SIGM_NOEU– \(\mathit{SIYY}\) en \(B\) |
NON_REGRESSION |
1,00E-006% |
SIGM_NOEU– \(\mathit{SIXX}\) en \(A\) |
NON_REGRESSION |
1,00E-006% |
VARI_NOEU– \(\mathit{V1}\) en \(B\) |
NON_REGRESSION |
1,00E-006% |
VARI_NOEU– \(\mathit{V2}\) en \(B\) |
NON_REGRESSION |
1,00E-006% |
VARI_NOEU– \(\mathit{V1}\) en \(A\) |
NON_REGRESSION |
1,00E-006% |
VARI_NOEU– \(\mathit{V2}\) en \(A\) |
NON_REGRESSION |
1,00E-006% |
On teste la valeur des composantes de contraintes et des variables internes pour le calcul élasto-plastique avec décharge (à l’instant correspond au déchargement final de \(0\) ):
Composante |
Type de référence |
Tolérance |
SIGM_NOEU– \(\mathit{SIYY}\) en \(B\) |
NON_REGRESSION |
1,00E-006% |
SIGM_NOEU– \(\mathit{SIXX}\) en \(A\) |
NON_REGRESSION |
1,00E-006% |
VARI_NOEU– \(\mathit{V1}\) en \(B\) |
NON_REGRESSION |
1,00E-006% |
VARI_NOEU– \(\mathit{V2}\) en \(B\) |
NON_REGRESSION |
1,00E-006% |
VARI_NOEU– \(\mathit{V1}\) en \(A\) |
NON_REGRESSION |
1,00E-006% |
VARI_NOEU– \(\mathit{V2}\) en \(A\) |
NON_REGRESSION |
1,00E-006% |
Modélisation C#
Caractéristiques de la modélisation#
On fait deux calculs sur un modèle en contraintes planes (C_PLAN) :
Calcul élasto-plastique: on charge jusqu’à \(p=243\mathit{MPa}\) ;
Calcul élasto-plastique: on fait un calcul par pilotage au-delà de la charge limite;
Caractéristiques du maillage#
On utilise le même maillage que la modélisation B qui comporte 315 TRIA6 et 686 nœuds.
Calcul avec charge limite#
Détection « manuelle » de la charge limite#
Lancer le module AsterStudy.
Puis en colonne gauche, cliquer surl’onglet Case View.
On définit le fichier de commandes du cas de calcul (cliquer à droite à CurrentCase et choisir Add Stage).
Nota:ajouter des commandes par Menu Commands→ Show All.
Le fichier de commande est presque identiqueà la modélisationprécédente.La seule différence est que nous proposons ici d’activer la gestion automatique de la liste d’instants, c’est-à-dire que la discrétisation temporelle est entièrement gérée par STAT_NON_LINE.
Lire le maillage au format MED: Commande LIRE_MAILLAGE.
Orienter la normale du bord sur lequel le chargement de traction sera appliqué: Catégorie Mesh / Commande MODI_MAILLAGE / ORIE_PEAU_2D en affectant le groupe haut dans GROUP_MA. On garde le même nom du maillage en utilisant reuse.
Définir les éléments finis utilisés: Commande AFFE_MODELE pouraffecter le phénomène MECANIQUEet la modélisation en contraintes planes 2D ( C_PLAN) à tous les éléments.
Lire la courbe de traction fournie dans le fichier forma03c.21 Commande LIRE_FONCTION.
Définir le matériau: Commande DEFI_MATERIAU /ELAS et TRACTION.
Affecter le matériau: Commande AFFE_MATERIAU.
Affecter les conditions aux limites cinématiques: Commande AFFE_CHAR_CINE/ MECA_IMPO pour la symétrie sur le quart de plaque (les groupes gauche et bas).
Affecter le chargement : Commande AFFE_CHAR_ MECA / FORCE_CONTOUR pour la force répartie sur le haut de la plaque. Le plus simple est de définir une charge unitaire ( FY=1.0 ) , que l’on multipliera ensuite par une fonction rampe avec le temps.
Créer une fonction rampe linéaire \(f=t\) pour multiplier le chargement mécanique unitaire: Commande DEFI_FONCTION. Par exemple, elle varie entre \((\mathrm{0.,}0.)\) et \((\mathrm{1000.,}1000.)\)
Activer la gestion automatique du pas de temps avec la commande DEFI_LIST_INST / METHODE=”AUTO” et DEFI_LIST avec PAS_MINI=1.e-6, PAS_MAXI=100 et VALE=(0., 50., 243) qui donne les trois instants de passage obligatoire de la liste automatique.
Calculer l’évolution de l’élasto-plastique: CommandeSTAT_NON_LINE/ COMPORTEMENT / RELATION=”VMIS_ISOT_TRAC” avec la liste d’instant définie précédemment.
On propose d’abord de constater «à la main» la charge limite. La charge limite correspond au moment où un point de Gauss atteint la valeur équivalente de Von Mises d’environ \(\mathrm{270MPa}\) . Par le calcul analytique (voir § 2.2 ), on a déterminé une borne supérieure du chargement qui provoque une plastification à cette valeur limite de \(\mathrm{270MPa}\) (autour du trou).
Pour cela, procéder comme dans la modélisation précédente mais faites un chargement au delà de \(p=\mathrm{230MPa}\) . Dans un premier temps, une valeur de \(p=\mathrm{245MPa}\) est une bonne référence.
On voit qu’au delà d’un certain chargement, la matrice tangente devient singulière: c’est le signe qu’on a atteint la charge limite et donc une tangente horizontale sur la courbe de traction. Le code va tenter de découper le pas de temps pour aller au-delà de ce point limite. Par dichotomie (découpe du pas de temps), il va s’approcher de la valeur de chargement limite. Selon les maillages, on trouve que \({p}_{\lim}\approx \mathrm{243MPa}\) .
La cause de la convergence difficile est bien la proximité de la charge limite. C’est pourquoi il faut subdiviser le pas de temps. On peut s’en rendre compte par la valeur du chargement et par la courbe contrainte-déplacement en haut de la structure : on peut constater que pour \(p=240\mathit{MPa}\) la charge limite n’est pas complètement atteinte (pas d’asymptote horizontale) mais que l’on s’en rapproche. Les isovaleurs de \(p\) montrent une zone de concentration de déformation plastique (assimilable à une ligne de glissement) inclinée de \(53°\) environ par rapport à la verticale, allant du point \(B\) au bord droit. Ceci correspond assez bien à la théorie qui dit que les lignes de glissement sont inclinées de \(54,44°\) (voir [bib2]). On a ici bien sûr une approximation de la ligne de glissement qui est en théorie d’épaisseur nulle. On peut aussi relever le déplacement vertical maximum suivant \(Y\) du point \(G\) . Il est d’environ \(\mathrm{5.7mm}\) .
Calcul au-delà de la charge limite par pilotage#
La meilleure solution si on souhaite atteindre le chargement limite et même aller au-delà (par la résolution d’un problème élastoplastique incrémental) est d’utiliser le pilotage de la contrainte imposée par le déplacement d’un point. C’est ce qu’on propose ici.
On pourra utiliser par exemple le déplacement DY du point \(G\) pour piloter la contrainte \({\sigma}_{yy}\) imposée sur \(\mathit{FG}\) . On l’augmentera jusqu’à \(6\mathit{mm}\) par exemple (dans le calcul précédent, on a observé que le déplacement maximum suivant \(Y\) du point \(G\) était d’environ \(\mathrm{5.7mm}\) , \(\mathrm{6mm}\) est donc bien au-delà de la charge limite).
On prendra un coefficient égal à 1. On utilisera donc un temps fictif \(t\) tel que \(\Delta t=\Delta {U}_{Y}(G)\times 1\) . C’est-à-dire que le temps fictif varie ici entre \(0\) et \(6s\) (pour représenter un déplacement de pilotageDY du \(G\) entre \(0\) et \(6\mathrm{mm}\) ).
Remarque: en alternative à ce type de calcul, une méthode de calcul de la charge limite est disponible dans Code_Aster(en \(\mathrm{3D}\) , axisymétrique et déformation plane): elle utilise un matériau de Norton-Hoff, des éléments quasi-incompressibles et des méthodes directes d’analyse limite fournissant un encadrement de la charge limite (cf. les documents [U2.05.04] et [R7.07.01] ).
Lancer le module AsterStudy.
Puis en colonne gauche, cliquer surl’onglet Case View.
On définit le fichier de commandes du cas de calcul (cliquer à droite à CurrentCase et choisir Add Stage)
Nota:ajouter des commandes par Menu Commands→ Show All.
Les principales étapes pour la création et le lancement du cas de calcul sont les suivantes:
Lire le maillage au format MED: Commande LIRE_MAILLAGE.
Orienter la normale du bord sur lequel le chargement de traction sera appliqué: Catégorie Mesh / Commande MODI_MAILLAGE / ORIE_PEAU_2D en affectant le groupe haut dans GROUP_MA. On garde le même nom du maillage en utilisant reuse.
Définir les éléments finis utilisés: Commande AFFE_MODELE pouraffecter le phénomène MECANIQUEet la modélisation en contraintes planes 2D ( C_PLAN) à tous les éléments.
Lire la courbe de traction fournie dans le fichier forma03c.21 Commande LIRE_FONCTION.
Définir le matériau: Commande DEFI_MATERIAU /ELAS etTRACTION.
Affecter le matériau à tous les éléments: Commande AFFE_MATERIAU.
Affecter les conditions aux limites cinématiques: Commande AFFE_CHAR_CINE/ MECA_IMPO pour la symétrie sur le quart de plaque (les groupes gauche et bas).
Affecter le chargement : Commande AFFE_CHAR_ MECA / FORCE_CONTOUR pour la force répartie sur le haut de la plaque. On définit une charge unitaire (FY=1.0) qui sera pilotée dans STAT_NON_LINE.
Créez la discrétisation temporelle à l’aide des commandes DEFI_LIST_REEL ( 10 pas de temps jusqu’à \(t=4s\) ) et DEFI_LIST_INST en activant la découpe automatique du pas de temps: ECHEC / EVENEMENT=”ERREUR” et ACTION=”DECOUPE”) .
Modifierle type du chargement (force répartie) : dans la commande STAT_NON_LINE / EXCIT, utiliser TYPE_CHARGE=”FIXE_PILO” pour la force appliquée. Et activer le pilotage: mot-clef facteur PILOTAGE:
TYPE=”DDL_IMPO”.
COEF_MULT=1.0.
GROUP_NO=”G”.
NOM_CMP=”DY”.
Calculer l’évolution de l’élasto-plastique:Commande STAT_NON_LINE/ COMPORTEMENT / RELATION=”VMIS_ISOT_TRAC”.
Observez dans le fichier «message» la valeur du paramètre ETA_PILOTAGE, \(\eta ` . Pour connaître la charge issue du pilotage, il suffit de faire :math:`{F}_{\mathit{pilote}}=\eta \times {F}_{y}\) avec \({F}_{y}=1.0\) , charge unitaire qu’on impose. On obtient en principe une bonne approximation de la charge limite (par valeur supérieure) que l’on pourra comparer à la borne supérieure analytique (\(243\mathrm{MPa}\) ).
On pourra tracer, en poursuite, la courbe force résultante-déplacement en \(G\) en fonction du temps.
Extraction du paramètre de pilotage:
Extraire les valeurs du paramètre ETA_PILOTAGEdans le résultat: commande RECU_FONCTION/NOM_PARA_RESU surle paramètre ETA_PILOTAGE.
Imprimer la fonction \(\text{ETA\_PILOTAGE}=f(t)\) : Commande IMPR_FONCTION.
On peut visualiser avec Salomé, la déformée, les isovaleurs de contraintes \(\mathrm{SIYY}\) et de la déformation plastique équivalente cumulée. A l’instant 6, on peut remarquer sur les isovaleurs de déformation plastique cumulée, la localisation des déformations au voisinage de \(B\) .
Grandeurs testées et résultats#
On teste la valeur de la charge limite pour le calcul élastoplastique sans pilotage, par rapport à la solution analytique de la borne minimale:
Composante |
Type de référence |
Valeur |
Tolérance |
SIGM_NOEU– \(\mathit{SIYY}\) en \(G\) |
ANALYTIQUE |
\(243\mathit{MPa}\) |
1.00% |
On teste la valeur de la charge limite pour le calcul élastoplastique avec pilotage (déplacement imposé vertical de \(6\mathit{mm}\) au point \(G\) ), par rapport à la solution analytique de la borne minimale, de la solution sans pilotage:
Composante |
Type de référence |
Valeur |
Tolérance |
SIGM_NOEU– \(\mathit{SIYY}\) en \(G\) |
ANALYTIQUE |
\(243\mathit{MPa}\) |
1.00% |
SIGM_NOEU– \(\mathit{SIYY}\) en \(G\) |
AUTRE_ASTER |
\(243.05\mathit{MPa}\) |
0.30% |
ETA_PILOTAGEen INST=6 |
ANALYTIQUE |
\(243\mathit{MPa}\) |
1.00% |
Modélisation D#
Caractéristiques de la modélisation#
Identique à la modélisation C, cette modélisation utilise un chargement suiveur, SUIV_PILO au lieu de FIXE_PILO (type de chargement: PRES_REP au lieu de FORCE_CONTOUR).
Si FIXE_PILO, le chargement est toujours fixe (indépendant de la géométrie).
Si SUIV_PILO, Le chargement est dit «suiveur», c’est-à-dire qu’il dépend de la valeur des inconnues: par exemple, la pression, étant un chargement s’appliquant dans la direction normale à une structure, dépend de la géométrie actualisée de celle-ci, et donc des déplacements.
Caractéristiques du maillage#
On utilise le même maillage que la modélisation C qui comporte 315 TRIA6 et 686 nœuds.
Grandeurs testées et résultats#
On teste la valeur de la charge limite pour le calcul élastoplastique sans pilotage, par rapport à la solution analytique de la borne minimale:
Composante |
Type de référence |
Valeur |
Tolérance |
SIGM_NOEU– \(\mathit{SIYY}\) en \(G\) |
ANALYTIQUE |
\(243\mathit{MPa}\) |
1.00% |
On teste la valeur de la charge limite pour le calcul élastoplastique avec pilotage (déplacement imposé vertical de \(\mathrm{6mm}\) au point \(G\) ), par rapport à la solution analytique de la borne minimale de la solution sans pilotage:
Composante |
Type de référence |
Valeur |
Tolérance |
SIGM_NOEU– \(\mathit{SIYY}\) en \(G\) |
ANALYTIQUE |
\(243\mathit{MPa}\) |
1.00% |
SIGM_NOEU– \(\mathit{SIYY}\) en \(G\) |
AUTRE_ASTER |
\(243.05\mathit{MPa}\) |
0.30% |
ETA_PILOTAGEen INST=6 |
ANALYTIQUE |
\(246\mathit{MPa}\) |
1.00% |
Synthèse des résultats#
Ce test permet de montrer comment mener le calcul d’une structure élasto-plastique et son dépouillement, et en particulier de mettre en évidence le bénéfice à utiliser le pilotage pour un problème de charge limite. On peut retenir de ce test quelques idées:
Même en dehors d’un comportement élasto-plastique parfait, il peut exister une charge limite. C’est le cas avec toutes les courbes de traction réelles. Il faut alors adapter la méthode de résolution à la solution mécanique et par exemple utiliser le pilotage;
Le découpage en petits incréments de charge est souvent nécessaire pour intégrer correctement la relation de comportement. Cela peut aider aussi à la convergence, il est donc conseillé d’utiliser le redécoupage automatique du pas de temps;
La recherche linéaire peut être utilisée pour aider à la convergence, ainsi que la subdivision automatique des pas de temps. En cas de décharge, la prédiction élastique est une solution efficace.