r7.07.01 Calcul de charge limite par la méthode de Norton-Hoff-Friaâ, comportement NORTON_HOFF

Résumé :

L’analyse limite permet de déterminer les chargements admissibles d‘une structure, de géométrie fixe donnée, constituée d’un matériau ayant un critère de résistance. On considère le cas de chargements constitués de la somme d’une charge permanente et d’une autre paramétrée par le facteur de charge, dont on cherche la valeur extrême supportable.

Après un rappel de la formulation théorique, on présente l’approche cinématique régularisée appliquée au critère de résistance de VonMises (méthode de Norton-Hoff-Friaâ) et mise en œuvre dans Code_Aster . On pourra se reporter à [bib4] pour les diverses méthodes de régularisation possibles proposées dans la littérature. On expose ensuite le calcul des solutions de ce problème non linéaire et le post-traitement fournissant une estimation de la charge limite (une valeur par excès dans tous les cas, et quand il n’y a pas de chargement permanent, une valeur par défaut).

Table des Matières

Aspects numériques du calcul de la charge limite

Relation de comportement de Norton-Hoff

Le tenseur des contraintes \(\sigma (u)\) vérifie la relation de comportement de Norton-Hoff. Le déviateur des contraintes associé à la vitesse de déformation est:

\({\sigma}^{D}\left(\mathrm{u}\right)=A\left(m\right).{\left(\sqrt{{\epsilon}^{D}\left(\mathrm{u}\right).{\epsilon}^{D}\left(\mathrm{u}\right)}\right)}^{m-2}.{\epsilon}^{D}\left(\mathrm{u}\right)\) éq 2.1-1

avec \(\text{tr}\varepsilon (u)=0\) , et: \(A(m)={k}^{1-m}{(\frac{2}{3})}^{m/2}{\sigma}_{y}^{m}\) .

Ce comportement est intégré de la même façon que les comportements incrémentaux élastoplastiques de Von Mises [R5.03.02]. Remarquons toutefois que, en un point d’intégration, le calcul du tenseur des contraintes en fonction du tenseur des déformations est explicite, aucun schéma itératif n’étant utilisé. De plus, aucune variable interne n’est nécessaire à l’intégration de ce comportement.

Dans Code_Aster , le calcul de la charge limite étant indépendant des coefficients d’élasticité , on choisit \(k={\sigma}_{y}\) , d’où \(A(m)={\sigma}_{y}{(\frac{2}{3})}^{m/2}\) . On retrouve ainsi le problème élastique incompressible quand \(m=2\) , pour un module d’Young \(E={\sigma}_{y}\) .

De plus, la suite des scalaires \(m\) est directement déduite de la liste d’instants (fictifs) de calcul par:

\(m=1+{10}^{1-t}\) ,

de sorte que quand l’instant augmente, \(m\) tende vers 1, et le comportement se rapproche d’un comportement rigide-plastique parfait, voir en unidimensionnel les courbes [fig. 2.1-a]. En pratique, on choisit la suite des valeurs du [tab. 2.1-a].

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

	1	1.5	1.7	2	…
	2	1.3	1.2	1.1	…	1

Tab. 2.1-a. Suite des valeurs de l’instant \(t\) , et valeurs de \(m\) correspondantes.

L’opérateur tangent, utilisé dans l’option FULL_MECA de la méthode de Newton, s’écrit grâce à [éq 2.1-1]:

\({\begin{array}{c}\frac{d{\sigma}^{D}}{d{\varepsilon}^{D}}\end{array}\mid }_{\varepsilon (u)}=A(m){\parallel {\varepsilon}^{D}\parallel }^{m-2}(\text{Id}\otimes \text{Id}+(m-2){\parallel {\varepsilon}^{D}\parallel }^{-2}{\varepsilon}^{D}\otimes {\varepsilon}^{D})\) éq 2.1-2

avec \({\varepsilon}^{D}\) , \({\sigma}^{D}\) les vecteurs des déformations et contraintes déviatoriques écrits en notations vectorielles de WALPOLE-COWIN:

\({\varepsilon}^{D}=({\varepsilon}_{11}^{D},{\varepsilon}_{22}^{D},{\varepsilon}_{33}^{D},\sqrt{2}{\varepsilon}_{12}^{D},\sqrt{2}{\varepsilon}_{23}^{D},\sqrt{2}{\varepsilon}_{31}^{D})\)

