u2.05.04 Notice d’utilisation pour le calcul de charge limite

Résumé

L’objectif de cette note est de donner les informations nécessaires pour qu’un utilisateur puisse réaliser des calculs de charge limite avec code_aster .

La première partie rappelle les grandes lignes de la méthode et ses propriétés.

La deuxième partie présente les différentes étapes nécessaires à la mise en œuvre dans code_aster .

Finalement, la troisième partie présente quelques remarques de mise en œuvre en s’appuyant sur le cas-test Aster [V6.04.124] et sur le cas industriel visant à déterminer la pression limite pour un joint de cuve Canopy.

Les grandes lignes de la méthode

Présentation

Les objectifs de l’analyse limite sont:

• l’analyse de sécurité face à un comportement extrême (État Limite Ultime E.L.U.);

• le dimensionnement rapide sans chercher à décrire l’ensemble du processus de ruine;

• la caractérisation énergétique de la ruine et la compréhension des modes de ruine;

• l’obtention d’une information simple sur l’évolution non linéaire du matériau de la structure.

L’analyse limite est un problème que l’on peut traiter de deux manières, cf. [bib8], [bib9], [bib10]:

• calcul de ruine «plastique» de structures élastoplastiques avec plateau ductile. Le trajet de chargement et le modèle de comportement du matériau doivent être décrits entièrement;

• calcul de perte de potentialité d’équilibre à critère de résistance donné, pour une direction de chargement donnée. Il s’agit d’un problème d’optimisation (sous contrainte) du paramètre de charge \(\lambda\) . Cette approche est appelée «calcul à la rupture».

Pour les matériaux standards, ces deux méthodes donnent le même résultat.

Le «calcul à la rupture» ou «analyse limite» (vocable désignant le calcul à la rupture dans le cas d’un matériau élastoplastique à règle d’écoulement normal) vise à déterminer directement, de façon simplifiée et sans avoir recours à la description du trajet de chargement par un calcul incrémental élastoplastique onéreux, la frontière du domaine de ruine plastique (et par déduction le domaine des charges supportables) pour une structure \(\Omega\) , de géométrie et de limites de résistance des matériaux données, soumise à un chargement donné par sa direction \(\text{F}\) , et d’amplitude paramétrée par le réel positif \(\lambda\) . Un chargement «permanent» \({\text{F}}_{0}\) , comme par exemple la pesanteur, peut éventuellement être présent en plus (sans être amplifié par \(\lambda\) ).

Fig. 33 Ingrédients du calcul à la rupture

Remarque

On ne peut pas tenir compte d’un quelconque changement de géométrie par le calcul à la rupture, comme cela arrive lors d’une ruine par flambage, ou pour un solide très souple. Cela constitue les hypothèses de la méthode: la configuration du solide est celle de sa géométrie initiale, les liaisons de la structure sont supposées données et fixes jusqu’à la ruine; de même les chargements, en efforts uniquement, sont de directions fixées.

Deux approches du calcul à la rupture sont accessibles:

• l’approche statique qui estime la valeur de charge limite par l’intérieur et qui nécessite la construction des champs de contraintes statiquement admissibles, ce qui est délicat en général par éléments finis. Elle consiste à maximiser le (ou les) paramètre(s) de chargement à la condition que les équations (linéaires) de la statique restent vérifiées et que le critère en contraintes ne soit pas violé;

• l’approche cinématique, duale de la précédente, qui estime la valeur de charge limite par l’extérieur et qui nécessite une minimisation par une méthode de régularisation du (ou des) paramètre(s) de chargement sous la condition que la puissance des efforts extérieurs reste supérieure à la puissance résistante, (définie à partir du critère de résistance), d’une fonctionnelle non régulière, qui doit en conséquence être régularisée dans le cas général.

L’emploi combiné (cas idéal!) de ces deux approches fournit des encadrements de la charge limite.

Un exemple analytique [bib7]

On considère un système hyperstatique à trois barres (voir Fig. 34), les barres ont un critère de résistance identique, exprimé en terme d’effort normal (ou tension) \(N\) : \(g(N)=\mid N\mid \le \stackrel{ˉ}{N}\) .

Le point \(D\) est soumis à une force \(\overrightarrow{F}\) de composantes \(({F}_{x},{F}_{y})\) non nulles, amplifiée par un facteur multiplicatif \(\lambda\) .

Fig. 34 Système à trois barres

Fig. 35 Domaine des tensions supportables

Fig. 36 Domaine des chargements supportables

L’espace de solutions statiquement admissibles (les tensions \({T}_{i}\) dans les barres) est défini par les équations:

