u2.05.06 Réalisation de calculs d’endommagement en quasi-statique (rupture fragile)#
Résumé:
Ce document a pour but de conseiller l’utilisateur de Code_Aster pour la réalisation de simulation à l’aide de modèles d’endommagement et notamment des techniques à mettre en œuvre pour avoir une solution de bonne qualité et de manière robuste.
Choix du modèle d’endommagement#
Pour représenter la dégradation d’un matériau voire sa rupture, l’une des méthodes possibles (notamment lorsqu’on ne connaît pas le mode de ruine) est d’utiliser une loi de comportement «adoucissante», c’est-à-dire telle que, une fois passé un seuil (en contrainte ou en déformation), la contrainte diminue lorsque la déformation augmente. Ce type de comportement peut être obtenu :
en introduisant une variable d’endommagement \(D\) comprise entre \(0\) et \(1\) (= mécanique de l’endommagement telle qu’introduite par Kachanov ou Lemaitre): lois ENDO_FRAGILE, ENDO_ISOT_BETON, ENDO_ORTH_BETON,MAZARS dans Code_Aster;
en utilisant des modèles de plasticité avec un «écrouissage» négatif comme par exemple les lois BETON_DOUBLE_DP, DRUCK_PRAGER, mais aussi les différentes lois de sol disponibles dans Code_Aster ;
en introduisant une variable telle que la porosité couplée à de la plasticité pour la rupture ductile: loi de ROUSSELIER dans Code_Aster .
Quelle que soit la méthode employée, toutes ces lois adoucissantes «locales» conduisent à la même difficulté : il arrive un moment où le problème devient mal posé et la solution obtenue devient fortement dépendante du maillage. En effet, l’endommagement se localise dans une bande ayant pour épaisseur un seul élément (d’où une énergie de fissuration qui tend vers zéro quand on raffine le maillage) et il y a dépendance du «chemin» de fissuration à la topologie du maillage. La pertinence de ces calculs (bien qu’encore assez répandus dans la littérature) est donc très contestable, et le fait de faire dépendre l’énergie dissipée de la taille de maille n’apporte qu’une solution très partielle.
Pour contourner ce problème, il est couramment admis qu’il faut introduire dans le problème à résoudre une longueur caractéristique, qui va contrôler l’épaisseur de la zone endommagée indépendamment du maillage et permettre à nouveau d’avoir convergence de la solution lorsqu’on raffine le maillage.
Plusieurs méthodes de régularisation ont été mises au point ces dernières années qui ont chacune leurs avantages et leurs limites. Dans Code_Aster , les méthodes disponibles s’appuient toutes sur une notion de gradient (par opposition aux méthodes de type intégrales, largement répandues dans la littérature également).
On distingue les méthodes qui introduisent :
le gradient de la variable d’endommagement (modèle GRAD_VARI, confer documentation [R5.04.01]);
les modélisations second gradient et second gradient de dilatation qui introduisent une énergie dépendant totalement, ou en partie, des composantes du gradient de la déformation (modèles 2DGet DIL, cf. doc [R5.04.03]);
les modélisations à gradient de gonflement (modèle INCO_UPG, cf. doc [R3.06.08]).
Dans tous les cas, le principe est de pénaliser d’un point de vue énergétique la localisation de l’endommagement.
Le tableau ci-dessous récapitule pour les différentes lois de comportement, quelle modélisation régularisée est disponible (et valide) dans Code_Aster .
Loi de comportement |
Modélisation |
ENDO_FRAGILE |
GRAD_VARI / 2DG |
ENDO_ISOT_BETON |
GRAD_VARI / 2DG |
MAZARS |
2DG |
ENDO_ORTH_BETON |
2DG |
DRUCK_PRAGER |
2DG/DIL |
ROUSSELIER |
INCO_UPG |
Tableau 1: correspondance loi d’endommagement/modélisation non-locale
Remarques :
Il n’existe pas de version non-locale du modèle BETON_DOUBLE_DP. La version implantée inclut toutefois une régularisation de type Hillerborg pour éviter le problème de l’énergie qui tend vers zéro quand la taille des éléments tend vers zéro.
Toutes les lois de comportement peuvent être utilisées avec les modélisations 2DG et DIL. Cependant, la modélisation 2DG n’a pour l’instant été utilisée qu’avec les lois de sols et la modélisation DIL, n’a de sens que pour régulariser l” «endommagement volumique» : elle est donc bien adaptée pour traiter le cas des matériaux dilatants et donc des sols. Outre DRUCK_PRAGER, il est donc possible d’utiliser les lois CAM_CLAY (qui est un cas particulier de la loi de Hujeux) et HUJEUX (cf. thèse d’Alexandre Foucault et CR-AMA-09-154), VISC_DRUCK_PRAG (voir note H-T64-2009-03498) et LETK, BARCELONE, CJS, HOEK-BROWN ; mais le retour d’expérience avec ces lois est encore faible.
La modélisation INCO_UPG est pour l’instant mal maîtrisée, et doit faire l’objet d’études et de développement complémentaires. On déconseille son utilisation.
Le coût CPU des modélisations régularisées est important, d’une part parce qu’elles nécessitent de mailler assez finement (voir § 3 ) et d’autre parce qu’elles introduisent des degrés de liberté supplémentaires rendant les matrices à inverser beaucoup plus grandes et moins creuses.
Le maillage#
Pour un calcul mécanique avec endommagement local, on peut utiliser indifféremment des maillages linéaires ou quadratiques. Par contre, il est conseillé d’avoir des tailles de mailles assez homogènes pour que la dissipation le soit également (rappel: la dissipation en local est liée à la taille de la maille). Il est également conseillé d’avoir procédé à l’identification des paramètres post-pic sur des essais utilisant à peu près les mêmes tailles de maille.
La plupart des modèles non-locaux s’appuient sur des maillages quadratiques. Le tableau suivant récapitule les éléments disponibles pour chaque modélisation.
2D |
AXI |
3D |
|
GRAD_VARI |
TRIA6 QUAD8 |
TRIA6 QUAD8 |
TETRA10 HEXA20 PENTA15 PYRAM13 |
2DG |
TRIA7 QUAD9 |
||
DIL |
TRIA7 QUAD9 TRIA6 QUAD8 |
-[*]_ |
|
INCO_UPG |
TRIA6 QUAD8 |
TRIA6 QUAD8 |
TETRA10 HEXA20 |
Tableau 2: mailles disponibles en fonction des modélisations
Concernant la finesse du maillage, on rappelle que pour que la régularisation du problème soit effective, il faut au moins trois mailles voire dix mailles (pour la modélisation GRAD_VARI) pour une bande d’endommagement. Ainsi, si la taille de la zone endommagée fait \(1\mathit{cm}\) , il faut que dans cette zone, la taille des mailles soient inférieures à \(3\mathit{mm}\) voire \(1\mathit{mm}\) pour GRAD_VARI.
La résolution du problème avec STAT_NON_LINE#
Le but de ce chapitre est de donner des conseils sur la manière de résoudre les problèmes adoucissants avec STAT_NON_LINE. On propose ici une certaine gradation des outils en fonction de leur originalité, de leur efficacité, de leurs difficultés de mise en œuvre,… permettant à un nouvel utilisateur de se faire sa propre expérience. En fonction de la loi de comportement utilisée, de la modélisation, de la structure, il pourra être plus judicieux d’utiliser directement l’une ou l’autre des méthodes (notamment le pilotage ou la recherche linéaire mixte). Toutefois, il est possible aussi d’augmenter légèrement le critère de convergence (RESI_GLOBA_RELA = \({10}^{-6}\) par défaut) afin de réduire le temps de calcul et d’aider à la convergence pour des structures fortement endommagées. Il ne faut pas dépasser \({10}^{-4}\) .
Newton-Raphson simple (METHODE=”NEWTON”)#
Il s’agit de la méthode la plus simple à utiliser et qui garantit la qualité de la solution obtenue, à condition bien sûr d’utiliser un critère de convergence suffisamment petit. Cependant :
les problèmes étant adoucissants, il est recommandé de travailler en déplacement imposé plutôt qu’en effort imposé. Sinon, dès que l’effort global à appliquer tant à diminuer, il n’y aura plus aucun moyen de converger.
Dans le cas de matériau fragile, on observe souvent une instabilité de la solution, c’est-à-dire une rupture brutale du matériau (snap-back dans la réponse globale force-déplacement). Dans ce cas, il y a peu de chance que Newton seul franchisse cette instabilité (sauf parfois dans le cas de petites instabilités ou du passage à l’état totalement cassé), et il convient de basculer vers la recherche linéaire mixte ou surtout le pilotage.
Figure 1: exemple de réponse instable
On recommande en général, d’activer le découpage automatique du pas de temps (plusieurs niveaux, voir la commande DEFI_LIST_INST) et d’utiliser la matrice tangente (réactualisée pour les lois types ENDO_** , pas forcément réactualisée pour Mazars) mais de permettre éventuellement le basculement vers la matrice sécante lorsqu’on a redécoupé beaucoup le pas de temps (mot-clé PAS_MINI_ELAS). Une fois que la structure est fortement endommagée, la matrice sécante est parfois plus efficace. Signalons que pour les lois de sols, il est inutile de réactualiser la matrice sécante, car on utilise la matrice élastique (REAC_ITER_ELAS=0).
Avec la matrice tangente, on peut autoriser entre 10 et 30 itérations de Newton selon la loi de comportement utilisé (ITER_GLOB_MAXI). En revanche, avec la matrice sécante, la convergence est beaucoup plus lente et il faut autoriser beaucoup plus d’itérations (ITER_GLOB_ELAS100).
Newton-Raphson + recherche linéaire#
Dans certains cas, la recherche linéaire classique (mot-clé RECH_LINEAIRE, METHODE=”CORDE”) peut aider à la convergence. Mais la méthode MIXTE est nettement plus robuste.
Newton-Raphson + recherche linéaire méthode mixte#
La recherche linéaire méthode mixte (mot-clé RECH_LINEAIRE, METHODE=”MIXTE”) est capable d’aller chercher des solutions plus loin que le bassin d’attraction de Newton et va ainsi permettre de franchir des snap-back. Pour éviter d’aller chercher des solutions trop éloignées (les problèmes d’endommagement n’admettent pas une solution unique), on conseille :
d’utiliser de petits pas de temps,
de ne pas activer le redécoupage du pas de temps (ou très peu), car en présence d’une grande instabilité la taille du pas de temps importe peu
d’autoriser un grand nombre d’itérations de Newton (ITER_GLOB_MAXI=1000),
d’utiliser la matrice tangente réactualisée (REAC_ITER=1),
d’autoriser un grand nombre d’itérations de recherche linéaire pour laisser à la méthode l’opportunité de trouver la solution dans la direction de descente (ITER_LINE_MAXI=50)
C’est en général la méthode la plus robuste, mais elle peut s’avérer relativement coûteuse, puisqu’on autorise beaucoup d’itérations. Par contre, elle n’est en général plus très efficace quand la structure est totalement cassée. Dans ce cas, elle a tendance à faire beaucoup d’itérations et peut parfois augmenter artificiellement la taille de la bande de localisation.
En cas de grands snap-back, correspondant à une rupture brutale de l’éprouvette, identifiable sur la courbe force-déplacement par une diminution brutale de l’effort appliqué, il est conseillé d’utiliser le pilotage, qui va chercher une solution continue afin de s’assurer de la solution de la qualité.
Pilotage#
Le pilotage ou méthode de continuation est une méthode qui permet de suivre la réponse d’une structure en cas d’instabilité (mot-clé PILOTAGE, confer documentation [R5.03.80]). Pour ce faire, l’intensité du chargement devient une nouvelle inconnue du problème. En conséquence, cette méthode n’est pas applicable lorsque le problème dépend explicitement du temps (endommagement couplé avec du fluage, chargement thermique, etc.).
Deux modes de pilotage adaptés aux problèmes d’endommagement existent dans Code_Aster. Il s’agit d’une part du pilotage par incrément de déformation (TYPE=”DEFORMATION”), et d’autre part, du pilotage par prédiction élastique (TYPE=”PRED_ELAS”). Le pilotage par prédiction élastique est en général plus performant, mais il n’est pas générique (contrairement au pilotage par déformation); il est disponible pour les lois ENDO_FRAGILE, ENDO_ISOT_BETON et ENDO_ORTH_BETON.
Préliminaire valable pour les 2 types de pilotages#
Lorsque vous décidez de piloter un chargement, vous ne maîtrisez plus son évolution au cours du temps. En effet, l’intensité du chargement devient une inconnue du problème : elle peut augmenter ou diminuer au cours des incréments. Vous pouvez néanmoins fixer la rapidité avec laquelle il va évoluer, grâce aux incréments de temps et au paramètre COEF_MULT (voir plus loin).
Vous ne pouvez piloter qu’un chargement ou alors plusieurs chargements qui sont proportionnels et donc définis dans la même occurrence de AFFE_CHAR_MECA et multipliée par la même fonction multiplicatrice. (ex. traction + cisaillement simultanés). Ainsi, si les efforts sont décomposés en deux termes : les efforts fixes \({F}_{\mathrm{fixe}}\) et les efforts pilotés \({F}_{\mathrm{pilo}}\) , le chargement vaudra \(F={F}_{\mathit{fixe}}+\eta {F}_{\mathit{pilo}}\) où \(\eta\) sera l’intensité de l’effort piloté.
Admettons que le chargement soit composé des conditions aux limites et d’un effort de traction.
TRAC = AFFE_CHAR_MECA(DDL_IMPO=(_F(GROUP_MA=”HAUT”, DY = 0.0002)
FCT = DEFI_FONCTION( VALE = 0. , 0., 1., 1.)
Si TRAC est le chargement à piloter, dans STAT_NON_LINE, il suffit de remplacer :
EXCIT=_F(CHARGE= TRAC, FONC_MULT=FCT)
par
EXCIT=_F(CHARGE= TRAC, TYPE_CHARGE=”FIXE_PILO”).
Au cours des pas de temps, vous pouvez suivre l’intensité du chargement grâce à la valeur ETA_PILOTAGE indiquée dans les tableaux de convergence de STAT_NON_LINE (voir exemple ci-dessous). Dans notre cas, à la fin du pas de temps considéré, le déplacement vaudra:
\(\mathit{dy}=0,0002\times \text{eta\_pilotage}=0,0002\times 8,66245\times {10}^{-3}\)
Vous pouvez contrôler les bornes du chargement en fixant dans le mot-clé PILOTAGE, les paramètres ETA_PILO_MAXet ETA_PILO_MIN. Quand ces valeurs maximales ou minimales sont atteintes, le calcul s’arrêtera, quel que soit l’instant. ETA_PILO_MIN est surtout utile quand on travaille à effort imposé pour arrêter le calcul quand l’effort devient très petit (parce que la structure est cassée) et ETA_PILO_MAX quand on travaille à déplacement imposé, lorsqu’on a atteint le déplacement souhaité.
On peut également renseigner les valeurs ETA_PILO_R_MIN et ETA_PILO_R_MAX qui fixe les valeurs physiques pour ETA. ETA_PILO_R_MIN est notamment indispensable pour ENDO_FRAGILE qui présente le même comportement en traction et en compression sous peine de sauter d’une solution à l’autre.
Enfin, il peut être intéressant de préciser le GROUP_MA sur lequel effectuer le pilotage. En effet, si l’endommagement est localisé, cela peut permettre de limiter les temps de calculs. En revanche, c’est indispensable dans le cas où il y a plusieurs lois de comportement dans le modèle, dont certaines non «pilotables» (pour pilotage type PRED_ELAS).
En général, le calcul avec pilotage demande plus de pas de temps qu’un calcul direct, car on limite l’évolution de l’endommagement. Si le calcul s’est arrêté parce qu’il a atteint la fin de la liste d’instant, et pas le chargement, il suffit de «prolonger» la liste d’instant initiale au moment de la POURSUITE du calcul.
Conseils pour la mise en œuvre du pilotage PRED_ELAS.#
Pour les lois ENDO_FRAGILE et ENDO_ISOT_BETON, le coefficient COEF_MULT est relié au pas de temps et à la variation d’endommagement minimale recherchée :
\(\text{COEF\_MULT}=\frac{\Delta t}{\Delta D}\)
Ainsi, si vous souhaitez qu’il y ait au moins un point qui voit son endommagement croître de 20%, et que vos pas de temps valent \(\Delta t=1\) , vous devez définir \(\text{COEF\_MULT}=\frac{1.}{0.2}=5\) . Ainsi plus COEF_MULT est petit, plus l’endommagement progresse rapidement.
Pour les autres lois, le coefficient COEF_MULT assure que, au moins un point de Gauss sorte du seuil d’élasticité linéarisé \({f}_{\mathrm{pred}\text{\_}\mathrm{elas}}\) d’une quantité \(\frac{\Delta t}{\text{COEF\_MULT}}\) ; on a donc \(\text{COEF\_MULT}=\frac{\Delta t}{\underset{\text{pts de Gauss}}{\mathrm{MAX}}({f}_{\text{pred\_elas}})}\) . Dans ce cas, la valeur du paramètre à fournir n’est pas évidente et dépend des lois. On conseille de mettre 1 pour commencer. En fonction de la progression de l’endommagement observé, on pourra accélérer le phénomène en diminuant COEF_MULT ou le ralentir en augmentant COEF_MULT. On peut aussi préférer de manière équivalente modifier la liste des pas temps.
On conseille en général d’opter pour SELECTION=”RESIDU”,
Lorsqu’on pilote le chargement, il ne faut pas hésiter à autoriser beaucoup de redécoupage (en utilisant la commande DEFI_LIST_INST) car certains passages sont délicats.
Conseils pour la mise en œuvre du pilotage DEFORMATION.#
Pour le pilotage en déformation, il est obligatoire de faire au moins un incrément de charge sans le pilotage. On conseille en général de ne basculer avec le pilotage que lorsqu’on ne converge plus avec les autres méthodes .
L’idée ici est d’assurer qu’au moins un point de la structure voit sa déformation progresser de façon monotone :
\(\text{COEF\_MULT}\times \underset{\text{pts de Gauss}}{\mathit{MAX}}(\frac{{\varepsilon}^{-}}{\mathrm{\parallel }{\varepsilon}^{-}\mathrm{\parallel }}.\Delta \varepsilon )=\Delta t\) .
La valeur de COEF_MULT n’est pas toujours facile à calibrer, mais c’est en général des valeurs importantes puisqu’inversement proportionnelle à l’incrément de déformation (exemple pour des pas de \(0,1s\) et un accroissement de déformation de \({10}^{\text{-6}}\) , \(\text{COEF\_MULT}=\frac{0,1}{{10}^{\text{-6}}}=100000\) ). Après un premier test, on ajustera si besoin la valeur de ce coefficient en augmentant la valeur pour charger moins vite, et en la diminuant pour endommager plus vite. On peut aussi préférer agir sur la liste des pas de temps.
La dynamique#
En dernier recours, on peut également essayer de se lancer dans un calcul dynamique. Dans certains cas, cela peut apporter des solutions mais c’est à utiliser avec beaucoup de précautions. Une documentation U2 a spécialement été rédigée pour rassembler les conseils, il s’agit de la doc [U2.04.07], qu’il est indispensable de consulter avant de se lancer dans une simulation en dynamique (pour résoudre un problème quasi-statique).
Remarque :
Depuis la version 10 de Code_Aster, deux méthodes seront également disponibles pour l’utilisateur : la gestion automatique du pas de temps (qui peut se révéler une aide intéressante) et la méthode IMPLEX pour les lois de comportement ENDO_FRAGILE et ENDO_ISOT_BETON mais uniquement en version locale.