\({\sigma}^{D}=({\sigma}_{11}^{D},{\sigma}_{22}^{D},{\sigma}_{33}^{D},\sqrt{2}{\sigma}_{12}^{D},\sqrt{2}{\sigma}_{23}^{D},\sqrt{2}{\sigma}_{31}^{D})\)

Pilotage

Le problème s’écrit sous forme variationnelle de la façon suivante sur l’espace des champs incompressibles.

Pour \(m\) donné, donc à un instant \(t\) donné, connaissant la solution à l’instant précédent (notée \({u}^{-},{\lambda}^{-}\) ), trouver \((\Delta \lambda ,\Delta u)\in \mathrm{IR}\times \tilde{{V}_{a}}\) tels que:

\(\lbrace \begin{array}{}{\int}_{\Omega}\sigma ({u}^{-}+\Delta u).\varepsilon (v)d\Omega ={L}_{0}(v)+({\lambda}^{-}+\Delta \lambda )L(v)\forall v\in {\tilde{}V}_{0}\\ L(v)={\int}_{\Omega}f.vd\Omega +{\int}_{{\Gamma}_{f}}F.v\text{dS}=1\end{array}\) éq 2.2-1

• \({L}_{0}\) est le chargement permanent et \(L\) le chargement piloté par le paramètre \(\lambda\) , cf. [§1.1],

• \(\tilde{{V}_{0}}\) est un espace de fonctions discrétisées sur la base d’éléments finis incompressibles, et est donc défini par un vecteur \((U)\) de degrés de liberté.

Ce problème admet une solution unique pour tout \(1\le m\le 2\) (voir [bib4]). Pour \(m=2\) le problème est de type élasticité linéaire incompressible.

Le problème discrétisé à l’instant \(t\) (donc pour une valeur de \(m\) , cf. [tab. 2.1-a]) peut s’écrire (en omettant les conditions aux limites pour simplifier):

\(\lbrace \begin{array}{}{F}_{int}(\Delta U;{U}^{-};..)={F}_{0}^{\text{ext}}+\lambda {F}^{\text{ext}}\\ L({U}^{-}+\Delta U)=1\end{array}\)

La recherche de \(\lambda\) assurant la condition \(L(U)=1\) est assurée par un algorithme de pilotage [R5.03.80].

Brièvement, le principe est le suivant: par linéarisation des équations portant sur les forces intérieures, on obtient, pour l’itération \(n\) de l’algorithme de Newton, cf. [R5.03.01]:

\(\underset{{K}_{T}}{\underset{\underbrace{}}{\left[\frac{\partial {F}_{int}}{\partial U}(\Delta {U}^{n})\right]}}\left[\delta U\right]=\underset{{R}^{\text{cst}}}{\underset{\underbrace{}}{\left[{F}_{0}^{\text{ext}}-{F}_{int}(\Delta {U}^{n})\right]}}+\lambda \underset{{R}^{\text{pilo}}}{\underset{\underbrace{}}{\left[{F}^{\text{ext}}\right]}}\) éq 2.2-2

On peut maintenant exprimer les corrections de déplacements \(\delta U\) et de multiplicateurs de Lagrange \(\delta \lambda\) en fonction de \(\lambda\) moyennant la résolution de ce système linéaire:

\(\left[\delta U\right]=\left[\delta {U}^{\text{cst}}\right]+\lambda \left[\delta {U}^{\text{pilo}}\right]\text{où}\left[\delta {U}^{\text{cst}}\right]={K}_{{T}^{-1}}{R}^{\text{cst}}\text{et}\left[\delta {U}^{\text{pilo}}\right]={K}_{{T}^{-1}}{R}^{\text{pilo}}\) éq 2.2-3

On peut substituer la correction de déplacement \(\delta U\) en fonction de son expression [éq2.2-2] dans l’équation de contrôle du pilotage du système \(L(U)=1\) ; il en résulte une équation scalaire en \(\Delta \lambda\) :

\(L({U}^{-}+\Delta {U}^{n}+\delta {U}^{\text{cst}}+\lambda \delta {U}^{\text{pilo}})=1\) soit