(7)\[\begin{split}\left \lbrace \begin{array}{c}{T}_{1}=\frac{\lambda ({F}_{x}\cos\alpha +{F}_{y}\sin\alpha )-X\sin\alpha }{\sin2\alpha }\\ {T}_{2}=X\\ {T}_{3}=\frac{\lambda ({F}_{y}\sin\alpha -{F}_{x}\cos\alpha )-X\sin\alpha }{\sin2\alpha }\end{array}\right.\end{split}\]

Pour \(\lambda\) positif, on constate que c’est en \({T}_{1}={T}_{2}=\stackrel{ˉ}{N}\) que l’extremum est atteint, voir Fig. 35. On trouve ainsi que la valeur maximale supportable du facteur de charge \(\lambda\) (ou charge limite), obtenue par l’approche statique, est:

(8)\[\stackrel{ˉ}{\lambda}=\stackrel{ˉ}{N}\frac{(1+2\cos\alpha )\sin\alpha }{{F}_{x}\cos\alpha +{F}_{y}\sin\alpha }\]

L’approche cinématique donne le même résultat [bib6]: c’est bien la charge limite de ce problème.

Quelques propriétés utiles du calcul de charge limite

Dans le cas idéal, les bornes supérieure et inférieure de la charge limite doivent être égales à la valeur limite. Avec l’approche numérique, on aura toujours un écart et c’est la borne inférieure qui est la plus pénalisante. Il est cependant à noter que, pour les structures dont on peut conduire le calcul suffisamment loin (comme dans les cas tests), la borne supérieure est en pratique celle qui est la plus proche de la valeur exacte.

On choisit fréquemment comme seuil de résistance la limite d’élasticité: cela va dans le sens de la sécurité.

On rappelle ci-après quelques propriétés utiles du calcul de charge limite (voir [bib4], [bib8]):

• la charge limite est proportionnelle à la valeur de la limite de résistance ou seuil \({\sigma}_{y}\) dans un solide homogène. Elle ne dépend pas de l’histoire du chargement subi par la structure au préalable;

• comme le critère de résistance est convexe (critère de vonMises), le domaine des chargements supportables (donc la frontière des chargements limites) dans l’espace des chargements est convexe. On peut donc approcher le domaine des chargements admissibles par le polyèdre généralisé construit sur les sommets, correspondant chacun à une direction choisie dans l’espace des chargements;

• les conditions de Dirichlet (un déplacement imposé) qui sont appliquées sur la partie \({\Gamma}_{u}\) du bord \(d\Omega\) de la structure, ou une déformation anélastique initiale –thermique, plastique…– n’ont pas d’effet sur le domaine des charges admissibles, (la ruine étant l’impossibilité de satisfaction de l’équilibre, le mode de ruine correspond à une vitesse-direction d’écoulement);

• la charge limite ne dépend pas de la présence éventuelle d’un champ de contraintes auto-équilibrées (contraintes résiduelles);

• pour un solide bidimensionnel, le matériau et la direction de chargement étant donnés, une borne inférieure obtenue par l’approche statique en contraintes planes est nécessairement inférieure à la charge limite exacte obtenue en déformations planes : \({\lambda}_{{C}_{\text{PLAN}}}^{-}\le {\lambda}_{{D}_{\text{PLAN}}}^{\lim}\) ;

• réciproquement, une borne supérieure obtenue par l’approche cinématique en déformations planes est nécessairement inférieure à la charge limite exacte obtenue en contraintes planes : \({\lambda}_{{D}_{\text{PLAN}}}^{+}\ge {\lambda}_{{C}_{\text{PLAN}}}^{\lim}\) . Ce résultat fournit donc un majorant. Si l’on souhaite traiter un problème en contraintes planes, il est nécessaire alors de faire l’approche cinématique sur une modélisation tridimensionnelle d’une «tranche» de solide;

• les charges limites obtenues en déformations planes avec le critère de Tresca valent \(\sqrt{3}/2\) fois celles trouvées avec le critère de von Mises;

• à géométrie et direction de chargement données, si on remplace dans une zone donnée de la structure le matériau de domaine de résistance \({G}_{1}\) par un matériau de domaine de résistance \({G}_{2}\subset {G}_{1}\) (\({g}_{2}(\tau )\le {g}_{1}(\tau )\) ) (par exemple: critère de Tresca inclus dans celui de von Mises), alors les fonctions d’appui (puissances résistantes maximales) sont: \({\pi}_{2}(\varepsilon (v))\le {\pi}_{1}(\varepsilon (v))\) et donc: \({\lambda}_{2}^{\lim}\le {\lambda}_{1}^{\lim}\) ;

