r5.03.09 Relations de comportement non linéaires 1D

Résumé:

Ce document décrit les quantités calculées par l’opérateur STAT_NON_LINE nécessaires à la mise en œuvre de l’algorithme non linéaire quasi statique décrit en [R5.03.01] dans le cas des comportements élastoplastiques ou viscoplastiques monodimensionnels. Ces comportements, sauf mention contraire, sont applicables aux éléments de BARRE, aux éléments de poutre et poutres multifibres (direction axiales seulement) et aux éléments d’armature de béton (modélisation GRILLE).

Les comportements décrits dans ce document sont:

• le comportement de Von Mises à écrouissage isotrope linéaire: VMIS_ISOT_LINE, et quelconque VMIS_ISOT_TRAC,

• le comportement de Von Mises à écrouissage cinématique linéaire: VMIS_CINE_LINE,

• le comportement de Von Mises à écrouissage linéaire, non symétrique en traction et compression: avec restauration du centre du domaine élastique: VMIS_ASYM_LINE. Ce dernier a été développé pour modéliser l’action du sol sur les Câbles à Isolation Gazeuse,

• les comportements viscoplastiques avec effet de l’irradiation: VISC_IRRA_LOG, GRAN_IRRA_LOG.

• le comportement de MAZARS dans sa version \(\mathrm{1D}\) . La version \(\mathrm{1D}\) du modèle de mazars permet de rendre compte de la restauration de rigidité en cas de refermeture des fissures.

• le comportement pour modéliser la relaxation des câbles pré-contraints.

La résolution est faite cas par une méthode d’intégration implicite, sauf mention contraire, à partir de l’instant de calcul précédent, on calcule le champ de contraintes résultant d’un incrément de déformation, et le comportement tangent qui permet de construire les matrices tangentes.

On décrit enfin une méthode, similaire à la méthode due à R.deBorst [R5.03.03] permettant d’utiliser tous les comportements disponibles en 3D dans les éléments 1D.

Comportement de Von-Mises à écrouissage isotrope linéaire : VMIS_ISOT_LINE ou VMIS_ISOT_TRAC

Équations du modèle VMIS_ISOT_LINE

Elles sont la restriction du comportement 3D [R5.03.02] au cas uniaxial:

\((\begin{array}{c}\stackrel{ˉ}{\dot{{\varepsilon}^{p}}}=\frac{3}{2}\stackrel{ˉ}{\dot{p}}\frac{\tilde{\sigma}}{{\sigma}_{\mathrm{eq}}}=\stackrel{ˉ}{\dot{p}}\frac{\sigma}{\mid \sigma \mid }\\ \frac{\sigma}{E}=\varepsilon -{\varepsilon}^{p}-{\varepsilon}^{\mathrm{th}}\\ {\sigma}_{\mathrm{eq}}-R(p)=\mid \sigma \mid -R(p)\le 0\\ (\begin{array}{c}\stackrel{ˉ}{\dot{p}}=0\mathrm{si}{\sigma}_{\mathrm{eq}}-R(p)<0\\ \stackrel{ˉ}{\dot{p}}\ge 0\mathrm{si}{\sigma}_{\mathrm{eq}}-R(p)=0\end{array}\end{array}\)

avec:

\(\stackrel{ˉ}{\dot{{\varepsilon}^{p}}}\)	vitesse de déformation plastique,
\(p\)	déformation plastique cumulée,
\({\varepsilon}^{\mathit{th}}=\alpha (T-{T}_{\mathit{ref}})\)	déformation thermique,
\(R(p)=\frac{E{E}_{T}}{E-{E}_{T}}p+{\sigma}_{y}\)	fonction d’écrouissage linéaire isotrope, ou \(R(p)\) bien affine par morceaux, déduite de la courbe de traction.

Dans le cas VMIS_ISOT_LINE, les données des caractéristiques de matériaux sont celles fournies sous le mot clé facteur ECRO_LINE ou ECRO_LINE_FO de l’opérateur DEFI_MATERIAU [U4.43.01].

/ ECRO_LINE = ( D_SIGM_EPSI = \({E}_{T}\) , SY = \({\sigma}_{y}\) )

/ ECRO_LINE_FO = ( D_SIGM_EPSI = \({E}_{T}\) , SY = \({\sigma}_{y}\) )

Dans le cas VMIS_ISOT_TRAC, les données des caractéristiques des matériaux sont fournies sous le mot clé facteur TRACTION de l’opérateur DEFI_MATERIAU [U4.43.01].

TRACTION = _F (SIGM = courbe_traction)

courbe_traction représente la courbe de traction, point par point. Le premier point permet de définir la limite d’élasticité \({\sigma}_{y}\) et le module d’Young \(E\) [R5.03.02].

ECRO_LINE_FO correspond au cas où \({E}_{T}\) et \({\sigma}_{y}\) dépendent de la température et sont alors calculés pour la température du point de Gauss courant. Le module d’Young \(E\) et le coefficient de Poisson \(v\) sont ceux fournis sous les mots clés facteurs ELAS ou ELAS_FO. Dans ce cas la courbe de traction est la suivante:

\(\lbrace \begin{array}{c}\begin{array}{cc}{\sigma}_{L}=E{\varepsilon}_{L}& \mathrm{si}{\varepsilon}_{L}<\frac{{\sigma}_{y}}{E}\\ {\sigma}_{L}={\sigma}_{y}+{E}_{T}({\varepsilon}_{L}-\frac{{\sigma}_{y}}{E})& \mathrm{si}{\varepsilon}_{L}\ge \frac{{\sigma}_{y}}{E}\end{array}\end{array}\)

Lorsque le critère est atteint on a:

\({\sigma}_{L}-R(p)=0\) , donc \({\sigma}_{L}-R({\varepsilon}_{L}-\frac{{\sigma}_{L}}{E})=0\) d’où:

\(R(p)=\frac{{E}_{T}E}{E-{E}_{T}}p+{\sigma}_{y}=Hp+{\sigma}_{y}\)

Dans le cas d’une courbe de traction, la démarche est identique à [R5.03.01].

Intégration de la relation VMIS_ISOT_LINE

Par discrétisation implicite directe des relations de comportement, de façon analogue à l’intégration 3D [R5.03.02] on obtient:

\(\lbrace \begin{array}{c}∣{\sigma}^{-}+\Delta \sigma ∣-R({p}^{-}+\Delta p)\le 0\\ E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})-({\sigma}^{-}+\Delta \sigma )+\frac{E}{{E}^{-}}{\sigma}^{-}=E\Delta p\frac{{\sigma}^{-}+\Delta \sigma }{∣{\sigma}^{-}+\Delta \sigma ∣}\\ \Delta p\ge 0\mathrm{si}∣{\sigma}^{-}+\Delta \sigma ∣=R({p}^{-}+\Delta p)\\ \Delta p=0\mathrm{si}∣{\sigma}^{-}+\Delta \sigma ∣<R({p}^{-}+\Delta p)\end{array}\)

Deux cas se présentent:

• \(∣{\sigma}^{-}+\Delta \sigma ∣<R({p}^{-}+\Delta p)\)

dans ce cas \(\Delta p=0\) soit \(\sigma =E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})+\frac{E}{{E}^{-}}{\sigma}^{-}\)

donc \(∣{\sigma}^{-}\frac{E}{{E}^{-}}+E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})∣<R({p}^{-})\)

• \(∣{\sigma}^{-}+\Delta \sigma ∣=R({p}^{-}+\Delta p)\)

dans ce cas \(\Delta p\ge 0\)

donc \(∣\frac{{\sigma}^{-}E}{{E}^{-}}+E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})∣\ge R({p}^{-})\) .

On en déduit l’algorithme de résolution:

posons \({\sigma}^{e}=\frac{E{\sigma}^{-}}{{E}^{-}}+E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})\)

si \(∣{\sigma}^{e}∣\le R({p}^{-})\) alors \(\Delta p=0\) et \(\Delta \sigma =E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})\)

si \(∣{\sigma}^{e}∣>R({p}^{-})\) alors il faut résoudre:

\(\begin{array}{}{\sigma}^{e}={\sigma}^{-}+\Delta \sigma +E\Delta p\frac{{\sigma}^{-}+\Delta \sigma }{∣{\sigma}^{-}+\Delta \sigma ∣}\\ {\sigma}^{e}=(1+\frac{E\Delta p}{∣{\sigma}^{-}+\Delta \sigma ∣})({\sigma}^{-}+\Delta \sigma )\end{array}\)

donc en prenant la valeur absolue:

\(∣{\sigma}^{e}∣=(1+\frac{E\Delta p}{∣{\sigma}^{-}+\Delta \sigma ∣})({\sigma}^{-}+\Delta \sigma )\)

soit, en utilisant \(\begin{array}{}∣{\sigma}^{-}+\Delta \sigma ∣=R({p}^{-}+\Delta p)\\ ∣{\sigma}^{e}∣=R({p}^{-}+\Delta p)+E\Delta p\end{array}\) .

On en déduit donc:

• dans le cas d’un écrouissage linéaire: \(\Delta p=\frac{∣{\sigma}^{e}∣-({\sigma}_{y}+{\mathrm{Hp}}^{-})}{E+H}\)

• et dans le cas d’un écrouissage quelconque, la courbe \(R(p)\) étant affinée par morceaux, on résout directement l’équation en \(\Delta p\) : de la même façon \(E\Delta p+R({p}^{-}+\Delta p)=∣{\sigma}^{e}∣\) qu’en 3D [R5.03.02].

Remarquons au passage que: \(\frac{{\sigma}^{e}}{∣{\sigma}^{e}∣}=\frac{\sigma}{R(p)}\)

alors \(\sigma =({\sigma}^{-}+\Delta \sigma )=\frac{{\sigma}^{e}}{∣{\sigma}^{e}∣}R(p)=\frac{{\sigma}^{e}}{1+\frac{E\Delta p}{R(p)}}\)

De plus, l’option FULL_MECA permet de calculer la matrice tangente \({K}_{i}^{n}\) à chaque itération. L’opérateur tangent qui sert à la construire est calculé directement sur le système discrétisé précédent. On obtient directement:

si \(∣{\sigma}^{e}∣>R({p}^{-})\) \(\frac{\delta \sigma }{\delta \varepsilon }={E}_{T}\)

sinon \(\frac{\delta \sigma }{\delta \varepsilon }=E\)

Remarque:

L’option RIGI_MECA_TANG qui permet de calculer la matrice tangente \({K}_{i}^{0}\) utilisée dans la phase de prédiction de l’algorithme de Newton, tient compte de l’indicateur de plasticité à l’instant précédent:

si \(\chi =1\frac{\delta \sigma }{\delta \varepsilon }={E}_{T}\) si \(\chi =0\frac{\delta \sigma }{\delta \varepsilon }=E\)

Variables internes

La relation de comportement VMIS_ISOT_LINE produit deux variables internes: \(p\) et \(\chi\)

Comportement de Von Mises, écrouissage cinématique linéaire 1D : VMIS_CINE_LINE

Équation du modèle VMIS_CINE_LINE

Pour des raisons de performances la relation est écrite en 1D. Elles sont la restriction du comportement 3D ([R5.03.02] et [R5.03.16]) au cas uniaxial. Le comportement 3D s’écrit:

\(\sigma =K(\varepsilon -{\varepsilon}^{p}-{\varepsilon}^{\mathrm{th}})\) avec \(K\) opérateur d’élasticité

\(X=C{\varepsilon}^{p}\)

\(F(\sigma ,R,X)={(\tilde{\sigma}-X)}_{\mathrm{eq}}-{\sigma}_{y}\) avec \({A}_{\mathrm{eq}}=\sqrt{\frac{3}{2}\tilde{A}\cdot \tilde{A}}\)

\(\dot{{\varepsilon}^{p}}=\dot{p}\frac{\partial F}{\partial \sigma }=\frac{3}{2}\dot{p\frac{\tilde{\sigma}-X}{{(\tilde{\sigma}-X)}_{\mathrm{eq}}}}\) \(\lbrace \begin{array}{c}\mathrm{si}F<0\dot{p}=0\\ \mathrm{si}F=0\dot{p}\ge 0\end{array}\)

Dans le cas uniaxial, les tenseurs s’écrivent:

\(\tilde{\sigma}=\sigma DX=XD{\varepsilon}^{p}=\frac{3}{2}{\varepsilon}^{p}D\) avec \(D=\left[\begin{array}{ccc}2/3& \text{}& \text{}\\ \text{}& -1/3& \text{}\\ \text{}& \text{}& -1/3\end{array}\right]\)

Tant que le chargement est monotone, on obtient immédiatement les relations suivantes:

\(p={\varepsilon}^{p}X=\frac{3}{2}C{\varepsilon}^{p}\sigma =\frac{3}{2}C{\varepsilon}^{p}+{\sigma}_{y}\sigma =F(\varepsilon )={\sigma}_{y}+\frac{E\cdot {E}_{T}}{E-{E}_{T}}p\)

C est déterminée par: \(C=\frac{2}{3}\frac{E{E}_{T}}{E-{E}_{T}}\) . On pose: \(H=\frac{E{E}_{T}}{E-{E}_{T}}=\frac{3}{2}C\)

La relation de comportement 1D s’écrit alors:

\((\begin{array}{c}∣{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X∣\\ E\Delta {\varepsilon}^{p}=E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})-({\sigma}^{-}+\Delta \sigma )+\frac{E}{{E}^{-}}{\sigma}^{-}\\ X=\frac{3}{2}C{\varepsilon}^{p}=H{\varepsilon}^{p}\\ ∣\sigma -X∣-{\sigma}_{y}\le 0\\ (\begin{array}{c}\dot{p}=0\mathrm{si}∣\sigma -X∣-{\sigma}_{y}<0\\ \dot{p}\ge 0\mathrm{si}∣\sigma -X∣-{\sigma}_{y}=0\end{array}\end{array}\)

Les données des caractéristiques de matériaux sont celles fournies sous le mot clé facteur ECRO_LINE ou ECRO_LINE_FO de l’opérateur DEFI_MATERIAU [U4.43.01]:

/ ECRO_LINE = ( D_SIGM_EPSI = \({E}_{T}\) , SY = \({\sigma}_{y}\) )

/ ECRO_LINE_FO = ( D_SIGM_EPSI = \({E}_{T}\) , SY = \({\sigma}_{y}\) )

Intégration de la relation VMIS_CINE_LINE

Par discrétisation implicite directe des relations de comportement, de façon analogue à l’intégration 3D ([R5.03.02] et [R5.03.16]) on obtient:

\(\lbrace \begin{array}{c}∣{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X∣-{\sigma}_{y}\le 0\\ E\Delta {\varepsilon}^{p}=E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})-({\sigma}^{-}+\Delta \sigma )+\frac{E}{{E}^{-}}{\sigma}^{-}\\ \Delta {\varepsilon}^{p}=\Delta p\frac{{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X}{∣{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X∣}\\ \frac{X}{H}-\frac{{X}^{-}}{{H}^{-}}=\Delta {\varepsilon}^{p}\\ \Delta p\ge 0\mathrm{si}∣{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X∣={\sigma}_{y}\\ \Delta p=0\mathrm{si}∣{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X∣<{\sigma}_{y}\end{array}\)

avec \(\Delta {\varepsilon}^{\mathrm{th}}=\alpha (T-{T}_{\mathrm{ref}})-{\alpha}^{-}({T}^{-}-{T}_{\mathrm{ref}})\)

Deux cas se présentent:

• \(∣{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X∣<{\sigma}_{y}\) dans ce cas \(\Delta p=0\) soit \(\sigma =E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})+\frac{E}{{E}^{-}}{\sigma}^{-}-\frac{H}{{H}^{-}}{X}^{-}\) donc \(∣{\sigma}^{-}\frac{E}{{E}^{-}}-{X}^{-}\frac{H}{{H}^{-}}+E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})∣<R({p}^{-})\) ,

• sinon \(\Delta p\ge 0\) .

Pour simplifier les écritures on posera: \({\sigma}^{e}=\frac{E}{{E}^{-}}{\sigma}^{-}-\frac{H}{{H}^{-}}{X}^{-}+E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})\) .