\({\int}_{\Omega}f.({U}^{-}+\Delta {U}^{n}+\delta {U}^{\text{cst}}+\lambda \delta {U}^{\text{pilo}})d\Omega +{\int}_{{\Gamma}_{f}}F.({U}^{-}+\Delta {U}^{n}+\delta {U}^{\text{cst}}+\lambda \delta {U}^{\text{pilo}})\text{dS}=1\) éq 2.2-4

ce qui sous forme discrétisée revient à:

\(\begin{array}{}\\ \sum_{e}{F}^{\text{ext}}.({U}^{-}+\Delta {U}^{n}+\delta {U}^{\text{cst}}+\lambda \delta {U}^{\text{pilo}})=1\end{array}\)

ce qui conduit à:

\(\lambda =\frac{1-\sum_{e}{F}^{\text{ext}}.({U}^{-}+\Delta {U}^{n}+\delta {U}^{\text{cst}})}{\sum_{e}{F}^{\text{ext}}.\delta {U}^{\text{pilo}}}\) éq 2.2-5

Post-traitement du calcul de la charge limite

Ayant obtenu la solution \(({\lambda}_{m},{u}_{m})\) , pour chaque instant, donc chaque \(m\) donné, il reste à utiliser la suite des \({\lambda}_{m}\) pour construire l’approximation de la charge limite. Pour cela on exploite les propriétés [éq 1.3-2], [éq 1.3-3], le fait que \(A(m)\) est croissante et la propriété issue de la minimisation [éq 1.3-2] (voir [bib7]).

De ces deux dernières, avec \(1\le r\le s\) , on déduit que pour \({u}_{r}\) et \({u}_{s}\) solutions respectives (vérifiant aussi la condition d’incompressibilité et de normalisation) de [éq 1.3-2] pour \(m=r\) et

:

\(({\int}_{\Omega}A(r){(\varepsilon ({u}_{r}).\varepsilon ({u}_{r}))}^{r/2}d\Omega )\le ({\int}_{\Omega}A(s){(\varepsilon ({u}_{s}).\varepsilon ({u}_{s}))}^{s/2}d\Omega )\)

Associée à la propriété [éq 1.3-3], on tire pour \(1\le r\le s\) , en notant \({\parallel \Omega \parallel }_{r}={\int}_{\Omega}A(r)d\Omega\) :

\({\int}_{\Omega}A(r)\sqrt{\varepsilon ({u}_{r}).\varepsilon ({u}_{r})}d\Omega \le {\parallel \Omega \parallel }_{r}^{1-\frac{1}{r}}{({\int}_{\Omega}A(r){(\varepsilon ({u}_{r}).\varepsilon ({u}_{r}))}^{r/2}d\Omega )}^{\frac{1}{r}}\le {\parallel \Omega \parallel }_{s}^{1-\frac{1}{s}}{({\int}_{\Omega}A(s){(\varepsilon ({u}_{s}).\varepsilon ({u}_{s}))}^{s/2}d\Omega )}^{\frac{1}{s}}\) éq 2.3-1

En effet, cette propriété résulte de l’inclusion des espaces fonctionnels \({L}^{r}\supset {L}^{s}\) , pour \(1\le r\le s\) , soit aussi, où \(h(x)\) joue le rôle d’une mesure variable (conditionnée par la limite de résistance \({\sigma}_{y}\) ):

\({\parallel \Omega \parallel }_{r}^{-\frac{1}{r}}{({\int}_{\Omega}h(x){\mid f(x)\mid }^{r}d\Omega )}^{\frac{1}{r}}\le {\parallel \Omega \parallel }_{s}^{-\frac{1}{s}}{({\int}_{\Omega}h(x){\mid f(x)\mid }^{s}d\Omega )}^{\frac{1}{s}}\)