• en particulier si on remplace un défaut \(1\) (trou, fissure) présent dans la structure par le défaut \(2\) contenant le défaut \(1\) , alors on a: \({\lambda}_{2}^{\lim}\le {\lambda}_{1}^{\lim}\) . De même, si la structure est hétérogène, avec deux zones dont les limites de résistance sont \({\sigma}_{\mathrm{y1}}\le {\sigma}_{\mathrm{y2}}\) , la charge limite sera supérieure à celle de la même situation homogène pour le seuil \({\sigma}_{\mathrm{y1}}\) et inférieure à celle pour le seuil \({\sigma}_{\mathrm{y2}}\) ;

• en présence d’une direction de chargement \(f=\lambda (\alpha {f}_{1}+(1-\alpha ){f}_{2})\) combinant deux directions \({f}_{1}\) et \({f}_{2}\) , \(\alpha \in \left[0,1\right]\) , alors la charge limite exacte vérifie: \({\lambda}^{\lim}(\alpha )\ge \stackrel{ˉ}{\lambda}=\frac{{\lambda}_{1}^{\lim}{\lambda}_{2}^{\lim}}{(1-\alpha ){\lambda}_{1}^{\lim}+\alpha {\lambda}_{2}^{\lim}}\) . Ce résultat reste valable pour des approximations par l’intérieur des charges limites.

Alors que les situations tridimensionnelles générales sont inabordables analytiquement, en 2D déformations planes (D_PLAN) et contraintes planes (C_PLAN) il est possible de construire à la main et de calculer des solutions pour l’approche statique et l’approche cinématique, à l’aide de champs construits par blocs, qui donnent des encadrements de la charge limite: ceci se révèle utile pour corroborer un résultat obtenu par éléments finis.

La mise en œuvre de l’analyse limite

Les étapes de calcul

Pour réaliser une analyse limite, il faut:

• utiliser un maillage 2D (plan ou axis) ou 3D compatible avec les éléments finis incompressibles;

• définir le modèle avec les éléments finis incompressibles;

• définir la limite de résistance du matériau \({\sigma}_{y}\) ;

• définir le chargement permanent \({\text{F}}_{0}\) et le chargement variable \(\text{F}\) qui est piloté par \(\lambda\) ;

• définir la condition d’incompressibilité;

• réaliser un calcul non-linéaire avec le comportement incrémental et le pilotage prévus pour l’analyse limite;

• post-traiter le calcul pour obtenir les valeurs supérieure et inférieure de la charge limite estimée.

Maillage

Pour le choix du maillage et l’utilisation des éléments incompressibles, on renvoie à la documentation [r3.06.08]

Remarque: l’usage d’élément à trois champs UPG nécessitera d” imposer une condition aux limites sur le gonflement (GONF), les nœuds sommets doivent donc être parfaitement identifiés dans le maillage:

• soit dès la phase de conception et de réalisation du maillage;

• soit, plus simplement et a posteriori, en faisant appel dans à la commande DEFI_GROUP [u4.22.01] avec le mot-clé CREA_GROUP_NOen précisant le groupe de mailles d’éléments incompressibles et en utilisant l’option CRIT_NOEUD = “SOMMET”.

Modèle

Options de modélisation

Il faut sélectionner des éléments incompressibles, on renvoie donc à la commande AFFE_MODELE [U4.41.01].

Condition d’incompressibilité

Pour exprimer la condition d’incompressibilité pour les éléments incompressibles UPG à trois chams, on utilise la commande AFFE_CHAR_MECA [u4.44.01] avec le mot-clé DDL_IMPO pour imposer à la composante GONF du groupe des nœuds sommets des éléments incompressibles de rester nulle.

Matériau

Le matériau utilisé pour l’analyse limite avec les éléments incompressibles est un matériau avec un critère de von Mises, élastoplastique parfaitement plastique.

Les données des caractéristiques de matériaux sont fournies sous le mot clé facteur ECRO_LINE de la commande DEFI_MATERIAU [u4.43.01].

La pente de la courbe de traction (\({E}_{T}\) ) est choisie nulle (opérande D_SIGM_EPSI), la seule donnée nécessaire à fournir est donc la limite d’élasticité (c’est-à-dire le seuil de résistance dans notre cas) (opérande SY).

Chargement

Le chargement variable \(\text{F}\) qui est piloté par \(\lambda\) doit nécessairement être de type effort (force, pression, pesanteur) et déclaré dans la commande AFFE_CHAR_MECA[U4.44.01].

Si la structure est aussi soumise à un chargement permanent \({\text{F}}_{0}\) , il faut penser à le rappeler lors du post-traitement (voir u2.05.04-post-traitement).

Liste d’instants

La liste d’instants sert à contrôler la méthode de régularisation de Norton-Hoff, cf. [bib1], [bib5], par l’intermédiaire d’un coefficient \(m\) , et non pas l’évolution du chargement comme lors d’un calcul ordinaire:

(9)\[m=1+{10}^{1-t}\]

de sorte que, quand l’instant devient suffisamment grand, \(m\) tende vers 1, et le comportement se rapproche d’un comportement rigide-plastique parfait, voir la courbe uniaxiale Fig. 37.

Fig. 37 Courbe contrainte – déformation pour différentes valeurs de l’instant \(t\)

En pratique, on choisira au début une liste d’instants à pas constants (voir Tableau 2.6-a) avant de raffiner les pas de calcul en cas de non convergence.

\(t\)	1.0	1.5	2.0	2.5	3.0	…
\(m=1+{10}^{1-t}\)	2.00	1.30	1.10	1.03	1.01

Si le document d’utilisation de la commande POST_ELEM avec le mot-clé CHAR_LIMITE [U4.81.22] recommande, en pratique, de se limiter à des instants compris entre 1 et 2 pour ne pas avoir des calculs trop longs tout en permettant d’obtenir une borne supérieure de la charge limite suffisamment précise, nous observons que la charge limite inférieure nécessite au moins 2 à 3 itérations supplémentaires (instant supérieur à 3) pour converger vers les valeurs de référence lors des cas tests.

Calcul

La modélisation avec des éléments incompressibles doit nécessairement utiliser la commande STAT_NON_LINE [u4.51.03] et le mot-clé COMPORTEMENT.

Pour l’analyse limite, il est obligatoire de faire appel aux opérandes suivants de la commande STAT_NON_LINE:

• Dans COMPORTEMENT, on choisit RELATION=”NORTON_HOFF”. Cet opérande est utilisé pour décrire la relation de comportement de viscosité (indépendante de la température) dans le calcul de charges limites de structures, à seuil de von Mises [u4.51.11]

• Dans PILOTAGE, on choisit TYPE=”ANA_LIM”.Ce mode de pilotage est spécifique au calcul de charge limite (loi NORTON_HOFF) par approche cinématique. Il doit être seulement appliqué à la charge déclarée via le mot-clé facteur EXCIT, opérande TYPE_CHARGE=”FIXE_PILO”.

Le calcul de la charge limite peut requérir beaucoup d’itérations de recherche linéaire et d’itérations de Newton. Il est donc aussi fortement conseillé d’employer les différentes options de calcul de la commande STAT_NON_LINE pour améliorer la convergence, comme la recherche linéaire dont la pratique montre qu’il suffit d’avoir recours à 2 ou 3 itérations.

Remarque

si on amplifie l’intensité du chargement \(L\to \mathrm{\beta }L\) (alors que l’on ne considère pas de chargement permanent \({L}_{0}=0\) ), les solutions dépendent du facteur \(\mathrm{\beta }\) selon les relations suivantes:

(10)\[{u}_{m}(\beta )={\beta}^{-1}{u}_{m}(1);{\sigma}^{D}({u}_{m}(\beta ))={\beta}^{1-m}{\sigma}^{D}({u}_{m}(1))\]

À la convergence pour \(m\to {1}^{+}\) , le chargement limite donné par la solution \({u}_{m}(\mathrm{\beta })\) est bien le même que celui donné par \({u}_{m}(1)\) , puisque \({\lambda}_{\lim}(\beta )={\lambda}_{\lim}(1)/\beta\) .

Post-traitement

A partir du résultat du calcul non linéaire réalisé, l’opérateur POST_ELEM [u4.81.22] et le mot-clé CHAR_LIMITE produisent une table qui donne, pour chaque instant du calcul, l’estimation de la borne supérieure CHAR_LIMI_SUP (\({\stackrel{ˆ}{\lambda}}_{m}\) ) de la charge limite supportée par la structure; cette suite est monotone décroissante quand \(m\to {1}^{+}\) , c’est-à-dire quand \(t\to +\infty\) .

En outre, en l’absence de chargement permanent \({F}_{0}\) (opérande CHAR_CSTE=”NON” qui est l’option par défaut), la table contient également l’estimation CHAR_LIMI_ESTIM (\({\underline{\lambda}}_{m}\) ) de la borne inférieure de la charge limite. Cette valeur (approximation de la jauge du convexe de résistance) n’est calculée qu’aux points de Gauss des éléments finis. Aussi la valeur \({\underline{\lambda}}_{m}\) obtenue pour chaque \(m\) , inférieure à \({\stackrel{ˆ}{\lambda}}_{m}\) [bib6], ne peut être considérée que comme une indication (cette suite n’est pas nécessairement monotone). Elle permet, avec la valeur par excès \({\stackrel{ˆ}{\lambda}}_{m}\) , de fournir un encadrement de la charge limite du problème discrétisé.