On en déduit l’algorithme de résolution:

1. si \(∣{\sigma}^{e}∣\le {\sigma}_{y}\) alors \(\Delta p=0,X={X}^{-}\frac{H}{{H}^{-}},\sigma =E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})+\frac{E}{{E}^{-}}{\sigma}^{-}\)

2. sinon il faut résoudre:

\(\lbrace \begin{array}{c}E\Delta {\varepsilon}^{p}=E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})-\Delta \sigma ={\sigma}^{e}-({\sigma}^{-}+\Delta \sigma )+{X}^{-}\frac{H}{{H}^{-}}\\ \Delta {\varepsilon}^{p}=\Delta p\frac{{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X}{∣{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X∣}=\Delta p\frac{\sigma -X}{∣\sigma -X∣}\\ X-\frac{H}{{H}^{-}}{X}^{-}=H\Delta {\varepsilon}^{p}\\ ∣{\sigma}^{-}+\Delta \sigma -{X}^{-}-\Delta X∣-{\sigma}_{y}=0\end{array}\)

Remarquons que: \(\frac{H}{{H}^{-}}{X}^{-}=X-H\Delta {\varepsilon}^{p}\) .

On déduit alors de la première équation: \({\sigma}^{e}=\sigma -X+(E+H)\Delta {\varepsilon}^{p}\)

On obtient donc, en éliminant \(\sigma -X\) de la deuxième équation:

\(\Delta {\varepsilon}^{p}={\sigma}^{e}\frac{\Delta p}{(E+H)\Delta p+{\sigma}_{y}}\)

En remplaçant \(\Delta {\varepsilon}^{p}\) dans la relation entre \({\sigma}^{e}\) et \(\sigma -X\) , on obtient:

\(\sigma -X={\sigma}^{e}(\frac{{\sigma}_{y}}{(E+H)\Delta p+{\sigma}_{y}})\)

En prenant la valeur absolue des deux membres de l’équation précédente, on trouve \(\Delta p\) :

\((E+H)\Delta p+{\sigma}_{y}=∣{\sigma}^{e}∣\)

Une fois \(\Delta p\) déterminé, on peut calculer:

\(\begin{array}{}\Delta {\varepsilon}^{p}=\Delta p\frac{{\sigma}^{e}}{∣{\sigma}^{e}∣}\\ X={X}^{-}+\Delta X=\frac{H{X}^{-}}{{H}^{-}}+H\Delta p\frac{{\sigma}^{e}}{∣{\sigma}^{e}∣}\end{array}\)

et en utilisant: \(\frac{\sigma -X}{{\sigma}_{y}}=\frac{{\sigma}^{e}}{∣{\sigma}^{e}∣}\) , on obtient directement: \(\sigma ={\sigma}_{y}\frac{{\sigma}^{e}}{∣{\sigma}^{e}∣}+X\)

De plus, l’option FULL_MECA permet de calculer la matrice tangente \({K}_{i}^{n}\) à chaque itération. L’opérateur tangent qui sert à la construire est calculé directement sur le système discrétisé précédent. On obtient directement:

si \(∣{\sigma}^{e}∣>R({p}^{-})\) \(\frac{\delta \partial }{\delta \varepsilon }={E}_{T}\)

sinon \(\frac{\delta \sigma }{\delta \varepsilon }=E\)

L’option RIGI_MECA_TANG qui permet de calculer la matrice tangente \({K}_{i}^{0}\) utilisée dans la phase de prédiction de l’algorithme de Newton est obtenue à l’aide de l’indicateur de plasticité \({\chi }^{-}\) de l’instant précédent:

• si \({\chi }^{-}=1\) alors \(\frac{\delta \sigma }{\delta \varepsilon }={E}_{T}\)

• si \({\chi }^{-}=0\) alors \(\frac{\delta \sigma }{\delta \varepsilon }=E\)

Variables internes

La relation de comportement VMIS_CINE_LINE produit deux variables internes: \(X\) .et \(\chi\)

Comportement de Von Mises, écrouissage cinématique linéaire 1D : vmis_CINE_gc

Équation du modèle VMIS_CINE_GC

Pour des raisons de performances la relation est également écrite en 1D pour une utilisation avec des éléments finis de type poutre multifibre. Les équations sont issues de la restriction du comportement 3D ([R5.03.02] et [R5.03.16]) au cas uniaxial.

Les équations du modèle sont les mêmes que celles du § 3.1 .

Les données des matériaux sont celles fournies sous le mot clé facteur ECRO_LINE de l’opérateur DEFI_MATERIAU[U4.43.01]:

/ECRO_LINE=_F(

♦D_SIGM_EPSI = \({E}_{T}\) [Réel]

♦SY = \({\sigma}_{y}\) [Réel]

◊SIGM_LIM =sigmlim [Réel]

◊EPSI_LIM =epsilim [Réel]

)

Les opérandes sigm_lim et espi_lim permettent de définir les bornes qui correspondent aux états limites de service et ultime, classiquement utilisées lors d’étude en génie civil.

◊SIGM_LIM=sigmlim

Définition de la contrainte limite.

◊EPSI_LIM=epslim

Définition de la déformation limite.

Ces bornes sont obligatoires lorsque l’on utilise le comportement VMIS_CINE_GC (Cf. [U4.51.11] Comportements non-linéaires, [U4.42.07] DEFI_MATER_GC). Dans les autres cas elles ne sont pas prises en compte.

Intégration de la relation VMIS_CINE_GC

La méthode d’intégration est identique à celle présentées au § 3.2 .

Variables internes

La modélisation supportée est 1D, le nombres de variables internes est de 6.

• \(\mathit{V1}\) :Cette variable représente la contrainte divisée par la contrainte limite sigmlim.

• \(\mathit{V2}\) :Cette variable représente la déformation totale divisée par la déformation limite epslim.

• \(\mathit{V3}\) :Écrouissage cinématique: XCINXX. En 1D seul un scalaire est nécessaire.

• \(\mathit{V4}\) :Indicateur plastique: INDIPLAS. Indique si le matériau a dépassé le critère élastique.

• \(\mathit{V5}\) :dissipation non récupérable: DISSIP. Lors de calculs sismiques il peut être utile à l’utilisateur de connaître l’énergie dissipée non récupérable. La variable dissip représente le cumul d’énergie non récupérable. L’incrément d’énergie non récupérable s’écrit sous la forme:

\(\Delta \mathit{Eg}=\frac{1}{2}({E}^{+}\Delta \varepsilon –({\sigma}^{+}-{\sigma}^{-}))\Delta \varepsilon\)

• \(\mathit{V6}\) :dissipation thermodynamique: DISSTHER. L’incrément de dissipation thermodynamique s’écrit sous la forme: \(\Delta \mathit{Eg}={\sigma}_{y}\dot{p}\) .

Comportement de Von Mises à écrouissage linéaire asymétrique : VMIS_ASYM_LINE

Équations du modèle VMIS_ASYM_LINE

Comportement asymétrique en traction et en compression

C’est un comportement découplé en traction et compression, construit à partir de VMIS_ASYM_LINE, mais avec des limites d’élasticité et des modules d’écrouissage différents en traction et en compression. Nous adoptons un indice \(T\) pour la traction et \(C\) pour la compression. Le comportement élastique en traction et compression est identique et caractérisé par le même module d’Young. Il y a deux domaines d’écrouissage isotrope définis par \({R}_{T}\) et \({R}_{C}\) . Les deux domaines sont indépendants l’un de l’autre.

\(\mathrm{YT}\)	limite d’élasticité en traction. En valeur absolue.
\(\mathrm{YC}\)	limite d’élasticité en compression. En valeur absolue.
\({p}_{T}\)	Variable interne en traction. Valeur algébrique.
\({p}_{C}\)	Variable interne en compression. Valeur algébrique.
\({E}_{\mathrm{TT}}\)	Pente d’écrouissage en traction.
\({E}_{\mathrm{TC}}\)	Pente d‘écrouissage en compression.

Les équations du modèle de comportement sont:

\(\lbrace \begin{array}{c}{\dot{\varepsilon}}^{p}=\dot{\varepsilon}-{E}^{-1}\sigma -\dot{{\varepsilon}^{\mathrm{th}}}\\ {\dot{\varepsilon}}^{p}={\dot{\varepsilon}}_{C}^{p}+{\dot{\varepsilon}}_{T}^{p}\\ {\dot{\varepsilon}}_{C}^{p}={\dot{\varepsilon}}_{C}\frac{\sigma}{∣\sigma ∣}\\ {\dot{\varepsilon}}_{C}^{p}={\dot{\varepsilon}}_{T}\frac{\sigma}{∣\sigma ∣}\\ \sigma -{R}_{T}({p}_{T})\le 0\\ -\sigma -{R}_{C}({p}_{C})\le 0\end{array}\mathrm{avec}(\begin{array}{c}{\dot{p}}_{C}=0\mathrm{si}-\sigma -{R}_{C}({p}_{C})<0\\ {\dot{p}}_{C}\ge 0\mathrm{si}-\sigma ={R}_{C}({p}_{C})\\ {\dot{p}}_{T}=0\mathrm{si}\sigma -{R}_{T}({p}_{T})<0\\ {\dot{p}}_{T}\ge 0\mathrm{si}\sigma ={T}_{T}({p}_{T})\end{array}\)

\({\dot{\varepsilon}}_{C}^{p}\) : vitesse de déformation plastique en compression,

\({\dot{\varepsilon}}_{T}^{p}\) : vitesse de déformation plastique en traction,

\({\dot{\varepsilon}}^{\mathrm{th}}\) : déformation d’origine thermique: \({\varepsilon}^{\mathrm{th}}=\alpha (T-{T}_{\mathrm{ref}})\)

On remarque que l’on ne peut avoir simultanément plastification en traction et en compression: soit \({\dot{p}}_{C}=0\) , soit \({\dot{p}}_{C}=0\) , soit les deux sont nulles.

Les données des caractéristiques de matériaux sont celles fournies sous le mot clé facteur ECRO_ASYM_LINE de l’opérateur DEFI_MATERIAU [U4.43.01].

ECRO_ASYM_LINE = _F(

DT_SIGM_EPSI = \({E}_{\mathrm{TT}}\) , SY_T= \({\sigma}_{\mathrm{yT}}\) ,

DC_SIGM_EPSI = \({E}_{\mathrm{TC}}\) , SY_C= \({\sigma}_{\mathrm{yC}}\) ,)

Le module d’Young E est fourni sous les mots clés facteurs ELAS ou ELAS_FO.

On calcule les fonctions d’écrouissage par: \(\begin{array}{}{R}_{T}(p)=\frac{{E}_{\mathrm{TT}}E}{E-{E}_{\mathrm{TT}}}{p}_{T}+{\sigma}_{\mathrm{yT}}={H}_{T}\cdot {p}_{T}+{\sigma}_{\mathrm{yT}}\\ {R}_{C}(p)=\frac{{E}_{\mathrm{TC}}E}{E-{E}_{\mathrm{TC}}}{p}_{C}+{\sigma}_{\mathrm{yC}}={H}_{C}\cdot {p}_{C}+{\sigma}_{\mathrm{yC}}\end{array}\)

Intégration du comportement VMIS_ASYM_LINE

Par discrétisation implicite directe de la relation de comportement asymétrique, de façon analogue à la précédente, on obtient:

\(\lbrace \begin{array}{c}\Delta {\varepsilon}^{p}=\Delta {\varepsilon}_{T}^{p}+\Delta {\varepsilon}_{C}^{p}\\ \Delta {\varepsilon}^{p}=\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}}-\frac{\Delta \sigma }{E}\\ \Delta {\varepsilon}_{T}^{p}=\Delta {p}_{T}\frac{{\sigma}^{-}+\Delta \sigma }{∣{\sigma}^{-}+\Delta \sigma ∣}\\ \Delta {p}_{T}\ge 0\mathrm{si}({\sigma}^{-}+\Delta \sigma )-{R}_{T}(\stackrel{ˉ}{{p}_{T}}+\Delta \stackrel{ˉ}{{p}_{T}})\le 0\\ \Delta {p}_{T}=0\mathrm{si}({\sigma}^{-}+\Delta \sigma )-{R}_{T}(\stackrel{ˉ}{{p}_{T}}+\Delta \stackrel{ˉ}{{p}_{T}})<0\\ \text{}\\ \Delta {\varepsilon}_{C}^{p}=\Delta {p}_{C}\frac{{\sigma}^{-}+\Delta \sigma }{∣{\sigma}^{-}+\Delta \sigma ∣}\\ -({\sigma}^{-}+\Delta \sigma )-{R}_{C}(\stackrel{ˉ}{{p}_{C}}+\Delta \stackrel{ˉ}{{p}_{C}})\le 0\\ \Delta {p}_{C}\ge 0\mathrm{si}-({\sigma}^{-}+\Delta \sigma )-{R}_{C}(\stackrel{ˉ}{{p}_{C}}+\Delta \stackrel{ˉ}{{p}_{C}})=0\\ \Delta {p}_{C}=0\mathrm{si}-({\sigma}^{-}+\Delta \sigma )-{R}_{C}(\stackrel{ˉ}{{p}_{C}}+\Delta \stackrel{ˉ}{{p}_{C}})<0\end{array}\)