On note les termes \({\tilde{\lambda}}_{m}\) de la suite ci-dessous, que l’on calcule en pratique par post-traitement à l’aide de \({u}_{m}\) (la puissance extérieure étant unitaire:

\({\tilde{\lambda}}_{m}={\parallel \Omega \parallel }_{m}^{1-\frac{1}{m}}{({\int}_{\Omega}A(m){(\varepsilon ({u}_{m}).\varepsilon ({u}_{m}))}^{m/2}d\Omega )}^{\frac{1}{m}}\text{}-{L}_{0}({u}_{m})\) éq 2.3-2

Cette suite \({\tilde{\lambda}}_{m}\) est donc décroissante pour \(m\to {1}^{+}\) et on montre [bib7] qu’elle converge vers, ce qui \({\lambda}_{\lim}\) permet un bon contrôle. Comme on peut minorer (sachant que

) le premier terme de [éq 2.3-1]:

\({\lambda}_{\lim}\le {\int}_{\Omega}{\sigma}_{y}\sqrt{\frac{2}{3}\varepsilon ({u}_{m}).\varepsilon ({u}_{m})}d\Omega \text{}-{L}_{0}({u}_{m})\le {\tilde{\lambda}}_{m}\)

on calcule donc pour chaque valeur de

(donc de l’instant \(t\) ) la charge limite \({\tilde{\lambda}}_{m}\) par excèsconvergeant aussi vers \({\lambda}_{\lim}\) :

\({\lambda}_{\lim}\le {\stackrel{ˆ}{\lambda}}_{m}={\int}_{\Omega}{\sigma}_{y}\sqrt{\frac{2}{3}\varepsilon ({u}_{m}).\varepsilon ({u}_{m})}d\Omega \text{}-{L}_{0}({u}_{m})\) éq 2.3-3

On juge de la qualité de l’approximation de la charge limite \({\lambda}_{\lim}\) par comparaison des différentes valeursde \({\tilde{\lambda}}_{m}\) qui convergent vers \({\lambda}_{\lim}\) par excès (en \(m\to {1}^{+}\) ). Ces termes sont calculés par intégration numérique aux points de Gauss des éléments finis.

Une autre interprétation de l’intérêt que cette suite apporte réside dans le fait qu’elle exploite directement l’expression de la fonction d’appui du convexe de résistance, c’est-à-dire la puissance dissipable dans les modes de ruine potentielle, appliquée aux solutions incompressibles et normalisées calculées \({u}_{m}\) .

Dans le cas où le chargement permanent est nul: \({L}_{0}=0\) , on peut facilement exploiter le champ de contraintes (quasiment statiquement admissible) calculé avec la solution \({u}_{m}\) et obtenir une valeur par estimation de la charge limite, qui serait nécessairement une vraie borne inférieure si l’équilibre était vérifié exactement (voir [bib4]).

On calcule donc la suite \({\underline{\lambda}}_{m}\) , qui n’a en revanche pas de propriétés de monotonie:

\({\underline{\lambda}}_{m}={\int}_{\Omega}\frac{A(m)}{m}.{(\sqrt{\varepsilon ({u}_{m}).\varepsilon ({u}_{m})})}^{m}d\Omega .{(\underset{x\in \Omega }{\text{Sup}}(\frac{\sqrt{\frac{3}{2}{\sigma}^{D}({u}_{m}).{\sigma}^{D}({u}_{m})}}{{\sigma}_{y}}))}^{-1}\le {\stackrel{ˆ}{\lambda}}_{m}\) éq 2.3-4

Cette maximisation (de la fonction appelée jauge du convexe de résistance) n’est calculée qu’aux points de Gauss des éléments finis. Aussi la valeur obtenue, pour chaque

, inférieure à \({\tilde{\lambda}}_{m}\) [bib4], ne peut être considérée que comme une indication.

Par contre, toujours dans le cas où les charges permanentes sont nulles, elle permet, avec la valeur par excès \({\tilde{\lambda}}_{m}\) , de fournir un encadrement de la charge limite du problème discrétisé.

Fonctionnalités et vérification

Pour réaliser un calcul dans Code_Aster en analyse limite avec la méthode de régularisation de Norton-Hoff-Friaâ avec le critère de résistance de VonMises, il faut:

• définir le modèle 2D (plan ou axis) ou 3D avec les éléments finis quasi-incompressibles, modélisations 3D_INCO_UPG, D_PLAN_INCO_UPG, ou AXIS_INCO_UPG;

• assurer la condition d’incompressibilité: GONF=0 dans AFFE_CHAR_MECA;

• définir uniquement la caractéristique du matériau \({s}_{y}\) , la charge limite étant indépendante de

  

et \(\nu\) ,

• définir le chargement permanent et celui qui est paramétré par \(\lambda\) ;

• définir la discrétisation en temps, (en pratique entre \({t}_{min}=1\) et \({t}_{max}=2\) à \(5\) );

• réaliser un calcul non-linéaire avec la relation de comportement NORTON_HOFF avec la commande STAT_NON_LINE [U4.51.03], et le pilotage ANA_LIM. On peut en pratique utiliser la recherche linéaire pour améliorer la convergence, et la subdivision du pas de temps,

• post-traiter le calcul pour obtenir la charge limite avec la commande POST_ELEM [U4.81.22].

L’utilisation de ces commandes est détaillée dans le document [U2.05.04].

En ce qui concerne le post-traitement, l’opérateur POST_ELEM produit alors une table qui donne pour chaque instant du calcul, c’est-à-dire pour des régularisations de plus en plus faibles, 2 paramètres:

• le paramètre “CHAR_LIMI_SUP” contient une borne supérieure de la charge limite, par intégration sur chaque élément fini et une somme sur l’ensemble des éléments du modèle: \({\stackrel{ˆ}{\lambda}}_{m}={\int}_{\Omega}{\sigma}_{y}\sqrt{\frac{2}{3}\varepsilon ({u}_{m}).\varepsilon ({u}_{m})}d\Omega \text{}-{L}_{0}({u}_{m})\)

• et, en l’absence de chargement constant, (CHAR_CSTE =”NON”), le paramètre “CHAR_LIMI_ESTIMEE” contient une estimation d’une borne inférieure \({\underline{\lambda}}_{m}\) correspondant à: \({\underline{\lambda}}_{m}={\int}_{\Omega}\frac{A(m)}{m}.{(\sqrt{\varepsilon ({u}_{m}).\varepsilon ({u}_{m})})}^{m}d\Omega .{(\underset{x\in \Omega }{\text{Sup}}(\frac{\sqrt{\frac{3}{2}{\sigma}^{D}({u}_{m}).{\sigma}^{D}({u}_{m})}}{{\sigma}_{y}}))}^{-1}\le {\stackrel{ˆ}{\lambda}}_{m}\)

Si un chargement constant est présent, (renseigner alors impérativement CHAR_CSTE = “OUI”), le paramètre PUIS_CHAR_CSTE représente la puissance du chargement constant dans le champ de vitesse solution du problème.

Plusieurs tests de vérification sont disponibles, en particulier le test SSNV124 [V6.04.124]. Sur ce problème très simple, un calcul analytique permet d’obtenir la charge limite exacte dans la direction du chargement, ainsi que les estimations produites par la méthode de régularisation. Pour plus de détails on se reportera à [bib4] et [bib5].

D’autre part des validations complémentaires ont été effectuées dans le cadre d’études comparatives, comme le benchmark européen LISA [bib8, bib10]: sur des calculs de charges limites en 2D, 2D axis et 3D, la méthode cinématique régularisée présentée ici permet de gagner un facteur de 6 à 10 sur le temps calcul par rapport à un calcul élastoplastique incrémental, et permet d’obtenir un encadrement de la charge limite, contrairement aux méthodes des autres participants.

Bibliographie

1. ANGLES J., VOLDOIRE F., Modélisation et calcul de la charge limite d’un composant fissuré, CR-MMN 1522-07, sept. 96.

2. FRIAA A., Loi de Norton-Hoff généralisée en plasticité et viscoplasticité, Thèse de doctorat, 1979.

3. FRIAA A., FREMOND M., Les méthodes statique et cinématique en calcul à la rupture et en analyse limite, Journal de Mécanique théorique et appliquée, Vol. 11, No 5, 881-905, 1982.

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

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

6. MICHEL-PONNELLE S., LORENTZ E. Éléments finis traitant la quasi-incompressibilité, document R3.08.06D.

7. VOLDOIRE F., Mise en œuvre de la méthode de régularisation de Norton-Hoff-Friaâ pour l’analyse limite des structures, note EDF/DER HI-74/97/026.

8. 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.

9. [R3.06.08] S. MICHEL-PONNELLE, E. LORENTZ, Éléments finis traitant la quasi-incompressibilité, 2005.

10. 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).

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
5	F.VOLDOIRE EDF-R&D/AMA	Texte initial
8.4	E.LORENTZ, S.MICHEL-PONNELLE EDF-R&D/AMA	Modification de la gestion de l’exposant de la loi de Norton-Hoff