En revanche, si un chargement permanent \({F}_{0}\) est présent (opérande CHAR_CSTE=”OUI”), une telle estimation de la borne inférieure n’est plus disponible et la table indique alors la puissance PUIS_CHAR_CSTE du chargement constant dans le champ de vitesse solution du problème.

La visualisation du champ de déplacement obtenu pour une valeur du coefficient \(m\to {1}^{+}\) donne une «idée» du mode de ruine de la structure étudiée. On trouve plus de détails sur le post-traitement dans [r7.07.01].

Quelques remarques de mise en œuvre

Cas test représentatif : ssnv146a

Description

Ces remarques sont basées sur la mise en œuvre du cas test ssnv146a qui est tiré du benchmark du projet européen Brite EuRam BE97-4547 « LISA » [bib8]. Il s’agit du calcul de la charge limite d’un réservoir avec un fond tori-sphérique en 2D axisymétrique, sous pression interne, voir Figure 3.a. Le rayon interne de la partie cylindrique est: \(\mathrm{49 }\mathit{mm}\) , tandis que l’épaisseur est: \(\mathrm{2 }\mathit{mm}\) . Le rayon de la partie sphérique à l’apex est \(\mathrm{98 }\mathit{mm}\) , tandis que le rayon du tore de raccordement est de \(\mathrm{20 }\mathit{mm}\) .

Tableau 13 Données du test SSNV146 (benchmark LISA) [bib8]

Géométrie	Le réservoir axisymétrique à fond torisphérique (voir Figure 3.1-a) a les caractéristiques suivantes: * rayon interne de la partie cylindrique: \(49\mathit{mm}\) ; * épaisseur: \(2\mathit{mm}\) ; * rayon de la partie sphérique à l’apex: \(98\mathit{mm}\) ; * rayon du tore de raccordement: \(20\mathit{mm}\) *
Propriétés de matériau	Limite d’élasticité : \({\sigma}_{y}=100\mathit{MPa}\)
Conditions aux limites	Conditions de symétrie
Chargements	Pression interne de \(1\mathit{MPa}\)

Le tableau suivant récapitule les résultats obtenus par les participants au benchmark utilisant le même maillage qui contient 34 éléments QUAD8 (dont deux éléments dans l’épaisseur) et 141 nœuds.

Tableau 14 Résultats du benchmark LISA [bib8]

Modélisation				Valeur supérieureestimée	Valeur inférieureestimée
	\(m=1,0476\) | \(n=21\) | \(t=2,3222\)			3.9514	3,6049
EDF (1)	\(m=1,0322\) | \(n=31\) | \(t=2,4914\)			3.9456	3,7090
	\(m=1,0141\) | \(n=71\) | \(t=2,8513\)			3.9404	3,8372
	\(m=1,0099\) | \(n=101\) (1) | \(t=3,0043\) (2)			3.9396	3.8673
Univ. de Liège/LTAS				3,931	néant
Centre de recherche FZJ				néant	3.997

Nota (1): le coefficient de régularisation par la loi de Norton-Hoff \(n={(m-1)}^{-1}\) (ancienne méthode d’EDF)

Nota (2): \(m=1+{10}^{1-t}\) , le paramètre ‘instant’ \(t\) n’était pas explicite dans cette ancienne méthode, on ajoute ici pour pouvoir mieux comparer avec les résultats dans les sections 4,2 et 4,3.

La figure Fig. 39 présente l’évolution des termes des deux suites \({\underline{\lambda}}_{m}\) et \({\stackrel{ˆ}{\lambda}}_{m}\) en fonction du coefficient de régularisation \(n={(m-1)}^{-1}\) avec \(m=1+{10}^{1-t}\) , cf. [bib7]. Les bornes supérieure et inférieure de la charge limite estimée sont calculées avec la liste d’instants, jouant directement sur le coefficient de régularisation par la loi de Norton-Hoff. Ceci permet de réaliser directement cette convergence et simplifier l’utilisation.

Fig. 38 Maillage initial et déformé pour \(m=1,0322\)

Fig. 39 convergence des suites \({\underline{\lambda}}_{m}\) et \({\stackrel{ˆ}{\lambda}}_{m}\) en fonction du coefficient de régularisation, vers la pression limite exacte (benchmark LISA). |

Ce sont les résultats du calcul EDF effectué dans le cadre du benchmark LISA avec l’ancienne méthode.

Résultats du cas test SSNV146

Le maillage est celui utilisé pour le benchmark «LISA». Le calcul a été mené jusqu’à l’instant \(t=2s\) . La version utilisée est la version d’exploitation STA 8.3.