L’intégration est similaire à celle de VMIS_ISOT_LINE pour chacune des directions de traction et de compression. Il faut bien voir que les centres des domaines d’élasticité sont des données (calculées explicitement au pas précédent) pour le problème incrémental à résoudre.

Quatre cas se présentent:

• \(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}}>0\) on pose \({\sigma}_{T}^{e}={\sigma}^{-}+E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})\)

◦ si \({\sigma}_{T}^{e}<{R}_{T}(\stackrel{ˉ}{{p}_{T}})\) dans ce cas \(\Delta {p}_{T}=0\) donc \(\sigma ={\sigma}_{T}^{e}\) et \(\frac{\delta \sigma }{\delta \varepsilon }=E\)

◦ sinon: \(\Delta {p}_{T}=\frac{∣{\sigma}_{T}^{e}∣-({\sigma}_{\mathrm{yT}}+{H}_{T}{p}_{T}^{-})}{E+{H}_{T}}\) , \(\Delta {p}_{C}=0\) \(\begin{array}{}\sigma =\frac{{\sigma}_{T}^{e}}{1+\frac{E\Delta {p}_{T}}{{R}_{T}({p}_{T})}}=\frac{{\sigma}_{T}^{e}}{∣{\sigma}_{T}^{e}∣}{R}_{T}({p}_{T})\\ \frac{\delta \sigma }{\delta \varepsilon }={E}_{\mathrm{TT}}\end{array}\)

• \(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}}<0\) on pose \({\sigma}_{C}^{e}={\sigma}^{-}+E(\Delta \varepsilon -\Delta {\varepsilon}^{\mathrm{th}})\)

◦ si \(-{\sigma}_{C}^{e}<{R}_{C}(\stackrel{ˉ}{{p}_{C}})\) dans ce cas \(\Delta {p}_{C}=0\) donc \(\sigma ={\sigma}_{C}^{e}\) et \(\frac{\delta \sigma }{\delta \varepsilon }=E\)