—- TABLE: ECHL1 NOM_PARA: CHAR_LIMI_SUP REFERENCE: NON_DEFINI OK ECHL1 RELA 0.690 % VALE: 3.9581130295563D+00 CHAR_LIMI_SUP TOLE 1.000 % REFE: 3.9310000000000D+00

La valeur de référence déclarée (3.931) est une valeur estimée pour la borne supérieure et elle a été fournie par l’université de Liège pour le benchmark «LISA», voir Tab.3.1-b.

Calcul plus poussé

Pour améliorer l’estimation des valeurs supérieure et inférieure de la charge limite estimée, le calcul a été poussé plus loin jusqu’à l’instant t =2,45 s (rappelons que ce n’est pas un temps physique, voir u2.05.04-liste-instants) au-delà duquel le calcul ne converge plus.

Tableau 15 Évolution des bornes sup et inf estimées de la charge limite. Maillage initial

INST	M	CHAR_LIMI_SUP \({\stackrel{ˆ}{\lambda}}_{m}\)	CHAR_LIMI_ESTIM \({\underline{\lambda}}_{m}\)
\(t=1\)	\(m=2\)	4,30410	1,27329
\(t=2\)	\(m=1,1\)	3.94863	3, 27091
\(t=2,125\)	\(m=1,0750\)	3.93504	3, 40992
\(t=2,25\)	\(m=1,0562\)	3.92505	3, 52215
\(t=2,450\)	\(m=1,0355\)	3.91429	3.65061

Quand on compare les résultats obtenus (Tableau 15) avec celles fournies par EDF l’ancienne version (Tableau 14), on observe que l’équivalence « instant de calcul» et « coefficient Norton-Hoff» ne semble pas être identique si l’on se base sur la charge limite supérieure ou la charge limite inférieure.

En effet, si la tendance décroissante (resp. croissante) de la charge limite supérieure (resp. inférieure) est confirmée (Fig. 39), alors, à \(t=2,45s\) , le coefficient Norton-Hoff équivalent serait supérieur à \(n=101\) pour la valeur supérieure alors qu’il serait inférieur à \(n=21\) pour la valeur inférieure.

Influence de la finesse des maillages

Le calcul du cas test avec le maillage initial ne converge plus après \(t=2.45s\) . Comme pour tout autre calcul en non-linéaire, le maillage a été raffiné pour tenter d’améliorer la convergence.

Ce calcul permettrait surtout de mieux approcher la charge limite par une meilleure estimation de la valeur inférieure.

Avec une discrétisation deux fois plus fine, soit un maillage de 136 éléments QUAD8 (soit 4 éléments dans l’épaisseur), le calcul a été amené à convergence jusqu’à \(t=2,153s\) .

Tableau 16 Evolution des bornes sup et inf estimées de la charge limite. Maillage fin

INST ( s )		CHAR_LIMI_SUP \({\widehat{\lambda}}_{m}\)	CHAR_LIMI_ESTIM \({\underline{\lambda}}_{m}\)
1.00000	\(m=2,000\)	4.32038	1.18088
2.00000	\(m=1,100\)	3.97963	3.26240
2.0625	\(m=1,074\)	3.97297	3.34124
2.125	\(m=1,07499\)	3.96726	3.41252
2.1531	\(m=1,07029\)	3.96498	3.44220

L’écart relatif \(\dfrac{\lambda_{\lim}^{\sup} - \lambda_{\lim}^{\inf}}{\dfrac{1}{2} \cdot (\lambda_{\lim}^{\sup} + \lambda_{\lim}^{\inf})}\) passe de 6,97% (maillage initial) à 14,11% avec le maillage plus fin.

On note aussi qu’aux mêmes instants (\(t=\mathrm{1 }s\) et \(\mathrm{2 }s\) ), on n’a plus exactement les mêmes valeurs supérieure et inférieure de la charge limite estimée avec le maillage initial (Tableau 15).

L’influence de la finesse du maillage se manifeste donc à la fois dans la convergence globale et dans la précision des valeurs estimées aux instants sauvegardés. On en conclut que le calcul de charge limite par les éléments finis est sensible aux maillages. Le maillage plus fin ne donne pas forcément une limite plus poussée à cause d’un arrêt de calcul plus tôt. Ceci n’est néanmoins pas un problème, sachant que l’analyse de charge limite par la loi Norton-Hoff tend vers un problème singulier qui poserai des difficultés de convergence. L’objectif de l’analyse est de trouver une suite de valeurs de \({\widehat{\lambda}}_{m}\) convergeant vers une charge limite \({\lambda}_{\lim}\) .

Evolution de l’estimation avec la liste d’instants

Pour le cas test ssnv146a considéré, la valeur de référence n’est pas une valeur théorique mais issue elle-même d’un calcul numérique.

Le cas test ssnv124, décrivant une situation d’un solide 2D ou 3D chargé de manière homogène, permet de comparer les valeurs estimées par rapport à une solution analytique, voir [v6.04.124] (pour laquelle \({\stackrel{ˆ}{\lambda}}_{m}={\lambda}_{\lim},\forall m\) ).

Tableau 17 Cas test SSNV124 évolution de l’erreur relative pour les valeurs sup et inf estimées pour la charge limite

	erreur relative par rapport à la valeur de référence (en %)
INST (s)	Valeur sup	Valeur inf
1.00	0.00	50.00
2.00	0.99	0.99
3.00	0.00	0.99
4.00	0.00	0.10
5.00	0.00	0.01

Le Tableau 17 montre la convergence des valeurs supérieure et inférieure estimées vers la valeur analytique de référence. Cette convergence est lente pour la valeur inférieure estimée et elle est beaucoup plus rapide pour la valeur supérieure.

Ces résultats tendent à montrer que si l’on considère la moyenne des valeurs supérieure et inférieure on obtiendrait une valeur conservative de la charge limite estimée.

Problèmes présentant des symétries

Notons au passage une facilité importante dans la mise en œuvre. Bien que la charge limite \({\lambda}_{\lim}\) soit calculée sous forme d’une intégrale sur le domaine, il n’est pas nécessaire de multiplier la valeur obtenue si l’on fait le calcul sur une sous-partie du solide, cellule de symétrie du problème.

Nous l’illustrons avec une variante 3D de la modélisation C axisymétrique du cas test ssnv124.

Le cas représente ainsi un cylindre de rayon interne \(a=1\mathit{mm}\) et externe \(b=3\mathit{mm}\) soumis à une pression interne de \(1\mathit{MPa}\) en paroi interne.

Le résultat théorique \({\lambda}_{\lim}=\frac{2\sqrt{3}}{3}.{\sigma}_{y}.\ln\frac{b}{a}\) donne une charge limite de \(8.00377\mathit{MPa}\) .

Fig. 40 Cylindre sous pression interne (test SSNV124c)

Avec le maillage du ¼ de cylindre présenté Fig. 41, on peut lire dans le tableau ci-dessous directement les résultats des valeurs supérieure et inférieure de la charge limite estimée.

Tableau 18 Evolution des bornes sup et inf estimées de la charge limite. Cas test SSNV124c

#TABLE_SDASTER NUME_ORDRE INST CHAR_LIMI_SUP CHAR_LIMI_ESTIM 1 1.00000E+000 8.00361E+000 2.34638E+000 2 1.69897E+000 8.00360E+000 5.91927E+000 3 2.00000E+000 8.00360E+000 6.84900E+000 4 3.00000E+000 8.00360E+000 7.87601E+000 5 4.00000E+000 8.00360E+000 7.99071E+000 6 5.00000E+000 8.00360E+000 8.00231E+000 7 6.00000E+000 8.00360E+000 8.00347E+000 8 7.00000E+000 8.00360E+000 8.00359E+000 9 8.00000E+000 8.00360E+000 8.00360E+000

Fig. 41 Maillage du quart de cylindre considéré (SSNV124)

Comparaison analyse limite et calcul élastoplastique incrémental jusqu’à la ruine sur un exemple

On sait qu’en élastoplasticité parfaite, l’opérateur tangent possède des valeurs propres nulles à partir d’un certain niveau de chargement en effort: cela signifie que le solide élastoplastique a atteint la ruine plastique.

Il peut être intéressant de faire un calcul élastoplastique (avec critère de von Mises) incrémental jusqu’à la ruine, afin d’obtenir une «borne» inférieure de la charge limite. Comme l’algorithme de Newton utilisé pour résoudre l’équilibre statique non linéaire diverge pour ce niveau de chargement, il faut utiliser un pilotage par longueur d’arc [r5.03.80], ou sur une variable de déplacement, servant à contrôler le chargement par la déformation du solide subie avant la ruine.

Voici un exemple industriel: il s’agit du calcul de la pression interne limite d’un Joint de cuve Canopy. Le problème est axisymétrique; la pièce est bloquée sur sa frontière supérieure. La limite de résistance (forfaitaire) est fixée à \(100\mathit{MPa}\) ; on admet un critère de vonMises. On a réalisé deux calculs:

• un calcul d’analyse limite avec la méthode présentée dans les paragraphes précédents;

• un calcul incrémental élastoplastique parfait (sans et avec écrouissage), avec critère de vonMises en grandes transformations (Simo-Miehe, voir [r5.03.21]), afin de vérifier que le changement de géométrie ne modifie pas substantiellement la prédiction de pression limite.

Fig. 42 Maillage déformé (amplifié)