◦ sinon: \(\Delta {p}_{C}=\frac{∣{\sigma}_{C}^{e}∣-({\sigma}_{\mathrm{yC}}+{H}_{C}{p}_{C}\stackrel{ˉ}{\rbrace )}{E+{H}_{C}}\) , \(\Delta {p}_{T}=0\) \(\begin{array}{}\sigma =\frac{{\sigma}_{C}^{e}}{1+\frac{E\Delta {p}_{C}}{{R}_{C}({p}_{C})}}=\frac{{\sigma}_{C}^{e}}{∣{\sigma}_{C}^{e}∣}{R}_{C}({p}_{C})\\ \frac{\delta \sigma }{\delta \varepsilon }={E}_{\mathrm{TC}}\end{array}\)

Remarque:

La matrice tangente initiale (option RIGI_MECA_TANG)est prise égale à la matrice élastique.

Variables internes

La relation de comportement VMIS_ASYM_LINE produit 2 variables internes: \({p}_{C}\) \({p}_{T}\) .

Elle n’est pas utilisable pour les éléments de grille.

Comportements VISC_IRRA_LOG et GRAN_IRRA_LOG

Le modèle présenté dans ce chapitre décrit les comportements viscoplastiques 1D VISC_IRRA_LOG et GRAN_IRRA_LOG (fluage et grandissement sous irradiation des alliages M5 et Zircaloy-4) pour la modélisation des assemblages combustibles, et applicable aux éléments de barres et poutres multi-fibres.

Formulation du modèle

Les équations sont les suivantes:

\(\lbrace \begin{array}{c}{\dot{\epsilon}}^{\mathit{vp}}=\dot{\epsilon}\frac{\sigma}{\left|\sigma \right|}\\ \dot{\epsilon}=\left|\sigma \right|.\left({e}^{\frac{-Q}{T}}\right).\dot{\Phi}\left(\frac{A\omega }{1+\omega \Phi }+B-C\omega {e}^{-\omega t}\right)\\ \frac{\dot{\sigma}}{E}=\dot{\epsilon}-{\dot{\epsilon}}^{\mathit{vp}}-{\dot{\epsilon}}^{g}-\dot{{\epsilon}^{\mathit{th}}}\end{array}\)

Ces relations sont déduites des essais de fluage FLETAN et REFLET [8] pour diverses valeurs de flux neutronique.

Les coefficients sont fournis sous le mot clé VISC_IRRA_LOG ou GRAN_IRRA_LOGde DEFI_MATERIAU et \(\Phi\) est la fluence neutronique (intégrale du flux par rapport au temps).

\({\varepsilon}^{g}\) représente la déformation de grandissement sous flux. Elle n’est prise en compte que dans le comportement GRAN_IRRA_LOG et s’exprime sous la forme:

\({\varepsilon}^{g}(t)=f(T,{\Phi}_{t}(x,y,z))\)

Remarques:

1. La fluence neutronique \({\Phi}_{t}(x,y,z)\) s’exprime obligatoirement en \({10}^{20}n/{\mathit{cm}}^{2}\) *. Par convention dans* DEFI_MATERIAU [U4.43.01], si la valeur fournie sous le mot-clé FLUX_PHI est égale à 1, c’est le champ de fluence qui est utilisé pour le comportement. Dans le cas contraire, la valeur fournie dans DEFI_MATERIAU est utilisée comme flux neutronique constant.

2. C’est un champ aux nœuds défini comme variable de commande dans la commande AFFE_MATERIAU .

3. Attention : Le champ d’irradiation est incrémental et correspond à l’historique d’irradiation (stockée en variable interne – cf ci-dessous) auquel on ajoute l’incrément du champ de fluence venant de la variable de commande.

Variables internes

Trois variables internes:

• \(\mathit{V1}\) : la déformation viscoplastique cumulée: \({\varepsilon}_{p}\) ;

• \(\mathit{V2}\) : mémorisation de l’historique d’irradiation (fluence).

• \(\mathit{V3}\) : la déformation de grandissement: \({\varepsilon}^{g}\) .

Intégration implicite

Par discrétisation implicite directe des relations de comportement, on obtient:

\(\lbrace \begin{array}{c}\Delta {\epsilon}^{\mathit{vp}}=\Delta p\frac{\sigma \left({t}^{-}+\Delta t\right)}{\left|\sigma \left({t}^{-}+\Delta t\right)\right|}\\ \Delta p=\left|\sigma \left({t}^{-}+\Delta t\right)\right|\left({e}^{\frac{-Q}{T}}\right).\left(\frac{A\omega }{1+\omega \Phi ({t}^{-}+\Delta t)}+B-C\omega {e}^{-\omega t}\right)\Delta \Phi \\ \frac{\sigma}{E}-\frac{{\sigma}^{-}}{{E}^{-}}=\Delta \epsilon -\Delta {\epsilon}^{\mathit{vp}}-\Delta {\epsilon}^{g}-\Delta {\epsilon}^{\mathit{th}}\\ \text{avec}\\ \Delta {\epsilon}^{\mathit{th}}=\alpha (T)\left(T-{T}_{\mathit{ref}}\right)-\alpha ({T}^{-})\left(T-{T}_{\mathit{ref}}\right)\\ \Delta {\epsilon}^{g}=f({T}^{+},{\Phi}_{t}^{+})-f({T}^{-},{\Phi}_{t}^{-})\end{array}\)

On peut résoudre explicitement ces équations en posant: \({\sigma}^{e}=\frac{E}{{E}^{-}}{\sigma}^{-}+E(\Delta \varepsilon -\Delta {\varepsilon}^{g}-\Delta {\varepsilon}^{\mathrm{th}})\)

alors le système se réduit à: :math:`sigma ={sigma}^{e}-Esigma left({e}^{frac{-Q}{T}}right).left(frac{Aomega }{1+omega Phi }+B-Comega {e}^{-omega t}right)Delta Phi `

donc la solution s’obtient immédiatement: \(\sigma =\frac{{\sigma}^{e}}{1+E\left({e}^{\frac{-Q}{T}}\right).\left(\frac{A\omega }{1+\omega \Phi }+B-C\omega {e}^{-\omega t}\right)\Delta \Phi }\)

et l’opérateur tangent s’écrit: \(\frac{\partial \sigma }{\partial \epsilon }=\frac{E}{1+E\left({e}^{\frac{-Q}{T}}\right).\left(\frac{A\omega }{1+\omega \Phi }+B-C\omega {e}^{-\omega t}\right)\Delta \Phi }\)

Modèle de MAZARS en 1D

Équations du modèle

L’objectif de cette modélisation est de rendre compte de la refermeture des fissures. Ce modèle n’est utilisé qu’avec les poutres multifibres.Les équations présentées dans le document [R7.01.08] “Modèle d’endommagement de Mazars »sont repriseset réécrites en \(\mathrm{1D.}\)

\(\lbrace \begin{array}{c}{\sigma}_{xx}=(1-{D}_{t})E{\langle {\varepsilon}_{xx}^{e}\rangle }_{+}\\ {\sigma}_{xx}=(1-{D}_{c})E{\langle {\varepsilon}_{xx}^{e}\rangle }_{-}\end{array}\)

avec:

• \(E\) :module d’Young,

• \({D}_{t}\) :la variable d’endommagement en traction.

• \({D}_{c}\) :la variable d’endommagement en compression.

• \({\varepsilon}_{xx}^{e}\) :la déformation élastique \({\varepsilon}_{xx}^{e}=\varepsilon -{\varepsilon}^{\mathit{th}}\)

• \({\varepsilon}^{\mathit{th}}=\alpha (T-{T}_{\mathit{ref}})\) : la dilatation thermique

L’endommagement reste toujours piloté par les extensions.

L’endommagement de traction existe si \({\epsilon}_{\mathit{eq}}\ge {\epsilon}_{t0}\) , de compression existe si \({\epsilon}_{\mathit{eq}}\ge {\epsilon}_{c0}\) :

\({D}_{c}=1-\frac{{\epsilon}_{c0}\left(1-{A}_{c}\right)}{{\epsilon}_{\mathit{eq}}}-\frac{{A}_{c}}{\exp\left[{B}_{c}({\epsilon}_{\mathit{eq}}-{\epsilon}_{c0})\right]}\phantom{\rule{6em}{0ex}}{D}_{c}\in [0,1[\)

\({D}_{t}=1-\frac{{\epsilon}_{t0}\left(1-{A}_{t}\right)}{{\epsilon}_{\mathit{eq}}}-\frac{{A}_{t}}{\exp\left[{B}_{t}({\epsilon}_{\mathit{eq}}-{\epsilon}_{t0})\right]}\phantom{\rule{6em}{0ex}}{D}_{t}\in [0,1[\)

où \({A}_{c},{A}_{t},{B}_{c},{B}_{t},{\epsilon}_{c0},{\epsilon}_{t0}\) sont des paramètres matériaux à identifier.

L’endommagement est piloté par la déformation équivalente \({\varepsilon}_{\mathit{eq}}\) . Les extensions sont primordiales dans le phénomène de fissuration du béton, la déformation équivalente introduite est définie à partir des valeurs positives des déformations, soit:

\(\lbrace \begin{array}{c}\text{Si}{\epsilon}_{xx}^{e}\ge 0\text{alors}{\epsilon}_{\mathit{eq}}=|{\epsilon}_{xx}^{e}|\\ \text{Si}{\epsilon}_{xx}^{e}\le 0\text{alors}{\epsilon}_{\mathit{eq}}=\sqrt{2}\nu |{\epsilon}_{xx}^{e}|\end{array}\)

Remarque:

Dans le cas où \({\varepsilon}_{xx}^{e}\le 0\) , en \(\mathrm{1D}\) les déformations principales dans les autres directions sont \({\varepsilon}_{yy}^{e}={\varepsilon}_{zz}^{e}=-\nu {\varepsilon}_{xx}^{e}\) . En utilisant la formule \({\varepsilon}_{\text{eq}}=\sqrt{{\langle {\varepsilon}_{1}\rangle }_{+}^{2}+{\langle {\varepsilon}_{2}\rangle }_{+}^{2}+{\langle {\varepsilon}_{3}\rangle }_{+}^{2}}\) on obtient bien l’expression précédente.

La matrice tangente à pour expression: \(\frac{d{\sigma}_{xx}}{d{\varepsilon}_{xx}^{e}}=(1-\tilde{D})E-\frac{d\tilde{D}}{d{\varepsilon}_{xx}^{e}}E{\varepsilon}_{xx}^{e}\) avec:

\(\begin{array}{ccc}\text{Si}{\epsilon}_{xx}^{e}\ge 0\text{et}{\epsilon}_{\mathit{eq}}\ge {\epsilon}_{t0}& \frac{d\stackrel{~}{D}}{d{\epsilon}_{xx}^{e}}=\frac{d{D}_{t}}{d{\epsilon}_{xx}^{e}}=& \left(\frac{{\epsilon}_{t0}\left(1-{A}_{t}\right)}{{\epsilon}_{\mathit{eq}}^{2}}+\frac{{A}_{t}{B}_{t}}{\exp\left[{B}_{t}({\epsilon}_{\mathit{eq}}-{\epsilon}_{t0})\right]}\right)\\ \text{Si}{\epsilon}_{xx}^{e}<0\text{et}{\epsilon}_{\mathit{eq}}\ge {\epsilon}_{c0}& \frac{d\stackrel{~}{D}}{d{\epsilon}_{xx}^{e}}=\frac{d{D}_{c}}{d{\epsilon}_{xx}^{e}}=& -\sqrt{2}\nu \left(\frac{{\epsilon}_{c0}\left(1-{A}_{c}\right)}{{\epsilon}_{\mathit{eq}}^{2}}+\frac{{A}_{c}{B}_{c}}{\exp\left[{B}_{c}({\epsilon}_{\mathit{eq}}-{\epsilon}_{c0})\right]}\right)\end{array}\)

Les cas test [V6.02.120], [V6.02.119], [V5.02.130] mettent enœuvrela loi de comportement demazars dans sa version \(\mathrm{1D}\) .

Figure 8.1-a: Comportement de Mazars dans sa version \(\mathrm{1D}\) .

Variables internes

La loi de comportement est écrite en découplant les endommagements de traction et de compression, les deux endommagements ne sont plus des variables internes [R7.01.08].

Cette loi est dédiée aux calculs de génie civil. Pour faciliter les interprétations des résultats, deux variables sont créées pour décrire l’état « limite” du matériau béton, conformément à ce qui ce fait dans les règlements de calcul de béton armé aux états limites.

• La variable CRITSIG donne des informations par rapport à l’état de contrainte. Cette variable représente la contrainte divisée par la contrainte limite du béton donnée par l’utilisateur SIGM_LIM.

• La variable CRITEPS donne des informations par rapport à l’état de déformation. Cette variable représente la déformation équivalente \({\varepsilon}_{\mathit{eq}}\) divisée par la déformation limite donnée par l’utilisateur à l’aide du mot clef EPSI_LIM.

Les valeurs de la contrainte limite SIGM_LIM et de la déformation limite EPSI_LIM sont modifiables par l’utilisateur au moment de la définition du matériau: defi_materiau[U4.43.01], DEFI_MATE r_GC[U4.42.07].

L’écriture de la loi de mazars ne permet pas de calculer une dissipation intrinsèque au modèle. Mais, lors de calculs sismiques il peut être utile à l’utilisateur de connaître l’énergie dissipée non récupérable. La variable dissip représente le cumul d’énergie non récupérable. L’incrément d’énergie non récupérable s’écrit \(\Delta \mathit{Eg}=\frac{1}{2}(E(1-{D}^{+})\Delta \varepsilon –({\sigma}^{+}-{\sigma}^{-}))\Delta \varepsilon\) .

Matériau non-adoucissant

Matériau adoucissant

Les variables internes pour la loi de mazars en \(\mathrm{1D}\) :

\(\mathit{V1}\) CRITSIG : Critère en contrainte

\(\mathit{V2}\) CRITEPS : Critère en déformation.

\(\mathit{V3}\) ENDO : Endommagement [R7.01.08].

\(\mathit{V4}\) EPSEQT : Déformation équivalente de traction

\(\mathit{V5}\) EPSEQC : Déformation équivalente de compression

\(\mathit{V6}\) RSIGMA : Rapport de tri-axialité.

\(\mathit{V7}\) TEMP_MAX : Température maximale atteinte dans le matériau

\(\mathit{V8}\) DISSIP : Énergie non récupérable.

Loi de comportement RELAX_ACIER

Le phénomène de relaxation des aciers utilisés en précontrainte est réglementé. Les principaux règlements sont: BPEL83, NF-EN-1992-1-1 Octobre 2005, AFCEN-ETCC-2010,

On souhaite pouvoir modéliser des déformations qui vont varier lentement au court du temps, notamment pour la prise en compte du fluage du béton et des variations de température. On souhaite également prendre en compte l’influence de la température sur le phénomène de relaxation.

Réglementairement,il serait possible de tenir compte de l’effet du fluage du béton, de la déformation thermique en faisant une combinaison linéaire des différents phénomènes (Cf règlementspour plus de détails). Cette démarche est incompatible avec un calcul aux éléments finis.

Ces règlements permettent également de corriger la relaxation après un «saut» de déformation, pouvant être dû à une variation de la contrainte issue d’un chargement extérieur (par exemple lors d’un essai décennal sur une enceinte). En aucun cas ce «saut» ne peut évoluer dans le temps, ce ne peut être qu’un Dirac.La méthode utilisée dans ce cas est celle du temps équivalent qui, par construction, est incompatible avec un code aux éléments finis.

Formulation du modèle

Pour que la loi de relaxation soit utilisable dans un code aux éléments finis pour des calculs de structure avec des variations de chargements tels que: fluage du béton, reprise de tension des câbles, prise en compte de l’influence de la thermique,… elle doit être incrémentale et thermodynamiquement correcte.

Plusieurs formulations pour la modélisation de la déformation liée à la relaxation existent dans la littérature.

La formulation adoptée est basée sur celle proposée par J.Lemaitre:

\(\sigma =E.{\epsilon}^{e}\) \(\epsilon ={\epsilon}^{e}+{\epsilon}^{\mathit{an}}\)

\({{\dot{\epsilon}}^{\mathit{an}}=⟨\frac{\sigma –R\left({\epsilon}^{\mathit{an}}\right)}{{f}_{\mathit{prg}}.k}⟩}^{n}\) avec \(R\left({\epsilon}^{\mathit{an}}\right)=\frac{{f}_{\mathit{prg}}.c.{\epsilon}^{\mathit{an}}}{{\left(1+{\left(b.{\epsilon}^{\mathit{an}}\right)}^{\mathit{nr}}\right)}^{\frac{1}{\mathit{nr}}}}\) [éq9.1-1]

La loi de comportement est 1D, et uniquement disponible pour les modélisations de type barre, qui sont utilisées pour modéliser les câbles de précontrainte dans l’opérateur DEFI_CABLE_BP.

Pour tenir compte de l’influence de la température sur la relaxation, tous les coefficients de la loi peuvent être des fonctions de la température.

Remarque:

Les paramètres \(c\) , \(b\) , \(n\) et \(\mathit{nr}\) sont sans unité, et donc indépendant des unités utilisées pendant l’étude. \(k\) est adimensionné par rapport à \({f}_{\mathit{prg}}\) donc par rapport aux contraintes, mais pas par rapport au temps. En effet \({\dot{\epsilon}}^{\mathit{an}}\) est homogène à \({[s]}^{-1}\) , si l’unité de temps est en seconde. Donc si l’on connaît \(k\) pour une vitesse de déformation dans une unité de temps, il est nécessaire de convertir sa valeur par rapport à l’unité de temps utilisée lors de l’étude.

Variables internes

Deux variables internes:

• • \(V1\) : la déformation anélastique cumulée: \({\epsilon}^{\mathit{an}}\) .

• \(V2\) : mémorisation de la raideur tangente au comportement.

Intégration explicite

L’intégration se fait par une méthode de type Runge-Kutta d’ordre cinq à pas variables.

Identification des paramètres

Plusieurs solutions sont possibles:

• Un travail analytique permet de déterminer des relations entre les coefficients présents dans les formules réglementaires et celle de la loi de comportement: tangente au début de la relaxation, asymptote pour un temps «infini», …

• Simulation à l’aide des formules réglementaires de plusieurs «essais de relaxation». Pour un type de câble donné qui correspond à un jeu de paramètres (k1, k2, \({\rho}_{1000}\) , \({f}_{\mathit{prg}}\) ). Différentes valeur de :math:`mu ` vont permettre d’obtenir un faisceau «d’essais» sur lequel une identification des paramètres peut être réalisée.

• Étude bibliographique pour retrouver le travail expérimental et théorique réalisé sur la relaxation des câbles de précontrainte, afin d’utiliser directement les essais et d’identifier les paramètres.

Résultat d’une identification

Cette identification a été réalisée par rapport à la simulationd’un essai de relaxation obtenu avec les formules réglementaires.

On se place dans le cas \({\rho}_{1000}=2.5\text{\%}\) (\(k1=6.0E-03\) , \(k2=1.1\) ) avec \({f}_{\mathit{prg}}=1800\mathit{MPa}\) , \(E=190000\mathit{MPa}\) . La déformation initiale correspond à un taux de chargement de \(\mu =0.75\) . La simulation ainsi que l’identification sont réalisées sur [0h, 4000h].

Les figures ci-dessous comparent le résultat de l’identification avec la courbe réglementaire.

Figure 9.5-a: Évolution de la contrainte au cours du temps: référence et code_aster.

Figure 9.5-b: Zoom sur les 10 premières heures, de l’évolution de la contrainte au cours du temps.

Les observations vis-à-vis de l’identification:

• Très bonne correspondance entre la loi identifiée et la loi réglementaire.

• Très bonne identification de la pente à l’origine et de l’asymptote.

La loi proposée est donc apte à décrire correctement la relaxation d’un câble de précontrainte, si sa relaxation suit les courbes réglementaires.

Utilisation dans code_aster

La définition des paramètres matériaux se fait classiquement dans la commande DEFi_materiau, mot clef RELAX_ACIER, tous les paramètres peuvent dépendre de la température.

\({f}_{\mathit{prg}}\)	Contrainte à rupture du câble. Cette grandeur est facultative, car elle peut être également définie dans les matériaux BPEL_ACIERou ETCC_ACIERet dans ce cas \({f}_{\mathit{prg}}\) est une constante. La valeur/fonction donnée sous RELAX_ACIER_CABLest prioritaire.
ECOU_K ECOU_N	Correspond au coefficient k, dans l’équation . Correspond au coefficient n, dans l’équation .
ECRO_N ECRO_B ECRO_C	Correspond au coefficient nr, dans l’équation . Correspond au coefficient b, dans l’équation . Correspond au coefficient c, dans l’équation .
\({f}_{\mathit{prg}}\)	Contrainte à rupture du câble. Cette grandeur est facultative, car elle peut être également définie dans les matériaux BPEL_ACIER ou ETCC_ACIER et dans ce cas \({f}_{\mathit{prg}}\) est une constante. La valeur/fonction donnée sous RELAX_ACIER_CABL est prioritaire.
ECOU_K ECOU_N	Correspond au coefficient k, dans l’équation . Correspond au coefficient n, dans l’équation .
ECRO_N ECRO_B ECRO_C	Correspond au coefficient nr, dans l’équation . Correspond au coefficient b, dans l’équation . Correspond au coefficient c, dans l’équation .

Méthode pour utiliser en 1D tous les comportements 3D

Comme pour le traitement des contraintes planes [R5.03.03], il est possible de bénéficier pour les modélisations 1Ddes comportements disponibles en 3D. On étend pour cela la méthode due à R.deBorst au cas 1D, en traitant cette condition (champ de contraintes unidimensionnel) non pas au niveau de la loi de comportement mais au niveau de l’équilibre. On obtient ainsi au cours des itérations de l’algorithme de STAT_NON_LINEdes champs de contraintes qui tendent vers un champ unidirectionnel. On vérifie, à convergence des itérations de Newton globales, que les champs de contraintes sont effectivement unidirectionnels, à une précision près, sinon on continue les itérations. La méthode consiste à décomposer les champs de déformations et de contraintes en une partie purement unidirectionnelle (direction x) et une partie relative aux autres directions, et d’effectuer une condensation statique en écrivant que les composantes des contraintes relatives aux autres directions sont nulles. On ne considère dans les tenseurs (d’ordre 2) que les termes diagonaux, écrits sous forme de vecteurs à 3 composantes. La direction x correspond à la direction de l’élément (barre, poutre multifibre) ou à la direction des armatures de grille. A un instant quelconque de la résolution du comportement incrémental, l’opérateur tangent \(D\) relie l’accroissement de contraintes à l’accroissement de déformation par:

\(d\sigma =\left[\frac{\partial \sigma }{\partial \varepsilon }\right]d\varepsilon =Dd\varepsilon\) que l’on réécrit \(\left[\begin{array}{}d{\sigma}_{x}\\ d{\sigma}_{y}\\ d{\sigma}_{z}\end{array}\right]=\left[\begin{array}{ccc}{D}_{11}& {D}_{12}& {D}_{13}\\ {D}_{21}& {D}_{22}& {D}_{23}\\ {D}_{31}& {D}_{32}& {D}_{33}\end{array}\right]\left[\begin{array}{}d{\varepsilon}_{x}\\ d{\varepsilon}_{y}\\ d{\varepsilon}_{z}\end{array}\right]\) .

En écrivant ces accroissements comme la différence entre les itérations \(n\) et \(n+1\) de Newton, on obtient:

\(d\sigma ={\sigma}^{n+1}-{\sigma}^{n}=\Delta {\sigma}^{n+1}-\Delta {\sigma}^{n}\) , \(d\varepsilon =\Delta {\varepsilon}^{n+1}-\Delta {\varepsilon}^{n}\)

A convergence, cet écart doit tendre vers zéro.

En introduisant les conditions \({\sigma}_{y}^{n+1}=0\) et \({\sigma}_{z}^{n+1}\) (comportement unidirectionnel), on obtient, pour l’itération \(n+1\) :

\(\left[\begin{array}{}d{\sigma}_{x}\\ d{\sigma}_{y}\\ d{\sigma}_{z}\end{array}\right]=(\begin{array}{}{\sigma}_{x}^{n+1}-{\sigma}_{x}^{n}\\ {\sigma}_{y}^{n+1}-{\sigma}_{y}^{n}\\ {\sigma}_{z}^{n+1}-{\sigma}_{z}^{n}\end{array})=(\begin{array}{}{\sigma}_{x}^{n+1}-{\sigma}_{x}^{n}\\ -{\sigma}_{y}^{n}\\ -{\sigma}_{z}^{n}\end{array})=\left[\begin{array}{ccc}{D}_{11}& {D}_{12}& {D}_{13}\\ {D}_{21}& {D}_{22}& {D}_{23}\\ {D}_{31}& {D}_{32}& {D}_{33}\end{array}\right]\left[\begin{array}{}d{\varepsilon}_{x}\\ d{\varepsilon}_{y}\\ d{\varepsilon}_{z}\end{array}\right]\)

Les deux dernières équations permettent d’exprimer \(d{\varepsilon}_{y}\) et \(d{\varepsilon}_{z}\) en fonction de \(d{\varepsilon}_{x}\) :

\(\lbrace \begin{array}{c}d{\varepsilon}_{y}=\frac{1}{{D}_{22}}(-{\sigma}_{y}^{n}-{D}_{21}d{\varepsilon}_{x}-{D}_{23}d{\varepsilon}_{z})\\ d{\varepsilon}_{z}=\frac{1}{{D}_{33}}(-{\sigma}_{z}^{n}-{D}_{31}d{\varepsilon}_{x}-{D}_{32}d{\varepsilon}_{y})\end{array}\)

soit

\(\begin{array}{}d{\varepsilon}_{y}=\frac{1}{\Delta}(-{D}_{33}{\sigma}_{y}^{n}+{D}_{23}{\sigma}_{z}^{n}+{D}_{y}d{\varepsilon}_{x})\\ d{\varepsilon}_{z}=\frac{1}{\Delta}(-{D}_{32}{\sigma}_{y}^{n}+{D}_{22}{\sigma}_{z}^{n}+{D}_{z}d{\varepsilon}_{x})\end{array}\)

avec \(\Delta ={D}_{33}{D}_{22}-{D}_{23}{D}_{32,}{D}_{y}={D}_{23}{D}_{31}-{D}_{21}{D}_{33,}{D}_{z}={D}_{32}{D}_{21}-{D}_{31}{D}_{22}\)

en reportant ces expressions dans la première équation, on obtient:

\({\sigma}_{x}^{n+1}={\sigma}_{x}^{n}+({D}_{11}+\frac{{D}_{12}{D}_{y}+{D}_{13}{D}_{z}}{\Delta})d{\varepsilon}_{x}+\frac{{D}_{12}{D}_{23}-{D}_{22}{D}_{13}}{\Delta}{\sigma}_{z}^{n}+\frac{{D}_{12}{D}_{32}-{D}_{12}{D}_{33}}{\Delta}{\sigma}_{y}^{n}\)

L’équilibre à l’itération \(n+1\) s’écrit:

\(\begin{array}{}\int{D}^{T}{\sigma}^{n+1}\mathrm{dv}=\int{B}^{T}{\sigma}_{x}^{n+1}\mathrm{dv}=\int{B}^{T}({D}_{11}+\frac{{D}_{12}{D}_{y}+{D}_{13}{D}_{z}}{\Delta})d{\varepsilon}_{x}\\ +\int{B}^{T}({\sigma}_{x}^{n}+\frac{{D}_{12}{D}_{23}-{D}_{22}{D}_{13}}{\Delta}{\sigma}_{z}^{n}+\frac{{D}_{12}{D}_{32}-{D}_{12}{D}_{33}}{\Delta}{\sigma}_{y}^{n})\mathrm{dv}\\ ={K}^{n}{\mathrm{du}}^{n+1}+\int{B}^{T}({\sigma}_{x}^{n}+\frac{{D}_{12}{D}_{23}-{D}_{22}{D}_{13}}{\Delta}{\sigma}_{z}^{n}+\frac{{D}_{12}{D}_{32}-{D}_{12}{D}_{33}}{\Delta}{\sigma}_{y}^{n})\mathrm{dv}\end{array}\)

On constate que la prise en compte du comportement unidimensionnel intervient à deux niveaux:

• dans la matrice tangente, par le terme correctif:

\(\int{B}^{T}\frac{{D}_{12}{D}_{y}+{D}_{13}{D}_{z}}{\Delta}B\mathrm{dv}\)

• dans l’écriture du second membre, par le terme correctif:

\(\frac{\int{B}^{T}}{\Delta}(({D}_{12}{D}_{23}-{D}_{22}{D}_{13}){\sigma}_{z}^{n}+({D}_{12}{D}_{32}-{D}_{12}{D}_{33}){\sigma}_{y}^{n})\mathrm{dv}\)

Pour mettre en œuvre cette méthode, il suffit de calculer ces termes correctifs et de les ajouter aux contraintes et matrice tangente obtenue de la résolution 3D du comportement. Pour cela il est nécessaire de stocker des informations d’une itération de Newton à l’autre, par le biais de 4 variables internes supplémentaires. Les étapes de la résolution sont:

• à l’itération \(n+1\) , les données sont: \(\Delta {u}^{n+1},{\sigma}^{-},{\alpha}^{-}\) et les 4 variables internes (calculées à l’itération \(n\) ):

\(\begin{array}{}\mathrm{V1}=\Delta {\varepsilon}_{y}^{n}+\frac{1}{\Delta}({D}_{23}{\sigma}_{z}^{n}-{D}_{33}{\sigma}_{y}^{n}-{D}_{y}\Delta {\varepsilon}_{x}^{n}),\mathrm{V2}=\frac{{D}_{y}}{\Delta}\\ \\ \mathrm{V3}=\Delta {\varepsilon}_{z}^{n}+\frac{1}{\Delta}({D}_{32}{\sigma}_{z}^{n}-{D}_{22}{\sigma}_{z}^{n}-{D}_{z}\Delta {\varepsilon}_{x}^{n}),\mathrm{V4}=\frac{{D}_{z}}{\Delta}\end{array}\)

• avant d’effectuer l’intégration du comportement (effectué en axisymétrique) on calcule:

\(\begin{array}{}\Delta {\varepsilon}_{y}^{n+1}=\Delta {\varepsilon}_{y}^{n}+\frac{1}{\Delta}(-{D}_{33}{\sigma}_{y}^{n}+{D}_{23}{\sigma}_{z}^{n}+{D}_{y}d{\varepsilon}_{x})\\ \Delta {\varepsilon}_{z}^{n+1}=\Delta {\varepsilon}_{z}^{n}+\frac{1}{\Delta}(-{D}_{32}{\sigma}_{y}^{n}+{D}_{22}{\sigma}_{z}^{n}+{D}_{z}d{\varepsilon}_{x})\end{array}\)

• l’intégration du comportement fournit les contraintes \({\sigma}^{n+1}\) et l’opérateur tangent \(D\) ,

• on modifie le second membre et la matrice tangente comme indiqué ci-dessus,

• on stocke les nouvelles variables internes et on vérifie si \(∣{\sigma}_{z}^{n+1}∣<\eta\) et \(∣{\sigma}_{y}^{n+1}∣<\eta\) , avec \(\eta =\xi ∣{\sigma}_{x}^{n+1}∣\) , \(\xi =\) RESI_INTE_RELA

Remarque :

Les quatre variables internes supplémentaires sont ajoutées après les variables internes de la loi de comportement.

Bibliographie

[1] J.GUEDES, P.PEGON, P.E.PINTO: ’’A Fibre/Timoshenko Beam Element in Castem 2000’’, Joint Research Centre, European Commission, Institute for Safety Technology,1994.

[4] G.MONTI, C.NUTI: ’’Non-linear Cyclic Behaviour of Reinforcing Bars Including Buckling’’, Journal of Structural Engineering, Vol. 118, No 12, December 1992.

[5] P.DEBONNIERES, M.ZIDI: «Introduction de la viscoplasticité dans le modules de thermomécanique de CYRANO3: principe, description et validation» Note HI-71/8334.

[6] J.M.PROIX, B.QUINNEZ, P.MASSIN, P.LACLERGUE: «Assemblages combustibles sous irradiation. Étude de faisabilité». Note HI-75/97/017/0.

[7] I. LE PICHON, P. GEYER: «Modélisation du comportement viscoplastique anisotrope des tubes de gainage des crayons combustibles» Note HT-B2/95/018/A

[8] L. RANCŒUR: «Comportement en fluage axial sous irradiation des tubes guides en M5 et Zircaloy-4 recristallisés» note MMC HT25-C2004-192/LRC du 25/10/2004.