La Fig. 42 montre que les déformées calculées sont très voisines entre les deux méthodes et permettent de prédire le mode de ruine.

La Fig. 43 montre la convergence des suites des bornes de pression limite \({\underline{\lambda}}_{m}\) et \({\stackrel{ˆ}{\lambda}}_{m}\) en fonction du coefficient de régularisation de la méthode de calcul d’analyse limite, vers la pression limite exacte. La moyenne arithmétique des deux bornes semble constituer une bonne estimation de la pression limite. Pour la dernière valeur du coefficient de régularisation choisie, on a obtenu les valeurs reportées au Tableau 19. On les compare avec le calcul incrémental élastoplastique parfait en grandes transformations. On constate que les valeurs sont très proches.

Tableau 19 Pressions limites

analyse limite: borne inférieure	analyse limite: borne supérieure	calcul incrémental élastoplastique parfait
\(18,72\mathit{MPa}\)	\(23,84\mathit{MPa}\)	\(23,25\mathit{MPa}\)

Fig. 43 Calcul d’analyse limite : convergence des suites \({\underline{\lambda}}_{m}\) et \({\stackrel{ˆ}{\lambda}}_{m}\) en fonction du coefficient de régularisation, vers la pression limite exacte.

Le calcul incrémental élastoplastique parfait en grandes transformations a été réalisé avec deux types de pilotage:

• par un déplacement radial d’un point particulier;

• par longueur d’arc (ce dernier pilotage ayant été plus performant en temps calcul);

et on l’a comparé les résultats avec ceux obtenus par un calcul en petites transformations.

Une des questions qui se pose est enfin l’effet de l’écrouissage du matériau, et donc du choix du seuil de résistance.

C’est pourquoi on a réalisé un calcul élastoplastique incrémental en tenant compte de la courbe de traction du matériau considéré (les paramètres principaux étant: \({\sigma}_{y}=195\mathit{MPa}\) , \({\sigma}_{u}=520\mathit{MPa}\) ).

La Fig. 44 permet de constater que sur ce cas les deux pilotages donnent une solution identique en grandes transformations, qui est plus «souple» qu’en petites transformations. La pression de ruine trouvée est supérieure à \(56,4\mathit{MPa}\) . Rapportée pour une valeur \({\sigma}_{y}=100\mathit{MPa}\) , on aurait «trouvé»: \(28,9\mathit{MPa}\) .

Cet exemple permet donc de jauger l’effet conservatif de la méthode d’analyse limite en ayant pris comme seuil de résistance la limite d’élasticité \({\sigma}_{y}\) .Par contre, prendre la limite ultime \({\sigma}_{u}\) comme seuil semble non conservatif, la structure subissant des changements de géométrie substantiels dès qu’elle plastifie, avant la ruine, dans le trajet considéré pour le calcul incrémental.

Fig. 44 Calcul incrémental élastoplastique en grandes transformations (Simo-Miehe) avec deux types de pilotage, comparé avec un calcul en petites transformations.

Bibliographie

[bib1]

[R7.07.01] Calcul de charge limite par la méthode de Norton-Hoff-Friaâ , 2007.

[bib2] [R3.06.08] Éléments finis traitant la quasi-incompressibilité , 2005.

[bib3] [R5.03.80] Méthodes de pilotage du chargement , 2001.

[bib4]

[V6.04.124] SSNV124 – Analyse limite régularisée. Loi de Norton-Hoff , 2007.

[bib5]

VOLDOIRE F. – Analyse limite des structures fissurées et critères de résistance , note EDF/DER HI-74/95/026, 1995.

[bib6] (1,2)

VOLDOIRE F. – Limit analysis by the Norton-Hoff-Friaâ regularising method . In M. Heitzer, M. Staat, LISA project report 2001, publication du John von Neumann Institute for Computing (2003).

[bib7] (1,2)

VOLDOIRE F. – Calcul à la rupture et analyse limite des structures , note EDF HI-74/93/082.

[bib8] (1,2,3,4,5)

LAHOUSSE A., VOLDOIRE F. – Calcul de charge limite et benchmark du projet européen Brite EuRam «LISA» , Note EDF/DER HI-74/98/026/A.

[bib9]

SALENçON J. – Cours des structures anélastiques, calcul à la rupture et analyse limite – Presses ENPC, 1983.

[bib10]

SALENçON J. – De l’élastoplasticité au calcul à la rupture – École Polytechnique, 2002.

[bib11] SAVE M.A., MASSONNET C.E., De SAXCéG. – Plastic limit analysis of plates, shells and disks – North Holland Series in Applied Math. & Mech., Elsevier, 1997.

[bib12] [R5.03.21] Modélisation élasto(visco)plastique avec écrouissage isotrope en grandes déformations , 2005.